楊 冬 徐 勁 杜建麗

（1 中國科學(xué)院紫金山天文臺 南京 210023）

（2 中國科學(xué)院空間目標與碎片觀測重點實驗室 南京 210023）

1 引言

彈道導(dǎo)彈預(yù)警在航天技術(shù)和戰(zhàn)略預(yù)警中具有廣闊的應(yīng)用前景, 正日趨成為國內(nèi)外研究的熱點課題. 彈道導(dǎo)彈的整運行過程分為助推段、中段和再入段三部分. 其中中段又叫自由飛行段, 是彈道導(dǎo)彈關(guān)閉發(fā)動機后在大氣層外飛行的過程, 這個階段飛行時間最長, 占整個彈道飛行時間的80%以上,且僅受自然力影響, 其軌道特征容易被掌握. 因此中段是各種無線測量設(shè)備進行捕獲、跟蹤、身份鑒別、計算預(yù)警信息（如發(fā)、落點時間、位置）的主要階段.

在此前的研究中, 彈道導(dǎo)彈中段的身份鑒別和發(fā)、落點信息推算是兩個不同的課題. 文獻[1–3]等都基于彈道導(dǎo)彈彈頭在雷達回波下反映出的物理特征, 例如彈頭的微動（如進動、章動、雷達散射截面（Radar Cross Section, RCS）的變化）存在周期性特點, 可通過雷達回波呈現(xiàn)的微多普勒效應(yīng)進行分析和識別; 又例如雷達對飛行目標進行的高分辨率一維成像（High Resolution Range Profile, HRRP）具有成像速度快、分辨率高的特點,且HRRP能提供沿雷達視線方向強散射中心的分布信息, 具有目標的結(jié)構(gòu)特征, 根據(jù)獲取的特征和目標特征庫的模板對比, 即可完成目標身份鑒別.而關(guān)于發(fā)、落點信息推算的研究都是在已知目標為彈道導(dǎo)彈的前提下, 利用彈道與地球有交點這一特點進行推算. 推算方法多是如文獻[4]介紹的使用卡爾曼濾波或改進的卡爾曼濾波對目標運動方程進行外推的方法進行計算. 這類外推方法本質(zhì)上是一種試探法, 即每一次外推結(jié)束后需要判斷是否進入地球內(nèi)部, 如果進入地球內(nèi)部則需采用二分法反復(fù)外推最終求得結(jié)果. 其計算效率與濾波步長的選取密切相關(guān), 最合適的濾波時間步長又與具體的彈道有關(guān), 導(dǎo)致很難有普適的最優(yōu)步長可供選擇. 此外這種方法未考慮目標在飛行時受到的攝動力, 也沒有加入高程信息的支持. 也有文獻提出把導(dǎo)彈軌道簡化為二體模型下的橢圓軌道后通過幾何關(guān)系求解彈道與地球交點[5], 但這種方法同樣沒有考慮攝動因素和高程的影響. 文獻[6]則另辟蹊徑, 根據(jù)導(dǎo)彈軌道與地球有交點, 而衛(wèi)星軌道與地球無交點這一特征, 提出了一種同時進行快速鑒別目標身份和生成預(yù)警信息的方法. 該方法計算速度快、處理過程簡單穩(wěn)健、計算效率高, 且考慮了地球?qū)δ繕说闹饕獢z動項. 但在使用該方法生成預(yù)警信息時,地球被粗略地簡化為一個參考橢球體, 高程信息依然被略去. 這樣雖然減少了計算量、提高了計算速度, 但是如果導(dǎo)彈的發(fā)點或落點位于高山、高原等地表與參考橢球面距離較大的地區(qū), 使用該方法會降低預(yù)警信息的精度, 且高程值越大, 精度降低得顯然也就越多. 因此有必要改進文獻[6]介紹的預(yù)警信息推算方法, 使新的方法在保持較高計算效率和穩(wěn)健、穩(wěn)定度的前提下盡量消除高程因素引起的誤差, 提高預(yù)警信息的精度. 本文將針對這一問題展開理論分析和探討.

用于支持發(fā)、落點推算的高程信息庫可使用全球數(shù)字高程模型（Digital Elevation Model,DEM）, 它是對地球表面地形地貌的一種離散數(shù)學(xué)表達. DEM的原理是將一片區(qū)域劃分為b行c列（b、c均為正整數(shù)）的四邊形, 計算每個四邊形的平均高程, 然后以三維矩陣的方式存儲這些平均高程. 用函數(shù)的形式可描述為:

