基于改進(jìn)聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)的AUV定位算法

李鑫濱 馬迎港 閆 磊 韓 松

（1.燕山大學(xué)，河北省工業(yè)計(jì)算機(jī)控制工程重點(diǎn)實(shí)驗(yàn)室，河北秦皇島 066004；2.東北大學(xué)秦皇島分校計(jì)算機(jī)與通信工程學(xué)院，河北秦皇島 066004）

1 引言

海洋中蘊(yùn)含著大量的礦產(chǎn)資源、海洋生物資源等，海洋探索與開(kāi)發(fā)成為國(guó)際各領(lǐng)域炙手可熱的研究方向［1］，水下自主機(jī)器人因?yàn)槠湓趫?zhí)行水下任務(wù)時(shí)的靈活性以及其經(jīng)濟(jì)性，已廣泛應(yīng)用于許多的海洋應(yīng)用之中，如水下圖像數(shù)據(jù)采集［2］、水下軍事檢測(cè)［3］、海洋突發(fā)環(huán)境應(yīng)急［4］等，但水下自主航行器（Autonomous Underwater Vehicle，AUV）在這些領(lǐng)域中的應(yīng)用必須依靠準(zhǔn)確的定位才能有效。水下無(wú)線傳感器網(wǎng)絡(luò)（Underwater Wireless Sensor Networks，UWSN）利用聲學(xué)信號(hào)作為通信介質(zhì)，使用基于距離的定位方法，為AUV的水下環(huán)境的定位通信提供了一個(gè)良好的平臺(tái)。在基于距離的定位方法之中，基于到達(dá)時(shí)間（Time of Arrival，TOA）由于有更好的魯棒性［5］，在UWSN中已經(jīng)得到了廣泛使用。一階線性化的拓展卡爾曼濾波（Extended Kalman Filter，EKF）通過(guò)使用TOA 的方法獲取信標(biāo)節(jié)點(diǎn)和AUV 之間的距離測(cè)量并實(shí)現(xiàn)AUV 的定位，由于其簡(jiǎn)單、計(jì)算復(fù)雜度低等優(yōu)點(diǎn)［6］，是目前估計(jì)AUV 位置最常用的方法。

傳統(tǒng)的EKF 算法的應(yīng)用是在已知過(guò)程噪聲和測(cè)量噪聲協(xié)方差的前提下進(jìn)行的，然而在實(shí)際的水下環(huán)境中由于外部干擾和環(huán)境變化等因素［7］，導(dǎo)致準(zhǔn)確的噪聲協(xié)方差矩陣很難確定，而且可能是時(shí)變的，使用錯(cuò)誤或不準(zhǔn)確的噪聲協(xié)方差矩陣會(huì)導(dǎo)致大量的累積估計(jì)誤差，甚至濾波發(fā)散。針對(duì)這個(gè)問(wèn)題，文獻(xiàn)［8］提出了自適應(yīng)調(diào)節(jié)先驗(yàn)誤差協(xié)方差矩陣的強(qiáng)跟蹤濾波器，文獻(xiàn)［9］提出了用于未知測(cè)量噪聲協(xié)方差矩陣的自適應(yīng)Sage-Husa 濾波算法，實(shí)時(shí)估計(jì)和修正噪聲統(tǒng)計(jì)量來(lái)自適應(yīng)跟蹤過(guò)程中的噪聲變化，然而當(dāng)同時(shí)估計(jì)過(guò)程噪聲協(xié)方差矩陣和測(cè)量噪聲協(xié)方差矩陣時(shí)，估計(jì)值與真實(shí)值之間存在著一個(gè)常值誤差，即不可能同時(shí)估計(jì)出真實(shí)的過(guò)程噪聲協(xié)方差矩陣和測(cè)量噪聲協(xié)方差矩陣。在文獻(xiàn)［10］中提出一種迭代拓展卡爾曼濾波AEKF，使用期望最大化算法（Expectation-Maximum，EM）同時(shí)自適應(yīng)估計(jì)預(yù)測(cè)誤差協(xié)方差矩陣和測(cè)量噪聲協(xié)方差矩陣，解決了上面的問(wèn)題，期望最大化算法［11］是一種通過(guò)迭代進(jìn)行極大似然估計(jì)的優(yōu)化算法，在對(duì)模型的參數(shù)進(jìn)行迭代當(dāng)中其似然函數(shù)值一般是遞增的，但EM 算法在優(yōu)化過(guò)程中的重要前提是需要合理的模型參數(shù)迭代初值，如果所設(shè)置的迭代參數(shù)初值與現(xiàn)實(shí)所對(duì)應(yīng)的參數(shù)差距很大，則算法很有可能迭代出錯(cuò)誤的模型參數(shù)，陷入到局部最優(yōu)解中［12］。所以，該算法在應(yīng)對(duì)水下突發(fā)意外環(huán)境中的定位時(shí)的魯棒性會(huì)大大下降。

