李鑫冉 趙海斌

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

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

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

引言

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

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

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

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

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

1 進化算法

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

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

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

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

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

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

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

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

(1)DE/rand/1

(2)DE/best/1

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

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

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

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

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

3 (a,e,M)優(yōu)選法

3.1 變量選擇

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

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

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

3.3 適值函數(shù)

令已知一組觀測量{ti,αi,δi,i=1,2,...,n},(αi,δi)代表ti時刻的赤經(jīng)和赤緯,則由歷元時刻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得到.由于地球與小行星的為相對位置分地內(nèi)和地外兩種情況,因此依據(jù)R和L的夾角對觀測幾何進行分類討論.

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

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

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

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

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

(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),此時有適值函數(shù)

可以看出,適值越小表示個體越優(yōu).

圖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,?,ω)取多組結(jié)果的中值作為最后結(jié)果.

4 數(shù)值試驗

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

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

表1 小行星實測數(shù)據(jù)與模擬數(shù)據(jù)定軌結(jié)果Table 1 The results of orbit determination with measured data and simulated data

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

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

表2 小行星軌道根數(shù)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 組數(shù)據(jù)進行進行分類計算:當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].為避免隨機數(shù)對結(jié)果的影響,采用不同隨機數(shù)對每組數(shù)據(jù)重復計算300 次.計算發(fā)現(xiàn),與在空間碎片的計算結(jié)果不同,優(yōu)化結(jié)果的適值差異明顯,并不像近圓軌道的解那樣彼此接近,因此適值的差異性同樣需要關(guān)注.表3 列出了2001 FR85,2006 SV189,2019 UJ10 各自適值最小的前5 組優(yōu)選結(jié)果.算法雖然不同于空間碎片近圓軌道下可以迅速準確找到軌道信息的表現(xiàn),但從適值大小分析,真實解在適值上仍具有較為明顯的優(yōu)勢,且小偏心率的適值優(yōu)勢比大偏心率更加突出.僅依靠1-3 天的觀測數(shù)據(jù),得到的最優(yōu)解與MPC 中給的軌道根數(shù)基本一致.

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

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

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

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

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

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

表4 小行星2011 FS2 軌道根數(shù)Table 4 Orbital Elements of 2011 FS2

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

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

表5 加入誤差的定軌結(jié)果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 結(jié)論

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

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

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