基于稀疏互質(zhì)L型陣列的二維測(cè)向算法

邵 華 蘇衛(wèi)民 顧 紅 王 燦

（南京理工大學(xué)電子工程與光電技術(shù)學(xué)院，江蘇 南京 210094）

引 言

在雷達(dá)、聲納和通信等領(lǐng)域，對(duì)輻射源進(jìn)行定位［1－3］，常需要先確定入射信號(hào)的二維波達(dá)角（2DDOA，包括方位角和俯仰角）。針對(duì)2D－DOA估計(jì)問題，陣列信號(hào)處理中已有許多行之有效的方法被提出，例如二維多重信號(hào)分類法［4］（2D－MUSIC），但大多數(shù)測(cè)向方法只適合于等間距面陣。與等間距面陣相比，稀疏面陣能在不損失陣列孔徑的前提下減小所需的陣元數(shù)，因此，研究稀疏面陣結(jié)構(gòu)及其相應(yīng)的波達(dá)角（DOA）估計(jì)算法具有極大的意義。

其中，文獻(xiàn)［5］提出二維旋轉(zhuǎn)不變（2D－ESPRIT）算法，它通過子陣間x軸和y軸的旋轉(zhuǎn)不變關(guān)系估計(jì)方位角和俯仰角，缺點(diǎn)是需要3個(gè)結(jié)構(gòu)相同的子陣，硬件成本高；基于此缺點(diǎn)，文獻(xiàn)［6］、［7］利用四階累積量［8－10］的陣列擴(kuò)展特性，分別沿x軸和y軸構(gòu)造與原陣列結(jié)構(gòu)相同的虛擬陣列，然后再由2D－ESPRIT算法進(jìn)行2D－DOA估計(jì)，即二維虛擬旋轉(zhuǎn)不變（2D－VESPA）算法，然而其要求參考陣元間距（原陣列參考陣元和虛擬陣列參考陣元之間的間距）不大于信號(hào)半波長(zhǎng)0.5λ；為突破陣元間距不大于0.5λ的限制條件，文獻(xiàn)［11］、［12］提出嵌套陣列，文獻(xiàn)［13］、［14］提出互質(zhì)陣列，利用子陣間的嵌套或互質(zhì)關(guān)系形成陣元間距不大于0.5λ的均勻矩形虛擬陣列，而允許物理陣元間距大于0.5λ.雖然這兩種算法是利用虛擬陣元消除了測(cè)向模糊，但是其實(shí)質(zhì)是僅采用一步就完成了解模糊。這導(dǎo)致正確解模糊時(shí)，系統(tǒng)允許兩個(gè)陣元接收信號(hào)間相位差的測(cè)量誤差（簡(jiǎn)稱系統(tǒng)容差）較小。

針對(duì)以上問題，本文基于模轉(zhuǎn)換［15］逐步解模糊的思想，首先構(gòu)造了一類由兩個(gè)具有相同互質(zhì)結(jié)構(gòu)的線陣組成的L型陣列，其中線陣1沿x軸放置，線陣2沿y軸放置；其次，利用一維虛擬旋轉(zhuǎn)不變（1DVESPA）算法估計(jì)各基線所對(duì)應(yīng)的模糊相位差；然后，結(jié)合基線間的互質(zhì)關(guān)系逐步解模糊，分別得到x軸和y軸中最長(zhǎng)基線對(duì)應(yīng)的無模糊方向余弦；最后，采用x軸和y軸特征向量（估計(jì)不同坐標(biāo)軸方向余弦時(shí)分別進(jìn)行特征分解）間的對(duì)應(yīng)關(guān)系，實(shí)現(xiàn)不同坐標(biāo)軸方向余弦的匹配，最終得到方位角和俯仰角的估計(jì)值。與2D－VEPSA算法相比，本算法最大的優(yōu)勢(shì)在于利用基線間的互質(zhì)關(guān)系突破0.5λ的限制條件，提高測(cè)向精度的同時(shí)保證較大的系統(tǒng)容差。

