鐘倫瓏，劉永玉，李雪艷

1.中國(guó)民航大學(xué)智能信號(hào)與圖像處理天津市重點(diǎn)實(shí)驗(yàn)室，天津 300300

2.中國(guó)民航大學(xué)，天津 300300

全球?qū)Ш叫l(wèi)星系統(tǒng)（GNSS）為航空領(lǐng)域提供覆蓋全球的高性能導(dǎo)航信息。國(guó)際民航組織（ICAO）將其選作新一代區(qū)域?qū)Ш较到y(tǒng)的主用系統(tǒng)，規(guī)定使用民用GNSS 信號(hào)進(jìn)行導(dǎo)航。但是，由于民用GNSS 信號(hào)的碼字信息和信號(hào)結(jié)構(gòu)是公開(kāi)的，GNSS 導(dǎo)航容易受到蓄意的欺騙式干擾影響[1-2]。及時(shí)、有效的欺騙式干擾檢測(cè)技術(shù)對(duì)民航飛行安全至關(guān)重要。根據(jù)檢測(cè)數(shù)據(jù)的來(lái)源不同，機(jī)載欺騙式干擾檢測(cè)技術(shù)可分為基于信號(hào)處理的檢測(cè)技術(shù)[3-4]和基于信息解算[5-6]的檢測(cè)技術(shù)。與基于信號(hào)處理的檢測(cè)技術(shù)相比，基于信息解算的檢測(cè)技術(shù)無(wú)須改變接收機(jī)的結(jié)構(gòu)[7]，可融合其他無(wú)干擾的傳感器信息，提高算法檢測(cè)性能，且可嵌入機(jī)載組合導(dǎo)航算法中，便于在區(qū)域?qū)Ш接?jì)算機(jī)中實(shí)現(xiàn)。

典型的基于信息解算的欺騙式干擾檢測(cè)技術(shù)對(duì)衛(wèi)星導(dǎo)航和其他傳感器信息進(jìn)行一致性檢測(cè)[8]，如基于GNSS和慣性導(dǎo)航系統(tǒng)（INS）緊組合導(dǎo)航的欺騙式干擾檢測(cè)技術(shù)[9]，利用INS 不受欺騙式干擾影響的特點(diǎn)，使用卡爾曼濾波等方式融合GNSS 測(cè)量信息和INS 導(dǎo)航信息，再結(jié)合適當(dāng)?shù)臋z測(cè)方法判斷GNSS 測(cè)量信息和INS 導(dǎo)航信息的一致性，實(shí)現(xiàn)欺騙式干擾檢測(cè)。傳統(tǒng)的新息序列卡方檢測(cè)方法直接以組合導(dǎo)航濾波器新息構(gòu)建卡方檢測(cè)統(tǒng)計(jì)量，但卡方檢測(cè)容易擴(kuò)大檢測(cè)區(qū)域，導(dǎo)致對(duì)微小欺騙量檢測(cè)延時(shí)較大，漏警率較高。為提高算法性能，Bhatti[10]提出將新息的變化率作為檢測(cè)統(tǒng)計(jì)量，張超[6]提出將抗差估計(jì)與新息卡方檢測(cè)算法相結(jié)合。兩種方法重新構(gòu)建基于新息的欺騙檢測(cè)統(tǒng)計(jì)量，提高了對(duì)欺騙的檢測(cè)效果，但仍不滿足ICAO 對(duì)巡航階段的告警時(shí)間要求。

KL散度是基于數(shù)理統(tǒng)計(jì)方法的重要概念之一，能夠有效描述樣本數(shù)據(jù)可能存在的兩個(gè)概率密度分布之間的差異。曹玉蘋[11]等提出基于KL 散度的無(wú)跡卡爾曼濾波的過(guò)程故障檢測(cè)方法，能夠及時(shí)檢測(cè)出故障的發(fā)生?；诠收蠙z測(cè)思想[12]，針對(duì)GNSS 欺騙式干擾威脅民航飛行安全問(wèn)題，本文引入KL散度算法，基于緊組合導(dǎo)航濾波過(guò)程判斷GNSS是否受到欺騙式干擾。首先使用擴(kuò)展卡爾曼濾波器進(jìn)行GNSS的觀測(cè)數(shù)據(jù)和INS的導(dǎo)航數(shù)據(jù)融合。然后在緊組合系統(tǒng)濾波器輸出新息的概率密度分布基礎(chǔ)上，構(gòu)建基于KL 散度算法的檢測(cè)統(tǒng)計(jì)量，進(jìn)行GNSS 欺騙式干擾檢測(cè)。最后進(jìn)行仿真試驗(yàn)驗(yàn)證本文所提方法的有效性。

