周童，張雅聲，周?？?，漢京濱

(航天工程大學(xué)，北京 101400)

0 引言

彈道導(dǎo)彈的早期預(yù)警設(shè)備包括高軌預(yù)警衛(wèi)星和遠(yuǎn)程預(yù)警雷達(dá)，高軌預(yù)警衛(wèi)星由于位置較高將率先發(fā)現(xiàn)處于主動段的彈道導(dǎo)彈，并利用導(dǎo)彈尾焰的紅外輻射對導(dǎo)彈進(jìn)行跟蹤。遠(yuǎn)程預(yù)警雷達(dá)可以根據(jù)高軌預(yù)警衛(wèi)星的指引在導(dǎo)彈關(guān)機(jī)后對導(dǎo)彈進(jìn)行接力跟蹤[1-2]。對遠(yuǎn)程彈道導(dǎo)彈來說，高軌預(yù)警衛(wèi)星探測到導(dǎo)彈關(guān)機(jī)后，遠(yuǎn)程預(yù)警雷達(dá)由于探測距離限制并不能及時捕獲到導(dǎo)彈目標(biāo)，在這種情況下就需要對高軌預(yù)警衛(wèi)星的探測數(shù)據(jù)進(jìn)行處理，生成預(yù)警信息，對遠(yuǎn)程預(yù)警雷達(dá)的搜索與捕獲過程進(jìn)行指引[3-5]。

高軌預(yù)警衛(wèi)星與遠(yuǎn)程預(yù)警雷達(dá)的交接問題主要涉及2個過程：一是預(yù)警衛(wèi)星探測到導(dǎo)彈關(guān)機(jī)后的彈道及誤差預(yù)報問題；二是2種設(shè)備交接班利用的預(yù)警信息的表示方法。在彈道預(yù)報方面，當(dāng)前的彈道預(yù)報方法主要有解析幾何法和數(shù)值積分法[6]。解析幾何法的依據(jù)是橢圓彈道理論，將導(dǎo)彈自由段的飛行軌跡近似為橢圓彈道，利用解析幾何法實現(xiàn)彈道的預(yù)報[7-9]。數(shù)值積分法的依據(jù)是導(dǎo)彈的動力學(xué)方程，利用數(shù)值逼近的方法計算導(dǎo)彈彈道[10-11]。解析幾何法計算過程簡單，速度較快，但計算精度較低；數(shù)值積分法計算過程復(fù)雜，速度較慢，但計算精度較高?？紤]到計算精度要求，數(shù)值積分法在實際工程中應(yīng)用更為廣泛。在預(yù)警信息表示方面，相關(guān)研究還較少。文獻(xiàn)[12]提出了一種導(dǎo)彈預(yù)警衛(wèi)星引導(dǎo)能力評價方法，在預(yù)警信息含義方面沒有詳細(xì)表述。文獻(xiàn)[13]給出了粗略的預(yù)警信息生成模型，但缺乏理論依據(jù)且空域的劃分比較粗略。

針對以上問題，本文首先給出了基于數(shù)值積分法的彈道及誤差預(yù)報過程，然后將彈道預(yù)報的結(jié)果與概率統(tǒng)計理論中的高斯概率誤差分布相結(jié)合提出了4種預(yù)警信息的表示方法，最后通過仿真驗證了預(yù)警信息生成方法的正確性。

1 基于數(shù)值積分法的彈道預(yù)報

對于遠(yuǎn)程彈道導(dǎo)彈來說，發(fā)動機(jī)關(guān)機(jī)后導(dǎo)彈將在大氣層外無動力飛行，此時可以利用考慮多種攝動因素的二體運(yùn)動模型對導(dǎo)彈運(yùn)動進(jìn)行描述，并以此建立導(dǎo)彈的運(yùn)動學(xué)方程。由于導(dǎo)彈的動力學(xué)方程為變系數(shù)非線性常微分方程，因此可以采用數(shù)值積分的方法對方程進(jìn)行求解。數(shù)值積分法計算彈道的一般過程為首先建立動力學(xué)方程，而后選取適當(dāng)?shù)姆e分方法按照一定的積分步長進(jìn)行計算，在給定初始狀態(tài)的情況下，可以得到任意時刻導(dǎo)彈的運(yùn)動狀態(tài)。

