李海森 ,李 珊 ,周 天

(1.哈爾濱工程大學(xué) 水聲工程學(xué)院，哈爾濱 150001;2.哈爾濱工程大學(xué) 水聲技術(shù)重點實驗室，哈爾濱 150001)

多波束測深聲吶是水下地形地貌測量最主要的儀器之一，在軍事和民用方面都有廣泛的需求。在測量時，先向水下發(fā)射窄帶脈沖聲信號，并接收海底及水體中散射體的反向散射信號，通過對回波到達角度和回波到達時間進行估計，獲得海底的深度信息。

目前多波束測深系統(tǒng)采用的方位估計方法都是建立在點信源模型之上的[1-3]，在水下環(huán)境及海底地形較復(fù)雜的情況下，海底反向散射信號往往具有空間角度擴散特性，基于點源假設(shè)的方位估計算法性能嚴(yán)重惡化，降低了多波束測深系統(tǒng)的精度[4-6]。很多學(xué)者針對分布源進行了研究，但現(xiàn)有方法不能有效的區(qū)分兩個相干的分布源，因而常假設(shè)兩個分布源之間不相干。文獻[7]提出采用Toeplitz方法解決這一問題并進行了驗證，但該方法無法估計角度擴展參數(shù)并且精度較低。

針對以上問題，為解決多波束測深聲吶相干分布源的方位估計問題，提出采用空間平滑的分布源廣義MUSIC算法，經(jīng)過公式推導(dǎo)證明了該方法具有理論依據(jù)，通過仿真驗證了算法解相干的有效性，分析其方位估計精度以及不同信噪比條件下的性能，并采用多波束測深系統(tǒng)的實驗數(shù)據(jù)驗證了算法的有效性。

1 分布源信號模型

1.1 點源模型

N個遠(yuǎn)場窄帶信號入射到由M個陣元組成的空間均勻直線陣列上，在窄帶遠(yuǎn)場信號源的假設(shè)下有：

X(t)=AS(t)+N(t)

(1)

其中X(t)為陣列的M×1維快拍數(shù)據(jù)矢量，N(t)為陣列的M×1維噪聲數(shù)據(jù)矢量，S(t)為空間信號的N×1維矢量，A為空間陣列的M×N維流型矩陣，且有：

A=[a(β1)a(β2) …a(βN)]

(2)

其中a(β)=[1 exp(-jβ) … exp(-j(M-1)β)]T，d為陣元間距，λ為波長，β=2πdsinθ/λ。

1.2 分布源模型

Valaee等[8]提出了兩種分布式目標(biāo)信號源模型：基于確定的角信號密度函數(shù)的相干分布式目標(biāo)信號源和基于確定的角功率密度函數(shù)的非相干分布式目標(biāo)信號源。多波束測深聲吶的回波信號是相干信號，可以用相干分布源模型描述。

在加性噪聲背景下，N個窄帶分布源信號到達接收陣列，陣元間距為半波長，接收的數(shù)據(jù)矢量可以表示為：

(3)

其中，βi和si(β-βi,t)分別為第i個分布源的中心波達角度和t時刻的角信號密度函數(shù)，式中的積分限根據(jù)β=2πdsinθ/λ=πsinθ，取-π≤β≤π。

相干信號的角信號密度函數(shù)可寫為：

si(β-βi,t)=si(t)gi(β-βi)

(4)

其中:gi(β-βi)是一個以βi為對稱中心的確定性函數(shù)，且滿足：

(5)

相干分布源信號模型可以進一步簡化為：

X(t)=BS(t)+N(t)

(6)

其中B=[b(β1)b(β2) …b(βn)]，且有：

(7)

當(dāng)角信號分布函數(shù)符合確定分布時，由式(7)可得到陣列流型矢量的閉式解，例如以中心波達角度βi為中心，在角度擴展Δi范圍內(nèi)符合均勻分布時，角信號密度函數(shù)為：

(8)

陣列流型矢量為：

b(βi)=

(9)

2 基于空間平滑的廣義MUSIC算法