1 信號(hào)模型

如圖1所示，L型陣列由x軸子陣（M＋1個(gè)全向陣元）和y軸子陣（M＋1個(gè)全向陣元）組成，兩個(gè)子陣在坐標(biāo)原點(diǎn)處共用陣元0，并且兩子陣內(nèi)相鄰陣元間的基線長(zhǎng)度分別為Dx，1，…，Dx，M和Dy，1，…，Dy，M，其中Dx，m＝Dy，m＝Dm，m＝1，…，M，基線D1最長(zhǎng)，各基線之間滿足如下比例關(guān)系

式中，Pm和Qm為互質(zhì)的正整數(shù)，并且Q2＜…＜QM.

圖1 L型陣列結(jié)構(gòu)示意圖

假設(shè)K個(gè)入射方向?yàn)椋é?，φ1），…，（θK，φK）的非高斯獨(dú)立窄帶遠(yuǎn)場(chǎng)信號(hào)源入射到上述陣元總數(shù)為2M＋1的L型陣列天線上，其中0≤φk＜2π和0≤θk＜π分別表示第k個(gè)入射信號(hào)的方位角和俯仰角。那么，L型陣列在t時(shí)刻接收的信號(hào)矢量為

式中：yx，m（t）和yy，m（t）分別表示在x軸和y軸上第m個(gè)陣元的接收信號(hào)；s（t）＝ ［s1（t），…，sK（t）］T為信號(hào)向量；w（t）＝ ［w0（t），…，w2M（t）］T表示陣列接收噪聲，噪聲為平穩(wěn)、時(shí)間和空間都互不相關(guān)的高斯白噪聲，且與信號(hào)相互獨(dú)立；A＝ ［a1，…，aK］是（2M＋1）×K維的陣列導(dǎo)向矩陣；ak＝］T，k＝1，…，K，其中

式中：qx，k，m＝exp（－jψx，k，m）和qy，k，m＝exp（－jψy，k，m），m＝1，…，M，其中ψx，k，m＝／λ和＝2πvkDy，m／λ，uk＝sinθkcosφk和vk＝sinθksinφk分別表示x軸和y軸方向余弦。當(dāng)基線Dx，m＝Dy，m＞λ／2時(shí)，測(cè)量相位差φx，k，m∈ ［0，2π）和φy，k，m∈ ［0，2π）具有周期性模糊，即

式中，Kx，k，m∈Z和Ky，k，m∈Z分別表示φx，k，m和φy，k，m的模糊數(shù)。這導(dǎo)致x軸和y軸測(cè)量方向余弦也具有模糊性。利用L型陣列輸出的N個(gè)快拍數(shù)據(jù)進(jìn)行無模糊2D－DOA估計(jì)。另外，為表述簡(jiǎn)潔，下文均省略時(shí)間t.

2 算法描述

2.1 四階累積量矩陣的構(gòu)造

利用圖1中L型陣列的特殊結(jié)構(gòu)，并結(jié)合四階累積量的定義［7－10］及空間接收信號(hào)的數(shù)學(xué)模型，構(gòu)造如下四階累積量矩陣

式中，m＝1，…，M，且

式中，γ4，sk表示第k個(gè)信號(hào)的四階累積量。由四階累積量的陣列擴(kuò)展特性可知，式（8）表示的物理含義為：原陣列沿x軸平移xm距離后產(chǎn)生的虛擬陣列，如圖2（a）所示，其中陣元m為虛擬陣列的參考陣元，稱這類虛擬陣列為x軸虛擬陣列組；式（9）則表示原陣列沿y軸平移ym距離后產(chǎn)生的虛擬陣列，如圖2（b）所示，其中陣元M＋m為虛擬陣列的參考陣元，類似的稱其為y軸虛擬陣列組。為了估計(jì)所有相鄰陣元間的相位差，選擇陣元0作為原陣列的參考陣元，其他2M個(gè)陣元分別為2M個(gè)虛擬陣列的參考陣元，那么可以定義如下矩陣

