CN104596520A - 一种基于半长轴滤波值连续性的航天器轨道机动检测方法 - Google Patents

Abstract

本发明公开一种基于半长轴滤波值连续性的航天器轨道机动检测方法，包括：基于运动学模型的轨道半长轴估计；基于轨道半长轴连续性的脉冲机动检测；基于轨道半长轴变化率的连续推力机动检测。采用本发明的技术方案，可以实现对轨道机动的快速准确检测与辨识，提高航天器轨道估计精度和可靠性。

Description

一种基于半长轴滤波值连续性的航天器轨道机动检测方法

技术领域

本发明涉及航空航天技术领域，尤其涉及一种基于半长轴滤波值连续性的航天器轨道机动检测方法。

技术背景

随着空间态势感知、攻防对抗和深空探测等领域任务的拓展，新型航天器的高机动、大范围运动特点使得对其进行实时高精度跟踪面临着更加严峻的挑战，尤其在空间攻防对抗等非合作目标跟踪中，未知的轨道机动极大的限制了航天器实时跟踪的精度，迫切需要有效的轨道机动检测和辨识技术。

非线性滤波是航天器实时轨道确定的主要方法，航天器轨道机动可以看作是其运动状态的突变，表现为一种系统性的偏差。因此轨道机动检测就是在滤波过程中或滤波后，引入必要的方法来检测或辨识轨道是否产生了偏离。其检测过程可以参考一般性的运动目标机动检测方法，又有其特殊性。航天器常见机动模式有脉冲机动、连续推力机动和bangbang机动，其中bangbang机动可以看作脉冲机动与连续推力机动的随机组合，而且目前应用相对较少，因此研究的热点主要集中于前两种机动的检测。检测算法实现的关键是检测统计量的构造，通常有两类方法，一是从测量信息或其某种特征的变化规律判断，二是基于滤波信息或输入估计的假设检验问题。后者利用了目标运动先验信息，且可以实现测量数据的融合，应用较广，但算法的检测和跟踪精度均依赖于模型的选择。鉴于航天器的空间轨道运动特征，有动力学和运动学两类模型可以利用。前者的约束力较强，常用于高精度的跟踪滤波器设计，且基于该滤波器信息设计的检测算法对微小机动检测能力较强，但由于对检测参数的依赖性也较大，检测滞后和虚警之间存在矛盾；后者对测量信息变化敏感，由于轨道机动会引起测量数据某些特征的变化(测量数据是轨道在测量空间的投影)，因此更利于航天器机动的检测，但其滤波结果的随机误差较大，不利于微小机动的检测与辨识。

常用机动检测和估计方法如等效噪声、输入检测与估计、模型转换等，但对于实时性和稳定性要求较高的实时轨道确定问题，上述方法对系统改造较大，检测过程所需数据积累和计算量相对较多，且适用性因实际系统状态噪声的不同而产生较大差异。此外，现有系统多以航天器位置与速度为滤波状态，因此机动检测手段也经常直接针对位置与速度进行，受测量异常值，特别是成片野值的影响较大，可靠性难以保证。

发明内容：

本发明要解决的技术问题是，提供一种基于轨道半长轴连续性变化监测的轨道机动检测方法，实现对轨道机动的快速准确检测与辨识。

为解决上述问题，本发明采用如下的技术方案：

一种基于半长轴滤波值连续性的航天器轨道机动检测方法包括：

步骤S1、基于运动学模型的轨道半长轴估计

在每个采样时刻，航天器跟踪系统利用对轨道运动学模型进行滤波，实时解算航天器位置与速度，进而获得该时刻的半长轴估计值；

步骤S2、基于轨道半长轴连续性的脉冲机动检测

利用窗口累积的轨道半长轴估值进行航天器脉冲机动检测，具体以下步骤：

步骤S2.1、利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴的统计特征；