研究人員以點源MUSIC算法為基礎(chǔ)，將其推廣到分布源參數(shù)估計中，即廣義MUSIC方法[8-9]，該方法中的相干分布源指同一分布源的各分量之間是相干的，而估計多個相干分布源時，假定不同分布源之間是不相干的。由于多波束測深系統(tǒng)接收到的不同方位的信號之間是相干的，因此，必須進行解相干才能獲得正確的測量結(jié)果。文獻[7]提出采用Toeplitz方法進行分布源解相干，但通過該方法獲得的二維空間譜僅能估計中心波達角度，無法對角度擴展參數(shù)進行估計，且精度較差。本文通過推導(dǎo)證明空間平滑方法能夠有效的對分布源信號進行解相干，獲得的二維空間譜可以估計分布源中心波達角度和角度擴展參數(shù)，從而提出了基于空間平滑解相干處理的廣義MUSIC算法。

2.1 廣義MUSIC算法

考慮加性噪聲和N個窄帶相干分布源的情況，根據(jù)式(6)的相干分布源模型，可得數(shù)據(jù)的協(xié)方差矩陣為：

(10)

(11)

通過二維譜峰搜索，找出極大值點對應(yīng)的中心波達角度和角度擴展參數(shù)。

2.2 空間平滑算法

空間平滑算法是針對點信源相干時，陣列接收的數(shù)據(jù)協(xié)方差矩陣的秩降低，信號子空間維數(shù)小于信源數(shù)而提出的[10]。對于均勻線陣，M為陣元數(shù)，N為信號源數(shù)?？臻g前向平滑技術(shù)原理如圖1所示，假設(shè)將M元直線陣分為相互交錯的p個子陣，每個子陣陣元數(shù)為m，則有M=p+m-1。

圖1 前向空間平滑算法原理

當(dāng)接收信號為相干分布源且角度擴展較小時，對相鄰兩子陣廣義方向矢量做泰勒近似可以獲得關(guān)于中心DOA的旋轉(zhuǎn)矩陣[11]。相鄰兩個子陣接收到的數(shù)據(jù)可用向量表示為：

X=B(μ)S+Ns

(12)

Y=C(μ)S+Ny

(13)

且有，

其中τ(θ)為兩子陣相同位置陣元上信號的傳輸時延。

相干分布式信源中心波達方向定義為：

(14)

在θ=θ0i處，對a(θ)進行泰勒級數(shù)展開有：

(15)

由于gi(θ;μi)對稱，則有：

(16)

設(shè)f(θ)=2π(d/λ)sinθ，當(dāng)d?λ時，可忽略f′(θ)=2π(d/λ)cosθ，同理可得：

(17)

顯然c(μi)≈b(μi)exp(j2πdsinθ0i/λ)。

矩陣表達形式為：

C(μ)≈B(μ)Φ

(18)

如圖1所示，以左側(cè)第一個子陣為參考子陣，第k個子陣接收的數(shù)據(jù)矩陣為：

Xk=BΦk-1S+N

(19)

其協(xié)方差矩陣為：

Rk=BΦk-1Rs(Φk-1)HBH+σ2I

(20)

對各子陣協(xié)方差矩陣求均值得到修正協(xié)方差矩陣：

(21)

(22)

可以簡化為：

(23)

其中，G是N×pN維矩陣：

G=[EΦE…Φp-1E]

(24)

其中，EEH=Rs/p。

(25)

其中，eij是矩陣E的第i行第j列的元素，且有：

(26)

由式(25)可見，要證明矩陣G的秩是N，即行滿秩，矩陣E的每行至少應(yīng)該有一個非零元素，且向量{Ψ1,…,ΨN}是線性無關(guān)的。由EEH=Rs/p可知，如果矩陣E有一行的元素全為零，那么這一行對應(yīng)的信號能量就是零，顯然不符合條件。因此，第一個條件成立。由以上條件可以知道，向量{Ψ1,…,ΨN}是線性無關(guān)的。

采用空間平滑解相干的代價是減小基陣的有效孔徑。同理，如果采用圖2所示的后向空間平滑算法，也能得到與上面類似的結(jié)論。

2.1 基本情況 396例患者中不符合入組60例，失訪38例，最終例符合條件298例。298例中，男215例，女83例；年齡（59.6±15.6）歲；行急診住院為177例，留院觀察121例；急性肺栓塞2例，主動脈夾層13例，重癥心肌炎1例；30 d死亡10例。

圖2 后向空間平滑算法原理

3 計算機仿真

