申敬松

(北京航空航天大學 宇航學院，北京100191)

劉 也* 胡松杰

(北京航天飛行控制中心 航天飛行動力學技術重點實驗室，北京100094)

韓 潮

(北京航空航天大學 宇航學院，北京100191)

人類開發(fā)太空的腳步給空間帶來了大量的碎片累積，這些碎片主要分布在近地軌道和同步軌道，具有高速運動和分布隨機化的特征，給在軌正常運行的航天器，特別是長期運行的載人航天器帶來了極大的安全隱患［1-2］.因此，對于需長期運行于低地軌道的載人航天器來說，為評估碎片碰撞帶來的安全風險并采取有效的預防措施，需要有效描述碎片分布和演化規(guī)律的環(huán)境建模技術的支持［2-3］.

空間碎片環(huán)境建模的方法主要是利用已有觀測數(shù)據(jù)，綜合考慮航天器區(qū)域分布和碎片的尺寸、速度、碰撞角和來源等因素，并在一定簡化假設下，借助數(shù)據(jù)擬合和概率統(tǒng)計等建立的［4］.已公開提供的模型按計算方式可以分為演化模式和工程模式，前者立足于采用分析手段刻畫空間碎片分布的演化規(guī)律，但其建模精度依賴于碎片分布影響因素的有效認識，這在目前理論知識下還是很難實現(xiàn)的.工程模式則依據(jù)碎片測量數(shù)據(jù)建模，可以給出比較準確的碎片分布和碰撞風險情況，適合分析空間環(huán)境的綜合影響，已廣泛服務于航天器設計、操控和為碎片觀測計劃制定等工程任務［1，3］，如美國國家宇航局(NASA)、歐洲空間局(ESA)、俄羅斯國家宇航局(RASA)的ORDEM，MASTER和SDPA等系列模型［5-6］.其中 NASA的 ORDEM 模型已更新到ORDEM2010，并在官方網(wǎng)站公開提供了ORDEM2000的免費版本［7-8］.

國內方面，一些單位也開展了空間碎片工程模式的分析和建模工作［3-5，9-10］，但由于可用數(shù)據(jù)限制和研究起步較晚的緣故，許多科學研究和工程實踐還依賴于國外公開的環(huán)境模型.為此，本文針對低軌載人航天器長期在軌運行需求，試圖通過碎片軌道測量數(shù)據(jù)統(tǒng)計分析，建立一種實用性較強的低軌空間碎片環(huán)境建模方法，并通過建模實例進行碎片時空分布規(guī)律的分析，以期為航天器碰撞風險評估與防護設計等工作提供借鑒.

1 碎片環(huán)境建模流程

1.1 碎片初選

碎片環(huán)境模型是描述近地軌道空間碎片三維分布特性的，且需要考慮其隨時間的變化情況.為便于統(tǒng)計分析，需要進行時空的離散化，并進行區(qū)域內碎片(也包含正常飛行目標)的遴選.

碎片的數(shù)量和分布是受航天發(fā)射、在軌解體與碰撞、大氣阻力和人為緩解(如回收)等因素綜合影響的，其中航天發(fā)射和大氣阻力等都是隨季節(jié)變化的，因此，鑒于統(tǒng)計時間與模型粒度和準確性的關系，這里以季度為基本時間單元統(tǒng)計碎片分布狀況.

空間的劃分相對簡單，對于近地空間可以按照距地面高度h、緯度φ和經度λ做出三維的網(wǎng)格劃分.為了處理球殼兩極附近的網(wǎng)格，引入最大緯度φ0，將每個高度層內大于φ0的網(wǎng)格合并為同一網(wǎng)格.網(wǎng)格劃分要兼顧計算量、存儲量和精度等要素，也可以引入相鄰高度層的網(wǎng)格比例，以調節(jié)不同高度層網(wǎng)格的大小.

碎片初選是為了剔除明顯不會飛入所考察空域的碎片.對于高度區(qū)域［hmin，hmax］，當碎片軌道近地點高度大于hmax或遠地點高度小于hmin，碎片才可能進入該高度層.例如TLE(Two Line Elements)編目數(shù)據(jù)可提供目標的偏心率e和角速度n，記地球物理常數(shù) μ，平均半徑r0，則可得穿越［hmin，hmax］的碎片需滿足的角速度取值范圍［nmin，nmax］:

其中Δn是為了補償球形地球和瞬時圓軌道假設的誤差，保證初選碎片不出現(xiàn)漏警.

1.2 建模流程

在模型覆蓋時空多維離散化的基礎上，利用經過初選的編目碎片數(shù)據(jù)，設計了如圖1所示的建模流程.基本過程分為3個步驟:

1)根據(jù)數(shù)據(jù)更新頻率和精度進行軌道預報，并計算與該段軌道交會的網(wǎng)格及交會起止時間;

2)在基本時間單元內(季度)，統(tǒng)計碎片在各網(wǎng)格內的駐留(飛行)時間和通量，獲得工程模式所需的基礎數(shù)據(jù)源;

3)根據(jù)應用需求，內嵌一系列數(shù)據(jù)插值和時間序列分析等方法，輸出應用模型所需的結果.

圖1 空間碎片環(huán)境建模流程Fig.1 Modeling diagram for space debris environment

2 碎片環(huán)境建模方法

2.1 碎片分布

下面采用搜索方式，依次在高度維、緯度維和經度維尋找與碎片交會的網(wǎng)格及起止時刻.

對碎片j第s段預報軌道，根據(jù)軌道時空連續(xù)性特點，可將其劃分為若干擬合區(qū)間［ti，s，ti+1，s］.對第 i 個 區(qū) 間 的 高 度 數(shù) 據(jù)，記 ai，sh3+bi，sh2+ci，sh+di，s=0 為其三次多項式擬合結果.為便于時間統(tǒng)計，需要將擬合區(qū)間劃分為幾個單調區(qū)間.導數(shù)方程求根得極值時刻ti1，s=［-2b+(4b2-12ac)1/2］/6a 與 ti2，s=［-2b-(4b2-12ac)1/2］/6a.進而可分為如下4種情況:

1)ti1，s與 ti2，s均在區(qū)間內，此時軌道按高度可劃分 單 調 區(qū) 間 ［ti，s，ti1，s］，［ti1，s，ti2，s］與 ［ti2，s，ti+1，s］;

2)ti1，s不在區(qū)間內，此時軌道按高度可劃分單調區(qū)間［ti，s，ti2，s］與［ti2，s，ti+1，s］;

3)ti2，s不在區(qū)間內，此時軌道按高度可劃分單調區(qū)間［ti，s，ti1，s］與［ti1，s，ti+1，s］;

4)ti1，s與 ti2，s均不在區(qū)間內，此時軌道按高度為一個單調區(qū)間［ti，j，ti+1，j］.

對每個單調區(qū)間［tb，te］對應的［hb，he］，根據(jù)設定的高度劃分Δh來搜索與其交會的高度層.若高度層［hm，hm+Δh］?［hb，he］，則根據(jù)三次方程求根公式［11］和根的范圍可計算交會時間［th，m，th，m+1］.同樣的，對［th，m，th，m+1］內的緯度數(shù)據(jù)采用類似的擬合、劃分單調區(qū)間和求根步驟，可獲得碎片交會的緯度維編號 n及起止時間［tφ，n，tφ，n+1］.進一步，同樣可計算［tφ，n，tφ，n+1］內碎片交會的經度維編號 k 及起止時間［tλ，k，tλ，k+1］.

綜上，就搜索到了與碎片交會的網(wǎng)格在高度、緯度與經度編號(m，n，k)，可計算駐留時間為

記該時間單元(季度)總時間T，對s，j和i求和，可得碎片駐留時間比:

航天器在空域中的碰撞風險分析常用到碎片通量的概念，它可以定義為該處碎片數(shù)量密度與速度的乘積［5］.通量方向由航天器和碎片相對速度方向決定，通常慣性空間一點的通量分為單方向通量和各方向全通量.當不考慮特定航天器時，可以采用全通量方式(即各方向通量之和)刻畫碎片特征，此時若假設每個網(wǎng)格內碎片速度是常值，記網(wǎng)格體積為V(m，n，k)，則可計算碎片的全通量:

式(4)若不考慮網(wǎng)格體積，則為網(wǎng)格內的碎片通量.當分析具體航天器的軌道時，在精度要求不高的情況下，在近地軌道區(qū)域可以認為單個位置的空間碎片通量全部集中在該位置的當?shù)厮矫鎯?，?ORDEM2000 工程模式［3］.

2.2 模型數(shù)據(jù)

模型的基礎數(shù)據(jù)源就是統(tǒng)計的碎片駐留時間和通量，可以將每個時間單元的統(tǒng)計結果存儲為一個基本文件.為了節(jié)省存儲空間，文件內除包含必要的時空網(wǎng)格參數(shù)，只需按順序存儲量化后的各網(wǎng)格的駐留時間比和碎片通量.順序存儲的數(shù)據(jù)也便于開發(fā)輔助工具進行可視化等處理.

進一步，為適應精簡模型的應用需求，在參數(shù)合理設定的基礎上，對碎片駐留時間比和通量進行分段的三維多項式擬合，并僅存儲多項式系數(shù).這樣就大大減少了模型的存儲量和預處理時間，而應用時可利用這些參數(shù)快速恢復駐留時間比與通量等數(shù)據(jù)，或直接進行其他參數(shù)的計算.

2.3 模型分析

駐留時間比 r(m，n，k)刻畫了網(wǎng)格(m，n，k)出現(xiàn)碎片的概率，記網(wǎng)格體積為V(m，n，k)，可得空間密度:

多數(shù)的空間碎片處于無主動機動狀態(tài)，長期的軌道進動會使碎片在經度上趨于均勻分布.此外，航天發(fā)射場的緯度分布導致了碎片軌道傾角的分布具有峰值特性，這會表現(xiàn)在碎片空域的緯度分布上［12］，故對于刻畫長期演化規(guī)律，可計算高度-緯度帶內的平均空間密度:

低軌載人航天器通常飛行在高度相對固定的近圓軌道，空間碎片環(huán)境分析時也通常給出各個高度層內的碎片分布，如ORDEM模型［7-8］.記高度層［hb，he］的體積為 V(b，e)，可以統(tǒng)計該高度層內的平均空間密度:類似地，為航天器長期運行的碰撞風險綜合評估提供參考，可以統(tǒng)計各個層次的碎片平均速度和通量等.特別地，給定航天器軌道oL和周期TL(L代表近地軌道衛(wèi)星)，航天器軌道周期內的平均碎片通量計算如下:

根據(jù)所存儲的分季度網(wǎng)格數(shù)據(jù)，在空間上采用三次多項式插值，可以給出該季度模型覆蓋空域內任意空間位置的碎片駐留時間，也就可以給出任意位置的碎片密度、平均速度和通量等.這與ORDEM模型的地基觀測視角［7］類似，只是時空劃分更加細致，但由于數(shù)據(jù)源限制，未給出分布與碎片尺寸的關系.

模型將碎片空間分布按季度存儲為單獨的文件，隨著數(shù)據(jù)積累就形成了一個時間序列.因此，通過分析碎片歷史數(shù)據(jù)的時間變化規(guī)律，可以分析碎片的時間演化特征.除阻力消減外，由于碎片發(fā)射頻率、在軌分裂與碰撞、人工緩解等因素的精確建模都很困難且嚴重依賴于數(shù)據(jù)積累，因此這里不考慮碎片演化的物理因素，直接采用時間序列方法分析歷史數(shù)據(jù)的變化規(guī)律，并據(jù)此進行碎片演化預報.以空間密度為例，若在無碎片驟變的平穩(wěn)段落，采用AR模型擬合與預報，計算公式為

其中，ε是噪聲序列;模型系數(shù)φ在擬合階段通過Yule-Walker方程求解.若存在碎片解體等突發(fā)事件時，考慮非平穩(wěn)的時間序列分析方法［13］.

目前主要給出兩類預報結果:一是高度層內空間碎片的分布，二是航天器軌道內平均流量，它們均以季度為時間粒度.當然，根據(jù)研究需要，也可利用駐留時間比文件進行其他參數(shù)的統(tǒng)計、擬合和預報分析，這將在以后研究中予以補充.

