基于自適應(yīng)噪聲方差的衛(wèi)星定位故障檢測法

陳含智，孫蕊，*，邱明，毛繼志，胡浩亮，張立東

(1.南京航空航天大學(xué)民航學(xué)院，南京 211106；2.中國航空無線電電子研究所民航空管航空電子技術(shù)實(shí)驗(yàn)室，上海 201109)

全球衛(wèi)星導(dǎo)航系統(tǒng)(global navigationsatellite system，GNSS)有著廣泛的用途，在軍事和民用領(lǐng)域發(fā)揮著重要作用。民航、自動駕駛等交通領(lǐng)域?qū)NSS觀測數(shù)據(jù)的質(zhì)量有較高的要求[1]。然而，觀測環(huán)境容易造成含有較大誤差的觀測數(shù)據(jù)，即故障數(shù)據(jù)，如環(huán)境反射的衛(wèi)星信號造成的多路徑效應(yīng)[2]。GNSS質(zhì)量控制旨在保證復(fù)雜環(huán)境下GNSS 的觀測數(shù)據(jù)質(zhì)量。故障檢測算法作為GNSS 質(zhì)量控制的重要組成部分，通過檢測和排除故障觀測數(shù)據(jù)獲得高精度可靠的定位結(jié)果，保障GNSS 安全地應(yīng)用于民航、自動駕駛等交通領(lǐng)域。目前為止，國內(nèi)外學(xué)者已提出多種衛(wèi)星定位故障檢測算法。其中，Lee[3]提出的偽距比較法、Parkinson 和Axelrad[4]提出的最小二乘殘差法和Sturza[5]提出的奇偶矢量法都是對當(dāng)前觀測值進(jìn)行一致性檢驗(yàn)，被稱為“快照法”[6]。一些學(xué)者又提出基于歷史信息和當(dāng)前信息一致性檢驗(yàn)的“濾波法”。例如，針對“快照法”難以檢測微小誤差[7]，沙海等[8]提出一種基于滑動窗內(nèi)多個歷元?dú)埐钍噶坷鄯e的抗差卡爾曼濾波方法。Bhattacharyya 和Gebre-Egziabher[9]提出基于卡爾曼濾波的故障檢測與保護(hù)水平計算方法。針對卡爾曼濾波不適用于非高斯問題[10]，王爾申等[11-12]又提出基于粒子濾波的故障檢測法。劉江等[13]提出了適用于非線性觀測特征的基于容積卡爾曼濾波的故障檢測法。吳云[14]比較了“快照法”和“濾波法”，通過基于靜態(tài)和動態(tài)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果得出，“濾波法”的檢測和識別性能優(yōu)于“快照法”的結(jié)論。

基于卡爾曼濾波的故障檢測法中，觀測噪聲方差矩陣影響故障檢測量和卡爾曼增益矩陣的確定，進(jìn)一步影響故障檢測效果和濾波結(jié)果的精度，因此，準(zhǔn)確確定觀測噪聲方差矩陣具有較大的意義。不同觀測條件下的觀測噪聲可能有較大差別，例如高速無人機(jī)上搭載接收機(jī)的觀測噪聲應(yīng)明顯大于靜態(tài)接收機(jī)的觀測噪聲。因此，按照經(jīng)驗(yàn)確定觀測噪聲方差矩陣并將其固定的做法并不合理，有可能造成濾波結(jié)果精度下降，甚至故障漏檢和誤檢。

為改進(jìn)基于卡爾曼濾波的故障檢測法將觀測噪聲方差矩陣固定的策略，本文提出基于自適應(yīng)噪聲方差的衛(wèi)星定位故障檢測算法，利用滑動窗內(nèi)的歷史新息序列實(shí)時估計觀測噪聲方差矩陣，再重構(gòu)故障檢測量與識別量，與根據(jù)誤警率確定的閾值比較來檢測和識別故障，完成故障檢測后計算卡爾曼增益矩陣，最后利用排除故障后的新息數(shù)據(jù)進(jìn)行測量更新并輸出濾波的狀態(tài)向量作為定位結(jié)果。

1 本文算法

1.1 基于自適應(yīng)噪聲方差的衛(wèi)星定位故障檢測算法框架

本文所提算法的框架如圖1 所示：首先初始化濾波參數(shù)，預(yù)測下一歷元的狀態(tài)向量，再計算新息和觀測噪聲方差矩陣，接著計算故障檢測量并通過與故障檢測閾值比較以判斷是否存在故障數(shù)據(jù)，若存在故障數(shù)據(jù)，則進(jìn)行故障識別并排除故障數(shù)據(jù)，完成故障檢測與識別后，計算該歷元狀態(tài)向量的最終估計值，輸出該歷元的定位結(jié)果。

圖1 本文算法框架Fig.1 Framework of the proposed algorithm

1.2 狀態(tài)方程和測量方程

設(shè)k歷元狀態(tài)向量為

