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

?

基于RJMCMC方法的水下被動目標聲源數(shù)和方位聯(lián)合估計

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

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

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

被動目標方位估計; 時域?qū)拵Ｐ? 貝葉斯高分辨; 模型階數(shù)檢測; 可逆跳變馬爾可夫鏈蒙特卡羅(RJMCMC)

0 引言

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

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

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

本文將文獻[4]的時域?qū)拵盘柲Ｐ秃蚏JMCMC的方法結(jié)合起來用于水下被動目標的聲源數(shù)和方位的聯(lián)合估計。在較少的快拍情況下可以取得較好的估計精度, 且在估計出目標方位(direction of arrival, DOA)的同時可以進行聲源數(shù)檢測和聲源幅度恢復。仿真結(jié)果表明了此方法的有效性。

1 基于時域?qū)拵盘柲Ｐ偷呢惾~斯高分辨DOA估計方法

本節(jié)簡述文獻[4]提出的時域?qū)拵盘柲Ｐ汀＿@一模型采用插值函數(shù)表示陣列接收的寬帶信號, 使得算法只需在時域內(nèi)進行處理, 避免了傳統(tǒng)的聚焦類方法在頻域內(nèi)處理的缺點, 從而無需對陣列信號進行傅里葉變換。

1.1 時域?qū)拵盘柲Ｐ?/p>

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

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

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

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

(4)

(6)

式中,f是系統(tǒng)的采樣頻率, 滿足

陣元間距滿足式(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)

重新定義插值矩陣為

重寫時域?qū)拵盘柲Ｐ?/p>

(18)

1.2 貝葉斯后驗概率密度函數(shù)

對于所有的快拍, 全似然函數(shù)為

(20)

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

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

(22)

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

(24)

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

其中

(26)

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

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

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

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

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

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

(28)

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

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

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

(30)

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

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

(32)

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

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

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

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

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

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

3 水下目標輻射噪聲模型

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

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

4 計算機仿真

采用如下仿真條件。使用=10元ULA, 相應波束寬度為=10.09°。信號快拍數(shù)為=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。系統(tǒng)采樣頻率f=8 000 Hz。水下聲速為=1 500 m/s, 陣元間距為f= 0.187 5 m。仿真參數(shù)詳見表1。

表1 仿真采用的參數(shù)值

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

圖4 單次運行時算法的聲源數(shù)估計結(jié)果

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

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

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

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

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

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

圖7 進行100次獨立蒙特卡洛仿真的聲源數(shù)估計結(jié)果

圖8 聲源數(shù)估計值為2時的方位估計結(jié)果

圖9 聲源數(shù)估計值為3時的方位估計結(jié)果

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

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

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

表2 不同信噪比下聲源數(shù)估計結(jié)果及方位估計均方根誤差

5 結(jié)束語

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

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

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

[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), 西北工業(yè)大學基礎研究基金(NPU-FFR-JC201004).

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