基于貝葉斯估計的淺海聲源被動定位方法

石海杰, 李京華, 劉麗麗, 常虹

(1.西北工業大學 電子信息學院, 陜西 西安 710129; 2.西安郵電大學 通信與信息工程學院, 陜西 西安 710121)

0 引言

在淺海環境中,多途干涉和環境參數易變不利于目標定位,需要充分利用聲場傳輸模型和聲源先驗信息改善定位效果。匹配場定位是基于聲場傳輸模型的方法,自誕生以來,相關研究不斷涌現[1-7],但該方法面臨的最大問題仍然是模型失配問題,特別是在淺海條件下,建立確定性意義上的聲傳播模型并不容易[8]。同時,為了提高空間信息獲取能力,匹配場定位常采用陣列形式接收信號,而在隱蔽性要求較高的情況下,需要的是一種利用單個水聽器就能實現定位的方法。Warping變換能夠提取寬帶脈沖不同階簡正波信號,根據群延遲理論可以實現單個水聽器條件下的目標定位[9],但實際上,窄帶信號具有更好的傳輸特性。Zhang等[10]提出一種基于自相關函數峰值提取的無源定位方法,該方法不需要估計特征聲線的相對延遲,利用單個水聽器就能實現定位,但無法適應非穩定環境。 Warner等[11]利用單個水聽器采用基于模型的貝葉斯估計方法實現了水聲參數的反演。Ren等[12]采用貝葉斯濾波方法,跟蹤海底環境參數,高飛等[13]利用貝葉斯方法進行地聲參數反演,該類方法能有效地解決環境參數非平穩變化的問題。Guo等[14]提出一種基于稀疏貝葉斯學習的魯棒水聲信道估計算法,用于冰下信道參數估計,仿真結果表明,所提出的迭代算法具有較高的估計精度。以上基于貝葉斯估計的水聲信號處理方法,雖未涉及定位方面的研究,但該理論卻是適合于水聲定位的。貝葉斯估計理論是一種基于統計模型的參數估計方法,由于可以通過時間累計效應來拓展空間信息的缺失,可以采用單個水聽器對目標進行定位。Ogiso等[15]利用遞歸貝葉斯濾波方法在室內環境條件下,利用方向標輔助方法,實現移動機器人自主定位。王彪等[16-17]將貝葉斯理論應用到多目標方向估計問題中,利用多矢量稀疏貝葉斯學習算法重構信號空間譜,建立未知稀疏源信號的壓縮感知模型,仿真結果表明該方法具有較高的空間分辨率和估計精度。

總之,在淺海環境條件下,聲信號在傳播過程中受到環境的影響尤為嚴重,多種狀態的海面和海底參數的相互作用導致了具有復雜空間干涉結構的聲場。而當低頻窄帶聲源移動時,聲場的空間復雜調制特性轉變成固定點接收信號的時間幅度調制特性。貝葉斯估計理論是一種利用概率密度函數描述聲場模型的方法,參數估計結果采用后驗概率的形式給出,能夠有效克服隨機環境對定位結果的影響,并通過時間累計效應來拓展空間信息的缺失,實現單個水聽器對目標進行定位。因此,基于貝葉斯估計的水聲目標定位方法的建模與求解是本文的主要研究內容。

1 淺海水下隨機聲場建模

1.1 信號傳輸幾何模型

信號傳輸幾何模型如圖1所示。圖1中,H為海水深度,d為水聽器布放深度,ρw、cw為海水密度和聲速,ρb、cb為海底密度和聲速,聲源S以徑向速度v[n]運動,n為離散時間,z[n]、r[n]分別為聲源深度和距離,a[n]為聲源發射的聲信號,y[n]為水聽器接收到的信號,Q為靜止不動的水聽器。

圖1 信號傳輸幾何模型Fig.1 Geometric model of signal transmission

聲源發射的聲信號為

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

(1)

式中:As[n]為聲源聲壓幅度;f0為原始頻率;φ0為初相。聲源狀態矢量s[n]包含四個分量:r[n]、v[n]、z[n]、As[n],目標定位通過估計r[n]、z[n]實現。水聽器接收的信號為

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

(2)