式中：xk、yk和z k分 別為k歷元時接收機(jī)在地心地固坐標(biāo)系中的x軸、y軸和z軸坐標(biāo)；rk為k歷元的接收機(jī)鐘誤差與光速的乘積；和分別為k歷元時接收機(jī)在地心地固坐標(biāo)系中的x軸、y軸和z軸坐標(biāo)對時間的導(dǎo)數(shù)，即接收機(jī)k歷元時在x軸、y軸和z軸方向上的速度分量；為k歷元時接收機(jī)鐘誤差與光速的乘積對時間的導(dǎo)數(shù)。

建立狀態(tài)方程為

式中：X k+1為k+1歷 元的狀態(tài)向量；W k為過程噪聲，假設(shè)W k為高斯白噪聲，其協(xié)方差矩陣為Q； Φk+1|k為狀態(tài)轉(zhuǎn)移矩陣，其式為

式中：I4×4為 4 階單位矩陣；t為k歷元與k+1歷元時間間隔； 04×4為4 階零矩陣。

建立觀測方程如下：

式中：Z k+1為 觀測向量；H k+1為 觀測矩陣；V k+1為偽距觀測噪聲，假設(shè)V k+1為高斯白噪聲，其協(xié)方差矩陣為R k+1。

觀測向量為

式 中：Nk+1為k+1歷 元 可 見 星 的 數(shù) 量；ρk+1,i(i=1,2,···,Nk+1)為k+1歷 元可見衛(wèi)星i經(jīng)過誤差改正后的偽距； ρ?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 觀測噪聲方差矩陣

自適應(yīng)濾波的一個目的是利用觀測信息，實(shí)時估計和修正噪聲統(tǒng)計特性或增益矩陣[15]。本文通過歷史新息序列實(shí)時估計觀測噪聲方差矩陣。

根據(jù)狀態(tài)方程，預(yù)測k+1歷元的狀態(tài)向量為

式中：P k為估計值的誤差協(xié)方差矩陣。

計算新息如下：

理論上新息的協(xié)方差矩陣為

根據(jù)極大似然準(zhǔn)則，Cηk+1在長度為L的滑動窗（假設(shè)滑動窗內(nèi)觀測的可見星相同）內(nèi)的最優(yōu)估計值為[16]

式中：j0=k?L+1。

為使滑動窗內(nèi)最新的新息對新息協(xié)方差的估計起到更大的影響，對式(11)做如下調(diào)整：

于是，R k+1的估計值為

假設(shè)各顆衛(wèi)星的偽距觀測噪聲互不相關(guān)，同時為了避免估計結(jié)果超出合理范圍，對估計結(jié)果進(jìn)行如下約束，得到修正后的偽距觀測噪聲協(xié)方差矩陣為

1.4 故障檢測

構(gòu)建故障檢測量為

由統(tǒng)計理論可知，當(dāng)k+1歷元各偽距觀測值均不含故障時：

式中： χ2(Nk+1?4) 為自由度為Nk+1?4的卡方分布。

給定誤警率PFA，故障檢測閾值TD可根據(jù)式（18）確定：

式中：fχ2(Nk+1?4)(x)為Nk+1?4自由度的卡方分布的概率密度函數(shù)。

于是可以得到故障檢測準(zhǔn)則：若 λk+1≤TD，則該歷元無故障發(fā)生；若λk+1＆gt;TD，則該歷元有故障發(fā)生。

1.5 故障識別

若檢測到有故障，則進(jìn)一步識別故障。構(gòu)建故障識別量為

式中： ηk+1(i)為 新息向量 ηk+1的 第i個 衛(wèi)星；(i,i)為新息協(xié)方差矩陣的第i行 第i列 衛(wèi)星，i=1,2,···,Nk+1。

統(tǒng)計量 φk+1,i在衛(wèi)星i的觀測數(shù)據(jù)無故障時服從標(biāo)準(zhǔn)正態(tài)分布。于是可以得到如下故障識別準(zhǔn)則：若φk+1,i≤TR， 則衛(wèi)星i的 偽距不含故障；若φk+1,i＆gt;TR, 則衛(wèi)星i的 偽距含故障。其中，TR為 統(tǒng)計量φk+1,i的閾值，TR可根據(jù)式（20）確定：

式中： e 為自然底數(shù)。

1.6 算法步驟

本節(jié)給出本文所提算法的詳細(xì)步驟。

步驟 1對濾波參數(shù)進(jìn)行初始化：令k=0，確定初始狀態(tài)向量， 根據(jù)經(jīng)驗(yàn)確定P0、Q、L、PFA、、β1、β2。

步驟 2根據(jù)式(7)預(yù)測k+1歷元的狀態(tài)向量，根據(jù)式(8)計算P k+1|k。

步驟 3根據(jù)式(9)計算新息ηk+1。

步驟 4確定觀測噪聲方差矩陣。若k+1歷元的所有可見星對應(yīng)滑動窗內(nèi)的新息序列的長度均等于L，根據(jù)式(14)計算， 否則根據(jù)式(15)計算。