對(duì)式（13）進(jìn)行奇異值分解（SVD）后，可得信號(hào)子空間Es＝ ［，…，，…，］T.

2.2 估計(jì)x軸基線對(duì)應(yīng)的測(cè)量相位差

抽取A的相應(yīng)行可以構(gòu)造原陣列和x軸虛擬陣列組對(duì)應(yīng)的導(dǎo)向矩陣Ax＝ ［AT，ATΛx，1，…，ATΛx，M］T.同樣由Es可得它們對(duì)應(yīng)的信號(hào)子空間Es，x＝ ［，，…，］T.因?yàn)镋s，x和Ax都表示信號(hào)子空間，所以它們之間可以相互轉(zhuǎn)換，令可逆矩陣T為轉(zhuǎn)換矩陣，則Es

為了估計(jì)基線對(duì)應(yīng)的測(cè)量相位差，令Fx，m＝（Es，x，m－1）?Es，x，m，m＝1，…，M，其中Es，x，0＝Es，0，并結(jié)合式（16）得

其中：（·）?表示矩陣偽逆；Λx，0是元素為0的K維方陣。對(duì)Fs，x，m進(jìn)行特征分解后得Λx，m－Λx，m－1，求其對(duì)角元素的相位可得測(cè)量相位差的估計(jì)值＝，…］和特征向量Tx，m.

2.3 x軸方向余弦的估計(jì)

由于基線越長(zhǎng)，測(cè)向精度越高。因此，為得到高精度x軸方向余弦估計(jì)，必須先通過模轉(zhuǎn)換解模糊算法估計(jì)最長(zhǎng)基線Dx，1對(duì)應(yīng)的無模糊相位差。然而，由模轉(zhuǎn)換解模糊算法的原理可知，它需要結(jié)合相同信號(hào)、不同基線所對(duì)應(yīng)的測(cè)量相位差才能正確解模糊。又由x軸基線對(duì)應(yīng)的測(cè)量相位差的求解過程可知，和是分別進(jìn)行特征分解后得到的，其中m≠n，m∈1［］M，n∈［1M］，故和信號(hào)源三者之間的順序具有任意性，所以對(duì)它們的配對(duì)是正確解模糊的前提。

根據(jù)矩陣的特征分解可知，特征值與特征向量是一一對(duì)應(yīng)的，即和的配對(duì)問題轉(zhuǎn)化成Tx，m和Tx，n的配對(duì)問題。假設(shè)表示Tx，m中第k1列向量表示Tx，n中第k2列向量，當(dāng)它們對(duì)應(yīng) 相 同 信 號(hào) 源 時(shí)，＝ 1，否 則?1.基于以上分析，對(duì)Tx，m和Tx，n做如下操作其中，矩陣G中1的位置反映和間的配對(duì)關(guān)系。不失一般性，假設(shè)以的順序?yàn)閰⒖?，利用上述方法分別對(duì)｛，m＝1，…，M｝進(jìn)行配對(duì)后得

結(jié)合基線Dx，m和配對(duì)后的測(cè)量相位差估計(jì)值，m＝1，…，M，k＝1，…，K，利用模轉(zhuǎn)換解模糊算法解最長(zhǎng)基線D1，x對(duì)應(yīng)測(cè)量相位差估計(jì)值φx，k，1的模糊，得到Dx，1對(duì)應(yīng)的無模糊測(cè)量相位差，再經(jīng)簡(jiǎn)單的計(jì)算后得高精度、無模糊x軸方向余弦估計(jì)值經(jīng)過逐步解模糊得出，與單步解模糊相比，它能保證較大的系統(tǒng)容差較［15］。

2.4 y軸方向余弦的估計(jì)

