近地小行星極短弧定軌的進化算法研究1)

李鑫冉 趙海斌

?(中國科學院行星科學重點實驗室,紫金山天文臺南京 210034)

?(月球與行星科學國家重點實驗室,澳門科技大學澳門 999078)

??(中國科學院比較行星學卓越創新中心合肥 230026)

引言

近地小行星是太陽系內一類特殊的天體,部分近地小行星軌道可能與地球相交,對地球安全和人類生存環境構成潛在威脅,如6500 萬年的全球物種大滅絕[1]、2013 年俄羅斯的車里雅賓斯克隕石墜落事件[2-3]、發現于2004 年6 月著名的危險小行星(99942)Apophis[4].因此,發現、監測近地小行星并計算其與地球的碰撞概率、開展危險程度的評估等相關研究是十分重要的,而其中能利用較少數據盡早確定近地小行星的軌道參數尤為關鍵,小行星參數的精確度對碰撞模型和危害評估的結果有很大影響[5],準確的軌道參數可以為后續的預警工作提供可靠的輸入,而這就涉及到小行星的定軌問題.

NASA 在 2005 年提出對至少 90%的直徑超過 140 m 的近地天體進行編目和特性獲取[6],Pan-STARRS[7]、Catalina[8]、NEOWISE[9]及LSST[10]、NEOCam[11]等大量的近地小行星大視場巡天項目的開展使得巡天能力不斷增強,得到了大量的觀測數據.但同時,新的觀測方式也使得無法對巡天中探測到的每一個目標進行后續的跟蹤觀測,因此獲得的弧長都很短,通常只有一個晚上的拍攝[12].為了提高巡天效率,未來采集的數據將會更為稀疏,并且軌道參數分布范圍很廣,包含眾多大偏心率軌道,這些短而稀疏的數據給軌道確定以及識別帶來了很大困難.對于這些過短的觀測弧段尤其是大偏心率極短弧段,利用傳統的Laplace 和Gauss 方法無法進行定軌,加之短弧定軌本身具有的病態性[12-16],使得定軌難度大幅度增加.由此,如何有效利用這些數據對小行星進行極短弧定軌,對巡天項目的充分利用及小行星的探測研究都有著重要意義.近年來針對這類問題,極短弧定軌的概念被明確提出并成為研究熱點.極短弧的具體弧長目前尚無嚴格的定義,通常無法用經典方法得到合理定軌結果的觀測弧段即稱為極短弧,以區別于傳統意義上的短弧定軌[12-13,17-18].

除經典計算方法外,優選法也可被利用來解決定軌問題,優選法克服了經典方法中迭代不收斂的現象,但更適合于解決一維的優選問題.對于多維情況,計算過程過于復雜.此外,方法對于初值的要求較高,而初值選取本身就是一個初軌計算問題,對于極短弧軌道計算問題也不適用.

Ansalone 和Curti[18]針對極短弧定軌問題下天基的模擬資料,應用遺傳算法(genetic algorithm,GA),將觀測首末時刻的斜距作為優選變量,使定軌問題轉換為一個優化問題,但采用的參數與通常選法相差較多.王志勝等[19]傳算法運用到短弧定軌的問題上來,研究基于測角資料的衛星短弧定軌.劉磊等[20]將遺傳算法應用于天基的短弧定軌問題,采用雙ρ 迭代模型,對稀疏數據進行定軌.李鑫冉和王歆[21-22]對參數調整后將遺傳算法應用到了極短弧定軌問題中,在地基的空間目標定軌問題中得到了較好的應用.進化算法在優選問題中可以將生物進化中優勝劣汰的現象應用到對最優解的搜尋中,在探索過程中通過積累經驗,啟發式地尋找最終解.算法已有較為成熟的理論基礎,在多個領域都有研究和應用[23-24].算法對先驗信息的依賴性較小,受數據中的噪聲影響也較小,采用進化算法研究定軌問題已成為新的趨勢.除遺傳算法外,進化算法中還包含多種不同進化機制的算法,算法各有特點和優勢,粒子群算法(particle swarm optimization,PSO)、差分進化算法(differential evolution,DE),及基于統計學思想的分布估計法(estimation of distribution algorithm,EDA)等已被應用于解決空間目標的短弧定軌問題,并在近圓軌道下有較好表現.

