近地小行星碰撞概率計算方法

李鑫冉，趙海斌,2，唐玉華，于喜雙，王秀海,5，李 彬,5

（1.中國科學院 行星科學重點實驗室，紫金山天文臺，南京 210034；2.中國科學院 比較行星學卓越創新中心，合肥 230026；3.探月與航天工程中心，北京 100195；4.國家國防科技工業局 重大專項工程中心，北京 100101；5.中國科學技術大學 天文與空間科學學院，合肥 230026）

引 言

近日點距離q≤1.3 AU 的小天體被稱為近地天體，其中大多數為日心軌道上的近地小行星。近地小行星與地球軌道發生相交或相切會給地球帶來潛在風險，其撞擊地球的過程可能會釋放數十Mt（直徑50 m）至數百萬Mt（直徑數km）的能量，對人類生命財產安全構成威脅。直徑在140 m以上與地球最小軌道距離在750萬 km以內的小行星被定義為潛在威脅小行星（Potentially Hazardous Asteroid，PHA），總數估計在幾萬顆以上，然而絕大多數還沒有被發現。這些沒有被發現的潛在威脅小行星隨時可能突然出現并與地球發生碰撞，因此需要對其監測預警并建立防御系統，以應對近地小行星撞擊的威脅。國際上早已開展關于近地小行星監測預警及防御的相關研究，并成立國際組織合作建立近地小行星的全球監測網，防范可能來自太空的威脅。

近地小行星撞擊監測作為一個重要的研究領域，碰撞概率（Impact Probability）計算是其中的關鍵因素和前提，在20多年的時間里被逐步建立起來并不斷改進和完善[1]，以提高碰撞概率的計算精度，并在此基礎上產生了多個評估系統。意大利比薩大學（University of Pisa）和空間動力學服務有限公司（SpaceDyS）研發了CLOMON和CLOMON2系統[2]。CLOMON2在近地小行星動力學網站（Near Earth Objects-Dynamic Site，NEODys）①https://newton.spacedys.com/neodys。在線發布了潛在威脅小行星在前后100年與地球發生密近交會的時間、距離及概率。美國國家航空航天局（National Aeronautics and Space Administration，NASA）的噴氣推進實驗室（Jet Propulsion Laboratory，JPL）開發了SENTRY和SENTRY-Ⅱ系統，該系統在近地天體研究中心（Center for Near Earth Object Studies，CNEOS）網站②https://cneos.jpl.nasa.gov/sentry。上給出了未來發生碰撞的概率和時間，以及后續發生密近交會的時間、距離，時間精度在1 min，距離精確到10–7AU。歐洲航天局（European Space Agency，ESA）AstOD系統的結果發表在近地天體協調中心（Near-Earth Objects Coordination Centre，NEOCC）網站①https://neo.ssa.esa.int。，給出了與地球發生密近交會的時間和距離，及發生碰撞的時間和概率，精度與SENTRY-Ⅱ相當。此外，對于臨近小行星進行碰撞監測的系統也陸續被開發出來，如JPL的SCOUT[3]、赫爾辛基大學（University of Helsinki）的NEORANGER[4]和比薩大學/SpaceDyS公司的NEOScan[5]。

近年來，中國也逐步展開了小行星監測預警、在軌處置等方面的研究。2021年航天日上提出了將論證建設近地小行星防御系統，2022年航天日上提出要完善建立地基天基對小行星的監測預警系統，組織編制近地小行星防御發展規劃，開發近地小行星防御仿真推演軟件并組織開展基本流程推演。

本文將對小行星碰撞概率計算方法研究的發展歷程進行簡要綜述，同時提出對未來該研究方向的關鍵問題和未來發展趨勢的思考。

1 碰撞概率理論基礎

1.1 碰撞概率理論基礎

