基于輔助陣元的幅相誤差和DOA 同時估計算法 *

趙海洲，王小哲，孫晨，郭藝奪

（空軍工程大學 防空反導學院，陜西 西安 710051）

0 引言

陣列擾動［1-3］作為在實際中常見的一種非理想條件，會嚴重影響信號處理的性能。考慮到陣列陣元在加工等流程以及實際中的應用，可將陣列擾動總結為3 類：陣元互耦［4-5］、陣元幅相誤差［6-7］和陣元位置誤差［8-9］。其中，陣元幅相誤差在現實陣列中是普遍存在的，也是各國學者研究波達方向（direction of arrival，DOA）估計時考慮的一個熱點問題。

對于陣元幅相誤差的校正，初期提出的算法均需要校準信源作為先驗信息［10-11］。為了實現陣元幅相誤差的自校正，文獻［12-13］利用子空間類算法循環估計信號DOA 和幅相誤差，通過不斷地迭代最終估計出了幅相誤差。該算法雖然實現了陣元幅相誤差的自校正，但多次循環迭代帶來的高復雜度使得算法并不實用。為了避免出現迭代過程，Paulraj 和Kailath 利用ULA 接收信號協方差矩陣具有Toeplitz 矩陣結構的特征，提出了一種基于子空間變換的幅相誤差估計算法［14］，但是該算法需要對協方差矩陣上三角的每一條對角線進行標定，計算過程煩瑣。針對這一問題，文獻［15］提出了一種基于協方差矩陣不同對角線的簡化標定算法，實現了低復雜度的幅相誤差校正。但是這2 類算法只適用于均勻線陣，為了解決受陣列結構限制的問題，文獻［7］中針對二維陣列提出了一種新型的基于子空間的幅相誤差校正算法。但是現有算法普遍計算復雜度偏高，且大多不能同時實現幅相誤差和DOA 的同時估計。

為了解決以上問題，本文提出了一種基于少量精確校正輔助陣元的幅相誤差和DOA 同時估計算法。該算法能夠實現信源功率存在較大差異時，高精度的誤差和角度的同時估計。

1 信號模型

1.1 均勻線陣信號模型

通常，理想陣列可認為是陣元間距為半波長的均勻線陣（uniform linear array，ULA），即陣元間距d滿足d=λ2，其中，λ為波長。假設L個遠場窄帶信號入射到陣元數為M的ULA 上，如圖1 所示。

在圖1 中，θl為第l個入射信號的角度。此時，可以寫出該陣列在θl方向上的導向矢量為

式中：?l=2πdsin(θl)λ，l= 1，2，…，L。

假設第l個信號在k時刻輻射到陣列上的波形為sl(k)，每個信號相互獨立。在不考慮陣元之間存在誤差的情況下，在k時刻陣列接收到的數據可以表示為

式中：n(k) ∈CM×1表示噪聲向量，一般噪聲設為加性高斯白噪聲，且與入射信號相互獨立。

由式（2）可得接收信號的協方差矩陣為

但在實際中，協方差矩陣往往都難以準確獲得，一般都用采樣協方差矩陣（sampling covariance matrix，SCM）來代替，其表達式為

式中：K為采樣快拍總數。

1.2 陣元幅相誤差下的信號模型

假設陣列中第m個陣元的幅度誤差和相位誤差分別為αm，φm，則陣列的幅度誤差矩陣和相位誤差矩陣分別為

式中：diag{·}表示對角化運算符；一般以第1 個陣元作為參考陣元，即有α1= 1，φ1= 0。

在考慮陣元幅相誤差時，模型（2）改寫為

此時，受誤差影響下的導向矢量為

2 幅相誤差和DOA 同時估計算法

2.1 算法原理

根據矩陣運算，式（8）可以進一步寫成

式中：Ψ=ΦΓ為幅相誤差矩陣；?=diag(Ψ)為幅相 誤 差 向 量；T?(θ) ∈CM×M的 第m列 表 示 為Ema(θ)，Em∈CM×M表示為

由子空間基本原理［16］知，在存在幅相誤差時可推出：

由于Un秩為M-L，因此對于任意角度θ都有Q?(θ)為奇異矩陣。若幅相誤差矩陣中的元素出現若干相同元素，則在滿足一定條件時獲得的對應于非入射信號方向的過渡矩陣將會是滿秩矩陣，此時便可利用這一性質來估計幅相誤差。

對于幅相誤差矩陣中出現若干相同元素這一假設，可通過在陣列的一端放置少量已校正陣元或幅相誤差近似的陣元來實現。當陣列中增加的已校正陣元數為P-1 時，陣列的幅度誤差和相位誤差系數有

此時，Γ和Φ可重寫為

