一种高斯混合协方差演化的碰撞概率计算方法及系统

技术领域

本发明属于航天技术领域，涉及一种高斯混合协方差演化的碰撞概率计算方法及系统。

背景技术

在碰撞预警工程中，需要预测未来几天内两个空间物体最接近的时间，预报轨道状态和协方差，计算碰撞概率。各种碰撞概率计算方法均假设轨道误差服从高斯分布，所以一般采用线性方法预报轨道协方差，该方法对非线性轨道动力系统采用一阶线性近似，忽略了高阶项的影响。事实上，当初始轨道不确定性较大或轨道预报周期较长时，预报的轨道协方差将不再服从高斯分布。因此，亟需一种针对初始轨道协方差较大或轨道预报周期较长时轨道协方差演化及碰撞概率计算方法，提高碰撞概率计算精度，降低碰撞预警虚警率。

发明内容

本发明的目的在于克服现有技术缺陷，提出了一种高斯混合协方差演化的碰撞概率计算方法。

为了实现上述目的，本发明提出了一种高斯混合协方差演化的碰撞概率计算方法，所述方法包括：

对待预测的两个空间物体分别基于高斯混合算法进行初始轨道协方差拟合，得到每个空间物体关于高斯单元的高斯混合概率密度函数；

对每个高斯单元进行协方差演化得到对应的高斯混合表示的轨道演化误差分布；

针对高斯混合表示的轨道演化误差分布通过高斯单元之间相互计算碰撞概率，并进行加权求和得到总的碰撞概率。

作为上述方法的一种改进，所述对待预测的两个空间物体分别基于高斯混合算法进行初始轨道协方差拟合，得到每个空间物体关于高斯单元的高斯混合概率密度函数；具体包括：

对每个空间物体基于划分的高斯单元，分别计算关于初始轨道的高斯概率密度；基于高斯混合算法对协方差进行拟合，得到该空间物体的高斯混合概率密度函数P(x)为：

其中，N为高斯单元的总数，p

作为上述方法的一种改进，所述对每个高斯单元进行协方差演化得到对应的高斯混合表示的轨道演化误差分布；具体包括：

对每个空间物体的每个高斯单元采用Unscented变换方法进行协方差演化，以第0时刻的初始轨道的状态量x的估计值和协方差为初值，以离散化非线性方程为模型，进行Unscented变换，直至得到第T时刻的估计值和协方差，实现高斯单元的均值和协方差传播获得高斯混合表示的轨道演化误差分布。

作为上述方法的一种改进，所述以第0时刻的初始轨道的状态量x的估计值和协方差为初值，以离散化非线性方程为模型，进行Unscented变换，直至得到第T时刻的估计值和协方差；具体包括：

根据第k-1时刻第i个高斯单元初始轨道的状态量x的估计值

构造Sigma点，得到下式：

其中，σ

根据计算所得Sigma点直接按照离散化非线性模型进行状态和观测预测，得到预测的对称分布Sigma样本点χ

χ

其中，f[·]表示与时间相关的非线性轨道预报函数，k表示第k时刻；

根据Sigma样本点χ

其中，χ

作为上述方法的一种改进，所述针对高斯混合表示的轨道演化误差分布通过高斯单元之间相互计算碰撞概率，并进行加权求和得到总的碰撞概率；具体包括：

根据下式得到在时间区间t

f

其中，P

P

其中，p

P

其中，θ,φ分别表示球形积分的方位角和高度角，v(r’,t)表示为相对积分速度，r’表示积分半径，满足下式：

其中，

σ

其中，P

对于球面积分采用Lebedev积分法，随时间变化的积分采用标准的积分公式；

根据下式得到总的碰撞概率P

其中，f

一种高斯混合协方差演化的碰撞概率计算系统，所述系统包括：高斯混合法轨道协方差拟合模块、协方差演化模块和碰撞概率计算模块；其中，

所述高斯混合法轨道协方差拟合模块，用于对待预测的两个空间物体分别基于高斯混合算法进行初始轨道协方差拟合，得到每个空间物体关于高斯单元的高斯混合概率密度函数；

所述协方差演化模块，用于对每个高斯单元进行协方差演化得到对应的高斯混合表示的轨道演化误差分布；

所述碰撞概率计算模块，用于针对高斯混合表示的轨道演化误差分布通过高斯单元之间相互计算碰撞概率，并进行加权求和得到总的碰撞概率。

一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述的方法。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行上述的方法。

与现有技术相比，本发明的优势在于：

空间目标轨道状态及不确定度演化是一个非线性动力学过程，对于初始轨道状态不确定度较大或预报期较长时，轨道误差分布会呈现非高斯性，本发明方法，对比已有技术，无论轨道不确定度是呈现高斯型分布还是非高斯型分布，采用高斯混合算法计算的碰撞概率达到蒙特卡洛计算的碰撞概率精度，高斯混合协方差计算方法提高碰撞概率计算精度。

附图说明

图1是本发明实施方式提供的一种高斯混合协方差演化的碰撞概率计算方法流程图；

图2是蒙特卡罗样本的径向沿迹平面投影及线性法、无迹卡尔曼滤波法和高斯混合方法的3σ误差等值线图。

具体实施方式

本发明方法是通过下述技术方案实现的：