步骤S2.2、采用逐点滑动的方式更新窗口数据，并通过检测新采样时刻的半长轴数据是否符合半长轴的统计特征征，进行脉冲机动检测；

步骤S2.3、若未发生机动，则利用新窗口数据计算半长轴的统计特征。

步骤S3、基于轨道半长轴变化率的连续推力机动检测

在脉冲机动检测基础上，利用窗口累积的轨道半长轴变化率估值进行航天器连续推力机动检测，具体包括以下步骤：

步骤S3.1、若未发生脉冲机动，则利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴变化率的统计特征；

步骤S3.2、采用逐点滑动的方式更新窗口数据，将窗口按时间分成前后两段，分别计算半长轴变化率，并通过检测两个变化率的差异是否符合半长轴变化率的统计特征，进行连续推力机动检测；

步骤S3.3、若未发生机动，则利用新窗口数据计算半长轴变化率的统计特征。

本发明提出利用轨道半长轴滤波值连续性的机动检测方法，具有实现简单，检测精度高、速度快和稳健性强特点，具体优点如下：

1、仅需在跟踪滤波器中引入简单的半长轴统计特征计算和连续性变化检测环节，即可分别实现脉冲机动和连续推力机动的检测，对现有跟踪系统的改造较小，方法实现方便且适应性较强；

2、采用逐点滑动窗口方式，不存在累积数据自身引入的解算滞后，且仅使用了简单的多点差分算法，以及半长轴与半长轴变化率的统计特征计算，检测过程所需计算量很小，保证了机动检测的速度；

3、检测过程基于窗口内的轨道半长轴估值进行，充分利用了多个时刻的测量数据，可以降低测量异常的影响，提高检测精度的同时，保障机动检测的可靠性；

附图说明

图1是本发明的算法流程示意图；

图2是实施算例中航天器轨道的脉冲机动检测示意图；

图3是实施算例中航天器轨道的连续推力机动检测示意图。

具体实施方式

如图1所示，本发明公开一种基于半长轴滤波值连续性的航天器轨道机动检测方法，包括以下步骤：

步骤S1、基于运动学模型的航天器轨道半长轴估计

在每个采样时刻t_k，航天器跟踪系统利用对轨道运动学模型进行滤波，实时解算航天器位置与速度，进而获得该时刻的半长轴估计值a_k。

从计算效率和精度出发，所述运动学模型为当前统计模型(Current Statistic Model，CS)。记t_k时刻航天器惯性系位置为X_k＝(x_k，y_k，z_k)^T，x_k，y_k，z_k分别为在惯性系x，y，z轴的分量。当前统计模型认为机动加速度在各个方向是解耦的，并且在时间轴上是一阶时间相关的随机过程。记 ${\tilde{X}}_k={[x_k,{\stackrel{\·}{x}}_k,{\stackrel{\·\·}{x}}_k]}^T,{\tilde{T}}_k={[y_k,{\stackrel{\·}{y}}_k,{\stackrel{\·\·}{y}}_k]}^T,{\tilde{Z}}_k={[z_k,{\stackrel{\·}{z}}_k,{\stackrel{\·\·}{z}}_k]}^T,$ 则状态转移方程

$\left[\begin{array}{c}{\tilde{X}}_{k+1}\\ {\tilde{Y}}_{k+1}\\ {\tilde{Z}}_{k+1}\end{array}\right]=\left[\begin{array}{ccc}{B}_k& 0_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& B_k& 0_{3\×3}\\ 0_{3\×3}& 0_{3\×3}& B_k\end{array}\right]\left[\begin{array}{c}{\tilde{X}}_k\\ {\tilde{Y}}_k\\ {\tilde{Z}}_k\end{array}\right]+\left[\begin{array}{ccc}{\mathrm{\Γ}}_k& 0_{1\×3}& 0_{1\×3}\\ 0_{1\×3}& {\mathrm{\Γ}}_k& 0_{1\×3}\\ 0_{1\×3}& 0_{1\×3}& {\mathrm{\Γ}}_k\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]---\left(1\right)$