式中:A(s[n])為接收信號聲壓幅度;fd為接收信號頻率;w[n]為環境噪聲。傳輸衰減用(傳輸模型)x[n]表示,則對于已知固定位置的水聽器:

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

(3)

1.2 聲場概率密度模型

利用信道的統計特性,以概率密度函數描述聲場傳輸模型,是解決環境參數非穩定的有效方法。概率密度函數估計有參數化方法和非參數化方法,核密度估計(KDE)是典型的非參數化方法。KDE法不事先設定隨機過程的分布規律,是一種從樣本本身出發研究數據分布特征的方法,因此理論上KDE法適用于任何復雜隨機過程的概率密度函數的求解。假設利用蒙特卡洛法,根據樣本數據統計獲得的未知概率密度函數為f(x)(樣本分布直方圖),KDE方法的目的就是按式(4)構造一個概率密度函數,盡量去接近f(x)。

(4)

式中:m為核函數個數;κ(x-xi;w),i=1,2,…,m為核函數,xi為隨機樣本,w為核函數寬度。核函數可以是均勻核函數、三角核函數、高斯核函數等。理論上任何平滑形式的峰值函數均可作為KDE的核函數,只要滿足函數曲線下方的積分面積為1即可。考慮到函數在波形合成計算上的易用性,也考慮高斯分布的普遍適用性,因此選擇使用高斯曲線作為KDE的核函數。高斯核函數可表示為

(5)

核函數都是以樣本點xi為均值、以w為方差的高斯分布。核函數寬度w能夠控制結果曲線的整體平坦程度,也就是樣本點數據在概率密度曲線形成過程中所占的比重。w越寬,樣本點在最終形成的曲線形狀中作用范圍越大,曲線就越平坦。反之,w越窄,樣本點在最終形成的曲線形狀中作用范圍越小,結果曲線就越陡峭。實際上,這也表明估計得到的概率密度是許多樣本點分布直方圖平滑的結果,這也是核密度估計方法適用廣泛的原因。

核函數寬度w按式(6)確定[18]:

(6)

式中:σx為隨機變量x的標準差。

2 基于貝葉斯估計的目標定位方法

2.1 貝葉斯估計水聲目標定位方法建模

基于貝葉斯估計的水聲目標定位方法,就是在已知水聽器位置的條件下,利用聲場的概率密度模型和接收信號y[n],確定聲源狀態矢量s[n]。根據聲場概率密度模型可以確定以聲源狀態矢量為條件的概率密度f(x[n]|s[n]),即先驗概率密度,再結合接收信號y[n]提取的信號幅度A(s[n]),利用貝葉斯濾波方法求后驗概率f(s[n]|x[1],…,x[n])。 采用后驗概率最大化函數作為參數估計的求解條件,則聲源狀態矢量估計值為

(7)

根據貝葉斯公式,有

(8)

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

(9)

式中各函數含義如下:

①f(s[n]|x[1],…,x[n])為聲源狀態矢量后驗概率,表示接收到信號為x[1],…,x[n]時聲源狀態矢量s[n]的概率密度,概率密度最大位置就是聲源位置。該函數就是貝葉斯估計方法最終要求的結果,反映的實質是利用接收信號確定聲源狀態。

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

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

(10)

式中:f(x[n]|s[n])稱為先驗概率,通過聲場概率密度模型確定。該函數反映的實質是在不同聲源狀態條件下,能夠接收到x[n]的可能性有多大。

③f(s[n]|x[1],…,x[n-1])表示當接收到信號為x[1],…,x[n-1]時,聲源狀態矢量s[n]的概率密度,實際上就是前一時刻的后驗概率。在貝葉斯估計理論中,該函數即表示前一迭代過程的結果,又是后一迭代過程的輸入,反映了估計結果不斷更新校正的過程。

④f(x[n]|x[1],…,x[n-1])表示不考慮聲源狀態情況下,n-1時刻以前聲場中某點已經接收到x[1],…,x[n-1]條件下,n時刻再次接收x[n]的概率密度。由于該函數與聲源狀態無關,根據無源波動方程(齊次)可知該函數應為一個常數,可以假設

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

(11)

式中:η為一個常數,滿足后驗概率對狀態矢量的積分為1。在進行后驗概率的求解過程中,可以先假設為1,求出后驗概率后,后驗概率歸一化常數的倒數就是η的值。將式(10)、式(11)代入式(9),整理得

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

