本發(fā)明屬于航天動力學技術和計算機技術領域，具體涉及一種大規(guī)?？臻g碎片分布演化的數(shù)值計算方法及系統(tǒng)。

背景技術：

隨著人類空間活動的日益增加，人為發(fā)射進入空間的物體不斷增多，使得地球外部空間的物體數(shù)量迅速增加，在軌航天器受到的碰撞威脅也愈加嚴重。

據(jù)nasa空間監(jiān)視網(wǎng)(ssn)的統(tǒng)計，地球同步軌道高度以下的大于10cm尺寸的空間物體數(shù)量達到15000多個，小于10cm尺寸的空間物體數(shù)量更多，粗略估計至少為400000個，致使地球外部空間達到了前所未有的擁擠狀態(tài)。這些空間物體中，只有很少的一部分是在軌正常工作的航天器，大部分是由于發(fā)射任務、在軌航天器爆炸分解、在軌航天器碰撞分解產(chǎn)生的空間碎片。因此，這些在軌道動力學的約束下沿軌道繞地球自由運動的空間碎片，會擠占有限的空間資源，使得諸如近地軌道和靜止軌道等區(qū)域可用的資源減少。同時，在大氣阻力、地球非球形攝動以及太陽、月球三體引力等攝動力的作用下，自由無控的空間碎片軌道會隨時間不斷變化，當空間碎片軌道與正常工作航天器的軌道存在交點時，就會存在與航天器發(fā)生碰撞的可能，對正常工作航天器產(chǎn)生威脅。另外，航天器發(fā)生碰撞后，會進一步使空間碎片的數(shù)量增加，從而使得產(chǎn)生新的碰撞的可能性增加，進而產(chǎn)生更多碎片，對正常運行航天器產(chǎn)生更大威脅，這種鏈式的雪崩式效應，最終會導致航天器在軌運行的成本迅速增加，甚至會出現(xiàn)近地或同步軌道范圍內(nèi)航天器均無法生存的情況。

由此可見，大量空間碎片對航天器造成的碰撞威脅，將是人類空間活動面臨的主要安全問題，如不加以管理，將導致航天器運行成本大大增加，甚至會使空間資源無法利用。有效的對當前空間碎片的運動分布狀態(tài)進行預測分析，是航天器在軌運行碰撞規(guī)避、空間碎片減緩措施制定以及確?？臻g資源可持續(xù)利用的重要的技術保障?，F(xiàn)有技術中，尚未出現(xiàn)有效的對當前空間碎片的運動分布狀態(tài)進行預測分析的相關內(nèi)容。

技術實現(xiàn)要素：

針對現(xiàn)有技術存在的缺陷，本發(fā)明提供一種大規(guī)?？臻g碎片分布演化的數(shù)值計算方法及系統(tǒng)，可有效解決上述問題。

本發(fā)明采用的技術方案如下：

本發(fā)明提供一種大規(guī)?？臻g碎片分布演化的數(shù)值計算方法，包括以下步驟：

s1，空間碎片狀態(tài)推演子系統(tǒng)獲取空間碎片的初始運動狀態(tài)數(shù)據(jù)以及空間碎片受到的空間攝動力，并將所述初始運動狀態(tài)數(shù)據(jù)以及所述空間攝動力輸入到數(shù)值積分模型中；所述數(shù)值積分模型按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)；然后，采用所述空間碎片未來運動狀態(tài)數(shù)據(jù)更新空間碎片當前運動狀態(tài)數(shù)據(jù)；

s2，判斷是否達到演化時間，如果未達到，將狀態(tài)更新后的空間碎片運動狀態(tài)和屬性數(shù)據(jù)傳輸給空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)；

s3，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)將所述狀態(tài)更新后的空間碎片運動狀態(tài)和屬性數(shù)據(jù)輸入到碎片解體模型中；所述碎片解體模型進行解體計算，得到空間碎片解體后產(chǎn)生新空間解體碎片的屬性數(shù)據(jù)；并將所述新空間解體碎片的屬性數(shù)據(jù)反饋傳輸?shù)剿隹臻g碎片狀態(tài)推演子系統(tǒng)；

s4，所述空間碎片狀態(tài)推演子系統(tǒng)將所述新空間解體碎片的屬性數(shù)據(jù)以及已有的空間碎片的運動狀態(tài)數(shù)據(jù)輸入到數(shù)值積分模型中，再次按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)，如此不斷循環(huán)s2-s4；直到達到演化時間，結束流程。

優(yōu)選的，s1中，所述空間攝動力包括大氣阻力攝動力、地球非球形引力攝動力、太陽光壓攝動力以及第三體引力攝動力；所述數(shù)值積分模型為4階阿達姆斯預估矯正積分模型。