一种高斯混合协方差演化的碰撞概率计算方法，其基本实施过程如下：

首先采用高斯混合算法对轨道初始轨道协方差拟合，将初始轨道协方差分解成较小的高斯单元；然后对每个高斯混合单元采用Unscented变换方法进行协方差演化，获得高斯混合表示的轨道演化误差分布；最后针对高斯混合表示的轨道误差采用高斯混合单元之间相互计算碰撞概率并加权求和来计算总的碰撞概率。

下面结合附图和实施例对本发明的技术方案进行详细的说明。

实施例1

如图1所示，本发明的实施例1提出了一种高斯混合协方差演化的碰撞概率计算方法，其具体步骤包括：

步骤一、通过有限个高斯概率密度函数加权求和来近似非线性概率密度函数，不同高斯单元的加权值由数值优化方法来确定。高斯混合概率密度函数表示为：

其中N为高斯概率密度单元函数的总数，α

把高斯混合函数中的每一项的均值和协方差演化来对非线性系统的误差进行量化。高斯混合与其它误差传播方法结合在一起使用，均值和协方差的演化通过UT转换方法实现。

具体为：

高斯混合的协方差演化时，首先找到一个平方根因子S，使SS

其中

通过选择特定的特征值和特征向量组成的平方根因子。因此，协方差矩阵P的分解为：

P＝VΛV

其平方根因子确定为S＝VΛ

P

其中v

步骤二、采用Unscented变换以k-1时刻系统状态x的估计值

具体为：

(1)构造Sigma点

式中，

(2)利用Sigma点进行非线性传播

根据计算所得Sigma点直接按照离散化非线性模型进行状态和观测预测，对应产生的Sigma样本点为：

χ

式中χ

(3)根据样本点计算状态均值和方差

步骤三、高斯混合的方法计算碰撞概率是计算两目标加权的高斯之间的碰撞概率

具体为：

高斯混合的方法计算碰撞概率基本公式是：

f

通过在两目标联合误差椭球面上积分计算瞬时碰撞概率变化率，并在两目标交会时间区间进行数值积分计算碰撞概率f

x

r

v

P

其在时间区间t

f

P

p

其中

σ

其中球面积分采用Lebedev积分法，随时间变化的积分可以使用标准的积分公式，如Newton-Cotes、Gauss-Konrod和Clenshaw-Curtis。

自此，就完成了高斯混合协方差演化的碰撞概率计算过程。

实施例2

本发明的实施例2提出了一种高斯混合协方差演化的碰撞概率计算系统，具体处理同实施例1的方法，所述系统包括：高斯混合法轨道协方差拟合模块、协方差演化模块和碰撞概率计算模块；其中，

所述高斯混合法轨道协方差拟合模块，用于对待预测的两个空间物体分别基于高斯混合算法进行初始轨道协方差拟合，得到每个空间物体关于高斯单元的高斯混合概率密度函数；

所述协方差演化模块，用于对每个高斯单元进行协方差演化得到对应的高斯混合表示的轨道演化误差分布；

所述碰撞概率计算模块，用于针对高斯混合表示的轨道演化误差分布通过高斯单元之间相互计算碰撞概率，并进行加权求和得到总的碰撞概率。

实施例3

一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现实施例1的方法。

实施例4

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序当被处理器执行时使所述处理器执行实施例1的方法。

仿真实例

仿真分析引入一个案例来验证高斯混合方法预报协方差计算碰撞概率，案例的初始协方差较大，案例中分别用线性化方法、无迹卡尔曼滤波方法、高斯混合方法+UT化方法、蒙特卡洛模拟方法对轨道状态及协方差进行预报。预报间隔为2天。其中蒙特卡洛方法根据初始轨道协方差模拟一百万个含有随机误差样本点，每一个样本点向前预报2天，得到预报时刻的轨道样本分布图形。蒙特卡洛方法是最准确的协方差预报方法，作为验证其他方法精度的依据。把以上预报模拟方法计算的协方差应用到两目标碰撞概率计算，得到不同协方差预报方法下的碰撞概率。

图2是蒙特卡罗样本的径向沿迹平面投影及线性法、无迹卡尔曼滤波法和高斯混合方法的3σ误差等值线。初始轨道协方差演化2天后，蒙特卡洛采样点分布已经不再服从高斯分布，而是在径向呈现月牙形分布。线性化方法和无迹卡尔曼滤波方法依然假设轨道状态分布为高斯型。无迹卡尔曼滤波方法可以描述轨道状态的总体分布区域，但是无法表示出具体分布形状。线性化方法甚至不能表示出轨道状态的分布区域。因此线性化方法和无迹卡尔曼滤波演化的协方差已经不能捕捉到这种轨道协方差分布特性。而高斯混合分布与蒙特卡洛样本点的分布拟和效果较好。表1是不同方法的碰撞概率计算结果。线性化方法、无迹卡尔曼滤波方法计算的碰撞概率小于蒙特卡洛方法计算的碰撞概率，其计算精度不能作为衡量空间目标碰撞预警的标准。高斯混合计算的碰撞概率结果可以看出高斯混合协方差计算方法使碰撞概率计算精度提高一个量级。

表1不同方法的碰撞概率计算结果

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。