(12)

用全概率公式展開f(s[n]|x[1],…,x[n-1]),根據狀態更新過程的馬爾可夫性,得

(13)

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

2.2 貝葉斯濾波問題的求解

觀察式(13),可以發現貝葉斯濾波最大的問題在于積分的求解,由于存在時間上的迭代,很難解得封閉形式解。對聲源狀態矢量離散化,通過求和解決積分問題是很好的思路,將聲源狀態矢量離散化成多維網格,認為單一固定網格內的狀態矢量的概率密度為常數。相當于利用離散的概率直方圖代替真實的狀態矢量概率密度函數,因此也稱為直方圖濾波。假設將聲源狀態各分量離散成Kr、Kv、Kz和Ka點,下標r、v、z、a分別對應聲源距離r[n]、速度v[n]、深度z[n]和聲壓幅度As[n],則總網格數為K=KrKvKzKa。設第k個網格的狀態矢量為sk,式(13)變成:

(14)

式中各函數含義如下:

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

②P(x[n]|sk[n])為離散形式先驗概率,表示在n時刻,聲源狀態矢量為sk條件下,聲場中對應點處接收到x[n]的概率。

③P(sk[n]|sj[n-1])為狀態轉移概率,運動目標的狀態更新過程可用如下方程表示:

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

(15)

式中:An為狀態轉移矩陣,

(16)

(17)

當n=1時狀態轉移概率變成聲源狀態矢量的初始分布,可以用均勻分布表示,因此

P(sj[0])=1/K,j=1,…,K

(18)

④P(sj[n-1]|x[1],…,x[n-1])是前一時刻的后驗概率,也是后一時刻迭代的輸入。

⑤η1為一個常數,滿足后驗概率對狀態矢量的求和為1。

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

3 實測驗證與仿真分析

3.1 SWellEx-96實驗說明

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

3.2 聲場概率密度函數估計結果

SWellEx-96實驗提供了51組聲速剖面數據,假設海水深度服從范圍180～200 m、間隔為1 m的均勻分布,聲速剖面和深度任意組合,構成包含1 071個樣本的隨機環境參數,利用Kraken程序仿真 1 071個對應的聲場樣本,KDE方法根據這些樣本可以估計聲場概率密度模型。圖2所示為假設將水聽器固定在(0 km,198 m)處(對應S5實驗中SHLA陣列布放位置),聲源在聲場中不同位置時仿真的聲場結果,聲源幅度為1,聲場聲壓對功率歸一化。其中,圖2(a)為1 071個樣本聲場中的任意一個,圖2(b)、圖2(c)為對1 071個樣本求均值和方差而得到的平均聲場和方差聲場。由仿真結果可以看出聲場中存在由干涉造成的明暗相間的條紋。仿真過程中發現不同環境參數對應的干涉結構不同,而干涉結構的存在會影響到定位結果的精度。從圖2(a)、圖2(b)、圖2(c)的對比中可以看出,平均聲場比樣本聲場要平滑許多,因為單個樣本是有隨機性的。據此現象可以想到,利用測量數據進行目標定位時,獲取的實測數據相當于隨機過程的一個樣本,一定存在隨機性,因此需要一種能夠利用模型的統計特性輔助定位方法來克服這種隨機性;從仿真的方差聲場中可以看出聲源靠近水聽器時所產生的聲場的離散性較大,隨機性較強,這是因為聲源與水聽器距離越近,干涉效應越強,較小的環境參數變化也會造成較明顯的干涉結果,表明兩者距離靠近并不利于提高定位精度。

圖2 歸一化聲壓場仿真結果Fig.2 Simulation results of normalized sound pressure field