直接通過小行星標稱軌道來計算其與地球的碰撞是不可靠的。由于小行星的軌道根數存在誤差，可通過小行星軌道根數的誤差協方差矩陣確定未來小行星所處的置信區域，評估其內部是否包含會與地球發生撞擊的小行星軌道簇，從而給出與地球的碰撞概率。?pik等[6-8]在20世紀50年代逐步建立了近地天體平均和長期碰撞概率的理論基礎。?pik[9]提出了第1個關于行星相遇的數學理論：假設小行星的運動可被處理為不同兩體問題解的組成，小行星相對于行星的相對速度定義了行星雙曲軌道入射漸近線的方向和速度，這個方向和速度是小行星日心軌道的半長徑a、偏心率e和傾角i的簡單函數，其中忽略了交會距離的項。密近交會時刻被計算為速度矢量在雙曲線軌道輸出漸近線方向的瞬時偏轉，忽略了太陽擾動以及小行星沿著連接2個漸近線的彎曲路徑行進實際花費的時間。根據?pik的理論，目標的軌道被認為是日心橢圓軌道，直到它進入由行星引力為主導的區域，此時軌跡變為雙曲線上的一支，當其離開行星影響范圍后進入一個新的日心橢圓軌道，其初始條件由雙曲三體問題的解給出。?pik的理論僅針對小行星和目標行星2個軌道發生實際接觸的情況，因此理論僅在最小接近距離為0時才是精確的。此外，理論沒有考慮到小行星隨后與同一顆行星（或另一顆）的相遇并不獨立于先前事件的發生。這一理論現在被稱為共振回歸，自20世紀70年代以來就被用于航天器導航，但直到20世紀90年代末才首次應用于小行星密近交會的研究[10-11]。針對以上?pik理論2個方面的不足，Valsecchi等[12]進行了發展，使得該理論可用于非0近距離交會事件的計算。

1.2 目標平面和坐標系

Kizner[13]在太陽系內星際航行中首次提出了小行星碰撞問題研究的主要工具—b平面，也稱為目標平面。每當小行星和地球之間有近距離接觸時，b平面即定義為以地球為原點且垂直于小行星關于地球相對速度矢量的平面，從而可通過軌道漸近線的位置b決定小行星與地球之間是否發生碰撞。實際應用時通常采用圖1所示的近似：垂直于近地小行星吻切軌道地心雙曲線的漸近線方向（又稱無攝相對速度方向）的平面。在密近交會中當小行星的相對速度大于地球的逃逸速度時都可以采用此近似方法。以地心為原點在b平面建立右旋坐標系(ξ,η,ζ)。圖1中，η為無攝相對速度方向v∞，即漸近線方向，垂直于b平面；ξ為地球日心速度v⊕在b平面上投影的反方向。令漸近線與b平面交點的坐標為(ξ,ζ)，其有明確的物理意義，其中 ξ為近地小行星和地球的離跡交會距離，如圖2所示，通過調整交會時間，可得到 ξ 坐標在b平面的最小值，即小行星和地球吻切軌道的最小距離，稱為地球最小的軌道距離（Minimum Orbit Intersection Distance，MOID）[14]。ζ為小行星和地球的沿跡交會距離，表示近地小行星提早或滯后于交會點的距離，間隔時間為

圖1 b平面示意圖[15]Fig.1 b-Plane[15]

圖2 b平面坐標系[14]Fig.2 b-plane coordinates[14]

其中：θ為圖1中v∞和v⊕的夾角。

其中：r⊕為地球半徑，當交點到地心的距離b

2 碰撞概率計算方法

2.1 線性碰撞概率計算方法

1993年，Muinonen等[16]發展了計算小行星碰撞概率的線性方法，將地球與小行星軌道根數誤差橢球間的距離用于評估碰撞概率，于1994年提出了潛在威脅小行星的概念，把與地球最小軌道交會距離小于0.05 AU、絕對星等小于22 mag（magnitude，星等）的小行星定義為潛在威脅小行星，確定了主要研究目標[17]。

JPL于90年代就開始研究，1994年Chodas等[18]采用線性方法預報了小行星、彗星與地球的密近交會事件，計算了軌道的不確定度和碰撞概率。根據小行星初始歷元狀態矢量的誤差協方差，通過線性過程求出目標時刻的誤差協方差從而獲得置信橢球體，置信橢球體與地球相交截面的概率密度即可用于刻畫小行星的碰撞概率。

小行星軌道根數的不確定性分析采用最小二乘法來描述，在不確定性分析中，只關注誤差協方差矩陣P。誤差協方差矩陣取決于測量誤差統計量和觀測參數的偏導數，在進行誤差協方差矩陣的狀態轉移計算P(t)=Φ(tk)P(tk)Φ(tk)T時，可采用序貫處理方法和偽序貫處理方法求取P(tk)。其中

其中：e1,e2,···,e6是歷時狀態的分量；x1,x2,···,x6是觀測時的狀態分量。Chodas[19]發現偽序貫處理方法比序貫處理方法更為穩健。在加入新的觀測結果獲得新的誤差協方差矩陣后，通過狀態轉移將誤差協方差矩陣映射到碰撞時刻。之后將由誤差協方差矩陣得到的置信橢球體進行坐標轉換，并投影到目標平面，得到置信橢圓。此時矩陣P(t)變換成C(t)，C(t)實際就是P(t)左上角的（2 × 2）子矩陣，使C(t)對角化可得

其中：σ1、σ2分別為橢圓的半長軸和半短軸。橢圓與地球在目標平面截面內二維高斯密度函數的積分就是碰撞概率為

