CN110161543A - 一种基于卡方检验的部分粗差抗差自适应滤波方法 - Google Patents

Abstract

本发明公开了一种基于卡方检验的部分粗差抗差自适应滤波方法，首先，基于观测模型异常检验量，分析了观测值之间的相关性，并针对由于观测值之间的相关性所导致的粗差误判问题，提出了部分粗差抗差方法；然后根据假设检验理论，构造了滤波模型整体检验量，基于卡方检验判断整体模型是否存在异常。当判定出模型存在故障时，本发明就采用部分粗差抗差自适应方法对异常的位置进行定位，并通过放大协方差，保障定位的精度和鲁棒性；最后设计了两组实验，采用三种方案进行对比分析，以验证本发明所提方法的性能。实验结果表明，该方法极大地消弱了观测值之间相关性的影响，能准确的识别粗差位置，明显降低了粗差探测的误警率，保证了定位的鲁棒性。

Description

一种基于卡方检验的部分粗差抗差自适应滤波方法

技术领域

本发明属于全球导航卫星系统粗差探测和识别技术领域，尤其涉及一种基于卡方检验的部分粗差抗差自适应滤波方法。

背景技术

GNSS多系统融合由于引入了更多的可视卫星，为高精度定位提供了更可靠的保证，但同时也应该注意到，在城市峡谷等复杂情况下，由于观测环境比较复杂，加之仪器振荡而表现出的不稳定特性，观测维数的大幅增加必然也会导致粗差出现的概率及复杂度成倍增加，给数据处理带来了极大地挑战。

目前，粗差探测和识别理论根据思想不同可以总结为两大类：一类是将粗差纳入函数模型，认为粗差观测值与正常观测值方差相同而期望不同，即数据探测法；另一类是将粗差纳入随机模型，认为粗差观测值与正常观测值期望相同而方差不同，即抗差估计理论。两者都极大的推动了测量数据处理理论的发展，相较于前者而言，抗差估计理论因其独特优势而得到更加广泛的应用。2016年，有学者就比较了上述两种粗差探测方法抵御粗差的效果差异，结果表明抗差估计理论更易实现粗差的探测和定位。

抗差估计理论最早是由丹麦学者Krarup提出，并将其引入测量界；随后Caspary对该理论进行完善，并作出了一系列研究；与此同时，国内学者周江文等人也开始对抗差估计理论展开研究，提出了IGG-I抗差估计方案；而在动态数据处理过程中，Kalman 滤波是GNSS定位与导航普遍采用的参数估计策略，有学者就基于贝叶斯推断理论，构建了抗差Kalman滤波方法，进一步保障了动态定位的精度和鲁棒性。

发明内容

发明目的：因为抗差Kalman滤波是通过构造观测模型异常检验量来进行粗差探测和定位。由于真误差不完全可知，因此粗差探测可能发生误警和漏检，特别是当系统中观测值异常检验量存在强相关时，由于观测值之间的相关性，可能会造成部分粗差被分配到其他正常观测值中，从而降低正常观测值对参数估计的贡献，导致粗差误判的概率明显升高，进而给参数估计带来显著偏移。本发明从实用价值方面出发，为避免由于观测值之间的相关性而导致粗差转移的问题，本发明提出了一种基于卡方检验的部分粗差抗差自适应滤波方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于卡方检验的部分粗差抗差自适应滤波方法，该方法包括如下步骤：

(1)使用大地测量型接收机采集卫星数据，通过对原始观测值进行处理，从而构造观测方程；

(2)根据假设检验理论，构造kalman滤波模型整体检验量，采用卡方检验判断整体模型是否存在异常，若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常，T_k指的是整体假设检验统计量，T_l指的是整体假设检验统计量阈值；

(3)当判定模型存在异常时，应首先进行观测噪声自适应，对由于检验量之间的相关性所导致的粗差误判进行部分观测值粗差抗差；

(4)当动力学模型存在异常扰动时，对系统噪声进行自适应，以消除动力学模型预报信息与动态载体运行轨迹之间的差异。

