陈含智，孙蕊，*，邱明，毛继志，胡浩亮，张立东
(1.南京航空航天大学民航学院，南京 211106；
2.中国航空无线电电子研究所民航空管航空电子技术实验室，上海 201109)
全球卫星导航系统(global navigationsatellite system，GNSS)有着广泛的用途，在军事和民用领域发挥着重要作用。民航、自动驾驶等交通领域对GNSS观测数据的质量有较高的要求[1]。然而，观测环境容易造成含有较大误差的观测数据，即故障数据，如环境反射的卫星信号造成的多路径效应[2]。GNSS质量控制旨在保证复杂环境下GNSS 的观测数据质量。故障检测算法作为GNSS 质量控制的重要组成部分，通过检测和排除故障观测数据获得高精度可靠的定位结果，保障GNSS 安全地应用于民航、自动驾驶等交通领域。目前为止，国内外学者已提出多种卫星定位故障检测算法。其中，Lee[3]提出的伪距比较法、Parkinson 和Axelrad[4]提出的最小二乘残差法和Sturza[5]提出的奇偶矢量法都是对当前观测值进行一致性检验，被称为“快照法”[6]。一些学者又提出基于历史信息和当前信息一致性检验的“滤波法”。例如，针对“快照法”难以检测微小误差[7]，沙海等[8]提出一种基于滑动窗内多个历元残差矢量累积的抗差卡尔曼滤波方法。Bhattacharyya 和Gebre-Egziabher[9]提出基于卡尔曼滤波的故障检测与保护水平计算方法。针对卡尔曼滤波不适用于非高斯问题[10]，王尔申等[11-12]又提出基于粒子滤波的故障检测法。刘江等[13]提出了适用于非线性观测特征的基于容积卡尔曼滤波的故障检测法。吴云[14]比较了“快照法”和“滤波法”，通过基于静态和动态数据的实验结果得出，“滤波法”的检测和识别性能优于“快照法”的结论。
基于卡尔曼滤波的故障检测法中，观测噪声方差矩阵影响故障检测量和卡尔曼增益矩阵的确定，进一步影响故障检测效果和滤波结果的精度，因此，准确确定观测噪声方差矩阵具有较大的意义。不同观测条件下的观测噪声可能有较大差别，例如高速无人机上搭载接收机的观测噪声应明显大于静态接收机的观测噪声。因此，按照经验确定观测噪声方差矩阵并将其固定的做法并不合理，有可能造成滤波结果精度下降，甚至故障漏检和误检。
为改进基于卡尔曼滤波的故障检测法将观测噪声方差矩阵固定的策略，本文提出基于自适应噪声方差的卫星定位故障检测算法，利用滑动窗内的历史新息序列实时估计观测噪声方差矩阵，再重构故障检测量与识别量，与根据误警率确定的阈值比较来检测和识别故障，完成故障检测后计算卡尔曼增益矩阵，最后利用排除故障后的新息数据进行测量更新并输出滤波的状态向量作为定位结果。
1.1 基于自适应噪声方差的卫星定位故障检测算法框架
本文所提算法的框架如图1 所示：首先初始化滤波参数，预测下一历元的状态向量，再计算新息和观测噪声方差矩阵，接着计算故障检测量并通过与故障检测阈值比较以判断是否存在故障数据，若存在故障数据，则进行故障识别并排除故障数据，完成故障检测与识别后，计算该历元状态向量的最终估计值，输出该历元的定位结果。
图1 本文算法框架Fig.1 Framework of the proposed algorithm
1.2 状态方程和测量方程
设k历元状态向量为
式中：xk、yk和z k分 别为k历元时接收机在地心地固坐标系中的x轴、y轴和z轴坐标；
rk为k历元的接收机钟误差与光速的乘积；
和分别为k历元时接收机在地心地固坐标系中的x轴、y轴和z轴坐标对时间的导数，即接收机k历元时在x轴、y轴和z轴方向上的速度分量；
为k历元时接收机钟误差与光速的乘积对时间的导数。
建立状态方程为
式中：X k+1为k+1历 元的状态向量；
W k为过程噪声，假设W k为高斯白噪声，其协方差矩阵为Q；
Φk+1|k为状态转移矩阵，其式为
式中：I4×4为 4 阶单位矩阵；
t为k历元与k+1历元时间间隔；
04×4为4 阶零矩阵。
建立观测方程如下：
式中：Z k+1为 观测向量；
H k+1为 观测矩阵；
V k+1为伪距观测噪声，假设V k+1为高斯白噪声，其协方差矩阵为R k+1。
观测向量为
式 中：Nk+1为k+1历 元 可 见 星 的 数 量；
ρk+1,i(i=1,2,···,Nk+1)为k+1历 元可见卫星i经过误差改正后的伪距；
ρ˜k+1,i(i=1,2,···,Nk+1)为 (xk+1,yk+1,z k+1)的近似位置(xk,yk,z k)到k+1历 元的可见星i的距离；
mk+1,i、nk+1,i、lk+1,i(i=1,2,···,Nk+1) 分 别为 (xk,yk,z k)到k+1历元可见星i的方向余弦。
观测矩阵如下：
1.3 观测噪声方差矩阵
自适应滤波的一个目的是利用观测信息，实时估计和修正噪声统计特性或增益矩阵[15]。本文通过历史新息序列实时估计观测噪声方差矩阵。
根据状态方程，预测k+1历元的状态向量为
式中：P k为估计值的误差协方差矩阵。
计算新息如下：
理论上新息的协方差矩阵为
根据极大似然准则，Cηk+1在长度为L的滑动窗（假设滑动窗内观测的可见星相同）内的最优估计值为[16]
式中：j0=k−L+1。
为使滑动窗内最新的新息对新息协方差的估计起到更大的影响，对式(11)做如下调整：
于是，R k+1的估计值为
假设各颗卫星的伪距观测噪声互不相关，同时为了避免估计结果超出合理范围，对估计结果进行如下约束，得到修正后的伪距观测噪声协方差矩阵为
1.4 故障检测
构建故障检测量为
由统计理论可知，当k+1历元各伪距观测值均不含故障时：
式中：
χ2(Nk+1−4) 为自由度为Nk+1−4的卡方分布。
给定误警率PFA，故障检测阈值TD可根据式（18）确定：
式中：fχ2(Nk+1−4)(x)为Nk+1−4自由度的卡方分布的概率密度函数。
于是可以得到故障检测准则：若 λk+1≤TD，则该历元无故障发生；
若λk+1＆gt;TD，则该历元有故障发生。
1.5 故障识别
若检测到有故障，则进一步识别故障。构建故障识别量为
式中：
ηk+1(i)为 新息向量 ηk+1的 第i个 卫星；
(i,i)为新息协方差矩阵的第i行 第i列 卫星，i=1,2,···,Nk+1。
统计量 φk+1,i在卫星i的观测数据无故障时服从标准正态分布。于是可以得到如下故障识别准则：若φk+1,i≤TR， 则卫星i的 伪距不含故障；
若φk+1,i＆gt;TR, 则卫星i的 伪距含故障。其中，TR为 统计量φk+1,i的阈值，TR可根据式（20）确定：
式中：
e 为自然底数。
1.6 算法步骤
本节给出本文所提算法的详细步骤。
步骤 1对滤波参数进行初始化：令k=0，确定初始状态向量， 根据经验确定P0、Q、L、PFA、、β1、β2。
步骤 2根据式(7)预测k+1历元的状态向量，根据式(8)计算P k+1|k。
步骤 3根据式(9)计算新息ηk+1。
步骤 4确定观测噪声方差矩阵。若k+1历元的所有可见星对应滑动窗内的新息序列的长度均等于L，根据式(14)计算， 否则根据式(15)计算。
步骤 5根据式(16)计算故障检测量，并根据故障检测准则判断该历元是否存在故障数据。
步骤 6若判断该历元不存在故障数据，则直接进入步骤8；
否则根据式(19)计算各颗卫星的故障识别量，并根据故障识别准则识别故障卫星。
步骤 7 对 ηk+1进行如下调整：
步骤 8更新k+1历 元的状态向量，输出定位结果：
增益矩阵为
于是测量更新后的k+1历元状态向量为
步骤 9k增加1 重复步骤2。
本文设计了静态和动态2 个实验，通过向采集的真实静态和动态数据中加入相应的伪距粗差，评估本文所提算法的故障检测性能。
本文定义了评价算法性能的指标。
故障检测比率（fault detectionrate,FDR）为
式中：nF为 一次实验中故障发生的历元数；
nFD为一次实验中故障发生且算法成功地检测到故障的历元数。
故障识别比率（fault identificationrate,FIR）为
式中：nFI为一次实验中故障发生且算法正确地识别故障数据的历元数。
虚警比率（false alarmrate,FAR）为
式中：nN为一次实验中没有故障发生的历元数；
nFA为一次实验中没有故障发生，但算法检测到故障的历元数。
漏检比率（m issed detectionrate,MDR）为
式中：nMD为一次实验中故障发生但是算法未检测到故障的历元数。易知RFD+RMD=100%。
故障检测算法及定位结果的评估法均以最小二乘残差法和基于卡尔曼滤波的故障检测法为对照算法。其观测噪声方差矩阵均设置为I Nk+1×Nk+1。
阶跃故障和斜坡故障是2 种常见的故障类型。根据同时存在故障数量，又可分为单故障和多故障。为了充分验证本文所提算法在不同使用环境（静态观测和动态观测）和不同故障类型下的性能，本文设计如表1 所示的实验场景。
表1 本文的实验场景Table 1 The designed scenarios
2.1 基于静态数据的仿真实验
本文采用的静态数据于2019 年12 月29 日在南京航空航天大学将军路校区图书馆附近观测获得。图2 显示了实验所用的GNSS 接收机的位置和观测环境。接收机的GNSS 板卡型号为NovAtel OEM 7500，采样率为10Hz。选取了其中上午10:45—10:50 共300 s 的GPS 观测数据进行仿真实验。
图2 静态观测实验的环境Fig.2 The environment of static observation experiment
接收机的参考位置通过千寻提供的网络实时动 态 差 分 技 术（network real timekinematic,NRTK）对接收机持续1 h 的GNSS 观测数据进行解算后得到。本文所提算法对首历元的数据利用最小二乘法进行解算得到初始状态向量，设置初始滤波参数P0=diag(1,1,1,1,0,0,0,1)（ d iag表示对角矩阵，()中的元素为对角线元素），Q=0.1I8×8，L=100，PFA=0.01%，=1m2， β1=5 ， β2=0.1。
2.1.1 原始静态数据的实验
3 种故障检测法对原始静态数据进行故障检测时的故障检测量如图3 所示。实验中，他们的故障检测量均未超过阈值，虚警比率均为0。
图3 原始静态数据的故障检测量Fig.3 Fault detection statistics based on the original static data
图4 给出了本文所提算法对原始数据进行解算后定位误差的变化情况。表2 对比了3 种算法的定位精度指标，可知本文所提算法的水平均方根误差为0.191m，水平方向上的定位精度相较于最小二乘残差法、基于卡尔曼滤波的故障检测法都提高了13.96%；
所提算法的3D 均方根误差为0.685m，3D 定位精度与前2 种算法相比分别提高了26.11%和25.46%。
图4 原始静态数据下所提算法的定位误差Fig.4 The positioning error of the proposed algorithm based on the original static data
表2 原始静态数据的定位精度指标Table 2 The positioning accuracy for the candidate algorithms based on the original static data m
2.1.2 含单阶跃故障静态数据的实验
往可见星G09 在观测时刻150～200s（含1 5 0 s和200s）内的伪距添加5 m 的阶跃故障后，使用本文所提算法和2 种对照算法进行故障检测。故障检测量如图5 所示，150 s 时所提算法的故障检测量阶跃至135m 左右超过阈值，故障检测量超出其阈值一直持续到200 s。基于卡尔曼滤波的故障检测法与最小二乘残差法的故障检测量虽然也在故障发生期间阶跃升高，但却未超过故障检测阈值，因此未检测到故障。由图6 可知，150～200 s 内所提算法在检测到故障后，G09 的故障识别量超过其阈值，被正确地标记为观测数据含故障的卫星。如表3所示，此时所提算法的故障检测比率和故障识别比率均为100%，虚警比率为0；
而最小二乘残差法和基于卡尔曼滤波的故障检测法未能检测到故障。
表3 静态数据含5 m 阶跃故障下的故障检测指标Table 3 Fault detection perform ance for a 5 m step error added on the static data %
图5 静态数据含5 m阶跃故障的故障检测量Fig.5 Fault detection statistics with a 5 m step error added on observations in the static mode
图6 静态数据含5 m阶跃故障下所提算法的故障识别量Fig.6 Fault identification statistics of the proposed algorithm with a 5 m step error added on the observations in the static mode
由表4 可知，所提算法此时的水平均方根误差为0.192m，3D 均方根误差为0.6 8 7m。对比对照算法的精度指标知，相较于最小二乘残差法、基于卡尔曼滤波的故障检测法，所提算法的水平定位精度分别提高82.14%和82.17%，3D 定位精度分别提高71.67%和71.68%，其中，本文的定位精度均根据均方根误差计算得到。
表4 静态数据含5 m 阶跃故障下的定位精度指标Table 4 The positioning accuracy for the candidate algorithm s based on the static data w ith a 5 m step error m
若向可见星G09 从观测时刻150～200s（含150 s和200 s）内的伪距依次添加 0,1,···,10 m的阶跃故障，3 种算法的故障检测比率、故障识别比率、3D 均方根误差如图7 所示。3 种算法的虚警比率始终为0，漏检比率因其和故障检测比率之和为1。由图7 可知所提算法和基于卡尔曼滤波的故障检测法分别自单星伪距故障超过3m 和8 m 后达到100%的故障检测比率和故障识别比率。最小二乘残差法在单星伪距故障为8m 时，故障检测比率达到100%；
在单星伪距故障为9m 时，故障识别比率达到100%。所提算法在故障小于2m 时3D 均方根误差随故障的增大而增大，当故障超过3m 后，由于所提算法能100%地检测并识别故障，3D 均方根误差保持在0.69m。基于卡尔曼滤波的故障检测法的3D 均方根误差随着故障增加至6m 而逐渐上升至3m 左右，在故障超过7m 后3D 均方根误差下降至1m 左右。最小二乘残差法的3D 均方根误差在7m 故障时超过3m，随后随着故障的增加而降低。所提算法的3D 均方根误差曲线始终在另外2 种算法的下方，因此，所提算法在不同的故障条件下均具有最高的定位精度。
图7 静态数据含不同大小的阶跃故障时的故障检测比率、故障识别比率和3D均方根误差Fig.7 The FDR, FIR and 3D accuracy (RMSE) for the candidate algorithms based on the static data with various step errors
2.1.3 含单斜坡故障的静态数据下的实验
若在观测时刻100～200s（不含1 0 0 s，含2 0 0 s）内往可见星G09 的伪距中添加以0.2m/s 的速率增大的斜坡故障，3 种算法的故障检测量如图8 所示（TSE为 斜坡故障的持续时间；
VSE为斜坡故障的速率）。故障发生期间，3 种算法的故障检测量均逐渐上升，并在150 s 之前超过阈值。所提算法开始检测到斜坡故障的时间要略微早于另外2 种算法。
图8 静态数据在T SE=100 s，V SE=0.2 m/s故障下的故障检测量Fig.8 Fault detection statistics for the static data with a 0.2 m/s ramp error added for 100 s
图9 中所提算法计算得到的G09 卫星的故障识别量和其故障检测量一样于观测时刻150 s 之前超过阈值，其伪距被识别出含有故障。图10 说明基于卡尔曼滤波的故障检测法大约在观测时刻160 s后正确识别故障；
识别出故障后，G09 的故障识别量经过了短时间的快速上升。图11 说明最小二乘残差法仅在观测时刻200s 前的大约2 0 s 时间内正确识别了故障。
图9 静态数据在T SE=100 s，V SE=0.2 m/s故障下所提算法的故障识别量Fig.9 Fault identification statistics of the proposed algorithm with a 0.2 m/s ramp error added for 100 s in the static mode
图10 静态数据在T SE=100 s，V SE=0.2 m/s故障下基于卡尔曼滤波的故障检测法的故障识别量Fig.10 Fault identification statistics of the KF-based method with a 0.2 m/s ramp error added for 100 s in the static mode
图11 静态数据在T SE=100 s，V SE=0.2 m/s故障下最小二乘残差法的故障识别量Fig.11 Fault identification statistics of the least square residual based method with a 0.2 m/s ramp error added for 100 s in the static mode
表5 计算了3 种算法的故障检测指标。最小二乘残差法和基于卡尔曼滤波的故障检测法的故障检测比率均为57.8%。虽然所提算法的故障检测比率为61.6%较这2 种算法没有明显地提高，但是本文所提算法的故障识别比率为51.4%，明显地高于其他2 种算法20%和35.2%的故障识别率。
表5 静态数据在T SE = 100 s，V SE = 0.2 m/s故障下的故障检测指标Table 5 Fau lt detection performance for the candidate algorithm s with a 0.2 m/s ram p error added for 100 s in the static mode %
表6 对比了该种故障下3 种算法的定位精度指标。所提算法的3D 均方根误差为0.712m，定位精度分别比2 种对照算法提高72.2%和64.4%。
表6 静态数据在T SE = 100 s，V SE = 0.2 m/s故障下的定位精度指标Table 6 The positioning accuracy for the candidate algorithm s w ith a 0.2 m/s ram p error added for 100 s in the static mode m
图12 显示了当卫星G09 的观测时刻100～200 s（不含100s，含2 0 0 s）内的斜坡故障速率从0 m/s增加到0.5m/s 时的实验结果。3 种算法的故障检测比率比较接近，但所提算法的故障识别比率要明显高于其他2 种算法。2 种对照算法的3D 均方根误差先随斜坡故障的速率的增加而增加，然后降低，最大值超过了2m，但是所提算法的3D 均方根误差始终稳定地保持小于1m。
图12 静态数据含不同大小的斜坡故障时的故障检测比率、故障识别比率和3D均方根误差Fig.12 The FDR, FIR and 3D accuracy (RMSE) for candidate algorithms with various ramp errors in the static mode
2.1.4 含多故障的静态数据的实验
最小二乘残差法基于单故障假设，每个历元只能排除一个故障，因此，2.1.4 节和2.2.4 节中的实验只将基于卡尔曼滤波的故障检测法作为对照组，验证所提算法对多故障的检测性能。
观测时刻150～200s（含150s 和200s）内向卫星G09 和G16 的伪距中添加5m 的阶跃故障，2 种算法的故障检测量如图13 所示。故障发生期间，所提算法的故障检测量始终高于阈值，但是基于卡尔曼滤波的故障检测法未检测出故障。
图13 静态数据中卫星G09和G16的伪距同时含5 m的故障下的故障检测量Fig.13 Fault detection statistics for 5 m step errors added on the pseudorange of G09 and G16 in the static mode
如图14 所示，故障发生期间，所提算法的卫星G09 和G16 的故障识别量均超过了阈值。
图14 静态数据中卫星G09和G16的伪距同时含有5 m的阶跃故障下所提算法的故障识别量Fig.14 Fault identification statistics for 5 m step errors added on pseudorange of G09 and G16 in the static mode
由表7 可知，所提算法能100%地检测和识别2 颗卫星的伪距同时存在的5m 的阶跃故障。而基于卡尔曼滤波的故障检测法的故障检测检测比率和识别比率均为0。
表7 静态数据中卫星G09 和G16 的伪距同时有5 m 的阶跃故障下的故障检测指标Table 7 Fault detection performance for candidate algorithm s w ith 5 m step errors added on pseudorange of G09 and G16 in static mode %
2 种算法的定位精度指标如表8 所示。此时基于卡尔曼滤波的故障检测法的3D 均方根误差为1.773m；
所提算法的3 D 均方根误差为0.7 0 4 m，定位精度大约均提高了60.3%
表8 静态数据中卫星G09 和G16 的伪距同时含有5 m的阶跃故障下的定位精度指标Table 8 Positioning accuracy performance for candidate algorithm s w ith 5 m step errors added on pseudorange of G09 and G16 in static mode m
图15 为向卫星G10 和G12 的观测时刻150～200s（含150 s 和200 s）内的伪距添加的阶跃故障逐渐增大后的故障检测和定位情况。不同大小的多故障下，2 种算法的虚警比率均为0。基于卡尔曼滤波的故障检测法能100%地检测超过8m 的多故障，能100%地识别超过10m 的多故障。所提算法分别在多故障达到3m 和4m 时能100%地检测和识别多故障。多故障低于2 m 时，2 种算法的3D 均方根误差比较接近，但是随着多故障的进一步增加，基于卡尔曼滤波的故障检测法的3D 均方根误差在多故障低于8m 时一直上升至高于2 m；
而所提算法的3D 均方根始终大致地保持在低于1m 的水平。
图15 静态数据含不同大小的多故障时故障检测比率、故障识别比率和3D均方根误差Fig.15 FDR, FIR and 3D accuracy (RMSE) for candidate algorithms with various step errors on two satellites in the static mode
2.2 基于动态数据的仿真实验
本文采用的动态数据来自2018 年1 月31 日于台湾南投进行的无人机动态观测实验。无人机的飞行轨迹如图16 所示，机载GNSS 接收机型号为Trimble BD9 8 2，采样率为1 0 Hz。实验中使用的无人机的参考轨迹是通过机载VLP-16 Velodyne 激光雷达进行近景摄影测量获得的，具有cm 级的定位精度。
图16 无人机的飞行轨迹Fig.16 Flight trajectory of the unmanned aerial vehicle (UAV)
所提算法对第1、第2 个历元的数据利用最小二乘法进行解算得到初始状态向量，设置初始滤波参数P0=diag(10,10,10,10,2,2,2,2)，Q=0.1I8×8，L=100，PFA=0.01%，=5 m2，β1=2 ，β2=0.2。
2.2.1 原始动态数据下的实验
由图17 可知，原始动态数据下，3 种算法整个过程故障检测量均未超过故障检测阈值，虚警比率都为0。
图17 原始动态数据的故障检测量Fig.17 Fault detection statistics based on original dynamic data
所提算法对原始动态数据进行定位解算后各方向定位误差的变化情况如图18 所示。定位精度指标如表9 所示，表9 中所提算法的水平均方根误差和3D 均方根误差分别为1.567m 和2.233m，对比最小二乘残差法和基于卡尔曼滤波的故障检测法，所提算法水平方向上的定位精度分别提高1.14%和0.57%，3D 方向上的定位精度分别提高33.4%和32%。
表9 原始动态数据的定位精度指标Table 9 Positioning accuracy performance for candidate algorithm s based on original dynam ic data m
图18 原始动态数据下所提算法的定位误差Fig.18 Positioning error of the proposed algorithm based on original dynamic data
2.2.2 含单阶跃故障的动态数据下的实验
往可见卫星G10 在观测时刻300～400s（含300s和400s）内的伪距添加1 0 m 的阶跃故障后，使用本文所提算法和2 种对照算法进行故障检测。
如图19 所示，观测时刻300～400 s 内所提算法的故障检测量均高于故障检测阈值，其余时间内低于阈值。对照算法中，基于卡尔曼滤波的故障检测法在故障发生期间其故障检测量均低于阈值，故障检测比率为0，而最小二乘残差法具有较高的故障检测比率。如图20 所示，所提算法的卫星G10 的故障识别量在300～400 s 内均超过相应阈值，其他可见卫星的故障识别量在此段时间内均低于相应阈值，因此，所提算法能100%正确识别故障。图21表明最小二乘残差法在一半左右的故障发生的历元上未能正确识别故障。
图19 动态数据含10 m阶跃故障下的故障检测量Fig.19 Fault detection statistics with a 10 m step error added on observations in the dynamic mode
图20 动态数据含10 m阶跃故障下所提算法的故障识别量Fig.20 Fault identification statistics with a 10 m step error added on observations in the dynamic mode
图21 动态数据含10 m阶跃故障下最小二乘残差法的故障识别量Fig.21 Fault identification statistics of least square residual based method with a 10 m step error added on the observations in the dynamic mode
如表10 所示，所提算法在处理单星伪距含10m故障的动态数据时，故障检测效果较好，故障检测比率和故障识别比率均为100%，并且虚警比率为0。最小二乘残差法的故障检测比率虽然也达到了98.7%，但其故障识别比率却不到60%；
基于卡尔曼滤波的故障检测法未能检测到故障。
表10 动态数据含10 m 阶跃故障下的故障检测指标Table 10 Fault detection perform ance for candidate algorithm s w ith a 10 m step error added on the observations in the dynam ic mode %
表11 给出了所提算法与2 种对照算法在处理含故障的动态数据时的定位精度指标，此时所提算法定位结果的水平均方根误差和3D 均方根误差分别为1.592m 和3.049m，相较于最小二乘残差法，水平定位精度和3D 定位精度分别提高2.63%和39.87%；
相较于基于卡尔曼滤波的故障检测法，水平定位精度和3D 定位精度分别提高和3.28%和52.33%。
表11 动态数据含10 m 阶跃故障下的定位精度指标Table 11 Positioning accuracy performance for candidate algorithm s w ith a 10 m step error added on observations in dynam ic mode m
若向可见星G10 从300～400s（含300s 和400s）内的伪距依次添加 0,2,···,18 m的阶跃故障，3 种算法的故障检测效果如图22 所示。由图22 可知，为使故障检测比率和故障识别比率均达到100%，所提算法和基于卡尔曼滤波的故障检测法分别需要故障超过10m 和1 6m。最小二乘残差法在故障超过10m 时故障检测比率达到1 0 0%，在故障超过1 2 m时故障识别比率达到100%。所提算法的3D 均方根误差在故障增加到2m 的过程中经历了较为缓慢的上升，随后便基本稳定在3m。但是基于卡尔曼滤波的故障检测法的3D 均方根误差在故障达到14m 前一直都在上升，最高点的3D 均方根误差超过了7m，当故障超过了1 6 m，其3D 均方根误差下降至3.5m 左右。最小二乘残差法的3D 均方根误差在8m 故障时达到最大值6 m。故障从0 m 增加至18m 的过程中，所提算法的3D 均方根误差始终是3 种算法中最低的，所提算法在不同的故障条件下具有更高的定位精度。
图22 动态数据含不同大小的阶跃故障时的故障检测比率、故障识别比率、3D均方根误差Fig.22 The FDR, FIR and 3D accuracy (RMSE) for the candidate algorithms based on dynamic data with various step errors
2.2.3 含单斜坡故障的动态数据下的实验
本节实验中向卫星G10 的伪距添加从观测时刻400～6 0 0 s（不含4 0 0 s，含6 0 0 s）的速率为0.2m/s的斜坡故障。
3 种算法对含斜坡故障的动态数据的故障检测量如图23 所示。所提算法的故障检测量从400 s开始斜坡上升，并在450 s 左右超过阈值，所提算法开始检测到斜坡故障。对照算法和所提算法有类似的故障检测量的变化情况；
开始检测到斜坡故障的时间比较接近。
图23 动态数据在T SE=200 s，V SE=0.2 m/s故障下的故障检测量Fig.23 Fault detection statistics for the dynamic data with a 0.2 m/s ramp error added for 200 s
如图24 所示，所提算法计算的卫星G10 的故障识别量从400 s 开始斜坡上升，并大致在所提算法检测到故障的同一时间超过阈值，此后一直到故障不再发生，所提算法均能正确地识别故障。如图25所示为基于卡尔曼滤波的故障检测法识别出故障的时间要晚于所提算法。最小二乘残差法在斜坡故障增大到一定值后，多颗卫星的故障识别量超过阈值，此时将他们中的最大值对应的数据标记为故障数据。最小二乘残差法在故障发生期间的故障识别结果如图26 所示。最小二乘残差法在440 s左右首次正确识别故障，然后，一小段时间内在正确识别故障和错误识别故障间波动，过后便能一直正确识别含斜坡故障的数据。
图24 动态数据在T SE=200 s，V SE=0.2 m/s故障下所提算法的故障识别量Fig.24 Fault identification statistics of the proposed algorithm with a 0.2 m/s ramp error added for 200 s in the dynamic mode
图25 动态数据在T SE=200 s，V SE=0.2 m/s故障下基于卡尔曼滤波的故障检测法的故障识别量Fig.25 Fault identification statistics of the KF-based method with a 0.2 m/s ramp error added for 200 s in the dynamic mode
图26 动态数据在T SE=200 s，V SE=0.2 m/s故障下最小二乘残差法的故障识别结果Fig.26 Fault identification result of the least square residual based method with a 0.2 m/s ramp error added for 200 s in dynamic mode
由表12 可知，在该种斜坡故障条件下，所提算法的故障检测比率和识别比率均为66.25%。虽然最小二乘残差法和基于卡尔曼滤波的故障检测法的故障检测比率都略高于所提算法，但是这2 种算法的故障识别比率都明显低于所提算法。
表12 动态数据在 T SE= 200 s，V SE = 0.2 m/s故障下的故障检测指标Table 12 Fault detection perform ance for candidate algorithm w ith a 0.2 m/s ram p error added for 200 s in dynam ic m ode %
表13 计算了3 种算法的定位精度指标，2 种对照算法的3D 均方根误差分别为4.253m 和3.689m，然而所提算法的3D 均方根误差为2.552m，相较于2 种对照算法分别降低40%和30.2%。
表13 动态数据在 T SE= 200 s，V SE = 0.2 m/s故障下定位精度指标Table 13 Positioning accuracy for candidate algorithm s w ith a 0.2 m/s ramp error added for 200 s in dynam ic mode m
图27 中，若向卫星G10 的观测时刻400～600 s（不含400s，含6 0 0 s）伪距中添加速率逐渐增加的斜坡故障，最小二乘残差法的故障检测比率一直略高于所提算法，但是所提算法的故障识别比率始终高于另外2 种对照算法。此外所提算法的定位精度也始终高于2 种对照算法。
图27 动态数据含不同大小的斜坡故障时的故障检测比率、故障识别比率和3D均方根误差Fig.27 The FDR, FIR and 3D accuracy (RMSE) for candidate algorithms with various ramp errors in the dynamic mode
2.2.4 含多故障的动态数据下的实验
本节比较分析所提算法对多故障的检测效果。如图28 所示，当向卫星G10 和G12 在观测时刻400～5 0 0 s（含4 0 0 s 和5 0 0 s）的伪距中注入1 0m的阶跃故障时，基于卡尔曼滤波的故障检测法仅在十分少的历元检测到了故障；
而所提算法能100%地检测到多故障。
图28 动态数据中卫星G10和G12同时含有10 m的故障下的故障检测量Fig.28 Fault detection statistics for 10 m step errors added on pseudorange of G10 and G12 in the dynamic mode
所提算法和基于卡尔曼滤波算法的故障检测量计算的故障识别量分别如图29 和图30 所示。由图29 可知所提算法基本能正确识别多故障，只有在少量历元，卫星G12 的故障识别量没有超过阈值，没有被识别出来。但是在图30 中，只有十分少量的历元内，卫星G12 的故障识别量超过了阈值，即使是在这些历元内，多故障也不能被完全识别。
图29 动态数据中卫星G10和G12的伪距同时含有10 m的阶跃故障下所提算法的故障识别量Fig.29 Fault identification statistics of the proposed algorithm with 10 m step errors added on the pseudorange of G10 and G12 in dynamic mode
图30 动态数据中卫星G10和G12的伪距同时含有10 m的阶跃故障下基于卡尔曼滤波的故障检测法的故障识别量Fig.30 Fault identification statistics of KF-based method with 10 m step errors added on pseudorange of G10 and G12 in dynamic mode
由表14 可知基于卡尔曼滤波的故障检测法的故障检测比率仅为3.7%，故障识别比率为0。所提算法的故障检测比率和故障识别比率分别为100%和85.6%，显著优于基于卡尔曼滤波的故障检测法。
表14 动态数据中卫星G10 和G12 的伪距同时含有10 m阶跃故障下的故障检测指标Table 14 Fault detection performance for candidate algorithm s w ith 10 m step errors added on pseudorange of G10 and G12 in dynam ic mode %
表15 比较了2 种算法的定位精度指标。所提算法的3D 均方根误差为2.645m，而基于卡尔曼滤波的故障检测法的3D 均方根误差为6.618m。在动态数据含有10m 的多故障条件下，所提算法的定位精度相较对照算法提高了60.03%。
表15 动态数据中卫星G10 和G12 的伪距同时含有10 m 的故障下的定位精度指标Table 15 Positioning accuracy performance for candidate algorithm s w ith 10 m step errors added on pseudorange of G10 and G12 in dynam ic mode m
为了进一步比较2 种算法的性能，图31 给出了卫星G10 和G12 的观测时刻400～500s（含4 0 0 s和500 s）的伪距含有 0,2,···,16 m的阶跃故障时，2 种算法的故障检测比率、故障识别比率和3D均方根误差。所提算法的故障检测比率和故障识别比率分别在多故障超过8m 和1 2m 时达到了100%；
基于卡尔曼滤波的故障检测法的故障检测比率和故障识别比率分别在多故障超过12m 和1 4m 时达到了100%。随着多故障从0 增大到18m，由于所提算法的滤波采取了自适应的噪声方差策略并且有能检测到较小的多故障，所提算法的3D 均方根误差一直稳定在2m 左右，但是基于卡尔曼滤波的故障检测法在多故障达到了12m 时，其3D 均方根误差也达到了峰值8m。
图31 动态数据含不同大小的多故障时故障检测比率、故障识别比率和3D均方根误差Fig.31 The FDR, FIR and 3D accuracy RMSE for candidate algorithms with various step errors on two satellites
本文提出了基于自适应噪声方差的卫星定位故障检测算法，利用滑动窗内的新息序列实时估计观测噪声方差矩阵，降低了设置观测噪声矩阵的主观性，且能更好地适应实际观测噪声随观测环境而变化。在自适应噪声方差的基础上，所提算法的故障检测量和故障识别量能更好地满足对应的统计学概率分布，具有较高的故障检测比率和识别比率，特别是对于微小故障；
而且完成故障检测后，自适应噪声方差的策略也优化了卡尔曼增益矩阵，一定程度上提高了滤波的精度。实验结果表明：
1）若给定的静态数据中单星伪距的阶跃故障逐渐增大，所提算法在阶跃故障超过3m 时，其故障检测比率和故障识别比率同时达到100%；
而基于卡尔曼滤波的故障检测法在阶跃故障超过8m时，其故障检测比率和故障识别比率才均达到100%；
最小二乘残差法分别在阶跃故障超过8m 和超过9 m时，达到了100%的故障检测比率和故障识别比率。
2）若给定的静态数据中单星伪距的斜坡故障的速率逐渐增大，所提算法和对照算法的故障检测比率的变化情况接近，但所提算法的故障识别比率明显高于对照算法。
3）若给定的静态数据所含的多故障逐渐增大，所提算法分别在多故障超过3m 和超过4 m 时100%地检测和识别故障；
而基于卡尔曼滤波的故障检测法分别在多故障超过8m 和超过1 0 m 时100%地检测和识别故障。
4）若给定的动态数据中单星伪距的阶跃故障逐渐增加，要使故障检测比率和识别比率同时达到100%，所提算法需要故障超过10m；
而基于卡尔曼滤波的故障检测法需要故障超过16m；
虽然最小二乘残差法在故障超过10m 时，其故障检测比率也达到了100%，但故障超过12m 时其故障识别比率才能达到100%。
5）若给定的动态数据中单星伪距的斜坡故障的速率逐渐增加，虽然最小二乘残差法的故障检测比率略高于所提算法，但是所提算法的故障识别比率一直明显高于另外2 种对照算法。
6）若给定的动态数据所含的多故障逐渐增大，所提算法和基于卡尔曼滤波的故障检测法分别在故障超过8m 和1 2m 时，达到100%的故障检测比率；
这2 种算法分别在故障超过12m 和14m 时，达到了100%的故障识别比率。
在本文所有故障情景下，由于所提算法采取了根据滑动窗中的历史新息估计观测噪声方差的策略，提高了基于卡尔曼滤波的故障检测法的故障检测比率和识别比率，同时优化了测量更新中各颗卫星的权重，所以应用所提算法进行质量控制后的3D均方根误差基本不随故障的增大而发生变化，其定位精度一直高于2 种对照算法。
猜你喜欢动态数据伪距方根方根拓展探究中学生数理化·七年级数学人教版(2023年3期)2023-03-21我们爱把马鲛鱼叫鰆鯃飞天(2019年6期)2019-07-08云计算环境下动态数据聚集算法研究计算机测量与控制(2018年1期)2018-02-05北斗伪距观测值精度分析测绘科学与工程(2017年3期)2017-08-16均方根嵌入式容积粒子PHD 多目标跟踪方法自动化学报(2017年2期)2017-04-04颞下颌关节三维动态数据测量的初步研究中华老年口腔医学杂志(2016年3期)2017-01-15GNSS伪距粗差的开窗探测及修复测绘通报(2016年12期)2017-01-06联合码伪距和载波宽巷组合的相对定位技术研究导航定位学报(2015年2期)2015-06-05基于动态数据驱动的突发水污染事故仿真方法浙江大学学报（工学版）(2015年1期)2015-03-01数学魔术新高考·高二数学(2014年7期)2014-09-18