凱斯勒臨界線：以平行化 Cowell 公式模擬軌道碎片級聯

封存文章

本文保留一個以合成碎裂情境探索軌道碎片動力學的教學模型。圖中數值不是對真實碰撞事件的重建，亦未用正式 conjunction-assessment 資料驗證，因此不應視為實際碰撞風險預測。

在低地球軌道，細小物件也可具有極高相對動能。2009 年 Cosmos 2251 與 Iridium 33 相撞，是第一宗兩枚完整人造衞星之間的超高速碰撞，並產生數以千計可追蹤碎片。這類事件說明了凱斯勒效應的核心機制：物件密度上升令碰撞增加，而新碎片又進一步提高日後碰撞率。

這個專案提出一個較窄、可計算的問題：

地球的 $J_2$ 扁率、大氣阻力與太陽輻射壓，如何改變一次合成碎裂事件在短時間內的幾何形狀與軌道衰減？

模型由 Cartesian state 直接積分運動方程，並把每件碎片的傳播分派至不同程序。文章分為物理模型、計算架構、合成情境與結果判讀四部分。

一、物理引擎

理想二體問題產生 Kepler 橢圓，但真實低軌物件還受到地球非球形重力、稀薄大氣與太陽光壓影響。這裡使用 Cowell 公式，直接積分

\ddot{\mathbf r} = \mathbf a_{\mathrm{grav}} +\mathbf a_{J_2} +\mathbf a_{\mathrm{drag}} +\mathbf a_{\mathrm{SRP}}.

1. 中心重力

把地球視為點質量時，

\mathbf a_{\mathrm{grav}} =-\frac{\mu}{\lVert\mathbf r\rVert^3}\mathbf r,

其中 $\mu=GM_\oplus\approx3.986\times10^{14}\ {\rm m^3\,s^{-2}}$ 。只有這一項時，系統就是標準二體問題，並在連續模型中守恆能量與角動量。

2. $J_2$ 擾動

地球因自轉而呈扁球形。重力勢的首個主要非球形修正是第二帶諧係數 $J_2$ ：

U(r,\phi) =\frac{\mu}{r}\left[ 1-J_2\left(\frac{R_E}{r}\right)^2 \frac{3\sin^2\phi-1}{2} \right].

由勢函數取梯度可得程式所用的加速度。以 Cartesian 分量表示，

\mathbf a_{J_2} = \frac{3\mu J_2R_E^2}{2r^5} \begin{bmatrix} x(5z^2/r^2-1)\\ y(5z^2/r^2-1)\\ z(5z^2/r^2-3) \end{bmatrix}.

它會引起升交點進動；對近圓軌道，一階近似為

\dot\Omega \approx -\frac32nJ_2 \left(\frac{R_E}{a}\right)^2 \frac{\cos i}{(1-e^2)^2}.

碎裂後各碎片的半長軸、偏心率與傾角稍有不同，因而以不同速率進動。原本集中在一個位置的碎片雲會先沿軌道剪切，再逐漸在軌道面方向展寬。

3. 大氣阻力

低軌並非真空。相對於共同旋轉大氣的阻力可寫成

\mathbf a_{\mathrm{drag}} =-\frac12\rho(h)C_d\frac{A}{m} \lVert\mathbf v_{\mathrm{rel}}\rVert \mathbf v_{\mathrm{rel}},

其中

\mathbf v_{\mathrm{rel}} =\mathbf v_{\mathrm{ECI}} -\boldsymbol\omega_\oplus\times\mathbf r.

$C_dA/m$ 控制碎片對阻力的敏感度：面積質量比高的薄片會較快失去軌道能量，緻密碎片則維持較久。教學模型採用靜態指數密度

\rho(h)=\rho_0\exp\left[-\frac{h-h_0}{H}\right].

這是重要簡化。真實熱層密度會隨太陽通量、地磁活動、地方時與緯度改變；正式研究應採用如 NRLMSISE-00 的模型並傳播密度不確定性。

4. 太陽輻射壓

光子攜帶動量。在 1 AU 附近，太陽輻射壓約為 $4.56\times10^{-6}\ {\rm N\,m^{-2}}$ 。簡化模型使用

\mathbf a_{\mathrm{SRP}} =-P_\odot C_R\frac{A}{m} \hat{\mathbf r}_{\odot\rightarrow\mathrm{sat}}\nu,

其中 $\nu$ 是日照／地影函數。原實驗把 $\nu$ 固定為 1，即假設碎片全程受光；這不是實際軌道的忠實描述，只用來隔離高面積質量比碎片對 SRP 的敏感度。

二、數值積分與平行化

1. RK45

每件碎片以 scipy.integrate.solve_ivp(method="RK45") 傳播。 RK45 是具內嵌誤差估計的顯式 Runge–Kutta 5(4) 法，可根據局部截斷誤差自動調整步長。設定 rtol=1e-5、atol=1e-8 並不自動「保證物理精確」：仍需以更嚴格 tolerance 重跑、比較守恆量，以及檢查輸出取樣是否足以捕捉近地點變化。

2. 多程序 Monte Carlo