考慮由16個陣元組成的均勻線陣，陣元間距為半波長，在加性白噪聲條件下，兩個相干分布源的角信號分布函數(shù)采用均勻分布，波達方向分別是θ1=10°,θ2=-5°，角度擴展分別為μ1=2°,μ2=4°。信噪比為10 dB，快拍數(shù)50。采用各種解相干方法得到的廣義MUSIC二維空間譜如圖3所示，圖3(a)不進行解相干，圖3(b)采用文獻[7]中的Toeplitz方法，圖3(c)采用前向空間平滑進行解相干，取子陣數(shù)為3，子陣陣元數(shù)為14。

當(dāng)不考慮角度擴展參數(shù)的估計時，可將二維譜投影到波達角度方向上進行一維搜索，如圖4所示。

圖3 二維空間譜對比

圖4 空間譜在波達角度方向上的投影

通過以上仿真可以得出以下結(jié)論：① 如圖3(a)、圖4(a)所示，當(dāng)分布源相干時，必須進行解相干，否則將不能得到正確的結(jié)果；② 如圖3(b)、圖4(b)所示，Toeplitz方法能夠分辨兩個相干分布源，得到中心波達角度結(jié)果θ1=10.4°,θ2=-6.6°，但是角度擴展信息完全丟失；③ 如圖3(c)、圖4(c)所示，空間平滑算法可以對分布源進行解相干，并且能夠同時估計中心波達角度和角度擴展，結(jié)果為θ1=10°,θ2=-5.4°，μ1=2.2°,μ2=4.3°；④ 對比圖4(b)和圖4(c)的結(jié)果，采用空間平滑估計中心波達角度比Toeplitz方法誤差更小，結(jié)果更加精確。

由于多波束測深聲吶邊緣波束的反向散射強度較弱，回波信噪比較低，因此，有必要對算法在低信噪比下的性能進行分析。仿真條件與上面相同，在每個信噪比下進行50次獨立實驗，求得分布源方位估計標(biāo)準(zhǔn)差，如圖5所示，圖5(a)和圖5(b)分別為兩個目標(biāo)的估計標(biāo)準(zhǔn)差。

圖5 方位估計誤差隨信噪比的變化

由圖5可見：① 算法的方位估計誤差隨信噪比的變化不大，因此，當(dāng)邊緣波束回波的信噪比較低時，依然可以準(zhǔn)確的估計分布源的方位；② 角度擴散較小的目標(biāo)，其方位估計的精度較高。

4 多波束測深系統(tǒng)實驗數(shù)據(jù)分析

為驗證算法的有效性，本文對自主研制的國產(chǎn)高分辨多波束測深聲吶的湖試數(shù)據(jù)進行處理。該試驗于2011年1月，在湖北省長陽縣清江水庫進行。多波束測深系統(tǒng)頻率為300 kHz，采樣頻率為40 kHz，接收基陣采用陣元數(shù)為40的均勻直線陣，陣元間距半波長。

采用基于空間平滑的廣義MUSIC算法對多波束測深系統(tǒng)原始數(shù)據(jù)進行處理，有以下幾點說明：① 由于在未知分布源角信號密度分布情況的情況下，采用均勻分布往往能得到更為準(zhǔn)確的估計結(jié)果[12]，因此假設(shè)分布源的角密度函數(shù)符合均勻分布；② 算法采用前向空間平滑，將陣元數(shù)為40的均勻直線陣分為9個子陣，每個子陣陣元數(shù)為32，對各子陣的數(shù)據(jù)協(xié)方差矩陣求平均，得到修正的數(shù)據(jù)協(xié)方差矩陣；③ 參考點源信源數(shù)估計，在平坦區(qū)域信源數(shù)一般為2，為防止遺漏信號，在算法中取信源數(shù)為3；④ 譜峰搜索范圍取中心波達角度范圍-90°至90°，取角度擴展范圍為0°至6°。

如圖6所示為兩個時刻的分布源廣義MUSIC算法空間譜。

在圖6(a)中，從空間譜上可以清晰的分辨兩個方向的信號，譜峰所在位置處的角度擴展很小，此時信源符合點源模型，點源MUSIC算法能夠較好的逼近分布源廣義MUSIC方法，獲得良好的DOA估計結(jié)果；圖6(b)為另一時刻的廣義MUSIC空間譜，根據(jù)譜峰搜索結(jié)果，兩個目標(biāo)的角度擴展分別為2°和3.8°，此時，多波束測深系統(tǒng)接收的海底反向散射信號具備分布源特性。

