CN111985077B - 一种航天器外弹道跟踪数据斑点型野值识别和修正方法 - Google Patents
一种航天器外弹道跟踪数据斑点型野值识别和修正方法 Download PDFInfo
- Publication number
- CN111985077B CN111985077B CN202010664267.6A CN202010664267A CN111985077B CN 111985077 B CN111985077 B CN 111985077B CN 202010664267 A CN202010664267 A CN 202010664267A CN 111985077 B CN111985077 B CN 111985077B
- Authority
- CN
- China
- Prior art keywords
- spacecraft
- trajectory tracking
- tracking data
- formula
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000013499 data model Methods 0.000 claims abstract description 11
- 238000013178 mathematical model Methods 0.000 claims abstract description 8
- 238000005070 sampling Methods 0.000 claims description 34
- 238000004422 calculation algorithm Methods 0.000 claims description 31
- 238000005259 measurement Methods 0.000 claims description 23
- 239000011159 matrix material Substances 0.000 claims description 9
- 230000009897 systematic effect Effects 0.000 claims description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000009825 accumulation Methods 0.000 claims description 3
- 238000012937 correction Methods 0.000 abstract description 5
- 238000013213 extrapolation Methods 0.000 description 7
- 230000002159 abnormal effect Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000011217 control strategy Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000003909 pattern recognition Methods 0.000 description 2
- 235000015842 Hesperis Nutrition 0.000 description 1
- 235000012633 Iberis amara Nutrition 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Navigation (AREA)
Abstract
本发明公开了一种航天器外弹道跟踪数据斑点型野值识别和修正方法,具体为:首先,建立含有斑点型野值航天器外弹道跟踪数据时间序列数学模型;对航天器外弹道跟踪数据进行分割;之后求取各段航天器外弹道跟踪数据的标准差;对航天器外弹道跟踪数据模型的参数进行估计;识别和修正航天器各段外弹道跟踪数据中的斑点型野值;最后,将识别和修正野值后的各段航天器外弹道跟踪数据进行拼接。本发明的方法能够自动地识别出航天器外弹道跟踪数据中的斑点型野值,并进行高精度的修正,免去了专业工作人员通过人眼方式判断和修正斑点型野值的时间消耗和人员消耗。
Description
技术领域
本发明属于航空航天跟踪数据处理技术领域,具体涉及一种航天器外弹道跟踪数据斑点型野值识别和修正方法。
背景技术
由于跟踪环境、设备工况和操作过程中各种突发性异变因素的影响,在运载火箭、人造卫星和导弹等各类飞行器跟踪测量数据中,不可避免地存在着诸如异常数据、异常点(Outliers)等明显偏离大部分正常数据点(Inliers)所呈现变化趋势的小部分数据,这部分异常数据称之为野值。几十年航天发射测控的实践和国内外大量的理论研究、仿真分析、应用实践表明:野值对航天测控等领域广泛采用的经典处理方法有极大的破坏作用,轻则影响计算结果的可靠性,重则导致计算结果失真、控制策略错误,甚至整个算法崩溃,直接威胁到飞行器轨道确定的可靠性和控制策略的安全性,因此必须提前将它们进行识别。
外弹道跟踪测量数据中的野值根据其表现形式分为孤立型野值和斑点型野值两类,其中孤立型野表现为采样时间上不连续的一个孤立跳点。对于斑点型野值,因某种相关性影响,表现为采样时间上连续的多个连续跳点。传统的野值剔除方法(包括莱特准则、罗曼诺夫斯基准则、格拉布斯准则及肖维勒准则等)对于弹道数据的统计差特性具有约束性,导致这些方法在使用时受限;中值滤波差分法可有效剔除斑点型野值,但剔除效果不稳定,会出现数据失真或缺失的情况;外推拟合法可有效剔除孤立型和斑点型野值,但斑点型野值剔除后会出现数据缺失。目前,还没有一种有效的方法对外弹道跟踪测量数据中的斑点型野值进行识别和修正。
发明内容
本发明的目的是提供一种航天器外弹道跟踪数据斑点型野值识别和修正方法,能够利用抗野值识别算法从含有斑点型野值的航天器外弹道跟踪数据估计出正常数据的模型参数,能够设定动态阈值将不符合正常数据模型参数的野值进行识别和修正。
本发明所采用的技术方案是,一种航天器外弹道跟踪数据斑点型野值识别和修正方法,具体按照以下步骤实施:
步骤1,建立含有斑点型野值航天器外弹道跟踪数据时间序列数学模型;
步骤2,对航天器外弹道跟踪数据进行分割;
步骤3,求取各段航天器外弹道跟踪数据的标准差;
步骤4,对各段航天器外弹道跟踪数据模型的参数进行估计;
步骤5,识别和修正各段航天器外弹道跟踪数据中的斑点型野值;
步骤6,将识别和修正野值后的各段航天器外弹道跟踪数据进行拼接。
本发明的特点还在于,
步骤1中,具体为:
在时间T1,T2,…,Ti,…,Tn上采样得到一组航天器外弹道跟踪数据Y1,Y2,…,Yi,…,Yn,即,时刻Ti(i=1…n)的航天器外弹道跟踪采样数据为Yi(i=1…n),时刻Ti(i=1…n)航天器外弹道跟踪数据Yi(i=1…n)的数学模型函数,如式(1)所示:
Yi=Xi+Si+ωi+δi (1);
式(1)中,Si(i=1…n)为系统误差;ωi(i=1…n)为随机误差项,是标准差为σ的高斯白噪声,{εi}~N(0,σ),δi表示野值,当某个采样时刻出现野值时,δi的绝对值远大于高斯白噪声的标准差σ,满足σ<<|δi|,否则δi=0,Xi为外弹道跟踪测量数据时刻Ti(i=1…n)的真值,被表示为三次多项式,如式(2)所示;
Xi=θ3(Ti)3+θ2(Ti)2+θ1(Ti)1+θ0 (2);
式(2)中,θ0 θ1 θ2 θ3依次为三次多项式Xi的常数项、一次项系数、二次项系数和三次项系数。
步骤2中,具体为:
假定在时间t1,t2,…,ti,…,tn上采样得到一组航天器外弹道跟踪数据y1,y2,…,yi,…,yn,即,时刻ti(i=1…n)的航天器外弹道跟踪采样数据为yi(i=1…n)。将y1,y2,…,yi,…,yn分割成L段(L为整数),且满足分割后每段航天器外弹道跟踪数据序列的数据长度为l(l>4)的整数,则第k(1≤k≤L)段航天器外弹道跟踪数据向量mk如式(3)所示;
mk=[y(k-1)×l+1,y(k-1)×l+2,…,yk×l],1≤k≤L (3);
进一步,第k(1≤k≤L)段航天器外弹道跟踪数据mk中第j个测量数据如式(4)所示;
步骤3中,具体为:
第k(1≤k≤L)段航天器外弹道跟踪数据集合mk的第j个测量数据表示为式(1)的形式,令i=z,Yi=yz,Ti=tz,则式(1)变形为式(5);
其中,z=((k-1)×l+j),(1≤j≤l),tz代表第z个采样时刻,依次为第k段航天器外弹道跟踪数据mk对应的三次多项式的常数项、一次项系数、二次项系数和三次项系数,Sk/>依次是第k段航天器外弹道跟踪数据mk在时刻z的系统误差和随机误差。
第k(1≤k≤L)段航天器外弹道跟踪数据mk中一共有l个数据点y(k-1)×l+1,y(k-1)×l+2,…,yk×l,将式(5)全部化为方程组的形式,如式(6)所示;
将式(6)化为矩阵形式,如式(7)所示;
Mk=Tkθk+ωk (7);
式(7)中的各个矩阵和向量展开如式(8)所示;
式(7)中,Mk是第k段航天器外弹道跟踪数据mk的向量形式,Tk是第k段航天器外弹道跟踪数据对应的时间矩阵,θk是第k段航天器外弹道跟踪数据模型参数向量,ωk是第k段航天器外弹道跟踪数据模型随机误差向量形式;
定义式(9)为式(7)第k(1≤k≤L)段航天器外弹道跟踪数据最小二乘拟合的损失函数J(θk);
则式(9)中损失函数J(θk)的偏导数如式(10)所示,并令其为0向量;
求解式(10)得到第k段航天器外弹道跟踪数据模型参数向量θk的估计值如式(11)所示;
则第k段外弹道跟踪测量数据mk的随机误差的标准差的估计值,如式(12)所示;
则所有L段外弹道跟踪测量数据集合{m1,m2,…mk…,mL-1,mL}对应的随机误差的标准差估计值集合为
步骤4中,其具体步骤是:
已知第k(1≤k≤L)段共l个航天器外弹道跟踪数据采样点集合,如式(13)所示;
Φk={(t(k-1)×l+1,y(k-1)×l+1),(t(k-1)×l+2,y(k-1)×l+2),…,(tk×l,yk×l)},1≤k≤L(13);
建立抗野值估计模型,如式(14)所示;
φ(Θ)=Θ3(t)3+Θ2(t)2+Θ1t+Θ0+S (14);
其中,Θ3 Θ2 Θ1 Θ0 S依次是航天器外弹道跟踪数据抗野值估计模型的三次项系数,二次项系数,一次项系数,常数项系数和系统误差,t是航天器外弹道跟踪数据对应的采样时间变量;
样本子集个数:u=4,样本集合:U=Φk,一致性集合判别阈值:sigma=2×min(R),min(R)是集合随机误差的标注差估计值的集合的最小元素;
算法迭代次数:
步骤4.1:从样本集U中随机抽取u个样本点,作为子集ηc;
步骤4.2:将样本子集ηc中的u个样本点依次代入具体模型中,解方程组,得到模型参数
步骤4.3:计算出本次模型参数对应的支持度函数/>
对于第k(1≤k≤L)段航天器外弹道跟踪数据第c次执行后的支持度函数有如下定义:步骤4.2中第c次执行后得到模型参数/>时,支持度函数/>就是样本集U中所有l个样本到/>的距离小于一致性集合判断阈值sigma的样本点的个数,如式(15)所示;
其中,函数ρ(·)定义为积累函数,其表达形式如式(16)所示;
其中,Uc表示样本集合U中的第c个样本点,因此,第k(1≤k≤L)段航天器外弹道跟踪数据经过ξ次执行后,一共可以得到ξ个模型参数,即,它们一一对应ξ个支持度函数/>
步骤4.4,设中的最大值为/>则将与/>对应的模型参数作为利用抗野值模型估计迭代算法估计得到的第k(1≤k≤L)段航天器外弹道跟踪数据参数的估计值;
因此,第k(1≤k≤L)段航天器外弹道跟踪数据的估计值为式(17)所示;
步骤5中,具体为:
将第k(1<k≤L)段外弹道跟踪测量数据(t(k-1)×l+1,y(k-1)×l+1),(t(k-1)×l+2,y(k-1)×l+2),…,(tk×l,yk×l)代入公式(5)可得公式(18);
将式(17)和(18)作差并取绝对值,如式(19)所示;
定义:因此式(19)可化为式(20);
当时,将yz识别为野值点,并用/>修正它,否则将其识别为正常点。
步骤6中,具体为:设识别和修正野值后的第k(1≤k≤L)段航天器外弹道跟踪数据向量如式(21)所示;
对应采样时刻向量Tk=[t(k-1)×l+1,t(k-1)×l+2,…,tk×l];则完整弹道被表示为对应的采样时刻[T1,T2,…,Tk,…,TL]。
本发明的有益效果是,
本发明方法能够从含有野值(outliers)的数据中识别出正常数据(inliers)的模型参数,自适应地给出野值识别的动态阈值,利用动态阈值将不符合正常数据模型参数的野值进行识别和修正,解决了现有野值识别和修正算法无法有效识别和修正斑点型野值的缺陷。本发明的方法能够自动地识别出航天器外弹道跟踪数据中的斑点型野值,并进行高精度的修正,免去了专业工作人员通过人眼方式判断和修正斑点型野值的时间消耗和人员消耗。
附图说明
图1是本发明一种航天器外弹道跟踪数据斑点型野值识别和修正方法的流程图;
图2是本发明方法中对航天器外弹道跟踪数据模型的参数进行估计的流程图;
图3是四种斑点型野值识别算法拟合曲线对比图;
图4是四种斑点型野值识别算法误差结果对比图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种航天器外弹道跟踪数据斑点型野值识别和修正方法,如图1所示,具体按照以下步骤实施:
步骤1,建立含有斑点型野值航天器外弹道跟踪数据时间序列数学模型;
具体为:在时间T1,T2,…,Ti,…,Tn(时间序列长度5~10秒之间)上采样得到一组航天器外弹道跟踪数据Y1,Y2,…,Yi,…,Yn,即,时刻Ti(i=1…n)的航天器外弹道跟踪采样数据为Yi(i=1…n),由航空工程实践经验,时刻Ti(i=1…n)航天器外弹道跟踪数据Yi(i=1…n)的数学模型函数,如式(1)所示:
Yi=Xi+Si+ωi+δi (1);
式(1)中,Si(i=1…n)为系统误差,一般为常数;ωi(i=1…n)为随机误差项,是标准差为σ的高斯白噪声,{εi}~N(0,σ),δi表示野值,其出现的时刻不可预测,当某个采样时刻出现野值时,δi的绝对值远大于高斯白噪声的标准差σ,满足σ<<|δi|,否则δi=0,Xi为外弹道跟踪测量数据时刻Ti(i=1…n)的真值,一般被表示为三次多项式,如式(2)所示;
Xi=θ3(Ti)3+θ2(Ti)2+θ1(Ti)1+θ0 (2);
式(2)中,θ0 θ1 θ2 θ3依次为三次多项式Xi的常数项、一次项系数、二次项系数和三次项系数;
步骤2,对航天器外弹道跟踪数据进行分割;
假定在时间t1,t2,…,ti,…,tn上采样得到一组航天器外弹道跟踪数据y1,y2,…,yi,…,yn,即,时刻ti(i=1…n)的航天器外弹道跟踪采样数据为yi(i=1…n)。将y1,y2,…,yi,…,yn分割成L段(L为整数),且满足分割后每段航天器外弹道跟踪数据序列的数据长度为l(l>4)的整数,则第k(1≤k≤L)段航天器外弹道跟踪数据向量mk如式(3)所示;
mk=[y(k-1)×l+1,y(k-1)×l+2,…,yk×l],1≤k≤L (3);
进一步,第k(1≤k≤L)段航天器外弹道跟踪数据mk中第j个测量数据如式(4)所示;
步骤3,求取各段航天器外弹道跟踪数据的标准差;
具体为:根据工程经验,第k(1≤k≤L)段航天器外弹道跟踪数据集合mk的第j个测量数据可以表示为式(1)的形式,令i=z,Yi=yz,Ti=tz,则式(1)变形为式(5);
其中,z=((k-1)×l+j),(1≤j≤l),设δi=0,tz代表第z个采样时刻。依次为第k段航天器外弹道跟踪数据mk对应的三次多项式的常数项、一次项系数、二次项系数和三次项系数,Sk/>依次是第k段航天器外弹道跟踪数据mk在时刻z的系统误差(常量)和随机误差。
第k(1≤k≤L)段航天器外弹道跟踪数据mk中一共有l个数据点y(k-1)×l+1,y(k-1)×l+2,…,yk×l,将式(5)全部化为方程组的形式,如式(6)所示;
将式(6)化为矩阵形式,如式(7)所示;
Mk=Tkθk+ωk (7);
式(7)中的各个矩阵和向量展开如式(8)所示;
式(7)中,Mk是第k段航天器外弹道跟踪数据mk的向量形式(已知量),Tk是第k段航天器外弹道跟踪数据对应的时间矩阵(已知量),θk是第k段航天器外弹道跟踪数据模型参数向量(未知量),ωk是第k段航天器外弹道跟踪数据模型随机误差向量形式(未知量);
定义式(9)为式(7)第k(1≤k≤L)段航天器外弹道跟踪数据最小二乘拟合的损失函数J(θk);
式(9)中,z=((k-1)×l+j),tz为第z个采样数据对应的时间点;则式(9)中损失函数J(θk)的偏导数如式(10)所示,并令其为0向量;
求解式(10)可以得到第k段航天器外弹道跟踪数据模型参数向量θk的估计值如式(11)所示;
则第k段外弹道跟踪测量数据mk的随机误差的标准差的估计值,如式(12)所示;
则所有L段外弹道跟踪测量数据集合{m1,m2,…mk…,mL-1,mL}对应的随机误差的标准差估计值集合为
步骤4,对各段航天器外弹道跟踪数据模型的参数进行估计,如图2所示,其具体步骤是:
已知第k(1≤k≤L)段共l个航天器外弹道跟踪数据采样点集合,如式(13)所示;
Φk={(t(k-1)×l+1,y(k-1)×l+1),(t(k-1)×l+2,y(k-1)×l+2),…,(tk×l,yk×l)},1≤k≤L(13);
利用抗野值模型估计迭代算法估计第k(1≤k≤L)段航天器外弹道跟踪数据参数开始前初始化迭代算法参数如下:
建立抗野值估计模型,如式(14)所示;
φ(Θ)=Θ3(t)3+Θ2(t)2+Θ1t+Θ0+S (14);
其中,Θ3 Θ2 Θ1 Θ0 S依次是航天器外弹道跟踪数据抗野值估计模型的三次项系数,二次项系数,一次项系数,常数项系数和系统误差,t是航天器外弹道跟踪数据对应的采样时间变量;
样本子集个数:u=4,样本集合:U=Φk,一致性集合判别阈值:sigma=2×min(R),(min(R)是集合随机误差的标注差估计值的集合的最小元素);
算法迭代次数:
利用抗野值模型估计迭代算法估计第k(1≤k≤L)段航天器外弹道跟踪数据参数迭代算法如下:
对于第k(1≤k≤L)段航天器外弹道跟踪数据重复ξ次执行步骤(4.1)~(4.3),其中第c,1≤c≤ξ次执行后的结果是:
步骤4.1:从样本集U中随机抽取u个样本点,作为子集ηc;
步骤4.2:将样本子集ηc中的u个样本点依次代入具体模型中,解方程组,得到模型参数其中,(te,ye)是u个样本点中的其中一个,也即U中的一个随机样本点,e是时间角标;/>Sk依次是第k(1≤k≤L)段航天器外弹道跟踪数据抗野值估计模型第c次执行后对应的三次项系数,二次项系数,一次项系数,常数项系数和系统误差;
步骤4.3:计算出本次模型参数对应的支持度函数/>
对于第k(1≤k≤L)段航天器外弹道跟踪数据第c次执行后的支持度函数有如下定义:步骤4.2中第c次执行后得到模型参数/>时,支持度函数/>就是样本集U中所有l个样本到/>的距离小于一致性集合判断阈值sigma的样本点的个数,如式(15)所示;
其中函数ρ(·)被定义为积累函数,其表达形式如式(16)所示;
其中,Uc表示样本集合U中的第c个样本点,因此,第k(1≤k≤L)段航天器外弹道跟踪数据经过ξ次执行后,一共可以得到ξ个模型参数,即,它们一一对应ξ个支持度函数/>
步骤4.4,假定中的最大值为/>则将与/>对应的模型参数作为利用抗野值模型估计迭代算法估计得到的第k(1≤k≤L)段航天器外弹道跟踪数据参数的估计值;
因此,第k(1≤k≤L)段航天器外弹道跟踪数据的估计值为式(17)所示;
步骤5,识别和修正航天器各段外弹道跟踪数据中的斑点型野值;
具体为:将第k(1<k≤L)段外弹道跟踪测量数据(t(k-1)×l+1,y(k-1)×l+1),(t(k-1)×l+2,y(k-1)×l+2),…,(tk×l,yk×l)代入公式(5)可得公式(18);
将式(17)和(18)作差并取绝对值,如式(19)所示;
定义:因此式(19)可化为式(20);
当时,将yz识别为野值点,并用/>修正它,否则将其识别为正常点;
步骤6,将识别和修正野值后的各段航天器外弹道跟踪数据进行拼接;
设识别和修正野值后的第k(1≤k≤L)段航天器外弹道跟踪数据向量如式(21)所示;
对应采样时刻向量Tk=[t(k-1)×l+1,t(k-1)×l+2,…,tk×l];则完整弹道可以被表示为对应的采样时刻[T1,T2,…,Tk,…,TL]。
为了进一步验证本发明的算法对斑点型野值的识别效果,在给出仿真数据的基础上将本发明方法和多个经典方法进行对比分析。
将含有斑点型野值的航天器外弹道跟踪数据时间序列数学模型具体化为式(22)和式(23)得到含有斑点型野值的弹道跟踪测量仿真数据。ti(i=1,2,3)是时间变量,采样频率20个点/秒,时间变量t2∈[5,6)∪[16,17.5)对应两个斑点型野值,时间变量t1∈[0,5)∪[6,16)∪(17.5,20]对应正常数据点,时间变量t3=t1∪t2对应含有两类斑点型野值的弹道跟踪测量仿真数据;y代表弹道跟踪测量(位移,角度等),ω是均值为0,标准差为1的高斯白噪声,δ是野值,其绝对值远远大于标准差1,|δ|》1。
使用三个传统的野值识别算法(趋势外推算法、一阶差分趋势外推拟合算法和中值滤波差分算法)和本申请算法对由式(23)的数据模型产生的含有两类斑点型野值的弹道跟踪数据进行拟合(识别和修正)得到图3四种斑点型野值识别算法拟合曲线对比图。其中标准对比数据曲线由式(22)的数据模型产生,用来定量对比分析四种算法对斑点型野值的识别和修正效果。将图3中四种斑点型野值拟合曲线分别和标准对比曲线作差求绝对值得到图4四种斑点型野值识别算法误差结果对比图。
观察图3和图4可以发现趋势外推算法只是对斑点型野值起到了一定的抑制作用,根部无法完成对斑点型野值的识别和修正;由于中值滤波差分算法和一阶差分外推拟合算法对差分后的数据进行了修正,导致积分还原后的弹道数据出现了整体偏移的现象,随着斑点型野值的个数增多,数据整体偏移的程度加大;在两个斑点型野值区域(5-6秒和17-17.5秒之间),与一阶差分外推算法、中值滤波差分算法和趋势外推算法拟合(识别和修正)曲线相比,本申请所提算法和标准对比曲线的误差绝对值最小,对斑点型野值的识别和修正效果最优。
Claims (4)
1.一种航天器外弹道跟踪数据斑点型野值识别和修正方法,其特征在于,具体按照以下步骤实施:
步骤1,建立含有斑点型野值航天器外弹道跟踪数据时间序列数学模型;具体为:
在时间T1,T2,…,Ti,…,Tn上采样得到一组航天器外弹道跟踪数据Y1,Y2,…,Yi,…,Yn,即,时刻Ti(i=1…n)的航天器外弹道跟踪采样数据为Yi(i=1…n),时刻Ti(i=1…n)航天器外弹道跟踪数据Yi(i=1…n)的数学模型函数,如式(1)所示:
Yi=Xi+Si+ωi+δi (1);
式(1)中,Si(i=1…n)为系统误差;ωi(i=1…n)为随机误差项,是标准差为σ的高斯白噪声,{εi}~N(0,σ),δi表示野值,当某个采样时刻出现野值时,δi的绝对值远大于高斯白噪声的标准差σ,满足σ<<|δi|,否则δi=0,Xi为外弹道跟踪测量数据时刻Ti(i=1…n)的真值,被表示为三次多项式,如式(2)所示;
Xi=θ3(Ti)3+θ2(Ti)2+θ1(Ti)1+θ0 (2);
式(2)中,θ0θ1θ2θ3依次为三次多项式Xi的常数项、一次项系数、二次项系数和三次项系数;
步骤2,对航天器外弹道跟踪数据进行分割;具体为:
假定在时间t1,t2,…,ti,…,tn上采样得到一组航天器外弹道跟踪数据y1,y2,…,yi,…,yn,即,时刻ti(i=1…n)的航天器外弹道跟踪采样数据为yi(i=1…n);将y1,y2,…,yi,…,yn分割成L段,且满足分割后每段航天器外弹道跟踪数据序列的数据长度为l的整数,则第k(1≤k≤L)段航天器外弹道跟踪数据向量mk如式(3)所示;
mk=[y(k-1)×l+1,y(k-1)×l+2,…,yk×l],1≤k≤L (3);
进一步,第k(1≤k≤L)段航天器外弹道跟踪数据mk中第j个测量数据如式(4)所示;
步骤3,求取各段航天器外弹道跟踪数据的标准差;具体为:
第k(1≤k≤L)段航天器外弹道跟踪数据集合mk的第j个测量数据表示为式(1)的形式,令i=z,Yi=yz,Ti=tz,则式(1)变形为式(5);
其中,z=((k-1)×l+j),(1≤j≤l),tz代表第z个采样时刻,依次为第k段航天器外弹道跟踪数据mk对应的三次多项式的常数项、一次项系数、二次项系数和三次项系数,/>依次是第k段航天器外弹道跟踪数据mk在时刻z的系统误差和随机误差;
第k(1≤k≤L)段航天器外弹道跟踪数据mk中一共有l个数据点y(k-1)×l-1,y(k-1)×l-2,…,yk×l,将式(5)全部化为方程组的形式,如式(6)所示;
将式(6)化为矩阵形式,如式(7)所示;
Mk=Tkθk+ωk (7);
式(7)中的各个矩阵和向量展开如式(8)所示;
式(7)中,Mk是第k段航天器外弹道跟踪数据mk的向量形式,Tk是第k段航天器外弹道跟踪数据对应的时间矩阵,θk是第k段航天器外弹道跟踪数据模型参数向量,ωk是第k段航天器外弹道跟踪数据模型随机误差向量形式;
定义式(9)为式(7)第k(1≤k≤L)段航天器外弹道跟踪数据最小二乘拟合的损失函数J(θk);
则式(9)中损失函数J(θk)的偏导数如式(10)所示,并令其为0向量;
求解式(10)得到第k段航天器外弹道跟踪数据模型参数向量θk的估计值如式(11)所示;
则第k段外弹道跟踪测量数据mk的随机误差的标准差的估计值,如式(12)所示;
则所有L段外弹道跟踪测量数据集合{m1,m2,…mk…,mL-1,mL}对应的随机误差的标准差估计值集合为
步骤4,对各段航天器外弹道跟踪数据模型的参数进行估计;
步骤5,识别和修正航天器各段外弹道跟踪数据中的斑点型野值;
步骤6,将识别和修正野值后的各段航天器外弹道跟踪数据进行拼接。
2.根据权利要求1所述的一种航天器外弹道跟踪数据斑点型野值识别和修正方法,其特征在于,所述步骤4中,其具体步骤是:
已知第k(1≤k≤L)段共l个航天器外弹道跟踪数据采样点集合,如式(13)所示;
Φk={(t(k-1)×l+1,y(k-l)×l+1),(t(k-l)×l+2,y(k-l)×l+2),…,(tk×l,yk×l)},1≤k≤L (13);
建立抗野值估计模型,如式(14)所示;
φ(Θ)=Θ3(t)3+Θ2(t)2+Θ1t+Θ0+S (14);
其中,Θ3Θ2Θ1Θ0S依次是航天器外弹道跟踪数据抗野值估计模型的三次项系数,二次项系数,一次项系数,常数项系数和系统误差,t是航天器外弹道跟踪数据对应的采样时间变量;
样本子集个数:u=4,样本集合:U=Φk,一致性集合判别阈值:sigma=2×min(R),min(R)是集合随机误差的标注差估计值的集合的最小元素;
算法迭代次数:
步骤4.1:从样本集U中随机抽取u个样本点,作为子集ηc;
步骤4.2:将样本子集ηc中的u个样本点依次代入具体模型中,解方程组,得到模型参数其中,(te,ye)是u个样本点中的其中一个,也即U中的一个随机样本点,e是时间角标;
步骤4.3:计算出本次模型参数对应的支持度函数/>
对于第k(1≤k≤L)段航天器外弹道跟踪数据第c次执行后的支持度函数有如下定义:步骤4.2中第c次执行后得到模型参数/>时,支持度函数/>就是样本集U中所有l个样本到/>的距离小于一致性集合判断阈值sigma的样本点的个数,如式(15)所示;
其中,函数ρ(·)定义为积累函数,其表达形式如式(16)所示;
其中,Uc表示样本集合U中的第c个样本点,因此,第k(1≤k≤L)段航天器外弹道跟踪数据经过ξ次执行后,一共可以得到ξ个模型参数,即,它们一一对应ξ个支持度函数/>
步骤4.4,设中的最大值为/>则将与/>对应的模型参数作为利用抗野值模型估计迭代算法估计得到的第k(1≤k≤L)段航天器外弹道跟踪数据参数的估计值;
因此,第k(1≤k≤L)段航天器外弹道跟踪数据的估计值为式(17)所示;
3.根据权利要求2所述的一种航天器外弹道跟踪数据斑点型野值识别和修正方法,其特征在于,所述步骤5中,具体为:
将第k(1<k≤L)段外弹道跟踪测量数据(t(k-1)×l+1,y(k-1)×l+1),(t(k-1)×l+2,y(k-1)×l+2),…,(tk×l,yk×l)代入公式(5)可得公式(18);
将式(17)和(18)作差并取绝对值,如式(19)所示;
定义:因此式(19)可化为式(20);
当时,将yz识别为野值点,并用/>修正它,否则将其识别为正常点。
4.根据权利要求3所述的一种航天器外弹道跟踪数据斑点型野值识别和修正方法,其特征在于,所述步骤6中,具体为:设识别和修正野值后的第k(1≤k≤L)段航天器外弹道跟踪数据向量如式(21)所示;
对应采样时刻向量Tk=[t(k-1)×l+1,t(k-1)×l+2,…,tk×l];则完整弹道被表示为对应的采样时刻[T1,T2,…,Tk,…,TL]。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010664267.6A CN111985077B (zh) | 2020-07-10 | 2020-07-10 | 一种航天器外弹道跟踪数据斑点型野值识别和修正方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010664267.6A CN111985077B (zh) | 2020-07-10 | 2020-07-10 | 一种航天器外弹道跟踪数据斑点型野值识别和修正方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111985077A CN111985077A (zh) | 2020-11-24 |
CN111985077B true CN111985077B (zh) | 2024-03-22 |
Family
ID=73439029
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010664267.6A Active CN111985077B (zh) | 2020-07-10 | 2020-07-10 | 一种航天器外弹道跟踪数据斑点型野值识别和修正方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111985077B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112540974B (zh) * | 2020-12-24 | 2021-09-24 | 哈尔滨工业大学 | 一种基于二阶动量的航天器遥测数据野值点去除方法 |
CN113326620B (zh) * | 2021-06-03 | 2022-04-08 | 中国人民解放军32039部队 | 测距数据处理方法及装置 |
CN113589239B (zh) * | 2021-06-30 | 2023-05-16 | 中国西安卫星测控中心 | 一种雷达测量数据精度容错估计方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102322861A (zh) * | 2011-05-31 | 2012-01-18 | 电子科技大学 | 一种航迹融合方法 |
CN108304612A (zh) * | 2017-12-26 | 2018-07-20 | 南京邮电大学 | 基于噪声补偿的迭代平方根ckf的汽车雷达目标跟踪方法 |
CN108681331A (zh) * | 2018-05-21 | 2018-10-19 | 济南大学 | 一种近空间飞行器的姿态跟踪控制方法 |
CN109459705A (zh) * | 2018-10-24 | 2019-03-12 | 江苏理工学院 | 一种抗野值鲁棒无迹卡尔曼滤波的动力电池soc估计方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105814499B (zh) * | 2013-06-14 | 2021-04-02 | 华莱士·E·拉里莫尔 | 用于监测和控制具有可变结构或可变运行条件的动态机器的动态模型认证方法和系统 |
-
2020
- 2020-07-10 CN CN202010664267.6A patent/CN111985077B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102322861A (zh) * | 2011-05-31 | 2012-01-18 | 电子科技大学 | 一种航迹融合方法 |
CN108304612A (zh) * | 2017-12-26 | 2018-07-20 | 南京邮电大学 | 基于噪声补偿的迭代平方根ckf的汽车雷达目标跟踪方法 |
CN108681331A (zh) * | 2018-05-21 | 2018-10-19 | 济南大学 | 一种近空间飞行器的姿态跟踪控制方法 |
CN109459705A (zh) * | 2018-10-24 | 2019-03-12 | 江苏理工学院 | 一种抗野值鲁棒无迹卡尔曼滤波的动力电池soc估计方法 |
Non-Patent Citations (1)
Title |
---|
弹道跟踪数据野值剔除方法性能分析;侯博文;王炯琦;周萱影;李冬;何章鸣;;上海航天;20180825(第04期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN111985077A (zh) | 2020-11-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111985077B (zh) | 一种航天器外弹道跟踪数据斑点型野值识别和修正方法 | |
CN108829928B (zh) | 一种涡轴发动机自适应部件级仿真模型构建方法 | |
US6539783B1 (en) | Methods and apparatus for estimating engine health | |
EP1926001B1 (en) | Reducing gas turbine performance tracking estimation non-repeatability | |
US7788014B2 (en) | Process and methodology for root cause identification in gas turbine engine performance tracking | |
CN108960334B (zh) | 一种多传感器数据加权融合方法 | |
CN108205310B (zh) | 一种基于elm滤波算法的航空发动机包线内气路故障识别方法 | |
CN108573116B (zh) | 一种基于长短时记忆网络的航空发动机过渡态推力估计方法 | |
CN110659722A (zh) | 基于AdaBoost-CBP神经网络的电动汽车锂离子电池健康状态估算方法 | |
US7177785B2 (en) | Systems and methods for improved aircraft performance predictions | |
CN109800449B (zh) | 一种基于神经网络的航空发动机压缩部件特性修正方法 | |
CN110398942B (zh) | 一种用于工业生产过程控制的参数辨识方法 | |
CN113283511A (zh) | 一种基于权重预分配的多源信息融合方法 | |
CN107357176B (zh) | 一种航空发动机试车数据建模方法 | |
CN111144018B (zh) | 一种基于航后数据的航空发动机整机剩余性能提取方法 | |
CN116482579A (zh) | 一种稳健的噪声量自修正的电连接器剩余寿命预测方法 | |
Yoon et al. | New modeling algorithm for improving accuracy of weapon launch acceptability region | |
CN106934124B (zh) | 一种基于量测变化检测的自适应变划窗方法 | |
CN116468174A (zh) | 飞行参数预测及置信度评价方法 | |
CN110705132A (zh) | 一种基于多源异质数据的制导控制系统性能融合评估方法 | |
CN112069592B (zh) | 一种航天器外弹道跟踪测速数据特征点识别方法 | |
Raikar et al. | Denoising signals used in gas turbine diagnostics with ant colony optimized weighted recursive median filters | |
Sundaram et al. | Aircraft engine health monitoring using density modelling and extreme value statistics | |
WO2005050337A2 (en) | Intelligent database for performance predictions | |
CN114235005B (zh) | 一种适用于长期加电下快速更新六项陀螺仪诸元方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |