陳筱月, 馬曉川,*, 李 璇
(1. 中國科學院聲學研究所水下航行器信息技術重點實驗室, 北京 100190;2. 中國科學院大學電子電氣與通信工程學院, 北京 100049)
來波方位(direction of arrival, DOA)估計作為陣列信號處理的研究熱點之一,在目標跟蹤,聲源定位以及無線通信等多個領域當中具有廣泛的應用[1-4]。常規波束形成(conventional beamforming, CBF)[5]是最為經典的DOA估計算法,其通過陣列加權得到指向性波束,調整加權向量使波束指向不同方位從而得到方位譜,方位譜峰值對應的方位即真實來波方位。基于無失真約束與波束輸出噪聲最小準則,Capon[6]提出了最小方差無失真方法(minimum variance distortionless response, MVDR),利用接收信號協方差矩陣自適應的調節波束加權向量,提高了估計結果的準確性。子空間類算法[7-8]利用噪聲子空間與信號子空間正交的特性,實現了超分辨。近年來,一些基于稀疏優化、解卷積以及神經網絡的DOA算法被提出。稀疏優化類算法構造過完備基矩陣將DOA估計視作優化問題,并通過對解向量進行稀疏先驗約束從而獲得稀疏解[9-11];解卷積后處理算法基于逆濾波的原理,消除點擴散函數造成的方位譜擴展,從而得到理想方位譜[12-13];神經網絡類算法將各個角度分為不同的類別,然后利用深度神經網絡等方法進行分類,從而獲得估計結果[14-15]。
上述大多方法獲得高精度DOA估計結果的前提是進行密集的方位搜索。然而,密集的方位搜索會導致計算量的升高,特別是對于二維DOA來說[16],計算量會按平方急劇膨脹。而對于實際的聲納雷達系統,由于實時性要求及硬件系統的限制,難以通過密集的方位掃描來獲得高精度的DOA結果。因此,實際工程當中一般會采用幅度插值[17]的方法得到來波方位,該方法基于CBF,但與傳統密集譜峰搜索不同的是,該方法僅在少數方位上計算波束輸出結果,然后利用不同的函數對輸出結果進行方位譜擬合,插值得到來波方位。其中,最常用的擬合函數有線性函數、二次函數、高斯函數[18]。研究者針對幅度比較測向法的誤差成因以及估計精度的提高做了大量相關工作。顧敏劍[19]使用高斯函數對波束圖進行近似,理論分析了系統噪聲,通道幅度特性不一致以及波束軸角指向對估計結果的影響。在高斯噪聲下,文獻[20]利用泰勒級數展開得到了振幅比較單脈沖雷達的DOA估計均方誤差的解析表達式,避免了使用蒙特卡羅計算均方誤差帶來的高計算量。文獻[21]中自適應地使用相鄰或交替天線進行DOA測量,有效的降低了由于噪聲與來波方位偏離波束指向所造成的誤差。文獻[22]中利用神經網絡模型對一階泰勒展開引起的非線性誤差進行補償,得到最終的DOA結果,提高了估計精度。文獻[17]中考慮波束偏離0°時,左右波束的不對稱,利用波束圖-3 dB點作為求解條件,通過非對稱函數進行方位譜擬合,得到了更好的估計結果。文獻[23]中,提出了一種基于雷達模式高斯斜率修正的優化算法,取得了比傳統算法更好的精度。而文獻[24]直接在真實方向圖基礎上進行公式推導,進而對目標角度進行求解。上述工作大多均在一維進行,針對二維情況,文獻[25]通過蜂窩網結構分布的二維多波束,使用二元二次方程進行曲面擬合,得到了較準確的測向結果。
本文首先通過對二維平面陣方位譜的推導,將二維幅度插值DOA估計解耦為兩個相互獨立的一維幅度插值DOA估計,從而避免了文獻[25]中使用多元高次方程進行曲面擬合。對于解耦后的一維幅度插值DOA估計,考慮波束指向偏離中心時左右波束不對稱,若使用對稱函數進行插值會帶來誤差,所以使用了不同的非對稱函數進行幅度插值,并且在角度求解的過程中,區別于文獻[17]中使用-3 dB波束寬度作為求解條件,本文提出以左右波束斜率比作為求解條件,提高了算法的性能表現。針對二維幅度插值中出現的角度失配問題,提出了以迭代插值的方式對估計結果進行修正,進一步地提高了角度估計精度。最后,通過仿真實驗分析了不同插值方法的誤差曲線以及角度失配,信噪比對插值精度的影響,并通過對湖試數據的處理對算法進行了驗證,本文算法能夠應用于對實時性有要求的系統中進行高精度單目標方位估計。
本文所使用的三維坐標系統[26]如圖1所示,O為三維坐標系的原點,在此坐標系統中定義空間角為Ω=(θ,φ),其中θ為Ω與平面yOz的夾角,范圍為[-90°,90°],而φ為Ω與平面xOz的夾角,范圍為[-90°,90°]。在上述定義下,θ與φ實際上無法同時取到各自的邊界值,兩者之間存在著如下的限制關系:

圖1 三維坐標系示意圖Fig.1 Diagram of the three-dimensional coordinate system
|θ|+|φ|≤90°
(1)
即超出式限制的空間角Ω是不存在的。
二維平面陣位于xOy平面上,陣元個數為M,陣元坐標為pm=[pmx,pmy,pmz]T,則陣元坐標矩陣為P=[p1,p2,…,pM],在下文中,為了仿真的計算方便,使用的為四元矩形陣,且陣元間距d等于半波長,但是實際上本文方法對于任意二維平面陣均適用。現一窄帶平面波s(t)由遠場入射至二維平面陣上,其頻率為f,波長為λ,聲速為c,入射角度為Ω0=(θ0,φ0),其入射方向的單位方向矢量v0為
(2)
則陣列接收數據x(t)可以寫為
x(t)=a(θ0,φ0)s(t)+n(t)
(3)
式中:n(t)為加性高斯白噪聲;a(θ0,φ0)為(θ0,φ0)方向上的導向矢量,其表達式為

(4)
式(3)為連續形式,實際中會經過采樣得到其離散形式如下:
x(n)=a(θ0,φ0)s(n)+n(n),n=1,2,…,N
(5)
式中:N為快拍數;x(n)∈CM×N;a(θ0,φ0)∈CM×1;s(n)∈C1×N;n(n)∈CM×N。
二維平面陣通過對采樣得到的陣元接收信號進行加權求和,即可實現波束形成,對期望方向的信號進行輸出。常規波束形成加權向量表達式如下:
(6)
波束輸出信號為
y(θ,φ,n)=ω(θ,φ)Hx(n)
(7)
進而可以計算出其功率為
(8)

為了應對上述情況,實際系統中多采用幅度插值的方法進行角度估計。其原理如圖2所示,圖中的曲面是按照0.1°的間隔對整個空間進行掃描得到的方位幅度譜,可近似認為是一個連續的曲面(這里需要注意的是,盡管根據式(1),圖中某些角度并不存在,但是為了之后插值與理解的方便,仍然計算了這些角度的波束輸出值),紅色星號所對應的位置(20°,10°)為通過波束密集掃描估計的來波方位。幅度插值法僅在少數的空間角度計算波束輸出,如圖2中的品紅色圓點對應的角度,該操作相當于對真實的方位幅度譜進行空間采樣。這里的空間采樣間隔則與陣列在水平豎直方向上的波束寬度有關。對于矩形陣來講,其在水平豎直兩個方向上的分辨率相同,所以在這兩個方向上的采樣間隔也相同,如圖3(a)所示。但對于任意陣,如2×4的五元L形陣,如圖4(a)所示,其在水平豎直兩個方向上的分辨率不同,因此在這兩個方向上的采樣間隔也不一樣,在分辨率更高水平的方向上則需要更密的采樣間隔,以免波束主瓣落在兩個采樣點之間,一般要求采樣間隔小于主瓣寬度的一半。之后利用曲面函數對最大采樣點及其周圍多個采樣點進行曲面擬合,插值得到的曲面最大值所對應的角度則為估計的來波角度。該方法避免了對二維空間進行密集掃描,有效的降低了計算量。

圖2 二維幅度插值原理示意圖Fig.2 Two-dimensional amplitude interpolation principle diagram

圖3 二維方位幅度譜及其一維切片(矩形陣)Fig.3 Two-dimensional azimuth amplitude spectrum and its one-dimensional section (rectangular array)