優(yōu)選的，s3中，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)具體用于：

首先，根據(jù)能量定律確定空間碎片解體后產(chǎn)生新空間解體碎片的數(shù)量n和新空間解體碎片特征尺寸l的分布；

然后，以新空間解體碎片特征尺寸l作為獨立變量，利用概率分布模型確定新空間解體碎片的面質比a/md、新空間解體碎片的迎風截面積a、新空間解體碎片相對于解體前空間碎片的速度增量δv；

最后，根據(jù)a/md、a確定碎片質量md。

優(yōu)選的，s3中，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)具體用于：

s3.1，確定解體碎片數(shù)量：

其中，解體碎片區(qū)分為爆炸解體碎片和碰撞解體碎片；

1)對于爆炸解體碎片，根據(jù)能量守恒定律，解體碎片特征尺寸l大于被研究解體碎片最小特征尺寸lc的碎片數(shù)量nf滿足函數(shù)關系：

其中：為解體碎片的最小特征尺寸；

cs是修正系數(shù)，與爆炸類型有關，對于歷史爆炸事件，修正系數(shù)cs的表達式：

其中：

其中：h為爆炸解體碎片的高度；

2)對于碰撞解體碎片，碰撞解體碎片數(shù)量滿足函數(shù)關系：

式中：為兩碰撞目標的等效質量：

其中：

其中：ep為等效碰撞能量，為災難性碰撞臨界能量；mp為相互碰撞中較小質量；mt為相互碰撞中較大質量；vr是碰撞時兩碎片的相對速度，即碰撞速度；

s3.2，確定解體碎片尺寸分布律：

爆炸解體碎片尺寸分布律f：

碰撞解體碎片尺寸分布律f：

s3.3，確定解體碎片面質比參數(shù)：

本步驟中，爆炸解體碎片和碰撞解體碎片的碎片面質比參數(shù)滿足同樣的分布函數(shù)；

(1)對于尺寸大于11cm解體碎片，定義解體碎片面質比的對數(shù)值γ＝log(a/md)，解體碎片特征長度的對數(shù)值θ＝log(l)，解體碎片的面質比分布由如下二元正態(tài)分布決定：

p(γ,θ)＝ε(θ)p1(γ)+(1-ε(θ))p2(γ)(19)

其中：p1(γ)和p2(γ)均是正態(tài)分布概率密度函數(shù)：

其中：權重參數(shù)ε(θ)∈[0,1]，是θ的飽和函數(shù)；正態(tài)分布概率密度函數(shù)的均值μk和正態(tài)分布概率密度函數(shù)的方差σk均是θ的函數(shù)；

對于航天器解體碎片，其面質比分布函數(shù)中的參數(shù)由下式?jīng)Q定：

對于火箭上面級解體碎片，其面質比分布函數(shù)的參數(shù)由下式?jīng)Q定

(2)對于尺寸小于8cm的航天器解體碎片和尺寸小于1.7cm的火箭上面級解體碎片，碎片的面質比滿足正態(tài)分布律p1(γ)，即權重函數(shù)ε＝1，

p(γ,θ)＝p1(γ)(23)

此時，分布函數(shù)參數(shù)由下式統(tǒng)一決定：

(3)對于尺寸介于8cm～11cm之間的航天器解體碎片和尺寸介于1.7cm～11cm之間的火箭上面級碎片，碎片的面質比通過隨機采樣的方法確定，即：首先生成服從均勻分布且取值區(qū)間在[0.0,1.0]內(nèi)的隨機變量ζ，然后將ζ與利用(25)計算得到的進行對比，若則利用公式(19)計算面質比參數(shù)，反之則利用公式(23)計算面質比參數(shù)；

步驟3.4，確定解體碎片的速度增量：

解體碎片相對于解體前目標的速度增量滿足正態(tài)分布規(guī)律，定義v＝lg(δv)，γ＝lg(a/md)，則v的分布函數(shù)p(ν)滿足

其中：v為解體碎片的速度；δv為解體碎片的速度增量；

解體碎片的均值μν和解體碎片的方差σν通過下式計算：

σv＝0.4。

優(yōu)選的，還包括：

s5，空間碎片運動狀態(tài)三維顯示子系統(tǒng)將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)進行三維顯示，即：動態(tài)的顯示所有空間碎片在空間中的位置分布狀態(tài)；具體包括：

步驟s5.1，接收空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)；

步驟s5.2，調(diào)用三維顯示接口函數(shù)，將空間碎片的位置和形狀在三維模型中進行渲染；

步驟s5.3，驅動三維模型，對空間碎片運動狀態(tài)進行動態(tài)更新。

優(yōu)選的，還包括：

