基于直方圖濾波的淺海聲源測距算法研究

劉麗麗, 李京華, 馮曉毅, 石海杰, 張曉彪

(西北工業大學 電子信息學院, 陜西 西安 710072)

水下運動目標的測距一直以來都是水聲學的研究熱點問題,針對水下實現淺海聲源的測距問題[1]。近年來許多研究者利用濾波技術來實現對淺海聲源的跟蹤及定位[2],郭曉樂等[3]提出了一種利用集合卡爾曼濾波算法的聲速剖面跟蹤反演和移動聲源跟蹤定位的方法,實現對聲源的跟蹤定位。高飛等[4]基于單水聽器多途聲到達時差,利用擴展卡爾曼濾波算法對深海淺層移動聲源進行定位研究,以首達波與次達波的聲到達時差為觀測值,構建EKF聲源定位模型,得到移動聲源的深度、水平距離及速度等信息。上述算法都屬于線性高斯目標的狀態估計問題。而對于處理非線性非高斯目標的狀態估計問題,本文提出了利用直方圖濾波來實現對淺海運動聲源的定位。首先,在統計模擬方法的基礎上利用簡正波建模方法對聲源在水下聲傳播進行建模仿真,對由多個環境變量引起的傳播損失進行分析,捕捉環境因素的不確定性；其次研究了統計模擬方法模型產生的聲場預測集合的統計特性,采用概率密度估計方法對聲場在不同水平距離-深度處的幅度分布進行了分析實驗,得出在不同水平距離-深度上聲壓幅度的概率密度分布;最后,利用對聲場模型的概率估計結果以及接收信號等先驗知識利用直方圖濾波實現對聲源的距離估計,并利用實測數據進行實驗以驗證算法有效性。

1 基于統計模擬方法的聲場建模

為了驗證本文測距算法的有效性,本文的仿真實驗采用SwellEx-96[5]實驗中的S5事件的海洋環境參數[6]。海洋環境參數如圖1所示。

圖1 SwellEx-96實驗海洋環境參數示意圖

在圖1中,最上面是水層;水層下方為沉積層;沉積層下方是巖層。在SwellEx-96實驗周圍海域中,在不同的時間段內進行了51次聲速剖面測量,測量結果如圖2所示。

圖2 SwellEx-96實驗期間測量的51次聲速剖面示意圖

如果確切地知道聲速剖面,海水深度和沉積物特征,就可以得到單一環境參數的聲場。然而在實際應用中,建模用的環境參數存在著不確定性,選擇不確定的參數值將導致建模的聲場與實際接收的聲場數據之間存在較大的差異。因此,為了分析環境參數的不確定性對聲場建模產生的影響,本文在單一環境參數簡正波建模的基礎上,使用統計模擬方法進行聲場建模,用SwellEx-96的51種聲速剖面以及不同接收器深度的環境參數條件進行多環境參數建模,從而得到不同環境參數下的聲場模型。

首先,利用簡正波模型實現單一環境參數下的聲場建模,聲源強度S在深度zs處的聲場p(r,z,t,zs)為[7]

(1)

式中,am(Km)是第m階簡正波的水平波數的模式激勵,(1)式表示簡正波建模產生的聲場是m階不同簡正波模式的疊加。依賴聲速曲線和第m階簡正波的激發系數Zm(zs),可得到聲場的形式為

(2)

單一環境參數建模時聲速剖面選取51組聲速剖面的第一組。圖3所示為頻率109 Hz的淺海聲源在水下9 m處的聲場傳播衰減圖。可以看出,聲源在正下方的傳播衰減最小,傳播過程中隨著水平距離和海水深度的增大,傳播衰減也逐漸增大。這是因為海底泥沙等引起的聲反射對其傳播所產生的干擾所導致。

圖3 淺海9 m處聲源的傳播衰減圖(分貝)