其中：(ξI，ζI)為橢圓中心的坐標。積分可通過計算機或采用Monte Carlo統計方法近似計算。大多數情況下橢圓區域遠大于地球截面，截面面積事實上就是地球的橫截面積。因此，碰撞概率可以近似考慮為以地球為中心的正態密度函數與地球截面面積的乘積，從而簡化計算過程。

1998年，Marsden[10]采用線性方法對近地小行星1997 XF11未來與地球發生碰撞的可能性進行了評估，認為將于2028年10月26日相撞，隨后基于1990—1997年的累計觀測數據進行重新計算排除了這一可能性。但研究發現，1997 XF11與地球在2028年發生密近交會后，軌道可能會受影響從而呈現非線性特征，使得線性的碰撞概率計算方法無法有效預測未來的碰撞概率。關于1997XF11碰撞概率的研究推動了撞擊危險評估理論和算法的發展，出現了非線性碰撞概率的計算方法，同時提出了自動化碰撞監測預警系統的概念。

2.2 非線性碰撞概率計算方法

線性方法存在局限性，當初始狀態矢量的置信橢球過于細長，或外推過程中出現非線性影響時，狀態矢量映射到目標平面上的過程就會出現非線性化，導致初始歷元狀態矢量置信橢球體本身拉長或在目標平面上的投影橢圓非常狹長，從而引起線性方法的失效。另一方面，觀測弧段過短、軌道誤差較大的目標在計算時的過大誤差，也會導致無法使用線性方法。由此產生了對非線性方法計算碰撞概率的研究。

2.2.1 Monte Carlo方法

1996年，Chodas等[20]利用非線性方法研究了Shoemaker-Levy彗星軌道的演化歷史，是非線性方法的最早應用，在軌道元素空間生成1 000個隨機點填充6維不確定性的橢球體，以獲得與軌道解的協方差矩陣一致的初始條件集合，計算發現彗星可能在1929 ±9年被木星俘獲，同時確定該彗星碰撞的可能性高達95%，并給出了碎片撞擊的時間。1999年Chodas等[21]用Monte Carlo方法取樣觀測歷元軌道根數的線性6維置信區間，然后采用非線性方法數值積分到碰撞時刻以計算碰撞概率，通過將初始歷元置信橢球體劃分為若干個子區域，對子區域進行線性化處理。Monte Carlo方法的主要思想是利用頻次來估計概率，方法直接采用最小二乘法原理，對軌道元素空間中的概率分布進行采樣，方法可對置信區域（Confidence Region）充分采樣，并且計算思路簡單，能對非線性方法進行完整的體現，但樣本數量過大會帶來計算成本的負擔。

2.2.2 變化線（Line Of Variation，LOV）方法

由于軌道的不確定性過大，需要考慮置信區域中一組殘差在固定閾值范圍內的軌道，通過抽樣獲得虛擬小行星（Virtual Asteroid）的軌道[22]，判斷它們在軌道外推時是否會與地球發生撞擊。發生碰撞的虛擬小行星被稱為虛擬碰撞體（Virtual Impactor），一旦確定發生碰撞的虛擬碰撞體，即獲得了導致碰撞的初始條件，可以進一步計算碰撞概率。但虛擬碰撞體的碰撞概率與虛擬碰撞體在軌道元素空間中的體積成比例，如果碰撞概率較小，那么采樣就會變得非常密集。因此需平衡搜索的完備性和計算成本。

意大利比薩大學Milani等[11]采用多重解方法研究了小行星1999 AN10，首次證認了碰撞解和觀測數據之間的相關性。同時還提出了一維抽樣的變化線（Line Of Variation，LOV）方法：在軌道根數空間中沿著一條可微的曲線進行采樣{x(σi)},i=–M,···,M，對應于LOV的參數{σi},i=–M,···,M，在某些情況下該曲線可以代表置信區域的主軸，外推時軌道的不確定性沿著軌道延伸。沿著曲線計算對應于協方差矩陣最大特征值的單位特征向量λ1=(σ1)2，在特征空間中以h=σ/p步長前進，其中p為整數。由于是在曲線上采點，點x的下一步x+h×σ1V1可能不再屬于該曲線，但相距不遠，需要對其進行微分修正[23]。在LOV變化線上進行采樣的實質就是獲得虛擬小行星集，也就是計算多個解，多解法在觀測歷元的軌道根數非線性置信區域的變量中心軸上取樣，然后數值積分到碰撞時刻來估計碰撞概率[23-24]。在對虛擬小行星外推時，目標是找到最小的交會距離。以交會距離的平方r2為LOV參數σ的函數。r2(σ) 是可微的，可以搜索使f(σ)=為0的項。如果有閉合區間[σ1,σ2]，其中的σ對于所有初始條件都會發生密近交會，并且導數在2個極值的值為f(σ1)<0，f(σ2)>0，則至少有1個r2(σ)的最小值在區間內。Milani等[25]對目標平面進行了修正，定義垂直于最接近時小行星的地心速度矢量的平面為目標平面。由于在目標平面上的分析是局部的，因此可在每個虛擬小行星的附近進行線性近似，通過對誤差協方差矩陣的計算獲得投影橢圓的半長軸和半短軸，即LOV上軌跡的拉伸長度和寬度。令=是LOV軌道在目標平面上的軌跡。偏導數Df將置信橢球體ZX映射到目標平面上的置信橢圓ZY上。利用協方差進行傳播，定義橢圓ZY