其中，分别为航天器t_k时刻x方向的速度和加速度。分别为航天器t_k时刻y方向的速度和加速度，分别为航天器k时刻z方向的速度和加速度 $B_k=\left[\begin{array}{ccc}1& \mathrm{dt}& {0.5\mathrm{dt}}^2\\ & 1& \mathrm{dt}\\ & & 1\end{array}\right],$ ${\mathrm{\Γ}}_k=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right],$ dt为采样时间间隔， $\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]$ 为机动加速度引入导致的状态噪声向量，其方差由机动加速度的概率密度函数确定。当前统计模型假设其满足修正的瑞利分布，三个方向相对独立，以x方向为例，t_k时刻ω_x的方差为

${\mathrm{\&sigma;}}_{x,k}^2=\left\{\begin{array}{cc}{({\stackrel{\&CenterDot;\&CenterDot;}{x}}_{\max }-{\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}}_{k-1})}^2/\mathrm{\&pi;}& {\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}}_{k-1}>0\\ (4-\mathrm{\&pi;}){({\stackrel{\&CenterDot;\&CenterDot;}{x}}_{-\max }-{\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}}_{k-1})}^2/\mathrm{\&pi;}& {\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}}_{k-1}<0\end{array}\right.---\left(2\right)$

其中，表示t_k-1时刻加速度估计值，与分别为加速度正向和负向最大值，一般根据目标运动规律先验知识获取。其它方向类似。获得轨道测量数据后，根据测量的非线性性，这里采用无迹滤波(the unscented filter，UKF)进行位置X_k＝(x_k，y_k，z_k)^T与速度的实时估计，具体可过程可参考文献《A new extension of the Kalmanfilter t0nonlinear systems》(Julier S J and Uhlmann J K.Proc AeroSense，Orlando，1997：54-65.)。进而由航天器位置X_k＝(x_k，y_k，z_k)^T与速度的估值计算其轨道半长轴a_k如下

$a_k={[\frac{2}{{(x_k+y_k^2+z_k^2)}^{1/2}}-\frac{({\stackrel{\&CenterDot;}{x}}_k^2+{\stackrel{\&CenterDot;}{y}}_k^2+{\stackrel{\&CenterDot;}{z}}_k^2)}{\mathrm{\&mu;}}]}^{-1}---\left(3\right)$

其中，μ为地球引力常数，μ＝398600.4418e⁹。

步骤S2、基于半长轴连续性的航天器轨道脉冲机动检测

利用窗口累积的轨道半长轴估值进行航天器脉冲机动检测，其检测过程包括如下步骤：

步骤S2.1、利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴的统计特征；

记t_i时刻的半长轴估值为a_i，若已连续获得了N个采样时刻的轨道估值{a_i，a_i+1，…，a_i+N-1}，N为积累窗口长度，采用序列差分检测来判断这些点是否满足初始化条件，计算公式如下：

$|a_{k-4}-{4a}_{k-3}+{6a}_{k-2}-{4a}_{k-1}+a_k|>\mathrm{\&beta;}\&CenterDot;{\mathrm{\&sigma;}}_a^0,i+4\&le;k\&le;i+N-1---\left(4\right)$

其中，β为第一阈值，对于圆轨道可取值1.5～3，对于偏心率较大的轨道可取值3～6；为半长轴的先验方差。

当窗口内任意相邻五点均满足公式(4)，则初始积累完成，即认为这连续N个采样时刻对应的轨道比较平稳，无机动存在，可以作为后续检测的基础。

若初始积累完成，在计算窗口{a_i，a_i+1，…，a_i+N-1}内，按照公式(5)和公式(6)计算窗口数据{a_i，a_i+1，…，a_i+N-1}的中值和标准差中值和标准差为半长轴的统计特征，并转入Step 2.2；否则，窗口向前滑动一个采样点，即取{a_i+1，a_i+2，…，a_i+N}，重新初始积累。

