聲矢量陣快速子空間方位估計算法

梁國龍，張 柯，安少軍，范 展

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

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

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

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

1.1 聲矢量陣輸出模型

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

式中: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個信源的 導(dǎo) 向 矢 量;βk=2πdsin(θk)/λ;S(t)=［s1(t)，s2(t)，…，sK(t)］T為陣列接收的信號矢量;Np(t)為聲壓傳感器接收到的零均值高斯白噪聲，Nvx(t)為x方向振速傳感器接收到的零均值高斯白噪聲，Nvy(t)為y方向振速傳感器接收到的零均值高斯白噪聲，各個方向接收到的噪聲相互獨(dú)立;Φvx為x方向振速傳感器輸出的系數(shù)矩陣;Φvy為y方向振速傳感器輸出的系數(shù)矩陣，它們可表示為

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

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

其中

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

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

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

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

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

2 聲矢量陣快速子空間估計算法

2.1 MSWF算法求解信號子空間

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

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

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

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

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

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

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

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

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

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

令R's=ΦvcRs，則有

式中:R's為對角矩陣為入射信號功率．

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

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

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

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

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

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

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

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

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

2.3 算法性能分析

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

3 仿真與分析

3.1 計算量比較

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

圖1 計算量對比

3.2 空間譜估計性能比較

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

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

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

3.3 DOA估計性能比較

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

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

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

圖3 DOA估計的均方根誤差曲線

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

4 結(jié)語

基于聲壓振速聯(lián)合信息處理，本文選擇參考陣元的電子旋轉(zhuǎn)矢量作為MSWF算法的期望信號，證明了由聲矢量陣觀測信號的互協(xié)方差矩陣及觀測信號與期望信號的互相關(guān)矢量構(gòu)成的Krylov子空間等價于信號子空間，并運(yùn)用MSWF算法快速求解 Krylov子空間的基.與常規(guī)的MSWF算法相比，本文算法在獲得更高的DOA估計性能的同時，沒有增加任何計算量;另外，本文算法的DOA估計性能稍遜于傳統(tǒng)子空間類聲矢量陣高分辨DOA估計算法，但是本文算法的計算量卻大為減小，這更有利于工程實(shí)現(xiàn)，尤其信噪比較高或陣元數(shù)目較多時，本文算法優(yōu)勢明顯.綜上所述，本文算法在實(shí)現(xiàn)聲矢量陣DOA估計中具有一定的工程應(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］姚直象，胡金華，姚東明.基于多重信號分類法的一種聲矢量陣方位估計算法［J］.聲學(xué)學(xué)報，2008，33(4):305-309.

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

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

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

［6］黃磊.快速子空間估計方法研究及其在陣列信號處理中的應(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)及其在多波束測深聲納中的應(yīng)用［J］.系統(tǒng)工程與電子技術(shù)，2011(12)，33(12):2613-2617.

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

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

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

［12］劉紅明，何子述，夏威，等.無參考信號條件下基于MSWF的 DOA 估計算法［J］.電子學(xué)報，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.