與x軸方向余弦的估計(jì)類似，y軸方向余弦的估計(jì)分為以下幾步：

1）由Es和A分別構(gòu)造原陣列和y軸虛擬陣列組對(duì)應(yīng)的信號(hào)子空間Es，y＝ ［，，…，］T和導(dǎo)向矩陣Ay＝ ［AT，ATΛy，1，…，ATΛy，M］T；

2）為了估計(jì)基線對(duì)應(yīng)的測(cè)量相位差，m＝1，…，M，令

其中，Es，y，0＝Es，0，Λy，0是元素為0的K維方陣；

3）對(duì)Fs，y，m進(jìn)行特征分解后得基線Dy，m對(duì)應(yīng)測(cè)量相位差的估計(jì)＝［，…，］和特征向量Ty，m，m＝1，…，M；

4）以為基準(zhǔn)，利用G＝，m，m＝2，…，M中1的位置對(duì)和進(jìn)行配對(duì)，最后得到配對(duì)后y軸各基線所對(duì)應(yīng)的測(cè)量相位差

5）對(duì)于第k個(gè)信號(hào)，利用模轉(zhuǎn)換算法解模糊相位差的模糊數(shù)，得到無模糊測(cè)量相位差再經(jīng)簡(jiǎn)單的計(jì)算后得高精度、無模糊y軸方向余弦估計(jì)

2.5 方位角和俯仰角的估計(jì)

由于和間并不一定對(duì)應(yīng)同一信號(hào)源，因此，先對(duì)它們進(jìn)行配對(duì)。與前面的配對(duì)類似，采用G＝中1的位置實(shí)現(xiàn)和的匹配，記匹配后的x軸和y軸方向余弦分別為和.根據(jù)方向余弦的定義，可得第k個(gè)入射信號(hào)方位角和俯仰角的估計(jì)值

3 計(jì)算機(jī)仿真試驗(yàn)

通過計(jì)算機(jī)仿真來比較本文算法和2D－VEPSA算法的測(cè)向性能。由于本文算法要求各陣元間距滿足一定的互質(zhì)關(guān)系，而2D－VEPSA算法要求參考陣元間距不大于信號(hào)半波長(zhǎng)，其他陣元位置任意，即它們對(duì)陣列結(jié)構(gòu)的要求不同。因此，在仿真中，兩種算法分別采用陣元數(shù)相同、結(jié)構(gòu)不同的天線a和天線b來比較它們的測(cè)向性能，其中天線a是陣元間距滿足互質(zhì)關(guān)系的9陣元L型陣，各陣元都為全向陣，并且為了減小天線a的體積，已對(duì)輔助陣元的位置做了調(diào)整，但并不影響基線間的比例關(guān)系，如圖3所示；天線b也是9陣元L型陣，其陣元1和陣元5的位置為 （λ／2，0，0）和（0，λ／2，0），其他陣元分布與天線a的相同，它作為天線a對(duì)比天線。

圖3 仿真中本文算法采用的天線a

必須說明的是：在仿真中，本文算法首先以陣元0到陣元｛m，m＝1，…，4｝為參考陣元，形成與原陣列相距 ｛（Dm，0，0），m＝1，…，4｝的x軸虛擬陣列組，然后以陣元0到陣元｛m，m＝5，…，8｝為參考陣元，形成與原陣列相距 ｛（0，Dm，0），m＝1，…，4｝的y軸虛擬陣列組，利用虛擬陣列與原陣列間的旋轉(zhuǎn)不變關(guān)系估計(jì)兩坐標(biāo)軸各基線對(duì)應(yīng)的模糊相位差，經(jīng)解模糊得到兩坐標(biāo)軸最長(zhǎng)基線（D1）對(duì)應(yīng)的無模糊方向余弦，簡(jiǎn)單計(jì)算后得三個(gè)信號(hào)的2D－DOA估計(jì)值；2D－VEPSA算法以陣元0、1和5為參考陣元，形成與原陣列相距 （λ／2，0，0）和（0，λ／2，0）的虛擬陣列組，利用虛擬陣列組與原陣列的旋轉(zhuǎn)不變關(guān)系估計(jì)陣元0和陣元1及陣元0和陣元5之間對(duì)應(yīng)的無模糊方向余弦，從而估計(jì)三個(gè)信號(hào)的2D－DOA.另外，仿真實(shí)驗(yàn)中使用均方根誤差（RMSE）對(duì)算法測(cè)向性能進(jìn)行評(píng)估，其中均方根誤差定義為