1.1 建立目標(biāo)運(yùn)動學(xué)方程

(1)

式(1)為時變的非線性微分方程，無法獲得解析解，通常采用數(shù)值積分的方法獲得其數(shù)值解。數(shù)值積分的過程可以描述為：已知目標(biāo)在k-1時刻的狀態(tài)xk-1，那么可以通過積分獲得目標(biāo)在下一時刻k的運(yùn)動狀態(tài)為

(2)

利用式(2)，對導(dǎo)彈動力學(xué)方程采用合適的數(shù)值積分方法，就可以計算得到導(dǎo)彈彈道[14]。

1.2 彈道誤差傳播分析

在利用數(shù)值積分法計算彈道的過程中，通常利用協(xié)方差分析描述函數(shù)法(CADET)對彈道誤差的傳播進(jìn)行分析，CADET的一般過程為首先對系統(tǒng)進(jìn)行擬線性化處理，然后利用線性化系統(tǒng)協(xié)方差分析原理計算隨機(jī)狀態(tài)變量的協(xié)方差。

(3)

式中：m為位置矢量r的均值。

因此擬線性系統(tǒng)的動態(tài)矩陣描述函數(shù)N可以表示為

(4)

式中：I為3×3的單位矩陣；

利用CADET方法可得導(dǎo)彈的狀態(tài)均值和協(xié)方差的傳播方程：

(5)

2 基于統(tǒng)計理論的狀態(tài)預(yù)測誤差概率分布

彈道導(dǎo)彈在發(fā)動機(jī)關(guān)機(jī)后進(jìn)入中段飛行，由于中段飛行環(huán)境的復(fù)雜性、彈道目標(biāo)和衛(wèi)星運(yùn)動的非線性以及預(yù)報算法模型的建模精度等因素的影響，對目標(biāo)運(yùn)動狀態(tài)預(yù)測誤差分布進(jìn)行精確建模分析非常困難。通常的做法是將這些誤差源帶來的隨機(jī)誤差統(tǒng)一處理并建模為服從零均值的高斯分布，并將傳感器視線觀測噪聲建模為零均值高斯噪聲。根據(jù)概率統(tǒng)計理論中心極限定理，可以將目標(biāo)運(yùn)動狀態(tài)預(yù)測誤差建模為高斯分布。

基于統(tǒng)計理論的高斯分布概率模型給出了狀態(tài)估計量在高斯假設(shè)條件下的概率分布計算[15]，下面對其進(jìn)行分析。

設(shè)ξ為狀態(tài)量m的估計，ξ為隨機(jī)變量且滿足高斯分布，則ξ的概率分布函數(shù)為

(6)

式中:N為狀態(tài)矢量ξ的維數(shù)；m為狀態(tài)矢量真實值；P為狀態(tài)估計量ξ的狀態(tài)誤差協(xié)方差矩陣，即

P=E[(ξ-m)(ξ-m)T].

(7)

若P為對稱正定矩陣，則存在逆陣P-1，且P-1為對稱正定矩陣。令

κ=(ξ-m)P-1(ξ-m)T.

(8)

(9)

(10)

可知，式(9)和式(10)不但描述了誤差概率密度隨誤差絕對值的變化，而且描述了誤差概率密度隨誤差方向的變化規(guī)律。

對位置誤差分布概率展開研究,令

(11)

(12)

根據(jù)矩陣分析理論，存在且唯一存在正交矩陣A，有AT=A(AT為A的轉(zhuǎn)置)，且使得

(13)

式中:λ1,λ2,λ3為協(xié)方差矩陣Pr的特征值。

令ζ=AT,γ，并將式(13)代入式(12)中，整理得

(14)

(15)

(16)