3 基于TLE的建模實例

前述的統(tǒng)計建模過程需要輸入碎片的軌道數(shù)據(jù)，并不嚴格限制數(shù)據(jù)來源，所用處理方法有一定通用性.本節(jié)進行模型實例化與分析，數(shù)據(jù)來源是北美防空司令部利用美國空間監(jiān)視網(wǎng)發(fā)布的TLE數(shù)據(jù)庫，這是目前可公開獲得的最完整且更新性最好的空間碎片數(shù)據(jù)源［1，7］.對其他編目數(shù)據(jù)可類似處理.

目前，TLE數(shù)據(jù)可提供10cm以上和部分1 cm以上的空間碎片編目數(shù)據(jù)，基本保證1～2 d更新一次，且提供了比較精確的軌道預報模型［14-15］.為了節(jié)省計算時間，采用120 s的采樣間隔，對每個碎片根據(jù)其軌道高度情況分別用 SDP4和SGP4模型［14］預報各編目數(shù)據(jù)更新周期(通常為1 d)內的軌道，并轉化為軌道高度、緯度和經度，進而采用前述的方法和流程進行碎片環(huán)境統(tǒng)計建模.以下簡稱獲得的模型為駐留時間統(tǒng)計模型(RTS，Residential Time Statistical model).

ORDEM2000模型可以按年度和尺度給出碎片的空間密度分布，其中大于10 cm的碎片與TLE有類似的測量數(shù)據(jù)來源，因此，將RTS與其進行對比驗證.

圖2給出了2011年和2014年碎片空間密度的計算結果，可見ORDEM2000與RTS的曲線變化趨勢基本一致，即碎片空間密度均隨高度降低而減小，這種現(xiàn)象的一個重要原因是低軌大氣阻力加大導致碎片快速下降所致.RTS曲線相對平滑，這是由于其采用三次多項式插值而ORDEM2000采用線性插值所致.此外，RTS的密度幅度略高.這是因為此處的ORDEM2000僅統(tǒng)計了10 cm以上的碎片且數(shù)據(jù)時間較早，而RTS所用新的TLE數(shù)據(jù)包含了10 cm以下碎片.

圖2 本文模型(RTS)與ORDEM2000模型的空間密度對比Fig.2 Space debris density according to RTS and ORDEM2000

圖3 駐留時間比的經緯度分布情況Fig.3 Geodetic distribution of the residual time ration

圖4 駐留時間比的緯度分布情況Fig.4 Latitude distribution of the residual time ration

圖3是2010年第一季度300km和420km高度層的碎片駐留時間比(以dB形式給出以便于顯示)，圖4則給出了它們在緯度帶內的分布情況.可見緯度分布上存在分布峰值，且兩極地區(qū)碎片密度較少;420 km的碎片密度明顯多于300 km;300 km的碎片分布在個別軌道上比較密集，而420 km的碎片在經度分布上則已比較均勻.進一步，通過綜合分析其他年度不同高度層碎片分布情況，表明:低軌碎片緯度分布主要集中在-60°～60°之間，基本上呈現(xiàn)南北對稱特性，但有一定的分布差異，其中20°，40°和60°緯度帶在一些高度層，特別是400～700 km高度，出現(xiàn)了局部分布密度較大的帶狀峰值;350 km高度以下碎片的經度分布有一定的不均勻性，但隨著高度增加碎片在經度上分布趨于均勻.

圖5是碎片空間密度隨時間的變化曲線，可見各高度層內的碎片空間密度相對比較穩(wěn)定，這些結論與ORDEM和MASTER等一些公開模型的分析結果是基本一致的［4-7］.

圖5 碎片空間密度隨時間的變化情況Fig.5 Space debris density changes with time

4 結論

