羅浚 柏祥華 安良

（1.東南大學(xué)水聲信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室，南京，210096）

（2.海裝上海局駐南京地區(qū)第一軍事代表室，南京，210000）

水聲被動(dòng)測(cè)距是水聲信號(hào)處理和水聲對(duì)抗領(lǐng)域的難點(diǎn)之一[1]，面臨著信號(hào)形式、頻率、DOA、時(shí)間等先驗(yàn)信息缺失問題。尤其是水聲脈沖信號(hào)的被動(dòng)定位，由于水聲脈沖信號(hào)是間歇式發(fā)射的非合作水聲信號(hào)，無法通過長(zhǎng)時(shí)間的積分獲得時(shí)間增益。信號(hào)占空比較大，導(dǎo)致所獲取的數(shù)據(jù)量相對(duì)于輻射噪聲等信號(hào)較少，這也加大了脈沖信號(hào)被動(dòng)測(cè)距的難度。目前，未知參量脈沖信號(hào)被動(dòng)測(cè)距有三類方法：三元陣測(cè)距、匹配場(chǎng)處理和目標(biāo)運(yùn)動(dòng)分析（Target Motion Analysis，TMA）[2]。對(duì)于三元陣測(cè)距，水聲信道的復(fù)雜性導(dǎo)致傳播時(shí)延起伏，實(shí)現(xiàn)多水聽器布陣也較為困難。匹配場(chǎng)處理的測(cè)距性能對(duì)海洋聲學(xué)參數(shù)的預(yù)估精度要求很高，如果參數(shù)失配則無法實(shí)現(xiàn)定位。因此，在已知系統(tǒng)的狀態(tài)模型與觀測(cè)模型的前提下，采用TMA是較為有效的被動(dòng)定位方法[3]。

經(jīng)典的TMA被動(dòng)定位方法分為純方位TMA、方位-信號(hào)到達(dá)時(shí)間差TMA等。前者需要觀測(cè)站機(jī)動(dòng)，因此具有很大的局限性[4]；后者充分利用DOA與 TOA作為觀測(cè)量建立模型，使得該方法不需要觀測(cè)站機(jī)動(dòng)。在非線性系統(tǒng)中，經(jīng)常使用EKF方法進(jìn)行跟蹤，通過保留系統(tǒng)泰勒展開式的一階項(xiàng)、省略高階項(xiàng)的方法得到近似線性化的系統(tǒng)，該方法必定會(huì)產(chǎn)生誤差，且誤差會(huì)隨著樣本數(shù)量累積，可能會(huì)導(dǎo)致濾波器發(fā)散[5]?；赨nscented變換的UKF方法沒有忽略系統(tǒng)高階項(xiàng)，而是通過一系列的采樣點(diǎn)獲取隨機(jī)變量的均值和方差[6]。當(dāng)狀態(tài)變量通過實(shí)際的非線性系統(tǒng)之后，后驗(yàn)均值與協(xié)方差估計(jì)值可精確至三階，進(jìn)而改善系統(tǒng)的計(jì)算精度與非線性效果[7-8]。

海洋水聲環(huán)境十分復(fù)雜，多途傳播、陣形畸變都是不可避免的，這些對(duì)聲信號(hào)傳播的影響很大[9]。這使得用于聲源定位的方位、信號(hào)到達(dá)時(shí)間差等參數(shù)估計(jì)誤差較電磁信號(hào)源大，導(dǎo)致用于電磁信號(hào)源定位的傳統(tǒng)TMA方法進(jìn)行水聲信號(hào)跟蹤時(shí)，結(jié)果容易發(fā)散、系統(tǒng)性能下降。

本文提出基于DOA與TOA的MAF-PLE-UKF聯(lián)合方法，首先使用MAF方法對(duì)用于初值估計(jì)的DOA和TOA觀測(cè)數(shù)據(jù)進(jìn)行平滑預(yù)處理，之后使用PLE方法估計(jì)目標(biāo)狀態(tài)初值，再使用UKF對(duì)目標(biāo)信息進(jìn)行跟蹤，以提高目標(biāo)跟蹤的收斂性能和定位精度。

1 定位模型建立

1.1 目標(biāo)狀態(tài)模型