式（7）為目標平面上（2 × 2）的矩陣，其特征值的平方根是橢圓半長軸和半長軸的長度，也就是LOV軌道處拉伸的長度和寬度。

在檢測每個虛擬小行星未來與地球的碰撞事件時，有些虛擬小行星完全不會發生碰撞事件，有些則可能發生多次。Milani等[26]首先將密近交會的軌道按時間分類，然后每一類被進一步劃分為連續的LOV段，通過分析LOV段搜索虛擬碰撞體。當給定的LOV段有多個點位于目標平面時，那么此LOV段的拉伸長度較短，可采用局部的線性方法進一步計算。相反，在強非線性情況下，LOV段拉伸很長，并且兩點之間的變化很快，在這種情況下需要在每個虛擬小行星的鄰域中進行局部分析。

由于虛擬小行星是一組對平滑曲線采樣獲得的點，Milani等[26]和Tommei[27]假設2個連續的虛擬小行星在目標平面上的軌跡連線穿越地球截面，如果軌跡的幾何形狀簡單，則可采用插值方法獲得位于地球截面內的點，即假設有2個連續的虛擬小行星xi和xi+1在目標平面上的軌跡點yi和yi+1的連線穿越地球截面，則插值方法可提供xi+δ，0 <δ<1，使yi+δ在地球撞擊截面內，而xi+δ附近有一個虛擬碰撞體，通過計算以該點為中心的概率密度函數可得到碰撞概率。

Milani等[24-26]根據分析，計算2004 FU4在2010年的碰撞概率為4 × 10–8；1997 AE12與地球的密近交會距離0.09 AU，隨著進一步觀測之后排除了對地球的威脅；同時計算了1997 XF11在2028年的碰撞概率，發現可能性極小。Tommei[27]討論了小行星Apophis在2029年和2036年的碰撞概率，其中2029年的碰撞概率高達2%。

LOV變化線方法較為復雜，但相對Monte Carlo方法減少了計算量。在采樣方式上，LOV變化線上的概率密度是由高斯方法定義，但LOV曲線中間的概率密度相比LOV兩端要高得多，采用固定步長進行采樣可能導致虛擬小行星的遺漏，因此采取與概率密度成反比的步長進行采樣是更優的選擇[28]，但該采樣方式也會導致計算成本的增加。當置信區形狀不穩定，出現多個方向時，一維采樣方式可能會出現遺漏虛擬小行星的情況。同時小行星在與地球密近交會后可能導致置信橢球體映射過程中的強非線性，從而導致置信區域拉伸較大且變化迅速或LOV采樣點處的拉伸非常大，使得LOV變化線變得復雜，帶來虛擬碰撞體搜索失敗的問題。此時需適當地增加采樣密度，使得給定的LOV段在目標平面上有多個點，從而進行連續虛擬小行星的插值計算[29]。

3 碰撞概率計算方法的改進

隨著對小行星探測技術的提升和觀測數據的增加，發現了越來越多對地球具有潛在危險的小行星。而小行星與地球發生碰撞事件的場景常常是較為復雜的，僅僅依靠經典的線性和非線性碰撞概率計算方法，無法滿足對碰撞概率高精度快速預報的要求，例如小行星飛掠地球的速度過快或者過慢，或已經發生過密近交會的小行星再次飛掠地球，都會使經典碰撞概率計算方法的結果產生偏差。因此根據實際需求，對碰撞概率計算方法進行改善的研究不斷涌現。

3.1 變化流形（Manifold Of Variations，MOV）方法

