聲矢量陣快速子空間方位估計(jì)算法

梁國(guó)龍，張 柯，安少軍，范 展

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

聲矢量傳感器由聲壓傳感器與軸向正交的振速傳感器組成，可同時(shí)、共點(diǎn)測(cè)量聲場(chǎng)的聲壓和質(zhì)點(diǎn)振速信息，與傳統(tǒng)的聲壓傳感器相比，其獲得的信息量大為增加.近十年來(lái)，許多文獻(xiàn)對(duì)聲矢量陣的高分辨DOA算法進(jìn)行了廣泛研究［1-5］.在文獻(xiàn)［1］中，基于Nehorai處理框架的聲矢量陣信號(hào)處理方法，僅僅把振速信息看成與聲壓相同的獨(dú)立信息來(lái)處理，并沒(méi)有利用聲矢量陣中聲壓和振速的相干性，即抗各向同性噪聲的能力.基于聲壓振速聯(lián)合處理的抗干擾能力，文獻(xiàn)［2-5］提出了一系列基于傳統(tǒng)子空間算法(諸如MUSIC、ESPRIT等)的聲壓與振速聯(lián)合處理的聲矢量陣高分辨DOA估計(jì)方法，與Nehorai框架的算法相比，其在多目標(biāo)方位估計(jì)精度、分辨信噪比門限、分辨成功概率等方面都具有更好的性能，但是，以上算法在提高聲矢量陣DOA估計(jì)性能的同時(shí)，其繁重的計(jì)算量卻沒(méi)有得到改善，這使得工程應(yīng)用受到限制.

在眾多的子空間快速估計(jì)算法中，多級(jí)維納濾波器(MSWF)［6-12］因無(wú)需估計(jì)協(xié)方差矩陣從而使其可以應(yīng)用在小樣本支撐的信號(hào)環(huán)境中，而且收斂速度快，能夠?qū)r(shí)變信號(hào)進(jìn)行快速跟蹤，更重要的是其無(wú)需特征值分解運(yùn)算，大大降低了運(yùn)算量.文獻(xiàn)［6］將多級(jí)維納濾波器應(yīng)用到標(biāo)量陣的子空間分解中，指出若取多級(jí)維納濾波器的期望信號(hào)為參考陣元的接收信號(hào)，則經(jīng)前向遞推得到的多級(jí)維納濾波器的匹配濾波器可作為信號(hào)子空間基的估計(jì)值.文獻(xiàn)［7］將參考陣元的輸出延時(shí)后作為多級(jí)維納濾波器的期望信號(hào)，在空時(shí)白噪聲條件下，提高了文獻(xiàn)［6］算法的DOA估計(jì)性能，但是其期望信號(hào)選擇方式影響了MSWF算法的實(shí)時(shí)性，且當(dāng)噪聲具有時(shí)間相關(guān)性時(shí)，該算法性能將與文獻(xiàn)［6］算法相同.文獻(xiàn)［8-9］在硬件平臺(tái)上實(shí)現(xiàn)了MSWF算法對(duì)信號(hào)子空間的快速估計(jì)，并取得了良好的效果.文獻(xiàn)［10］將多級(jí)維納濾波器應(yīng)用到MMUSIC(Modified MUSIC)算法中，實(shí)現(xiàn)了基于Nehorai框架的聲矢量陣相干源的快速DOA估計(jì).基于聲壓振速聯(lián)合處理的抗各向同性干擾特性，本文提出一種新的期望信號(hào)選擇方法，不影響MSWF算法的實(shí)時(shí)性，且使MSWF算法在保持高精度DOA估計(jì)的同時(shí)，大大減小了計(jì)算量.

1 矢量陣聲壓振速組合MUSIC算法

1.1 聲矢量陣輸出模型

考慮二維平面情況下，M個(gè)矢量傳感器等間距排列成一聲矢量陣線陣，假設(shè)均為各向同性陣元，放置于各向同性的均勻流體介質(zhì)中，陣列遠(yuǎn)場(chǎng)中在以線陣軸線的法線為參考的θk(k=1，2，…，K)處有K個(gè)波長(zhǎng)為λ的窄帶平面波入射，則聲矢量陣的輸出模型為

