單水聽器淺海聲源被動定位方法研究

石海杰,李京華,常虹

(1．西北工業(yè)大學電子信息學院， 陜西西安 710072； 2．中北大學信息與通信工程學院， 山西太原 030051；3．西安郵電大學通信與信息工程學院， 陜西西安 710121)

自冷戰(zhàn)結(jié)束以后,對海洋的關(guān)注更多轉(zhuǎn)向近岸淺海,由于干涉結(jié)構(gòu)和環(huán)境參數(shù)的非穩(wěn)定性,很難直接利用接收的聲信號進行目標定位,需要充分考慮聲場模型和聲源的先驗信息。匹配場定位是典型的基于模型的方法,自誕生以來相關(guān)方法不斷涌現(xiàn)[1-4],但該方法面臨的最大挑戰(zhàn)是模型失配問題,特別是在淺海環(huán)境下,建立確定意義上的聲傳播模型通常是很困難的[5]。同時,為了獲取更多空間信息,匹配場定位往往需要采用陣列接收信號,而在隱蔽低耗平臺條件下對運動目標定位時,需要的是一種利用單個水聽器就能實現(xiàn)定位的方法。李曉曼等[6]利用warping變換分離寬帶脈沖信號的簡正波,根據(jù)群延遲理論實現(xiàn)單個水聽器條件下的目標定位,但實際上,窄帶信號具有更好的傳輸特性。謝亮等[7]分析典型深海環(huán)境中脈沖信號的到達結(jié)構(gòu)特征,提出通過特征參數(shù)匹配搜索進行目標定位的方法,單個水聽器海上實驗的距離估計誤差不大于8%。Zhang等[8]提出一種基于自相關(guān)函數(shù)峰值提取的無源定位方法,利用單個水聽器就能實現(xiàn)定位,但該方法無法適應非穩(wěn)定環(huán)境。Warner等[9]利用單個水聽器，采用基于模型的貝葉斯估計方法實現(xiàn)了水聲參數(shù)和海底參數(shù)估計,但并未涉及定位方法的研究。貝葉斯估計理論是一種基于概率密度模型的方法,參數(shù)估計時以后驗概率的形式給出,由于可以通過時間累計效應來拓展空間信息的缺失,因此可以采用單個水聽器對目標進行定位。Ogiso等[10]利用遞歸貝葉斯濾波方法在室內(nèi)環(huán)境條件下,實現(xiàn)移動機器人自主定位,該方法屬于高精度定位,但并不適用于水中被動定位情況。王彪等[11]將貝葉斯理論應用到多目標方向估計問題當中,利用多矢量稀疏貝葉斯學習算法重構(gòu)信號空間譜,仿真結(jié)果表明該方法相對于傳統(tǒng)的算法具有更高的空間分辨率和估計精度。Ren等[12]采用貝葉斯濾波方法,利用矢量傳感器測得的船舶噪聲數(shù)據(jù)跟蹤海底環(huán)境參數(shù),該方法可以有效解決環(huán)境參數(shù)非平穩(wěn)變化的問題。Chu等[13]提出一種基于非同步陣列測量的貝葉斯推斷方法細化點陣函數(shù),實驗結(jié)果表明該方法即使在拉普拉斯噪聲干擾下也有較高的定位精度。Guo等[14]將貝葉斯方法用于冰下信道估計,利用估計的信道對脈沖噪聲進行迭代估計和消除,仿真結(jié)果表明,所提出的迭代算法具有較高的估計精度。Sun等[15]將貝葉斯技術(shù)應用于速度跟蹤,采用統(tǒng)計推斷方法對未知非平穩(wěn)噪聲的統(tǒng)計特性進行動態(tài)跟蹤,提出的算法可以有效地解決卡爾曼濾波對初始測量噪聲方差過于依賴的問題。文獻[16-18]應用貝葉斯理論對環(huán)境聚焦以期解決匹配場定位時的失配問題,其實質(zhì)是將未知環(huán)境參數(shù)和聲源位置進行聯(lián)合優(yōu)化估計,該類方法解決了環(huán)境失配問題卻增加了算法復雜程度。總之,在淺海環(huán)境條件下,聲信號在傳播過程中受到環(huán)境的影響尤為嚴重,多種狀態(tài)的海面和海底參數(shù)的相互作用導致了具有復雜空間干涉結(jié)構(gòu)的聲場。而當?shù)皖l線譜聲源移動時,聲場的空間復雜調(diào)制特性轉(zhuǎn)變成固定點接收信號的時間幅度調(diào)制特性。