值得注意的是，水下節(jié)點(diǎn)通常部署在無(wú)人值守的開(kāi)放聲環(huán)境中，這使外來(lái)的惡意節(jié)點(diǎn)可通過(guò)偽裝成為合法的鄰居信標(biāo)節(jié)點(diǎn)，然后在信標(biāo)節(jié)點(diǎn)與AUV的距離測(cè)量數(shù)據(jù)傳輸中，惡意節(jié)點(diǎn)在攔截消息后通過(guò)延遲響應(yīng)［13］，可以很容易的在傳輸數(shù)據(jù)包中的時(shí)間戳和在測(cè)距數(shù)據(jù)中注入惡意噪聲，使得接收節(jié)點(diǎn)得到錯(cuò)誤的距離測(cè)量，危及無(wú)線傳感器網(wǎng)絡(luò)中待定位的AUV位置準(zhǔn)確估計(jì)。因此，在水下定位系統(tǒng)中進(jìn)行惡意噪聲檢測(cè)以保證整個(gè)網(wǎng)絡(luò)的可用性［14］是十分有必要的。目前，研究人員已經(jīng)提出了一些惡意噪聲檢測(cè)的方法，文獻(xiàn)［15］中基于估計(jì)水下聲學(xué)信號(hào)的平均傳播速度來(lái)判斷和刪除惡意信標(biāo)節(jié)點(diǎn)信息。文獻(xiàn)［16］基于Cayley-Menger 行列式和l1最小化解碼算法的幾何約束，構(gòu)造了一個(gè)距離測(cè)量的離群值檢測(cè)模型，通過(guò)丟棄異常值、只使用可靠的數(shù)據(jù)來(lái)估計(jì)未知節(jié)點(diǎn)的位置。文獻(xiàn)［17］使用基于密度的空間聚類來(lái)識(shí)別惡意噪聲，但基于密度的識(shí)別方法存在參數(shù)設(shè)置敏感的問(wèn)題，如果惡意噪聲的攻擊的頻率比較高時(shí)，產(chǎn)生的測(cè)量點(diǎn)對(duì)應(yīng)的密度就會(huì)比較分散，此時(shí)針對(duì)識(shí)別惡意噪聲所設(shè)置閾值就會(huì)比較困難。此外，在篩選出未受惡意噪聲污染的測(cè)量點(diǎn)后，還存在需要對(duì)這些測(cè)量點(diǎn)進(jìn)行數(shù)據(jù)融合的問(wèn)題，傳統(tǒng)的聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)（Joint Probabilistic Data Association，JPDA）算法是處理一批測(cè)量點(diǎn)并對(duì)目標(biāo)進(jìn)行定位更新的有效算法［18］，通過(guò)對(duì)落入其以預(yù)測(cè)點(diǎn)為中心的關(guān)聯(lián)門的測(cè)量點(diǎn)進(jìn)行數(shù)據(jù)關(guān)聯(lián)進(jìn)行定位更新。在文獻(xiàn)［19］中對(duì)傳統(tǒng)JPDA 進(jìn)行了改進(jìn)，通過(guò)在預(yù)測(cè)點(diǎn)的周圍生成服從高斯分布的虛擬測(cè)量點(diǎn)以提高定位精度。但是數(shù)據(jù)關(guān)聯(lián)在當(dāng)前時(shí)刻的后驗(yàn)位置估計(jì)要使用所有的歷史時(shí)刻的測(cè)量點(diǎn)，隨著時(shí)刻增加，貝葉斯估計(jì)中的歷史測(cè)量占比也會(huì)隨之增加，而當(dāng)前獲取的測(cè)量所占的比重太小，當(dāng)系統(tǒng)存在模型誤差和噪聲誤差時(shí)，會(huì)減小當(dāng)前時(shí)刻獲取的測(cè)量點(diǎn)對(duì)位置估計(jì)的作用。

在本文中，首先我們?cè)贓KF 的更新過(guò)程中進(jìn)行測(cè)量噪聲誤差協(xié)方差迭代，然后基于計(jì)算噪聲誤差協(xié)方差矩陣的相似度的方法，檢測(cè)信標(biāo)節(jié)點(diǎn)與AUV節(jié)點(diǎn)之間的距離測(cè)量是否受到了惡意攻擊的污染，接著，為了豐富輸入測(cè)量點(diǎn)的樣本以提高定位的精度，我們?cè)跊](méi)有受到惡意噪聲污染的測(cè)量點(diǎn)以及預(yù)測(cè)點(diǎn)的周圍生成服從高斯分布的虛擬測(cè)量點(diǎn)，最后，對(duì)于沒(méi)有受到惡意噪聲污染的測(cè)量點(diǎn)、生成的虛擬粒子進(jìn)行數(shù)據(jù)關(guān)聯(lián)和更新，輸出當(dāng)前時(shí)刻最終的后驗(yàn)位置估計(jì)和后驗(yàn)誤差協(xié)方差，并且會(huì)作為下一時(shí)刻EKF 的先驗(yàn)估計(jì)的輸入，由于消除了惡意噪聲的影響，在提高定位精度的同時(shí)，在下一時(shí)刻惡意噪聲檢測(cè)階段能得到正確的參數(shù)迭代值，從而保證了MJPDA算法正常運(yùn)作。仿真結(jié)果表明，該方法能有效地檢測(cè)出異常值，并且與現(xiàn)有工作相比，所提出的定位方法具有較高的精度。

2 網(wǎng)絡(luò)架構(gòu)與問(wèn)題描述

2.1 網(wǎng)絡(luò)架構(gòu)

為了進(jìn)一步說(shuō)明所提出的水下多AUV 定位算法，我們考慮了一種包含三種不同類型節(jié)點(diǎn)的水聲無(wú)線傳感器網(wǎng)絡(luò)架構(gòu)，如圖1所示。

圖1 水聲無(wú)線傳感器網(wǎng)絡(luò)場(chǎng)景圖Fig.1 Illustration of the UWSN

1）水面浮標(biāo)節(jié)點(diǎn)：水面浮標(biāo)節(jié)點(diǎn)安裝有全球定位系統(tǒng)（Global Positioning System，GPS），通過(guò)與衛(wèi)星之間的電磁波通信獲取自身的準(zhǔn)確實(shí)時(shí)位置。水面浮標(biāo)節(jié)點(diǎn)的作用是為水下信標(biāo)節(jié)點(diǎn)提供自定位服務(wù)。