其次,進行多環境參數下的簡正波建模,設影響聲場建模的環境參數為se,使用聯合概率密度f(se)表示參數的不確定性。采用統計模擬方法從f(se)中隨機抽取環境參數,實現聲場建模。利用不同的環境參數重復該過程M次,則得到多種環境參數的聲場集合,如表1所示。

表1 基于統計模擬方法的聲場建模算法描述

統計模擬方法描述了不同位置處的聲壓值集合f(p|r,z),聲傳播過程可以看作是從環境參數到接收信號聲壓的映射p(r,z,se),表示為

(3)

若環境參數中有一個是變化的,則聲場的變化也取決于該參數的變化,設環境參數的先驗信息表示為

(4)

式中：μ為均值；σ2為方差；η表示幅度；f(se)表示關于環境參數的先驗信息,聯立(3)至(4)式得

f(p|r,z)=

(5)

本文針對同一個目標聲源,分析了聲速剖面和水聽器深度變化對建模的影響,其中,聲速剖面為圖2中51組聲速剖面曲線,水聽器深度從80～100 m深度范圍(步長為1)變化,使用枚舉法對聲速剖面和水聽器深度的1 071(51×21)種組合進行聲場仿真,結果如圖4所示,得到不同海洋環境參數下聲場的集合。

圖4 統計模擬方法下的簡正波聲學建模結果示意圖

2 基于概率密度估計的接收信號仿真

2.1 水下接收信號模型

水聽器接收聲源信號的示意圖如圖5所示。其中水下聲源為黑色橢圓,水聽器為黑色圓點。聲源到水聽器的水平距離為r[n],聲源深度為z[n]。聲源輻射的信號記為a[n],接收的信號記為x[n],n是離散時間,聲源接觸到海底會發生反射,聲信號的傳播方向用帶箭頭的直線表示。

圖5 水聽器接收聲源信號的示意圖

由聲源輻射的信號a[n]表示為

a[n]=ASL[n]cos(2πf0+θ)

(6)

式中：ASL[n]是聲源幅度;f0是聲源信號的頻率;θ是相位。

當信號在淺水環境中傳播時,與底部和表面接觸多次,從而導致多徑傳播現象[8],因此,聲信號傳播到水聽器會發生一些變化。單個水聽器上接收到的聲源信號x[n]表示為

x[n]=A(s[n])cos(2πfdn+φ)+W[n]

(7)

式中：W[n]是噪聲;fd是多普勒頻移[9];A(s[n])是接收信號幅度;s[n]是聲源位置的函數,稱為聲源狀態向量,它取決于聲源的水平距離、深度、徑向速度和幅度。

(7)式中的接收信號幅度A(s[n])可以分解為聲源幅度ASL[n],以及由傳播過程產生的幅度c[n],即不同水平距離-深度處信號的幅度變化,表示為

A(s[n])=ASL[n]c[n]

(8)

如果準確地知道環境參數se,那么就可精確得到由聲傳播引起的幅度衰減c[n]。但實際情況由于環境參數是變化的,因此,本文在統計模擬方法建模的基礎上,用概率密度估計方法計算出由環境參數引起的c[n]幅度值的變化,得到不同時刻(即不同水平距離-深度處)的幅度分布,由(7)、(8)式計算得出接收信號x[n]。

2.2 概率密度估計方法原理及仿真

采用密度估計方法對不同水平距離-深度處的幅度變化進行概率密度分布估計,從而實現接收信號的建模。

(9)

式中：ω(p-pi;h)是核函數;h為核寬度,核函數可以是任意的形式,這里定義是高斯核函數,表示為

(10)

(11)

進一步可簡化為

(12)

式中

(13)

(14)

聯立(13)式和(14)式可將(12)式簡化為

(15)

(15)式不是直接適用的,因為需要基礎概率密度函數f(p)的真實形式,本文中的f(p)為高斯核密度形式,因此,(15)式可以進一步簡化為

(16)

σ是基礎概率密度函數f(p)的標準差。

