基于RJMCMC方法的水下被動目標聲源數和方位聯合估計

陳 釗1, 劉鄭國2, 王海燕1, 申曉紅1, 和 斌1

?

基于RJMCMC方法的水下被動目標聲源數和方位聯合估計

陳 釗, 劉鄭國, 王海燕, 申曉紅, 和 斌

(1. 西北工業大學航海學院, 陜西西安, 710072; 2. 中國船舶重工集團公司, 北京, 100097)

針對傳統的子空間類寬帶信號處理方法需要較多的信號快拍數才能取得較好的估計效果這一問題, 本文將一種基于時域寬帶信號模型的貝葉斯高分辨估計方法用于水下被動目標的方位估計(DOA)。該方法使用寬帶信號的時域模型, 依據貝葉斯準則構建所需參數的后驗概率密度函數, 并采用可逆跳變馬爾科夫鏈蒙特卡羅(RJMCMC)方法執行貝葉斯計算, 對后驗概率密度函數進行峰值搜索。由于其能夠在不同維的參數空間中跳轉, 因此可以進行模型階數(聲源數)和目標方位的聯合估計。與傳統的子空間類方法相比, 這一方法能直接在時域進行信號處理且僅需顯著較少的信號快拍數。仿真結果表明, RJMCMC算法能較好地估計出聲源數和目標方位。

被動目標方位估計; 時域寬帶模型; 貝葉斯高分辨; 模型階數檢測; 可逆跳變馬爾可夫鏈蒙特卡羅(RJMCMC)

0 引言

精確定位遠場中相距很近的多個聲源是水聲信號處理中的關鍵問題之一。傳統的波束形成法受到瑞利限的制約, 分辨率很低。最大熵譜法和最小方差法是對常規波束形成法的簡單修正, 但兩者在估計性能上沒有大的突破。子空間類方法極大地突破了瑞利限的制約, 使高分辨目標方位估計方法向前邁進了一大步, 但是在低信噪比, 少快拍, 信號源相關以及存在系統誤差時, 此類方法的性能會受到很大的影響。最近十幾年, 貝葉斯高分辨方位估計方法受到普遍關注并被廣泛地研究。貝葉斯高分辨方位估計技術是利用信號和噪聲的聯合后驗概率密度函數對信號進行譜估計的一種方法。Lasenby和Fitzgerald最早在文獻[1]中提出了這種方法。它利用貝葉斯定理把待估計量看作是隨機變量, 引入被估計量的先驗概率, 因而改善了估計的性能。之后不同學者對貝葉斯高分辨方位估計方法進行了探索和改進, 提出了許多新的改進。

對于水下目標的檢測、方位估計和跟蹤問題, 與主動探測方式相比, 被動方式具有更好的隱蔽性。被動探測方式利用水下目標的輻射噪聲進行處理, 通常涉及寬帶信號處理問題。傳統的基于子空間類的高分辨方位估計方法常需要較多的快拍數且在頻域內進行處理。被廣泛使用的相干信號子空間方法(coherent signal subspace method, CSSM)需要方位預估, 其準確性嚴重影響到了目標方位估計的精度。

本文引入文獻[4]中提出的一種時域寬帶信號模型, 利用水下被動目標的輻射噪聲進行目標方位估計。首先得到基于貝葉斯模型的聲源數和目標參數的聯合后驗概率密度函數。貝葉斯高分辨方位估計方法需要對此概率密度函數進行多維搜索, 隨著待估計參數維數的增加, 計算量呈指數增長, 為此許多文獻提出使用可逆跳變馬爾可夫鏈蒙特卡羅(reversible jump Markov chain Monte Carlo, RJMCMC)方法進行后驗概率密度函數的峰值搜索。同時采用文獻[4]中的最大后驗(maximum aposterior, MAP)聲源恢復方法可以求得各個聲源的幅度。