式中:Yp(t)為聲壓傳感器的輸出矢量;Yvx(t)為x方向振速傳感器的輸出矢量;Yvy(t)為y方向振速傳感器的輸出矢量;A(θ)=［a(θ1)，a(θ2)，…，a(θK)］為聲矢量陣中M維聲壓陣的陣列流形矢量．其中，a(θk)=［1，e-jβk，…，e-j(M-1)βk］T為第k個(gè)信源的 導(dǎo) 向 矢 量;βk=2πdsin(θk)/λ;S(t)=［s1(t)，s2(t)，…，sK(t)］T為陣列接收的信號(hào)矢量;Np(t)為聲壓傳感器接收到的零均值高斯白噪聲，Nvx(t)為x方向振速傳感器接收到的零均值高斯白噪聲，Nvy(t)為y方向振速傳感器接收到的零均值高斯白噪聲，各個(gè)方向接收到的噪聲相互獨(dú)立;Φvx為x方向振速傳感器輸出的系數(shù)矩陣;Φvy為y方向振速傳感器輸出的系數(shù)矩陣，它們可表示為

1.2 聲壓振速聯(lián)合處理的測(cè)向方法

文獻(xiàn)［2］對(duì)矢量陣兩軸向的振速輸出進(jìn)行組合，可得到振速在某個(gè)參考方向θr上的電子旋轉(zhuǎn)矢量

其中

定義聲壓振速聯(lián)合信息處理的聲矢量陣P-V協(xié)方差矩陣為

式中:Rs=E［(S(t)SH(t))］，Rn為噪聲的協(xié)方差矩陣，Ⅰ為單位陣．

對(duì)R進(jìn)行特征值分解并將其特征向量按照特征值的大小降序排列可得

式中Un為矢量陣列噪聲子空間．則聲矢量陣聲壓振速組合MUSIC算法的空間譜為

上述算法既進(jìn)行了聲矢量陣聲壓振速聯(lián)合信息處理，又利用了矢量陣組合(P+Vc)Vc的指向性增益，因此獲得了較高的DOA估計(jì)性能.然而，由于R不在具有Toeplitz結(jié)構(gòu)，需通過(guò)一定的方法［3］對(duì)其進(jìn)行Toeplitz重構(gòu).

2 聲矢量陣快速子空間估計(jì)算法

2.1 MSWF算法求解信號(hào)子空間

文獻(xiàn)［6］指出，多級(jí)維納濾波器是一種有效的降維濾波技術(shù)，其在最小均方誤差意義下可得到維納霍夫方程的漸近最優(yōu)解而無(wú)需協(xié)方差矩陣的求逆運(yùn)算.對(duì)于標(biāo)量陣，其協(xié)方差矩陣為

而Wiener-Hoof方程RxWwf=rxd的漸進(jìn)最優(yōu)解為

基于相關(guān)相減的MSWF是一種有效的多級(jí)維納濾波結(jié)構(gòu)，其迭代過(guò)程如下:

步驟1 初始化d0(t)和x0(t)，

步驟2 前向遞推:令i=1，2，…，M;hi=Edi(k)==xi-1(t)-hidi(t)．

步驟3 后向遞推:eM(t)=dM(t)，

令i=M，M-1，…，1;

式中d0(t)為 MSWF算法中的期望信號(hào).文獻(xiàn)［6］指出，當(dāng)入射信號(hào)波形未知時(shí)，期望信號(hào)可取為參考陣元的輸出.令則MSWF的遞推過(guò)程等價(jià)于在Krylov子空間κ(M)(Rx，rxd)求解 Wiener Hopf方程［13］，經(jīng)M級(jí)遞推得到的各級(jí)濾波器的匹配濾波器h1，h2，…，hM構(gòu)成了M維Krylov子空間κ(M)(Rx，rxd)的一組基．文獻(xiàn)［7］指出，若rxd可表示為信號(hào)子空間對(duì)應(yīng)的所有特征向量的線性組合，則K維Krylo

v子空間κ(M)(Rx，rxd)等價(jià)于信號(hào)子空間．

2.2 新算法的原理與實(shí)現(xiàn)

基于聲壓振速聯(lián)合信息處理，取聲矢量陣互協(xié)方差矩陣為

給出期望信號(hào)d0(t)一種新的取法

定理 當(dāng)期望信號(hào)d0(t)取為聲矢量陣的電子旋轉(zhuǎn)矢量時(shí)，rxd可以表示為信號(hào)子空間對(duì)應(yīng)的所有特征向量的線性組合．

證明 期望信號(hào)和觀測(cè)數(shù)據(jù)的互相關(guān)函數(shù)為