式中:Γ()為gamma函數(shù);erf()為誤差函數(shù),表達(dá)形式為

(17)

3 預(yù)警信息生成模型

不同的交接班任務(wù)對預(yù)警信息的要求是不一樣的。較寬視場的接班設(shè)備覆蓋空域大，只需預(yù)警信息提供目標(biāo)的概略空域分布；接班設(shè)備視場再小一些，就需要縮小接班傳感器搜索的空域，此時需要預(yù)警信息提供目標(biāo)較為精確的空域分布；對于采用窄視場搜索的接班設(shè)備來說，為保證搜索數(shù)據(jù)率，快速發(fā)現(xiàn)目標(biāo)，則需要預(yù)警信息提供盡可能精確的目標(biāo)空域分布，同時還能指示目標(biāo)在空域中出現(xiàn)的概率。本節(jié)將按照上述思路，由精細(xì)到粗略分別研究建立概率橢球模型、概率球模型、均勻橢球模型和均勻球模型4種粒度的預(yù)警信息模型。

3.1 概率橢球模型

概率橢球模型是以目標(biāo)預(yù)測狀態(tài)的分布概率密度為基礎(chǔ)，用誤差橢球精細(xì)描述目標(biāo)狀態(tài)的分布概率。根據(jù)式(9)和式(10)可知，第k時刻目標(biāo)位置與速度狀態(tài)預(yù)測值的概率分布密度函數(shù)分別為

(18)

式中:Pr,k和Pv,k都是3×3的非負(fù)正定Hermitian矩陣，因此存在3×3的正定矩陣A和B，有AAT=I，BBT=I，且使得

(19)

(20)

式中：λr,i，λv,i(i=1,2,3)分別為協(xié)方差矩陣Pr,k和Pv,k的特征值，且有如下關(guān)系：

(21)

(22)

由式(22)可知，對于給定的概率Pr和Pv，相應(yīng)的κ的取值與狀態(tài)變量無關(guān)，可以直接求解κr和κv的值。下面不加區(qū)分的給出概率P取不同值時，利用Matlab符號工具箱求解得到的對應(yīng)的κ值，即

(23)

概率橢球模型是以預(yù)測值為原點(diǎn)、等概率面為橢球面的三維聯(lián)合高斯分布模型。它不僅描述了目標(biāo)狀態(tài)的分布概率密度隨偏離預(yù)測值距離的變化規(guī)律和變化方向，而且描述了目標(biāo)狀態(tài)位于橢球體內(nèi)的實際分布概率大小。

3.2 概率球模型

概率球模型以目標(biāo)預(yù)測狀態(tài)的分布概率密度為基礎(chǔ)，用誤差橢球最大誤差半徑所在軸方向的概率分布密度來描述狀態(tài)的概率分布密度，用誤差球來描述目標(biāo)狀態(tài)的分布概率。概率球模型的概率分布密度仍是高斯型概率分布密度函數(shù)的形式，但是，概率分布的大小以球體來描述。

假設(shè)已由式(19)和式(20)得到第k時刻目標(biāo)位置與速度狀態(tài)估計誤差協(xié)方差矩陣的特征值，并且特征值滿足如式(21)所示，則第k時刻目標(biāo)位置與速度狀態(tài)預(yù)測值的概率分布密度分別表示為

(24)

概率球在概率橢球的基礎(chǔ)上，將原來概率橢球各軸向不同的概率分布密度統(tǒng)一以同一個概率分布密度來描述目標(biāo)預(yù)測狀態(tài)的概率分布密度。其概率積分本質(zhì)上還是對高斯概率分布密度積分，利用前文推導(dǎo)的結(jié)果，直接給出以下結(jié)論。