巡天觀測的發展使得小行星的觀測數據大大增加，但同時巡天觀測獲得的弧段過短，給小行星軌道確定和危險評估帶來了很大的問題[30]。多解法是在置信區域內進行采樣，當觀測到的弧段足夠長時，置信區域的形狀是拉長狀的，因此可采用LOV方式以一維曲線作為該區域的主軸進行近似處理。但當觀測到的弧段很短時，置信區域變成平坦的盤面，在單一方向上不再具有明顯優勢，而是在2個方向上都有一定的寬度，導致一維LOV變化線無法完全代表置信區域的分布，從而無法獲得有代表性的采樣點。

對于由于弧段太短而無法采用傳統方法確定完整軌道的觀測弧段，稱之為極短?。═oo Short Arc，TSA）[31]。Milani等[31]通過線性回歸和其它擬合程序來獲得極短弧的直線弧段，將極短弧由可歸屬項（Attributable）來表示，可歸屬項由參考時間（觀測時間的平均值）下的2個平均角坐標和2個相應的角速率構成。令r和q分別表示目標和觀測站為時刻t在地球上的日心位置向量。利用少量的觀測數據計算赤經α、赤緯δ，以及其導數和，通過多項式擬合可得到可歸屬項（α，δ，，）[32]?？蓺w屬項中包含的信息無法求得目標相對太陽的距離r和徑向速度，由此引入容許區域（Admissible Region）的概念，在受約束的范圍內隨機生成距離r和徑向速度，與可歸屬項共同計算小行星的軌道根數。約束條件包括日心兩體能量、地心兩體能量、地球影響范圍半徑和地球半徑?？杀硎緸?/p>

其中：A=[α,δ,,]，B=[r,]。由于與通常情況不同，（6 × 6）的協方差矩陣不可用，不確定性需要描述為

其中：μ>0是一個參數；A0為可歸屬項角坐標的標稱值；CA0是對應的正規矩陣。此集合不是笛卡爾乘積，盡管在許多情況下它可以近似為A空間中置信橢球與由A0計算得到的容許區域D(A)的乘積

由于容許區域是緊致的，可以用有限數量的點對其進行采樣[5]。如果存在標稱解，則在標稱解的附近進行蛛網采樣[27]，否則采用（r,）上的矩形網格在容許區域進行掃描[33]。

蛛網采樣利用最小化觀測殘差RMS目標函數二次近似的水平曲線來獲得。目標函數的水平曲線是6維軌道根數空間中的同心5維橢球體，因此（r,）空間上的水平曲線由5維的邊緣橢球表示，然后對該曲線進行采樣。

矩形網格的采樣分為兩步，首先根據日心能量和容許區域中連接部分的數量劃分采樣區域并進行采樣；然后根據已有的網格應用雙約束微分校正，獲得完整的軌道樣本，并找到（r,）中的最小值和最大值，以此為約束在新的網格中進一步采樣。

極短弧誤差協方差矩陣的特征值往往存在有2個特征值比其它特征值大得多的情況，表明不確定性存在于2個方向上，因此置信區域具有二維結構[1]?？梢栽谌菰S區域上定義二維流形，即變化流形MOV，進行二維采樣，并采用條件概率的方法將殘差空間的概率密度函數非線性傳播到采樣空間。采用蛛網采樣的方式，利用最小化觀測殘差函數曲線來構建覆蓋平面（r,）內置信區域的蛛網。在每條水平曲線上選擇與固定方向相對應的點作為虛擬小行星，在橢圓極坐標系（R,θ）中建立正則網格，0≤θ<2π，0 ≤R≤MRMS（MRMS是定義RMS最大值的參數）。根據初始軌道的協方差矩陣和軌道本身，對網格的每個點進行變換。

1）將（R,θ）空間的網格映射到距離速度平面

其中：λ1和λ2是變量（r,）的（2 × 2）協方差矩陣的特征值；V1是對應于較大特征值λ1的特征向量。

2）由于標稱解在距離速度平面的原點，需要利用（rnorm,norm）向量對所有點進行移動，表示用最小二乘擬合計算出的初始軌道的真實位置

隨后通過微分修正的方法進行后續校正。

（1）由4維向量A0和（r,）平面上的（r0,）構成6維向量；

（2）考慮固定方向上的下一個點（r1,），構成（A0,r1,）；

（3）對A0進行微分修正得到，從而得到（r1,）；

（4）重復前面的步驟，分析沿著h方向的所有點。

Vigna[34]對2014 AA進行了計算，根據有3個觀測數據的短弧發現其有3%的碰撞概率，結合后續7個觀測數據，其碰撞概率提升為100%。

3.2 軌道測距方法（Ranging Methods）

對于在即將撞擊前幾天或幾小時才檢測到的小行星，經典碰撞概率算法不再適用。除了可采用容許區域的方法，還可利用軌道測距方法來計算，測距方法主要分為2類，軌道統計測距（Statistical Methods）和軌道系統測距（Systematic Methods）。

