基于輔助陣元的幅相誤差和DOA 同時(shí)估計(jì)算法 *

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

（空軍工程大學(xué) 防空反導(dǎo)學(xué)院，陜西 西安 710051）

0 引言

陣列擾動(dòng)［1-3］作為在實(shí)際中常見(jiàn)的一種非理想條件，會(huì)嚴(yán)重影響信號(hào)處理的性能?？紤]到陣列陣元在加工等流程以及實(shí)際中的應(yīng)用，可將陣列擾動(dòng)總結(jié)為3 類：陣元互耦［4-5］、陣元幅相誤差［6-7］和陣元位置誤差［8-9］。其中，陣元幅相誤差在現(xiàn)實(shí)陣列中是普遍存在的，也是各國(guó)學(xué)者研究波達(dá)方向（direction of arrival，DOA）估計(jì)時(shí)考慮的一個(gè)熱點(diǎn)問(wèn)題。

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

為了解決以上問(wèn)題，本文提出了一種基于少量精確校正輔助陣元的幅相誤差和DOA 同時(shí)估計(jì)算法。該算法能夠?qū)崿F(xiàn)信源功率存在較大差異時(shí)，高精度的誤差和角度的同時(shí)估計(jì)。

1 信號(hào)模型

1.1 均勻線陣信號(hào)模型

通常，理想陣列可認(rèn)為是陣元間距為半波長(zhǎng)的均勻線陣（uniform linear array，ULA），即陣元間距d滿足d=λ2，其中，λ為波長(zhǎng)。假設(shè)L個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)入射到陣元數(shù)為M的ULA 上，如圖1 所示。

在圖1 中，θl為第l個(gè)入射信號(hào)的角度。此時(shí)，可以寫(xiě)出該陣列在θl方向上的導(dǎo)向矢量為

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

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

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

由式（2）可得接收信號(hào)的協(xié)方差矩陣為

但在實(shí)際中，協(xié)方差矩陣往往都難以準(zhǔn)確獲得，一般都用采樣協(xié)方差矩陣（sampling covariance matrix，SCM）來(lái)代替，其表達(dá)式為

式中：K為采樣快拍總數(shù)。

1.2 陣元幅相誤差下的信號(hào)模型

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

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

在考慮陣元幅相誤差時(shí)，模型（2）改寫(xiě)為

此時(shí)，受誤差影響下的導(dǎo)向矢量為

2 幅相誤差和DOA 同時(shí)估計(jì)算法

2.1 算法原理

根據(jù)矩陣運(yùn)算，式（8）可以進(jìn)一步寫(xiě)成

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

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

由于Un秩為M-L，因此對(duì)于任意角度θ都有Q?(θ)為奇異矩陣。若幅相誤差矩陣中的元素出現(xiàn)若干相同元素，則在滿足一定條件時(shí)獲得的對(duì)應(yīng)于非入射信號(hào)方向的過(guò)渡矩陣將會(huì)是滿秩矩陣，此時(shí)便可利用這一性質(zhì)來(lái)估計(jì)幅相誤差。

對(duì)于幅相誤差矩陣中出現(xiàn)若干相同元素這一假設(shè)，可通過(guò)在陣列的一端放置少量已校正陣元或幅相誤差近似的陣元來(lái)實(shí)現(xiàn)。當(dāng)陣列中增加的已校正陣元數(shù)為P-1 時(shí)，陣列的幅度誤差和相位誤差系數(shù)有

此時(shí)，Γ和Φ可重寫(xiě)為

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

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

重新定義式（10）為

式中：G=M-P+ 1。

可得

因此，式（11）可重寫(xiě)為

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

當(dāng)G≤M-L，即L≤P- 1 時(shí)，一般角度θ時(shí)可以認(rèn)為Q?(θ)是一個(gè)滿秩矩陣。但是當(dāng)θ為入射信號(hào)的DOA 時(shí)，也即θ=θl(l= 1，2，…，L)，根據(jù)第1 節(jié)的分析，Q?(θ)則會(huì)成為一個(gè)奇異矩陣。此時(shí)由式（20）可知，?P是Q?(θ)中0 特征值對(duì)應(yīng)的特征向量。

因此，可構(gòu)造譜函數(shù)為

對(duì)式（21）或式（22）進(jìn)行譜峰搜索，得到前L個(gè)最大譜峰對(duì)應(yīng)的位置即可得到入射信號(hào)的角度估計(jì)值，并可計(jì)算出對(duì)應(yīng)方向?的過(guò)渡矩陣Q?(?)。由于幅相誤差矩陣是與方向無(wú)關(guān)的，因此?P可以估計(jì)為

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