令R's=ΦvcRs，則有

式中:R's為對(duì)角矩陣為入射信號(hào)功率．

由于聲矢量陣聲壓傳感器與振速傳感器接收到的噪聲互不相關(guān)，可得

將式(12)代入式(10)得

由上式可以看出，rxd為所有信號(hào)特征向量的線性組合，命題成立.

上述定理表明，當(dāng)取參考信號(hào)d0(t)=eTYvc(t)時(shí)，Krylov子空間 κ(K)(Rpv，rxd)等價(jià)于信號(hào)子空間．則經(jīng)K級(jí)遞推得到的各級(jí)濾波器的匹配濾波器h1，h2，…，hK即為K維Krylov子空間κ(K)(Rpv，rxd)的標(biāo)準(zhǔn)基，運(yùn)用MSWF快速子空間估計(jì)算法獲得聲矢量陣入射信號(hào)的子空間．聲矢量陣快速子空間方位估計(jì)算法的具體步驟如下:

步驟1 由式(3)得到參考方向θr上的電子旋轉(zhuǎn)矢量Yvc(t)．

步驟 2 取d0(t)=eTYvc(t)，x0(t)=Yp(t)．

步驟3 令i=1，2，…，K;hi=(t)］/

上式得到的各級(jí)濾波器系數(shù)H=［h1，h2，…，hK］張成信號(hào)子空間．

步驟4 由下式可得到聲矢量陣的空間譜估計(jì)為

2.3 算法性能分析

與文獻(xiàn)［2］相比本文的主要區(qū)別在于獲取噪聲子空間不同，文獻(xiàn)［2］算法通過(guò)計(jì)算陣列的協(xié)方差矩陣并對(duì)其進(jìn)行特征值分解來(lái)獲取噪聲子空間，需要的計(jì)算量為O(M2N+4M3/3)，而本文算法是通過(guò)MSWF遞推來(lái)獲得噪聲子空間，所需計(jì)算量?jī)H為O(KMN)，實(shí)際應(yīng)用中，入射信源數(shù)K遠(yuǎn)遠(yuǎn)小于陣元數(shù)M，故本文算法的計(jì)算量遠(yuǎn)小于文獻(xiàn)［2］算法．另外，由式(10)～(13)可以看出，本文提出的算法有效地抑制了各向同性噪聲，充分利用了聲壓振速信息聯(lián)合處理的良好抗噪能力，提高了聲矢量陣的處理增益．式(13)中，當(dāng)存在cos(θk-θr)≤0的入射信號(hào)時(shí)，R's將分別變?yōu)榉菨M秩矩陣和不定陣，文獻(xiàn)［3］進(jìn)行了深入的討論．

3 仿真與分析

3.1 計(jì)算量比較

假設(shè)2個(gè)信源入射到陣列，圖1(a)表示快拍數(shù)N=200時(shí)，文獻(xiàn)［2］算法與本文算法的計(jì)算量隨陣元數(shù)的變化曲線，隨著陣元數(shù)的增加，本文算法的計(jì)算量數(shù)值在一條斜率很小的直線上，計(jì)算量增加并不明顯，而文獻(xiàn)［2］算法的計(jì)算量呈指數(shù)增長(zhǎng)趨勢(shì)，陣元越多，計(jì)算量增加越明顯.在圖1(b)中，陣元數(shù)M=16，文獻(xiàn)［2］算法與本文算法的計(jì)算量隨快拍數(shù)的增加呈直線增長(zhǎng)趨勢(shì)，與文獻(xiàn)［2］算法相比，代表本文算法計(jì)算量直線的斜率很小，這說(shuō)明隨著快拍數(shù)的增加，文獻(xiàn)［2］算法計(jì)算量的增幅遠(yuǎn)大于本文算法.

圖1 計(jì)算量對(duì)比

3.2 空間譜估計(jì)性能比較

假設(shè)16個(gè)聲矢量傳感器沿x軸以d=λ/2等間距布放，其中λ為入射信號(hào)源在水中的波長(zhǎng)，3個(gè)互不相關(guān)的等功率高斯窄帶信號(hào)分別從10°，15°和40°方向入射，快拍數(shù)為1 000.