如圖7所示，為算法對整ping多波束測深數(shù)據(jù)進行處理得到DOA-TOA曲線。從圖7可以看出：① 空間平滑可以對相干分布源進行解相干是有效的；② 算法對點源和分布源都有效；③ 該方法在小角度范圍內(nèi)的測量結(jié)果較為發(fā)散，外側(cè)的性能遠(yuǎn)優(yōu)于內(nèi)側(cè)，說明該方法適應(yīng)性尚有一定的局限性。主要原因是：第一，小角度范圍內(nèi)回波持續(xù)樣本點數(shù)較少(正下方水深約56 m，在-10°至10°內(nèi)約有50個采樣點)，算法選取快拍數(shù)為20(兩倍脈寬，脈寬為0.1 ms)，小快拍數(shù)使算法的性能受到影響；第二，DOA隨時間的變化率較大，因此測量精度降低。

圖6 廣義MUSIC空間譜

5 結(jié) 論

針對復(fù)雜水聲環(huán)境中點源方位估計算法性能較差的問題，本文提出基于空間平滑的廣義MUSIC算法對多波束測深聲吶相干分布源的方位進行估計。公式推導(dǎo)證明了算法具有理論依據(jù)；計算機仿真證明了算法能夠獲得比Toeplitz方法更高的方位估計精度，且能對角度擴展進行估計，分析了算法在不同信噪比下的方位估計性能以保證邊緣波束信噪比較低時算法的穩(wěn)健性；最后采用多波束測深系統(tǒng)的實驗數(shù)據(jù)對算法進行了驗證。

參 考 文 獻

[1]李海森，陳寶偉，么 彬，等.多子陣高分辨海底地形探測算法及其FPGA和DSP陣列實現(xiàn)[J].儀器儀表學(xué)報,2010,31(2):281-286.

LI Hai-sen,CHEN Bao-wei,YAO Bin,et al.Implementation of high resolution sea bottom terrain detection method based on FPGA and DSP array[J].Chinese Journal of Scientific Instrument.2010,31(2):281-286.

[2]周 天.超寬覆蓋海底地形地貌高分辨探測技術(shù)研究[D].哈爾濱: 哈爾濱工程大學(xué),2005.

[3]么 彬，李海森，周 天，等．多子陣超寬覆蓋海底地形探測方法試驗研究[J]．哈爾濱工程大學(xué)學(xué)報，2008，29(10)：1076-1081．

YAO Bin，LI Hai-sen，ZHOU Tian,et al．Test study of multiple sub-array ultra-wide coverage seabed terrain detection method[J]．Journal of Harbin Engineering University,2008.29(10):1076-1081.

[4]Astely D,Ottersten B.The effect of local scattering on direction of arrival estimation with MUSIC [J].IEEE Trans.Signal Processing,1999,47 (12): 3220-3234.

[5]Tabrikian J,Messer H.Robust localization of scattered sources[J].Proceedings of the Tenth IEEE Workshop on Statistical Signal and Array Processing,Pocono Manor,Paname: IEEE,2000,453-457.

[6]Raich R,Goldberg J and Messer H.Bearing estimation for a distributed source via the conventional beamformer[J].Proceedings of Ninth IEEE SP Workshop on Statistical Signal and Array Processing,Portland OR,USA: IEEE,1998,5-8.

[7]馬培峰，嚴(yán)勝剛.相干分布源的方位估計[J].聲學(xué)技術(shù)，2007,26(5):817-821.

MA Pei-feng,YAN Sheng-gang.Estimation of the DOA of coherent distributed sources[J].Technical Acoustics,2007,26(5):817-821.

[8]Valaee S,Champagne B,Kabal P.Parametric localization of distributed sources[J].IEEE Trans.SP,1995,43(9): 2144-2153.

[9]Lee Y U,Choi J,Song I,et al.Distributed source modeling and direction-of-arrival estimation techniques[J].IEEE Trans.SP,1997,45(4):960-969.

[10]Shan T J,Wax M,Kailath T.On spatial smoothing for estimation of coherent signals[J].IEEE Trans.On ASSP,1985,33(4): 806-811.

[11]鄭 植，李廣軍.低復(fù)雜度相干分布源中心DOA估計方法[J].信號處理，2010,26(10):1516-1520.

ZHENG Zhi,LI Guang-jun.Low-complexity method for the central DOA estimation of coherently distributed source [J].Signal Processing,2010,26(10):1516-1520.

[12]王永良，陳 輝，彭應(yīng)寧，等.空間譜估計理論與算法[M].北京：清華大學(xué)出版社，2004：320.