近地天體與空間目標初軌確定的多解問題

趙柯昕, 甘慶波,*, 楊志濤, 劉 靜

(1. 中國科學院國家天文臺, 北京 100012; 2. 國家航天局空間碎片監測與應用中心, 北京 100012;3. 中國科學院大學, 北京 100049)

0 引 言

利用測角數據進行初軌確定的問題從提出至今,已經有200多年的歷史。在此期間,提出了許多方法來解決該問題。最早提出的方法是Laplace法和Gauss法。這兩種方法均使用三次測角數據，且至今仍有廣泛應用。Laplace法通過對觀測幾何構型進行微分,構建了關于中間觀測時刻目標地心距離的八次方程。Gauss法通過假設三次觀測時刻的目標地心位置矢量位于同一平面,同樣構建了關于中間觀測時刻目標地心距離的八次方程。雖然兩種方法的八次方程具有相同形式,但系數的具體定義不同。而在初軌確定的過程中,求解這兩種形式的八次方程都需要面對從多個根選擇出符合觀測數據的正確解的問題,若無法解決該問題，會導致后續軌道參數計算的失敗。由于計算機技術的發展,Escobal于1965年提出了雙R迭代法,但是該方法的魯棒性較差。Gooding于1993年提出了基于Lambert方程的初軌計算方法,在將其應用于天基空間目標監測時也相當穩健。Karimi于2011年提出了一種使用拉格朗日系數的迭代計算方法,該方法在共面情況下不存在奇異性。桑吉章等人將利用測角數據的初軌確定問題轉換為基于兩個位置矢量的初軌確定問題,使用仿真和實測數據發現該方法的成功率和精度都優于Gauss法和Gooding法,但該方法對于初值選擇較為苛刻。Wishnek等人于2021年提出了一種使用角度測量值的允許域的初軌確定新方法,發現該方法魯棒性較好。Gronchi等人推廣了Mossotti提出的利用多次觀測和線性方程進行初軌確定的方法,將其應用于地球衛星的軌道確定中。隨著優化算法和神經網絡的興起,Ansalone、李鑫冉、Hu等人將遺傳算法、粒子群算法、進化算法等方法應用于地基、天基空間目標初軌確定問題和小行星初軌確定問題，對于近圓軌道獲得了較好的效果。Handley等人提出了一種基于自適應Nelder-Mead最優化方法的天基空間目標確定方法,尤其在高軌對高軌目標時有較強性能。

針對Gauss法和Laplace法在初軌確定過程中求解八次方程時遇到的根的個數問題，已經有了較為詳細的研究,但是對于根的選擇問題研究較少。Charlier針對近地天體定軌問題,對觀測目標與太陽的距離量和觀測目標與地球之間的斜距量進行歸一化，重構了八次方程,發現八次方程共有8個根,其中包含3個正實根,1個負實根和4個虛根。他將非偽解定義為使得斜距量具有意義的正實根。同時給出了判別非偽解個數的幾何條件,認為非偽解的個數只取決于中心天體、觀測平臺和空間目標三者在某一時刻平面上的相對位置。1962年,Danby等人使用正弦定理重新表示了觀測幾何構型,將求解八次方程根的過程轉化為求解兩個正弦曲線交點的過程,但是該方法仍面臨多個非偽解無法選擇的問題。2009年,Gronchi將Charlier的理論進行了推廣,提出了Charlier極限概念,但沒有考慮多解問題。2012年,Der將可見的空間目標應用于觀測平臺所處的地平面之上,作為判別條件。若計算得到的根對應的觀測視線穿過地球,則該根為錯解予以舍棄,但對于有些觀測數據，仍需要其他信息才能進行非偽解的選擇。2016年,Wie將八次方程分為兩種情況進行簡單近似,認為可以利用近似形式對非偽解進行選擇,顯然這是一個不準確的情況。

