劉新建，胡佳鑫，楊樂平

（國防科技大學(xué)空天科學(xué)學(xué)院，長沙410073）

0 引言

彈道落點(diǎn)預(yù)報是彈道導(dǎo)彈預(yù)警和反導(dǎo)中的重要任務(wù)，但彈道導(dǎo)彈發(fā)射陣地的位置估計在彈道導(dǎo)彈防御中也具有重要作用。對于戰(zhàn)略彈道導(dǎo)彈，要準(zhǔn)確判定是哪個敵對國家進(jìn)行的攻擊，何時何地發(fā)射的，以便迅速進(jìn)行反擊；對于戰(zhàn)術(shù)彈道導(dǎo)彈，由于發(fā)射也需要調(diào)平和方位瞄準(zhǔn)，一般是事先安排了有陣地依托的發(fā)射，也要盡可能準(zhǔn)確知道敵方是何時何地發(fā)射的，以便告知偵察衛(wèi)星對該區(qū)域掃描拍照，通過衛(wèi)星圖像進(jìn)一步精確鎖定坐標(biāo)，以便轟炸發(fā)射陣地或彈道導(dǎo)彈貯存基地。

Danish［1］較早研究了天基測量戰(zhàn)術(shù)彈道導(dǎo)彈（TBM）主動段彈道的發(fā)射點(diǎn)參數(shù)估計問題；Isaacson等［2］研究了戰(zhàn)術(shù)彈道導(dǎo)彈的狀態(tài)估計和軌跡預(yù)報。國內(nèi)李盾等［3］基于同步軌道預(yù)警衛(wèi)星對助推段彈道導(dǎo)彈的測量信息，利用改進(jìn)的高斯－牛頓法探討了發(fā)射點(diǎn)參數(shù)的最大似然估計，但必須利用主動段目標(biāo)彈道曲線的先驗(yàn)信息；宋偉等［4］、陽曙光等［5］研究了雙星交匯測量戰(zhàn)術(shù)彈道導(dǎo)彈主動段的彈道參數(shù)，并利用多項(xiàng)式擬合進(jìn)行彈道外推估計發(fā)射點(diǎn)參數(shù)；王瑞等［6］研究了天基紅外低軌衛(wèi)星系統(tǒng)對自由段空間目標(biāo)定位的最小二乘、卡爾曼跟蹤濾波及自適應(yīng)濾波的方法。但鮮見基于對大氣層外被動段彈道的短弧測量對發(fā)射點(diǎn)位置進(jìn)行估計的研究，而且這種多項(xiàng)式彈道擬合外推發(fā)射點(diǎn)的方法較為簡單，對TBM主動段可能有效，但對中遠(yuǎn)程、洲際彈道導(dǎo)彈的被動段測量信息反推發(fā)射點(diǎn)未必有效。

發(fā)射點(diǎn)的位置確定比彈道的落點(diǎn)預(yù)報更難，前者是逆問題；后者是正問題，相對要簡單。

本文基于天基或地基設(shè)備的被動段短弧測量數(shù)據(jù)，假設(shè)在WGS84系中，得到了來襲彈道導(dǎo)彈大氣層外的部分彈道，提出發(fā)射陣地位置估計方法。

1 濾波定軌

（1）坐標(biāo)參數(shù)轉(zhuǎn)換

假設(shè)ti∶xi，yi，zi，vxi，vyi，vzi（i＝1，2，…，n）為彈道導(dǎo)彈在真空中的一段短弧軌跡測量數(shù)據(jù)，分別從WGS84系轉(zhuǎn)換到某一地心慣性系，如J2000 系，得到ti∶Xi，Yi，Zi，Vxi，Vyi，Vzi。

其中，ti為飛行時刻，Xi、Yi、Zi為彈道導(dǎo)彈在該時刻的3個直角坐標(biāo)位置分量，Vxi、Vyi、Vzi為該時刻的3個速度分量。

（2）均值濾波去噪

由于每個測量值含有噪聲，并非是彈道導(dǎo)彈的真實(shí)狀態(tài)，因此需要去噪濾波。工程上有多種數(shù)字濾波方法，如卡爾曼濾波、最小二乘濾波、自適應(yīng)濾波、平均值濾波。特別對于地基雷達(dá)，由于地球曲率的影響，可見弧段極短，而彈道導(dǎo)彈飛行速度極快，獲得的測量數(shù)據(jù)不會很多，因此可采用單點(diǎn)均值濾波。均值濾波最簡單，沒有延遲，對少量的短弧測量數(shù)據(jù)魯棒可靠。