圖4 二維方位幅度譜及其一維切片(2×4 L形陣)Fig.4 Two-dimensional azimuth amplitude wpectrum and its one-dimensional section (2×4 L-shaped array)
然而,進行二維曲面擬合會面臨復雜的多元高次方程組[25],求解較為困難,且面對任意陣型時,并不一定能夠找到合適的擬合函數,從而會導致估計結果誤差增加。因此,本文提出將二維幅度插值解耦為兩個不同方向(θ方向,φ方向)的一維幅度插值,下面對該方式的可行性進行驗證。對于式(8),在忽略噪聲項后得到的幅度表達式如下:
(9)
進一步將式(6)代入式(9)可以得到
(10)

然而,θ0與φ0是所需要估計的量,實際上其具體值事先未知,僅知道在方位幅度譜上的等間隔采樣值(即圖3(a)與圖4(a)中的黑色下三角點與品紅色圓點處),而其中黑色下三角點的交叉處則是這些采樣點中的最大幅度所對應的角度,記作(θm,φm)。因此,將二維曲面插值分解為(θ,φm)與(θm,φ)兩個切面(黑色下三角點指示)上的一維曲線插值。圖3(b)、圖3(c)與圖4(b)、圖4(c)分別給出了四元矩形陣與2×4的五元L形陣的方位幅度譜及其在(θ0,φ0)的兩個方向上(紅色上三角)與(θm,φm)的兩個方向上(黑色下三角)的一維切片。可以看出,對于矩形陣而言,其方位幅度譜在(θ0,φ0)的兩個方向上與(θm,φm)的兩個方向上的一維切片雖然存在幅度上的差異,但是最大值均所對應的角度相同。并且對于矩形陣來說,在任意角度的兩個方向上的一維切片均滿足此條性質。然而,對于其他任意二維平面陣,該條性質并不一定滿足,如圖4所示,L形陣方位幅度譜在(θ0,φ0)的兩個方向上與(θm,φm)的兩個方向上的一維切片極值點并不相同,若利用(θm,φm)的兩個方向上采樣點進行插值,則會產生誤差,在本文中這種誤差被稱作角度失配誤差。為了使算法能適用于任意二維平面陣,后續章節提出了使用多次迭代插值的方法來降低該項誤差。


圖5 φ=φm的一維切面Fig.5 One-dimensional section of φ=φm
3.1.1 對稱線性插值
(11)

(12)

(13)
3.1.2 對稱二次插值
二次插值則是用二次函數對方位幅度譜進行擬合,二次函數的表達式如下:
(14)

(15)
3.1.3 對稱高斯插值
高斯插值則是用高斯函數對方位幅度譜進行擬合,其表達式如下:
(16)
(17)
上述方法均認為對真實二維方位幅度譜進行切片得到的一維方位幅度譜是關于其極大值左右對稱的。然而,事實上,這一點僅當一維方位幅度譜的極大值位于0°的時候滿足,當其極大值偏離0°的時候,左右兩側將不再對稱。如圖6所示,為真實來波位于不同方位時的一維方位幅度譜,可以看出,隨著來波角度偏離0°越嚴重,方位幅度譜的左右不對稱性就越來越嚴重。因此,如果用對稱函數進行插值擬合,會給估計結果帶來較大的誤差,為了解決這個問題,下面提出使用非對稱函數進行插值擬合,以提高估計精度。

圖6 不同來波方位下的方位幅度譜Fig.6 Azimuth amplitude spectrum in different wave directions
3.2.1 非對稱線性插值
非對稱線性插值認為一維方位幅度譜在其極大值左右兩側的斜率不同,因此線性擬合函數如下:
(18)
將α、β、γ與相應的波束幅度響應Uα、Uβ、Uγ(Uα≥Uγ)代入式(18)可以得到方程組:
(19)
(20)

(21)

(22)

(23)
3.2.2 非對稱二次插值
非對稱二次插值使用兩個不同的二次函數對一維方位幅度譜極值點左右兩側進行擬合,兩個二次函數僅在二次項系數上存在差別,其插值函數如下:
(24)
將α、β、γ與相應的波束幅度響應Uα、Uβ、Uγ(Uα≥Uγ)代入式(24)可以得到方程組:
(25)