在此模型中，碎片之間不計互相引力或碰撞，因此每條軌跡可獨立傳播，是典型的 embarrassingly parallel 問題。Python thread 受 GIL 限制，但多程序可以讓每個 worker 擁有獨立 interpreter：

import multiprocessing as mp

def worker(initial_state):
    return propagate(initial_state)

if __name__ == "__main__":
    debris = generate_debris(n=1_000)
    with mp.Pool(processes=mp.cpu_count()) as pool:
        trajectories = pool.map(worker, debris)

實際加速不會必然等於 CPU 核心數：序列化、程序啟動、記憶體頻寬與每條軌跡的不均勻工作量都會造成 overhead。因此原文所述「十二核心線性加速」應視為特定電腦上的觀察，而非一般保證。

三、合成碎裂情境

實驗設定一件位於近太陽同步軌道的母體：

質量：2,000 kg；
高度：800 km；
傾角： $98^\circ$ 。

另設一個高度 780 km、傾角 $45^\circ$ 的合成目標，讓它與碎片軌道面交會。碎裂模型只借用 NASA Standard Breakup Model 的若干概念，並不是完整實作：

速度增量 $\Delta v$ 以零均值 Gaussian 分布取樣，標準差為 $150\ {\rm m\,s^{-1}}$ ；
方向在球面上均勻取樣；
面積質量比 $A/m$ 取 log-normal 分布。

這些假設可製造低、高能量和不同阻力敏感度的碎片，但沒有由母體材質、撞擊能量或實測 fragment-size distribution 校準。

四、結果如何判讀

1. 碎片雲的剪切

碎片雲隨時間演化

圖 1　合成碎片雲的短期演化；畫面主要呈現不同軌道週期造成的相位剪切。

碎裂後，沿飛行方向獲得正向 $\Delta v$ 的碎片通常進入較高能量、較長週期的軌道；獲得反向 $\Delta v$ 的碎片則可能進入較短週期軌道。不同平均運動令雲團沿軌道拉長，最後繞成環狀。 $J_2$ 造成的差異進動在更長時段把薄環擴成較厚的三維帶。

不過 24 小時模擬只足以展示初期剪切；要聲稱形成全球殼層，需要更長時間積分、更完整擾動模型與數值誤差檢查。

2. Gabbard 圖

軌道碎片分析儀表板

圖 2　合成結果儀表板；左下為軌道週期對近／遠地點高度的 Gabbard 表示。

Gabbard 圖通常把軌道週期與近地點、遠地點高度配對。高能碎裂會產生近似 X 形分支：

頂點接近母體原有週期與高度；
失去軌道能量的碎片具有較短週期，部分近地點下降；
獲得能量的碎片具有較長週期，遠地點上升。

低近地點碎片受較強阻力，長期會較快衰減，但在一天尺度內要把圖形分支的缺口直接歸因於 de-orbit 仍需逐件核對再入條件，而不能只靠視覺判斷。

3. 「警戒區」並不等於碰撞概率

儀表板亦計算碎片與合成目標的最近距離分布。設定 50 km 門檻可用來篩選近接事件，但它不是 hard-body radius，也不代表碰撞。簡單 flux 模型

P_c=1-\exp(-FAT)

只有在通量 $F$ 、有效截面 $A$ 、暴露時間 $T$ 與 Poisson 假設皆合理時才有意義。正式 conjunction assessment 需要相對狀態、兩物體尺寸、位置協方差與誤差橢球；本實驗沒有這些資料，所以不能支持「撞擊只是時間問題」之類斷言。

4. 軌道衰減

軌道高度包絡線

圖 3　部分碎片的高度時間序列，用以比較面積質量比與阻力敏感度。

高度的週期波動主要反映近地點與遠地點。若平均軌道能量隨時間下降，波動包絡線會下移；高 $A/m$ 碎片一般較明顯。但指數大氣加上短時段模擬只能示範機制，不能可靠估計數十年至數百年的軌道壽命。

結論

這個合成實驗把四個概念連成一條可計算鏈：

$\Delta v$ 改變碎片的軌道能量與週期；
Keplerian shear 把局部雲團沿軌道拉長；
$J_2$ 差異進動令軌道面逐漸分散；
大氣阻力按面積質量比選擇性移除低近地點碎片。

它適合用來學習 Cowell 傳播、擾動力、Monte Carlo 初始條件與平行計算，但仍缺少高保真碎裂分布、日地幾何、地影、時變大氣、碰撞協方差與長期人口演化。因此，最穩妥的結論不是預測某顆衞星何時被撞，而是展示單一碎裂事件如何把局部初始條件轉化為分布廣泛、壽命差異很大的軌道風險來源。

參考資料

Kessler, D. J., & Cour-Palais, B. G. (1978). Collision frequency of artificial satellites: The creation of a debris belt. Journal of Geophysical Research, 83(A6), 2637–2646.
Vallado, D. A. (2013). Fundamentals of Astrodynamics and Applications (4th ed.). Microcosm Press.
NASA Orbital Debris Program Office. NASA Standard Breakup Model and orbital-debris technical resources.
Montenbruck, O., & Gill, E. (2000). Satellite Orbits: Models, Methods and Applications. Springer.