1 基于偽距測(cè)量的欺騙式干擾模型

GNSS欺騙式干擾的根本目的是對(duì)用戶接收機(jī)造成惡意攻擊，使用戶接收機(jī)捕獲到虛假的干擾信號(hào)，形成錯(cuò)誤的測(cè)量偽距，導(dǎo)致接收機(jī)定位結(jié)果偏離真實(shí)的用戶目標(biāo)位置[13]。根據(jù)GNSS 基本原理[14]，用戶接收機(jī)R與第j顆可見(jiàn)衛(wèi)星的偽距測(cè)量值ρj為

式中，rj為用戶接收機(jī)R到第j顆可見(jiàn)衛(wèi)星的真實(shí)距離；δtu為GNSS接收機(jī)時(shí)鐘誤差的等效距離；vρ是接收機(jī)內(nèi)部噪聲與GNSS 測(cè)量噪聲總和，vρ～N(0,)，σG為偽距測(cè)量噪聲標(biāo)準(zhǔn)差。

欺騙式干擾使用戶接收機(jī)定位結(jié)果偏離真實(shí)位置的關(guān)鍵是在偽距測(cè)量值上引入了附加的偽距欺騙量。當(dāng)存在欺騙式干擾時(shí)，受附加的偽距欺騙量的影響，第j顆可見(jiàn)衛(wèi)星的偽距測(cè)量值變?yōu)?/p>

式中，Δρj為欺騙式干擾引入的附加偽距。當(dāng)Δρj=a（a為常值）時(shí)，存在突變式欺騙干擾；當(dāng)Δρj=v(t-ts)時(shí)，存在緩變式欺騙干擾，v為偽距欺騙速率，t為導(dǎo)航時(shí)間，ts為欺騙發(fā)生時(shí)刻。

2 基于緊組合導(dǎo)航的欺騙檢測(cè)系統(tǒng)模型

GNSS和INS緊組合導(dǎo)航是基于偽距、偽距率的組合導(dǎo)航方式[15]。相對(duì)于其他組合導(dǎo)航方式，緊組合導(dǎo)航有更高的精度和更強(qiáng)的抗干擾能力。在GNSS/INS緊組合導(dǎo)航系統(tǒng)中嵌入欺騙式干擾檢測(cè)算法，利用卡爾曼濾波器新息的概率密度分布，構(gòu)建基于KL散度算法的檢測(cè)統(tǒng)計(jì)量，實(shí)現(xiàn)欺騙式干擾檢測(cè)功能。檢測(cè)無(wú)欺騙時(shí)系統(tǒng)校正INS進(jìn)行組合導(dǎo)航輸出，否則系統(tǒng)進(jìn)行欺騙異常告警。設(shè)計(jì)的基于緊組合導(dǎo)航的欺騙檢測(cè)系統(tǒng)模型如圖1所示。

2.1 緊組合導(dǎo)航系統(tǒng)狀態(tài)方程

GNSS/INS 緊組合導(dǎo)航系統(tǒng)的狀態(tài)量由INS 的誤差狀態(tài)和GNSS 的誤差狀態(tài)構(gòu)成，包括18 維的INS 誤差狀態(tài)和二維的GNSS 誤差狀態(tài)。由緊組合導(dǎo)航系統(tǒng)理論可知[16]，系統(tǒng)狀態(tài)方程為

式中，F(xiàn)為狀態(tài)轉(zhuǎn)移矩陣；G為系統(tǒng)噪聲矩陣；W為系統(tǒng)噪聲矢量；X為系統(tǒng)狀態(tài)量。

圖1 中，卡爾曼濾波器新息為緊組合導(dǎo)航系統(tǒng)實(shí)際量測(cè)值與預(yù)測(cè)值之差。離散化式（3）所示的狀態(tài)方程和式（5）所示的觀測(cè)方程，可得