(26)
A1=ηq(β)-1-M1ηq(β)+M1
A2=1-ηq(β)+M2ηq(β)-M2
B1=2β-2ηq(β)α+2M1ηq(β)α-2M1γ
B2=2ηq(β)α-2γ-2M2ηq(β)β+2M2γ
C1=ηq(β)α2-β2-M1ηq(β)α2+M1γ2
C2=γ2-ηq(β)α2+M2ηq(β)β2-M2γ2

(27)

(28)
3.2.3 非對稱高斯插值
非對稱高斯插值使用兩個不同的高斯函數對一維方位幅度譜極值點左右兩側進行擬合,使用的非對稱高斯插值函數如下:
(29)
將α、β、γ與相應的波束幅度響應Uα、Uβ、Uγ(Uα≥Uγ)代入式(29)可以得到方程組:
(30)

(31)
M3=[ln(Uα/Uβ)]/[ln(Uα/Uγ)]
M4=ln(Uγ/Uα)/ln(Uγ/Uβ)
A3=-ηg(β)+1+M3ηg(β)-M3
A4=-1+ηg(β)-M4ηg(β)+M4
B3=-2β+2ηg(β)α-2M3ηg(β)α+2M3γ
B4=-2ηg(β)α+2γ+2M4ηg(β)β-2M4γ
C3=-ηg(β)α2+β2+M3ηg(β)α2-M3γ2
C4=-γ2+ηg(β)α2-M4ηg(β)β2+M4γ2

(32)

(33)

(34)
下面將對二維幅度插值方法的估計誤差進行討論,并分析不同因素對算法估計結果精度的影響。
仿真使用的陣列為四元矩形陣并按照半波長布陣,來波信噪比為20 dB。在與z軸夾角小于50°的區域當中,按照10°的間隔進行幅度采樣,使用上面提出的6種插值方法分別進行角度估計,并通過100次蒙特卡羅實驗得到該區域內的二維誤差曲面。圖7給出了非對稱高斯插值的估計誤差曲面,為了方便對不同方法的估計精度進行觀察與對比,之后不依次展示各個方法的估計誤差曲面,而是展示各方法二維誤差曲面在φ=0°處的切面,如圖8所示。

圖7 非對稱高斯插值二維誤差曲面Fig.7 Two-dimensional error surface of asymmetric Gaussian interpolation

圖8 不同插值方法方位估計誤差曲線(φ=0°切面)Fig.8 Azimuth estimation error curves of different interpolation methods (φ=0° section)
可以看出,隨著入射角度逐漸偏離0°,估計誤差呈現振蕩上升。對稱插值方法中,由于線性函數與方位譜的相似程度較差,其估計結果誤差最大,而二次函數,高斯函數與方位譜具有較好的相似度,其擬合結果更加準確。而非對稱插值相對對稱插值,線性函數的整體估計誤差并沒有明顯的改變,只是其誤差曲線的極值點出現了移動,二次函數與高斯函數的誤差出現了較為明顯的下降,得到了更加精確的估計結果。進一步可以發現,非對稱二次插值與非對稱高斯插值的誤差曲線的極小值與幅度采樣角度一致,而誤差曲線的極大值出現在兩個相鄰采樣角度中間,即當入射信號方位在幅度采樣角度附近的時候可以得到更小的估計誤差。

由圖8(b)可以看出,非對稱二次插值與非對稱高斯插值的誤差曲線極大值點出現在兩個相鄰的幅度采樣角度中點處,極小值點則剛好出現在幅度采樣角度上。這是因為當真實角度剛好位于兩個相鄰幅度采樣點中點的時候,真實角度離兩側的幅度采樣點距離最遠,所產生的角度失配最嚴重,誤差較高;而入射角度剛好位于幅度采樣點處時,幾乎沒有角度失配產生,所以誤差很小。于是本文提出多次迭代插值,在每次迭代插值中,以上一次的插值結果為基礎建立新的采樣網格,這樣,真實來波方位距離新的采樣網格點更近,因此產生的失配誤差就越小,下面分別對矩形陣與任意二維陣的迭代插值方式進行介紹。