本文提出的基于貝葉斯估計的目標定位方法是一種基于模型的方法,不直接利用聲場環(huán)境參數(shù)進行匹配,也不聚焦估計聲場環(huán)境參數(shù),而是利用聲場概率密度模型消除非穩(wěn)定因素,參數(shù)估計結(jié)果采用后驗概率的形式給出,能夠有效克服隨機環(huán)境對定位結(jié)果的影響,而且能夠通過時間累計效應拓展空間信息的缺失,實現(xiàn)單個水聽器對目標進行定位,因此基于貝葉斯估計水聲目標定位方法是本文主要研究內(nèi)容。

1 隨機聲場統(tǒng)計模型

為了解決非穩(wěn)定性環(huán)境參數(shù)導致的模型失配問題,一種有效的方法是利用信道的統(tǒng)計特性,以概率密度函數(shù)描述聲場傳輸模型。概率密度函數(shù)估計有參數(shù)化方法和非參數(shù)化方法,核密度估計(kernel density estimation,KDE)是典型的非參數(shù)化方法。KDE法不事先設定隨機過程的分布規(guī)律,而是從樣本本身出發(fā)研究數(shù)據(jù)分布特征,所以理論上KDE法適用于任何隨機過程概率密度函數(shù)的求解。

(1)

式中：κ(x-xi;w),i=1,2,…,n為核函數(shù),w為核函數(shù)寬度；n為核函數(shù)個數(shù)。核函數(shù)可以是均勻核函數(shù)、三角核函數(shù)、高斯核函數(shù)等。理論上任何平滑形式的峰值函數(shù)均可作為KDE的核函數(shù)來使用,只要滿足函數(shù)曲線下方的積分面積等于1即可。根據(jù)(1)式可知,估計結(jié)果實際上是多個樣本數(shù)據(jù)點的核函數(shù)曲線的疊加,由于鄰近樣本點之間會發(fā)生波形疊加,因此最終所形成的曲線形狀與選擇的核函數(shù)關(guān)系并不密切(類比任意波形信號都可以很好地由正弦波、方波或者三角波正交基合成)。考慮到函數(shù)在波形合成計算上的易用性以及高斯分布的普遍適用性,選擇使用高斯曲線作為KDE的核函數(shù)。高斯核函數(shù)可表示為

(2)

每個核函數(shù)都是以樣本點xi為均值,以w為方差的高斯分布。核函數(shù)寬度w能夠控制結(jié)果曲線的整體平坦程度,也就是樣本點數(shù)據(jù)在概率密度曲線形成過程中所占比重。w越寬,樣本點在最終形成的曲線形狀中作用范圍越大,曲線越平坦。反之w越窄,樣本點在最終形成的曲線形狀中作用范圍越小,結(jié)果曲線越陡峭。實際上,這也說明估計得到的概率密度是許多樣本點分布直方圖平滑的結(jié)果,這就是核密度估計方法確定概率密度函數(shù)的核心思想。

核函數(shù)寬度w按(3)式確定[19]

(3)

式中,σx為隨機變量x的標準差,通過(4)式確定

(4)

(5)

2 貝葉斯估計定位方法建模與求解

2.1 淺海聲場信號傳輸幾何模型

如圖1所示,聲源S以徑向速v[n]運動,n為離散時間,z[n],r[n]分別為聲源深度和距離,Q為靜止不動的水聽器。聲源發(fā)射的聲信號為

a[n]=As[n]cos(2πf0n+φ0)

(6)

式中,As[n]為聲源聲壓幅度。聲源狀態(tài)矢量s[n]包含4個分量:r[n],v[n],z[n],As[n],目標定位通過估計r[n],z[n]實現(xiàn)。水聽器接收的信號為

y[n]=A(s[n])cos(2πfdn+φo)+w[n]

(7)

其中,A(s[n])為接收信號聲壓幅度。傳輸衰減用x[n]表示,則對于已知固定位置的水聽器,x[n]由(8)式?jīng)Q定