式中：和分別表示第k個(gè)目標(biāo)俯仰角和方位角估計(jì)值；N表示Monte－Carlo試驗(yàn)次數(shù)。

仿真A：假設(shè)存在3個(gè)非高斯、獨(dú)立、窄帶遠(yuǎn)場(chǎng)信號(hào)源投射到這兩個(gè)天線上，3個(gè)信源方位角分別為10°、20°和30°，俯仰角分別為15°、25°和35°，信源間相互獨(dú)立，且與噪聲互不相關(guān)，噪聲為平穩(wěn)、時(shí)間和空間都不相關(guān)的高斯白噪聲。對(duì)于測(cè)量噪聲，三個(gè)信號(hào)源有相同的信噪比（SNR），SNR取值范圍為－10～20dB，取值步長(zhǎng)為5dB，每次獨(dú)立仿真采用2 000次快拍數(shù)據(jù)，獨(dú)立重復(fù)200次Monte－Carlo試驗(yàn)。

如圖4給出了兩種算法的2D－DOA估計(jì)的RMSE隨SNR的變化曲線。圖5是本文算法中正確解模糊的概率隨SNR的變化曲線。結(jié)合圖4和圖5可以看出，在正確解模糊的情況下（RSN≥0 dB），本文算法相對(duì)于2D－VEPSA算法在測(cè)向精度上具有較大改善，這主要是因?yàn)楸疚乃惴ɡ媒饽：惴ㄍ黄屏?D－VEPSA算法中參考陣元間距不能大于半波長(zhǎng)的限制，而測(cè)向精度與參考陣元的間距成正比，因此，本文算法能達(dá)到更高的測(cè)向精度。但是，當(dāng)RSN＜0dB，本文算法的測(cè)向精度逐漸劣于2D－VEPSA算法，其原因是更低的信噪比超出了正確解模糊的系統(tǒng)容差，從而存在解模糊錯(cuò)誤，導(dǎo)致本文算法測(cè)向性能急劇下降。

仿真B：假設(shè)兩個(gè)陣列接收信號(hào)的SNR為10 dB，當(dāng)快拍數(shù)從200變化到2 000，變化步長(zhǎng)為200時(shí)，估計(jì)兩種算法的測(cè)向性能隨快拍數(shù)的變化，其他條件與仿真A相同。

如圖6是兩種算法的2D－DOA估計(jì)的RMSE隨快拍數(shù)的變化曲線。很明顯，兩種算法的測(cè)向性能隨快拍數(shù)的增大而提高，但變化的陡峭程度不如SNR.另外，本文算法測(cè)向 RMSE要小于2DVEPSA算法，其原因也是本文算法構(gòu)造稀疏互質(zhì)L型陣列，利用陣元間的互質(zhì)關(guān)系解模糊，擴(kuò)大陣列的孔徑，提高測(cè)向精度。

圖4 兩種算法的2D－DOA估計(jì)的RMSE隨SNR的變化曲線

圖5 本文算法中正確解模糊概率隨SNR的變化

圖6 兩種算法的2D－DOA估計(jì)的RMSE隨快拍數(shù)的變化曲線

4 結(jié) 論