假設(shè)觀測(cè)站固定不動(dòng)，以觀測(cè)站為原點(diǎn)，正東方向?yàn)閤軸、正北方向?yàn)閥軸建立直角坐標(biāo)系。目標(biāo)按照某一未知方向做勻速直線運(yùn)動(dòng)，并發(fā)射周期恒定的脈沖信號(hào)，如圖1所示。

圖1 目標(biāo)與觀測(cè)站的運(yùn)動(dòng)幾何關(guān)系的模型

記初始時(shí)刻為0時(shí)刻，令tk表示第k時(shí)刻，則目標(biāo)在第k時(shí)刻的狀態(tài)矢量為

式中，xk、yk為目標(biāo)第k時(shí)刻的位置，Tr為目標(biāo)發(fā)射脈沖的周期，dx、dy為目標(biāo)在Tr時(shí)間內(nèi)沿x、y軸方向的位移。系統(tǒng)下一時(shí)刻的狀態(tài)方程為

式中，Wk表示系統(tǒng)噪聲，假設(shè)其服從高斯分布，其協(xié)方差矩陣記為Qk。Φ(k+1/k)表示系統(tǒng)狀態(tài)矩陣：

1.2 目標(biāo)觀測(cè)模型

觀測(cè)站進(jìn)行第k次觀測(cè)時(shí)，目標(biāo)的角度βk與相鄰脈沖的到達(dá)時(shí)間差ΔTk分別為

根據(jù)式（3）、（4），可將觀測(cè)模型表達(dá)式記為

式中，nk為觀測(cè)誤差矩陣，其協(xié)方差矩陣記為Rk。

2 定位方法原理

2.1 滑動(dòng)平均濾波法進(jìn)行數(shù)據(jù)預(yù)處理

針對(duì)水聲脈沖參數(shù)估計(jì)誤差較大的問題，首先使用滑動(dòng)平均濾波法對(duì)用于狀態(tài)初值估計(jì)的脈沖DOA與TOA觀測(cè)數(shù)據(jù)進(jìn)行平滑處理。該方法的基本思想是：一個(gè)時(shí)間點(diǎn)的取值由當(dāng)前時(shí)間點(diǎn)及其鄰域內(nèi)時(shí)間點(diǎn)的平均值來代替[10]。通常是使用一個(gè)濾波窗實(shí)現(xiàn)的，設(shè)窗長(zhǎng)為M。在時(shí)間點(diǎn)tk時(shí)，窗口選中的時(shí)間區(qū)域?yàn)镾t，輸入數(shù)據(jù)為g(n)，則該時(shí)間點(diǎn)濾波器的輸出f(tk)為

選取窗長(zhǎng)為5個(gè)時(shí)間點(diǎn)，使用滑動(dòng)平均濾波法，對(duì)DOA和TOA觀測(cè)數(shù)據(jù)進(jìn)行預(yù)處理。

2.2 偽線性方法估計(jì)初值

實(shí)際應(yīng)用中，UKF方法對(duì)初始值的準(zhǔn)確性較為敏感。在沒有水聲脈沖參數(shù)先驗(yàn)信息的情況下，本模型采用偽線性方法對(duì)初值進(jìn)行估計(jì)。對(duì)于第k次觀測(cè)，根據(jù)坐標(biāo)的幾何關(guān)系，可將（3）、（4）式轉(zhuǎn)化為線性表示形式：

假設(shè)含有噪聲的觀測(cè)量為βmk、ΔTmk，它們與真實(shí)值βk、ΔTk的關(guān)系為

其中，nβk與nΔTk分別表示方位、脈沖到達(dá)時(shí)間差的估計(jì)誤差，并服從均值為零的高斯分布。將觀測(cè)關(guān)系式（8）、（9）代入線性化的關(guān)系式（6）、（7）中。在噪聲較小，并且觀測(cè)目標(biāo)距離觀測(cè)站較遠(yuǎn)時(shí)，可做式（10）～（13）的近似：

對(duì)公式進(jìn)行計(jì)算與整理可得到，系統(tǒng)線性觀測(cè)方程的等效角度與脈沖到達(dá)時(shí)間差的測(cè)量噪聲為rknβk與nΔTk（rk表示k時(shí)刻目標(biāo)與觀測(cè)站間的距離）。記系統(tǒng)觀測(cè)方程的等效噪聲Nk為

由此便可將觀測(cè)模型表達(dá)式（5）表示為線性形式：

對(duì)于每一時(shí)刻tk，根據(jù)式（14），觀測(cè)方程都可表示為矩陣形式：