本文將進化算法引入小行星的極短弧定軌問題,構建計算框架,以差分進化算法為代表采用模擬資料進行計算驗證,并比較算法在不同偏心率下短弧定軌問題中的表現,探討大偏心率下算法的特征.

1 進化算法

20 世紀60 年代進化算法基于模擬自然進化的方法首次被提出,70 年代出現了相關的理論研究,直至21 世紀基本成熟[25].

這其中最具代表性的就是GA 算法,算法通過模擬自然進化中優勝劣汰的過程搜索最優解,基于適應度來選擇父代進行雜交,GA 算法產生的子代有概率發生變異,從而在進化的同時尋找新的可能性.從20世紀70 年代被Holland 和De Jong 提出以來被廣泛應用,收斂性和全局搜索能力已得到證明[25].

PSO 算法于1995 年被Kennedy 和Eberhart[26]提出,源于對鳥類捕食行為的研究,即鳥類找到食物最簡單有效的方法就是搜尋當前距離食物最近的鳥的附近區域.與GA 算法不同,它基于群體智能而不是遺傳操作,利用群體中個體對信息的分享,使整個群體的運動在問題求解空間中產生從無序到有序的演化,最終獲得最優解.

DE 方法在1996 年由Storn 和Price[27]提出,是目前最有效的隨機參數優選算法之一,它模擬生物進化,將初始種群中兩個個體的向量差作為變異方向,疊加到第三個個體上,以此產生新個體,反復迭代使得適應環境的個體被保留下來.不同于遺傳算法原本采用二進制編碼適用于離散問題的求解,它適用于求解連續變量的優化問題.算法構造思想借鑒了GA 算法和PSO 算法,算法保留了GA 算法的進化過程,同時以類似于PSO 算法中的更新方法替代GA 算法中的遺傳操作,因此參數和算子較少,計算復雜性降低.

進化算法的進化機制多種多樣,但通常的流程是相同的.一般將需要被優選的變量稱為個體,通過優選方法隨機生成一定數量的個體組成初始種群.種群內個體數目稱為種群數,種群數越大搜索能力也越強,但計算效率隨之降低.適值函數用于評估個體優劣,即優選法中的目標函數,對于最小化問題,個體適值越小則越優秀.通過進化在滿足預定條件時終止即得到最優解.具體流程如圖1.

圖1 進化算法流程圖Fig.1 The flowchart of Evolutionary Algorithm

2 初軌計算的差分進化算法

DE 算法的主要步驟與GA 算法類似,主要包括變異(mutation)、交叉(crossover)、選擇(selection)三種操作,但次序不同.算法隨機生成初始種群X={x1,x2,...,xNP},其中NP為種群數,xi={xi1,xi2,...,xiD}為D維向量,D為優選變量的維數.DE 算法先進行變異操作,對每個個體xi變異得到個體Vi,變異方式較GA 算法大為簡化,常見的變異方式有以下三種

(1)DE/rand/1

(2)DE/best/1

(3)DE/target-to-best/1

