一種用于水下小尺度運動陣列的目標波達方向估計方法

郭拓， 王英民， 張立琛

(西北工業大學 航海學院， 陜西 西安 710072)

一種用于水下小尺度運動陣列的目標波達方向估計方法

郭拓， 王英民， 張立琛

(西北工業大學 航海學院， 陜西 西安 710072)

為解決水下小尺度運動陣列對目標定位的孔徑和樣本不足問題，采用合成孔徑方法實現孔徑的擴展。通過分析相干信源情況下不同陣元數與樣本數之比時樣本協方差矩陣的譜分離特性，提出采用小快拍的主特征空間目標波達方向估計方法。研究結果表明：仿真實驗中，陣元數與樣本數之比為5時，該方法依然可以正確分辨4個目標；水池試驗中，陣元數與樣本數之比為4.8時，該方法同樣可以準確分辨3個相鄰目標。該方法不需要先驗信息如信源個數等，解決了小尺度運動陣列在實際應用時所暴露的孔徑和樣本方面的固有局限，可滿足實際應用需求。

聲學； 小尺度運動陣； 主特征空間； 波達方向估計； 小快拍

Abstract: The underwater small-scale moving array is limited by its small aperture and sample size for direction of arrival (DOA) estimation of coherent target. The passive synthetic aperture technique is used to solve the problem of insufficient aperture. The spectral separation of sample covariance matrix is stu-died in the case of coherent target. A DOA estimation method for coherent target based on main feature space is proposed, which uses small snapshots. In the simulation, the proposed method can still distinguish the four targets correctly when the ratio of the number of sensors to the number of samples is 5; in the water tank experiment, the above ratio is 4.8, and it can identify three adjacent targets clearly. The proposed method can achieve better resolution for coherent targets in small samples, which meets the application needs of small-scale motions array of targeting, and requires no priori information of signal source number.

Key words: acoustics; small-scale moving array; main feature space; DOA estimation; small snapshot

0 引言

小尺度聲納基陣由于其體積小、成本低、基陣校準和使用方便，廣泛應用于魚雷、水雷、水下無人航行器等小尺度安裝平臺上。但該聲納平臺的運動特性導致在信號的局部平穩過程中所獲得的樣本數較小，小尺度陣即陣列孔徑較小，因此基于小尺度運動陣列的目標波達方向(DOA)估計成為一個關鍵問題。

DOA估計作為陣列信號處理的一個重要分支，已經廣泛用于雷達、聲納、通信、無源定位等多個領域[1-3]。陣列DOA估計經過幾十年的發展，已經形成了兩類技術，即基于波束形成技術和基于子空間分解技術。其中：基于波束形成技術的典型方法有常規波束形成器(CBF)和最小方差無畸變響應(MVDR)波束形成器；基于子空間分解技術的經典方法[4]是多重信號分類(MUSIC)算法。

上述傳統DOA估計算法大都基于孔徑較大的陣列，其估計性能隨著陣列孔徑的減小會迅速下降，且都需要估計基陣的統計協方差矩陣。在實際數據處理中都以樣本協方差矩陣代替統計協方差矩陣，在小樣本情況下如此估計出的統計協方差矩陣誤差很大，會導致DOA估計算法失效。特別是對于相干目標，經典的高分辨方法如MUSIC隨著信號相干程度的增加，空間分別率逐漸下降，甚至完全失效。因此針對水下小尺度運動陣列對相干目標DOA估計成為亟待解決的問題，以實現在小孔徑情況下采用小樣本進行方位估計。