式中，對(duì)角度βk與脈沖到達(dá)時(shí)間差ΔTk的測(cè)量是相互獨(dú)立的，方差為。由此可得Nk的協(xié)方差矩陣Rk為

經(jīng)過N次觀測(cè)后，線性觀測(cè)方程表示為

對(duì)初始狀態(tài)量X0的加權(quán)最小二乘估計(jì)[11]：

式中，該算法需要預(yù)先設(shè)定矩陣R的值進(jìn)行計(jì)算。0為偽線性估計(jì)的系統(tǒng)狀態(tài)初值。

2.3 無跡卡爾曼濾波方法定位與跟蹤

獲得估計(jì)的初值后，使用無跡卡爾曼濾波方法對(duì)目標(biāo)狀態(tài)進(jìn)行跟蹤。UKF以Unscented變換為基礎(chǔ)[12]。Unscented變換根據(jù)系統(tǒng)狀態(tài)向量x的維數(shù)n，選取2n+1個(gè)近似高斯分布的Sigma點(diǎn)，均值為系統(tǒng)狀態(tài)量均值，方差為系統(tǒng)協(xié)方差P[13]。通過這組Sigma點(diǎn)對(duì)系統(tǒng)狀態(tài)向量的概率密度函數(shù)進(jìn)行近似化，而不是對(duì)系統(tǒng)的非線性函數(shù)進(jìn)行近似線性化，從而克服了線性化所帶來的濾波誤差[14]。

在第k+1次觀測(cè)時(shí)，基于DOA與TOA的UKF算法的具體計(jì)算步驟如下：

（1）首先根據(jù)k時(shí)刻的狀態(tài)估計(jì)量(k/k)、系統(tǒng)估計(jì)協(xié)方差P(k/k)以及Unscented變換的原理，求得最初的2n+1個(gè)Sigma點(diǎn)。將Sigma點(diǎn)組成的向量記為Ei(k/k)，相應(yīng)的權(quán)值記為和，上標(biāo)m與均值相對(duì)應(yīng)，上標(biāo)c與方差相對(duì)應(yīng)。其中i為從1～2n+1的值。Sigma向量和采樣權(quán)重的取值分別為

（2）將系統(tǒng)的狀態(tài)方程（1）式記為f，求得Sigma點(diǎn)的預(yù)測(cè)值Ei(k+1/k)，預(yù)測(cè)值的估計(jì)(k+1/k)和預(yù)測(cè)值的協(xié)方差P(k+1/k)為

（3）根據(jù)系統(tǒng)的觀測(cè)方程（5）式，計(jì)算觀測(cè)值Zi(k+1/k)，并計(jì)算觀測(cè)值的估計(jì)(k+1/k)和協(xié)方差Pzz(k+1/k)：

（4）計(jì)算Ei(k+1/k)、Zi(k+1/k)的協(xié)方差Pxz(k+1/k)與系統(tǒng)的卡爾曼增益G(k+1)：

（5）進(jìn)行協(xié)方差P(k+1/k+1)和系統(tǒng)狀態(tài)(k+1/k+1)的更新：

由此便完成對(duì)k+1時(shí)刻狀態(tài)方程的估計(jì)。

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

假設(shè)觀測(cè)站靜止于坐標(biāo)原點(diǎn)，目標(biāo)真實(shí)的初始狀態(tài)位置x0為22 km，y0為16 km。目標(biāo)發(fā)射周期Tr為30 s的脈沖信號(hào)；其在一個(gè)脈沖周期內(nèi)的位移dx為-150 m，dy為60 m，其中系統(tǒng)位置噪聲標(biāo)準(zhǔn)差σrx與σry為1 m，一個(gè)脈沖周期內(nèi)系統(tǒng)位移的噪聲標(biāo)準(zhǔn)差σdx與σdy為0.1 m。觀測(cè)器每隔1個(gè)脈沖周期測(cè)得一組角度與相鄰脈沖到達(dá)時(shí)間差的數(shù)據(jù)。使用300組觀測(cè)數(shù)據(jù)進(jìn)行仿真。

