0 引 言
近年来，随着各类低轨星座开始大规模部署，低轨卫星数量急剧增长。星座卫星失效后若长期留轨，存在与其他失效卫星、箭体等撞击而解体的严重风险。卫星在轨解体将产生大量碎片，并在动力学作用下不断演化[1]：早期，解体碎片之间的相对距离较近，但由于轨道周期差异，碎片的平近点角逐渐分散，沿轨迹方向不断拉伸，最终形成碎片环；中期，主要在地球非球形引力作用下，各碎片的升交点、拱点逐渐分散；当不同碎片的升交点、拱点均可视为随机均匀分布时，碎片云进入长期演化阶段。由于解体事件具有突发性，对解体碎片的探测识别和定轨编目一般会滞后于解体时刻[2]。同时，受监测装备的能力限制，低轨区域通常只能实现10 cm以上碎片的稳定探测编目。因此，解体初期碎片云的监测预警难度极大，准确评估其对其他星座卫星的短期碰撞风险，可为突发解体事件的处置应对提供重要支撑。
在解体碎片云的危害评估问题研究中，由于解体碎片的数量规模极大且大部分无法被编目，相关方法一般借鉴气体动力学理论，引入空间密度的概念来描述碎片云分布，并以碰撞通量等进行碰撞风险量化评估。按照空间密度计算方法的不同，现有方法主要分为以下三类：
（1）离散方法。其基本思想是将解体碎片视为单独的个体，分别对所有碎片进行轨道预报，以获得碎片云在不同时刻的空间分布，进而评估其与航天器的碰撞风险。计算所需的解体碎片通过对解体模型进行随机抽样来生成。目前，国外的ORDEM、MASTER及国内的SDEEM系列空间碎片环境工程模型均采用离散方法进行解体碎片云的演化计算。例如，文献[3-4]分别利用SDEEM模型及MASTER模型开展了星座卫星解体事件对其他卫星碰撞风险分析；文献[5]则直接利用离散方法分析了巨型星座中处于升轨期的卫星突发解体事件造成的碰撞风险等。由于离散方法需通过多次解体碎片的蒙特卡罗抽样生成和预报计算以获得可信评估结果，将耗费大量计算时间，极大限制该类方法的应用范围。
（2）演化方法。其基本思想是将碎片云视为整体，利用概率密度函数描述其空间分布状态，并通过建立演化方程实现对碎片云空间密度分布的演化计算，在此基础上评估不同时刻航天器的碰撞风险。1976年，Heard[6]首次提出可以将解体碎片云视为由大量互不接触的粒子构成的连续体，并给出了圆轨道、只考虑二体运动的碎片云空间密度演化函数。随后，相关学者开展了一系列拓展研究，如McInnes[7]利用流体连续性方程对解体碎片云的密度演化进行建模描述，Letizia等[8]在碎片云的密度演化方程中进一步考虑了指数大气密度模型的影响，并将其应用于解体参数对碰撞风险评估结果的敏感性分析中，Frey等[9]提出一种适用于不同碎片云密度演化计算方法的通用碰撞风险评估框架，Wen等[10]和金紫涵等[11]则针对碎片云中期演化的动力学特征，建立了轨道根数形式的碎片云密度演化计算方法，并开展了低轨星座中卫星解体数月内对星座卫星的碰撞风险分析等。但是，上述方法一般需要对碎片分布和轨道动力学模型进行一定简化。解体初期，碎片云的密度分布极不均匀且变化剧烈，简化处理将极大影响演化计算与危害评估结果的可信度。
目前，在以星座为分析对象的解体事件危害评估研究中，相关工作主要集中于解体碎片云的中、长期演化规律与撞击风险分析，短期危害评估相关工作尚不多见。为此，本文以低轨星座中某卫星受碎片撞击而发生突发解体事件为背景，首先建立基于边值问题的解体碎片云短期碰撞风险评估算法，并利用蒙特卡罗方法对算法进行验证。在此基础上，针对给定的星座卫星解体场景，开展解体碎片云的短期碰撞风险评估，获得不同星座卫星受不同尺寸解体碎片的短期碰撞风险变化规律。
1 碰撞风险评估算法
如图1 所示，记解体事件发生时刻为${t_1}$ ，解体点的地心矢径为${{\boldsymbol{r}}_1}$ 。在${t_2}$ 时刻，待评估卫星的地心矢径为${{\boldsymbol{r}}_2}$ 。
危害评估的关键在于确定在该时刻能与卫星相撞的解体碎片，本文将这类碎片定义为危险碎片。这类碎片可视为从解体点出发，经过${t_2} - {t_1}$ 时间后到达卫星所在位置。由于起始位置、终端位置、飞行时间均为已知条件，通过求解两点边值问题（兰伯特问题），可以精确计算出危险碎片在解体时刻所需具备的飞散速度。根据解体碎片飞散速度的概率分布，即可得到解体事件产生危险碎片的概率，并进一步推导得到解体事件对卫星的碰撞风险。碰撞风险评估算法流程如图2 所示。
1.1 危险解体碎片的飞散速度求解
本文针对任意给定时刻${t_2}$ ，危险碎片均可被视为在${t_1}$ 时刻从解体位置${{\boldsymbol{r}}_1}$ 出发，在${t_2}$ 时刻飞行至卫星所在位置${{\boldsymbol{r}}_2}$ 。通过求解兰伯特问题，可以确定危险碎片的出发速度，该速度相对于解体时刻母体速度的增量即为碎片应该具备的飞散速度。
兰伯特问题的求解是轨道力学中的经典问题[19]。在二体轨道假设下，若某空间目标经转移时间$\Delta t$ 从位置矢量${{\boldsymbol{r}}_1}$ 转移到${{\boldsymbol{r}}_2}$ ，则转移时间与转移轨道半长轴$a$ 、距离和${{\boldsymbol{r}}_1} + {{\boldsymbol{r}}_2}$ 、弦长$c$ 之间满足兰伯特飞行时间方程。通过求解该方程，即可确定转移轨道。通常情况下，当$\Delta t$ 较小时，转移轨道的解具有唯一性。当$\Delta t$ 较大时，转移轨道可能包含多个完整圈次，对应不同的转移轨道解。针对给定的转移时间$\Delta t$ ，可能的最大飞行圈数为[20]：
式中，[·]表示取整，$ c $ 为弦长。满足二体运动规律的可能转移轨道有$2{N_{\max }} + 1$ 个。
由于本文只关注碎片云的短期撞击风险，碎片轨道转移时间相对较短。因此，求解转移轨道时，只需考虑J2项摄动这一主要摄动力的影响。在进行危险碎片的飞散速度求解时，针对任意可能的转移圈次，首先利用Izzo[24]提出的二体兰伯特算法求解碎片在${t_1}$ 时刻的出发速度。随后，以此速度为初值，利用考虑J2项摄动的碎片运动方程和牛顿迭代方法，对碎片的出发速度进行迭代修正。记第k次迭代得到的碎片出发速度为${\boldsymbol{v}}_1^k$ ，摄动运动模型给出的碎片终端实际位置为${\boldsymbol{r}}_2^k$ ，则$ {\boldsymbol{v}}_1^{k + 1} $ 由式（2）迭代更新：
式中，$ {\dfrac{{\partial {\boldsymbol{r}}_2^k}}{{\partial {\boldsymbol{v}}_1^k}}} $ 利用有限差分方法求解，当$ \left\| {{{\boldsymbol{r}}_2} - {\boldsymbol{r}}_2^k} \right\| $ 足够小时，迭代终止。
1.2 特定尺寸危险解体碎片的概率密度计算
在CSBM模型中，以解体碎片的质量作为独立变量，碎片质量分布为：
式中，$ C{N_{{m_f}}} $ 为质量≥${m_f}$ 的解体碎片总数；$ {t_i}(i = 0,\; \cdots ,\;3) $ 为常数；$ M $ 为撞击通道总质量；$\varphi $ 为由撞击通道内能量密度决定的参数；$V$ 为相对撞击速度；$\rho $ 为解体航天器名义密度。
当解体产生的碎片质量下限${m_{{\mathrm{lb}}}}$ 和上限${m_{{\mathrm{ub}}}}$ 给定时，解体碎片总数可由式(3)确定，进而可得解体碎片质量分布的概率密度函数为：
式中，$ c = {t_2}\lg V + {t_3}\lg \rho $ 。
为便于后续推导，利用一维随机变量X与$Y = g(X)$ 之间的概率密度换元公式：
可得解体碎片对数质量$\xi = \lg ({m_f})$ 的概率密度函数为：
对于质量为$\xi $ 的解体碎片，CSBM模型认为其对数形式的飞散速度值$\upsilon = \lg (\Delta v)$ 服从正态分布，飞散速度的方向则满足完全随机分布，表示为条件概率密度的形式为：
式中，$ \mu = - 0.030\;49{\xi ^2} - 0.185\;65\xi + 2.626\;82 $ ，$ \sigma = 0.4 $ 。
解体碎片的质量$\xi $ 给定时，其对数形式的面质比$\chi = \lg {A_M}$ 则服从一个三峰正态分布，表示为条件概率密度形式为：
式中，${\alpha _i}$ 、${\mu _i}$ 和${\sigma _i}$ 均为由$\xi $ 决定的参数。
利用CSBM模型中碎片特征尺寸$\lambda = \lg L$ 与面质比$\chi $ 之间的变换公式：
可得碎片特征尺寸的条件概率密度函数为：
结合式（6）~式（7）和式（10），可得碎片飞散速度值与特征尺寸的联合概率密度函数为：
利用概率密度换元公式，可得飞散速度值$\Delta {\boldsymbol{v}}$ 与特征尺寸${L_c}$ 的联合概率密度函数为：
由飞散速度方向的全向均匀分布假设，可以进一步得到飞散速度矢量$ \Delta {\boldsymbol{v}} $ 与特征尺寸$L$ 的联合概率密度为：
根据式（13），即可利用特定尺寸危险碎片应该具备的飞散速度矢量，求出其对应的概率密度。
1.3 碎片云对卫星的碰撞风险描述
在空间碎片环境的碰撞风险评估中，常用通量的概念描述碎片环境对卫星撞击风险。对于已知撞击截面积为${A_c}$ 的空间物体，以恒定速率v在空间碎片密度为$\rho $ 的轨道上运行，其在有限时间$\Delta t$ 内的平均碰撞次数c为[29]：
对于在${t_1}$ 时刻具有飞散速度矢量$\Delta {\boldsymbol{v}}$ 和特征尺寸$L$ 的危险解体碎片，式(13)给出的概率密度可视为解体碎片在飞散速度空间中的概率密度。其对应于在${t_2}$ 时刻的解体碎片位置分布中，位置${{\boldsymbol{r}}_2}$ 处的特征尺寸为$L$ 的碎片概率密度表示为：
位置${{\boldsymbol{r}}_2}$ 处的特征尺寸在$[{L_0},{L_1}]$ 区间的碎片概率密度为：
将卫星外形简化为半径R的球形，则根据式(14)，卫星在${t_2}$ 时刻受到的在$[{L_0},{L_1}]$ 区间的任一尺寸解体碎片瞬时撞击率可表示为：
式中，N为从解体点${{\boldsymbol{r}}_1}$ 转移到位置${{\boldsymbol{r}}_2}$ 的两点边值问题的可行解个数，$ {A_s} = {\text{π }}{R^2} $ 为截面积；$ v_{{\mathrm{sf}}}^i $ 为${t_2}$ 时刻解体碎片与卫星的相对速度，在两点边值问题求解时一并给出。
记解体事件产生的碎片总数为${N_s}$ ，则卫星在${t_2}$ 时刻受到的尺寸在$[{L_0},{L_1}]$ 区间的解体碎片瞬时撞击数可表示为：
对式(18)从解体时刻${t_1}$ 积分至${t_2}$ ，即可得到解体碎片在${t_2}$ 时刻对卫星的累积撞击数：
由式（17）~式（19），即可对碎片云与卫星的碰撞风险进行量化描述。
2 算法验证
利用蒙特卡罗方法可以对本文提出的碰撞风险评估算法进行验证。蒙特卡罗方法的基本思想是通过对造成非确定性事件的随机参数进行依概率的随机采样，生成大量随机参数样本，再通过对这些样本对应的事件结果进行统计分析，从而对随机参数引起的不确定性进行量化。理论上，当样本量足够大时，蒙特卡罗方法的统计结果可以很好地逼近真实情况。
用于算法验证的评估场景设定如下：某卫星在${T_0}$ 时刻发生解体，并产生大量解体碎片，解体后的第8 580 s对某航天器造成碰撞风险。解体卫星及航天器的轨道参数如表1 所示。
表 1 解体卫星及航天器的轨道参数（J2000）Table 1 Orbital parameters of the breakup satellite and the spacecraft （J2000） |
| 对象 | 时间 | 轨道参数/（km、km/s） |
| 解体卫星 | T0 | 位置：[−765.4, −7 209.1, −1 069.8] |
| 速度：[1.27, 0.93, −7.20] | ||
| 航天器 | T0+8 580 s | 位置：[1 631.4, 5 731.1, −4 265.5] |
| 速度：[−2.52, −3.66, −5.88] |
考虑到在蒙特卡罗方法中，当航天器与碎片的碰撞概率较低时，需采集大量碎片样本方能获得较为准确的计算结果。为便于验证并提高计算效率，本文暂不考虑解体碎片的尺寸分布特性，并假设解体碎片飞散速度的3个分量满足式(20)所示的正态分布，以保证碎片与航天器的碰撞概率达到足够高的水平。
针对上述场景，分别利用本文算法及蒙特卡罗方法，计算任意解体碎片在T0+8580 s时对航天器的瞬时撞击数，并对结果进行对比分析。在计算过程中，假设航天器为半径5 m的球形。
蒙特卡罗方法的具体计算流程如下：首先，根据碎片飞散速度的概率密度函数，通过随机抽样产生${N_m}$ 颗碎片；接着，利用J2项摄动运动方程对每个碎片进行轨道外推，并判断其在T0+8 580 s时是否击中航天器；最后，根据击中航天器的各碎片交会参数，计算单颗碎片与航天器的瞬时撞击数。设击中航天器的碎片总数为${N_{{\text{hit}}}}$ ，其中第$i$ 颗碎片与航天器交会的相对速度为$v_{{\text{sf}}}^i$ ，则蒙特卡罗方法计算单颗解体碎片的瞬时撞击数为：
蒙特卡罗方法和本文算法的计算结果如图3 所示：
3 星座解体事件的短期危害评估
3.1 解体场景建模
在本文设定的星座卫星解体场景中，低轨星座采用Walker星座构型。其主要参数包括卫星总数N、轨道面个数P、相位因子F、轨道高度h、轨道倾角i，具体取值如表2 所示。
表 2 星座参数Table 2 Constellation parameters |
| 参数 | N | P | F | h/km | i/(°) |
| 取值 | 84 | 12 | 1 | 950 | 81 |
星座各轨道面按升交点均匀分布，轨道面内各卫星相位等间隔分布，且相邻轨道面对应序号的卫星之间相位差为：
记${\text{Sa}}{{\text{t}}_{i - j}}$ 为星座中第i个轨道面的第j颗卫星，${\text{Sa}}{{\text{t}}_{1 - 1}}$ 为基准卫星，其在${T_0}$ 时刻的轨道参数如表3 所示：
表 3 基准卫星初始轨道参数Table 3 Initial orbital parameters of the baseline satellite |
| 参数 | a/km | e | i/(°) | f/(°) | ||
| 取值 | 0 | 81.002 | 82.601 | 153.436 | 35.064 |
为便于考察解体事件对星座卫星的影响，假设星座构型保持稳定，各卫星的轨道运动采用二体轨道预报模型进行仿真计算。
在${T_0}$ 时刻，假设星座中的已报废卫星${\text{Sa}}{{\text{t}}_{1 - 1}}$ 受空间碎片撞击而发生解体，解体相关参数如表4 所示。
表 4 解体相关参数Table 4 Parameters related to the breakup event |
| 参数/单位 | 取值 |
| 解体卫星质量/kg | 730 |
| 空间碎片质量/kg | 3 |
| 卫星名义密度/g·cm−3 | 0.08 |
| 相对撞击速度/km·s−1 | 10.5 |
| 卫星的撞击通道质量比 | 0.16 |
利用CSBM模型对解体事件进行建模，令式(3)中的$ C{N_{{{{m}}_{{f}}}}} = 1 $ ，得到解体碎片的质量上限为${m_{{\text{ub}}}} = 20.19 \;{\text{kg}}$ 。同时，考虑到毫米级以下碎片对卫星的碰撞危害相对有限，参照φ为1 mm的铝球质量，将解体碎片的质量下限设置为${m_{{\text{lb}}}} = 1 \times {10^{ - 6}}\;{\text{kg}}$ ，据此计算，解体产生的碎片总数约为299.06万个。碎片飞散速度和特征尺寸的联合概率密度分布如图4 所示。
针对上述解体场景，利用本文算法开展解体事件的短期危害评估，并分析星座卫星的撞击风险变化规律。考虑到低轨卫星解体碎片云的短期演化阶段通常持续约10个轨道周期[30]，评估时间设置为解体后的12个轨道周期内。在本文算法中，所有卫星外形均等效处理为半径5 m的球形。
3.2 共面卫星的碰撞风险分析
以与解体卫星相邻的${\text{Sa}}{{\text{t}}_{1 - 2}}$ 卫星为例，其受不同尺寸级别碎片撞击风险情况如图5 所示。
由星座构型特征可知，${\text{Sa}}{{\text{t}}_{1 - {{j}}}}({{j}} = 2, \cdots ,7)$ 在同一个轨道面上按等相位差分布。因此，${\text{Sa}}{{\text{t}}_{1 - {{j}}}}({{j}} = 7, \cdots ,2)$ 将依次穿越${\text{Sa}}{{\text{t}}_{1 - 1}}$ 产生的碎片云，导致其受碎片的瞬时碰撞风险依次达到最大，如图8 所示。
3.3 异面卫星的碰撞风险分析
首先，以与解体卫星异面的3颗卫星${\text{Sa}}{{\text{t}}_{2 - 1}}$ 、${\text{Sa}}{{\text{t}}_{4 - 5}}$ 、${\text{Sa}}{{\text{t}}_{6 - 7}}$ 为例，其在解体事件发生后短期内受解体碎片的瞬时撞击风险变化情况如图9 所示。图9 中，黑色散点表示解体碎片对各卫星的瞬时撞击数，对应于左侧纵坐标；蓝色实线表示各卫星与解体点的伪相对距离，以地球赤道平均半径${R_E} = 6\;378.137\;{\text{km}}$ 为度量单位，对应于右侧纵坐标。需要说明的是，该相对距离是在假设解体点仍沿着解体卫星轨道继续运动的情况下计算得到的，用于评估卫星与解体碎片云间的接近程度。
由图9 可知，${\text{Sa}}{{\text{t}}_{2 - 1}}$ 和${\text{Sa}}{{\text{t}}_{6 - 7}}$ 受解体碎片的最大瞬时碰撞数接近10−6量级，与解体点的最小伪相对距离小于0.3RE，且瞬时碰撞数呈周期性变化。在每个变化周期内，瞬时碰撞数的最大值对应伪相对距离的最小值，而后续变化周期中最大瞬时碰撞数呈逐渐变小的趋势。${\text{Sa}}{{\text{t}}_{4 - 5}}$ 受解体碎片的最大瞬时碰撞数则在10−8量级，与解体点的伪相对距离大于1.9RE，瞬时碰撞数虽也呈周期性变化特征，但在每个周期中最大瞬时撞击数呈逐渐增大的趋势。
进一步以${\text{Sa}}{{\text{t}}_{2 - 1}}$ 为例，分析瞬时碰撞数与伪相对速度的变化关系，如图10 所示。可以看出，瞬时碰撞率的最大值与伪相对速度的最大值相对应。
上述分析表明，对于${\text{Sa}}{{\text{t}}_{2 - 1}}$ 和${\text{Sa}}{{\text{t}}_{6 - 7}}$ 这类与解体卫星存在高速接近的异面卫星而言，其受解体碎片的最大瞬时撞击风险出现在卫星和碎片云的周期性高速接近的时刻，且每次接近所造成的最大撞击风险逐渐降低。对于${\text{Sa}}{{\text{t}}_{4 - 5}}$ 这类始终与解体卫星保持较远相对距离的异面卫星而言，其主要受到碎片云逐渐扩散带来的撞击风险，虽然最大撞击风险随时间推移而逐渐升高，但整体风险等级明显较低。
总体而言，碎片云的撞击风险主要集中在共面卫星及与其高速交会的少量异面卫星上，且在解体后迅速达到10−4量级。相比之下，其他星座卫星的撞击风险则随时间推移而逐渐累积，但始终保持在10−4量级以下。这表明，在星座卫星解体事件的应急处理中，对解体事件的及时感知对于其他星座卫星的风险规避至关重要。
4 结 语
本文以低轨星座中某卫星在轨解体事件为背景，首先构建了基于边值问题的解体碎片云短期碰撞风险评估算法，并利用蒙特卡罗方法验证了该算法的正确性；随后，针对给定的星座卫星解体场景开展了解体事件的短期危害评估。结果表明，星座卫星解体事件的短期危害具有以下四个特征：
（1）对于共面卫星，碎片云沿迹向扩散是造成碰撞风险的主要因素。由于碎片云持续扩散，卫星每次穿越碎片云时，受到的最大瞬时撞击风险逐渐降低。
（2）对于与解体卫星存在高速接近的异面卫星，与碎片云的高速交会是造成碰撞风险的主要因素，且每次交会所造成的最大撞击风险也逐渐降低。
（3）对于与解体卫星不存在高速接近的异面卫星，碎片云扩散仍会引发碰撞风险，虽然其最大撞击风险随时间推移而逐渐升高，但整体风险等级相对较低。
（4）解体碎片云的撞击风险高度集中在与其高速交会的少量异面卫星上，且在解体后迅速达到最大瞬时撞击风险，因此，对解体事件的及时感知与监测是有效规避风险的关键。
本文算法和相关结论，可为星座卫星解体事件的早期危害快速评估提供有力支撑。