分塊迭代導(dǎo)向最小方差算法及其在線譜目標檢測中的應(yīng)用

郭浩泉， 朱代柱,2

(1.上海船舶電子設(shè)備研究所， 上海 201108； 2.水聲對抗技術(shù)重點實驗室， 上海 201108)

高方位分辨力的波束形成在目標檢測、參數(shù)估計與識別中具有重要的意義，在多目標和強弱目標同時存在的情況下，有效抑制多目標相干干擾是后續(xù)信號處理的基礎(chǔ)。導(dǎo)向最小方差(steered minimum variance, STMV)波束形成算法以其高方位分辨力、無旁瓣以及良好穩(wěn)健性等優(yōu)越性能而受到廣泛關(guān)注，其中最主要的運算量取決于協(xié)方差矩陣的求逆過程，對陣元數(shù)為M的基陣進行經(jīng)典STMV分析，每個導(dǎo)向方位的每個頻率點都需要進行O(M3)次乘法運算[1-2]，這對大規(guī)模基陣的寬帶信號處理是難以承受的，必須尋求更低復(fù)雜度的簡化算法。

為了降低STMV算法中矩陣求逆過程的復(fù)雜度，許多學(xué)者做了大量研究，基本可以歸納為降維處理、遞歸迭代、矩陣分解和矩陣分塊4大類方向：1)在降維處理方面，文獻[3-7]主要利用盲分離的自適應(yīng)波束形成、對導(dǎo)向矢量進行處理以降低干擾的自適應(yīng)波束形成、結(jié)合穩(wěn)健Capon波束形成器和Krylov子空間、基于陣形在線估計聚焦和Krylov子空間等降維技術(shù)。2)在遞歸迭代方面，文獻[8-9]主要利用基于統(tǒng)計優(yōu)化準則的迭代MVDR算法、基于矩陣求逆引理的遞歸MVDR算法進行迭代運算減少計算量，文獻[10]提出了基于一階遞歸濾波器和矩陣求逆引理的迭代STMV算法，將計算量降低至經(jīng)典STMV算法的2/M。3)在矩陣分解方面，文獻[11] 主要利用基于QR分解的自適應(yīng)波束形成算法將求解權(quán)向量的問題轉(zhuǎn)化為求解三角線性方程的問題，避免了自相關(guān)矩陣的估計和求逆運算；文獻[12-14]提出了子陣STMV波束形成算法，大幅度降低了運算量。4)在矩陣分塊方面，文獻[15-17]主要利用協(xié)方差矩陣為正定Hermite陣的對稱性并結(jié)合Strassen矩陣求逆的高效性降低協(xié)方差矩陣求逆的復(fù)雜度。

這些研究將全陣平均分為多個子陣的矩陣分解方法計算量最小，但由于削弱了子陣之間的信號相關(guān)性，在本質(zhì)上是一種近似算法，所以與全陣STMV算法相比性能略有損失；并且，如何分配子陣數(shù)目和子陣陣元數(shù)成為一個繞不開的問題，缺乏尋求最優(yōu)解的理論依據(jù)。其他降維處理、遞歸迭代和矩陣分塊研究方面，計算量最小的是文獻[10]在保持STMV性能的同時將計算量降低至經(jīng)典STMV算法的2/M。

本文針對經(jīng)典STMV的大運算量問題，提出了一種低復(fù)雜度的BISTMV算法，通過對矩陣求逆公式的層層分解與分塊迭代，可以大幅度降低協(xié)方差矩陣求逆的復(fù)雜度，更加有利于工程應(yīng)用。

1 導(dǎo)向最小方差波束形成

1.1 STMV的基本原理

基于導(dǎo)向協(xié)方差矩陣的導(dǎo)向最小方差波束形成算法的基本原理是通過波束旋轉(zhuǎn)，使指定方向的響應(yīng)為1，而陣列輸出功率最小。

假設(shè)各陣元時域輸出為x(t)=[x1(t)…xM(t)]T，相應(yīng)的頻域輸出為X(fk)，各陣元的權(quán)向量為W=[w1w2…wM]T，則導(dǎo)向矩陣為：

(1)

式中：M為陣元數(shù)；θ為方位角；d為陣元間距；c為聲速。

陣列的輸出為：

(2)

式中：h為頻率點上限；l為頻率點下限。

陣列輸出的平均功率為：

P(θ)=E[y(t,θ)y(t,θ)H]

(3)

STMV波束形成器的最佳權(quán)向量w是使指定方向的響應(yīng)為1，而陣列輸出功率最小，故：

