王 彪,吳 楠,孟凡坤,王建濤

(1.信息工程大學,河南鄭州 450000;2.中國人民解放軍95129部隊,河南開封 475000)

彈道導彈射程遠、威力大,已經成為大國戰(zhàn)略威懾的重要武器。為應對這種威脅,近年來，用于導彈防御的導彈預警系統(tǒng)得到了很大發(fā)展。由于彈道導彈在自由段的飛行時間和距離占全彈道的80%以上[1],因此，自由段成為導彈預警研究的主要過程。導彈在自由段速度快、距離遠,而單部雷達探測范圍有限,往往需要多部雷達進行預警跟蹤。在多雷達預警中,當前一部雷達與后一部雷達對導彈的探測范圍不重疊時,需要進行雷達交接時刻的預報[2-3],來提高整體預警能力,使用有效的彈道預報方法可以更好地提高雷達交接時刻的預報性能[4-5]。

傳統(tǒng)的彈道預報方法包括以求解軌道根數(shù)為核心的解析幾何法和以求解目標動力學微分方程為核心的數(shù)值積分法[6]。文獻[6]分析了兩種方法的預報特性。解析幾何法將導彈視為二體運動,未考慮運動中的攝動因素,計算簡單，預報快捷,但是預報精度不高。數(shù)值積分法考慮了導彈運動中的攝動因素,建立了復雜的動力學微分方程,預報精度高，但計算量大。

恩克法由德國科學家恩克(Johann F. Encke)提出,曾應用在計算短周期彗星和小行星軌道[7]。本文改進恩克法，并將其應用在雷達交接時刻預報中,首先，介紹了恩克法彈道預報算法的具體步驟,其次，對雷達交接時刻計算流程進行了詳細分析,并基于無跡變換的方法進行了誤差分析[8]。最后，對算法進行了仿真實驗,并與傳統(tǒng)方法進行了對比,驗證了算法的有效性。

1 雷達交接時刻預報算法

1.1 雷達探測模型

從雷達站觀測導彈飛行的觀測數(shù)據中獲得導彈的狀態(tài)值,首先，需要建立雷達的觀測模型。雷達站的位置通常使用經度、緯度、高度(L,B,H)來描述,而導彈目標的位置[X,Y,Z]T通常在地固系(ECEF)下進行描述,雷達對導彈目標位置的探測通常使用斜距、方位角、仰角(R,A,E)來確定,因此，需要進行坐標系間的轉換,其轉換過程如下。

1)將雷達站位置(L,B,H)轉換為地固系坐標

(1)

2)將地固系下導彈對雷達的相對位置轉換為雷達直角坐標系ENU下的坐標

(2)

(3)

3)將導彈在雷達直角坐標系的坐標轉化為雷達探測數(shù)據(R,A,E)

(4)

1.2 恩克法彈道預報方法

雷達交接時刻預報的重要環(huán)節(jié)是如何對導彈進行彈道預報。本文使用的恩克法彈道預報主要分為三個方面:利用二體力學模型建立導彈運動的基準軌道，求解基準軌道與實際軌道的偏差值，導彈狀態(tài)在不同坐標系間的轉換,具體流程如下。

1)求解基準軌道根數(shù)與t時刻的位置ρ

基準軌道是理想的二體軌道,可由6個獨立的軌道根數(shù)(半長軸a、偏心率e、軌道傾角i、升交點赤經Ω、近地點幅角ω、真近點角f)來確定,軌道根數(shù)的求解則可由初始時刻的r0、v0算得[7]。

軌道確定后,軌道的偏近點角E,平近點角M,平均角速度n及過近地點時刻τ的計算如下

(5)

(6)

2)求解t時刻真實軌道與基準軌道的偏差量δr

本文分析導彈的受力情況,可得實際軌道和基準軌道的動力學微分方程分別為

(7)

(8)

式中,ap為攝動力加速度表達式。導彈在自由段受到的主要攝動力為地球非球形引力,該攝動力中J2攝動影響最大,為 10-3量級,而J4項攝動與其他攝動力均為 10-6量級以下,因此通常只考慮J2攝動。

由基準軌道的定義可知δr=r-ρ,可得

(9)

(10)

因式(10)為非線性微分方程,無法獲得解析解,故采用工程上常用的四階龍格-庫塔(R-K)方法計算其數(shù)值解。

3)計算t時刻導彈的狀態(tài)量r

在獲得t時刻導彈基準軌道的位置、速度及其偏差量后,t時刻真實軌道的位置、速度可由兩者之和得到,具體表達式分別為：

(11)

(12)

由于r(t)、v(t)是在地慣系下表示的,需要將其再統(tǒng)一到地固系下

(13)

(14)