圖1 基于緊組合導(dǎo)航的欺騙檢測(cè)系統(tǒng)模型Fig.1 Diagram of tightly-coupled navigation based spoofing interference detection system

式中，φ為東北天方向的姿態(tài)角誤差；δv為東北天方向的速度誤差；δ為緯經(jīng)高方向的位置誤差；ε為INS陀螺儀誤差；?為INS 加速度計(jì)零偏；δtu為GNSS 接收機(jī)時(shí)鐘誤差的等效距離；δtru為GNSS接收機(jī)時(shí)鐘頻率誤差的等效距離率。

2.2 緊組合導(dǎo)航系統(tǒng)量測(cè)方程

GNSS/INS緊組合導(dǎo)航系統(tǒng)量測(cè)的實(shí)際觀測(cè)量由偽距、偽距率量測(cè)值構(gòu)成。偽距、偽距率量測(cè)值為INS 等效的偽距、偽距率與GNSS 接收機(jī)觀測(cè)的偽距、偽距率之差。INS等效的偽距、偽距率由INS推算的載體位置、速度相應(yīng)與星歷中衛(wèi)星的位置、速度計(jì)算得到。

在地心地固坐標(biāo)系下[16]，GNSS/INS緊組合導(dǎo)航系統(tǒng)量測(cè)方程為

式中，H為量測(cè)矩陣；V為量測(cè)噪聲矢量；X為系統(tǒng)狀態(tài)量；Z為系統(tǒng)觀測(cè)量，由n顆可見(jiàn)衛(wèi)星的偽距、偽距率量測(cè)值構(gòu)成，即

式中，ρI,j、ρ?I,j(j=1，…，n)為INS 與第j顆可見(jiàn)衛(wèi)星等效的偽距、偽距率；ρG,j、ρ?G,j為GNSS 接收機(jī)觀測(cè)的第j顆可見(jiàn)衛(wèi)星星歷中原始的偽距、偽距率。

3 基于KL散度的欺騙式干擾檢測(cè)方法

式中，k表示當(dāng)前時(shí)刻；Zk為系統(tǒng)觀測(cè)量；Xk為系統(tǒng)狀態(tài)量；Hk為系統(tǒng)觀測(cè)方程；Φk為狀態(tài)轉(zhuǎn)移矩陣；Γk為噪聲分配矩陣；Wk-1與Vk是互不相關(guān)的高斯白噪聲。

新息定義為

式中，δZk為濾波器新息矢量；Zk為系統(tǒng)觀測(cè)量；Zk|k-1為系統(tǒng)預(yù)測(cè)值；Xk|k-1為狀態(tài)一步預(yù)測(cè)值；Hk為系統(tǒng)觀測(cè)方程。

緊組合導(dǎo)航系統(tǒng)的量測(cè)噪聲和系統(tǒng)噪聲均為互不相關(guān)的高斯白噪聲，因此衛(wèi)星未受欺騙時(shí)新息服從零均值高斯分布，衛(wèi)星受欺騙時(shí)引入的附加偽距會(huì)改變新息值的均值，但新息值依然服從高斯分布。

3.1 傳統(tǒng)的新息序列卡方檢測(cè)方法

傳統(tǒng)的新息序列卡方檢測(cè)方法利用欺騙式干擾后新息的變化來(lái)檢測(cè)干擾。若當(dāng)前時(shí)刻GNSS 受到欺騙式干擾，則GNSS 接收機(jī)觀測(cè)的衛(wèi)星偽距增大，以致濾波器實(shí)際量測(cè)值增大。又因?yàn)V波器的預(yù)測(cè)值未受當(dāng)前時(shí)刻欺騙式干擾的影響，從而濾波器新息值增大并反映出欺騙信息。檢測(cè)統(tǒng)計(jì)量定義為

式中，St(t= 1,2,…,k)為各濾波時(shí)刻新息矢量的協(xié)方差矩陣。

根據(jù)卡方分布分位數(shù)的計(jì)算原理，可由虛警率Pf確定檢測(cè)閾值T。

3.2 基于KL散度的欺騙式干擾檢測(cè)方法

