雙平行線陣的水下二維DOA估計算法

姚 澤，張 歆,2

（1.西北工業(yè)大學(xué)航海學(xué)院，陜西西安 710072；2.鵬城實驗室，廣東深圳 518000）

關(guān)鍵字：雙平行線陣；傳播算子算法；二維波達方向估計；自主式水下航行器

0 引 言

自主水下航行器（Autonomous Underwater Vehicles,AUV）在進行水下搜救、環(huán)境監(jiān)測、數(shù)據(jù)采集等領(lǐng)域的需求越來越多，其工作也越來越復(fù)雜[1]。單個AUV所實現(xiàn)的工作是有限的，而集群AUV協(xié)同系統(tǒng)具有空間、時間、功能的分布性，能擴展單AUV的功能，但是同時也提高了工作的復(fù)雜度[2]。在集群AUV進行水下作業(yè)時，協(xié)同作業(yè)成為了其完成復(fù)雜任務(wù)的重要環(huán)節(jié)，而多個AUV水下協(xié)同工作時，彼此之間的位置信息就十分重要。對多個AUV目標(biāo)的二維波達方向（Direction of Arrival,DOA）估計是AUV水下協(xié)同作業(yè)必須要解決的問題。目前水下AUV一般采用通信、定位一體化方法進行位置信息的確定，利用通信信號時延進行距離估計，再進一步對其進行二維DOA估計，通常使用長、短基線對水下目標(biāo)進行估計。在小型AUV上，基陣陣列長度等受載體體積和換能器尺寸的限制；小尺寸、低功耗的工作環(huán)境條件對于水下DOA估計提出了更高的要求。實際應(yīng)用時，所采用的算法和尺寸需要結(jié)合上述條件進行選擇，還要保證一定的估計精度。

DOA估計應(yīng)用于雷達、生物、聲吶等領(lǐng)域[3-4]，其中，二維DOA估計得到了廣泛關(guān)注[5-7]。要獲得高精度的DOA估計需要合適的陣型和DOA估計算法。在二維DOA估計中，一般采用直線陣、L型陣、平面陣或矢量傳感器，并且按照陣元排布也可以分為均勻陣、非均勻陣。近年來，非均勻陣中的互質(zhì)陣[8-9]得到了關(guān)注，其通過互質(zhì)關(guān)系對陣列進行虛擬重構(gòu)的特性使得陣列測向自由度得到提升，但其重構(gòu)造成計算量加大，導(dǎo)致整體處理過程的復(fù)雜度提高。在常用平面陣列結(jié)構(gòu)中，由等距線陣構(gòu)成的衍生陣近年來由于結(jié)構(gòu)簡單、易于實施、估計效果良好等受到人們的廣泛重視，平行陣就是其中的一種。

二維多重信號分類算法（Multiple Signal Classification,MUSIC）[10]和二維旋轉(zhuǎn)不變子空間算法（Estimation of Signal Parameters via Rotational Invariance Techniques,ESPRIT）[11]是DOA估計的經(jīng)典算法，并且都屬于子空間算法。子空間算法的核心就是對接收信號的協(xié)方差矩陣進行特征值分解或奇異值分解，從而獲得信號子空間或者噪聲子空間，進行二維譜峰搜索，這類計算較為復(fù)雜，往往涉及到龐大的計算量。文獻[12]中指出，子空間算法在大規(guī)模多輸入多輸出（Multiple-Input Multiple-Output,MIMO）系統(tǒng)中，這種劣勢尤為明顯，這在多源信號估計的條件下同樣存在。文獻[13]中提出了最大似然法，根據(jù)最大似然準(zhǔn)則，對輸出數(shù)據(jù)進行時空二維處理來獲得估計結(jié)果。之后又提出了子空間算法的改進算法，但算法計算復(fù)雜度較高。復(fù)雜度過高不利于算法在處理能力較弱、或?qū)挠邢拗频膱鼍皯?yīng)用，AUV作為低功耗水下平臺，其分配在定位方面的處理資源有限，無法進行高復(fù)雜度的定位處理計算。二維傳播算子（Propagator Method,PM）算法不需要互相關(guān)矩陣的特征值分解和接收數(shù)據(jù)的奇異值分解就可以快速得到噪聲子空間[14]，計算復(fù)雜度大大降低，陣列信號的處理性能也得到了極大地提高，更多地應(yīng)用于低功耗系統(tǒng)中。近些年來PM算法的應(yīng)用越來越多，基于PM算法的二維DOA估計方法得到了廣泛關(guān)注。

