王偉東， 李向水， 譚偉杰， 鄒波蓉， 李 輝

(1. 河南理工大學(xué) 物理與電子信息學(xué)院，河南 焦作 454000；2.貴州大學(xué) 大數(shù)據(jù)產(chǎn)業(yè)發(fā)展應(yīng)用研究院 公共大數(shù)據(jù)國家重點實驗室，貴陽 550025)

隨著隱身技術(shù)的不斷發(fā)展，借助消聲瓦新型材料，武器裝備，例如魚雷、水雷、水面艦等，輻射噪聲的頻率與強度逐年降低，給被動探測設(shè)備帶來了困難[1-2]。為了提高遠程低頻目標的方位估計性能，通常需要增加聲壓傳感器陣列的孔徑。然而，對于空間受限的探測平臺，想要增加陣列的孔徑通常是比較困難的。矢量傳感器陣列是解決上述難題的理想方案。較之僅能測量聲場中聲壓信息的聲壓傳感器，矢量傳感器能共點、同步、獨立地測量聲場中的聲壓和振速信息，通過測量的聲壓和振速信息可實現(xiàn)全空間范圍內(nèi)目標方位角的無左右舷模糊估計[3-4]，為空間受限平臺探測系統(tǒng)進行遠程低頻目標的探測奠定了基礎(chǔ)。

Nehorai等[5]利用克羅內(nèi)克乘積將矢量傳感器及其陣列信號模型與傳統(tǒng)的聲壓傳感器陣列信號模型統(tǒng)一以來，基于已日趨成熟的聲壓傳感器陣列方位估計方法，如聲能流方法、波束形成方法[6]、高分辨子空間類方法[7-8]等，廣泛被應(yīng)用矢量傳感器陣列的目標方位估計中。聲能流方法及波束形成方法具有結(jié)構(gòu)簡單、計算量小的優(yōu)勢，但是兩種方法的分辨性能較差。子空間類方法，諸如多重信號分類(multiple signal classification, MUSIC)方法[9]和旋轉(zhuǎn)不變技術(shù)估計信號參數(shù)方法[10]及它們的改進方法[11-12]，通過利用信號子空間和噪聲子空間的正交性突破了波束形成方法無法實現(xiàn)一個波束范圍內(nèi)兩目標的分辨這一弊端，實現(xiàn)了目標的高分辨方位估計。然而，該類方法對信號子空間與噪聲子空間的正交性具有較強的敏感性使其處理相干信號時方位估計性能急劇下降。在傳統(tǒng)的聲壓傳感器陣列信號處理中，空間平滑技術(shù)[13]是解決相干信號方位估計的有效方法之一。在陣列孔徑較大時，該方法能夠被應(yīng)用于矢量傳感器陣列的相干信號處理中。但是，空間平滑技術(shù)是以損失陣列的孔徑為代價來提高相干信號的方位估計性能。為了解決此難題，文獻[14-16]利用矢量傳感器固有的指向性，相繼提出了振速域平滑方法、振速域差分平滑方法及無左右舷模糊的平滑方法。相較于空間平滑技術(shù)，盡管振速類平滑方法在未損失陣列孔徑的前提下提高了相干信號源的方位估計性能，但是，矢量傳感器陣列面臨損失自由度的難題。