(4)

1.2 迭代STMV算法原理

STMV波束形成時，對N次快拍的頻域輸出Xn(fk)進行求和，得到平均頻域輸出為：

(5)

則導(dǎo)向協(xié)方差矩陣Rstmv(θ)可表示為：

(6)

(7)

(8)

式中：α為一階遞歸濾波器的系數(shù)；β=1-α。

由矩陣求逆引理：

(9)

(10)

(11)

1.3 分塊迭代STMV算法原理

1.3.1 分塊矩陣求逆算法

關(guān)于2階分塊矩陣的求逆公式，文獻[18]已詳細闡述。

(12)

(13)

把式(13)命名為“組合求逆公式”。

證明考慮到：

則

(14)

式中E-1可由式(12)求得。

(15)

1.3.2 算法原理

假設(shè)U被等分成K段，記為Ui，i=1,2，…,K，每段的長度為k=M/K。則導(dǎo)向協(xié)方差矩陣Rstmv(θ)被等分成K×K塊，每塊矩陣的階數(shù)為k?？紤]到Rstmv(θ)為埃爾米特矩陣，則式(7)可表示為：

(16)

(17)

2) 構(gòu)造矩陣

(18)

(19)

3) 應(yīng)用分塊矩陣求逆公式遞推求解E-1；

2 計算量分析

2.1 經(jīng)典STMV算法計算量

假設(shè)陣元數(shù)為M，信號帶寬為B，利用N次快拍的頻域輸出估計互譜密度矩陣，注意到乘法運算的計算時間為加法運算的數(shù)十倍(約40倍)，因此本文僅考慮乘法運算的計算量。其計算量主要集中在2部分：

1)估計互譜密度矩陣時，需在每個頻點處需進行NM2次乘法運算，總共需進行NM2B次乘法運算；

2) 估計空間譜估計時，需在每個波束進行如下計算：

①估計導(dǎo)向協(xié)方差矩陣，在每個頻點處需進行2M2次乘法運算，總共需進行2M2B次乘法運算；

②對導(dǎo)向協(xié)方差矩陣求逆需進行M3次乘法運算，然后利用式(4)進行空間譜估計，需進行M2+M次乘法運算。

綜上所述，假設(shè)對L個波束進行空間譜估計，所需的計算量為：

NM2B+[M3(2B+1)+M2+M]L

(20)

2.2 迭代STMV算法計算量

假設(shè)m階矩陣求逆的計算量為C(m)，則根據(jù)式(12)、(14)、(15)，2×2、3×3、4×4分塊矩陣求逆的計算量依次為：

x2=2C(m)+4m3

(21)

x3=x2+C(m)+11m3

(22)

x4=x3+C(m)+21m3

(23)

綜上所述，考慮n×n分塊矩陣求逆的計算量的通項公式。n×n分塊矩陣求逆的計算量為：

(24)

結(jié)合式(21)～(23)，可得通項公式為：

3(n-1))m3

(25)

假設(shè)陣元數(shù)為M，信號帶寬為B。其計算量主要集中在2部分：

綜上所述，假設(shè)對L個波束進行空間譜估計，所需的計算量為：

[M2(4B+1)+M(B+1)]L

(26)

2.3 分塊迭代STMV算法計算量

假設(shè)陣元數(shù)為M，信號帶寬為B，U被等分成K段，記為Ui，i=1,2，…,K。每段的長度為k=M/K。其計算量主要集中在5部分：

4) 利用分塊矩陣求逆公式求解E-1，需進行(K-1)C(k)+((K-2)3+4(K-2)2+3(K-2))k3/2次乘法運算；

綜上所述，假設(shè)對L個波束進行空間譜估計，所需的計算量為：

(M2B+(3k2+k)B+(K(K+1)/2-1)k2B+

(K-1)C(k))L+((K-2)3+4(K-2)2+

3(K-2))k3L/2+((K-1)2k3+(K-1)k3+

M2+M)L

(27)

式中：C(k)為k階矩陣求逆的計算量，取C(k)=k3。同時考慮到k=M/K，式(27)可化簡為：

(28)

(29)

假設(shè)信號帶寬B=5 kHz，在0°～180°方位內(nèi)每隔2°形成一個預(yù)成波束，即對L=91個波束進行空間譜估計，快拍數(shù)N=4，K=16，STMV算法與分塊迭代STMV算法計算量的比值隨陣元數(shù)M的變化如圖 1 所示。

圖1 2種算法計算量比值

