三平行線陣改進傳播算子算法的波達方向估計

胡偉偉 王昌明 陸建山 孟翔飛

（南京理工大學機械工程學院，江蘇 南京210094）

引 言

波達方向估計能夠定位信號源，是陣列信號處理的重要研究內容之一，廣泛應用于雷達、聲納、移動通訊、地震等領域［1-2］.自20世紀70年代末開始，出現以多重信號分類算法（Multiple Signal Classifi-cation，MUSIC）和旋轉不變子空間算法（Estimating Signal Parameters Via Rotational Invariance Techniques，ESPRIT）為代表的子空間分解類算法，標志著波達方向估計向現代超分辨測向技術飛躍［3-4］.但是由于MUSIC和ESPRIT算法需要對數據的協方差矩陣進行特征分解，計算量很大，不利于對信源的快速定位.為了解決這一問題，出現了PM算法，PM算法在求解信號與噪聲子空間時，以線性變換代替特征分解，大大降低運算量［5］.文獻［6］采用三平行立體陣，把傳播算子原理引入到ESPRIT算法中，以線性變換代替特征分解求解旋轉不變關系矩陣，并可以實現自動配對；文獻［7］采用最小冗余線陣，降低陣列冗余度，最大限度地減少陣元數目，利用傳播算子算法，降低了運算量；文獻［8］在對相干分布式信源進行DOA估計時，采用三平行立體陣，利用均勻線陣導向矢量間存在的旋轉不變關系，把對導向矢量的求解轉化為對旋轉不變關系矩陣的求解過程，在其中引入傳播算子算法，提高DOA估計精度.

上述文獻與傳統子空間分解類算法相比大大降低了運算量.但是在估算傳播算子時，信號的協方差矩陣中存在較大的冗余度，可以剔除其中的冗余數據，降低協方差矩陣的維數，通過這樣的方法進一步降低運算量，同時構造新的傳播算子陣，以補償孔徑損失；由三平行立體陣可以建立一個三維方程組，而估計的參數只有方位角和俯仰角這兩項，上述文獻未對該超定方程組的處理給出方法，本文引入非線性最小二乘法，以解決這一問題.

1 陣列結構與數據模型

陣列結構如圖1所示，該結構由三個分別位于x軸，x-y平面和x-z平面的均勻線陣組成，陣元間距為d.其中位于x軸上有M＋1個陣元，分為X和Y兩個子陣，X陣由前M 個陣元構成，Y陣由后M個陣元構成；位于x-z和x-y平面上都有M個陣元，分別用子陣Z和W 表示.假設有K個遠場窄帶信源入射到該陣列，θi和φi分別表示第i個信源的方位角和俯仰角.

以坐標原點為參考點，則四個子陣的接收信號可表示為：

圖1 三平行立體線陣結構

式中：

X（t）＝［x1（t），x2（t），…，xM（t）］T；

Y（t）＝［x2（t），x3（t），…，xM＋1（t）］T；

Z（t）＝［z1（t），z2（t），…，zM（t）］T；

W（t）＝［w1（t），w2（t），…，wM（t）］T分別是四個子陣接收數據矩陣；

A（Θ）＝［α（θ1，φ1），α（θ2，φ2），…，α（θK，φK）］是陣列流型矩陣；為導向向量；λ是入射信號波長；

S（t）＝［s1（t），s2（t），…，sK（t）］T是源信號向量；Nx（t），Ny（t），Nz（t），Nw（t）均為期望為0，方差為σ2的加性高斯白噪聲；

分別是子陣Y，Z，W 相對于子陣X的旋轉不變關系矩陣.

2 去冗余PM算法實現

文獻［6］提出的ESPRIT-PM算法構造的陣列流型矩陣是由4個子陣的流型矩陣直接組成，子陣間有重疊陣元，估算傳播算子的陣列協方差矩陣含有冗余數據，考慮剔除其冗余數據，以降低運算量.

假設總陣列的廣義方向為B，其前K行線性獨立，則B可以分解為

其中B1取B的前K 行，B2為B的后3 M＋1-K行.

根據傳播算子定義可以得到

式中，P為傳播算子.

傳播算子可以通過整個陣列的協方差矩陣計算得到，計算步驟如下：

①估計陣列的協方差矩陣

②對R進行如下分塊

式中，R1和R2分別是（3 M＋1）×K和（3 M＋1）×（3 M＋1-K）維矩陣.

③ 如果無噪聲干擾，則矩陣R1與R2存在如下關系

但是實際有噪聲存在，故式（9）等號不成立，因此，得到傳播算子的最小二乘解

得到傳播算子P后，構造矩陣

式中：I為K 維單位陣；P′1，P′2，P′3和P′4為P′均勻劃分的四個子陣，每個子陣都是M行，則得

根據式（11）可以得到

由式（12）可得

式中，P′＋1代表取該矩陣的廣義逆.

從式（13）可看出，Φ1（Θ），Φ2（Θ），Φ3（Θ）對角元素分別對應P′＋1P′2，P′＋1P′3，P′＋1P′4的特征值［λ11，λ12，…，λ1K］， ［λ21，λ22，…，λ2K］， ［λ31，λ32，…，λ3K］.如果存在多個信源，則需要對特征值進行配對，否則不能得到正確結果，根據陣列結構特點，可以通過文獻［9］完成特征值配對過程，在此不再贅述.

由此可以得到與第i個信源方位角θi和俯仰角φi相關的方程組

從以上分析可以看出，利用本文算法在構建協方差矩陣與傳播算子估計時，都只需要對（3 M＋1）×（3 M＋1）維矩陣進行復乘與線性運算，而ESPRIT-PM算法在執行相應操作時，需要對（4 M）×（4 M）維矩陣進行復乘與線性運算，在陣元數較大的場合，計算量節約很明顯；同時，構造出P′矩陣，可以補償陣列孔徑損失.