圖3所示為聲場概率密度估計結果,其中包含聲源靠近水聽器、聲源在海面遠場這些特殊位置的估計結果,也包括典型的一般位置(聲源在中間深度,中間距離)的估計結果。從樣本的綠色柱狀直方圖對比可以看出:1)不同位置的直方圖是不同的,同一區域的直方圖也不同,特殊位置與一般位置的直方圖都是隨機的,沒有可觀察到的一般規律可遵循;2)近場、特別是聲源靠近水聽器位置處的直方圖更不規律,更容易出現多峰情況,表明這些位置聲場信號的隨機性更強;3)在仿真過程中發現臨近位置的直方圖區別可能也很大,這是由聲場干涉結構決定的,干涉結構的紋理越細,聲場隨機性越強,這也再次表明只有基于統計特性的信號處理方法才能適用于此處。圖3中的紅色劃線表示核密度方法估計的概率密度函數,與高斯模型和最大熵估計得到的概率密度模型對比可以看出:高斯模型法得出的概率密度函數始終保持高斯函數固有的鐘形脈沖形式;最大熵法好于高斯模型法,但對多峰分布的適應能力也不強;核函數方法估計出的概率密度函數在不同情況下都有很好的適應性,即使樣本直方圖分布具有多個峰值,估計結果與樣本分布也匹配得很好,因此本文采用KDE方法進行聲場概率密度函數估計。

圖3 概率密度函數估計結果Fig.3 Estimation results of probability density function

3.3 SWellEx-96實驗數據預處理

為驗證定位算法,以SHLA的陣元01、陣元17和陣元32作為3個獨立水聽器,提取它們的數據進行前置濾波(抗混疊)、降采樣、混頻、濾波處理,最終提取的經過傳輸調制的幅度值,作為A(s[n])的估計值。圖4所示為以陣元01為代表的預處理結果圖,原始數據有50 min時長,選取5 min進行展示,原始數據頻域圖中能夠很清晰地看到有109 Hz和112 Hz信號的存在,由于聲源是移動目標,線譜被展寬了。通過處理后時域、頻域信號可以看出109 Hz信號有效地被提取出來。

圖4 實測數據預處理Fig.4 Pre-processing of measured data

圖5所示為仿真的模型聲場(Kraken)信號與陣元01實測數據聲場信號隨距離衰減的對比圖。由圖5可以看出,實測的真實值相當于隨機過程的一個樣本,與均值相比信號的隨機起伏明顯更大,但都是在最大值和最小值之間變化的,并且能夠看出實測信號與仿真信號衰減的總體趨勢基本吻合,都是由近及遠衰減變大,同時伴有由干涉現象引起的幅度起伏。實測信號隨機性更強還有實際測量環境噪聲的影響,仿真信號并不考慮噪聲,其信號幅度變化會更平順一些。另外可以看出,實測信號的干涉結構和信號幅度與仿真信號并不匹配,特別是遠場信號幅度差別更大。信號幅度的差別主要是因為實測環境的背景噪聲造成的,遠場時信噪比降低,幅度差異更大。正如前文所述,干涉結構的差異是由單次測量的隨機性造成的,這也正是本文采用基于統計理論的貝葉斯估計方法進行目標跟蹤和定位的原因。此仿真結果,一方面表明了聲場模型的有效性,另一方面表明了采用概率密度模型描述聲場的必要性。

圖5 陣元01實測數據傳輸損失Fig.5 Transmission loss of measured data from array element 01

3.4 基于貝葉斯估計的水聲目標定位結果

聲源狀態矢量按表1參數離散化,聲源幅度為歸一化聲壓幅度。信號預處理時將實測數據采樣率降為1 000 Hz,因此T=0.001 s,狀態更新噪聲參數按照表2取值。根據式(17)和式(18)可以確定狀態轉移概率和初始狀態分布,根據概率密度函數可以確定先驗概率,根據式(14)可以進行后驗概率迭代,根據式(7)可以確定目標位置。

表1 聲源狀態矢量參數

表2 狀態更新參數表

定位結果如圖6和圖7所示,圖6為三陣元分別在3個時刻的定位結果,圖中藍色圓圈中心位置代表聲源的真實位置,紅色十字中心位置是該時刻的定位結果,灰度圖表示該時刻目標在定位區域內各點可能出現的后驗概率。圖7為定位結果相對誤差隨時間的變化情況。綜合觀察圖6和圖7可以得出如下結論:1)表示概率分布的灰度圖以一個具有較大值區域(灰度較暗的區域)的方式跟進目標真實位置,最大概率值點與該區域有較高的一致性,而非某些突兀的概率奇點作為定位結果,能表明本文方法是一種具有較好可靠性的方法;2)開始時刻的定位精度沒有后續時刻定位精度高,這是因為開始時刻聲源狀態矢量采用均勻分布方式表示,后續時刻聲源狀態矢量用前一時刻后驗概率表示,隨著時間的推進,定位結果將得到不斷糾正,越來越接近真實目標位置;3)距離定位結果的穩定性好于深度定位結果,這是因為聲場模型干涉結構在距離方向的尺度遠大于深度方向的尺度,尺度越小干涉越復雜,定位越困難,穩定性越差。