步驟 5根據(jù)式(16)計算故障檢測量，并根據(jù)故障檢測準(zhǔn)則判斷該歷元是否存在故障數(shù)據(jù)。

步驟 6若判斷該歷元不存在故障數(shù)據(jù)，則直接進(jìn)入步驟8；否則根據(jù)式(19)計算各顆衛(wèi)星的故障識別量，并根據(jù)故障識別準(zhǔn)則識別故障衛(wèi)星。

步驟 7 對 ηk+1進(jìn)行如下調(diào)整：

步驟 8更新k+1歷 元的狀態(tài)向量，輸出定位結(jié)果：

增益矩陣為

于是測量更新后的k+1歷元狀態(tài)向量為

步驟 9k增加1 重復(fù)步驟2。

2 仿真驗(yàn)證與分析

本文設(shè)計了靜態(tài)和動態(tài)2 個實(shí)驗(yàn)，通過向采集的真實(shí)靜態(tài)和動態(tài)數(shù)據(jù)中加入相應(yīng)的偽距粗差，評估本文所提算法的故障檢測性能。

本文定義了評價算法性能的指標(biāo)。

故障檢測比率（fault detectionrate,FDR）為

式中：nF為 一次實(shí)驗(yàn)中故障發(fā)生的歷元數(shù)；nFD為一次實(shí)驗(yàn)中故障發(fā)生且算法成功地檢測到故障的歷元數(shù)。

故障識別比率（fault identificationrate,FIR）為

式中：nFI為一次實(shí)驗(yàn)中故障發(fā)生且算法正確地識別故障數(shù)據(jù)的歷元數(shù)。

虛警比率（false alarmrate,FAR）為

式中：nN為一次實(shí)驗(yàn)中沒有故障發(fā)生的歷元數(shù)；nFA為一次實(shí)驗(yàn)中沒有故障發(fā)生，但算法檢測到故障的歷元數(shù)。

漏檢比率（m issed detectionrate,MDR）為

式中：nMD為一次實(shí)驗(yàn)中故障發(fā)生但是算法未檢測到故障的歷元數(shù)。易知RFD+RMD=100%。

故障檢測算法及定位結(jié)果的評估法均以最小二乘殘差法和基于卡爾曼濾波的故障檢測法為對照算法。其觀測噪聲方差矩陣均設(shè)置為I Nk+1×Nk+1。

階躍故障和斜坡故障是2 種常見的故障類型。根據(jù)同時存在故障數(shù)量，又可分為單故障和多故障。為了充分驗(yàn)證本文所提算法在不同使用環(huán)境（靜態(tài)觀測和動態(tài)觀測）和不同故障類型下的性能，本文設(shè)計如表1 所示的實(shí)驗(yàn)場景。

表1 本文的實(shí)驗(yàn)場景Table 1 The designed scenarios

2.1 基于靜態(tài)數(shù)據(jù)的仿真實(shí)驗(yàn)

本文采用的靜態(tài)數(shù)據(jù)于2019 年12 月29 日在南京航空航天大學(xué)將軍路校區(qū)圖書館附近觀測獲得。圖2 顯示了實(shí)驗(yàn)所用的GNSS 接收機(jī)的位置和觀測環(huán)境。接收機(jī)的GNSS 板卡型號為NovAtel OEM 7500，采樣率為10Hz。選取了其中上午10:45—10:50 共300 s 的GPS 觀測數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn)。

圖2 靜態(tài)觀測實(shí)驗(yàn)的環(huán)境Fig.2 The environment of static observation experiment

接收機(jī)的參考位置通過千尋提供的網(wǎng)絡(luò)實(shí)時動 態(tài) 差 分 技 術(shù)（network real timekinematic,NRTK）對接收機(jī)持續(xù)1 h 的GNSS 觀測數(shù)據(jù)進(jìn)行解算后得到。本文所提算法對首歷元的數(shù)據(jù)利用最小二乘法進(jìn)行解算得到初始狀態(tài)向量，設(shè)置初始濾波參數(shù)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 原始靜態(tài)數(shù)據(jù)的實(shí)驗(yàn)

3 種故障檢測法對原始靜態(tài)數(shù)據(jù)進(jìn)行故障檢測時的故障檢測量如圖3 所示。實(shí)驗(yàn)中，他們的故障檢測量均未超過閾值，虛警比率均為0。

圖3 原始靜態(tài)數(shù)據(jù)的故障檢測量Fig.3 Fault detection statistics based on the original static data

圖4 給出了本文所提算法對原始數(shù)據(jù)進(jìn)行解算后定位誤差的變化情況。表2 對比了3 種算法的定位精度指標(biāo)，可知本文所提算法的水平均方根誤差為0.191m，水平方向上的定位精度相較于最小二乘殘差法、基于卡爾曼濾波的故障檢測法都提高了13.96%；所提算法的3D 均方根誤差為0.685m，3D 定位精度與前2 種算法相比分別提高了26.11%和25.46%。