${\stackrel{\&OverBar;}{a}}_i=\mathrm{med}(a_i,a_{i+1},\&CenterDot;\&CenterDot;\&CenterDot;,a_{i+N-1})---\left(5\right)$

${\mathrm{\&sigma;}}_a^i=\mathrm{std}(a_i,a_{i+1},\&CenterDot;\&CenterDot;\&CenterDot;,a_{i+N-1})---\left(6\right)$

其中，med(·)表示中值计算函数，std(·)表示标准差计算函数。

步骤S2.2、采用逐点滑动的方式更新窗口数据，并通过检测新采样时刻的半长轴数据是否符合半长轴的统计特征，进行脉冲机动检测；

当滤波获得窗口右侧采样时刻半长轴估值a_i+N后，若下式成立，则认为该f+N时刻航天器潜在机动

$|a_{i+N}-{\stackrel{\&OverBar;}{a}}_i|>\mathrm{\&alpha;}\&CenterDot;{\mathrm{\&sigma;}}_a^i---\left(7\right)$

其中，α为第二阈值，其反映了检测的灵敏度，通常取值为3～6。

为了减小测量异常的影响，保证检测可靠性，需要连续多个采样点超出阈值时才判断机动发生。这里选择四分之一的窗口数据，即当N/4个采样点满足公式(7)时，认为航天器在该时段内发生了脉冲机动，并将第一个潜在机动的时刻i+N做为脉冲机动起始时刻，窗口向右滑动一个采样点，转Step 2.1，重新开始窗口初始积累，并将初始积累结束的时刻作为脉冲机动的结束时刻；否则，窗口向右滑动一个采样点，转Step2.3。

步骤S2.3、若未发生机动，则利用新窗口数据计算半长轴的统计特征。

此时认为航天器并未处于脉冲机动段落，需要重新计算半长轴的统计特征，即中值和标准差，并转入Step 2.2，开始下一时刻的机动检测。为了更有效的抑制随机误差和异常值的影响，中值计算仍采用公式(5)，即但标准差的计算引入如下的平滑方法：

${\mathrm{\&sigma;}}_a^{i+1}=\mathrm{\&gamma;}\&CenterDot;{\mathrm{\&sigma;}}_a^i+(1-\mathrm{\&gamma;})\&CenterDot;\mathrm{std}(a_{i+1},a_{i+2},\&CenterDot;\&CenterDot;\&CenterDot;,a_{i+N})---\left(8\right)$

其中，γ为渐消因子，取值越大则平滑效果越好，但会降低对轨道半长轴变化的敏感性，通常计算时取值0.8～0.95。

步骤S3、基于半长轴变化率连续性的航天器轨道连续推力机动检测

在脉冲机动检测基础上，利用窗口累积的轨道半长轴变化率估值进行航天器连续推力机动检测，其检测过程包括如下步骤：

步骤S3.1、在脉冲机动检测基础上，若未发生脉冲机动，则利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴变化率的统计特征；

设积累窗口的长度为N，采用二次多项式拟合方式计算轨道半长轴在窗口{a_i，a_i+1，…，a_i+N-1}内的变化率。具体的，设采样时刻为T＝[t_i，t_i+1，…，t_i+N-1]^T，则拟合系数为C＝(D^TD)^-1D^TL，变化率的计算方法为

$\stackrel{\&CenterDot;}{L}=\mathrm{EC}=E{\left(D^{T}D\right)}^{-1}{D}^{T}L---\left(9\right)$

其中，表示a_i的变化率， $D=\left[\begin{array}{ccc}{t}_i^2& t_i& 1\\ .& .& .\\ .& .& .\\ .& .& .\\ t_{i+N-1}^2& t_{i+N-1}& 1\end{array}\right],E=\left[\begin{array}{ccc}{2t}_i& 1& 0\\ .& .& .\\ .& .& .\\ .& .& .\\ {2t}_{i+N-1}& 1& 0\end{array}\right].$

进一步，采用序列差分检测方法，判断是否满足初始化条件，差分检测公式为：

$|{\stackrel{\&CenterDot;}{a}}_{k-4}-4{\stackrel{\&CenterDot;}{a}}_{k-3}+6{\stackrel{\&CenterDot;}{a}}_{k-2}-4{\stackrel{\&CenterDot;}{a}}_{k-1}+{\stackrel{\&CenterDot;}{a}}_k|>\mathrm{\&beta;}\&CenterDot;{\mathrm{\&sigma;}}_{\stackrel{\&CenterDot;}{a}}^0,i+4\&le;k\&le;i+N---\left(10\right)$

其中，β为第一阈值，对于圆轨道可取值1.5～3，对于偏心率较大的轨道可取值3～6；为半长轴变化率的先验方差。

当窗口内任意相邻五点均满足公式(10)，则初始积累完成，即认为这连续N个采样时刻对应的轨道比较平稳，无连续推力机动存在，可以作为后续检测的基础。

初始积累完成后，计算半窗口数据的标准差标准差为半长轴变化率的统计特征，其中，[·]表示取整运算，并转入Step 3.2。

标准差计算公式为：

${\mathrm{\&sigma;}}_{\stackrel{\&CenterDot;}{a}}^i=\mathrm{std}({\stackrel{\&CenterDot;}{a}}_i,{\stackrel{\&CenterDot;}{a}}_{i+1},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&CenterDot;}{a}}_{i+N-1})---\left(11\right)$

需要指出的是，满足积累条件后只能说明目前航天运动特征相对平稳，其可能处于机动段落，也可能处于连续常值推力机动的段落。

步骤S3.2、采用逐点滑动的方式更新窗口数据，将窗口按时间分成前后两段，分别计算半长轴变化率，并通过检测两个变化率的差异是否符合半长轴变化率的统计特征，进行连续推力机动检测；

当滤波获得窗口右侧采样时刻半长轴估值a_i+N后，将数据分为 $A=\{{\stackrel{\&CenterDot;}{a}}_i,{\stackrel{\&CenterDot;}{a}}_{i+1},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&CenterDot;}{a}}_{i+\left[0.5N\right]}\}$ 和 $B=\{{\stackrel{\&CenterDot;}{a}}_{i+\left[0.5N\right]+1},{\stackrel{\&CenterDot;}{a}}_{i+\left[0.5N\right]+2},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&CenterDot;}{a}}_{i+N}\}$ 两部分，若下式成立，则认为时刻f+[0.5N]航天器的机动特性发生变化(连续机动开始或结束)

$|\mathrm{med}\left(A\right)-\mathrm{med}\left(B\right)|>\mathrm{\&alpha;}\&CenterDot;{\mathrm{\&sigma;}}_{\stackrel{\&CenterDot;}{a}}^i---\left(12\right)$

其中，[·]表示取整运算，α为第二阈值，反映了检测的灵敏度，而为了兼顾测量异常的影响，通常取值为3～6。

公式(12)只能给出机动特性变化的判断，进一步，可以通过前后窗口比对判断时刻f+[0.5N]是机动起始还是终止。具体的如下

为了保证检测可靠性，仍采用步骤S2中Step 2.2的多点判读方式。当连续有四分之一的窗口数据(即N/4个采样点)满足公式(12)时，则认为航天器在该时段内的机动特性发生了改变，将第一个时刻f+[0.5N]作为变化起始时刻，窗口向右滑动一个采样点，转Step3.1，重新开始窗口初始积累；否则，窗口向右滑动一个采样点，转Step3.3。

步骤S3.3、若未发生机动，则利用新窗口数据计算半长轴变化率的统计特征。

当判断航天器未发生机动后，重新计算半长轴变化率计算方法见公式(9)，并采用下(公式14)计算半长轴变化率标准差