本文設(shè)計適合了AUV平臺的水下二維DOA估計算法，選取了雙平行線陣。雙平行線陣由兩組平行線陣組成，兩者相互垂直構(gòu)成陣列模型，具有平移不變性，有較好的陣列分辨率，易于構(gòu)造協(xié)方差矩陣。且本文采用的PM算法可以通過協(xié)方差矩陣快速得到噪聲子空間，不需要對協(xié)方差矩陣的特征分解計算。輸出DOA估計結(jié)果時，避免了二維譜峰搜索，可以得到自動配對的二維DOA估計結(jié)果，具有較好的二維DOA估計性能，減少了計算量，計算方法復(fù)雜度相對較低。這種復(fù)雜度較低的算法，對于整個系統(tǒng)的功耗控制較好，利用通信信號與聲學(xué)定位一體化方法，將高頻通信信號作為AUV水下信號源進行DOA定位處理，高頻信號的短波長特性，使得小尺寸陣列在水下DOA估計中得以使用，整體估計方法滿足了AUV水下DOA估計平臺對于小尺寸、低功耗的要求，并且有較好的估計精度。

1 算法模型

1.1 雙平行線陣陣列模型

陣列幾何結(jié)構(gòu)如圖1所示，在平行線陣基礎(chǔ)上，進行垂直方向的擴展，得到雙平行線陣。在x-y平面中，有Lx1，Lx2，Ly1，Ly2四條均勻線陣，其中Lx1和Ly1由相鄰陣元間距為d的N+1個陣元組成；Lx2和Ly2由陣元間距為d的N個陣元組成；Lx1和Lx2平行于x軸放置，彼此平行間距為d，構(gòu)成沿x軸方向的平行線陣，Ly1和Ly2平行于y軸放置，彼此平行間距為d，構(gòu)成沿y方向的平行線陣。

圖1 雙平行線陣陣列結(jié)構(gòu)圖Fig.1 Configuration of double parallel linear array

式中：λ是信號波長，Lx2的導(dǎo)向矢量矩陣Ax2是Ax1的前N行，定義：

其中，diag（·）表示對角矩陣。同理可得與y軸平行的平行線陣Ly1和Ly2的陣列輸出信號yy1和yy2，表示為

Ly2的導(dǎo)向矢量矩陣Ay2是Ay1的前N行，定義：

兩種平行線陣的信號形式相似，本文以與x軸平行的平行線陣為主進行二維DOA算法介紹，同理可得與y軸平行的線陣結(jié)果。

1.2 二維DOA估計算法

1.2.1 陣列輸出重構(gòu)

將與x軸平行的平行線陣的陣列輸出合并表示為

其中，Ax為與x軸平行的平行線陣的矢量導(dǎo)向矩陣

1.2.2 基于PM算法的二維DOA估計

將上一步重構(gòu)得到的x軸方向平行線陣的矢量導(dǎo)向矩陣Ax進行分塊：

其中：右上角表示矩陣的行列數(shù)，且A1是非奇異矩陣，這樣分解得到的A1、A2是具有線性關(guān)系的兩個子矩陣，這種關(guān)系可以表示為

這里的PH是一個行列數(shù)為K×（2N+1-K）的矩陣，定義為傳播算子矩陣，之后的算法圍繞這里的PH矩陣進行。考慮接收數(shù)據(jù)的協(xié)方差矩陣為

由于實際應(yīng)用中，快拍數(shù)為有限值，所以無法精確得到式（14）中的Rx理論值，一般使用離散形式的計算進行近似，即：

同樣，對此協(xié)方差矩陣進行分塊處理：

理論上，這種分塊處理得到的R1、R2也是具有線性關(guān)系的兩個子矩陣，在無噪聲情況下可表示為

這里的P就是傳播算子矩陣。但是實際應(yīng)用中要考慮到噪聲的情況，無法得到上述精確的對應(yīng)關(guān)系，一般使用近似值來估計傳播算子矩陣P：

得到傳播算子的估計值后，對其進行重構(gòu)：

取Pxx的前N行得到矩陣Px1，由式（12）和式（20）可得[15]：

分別取Pxx的后N行、Pxy的前（N-1）行和Pxy的后（N-1）行得到矩陣Px2、Px3、Px4。定義矩陣：

其中，amp（·）表示取幅度值，angle（·）表示取相角值。同理可得y軸方向平行線陣的二維DOA估計結(jié)果。

將兩個方向的平行線陣進行并行處理方式，獨立進行DOA估計，當(dāng)兩個獨立估計結(jié)果得出后，進行融合操作，提高DOA估計精度。

融合規(guī)則：第m次蒙特卡洛實驗得到的估計結(jié)果，將兩個方向方位角和俯仰角的估計值進行均值處理，得到最終DOA角度估計結(jié)果。