s6，在固定時刻或者設定的時刻，空間碎片運動狀態(tài)輸出子系統(tǒng)將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)輸出；其中，所述空間碎片未來運動狀態(tài)包括：空間碎片的位置、速度矢量和面值比、質量、來源屬性參數(shù)；此外，空間碎片運動狀態(tài)輸出子系統(tǒng)還輸出不同高度層內(nèi)空間碎片狀態(tài)數(shù)據(jù)。

優(yōu)選的，還包括：

s7，空間碎片演化結果處理子系統(tǒng)根據(jù)空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)，計算空間碎片狀態(tài)量，進一步計算到空間碎片密度；根據(jù)空間碎片密度得到空間碎片沿軌道高度方向和緯度方向上的分布狀態(tài)。

優(yōu)選的，所述空間碎片密度通過以下步驟計算得到：

步驟s7.1，確定空間碎片在高度區(qū)間[r,r′]內(nèi)的平均空間密度分布其中，r和r′分別為高度區(qū)間的最小地心距和最大地心距；包括以下步驟：

步驟s7.1.1，當目標軌道為橢圓軌道，近地點地心距和遠地點地心距分別為ra和rp，當r′≤ra且r≥rp時，通過下式計算：

其中：a為目標軌道半長軸；

步驟s7.1.2，當目標軌道為橢圓軌道，且r′＞ra或r＜rp時，則令r′＝ra或r＝rp，帶入(28)進行計算，得到

步驟s7.1.3，當目標軌道為圓軌或近圓軌道，即rp→ra∈[r,r′]時，則令r′＝ra或r＝rp，公式(28)變?yōu)楣?29)，通過公式(29)計算

步驟s7.1.4，其他情況，

步驟s7.2，確定空間碎片在緯度區(qū)間[β,β′]內(nèi)的平均空間密度分布函數(shù)其中，β和β′分別是緯度區(qū)間的最小值和最大值；

具體的，若β′＞0，則直接執(zhí)行步驟s7.2.1；

若β′≤0，則令β′＝-β，β＝-β′，均轉換為關于原點對稱的正值區(qū)間后，再執(zhí)行步驟s7.2.1；

步驟s7.2.1，碎片軌道傾角記為i，當β′≤i且β≥-i時

步驟s7.2.2，當β≤-i且β′≤i時

步驟s7.2.3，當β′≥i且-i≤β≤i時

步驟s7.2.4，當β′≥i且β≤-i時

步驟s7.2.5，當β＞i時，

本發(fā)明還提供一種大規(guī)?？臻g碎片分布演化的數(shù)值計算系統(tǒng)，包括：

空間碎片狀態(tài)推演子系統(tǒng)，用于獲取空間碎片的初始運動狀態(tài)數(shù)據(jù)以及空間碎片受到的空間攝動力，并將所述初始運動狀態(tài)數(shù)據(jù)以及所述空間攝動力輸入到數(shù)值積分模型中；所述數(shù)值積分模型按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)；然后，采用所述空間碎片未來運動狀態(tài)數(shù)據(jù)更新空間碎片當前運動狀態(tài)數(shù)據(jù)，得到空間碎片新運動狀態(tài)數(shù)據(jù)；

空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)，用于將所述空間碎片新運動狀態(tài)數(shù)據(jù)輸入到碎片解體模型中；所述碎片解體模型進行解體計算，得到空間碎片解體后產(chǎn)生新空間解體碎片的屬性數(shù)據(jù)；并將所述新空間解體碎片的屬性數(shù)據(jù)反饋傳輸?shù)剿隹臻g碎片狀態(tài)推演子系統(tǒng)；

空間碎片運動狀態(tài)三維顯示子系統(tǒng)，用于將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)進行三維顯示，即：動態(tài)的顯示所有空間碎片在空間中的位置分布狀態(tài)；

空間碎片運動狀態(tài)輸出子系統(tǒng)，用于將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)輸出；其中，所述空間碎片未來運動狀態(tài)包括：空間碎片的位置、速度矢量和面值比；此外，空間碎片運動狀態(tài)輸出子系統(tǒng)還輸出不同高度層內(nèi)空間碎片狀態(tài)數(shù)據(jù)；

空間碎片演化結果處理子系統(tǒng)，用于根據(jù)空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)，計算空間碎片狀態(tài)量，進一步計算到空間碎片密度；根據(jù)空間碎片密度得到空間碎片沿軌道高度方向和緯度方向上的分布狀態(tài)。

本發(fā)明提供的大規(guī)?？臻g碎片分布演化的數(shù)值計算方法及系統(tǒng)具有以下優(yōu)點：