进一步的，步骤(2)的具体方法如下：

设置Kalman滤波中新息向量及其协方差矩阵表示：

V_k,k-1＝A_kX_k,k-1-L_k (1)

式中，V_k,k-1表示新息向量；A_k为观测设计矩阵；X_k,k-1为状态预报向量；L_k为观测向量；表示新息向量协方差矩阵；R_k表示观测向量协方差矩阵；Q_k,k-1为状态预报向量协方差；

当Kalman滤波整体模型未发生异常时，新息向量满足零均值的高斯分布；当异常发生时，异常不会改变观测值的分布类型，但使其概率分布产生一定的偏移，即如下假设检验问题：

式中，H₀是原假设，表示整体模型无异常；H₁是备择假设，表示整体模型存在异常；λ为概率分布偏移量，统计学上称为非中心化参数；N表示正态分布；

若已知先验单位权方差为σ²，取整体假设检验统计量为：

其中，T_k满足自由度为t的χ²分布；

根据给定的显著水平α、自由度t和分布类型可以确定阈值：

若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常。

进一步的，步骤(3)的具体方法如下：

设第i颗卫星观测模型异常检验量为：

式中，n表示观测值维数，e_i表示第i个元素为1，其它元素为0的n_×1矩阵。

同理，第j颗卫星观测模型异常检验量为：

根据误差传播定律可得两个检验量之间的相关系数为：

在观测噪声自适应过程中，由于观测值存在相关性ρ_i,j，所以对不同类型观测值进行分类，并且每次只对粗差最大的观测值进行抗差。

进一步的，步骤(4)的具体方法如下：

采用自适应滤波方法，基于状态不符值构造自适应因子动态调节状态预报向量和观测向量的权比，消除动力学模型预报信息与动态载体运行轨迹之间的差异，具体实施方法如下：

预测状态向量不符值可表示为：

ΔX＝X_k,k-X_k,k-1 (10)

式中，ΔX为状态向量不符值，X_k,k为抗差Kalman滤波解，X_k,k-1为状态预报向量；

在自适应滤波过程中，状态预测向量X_k,k-1的协方差矩阵应该等于实际预测向量的偏离程度，即：

Q_ΔX＝Q_k,k-1 (11)

式中，Q_ΔX为ΔX的方差-协方差向量，Q_k,k-1为状态预报向量的协方差；

取：

Q_ΔX＝ΔXΔX^T (12)

由式(11)和式(12)即可求解自适应因子α，即：

式(13)即为理论自适应因子的估计值，通过使得Q_k,k-1＝αQ_ΔX，旨在保证滤波器输出噪声的不确定度与理论噪声不确定度相当，以消除动力学模型预报信息与动态载体运行轨迹之间的差异。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

采用本发明所提出的技术手段，通过每次只对不同观测类型粗差最大的观测值进行抗差，可以有效避免由于粗差的转移而对正常观测值产生误判的问题，可以明显消弱观测值之间相关性的影响，准确的识别粗差的位置，从而保证定位的精度和鲁棒性。

附图说明

图1GPS伪距观测值之间的相关性；

图2GPS伪距+多普勒观测值之间的相关性；

图3基于卡方检验的部分粗差抗差自适应滤波方法框架；

图4伪距实验N方向定位偏差；

图5伪距实验E方向定位偏差；

图6伪距实验U方向偏差；

图7伪距实验卫星数；

图8伪距+多普勒实验N、E、U方向定位偏差；

图9伪距+多普勒实验N、E、U方向测速偏差。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本文基于大地测量型接收机采集的数据，通过对原始观测值进行处理,包括卫星位置计算，大气延迟处理等，提出了一种基于卡方检验的部分粗差抗差自适应滤波方法，该方法包括如下步骤：

(1)根据假设检验理论，构造kalman滤波模型整体检验量，采用卡方检验判断整体模型是否存在异常，若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常，T_k指的是整体假设检验统计量，T_l指的是整体假设检验统计量阈值；

