半解析卫星轨道不确定性快速估计方法

于彦君，岳程斐，李化义，陈雪芹，刘培

（1.哈尔滨工业大学卫星技术研究所，哈尔滨 150001；2.哈尔滨工业大学（深圳）空间科学与应用技术研究院，深圳 518055；3.上海卫星工程研究所，上海 710119）

0 引言

轨道不确定性传播是指在给定轨道初始状态及其相关不确定性的基础上，预测未来的轨道状态及其统计特性。轨道不确定性传播在空间态势感知与决策等相关任务中起着重要作用，例如轨迹跟踪、碰撞概率分析、传感器资源管理和异常检测等［1］。对于航天器而言，其轨道状态矢量的真值是无法确切知道的，当存在测量信息时，轨道状态误差可以通过滤波算法［2］进行估计，该过程通常是收敛的；当不存在测量信息时，可利用滤波算法的预测步骤进行不确定性估计，即进行轨道不确定性传播，该过程通常是发散的。在轨道不确定性传播问题中，初始轨道不确定性在任何坐标系中都不是绝对高斯的，但目前通常采用高斯假设。高斯误差超椭圆体随着传播时间的增加而放大和旋转，并且由于动力学的非线性性质而逐渐变为非高斯的。在笛卡尔坐标系中，非高斯轨道不确定性通常位于“香蕉型”弧线上，该弧线沿航天器的标称轨迹弯曲［3］。轨道不确定性的典型传播过程如图1所示。

图1 轨道不确定性传播示意图［1］Fig.1 Illustration of uncertainty propagation process［1］

轨道不确定性传播方法包含局部线性化方法和蒙特卡洛法（Monte Carlo method，MC）等［4-5］。然而，由于轨道动力学的非线性特质，局部线性化方法过低的精度和蒙特卡洛方法的高计算量极大地阻碍了这2种方法的应用。为了在降低计算量的同时保证合适的精度，针对不同场景，学者们提出了多种适用于非线性轨道不确定性传播的方法，例如高斯和模型（Gaussian mixture model，GMM）［6-7］、无迹变换法（Unscented transformation，UT）［8］、状态转移张量法［9］、微分代数法［10］、多项式展开法［11］、坐标转换法［12-13］以及将上述方法相互结合的混合方法等［14-17］。在这些方法中，UT法［18］解决了线性化方法和蒙特卡洛法不能兼顾计算效率和精度这一难题。其通过非线性积分从初始分布中确定性地选择样本来近似未来时间的概率分布以取代复杂的高阶分析。这一过程中，通过近似概率分布代替近似非线性变换，提高了计算精度；通过较少的样本计算代替MC 法中的大样本计算，提高了计算效率。坐标转换法则通过选择在非线性较弱的坐标系中建立动力学模型，在相同计算量条件下，达到大幅延长不确定性真实度（Uncertainty realism，UR）［19］保持时间的效果。坐标转换法与UT 法相结合可用于轨道不确定性的传播，如文献［20］和［21］。

随着低轨巨型卫星星座的发展，卫星间碰撞概率显著增加，尤其当单颗卫星同时面临与多个卫星碰撞的风险时，快速估计分析轨道不确定性传播结果并评估碰撞风险成为亟待突破的问题。特别是考虑星上有限计算能力条件下，在确保不确定性估计精度的基础上，快速获得某一特定时刻或时段内的轨道状态及其统计特性具有重要的应用价值。

针对这一问题，本文提出了一种基于半解析法的轨道不确定性估计方法。该方法是一种混合方法，具体特点如下：首先选择在轨道要素坐标系中进行不确定性传播，延长UR 的保持时间，跟踪笛卡尔坐标系中非高斯形式的轨道不确定性。其次，采用解析模型结合迭代法实现对任意时刻的轨道要素的快速计算，并将该方法应用于轨道不确定性估计，大幅提升计算效率。此外，采用球形单边采样策略而非目前常用的对称采样策略，进一步降低由采样点数量带来的计算量，并结合球形单边采样策略特性和高斯摄动方程对采样点进行修正，再估计传播后的轨道不确定性，以弥补采用解析模型传播样本点导致的精度降低。最后，将提出的算法与高斯和模型相结合，进一步改善轨道不确定性估计精度。