2）水下信標(biāo)節(jié)點(diǎn)：水下已固定好的各個(gè)信標(biāo)節(jié)點(diǎn)通過(guò)與水面浮標(biāo)節(jié)點(diǎn)的直接通信，并獲得其自身位置。信標(biāo)節(jié)點(diǎn)向待定位的水下多AUV 節(jié)點(diǎn)發(fā)送水聲信號(hào)，AUV在接收到信號(hào)后基于TOA方法獲取到各個(gè)信標(biāo)節(jié)點(diǎn)與AUV 之間準(zhǔn)確的距離測(cè)量。然而，外來(lái)的惡意節(jié)點(diǎn)可以捕獲水聲網(wǎng)絡(luò)中的水下信標(biāo)節(jié)點(diǎn)并偽造其信號(hào)特征，通過(guò)向傳輸?shù)木嚯x測(cè)量數(shù)據(jù)包中注入惡意噪聲，使得AUV接收到來(lái)自被惡意攻擊后的錯(cuò)誤測(cè)距信息，從而影響AUV目標(biāo)節(jié)點(diǎn)的定位精度。

3）待定位的多AUV節(jié)點(diǎn)：在水下待定位的多個(gè)AUV 移動(dòng)目標(biāo)節(jié)點(diǎn)配備了邊緣計(jì)算系統(tǒng)，在得到當(dāng)前時(shí)刻來(lái)自各個(gè)信標(biāo)節(jié)點(diǎn)的測(cè)距信息后，使用本文提出的改進(jìn)聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)算法對(duì)受到污染的距離測(cè)量進(jìn)行篩選，并對(duì)其自身進(jìn)行位置估計(jì)。

2.2 AUV運(yùn)動(dòng)系統(tǒng)模型

在實(shí)際的水下AUV定位中，由于通過(guò)AUV配備的壓力傳感器可以獲得精確的深度，所以AUV的深度獲取和平面的位置估計(jì)是相互獨(dú)立的。此外，由于復(fù)雜的水下環(huán)境，AUV 的垂直速度存在顯著的測(cè)量誤差，這可能會(huì)大大降低深度估計(jì)的運(yùn)動(dòng)模型的準(zhǔn)確性。因此，我們?cè)谟芯_的深度測(cè)量時(shí)，不必建立狀態(tài)變量中含有深度信息的AUV運(yùn)動(dòng)狀態(tài)模型。所以在本文中，待定位的AUV節(jié)點(diǎn)的深度并沒(méi)有被包含在狀態(tài)變量之中，三維定位問(wèn)題被簡(jiǎn)化為二維定位問(wèn)題。狀態(tài)變量隨著時(shí)刻k變化的模型為：

其中，(xAUV，n(k)，yAUV，n(k))，n=1，…，N為N個(gè)待定位AUV 在時(shí)刻k的位置坐標(biāo)，Δt為單位采樣時(shí)間，分別為第n個(gè)AUV 上配備的DVL 提供的其右舷和前向速度，是第n個(gè)AUV 使用羅盤測(cè)量的航向，wk-1，n=[wx，k-1，n，wy，k-1，n]T表示服從高斯分布的過(guò)程噪聲，即wk-1，n～N(0，Qk-1，n)，其中Qk-1，n是過(guò)程噪聲協(xié)方差，T為轉(zhuǎn)置運(yùn)算。

本文使用基于TOA 方法獲取M個(gè)信標(biāo)節(jié)點(diǎn)與N個(gè)AUV 目標(biāo)節(jié)點(diǎn)之間的距離測(cè)量。然而，水聲通信環(huán)境中影響TOA 測(cè)量模型的因素十分復(fù)雜，若不考慮水聲環(huán)境因素（如時(shí)變多普勒、多徑效應(yīng)、分層效應(yīng)等）直接使用TOA 測(cè)距模型，一定會(huì)影響所得到測(cè)距信息的準(zhǔn)確性。最近關(guān)于水聲定位的一些文獻(xiàn)［20-22］中對(duì)于在水聲環(huán)境中使用TOA 的方法獲取準(zhǔn)確測(cè)距信息已經(jīng)得到了充分的研究，為了簡(jiǎn)化，我們?cè)谝酝芯康幕A(chǔ)上假設(shè)了在水下環(huán)境中使用TOA 方法來(lái)獲取準(zhǔn)確的距離信息的合理性，以此更加貼近真實(shí)的水下環(huán)境。

在第k個(gè)時(shí)刻，定義第m個(gè)信標(biāo)節(jié)點(diǎn)與n個(gè)AUV之間的相對(duì)距離測(cè)量為：

其中，AUV 在k時(shí)刻的位置定義為(xAUV，n(k)，yAUV，n(k)，hAUV，n(k))，n=1，…，N，hAUV，n(k)和hBN，m(k)分別為通過(guò)壓力傳感器準(zhǔn)確測(cè)量得到AUV 和信標(biāo)節(jié)點(diǎn)的深度。

由于AUV 和信標(biāo)節(jié)點(diǎn)的深度是由壓力傳感器精確測(cè)量得到的，所以我們利用TOA 的相對(duì)距離以及深度信息進(jìn)行反算可以得到AUV 和信標(biāo)節(jié)點(diǎn)之間的水平距離，因此這里的三維距離測(cè)量lm，n(k)可以被轉(zhuǎn)換為二維距離測(cè)量dm，n(k)

考慮TOA方法的距離測(cè)量噪聲，在第k個(gè)時(shí)刻第m個(gè)信標(biāo)節(jié)點(diǎn)n個(gè)AUV之間的二維測(cè)量方程可以被表示為：

其中，M個(gè)信標(biāo)節(jié)點(diǎn)的位置坐標(biāo)(xBN，m(k)，yBN，m(k))，m=1，…，M假設(shè)通過(guò)浮標(biāo)節(jié)點(diǎn)提供的自定位服務(wù)獲取為已知，并且通過(guò)水聲通信廣播給所有的水下AUV 目標(biāo)節(jié)點(diǎn)；δm，n，k為未遭受惡意噪聲攻擊或者存在惡意噪聲攻擊下的距離測(cè)量噪聲。在沒(méi)有遭受攻擊的情況下，該信標(biāo)節(jié)點(diǎn)的測(cè)量噪聲δm，n，k=ρ被建模為高斯白噪聲ρ～N(0，Rm，n，k)，其中Rm，n，k為測(cè)量噪聲協(xié)方差；在本文中我們考慮了與文獻(xiàn)［15］中類似的惡意噪聲攻擊模型，存在惡意噪聲攻擊的情況下，測(cè)量噪聲δm，n，k=ρ+γc，其中為惡意噪聲，μγ和分別為惡意噪聲的均值和協(xié)方差。

