基于空間平滑的多波束測深聲吶相干分布源方位估計

李海森 ,李 珊 ,周 天

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

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

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

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

1 分布源信號模型

1.1 點源模型

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

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

(1)

其中X(t)為陣列的M×1維快拍數據矢量，N(t)為陣列的M×1維噪聲數據矢量，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]提出了兩種分布式目標信號源模型：基于確定的角信號密度函數的相干分布式目標信號源和基于確定的角功率密度函數的非相干分布式目標信號源。多波束測深聲吶的回波信號是相干信號，可以用相干分布源模型描述。

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

(3)

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

相干信號的角信號密度函數可寫為：

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

(4)

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

(5)

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

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

(6)

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

(7)

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

(8)

陣列流型矢量為：

b(βi)=

(9)

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

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

2.1 廣義MUSIC算法

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

(10)

(11)

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

2.2 空間平滑算法

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

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

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

X=B(μ)S+Ns

(12)

Y=C(μ)S+Ny

(13)

且有，

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

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

(14)

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

(15)

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

(16)

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

(17)

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

矩陣表達形式為：

C(μ)≈B(μ)Φ

(18)

如圖1所示，以左側第一個子陣為參考子陣，第k個子陣接收的數據矩陣為：

Xk=BΦk-1S+N

(19)

其協方差矩陣為：

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

(20)

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

(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的每行至少應該有一個非零元素，且向量{Ψ1,…,ΨN}是線性無關的。由EEH=Rs/p可知，如果矩陣E有一行的元素全為零，那么這一行對應的信號能量就是零，顯然不符合條件。因此，第一個條件成立。由以上條件可以知道，向量{Ψ1,…,ΨN}是線性無關的。

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

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

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

3 計算機仿真

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

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

圖3 二維空間譜對比

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

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

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

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

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

4 多波束測深系統實驗數據分析

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

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

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

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

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

圖6 廣義MUSIC空間譜

5 結 論

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

參 考 文 獻

[1]李海森，陳寶偉，么 彬，等.多子陣高分辨海底地形探測算法及其FPGA和DSP陣列實現[J].儀器儀表學報,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]周 天.超寬覆蓋海底地形地貌高分辨探測技術研究[D].哈爾濱: 哈爾濱工程大學,2005.

[3]么 彬，李海森，周 天，等．多子陣超寬覆蓋海底地形探測方法試驗研究[J]．哈爾濱工程大學學報，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]馬培峰，嚴勝剛.相干分布源的方位估計[J].聲學技術，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]鄭 植，李廣軍.低復雜度相干分布源中心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]王永良，陳 輝，彭應寧，等.空間譜估計理論與算法[M].北京：清華大學出版社，2004：320.