算法具體實現(xiàn)步驟總結(jié)如下：

（1）根據(jù)式（15）構(gòu)建協(xié)方差矩陣；

（2）式（18）和（19）利用協(xié)方差矩陣分解重構(gòu)得到估計值P和Pc；

（3）對矩陣pc分解重構(gòu)得到Px1、Px2、Px3、Px4，并且得到bk和A′；

（5）同理進行y軸方向平行陣的參數(shù)計算，基于y軸方向平行線陣進行二維DOA估計；

（6）結(jié)果融合得到雙平行線陣的DOA估計結(jié)果。

2 仿真分析

為驗證本文算法的性能指標(biāo)，使用Matlab仿真軟件對算法進行仿真。仿真中，所采用的雙平行線陣的兩個子陣陣元數(shù)分別為4和3，仿真使用水下通信定位信號的頻率為37.5 kHz，對應(yīng)的陣列陣元間距為1 cm，進行算法性能測試。采用求均方根誤差（Root Mean Square Error,RMSE）作為衡量算法DOA估計性能的標(biāo)準(zhǔn)，定義為

式中：M表示蒙特卡洛實驗次數(shù)，θmk和φmk分別表示第m次實驗中、第k個信號的估計值。仿真結(jié)果中誤差單位均為（°）。

2.1 估計效果測試

圖2給出了在信噪比（Signal to Noise Ratio,SNR）為15 dB的情況下，進行200次蒙特卡洛實驗，快拍數(shù)S=200，信源數(shù)K=3，角度分別為（?25°,45°）、（5°,15°）和（30°,60°）時算法的仿真情況。

從圖2可以看出，本文采用的算法在設(shè)置條件下可以較好地對目標(biāo)方位進行估計，得到的DOA估計結(jié)果基本和目標(biāo)位置吻合，驗證了本文所使用高頻通信信號配合小尺寸陣列進行水下DOA估計的可行性。

圖2 SNR為15dB情況下雙平行線陣及傳播算子算法的二維DOA估計結(jié)果（S=200,K=3）Fig.2 Two-dimensional DOA estimation result of double parallel linear arrays and propagator method when SNR is 15 dB （S=200,K=3）

2.2 估計誤差隨參數(shù)變化情況

下面測試算法的多源信號估計性能，進行 200次蒙特卡洛實驗，快拍數(shù)為200，計算均方根誤差，仿真結(jié)果如圖3所示。

圖3 SNR為15dB情況下DOA估計的均方根誤差隨信源數(shù)變化情況（S=200）Fig.3 Variation of RMSE of DOA estimation with number of sources when SNR is 15 dB （S=200）

從圖3可以看出信源數(shù)對本文所用水下DOA估計算法的影響，由于本文采用小尺寸陣列，子陣陣元數(shù)分別為4和3，進而當(dāng)信源數(shù)為4時，由于陣元個數(shù)的限制，本文算法精度下降。當(dāng)信源數(shù)小于4時，本文算法在考慮小尺寸陣列的情況下有較好的估計效果。圖4、5測試了算法估計誤差隨快拍數(shù)和信噪比的變化情況。

從圖4～5可以看出，與圖3相對應(yīng)，在SNR為15 dB、S=200的條件下，對于3個目標(biāo)的估計RESM最大在1°左右，有較好的估計效果。

圖4 SNR為15 dB情況下DOA估計的均方根誤差隨快拍數(shù)變化情況（K=3）Fig.4 Variation of RMSE of DOA estimation with number of snapshots when SNR is 15 dB （K=3）

圖5 DOA估計的均方根誤差隨SNR變化情況（K=3,S=200）Fig.5 Variation of RMSE of DOA estimation with SNR（K=3,S=200）

3 結(jié) 論

本文為解決集群AUV水下協(xié)同工作所遇到的定位問題，并且滿足AUV水下平臺對于小尺寸、低功耗的要求，提出了一種雙平行線陣結(jié)合PM方法的二維水下DOA估計算法，利用雙平行陣列融合提高了一定的估計精度；使用PM算法，避免二維譜峰搜索造成的大量計算，實現(xiàn)了二維角度自動配對，降低了復(fù)雜度，從而達到低功耗要求；采用通信與水下定位一體化方法，利用高頻通信信號作為定位信號源，使用小尺寸基陣，開展集群 AUV水下定位方法，對多源信號估計精度較好，滿足了集群AUV協(xié)同工作時的定位要求。本文所采用的陣列結(jié)構(gòu)簡單，算法復(fù)雜度較低，對整個系統(tǒng)功耗控制較好，在載體限制陣列尺寸以及平臺低功耗的要求下，有較好的角度估計性能。