稀疏信號處理方法是最近十幾年興起的方位估計方法。早期，多數(shù)稀疏表示方位估計方法[17-18]是基于聲壓傳感器陣列進行的研究。該類算法對信號之間相關(guān)性的敏感度較低，被應(yīng)用于矢量傳感器陣列方位估計時能解決矢量傳感器陣列損失自由度的難題。國內(nèi)外學(xué)者對矢量傳感器陣列稀疏表示方位估計方法的研究成果并不多見。針對低信噪比和小快拍下高速運動目標的穩(wěn)健高分辨方位估計問題，文獻[19]提出了基于矢量傳感器陣列的穩(wěn)健高分辨方位估計方法。然而，該方法的方位估計性能受噪聲功率門限值影響較大。此外，該方法涉及到聲壓與振速的互相關(guān)操作，導(dǎo)致矢量傳感器陣列損失自由度的難題。針對各向同性環(huán)境噪聲下矢量傳感器陣列相干信號方位估計性能下降問題，文獻[20]提出了基于增廣互協(xié)方差矩陣矢量稀疏表示的方位估計方法。盡管該方法的估計精度有所提升，但是矢量傳感器陣列損失自由度的難題依然存在。為了解決小孔徑矢量傳感器陣列損失自由度的難題，文獻[21]提出了稀疏信號功率迭代補償?shù)氖噶總鞲衅麝嚵蟹轿还烙嫹椒?。在該方法中，當前一次迭代估計的稀疏信號功率是基于前一次稀疏信號功率的估計結(jié)果進行的補償，提高了矢量傳感器陣列在信噪比較低和角度間隔較小時的方位估計性能。

然而，上述方法都是基于陣列中振速軸指向同一方向的理想假設(shè)，對矢量傳感器陣列的方位估計性能進行的研究。由于矢量傳感器具有偶極子指向性，在安裝矢量傳感器陣列時，難以保障陣列中的振速軸指向同一方向，簡稱為振速軸向不一致，使矢量傳感器陣列的方位估計性能下降。文獻[22]在矢量傳感器陣列信號模型中引入了一個軸向角度偏差參數(shù)，在迭代自適應(yīng)方法(iterative adaptive approach，IAA)的基礎(chǔ)上建立了關(guān)系稀疏信號和軸向偏差矩陣的代價函數(shù)。然后提出了一種交替迭代自適應(yīng)方法(alternating iterative adaptive approach，AIAA)，該方法基于貝葉斯信息準則對信號進行補償，進一步提高了軸向偏差參數(shù)偏差模型下矢量傳感器陣列的方位估計精度。然而，該方法是基于矢量傳感器陣列輸出的一階統(tǒng)計量進行的研究。相較于一階統(tǒng)計量，矢量傳感器陣列輸出的二階統(tǒng)計量具有更多的自由度和輸出信噪比。

因此，針對振速軸向不一致時矢量傳感器陣列方位估計性能惡化問題，提出了一種兩步加權(quán)交替迭代方法(two-step weighted alternating iterative approach, TWAIA)。首先，基于矢量傳感器陣列的二階統(tǒng)計量，采用正則化加權(quán)協(xié)方差矩陣擬合方法估計稀疏信號功率。其次，采用正則化加權(quán)最小二乘估計軸向角度偏差矩陣。在每次迭代中，為了提高稀疏信號功率在空域的稀疏性，稀疏信號功率補償項被約束；為了消除軸向角度偏差對方位估計性能的影響，基于矢量傳感器固有的指向性，期望的軸向角度偏差矩陣被重構(gòu)。仿真結(jié)果表明，與現(xiàn)有方法相比，該方法提高了振速軸向不一致時矢量傳感器陣列的方位估計精度。

1 問題描述

1.1 理想矢量傳感器陣列數(shù)據(jù)模型

h(θk)=[1 cosθksinθk]T

(1)

式中，θk∈(-180°,180°]為第k個信號源的水平方位角。

圖1 理想矢量傳感器均勻線陣模型Fig.1 The ideal uniform vector sensor array model

于是，由M個矢量傳感器組成均勻線列陣的陣列流形矩陣為

A(θ)=[a(θ1),…,a(θK)]=

[ap(θ1)?h(θ1),…,ap(θK)?h(θK)]=

Ap(θ)⊙H(θ)

(2)

