宇航学报
JOURNAL OF ASTRONAUTICS
1999年 第2期　No.2 1999



靶场外测数据野值点的统计诊断技术
胡绍林　孙国基
　　摘　要　本文修正了野值点的数学描述性定义,建立了野值的分类模型,创立了多个野值点及野值斑点的“稳健-似然比”检验技术,并进行了仿真计算。
　　主题词　统计诊断　野值　建模
STATISTICAL DIAGNOSIS METHOD FOR 
OUTLIERS FROM SPACECRAFT TRACKING DATA
Hu Shaolin Sun Guoji
(System Engineering Institute,Xi'an Jiaotong University.Xi'an.710049)
　　Abstract　A modified definition and some classified models are built for outliers,a high efficient “robust-LR (LR:likelihood ratio) test”method is established to detect and to diagnose multi-outliers or patch-types outliers from measurement data set,and a simulation is done to show the availability of this new diagnosis approach.
　　Key words　Statistical diagnosis Outliers Modeling
1　引　言
　　众所周知,无论是靶场测量数据高精度处理中广泛采用的多项式滤波技术、平滑微分技术、功率谱分析技术,还是基于Kalman滤波理论的实时预报和控制技术,从应用统计学的角度都可以看作为:在一定的公式或模型假定下,依据对动态目标的一组跟踪测量数据所形成的数据集合D,采用统计的或信号处理的方法分析、控制或预报目标运动状态及特性。由此可见,数据集合D是进行统计推断或工程决策的起点,集合D中各数据点的质量如何,直接影响着我们做出判决的精确度和可靠性。
　　多年来靶场外测数据处理工作的实践告诉我们,即使是高精度的测量设备,由于多种偶然因素的综合影响或作用,采样数据集合中往往包含1～2%有时甚至多达10～20%(例如,雷达进行高仰角跟踪)的数据点严重偏离目标真值。工程数据处理领域称这部分异常数据为野值。
　　野值对高精度数据处理工作有着十分不利的影响。近几十年来,国际统计界大量的研究结果表明,无论是基于最小二乘理论的模型参数最优估计、最优线性滤波算法(包括多项式滤波与递推Kalman滤波等)还是频谱分析中著名的极大熵谱估计,对采样数据中包含的野值点反应都极为敏感［1,2,5,7,8］。工程数据处理的实践也证实,在航天靶场试验的数据事后处理工作中,采样数据合理性检验是数据预处理工作的重要一环。它对改进处理结果的精度、提高处理质量有着相当的重要性。
　　工程数据野值点的识别方法和处理技术是一个涉及到多学科知识综合运用的应用性研究课题。对于航天工程中各类长弧段测量数据集合,要准确识别出所含有的全部野值点是一件十分不容易的事。特别是随着工程对象复杂化和高速计算机迅速普及,人们碰到的数据集合越来越大、数据维数越来载高,检测采样数据质量和诊断野值点的工作变得越来愈困难。因此,建立高效、可靠的野值点检测与诊断方法,其必要性日益为工程数据处理工作者所重视。
　　本文,从野值点的模型分析入手,构造一套检测和诊断野值点的统计算法。
2　野值点的建模
　　什么是野值?在探索性数据分析和数据处理领域中,野值又称异常值(Outliers)。其定义一直不很明确,持不同态度的应用统计学家对其定义也不一样［5］。作者在文［5］中对几种常用的定义进行过较为详尽的剖析和比较。
　　本文建议采用野值点的定义为:集合D中严重偏离大部分数据所呈现趋势的小部分数据点。这一定义强调主体数据所呈现的“趋势”,以偏离数据集合主体的变化趋势为判别异常数据依据,并明确指出野值在集合D中只占小部分(即最多不超过一半),这从直观上是合理的。本定义不但可以覆盖简单随机抽样情形,包括Edeworth(1887)的定义和文［6,7］的定义,而且也覆盖了时间序列及随机系统中的野值数据情形,具有较为广泛的适应性。
2.1　野值点的成因分析
　　采样数据集合中出现野值点的原因很多。就航天靶场试验的外测数据而言,产生野值主要有如下几个方面的原因:
　　(1)操作和记录时的过失,以及数据复制和计算处理时所出现的过失性错误。由此产生的误差称为过失误差(Gross Error)。
　　(2)采样环境的变化。取样母体的突然改变使得部分数据与原先样本的模型不符合,例如大地测量时地震或其它突变现象的出现、雷达跟踪时应答机工作状态的不稳定、靶场试验时飞行目标的变轨与各级关机点附近的取样,等等。
　　(3)实际采样数据中也可能出现另一类异常数据,它既不是来自操作和处理的过失也不是由突发性强影响因素导致的,而是某些服从长尾分布(long-tailed distribution)的随机变量(例如,服从t分布的随机变量)作用的结果。
2.2　野值点的分类
　　工程采样数据中出现的野值点,比较常见的有如下几种类型:
　　(1)离群点　又称跳点,是统计诊断领域中研究得最多也比较成熟的一类［7］。简单随机抽样情形下,野值点主要表现为离群点。
　　(2)孤立型野值　它的基本特点是,某一采样时刻ti处的测量数据是否为野值与ti-1及ti+1时刻数据的质量无必然联系。而且,比较常见的是当ti时刻的测量数据呈现异常时,在ti时刻的一个邻域内(如等时刻)采样数据质量是好的,即野值点的出现是孤立的。母体不相依情形下,动态测量数据中孤立(isolated)异常值的出现也是比较普遍的情形之一。
　　(3)野值斑点　简称斑点(Patchy),是指成片出现的异常数据。它的基本特征是,在ti-p时刻出现的野值,也可能带动y(ti-p+1),…,y(ti)均严重偏离真值。雷达跟踪高仰角目标的测量数据序列中,野值斑点的出现是比较常见的故障现象。
2.3　野值点的模型描述
　　假定某测站对空间飞行目标进行跟踪观测,记测量数据集合为
　　　　　　D={y(t1),…y(tN)}　　　　　　　　　　　　　　(1)
　　(1)测量数据的“加性(additive)”分解:由于目标运动的连续性,我们可以假定测量对象y　(t)可合理地分解成四部分:
　　　　　　y(t)=ytr(t)+yp(t)+εs(t)+ε0(t)　　　　　　　(2)
其中：①ytr(t)称为趋势分量,描述的是目标的趋性运动;
　　　②yp(t)称为周期分量,描述的是目标周期性变化,即目标在趋性“进动”的同时,受周期性变化因素作用,每隔一段时间后运动状况出现一定的相似性;
　　　③εs(t)是测量数据随机误差分量,反映的是测量设备状态随机误差和目标飞行环境中多种不确定性扰动共同作用的结果。
　　　④ε0(t)称为过程的污染分量或突变性分量,其作用结果是使测量数据发生严重偏离,它的量级要么为零,要么显著大于随机误差标准差,即
　　　　　　　　　　　　　　　　　　　(3)
式中,γ为异常数据出现的几率。大量数据分析表明,γ通常介于1%～2.5%之间。
　　(2)加性异常数据(Additive Outliers,AO)模型:在某些采样时刻ti0,由于测量环境突然发生的强影响因素影响,以及记录与操作的失误,使得测量数据严重偏离大部分数据所呈现的趋势,用模型描述为:
　　　　　　ε0(t)=λδ(i,i0)　　　　　　　　　　　　　　(4)
式中,λ为测量数据异常幅度,δ(i,i0)为Kronecker delta函数。
　　(3)新息异常数据(Innovation Outliers,IO)模型［4］:对于观测值模型(2),当随机过程ε0　(t)服从长尾分布且具有广义平稳性时,利用平稳随机过程及现代时间序列分析理论,可进一步假定ε0(t)服从ARMA(p,q)模型。
3　野值点的统计诊断
　　对于简单随机抽样情况下野值点的识别,文献［1,6,7］及其所附参考资料中都有详细的讨论;孤立野值的诊断问题,也已经有大量的研究成果可供采用［4］;对于多个野值点(特别是斑点型野值)的诊断与处理,因涉及到误差的前后相依性,见诸文献的成功的处理方法或程序尚不多见。
　　考虑到斑点型野值点统计诊断的技术难度和工程实用价值,本文将对此进行一些初步的探索和研究。本文的基本想法是,将长弧段测量数据集合中野值点的诊断和处理工作分成三个阶段:
　　①容错处理:构造一条稳健的“柱心”轨线;
　　②可疑点的筛选:通过实测数据与该“柱心”轨线的比较,初步确定出可能包含野值点个数s及其所在位置(ti1,…,tis);
　　③统计诊断:构造检验统计量和检验门限,进一步核实②中所确定的s个测量数据确实为异常数据。
　　下面,从工程应用的角度建立起实现这一方法的基本步骤及具体算法。
3.1　“柱心”的稳健多项式拟合
　　根据飞行器运动学和轨道动力学理论,可以假定:在有限长度的被观测弧段上,目标的运动轨迹y(t)是几乎处处连续、按段光滑的。由Weierstrass逼近定理,测量对象y(t)可以被代数多项式与三角多项式的线性组合一致逼近,即
　　　　　　(5)
其中,T为周期性变化分量ype　(t)的演变周期。
　　基于式(1)给出的测量数据集合D={y(ti):i=1,2,…,N},拟合模型(5)可以改写为线性回归模型的形式:
　　　　　　　　　　　　　　　　　　(6)
式中，1)　
　　　2)　β=(a0,a1,…,a1,b1,…bm,c1,…,cn)τ,ε(t)=εs(t)+ε0(t)
　　构造线性拟合模型参数向量β及随机误差根方差σ的M-型稳健估计方程:
　　　　　　　　　　　　　　　　　　　　(7)
其中，1)　
　　　2)　　　(门限c1和c2为适当选取的非负常数)
　　显然,当选取门限c1=+∞时,ρ(ω)=0.5ω2和φ(ω)=ω均为无界函数,式(7)为最小二乘估计方程;当适当选取门限c1和c2为有限非负常数时,ρ(ω)和φ(ω)均为有界函数,此时式(7)确定的参数估计(,)为稳健-抗扰性估计,具有良好的抗野值能力。
　　下面,给出模型参数M-型估计方程(7)的一组迭代求解算法:
　　　　(8)
并置初值为:
　　　　　　　　　　　　　　　　(9)
　　利用上述的稳健-抗扰估计值,不难给出测量对象y　(t)的具有抗野值干扰能力的稳健拟合曲线:
　　　　　　　　　　　　　(10)
注记:在工程数据预处理时,为了简化计算,也可适当采用如下形式的Huber型或Hampel型φ函数:
　　(1)　　　　　　　　　　　　　　　　　(11)
　　(2)　　　　　　　　　　　　　　　　(12)
式中,c为根据数据和经验适当选取的非负常数。
3.2　野值点位置的确定
　　计算跟踪测量数据集合D中各数据点的规范化拟合残差:
　　　　　　　　　　　　　(13)
绘制拟合残差集合散点图,并记录所有使的时间点及数据个数s。利用“3σ准则”可初步认为这s个点最有可能是野值点。由此得出可疑为野值的数据点个数s,及其所在的位置{ti1,…,tis}。
3.3　野值点的似然比检验
　　假定模型(2)的随机误差分量为独立同分布(iid)的Gauss白噪声。当测量数据集合D中不含野值时,可以证明［7］:关于可疑点{y(ti1),…,y(tis)}是否为野值点的似然比统计量
　　　　　　　　　　　　　　　　(14)
式中,Rss和Rss分别为删除可疑点{y(ti1),…,y(tis)}前/后最小二乘拟合的残差平方和。
　　由式(14)可以看出,似然比函数LR(s)是的严格增函数,并且在上述正态性假定下可以证明［7］统计量服从参数为的β-分布。因此,构造检验统计量:
　　　　　　　　　　　　　　　　　　　(15)
为数据合理性检验的统计检测指标。
　　简记上述β-分布的α下分位点为Cα,即(α=0.05或0.025)。当统计量LR(s)＞1-Cα时,可以100×(1-α)%的概率认为集合D中可疑点{y(ti1),…,y(tis)}是真正的故障数据。
4　仿真计算
　　图1是计算机仿真生成的某地基雷达跟踪空间目标“测量”的20秒数据,采样频率为每秒10点。图1中第50、91～110和第200点为采用外部叠加方法形成的22个野值点,仿真数学模型为带随机分量的4阶代数多项式,模型设定参数见(表1);孤立野值点偏置量:ε0(50)=-10,ε0(150)=10;斑点型野值的数学模型为:ε0(i+1)=1-ε0(i),ε0(90)=-5(i=90,…,109)。
　　基于上述仿真数据,下表1中同时给出了模型系数的(全点)最小二乘估计、按式(8)～(9)计算的M-型稳健估计和(剔点)最小二乘估计,以及各种拟合的残差平方和。
　　　
图1　含22个野值点仿真数据散点图　　　图2　基于M-型估计的拟合曲线　
　　　　表1　仿真模型参数的设定值及统计估计值

