0 引 言
近年来，国内外学者和机构针对碎片云的产生开展了大量研究，并建立了多种解体模型，如美国国家航空航天局(National Aeronautics and Space Administration，NASA)的标准解体模型(Standard Breakup Model, SBM) [2]、欧空局的Battelle模型[4]，以及中国空气动力研究与发展中心的侵彻与剥落解体模型[5]等。Cimmino等[6]和Diserens[7]分别从参数校准与适应性角度，对NASA SBM进行了改进分析。在碎片云碰撞概率计算方面，以蒙特卡洛仿真[8]为代表的传统方法，通过对初始状态大量采样并逐条传播轨道来统计碰撞事件，其精度随采样数量增加而收敛；然而，对于概率低至10−8量级的短期碰撞事件，所需采样数量极大，计算代价难以承受。为提高计算效率，基于连续介质假设的密度演化方法被提出：该方法在轨道根数空间中传播空间密度函数[9]，计算成本与碎片数量无关。Frey和Colombo[10]进一步将这一思想与NASA SBM衔接，通过变量代换技术直接从NASA SBM导出碎片初始密度的解析表达式。不过，密度方法通常需引入动力学简化假设，在碎片云产生后数小时内的强非均匀、快时变阶段，其精度存在局限。
在BVP框架下，存在若干具有代表性的实现方案。Parigini等[11]采用微分代数技术(Differential Algebra, DA)，在兰伯特解和状态转移矩阵上构造高阶泰勒展开，并结合自动域分割(Automatic Domain Splitting, ADS)在被积函数强非线性区域递归细分积分域。该方法在开普勒动力学和J2摄动条件下均得到验证，但其分割过程需在子域内反复重构局部泰勒展开，未通过误差检验的展开难以直接参与最终积分，因而在高非线性时间窗口内可能带来较高的计算需求与耗时开销。文献[12, 15]则从碎片初始分布的连续概率密度函数出发，通过多圈兰伯特解联合传播碎片的位置-速度联合分布，推导碰撞概率的解析表达式。Trombetta等[3]在德国空间运行中心的工程实践中，以均匀时间步进配合梯形积分验证了BVP框架的工程可行性，但时间离散策略在精度与效率方面仍有优化空间。
本文在保留Parigini等[11]提出的边值问题框架及物理模型的基础上，引入洛巴托伪谱离散方法(Lobatto Pseudospectral Method, LPM)以替代原有的微分代数法与自适应网格细化策略，提出一种由LPM驱动的短期碎片云碰撞概率计算方法。在时间窗口内，通过LGL节点直接采样兰伯特解并构造拉格朗日插值，利用洛巴托求积权重完成概率积分；同时引入基于导数极值检测、曲率分析与局部误差评估的细化策略，构造局部LGL离散，以最少的子区间数量实现预设精度要求。
1 问题描述
设在t0时刻空间位置r1处产生新的碎片云，通过NASA SBM可得该碎片云包含Nd个碎片，其可能对位于位置r2、速度为v2的重要卫星构成碰撞风险。如图1 所示，为第i个碎片在r2处与卫星发生碰撞的示意，其中碎片云的速度矢量定义为vi = v1 + Δvi（v1为碎片生成前的速度，Δvi为碎片生成时的速度增量）。本文聚焦于碎片云产生后短期内与卫星碰撞概率的快速计算问题。
根据Jenkin[8]的观点，第i个碎片的碰撞概率可基于合理的碎片分布假设，在速度空间中进行分析。该计算需沿卫星截面积路径，对速度空间内的概率密度函数进行积分，即：
式中，ξ为卫星在速度空间中的运动路径，fi(Δv)为速度空间内的概率密度函数，dV为卫星截面积在dt时间内扫过的微分体积，其表达式可写为：
将式(2)代入式(1)，可得：
得到每一个碎片的碰撞概率Pi后，整个碎片云的碰撞概率PC即可写为：
由于直接对碎片云开展碰撞概率计算存在计算复杂、计算量大等问题，因此首先筛选碎片云与卫星可能发生交会的时间和位置窗口，其次采用LPM离散方法对各窗口进行初始离散并求解碰撞概率，最终得到整个计算时间内的累计碰撞概率。基于此，式(3)可改写为：
式中，k表示初步筛选出的时间窗口数量。
2 碎片云碰撞概率计算
2.1 时空初筛
为避免在整个时间域上盲目计算造成的算力浪费，首先通过初步交会评估确定可能发生碰撞的位置区间与时间区间。由于碰撞仅可能发生在卫星轨道上，结合碎片云生成点与卫星位置的几何关系，利用Battin[21]自由时间最小脉冲方法求解碎片到达卫星位置所需的速度增量Δv，并依据速度增量阈值初步筛选位置窗口；随后通过求解兰伯特问题确定所需时间窗口。本文设定1 200 m/s为速度增量上限，超过该阈值的碎片速度增量出现概率极低，对碰撞概率的影响可忽略不计。初步评估问题可转化为一个可判定问题，即对于卫星轨道上的任意位置r2，是否存在满足||Δv|| ≤ Δvmax的单脉冲转移，使得碎片能从生成点r1到达该位置。
在转移平面内，碎片速度可分解为弦向速度vc与径向速度vρ。根据Battin[21]自由时间最小脉冲理论，所有能从r1到达r2的平面内单脉冲速度解满足：
式中，θ为矢量r1与r2之间的转移角，c = || r2−r1||为弦长，r1 = || r1||、r2 = || r2||分别为碎片生成点和卫星位置的地心距，μ为中心天体引力常数。式(6)表明，在不限定转移时间的情况下，满足两点位置约束的径向速度与弦向速度组合位于vρ–vc平面内的一条双曲线上。
由图2 几何关系可知，总速度增量为Δv = vn + Δvπ，且vn与Δvπ正交，有||Δv||2 = ||vn||2 + ||Δvπ||2，故总速度增量约束可写为：
式中，vn = (${\boldsymbol{V}}^{\mathrm{T}}_1 $ un)un为交会前速度的法向分量，vπ = v1−vn为交会前速度在转移平面内的投影。若vn超过阈值上限，则无论如何调整平面内速度，都无法构造满足约束的单脉冲交会，因此该位置将被剔除；仅当vn ≤ Δvmax时，才有必要继续检查平面内可达性。将式(6)改写为vρ关于vc的函数，再代入速度增量达到边界时的式(7)~(8)，即可将该问题转化为关于弦向速度vc的四次代数方程$v_c^4 + a_3v_c^3+ a_2v_c^2 $ +a1vc+a0 = 0[21]。对于上述一元四次方程，本文采用伴随矩阵特征值法进行数值求解，具体方法是将该四次多项式构造为伴随矩阵：
则矩阵C的特征值即为该四次方程的全部根，求得根后，仅保留虚部小于给定容差的有限实根。若有效实根集为空，则认为在给定速度增量上限下，该卫星的真近点角位置不可达；若有效实根集非空，则说明存在至少一组平面内速度分量可使碎片由r1到达r2，此时取其实根范围[vc,min, vc,max]作为弦向速度的可行区间，并保留对应的卫星真近点角ftarget，将其归入可达域集合F。
得到ftarget∈F后，还需将几何位置可达域转化为时间交会区间。对此，本文进一步在保留的几何位置上执行多圈时间匹配：对于集合F中的每一个真近点角，分别写出碎片和卫星从生成时刻到该位置所需的飞行时间——两者均为各自平近点角差值与转移圈数的函数，令二者相等即可构成关于Δv和t的联合求解问题，进而完成交会时间区间的筛选。这一步可利用Izzo的多圈兰伯特求解器[23]实现。需要强调的是，这些时间窗口可能存在重叠，但只要对应转移圈数或分支不同，即代表不同的转移路径。
2.2 时间伪谱离散方法
对于任意时间窗口Tk中的时间t，碎片从生成位置r1转移到卫星位置r2(t)的多圈兰伯特问题可写为：
式中，μ为中心天体引力常数。由此得到每个窗口上的速度增量为：
若在某一时间区间Tk = [tk,min, tk,max]内存在一族速度增量vk(t)可使碎片与卫星交会，则下一步需据此计算第i个碎片在该时间区间内的碰撞概率积分。
时间窗口内碰撞概率积分的难点在于，Δvk(t)并非显式函数，而是通过离散时间窗口、在离散时刻反复求解多圈兰伯特问题得到。若直接在时间窗口Tk上采用均匀加密采样，一方面会显著增加兰伯特问题的求解次数，另一方面又难以兼顾曲线平缓区域与快速变化区域的近似效率，因此引入伪谱离散方法以平衡计算效率与精度。
由于插值节点τ仅在区间[−1, 1]内部生成，为此需先将窗口时间归一化至标准区间。对于第k个时间窗口Tk = [t0, tf]，引入时间变量τ，定义为：
在标准区间[−1, 1]上取n + 1个LGL积分点，后续关于插值、微分和求积的离散关系均先在标准空间中建立，再映射到物理时间t。考虑n次勒让德正交多项式序列Ln(τ)，则洛巴托节点生成多项式定义为：
则n + 1个LGL节点由区间两端点与内部导数零点组成，即：
在节点集{τj}n j=0上，引入拉格朗日基函数：
因而时间窗口内任意足够光滑的函数都可以在节点上通过全局多项式进行拟合。根据式(11)可知，节点速度增量可记为：
则其n次拉格朗日插值多项式表达式为：
式中，各节点上的Δvk,j由多圈兰伯特问题求解得到，节点之外的连续曲线则由插值多项式表示。定义微分矩阵D：
矩阵D是(n+1)×(n+1)维矩阵，因此速度增量对时间τ的导数为：
从而推导得到：
矩阵D的每一行对应某个配置点上的导数信息，因此节点处速度增量的导数可由节点函数值的线性组合得到。至此，全局LGL近似已构建完成，但该近似无法直接用于最终积分——仅依赖单一全局多项式时，在局部曲率较大或导数变化较快的区域，近似精度可能会受到损失。为此，可先通过全局近似刻画窗口内速度增量曲线的整体变化趋势，再根据速度增量轨迹的导数极值、曲率峰值及误差等信息，筛选出需要进行局部离散的区域。
设窗口上的全局伪谱重构为$ \Delta \boldsymbol{v}_{k}^{p}\left(t\right) $ ，并定义：
式中，$ \tilde{\kappa } $ 用于刻画速度增量曲线的局部弯曲强度。极值区域集为：
高曲率区域集为：
式中，Qq(·)为分位阈值算子，用于从窗口内曲率序列中提取高曲率阈值。两个集合与时间区间端点共同形成初始断点集：
在断点处重新生成子区间，设第k个窗口经局部细化后被分为若干个子区间：
针对某一个子区间，引入局部归一化映射并再次进行伪谱离散构造拉格朗日插值多项式和微分矩阵：
对第s个子区间，局部相对误差为：
若某个局部子区间的相对误差超过阈值，则在误差峰值附近进行二分处理，并通过最小区间约束避免产生过窄子区间，从而形成以局部误差为驱动的递归时间细化策略。在局部子区间Tk,s上，若记任意被积函数为gk,s(t)，则由式(26)可得：
利用LGL求积权重ωj，可将式(28)写为：
由此，通过LPM离散并重构得到Δvk及微分矩阵D等信息，用于后续碰撞概率数值积分计算。
2.3 卫星碰撞截面积计算
确定时间窗口内各节点兰伯特解后，需将卫星映射到速度空间中，以计算碰撞截面积。构造初始时刻碎片状态为：
在二体动力学模型下传播至交会时刻，得到状态转移矩阵：
令δr0 = 0，通过数学推导可得初始速度扰动到末端位置扰动的映射关系：
设卫星为半径为R的球体，则交会时刻卫星末端位置的空间域可写为Br(R) = {δrf ∈ R3: ||δrf || ≤ R}。为在速度空间中计算碰撞截面积，将该球域逆映射到初始速度空间，令J = ${\boldsymbol{\varPhi}}^{-1}_{\boldsymbol{r}\boldsymbol{v}} $ ，则该球域在线性映射下形成椭球。速度增量在t时刻的切向单位向量为：
碰撞概率积分仅与该椭球在法平面Πt={δv0∈ R 3: ${\boldsymbol{u}}^{\mathrm{T}}_t $ δv0 = 0}上的中心截面积相关。通过对几何关系的进一步推导可知，速度空间中以ut为法向的椭球截面，对应末端位置空间中以nr为法向的球体中心截面：
经几何推导，该截面的面积由线性映射的雅可比矩阵行列式与切向投影关系给出，即可得到时间区间上卫星的有效碰撞截面积
2.4 碰撞概率计算
在NASA的SBM模型中，速度增量与特定的面质比A/M相关，其取值来源于速度增量概率密度函数fi(Δvk)，且该速度增量符合如下正态分布：
式中，均值μ是关于χi=lg(A/M)的函数，标准差σ=0.4。假设速度增量各向同性，则速度增量的概率密度函数可写为：
将式(20)、式(35)和式(37)代入式(5)，可得碰撞概率计算表达式为：
式中，Ik(τ, χ) = fk(τ, χ)Ak(τ)||dΔvk/dτ||。对于NASA SBM生成的Nd个碎片，将具有相同面质比的碎片归为同一类别；针对面质比为χi的类别，第k个时间窗口的碰撞概率可表示为：
这一步可通过式(29)计算，最终利用式(4)计算得到整个碎片云的碰撞概率。
3 仿真分析
考虑新生成碎片云与一颗正常工作卫星之间的碰撞风险，相关轨道参数如表1 所示[11]，该卫星半径R为1 m。通过NASA SBM可产生约10 000个特征长度大于1 cm的碎片，且碎片面质比范围为10−3至102。
表 1 轨道参数设置Table 1 Orbital parameters setting |
| 物体 | A/km | e | i/(°) | Ω/(°) | ω/(°) | f/(°) |
| 生成碎片云的物体卫星 | 7 223 | 0.003 8 | 98.73 | 320.97 | 354.36 | 344.92 |
| 工作卫星 | 7 038 | 0.001 4 | 98.03 | 0 | 0 | 342.79 |
本文列出LPM初始离散与局部细化过程中采用的主要参数，如表2 所示。
表 2 LPM离散与局部细化参数设置Table 2 Parameters of LPM discretization and local refinement |
| 参数 | 数值 |
| 全局LGL阶数Ng | 25 |
| 局部LGL阶数Nl | 12 |
| 曲率分位数q | 0.88 |
| 误差阈值 | 5×10−3 |
| 最大细化层数 | 3 |
| 初始断点最小间距 | 10 s |
| 最小子区间长度 | 8 s |
对于任意时间区间Tk，均需满足法向速度vn≤vmax，其与真近点角的关系如图3 所示。其中，f表示卫星轨道上的真近点角，即前文可达域判定中的ftarget。
通过对空间位置的初步筛选，可得到碎片云与卫星的交会可能发生在真近点角f$\in $ [79°, 107°]和f$\in $ [258°, 287°]。在明确空间交会可行域后，结合碎片与卫星的轨道周期开展时间维度的匹配筛选，具体如图3 和图4 所示，从而确定碎片与卫星可能发生碰撞的时间与空间区域。
本文以第一个时间窗口T1的计算为例。若碎片要在该短期窗口内与卫星发生碰撞，其所需的速度增量Δv随交会时间t呈现出强烈的非线性变化。这种剧烈的非线性特性导致直接通过数值积分求解碰撞概率的过程极其耗时。为此，本文采用LPM伪谱离散方法对速度增量进行离散重构：在该时间窗口内进行高密度时间采样，并在每个采样时刻直接求解多圈兰伯特问题以获取速度增量的真值，结果如图5 所示。可以看到，全局单次LGL拟合已能大致捕捉Δv曲线的基本形态，但在曲线中部非线性较强的区域出现了明显偏离；经局部细化后，分段重构曲线能够准确跟踪窗口内速度增量的局部变化特征。该窗口全局单次拟合的最大相对误差约为1.57×10−2，经LPM分段处理后，误差降至2.18×10−4，精度提升约两个数量级，且仅需划分2个子区间即可实现。
此外，还对比了单段全局LPM离散与多次局部LPM离散两种方式下的碰撞概率计算结果。在碎片云生成初期，单段全局离散方法存在较大偏差，导致PC计算结果明显偏高，无法精确拟合12 h内的碰撞概率变化趋势；而采用局部离散方法得到的碰撞概率在整个仿真时间段内均表现出更高的计算精度，其曲线与微分代数方法的结果更为贴合。
为进一步验证本方法对离散参数的稳定性，针对局部伪谱阶数开展了敏感性分析（如表3 所示），并以高阶局部LGL方法(Nl = 50）对应的碰撞概率作为参考值。
表 3 局部LGL阶数敏感性分析Table 3 Sensitivity analysis of the local LGL order |
| 局部LGL 阶数Nl | 子区间 数量 | 计算 时间/s | 碰撞概率PC | 相对 误差% |
| 6 | 71 | 2.094 | 9.916 5×10−8 | 0.017 4 |
| 10 | 50 | 1.818 | 9.936 7×10−8 | 0.186 0 |
| 12 | 50 | 1.836 | 9.915 2×10−8 | 0.030 6 |
| 25 | 50 | 2.132 | 9.918 3×10−8 | 0.000 2 |
| 50 | 50 | 2.582 | 9.918 2×10−8 |
结果表明，提高局部伪谱阶数可减少所需的时间分段数量，但对最终累计碰撞概率的影响较小，碰撞概率相对变化小于0.2%，计算时间未见显著变化。
4 结 语
本文针对碎片云生成后短期内碎片云与卫星的碰撞风险快速评估问题，在Parigini等[11]提出的边值问题短期碎片云碰撞概率分析框架基础上，提出一种基于LPM的时间离散与窗口积分概率计算方法。在不改变原有物理建模框架的前提下，将基于LGL节点的拉格朗日插值、伪谱微分矩阵与洛巴托求积的时间离散体系引入窗口积分，以替代原有的DA泰勒展开与ADS策略。仿真结果显示，采用该方法计算的碰撞概率误差约为1.7%，计算时长由201.124 s缩短至1.836 s。
总体而言，本文工作表明，在短期碎片云碰撞概率问题中，将LPM引入时间窗口的时间离散，是一种兼具物理一致性与数值稳定性的实现方式。本文方法主要适用于碎片云生成后短期内、碎片轨道尚未充分扩散且时间窗口可由多圈兰伯特问题描述的场景，当前算例主要在二体动力学框架下验证，复杂摄动、非球形卫星姿态变化、碎片速度方向非各向同性以及实际观测数据等因素仍有待进一步研究。