圖4 原始靜態(tài)數(shù)據(jù)下所提算法的定位誤差Fig.4 The positioning error of the proposed algorithm based on the original static data

表2 原始靜態(tài)數(shù)據(jù)的定位精度指標(biāo)Table 2 The positioning accuracy for the candidate algorithms based on the original static data m

2.1.2 含單階躍故障靜態(tài)數(shù)據(jù)的實(shí)驗(yàn)

往可見星G09 在觀測時刻150～200s（含1 5 0 s和200s）內(nèi)的偽距添加5 m 的階躍故障后，使用本文所提算法和2 種對照算法進(jìn)行故障檢測。故障檢測量如圖5 所示，150 s 時所提算法的故障檢測量階躍至135m 左右超過閾值，故障檢測量超出其閾值一直持續(xù)到200 s?；诳柭鼮V波的故障檢測法與最小二乘殘差法的故障檢測量雖然也在故障發(fā)生期間階躍升高，但卻未超過故障檢測閾值，因此未檢測到故障。由圖6 可知，150～200 s 內(nèi)所提算法在檢測到故障后，G09 的故障識別量超過其閾值，被正確地標(biāo)記為觀測數(shù)據(jù)含故障的衛(wèi)星。如表3所示，此時所提算法的故障檢測比率和故障識別比率均為100%，虛警比率為0；而最小二乘殘差法和基于卡爾曼濾波的故障檢測法未能檢測到故障。

表3 靜態(tài)數(shù)據(jù)含5 m 階躍故障下的故障檢測指標(biāo)Table 3 Fault detection perform ance for a 5 m step error added on the static data %

圖5 靜態(tài)數(shù)據(jù)含5 m階躍故障的故障檢測量Fig.5 Fault detection statistics with a 5 m step error added on observations in the static mode

圖6 靜態(tài)數(shù)據(jù)含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。對比對照算法的精度指標(biāo)知，相較于最小二乘殘差法、基于卡爾曼濾波的故障檢測法，所提算法的水平定位精度分別提高82.14%和82.17%，3D 定位精度分別提高71.67%和71.68%，其中，本文的定位精度均根據(jù)均方根誤差計算得到。

表4 靜態(tài)數(shù)據(jù)含5 m 階躍故障下的定位精度指標(biāo)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）內(nèi)的偽距依次添加 0,1,···,10 m的階躍故障，3 種算法的故障檢測比率、故障識別比率、3D 均方根誤差如圖7 所示。3 種算法的虛警比率始終為0，漏檢比率因其和故障檢測比率之和為1。由圖7 可知所提算法和基于卡爾曼濾波的故障檢測法分別自單星偽距故障超過3m 和8 m 后達(dá)到100%的故障檢測比率和故障識別比率。最小二乘殘差法在單星偽距故障為8m 時，故障檢測比率達(dá)到100%；在單星偽距故障為9m 時，故障識別比率達(dá)到100%。所提算法在故障小于2m 時3D 均方根誤差隨故障的增大而增大，當(dāng)故障超過3m 后，由于所提算法能100%地檢測并識別故障，3D 均方根誤差保持在0.69m。基于卡爾曼濾波的故障檢測法的3D 均方根誤差隨著故障增加至6m 而逐漸上升至3m 左右，在故障超過7m 后3D 均方根誤差下降至1m 左右。最小二乘殘差法的3D 均方根誤差在7m 故障時超過3m，隨后隨著故障的增加而降低。所提算法的3D 均方根誤差曲線始終在另外2 種算法的下方，因此，所提算法在不同的故障條件下均具有最高的定位精度。

圖7 靜態(tài)數(shù)據(jù)含不同大小的階躍故障時的故障檢測比率、故障識別比率和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 含單斜坡故障的靜態(tài)數(shù)據(jù)下的實(shí)驗(yàn)

若在觀測時刻100～200s（不含1 0 0 s，含2 0 0 s）內(nèi)往可見星G09 的偽距中添加以0.2m/s 的速率增大的斜坡故障，3 種算法的故障檢測量如圖8 所示（TSE為 斜坡故障的持續(xù)時間；VSE為斜坡故障的速率）。故障發(fā)生期間，3 種算法的故障檢測量均逐漸上升，并在150 s 之前超過閾值。所提算法開始檢測到斜坡故障的時間要略微早于另外2 種算法。

圖8 靜態(tài)數(shù)據(jù)在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 衛(wèi)星的故障識別量和其故障檢測量一樣于觀測時刻150 s 之前超過閾值，其偽距被識別出含有故障。圖10 說明基于卡爾曼濾波的故障檢測法大約在觀測時刻160 s后正確識別故障；識別出故障后，G09 的故障識別量經(jīng)過了短時間的快速上升。圖11 說明最小二乘殘差法僅在觀測時刻200s 前的大約2 0 s 時間內(nèi)正確識別了故障。