當(dāng)入射信號(hào)的SNR 較低時(shí)，估計(jì)得到的角度會(huì)有較大誤差，此時(shí)利用對(duì)應(yīng)的過(guò)渡矩陣Q?(θ)估計(jì)?P時(shí)會(huì)有較大的誤差。因此，當(dāng)多個(gè)入射信號(hào)之間的功率相差過(guò)大時(shí)，在最終估計(jì)?P時(shí)對(duì)得到的過(guò)渡矩陣進(jìn)行篩選，選擇較大功率信號(hào)對(duì)應(yīng)的過(guò)渡矩陣Q?(θ)來(lái)估計(jì)?P會(huì)提高估計(jì)精度。此時(shí)有2 種策略：

（1） 選擇最大譜峰對(duì)應(yīng)的過(guò)渡矩陣來(lái)估計(jì)?P

通過(guò)譜峰搜索獲得的最大譜峰對(duì)應(yīng)著功率最強(qiáng)的入射信號(hào)，選擇此譜峰對(duì)應(yīng)的過(guò)渡矩陣Q?max(θ)用來(lái)估計(jì)?P，此時(shí)式（23）重寫(xiě)為

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

（2） 剔除掉對(duì)應(yīng)最小譜峰的過(guò)渡矩陣來(lái)估計(jì)?P

假設(shè)獲得的最小譜峰對(duì)應(yīng)的過(guò)渡矩陣為Q?(?)，在后續(xù) 估 計(jì)?P時(shí) 忽略矩陣Q?(?)，此 時(shí) 式（23）重寫(xiě)為

2.2 算法步驟

綜上，對(duì)所提的陣元幅相誤差和DOA 同時(shí)估計(jì)算法總結(jié)如下：

第1 步：估計(jì)協(xié)方差矩陣，計(jì)算噪聲子空間Un；

第2 步：計(jì) 算 變 換 矩 陣T?(θ)，構(gòu) 造 譜 函 數(shù)Pdet(θ)；

第3 步：譜峰搜索，找到L個(gè)最大譜峰的位置，估計(jì)角度，計(jì)算? 對(duì)應(yīng)的過(guò)渡矩陣Q?(?)；

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

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

3 計(jì)算機(jī)仿真

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

實(shí)驗(yàn)1： 陣元幅相誤差和DOA 同時(shí)估計(jì)算法的有效性驗(yàn)證

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

式中：δm和ηm均為相互獨(dú)立且均勻分布于[-0.5，0.5]中的隨機(jī)變量；σα和σφ分別為δm和ηm的標(biāo)準(zhǔn)偏差，分別設(shè)置為0.1 和15°。入射信號(hào)角度分別為-30°，-10°，40°，對(duì)應(yīng)的SNR 均為10 dB，采樣快拍數(shù)為200。仿真得到的空間譜和幅相誤差估計(jì)值分別如圖2 和表1 所示。

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

圖2 空間譜Fig. 2 Spatial spectrum

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

實(shí)驗(yàn)2： 陣元幅相誤差和DOA 同時(shí)估計(jì)算法的DOA 估計(jì)性能

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

圖3 DOA 估計(jì)值的RMSE 與SNR 的關(guān)系Fig. 3 Relationship between RMSE and SNR of estimated DOA

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

實(shí)驗(yàn)3： 陣元幅相誤差和DOA 同時(shí)估計(jì)算法陣元幅相誤差估計(jì)性能

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

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

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

實(shí)驗(yàn)4： 算法運(yùn)行時(shí)間比較

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

圖5 2 種算法的運(yùn)行時(shí)間對(duì)比Fig. 5 Comparison of running time between two algorithms

4 結(jié)束語(yǔ)

本文通過(guò)設(shè)置少量的輔助陣元，利用陣元幅相誤差矩陣為對(duì)角矩陣的結(jié)構(gòu)特點(diǎn)，并基于子空間類算法的原理，提出了一種陣元幅相誤差和DOA 同時(shí)估計(jì)算法。針對(duì)信號(hào)功率相差較大問(wèn)題，給出了2 種不同的解決策略，并取得了較好的效果。計(jì)算機(jī)仿真結(jié)果表明：所提算法能夠以較高的精度同時(shí)估計(jì)出幅相誤差和目標(biāo)的DOA，尤其在存在功率差異較大信號(hào)時(shí)，能夠顯著提高對(duì)誤差的估計(jì)精度。