王偉東, 李向水, 譚偉杰, 鄒波蓉, 李 輝
(1. 河南理工大學(xué) 物理與電子信息學(xué)院,河南 焦作 454000;2.貴州大學(xué) 大數(shù)據(jù)產(chǎn)業(yè)發(fā)展應(yīng)用研究院 公共大數(shù)據(jù)國家重點(diǎn)實(shí)驗(yàn)室,貴陽 550025)
隨著隱身技術(shù)的不斷發(fā)展,借助消聲瓦新型材料,武器裝備,例如魚雷、水雷、水面艦等,輻射噪聲的頻率與強(qiáng)度逐年降低,給被動(dòng)探測設(shè)備帶來了困難[1-2]。為了提高遠(yuǎn)程低頻目標(biāo)的方位估計(jì)性能,通常需要增加聲壓傳感器陣列的孔徑。然而,對于空間受限的探測平臺,想要增加陣列的孔徑通常是比較困難的。矢量傳感器陣列是解決上述難題的理想方案。較之僅能測量聲場中聲壓信息的聲壓傳感器,矢量傳感器能共點(diǎn)、同步、獨(dú)立地測量聲場中的聲壓和振速信息,通過測量的聲壓和振速信息可實(shí)現(xiàn)全空間范圍內(nèi)目標(biāo)方位角的無左右舷模糊估計(jì)[3-4],為空間受限平臺探測系統(tǒng)進(jìn)行遠(yuǎn)程低頻目標(biāo)的探測奠定了基礎(chǔ)。
Nehorai等[5]利用克羅內(nèi)克乘積將矢量傳感器及其陣列信號模型與傳統(tǒng)的聲壓傳感器陣列信號模型統(tǒng)一以來,基于已日趨成熟的聲壓傳感器陣列方位估計(jì)方法,如聲能流方法、波束形成方法[6]、高分辨子空間類方法[7-8]等,廣泛被應(yīng)用矢量傳感器陣列的目標(biāo)方位估計(jì)中。聲能流方法及波束形成方法具有結(jié)構(gòu)簡單、計(jì)算量小的優(yōu)勢,但是兩種方法的分辨性能較差。子空間類方法,諸如多重信號分類(multiple signal classification, MUSIC)方法[9]和旋轉(zhuǎn)不變技術(shù)估計(jì)信號參數(shù)方法[10]及它們的改進(jìn)方法[11-12],通過利用信號子空間和噪聲子空間的正交性突破了波束形成方法無法實(shí)現(xiàn)一個(gè)波束范圍內(nèi)兩目標(biāo)的分辨這一弊端,實(shí)現(xiàn)了目標(biāo)的高分辨方位估計(jì)。然而,該類方法對信號子空間與噪聲子空間的正交性具有較強(qiáng)的敏感性使其處理相干信號時(shí)方位估計(jì)性能急劇下降。在傳統(tǒng)的聲壓傳感器陣列信號處理中,空間平滑技術(shù)[13]是解決相干信號方位估計(jì)的有效方法之一。在陣列孔徑較大時(shí),該方法能夠被應(yīng)用于矢量傳感器陣列的相干信號處理中。但是,空間平滑技術(shù)是以損失陣列的孔徑為代價(jià)來提高相干信號的方位估計(jì)性能。為了解決此難題,文獻(xiàn)[14-16]利用矢量傳感器固有的指向性,相繼提出了振速域平滑方法、振速域差分平滑方法及無左右舷模糊的平滑方法。相較于空間平滑技術(shù),盡管振速類平滑方法在未損失陣列孔徑的前提下提高了相干信號源的方位估計(jì)性能,但是,矢量傳感器陣列面臨損失自由度的難題。
稀疏信號處理方法是最近十幾年興起的方位估計(jì)方法。早期,多數(shù)稀疏表示方位估計(jì)方法[17-18]是基于聲壓傳感器陣列進(jìn)行的研究。該類算法對信號之間相關(guān)性的敏感度較低,被應(yīng)用于矢量傳感器陣列方位估計(jì)時(shí)能解決矢量傳感器陣列損失自由度的難題。國內(nèi)外學(xué)者對矢量傳感器陣列稀疏表示方位估計(jì)方法的研究成果并不多見。針對低信噪比和小快拍下高速運(yùn)動(dòng)目標(biāo)的穩(wěn)健高分辨方位估計(jì)問題,文獻(xiàn)[19]提出了基于矢量傳感器陣列的穩(wěn)健高分辨方位估計(jì)方法。然而,該方法的方位估計(jì)性能受噪聲功率門限值影響較大。此外,該方法涉及到聲壓與振速的互相關(guān)操作,導(dǎo)致矢量傳感器陣列損失自由度的難題。針對各向同性環(huán)境噪聲下矢量傳感器陣列相干信號方位估計(jì)性能下降問題,文獻(xiàn)[20]提出了基于增廣互協(xié)方差矩陣矢量稀疏表示的方位估計(jì)方法。盡管該方法的估計(jì)精度有所提升,但是矢量傳感器陣列損失自由度的難題依然存在。為了解決小孔徑矢量傳感器陣列損失自由度的難題,文獻(xiàn)[21]提出了稀疏信號功率迭代補(bǔ)償?shù)氖噶總鞲衅麝嚵蟹轿还烙?jì)方法。在該方法中,當(dāng)前一次迭代估計(jì)的稀疏信號功率是基于前一次稀疏信號功率的估計(jì)結(jié)果進(jìn)行的補(bǔ)償,提高了矢量傳感器陣列在信噪比較低和角度間隔較小時(shí)的方位估計(jì)性能。
然而,上述方法都是基于陣列中振速軸指向同一方向的理想假設(shè),對矢量傳感器陣列的方位估計(jì)性能進(jìn)行的研究。由于矢量傳感器具有偶極子指向性,在安裝矢量傳感器陣列時(shí),難以保障陣列中的振速軸指向同一方向,簡稱為振速軸向不一致,使矢量傳感器陣列的方位估計(jì)性能下降。文獻(xiàn)[22]在矢量傳感器陣列信號模型中引入了一個(gè)軸向角度偏差參數(shù),在迭代自適應(yīng)方法(iterative adaptive approach,IAA)的基礎(chǔ)上建立了關(guān)系稀疏信號和軸向偏差矩陣的代價(jià)函數(shù)。然后提出了一種交替迭代自適應(yīng)方法(alternating iterative adaptive approach,AIAA),該方法基于貝葉斯信息準(zhǔn)則對信號進(jìn)行補(bǔ)償,進(jìn)一步提高了軸向偏差參數(shù)偏差模型下矢量傳感器陣列的方位估計(jì)精度。然而,該方法是基于矢量傳感器陣列輸出的一階統(tǒng)計(jì)量進(jìn)行的研究。相較于一階統(tǒng)計(jì)量,矢量傳感器陣列輸出的二階統(tǒng)計(jì)量具有更多的自由度和輸出信噪比。
因此,針對振速軸向不一致時(shí)矢量傳感器陣列方位估計(jì)性能惡化問題,提出了一種兩步加權(quán)交替迭代方法(two-step weighted alternating iterative approach, TWAIA)。首先,基于矢量傳感器陣列的二階統(tǒng)計(jì)量,采用正則化加權(quán)協(xié)方差矩陣擬合方法估計(jì)稀疏信號功率。其次,采用正則化加權(quán)最小二乘估計(jì)軸向角度偏差矩陣。在每次迭代中,為了提高稀疏信號功率在空域的稀疏性,稀疏信號功率補(bǔ)償項(xiàng)被約束;為了消除軸向角度偏差對方位估計(jì)性能的影響,基于矢量傳感器固有的指向性,期望的軸向角度偏差矩陣被重構(gòu)。仿真結(jié)果表明,與現(xiàn)有方法相比,該方法提高了振速軸向不一致時(shí)矢量傳感器陣列的方位估計(jì)精度。
h(θk)=[1 cosθksinθk]T
(1)
式中,θk∈(-180°,180°]為第k個(gè)信號源的水平方位角。
圖1 理想矢量傳感器均勻線陣模型Fig.1 The ideal uniform vector sensor array model
于是,由M個(gè)矢量傳感器組成均勻線列陣的陣列流形矩陣為
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個(gè)聲源在聲壓傳感器陣列上的響應(yīng);a(θk)為第k個(gè)聲源在矢量傳感器陣列上的響應(yīng);λ為聲波的波長;?為Kronecker積;⊙為Khatri-Rao積。于是,矢量傳感器陣列在t時(shí)刻的輸出矢量為
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)
由于矢量傳感器本身固有的偶極子指向性,在安裝矢量傳感器陣列時(shí)通常會存在軸向角度偏差,即陣列中的振速軸沒有指向同一個(gè)方向,如圖2所示,其中,βm為第m個(gè)矢量傳感器的軸向角度偏差,β=[β1,β2,…,βM]T為矢量傳感器陣列的振速軸向偏差矢量。
圖2 振速軸向不一致下矢量傳感器均勻線列陣模型Fig.2 The uniform vector sensor array model with velocity axial inconsistency
因此,當(dāng)?shù)趉個(gè)聲源入射在第m個(gè)矢量傳感器上時(shí),式(1)可修正為
h(θk,βm)=[1·cos(θk-βm) sin(θk-βm)]T=
gmh(θk)
(6)
式中,gm為第m個(gè)矢量傳感器的軸向角度偏差矩陣
(7)
振速軸向不一致下第k個(gè)聲源在矢量傳感器陣列上可表示為
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)
根據(jù)式(10),在無噪聲情況下,第n個(gè)網(wǎng)格相應(yīng)的信號協(xié)方差矩陣可表示為
(11)
(12)
于是,相應(yīng)于第n個(gè)網(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)系。
假設(shè)F(x)=xq是關(guān)于變量x的函數(shù)。不難看出,F(xiàn)(x)也是關(guān)于變量x的凹函數(shù)。當(dāng)變量x>0時(shí),根據(jù)文獻(xiàn)[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á)式和所提算法的收斂性進(jìn)行推導(dǎo)。
(1)式(27)的證明
(34)
其中
(37)
(39)
其中
(40)
從式(39)中可以看出,在計(jì)算N個(gè)離散網(wǎng)格上對應(yīng)的稀疏信號功率時(shí),共需計(jì)算N次干擾加噪聲協(xié)方差矩陣Σn的逆,因此,計(jì)算量較大。為了減少干擾加噪聲協(xié)方差矩陣Σn的逆運(yùn)算,根據(jù)Woodbury公式,可得
(41)
其中
(42)
把式(41)代入式(39),可得
(43)
其中
(44)
證畢。
(2)收斂性證明
易看出,若是所提方法具有收斂性,需要滿足如下條件
(45)
式(45)的證明過程如下。
(46)
根據(jù)式(16)、式(24)及式(25),可知
(47)
(48)
根據(jù)式(46)、式(47)和式(48),可得式(45)是成立的。因此,所提算法具有收斂性。
證畢。
類似于式(15)到式(26)的轉(zhuǎn)化,式(16)同樣可以轉(zhuǎn)化為關(guān)于軸向角度偏差矩陣G的線性函數(shù),即最小化式(16)關(guān)于軸向角度偏差矩陣G等價(jià)于最小化下式,
G=arg minH(G)
s.t.GGH=GHG=I
(49)
其中
(50)
為了解決式(49)的優(yōu)化問題,首先通過最小化式(50)獲得軸向角度偏差矩陣G的粗略估計(jì)。然后,利用矢量傳感器固有的方向性,對粗略估計(jì)的軸向角度偏差矩陣進(jìn)行重構(gòu)使其滿足期望的軸向角度偏差矩陣??梢钥闯?,通過最小化式(50)求解的軸向角度偏差矩陣粗略估計(jì)是與所有網(wǎng)格上相應(yīng)的來波方向有關(guān)系的,為了進(jìn)一步提高軸向角度偏差矩陣的估計(jì)精度,僅利用與K個(gè)目標(biāo)方位對應(yīng)的信號求解軸向角度偏差矩陣。則式(50)進(jìn)一步轉(zhuǎn)化為
(51)
(52)
其中
(53)
(54)
(55)
式中,η=3m-1。假設(shè)第m個(gè)矢量傳感器的振速軸向作為參考方向,可求得相對軸向角度偏差矩陣為
(56)
則估計(jì)的第m個(gè)矢量傳感器的指向性角度偏差可表示為
(57)
(58)
則矢量傳感器陣列的軸向角度偏差矩陣可表示為
(59)
式中,軸向角度偏差矩陣初始化為
(60)
所提算法總結(jié)如下。
第三步:根據(jù)式(30),更新R(j)。
對所提TWAIA算法計(jì)算復(fù)雜度進(jìn)行分析,這里假設(shè)迭代次數(shù)為Titer(Titer≤Tmax)。所提算法初始化部分需MN+(2L+1)N次乘法與MN+(2L-1)N次加法,計(jì)算協(xié)方差矩陣R的逆需M2N+MN+21M3次乘法與M2N+21M3次加法,計(jì)算網(wǎng)格上對應(yīng)稀疏信號功率需7M2N+5MN+3N次乘法與6M2N+2MN+2N次加法,估計(jì)和重構(gòu)軸向角度偏差矩陣需12M2N+3MN+KN次乘法與4M2N+MN次加法。所提算法的計(jì)算量為O{31TiterM2N+(12Titer+2)MN+(5Titer+4L)N+TiterKN+42M3}。MUSIC不需要進(jìn)行逆矩陣和迭代運(yùn)算,計(jì)算量為O{M2L+M3N}。SAMV-1和AIAA的計(jì)算量分別為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算法具有較高的計(jì)算量,但是在比較的幾種方法中,TWAIA算法的方位估計(jì)精度是最高的。
在本章中,將所提TWAIA方法與經(jīng)典的 MUSIC方法、SAMV-1方法及AIAA方法進(jìn)行比較,評估所提方法的有效性與穩(wěn)健性,同時(shí)給出了矢量傳感器陣列的克拉美羅下界(Cramer-Rao lower bound, CRLB)曲線,以此作為比較的基準(zhǔn)。所有方法都是基于一個(gè)由M元矢量傳感器組成的均勻線列陣進(jìn)行的研究,并且相鄰陣元間的間距為一個(gè)波長。除了位于坐標(biāo)原點(diǎn)處矢量傳感器的振速軸方向作為參考之外,在每次蒙特卡洛仿真中,其它矢量傳感器的振速軸向相較于參考陣元的軸向角度偏差滿足均值為β,方差為ρ2的均勻分布。等功率遠(yuǎn)場窄帶不相干的信號為復(fù)指數(shù)信號,其中心頻率f=711 Hz,系統(tǒng)采樣頻率為fs=8 kHz,聲源的角度掃描范圍為[-180°∶2°∶180°]。在所提方法中,τ=1×10-3,Tmax=20,根據(jù)Wang等的研究和文獻(xiàn)[27]的分析,υ=45,λs=λg=0.625,q=0.5。下文仿真中把均方根誤差作為估計(jì)精度的衡量標(biāo)準(zhǔn),并定義均方根誤差(root mean square error, RMSE)為
(61)
試驗(yàn)一:均方根誤差與信噪比的關(guān)系
假設(shè)兩個(gè)等功率遠(yuǎn)場不相干的信號源入射在矢量傳感器陣列上,圖3給出了兩個(gè)信號源的方位角分別位于θ1=-14°,θ2=38°,陣元個(gè)數(shù)為4,快拍數(shù)為300,β=15,ρ=4時(shí),均方根誤差與信噪比的關(guān)系。從圖3可以看出,隨著信噪比的增加,MUSIC的均方根誤差相較于CRLB具有較大的間隔,且MUISC的均方根誤差變化趨勢較緩慢,這主要是由于MUSIC未考慮軸向角度偏差對矢量傳感器陣列方位估計(jì)的影響。盡管SAMV-1方法也未考慮軸向角度偏差對矢量傳感器陣列方位估計(jì)的影響,但是相較于MUSIC方法,SAMV-1方法與CRLB之間的間隔有所減小,這表明SAMV-1方法對矢量傳感器陣列軸向角度偏差具有較低的敏感性。相較于MUAIC方法和SAMV-1方法,盡管AIAA通過估計(jì)軸向角度偏差矩陣對觀測數(shù)據(jù)進(jìn)行了校準(zhǔn),減小了它對目標(biāo)方位角度估計(jì)性能的影響。顯著的是,AIAA方法在信噪比較高時(shí),它的估計(jì)性能接近于CRLB,但在信噪比較低時(shí),它的估計(jì)性能仍較差。這主要是由于AIAA方法在估計(jì)稀疏信號時(shí)是根據(jù)陣列的一階統(tǒng)計(jì)特性進(jìn)行的估計(jì),加之軸向角度偏差矩陣采用最小二乘估計(jì)未考慮軸向角度偏差矩陣與殘差之間的關(guān)系導(dǎo)致AIAA方法在低信噪比時(shí)方位估計(jì)精度不高。相反,所提方法的均方根誤差與CRLB之間具有較小的間隔,這主要?dú)w因于所提方法是依據(jù)矢量傳感器陣列輸出的二階統(tǒng)計(jì)量估計(jì)稀疏信號功率,并且采用正則化加權(quán)最小二乘估計(jì)軸向角度偏差矩陣,并通過戶用參數(shù)較好的平衡了信號在空域的稀疏性及軸向角度偏差矩陣與殘差之間的關(guān)系,使得所提方法在低信噪比下矢量傳感器陣列的方位估計(jì)精度得以提高。
圖3 均方根誤差與信噪比的關(guān)系Fig.3 RMSE versus SNR
試驗(yàn)二:均方根誤差與快拍數(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方法的估計(jì)精度有很大提升,但是與CRLB之間仍有較大的間隔。在快拍數(shù)小于40時(shí),AIAA方法方位估計(jì)性能惡化,這主要是由于在快拍數(shù)較小時(shí),恢復(fù)的稀疏信號不夠精確,導(dǎo)致最終估計(jì)的稀疏信號功率誤差較大,從而致使方位估計(jì)誤差較大。所提方法依據(jù)矢量傳感器陣列輸出的二階統(tǒng)計(jì)量估計(jì)稀疏信號功率,并采用正則化參數(shù)及用戶參數(shù)進(jìn)行約束有效的提高了快拍數(shù)不足時(shí)AIAA方法存在的弊端,并且當(dāng)快拍數(shù)充足時(shí)所提方法的均方根誤差曲線近似于CRLB,表明了所提方法具有較好的方位估計(jì)性能。
圖4 均方根誤差與快拍數(shù)的關(guān)系Fig.4 RMSE versus the number of snapshot
試驗(yàn)三:均方根誤差與陣元個(gè)數(shù)的關(guān)系
圖5比較了所提方法與現(xiàn)有方位估計(jì)方法的均方根誤差與矢量傳感器個(gè)數(shù)的關(guān)系,其中,信噪比為8 dB,陣元個(gè)數(shù)從2到7變化,其他仿真參數(shù)與圖3保持一致??梢钥吹剑S著陣元數(shù)的增加,MUSIC方法、SAMV-1方法、AIAA方法及所提方法的估計(jì)性能逐漸提高。MUSIC方法在陣元個(gè)數(shù)為8時(shí)無法達(dá)到無偏估計(jì)。SAMV-1方法在陣元個(gè)數(shù)為7時(shí)可實(shí)現(xiàn)無偏估計(jì)。AIAA方法及所提方法在陣元個(gè)數(shù)為5時(shí)即可實(shí)現(xiàn)無偏估計(jì)。相較于未校準(zhǔn)的MUSIC方法和SAMV-1方法,在無偏估計(jì)性能的情況下,所提方法與AIAA方法需要的陣列孔徑最小。此外,在具有相同孔徑的條件下,相較于AIAA方法所提方法具有更高的估計(jì)精度,這充分展現(xiàn)了所提方法在小孔徑時(shí)優(yōu)越的方位估計(jì)性能。
圖5 均方根誤差與陣元個(gè)數(shù)的關(guān)系Fig.5 RMSE versus the number of snapshot
試驗(yàn)四:均方根誤差與偏差均值的關(guān)系
圖6比較了所提方法與現(xiàn)有方位估計(jì)方法的均方根誤差與矢量傳感器陣列軸向角度偏差均值之間的關(guān)系,其中,信噪比為8 dB,β從5到35變化,其他仿真參數(shù)與圖3保持一致??梢钥闯?,幾種方法的均方根誤差都會隨著β的增加而增加。相較于未考慮陣列軸向角度偏差參數(shù)的MUSIC方法和SAMV-1方法,AIAA方法和所提方法的均方根誤差隨軸向角度偏差均值增大的速率較大。盡管AIAA方法和所提方法的均方根誤差隨軸向角度偏差均值增大的速率較小,但是較于AIAA方法,所提方法的均方根誤差更小,這展示出所提方法具有更好的方位估計(jì)性能。
圖6 均方根誤差與軸向角度偏差均值之間的關(guān)系Fig.6 RMSE versus the mean value of axial angle bias
為了提高低信噪比振速軸向不一致時(shí)矢量傳感器陣列的方位估計(jì)精度,提出了一種兩步加權(quán)交替迭代自適應(yīng)的矢量傳感器陣列方位估計(jì)方法。該方法采用正則化加權(quán)稀疏協(xié)方差矩陣擬合方法估計(jì)稀疏信號功率,通過對稀疏信號功率進(jìn)行補(bǔ)償提高了其在空域的稀疏性;利用加權(quán)最小二乘估計(jì)軸向角度偏差矩陣,并根據(jù)矩陣中偏差的分布特性,重構(gòu)了期望的軸向角度偏差矩陣,以此交替的方式迭代更新稀疏信號功率和軸向角度偏差矩陣,使得最終估計(jì)的稀疏信號功率更加精確,從而提高了振速軸向不一致時(shí)矢量傳感器陣列的方位估計(jì)精度。仿真結(jié)果表明,在低信噪比、少快拍、小孔徑下,相較于MUSIC方法、SAMV-1方法和AIAA方法,所提方法提高了振速軸向不一致時(shí)矢量傳感器陣列的方位估計(jì)精度。