強相干干擾下微弱信號波達方向估計

朱 偉 陳伯孝

（西安電子科技大學 雷達信號處理重點實驗室，陜西 西安710071）

引 言

微弱信號的檢測和參數估計廣泛應用于雷達、通信、聲納和地震探測等領域.當存在強干擾源時，系統對微弱信號波達方向的估計精度會受到很大的影響.在陣列信號處理中通常的處理方法是使用超分辨算法對干擾和待檢測信號同時處理，或者根據信號和干擾在空域、頻域和時域的差異，在不同的域中進行分離、抑制，對待檢測信號進行波達方向估計.前一種方法在一定程度上可以同時估計微弱信號和強干擾的波達方向，但是當干擾和信號相干或者位于一個波束寬度內時，干擾源與信號源的波達方向（Direction Of Arrival，DOA）估計效果較差，且估計信噪比門限較高；后一種方法主要有方向圖零點綜合、信號分離和數字波束形成技術等.

近年來，學者們提出了多種算法來進行強干擾環境下的弱信號DOA估計［1-14］.文獻［1］提出了使用粒子群等優化算法進行陣列方向圖綜合，以降低副瓣電平并在干擾方向形成零點.文獻［2］提出的松弛算法、文獻［3］提出的潔凈技術以及文獻［4］提出的盲信號分離法，均是采用信號分離的方法使包含多個信號的陣列輸出數據分離成幾個數據塊，而具體的某一個數據塊只包含某一個信號的信息，從而分離出強干擾，達到干擾抑制的目的.數字自適應波束形成是通過控制波束的零點去對準干擾的來向，使用穩健的波束形成算法來形成零點，主要有對角加載算法［5］、最壞條件最優法［6］和迭代自適應最小方差法［7］等方法，但是當信號和干擾相干時，現有的大多數魯棒波束形成算法性能都會嚴重下降.在以上算法基礎上，文獻［8］提出利用陣列環境中強干擾的先驗知識構造阻塞矩陣來抑制干擾，從而實現對特定區域內低信噪比信號的DOA估計.文獻［9］提出將接收的數據矢量劃分成兩個子矢量之后再重新構造新的數據矢量，在某些不確定集的約束下進行波束形成抑制相干干擾.文獻［10］提出在已知干擾源個數和入射方向的前提下，將陣列劃分成若干子陣，在子陣上進行波束形成以抗干擾，再對加權后的子陣進行微弱信號的DOA估計，該算法在信號與干擾不在一個波束內且互不相干時有良好性能.文獻［11］提出一種自適應加權空間平滑解相干，然后利用線性約束最小方差準則得到子陣波束形成器最佳權矢量，再利用子陣間的相位關系對全陣進行波束形成的方法.文獻［12］提出將陣列劃分為兩個虛擬子陣，分別對子陣進行波束形成來抑制干擾，然后利用子陣間的相位偏移來對弱信號進行DOA估計，該算法能顯著消除同信道干擾的影響.文獻［13］提出了使用迭代自適應波束形成來抑制干擾.文獻［14］提出了一種迭代超分辨處理方法.以上算法在信號與干擾相干或者干擾位于主瓣內時性能較差.

針對空間干擾與待檢測信號相干且空間間隔較小的情況，本文提出了一種強相干干擾下微弱信號DOA估計的新方法.該算法首先對強干擾源的個數和入射方向進行精確估計，再估計微弱信號源個數和參考入射方向，然后將陣列劃分成兩個虛擬子陣，使用波束形成來進行干擾抑制并提高待估計信號的信噪比，最后利用子陣間的相位偏移使用比相單脈沖進行DOA估計，并對子陣波束形成和DOA估計進行迭代運算以提高估計精度.

1 理論分析

1.1 信號模型

假設空間陣列是由M個陣元組成的均勻線陣，陣元間距為d，有P個相干微弱信號和K個相干強干擾以平面波形式入射到陣列上，P＋K＜M，入射方向分別為θSp（p＝1，…，P）和θJk（k＝1，…，K）.則陣列接收數據可表示為

式中：X（t）為M×1維陣列接收數據矢量；a（θ）＝［1，ej2πdsin（θ）/λ，…，ej2π（M-1）dsin（θ）/λ］T為θ 方 向 的 導 向矢量，其中λ為信號波長，上標T表示矩陣轉置；sSp（t）和sJk（t）分別表示第p個信號和第k個干擾，信號和干擾相干，sSp（t）＝cSps（t），sJk（t）＝cJks（t），cSp和cJk為復常數；N（t）為M×1維噪聲矢量，各陣元噪聲統計獨立，服從零均值、方差為σ2的復高斯分布.將式（1）右邊兩個求和項用矩陣表示為