(2)当判定模型存在异常时，应首先进行观测噪声自适应，因此，本发明分析了观测模型异常检验量之间的相关性，针对由于检验量之间的相关性所导致的粗差误判问题，提出了部分粗差抗差方法；

(3)当动力学模型存在异常扰动时，此时仅依靠观测噪声自适应过程仍无法保证最终定位结果的鲁棒性。因此，就需对系统噪声进行自适应，以消除动力学模型预报信息与动态载体运行轨迹之间的差异。

步骤(1)，整体检验量构造方法如下：

Kalman滤波中新息向量及其协方差矩阵表示：

V_k,k-1＝A_kX_k,k-1-L_k (1)

式中，V_k,k-1表示新息向量；A_k为观测设计矩阵；X_k,k-1为状态预报向量；L_k为观测向量；表示新息向量协方差矩阵；R_k表示观测向量协方差矩阵；Q_k,k-1为状态预报向量协方差。

当Kalman滤波整体模型未发生异常时，新息向量满足零均值的高斯分布；当异常发生时，异常不会改变观测值的分布类型，但使其概率分布产生一定的偏移，即如下假设检验问题：

式中，H₀是原假设，表示整体模型无异常；H₁是备择假设，表示整体模型存在异常；λ为概率分布偏移量，统计学上称为非中心化参数；N表示正态分布。上述的假设检验问题是卡方检验的基础。

从式(1)可以看出，新息向量不仅包含了当前历元的观测模型信息，还包含了动力学模型信息，因此可以将公式(1)作为判断Kalman滤波整体模型是否存在异常的一个重要指标。若已知先验单位权方差为σ²，本实施例设为1，取整体假设检验统计量为：

其中，T_k满足自由度为t的χ²分布。

根据给定的显著水平α、自由度t和分布类型可以确定阈值，本实施例α设置为0.001、其值：

整体检验的目的是探测Kalman滤波的整体模型是否存在异常，若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常。

步骤(2)，部分粗差抗差算法，公式推导如下：

当判定模型存在异常时，应首先进行观测噪声自适应。而在传统的观测噪声自适应抗差迭代过程中，所采用的方法是对所有的观测值进行抗差，这并未考虑到观测值之间的相关性，可能会造成部分粗差被分配到其他正常观测值中，从而降低正常观测值对参数估计的贡献。因此该部分内容重点分析了观测值之间的相关性，主要是基于原始观测值构造观测模型异常检验量进行分析，旨在削弱观测值之间相关性的影响，提出一种部分粗差抗差算法。需要说明的是，观测值之间的相关性分析是部分粗差抗差算法的理论基础，而并非是观测噪声自适应过程中的所必须的步骤，也就是说，观测噪声自适应过程仅是利用观测模型异常检验量来识别粗差，识别粗差后，通过放大观测协方差矩阵来实时动态消除粗差的影响，并不涉及观测值相关性分析。

观测模型异常检验量是在观测值的基础上构造的，主要用于判断观测值的粗差大小。第i颗卫星观测模型异常检验量为：

式中，n表示观测值维数，e_i表示第i个元素为1，其它元素为0的n_×1矩阵。

同理，第j颗卫星观测模型异常检验量为：

根据误差传播定律可得两个检验量之间的相关系数为：

公式9的推导主要是用于观测值之间相关性的分析，旨在提出部分粗差抗差算法。相关系数表征两个变量之间的相关程度，相关系数绝对值越大，相关性越强，反之，相关系数越接近于0，相关度越弱。相关强度可以通过表1进行判断。

表1相关系数表

考虑到非差观测值形式简单且不受参考星异常观测的影响，下面基于非差伪距原始观测值以及非差伪距+多普勒原始观测值两种观测模型利用公式6-9对观测值之间的相关性进行分析。选取南京市CORS站两组数据进行实验分析，具体实验信息如表2所示。

表2实验信息表

表3给出了某一历元GPS卫星伪距观测值之间的相关系数值，由于呈现对称关系，故只列出上三角部分，图1更为形象的展示了表3中的信息，忽略了相关系数的正负，同时将相同观测值之间的相关系数置为0。图2给出了某历元GPS卫星伪距和多普勒观测值之间的相关性，1-10伪距观测值，11-20为多普勒观测值，由于观测值维数较大，不再列出相应的表格。

