雙層三階多項式擬合的彈道導彈發射點估計方法*

邱明劼,周垂紅,汪圣利

(1.南京電子技術研究所,南京 210039;2.中國電子科學研究院,北京 100041)

0 引 言

傳統炮位偵校雷達的作戰對象以迫擊炮、榴彈炮等無動力炮彈為主,此類目標運動軌跡符合拋物線彈道規律,發射點較好預測。但隨著作戰縱深的日益增加,中近程戰術彈道導彈逐步成為炮位偵校雷達的新型作戰對象。

與無動力的炮彈不同,彈道導彈的飛行過程可分為主動段和被動段兩個部分[1]。由于彈道導彈主動段和被動段運動存在較大差異,為了實現對發射點的精確估計,炮位偵校雷達需針對主動段軌跡而非被動段軌跡進行建?；驍M合,從而外推得到發射點位置。然而,由于導彈在主動段的受力情況十分復雜,對于非合作目標,在無先驗信息情況下難以準確推斷導彈的運動模型,僅能根據理論方法構建近似模型[2-3]。目前對主動段彈道近似建模的方法可以分為兩類:一是利用主動段彈道輪廓先驗知識進行建模[4]。該類方法通常需要預先積累大量先驗數據,難度大且通常只能對已知型號導彈起到較好效果,對未知型號的估計效果差。二是利用通用數學模型對主動段彈道進行擬合[5-7]。由于該類方法不需要積累先驗知識,具有可行性強、適用性廣的優勢特點,因而是目前研究的主要方向。

現有的傳統彈道外推算法大多使用各種濾波(如粒子濾波)或平滑算法(如最小二乘)處理雷達測量的位置和速度數據后,直接代入質點彈道方程計算,但在長距離情況下外推誤差較大,難以適應彈道導彈的遠距離外推需求[8]。文獻[9]和[10]分別采用了簡單線性模型和重力轉彎模型進行彈道擬合,但均無法準確反映主動段運動特性,導致估計精度大幅下降。文獻[11]提出了基于時間逆轉的反向卡爾曼濾波方法,對于傳統炮彈發射點估計效果較好,但仍無法適應彈道導彈主動段運動特征。

本文提取彈道導彈主動段在垂向和射向兩個方向的運動特征,將發射點零時刻估計作為發射點估計的中間過程,提出了一種基于雙層擬合的發射點估計方法。該方法首先利用發射時間和導彈垂直位移之間的關系,利用三次多項式迭代估計導彈發射零時,然后利用時間與射向之間的物理關系,估計出導彈發射點位置。本文所提方法估計參數數量少,計算簡單,定位誤差較傳統方法降低了50%以上,具有重要軍事應用價值。

1 彈道導彈運動軌跡仿真

通過正向設計,精確仿真彈道導彈軌跡是進行彈道運動軌跡研究和發射點估計的重要前提。表1給出了彈道導彈的飛行階段及其受力情況,以此為理論基礎進行軌跡仿真[12]。

表1 彈道模型的分段及其受力情況

假設彈道導彈關機點以后僅受到重力的作用,導彈目標的位置為p=[px,py,pz]T,速度為v=[vx,vy,vz]T,則依據圓球體地球模型建立的彈道導彈運動軌跡可以表示為

(1)

式中:aT為推進力產生的加速度;aD為空氣阻力產生的加速度;aG為地心引力產生的加速度;aC為科氏力和離心力產生的加速度。aT,aD,aG,aC構成了彈道目標的加速度[13]。

(2)

(3)

設ρ(h)是空氣密度函數,ρ0=1.22 kg/m3,k=0.141 41×10-3m-1,h(t)為t時刻目標海拔高度,v(t)為t時刻目標速度,uv是單位向量,m(t)為t時刻目標質量,cD(v)是阻力系數,S是目標體與其速度方向的截面積,則空氣阻力加速度可近似表示為

(4)

ρ(h(t))=ρ0e-kh(t)

(5)

(6)

設ω=[0,0,ω]′是地球自轉角速度矢量,則科氏力和離心力產生的加速度為

(7)

綜上所述,根據式(1)～(7),可得到彈道導彈主動段運動模型為

(2ωvy+ω2px)

(-2ωvx+ω2py)

(8)

在被動段沒有了推力,運動模型變為

(2ωvy+ω2px)

(-2ωvx+ω2py)

(9)

根據式(8)和式(9)即可通過正向設計的方式構建出彈道導彈的仿真運動軌跡,用于對本文所提方法進行仿真驗證。

2 主動段雙層擬合方法

圖1為主動段雙層擬合方法的流程圖,主要包括擬合坐標系建立、發射零時刻估計、主動段水平位移擬合和發射點估計及誤差計算4個步驟。

圖1 主動段雙層擬合方法流程

2.1 擬合坐標系建立

