申敬松

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

劉 也* 胡松杰

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

韓 潮

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

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

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

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

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

1.1 碎片初選

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

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

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

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

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

1.2 建模流程

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

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

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

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

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

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

2.1 碎片分布

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

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

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

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

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

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

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

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

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

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

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

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

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

進一步，為適應(yīng)精簡模型的應(yīng)用需求，在參數(shù)合理設(shè)定的基礎(chǔ)上，對碎片駐留時間比和通量進行分段的三維多項式擬合，并僅存儲多項式系數(shù).這樣就大大減少了模型的存儲量和預(yù)處理時間，而應(yīng)用時可利用這些參數(shù)快速恢復(fù)駐留時間比與通量等數(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)，長期的軌道進動會使碎片在經(jīng)度上趨于均勻分布.此外，航天發(fā)射場的緯度分布導(dǎo)致了碎片軌道傾角的分布具有峰值特性，這會表現(xiàn)在碎片空域的緯度分布上［12］，故對于刻畫長期演化規(guī)律，可計算高度-緯度帶內(nèi)的平均空間密度:

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

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

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

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

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

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

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

圖2給出了2011年和2014年碎片空間密度的計算結(jié)果，可見ORDEM2000與RTS的曲線變化趨勢基本一致，即碎片空間密度均隨高度降低而減小，這種現(xiàn)象的一個重要原因是低軌大氣阻力加大導(dǎo)致碎片快速下降所致.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 駐留時間比的經(jīng)緯度分布情況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則給出了它們在緯度帶內(nèi)的分布情況.可見緯度分布上存在分布峰值，且兩極地區(qū)碎片密度較少;420 km的碎片密度明顯多于300 km;300 km的碎片分布在個別軌道上比較密集，而420 km的碎片在經(jīng)度分布上則已比較均勻.進一步，通過綜合分析其他年度不同高度層碎片分布情況，表明:低軌碎片緯度分布主要集中在-60°～60°之間，基本上呈現(xiàn)南北對稱特性，但有一定的分布差異，其中20°，40°和60°緯度帶在一些高度層，特別是400～700 km高度，出現(xiàn)了局部分布密度較大的帶狀峰值;350 km高度以下碎片的經(jīng)度分布有一定的不均勻性，但隨著高度增加碎片在經(jīng)度上分布趨于均勻.

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

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

4 結(jié)論

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

References)

［1］王海福，馮順山，劉有英.空間碎片導(dǎo)論［M］.北京:科學(xué)出版社，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è)大學(xué)，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è)大學(xué)，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］.中國空間科學(xué)技術(shù)，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］王海福，余慶波，劉有英.空間碎片碰撞風(fēng)險評估系統(tǒng)［J］.北京理工大學(xué)學(xué)報，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ù)學(xué)手冊》編寫組.數(shù)學(xué)手冊［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］.長沙:國防科技大學(xué)出版社，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