${\mathrm{\&sigma;}}_{\stackrel{\&CenterDot;}{a}}^{i+1}=\min \{\mathrm{std}({\stackrel{\&CenterDot;}{a}}_{i+1},{\stackrel{\&CenterDot;}{a}}_{i+2},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&CenterDot;}{a}}_{i+\left[0.5N\right]+1}),\mathrm{std}({\stackrel{\&CenterDot;}{a}}_{i+\left[0.5N\right]+3},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&CenterDot;}{a}}_{i+N})\}---\left(14\right)$

统计特征更新后，转入Step 3.2，开始下一时刻的机动检测。

本发明采用的方法可以实现航天器脉冲和连续推力机动的检测与辨识，提高航天器跟踪的精度和可靠性。以某典型卫星轨道为例，给出本发明方法的一个验证实例。卫星的初始轨道根数和机动特性如表1所示，图2、3给出了航天器机动检测的示意图，检测结果如表2所示。可见在本算例中，方法的机动检测率达到了100％，其中脉冲机动开始时刻的检测并无明显滞后，机动阶数时刻相差在几秒以内；连续推力机动的开始和结束时刻滞后小于100s，远低于机动持续时间(1/10以下)。因此，算例验证了本发明方法的可行性和有效性。

表1：

表2：

说明书中描述的只是该发明的具体实施方式。虽然结合附图描述了本发明的实施方式，但是本领域内熟练的技术人员可以在所附权利要求的范围内做出各种变形或修改。

Claims (2)

1.一种基于半长轴滤波值连续性的航天器轨道机动检测方法，其特征在于，包括：

步骤S1、基于运动学模型的轨道半长轴估计

在每个采样时刻，航天器跟踪系统利用对轨道运动学模型进行滤波，实时解算航天器位置与速度，进而获得该时刻的半长轴估计值；

步骤S2、基于轨道半长轴连续性的脉冲机动检测

利用窗口累积的轨道半长轴估值进行航天器脉冲机动检测，具体以下步骤：

步骤S2.1、利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴的统计特征；

步骤S2.2、采用逐点滑动的方式更新窗口数据，并通过检测新采样时刻的半长轴数据是否符合半长轴的统计特征，进行脉冲机动检测；

步骤S2.3、若未发生机动，则利用新窗口数据计算半长轴的统计特征。

步骤S3、基于轨道半长轴变化率的连续推力机动检测

在脉冲机动检测基础上，利用窗口累积的轨道半长轴变化率估值进行航天器连续推力机动检测，具体包括以下步骤：

步骤S3.1、若未发生脉冲机动，则利用差分检测的方法实现窗口半长轴估值数据的累积，并计算半长轴变化率的统计特征；

步骤S3.2、采用逐点滑动的方式更新窗口数据，将窗口按时间分成前后两段，分别计算半长轴变化率，并通过检测两个变化率的差异是否符合半长轴变化率的统计特征，进行连续推力机动检测；

步骤S3.3、若未发生机动，则利用新窗口数据计算半长轴变化率的统计特征。

2.如权利要求1所述的一种基于半长轴滤波值连续性的航天器轨道机动检测方法，其特征在于，步骤S1中，设t_k时刻航天器惯性系位置为X_k＝(x_k，y_k，z_k)^T，其中，x_k，y_k，z_k分别为在惯性系x，y，z轴的分量，设t_k时刻航天器惯性系速度其中，分别为航天器t_k时刻x方向的速度，分别为航天器t_k时刻y方向的速度，分别为航天器k时刻z方向的速度；半长轴估计值a_k计算公式如下：

$a_k={[\frac{2}{{(x+y_k^2+z_k^2)}^{1/2}}-\frac{({\stackrel{.}{x}}_k^2+{\stackrel{.}{y}}_k^2+{\stackrel{.}{z}}_k^2)}{\mathrm{\&mu;}}]}^{-1}$

其中，μ为地球引力常数，μ＝398600.4418e⁹。