通過(guò)上述分析可知，衛(wèi)星受欺騙前后新息的概率密度分布是不同的。利用新息可能存在的兩種概率密度分布，本文基于KL散度構(gòu)建欺騙檢測(cè)統(tǒng)計(jì)量，將欺騙誤差信息反映為兩種概率密度分布之間的“距離”大小。與上述的傳統(tǒng)新息序列卡方檢測(cè)方法相比，此法可避免卡方檢測(cè)統(tǒng)計(jì)量容易擴(kuò)大檢測(cè)區(qū)域，導(dǎo)致對(duì)緩變式欺騙和微小的突變式欺騙檢測(cè)易發(fā)生漏檢的問(wèn)題。

KL 散度是一種衡量樣本數(shù)據(jù)可能存在的兩種概率密度分布之間相似度大小的方法，被廣泛應(yīng)用在基于數(shù)理統(tǒng)計(jì)研究方法的信息論領(lǐng)域[17]。假設(shè)x存在兩種可能的概率密度函數(shù)p(x)和q(x)，則KL散度為

KL 散度描述了概率密度函數(shù)p(x)和q(x)之間的“距離”，滿足非負(fù)特性和非對(duì)稱性，當(dāng)且僅當(dāng)描述的兩個(gè)概率密度函數(shù)完全相同時(shí)，KL散度取到最小值0。

對(duì)第j顆可見(jiàn)星構(gòu)建基于KL散度算法的欺騙檢測(cè)統(tǒng)計(jì)量。式（8）中，當(dāng)前時(shí)刻新息矢量δZk的第j行表示第j顆衛(wèi)星的新息值。以歷史時(shí)刻到當(dāng)前時(shí)刻第j顆衛(wèi)星的新息序列{δZt,j|t= 1,2,…,k}作 為KL 散 度 算 法 檢 測(cè) 樣 本xi(i=1,2,…,k)，當(dāng)前時(shí)刻的樣本均值和方差定義為

基于二元假設(shè)檢驗(yàn)[18]，定義原假設(shè)H0:GNSS 未受欺騙，= 0；備擇假設(shè)H1:GNSS受欺騙≠0。在二元假設(shè)條件下，當(dāng)前時(shí)刻觀測(cè)樣本的概率密度函數(shù)可表示為

對(duì)KL散度進(jìn)行簡(jiǎn)化[11]，定義檢測(cè)統(tǒng)計(jì)量為

將式（13）、式（14）代入式（15），化簡(jiǎn)得到檢測(cè)統(tǒng)計(jì)量

式中，p(xk|H0)為當(dāng)前時(shí)刻GNSS未受到欺騙式干擾時(shí)新息的概率密度函數(shù)；p(xk|H1)為當(dāng)前時(shí)刻GNSS受到欺騙式干擾時(shí)新息的概率密度函數(shù)；Dk≥0，當(dāng)且僅當(dāng)p(xk|H1)=p(xk|H0)時(shí)，Dk= 0。

在欺騙發(fā)生后的告警時(shí)間內(nèi)，若Dk＞Td，則判斷第j顆衛(wèi)星受到了欺騙式干擾，否則判斷第j顆衛(wèi)星未受到欺騙式干擾。檢測(cè)門限Td由虛警率Pf和漏警率Pm計(jì)算得到。

4 仿真結(jié)果與分析

4.1 試驗(yàn)設(shè)置和評(píng)價(jià)指標(biāo)

模擬飛機(jī)巡航階段中衛(wèi)星受到不同的欺騙式干擾，對(duì)比本文KL 散度檢測(cè)方法與傳統(tǒng)新息序列卡方檢測(cè)方法的檢測(cè)性能。設(shè)置三組仿真試驗(yàn):試驗(yàn)1 驗(yàn)證在突變式欺騙條件下兩種方法的檢測(cè)時(shí)間；試驗(yàn)2 驗(yàn)證在緩變式欺騙條件下兩種方法的檢測(cè)時(shí)間；試驗(yàn)3 驗(yàn)證在不同的欺騙方式下兩種方法的欺騙檢測(cè)率。按照ICAO 附件10 標(biāo)準(zhǔn)規(guī)定[19]，仿真實(shí)驗(yàn)相關(guān)參數(shù)設(shè)置見(jiàn)表1。

表1 仿真參數(shù)Table 1 Simulation parameters