低地軌道空間碎片的增多已經嚴重影響了航天器的正常工作和安全飛行，因此，建立空間環(huán)境模型以評估航天器在軌運行的碰撞風險是十分必要和迫切的.本文借鑒國內外空間碎片環(huán)境工程模式的建模思想，綜合利用軌道預報和多項式擬合及求根方法，進行了空間碎片分布統(tǒng)計，并在此基礎上給出了基于多項式插值和時間序列分析的碎片環(huán)境模型描述方法，該建模方法可以適用于各種編目軌道數(shù)據(jù).以TLE數(shù)據(jù)為例進行了模型驗證和空間碎片分布規(guī)律的研究，獲得了與ORDEM2000模型類似但尺度更加細致的計算結果.文中所述碎片環(huán)境建模方法和分析結論可為低軌道航天器，特別是近地載人航天器的碰撞風險評估和防護設計等工作提供參考.后續(xù)工作將集中于模型長時間預報性能的改進，以及模型在航天器軌道設計和碰撞規(guī)避能力設計等工作中的具體應用.

References)

［1］王海福，馮順山，劉有英.空間碎片導論［M］.北京:科學出版社，2010 Wang Haifu，F(xiàn)eng Shunshan，Liu Youying.Introduction of space debris［M］.Beijing:Science Press，2010(in Chinese)

［2］ Johnson N L.History of on-orbit satellite fragmentations［R］.NASA/TM-2008-214779，2008

［3］彭科科.空間碎片環(huán)境探測數(shù)據(jù)處理方法及工程模型建模方法研究［D］.哈爾濱:哈爾濱工業(yè)大學，2010 Peng Keke.Research on sapce debris environment exploration data processing measure and engineering modeling method［D］.Harbin:Harbin Institute of Technology，2010(in Chinese)

［4］張平平.近地軌道空間碎片軌道參數(shù)分布規(guī)律研究［D］.哈爾濱:哈爾濱工業(yè)大學，2006 Zhang Pingping.Research on space debris orbital parameter distribution uule in LEO［D］.Harbin:Harbin Institute of Technology，2006(in Chinese)

［5］唐頎，龐寶君，張偉.空間碎片環(huán)境工程模式參數(shù)分析［J］.中國空間科學技術，2004，24(5):22-27 Tang Qi，Pang Baojun，Zhang Wei.Analysis of the parameters in space debris environment engineering model［J］.Chinese Space Science and Technology，2004，24(5):22-27(in Chinese)

［6］ Fukushige S，Akahoshi Y，Kitazawa Y.Comparison of debris environment models:ORDEM2000，MASTER2001 and MASTER2005［J］.IHI Engineering Review，2007，40(1):31-41

［7］ Liou J C.The new NASA orbital debris engineering model［R］.ORDEM2000，NASA/TP-2002-210780，2002

［8］ Liou J C.NASA's ORDEM2010 status［C］//28th Inter-agency Debris Coordination Committee Meeting(IADC).Trivandrum，India:NASA，2010:1-20

［9］王海福，余慶波，劉有英.空間碎片碰撞風險評估系統(tǒng)［J］.北京理工大學學報，2008，28(12):1039-1042 Wang Haifu，Yu Qingbo，Liu Youying.Orbital debris risk assessment system［J］.Transactions of Beijing Institute of Technology，2008，28(12):1039-1042(in Chinese)

［10］ Xu X L，Xiong Y Q.A research on the collision probability calculation of space debris for nonlinear relative motions［J］.Chinese Astronomy and Astrophysics，2011(35):304-317

［11］《數(shù)學手冊》編寫組.數(shù)學手冊［M］.北京:高等教育出版社，2010:87-89“Handbook of Mathematics”Writing Group.Handbook of mathematics［M］.Beijing:Higher Education Press，2010:87-89(in Chinese)

［12］ Klinkrad H.Space debris-models and risk analysis［M］.New York:Springer＆ Praxis，2006

［13］王正明，易東云.測量數(shù)據(jù)建模與參數(shù)估計［M］.長沙:國防科技大學出版社，1996 Wang Zhengming，Yi Dongyun.Instrumental data modeling and parameter estimation［M］.Changsha:University of National Defence Technology Press，1996(in Chinese)

［14］ Felix R，Ronald L.Space track report No.3:models for propagation of NORAD element sets［R］.Peterson:Aerospace Defense Command，United States Air Force，1980

［15］ Levit C，Marshall W.Improved orbit predictions using two-line elements［J］.Advances in Space Research，2011，47(7):1107-1115