3 超定非線性方程組求解

由于有噪聲干擾存在，式（14）是一個超定非線性方程組，根據非線性最小二乘法可以求得θi，φi.由最小二乘法可以得到以下優化函數

式（15）中f1，f2，f3作如下定義：

針對該非線性最小二乘問題，可以用Gauss-Newton法求解，具體步驟如下：

① 給定解的初始估計 ［θi，φi］（1），置k＝1；

② 如果 ［θi，φi］（k）滿足精度要求，停止迭代；如果不滿足精度要求，執行③；

③ 解方程組ATkAkδ（k）＝-ATkrk得δ（k）；

④ 置 ［θi，φi］（k＋1）＝［θi，φi］（k）＋δ（k），k：＝k＋1后轉②.

式中：rk是向量函數r＝［f1f2f3］T第k次迭代時的值；A是向量函數r的Jacobian矩陣；Ak是A第k次迭代時的值.

4 仿真實驗與分析

為了驗證方法的正確性與有效性，采用圖1所示的陣列結構進行計算機仿真實驗.取子陣陣元數M＝8，子陣間距為0.5λ，兩個等功率信源分別位于（50°，40°）和（70°，60°），Gauss-Newton法以 ESPRITPM算法結果為迭代初始值，迭代精度為10-3.定義均方根誤差為

實驗1：在信噪比10dB，快拍數300的條件下，分別對ESPRIT-PM算法和本文算法進行100次獨立仿真實驗，可以得出兩種算法星座圖，如圖2所示.從圖中可以看出，本文算法有較高的估計準確度和穩健性.

實驗2：本例對ESPRIT-PM算法和本文提出的方法進行對比.快拍數N＝300，信噪比從0dB變化到30dB，對每個信噪比做100次的獨立仿真實驗，如圖3所示.從圖中可以看出，本文所用算法在低信噪比條件下明顯優于ESPRIT-PM算法，可見本文算法在低信噪比情況下具有更強的適應能力.

圖3 均方根誤差隨信噪比變化曲線

實驗3：本例對ESPRIT-PM算法和本文提出的方法進行對比.在信噪比20dB條件下，快拍數從100變化到1 500，對每個快拍數做100次的獨立仿真實驗，如圖4所示.從圖中可以看出，隨著快拍數的增加，兩信源的均方根誤差皆減小并且趨于接近，本文所用算法得到的均方根誤差明顯低于ESPRIT-PM算法得到的均方根誤差，可見本文算法具有更高的測量精度.

圖4 均方根誤差隨快拍數變化曲線

5 結 論

在三平行立體線陣中，劃分的子陣間存在部分重疊陣元，其直接構造的協方差矩陣含有冗余數據，這樣估算傳播算子時會增加較大計算量.本文采用子陣合并的方式，使協方差矩陣由（4 M）×（4 M）維降至（3 M＋1）×（3 M＋1）維，減少傳播算子估計時的計算量；對傳播算子陣進行重構，補償陣列孔徑損失，提高二維角估計精度；引入非線性最小二乘法，解決由三平行立體線陣構建的三維超定方程組求解二維角度問題.仿真結果表明，本文算法在低信噪比以及相同快拍數條件下，估計性能良好.該算法保持傳統傳播算子算法計算量小的優勢，通過剔除協方差矩陣冗余數據進一步降低運算量，同時具備較高的估計精度和穩健性，因此有廣闊的應用前景和實用價值.

［1］HAYKIN S，REILLY J P，KEZYS V，et al.Some aspects of array signal processing［J］.IEEE Proceedings-F，1992，139（l）：1-26.

［2］KRIM H，VIBERG M.Two decades of array signal processing research：the parametric approach［J］.IEEE Signal Processing Magazine，1996，13（4）：67-94.

［3］SCHMIDT R O.Multiple Emitter Location and Signal Parameter Estimation［J］.IEEE Trans AP，1986，34（3）：276-280.

［4］PAUIRAJ A，ROY R，KAILATH T.A subspace rotation approach to signal parameter estimation［J］.Proc of IEEE，1986，74（7）：1044-1046.

［5］MUNIER J，DELISLE G Y.Spatial analysis using new properties of the cross-spectral matrix［J］.IEEE Transactions on Signal processing，1991，39（3）：746-749.

［6］蘇淑靖，顏景龍，馬維賢，等.基于傳播算子的二維波達方向估計新算法［J］.北京理工大學學報，2008，28（7）：602-605.SU Shujing，YAN Jinglong，MA Weixian，et al.A new algorithm for 2DDOA estimation based on propagator method［J］.Transactions of Beijing Institute of Technology，2008，28（7）：602-605.（in Chinese）

［7］趙大勇，陳 超，刁 鳴.基于最小冗余線陣的二維傳播算子DOA估計［J］.系統工程與電子技術，2011，33（4）：724-727.ZHAO Dayong，CHEN Chao，DIAO Ming.Propagator method for two-dimensional DOA estimation based on minimum redundancy linear array［J］.Systems Engineering and Electronics，2011，33（4）：724-727.（in Chinese）

［8］ZHENG Zhi，LI Guangjun，TENG Yunlong.2DDOA Estimator for multiple coherently distributed sources using modified propagator［J］.Circuits System Signal Process，2012，31：255-270.

［9］刁 鳴，繆善林.一種二維ESPRIT算法參數配對新方法［J］.系統工程與電子技術，2007，29（8）：1226-1229.DIAO Ming，MIAO Shanlin.New method of parameter matching for 2DESPRIT algorithms［J］.Systems Engineering and Electronics，2007，29（8）：1226-1229.（in Chinese）