其中r為互不相同的均勻分布的隨機整數,ri∈[1,NP]且ri≠i;S為縮放因子,一般在區間[0,1]中取值,但多數文獻建議取較大的值,綜合文獻[27-30],范圍在[0.4,1]比較合適,S取值決定了算法的全局搜索能力,越大的取值全局搜索能力越強;xbest為由適值函數所確定的當代最優的個體,DE/*/*是DE 算法變異的表達方式,兩個* 依次表示變異基和差分數量.經變異所得的種群V={V1,V2,...,VNP}和原種群X交叉操作,得到新種群U.具體如下

其中r為[0,1]區間均勻分布的隨機數,rand為[1,D]上均勻分布的隨機整數,CR為交叉概率.此操作使得Ui以一定概率接受變異個體的分量,但確保至少有一個分量來自變異個體,CR決定了種群的多樣性,文獻中建議CR取0.1 或0.9[27-30]作為初始嘗試值.最后進行選擇操作,DE 算法采用了貪婪操作,如果新個體優于初始個體,則取而代之,否則初始個體保留下來,進入下一次的進化

這里k表示進化代數,表示xi(k)進化到第k代的個體,函數F(·)表示求解個體的適值.從選擇方式可看出最優個體一定會進入下一代,每一代種群不會劣于前一代.

通過上述三個操作完成了一次種群的進化,并通過不斷迭代求解出最優解.算法常用操作中的選擇和交叉操作都只有一種方式,而變異操作選擇也比較少,且基本形式是相同的.

3 (a,e,M)優選法

3.1 變量選擇

對于優化問題,優化變量過高會帶來求解困難的問題,即使如今計算能力已有了大幅度的提升.因此本文采用了3 個Kepler 根數,即歷元時刻t0的(a,e,M0)作為優化變量,與Ansalone 和Curti[18]采用首尾觀測時刻的斜距作為優化變量的方法不同,在只增加一維的情況下,使得優化結果不再需要依賴觀測量就可以得到完整的解,便于資料處理.

3.2 初始種群生成生成和終止條件

根據先驗信息定義優選變量的值域,由于進化算法對初值要求較低,無確切信息時可將范圍取的大一些:a∈[al,au],e∈[el,eu],M0∈[Ml,Mu].初始種群中每個個體的每個變量都在取值范圍內隨機選取,重復NP次即得到整個種群{xNP}.終止條件選取較為普通的迭代次數達最大進化代數G終止或連續C代沒有進化.

3.3 適值函數

令已知一組觀測量{ti,αi,δi,i=1,2,...,n},(αi,δi)代表ti時刻的赤經和赤緯,則由歷元時刻t0的(a,e,M0)可得ti時刻黃道坐標下的近點角Mi,fi和Ei,進一步可得

其中Li=(cos δicos αi,cos δisin αi,sin δi)T,ri為ti時刻目標的日心位置矢量,ri=|ri|,ρi為目標斜距,Ri為測站的日心位置矢量可由測站地心位置矢量Re和地日位置矢量RS得到.由于地球與小行星的為相對位置分地內和地外兩種情況,因此依據R和L的夾角對觀測幾何進行分類討論.

當R·L<0 時,即圖2 所示,此時觀測目標的所在位置有A,B,C 三種情況:

(2)當|r| >|R|時,即小行星的軌道高于地球,處于位置C;

(3)其他,此時目標位置有兩種可能A 或B,無法根據已有條件確定其具體位置,需分別計算.但這一計算非常容易,不會造成過多負擔.

當R·L>0 時,即圖3 所示,有2 種情況.

(1)|r| <|R| 時,目標軌道低于地球,但考慮到R和L的夾角,此情況不可能發生;

(2)其他,此時目標只可能位于A,即軌道在地球之上.

圖2 R·L<0 時地球與小行星的位置Fig.2 The locations of Earth and asteroid when R·L<0

基于上述分析,則可得式(1).當目標位于圖2 位置B、C,及圖3 位置A 時,式中取“+”,其他情況取“-”,則任意一對觀測時刻(tk,tj)且tk>tj,可得對應的(rk,rj)、(fk,fj),此時有適值函數

可以看出,適值越小表示個體越優.

圖3 R·L>0 時地球與小行星的位置Fig.3 Locations of Earth and asteroid when R·L>0

3.4 (i,?,ω)的求解

通過以上計算已可得t0的(a,e,M0),從而得到每個觀測時刻的位置矢量ri,從而可得(i,?,ω),考慮到計算精度,可由每對(ri,rj)得到的(i,?,ω)取多組結果的中值作為最后結果.

4 數值試驗

基于DE 算法采用MATLAB 編寫程序,算法參數選擇NP=300,S=1.0,CR=0.9,變異方法選擇DE/rand/1,最大迭代次數G=200,連續迭代次數C=30 最優適值的相對變化小于10?12則提前結束計算.考慮到近地小行星,令值域選擇范圍為a∈[0.8,4.0],M∈[0,2π].

選取三組偏心率不同的軌道,分別計算其軌道根數,并與傳統的Laplace 方法進行比較.表1 給出了實測數據與模擬數據的定軌結果,模擬數據基于其觀測數據的初始時刻、觀測時刻和測站數據生成,同時保留了觀測的幾何構型,其中POD 代表已獲得的軌道根數,作為參考標準.

表1 小行星實測數據與模擬數據定軌結果Table 1 The results of orbit determination with measured data and simulated data

可以看出,模擬數據下兩種方法都可得到初軌結果,Laplace 方法更接近準確值.采用實測數據時,當偏心率較小,DE 算法的結果偏差稍大,當偏心率逐漸增大時,Laplace 方法的結果偏離程度增大,至e>0.6時,已得不到有效結果,而DE 算法雖然出現偏差,依然可以得到有效解為后續工作提供軌道范圍的參考信息.另一方面,Laplace 方法只能由單一解判斷軌道信息,當計算出現困難時,得到的結果完全無效,無法指導后續工作.而進化算法的結果并不僅僅是單一的解,有效范圍內的解都是有效的,可根據多組解的分布判斷結果的有效性,并提示其存在范圍.對于大偏心率極短弧軌道,DE 算法的適用性更廣.

極短弧定軌問題本身存在困難,當偏心率增大時變得更為復雜,稀疏數據中的誤差也可能對計算帶來很大影響,因此,為了重點關注算法的計算規律,采用模擬數據對問題進行簡化,主要探討DE 算法在極短弧下解的特征.增加不同偏心率的小行星進行比較,選取MPC 中小行星共9 組,已經確定的軌道如表2,偏心率覆蓋[0,0.7]的范圍,每組數據的時間跨度為1-3 天,數據點不足10 個.

表2 小行星軌道根數Table 2 Orbital elements of asteroids

圖4 給出了小行星2001 FR85 在一次完整計算過程中軌道半長徑a和適值的變化.其中e∈[0.0,0.3],NP=300.可以看出DE 算法的效率很高,收斂速度很快.

圖4 半長徑a 和適值F 的收斂過程Fig.4 The convergence process of the semi-major and fitness value

考慮到大偏心率定軌相較近圓軌道更為復雜,不易求解,試驗時偏心率e的取值范圍不直接擴大到[0,1],而是對9 組數據進行進行分類計算:當e∈[0.0,0.3]時,值域選擇范圍為e∈[0.0,0.3];當e∈[0.3,0.6]時,值域選擇范圍為e∈[0.3,0.6];當e∈[0.6,1)時,值域選擇范圍為e∈[0.6,0.9].為避免隨機數對結果的影響,采用不同隨機數對每組數據重復計算300 次.計算發現,與在空間碎片的計算結果不同,優化結果的適值差異明顯,并不像近圓軌道的解那樣彼此接近,因此適值的差異性同樣需要關注.表3 列出了2001 FR85,2006 SV189,2019 UJ10 各自適值最小的前5 組優選結果.算法雖然不同于空間碎片近圓軌道下可以迅速準確找到軌道信息的表現,但從適值大小分析,真實解在適值上仍具有較為明顯的優勢,且小偏心率的適值優勢比大偏心率更加突出.僅依靠1-3 天的觀測數據,得到的最優解與MPC 中給的軌道根數基本一致.

圖5 給出了9 條軌道的(a,e)概率密度分布圖,圖中顏色越淺表示聚集度越高.可以看到,當偏心率較小(e<0.1)時,最優解主要集中分布在真實軌道附近,且與適值最小的結果相吻合,DE 算法可以得到有效的結果.而當偏心率逐漸增大(e>0.3)時,求得解的分布區域發生偏離,或出現多個分布區域,且分布最集中的區域也并不是真實解的所在區域,分布不再明顯集中于真實軌道附近.2018 XB5 和2011 HT 有多個分布聚集的區域,真實解包含在其中的某個區域但不是聚集度最有優勢的區域.雖然DE 算法搜索到了真解,但它在獲得的多個解中并沒有明顯優勢,個體分布較少.小偏心率的軌道更為穩定從而更易被搜索得到解,而當偏心率增大時,進化過程對觀測數據的敏感性可能產生變化.

表3 各組小行星軌道根數Table 3 Orbital Elements of Different Asteroid

圖5 概率密度分布圖Fig.5 Probability Density

圖5 概率密度分布圖(續)Fig.5 Probability Density(continued)

因此,對于小偏心率軌道,算法易于直接尋找到最優解,而對于大偏心率軌道,需結合分布密度和解的適值進行選擇,偏心率的增大使軌道更為復雜,也導致算法在搜尋最優解的靈敏度上有所下降,在搜尋過程中容易發生最優解搜尋方向的偏離,雖然可以找到適值最優的解,但數量較少,適值不具有優勢的解會被大量搜索到,真實解的小范圍內的聚集現象,在整體分布上沒有明顯優勢.

2011 FS2 僅包含半天的觀測數據,而2019 UJ10偏心率較大,其在概率密度分布圖上的聚集更加難以顯現.因此將搜索的解空間進行縮小,提高算法的靈敏度,再次進行試驗.表4 中列出了2011 FS2 分別在e∈[0.0,0.3]和e∈[0.1,0.2]內搜尋的結果,均為適值最小的5 個解.可以看到隨著搜索區間的縮小,搜索效率得到提高,最優解的存在被凸顯出來.將2019 UJ10 的搜索空間縮小至e∈[0.6,0.7]得到了圖6 的概率密度分布圖.可知,雖然仍有另一個干擾解的存在,但算法明顯搜索到真實解的存在區域,并且在此范圍內呈現聚集狀.因此算法在計算大偏心率及過短弧段的軌道時,搜索空間中真實解是存在的,只是在大范圍搜索中不易搜到,聚集分布不明顯.計算時可以通過將區間進行約束劃分,分段計算最優解,來提高算法的靈敏度及搜索能力,提升聚集程度,同時,結果需結合解的分布聚集區域和適值最優的個體考慮.

表4 小行星2011 FS2 軌道根數Table 4 Orbital Elements of 2011 FS2

圖6 e ∈[0.6,0.7]概率密度分布圖Fig.6 Probability density of e ∈[0.6,0.7]

不同類型軌道受誤差的影響不同,對多組模擬數據加入隨機觀測誤差在DE 算法下進行比較.表5 給出了分別加入0.1′′,0.2′′誤差的2001 FR85,2006SV189 和2019UJ10 的定軌結果.可以看到在0.1′′的誤差下,仍可得到有效的定軌結果.當誤差擴大到0.2′′時,解的分布仍涵蓋真實解,可提示軌道信息的參考范圍,但隨著偏心率增大,受誤差影響也增大了,大偏心率的軌道定軌結果會受到明顯影響.圖7 給出了約束解空間后加入0.2′′誤差的2011 FS2 的概率密度分布圖,星號表示真實解所在位置,在加入誤差后解的分布仍覆蓋真實解.

表5 加入誤差的定軌結果Table 5 The results of orbit determination with error

圖7 加入0.2′′ 誤差2011 FS2 的概率密度分布圖Fig.7 Probability Density of 2011 FS2 with 0.2′′ error

5 結論

極短弧定軌問題采用經典方法存在很大困難,甚至無法得到有效解,包括DE 算法在內的進化算法將這一反問題正向處理,避免了經典方法固有的病態性,并且DE 算法參數較少,操作簡便,易于實現.進化算法計算框架基本一致,只是進化機制和側重不同,對于不同需求的問題改用不同進化算法時操作更為便捷,如EDA 方法更注重整體搜索,DE 算法更注重局部搜索,換用不同進化算法時計算框架可保持不變.

根據進化算法的特點,將其應用于近地小行星的極短弧定軌進行探索.對于小偏心率軌道,DE 算法利用少于3 天的少量數據得到的軌道信息與利用多天多站定軌下的信息一致,可為后續工作提供可參考的信息.對于復雜的大偏心率軌道和弧段更短的軌道,進化算法的表現不如小偏心率軌道下良好,搜索靈敏度降低,不易搜索到最優解,僅在局部有向最優解聚集的現象.因此,需要縮小搜索空間提高算法靈敏度,并結合分布區域和適值最優進行討論.加入誤差后,較小的誤差對定軌結果影響較小,隨著誤差增大,盡管適值最優解受到干擾,尤其大偏心率軌道的定軌受影響較大,解的分布仍涵蓋真實解所在區域.

小行星軌道較為多樣,大偏心率軌道數量較多,并且實際觀測中觀測位置不同、與地球的相對位置不同,也會對算法產生影響,尤其對于極短弧定軌問題,數據量較少,觀測數據差異和誤差的影響更不可忽視.在模擬數據的研究基礎上,未來需對觀測位置和時刻做進一步研究,在不同情況下分類計算,提高算法在大偏心率下的搜索效率,完善進化算法在小行星極短弧定軌方面的應用.