針對小樣本數問題，為了獲得與統計協方差矩陣更加相近的樣本協方差矩陣，通常情況下會對樣本協方差矩陣進行修正，目前已經發展了很多修正方法，如：1)對角加載方法[5]、收縮因子方法[6-8]，二者都是通過對樣本協方差矩陣的特征值進行加權操作以改進其性能；2)前后向平均[9]和空間光滑技術[10]，通過移動窗口虛擬出更多的快拍，以改進樣本協方差矩陣估計；3)偽隨機噪聲再采樣技術[11-13]，通過人為生成偽隨機噪聲信號加在接收到的信號上，以改進樣本協方差矩陣估計；4)子空間泄露分析方法[14]，通過MUSIC算法的估計改進所需的樣本協方差矩陣，以減少樣本協方差矩陣方位估計時信號子空間的泄漏。上述方法的思路都是設法改進樣本協方差矩陣，以使其與大樣本時的樣本協方差矩陣更加相近，但沒有考慮樣本協方差矩陣與樣本數的關系，以及樣本數不同時協方差矩陣本身所具有的不同特性。

本文首先采用被動合成孔徑技術進行陣列的孔徑擴展，解決陣列的孔徑不足問題；對于小樣本問題將打破上述通過改進樣本協方差矩陣來改進方位估計的思路，而是通過研究不同快拍時樣本協方差矩陣本身具有的一些特性如譜的分離特性，提出基于小樣本協方差矩陣的主特征空間DOA估計方法。

1 陣列接收信號模型

如圖1所示，水平線陣由N個陣元組成，假設有D個相互獨立的平面波入射，陣列第k次快拍X(k)為

X(k)=A(θ)S(k)+N(k),

(1)

式中：X(k)=[x1(k),x2(k),…,xN(k)]T；A(θ)為陣列流形矩陣，A(θ)=[a(θ1),a(θ2),…,a(θD)]，θj為第j個平面波的水平入射角，a(θj)=[1,ej2πdsin (θj)/λ,…,ej2π(N-1)dsin (θj)/λ]為第j個入射平面波的導向向量,d為陣元間距，其值等于波長λ的一半 ；S(k)=[s1(k),s2(k),…,sD(k)]T,sj(k)為第j個平面波信號的第k次快拍；N(k)=[n1(k),n2(k),…,nN(k)]為陣列接收到的加性高斯白噪聲，ni(k)為第i個陣元接收到的噪聲。

圖1 陣列接收信號模型Fig.1 Array signal model

當信號被假設具有零均值時，陣列接收信號的協方差矩陣R=E[X(k)XH(k)]，即協方差矩陣需要通過求集總平均得到。而實際應用中只能得到有限的采樣數據，如有M次觀測樣本X(1),X(2),…,X(M)，則用這M個樣本估計的采樣協方差矩陣為

(2)

2 基于主特征空間的DOA估計方法

2.1 樣本協方差矩陣特征的譜分離特性

令y=N/M，y表示陣元數與樣本數之比，文獻[15-17]采用隨機矩陣理論研究樣本協方差矩陣譜的極限分布，研究表明樣本協方差矩陣具有譜分離特性。文獻[18]詳細研究了樣本協方差矩陣的譜分離特性，并給出樣本協方差矩陣的譜隨y變化的特點，研究了y為1.00(樣本數等于陣元數)、y為0.01(樣本數是陣元數的100倍)以及在y由1.00逐漸變小至0.01的過程中的特征值分布表現形式，研究發現：特征值分布表現為一簇一簇的，每簇中有幾個大小相近的特征值，其中每個簇的中心漸近于統計協方差矩陣R的一個真實特征值。每個簇的大小取決于y，即估計R的樣本數，如果每個陣元所取的樣本數小，則每個簇將會很大，的許多特征值將匯集在一個簇中，無法分離；反之，如果每個陣元所取的樣本數足夠大，則每個簇只對應R的一個真實特征值，這樣R的譜就得到很好的分離。