圖2為不同信噪比時(shí)，文獻(xiàn)［2，6］算法及本文算法的空間譜估計(jì)，從圖2(a)中可以看出，當(dāng)信噪比為15 dB時(shí)，3種算法均能準(zhǔn)確分辨出3個(gè)入射信號(hào)的方位，其中文獻(xiàn)［6］算法的“譜峰”較鈍，本文算法與文獻(xiàn)［2］算法的“譜峰”較尖銳，這表明高信噪比條件下，本文算法與文獻(xiàn)［2］算法的分辨性能相近，而文獻(xiàn)［6］算法的分辨性能較差;圖2(b)中，信噪比為5 dB，本文算法的“譜峰”略鈍于文獻(xiàn)［2］算法，但都能準(zhǔn)確分辨角度相近的2個(gè)入射信號(hào)，而文獻(xiàn)［6］算法卻未能成功分辨出入射角為10°的信號(hào)，這表明信噪比較低時(shí)，本文算法的分辨性能稍遜于文獻(xiàn)［2］算法，卻能夠分辨入射角度相近的信源，而文獻(xiàn)［6］算法的分辨性能嚴(yán)重下降，不能分辨入射角度相近的信源.

圖2 3種算法不同信噪比的空間譜

3.3 DOA估計(jì)性能比較

假設(shè)2個(gè)不相關(guān)等功率窄帶信號(hào)的入射角度分別為10°和20°，除信噪比和快拍數(shù)外，其他仿真條件同圖2，考察3種算法的DOA估計(jì)均方根誤差(Root Mean Square Error，RMSE)隨信噪比和快拍數(shù)的變化，定義DOA估計(jì)的均方根誤差為

式中:K為信號(hào)數(shù)目表示每次運(yùn)算得到的入射信號(hào)的DOA估計(jì)值，Ω為蒙特卡羅次數(shù)．

圖3(a)、(b)分別表示入射信號(hào)DOA估計(jì)的均方根誤差隨信噪比和快拍數(shù)的變化曲線，其中帶圓圈標(biāo)記的為文獻(xiàn)［6］算法的仿真結(jié)果，帶方塊標(biāo)記的為文獻(xiàn)［2］算法的仿真結(jié)果，帶星號(hào)標(biāo)記的為本文算法的仿真結(jié)果.在圖3(a)中，快拍數(shù)為200，橫軸為信噪比，從0 dB，間隔2 dB，變化到20 dB，縱軸為DOA估計(jì)的均方根誤差，每一個(gè)信噪比數(shù)據(jù)進(jìn)行500次蒙特卡羅獨(dú)立統(tǒng)計(jì).如圖3(a)所示，與文獻(xiàn)［6］算法相比，本文算法的θRMSE更小，當(dāng)RSN＜14 dB時(shí)，兩種算法的θRMSE差值較大，最大差值為0.25°，與文獻(xiàn)［2］算法相比，算法θRMSE略大，兩者的最大差值為0.12°，當(dāng)RSN＞12 dB時(shí)，兩種算法的θRMSE大致相等，這表明不同信噪比條件下，本文算法的DOA估計(jì)性能明顯優(yōu)于文獻(xiàn)［6］算法，當(dāng)信噪比較高時(shí)與文獻(xiàn)［2］算法的DOA估計(jì)性能相當(dāng).在圖3(b)中，信噪比為12 dB，橫軸為快拍數(shù)，從50，間隔50，變化到500，縱軸為DOA估計(jì)的均方根誤差，每一個(gè)快拍數(shù)進(jìn)行500次蒙特卡羅獨(dú)立統(tǒng)計(jì).如圖3(b)所示，本文算法的θRMSE比文獻(xiàn)［6］算法小0.05°，而只比文獻(xiàn)［2］算法大0.015°，且隨著快拍數(shù)的增加，3種算法的θRMSE變化不大，與圖2(a)結(jié)論相似，在不同快拍數(shù)條件下，本文算法的DOA估計(jì)性能明顯優(yōu)于文獻(xiàn)［6］算法，與文獻(xiàn)［2］算法性能相近.

圖3 DOA估計(jì)的均方根誤差曲線

圖2、3結(jié)果表明，與常規(guī)的MSWF算法相比，本文算法將聲壓振速聯(lián)合信息處理思想應(yīng)用到MSWF算法中獲得了更高的DOA估計(jì)性能，且在高信噪比條件下，本文算法的DOA估計(jì)性能接近文獻(xiàn)［2］算法，其原因是聲壓振速聯(lián)合信息處理抑制了各向同性噪聲，提高了聲矢量陣的處理增益.相比MUSIC類算法，基于MSWF的DOA估計(jì)算法已在硬件中實(shí)現(xiàn)并擁有較好的DOA估計(jì)性能，故本文算法更易于工程實(shí)現(xiàn).