以上是恩克法預報的三個步驟,在步驟2)中，恩克法與傳統(tǒng)的數(shù)值積分法一樣，需要進行微分方程的求解,不同之處為：數(shù)值積分法積分求解的對象是導彈狀態(tài)矢量,需要進行小時長積分多次迭代計算才能有較高精度;而恩克法積分求解的對象是實際軌道與基準軌道的攝動偏差量,是一個微小的量,在大時長積分的情況下,對整體位置的精度影響不大,從而可以進行單步積分,減小了整體的運算量節(jié)約了運算時間。恩克法彈道預報流程如圖1所示。

圖1 恩克法彈道預報流程圖

1.3 雷達交接時刻預報流程

根據建立的雷達探測模型,使用彈道預報方法,就可以進行雷達交接時刻的預報。假設前置雷達穩(wěn)定探測到導彈某一時刻的狀態(tài)矢量(位置和速度),已知后置雷達的位置及威力空間γmax=(γmax,Amax,Emax),則雷達交接時刻的預報流程如下:

1)確定彈道預報的時間范圍[t0,tf],其中t0為初始時刻,tf為導彈落點預報時刻,落點預報時刻可通過求解導彈的二體運動方程與地球理想球面方程計算獲得;

2)將時間范圍[t0,tf]進行2n(初始時令n=1)等分,得到時刻點t1，t2，…，t2n+1,使用彈道預報方法分別獲得對應時間點的位置r1，r2，…，r2n+1,進而，可計算出對應時刻點下，導彈與后置雷達的相對位置(γ1γ2…γ2n+1);

3)設t時刻雷達對導彈的觀測量γ=(γ,A,E),規(guī)定:當且僅當γ<γmax,Aγmax;若存在ti時刻滿足γi<γmax,且ti-1時刻滿足γi-1>γmax,說明預報交接點在區(qū)間[ti-1,ti]內(記為[tA,tB]),轉入步驟4)；若不存在ti時刻滿足上述條件,令n=n+1,并轉到步驟2);

4)使用彈道預報方法獲得(tA+tB)/2時刻導彈的位置r,并轉換為雷達與導彈的斜距γ,如果γ<γmax,則說明預報交接點在區(qū)間[ti-1,(ti-1+ti)/2],并更新tB=(tA+tB)/2;如果γ>γmax,說明預報交接點在區(qū)間[(ti-1+ti)/2,ti],更新tA=(tA+tB)/2;

5)設ε為交接時刻計算的精確度,當|tA-tB|≥ε說明精確度還不夠,需轉到步驟4),當|tA-tB|<ε時,即滿足所要求的精確度時,計算終止,得出最后的預報交接時刻為(tA+tB)/2,再次使用預報方法，獲得交接時刻導彈的位置和速度。

雷達交接時刻的計算流程如圖2所示。

圖2 雷達交接時刻的預報流程圖

1.4 基于無跡變換的交接時刻誤差傳播分析

在對預報時刻的誤差傳播情況進行分析時,需要計算預報時刻的狀態(tài)協(xié)方差矩陣,現(xiàn)有的主要方法包括處理線性系統(tǒng)的求解雅克比矩陣方法、將非線性系統(tǒng)擬線性化的協(xié)方差分析描述函數(shù)法、基于無跡變換(Unscented Transform,UT)的協(xié)方差矩陣傳播分析法。本文采用的是基于UT的協(xié)方差傳播分析方法,該方法不需要求解偏導數(shù)矩陣,且不要求系統(tǒng)具有一階可微性,適用于各類線性及復雜的非線性系統(tǒng),與恩克法彈道預報融合后,其協(xié)方差傳播的流程如圖3所示。

1)按照如下樣點選取規(guī)則,構造2n+1個σ點和相應的權重

(15)

2)由恩克法彈道預報方法計算交接時刻下對應樣點的預報狀態(tài)矢量φ

φ(i)=f(σ(i)),i=0,1,…,2n

(16)

式中,函數(shù)f(x)表示使用預報方法求交接時刻對應的狀態(tài)矢量。

(17)

式中,Wi為權重值,與公式(15)相同。

由雷達交接時刻的均值和協(xié)方差矩陣,可以在三維坐標系中以誤差橢球和二維坐標系中以誤差橢圓的形式，對均值和誤差情況進行可視化直觀展示。誤差橢球(圓)的中心為預報均值,其長半軸、短半軸的大小由各個方向上的標準差和誤差描述倍數(shù)決定,其中心軸與坐標軸的夾角可由協(xié)方差矩陣求得,具體求解公式可參閱文獻[9]。

2 仿真分析

2.1 仿真參數(shù)設置

為驗證算法,本文使用了Matlab進行仿真,仿真平臺如下:計算機CPU為Intel Core i7-7700HQ(2.80 GHz),內存為8 GB,顯卡為NVIDIA GeForce GTX 1050,仿真軟件為Matlab2017b。