概率球模型是以預(yù)測值為原點(diǎn)、等概率面為球面的三維聯(lián)合高斯分布模型。它描述了目標(biāo)狀態(tài)的分布概率密度隨偏離預(yù)測值距離的變化規(guī)律，但沒有描述概率密度變化方向，即認(rèn)為在概率球內(nèi)目標(biāo)狀態(tài)分布概率密度是各向同性分布。概率球模型也描述了目標(biāo)狀態(tài)位于概率球內(nèi)的概率分布大小，但是，由于其弱化了目標(biāo)狀態(tài)的實際概率分布密度，所以目標(biāo)在概率球內(nèi)的實際概率分布大于理論得到的概率分布。

3.3 均勻橢球模型

均勻橢球模型是在概率橢球模型的基礎(chǔ)上，認(rèn)為在橢球內(nèi)部目標(biāo)狀態(tài)的概率分布密度是均勻的。因此，采用均勻橢球模型的預(yù)警信息首先需要基于概率橢球模型，然后根據(jù)給定的概率分布分別計算橢球體三軸半徑。

假設(shè)已由式(19)和式(20)得到第k時刻目標(biāo)位置與速度狀態(tài)估計誤差協(xié)方差矩陣的特征值，并且特征值滿足如式(21)所示，則第k時刻目標(biāo)位置與速度狀態(tài)預(yù)測值的概率分布密度分別表示為

(25)

均勻橢球模型弱化了目標(biāo)狀態(tài)的概率密度分布特性，將內(nèi)部每一點(diǎn)都視為均勻分布，但其描述了目標(biāo)狀態(tài)的概率分布的變化方向。因此，采用均勻橢球模型生成的預(yù)警信息能夠用橢球三軸半徑的大小來描述目標(biāo)狀態(tài)概率分布的變化方向，但是不能描述目標(biāo)狀態(tài)概率分布密度的變化趨勢。

3.4 均勻球模型

均勻球模型是在概率球的基礎(chǔ)上，認(rèn)為誤差球內(nèi)部目標(biāo)狀態(tài)的概率分布密度是均勻的。

假設(shè)已由式(19)和式(20)得到第k時刻目標(biāo)位置與速度狀態(tài)估計誤差協(xié)方差矩陣的特征值，則第k時刻目標(biāo)位置與速度狀態(tài)預(yù)測值的概率分布密度分別表示為

(26)

均勻球模型弱化了目標(biāo)狀態(tài)的概率密度分布特性，認(rèn)為目標(biāo)狀態(tài)概率分布各向同性，且內(nèi)部每一點(diǎn)都為均勻分布。因此，采用均勻球模型生成的預(yù)警信息粒度是最粗的，在4種模型中能夠提供的指示信息最少，精確最差，能力最弱。但是，該模型形式簡單且易于描述，在接班設(shè)備視場覆蓋能力較強(qiáng)、覆蓋空域較大的情況下，能夠達(dá)到較好的引導(dǎo)效率。

4 仿真校驗

本節(jié)將通過仿真，對上文提出的高軌預(yù)警衛(wèi)星預(yù)警信息生成過程進(jìn)行模擬分析，并利用Monte Carlo仿真驗證預(yù)警信息生成方法的正確性。

4.1 彈道誤差傳播與預(yù)警信息生成

設(shè)預(yù)警衛(wèi)星估計得到導(dǎo)彈在關(guān)機(jī)點(diǎn)時刻t0的狀態(tài)均值m0和協(xié)方差P0分別為m0=(2 155 145 m,6 002 997 m，-1 069 m，-426 m/s，5 867 m/s,1 843 m/s)，

利用數(shù)值積分法對導(dǎo)彈狀態(tài)均值和協(xié)方差進(jìn)行預(yù)報計算，得到t=600 s時的導(dǎo)彈狀態(tài)均值m和協(xié)方差P。數(shù)值積分采用四階龍格-庫塔方法，積分步長設(shè)為0.1 s。圖1給出了導(dǎo)彈狀態(tài)變量均值隨時間變化的仿真結(jié)果。