Muinonen等[34]采用軌道統計測距方法嚴格映射軌道根數的概率密度，預測了小行星1998 OX4的碰撞概率。在Virtanen等[35]利用軌道統計測距技術計算軌道根數概率密度的基礎上，Muinonen等不再從整個觀測集中隨機選擇兩對觀測值，而是使用一對觀測值，對赤經和赤緯隨機采樣生成樣本，同時樣本內觀測值和計算值之間的殘差在定義閾值內，并且改進了迭代過程，減少了計算成本。Muinonen等[34]首先生成大量無偏的軌道根數，將其分布作為映射軌道根數的概率密度；其次在每個軌道根數集附近，沿著最小二乘協方差矩陣的主特征向量計算一維區間；然后在每個一維區間中計算在給定時間內導致碰撞的根數的子區間，由此可以根據導致碰撞的子區間軌道根數的概率密度計算小行星的碰撞概率。則小行星的碰撞概率為

其中：N1為初始軌道區間 ?i的數目；N2為每個軌道區間的樣本軌道數；表示發生碰撞的軌道區間i；Λ為觀測誤差的協方差矩陣；?ψ(P)為天空平面的殘差；Pij和分別為軌道區間和碰撞軌道區間中的樣本軌道j。通常碰撞概率可以簡化計算為

Muinonen等[34]據此給出了1998 OX4在2014、2038、2044、2046年的碰撞概率分別為5 × 10–7、2 ×10–5、9 × 10–6、4 × 10–5。

Chesley[38]和Farnocchia等[39]引入了軌道系統測距，Farnocchia等[39]采用系統測距的方法在弱約束的地心距離和速率的空間中進行掃描，而位置和速率與觀測值直接相關。根據觀測數據可以獲得小行星在天球上位置和速度的4個分量，即可歸屬項（α,δ,,），再利用掃描的方式在均勻分布的網格中采樣獲取地心距離和速率兩個分量（r,），即可得到小行星在極坐標的完整坐標。對每個網格固定r=ri，=，找到使函數最小的可歸屬項Aij最佳擬合

其中：v為天體測量殘差的向量；W為權重矩陣。Q的最小值通過微分校正得到

起 始 的A為(α1,δ1，(α N–α1)/(tN–t1)，(δN–δ1)/(tN–t1))，N表示最后一個觀測值。受約束的最佳擬合解可以轉換為軌道，然后傳播該軌道以尋找未來的密近交會。由于每個網格點都可找到可歸屬項的最佳擬合解，并且包含了協方差矩陣，因此可以進行碰撞概率的計算。Farnocchia等[39]給出了2008 TC3、2014AA等多顆小行星的碰撞概率，根據觀測點數量的不同，撞擊概率均約為10–3～1.0。

3.3 “鎖眼”（Keyholes）

1998年關于小行星1997 XF11是否會在未來與地球相撞的討論，推動了對撞擊危險評估的重大研究，并導致了“鎖眼”概念的確定。若地球軌道周期和密近交會后小行星的軌道周期構成簡單公約，則未來再次在相同位置發生的密近交會稱為共振回歸。如果周期的比率很接近，則可能會發生后續的交會，但時間將比上一次早或晚[40]。導致共振碰撞的b平面坐標位于可預測的Valsecchi圓上，軌道不確定性區域和Valsecchi圓之間的交點稱為“鎖眼”，由Chodas[41]最早提出這一概念，它表示特定近距離相遇的b平面上的小區域，如果小行星經過其中的一個“鎖眼”，那么它將在隨后的回歸中撞擊地球?！版i眼”也可以用于指示b平面上的某個區域，該區域不一定導致之后的碰撞，但可以導致非常近的密近交會?！版i眼”與在給定日期發生的下一次碰撞的小行星的軌道半長徑相關聯。

“鎖眼”可用于危險小行星的軌道偏轉策略優化，如果“鎖眼”的尺寸較小，則在碰撞發生前使小行星發生小幅度的偏轉并移出“鎖眼”相比其它的偏轉方案會更加容易。

3.4 Yarkovsky效應

對于大多數小行星來說，影響碰撞預測的不確定性是由計算軌道的不確定性引起的。但對于某些小行星，不確定性的主要來源是非引力攝動，尤其是Yarkovsky效應[42]，主要表現在半長徑的長期變化。當軌道受到很好的約束，并且當小行星碰撞概率需要考慮長期傳播的影響時，Yarkovsky效應也會成為撞擊危險評估問題的關鍵考慮因素。Milani等[43]發現在模擬Bennu的軌跡和評估其撞擊危險時，需要考慮的關鍵因素是Yarkovsky效應，這種由各向異性發射的熱輻射引起的反沖加速度會導致半長軸的變化[44]。Yarkovsky效應引起的半長徑變化為

