一種光學短弧初軌確定的軌道偏心率判定方法

張鄭元, 李彬,2, 趙廣宇,4, 曾春平, 葉釗, 桑吉章,2

(1. 武漢大學測繪學院, 湖北 武漢, 430079; 2. 湖北珞珈實驗室, 湖北 武漢, 430079;(3. 航天東方紅衛星有限公司 北京 100094; 4. 中國科學院空天信息創新研究院 北京 100094)

1 引言

經典的空間目標初軌確定方法, 如Gauss 法、Laplace 法, 最初是用來解算小行星的軌道。 當其應用于近地衛星軌道的光學初始軌道確定時, 往往不收斂或得不到精確解[1,2]。 隨著計算機技術的發展, 如double-r法、 Gooding 法的數值算法的出現, 解決了經典算法針對近地目標光學初軌確定時存在的不易收斂的問題[3-6]。 然而, 由于光學短弧對軌道的幾何約束能力很差, 這些數值算法在僅利用3 方向觀測值進行短弧初軌解算時,估計精度仍然難以達到應用期望[7,8]。

隨著光學觀測技術的發展, 光學觀測的數十秒至數分鐘弧段可能有數十至上百個觀測方向,因此可以組成許多3 方向組合。 Henderson 首先使用Gooding 方法解算多組3 方向組合, 然后通過遺傳算法挑選出比較好的解算結果取平均作為最終解[9]。 這種做法增強了Gooding 法在光學短弧初軌解算時的魯棒性, 但由于即使在小偏心率初軌解算時, Gooding 初軌參數集的分布也并非高斯分布, 而在偏心軌道時, 軌道參數估計分布的均值與軌道參數真值之間往往存在顯著偏差, 因此最終的初軌參數經常也存在很大的誤差。

在處理短弧多方向初軌確定問題時, 章品對弧段的首尾方向的斜距進行網格搜索, 通過由假設斜距所得到的計算觀測值與真實觀測值之間的差異大小, 來確定多組軌道解, 然后用小偏心率約束, 從多組軌道解中挑選最終解[10]。 這種方法會導致偏心軌道解算結果的誤差較大。 因此, 目前迫切需要一種更準確的初軌確定方法, 其應充分利用短弧光學弧段的信息, 且能夠自適應各種偏心率下的軌道。

有鑒于此, 本文將提出一種短弧多方向初軌參數確定方法。 首先通過一定策略形成多個3 方向組合, 并用經典的Gooding 方法解算每個3 方向組合的初軌參數, 如果收斂, 則立刻得到一個候選軌道解。 然后, 利用眾多待選軌道解的“半長軸-偏心率” 的二維分布特征, 初步確定目標軌道的偏心率, 最后根據初定偏心率從待選軌道解集合中確定最終解。 本文注意到, Gooding 方法在進行光學短弧初軌確定時需要較為精確的初始距離估值, 否則容易造成不收斂的問題。 因此,本文也對距離搜索法的小偏心率約束進行了優化, 使其能夠在偏心軌道上也能提供準確的距離。

本文章節安排如下, 第2 節介紹一種改進的距離搜索法, 可以緩解該算法的小偏心率約束準則在偏心軌道上出現大估計偏差的情況。 第3 節討論Gooding 方法解算3 方向組合所得到待選解半長軸-偏心率的分布與目標軌道真偏心率之間的關系, 并給出一種選解策略。 第4 節開展大規模的仿真實驗驗證方法的有效性。 第5 節對本文研究進行總結。

2 改進的距離搜索法

假設在觀測時刻ti, 傳感器與目標之間的關系可用方程(1) 表示, 即:

式中,αi和δi分別表示觀測時刻ti對應的“傳感器-目標” 之間的赤經和赤緯。

由于距離搜索法采用小偏心率約束準則, 其在小偏心率軌道上具有良好的斜距初值估計精度。 但當其應用于偏心軌道的初始斜距估計時,會出現較大的偏差。 Gooding 方法在使用距離搜索法所提供的斜距初值進行初軌確定時, 在偏心軌道上極易收斂到錯誤解或不收斂。 因此, 本文根據距離搜索法待選解的偏心率分布特征提出一種修正方法, 以較好估計目標的偏心率大小, 并采用合適的策略選擇最終解。

距離搜索法對弧段首尾時刻“傳感器- 目標” 斜距ρfore和ρaft進行網格搜索, 通過由假設斜距所得到的計算觀測值與真實觀測值之間的差異大小, 判別并得到多個待選軌道解[10]。 其中, 一組首尾時刻斜距的網格搜索組合見下式:

由一組可選變量xjk, 便可計算首尾時刻目標對應的目標位置矢量r?fore和r?aft, 代入Gauss-Lambert 方程, 便可以得到一組軌道根數解[3]。通過軌道根數解和二體方程, 并考慮攝動項, 便可以計算得到目標在觀測時刻ti的位置矢量。

那么, 用統計量來描述這些觀測資料之間的偏差, 取統計量都滿足閾值σn時所對應的xjk計算得到待選的軌道參數解, 待選的軌道解個數表示為Nn, 其中n表示距離搜索法的循環次數。

需要注意的是, 解的個數N0與判別閾值σ0的選擇有極大的關系。 低軌目標相比于高軌目標而言, 需要更大的閾值, 才能保證相近的數量級。 其原因在于, 在相同的斜距網格搜索步長下, 當斜距發生一個搜索網格的變化時, 低軌目標所產生的角度觀測值變化更大。

為保證不同類型軌道解的個數都能夠滿足后續軌道特征的分析, 本文采取以下優化措施。 先使用大閾值σ0來確定待選解的數目N0, 如果N0過高, 就采用較小的閾值σ1再重新得到一組新待選解數目N1。 重復此步驟, 直至閾值σn使得待選解數目Nn滿足需求。 一般Nn處于200 ～500 之間較好,n的次數要小于3。

當得到距離搜索法的待選解后, 需要一個合適的選解策略。 原始的距離搜索法一般選擇偏心率最小的一批待選解, 然后計算它們的平均值作為最終解輸出。 為了能夠詳細地說明這種策略,本文假設了一個900km 高的太陽同步軌道觀測平臺, 該軌道高度基于2000 國家大地坐標系所采用的地球半長軸參數, 即6378.137km, 并使其觀測兩個仿真目標。 其中第一個目標的半長軸為7500km, 偏心率為0.001, 另外一個目標的半長軸為9000km, 其偏心率為0.1。 目標-1 的弧長設計為120s, 目標-2 的弧長設計為240s, 角度觀測誤差為2″, 只考慮J2長期項對軌道的影響。觀測平臺和兩個目標的軌道根數見表1。

表1 天基觀測平臺和仿真目標的軌道根數Table 1 Orbital elements of the space-based observation platform and the simulation target

兩個目標分別使用前述的閾值策略得到待選解, 這些待選解的“半長軸-偏心率” 分布見圖1。 需要注意, 所有待選解的軌道根數都以弧段第一個觀測歷元為參考時刻。

圖1 距離搜索法下待選解的半長軸-偏心率分布圖Fig.1 Semi-major axis-eccentricity distribution of the solution to be selected under the range search method

由圖1 可見, 無論目標軌道是否為近圓軌道,待選解的半長軸-偏心率分布都呈現多條曲線平行或交錯的形式。 顯然, 基于半長軸-偏心率分布并不能有效地確定目標軌道的偏心率, 這會給最終解的確定帶來影響。 而且, 這對于偏心軌道的影響更大一些。

不過, 分析待選軌道解的偏心率分布可以發現, 偏心率的頻率分布并不是標準的高斯分布,可用概率密度函數f(e) 表示。 如果按照一定的偏心率間隔(ei,ej) 來計算待選解偏心率出現的頻率直方圖, 目標軌道的真實偏心率所對應的直方圖頻率是最高的, 其中i、j由算法設計的偏心率間隔數目確定。 此過程可用式(5) 表示:

式中,emax表示最高的直方圖頻率所對應的偏心率區間, 可用“max-間隔” 表述。

以前述兩個目標的待選解偏心率分布直方圖為例(見圖2), 圖中的偏心率間隔為0.05。

圖2 距離搜索法下待選解偏心率的頻率分布直方圖Fig.2 The frequency distribution histogram of the deeccentricity to be selected under the range search method

在圖2 (a) 中, 待選解偏心率處于0 ～0.05之間的頻率為0.22, 處于0.05 ～0.1 之間的頻率為0.18, 它們顯著大于其他偏心率區間的直方圖頻率。 在圖2 (b) 中, 待選解偏心率處于0.05 ～0.1之間的頻率為0.14, 而0.1 ～0.15 之間的頻率為0.23, 也顯著大于其他偏心率區間的頻率。 由此可見, 目標軌道的真偏心率就處于待選解偏心率最大頻率所對應的偏心率間隔 (即max-間隔)。

在此基礎上, 可對距離搜索法進行優化。 當max-間隔處于0 ～0.1 之間時, 表明目標軌道是近圓軌道, 此時可用小偏心率約束策略來選擇最終的初軌解。 當max -間隔的最小值大于0.1 時,表明目標軌道是偏心軌道。 此時, 可以把滿足max -間隔的前后兩個間隔的待選解組合成一個新的待選解集合, 然后在這個集合中選擇半長軸和偏心率的中位值作為最終初軌解的半長軸和偏心率。

用前述改正算法和原始距離搜索法對目標-1/目標-2 的測角弧段進行200 次初軌解算, 每次解算時的角度誤差隨機添加。 初軌最終解的半長軸和首尾斜距估計誤差的均方根誤差(Root-Mean-Square, RMS) 見表2。

表2 距離搜索法和改正算法的半長軸和首尾斜距估計RMSTable 2 Semi-major axis and fore-tail slant distance estimate RMS by range search method and correction algorithm

由表2 可知, 由于使用小偏心率約束策略,原始的距離搜索法解算目標-1 的軌道半長軸、首斜距、 尾斜距的估計誤差分別為102.4km、117.3km、 127.1km, 要優于目標-2 的428.7km、344.0km、 270.7km。 利用優化算法, 兩個目標的軌道半長軸/斜距的估計精度都得到了較大提升。改正算法對于目標-1 的軌道半長軸/首斜距/尾斜距的估計精度為49.4km、 62.7km、 72.9km,對于目標-2 的軌道半長軸/首斜距/尾斜距的估計精度為175.9km、 138.8km、 103.5km, 都高于原始距離搜索法的估計精度。 以上實驗說明, 利用改進的距離搜索法, 可以為Gooding 法提供更加精確的斜距估計初值。

3 軌道偏心率判定方法

基于改進的距離搜索法所提供的斜距初值, 就可以基本確保Gooding 方法利用3 方向觀測值解算初軌參數的收斂。 現假設觀測時刻為ti(i =1,2,3) , 距離搜索法所提供的首尾時刻斜距值為ρ1、ρ3。 那么, 由式(2) 和Gooding-Lambert 方程, 可以得到目標位于中間ti時刻的一個解[4]。 若ρ1、ρ3為真值, 則與重合, 否則兩者之間有偏差, 該偏差可用圖3 表示[5-6]。

圖3 目標函數示意圖Fig.3 Objective function diagram

基于圖3 中的幾何關系, Gooding 方法用f和g表示偏差, 并作為首尾時刻斜距ρ1、ρ3的目標函數, 同時將ρ1、ρ3作為自變量, 并用x、y表示。 那么, Gooding 方法的函數關系可以表示為:

使用Newton-Raphson 迭代法, 自變量的迭代值應為:

式中,δx、δy為自變量的迭代值,fx、fy、gx、gy分別為目標函數f、g對自變量的x、y的偏導數, 可用數值微分方法得到。

經過式(6) -(7) 的反復迭代, 并結合Gooding-Lambert 方程, 即可得到一組由3 方向觀測值得到的初軌解。 而對于一個具有數十甚至上百方向觀測值的光學測軌弧段, 可以組成許多3 方向組合。

本文將一個光學弧段分為首部、 中部、 尾部三部分。 首部和尾部的弧段時長占總弧段時長的25%, 中部的弧段時長占50%。 每個3 方向組合的第一個方向來自首部, 第二個來自中部, 第三個來自尾部。 選擇合適的方向組合, 保證Gooding方法所解算的待選解總數在200 ～500 之間, 和上節所討論的距離搜索法的待選解數目一致。 同時, 將待選解數目在10 以下的弧段舍棄, 因為此時不能從待選解的分布中識別出有效的半長軸-偏心率分布特征。

與距離搜索法類似, 此處同樣存在如何從待選解集中確定最終解的問題, 其中的關鍵是判定目標軌道的偏心率。 為此, 我們分析待選解半長軸-偏心率的二維分布特征。 設表1 所列的天基觀測平臺分別觀測了表3 所列的6 個目標。 真軌道生成時只考慮了J2長期項攝動, 角度觀測誤差為2″, 觀測頻率為1Hz。 這些目標的Gooding 法待選解的半長軸-偏心率分布見圖4 和圖5。

圖4 類型A 的半長軸-偏心率分布圖Fig.4 Semi-major axis-eccentricity distribution of type A

圖5 類型B 和類型C 的半長軸-偏心率分布圖Fig.5 Semi-major axis-eccentricity distribution for type B and type C

表3 TLE 目標的軌道根數Table 3 TLE Number of orbital elements of the target

圖4 顯示了Gooding 方法待選解最常見的半長軸-偏心率分布, 約占第三節實驗弧段總數的92%。 與距離搜索法的多條半長軸-偏心率曲線不同, 盡管分布密度不均勻或線斷裂, 但每個分布都呈現為光滑連續的折線或曲線。 對于這種類型的分布, 本文將其統稱為類型-A, 從這種類型的待選解中, 來確定最終解相對比較簡單, 而且確定最終解的準確性通常很高。

類型A 的半長軸-偏心率分布可以用一個或兩個高次多項式方程進行擬合, 此時可假設待選解組合的偏心率為因變量ε, 半長軸為自變量α, 那么擬合方程為:

式中,ai(i =0,1,2) 為擬合多項式的系數。

在式(7) 基礎上, 還需要判斷擬合線段是否離散, 可用擬合優度指標R進行判定, 該指標的計算公式為:

其中,SSE表示殘差平方和, 即

式中,εi表示待選解的真實偏心率值,ε︿i為由式(7) 計算得到的回歸偏心率值。

其中,SSR表示回歸平方和, 即

圖4 (a) 為近圓軌道的分布。 它呈現為“V” 形, 底部的偏心率接近于0, 并且待選解的分布在底部附近密度最大。 這種V 型分布被稱為類型-A1, 在偏心率接近于零的軌道上是最典型的分布。 其一般可以擬合出兩條多項式方程, 且兩個多項式都滿足a2≈0,R ＞0.95 。 但隨著偏心率的增加, 它逐漸演變成一條帶弧度的曲線, 如圖4 (b) 所示。 將偏心率最小的待選解作為類型-A1 的最終解, 將會獲得較為準確的初軌參數估計值。 圖4 (b) 中的分布被記為類型-A2, 可以用二次多項式近似, 即在由式(7) 所擬合的多項式中,≥0,R ＞0.95。 真解位于曲線曲率最大的位置。 當軌道偏心率大于0.15 時, 分布曲線呈現為類型- A3 或類型-A4, 分別如圖4 (c)和4 (d) 所示。 其中, 類型-A3 最為常見, 其真解的半長軸和偏心率與待選解半長軸和偏心率的中位數接近。 其只能擬合出一條多項式方程,且滿足a2≈0,R ＞0.95 。 但是, 如果是分布更為廣泛的類型-A4 型, 則很難得出準確的最終解決方案。 類型-A4 的特征比較明顯, 在4 個亞型中半長軸和偏心率分布范圍最廣。 此時, 雖然可以擬合出多項式, 但是R ＜0.8 , 真解通常位于待選解分布的底部附近, 但仍然很難確定準確的最終解。 幸運的是, 類型-A4 僅占所有分布的4%左右。

圖5 (a) 顯示了半長軸-偏心率分布的另一種常見形狀, 稱為類型- B, 約占所有分布的6%。 在這種類型中, 待選解聚集在一個封閉的區域內。 這種類型的真解通常非常靠近區域的質心, 因此我們使用質心來確定最終解。

其余的半長軸-偏心率分布不太規則, 難以描述, 約占總數的2%。 這種類型稱為類型-C,如圖5 (b) 所示。 造成類型-C 半長軸-偏心分布的主要原因是距離搜索法的斜距初值估計不準確, 導致Gooding 方法收斂錯誤。 在這種情況下,沒有辦法確定最終解, 一般判定初軌解算失敗。

需要注意的是類型-B 和類型-C 的擬合優度指標R都較低, 需要加入專家判斷來輔助分類, 機器學習算法也可以提供一個初步的分類參考[11]。

綜上所述, 本方法通過從候選解集找到半長軸-偏心率分布中最密集的點來確定最終解。 因此, 如果半長軸-偏心率的分布范圍很廣, 則很難確定準確的最終解。

基于以上討論, 圖6 為根據半長軸-偏心率分布判定目標軌道偏心率的算法流程。

圖6 算法流程圖Fig.6 Algorithm flow chart

4 實驗結果

本文設計天基仿真實驗來評估上述最終解優選策略的有效性。 考慮到不同軌道高度目標與觀測平臺的相對速度的不同, 觀測場景和幾何特征大不相同, 因此本文區分了不同種類的軌道目標。 它們的軌道參數可見表4。

表4 不同軌道類型的軌道參數Table 4 Orbital parameters of different orbits

天基觀測平臺軌道信息見表1, 設計為900km 高的太陽同步軌道, 該軌道高度基于2000國家大地坐標系所采用的地球半長軸參數, 即6378.137km。 仿真目標的初始軌道信息來自于2023年6月1日-2023年6月3日的TLE 數據庫, 并轉化為TLE 根數對應時刻下的開普勒根數。 同時, 設計LEO1和LEO2的弧段長度為120s, MEO、 HEO 和GEO 目標的弧段長度為240s。 每個目標的軌道考慮J2的長期項影響,目標角度誤差為2″, 觀測頻率為1Hz。 最后捕捉到了9325 個目標, 假設每個目標只觀測到一個光學弧段。

將原始的Gooding 方法通過多組待選解求平均的策略, 記為算法-1。 同時, 將基于Gooding方法, 并使用優選策略進行軌道偏心率判定方法, 記為算法-2。 這兩種算法都使用改正的距離搜索法所提供的斜距初值進行初軌解算。 分別將兩種算法應用于各類型下的光學弧段初軌解算。需要注意的是, 由于優選算法只是針對于待選解進行優化的算法, 因此其收斂率與原始Gooding方法是相同的。 在本次實驗中都為94.7%, 即8834 弧段成功收斂。 兩個算法的半長軸與傾角誤差RMS 見表5。

表5 兩種算法的半長軸和傾角誤差RMSTable 5 The semi-major axis and inclination error RMS of the two algorithms

從表5 中可以看出, 相比于算法-1 而言,算法-2 半長軸、 軌道傾角的估計精度更加精確。而且相比于偏心軌道而言, 算法-2 在小偏心率軌道上的初軌估計精度的提高更顯著。 算法-2在LEO1、 MEO、 GEO 軌道上的半長軸估計精度為55.93km、 46.64km、 150.16km, 而算法-1 只有162.10km、 357.95km、 2153.03km, 說明新算法能夠通過待選解的半長軸-偏心率分布, 準確地區分出目標的軌道屬性, 并給出較準確的軌道估值。 算法-2 在LEO2、 HEO 軌道上的半長軸估計精度為86.90km、 419.47km, 相比于算法-1的120.16km、 467.53km, 也有較大提升, 也說明了新算法在偏心軌道上的作用。 圖7 為兩種算法在前述不同軌道類型下的半長軸、 軌道傾角的誤差分布箱線圖。

圖7 天基仿真情景下Gooding 法和I-Gooding 法的誤差箱線圖Fig.7 Errorbox plots of Gooding method and I-Gooding method under the background of space-based simulation

從圖7 (a) 中可以看出, 相比于算法-1 的半長軸誤差分布而言, 算法-2 在小偏心率目標上誤差分布更加集中。 這是由于新算法能夠根據Gooding 算法待選解軌道半長軸-偏心率分布, 準確地識別小偏心率軌道, 并據此給出精確解。 同時也可以看出, 算法-2 在偏心軌道上的誤差分布相比于算法-1 更集中。 從圖7 (b) 可以看出, 算法-2 的軌道傾角誤差分布相比于算法-1更靠近0, 說明新算法在軌道傾角的估計中也能獲得較高的精度。

以上實驗結果說明, 根據Gooding 法待選解的半長軸-偏心率分布, 本方法可準確識別目標軌道的偏心率, 并得到較精確的初軌。 這給今后的光學觀測初軌確定工作提供了新的思路。

5 結論

隨著光學觀測技術的發展, 一個數十秒至數分鐘長的短弧一般有數十甚至上百個觀測方向。為使經典的Gooding 初軌確定方法能夠有效地處理光學多向觀測弧段, 本文主要從以下兩個方面對現有算法進行了改進。 首先, 對于只能應用于小偏心率初軌計算的距離搜索法, 按照其待選軌道解的偏心率分布特征, 準確識別出目標軌道的偏心率, 從而使距離搜索法應用于偏心軌道時也可得到較為精確的初軌參數與斜距。 其次, 根據Gooding 法的3 方向待選解半長軸-偏心率分布,能夠準確識別目標軌道的偏心率, 并能夠根據這些特征選擇合適的初軌解。

利用天基仿真光學觀測數據的實驗驗證了改進距離搜索法和優選算法的有效性。 實驗結果表明, 對于不同軌道類型的光學短弧初軌確定, 所提出的軌道偏心率判定方法均可提升初軌估計結果的精度。 新算法在近圓軌道初軌解算中表現優異, 對于LEO1、 MEO、 GEO 短弧初軌確定, 半長軸誤差的 RMS 值為 55.93km、 46.64km、150.16km, 傾角誤差的RMS 為0.18°、 0.05°、0.14°。 對于偏心軌道LEO2、 HEO, 應用新算法的半長軸誤差RMS 為86.90km、 419.47km, 傾角誤差RMS 為0.16°、 0.90°。

因此, 相信本文所提出的空間目標光學短弧軌道偏心率判定方法有助于光學空間監視中的新目標識別和關聯編目, 從而為自主編目能力提升做出貢獻。