由表1的虛警率和漏警率計(jì)算式（10）和式（17）的欺騙檢測(cè)閾值。仿真中設(shè)置仰角截止角為30o，保持可見(jiàn)星數(shù)為5顆，對(duì)單星施加欺騙。

試驗(yàn)1和試驗(yàn)2中分別施加突變欺騙量和緩變欺騙量，具體的欺騙量設(shè)置見(jiàn)表2。

表2 欺騙量設(shè)置Table 2 Settings of spoofing amount

由于機(jī)載組合導(dǎo)航的卡爾曼濾波估計(jì)技術(shù)是一種最優(yōu)估計(jì)，很小的欺騙量對(duì)組合導(dǎo)航結(jié)果影響有限，因此設(shè)置的最小突變欺騙量為10m，最小緩變欺騙量為0.3m/s。

試驗(yàn)3 中不同欺騙方式下對(duì)每個(gè)偽距欺騙量進(jìn)行100次蒙特卡羅試驗(yàn)，計(jì)算欺騙檢測(cè)率P

式中，ts為欺騙發(fā)生時(shí)刻；Δt為欺騙告警時(shí)間；Num[Dk＞Td,ts≤t≤ts+Δt]為成功檢測(cè)的試驗(yàn)次數(shù)；N為蒙特卡羅試驗(yàn)次數(shù)。

參照ICAO 附件10 對(duì)巡航階段的完好性要求，告警時(shí)間Δt取最大值30s，且告警時(shí)間內(nèi)算法欺騙檢測(cè)的漏警率不能超過(guò)0.001，即在受到欺騙式干擾30s 內(nèi)檢測(cè)算法的欺騙檢測(cè)率需超過(guò)99.9%。

4.2 仿真結(jié)果分析

（1）試驗(yàn)1

分別施加10m、25m 和55m 的突變欺騙量，KL 散度檢測(cè)方法的檢測(cè)統(tǒng)計(jì)量隨時(shí)間的變化如圖2所示。

圖2 施加突變欺騙量時(shí)KL散度方法檢測(cè)時(shí)間Fig.2 Detection time of KL divergence method when applying abrupt spoofing amounts

由圖2（a）可以看出，在第200s施加最小10m的突變欺騙量時(shí)，KL 散度檢測(cè)方法的檢測(cè)統(tǒng)計(jì)量在第213.57s 超過(guò)檢測(cè)閾值，滿足ICAO對(duì)巡航階段的30s告警時(shí)間要求。對(duì)比圖2（a）、圖2（b）和圖2（c）可知，KL散度檢測(cè)方法的檢測(cè)統(tǒng)計(jì)量超過(guò)檢測(cè)閾值的時(shí)間隨著突變欺騙量的增大而逐漸減小，即算法的檢測(cè)速度隨著突變欺騙量的增大而加快。

對(duì)比KL 散度檢測(cè)方法與傳統(tǒng)新息序列卡方檢測(cè)方法對(duì)不同突變欺騙量的檢測(cè)時(shí)間，見(jiàn)表3。

表3 不同突變欺騙量的檢測(cè)時(shí)間對(duì)比Table 3 Comparison between detection time under different abrupt spoofing amounts

從表3 中可以看出，在相同的突變欺騙量條件下，KL散度檢測(cè)方法相較于傳統(tǒng)新息序列卡方檢測(cè)方法有效地縮短了欺騙檢測(cè)時(shí)間。特別是在10m的微小突變欺騙量情況下，傳統(tǒng)新息序列卡方檢測(cè)方法不能檢測(cè)出欺騙干擾，而KL散度檢測(cè)方法能以較快的速度檢測(cè)出來(lái)。

（2）試驗(yàn)2

分別施加0.3m/s、0.5m/s 和1m/s 的緩變欺騙量，KL 散度檢測(cè)方法的檢測(cè)統(tǒng)計(jì)量隨時(shí)間的變化如圖3所示。

由圖3（a）可看出，在第200s 施加最小0.3m/s 的緩變欺騙量，KL 散度檢測(cè)方法的檢測(cè)統(tǒng)計(jì)量在第225.86s 超過(guò)檢測(cè)閾值。雖然欺騙檢測(cè)延時(shí)較長(zhǎng)，但仍滿足ICAO 對(duì)巡航階段的30s告警時(shí)間要求。對(duì)比圖3（a）、圖3（b）和圖3（c），算法的檢測(cè)速度隨著緩變欺騙速率的增大而加快。