注：三維距離測(cè)量lm，n(k)被轉(zhuǎn)換為二維距離測(cè)量dm，n(k)的過(guò)程中距離測(cè)量噪聲的分布會(huì)偏差為長(zhǎng)尾的彩色測(cè)量噪聲，從而導(dǎo)致測(cè)量噪聲并不滿足高斯分布。然而，在本文中并不需要建立上述測(cè)量噪聲的精確模型，通過(guò)自適應(yīng)地估計(jì)預(yù)測(cè)協(xié)方差矩陣和測(cè)量噪聲協(xié)方差矩陣，緩解了具有不準(zhǔn)確的測(cè)量噪聲模型的負(fù)面影響；并且在本文中使用基于矩陣相似度的惡意噪聲識(shí)別機(jī)制，將測(cè)量噪聲的模型簡(jiǎn)化為高斯白噪聲不會(huì)影響到惡意噪聲的檢測(cè)部分。所以我們?cè)谌S的水下AUV 定位問(wèn)題轉(zhuǎn)換為二維定位問(wèn)題的過(guò)程中，將測(cè)量噪聲模型在本文中簡(jiǎn)化為一個(gè)高斯白噪聲是合理的設(shè)置。

由式（1）、（3），將建立的AUV 運(yùn)動(dòng)系統(tǒng)模型表述為一般形式：

其中，狀態(tài)轉(zhuǎn)移矩陣F=I2，I2為二維單位矩陣，過(guò)程噪聲wn，k-1與測(cè)量噪聲δm，n，k互不相關(guān)，即，其中E[?]=0為期望運(yùn)算。

3 改進(jìn)聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)算法

在上一節(jié)中，我們基于TOA 的方法獲取了假設(shè)位置已知的M個(gè)信標(biāo)節(jié)點(diǎn)與待定位的N個(gè)AUV目標(biāo)節(jié)點(diǎn)之間的距離測(cè)量，對(duì)于其中一個(gè)AUV目標(biāo)節(jié)點(diǎn)，在時(shí)刻k可以從M個(gè)信標(biāo)獲取對(duì)于該AUV的距離信息。在本節(jié)中我們對(duì)該AUV 目標(biāo)節(jié)點(diǎn)進(jìn)行EKF 的先驗(yàn)估計(jì)，得到該AUV 目標(biāo)節(jié)點(diǎn)的先驗(yàn)位置估計(jì)（也就是AUV 的預(yù)測(cè)點(diǎn)）和先驗(yàn)誤差協(xié)方差Pk|k-1，然而，在開(kāi)放的水聲環(huán)境中，外來(lái)的惡意節(jié)點(diǎn)很容易在所得到的M個(gè)距離測(cè)量中注入惡意噪聲，從而導(dǎo)致最終定位精度的下降。所以我們?cè)贏UV的初步定位之前通過(guò)惡意噪聲檢測(cè)以排除其影響。

首先，利用期望最大化算法對(duì)M個(gè)距離測(cè)量對(duì)應(yīng)的先驗(yàn)協(xié)方差Pm，k|k-1和測(cè)量噪聲誤差協(xié)方差Rm，k進(jìn)行在線迭代估計(jì)；接著，我們基于計(jì)算測(cè)量噪聲誤差協(xié)方差矩陣相似度來(lái)識(shí)別距離測(cè)量是否受到惡意噪聲污染，并且丟棄受到惡意噪聲污染的距離測(cè)量。對(duì)于保留的P個(gè)未被惡意噪聲污染的距離測(cè)量每三個(gè)一組進(jìn)行分組，利用最小二乘法對(duì)AUV目標(biāo)節(jié)點(diǎn)完成對(duì)AUV的初步定位，由此我們?cè)诋?dāng)前時(shí)刻可以獲得個(gè)AUV 目標(biāo)節(jié)點(diǎn)位置的最小二乘解，n=1，…，Ni（由初步定位得到的Ni個(gè)測(cè)量點(diǎn)）；為了豐富測(cè)量點(diǎn)的樣本以提高定位的精度，與傳統(tǒng)的JPDA 算法相比，我們?cè)跊](méi)有受到惡意噪聲污染的測(cè)量點(diǎn)以及預(yù)測(cè)點(diǎn)的周圍生成服從高斯分布的虛擬測(cè)量點(diǎn)；最后，對(duì)于沒(méi)有受到惡意噪聲污染的測(cè)量點(diǎn)以及生成的虛擬粒子進(jìn)行數(shù)據(jù)關(guān)聯(lián)和更新，完成當(dāng)前時(shí)刻AUV的精確定位。輸出當(dāng)前時(shí)刻最終的后驗(yàn)位置估計(jì)和后驗(yàn)誤差協(xié)方差，并且會(huì)作為下一時(shí)刻EKF 的先驗(yàn)估計(jì)的輸入。由于輸出的后驗(yàn)位置估計(jì)和后驗(yàn)誤差協(xié)方差是使用通過(guò)惡意噪聲檢測(cè)階段的測(cè)量計(jì)算得到的，所以能有效保證MJPDA算法在惡意噪聲環(huán)境下的魯棒性。如圖2為MJPDA的流程圖。

圖2 改進(jìn)的聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)算法流程圖Fig.2 Flowchart of MJPDA

3.1 EKF先驗(yàn)估計(jì)