本文將文獻[4]的時域寬帶信號模型和RJMCMC的方法結合起來用于水下被動目標的聲源數和方位的聯合估計。在較少的快拍情況下可以取得較好的估計精度, 且在估計出目標方位(direction of arrival, DOA)的同時可以進行聲源數檢測和聲源幅度恢復。仿真結果表明了此方法的有效性。

1 基于時域寬帶信號模型的貝葉斯高分辨DOA估計方法

本節簡述文獻[4]提出的時域寬帶信號模型。這一模型采用插值函數表示陣列接收的寬帶信號, 使得算法只需在時域內進行處理, 避免了傳統的聚焦類方法在頻域內處理的缺點, 從而無需對陣列信號進行傅里葉變換。

1.1 時域寬帶信號模型

假定有元均勻線列陣(uniform linear array, ULA),個遠場聲源s(t)分別從θ,=0,…,-1 入射, 聲源s(t)為有限帶寬的寬帶或窄帶信號, 其頻率響應和帶寬分別為

當聲源入射方位角≤π/2時, 陣元1首先接收到信號, (-1)τ時延之后陣元才接收到信號, 如圖1所示。聲源s()在第個陣元上產生的信號可以表示為

其中:τ為第個聲源信號的陣元間延遲(inter sensor delay, ISD), 定義如下

式中:為水中聲速;為均勻線列陣的陣元間距;(mτ)是對應于時延mτ的插值序列的第l個元素, 采用加權sinc函數作為插值序列, 其中(l-)為漢明窗,為歸一化時延值

(4)

(6)

式中,f是系統的采樣頻率, 滿足

陣元間距滿足式(8) (其中為可能的最大目標聲源頻率成分), 本文取=c/f

(8)

由式(3)和式(8)可知

(10)

寫成矩陣-向量形式為

因此

(12)

式中:()∈R^×是第l個插值矩陣;[,,…,τ]∈R為與聲源入射方位相對應的時延向量;∈R是次快拍時刻的聲源幅度向量;∈R是次快拍時刻服從(0,I)分布的高斯白噪聲。最終定義一個向量∈R

它可以看作是快拍與其基于插值逼近之間的誤差, 可寫作

(14)

對于入射方位角>π/2的情況, 陣元首先接收到入射信號, (-1)τ時延之后陣元1才接收到信號, 如圖2所示。

圖2 入射角φk>π/2時的陣列信號接收模型

式中:為一個交換矩陣, 它交換矩陣()的每一行

(16)

重新定義插值矩陣為

重寫時域寬帶信號模型

(18)

1.2 貝葉斯后驗概率密度函數

對于所有的快拍, 全似然函數為

(20)

式中,為真實聲源數的一個估計。由此通過貝葉斯理論, 聲源參數的后驗概率密度函數可以表示為全似然函數和聲源參數先驗函數的乘積

式中,是一個超參數, 它與信噪比有密切關系。聲源幅度先驗的選取同文獻[4],選為均勻先驗

(22)

的先驗選擇為期望值Λ為的泊松分布, 其中Λ是期望的聲源數

(24)

對于上述各參數, 人們感興趣的是聲源數和聲源時延。將式(22)~式(24)代入式(21), 并把“多余”參數和積分掉，可以得到

其中

(26)

式(25)即所需參數和的聯合后驗概率密度函數, 通常需要對后驗概率密度函數進行多維峰值搜索, 最大峰值所對應的時延和模型階數值便是入射聲源的DOA和數目。傳統峰值搜索方法計算量很大, 基于馬爾可夫鏈蒙特卡羅(MCMC)的方法可以大大地減少運算量。

2 基于RJMCMC算法的貝葉斯高分辨方位估計

研究的目標是在貝葉斯框架下進行聯合模型選擇(聲源數估計)和參數估計, 即通過最大化式(25)來求得聲源數和時延。通常采用的方法是RJMCMC算法。

傳統的RJMCMC方法由3種不同的運動組成: 更新運動(update move)、消失運動(death move)和產生運動(birth move)。

每次迭代執行上述3種運動的概率分別為u,d和b, 它們由式(27)確定, 其中為每次迭代模型階數變化的概率, 通常選取為=0.5,()為式(24)中服從期望為Λ的泊松分布