表3 GPS伪距观测值之间的相关系数表

通过表3可以很明显的发现G10和G18卫星之间的相关系数达到0.75，根据表1 可知，两者呈现强相关性，这就意味着如果利用公式6构造的G10卫星观测模型异常检验量存在10m的粗差，那么就会对应的引起G18卫星观测模型异常检验量产生7.5m的偏差，在观测值抗差过程中，会认为G18卫星存在粗差，从而进行相应的降权或淘汰，特别是在存在多粗差的情况下，所带来的影响更为严重。因此，本步骤的研究重点就是如何有效的减弱观测值相关性的影响，从而避免粗差的转移。

从图2中可以看出，不同类型观测值之间的相关系数接近于0，相关性较弱，基本可以忽略，即伪距粗差不会转移到多普勒观测值上，而相同类型观测值之间存在一定的相关性。

因此，在观测噪声自适应过程中，应该对不同类型观测值进行分类，并且每次只对粗差最大的观测值进行抗差，避免因为粗差的转移而对正常观测值产生误判。以伪距、多普勒观测值为例，在抗差循环过程中，应该每次只对粗差最大的伪距观测值和多普勒观测值进行抗差，最大粗差的选取是根据公式6分别计算伪距和载波观测值的观测模型异常检验量，然后进行排序得到，而非对所有探测到的粗差都进行抗差，这也是本发明所提到的部分粗差抗差方法。

步骤(3)，系统噪声自适应方法如下：

Kalman滤波过程中，当动力学模型存在异常扰动时，此时仅依靠观测噪声自适应过程仍无法保证最终定位结果的鲁棒性，这主要是因为系统噪声设置的不合理导致的，即系统噪声设置的过小，相当于对状态预测向量施加了一个紧约束，从而影响了最终的定位结果。为解决系统噪声设置的不合理，本发明采用自适应滤波方法，基于状态不符值构造自适应因子动态调节状态预报向量和观测向量的权比，消除动力学模型预报信息与动态载体运行轨迹之间的差异，提高定位的鲁棒性，具体实施方法如下：

预测状态向量不符值可表示为：

ΔX＝X_k,k-X_k,k-1 (10)

式中，ΔX为状态向量不符值，X_k,k为抗差Kalman滤波解，X_k,k-1为状态预报向量。

在自适应滤波过程中，状态预测向量X_k,k-1的协方差矩阵应该等于实际预测向量的偏离程度，即：

Q_ΔX＝Q_k,k-1 (11)

式中，Q_ΔX为ΔX的方差-协方差向量，Q_k,k-1为状态预报向量的协方差。

取：

Q_ΔX＝ΔXΔX^T (12)

由式(11)和式(12)即可求解自适应因子α，即：

式(13)即为理论自适应因子的估计值。通过使得Q_k,k-1＝αQ_ΔX，旨在保证滤波器输出噪声的不确定度与理论噪声不确定度相当，以消除消除动力学模型预报信息与动态载体运行轨迹之间的差异。在实际过程中，这是一个不断迭代的过程，以状态向量为三维位置和速度为例，此时，α即为6*6对角线矩阵：

式中，为Q_k,k-1第i个对角线元素，为Q_ΔX第i个对角线元素。

基于卡方检验的部分粗差抗差自适应滤波方法主要包含三个过程：整体检验，观测噪声自适应和系统噪声自适应，具体流程如图3所示，图中v_j,max表示任意第j类观测值标准化残差最大值，k₀为常量，一般取1.0～1.5之间。

本实施例中，从定位结果方面进行分析，突出说明观测值之间的相关性对粗差转移所造成的影响，主要做了如下两组实验：

(1)、伪距实验分析：

实验数据采用表2中第一组数据，使用伪距单点定位模型，观测值采用原始伪距观测值，采用如下三种方案进行对比:

方案一：普通Kalman滤波方法，未添加粗差。

方案二：本发明所提方法流程框架，抗差迭代过程中，每次只对粗差最大的伪距观测值进行抗差(即部分粗差抗差方法)。在1953个历元后，随机对三颗卫星分别人为添加12m，10m，7m的粗差。

方案三：本发明所提方法流程框架，抗差迭代过程中，每次对所有粗差观测值都进行抗差，即传统方法。在1953个历元后，随机对三颗卫星分别人为添加12m，10m，7m 的粗差。

图4-6分别给出N、E、U三个方向定位偏差，图7给出整个时段卫星数变化，表4对定位误差进行统计。首先，就部分粗差抗差方法(方案二)而言，在添加粗差前后，定位结果基本保持一致，这也从侧面说明此方法可以正确识别粗差，并对粗差进行相应的降权或淘汰，降低观测异常对定位的影响，从而保证定位的可靠性。而对于传统方法(方案三)而言，在添加粗差前后，N、E、U三个方向的定位结果均出现了明显的偏移，特别是在卫星数较少的情况下，偏差更加明显，这主要是因为观测值之间的相关性，导致观测粗差的转移，从而对正常观测值产生误判，进行相应的降权或淘汰，从而引起较大的定位偏移，这也说明此种方法不能正确地对粗差进行定位。

结合表4可以发现，方案二的定位统计结果略差于方案一，这主要是方案一未添加粗差，为定位提供了更多的冗余信息，但方案二的定位统计结果明显优于方案三，这也说明了所提方法在相同类型观测值抗差过程中的鲁棒性。值得注意的是，在7000个历元之后，方案三的定位结果得到了一定的保障，结合图7可知，这主要是因为可视卫星数增多导致的，降权或淘汰粗差卫星后，仍然可以保证定位的鲁棒性。

表4伪距实验N、E、U定位误差统计表

(2)、伪距实验分析：

实验数据采用表2中第二组数据，使用多普勒+伪距的单点测速定位模型，即采用原始伪距和多普勒观测值进行单点的位置和速度解算，采用如下三种方案进行对比：

方案一：普通Kalman滤波方法，未添加粗差。

方案二：本发明所提方法流程框架，抗差迭代过程中，每次对粗差最大的伪距观测值和多普勒观测值进行抗差，即部分粗差抗差方法。在2000个历元后，随机选取某两颗卫星，伪距部分分别添加25m，18m的粗差，多普勒部分添加0.5m/s，0.42m/s的粗差。

方案三：本发明所提方法流程框架，抗差迭代过程中，每次对所有粗差观测值都进行抗差，即传统方法。在2000个历元后，随机选取某两颗卫星，伪距部分分别添加25m ，18m的粗差，多普勒部分添加0.5m/s，0.42m/s的粗差。

图8分别给出了N、E、U三个方向的定位偏差，图9分别给出了N、E、U三个方向的测速偏差，表5对定位和测速偏差进行了统计。首先，就方案三而言，在添加粗差后，无论是定位误差还是测速误差均发生了明显偏移，严重影响了定位的可靠性，实验结果较差。而方案二由于考虑到了观测值之间的相关性，在添加粗差后，定位精度和测速精度都得到了很好的保证，且在抗差过程中，每次选取伪距和多普勒粗差最大值进行抗差，从测速图中可以看出，伪距的大粗差也并未对测速精度产生影响，即不同观测值之间相关性可以忽略。需要说明的是，从定位/测速图中可以看出，方案二中个别历元实验结果会存在微弱地跳跃现象，特别是N方向的测速偏差，这可能是由于这些历元添加粗差后冗余观测信息不足，导致实验结果存在些许偏差。

从表5中的统计结果可以看出，方案一无论是定位精度还是测速精度均表现最优，这主要是因为方案一未添加粗差同时原始数据质量较优。需要说明的是，相较于方案一，方案二中U方向的定位误差明显较优，这可能是由于动态调节了观测值之间的权重所导致的。总体而言，方案二定位和测速精度略劣于方案一，但明显优于方案三，这也说明了所提方法在不同类型观测值抗差过程中的鲁棒性。