　a0a1a2a3a4σRss
仿真设定参数-5.016.0-0.5-0.030.0010.1　
全点最小二乘-4.85715.960-0.5157-0.027630.0009241.83099653.746
M-型估计-4.93015.984-0.4998-0.029910.0009960.101392.00452
剔点最小二乘-4.94115.993-0.5017-0.029780.0009930.099741.71095

　　分别按M-型估计和全点最小二乘估计计算不同拟合的标准化残差,如图3～图4所示:
　　　　
图3　基于M-型估计的标准化拟合残差　　　图4　基于全点最小二乘估计的拟合残差
　　从表1和图2可以看出,M-型估计具有良好的抗野值影响能力;从图3还可以看出,拟合残差中不包含仿真模型的趋势分量,设定的22个野值点所对应的残差绝对值均远远超过3,其余残差的绝对值明显地小于3,按“3σ准则”误报可疑点的个数为0;从图4可以看出,野值点对最小二乘拟合的不利影响是十分明显的,不但残差中残余有明显的趋势分量,除22个野值点外还有11个点的残差绝对值超过3,按“3σ准则”导致误报可疑点11个。
　　仿真结果说明,本文的“稳健-似然比”检验方法对野值斑点具有较高的检测效率和可靠性。
　　致谢:本文完稿过程中得到李济生院士的指导和帮助,在此谨致谢忱。
　　作者单位：西安交通大学系统工程研究所.西安.710049
参考文献
　1　PJ.Huber FJ.Robust Statistics.Wiley,New York,1981
　2　Hampel,et.Robust Statistics-approach based on influence.John Wildy and Sons Press,1986
　3　胡绍林,范金城.周期性测量数据三角多项式滤波及预报.西安交通大学学报,1996,30(2):110-116
　4　Tiao GC.Autoregressive Moving Average Models.Inter.Problems and Outlier Detection in Time Series,Handbook of Stat.5
　5　范金城,胡绍林.动态测量数据抗扰性研究综述.数理统计与应用概率,1996,11(3):244-248
　6　Wang Xueren,et.Assessment of Local Inference in Multivariate Analysis.Technical Reports(1990-1991),Yunnal University Press,1991
　7　韦博成等.统计诊断引论.南京:东南大学出版社,1991
　8　胡峰,范金城.动态测量系统的有界影响滤波.控制理论及应用,1992,11(1):36-45
　　本文为西安交通大学研究生博士学位论文基金(DFXJU1998-10)资助
收稿日期:1997年4月21日,修回日期:1998年2月24日