彈道曲線是在發射坐標系O′X′Y′Z′下進行描述的,但由于雷達無法提前知悉非合作目標的發射點位置,僅能通過探測數據擬合發射坐標系,從而在這個坐標系下進行發射點估計。在本文中將雷達依據觀測數據建立的坐標系稱為擬合坐標系。

擬合坐標系是以雷達觀測到的第一個導彈位置點在水平面上的投影作為原點O,以雷達觀測到的主動段數據在水平面上的投影的最小二乘擬合結果作為X軸,垂直于水平面豎直向上為Y軸,如圖2中虛線坐標系OXYZ所示。

圖2 擬合坐標系OXYZ、發射坐標系OX′Y′Z′和臨時坐標系OtXtYtZt示意圖

為了建立擬合坐標系,需要首先依據雷達的部署位置,首先建立臨時坐標系OtXtYtZt,然后再通過坐標旋轉的方法,將其旋轉至擬合坐標系下。

下面具體闡述擬合坐標系建立的過程:

1) 以雷達觀測到的第一個彈道導彈軌跡點在地面上的投影為臨時坐標系原點Ot,其與雷達位置之間的連線為臨時坐標系Xt軸,垂直地面向上為臨時坐標系Yt軸,并根據右手法則確定臨時坐標系Zt軸,則觀測到的臨時坐標系下的軌跡點集合為Pt={[0,0,zt0]T,[xt1,yt1,zt1]T,[xt2,yt2,zt2]T…}。

2)對軌跡點集合Pt在XtOtYt平面上的投影進行一階函數最小二乘擬合,求得各軌跡點Yt軸坐標關于Xt軸坐標的關系,即ytn=f(xtn)=axtn+b,其中a和b是最小二乘擬合得到的擬合參數。

3)根據求得的參數a,得到臨時坐標系OtXtYtZt與擬合坐標系OXYZ之間的角度偏差φ=arctana,然后使用旋轉矩陣將軌跡點坐標進行轉換,即

(10)

得到在擬合坐標系下的軌跡點集合P={[x0,y0,z0]T,[x1,y1,z1]T,[x2,y2,z2]T…}。

通過上述步驟,即可將坐標由雷達初始探測得到的臨時坐標系轉換至擬合坐標系中。在后續的兩次三階多項式擬合過程中,均認為彈道僅在擬合坐標系的XOY平面內運動,在可以將導彈在三維空間中的位移情況簡化為在彈道平面內的二維運動,簡化彈道估計的同時完成了位移分量的分解,得到在X方向的水平位移和Y方向的垂直位移兩個分量。

但從圖2中可以看到,擬合坐標系Y軸與導彈發射坐標系的Y′軸方向相同,但X軸和Z軸均與實際發射坐標系的X′軸和Z′軸存在一個角度偏差θ,這是由于實際彈道主動段在Z'方向的運動以及雷達的觀測誤差影響所導致。在仿真條件下,可以通過對無誤差的仿真彈道在X′O′Z′平面投影的坐標進行一次多項式擬合,從而與擬合坐標系下的軌跡點集合擬合結果進行比較,根據兩者擬合后一次函數的斜率求得角度偏差θ。但在實際作戰條件下,將無法獲取這個角度,這將使得沿著擬合坐標系X軸進行外推的方向與實際彈道運動方向之間存在偏差,而這正是發射點估計的誤差來源之一。

2.2 發射零時刻估計

在2.1節建立的擬合坐標系中,當忽略彈道軌跡在Z軸方向上的運動時,可認為zTn≡0(n=1,2,3…),則導彈絕對發射時刻TS的彈道點坐標可表示為(xTS,yTS,zTS)=(xTS,0,0),Ti時刻的彈道點坐標為(xTi,yTi,0)。通過彈道導彈運動模型仿真結果和大量實測數據可以看出,導彈主動段在垂直方向的位移量yTi和水平方向的位移量xTi-xTS均為隨相對發射時間ti=Ti-TS(i=1,2,3…)單調增長的平滑曲線。

盡管TS和xTi-xTS均為未知量,但絕對發射時間Ti和在Y軸方向上的位移量yTi均直接已知,因此可以假設TS=0 s,則ti=Ti,通過三階函數最小二乘擬合的方法,擬合出時間ti與yti的函數關系yti=g(ti)?？紤]到發射的導彈在相對發射時間0 s以前,其高度方向位移應始終等于0,而在相對發射時間0 s之后,高度方向位移應當始終大于0,因此在進行函數擬合時,基于上述事實對擬合過程進行約束:當ti≤0時,主動段Y方向位移均等于0;當ti>0時,主動段Y方向位移均大于等于0,即

(11)

從而使得擬合得到的三次函數曲線系數快速收斂,且不會因為過擬合導致外推趨勢發散。

根據上述兩個事實可知,當擬合出的三次函數曲線最靠近初始觀測值的零點坐標越接近于0時,即t0越接近于0 s,所擬合出的函數應當越符合事實情況。因此,在擬合出時間-垂向位移函數曲線后,首先使用梯度下降法求出離第一個觀測點最近的擬合曲線零點坐標(tz,0),然后將發射零時刻加上tz與優化速率η乘積代表的時間修正量Δt,對TS進行修正。該過程也可等效為將所有觀測數據時間Ti減去由tz與優化速率η乘積代表的時間修正量Δt,對Ti進行修正。

(12)

之后重復利用新的相對時間ti與位移yti進行擬合、求零點、修正,使得逐次求得的三次函數零點tz逐步與0 s接近。當擬合曲線零點坐標tz小于所設定的閾值ttreshold,即可認為在滿足上述約束條件下,時間-垂向位移擬合曲線與真實彈道在Y方向的曲線吻合,從而確定了每個軌跡點的相對時刻ti。

發射零時刻估計的具體過程如下:

1) 假定TS=0,利用每個軌跡點的絕對時刻Ti和高度方向位移yti進行三次函數的最小二乘擬合,得到擬合函數

2) 使用梯度下降法,從第一個觀測到的軌跡點開始,求出離其最近的擬合曲線零點坐標(tz,0)。

3)若|tz|≤ttreshold,即認為此時的TS已經足夠接近真實時間,擬合得到的曲線與實際情況吻合,得到擬合結果g(ti);|tz|>ttreshold,則使用優化速率η乘以tz得到時間修正量Δt。

4) 將TS加上時間修正量Δt,重新求得每個軌跡點的相對時間ti=Ti-TS。

5)重復步驟1)～3),直至得到擬合結果g(ti)。

2.3 主動段水平位移擬合

在擬合坐標系中,發射點位置估計實際上已經被轉化為求取相對時間t0時刻水平X方向位移量的問題。在2.2節中已經將TS收斂至真實絕對發射時間,從而使得發射點的相對發射時間t0≈0 s,進而進行第二次三次函數的最小二乘擬合。此時將各軌跡點的相對發射時間ti作為自變量,以軌跡點的X方向位移作為應變量進行最小二乘擬合,即可得到導彈在水平面投影上的時間-水平位移擬合結果xTi=h(ti)。

在進行三次函數擬合之前,同樣需要考慮如下的約束條件。由于是采用三次多項式對X方向的位移進行擬合,而對位移求導即可得到對應方向上的運動速度,因此在X方向上的速度應為二次多項式形式?？紤]彈道導彈通常采用垂直發射,在X方向上初始速度為0,因此首先對時間-速度曲線構建無常數項的二次多項式擬合模型:

vti=k(ti)=avti2+bvti

(13)

接著再對上式進行積分即可得到X方向位移關于時間Ti的三次多項式模型:

xti=p(ti)=axti3+bxti2+cx

(14)

利用該三次多項式模型對數據進行最小二乘擬合,即可得到主動段水平位移關于發射相對時間的具體表達式。由于這次擬合沒有相關的事實約束條件,因此其擬合精度與t0≈0 s的準確性密切相關。

2.4 發射點估計及誤差計算

由于擬合坐標系中的XOY平面是對三維發射坐標系中X′O′Y′平面的近似估計,兩者之間存在角度偏差θ,為了計算發射點估計結果的誤差,需要將在發射坐標系中的真實發射點位置坐標根據角度θ變換至擬合坐標系下。

根據2.1節所述,利用一次多項式擬合求得擬合坐標系和真實發射坐標系之間的夾角θ后,設發射點真實位置在發射坐標系中的坐標為(xT′,yT′,zT′),雷達觀測首點在發射坐標系中的坐標為(xF′,yF′,zF′),則發射點真實位置在擬合坐標系中的坐標(xS,yS,zS)可表示為

(15)

從而得到最終的發射點定位誤差ε為

(16)

當雷達與導彈發射點之間距離為R時,定位精度P定義為定位誤差除以偵察距離:

(17)

3 仿真實驗及結果分析

3.1 仿真條件

仿真過程采用搭載主頻為3.00 GHz的i5處理器的計算機進行,仿真語言為Matlab語言。

仿真時,首先根據式(8)和式(9)表示的彈道導彈主動段和被動段運動模型,調整參數生成射程200～1 000 km的近程戰術彈道導彈軌跡,然后設置雷達位置并生成雷達觀測數據。結合實際作戰場景情況,假設炮位偵察雷達部署于導彈射向方向,且距離發射點位置300 km,雷達探測角度精度0.1°、距離精度100 m,探測仰角范圍為2°～25°,仿真場景如圖3所示。

圖3 仿真實驗場景設置圖

