高軌預警衛星預警信息生成方法研究*

周童，張雅聲，周海俊，漢京濱

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

0 引言

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

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

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

1 基于數值積分法的彈道預報

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

1.1 建立目標運動學方程

(1)

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

(2)

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

1.2 彈道誤差傳播分析

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

(3)

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

因此擬線性系統的動態矩陣描述函數N可以表示為

(4)

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

利用CADET方法可得導彈的狀態均值和協方差的傳播方程：

(5)

2 基于統計理論的狀態預測誤差概率分布

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

基于統計理論的高斯分布概率模型給出了狀態估計量在高斯假設條件下的概率分布計算[15]，下面對其進行分析。

設ξ為狀態量m的估計，ξ為隨機變量且滿足高斯分布，則ξ的概率分布函數為

(6)

式中:N為狀態矢量ξ的維數；m為狀態矢量真實值；P為狀態估計量ξ的狀態誤差協方差矩陣，即

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

(7)

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

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

(8)

(9)

(10)

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

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

(11)

(12)

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

(13)

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

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

(14)

(15)

(16)

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

(17)

3 預警信息生成模型

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

3.1 概率橢球模型

概率橢球模型是以目標預測狀態的分布概率密度為基礎，用誤差橢球精細描述目標狀態的分布概率。根據式(9)和式(10)可知，第k時刻目標位置與速度狀態預測值的概率分布密度函數分別為

(18)

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

(19)

(20)

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

(21)

(22)

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

(23)

概率橢球模型是以預測值為原點、等概率面為橢球面的三維聯合高斯分布模型。它不僅描述了目標狀態的分布概率密度隨偏離預測值距離的變化規律和變化方向，而且描述了目標狀態位于橢球體內的實際分布概率大小。

3.2 概率球模型

概率球模型以目標預測狀態的分布概率密度為基礎，用誤差橢球最大誤差半徑所在軸方向的概率分布密度來描述狀態的概率分布密度，用誤差球來描述目標狀態的分布概率。概率球模型的概率分布密度仍是高斯型概率分布密度函數的形式，但是，概率分布的大小以球體來描述。

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

(24)

概率球在概率橢球的基礎上，將原來概率橢球各軸向不同的概率分布密度統一以同一個概率分布密度來描述目標預測狀態的概率分布密度。其概率積分本質上還是對高斯概率分布密度積分，利用前文推導的結果，直接給出以下結論。

概率球模型是以預測值為原點、等概率面為球面的三維聯合高斯分布模型。它描述了目標狀態的分布概率密度隨偏離預測值距離的變化規律，但沒有描述概率密度變化方向，即認為在概率球內目標狀態分布概率密度是各向同性分布。概率球模型也描述了目標狀態位于概率球內的概率分布大小，但是，由于其弱化了目標狀態的實際概率分布密度，所以目標在概率球內的實際概率分布大于理論得到的概率分布。

3.3 均勻橢球模型

均勻橢球模型是在概率橢球模型的基礎上，認為在橢球內部目標狀態的概率分布密度是均勻的。因此，采用均勻橢球模型的預警信息首先需要基于概率橢球模型，然后根據給定的概率分布分別計算橢球體三軸半徑。

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

(25)

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

3.4 均勻球模型

均勻球模型是在概率球的基礎上，認為誤差球內部目標狀態的概率分布密度是均勻的。

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

(26)

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

4 仿真校驗

本節將通過仿真，對上文提出的高軌預警衛星預警信息生成過程進行模擬分析，并利用Monte Carlo仿真驗證預警信息生成方法的正確性。

4.1 彈道誤差傳播與預警信息生成

設預警衛星估計得到導彈在關機點時刻t0的狀態均值m0和協方差P0分別為m0=(2 155 145 m,6 002 997 m，-1 069 m，-426 m/s，5 867 m/s,1 843 m/s)，

利用數值積分法對導彈狀態均值和協方差進行預報計算，得到t=600 s時的導彈狀態均值m和協方差P。數值積分采用四階龍格-庫塔方法，積分步長設為0.1 s。圖1給出了導彈狀態變量均值隨時間變化的仿真結果。

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

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

4.2 預警信息生成方法正確性驗證

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

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

圖4為導彈速度預警信息生成方法正確性驗證結果，圖中坐標軸與各標識含義與圖3中相同。

5 結束語

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