1 基于半解析法的轨道要素快速计算方法

不确定性传播算法的有效性可以通过使用非线性程度较低的坐标系得到显著提高。例如，相较于笛卡尔坐标系中的卫星轨道动力学，采用轨道要素作为动力学模型能够吸收非线性动力学中最主要的项（即引力中的r-2项），使得轨道要素空间的高斯分布能够更准确地捕捉到在笛卡尔空间中经常观察到的非高斯行为，并能够更长时间地保持UR［22］。本文采用经典轨道要素进行轨道不确定性传播，针对低轨卫星给出了考虑J2 摄动与大气阻力摄动的轨道要素模型和解析解。

1.1 轨道要素长期项与周期项建模

瞬时轨道要素通常可分为长期项、长周期项和短周期项。而上述长期项和周期项可以基于高斯摄动方程［23］（Gauss variational equations，GVE）进行推导。本文采用NTW坐标系下的描述形式。NTW坐标系以航天器为中心，T轴与轨道相切，以速度矢量所指方向为正方向；W轴沿角动量矢量方向；N轴垂直于T轴和W轴，方向由右手定则确定。当0 ＜e＜1，i≠0，i≠180°时，对应GVE为：

令参数向量α≜[a，e，i，Ω，ω，M]T，则在考虑J2摄动和指数模型大气阻力的影响下，长期项的时间导数［24］为：