本發(fā)明提供的大規(guī)?？臻g碎片分布演化的數(shù)值計算方法及系統(tǒng)，能夠對大規(guī)?？臻g碎片分布狀態(tài)進行長期演化計算，能夠動態(tài)三維演示碎片分布狀態(tài)，提供了碎片分布狀態(tài)數(shù)據(jù)的輸出和處理分析功能，可以為航天器在軌運行空間碎片碰撞規(guī)避、空間碎片減緩措施制定以及確保空間資源可持續(xù)利用的重要的保障。

附圖說明

圖1為本發(fā)明提供的大規(guī)?？臻g碎片分布演化的數(shù)值計算方法的流程示意圖；

圖2為本發(fā)明提供的空間碎片狀態(tài)推演子系統(tǒng)的結構示意圖；

圖3是大規(guī)模空間碎片分布演化數(shù)值計算系統(tǒng)運行時，空間碎片初始分布狀態(tài)的三維顯示圖；

其中：在圖3中，積分步數(shù)為0步；當前空間碎片數(shù)量為17000個；已發(fā)生碰撞次數(shù)為0次；

圖4是大規(guī)?？臻g碎片分布演化數(shù)值計算系統(tǒng)運行過程中，產(chǎn)生新的解體碎片的分布狀態(tài)示意圖。

具體實施方式

為了使本發(fā)明所解決的技術問題、技術方案及有益效果更加清楚明白，以下結合附圖及實施例，對本發(fā)明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

本發(fā)明提供一種大規(guī)?？臻g碎片分布演化的數(shù)值計算方法，應用于由5個子系統(tǒng)構成的大規(guī)?？臻g碎片分布演化的數(shù)值計算系統(tǒng)中，大規(guī)?？臻g碎片分布演化的數(shù)值計算系統(tǒng)包括：空間碎片狀態(tài)推演子系統(tǒng)、空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)、空間碎片運動狀態(tài)三維顯示子系統(tǒng)、空間碎片運動狀態(tài)輸出子系統(tǒng)和空間碎片演化結果處理子系統(tǒng)。

大規(guī)?？臻g碎片分布演化的數(shù)值計算方法包括以下步驟：

s1，空間碎片狀態(tài)推演子系統(tǒng)獲取空間碎片的初始運動狀態(tài)數(shù)據(jù)以及空間碎片受到的空間攝動力，并將所述初始運動狀態(tài)數(shù)據(jù)以及所述空間攝動力輸入到數(shù)值積分模型中；所述數(shù)值積分模型按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)；然后，采用所述空間碎片未來運動狀態(tài)數(shù)據(jù)更新空間碎片當前運動狀態(tài)數(shù)據(jù)。

本步驟中，所述空間攝動力包括大氣阻力攝動力、地球非球形引力攝動力、太陽光壓攝動力以及第三體引力攝動力；所述數(shù)值積分模型為4階阿達姆斯(adams-bashforth)預估矯正積分模型。

空間碎片狀態(tài)推演子系統(tǒng)的結構如圖2所示，具體實現(xiàn)過程如下：

步驟s1.1，計算空間碎片受到的大氣阻力攝動力。

大氣阻力作為唯一清除在軌空間碎片的因素，其作用隨軌道高度變化較為明顯。計算表明，250km的典型空間碎片，在大氣阻力的作用下，2個月內(nèi)會再入到大氣層內(nèi)；而600km高度的空間碎片，再入到大氣層的時間達到15年，高于850km的空間碎片，再入到大氣層的時間將需要數(shù)個世紀之久，而地球同步軌道處的空間碎片在通常情況下永遠不會再入到大氣層中，除非空間環(huán)境發(fā)生了大規(guī)模變化。

在軌空間碎片受到的大氣阻力隨著大氣狀態(tài)的不同而變化，阻力加速度計算式為

其中：是大氣阻力攝動力；cd是空間碎片的阻力系數(shù)；a是空間碎片的迎風截面積；md是空間碎片的質量；ρ是空間碎片所在位置的大氣密度；是空間碎片與當?shù)卮髿獾南鄬λ俣龋?/p>

假設大氣相對地球靜止，即大氣相對于地球的速度其中為地球自轉角速度，為地心矢徑；若空間碎片的速度為則相對速度阻力加速度的計算式可進一步改寫為

步驟s1.2，計算空間碎片受到的地球非球形引力攝動力。

在地固坐標系中，地球引力位函數(shù)作為拉普拉斯方程的解，其非球形部分u為：

其中：和為勒讓德多項式相關項，定義為：

其中：

g為引力常數(shù)，me是地球質量，r為地心距；

u表示地球引力場位函數(shù)；

λ和分別表示單位質點在地固坐標系中的地心經(jīng)和緯度；

ae表示地球平均半徑；