仿真實驗采用30個陣元，4個相干目標，噪聲為高斯白噪聲，可得如圖2所示的特征值隨不同y的變化曲線。由圖2可知：當y很小時，即樣本數遠遠大于陣元數即y很小時，的特征值會在最大特征值序數處陡峭升高，即最大特征值遠大于其他特征值；當y較大時，特征值不再陡峭升高，緊鄰最大特征值處還有較大的特征值。圖2中特征值隨y變化的特征，是相干源情況下樣本協方差譜分離特性的具體體現。從子空間的角度來解釋上述相干源時的譜分離特性，當樣本數遠遠大于陣元數即y很小時，樣本協方差矩陣得到較好的估計，相干信號的信息全部集中在一個最大特征值所對應的子空間上；當y較大時，樣本協方差矩陣沒有得到很好的估計，此時信號子空間發生了泄漏，使得除了最大特征值之外還出現一些連續的較大特征值，但是總體而言所有相干信號的主要信息仍然集中于最大特征值所對應的子空間中。本文將這個最大特征值所對應的子空間稱為主特征空間，而不是MUSIC算法所假設的那樣，即信號子空間就是s個信源所對應的s個較大特征值所構成的空間。通過上述對相干源情況下不同樣本數時譜分離特性的仿真實驗與分析，顯然信號信息主要集中在主特征空間而不是MUSIC算法所使用的s個較大特征值構成的空間，從而很好地解釋了MUSIC算法對相干目標DOA估計時失效的原因。

圖2 相干目標的譜分離性Fig.2 Spectral separation of coherent target

綜上所述，無論快拍數大小，目標為相干信號時其樣本協方差矩陣具有很好的譜分離特性，下面將基于譜分離特性提出相干信號的目標DOA估計方法。

2.2 基于被動合成孔徑的主特征空間DOA估計方法

被動合成孔徑聲納有幾種經典算法，本文采用較實用的擴展拖曳陣測量(ETAM)方法[19]。

假設目標是窄帶聲源，f0為聲源頻率，則目標信號可表示為s(tg)=Aexp (j2πf0tg)，A為信號幅度，tg為第g個時刻點。tg時刻第i個陣元接收到θj的信號可表示為

(3)

式中：f為目標頻率；c為水下聲速；εi,g為噪聲。

由于陣列與目標之間的相對運動，接收到的目標頻率f為

f=f0(1±vsinθj/c)，

(4)

式中：v為陣列與目標之間的相對運動速度。

將(4)式代入(3)式，得

(5)

將(5)式化簡，略去sinθ的2階項，整理得

(6)

如果陣列不移動，則tg+τ時刻(τ為時間間隔)第i個陣元接收到的信號可表示為

(7)

如果選擇合適的間隔時間τ，使得vτ=qd，即運動的距離正好等于q個陣元間距d，則tg+τ時刻第i個陣元接收到的信號變換為

(8)

而tg時刻第i+q個陣元接收到的信號為

(9)

比較(8)式與(9)式發現，二式只相差一個相位補償項，說明可以通過第i個陣元的運動來虛擬出tg時刻的第i+q個陣元。當陣列勻速運動時，時刻tg的第N-q+p個陣元與時刻tg+τ的第p個陣元在空間中同一個位置，即不同的陣元在不同的時刻接收空間相同位置上的信號。N-q為兩個時刻重疊陣元數，利用N-q個重疊陣元計算兩個時刻的相位補償，q即為一次合成孔徑虛擬出的陣元數，N為陣列物理陣元個數，經過u次合成則可以虛擬出q×u個陣元，經過u次孔徑合成的總陣元數為N+q×u.

經過虛擬孔徑合成后，物理陣和虛擬陣列接收到的信號可表示為

Xe(k)=
[x1(k),x2(k),…,xN(k),xN+1(k)，…，xN+q×u(k)]T.

(10)

經ETAM算法進行孔徑擴展后的實際使用樣本協方差矩陣可以用(11)式估計：

(11)

(12)

式中：λi為e的特征值；ei為該特征值所對應的特征向量；E為以所有特征向量按列排所構成的特征矩陣。

2.1節已將樣本協方差矩陣的最大特征值所對應的空間定義為主特征空間，此處定義只包括最大特征值所對應的空間，并非主特征空間只包含最大特征值所對應的空間，也可以根據信源的特點擴展主特征空間。定義e可以劃分為