對(duì)于水聲脈沖信號(hào)的DOA估計(jì)，近代聲吶的測(cè)向精度為0.25°～1°[15]。對(duì)于水聲脈沖信號(hào)的TOA估計(jì)，可使用短時(shí)傅里葉變換等時(shí)頻分析方法對(duì)接收的信號(hào)進(jìn)行處理得到。該方法的估計(jì)精度為窗函數(shù)的移動(dòng)步長(zhǎng)，可選取為脈沖信號(hào)脈寬的百分之一[16]，通常水聲脈沖信號(hào)脈寬為毫秒級(jí)至秒級(jí)，因此TOA觀測(cè)噪聲參數(shù)可選取毫秒級(jí)進(jìn)行仿真。本次仿真選取角度觀測(cè)噪聲標(biāo)準(zhǔn)差σβ為1°，相鄰脈沖到達(dá)時(shí)間差噪聲標(biāo)準(zhǔn)差σΔT為5 ms。

在進(jìn)行PLE估計(jì)時(shí)，矩陣R初值的選擇將影響估計(jì)的性能。由于沒有先驗(yàn)信息，σβ、σΔT、rk這3個(gè)參數(shù)的真值是未知的，因此首先通過仿真分析矩陣R初值選擇對(duì)距離估計(jì)的影響。圖2～4為這3個(gè)參數(shù)取不同初值時(shí)PLE距離估計(jì)相對(duì)誤差曲線，仿真過程中σβ真值為 1°、σΔT真值為 5 ms、rk真值為27 km。在進(jìn)行PLE距離估計(jì)時(shí)，依次改變其中1個(gè)參數(shù)的取值，另外兩個(gè)參數(shù)均取真實(shí)值，構(gòu)造R矩陣。每個(gè)參數(shù)改變時(shí)的Monte-Carlo仿真次數(shù)均為100次。

對(duì)比圖2～4的結(jié)果可知，σβ與rk取不同初值時(shí)，距離估計(jì)相對(duì)誤差變化不大；而距離估計(jì)相對(duì)誤差對(duì)σΔT初值的選擇較為敏感。在之后的偽線性估計(jì)仿真中，協(xié)方差矩陣R的初值設(shè)定時(shí)，參數(shù)σβ取為1.5°；σΔT取為 6 ms；距離rk設(shè)定為 81 km。

圖2 不同DOA標(biāo)準(zhǔn)差時(shí)的PLE距離估計(jì)相對(duì)誤差

圖3 不同TOA標(biāo)準(zhǔn)差時(shí)的PLE距離估計(jì)相對(duì)誤差

圖4 不同距離初值時(shí)的PLE距離估計(jì)相對(duì)誤差

圖5和圖6為偽線性距離初值估計(jì)誤差百分比的區(qū)間分布。圖5為未使用滑動(dòng)平均濾波法對(duì)脈沖參數(shù)進(jìn)行預(yù)處理后的估計(jì)結(jié)果，圖6為使用參數(shù)預(yù)處理的結(jié)果，Monte-Carlo仿真次數(shù)為100次。

圖5 未進(jìn)行平滑預(yù)處理的預(yù)測(cè)初始距離與實(shí)際距離的比值分布情況

圖6 進(jìn)行平滑預(yù)處理后的預(yù)測(cè)初始距離與實(shí)際距離的比值分布情況

對(duì)比圖 5、6可見，未進(jìn)行參數(shù)預(yù)處理時(shí)，僅15%比例數(shù)據(jù)的距離估計(jì)誤差低于30%；進(jìn)行參數(shù)預(yù)處理后，66%比例數(shù)據(jù)的距離估計(jì)誤差低于30%。使用滑動(dòng)平均濾波法進(jìn)行數(shù)據(jù)預(yù)處理，可以減小偽線性初始距離估計(jì)的誤差。

改變用于PLE初值估計(jì)的組數(shù)，對(duì)比進(jìn)行數(shù)據(jù)平滑處理前后系統(tǒng)跟蹤的收斂率變化情況，如表 1所示。從表中可知，不進(jìn)行平滑預(yù)處理時(shí)，至少需要30組觀測(cè)數(shù)據(jù)進(jìn)行PLE估計(jì)才能達(dá)到90%以上的收斂率；而進(jìn)行平滑預(yù)處理后，利用 20組觀測(cè)數(shù)據(jù)即可達(dá)到 90%以上的收斂率。因此選取前 20組數(shù)據(jù)進(jìn)行PLE初值估計(jì)。

表1 系統(tǒng)跟蹤收斂率與初值選取組數(shù)的關(guān)系

若N為Monte-Carlo仿真次數(shù)，卡爾曼濾波器距離估計(jì)的均方根誤差（RMSE）為