式中,Lj、ψj是大地經(jīng)緯度坐標,Hj是經(jīng)緯度坐標（Lj,ψj）處的高程. 目前常見的高分辨率全球DEM的生成主要得益于全球衛(wèi)星遙感、衛(wèi)星測高、船載測深等資料的獲取[7]. 一般來說, DEM里小四邊形的邊長就是DEM的分辨率. 常見的分辨率有90 m、30 m、12.5 m、5 m等. 分辨率是描述地形精確程度的一個重要指標, 邊長越小同一片區(qū)域分割的四邊形越多, 這就意味著分辨率越高, 描述地形越精確.

2 彈道導(dǎo)彈預(yù)警的原始方法

為了后續(xù)闡述方便, 簡要介紹一下文獻[6]提出的方法. 當(dāng)導(dǎo)彈處于發(fā)點或落點時, 其所處的彈道位置與地球表面相交, 導(dǎo)彈地心距r（f）與其星下點地心距R（f）相等, 即

其中f為目標于t時刻的真近點角. 求解（1）式的方法就是通過軌道運動特點鑒別目標身份和推算預(yù)警信息的核心內(nèi)容,f有解就證明目標與地球有交點, 通過軌道根數(shù)就能夠推算出發(fā)、落點等預(yù)警信息. 直接求解（1）式難度較高, 因此使用兩個步驟迭代求解: 第一步迭代是在二體運動模型下求解初始值. 迭代初值是由初始根數(shù)求得的目標真近點角f0, 收斂需要滿足的判定門限為|rk-Rk|<ε, 其中|rk-Rk|是第k次迭代后,目標地心距與星下點地心距之差的絕對值,ε >0是收斂門限,k表示第k次迭代. 二體運動模型下,r可表達為:

其中a、e分別為開普勒根數(shù)中的軌道半長徑和偏心率, 而R可表達為:

其中i、ω分別為開普勒根數(shù)中的軌道傾角和近地點幅角,Rp為地球極半徑,δ為地球扁率.

之后第二步是引入一階長期項和一階短周期項的影響, 對第一步處理獲取的初值進行修正,得到精確的交點時刻. 這一步驟的收斂門限為|r′k - R′k| < ε, 其中k′表示第k′次迭代. 這一步迭代收斂后可得到精確的交點時刻T0, 通過交點時刻可進一步求得目標在地固坐標系的位置矢量以及軌道與地球表面精確交點的地理經(jīng)緯度L和ψ.

3 有DEM支持的精確交點改進算法

原方法得出軌道與參考橢球體交點信息后, 為了進一步求得DEM支持下更精確的交點, 本文提出在原方法計算結(jié)果的基礎(chǔ)上進行新一輪的迭代計算的新方法.

在新一輪的迭代計算中, 收斂的條件設(shè)置為|HD-HC| < ε2, 其中HC為計算求得的t時刻目標的高程,HD為通過DEM查詢t時刻目標位置的經(jīng)緯度得出的高程,ε2>0是給定的門限,與原方法的門限是不同的值. 本輪迭代是以高程差小于特定門限作為收斂標準, 因此該新一輪迭代過程可稱為“高程迭代”部分. 具體迭代方法如下:

（1）使用原方法計算得出的交點地理經(jīng)緯度L和ψ查詢DEM, 得出高程HD. 再通過地固坐標系至站心地平坐標系的旋轉(zhuǎn)平移轉(zhuǎn)換, 可求得目標t時刻在以目標自身為原點的地平坐標系下垂直方向的速度vhz, 求得時間修正量?t為:

（2）令T=t+?t, 加入一階長期項和一階短周期項的影響, 重新外推計算得到t時刻目標的瞬時根數(shù)σ, 并由該根數(shù)計算t時刻的目標地心距r和目標在軌道坐標系中的位置矢量→r和速度矢量˙→r.

（3）通過軌道坐標系至地固坐標系的坐標旋轉(zhuǎn)變換可得到t時刻目標在地固坐標系中的位置矢量, 由按以下文獻[8]計算目標軌道與地球表面精確交點的地理經(jīng)緯度L、ψ和高程HC:

其中Re是地球赤道半徑,x、y、z分別為目標在地固坐標系位置矢量的三個分量, 此外:

（4）使用L、ψ查詢DEM得出高程HD, 并通過坐標轉(zhuǎn)換求得目標t時刻在以自身為原點的地平坐標系下垂直方向的速度vhz. 由下式求t時刻的修正量:

（5）令T=t+ ?t回到本迭代方法的第二步繼續(xù)計算t和?t, 回到本步驟時, 如果滿足迭代條件|HD-HC| < ε2, 則終止迭代計算過程, 否則依然回到第二步繼續(xù)計算.

迭代收斂以后, 得到的計算結(jié)果就是目標位于彈道和地表的交點時更為精確的時刻t和發(fā)、落點經(jīng)緯度L和ψ, 以及根據(jù)該經(jīng)緯度查詢DEM得到的高程值HD.

正如本文開頭所述, DEM數(shù)據(jù)庫中儲存的高程點是離散的, 這就意味著目標軌道與地表交點區(qū)域的高程越高、坡度越大、DEM自身分辨率越高時, 可能會導(dǎo)致連續(xù)兩次查詢得到的高程數(shù)據(jù)HD的差異較大, 從而使得迭代難以收斂. 為了應(yīng)對這種不收斂的情況, 首先可以增大收斂門限ε2的值; 其次可以使用優(yōu)選策略: 設(shè)定最大迭代次數(shù)kmax, 當(dāng)?shù)螖?shù)達到kmax次時停止迭代, 選擇歷次迭代中|HD- HC|的值最小的那一次的迭代計算求得的L、ψ以及相應(yīng)查詢DEM獲得的高程HD作為最優(yōu)近似值, 視為最終計算結(jié)果.

4 數(shù)值實驗和分析

由于缺乏實測彈道數(shù)據(jù), 我們進行了仿真測試和分析來評估該方法的計算效能. 測試中使用的DEM數(shù)據(jù)來源于中科院網(wǎng)信中心地理數(shù)據(jù)云平臺的SRTMDEM （Shuttle Radar Topography Mission DEM, 航天飛機雷達地形測量任務(wù)DEM數(shù)據(jù)） 90 m分辨率高程數(shù)據(jù)1http://www.gscloud.cn/., 該數(shù)據(jù)是在WGS84（World Geodetic System 84）地固系參考橢球體上的投影. 首先利用已有的成熟彈道設(shè)計軟件設(shè)計得到了5組彈道, 再用試探法通過高精度數(shù)值外推得到了這5組彈道與地球表面的精確交點作為發(fā)、落點的真實值, 其軌道根數(shù)、射程及交點的高程列于表1. 本文中的射程是指慣性空間中, 導(dǎo)彈的軌道在地球外部那一部分的弧長. 使用的根數(shù)格式是衛(wèi)星定軌中經(jīng)常采用的第一類無奇點根數(shù):除a、i和?外, 另外的三個根數(shù)為ξ=ecosω、η=-esinω、λ=M+ω, 其中M為軌道平近點角.

表1 仿真實驗中使用的軌道根數(shù)、射程及發(fā)、落點高程數(shù)據(jù)Table 1 Orbital elements, range and elevation of launching and landing points used in simulation experiments

在實際應(yīng)用中, 目標身份鑒別及預(yù)警信息生成是雷達在跟蹤目標時進行實時處理的一部分, 因此模擬雷達實時處理過程, 從中考察該方法的計算效率和計算精度的變化情況是最好的選擇. 我們設(shè)計實驗對表1的這5組射程和高程各異的根數(shù)均進行了仿真計算. 參考文獻[6]的方法, 首先根據(jù)彈道根數(shù)及發(fā)、落點選擇合適的觀測站址, 其次生成使用該觀測站來跟蹤目標的仿真觀測數(shù)據(jù), 然后在這些仿真觀測數(shù)據(jù)中加入角度0.1?、距離50 m的隨機誤差. 之后使用實驗室已有的程序?qū)⑶?0 s的數(shù)據(jù)進行濾波初始化, 最后使用卡爾曼濾波程序?qū)?0 s后的每秒一點的仿真數(shù)據(jù)進行濾波來更新目標的根數(shù).

每次濾波更新根數(shù)后, 我們分別采用原方法和本文提出的新方法產(chǎn)生目標兩個交點的位置和時刻, 觀察兩種方法計算出的兩個交點時刻和位置的計算誤差變化. 其中原方法的收斂精度設(shè)為0.01 m,新方法的高程迭代收斂精度ε2設(shè)為1 m, 最大高程迭代次數(shù)kmax設(shè)為5次, 交點時刻位置誤差的計算方法為:

其中tc和ts分別是時刻的計算值和真實值,→Rc和→Rs是地固坐標系下的位置矢量的計算值和真實值.

至于計算效率方面的變化, 新方法的計算量變化主要來源于高程迭代所增加的迭代計算, 因此在新方法計算效率的評價方面, 主要觀察增加的高程迭代次數(shù)的多少. 假定濾波過程中有l(wèi)次預(yù)警信息生成時進行了m（m= 1,2,···）次高程迭代,m的值可以用來評價新方法相比于原方法在計算效率上的變化. 事實上, DEM的讀取速度也是影響計算效率的一個重要因素, 但其涉及的是數(shù)據(jù)庫讀取的相關(guān)技術(shù), 已超出本文的討論范圍, 因此僅假定讀取DEM的速度與讀取內(nèi)存變量速度一致, 不做進一步展開.

表2列出了各算例在生成預(yù)警信息時增加的高程迭代次數(shù)的情況, 由表2可以看出, 實驗的5個算例的發(fā)、落點各自都進行了270次預(yù)警信息生成計算. 其中大部分情況下高程迭代2次即可收斂, 僅算例1和算例2的發(fā)點推算時, 超過2次高程迭代的預(yù)警信息生成次數(shù)較多, 甚至未收斂的生成次數(shù)分別為72次與99次, 這是由于這兩個發(fā)點處于高程高、坡度大的山地中. 算例1的落點位于玻利維亞的烏尤尼鹽沼, 此處高程變化很小, 因此雖然該點高程達3799 m, 卻仍然能保證大體上迭代2次即收斂. 至于發(fā)、落點高程都很低的算例5, 全程高程迭代基本1次收斂.

表2 新方法生成預(yù)警信息時增加的高程迭代次數(shù)的情況Table 2 Increase of elevation iterations when the new method generates warning information

圖1、3、5、7、9展示了5個算例的彈道目標發(fā)、落點位置誤差隨雷達跟蹤時間的變化情況,圖2、4、6、8、10展示了5個算例彈道目標發(fā)、落點時刻誤差隨雷達跟蹤時間的變化情況. 由圖1至圖10可以看出, 在濾波初始階段, 由于定軌誤差較大, 因此新方法和原方法計算的時刻和位置誤差均很大, 兩種方法之間并無明顯差異. 交點時刻和位置誤差最大是算例1, 分別達上百秒和五百公里量級, 誤差最小的算例5也分別達到了十幾秒和百公里量級; 濾波過程進行了3.5–4 min后, 交點的時刻和位置誤差均保持穩(wěn)定, 此時可以看到, 對于前四個算例以及第五個算例的發(fā)點, 新方法計算得出的交點時刻和交點位置精度都好于原方法. 交點位置的高程值越大, 新方法帶來的精度提升就越明顯. 對于算例1的彈道交點時刻誤差, 新方法的精度為0.5 s左右, 交點位置誤差為3.8 km左右, 而原方法給出的兩種誤差分別為2 s左右和10 km左右.對于算例2–4, 新方法的交點位置和交點時刻誤差分別為百米量級和0.01 s左右, 均優(yōu)于原方法公里量級和0.02 s左右的誤差. 對于發(fā)、落點位置接近參考橢球的算例5, 兩種方法推算的交點位置和交點時刻差距很小, 對于高程75 m的發(fā)點, 新方法精度僅略微優(yōu)于原方法, 而對于高程僅-38 m的落點,圖9和圖10顯示兩種方法的計算結(jié)果曲線重合在一起, 意味著兩種方法的推算精度幾乎一致.

圖1 兩種方法對算例1的發(fā)、落點位置計算誤差隨跟蹤時間的變化情況Fig.1 The variation of the calculation error with the tracking time for the launch and landing distance of example 1 by the two methods

圖2 兩種方法對算例1的發(fā)、落點時刻計算誤差隨跟蹤時間的變化情況Fig.2 The variation of the calculation error with the tracking time for the launch and landing time of example 1 by the two methods

圖3 兩種方法對算例2的發(fā)、落點位置計算誤差隨跟蹤時間的變化情況Fig.3 The variation of the calculation error with the tracking time for the launch and landing distance of example 2 by the two methods