采用核密度估計的方法對聲源在水下傳播時聲壓隨距離衰減情況進行統計,得到接收信號幅度在不同水平距離-深度處的幅度衰減,其統計過程示意圖如圖6所示。圖6a)為統計模擬方法產生的聲場集合,圖6b)為對聲場集合中所有紅色方塊(表示水平距離-深度單元格)中聲壓幅度值統計的概率分布直方圖,是對聲場集合中1 071個傳播損失圖中每一個紅色方塊進行聲壓幅度的統計,其中紅色方塊的取值范圍是水平距離20 m,深度4 m。

圖6 某水平距離-深度位置處的聲壓幅度直方圖

圖7為對頻率109 Hz淺海聲源的聲場建模結果采用核密度估計的結果。由圖7可以看出核密度估計得出的概率密度函數與直方圖有很好的對應關系,對數據的擬合較好,因此選擇核密度方法實現聲源對聲源在不同位置處的幅度衰減估計,從而實現接收信號的仿真。

圖7 核密度估計結果示意圖

圖8為仿真接收信號與實測接收信號的時域與頻域波形對比圖,仿真接收信號時采用了與SwellEx-96海洋環境相同的參數,實際接收信號數據則是來自SwellEx-96水平線陣HLA South中的一個水聽器1在10～13 min時間段所接收的數據。

圖8 實測接收信號與仿真接收信號時域與頻域波形對比圖

由圖8可知,仿真的接收信號與實測接收信號之間除了存在幅度上的差異之外,時域波形起伏、頻譜分布等都很類似,因此表明仿真的接收信號是有效的,這里在原數據的基礎上進行了隔點采樣,采樣頻率降為5 000 Hz。

3 基于直方圖濾波的淺海聲源測距

3.1 直方圖濾波算法原理

直方圖濾波[11]是Bayes濾波[12]的一種基于網格的實現,將狀態空間分解成多個有限區域,劃分后利用直方圖濾波算法計算狀態向量,它給出狀態向量在每一個區域的概率,區域劃分越精細計算結果就越精確,本文將水平距離、深度劃分為20 m,4 m的區間,結合目標的先驗知識和狀態更新方程,利用直方圖進行逼近,實現對聲源的測距。

與大多數濾波器一樣,直方圖濾波也包含預測和更新兩步,下面將進行分析說明:

當聲源信號在水中傳播時,產生的壓力波動由一個水聽器接收,水聽器所接收的信號包含關于聲源位置和幅度的信息。

利用聲源運動參數和新的接收信號,可以在每個時間步長更新狀態向量,用后驗概率密度函數表示為f(sk[n]|x[1],…,x[n])。在前一時刻直方圖濾波的輸出表示為f(sk[n-1]|x[1],…,x[n-1]),使用直方圖濾波將其更新為當前時刻的輸出表示為

f(sk[n]|x[1],…,x[n])=ηf(x[n]|sk[n])×

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

(17)

式中，f(sj[n-1]|x[1,…,x[n-1]])表示直方圖濾波在n-1時刻產生的狀態向量,狀態向量s[n]由4個聲源參數組成,包括聲源到接收器的水平距離、深度、速度和幅度。