(13)

式中：λm為主特征空間的特征值；em為與其對應的特征向量；Em為主特征空間特征向量構成的矩陣；Eo為特征空間E中去除主特征后的特征信息矩陣。

根據子空間理論，信號子空間與噪聲子空間正交，且導向向量張成的空間與信號子空間張成的空間是同一空間，即Ena(θ)=0，Esa(θ)=1，其中En為噪聲子空間，Es為信號子空間。

主特征空間中包含的主要是信號信息，即主特征子空間Em中主要包含的是Es空間的信息，當樣本數較小時會發生子空間泄漏，有部分En的信息進入Em，使Ema(θ)=1，則arccos(Ema(θ))=0，故定義其倒數為主特征空間方法的方位譜估計，即

(14)

下面將合成孔徑后，采用主特征空間方法進行方位估計的方法簡稱為ETAM-MFS.

3 仿真實驗

通過數值仿真模擬水下小尺度運動陣列目標定位，首先采用被動合成孔徑算法ETAM進行孔徑擴展，然后采用主特征空間DOA估計方法估計目標方位。仿真采用一個半波長均勻分布的6元線陣模擬小尺度陣列，陣列運動速度為0.07 m/s，入射的信號為4個遠場窄帶信號，頻率為6.25 kHz，采樣頻率50 kHz，4個信號的入射角分別為-20°、-10°、0°、10°，信噪比皆為0 dB，噪聲為復高斯白噪聲，快拍數取6，ETAM算法經過8次合成得到24個虛擬陣元，總線陣為30個陣元，陣元數與樣本數之比y為5.

圖3為CBF算法給出的目標DOA估計，圖4為MVDR算法給出的目標DOA估計，圖5為MUSIC算法給出的目標DOA估計，圖6為ETAM-MFS方法給出的目標DOA估計。由圖3～圖6可以發現，在30個陣元、6快拍情形下，CBF方法、MVDR方法與MUSIC方法都不能正確估計出目標方位，而ETAM-MFS方法可以正確分辨出目標方位。由此可見，本文提出的方法可以在小快拍情況下很好地分辨相干目標，且不需要像MUSIC算法那樣預先擁有目標個數的先驗信息，這是相干源時樣本協方差矩陣的譜分離特性帶來的好處，即信號信息全部集中于主特征空間，只有很少的信號泄漏于其他空間。

圖3 CBF方法DOA估計Fig.3 DOA estimation of CBF

圖4 MVDR方法DOA估計Fig.4 DOA estimation of MVDR

圖5 MUSIC方法DOA估計Fig.5 DOA estimation of MUSIC

圖6 ETMA-MFS方法DOA估計Fig.6 DOA estimation of ETMA-MFS

4 水池試驗

為了驗證所提算法在小尺度運動陣列小樣本的情況下能夠分辨相干目標，在消聲水池(長20 m, 寬8 m，深7 m，6面消聲)進行了相關試驗。試驗采用由6個水聽器組成的均勻線陣，設計頻率為6.25 kHz，水平布放于水下2.5 m深處，陣元間距為半波長0.12 m，陣列勻速運動速度為0.07 m/s，使用電機勻速拉動；3個發射源發射頻率均為6.25 kHz的連續波信號，與接收陣在同一深度，3個源按遠場條件布放，與陣列形成的角度分別為24°、33°、44° ；采樣頻率為50 kHz. 首先采用被動合成孔徑算法ETAM合成6次，共得到24個陣元，快拍數取5，陣元數與樣本數之比y為4.8. 試驗系統如圖7所示，試驗現場如圖8所示。

圖7 水池試驗系統Fig.7 Water tank experimental system

圖8 水池試驗現場Fig.8 Water tank experimental site