為規(guī)格化的締合勒讓德多項式；

和為規(guī)格化的地球引力位系數(shù)；

n和m為多項式的階和次，n為取的最高階數(shù)。

步驟s2.3，計算空間碎片受到的太陽光壓攝動力。

太陽光壓攝動力加速度可以表示為

其中：為太陽光壓攝動力加速度；psp＝4.5605×10-6n/m2為地球附近太陽光壓常數(shù)；cr為太陽光壓系數(shù)；asr是與太陽垂直的目標橫截面積；au＝1.49597870×1011m是iau1976天文常數(shù)系統(tǒng)給出的天文單位；為太陽位置矢量。

步驟s2.4，計算空間碎片受到的第三體引力攝動力。

第三體引力攝動力加速度為

其中：為第三體引力攝動力加速度；s和l分別代表太陽和月球；mj為太陽或月球質量；代表太陽或月球的地心矢徑。太陽和月球的位置可以通過平根公式計算，也可以采用jpl星歷確定。

步驟s2.5，利用4階阿達姆斯預估矯正方法對空間碎片的運動狀態(tài)進行積分更新。該積分方法利用4階龍格-庫塔方法進行參數(shù)初始化，充分利用歷史計算結果，每步積分僅計算一次攝動力，減少積分計算量；采用4階積分方法，比高階積分具有跟高的計算穩(wěn)定性。

4階龍格-庫塔參數(shù)初始化方法為：

針對初值問題：

其中：t是時刻；y是狀態(tài)量，如空間碎片的位置和速度矢量的坐標；是狀態(tài)變量y的導數(shù)；f(t,y)是約束函數(shù)，如空間碎片受到大氣阻力加速度、地球非球形引力攝動力加速度以及太陽光壓攝動力加速度等；tn表示離散化后的不同時刻；yn是第n個時刻狀態(tài)變量的取值。

龍格-庫塔積分公式為：

其中：h＝tn+1-tn，為積分步長；k為積分式所取的階數(shù)，當取值為4的時候即為4階龍格-庫塔初始化方法；ci,ai,bij均為已知的常數(shù)項。

設已知4個時刻的函數(shù)值分別為fi-3,fi-2,fi-1,fi，則i+1時刻y的近似估計值為

s2，判斷是否達到演化時間，如果未達到，將狀態(tài)更新后的空間碎片運動狀態(tài)和屬性數(shù)據(jù)傳輸給空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)；

s3，s3，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)將所述狀態(tài)更新后的空間碎片運動狀態(tài)和屬性數(shù)據(jù)輸入到碎片解體模型中；所述碎片解體模型進行解體計算，得到空間碎片解體后產(chǎn)生新空間解體碎片的屬性數(shù)據(jù)；并將所述新空間解體碎片的屬性數(shù)據(jù)反饋傳輸?shù)剿隹臻g碎片狀態(tài)推演子系統(tǒng)；

其中，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)，具體用于：

首先，根據(jù)能量定律確定空間碎片解體后產(chǎn)生新空間解體碎片的數(shù)量n和新空間解體碎片特征尺寸l的分布；其中，碎片的特征尺寸l＝(lx+ly+lz)/3，即為碎片3個軸向尺寸的平均。

然后，以新空間解體碎片特征尺寸l作為獨立變量，利用概率分布模型確定新空間解體碎片的面質比a/md、新空間解體碎片的迎風截面積a、新空間解體碎片相對于解體前空間碎片的速度增量δv；

最后，根據(jù)a/md、a確定碎片質量md。

更具體的，所述空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)具體用于：

s3.1，確定解體碎片數(shù)量：

其中，解體碎片區(qū)分為爆炸解體碎片和碰撞解體碎片；

1)對于爆炸解體碎片，根據(jù)能量守恒定律，解體碎片特征尺寸l大于被研究解體碎片最小特征尺寸lc的碎片數(shù)量nf滿足函數(shù)關系：

其中：為解體碎片的最小特征尺寸；

cs是修正系數(shù)，與爆炸類型有關，對于歷史爆炸事件，修正系數(shù)cs的表達式：

其中：修正系數(shù)cs與目標探測直徑門限dssn直接相關，dssn與目標高度h的關系式如下：

其中：h為爆炸解體碎片的高度；

2)對于碰撞解體碎片，解體碎片的數(shù)量和分布與相互碰撞的能量有關。碰撞分為災難性碰撞和非災難性碰撞，兩種碰撞條件下解體碎片的數(shù)量和分布規(guī)律滿足不同的表達式。

一次碰撞是否為災難性碰撞可通過下式判斷：