式中：A（θ）＝［a（θS1），…，a（θSP），a（θJ1），…，a（θJK）］；

s（t）＝［cS1，…，cSP，cJ1，…，cJK］Ts（t）.

1.2 本文算法描述

本文算法的框圖如圖1所示，由干擾估計與目標粗估計、子陣波束形成和DOA估計三個模塊組成.首先進行強干擾估計，在陣列上抑制干擾之后再得到目標的參考方向；將陣列劃分成兩個虛擬子陣，在子陣上進行波束形成來抑制干擾并提高待估計信號的信噪比；利用兩個子陣的相位中心偏移來得到目標的精確DOA結果.

圖1 本文算法框圖

1.2.1 干擾估計與目標粗估計

對強干擾源的個數和入射方向進行精確估計.由于信號和干擾為相干信號，對數據協方差矩陣進行特征分解后，信號子空間退化為一維子空間，而噪聲子空間的維數擴大為M-1維，無法通過特征值分解直接得到空間干擾的個數.通常情況下，干擾遠強于待檢測信號和噪聲，相干信號源數目的估計方法［15］有平滑秩法、矩陣分解法和蓋氏源方法等，相干信號源DOA估計的超分辨算法［15］有解相干多重信號分類算法、最大似然算法和Toeplitz近似法等.在高信噪比下，以上算法都能得到精確結果，在此不進行贅述，得到干擾的DOA估計為Jk（k＝1，…，K）.

在實際的諸如無線通信系統［16］等應用中，可以獲得信號源的一些先驗知識（如信號源個數和參考入射方向等）.在雷達、聲納等實際應用中，可以先在陣列上進行波束形成以抑制干擾［10］，然后再對加權后的陣列使用幅相估計 （Amplitude and Phase Estimate，APES）算法［17］得到微弱信號的個數和粗估計，將粗估計結果作為弱信號的參考入射方向，滿足關系為

1.2.2 子陣波束形成

在子陣波束形成模塊中，將M個陣元劃分成兩個相同大小的虛擬子陣.為了充分利用陣列孔徑，通常有兩種劃分方法：最大重疊子陣法（Maximum O-verlapping Subarrays，MOSs）［18］和共軛子陣法（Conjugate Subarrays，CSs）［19］.

1）最大重疊子陣

最大重疊子陣包含兩個M-1個陣元組成的虛擬子陣，其中陣列的前M-1個陣元組成子陣A，陣列的后M-1個陣元組成子陣B，如圖2（a）所示.因此，子陣A和子陣B的接收數據分別為：

式中：

AA＿M（θ）＝［aA＿M（θS1），…，aA＿M（θSP），aA＿M（θJ1），…，

從式（4）和式（5）可以得出

式中，Φ＝diag ｛ej2πdsin（θS1）/λ， …，ej2πdsin（θSP）/λ，ej2πdsin（θJ1）/λ，…，ej2πdsin（θJK）/λ｝，diag表示對角陣.

2）共軛子陣

共軛子陣包含兩個M個陣元組成的虛擬子陣，子陣A由一個虛擬陣元和陣列的前M-1個陣元組成，子陣B為當前陣列，如圖2（b）所示.由于共軛子陣比最大重疊子陣多一個陣元，因此共軛子陣能獲得更優的波束形成結果，但是共軛子陣的前提條件為信號包絡為實包絡，即s（t）＝s＊（t），上標＊表示共軛，文獻［19-21］討論了實際中信號滿足實包絡的條件.因此，子陣A和子陣B的接收數據分別為：

式中：AA＿C（θ）＝［aA＿C（θS1），…，aA＿C（θSP），aA＿C（θJ1），…，aA＿C（θJK）］，

從式（7）和式（8）可以得出

式中：Φ＝diag｛ej2πdsin（θS1）/λ，…，ej2πdsin（θSP）/λ，ej2πdsin（θJ1）/λ，…，ej2πdsin（θJ K）/λ｝.

從式（6）和式（9）可以看出：共軛子陣和最大重疊子陣具有一樣的旋轉不變性，旋轉因子均為Φ，所不同的是最大重疊子陣的陣元數為M-1，而共軛子陣的陣元數為M，但共軛子陣中無孔徑損失子陣的獲得是基于信號包絡為實包絡，而最大重疊子陣則可運用復包絡的情況.若無干擾且各信號源之間相互獨立，在子陣A和子陣B間運用旋轉不變子空間（Estimating Signal Parameters Viarotational Invariance Techniques，ESPRIT）算法，可得到各信號源的波達方向.