圖9～圖12為采用被動合成孔徑ETAM算法得到24個虛擬陣元后，分別使用CBF方法、MVDR方法、MUSIC方法與MFS方法得到的目標DOA估計，圖中目標的真實方位用黑點標出。由圖9可知，CBF方法不能清晰分辨3個目標，且具有明顯的方位誤差，其中物理陣所對應的是實際物理6個陣元沒有進行合成孔徑得到的目標DOA估計。由圖10可知，MVDR方法只能估計出3個目標中的2個。由圖11可知，MUSIC方法也只能估計出3個目標中的2個。從圖12可以發現，ETAM-MFS算法可以很好地分辨3個目標，且沒有方位誤差。與其他3種方法相比，首先，本文所提出的方法不需要做協方差矩陣求逆運算，減少了計算量；其次，不需要先驗信息如信源個數；最后，最重要的是，本文所提出的方法在小樣本情形下表現優異，能夠準確分辨目標。

圖9 水池試驗合成孔徑前后DOA估計對比Fig.9 DOA estimations before and after synthetic aperture in water tank experiment

圖10 水池試驗MVDR方法DOA估計Fig.10 DOA estimation of MVDR in water tank experiment

圖11 水池試驗MUSIC算法DOA估計Fig.11 DOA estimation of MUSIC in water tank experiment

圖12 水池試驗ETAM-MFS方法DOA估計Fig.12 DOA estimation of ETAM-MFS in water tank experiment

5 結論

為解決水下小尺度運動陣列的目標DOA估計中的孔徑和樣本不足問題，本文首先采用被動合成孔徑算法ETAM彌補小尺度陣的孔徑不足，然后在對樣本協方差矩陣譜分離特性分析的基礎上提出了小快拍主特征空間DOA估計方法。在分析相干目標樣本協方差矩陣的譜分離特性時，還得出了MUSIC算法對相干目標分辨失效的原因。通過仿真實驗和水池試驗，將本文所提出的方法與傳統DOA估計方法(CBF算法、MVDR算法以及MUSIC算法)進行了對比，結果發現：當快拍數很小時，試驗中陣元數與快拍之比接近5時，該方法依然可以準確地分辨多個目標，而其他3種方法或是出現多個偽目標或是有目標不能估計出。

本文所提出的方法實現了水下小尺度運動陣列的目標方位估計，主要體現在：1)通過合成孔徑克服了小尺度陣列孔徑的不足；2)能夠實現小快拍相干目標的DOA估計，可以很好地彌補MUSIC算法的缺陷，且不需要先驗信息如信源個數。該方法適合于水下小尺度運動陣列目標定位的應用需求，期望能在無人水下航行器、魚雷等應用中有很好的應用前景。

References)

[1] Krim H, Viberg M. Two decades of array signal processing research: the parametric approach[J]. IEEE Signal Processing Magazine, 1996, 13(4):67-94.

[2] 吳振, 戴繼生, 朱湘臨, 等. 互耦效應下一種基于實值稀疏表示的波達方向估計算法[J]. 兵工學報, 2015, 36(2):294-298. WU Zhen,DAI Ji-sheng, ZHU Xiang-lin, et al. A real-valued sparse representation method for DOA estimation with unknown mutual coupling[J]. Acta Armamentarii, 2015, 36(2):294-298.(in Chinese)

[3] 鄒吉武, 孫大軍. 線陣雙基地聲納波束零點形成MUSIC算法[J]. 兵工學報, 2010, 31(3):364-368. ZOU Ji-wu, SUN Da-jun. MUSIC algorithm of beam null forming on linear array of bi-static sonar[J]. Acta Armamentarii, 2010, 31(3):364-368.(in Chinese)

[4] 王永良, 陳輝, 彭應寧, 等. 空間譜估計理論與算法[M]. 北京:清華大學出版社, 2004:2-5. WANG Yong-liang, CHEN Hui, PENG Ying-ning, et al. Theory and algorithm of spatial spectrum estimation[M]. Beijing: Tsinghua University Press, 2004:2-5.(in Chinese)