所謂均值濾波，就是依次利用測量數(shù)據(jù)ti：Xi，Yi，Zi，Vxi，Vyi，Vzi（i＝1，2，…，n），根據(jù)真空段自由飛行動力學(xué)模型式 （2），進(jìn)行數(shù)值積分，獲得測量數(shù)據(jù)段之后任選的定軌時刻tf的n個狀態(tài)數(shù)據(jù)。如果測量弧段的數(shù)據(jù)沒有噪聲，那么依次積分，預(yù)報tf時刻的六分量狀態(tài)應(yīng)該是完全相同的；但正因?yàn)榇嬖跍y量噪聲，所以依次積分獲得該時刻的狀態(tài)也存在噪聲，那么求平均值。

可是tf時刻的彈道導(dǎo)彈真實(shí)狀態(tài)只有一個，如果測量數(shù)據(jù)的噪聲遵循白噪聲分布，而且采樣點(diǎn)足夠多，則平均值逼近真值，就得到了彈道導(dǎo)彈真實(shí)狀態(tài)tf：Xf，Yf，Zf，Vxf，Vyf，Vzf。

（3）根據(jù)如下地慣系中空間飛行器的動力學(xué)

微分方程組［7］

其中，地球扁率J2＝1．08263×10－3， 地球引力常數(shù)μ ＝3．986005×1014m3／s2， 地球長半軸＝6378140m，地心矢徑再由上述定軌點(diǎn)的真實(shí)狀態(tài)，進(jìn)行彈道積分預(yù)報再入落點(diǎn)e的直角坐標(biāo)位置和速度Xe、Ye、Ze、Vxe、Vye、Vze。

2 發(fā)射點(diǎn)初始解析估計

為獲得初始估計，不考慮彈道導(dǎo)彈助推段彈道及再入段大氣阻力的影響，假設(shè)導(dǎo)彈在地面瞬間獲得足夠能量，滿足測量彈道所預(yù)報落點(diǎn)的位置和速度要求，不妨稱為脈沖發(fā)射點(diǎn)。脈沖發(fā)射點(diǎn)的位置估計是發(fā)射陣地位置的近似估計，是彈道重構(gòu)尋優(yōu)估計的基礎(chǔ)。

圖1 二體脈沖彈道示意圖Fig．1 Impulse ballistic sketch of two－body problem

1）如圖1所示，根據(jù)再入點(diǎn)e的直角坐標(biāo)和速度Xe、Ye、Ze、Vxe、Vye、Vze， 計算經(jīng)典軌道根數(shù)a、e、i、Ω、ω、M，并求其真近點(diǎn)角fe［8］。

2）求二體脈沖發(fā)射點(diǎn)的真近點(diǎn)角fk

于是不難求得其對應(yīng)的平近點(diǎn)角Mk和偏近點(diǎn)角Ek。

3）求二體飛行時間

其中，a是定軌得到的彈道長半軸，平近點(diǎn)角Me＝M。

4）求脈沖發(fā)射點(diǎn)的經(jīng)緯度

理論上，這是一個已知再入終點(diǎn)位置和飛行時間，反求初始點(diǎn)位置的逆向動力學(xué)問題。這里提出利用二體運(yùn)動軌道的封閉性特點(diǎn)，避開逆向積分問題。

根據(jù)二體運(yùn)動理論，二體脈沖發(fā)射點(diǎn)的直角坐標(biāo)有如下解析公式［9］

再將其從地慣系轉(zhuǎn)換到WGS84系的xk、yk、zk，由此計算脈沖發(fā)射點(diǎn)的經(jīng)緯度λk、фk，就是發(fā)射點(diǎn)的初步估計，誤差與主動段的彈道設(shè)計程序及參數(shù)有關(guān)，對于射程超過10000km的洲際彈道導(dǎo)彈，數(shù)值仿真表明本方法求得發(fā)射點(diǎn)最大誤差為200km。

初始解析估計的計算量很小，適合于彈道導(dǎo)彈防御系統(tǒng)實(shí)時快速預(yù)報，至少能準(zhǔn)確告知彈道導(dǎo)彈來自哪個敵對國家。

3 彈道重構(gòu)尋優(yōu)估計