文獻[10]提出的迭代算法的計算量大致是STMV算法計算量的2/M。為便于分析，同時考慮到L?N、B?M，式(20)和式(29)分別可近似化簡為：

(30)

當K足夠大時，本文提出算法的計算量大致是原STMV算法計算量的1/4M，約為文獻[10]迭代STMV算法計算量(2/M)的1/8。當M?5時，計算量顯著降低。

3 仿真分析

仿真條件如下：設(shè)48元均勻線陣接收，陣元間距d為信號上限頻率對應(yīng)的半波長，每快拍處理時間為3 s，采樣頻率為fs，處理帶寬為0.05～0.25fs，聲速為1 500 m/s。兩目標輻射噪聲信號及背景噪聲均為相互獨立的高斯白噪聲，目標1的方位為30°保持不變，陣元輸入端信噪比為6 dB；目標2的初始方位為50°，陣元輸入端信噪比為0 dB，在0.1fs處疊加一根線譜，該線譜相對于周圍連續(xù)譜的信噪比為15 dB，對總能量的貢獻可忽略不計，目標2每快拍的目標方位增加0.4°。取α=0.5，由STMV算法和分塊迭代BISTMV算法(K=6)估計的130 s空間譜歷程圖看不出明顯差異，截取不同時刻的空間譜如圖2所示，圖3為圖2(b)的局部放大圖。

圖2 不同時刻2種算法的空間譜

圖3 線譜目標局部放大圖

4 海試數(shù)據(jù)分析

為驗證算法的穩(wěn)健性，對海上試驗數(shù)據(jù)進行處理，試驗態(tài)勢如圖4所示。以拖船的航向為0°且直線航行，其后拖著一條均勻線列陣。目標1距離拖船較近，其輻射噪聲較強，目標4的初始方位為75°且朝著拖船方向航行。試驗時，航運繁忙且目標3、4相對于目標1、2距離拖船較遠，其輻射噪聲較弱。

圖4 試驗態(tài)勢

對接收信號進行預(yù)成波束形成，處理帶寬為0.025～0.1fs，由STMV算法和BISTMV算法(K=6)估計的空間譜歷程圖也看不出明顯差異，截取第88 s和141 s時刻的空間譜，如圖5所示，當存在強目標時(目標1、3)，仍可實現(xiàn)對弱目標(目標2、4)的有效檢測。

圖5 88 s和141 s 2種算法的空間譜

由圖5(a)和5(b)可以看出，BISTMV算法在同時存在多個強弱目標的情況下，對每個強弱目標的檢測效果都與STMV算法比較貼合，對其中的線譜目標還有1～4 dB(平均2.2 dB)的檢測增益,如圖6(a)和6(b)所示。究其原因，是由于BISTMV算法在式(8)對協(xié)方差矩陣進行求逆的過程中引入了一階遞歸濾波器。一階遞歸濾波器將前后快拍之間的信息進行了關(guān)聯(lián)，由于目標輻射噪聲信號是平穩(wěn)相關(guān)的，而背景噪聲是弱相關(guān)的，因此，這種關(guān)聯(lián)在大多數(shù)情況下對目標檢測是有利的：當處理寬帶信號時幾乎沒有差異，但當處理窄帶信號時效果會有改善，尤其是當目標輻射噪聲信號中包含有窄帶線譜信息時(如圖5(a)和5(b)中75°左右的目標4)；只有當目標方位快速變化時，這種關(guān)聯(lián)才會略微展寬該目標所在波束的寬度(如圖5(a)中127°左右的目標1)，造成輕微負面影響。但總體而言，BISTMV算法在繼承STMV算法的優(yōu)良特性的同時大幅減小了計算量，對窄帶線譜信息豐富的目標檢測還可以獲得一定程度的增益提升。

圖6 線譜目標局部放大圖

5 結(jié)論

1)本文提出了基于一階遞歸濾波器、矩陣求逆引理和矩陣分塊求逆的BISTMV算法，進行了計算量分析，大致是經(jīng)典STMV算法計算量的1/4M，是迭代STMV算法計算量的1/8，顯著減小了計算量，提高了其工程可實現(xiàn)性。

2)仿真和海試數(shù)據(jù)處理結(jié)果表明，該算法保持了STMV算法高方位分辨力和對微弱目標檢測的優(yōu)勢。

3)該算法對窄帶線譜信息豐富的目標檢測可以獲得平均2 dB左右的增益提升，對工程應(yīng)用，特別是線譜目標的檢測具有重要的意義。