在RJMCMC算法的每次迭代中, 從與上述運動對應的3種不同的建議分布中抽取候選樣本。假定第次迭代的樣本為(,)

(28)

1) 更新運動: 此時聲源數固定不變, 僅更新時延向量值。本文選取以為均值, 以為方差的高斯分布作為建議分布, 從中抽取候選樣本~(,)。

2) 產生運動: 當<時有效。其中,為最大允許的聲源數(本文選取為6)。此時從分布[T,T]中隨機抽取一個樣本′, 并且加入上次迭代的樣本中成為候選樣本。

3) 消失運動: 當>1時有效。此時從上次迭代的樣本中隨機地去除一個時延值, 然后將新的時延組合作為候選樣本, 假設第個時延值被隨機地去除。

(30)

于是候選樣本以如下概率被接受

其中:和分別表示上次迭代值和本次候選值,(|)為建議函數。在本文的仿真中, 為了簡化, 候選樣本以下列概率被接受

(32)

這樣不僅簡化了運算, 并且在RJMCMC算法的迭代過程中樣本總是盡量向著更高的概率區域跳轉。迭代結束之后以最后一個樣本值作為最終的聲源數和方位估計的結果。于是一個簡化的RJMCMC算法的可概括為以下步驟。

1) 初始化, 從均勻分布[1,]中隨機第選取初始模型階數; 從與對應的均勻分布[T,T]中抽取初始時延值。

2) 第次迭代, 從均勻分布[0,1]中抽取隨機數, 并且利用式(27)計算各個運動被執行的概率u, d和b。若u<u則執行更新運動, 若≤≤d+u則執行消失運動, 否則執行產生運動。

3) 使用產生的候選樣本和上次迭代樣本值, 依據式(29)計算樣本的接受概率。若樣本被接受, 則; 否則。

4) 如此迭代下去直到迭代終止。

5) 將最后一次迭代值作為聲源數和方位估計的結果。

3 水下目標輻射噪聲模型

對于諸如艦船, 潛艇和魚雷等主要的水下水面目標, 它們都具有許多往復運動的機械部件, 這些部件所產生的振動輻射聲波到海水中便形成目標的輻射噪聲。目標輻射噪聲具有特殊重要性, 被動聲納就是利用其特性從自噪聲或環境噪聲背景上把它區分出來。按照噪聲源的不同它可以分為3類: 機械噪聲, 螺旋槳噪聲和水動力噪聲。其中前2類在多數情況下是主要的噪聲源, 艦船噪聲就是這3類噪聲中連續譜成分和線譜成分的迭加。本文采用文獻[9]的方法產生陣列接收寬帶信號來仿真水下被動目標的輻射噪聲。如圖3給出了簡化的3個水下目標輻射噪聲的功率譜模型。

圖3 水下目標輻射噪聲功率譜模型

4 計算機仿真

采用如下仿真條件。使用=10元ULA, 相應波束寬度為=10.09°。信號快拍數為=64。考慮稍復雜的3個聲源的情況, 都具有相等的信噪比(signal-to-noise ratio, SNR)。聲源信號的入射角分別為=6.00°,=12.00°和=18.00°, 兩兩間的間距接近半個波束寬度。寬帶聲源的頻譜范圍為∈[900, 2 900] Hz,∈[1 300, 3 400] Hz和∈[1 100, 3 100]Hz。系統采樣頻率f=8 000 Hz。水下聲速為=1 500 m/s, 陣元間距為f= 0.187 5 m。仿真參數詳見表1。

表1 仿真采用的參數值

采用表1的參數進行仿真, 圖4、圖5和圖6給出了單次運行使用RJMCMC算法時貝葉斯高分辨水下被動目標方位估計方法(wide band Bayesian estimation method,WBBEM)的性能。為了程序運行效率, 本文中RJMCMC算法僅采用200次迭代。

圖4 單次運行時算法的聲源數估計結果