圖6 各陣元不同時刻定位結果Fig.6 Localization results of each array element at different times

圖7 定位結果隨時間變化情況Fig.7 Variation of localization results with time

以上結論中1)、2)比較容易理解,3)的理解可以參考圖8和圖9。圖8為將圖2的聲場樣本變換顯示比例,盡量放大距離方向顯示比例,并截取一部分顯示的結果,此時距離方向的縮小比例仍然是深度方向縮小比例的數倍,但從圖像的干涉條紋中已經能夠清晰地看出深度方向的干涉尺度遠小于距離方向的干涉尺度。干涉條紋尺度越小,表明干涉結構越復雜,也越不利于定位。再觀察圖9所示的二維歸一化聲壓傳輸損失,選100 m深度處歸一化聲壓隨距離的傳輸損失和5 km距離處歸一化聲壓隨深度傳輸損失情況,并繪制到同一坐標下,可以看出深度方向干涉結構尺度和距離方向干涉結構尺度的對比關系顯示得更加清晰。

圖8 三維歸一化聲壓傳輸損失Fig.8 Three-dimensional normalized sound pressure transmission loss

圖9 二維歸一化聲壓傳輸損失Fig.9 Two-dimensional normalized sound pressure transmission loss

基于以上分析,貝葉斯估計方法在定位水聲目標時,距離方向定位穩定性要好于深度方向定位穩定性。實際上,淺海環境目標的深度定位一直是一個難以解決的問題。從實測數據的定位結果來看,本文方法的可取之處在于距離方向定位的穩定性,更在于能夠較準確地估計目標的深度。

3.5 基于貝葉斯估計的水聲定位誤差

定義式(19)為距離估計相對誤差,定義式(20)為深度估計相對誤差:

(19)

(20)

式中:rmax-rmin為距離定位范圍,rmax-rmin=10 km;zmax-zmin為深度定位范圍,zmax-zmin=200 m;N為估計次數。根據式(19)和式(20)可以得到表3的結果,表3中均值是指3個陣元相對誤差的平均值,表明本文研究的定位方法,深度定位誤差為12.04%,距離定位定位誤差6.47%,深度誤差大于距離誤差,與之前分析的結果是一致的。由于分析的樣本包含開始階段誤差較大的估計結果,實際定位效果要好于上述結論。

表3 定位結果相對誤差表

4 結論

貝葉斯估計原理的實質是試圖用所有己知信息來構造系統狀態變量的后驗概率密度,即用系統模型預測狀態的先驗概率密度,再使用最近的量測值進行修正,得到后驗概率密度。這樣通過觀測數據的不斷更新來遞推計算狀態參數,由此獲得狀態的最優估計。本文提出的基于貝葉斯估計水聲目標定位方法,主要是用于解決非平穩過程描述和非平穩噪聲去除問題。研究的重要內容之一是獲得有效的先驗概率密度函數,準確地表示不確定性。利用貝葉斯信號處理技術,可以利用這種不確定的知識來改進被動定位性能。另一個重要內容是根據系統的體系結構建立所需的參數估計模型,包括確定未知參數與概率密度函數的關系,還有就是確保觀測量是概率密度函數所描述的隨機過程的樣本。最后一個內容就是后驗概率的求解,也就是迭代算法的設計,通常也稱為濾波算法。通過以上內容的研究,本文設計了一種有效的針對淺海環境中的運動目標定位方法。仿真和實驗結果表明:在淺海環境條件下,200 m×10 km的范圍內,聲源目標水平距離定位精度可達到6.47%,深度定位精度可達12.04%。由于采用單個水聽器就能實現運動目標定位,本文方法可用于隱蔽低耗平臺對目標的定位。