考慮地球曲率影響和雷達仰角探測范圍,從仿真生成的彈道中選擇出可被雷達觀測的彈道弧段,并添加雷達測量誤差。接著分別使用常規彈道外推的定軌方法[15],和本文提出的主動段雙層擬合方法進行發射點位置估計,并多次添加雷達誤差進行蒙特卡羅仿真,統計發射點精度結果。單次仿真流程如圖4所示。

圖4 單次仿真實驗流程

本文所提方法中提到的時間修正量閾值參數Ttreshold和優化速率參數η,在所有的仿真過程中均取2和0.5,發射時間TS初始值設置為-100 s,則相對時間t0=100 s。

3.2 仿真分析

根據彈道導彈運動模型生成的某條彈道軌跡如圖5所示,并以此彈道軌跡為例進行精度分析。該彈道在發射坐標系下各方向的位移與時間的關系曲線如圖6所示,其中紅色曲線部分即為彈道導彈主動段軌跡,黑色曲線部分為被動段軌跡。

圖5 某條仿真彈道軌跡

圖6 某條仿真彈道時間-位移分量曲線

選取添加測量誤差后的主動段數據,利用主動段雙層擬合方法則雷達可觀測到的彈道軌跡如圖7所示。

圖7 雷達可觀測彈道軌跡

(a)時間-位移擬合第1次迭代結果

將第一層三階擬合后得到的Ti用于和擬合坐標系中的X方向水平位移進行第二層三階擬合,得到擬合結果如圖9中綠色曲線所示。

圖9 主動段時間-水平位移三階函數擬合結果

從圖9中可以看出,在構建的擬合坐標系下,發射點位置估計結果為(-345.13,0,0)m。通過坐標變換,真實發射點位置轉換至擬合坐標系下的坐標值為(-70.13,0,-40.51)m,根據式(15)可得定位誤差為277.97 m,并根據式(16)計算得到本次仿真定位精度為277.97/300000=0.09%。

按照上述實驗方法,共仿真生成10條彈道,其中前7條彈道為單級助推導彈彈道,后3條為二級助推導彈彈道。每條彈道采用傳統彈道導彈定軌方法和本文方法分別進行100次蒙特卡羅實驗,并使用2.4節所提的精度計算方法進行統計,得到發射點定位誤差和精度情況如表2所示。從表2的結果可以看出,在仿真條件下,傳統利用定軌方法對具有一級和二級助推的仿真彈道發射點定位精度在0.7%～1.4%之間,而本文所提方法對發射點的估計精度可達到0.6%以內,定位誤差相比優于傳統方法降低50%以上。

表2 傳統彈道外推方法與本文方法定位誤差對比

接著進一步考慮不同雷達探測精度對本文算法的影響。對于現有雷達,角度精度通常在0.5°以內,且對于遠距離探測時角度精度帶來的測量誤差遠大于距離精度影響,因此僅分析雷達角度精度對定位精度的影響。

分別選取具有一級助推和二級助推的彈道1和彈道8作為雷達測量精度對定位精度影響分析的實驗彈道。假設雷達測角精度不超過0.5°,每0.01°設置一個精度采樣間隔,生成雷達測量誤差與理論彈道軌跡疊加,并使用本文方法在每個測角精度條件下進行100次蒙特卡羅實驗,統計對兩條彈道的發射點定位精度,結果如圖10所示。

從圖10可以看出,對于具有一級助推和二級助推的彈道導彈軌跡,其發射點定位精度隨著雷達精度變差而逐步變差,但即使在雷達測角精度達到0.5°的情況下,使用本文算法得到的定位誤差也僅與測角精度0.1°情況下使用傳統方法得到的定位精度相近,從而充分說明本文所提方法可有效指導我方部隊遂行精確對地打擊任務,具有明顯先進性。同時,本文算法在3.1節所述的仿真條件下,單次仿真運行時間約為0.03 s,滿足發射點解算的實時性要求,具備良好的工程可實現性。

4 結束語

本文提出了一種基于主動段雙層擬合的彈道導彈發射點估計方法,利用三階函數擬合彈道導彈軌跡的時間-高度曲線,較準確地估計出相對發射時刻后,再利用相對發射時間與水平距離位移進行第二次三階函數擬合,進而得到發射點的位置估計結果。在仿真實驗條件下,將本文方法與傳統彈道外推方法進行了對比,充分驗證了方法的有效性。

本文所提方法能夠良好處理具有一級和二級助推的彈道導彈發射點估計任務,但當主動段的運動狀態更加復雜時,其軌跡將難以使用簡單多項式函數進行擬合逼近,從而導致較大預測誤差。隨著機器學習和神經網絡的廣泛應用,其在軌跡預測方面已有許多學者開展研究。由于神經網絡具有對復雜非線性系統的良好逼近能力,將神經網絡的方法引入到彈道導彈發射點估計中提升定位精度可作為后續進一步研究的方向。