圖4 兩種方法對算例2的發(fā)、落點時刻計算誤差隨跟蹤時間的變化情況Fig.4 The variation of the calculation error with the tracking time for the launch and landing time of example 2 by the two methods

圖5 兩種方法對算例3的發(fā)、落點位置計算誤差隨跟蹤時間的變化情況Fig.5 The variation of the calculation error with the tracking time for the launch and landing distance of example 3 by the two methods

圖6 兩種方法對算例3的發(fā)、落點時刻計算誤差隨跟蹤時間的變化情況Fig.6 The variation of the calculation error with the tracking time for the launch and landing time of example 3 by the two methods

圖7 兩種方法對算例4的發(fā)、落點位置計算誤差隨跟蹤時間的變化情況Fig.7 The variation of the calculation error with the tracking time for the launch and landing distance of example 4 by the two methods

圖9 兩種方法對算例5的發(fā)、落點位置計算誤差隨跟蹤時間的變化情況Fig.9 The variation of the calculation error with the tracking time for the launch and landing distance of example 5 by the two methods

以上計算結(jié)果表明, 交點位置高程值越大, 使用新方法相較于原方法提升的交點位置和交點時刻的計算精度就越多. 對于交點高程很小的彈道,兩種方法的推算結(jié)果差異并不明顯. 在本實驗中的大部分情況下, 每次生成預(yù)警信息時也僅增加1–2次高程迭代的計算量. 在交點高程高、坡度大時,迭代次數(shù)會增加, 未收斂的情況會增多, 但是交點位置和時刻的精度依然有大幅提升. 可見新方法加入高程迭代后, 可在增加少量迭代次數(shù)的情況下,極大程度地消除高程帶來的影響.相比于其他算例,算例1采用新方法的計算精度也偏低. 這是因為其主要誤差來源是較長射程導(dǎo)致的更久的外推時間,進而使得測軌誤差有更長時間的放大. 而這一部分的計算誤差顯然是無法通過高程迭代消除的.

5 結(jié)論和討論

本文基于文獻[6]的研究提出了一種有高程信息庫DEM支持的預(yù)警信息生成方法, 仿真實驗結(jié)果表明, 該方法具有下列特點:

（1）該方法相較于原方法僅增加了少量高程迭代計算, 且對不能收斂的情況有截斷和優(yōu)選措施,因此能保持較高的計算速度、穩(wěn)健和穩(wěn)定度, 能夠快速提供彈道導(dǎo)彈預(yù)警信息, 以供導(dǎo)彈防御決策;

（2）在雷達跟蹤彈道導(dǎo)彈的實時處理過程中,該方法計算的預(yù)警信息相比原方法能夠提升交點位置和交點時刻的預(yù)報精度.交點位置的高程值越大,精度提升越明顯. 通常在濾波進行了3.5–4 min后,能夠提供比較準確的發(fā)、落點信息, 發(fā)、落點位置和時刻的提升幅度可分別達到6 km和1.5 s. 對于發(fā)、落點高程值越小的彈道目標, 其發(fā)、落點的計算精度的提升效果會降低.

此外, 本文設(shè)計的仿真實驗顯示, 在設(shè)備觀測的前2–3 min, 原方法和新方法計算的精度都不高,甚至出現(xiàn)了新方法計算的預(yù)警信息精度低于原方法的情況. 因此在實際應(yīng)用中, 可以考慮在濾波達到穩(wěn)態(tài)之前僅使用原方法進行目標身份鑒別和生成預(yù)警信息, 待濾波達到穩(wěn)態(tài)后再加入高程迭代計算, 這樣就可以在不影響總體預(yù)警信息精度的前提下, 減少多次高程迭代計算以及DEM數(shù)據(jù)庫讀取操作, 從而在觀測設(shè)備跟蹤過程中進一步提高計算速度和效率.

需要說明的是, 在實際情況下目標在主動段的受力過程難以獲知, 因此即使使用本文的方法計算的發(fā)點信息仍然會有較大誤差. 落點計算則不同,設(shè)計彈道一般要求在不考慮大氣動力影響的情況下準確到達落點. 在此階段的控制調(diào)整主要任務(wù)是使目標能準確擊中指定地點, 因此, 采用本文的方法來提升落點位置時刻的精度會有較大的應(yīng)用價值.