其中：ep為等效碰撞能量，為災難性碰撞臨界能量。mp為相互碰撞中較小質量；mt為相互碰撞中較大質量；vr是碰撞時兩碎片的相對速度，即碰撞速度；當上式成立時，碰撞為災難性碰撞，反之則為非災難性碰撞。

因此，給出下面的碰撞解體碎片數(shù)量滿足函數(shù)關系：

式中：為兩碰撞目標的等效質量：

其中：

s3.2，確定解體碎片尺寸分布律：

爆炸解體碎片尺寸分布律f：

碰撞解體碎片尺寸分布律f：

計算時，首先確定要研究解體碎片的最小尺寸lc，然后根據(jù)關系式(11)和(15)計算爆炸或解體產(chǎn)生大于lc尺寸碎片的總數(shù)目，再根據(jù)分布律(17)或(18)確定這些碎片具有的尺寸。利用分布律確定解體碎片尺寸時，需要對隨機變量抽樣，可采用反函數(shù)的方法進行隨機變量的抽樣。

s3.3，確定解體碎片面質比參數(shù)：

得到解體碎片的特征尺寸以后，碎片的面質比參數(shù)滿足確定的分布函數(shù)，利用隨機變量采樣法可以得到每一個碎片具有的面質比。

本步驟中，爆炸解體碎片和碰撞解體碎片的碎片面質比參數(shù)滿足同樣的分布函數(shù)；

(1)對于尺寸大于11cm解體碎片，定義解體碎片面質比的對數(shù)值γ＝log(a/md)，解體碎片特征長度的對數(shù)值θ＝log(l)，解體碎片的面質比分布由如下二元正態(tài)分布決定：

p(γ,θ)＝ε(θ)p1(γ)+(1-ε(θ))p2(γ)(19)

其中：p1(γ)和p2(γ)均是正態(tài)分布概率密度函數(shù)：

其中：權重參數(shù)ε(θ)∈[0,1]，是θ的飽和函數(shù)；正態(tài)分布概率密度函數(shù)的均值μi和正態(tài)分布概率密度函數(shù)的方差σi均是θ的函數(shù)；

對于航天器解體碎片，其面質比分布函數(shù)中的參數(shù)由下式?jīng)Q定：

對于火箭上面級解體碎片，其面質比分布函數(shù)的參數(shù)由下式?jīng)Q定

(2)對于尺寸小于8cm的航天器解體碎片和尺寸小于1.7cm的火箭上面級解體碎片，碎片的面質比滿足正態(tài)分布律p1(γ)，即權重函數(shù)ε＝1，

p(γ,θ)＝p1(γ)(23)

此時，分布函數(shù)參數(shù)由下式統(tǒng)一決定：

(3)對于尺寸介于8cm～11cm之間的航天器解體碎片和尺寸介于1.7cm～11cm之間的火箭上面級碎片，碎片的面質比通過隨機采樣的方法確定，即：首先生成服從均勻分布且取值區(qū)間在[0.0,1.0]內(nèi)的隨機變量ζ，然后將ζ與利用(25)計算得到的進行對比，若則利用公式(19)計算面質比參數(shù)，反之則利用公式(23)計算面質比參數(shù)；

步驟3.4，確定解體碎片的速度增量：

解體碎片相對于解體前目標的速度增量滿足正態(tài)分布規(guī)律，定義v＝lg(δv)，γ＝lg(a/md)，則v的分布函數(shù)p(ν)滿足

其中：v為解體碎片的速度；δv為解體碎片的速度增量；

解體碎片的均值μν和解體碎片的方差σν通過下式計算：

σv＝0.4。

s4，所述空間碎片狀態(tài)推演子系統(tǒng)將所述新空間解體碎片的屬性數(shù)據(jù)以及已有的空間碎片的運動狀態(tài)數(shù)據(jù)輸入到數(shù)值積分模型中，再次按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)，如此不斷循環(huán)s2-s4；直到達到演化時間，結束流程。

還包括：

s5，空間碎片運動狀態(tài)三維顯示子系統(tǒng)將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)進行三維顯示，即：利用directx等三維顯示技術，動態(tài)的顯示所有空間碎片在空間中的位置分布狀態(tài)；

具體包括：

步驟s5.1，接收空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)；

步驟s5.2，調(diào)用三維顯示接口函數(shù)(api)，將空間碎片的位置和形狀在三維模型中進行渲染；

步驟s5.3，驅動三維模型，對空間碎片運動狀態(tài)進行動態(tài)更新。

空間碎片的分布狀態(tài)如圖3和圖4所示。

還包括：

s6，在固定時刻或者設定的時刻，空間碎片運動狀態(tài)輸出子系統(tǒng)將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)輸出；其中，所述空間碎片未來運動狀態(tài)包括：空間碎片的位置、速度矢量和面值比、質量、來源屬性參數(shù)；此外，空間碎片運動狀態(tài)輸出子系統(tǒng)還輸出不同高度層內(nèi)空間碎片狀態(tài)數(shù)據(jù)。