由于在式（5）中的測(cè)量方程中的h(xn，k，xBN，m，n，yBN，m，n)是非線性的，所以我們使用EKF 的先驗(yàn)估計(jì)對(duì)進(jìn)行AUV 當(dāng)前位置的預(yù)測(cè)。利用k-1時(shí)刻通過(guò)MJPDA 得到的AUV 目標(biāo)節(jié)點(diǎn)后驗(yàn)位置估計(jì)和后驗(yàn)誤差協(xié)方差Pk-1|k-1，得到k時(shí)刻AUV 的先驗(yàn)狀態(tài)估計(jì)和先驗(yàn)誤差協(xié)方差Pk|k-1。

同時(shí)我們得到k時(shí)刻的AUV 位置的預(yù)測(cè)點(diǎn)

3.2 惡意噪聲污染檢測(cè)

在本節(jié)中我們?cè)O(shè)計(jì)了一種基于噪聲協(xié)方差矩陣相似度的惡意噪聲檢測(cè)機(jī)制，首先，針對(duì)水聲環(huán)境中由于外部干擾產(chǎn)生的時(shí)變?cè)肼晠f(xié)方差，將在線EM算法應(yīng)用于EKF中，基于歷史k個(gè)時(shí)刻信標(biāo)節(jié)點(diǎn)與AUV 之間的測(cè)距信息dm，1：k，m=1，…M，來(lái)在線估計(jì)先驗(yàn)誤差協(xié)方差Pm，k|k-1和測(cè)量噪聲協(xié)方差Rm，k，其中表示信標(biāo)節(jié)點(diǎn)m由時(shí)刻1 到時(shí)刻k與AUV 之間的距離測(cè)量；在得到迭代N次的測(cè)量噪聲協(xié)方差后，我們分別計(jì)算每個(gè)信標(biāo)對(duì)應(yīng)測(cè)量噪聲協(xié)方差矩陣與所設(shè)置的初代測(cè)量噪聲協(xié)方差矩陣Rk的相似性，用來(lái)判斷k時(shí)刻信標(biāo)節(jié)點(diǎn)m使用TOA 方法獲取的距離測(cè)量是否受到了惡意噪聲的污染。

使用EM 算法迭代參數(shù)Φm，k?{Pm，k|k-1，Rm，k}，主要分為兩步，即E步（Expectation）和M步（Maximization）迭代求解距離測(cè)量關(guān)于模型參數(shù)Φm，k的對(duì)數(shù)似然函數(shù)logp(xk，dm，1：k，Φm，k)的極大似然，從而得到參數(shù)Φm，k的估計(jì)值。

E 步：根據(jù)貝葉斯公式，距離測(cè)量dm，1：k和AUV在時(shí)刻k迭代估計(jì)的參數(shù)Φk完全數(shù)據(jù)聯(lián)合概率分布可以被定義為：

所以，用公式（9）～（11）可以計(jì)算出完全數(shù)據(jù)對(duì)數(shù)似然函數(shù)logpΦm，k(xk，dm，1：k)。

在l+1代，通過(guò)利用l代的后驗(yàn)狀態(tài)估計(jì)和后驗(yàn)協(xié)方差得到更精確的非線性測(cè)距測(cè)量方程的泰勒線性化近似h(xk，xBN，m，yBN，m)。

計(jì)算l+1迭代的卡爾曼增益：

計(jì)算l+1代的后驗(yàn)狀態(tài)估計(jì)：

計(jì)算l+1代的后驗(yàn)協(xié)方差：

迭代l次的參數(shù)下的后驗(yàn)密度函數(shù)可以被近似高斯分布，其均值和方差分別為l+1迭代中的后驗(yàn)估計(jì)和后驗(yàn)協(xié)方差。

M步：

在M 步中，極大化在E 步中得到的Q(Φm，k，)，然后求解等式獲取l+1 代的先驗(yàn)誤差協(xié)方差、測(cè)量噪聲協(xié)方差；并且在通過(guò)N次迭代后得到了每一個(gè)距離測(cè)量對(duì)應(yīng)的測(cè)量噪聲協(xié)方差矩陣。

當(dāng)UWSN 中部分信標(biāo)節(jié)點(diǎn)受到惡意噪聲的影響時(shí)，相應(yīng)迭代出的測(cè)量噪聲誤差協(xié)方差將產(chǎn)生異常值。借助圖像識(shí)別中的相似度原理，我們分別計(jì)算每個(gè)信標(biāo)對(duì)應(yīng)測(cè)量噪聲協(xié)方差矩陣與所設(shè)置的初代測(cè)量噪聲協(xié)方差矩陣Rk的相似性，定義和Rk的矩陣相似度為：

其中，相似度的計(jì)算結(jié)果αm∈[-1，1]，αm的值越接近1，說(shuō)明矩陣的相似度越高，越接近0 則說(shuō)明兩個(gè)矩陣的相似度越低。經(jīng)過(guò)惡意噪聲的污染迭代出的噪聲協(xié)方差矩陣與設(shè)置的初代噪聲協(xié)方差矩陣具有明顯差異，所以他們經(jīng)計(jì)算得到的相似性最低。

通過(guò)計(jì)算噪聲誤差協(xié)方差的相似度，我們將受到惡意噪聲污染的距離測(cè)量丟棄，保留未遭受惡意噪聲污染的距離測(cè)量。在識(shí)別檢測(cè)之前，在時(shí)刻k對(duì)于一個(gè)AUV 目標(biāo)節(jié)點(diǎn)我們有M個(gè)測(cè)距信息m=1，…M，在噪聲識(shí)別檢測(cè)以后保留有P個(gè)未遭受惡意噪聲污染的距離測(cè)量(P≤M)。這樣，在進(jìn)行初步定位估計(jì)之前，我們通過(guò)惡意噪聲污染檢測(cè)消除了其對(duì)AUV定位效果的影響。

3.3 分組和初步定位