圖5 單次運行時算法的方位估計結果

圖6 單次運行時算法的聲源幅度恢復結果

圖4給出了算法聲源數估計的收斂結果。由圖可知, 算法收斂相當迅速, 從第55次迭代算法便正確地估計出了入射聲源數。并且由多次獨立仿真可以看出, 大多數時, 算法在30次迭代之內便能正確地估計出聲源數為3。

圖5給出了單次運行時算法的方位估計結果。由圖可以看出, 在迭代算法結束時, 3個入射聲源方位都能被較為精確地估計出來。

圖6給出了采用MAP算法的聲源幅度恢復結果, 由圖可知, 3個聲源的幅度都被較好地恢復出來, 這能為后續的基于頻譜估計的目標識別打下基礎。

圖7~圖9給出了采用表1中的參數進行100次獨立蒙特卡羅(MC)仿真時RJMCMC算法的性能。圖7給出了進行100次獨立MC仿真時聲源數估計值的分布。由圖可以看出, 在100次獨立仿真中, 有83次算法正確地估計出入射聲源數為3階, 其余17次傾向于低估聲源數, 且估計結果均為2。

圖7 進行100次獨立蒙特卡洛仿真的聲源數估計結果

圖8 聲源數估計值為2時的方位估計結果

圖9 聲源數估計值為3時的方位估計結果

圖8給出了聲源數估計值為2時的DOA估計結果。由圖可以看出, 在這17次獨立蒙特卡羅仿真中, 算法總是大致地估計出了聲源1和3而遺漏了聲源2, 且聲源1的角度估計值傾向于偏大而聲源3的角度估計值傾向于偏小。由此可以看出具有相同信噪比的3個聲源中處于中間位置的聲源的估計受到了其兩側聲源的干擾, 同時兩側聲源的DOA估計結果同樣受到了中間聲源的影響而趨于向其靠近。圖中三角形的方位估計值中心相對于六角星的真實方位值明顯地向右下側偏移清晰地表現了這一趨勢。

圖9給出了正確地估計出聲源數3時的DOA估計結果。為了表述方便, 分別繪出聲源1和2, 聲源2和3的DOA估計結果。由圖可以看出單次運行的DOA估計值(星形)絕大多數都集中地分布在真實的DOA值(圖中用六角星表示)附近, 其中心位置(圖中用三角形表示)稍有偏離。聲源1和2的估計值中心約為[6.1°,11.9°]; 而聲源2和3的估計值中心約為[11.9°,18.2°]。由此可見算法很好地估計出了3個聲源的方位, 具有較好的性能。

表2給出了不同信噪比下, 采用RJMCMC算法的聲源數估計結果和方位估計的均方根誤差。本文出于程序運行的簡便采用的迭代次數偏少(僅200次), 因此對于聲源數的估計無法達到很理想的效果; 但是一旦算法采用了較多的迭代次數(例如1 000次以上), 則聲源數的估計性能將會得到改善。

表2 不同信噪比下聲源數估計結果及方位估計均方根誤差

5 結束語

最近10年貝葉斯高分辨技術越來越廣泛地應用到目標方位估計中, 在低信噪比, 少快拍和信源間隔很近的情況下取得了比子空間類方法高得多的方位估計精度。

本文將文獻[4]中提出的基于時域寬帶模型的貝葉斯高分辨方位估計方法引入水下被動目標方位估計當中, 在水下目標輻射噪聲的寬帶表征中采用了插值矩陣表示方法。與傳統的處理寬帶信號的聚焦類子空間方法相比, 本文方法具有顯著的優勢。首先, 它可以直接在時域內進行信號處理; 其次, 它在信號快拍數較少的情況下也能取得較好的估計結果, 而相比之下聚焦類寬帶方法在少快拍情況下性能會急劇惡化; 最后, 它可以在進行方位估計的同時進行模型階數(聲源數)的檢測和聲源幅度的恢復, 而沒有顯著地增加計算量, 其中信號幅度的恢復有助于基于譜估計的目標識別。