其中：ρ為小行星的密度；R為小行星的半徑；γ為小行星赤道平面相對于其軌道平面的傾角。當自轉方向發生變化時，半長徑變化的符號會改變。對于1999RQ36，計算得到半長徑的變化約為–12.5 ± 5 × 10–4AU/106a。Chesley等[45]發現了2 175—2196年的幾種可能的密近交會和撞擊，累積的撞擊概率為3.7 × 10–4。

Farnocchia等[46]通過統計分析的方法評估了小行星（99942）Apophis撞擊地球的風險，將軌道解和Yarkovsky效應的不確定性考慮在內。采用Monte Carlo方法模擬了Yarkovsky效應，充分考慮了物理特性的不確定性，尤其是自轉方向的不確定性。之后將不確定性信息映射到2029年的b平面上并識別了與后續撞擊相對應的“鎖眼”，評估了未來碰撞的危險等級。之后Farnocchia等[47]又對Bennu小行星在考慮Yarkovsky效應的情況下進行了碰撞危害的評估。Vokrouhlicky[48]結合Yarkovsky效應將不確定性映射到2029年的密近交會，并計算了導致共振回歸的已知“鎖眼”的權重，修改了21世紀下半葉Apophis的撞擊概率。Spoto等[49]在加入Yarkovsky效應后將軌道空間從6維擴展到7維，并采用LOV方法和Monte Carlo方法計算了2009 FD的碰撞概率。Tardioli等[50]在Farnocchia的基礎上利用“鎖眼”和Yarkovsky效應研究了小行星碰撞概率的上限。

3.5 碰撞概率計算方法效率和采樣完備性的改進

3.5.1 Monte Carlo方法的效率提升Romano等[51]采用2種基于Monte Carlo方法的采樣方式，線采樣[52]和子集模擬[53-55]，提高了碰撞概率經典的Monte Carlo方法的性能。兩者以不同的方式對初始不確定性區域進行采樣。線采樣方法通過使用線而不是點來搜索不確定性區域，然后沿著線的方向進行積分從而計算碰撞概率，提高了碰撞概率計算的準確性。在確定采樣方向后，在采樣區域內隨機產生θk（k=1,2,···,N）以生成線進行搜索，以函數Y(c)的值來評估沿著線搜索不確定性區域時線與碰撞區域的交點c

對于每條線 θk，和撞擊區域之間最多可以找到2個交點使得=0。則隨機數 θk對應的碰撞概率為

其中：?表示單位高斯概率，總的碰撞概率為

子集模擬方法將碰撞概率計算過程變為條件概率的乘積，通過逐步識別向碰撞事件趨近的條件，減少采樣所需的樣本總數。給定碰撞事件I，令I1?I2? ···?In=I為一系列的中間事件，則Ik=，碰撞概率可表示為。以計算最小地心距離的時間區間為區分不同事件的標準，一旦單個事件及其歷元時刻被確認，可以得到計算最小地心距離的樣本的時間區間，計算最小地心距離后可得到事件I1，采用MCMC方法在新區域生成樣本，從而得到事件I2，重復此過程則得到最終的碰撞概率。Romano等[51]對2010 RF12、2017 RH16和Apophis小行星的碰撞概率進行了計算，在計算結果精度相當的情況下提高了計算效率。

3.5.2 碰撞偽觀測值法（Impact Pseudo-Observation，IOBS）

Roa等[56]提出了在軌道確定的過程中，將碰撞條件作為觀測值尋找參數空間中導致碰撞區域的方法。撞擊偽觀測值的殘差是近距離接近時b平面的坐標，不確定性設為撞擊位置的地球橫截面弦的1/2，若與觀測值兼容，則進行軌道外推最終收斂到發生碰撞的解。觀測數據通常約束軌道不確定性分布的強方向，而撞擊偽觀測值約束弱方向。該方法外推軌道時沿著弱方向迭代，以最小化MOID的值飛向地球，在不簡化動力學假設的情況下，探索了軌道不確定性的多維分布空間。

搜索虛擬碰撞體的步驟分為2步：首先確定會發生密近交會的軌道，之后依據密近交會的時間對軌道進行分類尋找會發生碰撞的軌道。IOBS方法在對密近交會的軌跡進行分類后，對在b平面上的交點加入殘差作為偽觀測值，約束不確定性的弱方向

其中：bs(t?)表示歷元t?時刻相對于地球的b平面上的坐標，∥bs(t?)∥