圖9 靜態(tài)數(shù)據(jù)在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ài)數(shù)據(jù)在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ài)數(shù)據(jù)在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 種算法的故障檢測指標(biāo)。最小二乘殘差法和基于卡爾曼濾波的故障檢測法的故障檢測比率均為57.8%。雖然所提算法的故障檢測比率為61.6%較這2 種算法沒有明顯地提高，但是本文所提算法的故障識別比率為51.4%，明顯地高于其他2 種算法20%和35.2%的故障識別率。

表5 靜態(tài)數(shù)據(jù)在T SE = 100 s，V SE = 0.2 m/s故障下的故障檢測指標(biāo)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 種算法的定位精度指標(biāo)。所提算法的3D 均方根誤差為0.712m，定位精度分別比2 種對照算法提高72.2%和64.4%。

表6 靜態(tài)數(shù)據(jù)在T SE = 100 s，V SE = 0.2 m/s故障下的定位精度指標(biāo)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 顯示了當(dāng)衛(wèi)星G09 的觀測時刻100～200 s（不含100s，含2 0 0 s）內(nèi)的斜坡故障速率從0 m/s增加到0.5m/s 時的實(shí)驗(yàn)結(jié)果。3 種算法的故障檢測比率比較接近，但所提算法的故障識別比率要明顯高于其他2 種算法。2 種對照算法的3D 均方根誤差先隨斜坡故障的速率的增加而增加，然后降低，最大值超過了2m，但是所提算法的3D 均方根誤差始終穩(wěn)定地保持小于1m。

圖12 靜態(tài)數(shù)據(jù)含不同大小的斜坡故障時的故障檢測比率、故障識別比率和3D均方根誤差Fig.12 The FDR, FIR and 3D accuracy (RMSE) for candidate algorithms with various ramp errors in the static mode

2.1.4 含多故障的靜態(tài)數(shù)據(jù)的實(shí)驗(yàn)

最小二乘殘差法基于單故障假設(shè)，每個歷元只能排除一個故障，因此，2.1.4 節(jié)和2.2.4 節(jié)中的實(shí)驗(yàn)只將基于卡爾曼濾波的故障檢測法作為對照組，驗(yàn)證所提算法對多故障的檢測性能。

觀測時刻150～200s（含150s 和200s）內(nèi)向衛(wèi)星G09 和G16 的偽距中添加5m 的階躍故障，2 種算法的故障檢測量如圖13 所示。故障發(fā)生期間，所提算法的故障檢測量始終高于閾值，但是基于卡爾曼濾波的故障檢測法未檢測出故障。

圖13 靜態(tài)數(shù)據(jù)中衛(wèi)星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 所示，故障發(fā)生期間，所提算法的衛(wèi)星G09 和G16 的故障識別量均超過了閾值。

圖14 靜態(tài)數(shù)據(jù)中衛(wèi)星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 顆衛(wèi)星的偽距同時存在的5m 的階躍故障。而基于卡爾曼濾波的故障檢測法的故障檢測檢測比率和識別比率均為0。

表7 靜態(tài)數(shù)據(jù)中衛(wèi)星G09 和G16 的偽距同時有5 m 的階躍故障下的故障檢測指標(biāo)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 種算法的定位精度指標(biāo)如表8 所示。此時基于卡爾曼濾波的故障檢測法的3D 均方根誤差為1.773m；所提算法的3 D 均方根誤差為0.7 0 4 m，定位精度大約均提高了60.3%

表8 靜態(tài)數(shù)據(jù)中衛(wèi)星G09 和G16 的偽距同時含有5 m的階躍故障下的定位精度指標(biāo)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 為向衛(wèi)星G10 和G12 的觀測時刻150～200s（含150 s 和200 s）內(nèi)的偽距添加的階躍故障逐漸增大后的故障檢測和定位情況。不同大小的多故障下，2 種算法的虛警比率均為0?；诳柭鼮V波的故障檢測法能100%地檢測超過8m 的多故障，能100%地識別超過10m 的多故障。所提算法分別在多故障達(dá)到3m 和4m 時能100%地檢測和識別多故障。多故障低于2 m 時，2 種算法的3D 均方根誤差比較接近，但是隨著多故障的進(jìn)一步增加，基于卡爾曼濾波的故障檢測法的3D 均方根誤差在多故障低于8m 時一直上升至高于2 m；而所提算法的3D 均方根始終大致地保持在低于1m 的水平。

圖15 靜態(tài)數(shù)據(jù)含不同大小的多故障時故障檢測比率、故障識別比率和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 基于動態(tài)數(shù)據(jù)的仿真實(shí)驗(yàn)

本文采用的動態(tài)數(shù)據(jù)來自2018 年1 月31 日于臺灣南投進(jìn)行的無人機(jī)動態(tài)觀測實(shí)驗(yàn)。無人機(jī)的飛行軌跡如圖16 所示，機(jī)載GNSS 接收機(jī)型號為Trimble BD9 8 2，采樣率為1 0 Hz。實(shí)驗(yàn)中使用的無人機(jī)的參考軌跡是通過機(jī)載VLP-16 Velodyne 激光雷達(dá)進(jìn)行近景攝影測量獲得的，具有cm 級的定位精度。