圖3 施加緩變欺騙量時(shí)KL散度方法檢測(cè)時(shí)間Fig.3 Detection time of KL divergence method when applying slowly varying amounts of spoofing

對(duì)比KL 散度檢測(cè)方法與新息序列卡方檢測(cè)方法對(duì)不同緩變欺騙量的檢測(cè)時(shí)間，見(jiàn)表4。

從表4中可以看出，在不同的緩變欺騙量條件下，傳統(tǒng)新息序列卡方檢測(cè)方法的檢測(cè)時(shí)間均不滿足ICAO對(duì)巡航階段的30s 告警時(shí)間要求。而KL 散度檢測(cè)方法相對(duì)提高了算法對(duì)緩變欺騙量的檢測(cè)敏感度，加快了檢測(cè)速度。

表4 不同緩變欺騙量的檢測(cè)時(shí)間對(duì)比Table 4 Comparison of detection time under different slowly varying amounts of spoofing

（3）試驗(yàn)3

根據(jù)式（19）計(jì)算每個(gè)偽距欺騙量下的欺騙檢測(cè)率。如圖4所示，以5m為間隔施加突變欺騙量，對(duì)比KL散度檢測(cè)方法和新息序列卡方檢測(cè)方法的欺騙檢測(cè)率。新息序列卡方檢測(cè)方法的欺騙檢測(cè)率在突變欺騙量為25m 時(shí)收斂到100%。而KL散度檢測(cè)方法的欺騙檢測(cè)率在欺騙量為10m時(shí)即可收斂到100%，相對(duì)降低了算法對(duì)小于25m的微小欺騙量的檢測(cè)漏警率。

圖4 突變欺騙量條件下兩種方法的檢測(cè)率Fig.4 Detection rates of two methods under the condition of abrupt spoofing amounts

如圖5 所示，以0.1m/s 為間隔施加緩變欺騙量，對(duì)比KL散度檢測(cè)方法和新息序列卡方檢測(cè)方法的欺騙檢測(cè)率。在施加0.3m/s到1.4m/s的緩變欺騙量時(shí)，KL散度欺騙檢測(cè)方法的檢測(cè)率明顯優(yōu)于新息序列卡方檢測(cè)方法。當(dāng)緩變欺騙量達(dá)到0.3m/s時(shí)，本文KL散度檢測(cè)方法的欺騙檢測(cè)率即可收斂到100%。而新息序列卡方檢測(cè)方法的欺騙檢測(cè)率在緩變欺騙量達(dá)到1.4m/s時(shí)才收斂到100%。相比之下，本文方法有效地解決了傳統(tǒng)方法對(duì)緩變式欺騙干擾檢測(cè)能力差的問(wèn)題。

圖5 緩變欺騙量條件下兩種方法的檢測(cè)率Fig.5 Detection rates of the two methods under the condition of slowly varying amounts of spoofing

由以上三組仿真試驗(yàn)可知，本文方法相較于傳統(tǒng)方法對(duì)欺騙式干擾的檢測(cè)敏感度更高，檢測(cè)性能更好。對(duì)于突變式欺騙干擾，本文方法可實(shí)現(xiàn)10m 以上的突變欺騙量的檢測(cè)。對(duì)于緩變式欺騙干擾，本文方法可實(shí)現(xiàn)0.3m/s 以上的緩變欺騙量的檢測(cè)。

5 結(jié)束語(yǔ)

本文基于航空運(yùn)輸飛機(jī)的機(jī)載電子系統(tǒng)體制，利用機(jī)載GNSS 和INS 導(dǎo)航數(shù)據(jù)，結(jié)合KL 散度理論，提出一種GNSS 欺騙式干擾檢測(cè)方法。試驗(yàn)結(jié)果表明，相較于傳統(tǒng)方法，本文方法提高了對(duì)緩變式欺騙干擾和微小突變式欺騙干擾的檢測(cè)能力。同時(shí)，本文方法以軟件形式集中進(jìn)行信息處理，在不增加成本的前提下提高了算法的欺騙檢測(cè)率，降低了檢測(cè)漏警率。