“空间碎片”(Space debris)专指人类在太空活动中产生的废弃物及其衍生物，包括各类人造物体及其碎片和零部件。根据美国空间目标监测网(U.S. Space Surveillance Network, SSN)最新数据表示，目前太空中有约18 000个已编号可监控空间碎片目标^[1]。据估计，大约有370 000个空间碎片以高达22 000 m/h的速度环绕地球轨道飞行。因此，空间中正常运行或将要发射的航天器正面临着空间碎片的干扰与碰撞的巨大威胁。最著名的一次空间目标(Space objects)碰撞发生在2009年，铱星33 (Iridium 33)和宇宙2 251 (Cosmos 2251)通信卫星在西伯利亚北部相撞^[2]，这极大地引起了人们的关注，同时对于该次碰撞的研究也获得了一系列重要成果。因此，需要通过预测碎片碰撞的概率来评估碰撞风险，并通过对风险的评估来指导空间目标的机动决策。

空间目标碰撞概率的计算一般基于如下的假设：两目标在相遇期间在惯性系中的位置、速度矢量是已知的；两目标的相对运动是线性的；两个已知目标均等效于已知半径的球体。当两目标中心的距离小于等效半径之和时，认为碰撞事件发生^[3]。基于概率的解析方法在近年来也得到了长足的发展。其基本思想是将不确定性量，即位置误差看作是服从三维高斯分布的随机变量，碰撞概率考虑两目标的状态矢量、误差协方差矩阵以及等效球体的形状，将碰撞概率的计算转化为概率密度函数(PDF)在特定区域内的积分。Akella等将三维的概率密度函数投影到二维的碰撞平面上^[4]，白显宗等进一步将位置误差矩阵投影到相遇平面，并引入了一种压缩空间和坐标旋转的方法，得到PDF积分的一重积分的解析表达式^[5]。上述概率的计算基于位置误差矩阵的分析，同时也忽略了速度的不确定性。在假设条件不满足的情况下很难得到良好的结果。另外一类方法是将空间目标运行的不确定性归纳到模型传播中，针对输出进行多次采样来确定模型的性质。通常的方法包括蒙特卡洛(Monte-Carlo, MC)模拟。该方法将碰撞概率视作关于不确定性变量的函数，计算精确度较高，但需要大量的数值模拟，对目标在空间中的传播进行重复计算，对计算量要求较大。

代理模型充分考虑了不确定性因素在传播过程中对输出结果(即位置、速度矢量)的影响，同时由于其计算复杂度低、数学表达明晰等特点在仿真试验中得到了广泛的应用^[6-7]。给定试验设计及观测数据，代理模型通过回归或者插值等方法拟合观测数据，可给出相应的预测均值、方差等统计量。本文主要考虑两类代理模型，即混沌多项式逼近方法(Polynomial chaos expansion, PCE)及高斯过程回归方法(Gaussian processes regression, GPR)^[8]。

PCE的基本思想是将精确解在随机参数空间中进行(有限项)多项式展开，通过求解系数得到精确解关于随机参数多项式的线性组合，同时获得关于精确解的所有信息^[9-10]。与PCE模型指定基函数不同，GPR模型是一类典型的非线性模型，其建立在响应服从某一个确定的高斯过程的假设之上。一般的，可通过求解关于观测数据的似然函数得到关于模型相关超参数的极大似然估计，从而得到关于响应的估计。

1 基于蒙特卡洛方法的碰撞概率估计

设概率空间为(Ω, F, P)，其中Ω为事件空间，F为其σ代数，P是定义在F上的概率测度。设$ \mathit{\boldsymbol{x}} = {({x_1}, {x_2}, \ldots , {x_d})^{\rm{T}}}$为定义在有界区域$ \mathit{\boldsymbol{D}} \subset {{\bf{R}}^d}$上的d维独立随机变量，代表输入的不确定性。对于原轨道模型M，设t=[0, t_f]为时间变量，u₀为初始状态，则对于原模型满足如下条件

$ M\left( {t,\mathit{\boldsymbol{x}};{u_0}} \right) = f\left( {t,\mathit{\boldsymbol{x}}} \right),\;\;\;\;\;\left( {t,\mathit{\boldsymbol{x}}} \right) \in \left[ {0,{t_f}} \right] \times \mathit{\boldsymbol{D}} $