利用計算得到的t=600 s時的導(dǎo)彈狀態(tài)均值和協(xié)方差可以根據(jù)本文提出的預(yù)警信息生成模型得到導(dǎo)彈的預(yù)警信息。由前文分析可知，基于概率橢球模型生成的預(yù)警信息最具代表性，因此本文以概率橢球模型預(yù)警信息生成為例開展仿真。圖2為根據(jù)t=600s時的導(dǎo)彈狀態(tài)生成的地心慣性坐標(biāo)系下導(dǎo)彈位置和速度的預(yù)警信息。

圖中橢球體所包含的區(qū)域為基于本文預(yù)警信息生成方法計算得到的理論誤差分布，外圈橢球體以99.5%的概率涵蓋導(dǎo)彈位置和速度分布區(qū)域，內(nèi)圈橢球體以95%的概率涵蓋導(dǎo)彈位置和速度分布區(qū)域。值得注意的是誤差橢球的三軸方向并不一定平行于地心慣性系的3個坐標(biāo)軸，橢球的三軸方向為誤差變化量最大的3個方向，可由協(xié)方差矩陣的特征向量確定。

4.2 預(yù)警信息生成方法正確性驗證

本節(jié)采用的驗證方法如下：根據(jù)t0時刻的導(dǎo)彈關(guān)機(jī)點(diǎn)狀態(tài)估計值及誤差協(xié)方差矩陣作為狀態(tài)初值，根據(jù)數(shù)值積分法預(yù)報t=600 s時的導(dǎo)彈運(yùn)動狀態(tài)，并基于本文提出的預(yù)警信息生成方法計算t=600 s時的導(dǎo)彈狀態(tài)的理論誤差分布。同時，取t0時刻導(dǎo)彈真實彈道數(shù)據(jù)作為真值，以t0時刻導(dǎo)彈運(yùn)動狀態(tài)誤差為標(biāo)準(zhǔn)差，預(yù)報導(dǎo)彈在t=600 s時的運(yùn)動狀態(tài)，以500次Monte Carlo仿真結(jié)果作為導(dǎo)彈狀態(tài)預(yù)測的實際分布。最后，對比考察理論誤差分布和實際分布在統(tǒng)計意義上是否符合一致性。

圖3為導(dǎo)彈位置預(yù)警信息生成方法正確性驗證結(jié)果，三維圖中*為500次Monte Carlo生成的導(dǎo)彈實際位置，外圈橢球體以99.5%的概率涵蓋導(dǎo)彈位置分布區(qū)域，內(nèi)圈橢球體以95%的概率涵蓋導(dǎo)彈位置分布區(qū)域。為直觀比較，將導(dǎo)彈理論狀態(tài)與實際狀態(tài)分別投影到地心慣性坐標(biāo)系的3個主平面，橢圓中心為導(dǎo)彈狀態(tài)均值。其中‘—’代表誤差分布概率為99.5%，‘--’代表誤差分布概率為95%。

圖4為導(dǎo)彈速度預(yù)警信息生成方法正確性驗證結(jié)果，圖中坐標(biāo)軸與各標(biāo)識含義與圖3中相同。

5 結(jié)束語

本文以導(dǎo)彈防御為背景，研究了導(dǎo)彈在發(fā)動機(jī)關(guān)機(jī)后自由飛行過程中運(yùn)動誤差的傳播特性及誤差的描述方法。文章首先給出了基于數(shù)值積分法的彈道誤差傳播方法，而后基于數(shù)值積分法獲得的導(dǎo)彈狀態(tài)由細(xì)到粗提出了4種導(dǎo)彈預(yù)警信息模型，并通過仿真驗證了所提出模型的正確性。文章的研究內(nèi)容可以推廣到我國導(dǎo)彈預(yù)警系統(tǒng)的建設(shè)中，但在實際運(yùn)用中預(yù)警信息并非越精細(xì)越好，預(yù)警信息越詳細(xì)，交班設(shè)備生成預(yù)警信息的時間就越多，整個交接班過程效率反而可能達(dá)不到最優(yōu)，根據(jù)具體任務(wù)選擇合適的預(yù)警信息生成方法將會大大提高導(dǎo)彈預(yù)警系統(tǒng)的防御能力。