提出了一種基于稀疏互質(zhì)L型陣列的2D－DOA估計(jì)算法。與2D－VEPSA算法相比，利用陣元間的互質(zhì)關(guān)系，解決測(cè)向精度與測(cè)向模糊之間的矛盾，擴(kuò)大陣列孔徑，提高測(cè)向精度，同時(shí)利用逐步解模糊，保證了較大的系統(tǒng)容差。

［1］DOGANCAY K. Relationship between geometric translations and TLS estimation bias in bearings－only target localization［J］.IEEE Transactions on Signal Processing，2008，56（3）：1005－1017.

［2］DOGANCAY K.Exploiting geometric translations in TLS based robot localization from landmark bearings［C］／／17th European Signal Processing Conference.Glasgow，Scotland，24－28August，2009：95－99.

［3］RAO S K，RAJESWARI K R，LINGAMURTY K S.Unscented Kalman filter with application to bearingsonly target tracking［J］.IETE Journal of Research，2009，55（2）：63－67.

［4］CHAN A Y J，LITVA J.MUSIC and maximum techniques on two－dimensional DOA estimation with uniform circular array［J］.IEE Proceedings Radar，Sonar and Navigation，1995，142（3）：105－114.

［5］KEDIA V S，CHANDNA B.A new algorithm for 2－D DOA estimation［J］.Signal Process，1997，60（3）：325－332.

［6］LIU T H，MENDEL J M.Azimuth and elevation direction finding using arbitrary array geometries［J］.IEEE Transactions on Signal Processing，1998，46（7）：2061－2065.

［7］DOGAN M C，MENDEL J M.Applications of cumulants to array processing－part I：aperture extension and array calibration［J］.IEEE Transactions on Signal Processing，1995，43（5）：1200－1216.

［8］齊子森，郭 英，王布宏，等.基于四階累積量的共形陣列波達(dá)方向估計(jì)算法［J］.電波科學(xué)學(xué)報(bào)，2011，26（4）：735－744.

QI Zisen，GUO Ying，WANG Buhong，et al.DOA estimation algorithm for conformal array based on fourthorder cumulants［J］.Chinese Journal of Radio Science，2011，26（4）：735－744.（in Chinese）

［9］沈怡平，滕升華.基于四階累積量的穩(wěn)健盲波束形成算法［J］.電波科學(xué)學(xué)報(bào)，2008，23（6）：1056－1060.

SHEN Yiping，TENG Shenghua.Fourth－order cumulant based on robust blind beamforming algorithm［J］.Chinese Journal of Radio Science，2008，23（6）：1056－1060.（in Chinese）

［10］劉學(xué)斌，韋 崗，季 飛.基于四階累積量擴(kuò)展孔徑的線陣設(shè)計(jì)［J］.電波科學(xué)學(xué)報(bào)，2006，21（1）：126－130.

LIU Xuebin，WEI Gang，JI Fei.Design of linear array with augmented aperture based on fourth－order cumulant［J］.Chinese Journal of Radio Science，2006，21（1）：126－130.（in Chinese）

［11］PAL P，VAIDYANATHAN P P.Nested arrays：a novel approach to array processing with enhanced degrees of freedom［J］.IEEE Transactions on Signal Processing，2010，58（8）：4167－4181.

［12］PAL P，VAIDYANATHAN P P.Two dimensional nested arrays on lattices［C］／／IEEE ICASSP.Prague，Czech Republic，May 22－27，2011：2548－2551.

［13］PAL P，VAIDYANATHAN P P.Sparse sensing with co－prime samplers and arrays［J］.IEEE Transactions on Signal Processing，2011，59（2）：573－586.

［14］VAIDYANATHAN P P，PAL P.Theory of sparse coprime sensing in multiple dimensions［J］.IEEE Transactions on Signal Processing，2011，59（8）：3592－3608.

［15］WILLETT P K.Modulo conversion method for estimating the direction of arrival［J］.IEEE Transactions on Aerospace and Electronic Systems，2000，36（4）：1391－1396.