還包括：

s7，空間碎片演化結果處理子系統(tǒng)根據(jù)空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)，計算空間碎片狀態(tài)量，進一步計算到空間碎片密度；根據(jù)空間碎片密度得到空間碎片沿軌道高度方向和緯度方向上的分布狀態(tài)。

所述空間碎片密度通過以下步驟計算得到：

步驟s7.1，確定空間碎片在高度區(qū)間[r,r′]內(nèi)的平均空間密度分布其中，r和r′分別為高度區(qū)間的最小地心距和最大地心距；包括以下步驟：

步驟s7.1.1，當目標軌道為橢圓軌道，近地點地心距和遠地點地心距分別為ra和rp，當r′≤ra且r≥rp時，通過下式計算：

其中：a為目標軌道半長軸；

步驟s7.1.2，當目標軌道為橢圓軌道，且r′＞ra或r＜rp時，則令r′＝ra或r＝rp，帶入(28)進行計算，得到

步驟s7.1.3，當目標軌道為圓軌或近圓軌道，即rp→ra∈[r,r′]時，則令r′＝ra或r＝rp，公式(28)變?yōu)楣?29)，通過公式(29)計算

步驟s7.1.4，其他情況，

步驟s7.2，確定空間碎片在緯度區(qū)間[β,β′]內(nèi)的平均空間密度分布函數(shù)其中，β和β′分別是緯度區(qū)間的最小值和最大值；

具體的，若β′＞0，則直接執(zhí)行步驟s7.2.1；

若β′≤0，則令β′＝-β，β＝-β′，均轉換為關于原點對稱的正值區(qū)間后，再執(zhí)行步驟s7.2.1；

步驟s7.2.1，碎片軌道傾角記為i，當β′≤i且β≥-i時

步驟s7.2.2，當β≤-i且β′≤i時

步驟s7.2.3，當β′≥i且-i≤β≤i時

步驟s7.2.4，當β′≥i且β≤-i時

步驟s7.2.5，當β＞i時，

本發(fā)明還提供一種大規(guī)?？臻g碎片分布演化的數(shù)值計算系統(tǒng)，包括空間碎片狀態(tài)推演子系統(tǒng)、空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)、空間碎片運動狀態(tài)三維顯示子系統(tǒng)、空間碎片運動狀態(tài)輸出子系統(tǒng)和空間碎片演化結果處理子系統(tǒng)。下在對這5個子系統(tǒng)詳細介紹：

(1)空間碎片狀態(tài)推演子系統(tǒng)

空間碎片狀態(tài)推演子系統(tǒng)，用于獲取空間碎片的初始運動狀態(tài)數(shù)據(jù)以及空間碎片受到的空間攝動力，并將所述初始運動狀態(tài)數(shù)據(jù)以及所述空間攝動力輸入到數(shù)值積分模型中；所述數(shù)值積分模型按照給定的時間步長進行積分推演，得到空間碎片未來運動狀態(tài)數(shù)據(jù)；然后，采用所述空間碎片未來運動狀態(tài)數(shù)據(jù)更新空間碎片當前運動狀態(tài)數(shù)據(jù)，得到空間碎片新運動狀態(tài)數(shù)據(jù)。

也就是說，空間碎片狀態(tài)推演子系統(tǒng)以空間碎片的初始運動狀態(tài)為輸入，綜合考慮空間碎片受到的大氣阻力、非球形引力、太陽光壓、三體引力等攝動力的作用，采用高效穩(wěn)定的積分方法，積分推演碎片的未來的運動狀態(tài)。

其中，在計算大氣阻力的過程中，要首先確定高層大氣密度。太陽輻射作用使大氣密度呈現(xiàn)周日變化、季節(jié)變化等周期性變化規(guī)律；太陽活動會導致大氣密度發(fā)生劇烈變化，嚴重影響大氣阻力對目標的衰減作用。計算系統(tǒng)選取了harris-priester大氣密度模型，該模型能夠反映太陽活動對大氣密度的影響，模型復雜性低，對計算能力要求不高，適宜進行大規(guī)模的長期演化計算。

計算空間碎片受到的非球形攝動時，計算系統(tǒng)對其量級進行了分析，并進行了合理的取舍。地球扁率j2項是主要的攝動項，導致空間碎片半長軸和升交點赤經(jīng)周期性變化，影響著物體沿經(jīng)度方向上的分布規(guī)律；在扁率j2項和j2,2項影響下，geo軌道帶物體在地球赤道短軸上空，即東經(jīng)75°和西經(jīng)105°附近，呈現(xiàn)聚集分布規(guī)律；j3項和j4項的量級與j2,2項量級相當，均為10^-6，計算空間碎片受到的非球形攝動時也加以考慮。

