李海森 ,李 珊 ,周 天

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

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

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

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

1 分布源信號(hào)模型

1.1 點(diǎn)源模型

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

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

(1)

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

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

(2)

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

1.2 分布源模型

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

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

(3)

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

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

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

(4)

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

(5)

相干分布源信號(hào)模型可以進(jìn)一步簡(jiǎn)化為：

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

(6)

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

(7)

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

(8)

陣列流型矢量為：

b(βi)=

(9)

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

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

2.1 廣義MUSIC算法

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

(10)

(11)

通過(guò)二維譜峰搜索，找出極大值點(diǎn)對(duì)應(yīng)的中心波達(dá)角度和角度擴(kuò)展參數(shù)。

2.2 空間平滑算法

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

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

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

X=B(μ)S+Ns

(12)

Y=C(μ)S+Ny

(13)

且有，

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

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

(14)

在θ=θ0i處，對(duì)a(θ)進(jìn)行泰勒級(jí)數(shù)展開(kāi)有：

(15)

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

(16)

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

(17)

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

矩陣表達(dá)形式為：

C(μ)≈B(μ)Φ

(18)

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

Xk=BΦk-1S+N

(19)

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

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

(20)

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

(21)

(22)

可以簡(jiǎn)化為：

(23)

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

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

(24)

其中，EEH=Rs/p。

(25)

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

(26)

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

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

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

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

3 計(jì)算機(jī)仿真

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

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

圖3 二維空間譜對(duì)比

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

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

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

圖5 方位估計(jì)誤差隨信噪比的變化

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

4 多波束測(cè)深系統(tǒng)實(shí)驗(yàn)數(shù)據(jù)分析

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

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

如圖6所示為兩個(gè)時(shí)刻的分布源廣義MUSIC算法空間譜。

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

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

圖6 廣義MUSIC空間譜

5 結(jié) 論

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

參 考 文 獻(xiàn)

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

[3]么 彬，李海森，周 天，等．多子陣超寬覆蓋海底地形探測(cè)方法試驗(yàn)研究[J]．哈爾濱工程大學(xué)學(xué)報(bào)，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ì)[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ì)方法[J].信號(hào)處理，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)寧，等.空間譜估計(jì)理論與算法[M].北京：清華大學(xué)出版社，2004：320.