對于第p個待檢測信號，進行子陣波束形成時需要在第p個待檢測信號方向上得到最大的空間增益，而且抑制全部干擾，同時對待檢測信號之外的其他信號也有一定的抑制，即子陣波束形成權值矢量wp需要滿足以下線性方程組：

式中，上標H表示共軛轉置；

上述線性方程組無法直接計算，利用已估計的前p-1個信號的方向Sl（l＝1，…，p-1）、后P-p＋1個信號的參考方向Rl（l＝p，…，P）和已估計的干擾方向Jk（k＝1，…，K）得到以下線性方程組：

將方程組（11）寫為矩陣形式有

式中：

于是，可得子陣波束形成權值wp為

分別對虛擬子陣A和虛擬子陣B進行波束形成，波束形成輸出為

1.2.3 波達方向估計

通過波束形成后，干擾和除待檢測信號之外的其他信號均得到了不同程度的抑制，兩個子陣的相位中心間距為d，相應子陣A與子陣B間的相位中心偏移φp為

為了進一步提高信噪比，將yA（t）和yB（t）的L個快拍數據進行時域相參積累，得到SA和SB，SA和SB的幅度近似相等且相位差為φp，可以使用比相單脈沖［21］進行DOA估計，且有

利用相參積累結果SA和SB形成和信號S∑與差信號SΔ：

式中Im（）表示取虛部.

為了提高微弱信號的DOA估計精度，對子陣波束形成和DOA估計進行迭代運算，即用估計出的Sp代替式（13）中的Rp，并繼續利用式（15）重新計算子陣波束形成權值wp，再次進行子陣波束形成和DOA估計.設置一定的終止條件，可以對待檢測信號DOA進行較為準確的估計.

1.2.4 算法流程

綜上所述，算法流程總結如下：

1）使用相干信源估計方法［15］來估計干擾源個數，并采用解相干超分辨算法［15］得到強干擾的DOA估計結果Jk（k＝1，…，K），然后在陣列上進行波束形成以抑制干擾［10］，再對加權后的陣列使用幅度相位估計（Amplitude and Phase Estimation，APES）算法［17］進行微弱信號的粗估計，得到參考方向Rp（p＝1，…，P）；

2）將陣列劃分成兩個虛擬子陣，信號包絡為復包絡時按照最大重疊子陣劃分，信號包絡為實包絡時劃分成共軛子陣；

3）對第p個信號進行步驟4至步驟6的過程；

4）通過式（12）至式（15）得到子陣波束形成的初始權值w（1）p，w（1）p的上標表示第幾次迭代，然后利用式（16）分別對虛擬子陣A和虛擬子陣B進行波束形成；

5）對波束形成輸出結果進行相參積累，按照式（20）生成和信號與差信號，并通過式（22）得到第p個信號的 DOA 初始估計值的上標表示第幾次迭代；

為了進一步提高估計精度，可以將6）中得到的P個DOA估計值Sp（p＝1，…，P）作為參考方向Rp（p＝1，…，P），重復進行3）至6）的過程.

2 實驗結果分析

為了分析本文算法的性能，將本文算法與文獻［10］和文獻［12］中的算法進行比較.考慮一個由10個全向陣元組成的均勻線陣，陣元間距d為半波長，3 dB波束寬度為10.15°.假設接收信號包絡均為復包絡，將陣列按最大重疊子陣的方式構造虛擬子陣.

實驗1：多干擾下多目標的DOA估計性能比較

假設兩個相干干擾方向為-50°和16°，干噪比均為80dB，三個相干信號分別從-11°、0.5°和12°方向入射，12°方向的信號與16°方向的干擾處在同一個波束寬度內，信噪比為0dB，信號的參考入射方向為-10°、0°和10°，信號與干擾相干，快拍數為100.迭代終止條件為δ＝0.001°.10 000次 Monte-Carlo實驗DOA估計結果的統計直方圖如圖3所示.表1給出了三種算法的DOA估計均方根誤差.

表1 三種算法DOA估計均方根誤差

從以上統計直方圖可以清晰地看出：文獻［10］算法明顯無法分辨多干擾下的多目標，尤其是無法分辨一個波束寬度內的目標和相干干擾；文獻［12］算法能夠檢測多干擾下多目標，但是估計結果偏離真實角度較大，表明文獻［10］和文獻［12］算法對目標與相干干擾位于一個波束寬度內時不能進行DOA估計，而本文算法能很好地進行多相干干擾下多目標DOA估計，特別是在干擾與目標處于一個波束寬度內的情況下也能得到較好的結果.

實驗2：多干擾下單目標的DOA估計性能隨信噪比的變化