式 中 ：ΓP?diag{αP，αP+1，…，αM}；ΦP?diag{ }ejπφP，ejπφP+1，…，ejπφM。

此時，幅相誤差矩陣和幅相誤差向量分別可以表示為

重新定義式（10）為

式中：G=M-P+ 1。

可得

因此，式（11）可重寫為

式中：Q?(θ) ∈CG×G。

當G≤M-L，即L≤P- 1 時，一般角度θ時可以認為Q?(θ)是一個滿秩矩陣。但是當θ為入射信號的DOA 時，也即θ=θl(l= 1，2，…，L)，根據第1 節的分析，Q?(θ)則會成為一個奇異矩陣。此時由式（20）可知，?P是Q?(θ)中0 特征值對應的特征向量。

因此，可構造譜函數為

對式（21）或式（22）進行譜峰搜索，得到前L個最大譜峰對應的位置即可得到入射信號的角度估計值，并可計算出對應方向?的過渡矩陣Q?(?)。由于幅相誤差矩陣是與方向無關的，因此?P可以估計為

式中：ulmin為Q?(?)最小特征值對應的特征向量；(ulmin)1表 示ulmin的 第1 個 元 素。

當入射信號的SNR 較低時，估計得到的角度會有較大誤差，此時利用對應的過渡矩陣Q?(θ)估計?P時會有較大的誤差。因此，當多個入射信號之間的功率相差過大時，在最終估計?P時對得到的過渡矩陣進行篩選，選擇較大功率信號對應的過渡矩陣Q?(θ)來估計?P會提高估計精度。此時有2 種策略：

（1） 選擇最大譜峰對應的過渡矩陣來估計?P

通過譜峰搜索獲得的最大譜峰對應著功率最強的入射信號，選擇此譜峰對應的過渡矩陣Q?max(θ)用來估計?P，此時式（23）重寫為

式中：umax-min為Q?max(θ)最小特征值對應的特征向量。

（2） 剔除掉對應最小譜峰的過渡矩陣來估計?P

假設獲得的最小譜峰對應的過渡矩陣為Q?(?)，在后續 估 計?P時 忽略矩陣Q?(?)，此 時 式（23）重寫為

2.2 算法步驟

綜上，對所提的陣元幅相誤差和DOA 同時估計算法總結如下：

第1 步：估計協方差矩陣，計算噪聲子空間Un；

第2 步：計 算 變 換 矩 陣T?(θ)，構 造 譜 函 數Pdet(θ)；

第3 步：譜峰搜索，找到L個最大譜峰的位置，估計角度，計算? 對應的過渡矩陣Q?(?)；

第4 步：選擇合適的過渡矩陣進行特征分解，估計幅相誤差向量。

需要注意的是，這些是建立在假設L≤P- 1 的基礎上推導完成的。雖然此方法可以實現幅相誤差和DOA 的同時估計，但在應用中算法自由度受到輔助陣元數的限制，而且將多個已校正陣元或幅相誤差近似的陣元作為輔助陣元接收信號意義不大，因此，本方法還是更適合用于前期的陣列校正中。不難發現，所提算法對于陣列的排布并無特殊要求，也即該算法也可適用于非均勻線陣。同樣地，通過修改導向矢量的形式，該算法也可拓展到二維陣列中去。此外，為了保證DOA 估計無模糊，輔助陣列的布置通常設置為陣元間距半波長的ULA。

3 計算機仿真

此仿真主要分析所提的算法在陣元幅相誤差下DOA 和幅相誤差估計性能。仿真條件設置為：待校準陣列為最小陣元間距為0.5、陣元數為12 的ULA，波長λ= 1；輔助陣列是由4 個已校準陣元構成的ULA，放置在待校準陣列一端，且與待校準陣列共同組成陣元數為16 的ULA。蒙特卡羅實驗次數為100。

實驗1： 陣元幅相誤差和DOA 同時估計算法的有效性驗證

此仿真主要驗證所提陣元幅相誤差和DOA 同時估計算法的有效性。為了更普遍地驗證算法性能，仿真中設置幅度和相位誤差分別服從如下分布：

式中：δm和ηm均為相互獨立且均勻分布于[-0.5，0.5]中的隨機變量；σα和σφ分別為δm和ηm的標準偏差，分別設置為0.1 和15°。入射信號角度分別為-30°，-10°，40°，對應的SNR 均為10 dB，采樣快拍數為200。仿真得到的空間譜和幅相誤差估計值分別如圖2 和表1 所示。

表1 陣元幅相誤差估計值Table 1 Estimated sensor gain-phase errors

圖2 空間譜Fig. 2 Spatial spectrum

從圖2 中可以看出，所提出的陣元幅相誤差和DOA 同時估計算法獲得的空間譜可精準地估計出入射信號的DOA；從表1 中可看出，所提算法能夠精確地估計出幅相誤差，也即所提算法有效地實現了陣元幅相誤差和DOA 的同時估計。