4 結(jié)語(yǔ)

基于聲壓振速聯(lián)合信息處理，本文選擇參考陣元的電子旋轉(zhuǎn)矢量作為MSWF算法的期望信號(hào)，證明了由聲矢量陣觀測(cè)信號(hào)的互協(xié)方差矩陣及觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矢量構(gòu)成的Krylov子空間等價(jià)于信號(hào)子空間，并運(yùn)用MSWF算法快速求解 Krylov子空間的基.與常規(guī)的MSWF算法相比，本文算法在獲得更高的DOA估計(jì)性能的同時(shí)，沒(méi)有增加任何計(jì)算量;另外，本文算法的DOA估計(jì)性能稍遜于傳統(tǒng)子空間類聲矢量陣高分辨DOA估計(jì)算法，但是本文算法的計(jì)算量卻大為減小，這更有利于工程實(shí)現(xiàn)，尤其信噪比較高或陣元數(shù)目較多時(shí)，本文算法優(yōu)勢(shì)明顯.綜上所述，本文算法在實(shí)現(xiàn)聲矢量陣DOA估計(jì)中具有一定的工程應(yīng)用前景.

［1］JOHAM M，SUN Y，ZOLTOWSKI M D，et al.A new backward recursion for the multi-stage nested wiener filter employing Krylov subspace methods［J］. Military Communications Conference，2001，2(5):1210-1213.

［2］姚直象，胡金華，姚東明.基于多重信號(hào)分類法的一種聲矢量陣方位估計(jì)算法［J］.聲學(xué)學(xué)報(bào)，2008，33(4):305-309.

［3］白興宇，姜煜，趙春暉.基于聲壓振速聯(lián)合處理的聲矢量陣信源數(shù)檢測(cè)與方位估計(jì)［J］.聲學(xué)學(xué)報(bào)，2008(11)，33(1):56-61.

［4］姚直象，胡金華，姜可宇.矢量陣兩類陣處理方法研究［J］.兵工學(xué)報(bào)，2012(9)，33(9):1138-1142.

［5］姚直象，姜可宇，郭瑞，等.基于聲壓振速聯(lián)合處理的矢量旋轉(zhuǎn)不變子空間方位估計(jì)方法［J］.北京理工大學(xué)學(xué)報(bào)，2012(5)，32(5):513-517.

［6］黃磊.快速子空間估計(jì)方法研究及其在陣列信號(hào)處理中的應(yīng)用［D］.西安:西安電子科技大學(xué)，2005.

［7］安志娟，蘇洪濤，包志強(qiáng)，等.一種新的基于Krylov子空間的快速子空間分解［J］.系統(tǒng)工程與電子技術(shù)，2009(1)，31(1):29-33.

［8］周天，楊程，李海森，等.基于AccelDSP的 MM-MUSIC算法實(shí)現(xiàn)及其在多波束測(cè)深聲納中的應(yīng)用［J］.系統(tǒng)工程與電子技術(shù)，2011(12)，33(12):2613-2617.

［9］毛國(guó)華，王可人.一種基于DSP的快速M(fèi)USIC測(cè)向方法研究與實(shí)現(xiàn)［J］.電子信息對(duì)抗技術(shù)，2007(9)，22(5):10-13.

［10］王曉瑤，張國(guó)軍，王盼盼，等.聲矢量陣目標(biāo)估計(jì)的新方法［J］. 傳感器與微系統(tǒng)，2011(8)，30(8):73-76.

［11］莊學(xué)彬，陸名全，馮振明.一種數(shù)值穩(wěn)健且低復(fù)雜度的信號(hào)子空間估計(jì)新方法［J］.電子與信息學(xué)報(bào)，2011(1)，33(1):90-94.

［12］劉紅明，何子述，夏威，等.無(wú)參考信號(hào)條件下基于MSWF的 DOA 估計(jì)算法［J］.電子學(xué)報(bào)，2010(9)，38(9):1979-1983.

［13］JOHAM M，ZOLTOWSKI M D.Interpretation of multistage nested Wienerfilterin the Krylov subspace framework［J］.TechRepTR-ECE-00-51， Purdue University，USA，2000.