表5伪距+多普勒实验N、E、U定位和测速误差统计表

Claims (4)

1.一种基于卡方检验的部分粗差抗差自适应滤波方法，其特征在于，该方法包括如下步骤：

(1)使用大地测量型接收机采集卫星数据，通过对原始观测值进行处理，从而构造观测方程；

(2)根据假设检验理论，构造kalman滤波模型整体检验量，采用卡方检验判断整体模型是否存在异常，若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常，T_k指的是整体假设检验统计量，T_l指的是整体假设检验统计量阈值；

(3)当判定模型存在异常时，应首先进行观测噪声自适应，对由于检验量之间的相关性所导致的粗差误判进行部分观测值粗差抗差；

(4)当动力学模型存在异常扰动时，对系统噪声进行自适应，以消除动力学模型预报信息与动态载体运行轨迹之间的差异。

2.根据权利要求1所述的一种基于卡方检验的部分粗差抗差自适应滤波方法，其特征在于，步骤(2)的具体方法如下：

设置Kalman滤波中新息向量及其协方差矩阵表示：

V_k,k-1＝A_kX_k,k-1-L_k (1)

式中，V_k,k-1表示新息向量；A_k为观测设计矩阵；X_k,k-1为状态预报向量；L_k为观测向量；表示新息向量协方差矩阵；R_k表示观测向量协方差矩阵；Q_k,k-1为状态预报向量协方差；

当Kalman滤波整体模型未发生异常时，新息向量满足零均值的高斯分布；当异常发生时，异常不会改变观测值的分布类型，但使其概率分布产生一定的偏移，即如下假设检验问题：

式中，H₀是原假设，表示整体模型无异常；H₁是备择假设，表示整体模型存在异常；λ为概率分布偏移量，统计学上称为非中心化参数；N表示正态分布；

若已知先验单位权方差为σ²，取整体假设检验统计量为：

其中，T_k满足自由度为t的χ²分布；

根据给定的显著水平α、自由度t和分布类型可以确定阈值：

若T_k≤T_l，则认为无异常发生，否则应判定模型存在异常。

3.根据权利要求2所述的一种基于卡方检验的部分粗差抗差自适应滤波方法，其特征在于，步骤(3)的具体方法如下：

设第i颗卫星观测模型异常检验量为：

式中，n表示观测值维数，e_i表示第i个元素为1，其它元素为0的n_×1矩阵；

同理，第j颗卫星观测模型异常检验量为：

根据误差传播定律可得两个检验量之间的相关系数为：

在观测噪声自适应过程中，由于观测值存在相关性ρ_i,j，对不同类型观测值进行分类，并且每次只对粗差最大的观测值进行抗差。

4.根据权利要求3所述的一种基于卡方检验的部分粗差抗差自适应滤波方法，其特征在于，步骤(4)的具体方法如下：

采用自适应滤波方法，基于状态不符值构造自适应因子动态调节状态预报向量和观测向量的权比，消除动力学模型预报信息与动态载体运行轨迹之间的差异，具体实施方法如下：

预测状态向量不符值可表示为：

ΔX＝X_k,k-X_k,k-1 (10)

式中，ΔX为状态向量不符值，X_k,k为抗差Kalman滤波解，X_k,k-1为状态预报向量；

在自适应滤波过程中，状态预测向量X_k,k-1的协方差矩阵应该等于实际预测向量的偏离程度，即：

Q_ΔX＝Q_k,k-1 (11)

式中，Q_ΔX为ΔX的方差-协方差向量，Q_k,k-1为状态预报向量的协方差；

取：

Q_ΔX＝ΔXΔX^T (12)

由式(11)和式(12)即可求解自适应因子α，即：

式(13)即为理论自适应因子的估计值，通过使得Q_k,k-1＝αQ_ΔX，保证滤波器输出噪声的不确定度与理论噪声不确定度相当，以消除动力学模型预报信息与动态载体运行轨迹之间的差异。