上一節(jié)中得到了P個(gè)檢測(cè)為沒(méi)有受到惡意噪聲污染的距離測(cè)量，為了獲取到盡可能多的測(cè)量點(diǎn)以提高定位精度，我們直接將P個(gè)距離測(cè)量按照排列組合進(jìn)行分組，然后利用最小二乘法來(lái)獲取AUV初步定位，n=1，…，Ni，其中為P個(gè)測(cè)距信息按照三個(gè)一組被分為不同組的組數(shù)。

對(duì)于其中的一組，根據(jù)距離測(cè)量值可建立以下方程：

其中，dn，1(k)、dn，2(k)、dn，3(k)為第n組內(nèi)信標(biāo)節(jié)點(diǎn)的距離測(cè)量，(xBN，1：3，yBN，1：3)為分組內(nèi)對(duì)應(yīng)信標(biāo)節(jié)點(diǎn)的位置坐標(biāo)，為組內(nèi)信標(biāo)節(jié)點(diǎn)對(duì)應(yīng)經(jīng)迭代得到的噪聲誤差協(xié)方差矩陣。為了使用最小二乘法得到AUV 的初步定位，式（19）被表達(dá)為線性矩陣形式Y(jié)=AX：

為了求解AUV的初步定位，定義如下殘差向量方程：

使AUV的位置坐標(biāo)殘差的平方和最小，對(duì)殘差向量方程求導(dǎo)并令其為零：

得到AUV 位置的最小二乘解，完成AUV 的初步定位：

3.4 產(chǎn)生虛擬測(cè)量

考慮到信標(biāo)節(jié)點(diǎn)的部署成本，在同一片海域當(dāng)中不適宜部署過(guò)多的信標(biāo)節(jié)點(diǎn)。同時(shí)，在有惡意噪聲的環(huán)境下，MJPDA 算法丟棄了受惡意噪聲污染的測(cè)量點(diǎn)。針對(duì)測(cè)量樣本的稀疏性，與傳統(tǒng)的JPDA算法相比，MJPDA 為了提高定位精度，在通過(guò)AUV初步定位得到的良好測(cè)量點(diǎn)與預(yù)測(cè)點(diǎn)的附近產(chǎn)生了服從高斯分布的虛擬測(cè)量。

為了產(chǎn)生虛擬測(cè)量的范圍，我們首先確定未遭受惡意噪聲影響的信標(biāo)節(jié)點(diǎn)n當(dāng)前時(shí)刻對(duì)應(yīng)的新息協(xié)方差：

并且由新息協(xié)方差確定測(cè)量點(diǎn)對(duì)應(yīng)的關(guān)聯(lián)門的面積，在本文中關(guān)聯(lián)門是指以良好測(cè)量點(diǎn)與預(yù)測(cè)點(diǎn)為中心，用來(lái)限制產(chǎn)生虛擬測(cè)量的范圍。這個(gè)范圍的大小主要由當(dāng)前時(shí)刻該測(cè)量對(duì)應(yīng)的新息協(xié)方差和關(guān)聯(lián)門概率確定。落入到關(guān)聯(lián)門內(nèi)的測(cè)量被稱為候選測(cè)量點(diǎn)。通常，關(guān)聯(lián)門被認(rèn)為是一個(gè)“gsigma”橢圓面［23］，關(guān)聯(lián)門的面積為：

其中，α為MN 坐標(biāo)向量的維數(shù)，cα為α維關(guān)聯(lián)門的單位面積。本文中為二維定位場(chǎng)景，所以α=2，c2=π。同時(shí)，γ服從具有兩個(gè)自由度的卡方分布。然后，如果檢測(cè)到正確的輸入測(cè)量點(diǎn)的概率PG可以表示為：。

為了簡(jiǎn)化虛擬測(cè)量的生成過(guò)程，我們將橢圓關(guān)聯(lián)門簡(jiǎn)化為正方形，邊長(zhǎng)定義為l：

其中，(Nx(k)，Ny(k))為產(chǎn)生的所有虛擬測(cè)量的位置坐標(biāo)，N1(l，l/2)和N2(l，l/2)是均值為l、方差為l/2 的高斯分布。因此，最終輸入到數(shù)據(jù)關(guān)聯(lián)和更新階段中的測(cè)量點(diǎn)主要包含：經(jīng)過(guò)AUV初步定位得到的測(cè)量點(diǎn)、AUV 位置的預(yù)測(cè)點(diǎn)、以及在他們周圍產(chǎn)生的虛擬測(cè)量點(diǎn)，為了接下來(lái)推導(dǎo)方便，將它們統(tǒng)一定義為zm(k)。

3.5 數(shù)據(jù)關(guān)聯(lián)和更新

我們?cè)谏弦还?jié)當(dāng)中得到了要輸入到數(shù)據(jù)關(guān)聯(lián)算法的測(cè)量，則這測(cè)量關(guān)于先驗(yàn)估計(jì)的新息為：

在初步定位得到的量測(cè)點(diǎn)中，如果存在受到所選擇的信標(biāo)和AUV 的相對(duì)位置的影響導(dǎo)致定位精度極差的測(cè)量點(diǎn)，首先，我們會(huì)對(duì)其進(jìn)行篩選，將檢驗(yàn)統(tǒng)計(jì)量Tm(k)計(jì)算為：

為了驗(yàn)證輸入測(cè)量，定義了以下假設(shè)：

其中，H0，m如果為真，則輸入的測(cè)量值將被添加到有效測(cè)量值集合中，否則將被丟棄。注意，有效測(cè)量值為，i=1，2，…，Nv，其中Nv是有效測(cè)量集合中的有效測(cè)量點(diǎn)個(gè)數(shù)。所以在經(jīng)過(guò)以上的篩選之后，AUV 和信標(biāo)節(jié)點(diǎn)的相對(duì)位置較差導(dǎo)致所構(gòu)建的定位方程秩虧對(duì)最終的精確定位結(jié)果影響不大。

然后，在我們的推導(dǎo)當(dāng)中關(guān)聯(lián)矩陣可以被表示為：