本文以Laplace法八次方程為例,討論了非偽解個數與觀測幾何構型之間的關系,并提出了一種有效的對非偽解的判別方法。首先,本文構建了Laplace法的八次方程,討論了方程系數與根的個數之間的關系。利用觀測平臺地心距對經典八次方程進行歸一化,解決了平凡解的問題,并分析了非偽解的個數與地心、觀測平臺和空間目標三者相對位置的關系。隨后從計算半通徑的問題入手,引入了Shefer提出的通用方程,推導出了空間目標斜距量應滿足的約束條件,給出了空間目標斜距變化量與觀測時間間隔之間應滿足的非線性方程組的解析形式。利用該方程組作為判別方程,對多個非偽解進行選擇。最后，使用近地天體和天基衛星光學仿真數據驗證了本文提出的非偽解個數的確定方法與判別方法的正確性。

1 光學資料初軌確定問題

利用三次測角數據進行初軌確定是天體力學中的基本問題。本節使用3次角度觀測值分析了Laplace法中八次方程面臨的非偽解的個數與性質。首先定義3次觀測時刻與對應的赤經、赤緯測量值分別為,,,=1,2,3,觀測時刻的觀測平臺地心位置矢量為(=1,2,3),需要求解出空間目標的地心位置矢量(=1,2,3)。

1.1 光學觀測幾何構型

使用天基和地基光學觀測平臺對空間目標和近地天體進行初軌確定,單次觀測幾何構型如圖1所示。觀測近地天體時,可認為觀測平臺位于地心,同時觀測近地天體和天基空間目標觀測的動力學是類似的,后續僅考慮利用地基和天基光學觀測平臺對空間目標進行初軌確定的情況。

圖1 觀測空間目標和近地天體的幾何構型Fig.1 Geometrical configuration of observation of space object and near-Earth object

由圖1可知,空間目標的地心位置矢量可以表示為

=·+,=1,2,3

(1)

式中:·為觀測平臺到空間目標的距離矢量,是觀測平臺到空間目標視線的單位矢量,是斜距;下標代表3次不同的觀測時刻。平臺到空間目標視線的單位向量可以表示為

(2)

1.2 Laplace法八次方程

將式(1)點乘自身,可以得到關于空間目標地心距離的表達式:

(3)

對式(1)進行一階、二階微分,可以得到每個觀測時刻的空間目標的速度和加速度:

(4)

(5)

時間間隔較短時,可以假設空間目標為二體運動。聯立式(1)～式(5),并經過轉換,可得到Laplace形式的八次方程:

()=+++=0

(6)

式中:

各符號定義如下:

針對方程,可以總結系數的性質如下:

(1) 系數和是恒負的;

(2) 系數可能為負,也可能為正;

(3) 當系數<0時,總會有1個正實根,1個負實根和6個復根;

(4) 當系數>0時,存在至多3個正實根,1個負實根和復根。

初軌確定過程中所需要的真解是方程的1個正實根。Vallado指出,該方程可能存在多個使得斜距量有意義的根,真解的選擇較為困難。需將每個根與先驗信息進行比較或者利用多次觀測值分離出正確的根,但該過程是繁瑣且效率低下的。

1.3 平凡解與非偽解的性質

方程中的系數取決于地心、觀測平臺和空間目標的相對位置。使用觀測平臺的地心距對空間目標地心距和斜距進行歸一化,式(6)的系數變為

(7)

可得,=1是方程的一個平凡解。當利用天基平臺觀測空間目標時,該根代表了觀測平臺的地心距;當觀測近地天體時,該根代表了地球的日心距。=1會導致()=0,所以該根是方程的平凡解,應該被舍棄。

求解經過降階的七次方程時,會出現至多2個使得斜距有意義的根,即至多存在2個非偽解。由Charlier提出的非偽解個數的條件,進一步可以得到存在兩個非偽解的判別方程:

(8)

如圖2所示,當空間目標位于陰影區域和時,會存在2個非偽解,而當空間目標位于區域和時,只存在一個非偽解。這4個區域的分界線為=和=+23-53。

圖2 非偽解個數與地心、觀測平臺和空間目標三者相對位置的關系Fig.2 Relationship between the number of non-pseudo solutions and the geometrical configuration of Geocenter, observation platform and space target

初軌確定的典型問題是近地空間目標和近地天體的觀測定軌。其中,通常的典型場景有以下3種:

(1) 地基近地空間目標(如衛星、碎片等)觀測,幾何構型如圖1(a)所示。可以證明在初軌確定過程中不存在平凡解情況,本文研究暫不涉及此觀測場景。