f(sk[n]|sj[n-1])稱為狀態更新方程。對于具有K個可能值的狀態向量,f(sk[n]|sj[n-1])具有K2個可能性,因為對于j和k都有j,k=1,2,…,K。因此,若狀態向量在時間n-1處的第j個值已知,則該狀態在時間n處的第k個值的概率由f(sk[n]|sj[n-1]得出。狀態更新方程是一個概率密度函數,表示關于狀態向量值如何隨時間變化的先驗知識,如(18)式所示

(18)

式中：m為均值向量;Σ為協方差矩陣,m表示為

m=[μr,μv,μz,μSL]T

(19)

式中,下標r,v,z,SL分別表示聲源的水平距離、速度、深度以及幅度。

因為狀態向量是獨立更新的,所以協方差矩陣Σ表示為

(20)

(17)式Σf(sk[n]|sj[n-1])f(sj[n-1]|x[1],…,x[n-1])是狀態向量的總和,使用卷積來計算。在更新過程中,狀態向量f(sj[n-1]|x[1],…,x[n-1])通過與狀態更新方程進行卷積,得到的是不同的概率值,表示在已知n-1時刻所有數據的情況下,在時刻每個狀態向量值分布的可能性。

(17)式中的f(x[n]|sk[n])是對聲學建模結果用核密度方法估計的結果。η是確保(21)式輸出的是有效的概率函數

(21)

總體來說,(17)式表示的是直方圖濾波一次狀態更新的結果f(sk[n]|x[n]),它描述了單個接收信號的狀態向量sk[n]在時間n時的可能性。此時,后驗概率密度f(sk[n]|x[n])比單個接收信號吸收了更多的信息,因為直方圖濾波還使用狀態向量先驗概率密度函數f(sj[n-1])、狀態更新方程f(sk[n]|sj[n-1])和描述聲場的概率密度函數f(x[n]|sk[n])。比起貝葉斯濾波器,直方圖實現起來容易得多,因為積分已被求和取代,非常適合在計算機上進行數值計算。

上述介紹了直方圖濾波的狀態更新過程,它的更新是遞歸的,因此前一次更新產生的后驗概率密度函數是后一次更新的先驗,直方圖濾波的算法描述如表2所示。

表2 直方圖濾波測距算法描述

直方圖濾波得出的是聲源狀態向量的后驗概率密度,包括水平距離,深度,幅度和速度,可以將其任何一個進行邊緣化[13],得到不同時間內對該參數的估計結果,邊緣化深度、幅度和速度可表示為

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

?f(r[n],v[n],z[n],ASL[n]|x[1],…,x[n])

dASL[n]dv[n]dz[n]

(22)

類似地,也可以對其他參數進行邊緣化。

總而言之,直方圖濾波算法即使先驗知識為非高斯、多模態、相關或非參數概率密度函數的形式,仍然可以實現。當了解狀態向量的先驗概率密度函數、狀態更新、環境參數和聲場導致狀態向量后驗的變化時,可以利用先驗知識實現精確的測距。

3.2 仿真實驗及結果分析

3.2.1 算法相關參數設置

仿真時,將聲源的水平距離、深度、徑向速度以及幅度作為狀態向量,結合聲源運動過程中接收到的聲壓值以及狀態更新方程,在每個劃分區域內根據前一時刻估計得到的聲源狀態向量,估計出當前時刻的聲源狀態信息。狀態更新的目的是準確表示聲源狀態向量參數隨時間變化的先驗知識。其中狀態向量的網格參數(狀態空間)設置如表3所示,不同的步長表示不同的劃分區間。

表3 狀態向量網格參數

表4給出了實驗過程中(18)式所示的狀態更新方程的均值和協方差的數值。

表4 狀態更新的參數

3.2.2 仿真結果

采用直方圖濾波算法進行目標測距,接收信號頻率為109Hz和127Hz,聲源位于海平面下9m,沿水平方向勻速運動,速度為2.5m/s,聲源運動時長為50min,水聽器位于水下100m。圖9為在不同時間段對109Hz聲源的水平距離和深度的估計結果,其中127Hz聲源實驗結果也同理可得。

圖9中每個圖對應不同時間下的距離估計結果,聲源的實際位置用紅色O表示,聲源的估計位置用紅色X表示,灰度為概率密度的對數,灰色顏色越深表示目標在該位置處的可能性越大,反之則越小。因剛開始時先驗知識較少,所以測距結果灰色分布范圍較大,隨著時間的增大,狀態向量在不斷更新,聲源在真實位置的可能性比在其他位置的可能性變大,灰色顏色變深且越來越集中。

圖9 淺海109 Hz聲源的距離估計結果

圖10為109Hz聲源在邊緣化其他3個參數(聲源幅度、速度、深度)后在不同時間內對水平距離和深度的估計結果,灰色較深的區域表示目標最可能的位置,紅線是實際聲源水平距離和深度隨時間的變化結果。圖11為對109Hz聲源的速度和幅度估計結果,紅線是實際聲源速度和幅度與時間的函數關系。圖10至11是對聲源狀態向量的后驗概率密度進行邊緣化生成的,通過將每個參數的概率密度函數表示為時間的函數,可以清楚地看到隨時間變化的參數估計結果。在剛開始時,狀態向量參數存在明顯的不確定性。這有2個原因:①先驗概率密度函數是四維均勻分布的,排除掉不可能的狀態向量參數需要一些信息,例如,一開始很難判斷一個聲源的幅度是接近一個較低的還是較高的幅度;②傳播模型在很短的時間范圍內往往是不準確的,隨著測距算法狀態向量參數的不斷更新,狀態向量參數的估計結果會提高。從圖11b)可以看出對聲源幅度估計的結果很準確,因此可以將聲源幅度邊緣化,從而觀察聲源速度對水平距離和深度估計結果的影響,繪制聲源速度、水平距離以及深度的三維等概率輪廓。圖12為對109Hz聲源在第9分鐘時狀態向量的三維等概率面?？梢钥闯?聲源在第9分鐘時,盡管聲源速度有些不確定,但此時水平距離和深度幾乎沒有不確定性,表明了速度對于水平距離的估計影響較小,證明了距離估計的有效性。