其中，關(guān)聯(lián)矩陣的第一列都為1 表示測(cè)量值并非來(lái)自于AUV 目標(biāo)節(jié)點(diǎn)。因?yàn)樗械妮斎霚y(cè)量值都可以通過(guò)雜波或虛警來(lái)產(chǎn)生，關(guān)聯(lián)矩陣的第二列表示輸入的測(cè)量值是否與AUV相關(guān)聯(lián)，如果相關(guān)聯(lián)則對(duì)應(yīng)的元素為1，如果不關(guān)聯(lián)則對(duì)應(yīng)的元素為0。對(duì)關(guān)聯(lián)矩陣進(jìn)行拆分，我們得到互聯(lián)矩陣Ωi(k)，每一個(gè)互聯(lián)矩陣表示不同的互聯(lián)事件。

計(jì)算每一個(gè)互聯(lián)事件的概率，也就是測(cè)量點(diǎn)為目標(biāo)準(zhǔn)確定位的概率βm，與文獻(xiàn)［24］的隨機(jī)權(quán)重不同，本文所提算法的權(quán)重是確定性的。與隨機(jī)權(quán)重相比，確定性權(quán)重能夠使基于良好測(cè)量點(diǎn)的估計(jì)狀態(tài)的貢獻(xiàn)最大化，使惡意測(cè)量的估計(jì)狀態(tài)的貢獻(xiàn)最小化。

測(cè)量值并非來(lái)自于AUV目標(biāo)節(jié)點(diǎn)的概率β0：

最后，根據(jù)計(jì)算得到的互聯(lián)概率的值，我們得到后驗(yàn)狀態(tài)估計(jì)和后驗(yàn)誤差協(xié)方差。

卡爾曼增益矩陣：

后驗(yàn)狀態(tài)估計(jì)為：

后驗(yàn)誤差協(xié)方差計(jì)算為：

在得到當(dāng)前時(shí)刻最終的后驗(yàn)狀態(tài)估計(jì)和后驗(yàn)協(xié)方差，我們將它們放入到MJPDA的下一時(shí)刻的先驗(yàn)計(jì)算中使用，可以有效提高算法在惡意噪聲影響下的魯棒性。

4 仿真分析

在本節(jié)中，對(duì)于本文所提出的MJPDA 算法在MATLAB R2021a上進(jìn)行了仿真模擬，對(duì)二維場(chǎng)景下兩個(gè)同時(shí)做勻速運(yùn)動(dòng)的目標(biāo)場(chǎng)景進(jìn)行了測(cè)試，以驗(yàn)證所提定位算法的有效性。首先通過(guò)仿真實(shí)例證明在相同參數(shù)下所提定位算法優(yōu)于所對(duì)比的算法；然后通過(guò)進(jìn)行不同參數(shù)的變化對(duì)均方根誤差（Root Mean Square Error，RMSE）的影響，證明所提算法對(duì)于處理惡意噪聲提升定位性能是有意義的。

在圖3 所示的二維定位場(chǎng)景中，600 m×600 m的區(qū)域內(nèi)隨機(jī)部署了6 個(gè)信標(biāo)節(jié)點(diǎn)來(lái)定位2個(gè)AUV目標(biāo)節(jié)點(diǎn)。藍(lán)色實(shí)線和紅色實(shí)線分別代表兩AUV利用改進(jìn)的聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)算法在惡意噪聲影響下所估計(jì)得到的運(yùn)動(dòng)軌跡，藍(lán)色虛線和紅色虛線則分別表示兩AUV 沿著固定的軌跡進(jìn)行勻速直線運(yùn)動(dòng)，在軌跡附近的藍(lán)色點(diǎn)與紅色點(diǎn)表示為提高定位精度而產(chǎn)生的測(cè)量點(diǎn)。MJPDA 的門概率被設(shè)置為PG=0.99，檢測(cè)概率被設(shè)置為PD=0.95。不考慮信標(biāo)節(jié)點(diǎn)的功耗，假設(shè)所有的信標(biāo)節(jié)點(diǎn)總可以與AUV 節(jié)點(diǎn)進(jìn)行有效的通信。每個(gè)仿真實(shí)驗(yàn)結(jié)果通過(guò)500 次的蒙特卡洛模擬運(yùn)行得到，所對(duì)比的算法JPDA 和AEKF［9］的默認(rèn)參數(shù)設(shè)置與MJPDA 一致，在表1中給出了仿真中所使用的參數(shù)。

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

圖3 定位兩個(gè)AUV的位置估計(jì)Fig.3 Simulation results for two AUVs tracking

定義在每一次蒙特卡洛模擬中的平均定位誤差（Average Localization Error，ALE）：

根據(jù)ALE 我們可以得到三種算法平均誤差的累積分布函數(shù)（Cumulative Distribution Function，CDF），CDF是指在某個(gè)精度門限下的定位次數(shù)在總定位次數(shù)中所占比例，圖4 的仿真結(jié)果可以看出，MJPDA 的百分之八十的定位誤差都小于3.3 m，相比之下，JPDA 算法和AEKF 百分之八十的定位誤差分別小于5.1 m 和7.9 m。圖4 的仿真結(jié)果表明，總體來(lái)說(shuō)，定位系統(tǒng)處于存在較大惡意噪聲和累積估計(jì)誤差的影響情況下，MJPDA 的定位精度優(yōu)于現(xiàn)有算法JPDA和AEKF。

圖4 經(jīng)驗(yàn)分布函數(shù)與定位誤差Fig.4 CDF and Localization Error

接下來(lái)，在圖5 至圖7 的仿真結(jié)果中我們將顯示在不同參數(shù)的變化下（惡意噪聲出現(xiàn)的概率、惡意噪聲的均值、信標(biāo)節(jié)點(diǎn)數(shù)量）詳細(xì)討論MJPDA、JPDA 和AEKF 的三種算法的比較，利用RMSE 進(jìn)行算法性能的評(píng)價(jià)。