經過多次迭代插值,即可以一步步地降低角度失配誤差,使估計結果逼近真實值。下面通過仿真對提出的方法進行驗證,仿真條件與第4.1節所使用的仿真條件相同,圖9給出了非對稱插值方法經過一次迭代插值后得到的誤差曲線,*代表迭代插值結果。

圖9 單次非對稱插值誤差估計曲線與一次迭代非對稱插值誤差估計曲線對比Fig.9 Comparison of the error estimation curve of single asymmetric interpolation with that of one-iteration asymmetric interpolation
由仿真可以看出,通過以第一次插值結果為基準,重新進行幅度采樣并插值得到的估計誤差有了明顯的降低。實際上,對于矩形陣而言,只需要進行一次迭代插值就能夠得到理想的估計結果,而對于任意二維陣,其角度失配更加嚴重,需要經過多次迭代插值才能夠得到理想的估計結果。

從圖10可以發現,常規波束形成的方位譜僅在高信噪比的情況下與波束圖有較高的一致性,而隨著信噪比的逐漸降低,兩者之間的一致性也出現了明顯的下降。在第4.2節進行非對稱幅度插值的過程當中,為了對方程組進行求解,文獻[17]引入了波束圖的-3 dB點作為求解條件,本文則是定義了波束圖的左右波束斜率比作為求解條件。而幅度插值實際上是針對方位譜進行插值,在低信噪比下,由于真實方位譜會相對波束圖出現展寬,此時使用基于波束圖的額外求解條件進行求解會導致最終的插值結果出現誤差。盡管如此,由于無法事先知道方位譜,所以只能使用基于波束圖的求解條件。這里想要說明的是,本文所使用的求解條件波束圖斜率比對信噪比的變化并不敏感,其相對于-3 dB波束寬度更加穩定。

圖10 不同信噪比下二維方位譜切片(φ=0°)Fig.10 Section of two-dimensional azimuth amplitude spectrum at different signal to noise ratios (φ=0°)
圖11給出了不同來波方位下,方位譜的-3 dB寬度隨信噪比的變化,圖12給出了不同來波方位下,方位譜的斜率比隨信混比的變化曲線。

圖11 不同DOA下方位譜的-3 dB寬度隨信噪比的變化曲線Fig.11 Variation curve of -3 dB width of azimuth spectrum with signal to noise ratio under different wave directions

圖12 方位幅度譜的斜率比隨信噪比的變化曲線Fig.12 Curve of slope ratio of the azimuth spectrum varies with the signal to noise ratio
從圖11可以明顯看出,方位譜的-3 dB寬度隨著信噪比的降低逐漸變大,甚至當信噪比低到某個程度之后,方位譜的-3 dB寬度已經不存在了,并且對于小孔徑的陣列來說,當來波方位偏離0°之后,其方位譜的一側也不存在-3 dB寬度。因此,使用波束圖的-3 dB點作為求解條件進行非對稱插值,其估計結果受信噪比的影響較為嚴重,雖然文獻[17]中通過粗估信噪比,并根據信噪比對插值過程中左右波束寬度進行修正,但這樣仍然會存在較大的誤差。而本文定義的斜率比隨信噪比的降低變化較小,相對穩定,這代表雖然低信噪比下方位譜相對于波束圖有所展寬,但是在這個過程中斜率比幾乎不改變,因此使用波束圖的斜率比進行非對稱插值,能夠在低信噪比下獲得較高的插值精度。
接下來以10°作為采樣間隔,在不同的信噪比下計算了不同插值方法在與圖1中z軸之間空間夾角小于50°范圍內的平均估計誤差,經過100次蒙特卡羅計算得到的誤差曲線如圖13所示。