由于脈沖點(diǎn)的估計與發(fā)射陣地的真實(shí)位置存在誤差，需要進(jìn)一步修正。策略是在脈沖發(fā)射點(diǎn)附近進(jìn)行系列彈道重構(gòu)，仿真模擬彈道導(dǎo)彈從發(fā)射點(diǎn)垂直起飛、主動段程序轉(zhuǎn)彎、末級瞄準(zhǔn)與閉路控制關(guān)機(jī)，優(yōu)選出滿足定軌預(yù)報再入點(diǎn)約束的重構(gòu)彈道，輸出所對應(yīng)的發(fā)射點(diǎn)經(jīng)緯度。

3．1 尋優(yōu)方法

1）以脈沖發(fā)射點(diǎn)及對應(yīng)的自由段飛行時間為中心，設(shè)定經(jīng)緯度誤差邊界范圍，如±2°（對應(yīng)440km邊長的正方形），自由段時間偏差，如±100s，劃分三維網(wǎng)格如圖2所示，遍歷搜索可行解與最優(yōu)解，發(fā)射點(diǎn)位置精度與網(wǎng)格大小有關(guān)，顯然網(wǎng)格越細(xì)，搜索的時間越長。

2）根據(jù)脈沖點(diǎn)與再入落點(diǎn)的經(jīng)緯度，由球面三角公式估計彈道射程。

3）按照射程大小 （洲際、遠(yuǎn)程、中程、近程和短程），參考對應(yīng)的國內(nèi)外典型彈道導(dǎo)彈的參數(shù)，虛擬來襲彈道導(dǎo)彈的動力學(xué)仿真參數(shù)，如級數(shù)、各級的推力、質(zhì)量、秒流量、燃燒時間以及氣動阻力、升力系數(shù)。

4）按照彈道導(dǎo)彈主動段的彈道設(shè)計方法，設(shè)計轉(zhuǎn)彎程序如攻角函數(shù)、各級轉(zhuǎn)彎角速率［8］。

圖2 脈沖發(fā)射位置點(diǎn)附近的三維網(wǎng)格Fig．2 Three dimensional grid near impulse launch location

5）為確保再入位置約束，導(dǎo)彈末級采用Lambert制導(dǎo)方法［7］，根據(jù)末級當(dāng)前位置矢量，再入落點(diǎn)位置矢量和飛行時間要求，求取需要速度，利用速度增益制導(dǎo)構(gòu)成閉路關(guān)機(jī)控制，求得導(dǎo)彈末級飛行的俯仰角、偏航角和關(guān)機(jī)時間。

3．2 約束條件

以上給出了主要搜索變量的邊界約束和方法，但要模擬真實(shí)彈道，還需加入彈道重構(gòu)約束條件，即仿真彈道的再入點(diǎn)位置和速度六分量應(yīng)與彈道測量后的定軌值相同或逼近。

再入點(diǎn)位置誤差約束：Δλ≤ελ，Δφ≤εφ，ΔH ≤εH。

即經(jīng)度、緯度和高程三分量仿真值與定軌值的誤差小于設(shè)定值，因?yàn)槟┘壊捎昧薒ambert制導(dǎo)方法需要速度閉路制導(dǎo)，仿真會自動滿足，簡化了位置約束的求解。

再入角度誤差約束：ΔΘ≤εΘ，Δσ≤εσ，即速度的傾角和偏角誤差小于設(shè)定值

3．3 目標(biāo)函數(shù)

取再入速度模的誤差ΔVe為目標(biāo)函數(shù)，應(yīng)取最小。

4 仿真驗(yàn)證

這里以某三級彈道導(dǎo)彈A為例，給出數(shù)學(xué)仿真予以說明。

（1）導(dǎo)彈A的主要參數(shù)［10］

一級起飛質(zhì)量為34．5t、推力為912kN、工作時間為61s、推進(jìn)劑秒流量為360kg／s；

二級起飛質(zhì)量為12t、推力為270kN、工作時間為65s、推進(jìn)劑秒流量為108kg／s；

三級起飛質(zhì)量為5t、推力為155kN、工作時間為59s、推進(jìn)劑秒流量為62kg／s。

末修級工作時間由閉路關(guān)機(jī)控制指令決定。但彈道數(shù)值仿真時，為簡化起見，可不設(shè)末修級，第三級工作時間不固定，由閉路關(guān)機(jī)仿真時間決定，只要滿足自由段飛行時間約束和落點(diǎn)位置約束就可以。其他中短程的戰(zhàn)術(shù)彈道導(dǎo)彈類似，也可不設(shè)末修級。