式中:Ap(θ)=[ap(θ1),ap(θ2),…，ap(θK)]為聲壓傳感器陣列的陣列流形矩陣；H(θ)=[h(θ1),h(θ2),…，h(θK)]為各聲源在單矢量傳感器上的響應(yīng)；ap(θk)=[1,…,e-j2π(M-1)dcos θk/λ]T為第k個聲源在聲壓傳感器陣列上的響應(yīng)；a(θk)為第k個聲源在矢量傳感器陣列上的響應(yīng)；λ為聲波的波長；?為Kronecker積；⊙為Khatri-Rao積。于是，矢量傳感器陣列在t時刻的輸出矢量為

Y=A(θ)S+W

(4)

式中:Y=[y(1),y(2),…,y(L)]；S=[s(1),s(2),…,s(L)]；W=[w(1),w(2),…,w(L)]。

(5)

1.2 振速軸向不一致時矢量傳感器陣列模型

由于矢量傳感器本身固有的偶極子指向性，在安裝矢量傳感器陣列時通常會存在軸向角度偏差，即陣列中的振速軸沒有指向同一個方向，如圖2所示，其中，βm為第m個矢量傳感器的軸向角度偏差，β=[β1,β2,…,βM]T為矢量傳感器陣列的振速軸向偏差矢量。

圖2 振速軸向不一致下矢量傳感器均勻線列陣模型Fig.2 The uniform vector sensor array model with velocity axial inconsistency

因此，當?shù)趉個聲源入射在第m個矢量傳感器上時，式(1)可修正為

h(θk,βm)=[1·cos(θk-βm) sin(θk-βm)]T=

gmh(θk)

(6)

式中，gm為第m個矢量傳感器的軸向角度偏差矩陣

(7)

振速軸向不一致下第k個聲源在矢量傳感器陣列上可表示為

a(θk,β)=[ap(θk)?h(θk,β1),…,ap(θk)?h(θk,βM)]=

[a1p(θk)g1h(θk),…,aMp(θk)gMh(θk)]=

GAp(θk)F(θk)

(8)

式中:Ap(θk)=diag{apk?I3×1}；diag{·}為對角矩陣操作；apk=[a1k(θk),…,aMk(θk)]T；F(θk)=IM×1?h(θk)；Ii×1為i行的單位矢量；G為矢量傳感器陣列的軸向角度偏差矩陣

G=blkdiag{g1,g2,…,gM}

(9)

式中，blkdiag{·}為塊對角矩陣操作。

振速軸向不一致下，矢量傳感器陣列的稀疏矢量數(shù)據(jù)模型式(5)，可修正為

(10)

2 兩步加權(quán)交替迭代自適應(yīng)矢量傳感器陣列方位估計方法

根據(jù)式(10)，在無噪聲情況下，第n個網(wǎng)格相應(yīng)的信號協(xié)方差矩陣可表示為

(11)

(12)

于是，相應(yīng)于第n個網(wǎng)格上信號的干擾加噪聲協(xié)方差矩陣為

Σn=R-Rn

(13)

式中，R為信號加噪聲協(xié)方差矩陣，

(14)

(15)

其中

(16)

G(j+1)=arg minFG(G)

s.t.GGH=GHG=I

(17)

其中

式中：λg為正則化參數(shù)，平衡著軸向角度偏差矩陣與擬合誤差之間的關(guān)系。

2.1 估計稀疏信號功率

假設(shè)F(x)=xq是關(guān)于變量x的函數(shù)。不難看出，F(xiàn)(x)也是關(guān)于變量x的凹函數(shù)。當變量x>0時，根據(jù)文獻[24]可知，對于變量x0>0，可得

F(x)≤F(x0)+F′(x0)(x-x0)

(19)

(20)

(21)

(22)

(23)

(24)

其中

(25)

(26)

(27)

其中

(28)

(29)

(31)

(32)

(33)

對式(27)的解析表達式和所提算法的收斂性進行推導(dǎo)。

(1)式(27)的證明

(34)

其中

(37)

(39)

其中

(40)

從式(39)中可以看出，在計算N個離散網(wǎng)格上對應(yīng)的稀疏信號功率時，共需計算N次干擾加噪聲協(xié)方差矩陣Σn的逆，因此，計算量較大。為了減少干擾加噪聲協(xié)方差矩陣Σn的逆運算，根據(jù)Woodbury公式，可得

(41)

其中

(42)

把式(41)代入式(39)，可得

(43)

其中

(44)

證畢。

(2)收斂性證明

易看出，若是所提方法具有收斂性，需要滿足如下條件

(45)

式(45)的證明過程如下。

(46)

根據(jù)式(16)、式(24)及式(25)，可知

(47)

(48)

根據(jù)式(46)、式(47)和式(48)，可得式(45)是成立的。因此，所提算法具有收斂性。

證畢。

2.2 估計軸向角度偏差矩陣

類似于式(15)到式(26)的轉(zhuǎn)化，式(16)同樣可以轉(zhuǎn)化為關(guān)于軸向角度偏差矩陣G的線性函數(shù)，即最小化式(16)關(guān)于軸向角度偏差矩陣G等價于最小化下式，

G=arg minH(G)
s.t.GGH=GHG=I

(49)

其中

(50)

為了解決式(49)的優(yōu)化問題，首先通過最小化式(50)獲得軸向角度偏差矩陣G的粗略估計。然后，利用矢量傳感器固有的方向性，對粗略估計的軸向角度偏差矩陣進行重構(gòu)使其滿足期望的軸向角度偏差矩陣。可以看出，通過最小化式(50)求解的軸向角度偏差矩陣粗略估計是與所有網(wǎng)格上相應(yīng)的來波方向有關(guān)系的，為了進一步提高軸向角度偏差矩陣的估計精度，僅利用與K個目標方位對應(yīng)的信號求解軸向角度偏差矩陣。則式(50)進一步轉(zhuǎn)化為

(51)

(52)

其中

(53)

(54)

2.3 重構(gòu)期望的軸向角度偏差矩陣

(55)

式中，η=3m-1。假設(shè)第m個矢量傳感器的振速軸向作為參考方向，可求得相對軸向角度偏差矩陣為

(56)

則估計的第m個矢量傳感器的指向性角度偏差可表示為

(57)

(58)

則矢量傳感器陣列的軸向角度偏差矩陣可表示為

(59)

式中，軸向角度偏差矩陣初始化為

(60)

2.4 所提算法總結(jié)

所提算法總結(jié)如下。

第三步：根據(jù)式(30)，更新R(j)。

2.5 所提算法計算復(fù)雜度分析

對所提TWAIA算法計算復(fù)雜度進行分析，這里假設(shè)迭代次數(shù)為Titer(Titer≤Tmax)。所提算法初始化部分需MN+(2L+1)N次乘法與MN+(2L-1)N次加法，計算協(xié)方差矩陣R的逆需M2N+MN+21M3次乘法與M2N+21M3次加法，計算網(wǎng)格上對應(yīng)稀疏信號功率需7M2N+5MN+3N次乘法與6M2N+2MN+2N次加法，估計和重構(gòu)軸向角度偏差矩陣需12M2N+3MN+KN次乘法與4M2N+MN次加法。所提算法的計算量為O{31TiterM2N+(12Titer+2)MN+(5Titer+4L)N+TiterKN+42M3}。MUSIC不需要進行逆矩陣和迭代運算，計算量為O{M2L+M3N}。SAMV-1和AIAA的計算量分別為O{15TiterM2N+(5Titer+2)MN+(2Titer+4L)N+42M3}和O{(13Titer+2)M2N+(8Titer+7)MN+(3Titer+2L)MN+TiterKN+42M3}。因此，相較于MUSIC，IAA和AIAA算法，所提TWAIA算法具有較高的計算量，但是在比較的幾種方法中，TWAIA算法的方位估計精度是最高的。

3 仿真與性能分析

在本章中，將所提TWAIA方法與經(jīng)典的 MUSIC方法、SAMV-1方法及AIAA方法進行比較，評估所提方法的有效性與穩(wěn)健性，同時給出了矢量傳感器陣列的克拉美羅下界(Cramer-Rao lower bound, CRLB)曲線，以此作為比較的基準。所有方法都是基于一個由M元矢量傳感器組成的均勻線列陣進行的研究，并且相鄰陣元間的間距為一個波長。除了位于坐標原點處矢量傳感器的振速軸方向作為參考之外，在每次蒙特卡洛仿真中，其它矢量傳感器的振速軸向相較于參考陣元的軸向角度偏差滿足均值為β，方差為ρ2的均勻分布。等功率遠場窄帶不相干的信號為復(fù)指數(shù)信號，其中心頻率f=711 Hz，系統(tǒng)采樣頻率為fs=8 kHz，聲源的角度掃描范圍為[-180°∶2°∶180°]。在所提方法中，τ=1×10-3，Tmax=20，根據(jù)Wang等的研究和文獻[27]的分析，υ=45，λs=λg=0.625，q=0.5。下文仿真中把均方根誤差作為估計精度的衡量標準，并定義均方根誤差(root mean square error, RMSE)為

(61)

試驗一：均方根誤差與信噪比的關(guān)系

假設(shè)兩個等功率遠場不相干的信號源入射在矢量傳感器陣列上，圖3給出了兩個信號源的方位角分別位于θ1=-14°，θ2=38°，陣元個數(shù)為4，快拍數(shù)為300，β=15，ρ=4時，均方根誤差與信噪比的關(guān)系。從圖3可以看出，隨著信噪比的增加，MUSIC的均方根誤差相較于CRLB具有較大的間隔，且MUISC的均方根誤差變化趨勢較緩慢，這主要是由于MUSIC未考慮軸向角度偏差對矢量傳感器陣列方位估計的影響。盡管SAMV-1方法也未考慮軸向角度偏差對矢量傳感器陣列方位估計的影響，但是相較于MUSIC方法，SAMV-1方法與CRLB之間的間隔有所減小，這表明SAMV-1方法對矢量傳感器陣列軸向角度偏差具有較低的敏感性。相較于MUAIC方法和SAMV-1方法，盡管AIAA通過估計軸向角度偏差矩陣對觀測數(shù)據(jù)進行了校準，減小了它對目標方位角度估計性能的影響。顯著的是，AIAA方法在信噪比較高時，它的估計性能接近于CRLB，但在信噪比較低時，它的估計性能仍較差。這主要是由于AIAA方法在估計稀疏信號時是根據(jù)陣列的一階統(tǒng)計特性進行的估計，加之軸向角度偏差矩陣采用最小二乘估計未考慮軸向角度偏差矩陣與殘差之間的關(guān)系導(dǎo)致AIAA方法在低信噪比時方位估計精度不高。相反，所提方法的均方根誤差與CRLB之間具有較小的間隔，這主要歸因于所提方法是依據(jù)矢量傳感器陣列輸出的二階統(tǒng)計量估計稀疏信號功率，并且采用正則化加權(quán)最小二乘估計軸向角度偏差矩陣，并通過戶用參數(shù)較好的平衡了信號在空域的稀疏性及軸向角度偏差矩陣與殘差之間的關(guān)系，使得所提方法在低信噪比下矢量傳感器陣列的方位估計精度得以提高。

圖3 均方根誤差與信噪比的關(guān)系Fig.3 RMSE versus SNR

試驗二：均方根誤差與快拍數(shù)的關(guān)系

圖4比較了均方根誤差與快拍數(shù)的關(guān)系，除了信噪比為8 dB、快拍數(shù)在[20，40，60，80，100∶100∶600]內(nèi)變化之外，其他仿真參數(shù)的設(shè)置與圖3保持一致。從圖4中可以看出，不論快拍數(shù)較小還是較大，MUSIC方法的均方根誤差都較大。相較于MUSIC方法，盡管SAMV-1方法的估計精度有很大提升，但是與CRLB之間仍有較大的間隔。在快拍數(shù)小于40時，AIAA方法方位估計性能惡化，這主要是由于在快拍數(shù)較小時，恢復(fù)的稀疏信號不夠精確，導(dǎo)致最終估計的稀疏信號功率誤差較大，從而致使方位估計誤差較大。所提方法依據(jù)矢量傳感器陣列輸出的二階統(tǒng)計量估計稀疏信號功率，并采用正則化參數(shù)及用戶參數(shù)進行約束有效的提高了快拍數(shù)不足時AIAA方法存在的弊端，并且當快拍數(shù)充足時所提方法的均方根誤差曲線近似于CRLB，表明了所提方法具有較好的方位估計性能。

圖4 均方根誤差與快拍數(shù)的關(guān)系Fig.4 RMSE versus the number of snapshot

試驗三：均方根誤差與陣元個數(shù)的關(guān)系

圖5比較了所提方法與現(xiàn)有方位估計方法的均方根誤差與矢量傳感器個數(shù)的關(guān)系，其中，信噪比為8 dB，陣元個數(shù)從2到7變化，其他仿真參數(shù)與圖3保持一致?？梢钥吹剑S著陣元數(shù)的增加，MUSIC方法、SAMV-1方法、AIAA方法及所提方法的估計性能逐漸提高。MUSIC方法在陣元個數(shù)為8時無法達到無偏估計。SAMV-1方法在陣元個數(shù)為7時可實現(xiàn)無偏估計。AIAA方法及所提方法在陣元個數(shù)為5時即可實現(xiàn)無偏估計。相較于未校準的MUSIC方法和SAMV-1方法，在無偏估計性能的情況下，所提方法與AIAA方法需要的陣列孔徑最小。此外，在具有相同孔徑的條件下，相較于AIAA方法所提方法具有更高的估計精度，這充分展現(xiàn)了所提方法在小孔徑時優(yōu)越的方位估計性能。

圖5 均方根誤差與陣元個數(shù)的關(guān)系Fig.5 RMSE versus the number of snapshot

試驗四：均方根誤差與偏差均值的關(guān)系

圖6比較了所提方法與現(xiàn)有方位估計方法的均方根誤差與矢量傳感器陣列軸向角度偏差均值之間的關(guān)系，其中，信噪比為8 dB，β從5到35變化，其他仿真參數(shù)與圖3保持一致。可以看出，幾種方法的均方根誤差都會隨著β的增加而增加。相較于未考慮陣列軸向角度偏差參數(shù)的MUSIC方法和SAMV-1方法，AIAA方法和所提方法的均方根誤差隨軸向角度偏差均值增大的速率較大。盡管AIAA方法和所提方法的均方根誤差隨軸向角度偏差均值增大的速率較小，但是較于AIAA方法，所提方法的均方根誤差更小，這展示出所提方法具有更好的方位估計性能。

圖6 均方根誤差與軸向角度偏差均值之間的關(guān)系Fig.6 RMSE versus the mean value of axial angle bias

4 結(jié) 論

為了提高低信噪比振速軸向不一致時矢量傳感器陣列的方位估計精度，提出了一種兩步加權(quán)交替迭代自適應(yīng)的矢量傳感器陣列方位估計方法。該方法采用正則化加權(quán)稀疏協(xié)方差矩陣擬合方法估計稀疏信號功率，通過對稀疏信號功率進行補償提高了其在空域的稀疏性；利用加權(quán)最小二乘估計軸向角度偏差矩陣，并根據(jù)矩陣中偏差的分布特性，重構(gòu)了期望的軸向角度偏差矩陣，以此交替的方式迭代更新稀疏信號功率和軸向角度偏差矩陣，使得最終估計的稀疏信號功率更加精確，從而提高了振速軸向不一致時矢量傳感器陣列的方位估計精度。仿真結(jié)果表明，在低信噪比、少快拍、小孔徑下，相較于MUSIC方法、SAMV-1方法和AIAA方法，所提方法提高了振速軸向不一致時矢量傳感器陣列的方位估計精度。