空間目標雷達軌道改進原理及工程實現

黃曉斌， 張 燕， 魯 力， 肖 銳

(空軍預警學院， 湖北武漢 430019)

0 引言

雷達的優點是白天黑夜均能探測遠距離的目標，且不受霧、云和雨的阻擋，具有全天候、全天時的特點，并有一定的穿透能力。雷達已廣泛應用于社會經濟發展(如氣象預報、資源探測、環境監測等)和科學研究(天體研究、大氣物理、電離層結構研究等)中，但其應用最廣的還是軍事領域，如防空警戒、反導預警、空間監視、指揮引導等。

空間目標監視雷達是空間監視領域的骨干裝備，與常規防空雷達相比，主要增加了空間目標軌道確定功能，軌道改進實現該功能的核心模塊，軌道改進模塊的研制涉及非常復雜的航天動力學知識，這對從事雷達系統設計的非軌道力學專業的研究人員來說帶來極大困難。雖有一些開源定軌軟件可供借鑒，但這些軟件主要用于處理光學觀測數據，并不適合處理雷達觀測數據。為此，作者著手開發RadarOrbDet函數庫來輔助非軌道力學專業的人員進行空間目標監視雷達的系統設計，目前已完成坐標轉換模塊、初軌確定模塊和攝動分析模塊的開發，本文研究軌道改進模塊的開發。

軌道改進是空間目標軌道確定的重點和難點，它涉及的內容較多，包括攝動分析、偏微分方程的變分法分析、最小二乘優化等理論，本文聚焦空間目標雷達定軌中的批處理軌道改進問題，著重介紹了最小二乘軌道改進原理，講解了模型中的偏導數求解方法，描述了程序實現流程，最后用ODTK軟件驗證了本文軟件設計的有效性。

1 最小二乘軌道改進原理

最小二乘軌道確定的基本原理是找到一條軌道，使得理論觀測值和實際觀測值之間的殘差平方和(即損耗函數)最小，即求解，使得式(1)的值最小。

()=[-()][-()]

(1)

式中:=[();()]是衛星在時刻的位置和速度構成的狀態參量，由通過軌道預報可以得到一條參考軌道；=[;;…;]是次雷達實際觀測矢量，每次觀測矢量包含距離、方位和俯仰三個分量(,,)；()是由參考軌道計算得到的理論觀測值。

由于是一個非線性函數，式(1)描述的是一個非線性最小二乘問題，求解異常困難，通常采用泰勒展開將其轉化為一個線性最小二乘問題，然后通過迭代方式求解它的線性最小二乘近似解，表達式為

(2)

(3)

2 模型中的偏導數求解

式(2)的求解關鍵在于偏導數雅克比矩陣的計算，不失一般性，以某一時刻的雷達觀測量為例。為書寫方便，除時刻外，省略其他時刻索引。

利用求導鏈式法則，轉換為

(4)

式中，矩陣是當前時刻觀測量對當前狀態參量的偏導數矩陣，是當前時刻狀態參量對初始時刻狀態參量的偏導數矩陣，也稱為狀態誤差傳遞矩陣。

2.1 觀測偏導數矩陣求解

在忽略光行差等因素后，雷達觀測量只與衛星當前的位置有關，而與其速度無關。因此，矩陣可拆分成

(5)

式中，是衛星在慣性坐標系(ECI坐標系)下的位置矢量。

雷達觀測是基于站心地平坐標進行的(如南-東-天坐標系，SEZ)，直接計算觀測量對ECI坐標系下的位置矢量的偏導數比較困難，因此，矩陣的計算需進一步運用鏈式求導法則進行分解：

(6)

式中，是衛星在SEZ坐標系下的直角坐標，用(,,)表示；是雷達對衛星的觀測矢量在地固坐標系(ECEF坐標系)下的表示；是衛星在ECEF坐標系的位置矢量。

注意到

=-

(7)

式中是雷達在ECEF坐標系下的位置矢量，是個常矢量。

因此

(8)

根據坐標變換有

(9)

根據式(9)有

(10)

同理有

(11)

將式(8)、(10)和(11)代入式(6)有

(12)

在坐標系下，由直角坐標和極坐標的關系，易得

(13)

將式(13)代入式(12)，有

(14)

計算式(14)中的兩個矩陣時，都需要在當前觀測時刻進行。最后，將式(14)代入式(5)后，可得矩陣。

2.2 狀態誤差傳遞矩陣求解

設衛星的運動方程為

(15)

則

(16)

其中，矩陣為

(17)

由于采用迭代方式進行最小二乘問題求解，因此對矩陣的求解精度要求不高，從而對矩陣的求解精度也要求不高。因此，矩陣中的加速度可以只考慮地球中心引力和低階的非球形引力攝動，以加快計算機求解速度。

在只考慮地球引力(包含非球形攝動)情況下，加速度只與位置有關，而與速度無關，并且它的計算通常是在ECEF坐標系下進行，具體表達式可以參考文獻[18]，篇幅所限，這里不再贅述。

最后還需將ECEF坐標系下的偏導數矩陣轉換到ECI坐標系下。在忽略科里奧利力和離心力情況下，有如下關系式：

(18)

計算得到當前時刻的矩陣后，采用數值微分方程方法求解。

3 程序設計

圖1是最小二乘軌道改進算法的程序設計流程圖。

圖1 最小二乘軌道改進算法的程序設計流程圖

程序分為三層循環，第一層循環是數值微分方程求解的內部循環，輸出結果是第時刻的參考軌道和狀態誤差傳遞矩陣；第二層循環是遍歷每個觀測值，根據最小二乘原理求解單次軌道改進量；第三層循環是執行多次最小二乘軌道改進，使最終結果趨近于實際軌道狀態。

4 仿真實驗

4.1 接口調用示例

采用C語言按照第3部分的程序設計原理實現了最小二乘軌道改進的接口函數rodLeastSquare，并集成到RadarOrbDet庫中。rodAccelHarmonic接口函數的軟件調用方法如圖2所示。

圖2 最小二乘軌道改進接口函數說明

4.2 實驗比較

在STK11.2中仿真一顆衛星，其歷元軌道參數設置如表1所示，運動模式設置為HPOP(高精度軌道預報)，考慮21階非球形攝動。雷達站址的經緯高為(120°,30°,0)。選擇2021-09-07 19:15:01至2021-09-07 19:15:10間隔1秒的10點數據，對距離、方位和俯仰分別加上100 m、01°和01°的觀測誤差。

表1 衛星軌道根數

采用本文軟件和AGI公司開發的一款先進的軌道確定與分析軟件——Orbit Determination Tool Kit(ODTK)進行定軌比較，ODTK與本文軟件計算結果如表2所示，采用式(19)所描述的位置和速度相對誤差進行定量比較。

(19)

從表2可以看出，定軌的位置誤差為0.059 6%，速度誤差為0.35%，表明本文軟件算法是有效的，其誤差主要來源為計算程序的坐標與時間轉換、計算截斷誤差等。

表2 定軌結果對比表(ODTK與本文軟件)

5 結束語

本文針對空間目標雷達定軌中的批處理軌道改進問題，介紹了最小二乘軌道改進的原理，描述了模型中的偏導數求解，給出了程序設計思路，仿真表明了軟件設計的有效性。此外，本文開發的RadarOrbDet函數庫可以很好地輔助空間目標監視雷達系統設計人員進行定軌指標的快速設計和驗證。