（2）發(fā)射點(diǎn)和目標(biāo)點(diǎn)參數(shù)假設(shè)

發(fā)射點(diǎn)：經(jīng)度－118°，緯度30°；

目標(biāo)點(diǎn)：經(jīng)度112°，緯度34°；

兩點(diǎn)之間的射程為11643km。

（3）仿真過程

1）在脈沖估計的基礎(chǔ)上，編寫3自由度仿真優(yōu)化搜索程序。

2）設(shè)計一條從發(fā)射點(diǎn)到目標(biāo)點(diǎn)的標(biāo)稱攻擊彈道，獲得中段N個短弧采樣數(shù)據(jù)，例如50個軌跡點(diǎn)，每隔1秒一個采樣數(shù)據(jù)，適當(dāng)加入白噪聲，基于這N個數(shù)據(jù)，完成定軌計算。

3）基于以上方法，完成發(fā)射點(diǎn)的二體脈沖估計和彈道重構(gòu)尋優(yōu)估計。

（4）仿真結(jié)果

二體脈沖估計：

發(fā)射點(diǎn)位置誤差

經(jīng)度差：－0．127°（對應(yīng)14km）

緯度差：1．0015°（對應(yīng)110km）彈道重構(gòu)尋優(yōu)估計：

發(fā)射點(diǎn)位置誤差

經(jīng)度差：0．125°（對應(yīng)14km）

緯度差：0．0015° （對應(yīng)0．2km）

（5）約束條件逼近度

速度誤差：0．047m／s

當(dāng)?shù)厮俣葍A角誤差：0．85°

速度偏角誤差：0．027°

經(jīng)度差：0．0045° （對應(yīng)0．5km誤差）

緯度差：0．0044° （對應(yīng)0．5km誤差）

圖3給出了上述數(shù)值仿真算例的尋優(yōu)估計重構(gòu)彈道。

圖3 彈道重構(gòu)尋優(yōu)估計數(shù)值仿真軌跡示意圖Fig．3 Numerical simulation sketch of optimal estimation with trajectory regeneration

可見，對于本算例的二體快速脈沖估計，發(fā)射點(diǎn)的位置誤差為120km，這對于判斷來襲彈道導(dǎo)彈的敵對國家是夠用的。而經(jīng)過彈道重構(gòu)的尋優(yōu)估計，發(fā)射點(diǎn)位置誤差估計降到了15km左右，提高了可觀的估計精度，如果發(fā)射點(diǎn)尋優(yōu)網(wǎng)格進(jìn)一步細(xì)分，誤差還可能降低。

5 結(jié)論

提出的彈道導(dǎo)彈發(fā)射參數(shù)估計方法包括脈沖初始估計和彈道重構(gòu)仿真兩部分，前者實(shí)時性好，也是后者的基礎(chǔ)。后者以脈沖點(diǎn)為估計中心，對于戰(zhàn)略彈道導(dǎo)彈，以經(jīng)緯度±2°誤差為位置邊界，以二體飛行時間tke為中心，±100s為自由段的飛行時間邊界，劃分三維網(wǎng)格，并以再入點(diǎn)的位置三分量、當(dāng)?shù)厮俣葍A角、速度偏角共5個彈道重構(gòu)約束條件，以速度誤差為目標(biāo)函數(shù)，進(jìn)行三自由度彈道控制仿真與優(yōu)化搜索，求取發(fā)射陣地坐標(biāo)，這樣來減小發(fā)射點(diǎn)位置估計誤差。對于戰(zhàn)術(shù)彈道導(dǎo)彈，網(wǎng)格邊界可通過仿真進(jìn)一步縮小，方法類似。

這種彈道重構(gòu)方法，不需要完全一致的彈道導(dǎo)彈參數(shù)，只要大致接近就可，因?yàn)槭前凑諒椀缹?dǎo)彈攻擊目標(biāo)所采用的典型彈道規(guī)劃方法，基于動力學(xué)與控制仿真，并通過優(yōu)化搜索參數(shù)和Lambert末制導(dǎo)來滿足再入邊界條件的，使得數(shù)學(xué)仿真彈道能夠逼近實(shí)際飛行軌跡。

只要獲得了來襲彈道測量參數(shù)，本方法可以獲得較高置信度的估計，在彈道導(dǎo)彈預(yù)警和防御中有一定參考價值。