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

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

引 言

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

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

1 陣列結(jié)構(gòu)與數(shù)據(jù)模型

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

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

圖1 三平行立體線陣結(jié)構(gòu)

式中：

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分別是四個子陣接收數(shù)據(jù)矩陣；

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

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

分別是子陣Y，Z，W 相對于子陣X的旋轉(zhuǎn)不變關(guān)系矩陣.

2 去冗余PM算法實現(xiàn)

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

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

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

根據(jù)傳播算子定義可以得到

式中，P為傳播算子.

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

①估計陣列的協(xié)方差矩陣

②對R進行如下分塊

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

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

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

得到傳播算子P后，構(gòu)造矩陣

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

根據(jù)式（11）可以得到

由式（12）可得

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

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

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

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

3 超定非線性方程組求解

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

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

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

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

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

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

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

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

4 仿真實驗與分析

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

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

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

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

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

圖4 均方根誤差隨快拍數(shù)變化曲線

5 結(jié) 論

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

［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］.北京理工大學(xué)學(xué)報，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］.系統(tǒng)工程與電子技術(shù)，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算法參數(shù)配對新方法［J］.系統(tǒng)工程與電子技術(shù)，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）