(2) 天基近地空間目標觀測,幾何構型如圖1(b)所示。觀測平臺繞地球作軌道運動,分界線=代表以觀測時刻觀測平臺所在地心距為半徑、以地心為原點的圓,可觀測到位于區域和中的物體,此時需要考慮多個非偽解判別的問題。

(3) 太陽系近地小天體觀測,幾何構型如圖1(c)所示。引力中心為太陽,觀測平臺位于地球表面,地球自轉可忽略,分界線=代表地球公轉軌道。其動力學本質和天基近地空間目標觀測近似,同樣可以觀測到位于區域和中的物體,也需要考慮多個非偽解判別的問題。

2 非偽解的判別方法

本節從求解軌道半通徑的問題出發,引入Shefer的改進通用方程,推導出了空間目標斜距量應滿足的約束條件,給出了空間目標斜距變化量與觀測時間間隔之間應滿足的解析形式的非線性方程組。利用該方程組作為非偽解的判別方程,對多個非偽解進行選擇。

2.1 半通徑計算方程

當觀測時間間隔較小時,可以認為空間目標三次地心位置矢量位于同一平面:

·(×)=

(9)

給定3個觀測時刻的空間目標相對于觀測平臺的斜距量,由此可以利用觀測幾何構型計算出觀測時刻的空間目標地心位置矢量(,,),從而可以計算軌道的半通徑:

(10)

式中:符號±和?對應于叉乘中兩個位置向量夾角正弦值的正負號;符號×表示向量叉乘。本文考慮的是觀測時間間隔較小的情況,和,和,和之間的夾角都小于π,因此式(10)可以化簡為

(11)

現在考慮由第1、第3個觀測時刻的空間目標地心位置矢量(,;,)求解半通徑的問題。由Kepler定律,兩個位置矢量可以構成一個三角形和一個扇形,如圖3所示。圖中第一、第三觀測時刻的空間目標所在位置分別為點與點，真近點角分別為與。兩時刻的目標位置矢量和線段可以組成一個扇形面積為，而與線段可以組成一個三角形面積為。

圖3 地心位置矢量組成的三角形和扇形Fig.3 Triangle and sector composed of geocentric position vectors

引入輔助變量表示扇形面積和三角形面積的比值:

(12)

當計算出面積比后,可以較為方便地求出半通徑。為了求解引入下式:

(13)

-=

(14)

式中:

將式(13)代入式(14):

=1+(+)

(15)

針對橢圓軌道,函數表示為

(16)

(17)

式中:d=-,和為第1、第3觀測時刻目標的偏近點角。

式(13)和式(14)需要迭代求解,Escobal于1976年提出了一種計算方法。當地心位置矢量(，)已知時,可計算得到參數和。首次循環中,令=1,依次求出參數和函數的值,由此更新的值。重復上述過程直到達到預定精度。對于較短弧段,即兩時刻地心位置矢量(,)的夾角較小時,具有很快的收斂速度。但當夾角較大或夾角等于π時,則無法計算面積比。

2.2 Shefer方程

為了避免計算過程中出現奇點,Shefer將式(13)改寫為

(18)

式中:

位置矢量和重合會導致參數等于0,但該情況在實際觀測時是幾乎不會發生的。因此,忽略=0的情況,聯立式(13)和式(14),新方程被稱為Shefer方程:

()()=

(19)

式中:

()=+()()

(20)

(21)

針對橢圓軌道,函數可表示為

(22)

2.3 斜距約束與非偽解判別方法

定義3個觀測時間間隔=-和=-,考慮時間間隔內空間目標位置矢量的變化與時間間隔之間的關系。將斜距(,)代入式(19),可得斜距量需要滿足的約束條件:

(23)

式中:

由式(1)和式(11)可得函數和分別為

(24)

針對橢圓軌道,函數()可以使用Lerch方程形式表示:

(25)

式中:

式(23)是空間目標斜距需要滿足的約束方程。初軌確定過程中,若得到了多個非偽解,需進一步計算得到每個非偽解對應的第1、第3觀測時刻的斜距(,),并代入約束方程式(23)中。若使得約束方程和接近于0,即約束方程在(,)附近存在交點,則該解為真解。否則，該解可當作錯解舍棄。