A(s[n])=As[n]x[n]

(8)

圖1 聲場信號傳輸幾何模型

2.2 聲源狀態(tài)矢量貝葉斯估計建模

基于貝葉斯估計的水聲定位方法,就是在已知水聽器位置的條件下,利用聲場的概率密度模型和接收信號y[n],確定聲源狀態(tài)矢量s[n]。根據(jù)聲場概率密度模型可以確定以聲源狀態(tài)矢量為條件的概率密度f(x[n]|s[n]),也稱先驗概率密度,再結(jié)合接收信號y[n]提取的信號幅度A(s[n]),利用貝葉斯估計方法求解后驗概率,從而估計出聲源位置參數(shù)。若以后驗概率的最大化函數(shù)作為參數(shù)估計的求解條件,則聲源狀態(tài)矢量估計值為

(9)

根據(jù)貝葉斯公式,可得

f(s[n]|x[1],…,x[n])=

(10)

根據(jù)條件概率公式和乘法公式,(10)式變形為

(11)

(11)式中各函數(shù)含義如下:

1)f(s[n]|x[1],…,x[n])是要求解的后驗概率,表示當聲場中某一固定點接收到信號為x[1],…,x[n]時,聲源狀態(tài)矢量s[n]的概率密度。

2)f(x[n]|s[n],x[1],…,x[n-1])表示聲源狀態(tài)矢量為s[n],在聲場中某個固定點已經(jīng)接收到x[1],…,x[n-1]條件下,再次接收到x[n]的概率密度。聲場中某點接收的信號由聲源狀態(tài)矢量和傳輸模型決定,與該點之前接收的信號無關(guān),因此

f(x[n]|s[n],x[1],…,x[n-1])=f(x[n]|s[n])

(12)

3)f(s[n]|x[1],…,x[n-1])表示當聲場中某一固定點接收到信號為x[1],…,x[n-1]時,聲源狀態(tài)矢量s[n]的概率密度,實際上就是前一時刻的后驗概率。

4)f(x[n]|x[1],…,x[n-1])表示不考慮聲源狀態(tài)情況下,聲場中某點已經(jīng)接收到x[1],…,x[n-1]條件下,再次接收到x[n]的概率密度。由于該函數(shù)與聲源狀態(tài)無關(guān),因此可以假設

(13)

式中,η為常數(shù),滿足后驗概率對狀態(tài)矢量的積分為1。

將(12)～(13)式代入(11)式,整理得

(14)

用全概率公式展開f(s[n]|x[1],…,x[n-1]),得:

(15)

式中,f(s[n]|s[n-1])=f(s[n]|s[n-1],x[1],…,x[n-1])是狀態(tài)更新概率密度,狀態(tài)更新過程有馬爾科夫性。

至此,得到了遞歸形式的后驗概率模型,如(15)式所示,也是基于貝葉斯估計的水聲定位方法的模型,也稱為貝葉斯濾波公式,由于積分運算的存在,遞歸過程很難求得封閉形式解,需要研究有效的求解方法。

2.3 直方圖濾波法求解貝葉斯估計問題

直方圖濾波是貝葉斯濾波公式的離散化實現(xiàn)方法,將聲源狀態(tài)矢量離散化成多維網(wǎng)格,認為單一固定網(wǎng)格內(nèi)狀態(tài)矢量的概率密度為常數(shù)。相當于利用離散的概率直方圖代替連續(xù)的狀態(tài)矢量概率密度函數(shù),因此稱為直方圖濾波。

假設將聲源狀態(tài)各分量離散成Kr,Kv,Kz和KA點,聲源狀態(tài)矢量總網(wǎng)格數(shù)為K=KrKvKzKA。設第k個網(wǎng)格的狀態(tài)矢量為sk,則貝葉斯濾波公式(15)變成離散形式

(16)

(16)式中各函數(shù)含義如下:

1)P(sk[n]|x[1],…,x[n])表示當聲場中某一固定點接收到信號x[1],…,x[n]時,聲源狀態(tài)矢量為s[n]=sk的概率。

2)P(x[n]|sk[n])表示在n時刻,聲源狀態(tài)矢量為sk條件下,聲場中對應點處接收到x[n]的概率。