(1)

通过随机输入x在上述模型中的传播，得到关于输出状态$\mathit{\boldsymbol{u}}\left( {t, {\rm{ }}\mathit{\boldsymbol{x}}} \right) = {\left[ {\mathit{\boldsymbol{r}}\left( {t, {\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}} \right), {\rm{ }}\mathit{\boldsymbol{v}}\left( {t, {\rm{ }}\mathit{\boldsymbol{x}}} \right)} \right]^{\rm{T}}} $。给定两个空间目标(l = 1, 2)，则在t_k时刻，两目标之间的距离为

$ s\left( {{t_k},{x_i},{{x'}_i}} \right) = \left\| {{\mathit{\boldsymbol{r}}_1}\left( {{t_k},{x_i}} \right) - {\mathit{\boldsymbol{r}}_2}\left( {{t_k},{{x'}_i}} \right)} \right\| $

(2)

设在t_k时刻一共进行了M次模拟，则碰撞概率P_c为所有模拟中发生碰撞事件的频率，即

$ {P_{\rm{c}}} = \frac{{{\rm{\# }}\left\{ {s\left( {{t_k},{x_i},{{x'}_i}} \right) \le R} \right\}}}{M} $

(3)

式中R为两空间目标等效半径之和，i =1, 2, …, M。上述模型直接对原模型进行计算，一般情况下，为了得到收敛的概率，需要的进行大量的仿真实验。

2 基于代理模型的碰撞概率估计 2.1 PCE代理模型

在PCE中，假设u(t, x)是有限阶多项式序列的组合^[11]，形式如下

$ \mathit{\boldsymbol{u}}\left( {t,\mathit{\boldsymbol{x}}} \right) \approx \mathit{\boldsymbol{\bar u}}\left( {t,\mathit{\boldsymbol{x}}} \right) = \sum\limits_{\left| i \right| < N} {{c_i}\left( t \right){\phi _i}\left( \mathit{\boldsymbol{x}} \right)} $

(4)

式中$ \{ {\phi _i}\left( \mathit{\boldsymbol{x}} \right)\} $关于d维变量N阶多项式空间的PCE基函数集。设x的联合概率测度为ρ(x)，则有

$ \int_D {{\phi _m}\left( \mathit{\boldsymbol{x}} \right){\phi _n}\left( \mathit{\boldsymbol{x}} \right)\rho \left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = {\delta _{mn}} = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 0\\ 1 \end{array}&\begin{array}{l} m \ne n\\ m = n \end{array} \end{array}} \right. $

(5)

由方程(5)可以看出，基函数$\{ {\phi _i}\left( \mathit{\boldsymbol{x}} \right)\} $由ρ(x)所决定。通过传统的逼近理论可以得到

$ {\bf{E}}\left[ {{{\left( {\mathit{\boldsymbol{u}}\left( {t, \cdot } \right) - \mathit{\boldsymbol{\bar u}}\left( {t, \cdot } \right)} \right)}^2}} \right] \to 0,\;\;\;\;\;N \to \infty $

(6)

式(6)表明该逼近可以达到任意精度^[12]，保证了逼近的准确性。因此问题的关键在于系数c_i(t)的计算。本文采用非介入式投影方法对系数进行求解，具体如下

$ \begin{array}{l} {c_i}\left( t \right) = \mathit{\boldsymbol{E}}\left( {\mathit{\boldsymbol{u}}\left( {t,\mathit{\boldsymbol{x}}} \right){\phi _i}\left( \mathit{\boldsymbol{x}} \right)} \right) = \int {\mathit{\boldsymbol{u}}\left( {t,\mathit{\boldsymbol{x}}} \right){\phi _i}\left( \mathit{\boldsymbol{x}} \right)p\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \\ \;\;\;\;\;\;\;\;\;\sum\limits_{j = 1}^7 {{w_j}\mathit{\boldsymbol{u}}\left( {t,{X_j}} \right){\phi _k}\left( {{X_j}} \right)} \end{array} $

(7)

式中：第一个等式充分利用了基函数的正交性；第三个等式采用了高斯积分的方法；X_j与w_j分别为高斯积分点及权重。与传统最小二乘方法对比，该方法能够避免对矩阵求逆，计算的复杂度由样本点(即高斯积分点)数量决定。

2.2 GPR代理模型

在GPR中，假设u(t, x)服从某一个确定的高斯过程基础之上，即：$ \mathit{\boldsymbol{u}}\left( {t, \mathit{\boldsymbol{x}}} \right) \sim {\rm{GPR}}(0, k\left( {\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{x}}} \right){\rm{ }} + {\rm{ }}{\sigma ^2})$，其中k(x, x)为核函数。具体来说，给定试验设计(X, Y)，有如下似然函数

$ \mathit{\boldsymbol{u}}\left( {t,\mathit{\boldsymbol{X}}} \right)\left| {\mathit{\boldsymbol{X}} \sim N\left( {{\bf{0}},k\left( {\mathit{\boldsymbol{X}},\mathit{\boldsymbol{X}}} \right) + {\sigma ^2}I} \right)} \right. $

(8)

记K=k(X, X)∈R^N×N，由高斯过程假设以及式(8)，可得预测均值和方差分别为

$ \mathit{\boldsymbol{\bar u}}\left( {t,x} \right) = K_ * ^{\rm{T}}{\left[ {K + {\sigma ^2}I} \right]^{ - 1}}Y $

(9)

$ {\rm{Var}}\left[ {\mathit{\boldsymbol{\bar u}}\left( {t,\mathit{\boldsymbol{x}}} \right)} \right] = {\mathit{\boldsymbol{K}}_{ * * }} - \mathit{\boldsymbol{K}}_ * ^{\rm{T}}\left[ {K + {\sigma ^2}I} \right]{\mathit{\boldsymbol{K}}_ * } $

(10)

其中$ {\mathit{\boldsymbol{K}}_*} = k\left( {\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{x}}} \right) \in {\mathit{\boldsymbol{R}}^{N \times 1}}, {\mathit{\boldsymbol{K}}_{**}} = K\left( {\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{x}}} \right) \in \mathit{\boldsymbol{R}}$。实际上，高斯过程的预测能力由核函数所决定。本文共使用了线性核函数、高斯核函数以及Matern3-2核函数来进行计算。一般来说，针对似然函数(8)进行优化可以得到关于核函数参数的相关估计。

2.3 碰撞概率估计及精度分析

本文采用SGP4模型对空间目标进行仿真，该具体模型比较复杂，不在本文中阐述，可详细参考文献[1 3]。根据文献[14]，SGP4模型初始不确定性(Initial uncertainties)来源于重力加速度(Geopotential acceleration)、大气阻力(Atmospheric drag)、太阳光压(Solar radiation pressure)和三体引力(Third body gravity)。本文针对SGP4模型的特点，假设不确定性来源于地球半径误差、地球引力常数(影响重力加速度)以及BSTAR阻力系数，合计为x。设该三项不确定性量服从高斯分布，标准差分别为20 km，0.4 km³/s²以及10^-5。通过x在SGP4模型中的传播可以得到空间目标的状态$ \mathit{\boldsymbol{u}}\left( {t{\rm{ }}, {\rm{ }}\mathit{\boldsymbol{x}}} \right) = {\left[ {\mathit{\boldsymbol{r}}\left( {t, {\rm{ }}\mathit{\boldsymbol{x}}} \right), {\rm{ }}\mathit{\boldsymbol{v}}\left( {{\rm{ }}t, {\rm{ }}\mathit{\boldsymbol{x}}} \right)} \right]^{\rm{T}}}$，从而给定试验设计X可以得到相应的观测Y。本文经过SGP4模型计算得到的位置及速度矢量都是在地球固定坐标系(ECF坐标系)中的。注意到输入变量的随机性认为是服从高斯分布的，但经过模型的传播，最后的输出可能呈现任意的分布形式。

本文中碰撞概率的计算与MC方法类似，但是使用代理模型替代原模型，这样在计算量上能够得到很大的节约。基于代理模型的碰撞概率估计主要分为以下几步：

(1) 利用SGP4模型生成N个样本，并针对其得到PCE代理模型或GPR代理模型；

(2) 对代理模型进行M次模拟，即利用代理模型计算空间目标的状态；

(3) 利用(3)式计算碰撞概率。

显然，碰撞概率的精度由代理模型的精度来保证。本文针对PCE代理模型分别计算其留一交叉验证(Leave-one -out cross validation, LOOCV)误差和预测的标准差(或变异系数)。前者反映了代理模型的可靠性，后者反映了代理模型预测的离散程度。同样的，对于GPR模型本文将LOOCV误差作为其精度指标。

3 仿真校验

为验证上述方法的正确性，以发生在2005年07219号碎片(美国雷神火箭推进器)与26207号碎片(中国长征四号火箭第三级解体后的碎片)的碰撞事件^[15]为例，所选取两个空间目标的TLE(来源：https://www.space-track.org/)如表 1所示。

根据表 1，可得初始历元时刻为t₀= 2005.1.16. 47659.26。以07219号碎片为例，令时间步长为1 s，分别利用SGP4模型，PCE代理模型及GPR模型对上述数据进行1 d的模拟，预测其在${t_f} = {\rm{ }}{t_0} + 1{\rm{ d}} $时刻的位置(X, Y, Z)及速度(V_X, V_Y, V_Z)。对于GPR代理模型来说，本文共采用了线性核函数、高斯核函数以及Matern3-2核函数来对SGP4模型进行逼近。各模型的LOO误差及标准差(变异系数)如表 2所示。从表 2中可以看出，PCE模型的LOOCV误差远远比GPR模型的LOOCV误差要小，因此对于SGP4模型构建的PCE模型可靠程度更高。此外，PCE模型的变异系数均在1%以下，说明对于预测均值来说，其预测的方差较小，从另一个方面也能说明其精度较高。对于GPR模型来说，可以看出高斯核函数较其他两类函数的LOOCV误差较小，因此其对于SGP4模型的适应程度比线性核函数及Matern3-2函数要高。

此外，在一台主频3.2 GHz，内存32 GB的电脑上，针对t_f时刻本文分别用SGP4模型以及4类代理模型进行了10 000次的仿真，所需的计算时间如表 3所示。

由表 3可以看出，在计算时间上看，代理模型比原SGP4模型具有巨大的优势。PCE模型对比于GPR模型来说其计算复杂度要更低一个量级，不同核函数对GPR模型的计算时间影响并不大。针对这10 000次仿真，以X方向为例，对比其真实分布及通过代理模型逼近得到的经验累计分布图 1所示。

注意图 1中显示F(X) < 0. 1区域的局部放大图，可以看出PCE模型以及高斯核函数能够对SGP4模型进行较为精确的逼近，同时说明位置、速度等状态上的不确定性可以通过输入不确定性在模型中的传播表示出来。

下面针对碰撞概率进行计算。经查询，两个碎片的碰撞发生在t_c=2 005.1.17 :0214 GMT。通过SGP4模型仿真可得在t_c±50 min的空间轨迹如图 2所示，其相对距离如图 3所示。

利用代理模型，通过计算两空间碎片之间的最近距离，得到相应的碰撞时间点，记为t₀+t′。计算所得相关信息如表 4所示。

对t_c±100 s的时间间隔内进行仿真，针对每个代理模型，对每个时间节点进行10 000次仿真，设距离阈值为1 000 m，计算每个时间点的碰撞概率P_c，结果如图 4所示。

由于利用4种代理模型通过MC方法计算碰撞概率得到的结果基本上保持一致，因此在图 4中只用一条曲线表示碰撞概率。从图中可以看出，在大概15 s的时间范围内，碰撞概率急剧增大，并且会很快增长到0.1左右的阈值；在两个空间目标分离后，碰撞概率又急速下降，与真实情况相吻合。

4 结论

上述实验结果表明，对于具有不确定性输入变量的复杂模型，通过构造相应的代理模型，可以较大的降低计算量。不同的代理模型原模型的逼近效果不同，因此需要对各模型进行筛选。本文针对SGP4模型，PCE代理模型及基于高斯核函数的GPR模型具有非常高的精度，同时PCE代理模型本身的误差也能控制在较小的范围之内，具有较强的稳健性。与GPR模型相比较，PCE代理模型充分利用了输入变量的先验分布，并且引入了正交多项式基，在计算方面其效率要更高；针对空间中可能存在的碰撞事件，各代理模型均能够较为精确的对碰撞事件进行预测，PCE模型的预测效果最好。同时，代理模型计算量小的特性可保证其能应对需快速响应的紧急情况。