實驗2： 陣元幅相誤差和DOA 同時估計算法的DOA 估計性能

此仿真主要為了分析所提陣元幅相誤差和DOA 同時估計算法DOA 估計的性能，以同為子空間類算法的MUSIC（multiple signal classification）算法和DOA 估計值 的CRB（Cramér-Rao bound）作為對比。假設陣列接收到3 個信號，角度分別為-30°，-10°，40°，SNR 均為10 dB。在分析所提算法DOA估計性能與SNR 變化的關系時，快拍數設置為200，SNR 從-10 dB 到30 dB 變化，其余參數保持不變；在分析所提算法DOA 估計性能與快拍數變化的關系時，SNR 設置為10 dB，快拍數從50 到400 變化，其余參數保持不變。仿真結果如圖3 所示。

圖3 DOA 估計值的RMSE 與SNR 的關系Fig. 3 Relationship between RMSE and SNR of estimated DOA

由圖3 可知，隨著SNR 的增加，提出的陣元幅相誤差和DOA 同時估計算法的DOA 估計精度均有所增加；但在低SNR 時，所提算法的DOA 估計精度低于MUSIC 算法，這是由于所提算法在實現DOA 和陣元幅相誤差同時估計時，利用的是噪聲子空間構建的一個更小維度的過渡矩陣，在低SNR 時，難以準確估計出噪聲子空間，此時所提算法利用更低維度的矩陣便難以獲得更高精度的DOA 估計值。

實驗3： 陣元幅相誤差和DOA 同時估計算法陣元幅相誤差估計性能

此仿真主要為了分析所提陣元幅相誤差和DOA 同時估計算法對陣元幅相誤差的估計性能。以文獻［10］中提出的基于子空間的陣元幅相誤差估計算法作為對比算法，同時仿真中還給出了陣元幅相誤差估計值的CRB。為了方便，將采用式（21）和式（22）的所提算法分別稱為“所提算法1”“所提算法2”。為充分展示所提算法在不同入射信號功率下對陣元幅相誤差的估計性能，假設入射信號為1 個SNR 為10 dB、入射角度為-10°的期望信號和2個干噪比（interference-to-noise ratio，INR）為40 dB、入射角度分別為-30°和40°的干擾信號。在分析所提算法幅相誤差估計性能與SNR 變化的關系時，快拍數設置為200，SNR 從-10 dB 到50 dB 變化，其余參數保持不變；在分析所提算法幅相誤差估計性能與快拍數變化的關系時，SNR 設置為10 dB，快拍數從50 到400 變化，其余參數保持不變。蒙特卡羅實驗次數為100，仿真結果如圖4 所示。

圖4 陣元幅相誤差估計值的RMSE 和SNR 的關系，INR =40 dB，K = 200Fig. 4 Relationship between RMSE and SNR of estimated sensor gain-phase errors INR = 40 dB， K = 200

如圖4 所示，相較于對比算法，在存在一個低SNR 信號時，也即多個入射信號存在較大的功率差時，提出的2 種算法對幅度誤差和相位誤差的估計性能接近，雖然2 種算法均更加靠近幅度誤差和相位誤差估計值的CRB，但所提算法2 估計得到的幅相誤差精度更高。隨著信號SNR 的增加，所提的2種算法對幅相誤差的估計精度較對比算法均有所下降，這是由于此時所有信號功率接近，每一個信號均為可用信號，利用越多的信號來估計幅相誤差，估計的精度就越高。隨著SNR 的進一步增大，高功率信號對幅相誤差估計性能的改善更為明顯，因此所提算法1 的估計精度再次超越對比算法，更為接近CRB。

實驗4： 算法運行時間比較

為了直觀地看出本文所提算法在運算量方面的優勢，圖5 給出了P=2，L=100，SNR=10 dB 時，所提算法和文獻［10］中對比算法的運行時間隨陣元數變化曲線。從圖5 中可以看出，當陣元數較小時，所提算法的運行時間與對比算法的運行時間基本相當。但隨著陣元數的增大，所提算法在運算量上的優勢逐步得到體現。

圖5 2 種算法的運行時間對比Fig. 5 Comparison of running time between two algorithms

4 結束語

本文通過設置少量的輔助陣元，利用陣元幅相誤差矩陣為對角矩陣的結構特點，并基于子空間類算法的原理，提出了一種陣元幅相誤差和DOA 同時估計算法。針對信號功率相差較大問題，給出了2 種不同的解決策略，并取得了較好的效果。計算機仿真結果表明：所提算法能夠以較高的精度同時估計出幅相誤差和目標的DOA，尤其在存在功率差異較大信號時，能夠顯著提高對誤差的估計精度。