圖16 無人機(jī)的飛行軌跡Fig.16 Flight trajectory of the unmanned aerial vehicle (UAV)

所提算法對第1、第2 個歷元的數(shù)據(jù)利用最小二乘法進(jìn)行解算得到初始狀態(tài)向量，設(shè)置初始濾波參數(shù)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 原始動態(tài)數(shù)據(jù)下的實(shí)驗(yàn)

由圖17 可知，原始動態(tài)數(shù)據(jù)下，3 種算法整個過程故障檢測量均未超過故障檢測閾值，虛警比率都為0。

圖17 原始動態(tài)數(shù)據(jù)的故障檢測量Fig.17 Fault detection statistics based on original dynamic data

所提算法對原始動態(tài)數(shù)據(jù)進(jìn)行定位解算后各方向定位誤差的變化情況如圖18 所示。定位精度指標(biāo)如表9 所示，表9 中所提算法的水平均方根誤差和3D 均方根誤差分別為1.567m 和2.233m，對比最小二乘殘差法和基于卡爾曼濾波的故障檢測法，所提算法水平方向上的定位精度分別提高1.14%和0.57%，3D 方向上的定位精度分別提高33.4%和32%。

表9 原始動態(tài)數(shù)據(jù)的定位精度指標(biāo)Table 9 Positioning accuracy performance for candidate algorithm s based on original dynam ic data m

圖18 原始動態(tài)數(shù)據(jù)下所提算法的定位誤差Fig.18 Positioning error of the proposed algorithm based on original dynamic data

2.2.2 含單階躍故障的動態(tài)數(shù)據(jù)下的實(shí)驗(yàn)

往可見衛(wèi)星G10 在觀測時刻300～400s（含300s和400s）內(nèi)的偽距添加1 0 m 的階躍故障后，使用本文所提算法和2 種對照算法進(jìn)行故障檢測。

如圖19 所示，觀測時刻300～400 s 內(nèi)所提算法的故障檢測量均高于故障檢測閾值，其余時間內(nèi)低于閾值。對照算法中，基于卡爾曼濾波的故障檢測法在故障發(fā)生期間其故障檢測量均低于閾值，故障檢測比率為0，而最小二乘殘差法具有較高的故障檢測比率。如圖20 所示，所提算法的衛(wèi)星G10 的故障識別量在300～400 s 內(nèi)均超過相應(yīng)閾值，其他可見衛(wèi)星的故障識別量在此段時間內(nèi)均低于相應(yīng)閾值，因此，所提算法能100%正確識別故障。圖21表明最小二乘殘差法在一半左右的故障發(fā)生的歷元上未能正確識別故障。

圖19 動態(tài)數(shù)據(jù)含10 m階躍故障下的故障檢測量Fig.19 Fault detection statistics with a 10 m step error added on observations in the dynamic mode

圖20 動態(tài)數(shù)據(jù)含10 m階躍故障下所提算法的故障識別量Fig.20 Fault identification statistics with a 10 m step error added on observations in the dynamic mode

圖21 動態(tài)數(shù)據(jù)含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故障的動態(tài)數(shù)據(jù)時，故障檢測效果較好，故障檢測比率和故障識別比率均為100%，并且虛警比率為0。最小二乘殘差法的故障檢測比率雖然也達(dá)到了98.7%，但其故障識別比率卻不到60%；基于卡爾曼濾波的故障檢測法未能檢測到故障。

表10 動態(tài)數(shù)據(jù)含10 m 階躍故障下的故障檢測指標(biāo)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 種對照算法在處理含故障的動態(tài)數(shù)據(jù)時的定位精度指標(biāo)，此時所提算法定位結(jié)果的水平均方根誤差和3D 均方根誤差分別為1.592m 和3.049m，相較于最小二乘殘差法，水平定位精度和3D 定位精度分別提高2.63%和39.87%；相較于基于卡爾曼濾波的故障檢測法，水平定位精度和3D 定位精度分別提高和3.28%和52.33%。

表11 動態(tài)數(shù)據(jù)含10 m 階躍故障下的定位精度指標(biāo)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）內(nèi)的偽距依次添加 0,2,···,18 m的階躍故障，3 種算法的故障檢測效果如圖22 所示。由圖22 可知，為使故障檢測比率和故障識別比率均達(dá)到100%，所提算法和基于卡爾曼濾波的故障檢測法分別需要故障超過10m 和1 6m。最小二乘殘差法在故障超過10m 時故障檢測比率達(dá)到1 0 0%，在故障超過1 2 m時故障識別比率達(dá)到100%。所提算法的3D 均方根誤差在故障增加到2m 的過程中經(jīng)歷了較為緩慢的上升，隨后便基本穩(wěn)定在3m。但是基于卡爾曼濾波的故障檢測法的3D 均方根誤差在故障達(dá)到14m 前一直都在上升，最高點(diǎn)的3D 均方根誤差超過了7m，當(dāng)故障超過了1 6 m，其3D 均方根誤差下降至3.5m 左右。最小二乘殘差法的3D 均方根誤差在8m 故障時達(dá)到最大值6 m。故障從0 m 增加至18m 的過程中，所提算法的3D 均方根誤差始終是3 種算法中最低的，所提算法在不同的故障條件下具有更高的定位精度。