3 仿真驗證

3.1 近地天體光學仿真觀測數據驗證

使用阿波菲斯小行星的三次光學觀測的仿真數據進行驗證,仿真弧段長度為20天。觀測時刻使用簡約儒略日(modified julian date, MJD)表示。地球在J2000日心天球坐標系中的位置矢量由DE421星歷得到。地球的日心位置矢量單位為天文單位AU(1AU=1.496×10km)。角度測量值位于J2000地心天球坐標系下。日心、地球和小行星三者的相對位置幾何構型如圖4所示。地球位置矢量和光學仿真觀測數據如表1所示。

圖4 太陽、地球與阿波菲斯小行星的相對幾何構型Fig.4 Geometrical configuration of Sun, Earth and Apophis

表1 阿波菲斯小行星光學仿真觀測數據Table 1 Apophis optical simulation observation data

使用表1的數據構建Laplace形式的八次方程:

()=-2007 872 21+1241 891 36-0249 875 22

利用Aberth法求解該方程,此時得到的8個根分別為1112 631 01、0956 244 08、0705 489 55、-1538 699 52、-0198 278 99+0672 112 24i、-0.198 278 99-0.672 112 24i、-0.419 553 57+0.514 357 67i、-0.419 553 57-0.514 357 67i。

從系數和根的個數來看,系數和是負數,系數為正數,而此時具有3個正實根,1個負實根和4個虛根,該情況符合第12節中對解的性質的描述。

為了去除平凡解,使用第2個觀測時刻地球的日心距進行歸一化,得到七次方程:

求解方程,得到7個根分別為0696 823 31、1092 238 11、-1539 539 33、-0416 905 71+0515 330 21i、-0.416 905 71-0.515 330 21i、-0.201 046 69+0.672 035 89i、-0.201 046 69-0.672 035 89i。

可以得到兩個正實根,這兩個根對應的3個觀測時刻斜距如表2所示。

表2 正實根對應的斜距Table 2 Slant-ranges corresponding to the positive real roots AU

由此可知,雖然通過求解七次方程可獲得2個正實根,但是其中1個根所對應的斜距量為負值,因此只存在1個非偽解,不需要進行判別。而此時小行星位于區域,由相對位置關系可以判斷出僅存在1個非偽解,二者的判斷是一致的。可以直接計算出阿波菲斯小行星的Kepler軌道根數分別為=0923 022 94,=0190 446 69,=3321 212 1°,=204214 89°,=126611 01°,=172473 75°。其中,、、、、、分別為軌道半長軸、偏心率、軌道傾角、近地點俯角、升交點赤經和平近點角。

3.2 天基光學仿真觀測數據驗證

假設MJD=59 410166 667時刻下觀測平臺的軌道根數分別為=6 73814 km、=0001 48、=185°、=0°、=340°、=20°,空間目標的軌道根數分別為=7 17314 km、=0000 74、=943°、=340°、=630°、=10°。此時,地心、觀測平臺和空間目標歸一化后的相對幾何構型如圖5所示。得到的3個觀測時刻的角度測量值和觀測平臺的地心位置矢量如表3所示。角度測量值和地心位置矢量都位于J2000地心天球坐標系中。

圖5 地心、觀測平臺與空間目標的相對幾何構型Fig.5 Geometrical configuration of Geocentric, observation platform and space object

表3 低軌天基目標光學仿真觀測數據Table 3 Low Earthorbit object space-based optical simulation observation data

使用表3中的角度測量值與觀測平臺地心位置矢量(觀測平臺位置矢量轉化為以地球半徑=6 378.14 km為單位),按照式(6)構建Laplace形式的八次方程:

()=-13887 443+32004 384-20001 826=0

求解該方程,得到8個根為1058 812 10、1122 895 25、3637 456 07、-3804 624 64、-0599 388 75+0829 975 73i、-0.599 388 75-0.829 975 73i、-0.407 880 64+0.996 718 08i、-0.407 880 64-0.996 718 08i。

從系數和解的個數來看,系數和是負數,系數為正數,而此時具有3個正實根,1個負實根和4個虛根,該情況同樣符合第12節中對解的性質的描述。為了去除平凡解,使用測站平臺的地心距離=6 74812 km對表1中觀測平臺的位置矢量進行歸一化,構建七次方程:

求解該方程,得到7個根分別為1190 192 67、3634 720 84、-3806 610 26、-0574 040 97+0836 886 67i、-0.574 040 97-0.836 886 67i、-0.435 110 65+0.995 007 52i、-0.435 110 65-0.995 007 52i。

可以發現得到了2個正實根,這2個根都是使得斜距有意義的非偽解,需要進行判別。而此時空間目標位于區域，也可以判斷出存在2個非偽解需要判別,二者對非偽解個數的判斷是一致的。分別求解出2個非偽解對應3個觀測時刻的目標斜距,如表4所示,表中斜距單位是。

表4 非偽解對應的斜距Table 4 Slant-ranges corresponding to the non-pseudo solutions

將非偽解對應的第1和第3個觀測時刻斜距代入約束方程，可以發現非偽解1得到的方程與的值分別為0007 445 23和0005 868 73。而非偽解2得到的方程與的值分別為1634 331 25和1653 806 37。相比之下,解2對應的斜距附近不存在方程=0與=0的交點,而解1對應的斜距更接近方程的交點,因此可以判定解1為八次方程的真解。

為了更直觀地分析非偽解的選擇問題,將表3中的測量數據代入判別方程中,繪制出(,)=0和(,)=0隨和的變化曲線,如圖6和圖7所示。虛線表示(,)=0的曲線,實線表示(,)=0的曲線。圖6為斜距范圍在5倍地球半徑范圍內(,)∈[0,5]約束方程交點和Laplace形式非偽解對應點。圖7為圖6在斜距范圍在1倍地球半徑范圍內(,)∈[0,1]的放大情況。圖中點1和點2分別為非偽解1、非偽解2對應的斜距值的點。交點3,4,5為判別方程中兩個方程的交點。

圖6 5倍地球半徑范圍內約束方程的交點與Laplace形式非偽解對應的點Fig.6 Intersection points of constraint equations under five times Earth radius and the non-spurious solutions of Laplace form

圖7 1倍地球半徑范圍內約束方程的交點與Laplace形式非偽解對應的點Fig.7 Intersection points of constraint equations under Earth radius and the non-spurious solutions of Laplace form

由圖6與圖7可知,在5倍地球半徑范圍內,約束方程共有3個交點。Laplace形式的非偽解2對應的點2附近不存在交點,而非偽解對應的點1位于約束方程的交點3附近,由此可以認為非偽解1滿足本文提出的約束方程。通過進一步計算,得到非偽解1對應的目標Kepler軌道根數分別為=7 172698 km,=0000 795,=94299 305°,=62999 681°,=32456 546°,=154469 389°。

定軌結果與所選空間目標的標稱軌道根數較為符合,表明本文提出的多個非偽解判別方法的正確性。

4 結束語

八次方程的求解是近地天體和空間目標初軌確定的一個基本問題,本文從方程的性質、解和觀測構型的關系進行了解析,從數學上解釋了真解、非偽解、平凡解存在的本因,同時提出了真解的選擇方法。其中,非偽解選擇問題是利用Laplace法和Gauss法進行初軌確定需要解決的一大問題。

首先分析了初軌確定過程中,八次方程的系數和根的個數之間的關系。利用觀測平臺的地心距離進行歸一化,解決了平凡解問題。分析了非偽解的個數問題,證明了非偽解的個數只取決于地心、觀測平臺和空間目標的相對位置。隨后從求解軌道半通徑的方程出發,引入了Shefer改進方程,推導出了觀測時刻斜距應滿足的約束條件,得到了觀測時間間隔和斜距之間應滿足的非線性方程組的解析形式。利用該方程組,給出了非偽解的判別方法。最后分別利用近地天體和天基光學的短弧段觀測仿真數據進行了數值驗證,驗證了本文提出的非偽解個數的確定方法與非偽解的判別方法是正確的和有效的。

本文提出的非偽解判別方法以及建立的初軌確定方法與觀測數據組數無關,同樣適用于多資料的光學初軌確定情況。下一步工作中,將針對多組觀測資料的情況,考慮更豐富的觀測資料統計信息,形成工程實用的初軌確定方法。