圖5 惡意噪聲出現(xiàn)概率對(duì)均方誤差的影響Fig.5 RMSE and different probability of Malicious Noise

其中，K表示一次蒙特卡洛模擬中的總采樣時(shí)長(zhǎng)，Tn表示蒙特卡洛運(yùn)行的次數(shù)，且K=130，Tn=500。

在通常的水聲對(duì)抗情形下，噪聲干擾往往是飽和的，所以我們利用惡意噪聲出現(xiàn)概率來(lái)體現(xiàn)對(duì)抗環(huán)境中惡意噪聲干擾的飽和性。在圖5 中，橫軸為其他參數(shù)不變的情況下惡意噪聲出現(xiàn)的概率從0.1到1，縱軸為蒙特卡洛的定位誤差。隨著惡意噪聲對(duì)一個(gè)信標(biāo)節(jié)點(diǎn)與AUV 之間的距離測(cè)量產(chǎn)生影響概率的增加，對(duì)比的兩種算法的定位誤差都在增加，其中AEKF的定位誤差增加的較為顯著，這是因?yàn)檩^大的惡意噪聲破壞了測(cè)量噪聲誤差協(xié)方差的迭代環(huán)境，從而影響了定位效果。當(dāng)惡意噪聲出現(xiàn)的概率比較小的時(shí)候，MJPDA 和JPDA 算法的定位精度比較接近。隨著惡意噪聲出現(xiàn)的概率增加時(shí)，由于所提出的MJPDA 提前篩選出了受污染的距離測(cè)量，所以仍然可以保持較好的定位精度，仿真結(jié)果表明，MJPDA 在場(chǎng)景中經(jīng)常出現(xiàn)惡意噪聲時(shí)可以具有更好的定位精度。

在圖6中，保持了其他參數(shù)的不變的同時(shí)，惡意噪聲的均值在每一次蒙特卡洛中從3 m 增加到10 m，隨著惡意噪聲誤差均值的增加，AEKF的RMSE呈線性增加，JPDA 的RMSE 增加相對(duì)緩慢，MJPDA的定位誤差幾乎沒(méi)有增長(zhǎng)。當(dāng)惡意噪聲誤差的標(biāo)準(zhǔn)差從3 增加到10 時(shí)，仿真結(jié)果表明，JPDA 通過(guò)關(guān)聯(lián)門丟棄具有較大殘差的離群值來(lái)實(shí)現(xiàn)其魯棒性，但如果有多個(gè)離群值連續(xù)出現(xiàn)，其定位效果將會(huì)降低。MJPDA 受惡意噪聲誤差擾動(dòng)的影響較小，具有較強(qiáng)的魯棒性，這是因?yàn)镸JPDA在惡意噪聲的檢測(cè)環(huán)節(jié)中，丟棄了所有受到惡意噪聲污染的距離測(cè)量。

圖6 惡意噪聲的均值對(duì)均方誤差的影響Fig.6 RMSE and the mean value of Malicious Noise

為顯示信標(biāo)節(jié)點(diǎn)數(shù)量對(duì)于錨節(jié)點(diǎn)定位精度的影響，我們使信標(biāo)節(jié)點(diǎn)的數(shù)量從4 個(gè)增加到10 個(gè)，如圖7 所示。雖然增加了信標(biāo)節(jié)點(diǎn)的數(shù)量，但惡意噪聲對(duì)AEKF 的均方根誤差的影響仍然較大。同時(shí)，隨著信標(biāo)節(jié)點(diǎn)數(shù)量的增加，JPDA 和MJPDA 算法的定位誤差都在減小。在信標(biāo)節(jié)點(diǎn)數(shù)量較少時(shí)，由于MJPDA 在預(yù)測(cè)點(diǎn)和未受到惡意噪聲污染的測(cè)量點(diǎn)附近都生成了虛擬測(cè)量，比JPDA 算法擁有更豐富的樣本值，因此JPDA算法的效果不如MJPDA好。仿真結(jié)果表明，隨著信標(biāo)節(jié)點(diǎn)數(shù)量的增加，兩種算法的定位精度趨于一致，但在信標(biāo)節(jié)點(diǎn)數(shù)量較少時(shí)，MJPDA 擁有更好的定位精度，即MJPDA 在水下物聯(lián)網(wǎng)的節(jié)點(diǎn)部署稀疏時(shí)擁有更好的定位表現(xiàn)。

圖7 信標(biāo)節(jié)點(diǎn)數(shù)量對(duì)RMSE的影響Fig.7 RMSE and the number of beacon nodes

5 結(jié)論

本文對(duì)存在惡意噪聲攻擊情況下的AUV 定位進(jìn)行了研究，提出了一種改進(jìn)的聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)定位方法（MJPDA），考慮了水下環(huán)境中存在惡意噪聲攻擊情況下的AUV定位追蹤，設(shè)計(jì)了一種基于誤差協(xié)方差迭代算法的惡意噪聲攻擊檢測(cè)方法，通過(guò)基于計(jì)算當(dāng)前時(shí)刻迭代出的誤差協(xié)方差相似度的方式檢測(cè)出受到惡意噪聲污染的測(cè)量，并且完成AUV 的初步定位。與傳統(tǒng)的JPDA 相比，在預(yù)測(cè)點(diǎn)和未受到惡意噪聲污染的測(cè)量點(diǎn)附近都生成了虛擬測(cè)量點(diǎn)，對(duì)這些測(cè)量點(diǎn)進(jìn)行數(shù)據(jù)關(guān)聯(lián)處理與更新，完成對(duì)AUV目標(biāo)節(jié)點(diǎn)的精確定位。仿真結(jié)果證明，所提出的MJPDA定位算法能夠保證存在惡意噪聲攻擊情況的網(wǎng)絡(luò)環(huán)境下具有良好的性能，對(duì)研究如何安全定位水下機(jī)器人具有一定的意義。