其中：α為b平面上不確定性橢球的縮放因子，通常α=1。在加入偽觀測值后，外推過程經過不確定性區域向碰撞事件演變，而不是在不確定區域中尋找會發生碰撞的軌道。

最后根據軌道不確定性的初始分布利用Monte Carlo方法采樣可計算碰撞概率。當虛擬碰撞體在參數空間中的分布被概率密度函數完整描述時，則可以采用重點抽樣等方差縮減方法對虛擬碰撞體進行采樣后計算碰撞概率。引入碰撞概率函數pi(x)，則碰撞概率可以表示為

其中：Fi為參數空間中在某一日期撞擊地球的軌道所屬的區域；q(x)為觀測值置信橢球的概率密度分布。計算的難點在于找到適當的pi(x)對虛擬碰撞體的分布進行建模，而引入偽觀測值則可以解決這一問題，因為碰撞解的協方差對虛擬碰撞體分布進行了描述。

Roa等[56]采用此方法計算了2008 TS10、2017 LD等小行星的碰撞概率，與SENTRY系統相比，在10–7碰撞概率的情況下，二者計算結果相當，當碰撞概率減小到4 × 10–8后，該方法更加穩健。

4 碰撞概率計算方法總結和展望

4.1 碰撞概率計算方法框架

綜合以上對于碰撞概率計算方法的梳理，目前主要的碰撞概率計算流程的框架圖如圖3所示，在置信橢球體較為簡單時可直接采用線性方法快速計算碰撞概率。當線性方法失效時，可根據獲得的軌道根數是否為可靠的最小二乘解，選擇在置信區域或容許區域內進行采樣，依據計算成本和精度要求選擇不同采樣方式獲得虛擬小行星，對虛擬小行星進行外推獲得可與地球發生碰撞的虛擬碰撞體，最終得到小行星的碰撞概率。同時，在計算過程中可以加入IOBS等方法對計算的效率和穩定性進行改善。通常碰撞概率的計算軟件采用其中的一種方法進行計算，也會同時列出線性方法和非線性方法進行自主選擇。

圖3 碰撞概率計算流程圖Fig.3 Flow chart for impact probability

第1個被開發出來的碰撞監測系統CLOMON，及后續改進的Sentry 和 CLOMON2 系統均采用了置信區域下的LOV方法。針對臨近小行星進行碰撞監測的系統SCOUT，NEORANGER和NEOScan系統則采用了容許區域下的MOV方法。2021年7月Sentry-Ⅱ采用了新的IOBS方法進行碰撞概率的計算，在計算速度上Sentry-Ⅱ相比Sentry系統要慢，但搜索虛擬碰撞體的穩定性和準確性得到了大幅提升。

4.2 關鍵問題和未來展望

目前小行星碰撞概率計算方法在國外已較為成熟，后續改進主要集中于提高計算效率和提升搜索的完備性。

大量巡天觀測計劃的實施使得近地小行星觀測數據海量劇增，其中短弧數據越來越多，而基于短弧定軌對碰撞概率的快速評估仍然是危險評估的難點，因此有效利用短弧數據對近地小行星進行危險評估也成為未來研究的重點。另一方面，對于信息量缺乏的短弧數據，也可通過軌道關聯方法將多段短弧連接為足以進行精確計算的弧長，進而計算準確的碰撞概率。我國可以通過建立可靠的短弧段關聯的計算方法極大地促進小行星碰撞概率計算準確性的提高。

實際問題中小行星的軌道較為復雜，采用單一評估方法計算負荷高且可能導致撞擊概率的計算遺漏，后續可對不同觀測情況的小行星分類分階段處理，對擁有長期觀測弧段的小行星進行高精度的碰撞概率計算，對即將發生碰撞的小行星則進行快速預測，從而建立可應對多種場景的中國獨立自主的小行星監測預警系統。

小行星碰撞概率計算的準確性依賴于軌道精度，高精度的觀測資料可以很好地約束置信區域，減少計算成本，簡化計算方法。因此，采用多源數據觀測可提高觀測精度，從根本上提升小行星碰撞概率計算的準確性。中國目前正在建立的“中國復眼”，以及規劃中的大口徑高精度天體測量望遠鏡，對小行星監測預警系統的建立有著重要意義，準確的多源觀測資料將會有力支持小行星防御系統的建設。

5 結束語

小行星防御需要對小行星進行危險評估和潛在威脅小行星的編目管理，而碰撞概率的計算是評估近地小行星威脅程度最關鍵的因素。本文梳理了現有計算小行星碰撞概率的基本方法和框架流程，并對中國未來的發展方向進行了展望。鑒于中國目前在碰撞預警方面的能力與國外有較大差距，中國應找準突破口，從自身特色出發建立適用于中國的監測預警系統。