圖13 不同插值方法估計誤差隨信噪比的變化曲線Fig.13 Estimation error curves of different interpolation methods varies with the signal to noise ratio
在圖13(a)中,線性插值由于其與真實方位譜擬合較差,所以具有較高的估計誤差;而二次函數插值與高斯插值方法的性能相近,在各個信噪比下均低于線性插值的誤差,且在考慮真實方位譜的非對稱性后其誤差有明顯的下降。同時,對比圖13(b)可以發現,在經過一次迭代插值之后,對稱線性插值誤差存在略微下降,而非對稱線性插值誤差下降更加明顯;對稱高斯插值與對稱二次插值誤差反而有所升高,這是因為圖8(a)中這兩種方式的誤差曲線的極小值點并不是在采樣點附近,而是在兩個采樣點之間;而非對稱二次插值與非對稱高斯插值在低信噪比下誤差相對單次插值改變不明顯,但是在較高信噪比下(10 dB以上),出現了下降,取得了上述所有方法中最優秀的估計效果。
為了進一步對本文提出的非對稱幅度插值方位估計算法的有效性與估計精度進行驗證,下面將對實際的湖試采集數據進行處理。數據于杭州千島湖采集,千島湖是一個人造湖,水底存在一定起伏,水面較為平靜。實驗區域湖深約80 m,無水底山體干擾,實驗時接收陣列與目標均布放在水深20 m附近。其中,接收陣列為均勻半波長布陣四元矩形陣,陣列朝向斜向下,由于一些其他原因,四元矩形陣的右下角陣元沒有采集到有效數據,因此實際上使用的陣列為三元的L形陣列。目標為一人工聲源,發射信號為10 kHz的連續波(continuous wave, CW)信號。
對接收信號按照0.01°的間隔在整個空間進行密集波束形成掃描,得到二維方位譜如圖14所示,可以發現峰值出現在 (-17.85°,43.51°)的位置,該處即為密集掃描估計得到的來波方位。而在二維方位譜的左下角也存在峰值,需要說明的是,根據圖2坐標系的定義,該峰值所對應角度|θ|+|φ|>90°,實際上這個空間角并不存在,因此將其排除。以二維方位譜最大值所對應的空間角作為基準,采樣間隔為10°,在表1中給出了不同的幅度插值方法的估計結果進行對比。表1中,密集波束掃描的掃描間隔為0.01°,*代表進行了4次迭代插值。

表1 不同插值方法的估計結果Table 1 Estimated results of different interpolation methods (°)

圖14 密集掃描方位幅度譜Fig.14 Azimuth amplitude spectrum with intensive scanning
從表1的插值結果可以看出,在只進行單次插值的結果中,非對稱插值的估計結果雖然略微優于對稱插值,但是整體的估計結果與密集波束掃描得到的結果相差較大。這是因為采集數據使用的陣列為三元L形陣列而非矩形陣,而對于任意陣來說,在幅度采樣最大值處將其二維幅度插值分解為兩個方向上獨立的一維幅度插值會因為角度失配產生誤差,這一點已經在第3節中進行了詳細論述。為了降低由角度失配造成的誤差,使用了在第4.2節中提出的多次插值的方法,經過4次迭代插值后,可以發現各種插值方法的估計精度有了較為明顯的提升,其中非對稱二次插值與非對稱高斯插值的估計精度提升最為明顯,其估計結果已經相當的接近密集波束掃描的結果。而對稱插值法在經過多次插值后,其估計結果精度依然較差,事實上繼續增加迭代插值次數,對稱插值算法的估計精度也不會有較好的提升,這是由于這類插值方法沒有充分的考慮方位譜的左右不對稱性所導致。另外需要說明的是,雖然多次插值會給計算量帶來一定的提升,但是每次迭代插值只需要多計算5次波束輸出,而這所需的計算量是相當小的,遠小于以0.01°間隔進行空間掃描的計算量。
本文提出了一種任意陣列的二維非對稱幅度插值方位估計算法。算法將二維幅度插值解耦為兩個獨立的一維幅度插值,然后利用左右波束斜率比作為求解條件,通過非對稱函數對少量離散波束輸出進行幅度插值得到來波角度。進一步的,針對插值過程中的角度失配誤差提出迭代插值的方法進行修正。通過仿真,本文提出的算法相對其他插值方法估計誤差更小,湖試數據的處理結果說明,本文方法能夠以較小的計算量取得與密集掃描常規波束形成類似的估計精度。本文算法具有如下優點:一是采用插值的方法進行方位估計,計算量很小;二是對陣型沒有特殊要求,且能夠在陣列孔徑較小的情況下取得精確的估計結果。因此,本文算法可以應用于一些對實時性有要求的小型航行器上,實現對單目標的高精度實時方位估計。但由于本文算法是基于常規波束形成進行插值,目前只局限于對單目標的估計,如何將算法擴展至多目標,是下一步的研究目標。