計算系統(tǒng)中的積分方法是經(jīng)過優(yōu)選確定的。對空間碎片的運動進行動力學積分時，要求積分算法精度高且具有有一定穩(wěn)定性，避免長期計算時結果發(fā)散；由于空間碎片數(shù)量規(guī)模大，需要開發(fā)高效的并行計算方法。系統(tǒng)采用4階阿達姆斯預估矯正方法(adams-bashforth)進行積分，該積分方法利用4階龍格-庫塔方法進行參數(shù)初始化，充分利用歷史計算結果，每步積分僅計算一次攝動力，減少積分計算量；采用4階積分方法，比高階積分具有跟高的計算穩(wěn)定性；模型基于mpi的并行計算框架，利用多進程和多核處理技術，設計了高效的并行演化系統(tǒng)。

(2)空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)

空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)，用于將所述空間碎片新運動狀態(tài)數(shù)據(jù)輸入到碎片解體模型中；所述碎片解體模型進行解體計算，得到空間碎片解體后產(chǎn)生新空間解體碎片的屬性數(shù)據(jù)；并將所述新空間解體碎片的屬性數(shù)據(jù)反饋傳輸?shù)剿隹臻g碎片狀態(tài)推演子系統(tǒng)。

空間碎片碰撞/爆炸解體產(chǎn)生新碎片子系統(tǒng)的功能包括兩個方面：首先確定空間碎片之間發(fā)生碰撞的概率，然后根據(jù)碎片解體模型生成大量解體碎片。然后，根據(jù)能量定律確定解體后產(chǎn)生碎片的數(shù)量n和碎片特征尺寸l的分布；然后以特征尺寸l作為獨立變量，利用概率分布模型確定碎片的面質比a/md、有效橫截面積a和速度增量δv；最后根據(jù)a/md、a確定碎片質量md。

(3)空間碎片運動狀態(tài)三維顯示子系統(tǒng)

空間碎片運動狀態(tài)三維顯示子系統(tǒng)，用于利用directx等三維顯示技術，將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)數(shù)據(jù)進行三維顯示，子系統(tǒng)能夠根據(jù)物體的推演狀態(tài)，動態(tài)的顯示所有物體在空間中的位置分布狀態(tài)。

(4)空間碎片運動狀態(tài)輸出子系統(tǒng)

空間碎片運動狀態(tài)輸出子系統(tǒng)，用于在固定時刻或者設定的時刻，將所述空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)輸出，具體的，輸出包塊解體碎片在內(nèi)的當前所有碎片的位置、速度矢量和物體的面值比等屬性參數(shù)；此外，空間碎片運動狀態(tài)輸出子系統(tǒng)還輸出不同高度層內(nèi)空間碎片狀態(tài)數(shù)據(jù)；

(5)空間碎片演化結果處理子系統(tǒng)

空間碎片演化結果處理子系統(tǒng)，用于根據(jù)空間碎片狀態(tài)推演子系統(tǒng)推演到的空間碎片未來運動狀態(tài)，分析碎片沿軌道高度方向、緯度方向等空間維度方向上的變化情況，分析碎片規(guī)模隨時間的變化規(guī)律。計算空間碎片狀態(tài)量，進一步計算到空間碎片密度；根據(jù)空間碎片密度得到空間碎片沿軌道高度方向和緯度方向上的分布狀態(tài)。

本發(fā)明提供的大規(guī)?？臻g碎片分布演化數(shù)值計算方法及系統(tǒng)，能夠根據(jù)碎片的初始運動狀態(tài)，綜合考慮大氣阻力、地球非球形攝動力以及太陽/月球三體引力等復雜攝動力的作用，碎片相互碰撞產(chǎn)生新碎片的作用，對碎片的分布狀態(tài)進行長期演化計算。系統(tǒng)能夠三維動態(tài)的顯示空間碎片在軌分布狀態(tài)隨時間的變化，能夠對碎片分布狀態(tài)的演化計算結果進行統(tǒng)計分析，得到空間碎片在不同軌道高度、不同緯度等方向上的分布規(guī)律，從而為航天器在軌運行碰撞規(guī)避、空間碎片減緩措施制定以及確?？臻g資源可持續(xù)利用的重要的技術保障。

以上所述僅是本發(fā)明的優(yōu)選實施方式，應當指出，對于本技術領域的普通技術人員來說，在不脫離本發(fā)明原理的前提下，還可以做出若干改進和潤飾，這些改進和潤飾也應視本發(fā)明的保護范圍。