假設四個相干干擾方向為-50°、-30°、3°和40°，其中3°方向的干擾位于主瓣波束內，干噪比均為80dB，一個信號從-1°方向入射，信號的參考入射方向為1.5°，信號與干擾相干，快拍數為100，信噪比從-10dB變化至10dB.迭代終止條件為δ＝0.001°.1 000次 Monte-Carlo實驗統計的DOA估計均方根誤差如圖4所示.均方根誤差定義為

式中：P表示目標個數；N表示Monte-Carlo實驗次數；Spn為第p個目標的第n次獨立實驗估計結果；θSp為第p個目標的真實值.

從圖4結果曲線可看出：本文算法對多干擾下單目標的DOA估計性能優于文獻［10］和文獻［12］的算法，當信噪比大于8dB時，三種算法的估計性能趨于一致.圖5給出了輸入信噪比為10dB時子陣波束形成的方向圖，從中可以看出干擾能得到較好的抑制.

實驗3：多干擾下單目標的DOA估計性能隨快拍數的變化

信號和干擾的仿真條件與實驗2相同，信噪比為0dB，快拍數從10變化至1 000.迭代終止條件為δ＝0.001°.1 000次 Monte-Carlo實驗統計的DOA估計均方根誤差如圖6所示.

圖6 快拍數變化時的性能曲線

可見文獻［10］、文獻［12］和本文算法DOA估計精度隨快拍數增加而提高，其中本文算法估計性能最優.

實驗4：多干擾下單目標的DOA估計性能與導向矢量誤差的關系

信號和干擾的仿真條件與實驗2相同，信噪比為0dB，快拍數為100，信號的參考入射方向從-14.5°變化至1.5°.迭代終止條件為δ＝0.001°.1 000次Monte-Carlo實驗統計的DOA估計均方根誤差如圖7所示.

圖7 DOA估計性能與目標參考方向的關系

從圖7可看出：信號參考入射方向與真實入射方向相差大于一個波束時，性能急劇下降，算法無法收斂到目標的真實方向，但是參考入射方向與真實入射方向相差在一個波束寬度之內時，不影響目標的估計精度，但影響收斂速度，導向矢量誤差越大，迭代次數越多，反之迭代次數越少.

實驗5：空間鄰近目標的估計性能

假設四個相干干擾方向為-50°、-30°、30°和50°，干噪比均為80dB，兩相干信號從-2.5°方向和2.5°方向入射，兩個信號處于半個波束寬度內，信號與干擾相干，信號的參考入射方向為-10°和10°，快拍數為100，考察不同信噪比下本文算法對兩目標的估計性能，迭代終止條件為δ＝0.001°.由于兩個信號空間間距很近，為了提高估計精度，在此可以將步驟6中得到的結果作為參考方向，再進行一次步驟3至步驟6的過程，得到更精確的估計結果.信噪比分別為0dB和10dB時，10 000次Monte-Carlo實驗DOA估計結果的統計直方圖如圖8所示.圖9給出了信噪比從-6dB變化至10dB，10 000次Monte-Carlo實驗統計的DOA估計均方根誤差.

從圖8和圖9可以看出：信噪比越高，空間鄰近目標的分辨性能越好，DOA估計精度隨信噪比的增加顯著提高.

3 結 論

提出一種均勻線陣下強相干干擾環境中進行微弱信號DOA估計的算法，利用子陣波束形成來抑制干擾并提高待估計信號的信噪比，然后使用比相單脈沖進行DOA估計，并且對子陣波束形成和DOA估計進行迭代運算以提高DOA估計精度.對于強相干干擾環境下系統對微弱信號的DOA估計問題，特別是當信號和干擾處在同一波束寬度內時，本文算法具有良好的性能.計算機仿真結果表明，相比文獻［10］和文獻［12］算法，本文算法具有更低的估計偏差和估計方差，性能更優.

［1］ISMAIL T H，HAMICI Z M.Array pattern synthesis using digital phase control by quantized particle swarm optimization［J］.IEEE Trans on Antennas and Propagation，2010，58（6）：2142-2145.

［2］LI J，STOICA P.Efficient mixed-spectrum estimation with applications to target feature extraction［J］.IEEE Trans.on Signal Processing，1996，44（2）：281-295.

［3］AKHDAR O，MOUHAMADOU M，CARSENAT D，et al.A new CLEAN algorithm for angle of arrival denoising［J］.IEEE Antennas and Wireless Propagation Letters，2009，36（4）：478-481.

［4］冶繼民，張賢達，金海紅.超定盲信號分離的半參數統計方法［J］.電波科學學報，2006，21（3）：331-336.YE Jimin，ZHANG Xianda，JIN Haihong.Semi-parametric statistical approach for overdetermined blind source separation［J］.Chinese Journal of Radio Science，2006，21（3）：331-336.（in Chinese）