3)P(sk[n]|sj[n-1])表示狀態(tài)矢量轉(zhuǎn)移概率,在整個離散狀態(tài)矢量空間內(nèi),k和j的取值都可以是1～K,因此狀態(tài)轉(zhuǎn)移概率有K2個可能取值。

運動目標的狀態(tài)更新過程可用(17)式表示

s[n]=Ans[n-1]+εn

(17)

式中,An為狀態(tài)轉(zhuǎn)移矩陣

(18)

對于離散的狀態(tài)矢量,狀態(tài)轉(zhuǎn)移概率為

(19)

4)P(sj[n-1]|x[1],…,x[n-1])表示前一時刻聲源狀態(tài)矢量的后驗概率,是前一時刻迭代的輸出,也是后一時刻迭代的輸入。η1是一個常數(shù),滿足后驗概率對狀態(tài)矢量的求和為1。

后驗概率公式遞推至n=1時刻,需要聲源狀態(tài)矢量的初始分布作為迭代公式的輸入,均勻分布可以很好地表示狀態(tài)矢量的初始分布,因此

(20)

以上推導過程,將貝葉斯估計問題進行了離散化,(16)式就是貝葉斯濾波問題的數(shù)值解,該方法解決了狀態(tài)矢量積分難于求解的問題。

2.4 分級網(wǎng)格劃分直方圖濾波法

觀察直方圖濾波算法,不難發(fā)現(xiàn),網(wǎng)格的劃分方法是影響直方圖濾波效果的關(guān)鍵,網(wǎng)格劃分過大,估計精度必然不高,網(wǎng)格劃分過小,算法運算量大增,針對此問題,提出一種分級網(wǎng)格劃分的直方圖濾波方法,該方法在相同網(wǎng)格劃分精度情況下,可以大大降低算法的運算量。

圖2 網(wǎng)格劃分方法對比圖

以Kr=16為例,單級網(wǎng)格劃分方法需要16次后驗概率估計運算和16次搜索運算,而采用二級網(wǎng)格劃分方法則只需要8次后驗概率估計運算和8次搜索運算,在相同的預期估計精度情況下,改進方法效率提高一倍。以Kr=100為例,單級網(wǎng)格劃分方法需要100次后驗概率估計運算和100次搜索運算,而采用二級網(wǎng)格劃分方法則只需要20次后驗概率估計運算和20次搜索運算,在相同的預期估計精度情況下,改進方法效率提高5倍。隨著網(wǎng)格劃分數(shù)的增加,效率提高的效果將會更加顯著。

3 基于SWellEx-96數(shù)據(jù)的算法驗證

3.1 SWellEx-96實驗說明

SWellEx-96系列實驗是1996年5月在加利福尼亞圣地亞哥海岸附近進行的,其中的S5實驗條件為：拖曳深度為9 m的聲源以2.5 m/s速度運行,采用海底布放的水平陣列接收(south horizontal line array,SHLA)。由于實驗信息記錄完整,GPS數(shù)據(jù)準確,實驗數(shù)據(jù)質(zhì)量較高,因此非常適合于淺海環(huán)境條件下目標定位方法的驗證。實驗海域海深約200 m,SHLA布放深度198 m,本文利用SHLA中的3個陣元接收的數(shù)據(jù)進行驗證,聲源含多個頻率的線譜,本文利用最低頻率109 Hz信號,詳細的環(huán)境參數(shù)可參考http:∥swellex96.ucsd.edu/index.htm。

3.2 核密度估計法隨機聲場建模結(jié)果