[5] Carlson B D. Covariance matrix estimation errors and diagonal loading in adaptive arrays[J]. IEEE Transactions on Aerospace Electronic Systems, 1988, 24(4):397-401.

[6] Chen Y, Wiesel A, Eldar Y C, et al. Shrinkage algorithms for MMSE covariance estimation[J]. IEEE Transactions on Signal Processing, 2010, 58(10):5016-5029.

[7] Tong J, Schreier P J, Guo Q, et al. Shrinkage of covariance matrices for linear signal estimation using cross-validation[J]. IEEE Transactions on Signal Processing, 2016, 64(11):2965-2975.

[8] Auguin N, Morales-Jimenez D, Mckay M, et al. Robust shrinkage M-estimators of large covariance matrices[C]∥Proceedings of 2016 IEEE Statistical Signal Processing Workshop. Palma de Mallorca, Baleares, Spain:IEEE, 2016:1-4.

[9] Pillai S U, Kwon B H. Forward/backward spatial smoothing techniques for coherent signal identification[J]. IEEE Transactions on Acoustics Speech & Signal Processing, 1989, 37(1):8-15.

[10] Pal P, Vaidyanathan P P. A novel array structure for directions-of-arrival estimation with increased degrees of freedom[C]∥Proceedings of 2010 IEEE International Conference on Acoustics Speech and Signal Processing. Dallas, TX, US:IEEE, 2010:2606-2609.

[11] Gershman A B, Bohme J F. Improved DOA estimation via pseudo-random resampling of spatial spectrum[J]. IEEE Signal Processing Letters, 1997, 4(2):54-57.

[12] Vasylyshyn V. Removing the outliers in root-MUSIC via pseudo-noise resampling and conventional beamformer[J]. Signal Processing, 2013, 93(12):3423-3429.

[13] Qian C, Huang L, So H C. Improved unitary root-MUSIC for DOA estimation based on pseudo-noise resampling[J]. IEEE Signal Processing Letters, 2014, 21(2):140-144.

[14] Shaghaghi M, Vorobyov S A. Subspace leakage analysis and improved DOA estimation with small sample size[J]. IEEE Tran-sactions on Signal Processing, 2015, 63(12):3251-3265.

[15] Bai Z D, Silverstein J W. No eigenvalues outside the support of the limiting spectral distribution of large-dimensional sample covariance matrices[J]. Annals of Probability, 1998, 26(1):316-345.

[16] Bai Z D, Silverstein J W. Exact separation of eigenvalues of large dimensional sample covariance matrices[J]. Annals of Probability, 1999, 27(3):1536-1555.

[17] 李華, 白志東, 肖玉山. 大維隨機矩陣的漸進特征[J]. 東北師大學報:自然科學版, 2014, 46(4):1-8. LI Hua, BAI Zhi-dong, XIAO Yu-shan. The asymptotic properties of the large dimension random matrix[J]. Journal of Northeast Normal University:Natural Science Edition, 2014, 46(4): 1-8.(in Chinese)

[18] Mestre X, Lagunas M. Modified subspace algorithms for DoA estimation with large arrays[J]. IEEE Transactions on Signal Processing, 2008, 56(2): 598-614.

[19] Stergiopoulos S, Sullivan E J. Extended towed array processing by an overlap correlator[J]. Journal of the Acoustical Society of America, 1989, 86(1): 158-171.

NovelTargetDirection-of-arrivalEstimationMethodforUnderwaterSmall-scaleMovingArray

GUO Tuo, WANG Ying-min, ZHANG Li-chen

(School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, Shaanxi, China)

TB566

A

1000-1093(2017)09-1779-07

10.3969/j.issn.1000-1093.2017.09.015

2017-01-17

國家自然科學基金項目(61401362)

郭拓(1986—)，男，博士研究生。E-mail:guotuonwpu@126.com

王英民(1963—)，男，教授，博士生導師。E-mail:ywang@nwpu.edu.cn