［5］曾 操，廖桂生，楊志偉.一種加載量迭代搜索的穩健波束形成［J］.電波科學學報，2007，22（5）：779-784.ZENG Cao，LIAO Guisheng，YANG Zhiwei.Diagonal loading level estimation for robust beamforming［J］.Chinese Journal of Radio Science，2007，22（5）：779-784.（in Chinese）

［6］ELNASHAR A.Efficient implementation of robust adaptive beamforming based on worst-case performance optimisation［J］.IET Signal Processing，2008，2（4）：381-393.

［7］NAI S E，SER W，YU Z L，et al.Iterative robust minimum variance beamforming［J］.IEEE Trans on Signal Processing，2011，59（4）：1601-1611.

［8］陳 輝，蘇海軍.強干擾/信號背景下的DOA估計新方法［J］.電子學報，2006，34（3）：530-534.CHEN Hui，SU Haijun.A new approach to estimate DOA in presence of strong jamming/signal suppression［J］.Acta Electronica Sinica，2006，34（3）：530-534.（in Chinese）

［9］胡 航，景秀偉.基于近似理想方向圖的子陣級超分辨測向方法［J］.電波科學學報，2007，22（4）：646-658.HU Hang，JING Xiuwei.Subarray level superresolution direction finding approaches based on approximate ideal patterns［J］.Chinese Journal of Radio Science，2007，22（4）：646-658.（in Chinese）

［10］柴立功，羅景青.一種強干擾條件下微弱信號DOA估計的新方法［J］.電子與信息學報，2005，27（10）：1517-1520.CHAI Ligong，LUO Jingqing.A novel algorithm for weak signals DOA estimation under intensive interferences［J］.Journal of Electronics ＆Information Technology，2005，27（10）：1517-1520.（in Chinese）

［11］周 圍，張德民，吳 波，等.相干環境下LCMV自適應陣列抗干擾問題研究［J］.電子與信息學報，2007，29（7）：1604-1607.ZHOU Wei，ZHANG Demin，WU Bo，et al.Study on interference suppression for LCMV adaptive array in coherent environment［J］.Journal of Electronics ＆Information Technology，2007，29（7）：1604-1607.（in Chinese）

［12］WANG N Y，AGATHOKLIS P，ANTONIOU A.A new DOA estimation technique based on subarray beamforming［J］.IEEE Trans.on Signal Processing，2006，54（9）：3279-3289.

［13］DOISY Y，DERUAZ L，BEEN R.Interference suppression of subarray adaptive beamforming in presence of sensor dispersions［J］.IEEE Trans on Signal Processing，2010，58（8）：4195-4212.

［14］BLUNT S D，CHAN T，GERLACH K.Robust DOA estimation：the reiterative superresolution（RISR）algorithm［J］.IEEE Trans on Aerospace and Electronic Systems，2011，47（1）：332-346.

［15］王永良，陳 輝，彭應寧等.空間譜估計理論與算法［M］.北京：清華大學出版社，2009.

［16］OJANPERA T，PRASAD R.Wideband CDMA for third generation mobile communications［M］.Norwood：Artech House，1998.

［17］KENNY S R，MOE S，AHMED A K.High accuracy peak location and amplitude spectral estimation via tuning APES method［J］.Digital Signal Processing，2010，20：552-560.

［18］SOON V C，HUANG Y F.An analysis of ESPRIT under random sensor uncertainties［J］.IEEE Trans.on Signal Processing，1992，40（9）：2353-2358.

［19］TAYEM N，KWON H M.Conjugate ESPRIT （CSPRIT）［J］.IEEE Trans on Antennas and Propagation，2004，52（10）：2618-2624.

［20］DELMAS J P.Comments on “Conjugate ESPRIT（C-SPRIT）”［J］.IEEE Trans.on Antennas and Propagation，2007，55（2）：511.

［21］TAYEM N，KWON H M.Reply to comments on“Conjugate ESPRIT （C-SPRIT）”［J］.IEEE Trans.on Antennas and Propagation，2007，52（2）：512-513.

［22］朱 偉，陳伯孝，周 琦.兩維數字陣列雷達的數字單脈沖測角方法［J］.系統工程與電子技術，2011，33（7）：1503-1509.ZHU Wei，CHEN Baixiao，ZHOU Qi.Angle measurement method with digital monopulse for 2-dimensional digital array radar［J］.Journal of Systems Engineering and Electronics，2011，33（7）：1503-1509.（in Chinese）