式中：ρ0为轨道半径r0处的大气密度；K2=(CDS2/m)，K1=(CDS1/m)Q；H为标高；Q=(1 -rPωe；S1和S2分别是垂直于T轴和N轴方向 的横截面积；rP和vP分别代表近地点半径和近地点速度；ωe=7.292 115 855 3 × 10-5rad·s-1为地球角速度，RE为地球半径。

类比文献［26］中的推导方式，将文献［24］中给出的大气阻力长周期项递推方程进一步推导，可得大气阻力摄动长周期项解析解αdl（α）为：

式中：ni和nω分别为轨道倾角和近心点角的平均变化率，近似为：

J2 摄动影响下的一阶长周期项αl（α）、短周期项αs（α）和大气阻力摄动影响下的短周期项αds（α）在参考文献［24-25］中已给出，在此不再赘述。

1.2 任意时刻轨道要素快速计算方法

假设初始时刻为t0，需要获取t时刻的瞬时轨道要素。一般情况下，可通过数值方法，如龙格库塔法对式（1）直接进行积分获取。该方式需要从t0时刻逐步积分到t时刻。然而，在实际情况下，通常只对未来某一时刻或时段内的轨道要素感兴趣。例如，在进行碰撞概率计算时，在整个轨道周期内计算碰撞概率显然是不必要的，通常需要首先进行碰撞时间区间估计，随后仅在该时间区间内进行碰撞分析［27］。

基于轨道要素长期项和周期项模型能够快速计算任意时刻轨道要素。在1.1节中将瞬时轨道要素分解为长期项，长周期项和短周期项。由于长期项是线性项，t时刻的轨道要素长期项，即t时刻平均轨道要素约为：

其中，上标表示对应时刻，αˉt0可通过下式获得：

对于轨道要素周期项的解析解，需要采用瞬时轨道要素作为自变量，而式（8）给出的是t时刻平均轨道要素。因此为了提高精度，t时刻瞬时轨道要素可通过迭代法求解，其步骤为：

1）依据t时刻平均轨道要素计算瞬时轨道要素初值：

2 轨道不确定性半解析快速估计方法

在第1节中给出了快速计算任意时刻轨道要素的方法。因此，在给定t0时刻瞬时轨道要素的均值和协方差后（假设为高斯噪声），即可选取特定的西格玛点，通过改进UT 变换进行轨道不确定性传播。方便起见，将上述算法称为改进半解析无迹变换（Improved semianalytical unscented transformation，ISUT）轨道不确定性估计方法；ISUT 算法与GMM 结合能够进一步提高精度，将该算法称为GMMISUT算法。

2.1 球形单边采样策略

目前采用UT 法进行轨道不确定性传播的文献通常采用传统对称采样策略，即对于一个ns维空间，需要选择2ns+1 个西格玛点。考虑到UT 的计算成本与使用的西格玛点数成正比，对于高维问题，最小化西格玛点的数量可以最小化计算成本。文献［18，28］中给出了球形单边采样策略，它仅需ns+2个西格玛点即可捕获与传统对称采样策略相当的分布信息，达到与截断二阶滤波器相同的预测能力。因此，本文采用球形单边采样策略进行西格玛点的选取。球形单边采样方法通过迭代的方式获得采样点，其步骤可总结如下：

1）选择位于原点（ι=0）和以原点为中心的超球体（ι≠0）上的采样点权重：

2）使用比例无迹变换调整权重：

式中：比例因子σ∈（0，1］。

3）定义各个采样点权重：

式中：β是影响第0采样点协方差权重的参数。

4）当维数大于1 维时，通过式（14）迭代获得采样点：

其中，各个向量序列初值为：

2.2 ISUT方法

根据2.1小节中的球形单边采样策略获得采样点后，选择初始西格玛点集：

采用1.2节中的任意时刻轨道要素快速计算方法传播初始西格玛点集，可获得t时刻的西格玛点集xι（t）。由于该方法计算西格玛点集所采用的动力学基于GVE 的解析解，存在一定的截断误差。为了使得计算结果更接近于GVE 的结果，需要对xι（t）的分布进行调整，即对采样点进行修正。

由式（14）和式（16）可知，位于超球体原点的第0 西格玛点x0(t0)=为不确定性分布的均值。而其余ns+1 个西格玛点位于以原点为中心的超球体上。因此，可通过轨道递推获得GVE 动力学对应的t时刻的第0 西格玛点，计算其与x0（t）的差值Δx，并依据轨道要素间的关系对解析方式获得的西格玛点集xι（t）进行修正，从而减小误差。

本文所采用的动力学模型中平近点角M为快变量，其余5 个轨道要素为慢变量。在进行不确定性传播时，M发散速度较快，其余5 变量发散速度较慢。因此，对于各西格玛点的前5个轨道要素，采用依据差值Δx进行平移的修正策略，即：

对于发散较快的M，设计如下修正策略进行补偿：

基于修正后的西格玛点集，t时刻的瞬时轨道要素均值和误差协方差矩阵为：

则ISUT 轨道不确定性估计算法的流程如图2所示。

图2 ISUT轨道不确定性估计算法框图Fig.2 Diagram of ISUT uncertainty propagation algorithm

最终输出的t时刻瞬时轨道要素均值和误差协方差矩阵将用于构建t时刻轨道不确定性的概率密度函数（Probability density function，PDF），即：

式中：N（x；μU，PU）表示正态分布；x，μU=和PU=Pt分别表示输入变量、均值和误差协方差矩阵。

2.3 GMMISUT方法

根据式（20）可知，ISUT 法传播的不确定性PDF是高斯形式的，即仅能够传播前2阶统计矩（均值与方差）。对于非线性系统，在进行不确定性传播时，不确定性通常会呈现出非线性特性，即包含了高阶统计矩，如峰度和偏度。因此，为了能够进一步提高半解析法的估计精度，可在初始时刻将高斯形式的不确定性近似为GMM。GMM 是对高斯PDF 的直接拓展，其形式为多个高斯PDF的加权总和：

式中：p（x）为待分解的PDF为p（x）的GMM 近似；ηG，μG，PG分别表示第G分量的权重、均值和协方差。NGM表示GMM 组件的数量。为了确保PDF 的非负性，并确保PDF 曲线下的总面积为1，权重必须是非负的，且总和为1，即：

采用GMM 方法对PDF 进行描述在保留高斯PDF 特性的同时扩展了高斯PDF 的适用性。GMM可以近似多种多样的PDF，文献［29］中证明了采用GMM 对PDF 进行近似时，当GMM 组件数量趋向无穷大时均匀收敛。这是一个非常直观的结果，因为随着组件协方差趋近于零，GMM 的每个组件都趋向于脉冲函数。因此，通过减小组件的协方差、增加组件的数量，并适当分布组件的均值，可以轻松地近似大部分PDF，且在进行不确定传播时能够同时传播低阶和高阶统计矩，进一步提升不确定性估计精度。

目前，拆分高斯PDF 有2 种方法，第1 种为在多个方向进行拆分，第2 种为使用单变量拆分库在一个特定的方向进行拆分［30］。采用第1种方法通常需要求解L2 最小化优化问题，会引入额外的计算量。因此，本文采用第2种方法，将多变量高斯PDF在其最大特征值方向进行拆分，步骤如下：

1）将标准正态分布p（x；0，1）拆分为NGM个均质高斯分布组件：

表1 3组件拆分库Table 1 Three-component splitting library

图3 3组件拆分库组件及其总和与标准高斯分布的比较Fig.3 Components of the 3-component splitting libraries and their sum as compared with the standard Gaussian distribution

2）采用谱分解将协方差矩阵PU分解为：

式中：Λ是PU的特征值矩阵，λk是第k大的特征值。选择沿第k个特征向量的方向分解p（x；μU，PU），则式（21）中的参数计算方法为：

式中：υk是PU的第k个特征向量。对于每个GMM 组件，均采用ISUT算法进行传播，即构成GMMISUT 方法。在应用GMMISUT 实现不确定性传播的过程中，应合理选择NGM从而平衡计算效率与精度。

3 数值仿真

3.1 轨道不确定性真实度的似然一致性度量

似然一致性度量L 描述了2 个PDF 之间的重叠程度，对于彼此更一致的PDF来说，L会更大。在将数据点样本与PDF进行比较的情况下，2个PDF之间的似然一致性度量［32］（Likelihood agreement measure，LAM）定义为：

由于需要分析轨道不确定性PDF 与蒙特卡洛法生成的样本的一致性，可假设q(x)由狄拉克模型（Dirac mixture model，DMM）给出：

式中：δ(x-μj)是以μj为中心的Dirac-δ分布，权重γj=K-1。若p（x）为高斯PDF，将式（20）和式（27）代入式（26），可得：

若p（x）为GMM的PDF，将式（21）和式（27）代入式（26），可得：

因此，通过DMM 给定1 组样本点，并与高斯/高斯和模型进行比较，DMM 与高斯/高斯和模型具有相同分布的可能性可以分别通过式（28）和式（29）计算。LAM的值越高，则给定的高斯/高斯和模型更有可能生成DMM。该方法允许将多个高斯/高斯和模型与1 组DMM 进行比较。其中最准确的高斯/高斯和模型具有最高的LAM值。

3.2 初始条件和参数设置

卫星的传感器通常能够提供卫星的瞬时位置、速度和时间信息，且无需处理原始伪距或载波信息。为了更贴近真实情况，在地心惯性坐标系中给出卫星初始测量误差协方差及卫星位置速度均值。假设地心惯性坐标系中卫星的三轴位置和速度误差标准差分别为100 m 和1 m/s。位置均值rECI=［2 505.35 -6 439.95 1 857.00］ km，速度均值vECI=［2.81 -0.96 -6.84］ km/s。选取样本点数目为MC=10 000，通过坐标转换［22］，可得初始瞬时轨道要素的均值和误差协方差如表2所示。

表2 瞬时轨道要素的初始均值与协方差Table 2 Initial mean and covariance of osculating orbital elements

假设初始大气密度ρ0=2.34 × 10-13kg/m3，密度标高H=68.7 km，阻力系数CD=2.2。航天器垂直于切线方向截面积S1=0.2 m2，垂直于次法线方向的截面积S2=0.2 m2，航天器质量m=10 kg。其他参数：σ=1，β=2，W0=0.25。

3.3 仿真结果与分析

为了比较UT、ISUT和GMMISUT法的计算效率，定义UT 法和ISUT 法的计算消耗时间之比εT（t）、UT法和GMMISUT法的计算消耗时间之比εTG（t）为：

式中：τUT（t）为采用UT算法获取t时刻轨道不确定性所需的时间，τISUT（t）为采用ISUT 法获取t时刻轨道不确定性所需的时间，τGMMISUT（t）为采用GMMISUT法获取t时刻轨道不确定性所需的时间。显然，当εT(t) ＞1时，ISUT 法的计算效率高于UT法，当εTG(t) ＞1时，GMMISUT法的计算效率高于UT法。

采用ODE45 方法执行UT 法的积分步骤，并与ISUT 法和GMMISUT 法进行比较，结果如图4 所示。从图4 中可以看出，当t＜1d时，所需获取的轨道不确定性所在时刻t越大，ISUT 算法和GMMISUT 算法的计算效率相较于UT算法越高。当t≥1d后，εT（t）和εTG（t）趋于稳定，ISUT 法的计算效率约为UT 法的7.9倍，GMMISUT法的计算效率约为UT法的2.6倍。本文中对比分析所采用的UT 法采用了单边球形采样策略，采样点数目为8个，均需进行积分；对于ISUT法，采样点数目为8个，需对其中1 个采样点进行积分；对于GMMISUT法，采样点数目为24个，需对其中3 个采样点进行积分。因此可知，上述算法的计算负担大小主要取决于需要进行积分的采样点数目。采用ISUT法和GMMISUT 法能够减小需要数值积分计算的点的数目，因此计算效率更高。

图4 εT（t）和εTG（t）变化情况Fig.4 The variation of εT（t）and εTG（t）

为了分析采用UT、SUT、ISUT 和GMMISUT 法进行轨道不确定性传播的精度，根据3.1 小节中给出的LAM 计算方法，分别计算4种算法的LAM 值。在对比分析UT法，SUT 法和ISUT 法时，以ISUT 法的LAM 值为标准进行归一化，则3 种算法的归一化LAM 值随时间变化情况如图5 所示。可以看出，没有进行修正的SUT 法的归一化LAM 值始终较低，对于轨道不确定性跟踪效果较差。在第1 天内，ISUT法的精度与UT 法相近，甚至整体略优于UT 法。而在第2天之后，ISUT法的精度相较于UT法有显著的降低。ISUT 算法的误差大于UT 算法主要是由计算解析解过程中的截断引起的。相较于SUT法，ISUT法的采样点修正过程弥补了一部分解析计算造成的精度降低。

图5 UT、SUT和ISUT法的归一化LAM随时间变化情况Fig.5 Normalized LAM of UT method，SUT method and ISUT method over time

在对比分析UT、ISUT 和GMMISUT 法时，以UT法的LAM 值为标准进行归一化，则3 种算法的归一化LAM 值随时间变化情况如图6 所示。可以看出，在第1 天内，ISUT 法与GMMISUT 法效果相近，GMMISUT 法的精度略微高于ISUT 法。第2 天之后，GMMISUT 法的精度相较于ISUT 法有显著提升，且逐渐优于UT 法。这是由于GMM 法的引入使得ISUT 法的精度得以提升。基于GMM 方法的性能与每个组件的协方差有关，协方差越小，结果越好。这一结论具有直观意义，因为当GMM 中每个分量的协方差趋于零时，每个分量就会变得越来越像脉冲函数。从而更轻松地近似各种PDF。结合图4 和图6 可知，GMMISUT 法在计算效率优于UT 法的同时，在大部分时段给出了优于UT法的结果。

图6 UT、ISUT和GMMISUT法的归一化LAM变化情况Fig.6 Normalized LAM of UT method，ISUT method and GMMISUT method over time

4 结论

对于轨道不确定性传播问题，提出基于半解析法和球形单边采样无迹变换的轨道不确定性快速估计方法。该方法的特色为采用球形单边采样无迹变换结合高斯摄动方程及其解析解进行采样点的传播和修正，从而实现轨道不确定性快速估计。其中解析解考虑了J2 摄动和大气阻力摄动的解析形式，被用于实现样本点的快速传播，加快计算效率；球形单边采样在保证二阶精度的同时进一步降低计算量；基于高斯摄动方程和采样点特性的采样点修正缓解了由采用解析解进行样本点传播导致的精度降低问题。将该方法与GMM 结合能够进一步提高精度。仿真分析表明，所提算法在保证合理精度的同时，大幅提升了不确定性估计的效率，适用于在轨不确定性快速传播，可为在轨避障分析等提供初始输入。