圖22 動態(tài)數(shù)據(jù)含不同大小的階躍故障時的故障檢測比率、故障識別比率、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 含單斜坡故障的動態(tài)數(shù)據(jù)下的實(shí)驗(yàn)

本節(jié)實(shí)驗(yàn)中向衛(wèi)星G10 的偽距添加從觀測時刻400～6 0 0 s（不含4 0 0 s，含6 0 0 s）的速率為0.2m/s的斜坡故障。

3 種算法對含斜坡故障的動態(tài)數(shù)據(jù)的故障檢測量如圖23 所示。所提算法的故障檢測量從400 s開始斜坡上升，并在450 s 左右超過閾值，所提算法開始檢測到斜坡故障。對照算法和所提算法有類似的故障檢測量的變化情況；開始檢測到斜坡故障的時間比較接近。

圖23 動態(tài)數(shù)據(jù)在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 所示，所提算法計算的衛(wèi)星G10 的故障識別量從400 s 開始斜坡上升，并大致在所提算法檢測到故障的同一時間超過閾值，此后一直到故障不再發(fā)生，所提算法均能正確地識別故障。如圖25所示為基于卡爾曼濾波的故障檢測法識別出故障的時間要晚于所提算法。最小二乘殘差法在斜坡故障增大到一定值后，多顆衛(wèi)星的故障識別量超過閾值，此時將他們中的最大值對應(yīng)的數(shù)據(jù)標(biāo)記為故障數(shù)據(jù)。最小二乘殘差法在故障發(fā)生期間的故障識別結(jié)果如圖26 所示。最小二乘殘差法在440 s左右首次正確識別故障，然后，一小段時間內(nèi)在正確識別故障和錯誤識別故障間波動，過后便能一直正確識別含斜坡故障的數(shù)據(jù)。

圖24 動態(tài)數(shù)據(jù)在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ài)數(shù)據(jù)在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ài)數(shù)據(jù)在T SE=200 s，V SE=0.2 m/s故障下最小二乘殘差法的故障識別結(jié)果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ài)數(shù)據(jù)在 T SE= 200 s，V SE = 0.2 m/s故障下的故障檢測指標(biāo)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 種算法的定位精度指標(biāo)，2 種對照算法的3D 均方根誤差分別為4.253m 和3.689m，然而所提算法的3D 均方根誤差為2.552m，相較于2 種對照算法分別降低40%和30.2%。

表13 動態(tài)數(shù)據(jù)在 T SE= 200 s，V SE = 0.2 m/s故障下定位精度指標(biāo)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 中，若向衛(wèi)星G10 的觀測時刻400～600 s（不含400s，含6 0 0 s）偽距中添加速率逐漸增加的斜坡故障，最小二乘殘差法的故障檢測比率一直略高于所提算法，但是所提算法的故障識別比率始終高于另外2 種對照算法。此外所提算法的定位精度也始終高于2 種對照算法。

圖27 動態(tài)數(shù)據(jù)含不同大小的斜坡故障時的故障檢測比率、故障識別比率和3D均方根誤差Fig.27 The FDR, FIR and 3D accuracy (RMSE) for candidate algorithms with various ramp errors in the dynamic mode

2.2.4 含多故障的動態(tài)數(shù)據(jù)下的實(shí)驗(yàn)

本節(jié)比較分析所提算法對多故障的檢測效果。如圖28 所示，當(dāng)向衛(wèi)星G10 和G12 在觀測時刻400～5 0 0 s（含4 0 0 s 和5 0 0 s）的偽距中注入1 0m的階躍故障時，基于卡爾曼濾波的故障檢測法僅在十分少的歷元檢測到了故障；而所提算法能100%地檢測到多故障。

圖28 動態(tài)數(shù)據(jù)中衛(wèi)星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 可知所提算法基本能正確識別多故障，只有在少量歷元，衛(wèi)星G12 的故障識別量沒有超過閾值，沒有被識別出來。但是在圖30 中，只有十分少量的歷元內(nèi)，衛(wèi)星G12 的故障識別量超過了閾值，即使是在這些歷元內(nèi)，多故障也不能被完全識別。

圖29 動態(tài)數(shù)據(jù)中衛(wèi)星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 動態(tài)數(shù)據(jù)中衛(wèi)星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%，顯著優(yōu)于基于卡爾曼濾波的故障檢測法。

表14 動態(tài)數(shù)據(jù)中衛(wèi)星G10 和G12 的偽距同時含有10 m階躍故障下的故障檢測指標(biāo)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 種算法的定位精度指標(biāo)。所提算法的3D 均方根誤差為2.645m，而基于卡爾曼濾波的故障檢測法的3D 均方根誤差為6.618m。在動態(tài)數(shù)據(jù)含有10m 的多故障條件下，所提算法的定位精度相較對照算法提高了60.03%。