圖10 在不同時間下對聲源水平距離和深度估計結果 圖11 在不同時間下對聲源速度和幅度估計結果

圖12 聲源在9分時的三維等概率面

下面將對測距的誤差進行分析,通過選取不同的誤差窗口大小,從而計算對聲源距離估計結果的正確率。表5給出了在選取不同的誤差窗口下對距離估計結果的正確率,其中水平距離和深度所取的誤差分析窗口有所不同。

表5 不同誤差窗口距離估計的正確率

由表5可以看出,對水平距離估計正確率比較好的情況是窗口大小為±10m左右,正確率達到90%左右;而深度估計正確率比較好的情況則是在±0.5m左右,正確率也達到90%左右。誤差來源可能是對狀態空間的區間劃分以及狀態更新方程的參數設定,理論上狀態空間劃分越小距離結果越精確,但隨之會導致較大的計算量,本文在滿足應用性的前提性選擇了實驗所設定的劃分區間和參數。

圖13為109Hz和127Hz頻率聲源在不同時間的距離估計結果與實際距離的偏差,計算公式如(23)式所示。

圖13 距離估計結果誤差示意圖

(23)

圖13可以看出不同時刻距離估計的偏差,此外，127Hz聲源比109Hz聲源的距離估計偏差更大,這是因為頻率越高,聲源傳播的距離越短,因此在本文這種低頻遠距離的傳播環境時,高頻聲源信號相比低頻聲源信號會產生更大一點的偏差。

在文獻[14]中,采用匹配場方法在Swellex-96數據庫測試得出定位結果水平距離誤差除在個別時間內最大誤差在0.5km左右,其他的都相對比較準確,本文在實驗中估計的聲源軌跡與GPS數據也比較吻合,誤差也相對較小,因此本文對于淺海聲源的水下測距是有效的,具有一定的參考意義。

4 結 論

本文針對海洋環境參數是時變的,提出了采用統計模擬的方法來捕捉海洋環境的不確定性,從而實現聲場的建模,將環境不確定性映射為聲場的不確定性,而并非某單一環境下的聲場,使得建模結果更加接近真實的環境聲場。為了對這一結果進行統計分析,得到傳播損失在不同水平距離-深度處幅度的概率密度分布,采用了密度估計方法對該建模的結果進行估計,實現接收信號建模。最后提出了直方圖濾波測距算法,它利用了統計模擬方法聲學模型結果、接收信號數據以及聲源運動參數等先驗知識實現狀態更新,實現對聲源的距離估計,為實現水下聲信號的處理提供了一定的基礎。