按SWellEx-96實驗提供的51組聲速剖面數(shù)據(jù),構(gòu)成一組輸入?yún)?shù),假設海水深度服從范圍180～200 m、間隔為1 m的隨機均勻分布,由21個深度組成另外一組輸入?yún)?shù),2組參數(shù)任意組合,構(gòu)成包含1 071個樣本的隨機環(huán)境參數(shù),利用Kraken程序仿真1 071個對應的聲場分布樣本,KDE方法根據(jù)這些聲場樣本可以估計聲場概率密度模型。圖3所示為聲場中一些點的概率密度函數(shù)結(jié)果展示,其中包含海面近場、靠近聲源(海底近場)、海面遠場和海底遠場的估計結(jié)果,這些是典型的特殊位置,也包括典型的一般位置的估計結(jié)果。從樣本的綠色柱狀直方圖對比可以看出:①不同位置的直方圖是不同的,同一區(qū)域的直方圖也不同,特殊位置與一般位置的直方圖都是隨機的,沒有觀察到可遵循的一般規(guī)律;②近場、特別是靠近聲源位置的直方圖更不規(guī)律,更容易出現(xiàn)多峰情況,說明這些位置的聲場信號的隨機性更強;③在仿真的過程中發(fā)現(xiàn)臨近位置的直方圖區(qū)別可能也很大,這是由聲場干涉結(jié)構(gòu)決定的,干涉結(jié)構(gòu)的紋理越細,聲場隨機性越強,這也再次說明只有基于統(tǒng)計特性的信號處理方法才能適用于此處。

圖3中的紅色劃線表示核密度方法估計的概率密度函數(shù),與高斯模型和最大熵估計得到的概率密度模型對比可以看出:高斯模型法得出的概率密度函數(shù)始終保持高斯函數(shù)固有的鐘形脈沖形式,只在樣本分布為單峰情況下勉強適用;最大熵法好于高斯模型法,但對多峰分布的適應能力也不強;核函數(shù)方法估計出的概率密度函數(shù)在不同情況下都有很好的適應性,即使在靠近聲源位置,直方圖分布具有多個峰值的情況下,估計結(jié)果與樣本分布也匹配得很好,所以,本文采用核密度方法估計聲場概率密度函數(shù)f(x)。將聲源狀態(tài)矢量參數(shù)帶入f(x),就可以求得先驗概率P(x[n]|s[n])。

圖3 聲場概率密度函數(shù)估計結(jié)果

3.3 SWellEx-96實驗數(shù)據(jù)預處理

為了驗證單個水聽器能夠?qū)崿F(xiàn)定位,以陣元1、陣元17和陣元32作為3個獨立水聽器,提取它們的數(shù)據(jù)進行前置濾波(抗混疊)、降采樣、混頻、濾波處理,將最終提取的經(jīng)過傳輸調(diào)制的幅度值,作為A(s[n])的估計值,用以為后驗概率估計時備用。圖4所示為預處理結(jié)果圖,原始數(shù)據(jù)有50 min時長,選取了5 min進行展示和預處理,可以看出1號陣元信號幅度最小,32號陣元信號幅度最大,這種規(guī)律是由陣元位置關(guān)系決定的。原始數(shù)據(jù)頻域圖中能夠很清晰地看到有109 Hz和112 Hz信號的存在,由于聲源是移動目標,因此線譜被展寬了。處理后時域、頻域信號可以看出109 Hz信號被有效地提取出來。

圖4 實測數(shù)據(jù)預處理結(jié)果

3.4 基于貝葉斯估計的水聲定位結(jié)果

聲源狀態(tài)矢量按表1取值,確定各聲源狀態(tài)矢量的范圍和網(wǎng)格劃分數(shù)K,聲源幅度為歸一化范圍,因此無單位。信號預處理時將實測數(shù)據(jù)采樣率降為1 000 Hz,因此T=0.001 s,狀態(tài)更新噪聲參數(shù)按照表2取值。根據(jù)(19)和(20)式可以確定狀態(tài)轉(zhuǎn)移概率和初始狀態(tài)分布。單次定位取20個數(shù)據(jù)進行迭代,每隔1 min定位1次。

表1 聲源狀態(tài)矢量參數(shù)表

表2 狀態(tài)更新參數(shù)表