表15 動態(tài)數(shù)據(jù)中衛(wèi)星G10 和G12 的偽距同時含有10 m 的故障下的定位精度指標(biāo)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

為了進(jìn)一步比較2 種算法的性能，圖31 給出了衛(wèi)星G10 和G12 的觀測時刻400～500s（含4 0 0 s和500 s）的偽距含有 0,2,···,16 m的階躍故障時，2 種算法的故障檢測比率、故障識別比率和3D均方根誤差。所提算法的故障檢測比率和故障識別比率分別在多故障超過8m 和1 2m 時達(dá)到了100%；基于卡爾曼濾波的故障檢測法的故障檢測比率和故障識別比率分別在多故障超過12m 和1 4m 時達(dá)到了100%。隨著多故障從0 增大到18m，由于所提算法的濾波采取了自適應(yīng)的噪聲方差策略并且有能檢測到較小的多故障，所提算法的3D 均方根誤差一直穩(wěn)定在2m 左右，但是基于卡爾曼濾波的故障檢測法在多故障達(dá)到了12m 時，其3D 均方根誤差也達(dá)到了峰值8m。

圖31 動態(tài)數(shù)據(jù)含不同大小的多故障時故障檢測比率、故障識別比率和3D均方根誤差Fig.31 The FDR, FIR and 3D accuracy RMSE for candidate algorithms with various step errors on two satellites

3 結(jié) 論

本文提出了基于自適應(yīng)噪聲方差的衛(wèi)星定位故障檢測算法，利用滑動窗內(nèi)的新息序列實(shí)時估計觀測噪聲方差矩陣，降低了設(shè)置觀測噪聲矩陣的主觀性，且能更好地適應(yīng)實(shí)際觀測噪聲隨觀測環(huán)境而變化。在自適應(yīng)噪聲方差的基礎(chǔ)上，所提算法的故障檢測量和故障識別量能更好地滿足對應(yīng)的統(tǒng)計學(xué)概率分布，具有較高的故障檢測比率和識別比率，特別是對于微小故障；而且完成故障檢測后，自適應(yīng)噪聲方差的策略也優(yōu)化了卡爾曼增益矩陣，一定程度上提高了濾波的精度。實(shí)驗(yàn)結(jié)果表明：

1）若給定的靜態(tài)數(shù)據(jù)中單星偽距的階躍故障逐漸增大，所提算法在階躍故障超過3m 時，其故障檢測比率和故障識別比率同時達(dá)到100%；而基于卡爾曼濾波的故障檢測法在階躍故障超過8m時，其故障檢測比率和故障識別比率才均達(dá)到100%；最小二乘殘差法分別在階躍故障超過8m 和超過9 m時，達(dá)到了100%的故障檢測比率和故障識別比率。

2）若給定的靜態(tài)數(shù)據(jù)中單星偽距的斜坡故障的速率逐漸增大，所提算法和對照算法的故障檢測比率的變化情況接近，但所提算法的故障識別比率明顯高于對照算法。

3）若給定的靜態(tài)數(shù)據(jù)所含的多故障逐漸增大，所提算法分別在多故障超過3m 和超過4 m 時100%地檢測和識別故障；而基于卡爾曼濾波的故障檢測法分別在多故障超過8m 和超過1 0 m 時100%地檢測和識別故障。

4）若給定的動態(tài)數(shù)據(jù)中單星偽距的階躍故障逐漸增加，要使故障檢測比率和識別比率同時達(dá)到100%，所提算法需要故障超過10m；而基于卡爾曼濾波的故障檢測法需要故障超過16m；雖然最小二乘殘差法在故障超過10m 時，其故障檢測比率也達(dá)到了100%，但故障超過12m 時其故障識別比率才能達(dá)到100%。

5）若給定的動態(tài)數(shù)據(jù)中單星偽距的斜坡故障的速率逐漸增加，雖然最小二乘殘差法的故障檢測比率略高于所提算法，但是所提算法的故障識別比率一直明顯高于另外2 種對照算法。

6）若給定的動態(tài)數(shù)據(jù)所含的多故障逐漸增大，所提算法和基于卡爾曼濾波的故障檢測法分別在故障超過8m 和1 2m 時，達(dá)到100%的故障檢測比率；這2 種算法分別在故障超過12m 和14m 時，達(dá)到了100%的故障識別比率。

在本文所有故障情景下，由于所提算法采取了根據(jù)滑動窗中的歷史新息估計觀測噪聲方差的策略，提高了基于卡爾曼濾波的故障檢測法的故障檢測比率和識別比率，同時優(yōu)化了測量更新中各顆衛(wèi)星的權(quán)重，所以應(yīng)用所提算法進(jìn)行質(zhì)量控制后的3D均方根誤差基本不隨故障的增大而發(fā)生變化，其定位精度一直高于2 種對照算法。