后置雷達的位置(緯度40°、西經80°、高度500 m),最大探測距離為1 600 km,仰角、方位角的范圍不做限制,利用以上已知條件,對雷達交接時刻進行預報。

在與傳統(tǒng)彈道預報方法計算的雷達交接時刻進行對比時,由于數(shù)值法與本文算法都使用了龍格-庫塔積分方法,兩者的計算時間、計算精度都與積分時長及積分次數(shù)有關,在預報時間相同的條件下,隨著積分步數(shù)的增多,計算精度越高,但計算時間也會越長。為有效檢驗本文算法的效能,分析不同方法計算的精度與運算耗時及誤差傳播情況,綜合考慮了積分步數(shù)和預報精度,恩克法采用一步積分,數(shù)值法采用三步積分。同時由于導彈運動速度高,將交接時刻計算的精確度設置到0.001 s。

根據真實數(shù)據,導彈進入后置雷達探測范圍時刻(雷達交接時刻)的真實值為1 333.831 s,此時導彈的狀態(tài)矢量xh=(-20 170.2-5 679 330.5 5 109 216.7 5 310.1-755.9-1 356.6)T。

2.2 仿真結果及分析

根據上述實驗條件設置,本文分別使用恩克法和傳統(tǒng)的解析法、數(shù)值法對雷達交接時刻及其相應的位置速度進行預報,并使用無跡變換的方法，計算交接時刻的狀態(tài)協(xié)方差矩陣,三種算法的運算時間采用運算100次取平均值獲得。

三種方法計算的雷達交接時刻值及與真實值的誤差如表1所示。對比可知本文恩克法預報的交接時刻的預報精度最高,誤差僅為0.115 s。解析法與數(shù)值法精度相當,誤差分別為0.388 s、0.389 s。

表1 三種方法預報雷達交接時刻值及誤差

三種方法計算交接時刻時的位置速度均值及真實值如表2所示。由表可知,Y軸方向為三個方向中偏離最大,其中解析法偏離大于3 km,數(shù)值法偏離大于1 km,恩克法偏離近1 km。說明解析法預報的狀態(tài)均值誤差較大,數(shù)值法與恩克法預報的位置精度較高且精度相當。

三種方法使用無跡變換計算交接時刻的位置速度標準差結果如表3所示。分析可知,三種方法在不同方向上計算的位置標準差都相差不大,其最大差值僅為4 m,說明三種方法的誤差傳播發(fā)散程度相當。

表2 三種方法預報的位置速度均值及與真實值對比

表3 三種方法在交接時刻下位置速度的標準差

三種方法的運算耗時如表4所示。對比可知,解析法耗時最短,恩克法次之,數(shù)值法耗時最長。原因在于解析法運算只進行了一次二體預報，其計算耗時短,恩克法較解析法多了一次積分運算,而數(shù)值法則需要進行三次積分計算其耗時最長。

表4 三種預報方法運算耗時

為直觀地展示三種方法在交接時刻的誤差傳播情況,本文以蒙特卡洛(Monte Carlo)打靶仿真1 000次獲得的交接時刻散點分布為基準,分別與三種方法計算的交接時刻三維誤差橢球及兩個方向的二維誤差橢圓進行對比,用來驗證算法的準確性和有效性,計算結果比較如圖4、圖5所示。圖中，圓圈表示1 000次蒙特卡洛打靶點,實線表示解析法的計算結果,虛線表示數(shù)值法的計算結果,點劃線表示恩克法的計算結果。由圖中可以看出,三種方法中,解析法誤差橢球(圓)離仿真打靶點最遠，位置精度最差,數(shù)值法、恩克法位置預報精度相當,與前文分析一致。

圖4 仿真打靶(1 000次)交接時刻位置與誤差橢球三維分布情況

圖5 仿真打靶(1 000次)交接時刻位置與誤差橢圓二維分布情況

綜上分析,解析法雖然耗時最短,但是，其預報的交接時刻及狀態(tài)均值的精度最低。數(shù)值法雖然預報的狀態(tài)均值精度與恩克法相當,但交接時刻的時間預報精度低,存在時間與空間不完全匹配的問題,同時運算耗時最長。恩克法對比解析法考慮了攝動因素影響,對比數(shù)值法能夠進行大步長積分,其運算耗時較短,交接時刻及狀態(tài)均值的計算精度最高。

3 結束語

本文介紹了導彈預警中的雷達探測模型及恩克法彈道預報的具體步驟,對雷達交接時刻計算的流程進行了詳細說明,使用無跡變換對預報交接時刻的誤差情況進行了分析。實驗表明,與傳統(tǒng)的解析幾何法和數(shù)值積分法相比,本文算法考慮了導彈運動中受到的攝動力影響,能夠進行大步長的積分運算,在相對較低的運算時間內,提高了雷達交接時刻的預報精度。