系統(tǒng)狀態(tài)量Xk跟蹤精度的克拉美羅下界（CRLB）的公式為

公式詳細(xì)推導(dǎo)參見文獻(xiàn)[11]。

為了對(duì)比不同方法距離估計(jì)的性能，根據(jù)對(duì)系統(tǒng)狀態(tài)量Xk估計(jì)的CRLB，可選取距離估計(jì)誤差參考值為[17]：

式中，CRLBk(1,1)和CRLBk(2,2)分別為系統(tǒng)狀態(tài)量xk、yk的CRLB值。

圖7為UKF與MAF-PLE-UKF方法距離估計(jì)的RMSE對(duì)比圖。UKF通過空間幾何三點(diǎn)法估計(jì)目標(biāo)初始位置[18]。100次Monte-Carlo仿真中，UKF收斂率為 74%，而 MAF-PLE-UKF方法收斂率為93%，較 UKF提升了 19%收斂率，且 MAF-PLEUKF方法的距離估計(jì)性能也明顯優(yōu)于UKF。

圖7 UKF與MAF-PLE-UKF方法距離估計(jì)的RMSE對(duì)比圖

圖8～9分別為MAF平滑預(yù)處理前后PLE-UKF與 PLE-EKF算法的距離估計(jì)性能對(duì)比圖。由于選取了前20組數(shù)據(jù)進(jìn)行PLE初值估計(jì)，因此從第21組進(jìn)行系統(tǒng)測(cè)距的性能分析。

圖8 未進(jìn)行數(shù)據(jù)平滑處理的PLE-EKF方法與PLE-UKF方法距離估計(jì)的RMSE對(duì)比圖

圖9 進(jìn)行數(shù)據(jù)平滑處理后的PLE-EKF方法與PLE-UKF方法距離估計(jì)的RMSE對(duì)比圖

圖8為未進(jìn)行數(shù)據(jù)平滑處理的情況。在第100組時(shí)，rek為 626.3 m，PLE-EKF求得的距離估計(jì)RMSE為 2048 m，而 PLE-UKF求得的距離估計(jì)RMSE為1111 m。在第200組時(shí)，rek為337.1 m，PLE-EKF求得的距離估計(jì)RMSE為 826.8 m，而PLE-UKF求得的距離估計(jì)RMSE為475.7 m。在收斂率方面，PLE-EKF收斂率為63%，而PLE-UKF收斂率為 77%。在未進(jìn)行數(shù)據(jù)平滑處理時(shí)，PLE-UKF收斂性能明顯優(yōu)于 PLE-EKF，有更好的寬容性。

圖9為進(jìn)行MAF平滑預(yù)處理后的情況。在第100組時(shí)，rek為626.3 m，PLE-EKF求得的距離估計(jì)RMSE為1318 m，而PLE-UKF求得的距離估計(jì)RMSE為938.6 m。在第200組時(shí)，rek為337.1 m，PLE-EKF求得的距離估計(jì)RMSE為 506.3 m，而PLE-UKF求得的距離估計(jì)RMSE為401.9 m。在收斂率方面，PLE-EKF收斂率為86%，而PLE-UKF收斂率為93%。因此，PLE-UKF較PLE-EKF有更快的收斂速度、更高的收斂率與更好的寬容性。

4 結(jié)束語

本文針對(duì)水聲脈沖信號(hào)被動(dòng)定位問題，提出了使用MAF-PLE-UKF聯(lián)合方法進(jìn)行位置距離估計(jì)。該方法解決了在參數(shù)估計(jì)誤差較大時(shí)，傳統(tǒng) TMA方法容易發(fā)散、收斂性能差等問題。經(jīng)過數(shù)據(jù)仿真實(shí)驗(yàn)，相比于傳統(tǒng)UKF方法，該方法提升了19%的收斂率，并有效降低距離估計(jì)誤差。而且，UKF算法在收斂速度、收斂率和寬容性上也明顯優(yōu)于EKF算法。本文的處理方法能夠有效提升水聲脈沖信號(hào)的被動(dòng)測(cè)距性能。利用PLE進(jìn)行距離初值估計(jì)時(shí)，系統(tǒng)協(xié)方差矩陣R中σΔT參數(shù)的設(shè)定對(duì)算法初始距離估計(jì)影響較大，需要進(jìn)一步的研究，以改善系統(tǒng)性能。