本文采用RJMCMC算法對求得的目標聲源數與目標參數的聯合后驗概率密度函數進行峰值搜索。并且選取了較復雜的3個目標輻射聲源的情況進行計算機仿真。仿真結果顯示WBBEM算法不僅能夠較準確地估計出目標輻射聲源數, 同時可以較為精確地估計出目標方位, 并且聲源幅度也能被很好地恢復出來, 取得了較好的估計效果。

[1] Lasenby J, Fitzgerald W J. A Bayesian Approach to High Resolution Beamforming[J]. Radar and Signal Processing, IEEE Proceedings F, 1991, 138(6): 539- 544.

[2] Huang J G, Sun Y. A New Gibbs Sampling DOA Estimator Based on Bayesian Method[J]. Acoustic, Speech and Processing, 2003(5): 229-232.

[3] Sun Y, Liu K W, Huang J G, et al. Bayesian DOA Estimator by Gibbs Sampling[J]. Computer Engineer- ing and Applications, 2002, 38(12): 36-38.

[4] William N, Reilly J P, Kirubarajan T, et al. Wideband Array Signal Processing Using MCMC Methods[J]. IEEE Transactions on Signal Processing, 2005, 53(2): 411-426.

[5] Andrieu C, Doucet A. Joint Bayesian Model Sele- ction and Estimation of Noisy Sinusoids Via Rever- sible Jump MCMC[J]. IEEE Transactions on Signal Processing, 1999, 47(10): 2667-2676.

[6] Green P J. Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determina- tion[J]. Biometrika, 1995, 82(4): 711-732.

[7] Larocque J R, Reilly J P. Reversible Jump MCMC for Joint Detection and Estimation of Sources in Colored Noise[J]. IEEE Transactions on Signal Proce- ssing, 2002, 50(2): 231-240.

[8] 于紅旗, 劉劍, 黃知濤, 等. 陣列接收寬帶信號的建模方法及仿真[J] . 電子對抗, 2006(4): 7-12.

[9] William N, Reilly J P, Kirubarajan T. Wideband Array Signal Processing Using Sequential Monte Carlo Methods[R]. Department of Electrical and Co- mputer Engineering, McMaster University, 2003: 1-28.

(責任編輯: 楊力軍)

Joint Estimation of Source Number and DOA for Underwater Passive Target Based on RJMCMC Method

CHEN Zhao, LIU Zheng-guo, WANG Hai-yan, SHEN Xiao-hong, HE Bin

(1.College of Marine Engineering, Northwestern Polytechnical University, Xi′an 710072, China; 2.China Shipbuilding Industry Corporation, Beijing 100097, China)

Conventional wideband signal processing method of subspace class wideband requires more snapshots to get a better estimation result. We apply a Bayesian high resolution direction of arrival(DOA) estimation method based on the time-domain wideband signal model to underwater passive target DOA estimation. Compared with conventional subspace class methods, the proposed method can process signal directly in time-domain and needs remarkably less snapshots. We use a time-domain wideband signal model to construct a posterior probability density function(PDF) of the parameters according to Bayes law, and then use the reversible jump Markov chain Monte Carlo(RJMCMC) method for Bayesian estimation to search peaks of the posterior PDF. This method can jump among parameter spaces with different orders, thus can implement joint model order (i.e. source number) and DOA estimation. Simulation results show that the RJMCMC method exhibits good performance in joint estimation of model order and DOA.

passive target direction of arrival (DOA) estimation; time-domain wideband model; Bayesian high resolution; model order detection; reversible jump Markov chain Monte Carlo (RJMCMC)

TJ630.34; TN911.7

A

1673-1948(2011)06-0415-08

2011-03-17;

2011-05-19.

國家自然科學基金(60972153), 博士點基金(20096102110038), 西北工業大學基礎研究基金(NPU-FFR-JC201004).

陳 釗(1983-), 男, 在讀博士, 研究方向為水下目標的識別、方位估計與跟蹤.