定位結(jié)果如圖5和圖6所示,圖5是三陣元分別在3個時刻的定位結(jié)果,圖中藍色圓圈中心位置代表聲源的真實位置,紅色十字中心位置是該時刻的定位結(jié)果,灰度圖表示該時刻目標在定位區(qū)域內(nèi)后驗概率分布。圖6是定位結(jié)果與真實值隨時間變化的對比情況。綜合觀察圖5和圖6的定位結(jié)果圖可以得出如下結(jié)論:①表示概率分布的灰度圖是以一個具有較大值區(qū)域(灰度較暗的區(qū)域)的方式跟進目標真實位置,最大概率值點與該區(qū)域有較高的一致性,而非某些突兀的概率奇點作為定位結(jié)果,這能說明本文方法是一種具有較高可靠性的方法;②開始時刻的定位精度沒有后續(xù)時刻定位精度高,這是由于開始時刻聲源狀態(tài)矢量采用均勻分布方式表示,后續(xù)時刻聲源狀態(tài)矢量是前一時刻的后驗概率,隨著時間的推進,定位結(jié)果將得到不斷糾正,越來越接近真實目標位置;③距離定位結(jié)果的穩(wěn)定性好于深度定位結(jié)果,這是由于聲場模型干涉結(jié)構(gòu)在距離方向的尺度遠大于深度方向的尺度,尺度越小干涉越復雜,定位越困難,穩(wěn)定性越差。實際上,淺海環(huán)境目標的深度定位一直是一個難以解決的問題。從實測數(shù)據(jù)的定位結(jié)果來看,本文方法的可取之處在于距離方向定位的穩(wěn)定性,更在于能夠較準確估計目標的深度。

圖5 各陣元不同時刻定位結(jié)果

圖6 定位結(jié)果隨時間變化情況

3.5 分級網(wǎng)格劃分直方圖濾波法效率仿真

圖7所示為單級網(wǎng)格劃分與分級網(wǎng)格劃分實測效率對比圖,仿真時以普通辦公用計算機為硬件平臺,以仿真工具Matlab為軟件平臺。網(wǎng)格化參數(shù)K取值:[1,25,100,400,2 500,4 900,6 400,10 000],迭代次數(shù)為50,圖中的橫坐標為K值取對數(shù),縱坐標為相對時間取對數(shù),ΔT為K值取1時單次迭代的時間開銷。觀察圖7的效率對比圖可以得出如下結(jié)論:①分級網(wǎng)格劃分直方圖濾波法對算法效率改善是非常明顯的,網(wǎng)格化參數(shù)K取值越大,效率提高效果越明顯,這一點與理論分析是一致的;②當網(wǎng)格化參數(shù)K取值較小時(小于25次),算法總的開銷時間還會受到其他因素 (包括臨時變量初始化、概率密度函數(shù)裝載、后驗概率空間分配等)影響,隨著K值的增加,算法開銷時間主要決定于后驗概率的迭代求解過程,K值越大,這種決定性越強。

圖7 分級網(wǎng)格劃分法效率對比

3.6 基于貝葉斯估計目標定位方法的精度

定義(21)式為距離均方根誤差,定義(22)式為深度均方根誤差

(21)

(22)

可以得到表3的結(jié)果,表中均值是指3個陣元均方根誤差的平均值。說明本文研究的定位方法,深度定位精度可達到35 m,距離定位精度可以達到0.69 km。考慮做均方根誤差的樣本包含定位開始階段偏差較大的定位結(jié)果,因此實際定位效果要更好。r[n]范圍為10 km,直觀估計目標距離估計誤差在7%以內(nèi)。與文獻[18](6%)、文獻[7](8%)、文獻[6](10%)方法相比,本文方法距離定位結(jié)果總體略好,并且效率高出許多。由于合適的深度定位精度對比的文獻不多,所以沒有進行深度定位精度的橫向比較。

表3 定位結(jié)果均方根誤差表

4 結(jié) 論

貝葉斯濾波原理的實質(zhì)是試圖用所有已知信息來構(gòu)造系統(tǒng)狀態(tài)變量的后驗概率密度,即用系統(tǒng)模型預測狀態(tài)的先驗概率密度,再使用最近的量測值進行修正,得到后驗概率密度。這樣通過觀測數(shù)據(jù)的不斷更新來遞推計算狀態(tài)參數(shù),由此獲得狀態(tài)的最優(yōu)估計。本文將貝葉斯估計理論應用于水聲信號的跟蹤與定位算法中,根據(jù)淺海條件和單個水聽器接收條件,利用聲場環(huán)境先驗信息和水聽器接收的聲場數(shù)據(jù),通過統(tǒng)計算法確定聲源狀態(tài)矢量后驗概率,充分利用環(huán)境的不確定性和水聲傳輸模型的先驗知識較好地實現(xiàn)了聲源目標被動跟蹤和定位。