基于恩克法的雷達(dá)交接時(shí)刻預(yù)報(bào)與誤差分析

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

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

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

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

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

1 雷達(dá)交接時(shí)刻預(yù)報(bào)算法

1.1 雷達(dá)探測(cè)模型

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

1)將雷達(dá)站位置(L,B,H)轉(zhuǎn)換為地固系坐標(biāo)

(1)

2)將地固系下導(dǎo)彈對(duì)雷達(dá)的相對(duì)位置轉(zhuǎn)換為雷達(dá)直角坐標(biāo)系ENU下的坐標(biāo)

(2)

(3)

3)將導(dǎo)彈在雷達(dá)直角坐標(biāo)系的坐標(biāo)轉(zhuǎn)化為雷達(dá)探測(cè)數(shù)據(jù)(R,A,E)

(4)

1.2 恩克法彈道預(yù)報(bào)方法

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

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

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

軌道確定后,軌道的偏近點(diǎn)角E,平近點(diǎn)角M,平均角速度n及過近地點(diǎn)時(shí)刻τ的計(jì)算如下

(5)

(6)

2)求解t時(shí)刻真實(shí)軌道與基準(zhǔn)軌道的偏差量δr

本文分析導(dǎo)彈的受力情況,可得實(shí)際軌道和基準(zhǔn)軌道的動(dòng)力學(xué)微分方程分別為

(7)

(8)

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

由基準(zhǔn)軌道的定義可知δr=r-ρ,可得

(9)

(10)

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

3)計(jì)算t時(shí)刻導(dǎo)彈的狀態(tài)量r

在獲得t時(shí)刻導(dǎo)彈基準(zhǔn)軌道的位置、速度及其偏差量后,t時(shí)刻真實(shí)軌道的位置、速度可由兩者之和得到,具體表達(dá)式分別為：

(11)

(12)

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

(13)

(14)

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

圖1 恩克法彈道預(yù)報(bào)流程圖

1.3 雷達(dá)交接時(shí)刻預(yù)報(bào)流程

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

1)確定彈道預(yù)報(bào)的時(shí)間范圍[t0,tf],其中t0為初始時(shí)刻,tf為導(dǎo)彈落點(diǎn)預(yù)報(bào)時(shí)刻,落點(diǎn)預(yù)報(bào)時(shí)刻可通過求解導(dǎo)彈的二體運(yùn)動(dòng)方程與地球理想球面方程計(jì)算獲得;

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

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

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

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

雷達(dá)交接時(shí)刻的計(jì)算流程如圖2所示。

圖2 雷達(dá)交接時(shí)刻的預(yù)報(bào)流程圖

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

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

1)按照如下樣點(diǎn)選取規(guī)則,構(gòu)造2n+1個(gè)σ點(diǎn)和相應(yīng)的權(quán)重

(15)

2)由恩克法彈道預(yù)報(bào)方法計(jì)算交接時(shí)刻下對(duì)應(yīng)樣點(diǎn)的預(yù)報(bào)狀態(tài)矢量φ

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

(16)

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

(17)

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

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

2 仿真分析

2.1 仿真參數(shù)設(shè)置

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

后置雷達(dá)的位置(緯度40°、西經(jīng)80°、高度500 m),最大探測(cè)距離為1 600 km,仰角、方位角的范圍不做限制,利用以上已知條件,對(duì)雷達(dá)交接時(shí)刻進(jìn)行預(yù)報(bào)。

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

根據(jù)真實(shí)數(shù)據(jù),導(dǎo)彈進(jìn)入后置雷達(dá)探測(cè)范圍時(shí)刻(雷達(dá)交接時(shí)刻)的真實(shí)值為1 333.831 s,此時(shí)導(dǎo)彈的狀態(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 仿真結(jié)果及分析

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

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

表1 三種方法預(yù)報(bào)雷達(dá)交接時(shí)刻值及誤差

三種方法計(jì)算交接時(shí)刻時(shí)的位置速度均值及真實(shí)值如表2所示。由表可知,Y軸方向?yàn)槿齻€(gè)方向中偏離最大,其中解析法偏離大于3 km,數(shù)值法偏離大于1 km,恩克法偏離近1 km。說明解析法預(yù)報(bào)的狀態(tài)均值誤差較大,數(shù)值法與恩克法預(yù)報(bào)的位置精度較高且精度相當(dāng)。

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

表2 三種方法預(yù)報(bào)的位置速度均值及與真實(shí)值對(duì)比

表3 三種方法在交接時(shí)刻下位置速度的標(biāo)準(zhǔn)差

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

表4 三種預(yù)報(bào)方法運(yùn)算耗時(shí)

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

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

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

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

3 結(jié)束語

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