基于子空間跟蹤的遞推穩(wěn)健自適應(yīng)波束形成

摘 要:為了解決實(shí)時(shí)系統(tǒng)應(yīng)用中，通用秩信號(hào)模型的穩(wěn)健自適應(yīng)波束形成算法受計(jì)算復(fù)雜度高的限制這一問(wèn)題，基于子空間跟蹤理論，運(yùn)用梯度下降法提出了一種遞推穩(wěn)健自適應(yīng)波束形成算法，有效降低了計(jì)算量，提高了系統(tǒng)性能，改善了陣列輸出的信干噪比。仿真結(jié)果表明，該算法具有快速收斂性與優(yōu)越的穩(wěn)健性。

關(guān)鍵詞:穩(wěn)健自適應(yīng)波束形成; 通用秩信號(hào)模型; 子空間跟蹤; 梯度下降; 計(jì)算復(fù)雜度

中圖法分類號(hào):TN911.7文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1001-3695(2010)06-2260-04

doi:10.3969/j.issn.1001-3695.2010.06.075

Recursive robust adaptive beamforming

WANG Zhiqiao， HUANG Hua

(School of Electrical Engineering Information， Sichuan University， Chengdu 610065， China)

Abstract:The robust adaptive beamforming， based on generalrank signal models， is restricted by the great computational complexity in the realtime systems. Based on subspace tracking theory， this paper proposed a recursive robust adaptive beamforming by gradient descent method.The proposed algorithm decreased computational complexity， developed online implementations， made the mean output array SINR. Simulations results demonstrate that the proposed method can converg rapidly and has excellent robust performance.

Key words:robust adaptive beamforming; generalrank signal models; subspace tracking; gradient descent; computational complexity

0 引言

波束形成是陣列信號(hào)處理的關(guān)鍵技術(shù)之一[1]。傳統(tǒng)自適應(yīng)波束形成算法假定訓(xùn)練序列中不含期望信號(hào)，當(dāng)系統(tǒng)有誤差時(shí)，期望信號(hào)會(huì)被當(dāng)做干擾抑制，系統(tǒng)性能嚴(yán)重下降[1，2]。為使系統(tǒng)存在誤差時(shí)仍能保持良好的性能，出現(xiàn)了大量穩(wěn)健自適應(yīng)波束形成算法，主要有線性約束波束形成算法(LCMV)、對(duì)角加載波束形成算法(LSMI)、特征空間波束形成算法。目前，基于最差情況下優(yōu)化和顯式系統(tǒng)誤差建模的穩(wěn)健自適應(yīng)波束形成算法得到了廣泛深入的研究[2~12]。文獻(xiàn)[5]中，考慮信道散射效應(yīng)和信號(hào)協(xié)方差陣與數(shù)據(jù)協(xié)方差陣存在誤差時(shí)，推導(dǎo)出基于通用秩信號(hào)模型的穩(wěn)健自適應(yīng)波束形成算法。

上述穩(wěn)健自適應(yīng)波束形成算法都是批處理算法，即每次數(shù)據(jù)快拍到來(lái)時(shí)，需計(jì)算波束形成器權(quán)向量，計(jì)算時(shí)要對(duì)數(shù)據(jù)協(xié)方差陣求逆或特征值分解，其計(jì)算復(fù)雜度為O(M3)。文獻(xiàn)[13]在LCMV波束形成器基礎(chǔ)上降秩，得到低計(jì)算成本的迭代波束形成器;文獻(xiàn)[14]用約束卡爾曼濾波器降低了文獻(xiàn)[4]中算法的復(fù)雜度。這兩種算法只適用于點(diǎn)信源(信號(hào)協(xié)方差陣秩為1)。為了得到一般信源情況下低計(jì)算復(fù)雜度的穩(wěn)健自適應(yīng)波束形成算法，本文運(yùn)用信號(hào)子空間跟蹤理論，推導(dǎo)出一種梯度下降的遞推穩(wěn)健自適應(yīng)波束形成算法。

1 理論基礎(chǔ)

1.1 信號(hào)模型

窄帶波束形成器輸出

y(n)=wHx(n) (1)

其中:n是時(shí)間下標(biāo);x(n)是M×1維陣列觀察數(shù)據(jù)的復(fù)向量;w是M×1維波束形成器的復(fù)值權(quán)向量，最優(yōu)權(quán)向量由最大化信干噪比(SINR)獲得;M是陣元個(gè)數(shù);(#8226;)T與(#8226;)H分別代表轉(zhuǎn)置與共軛轉(zhuǎn)置。陣列觀察數(shù)據(jù)向量為

x(n)=s(n)+i(n)+n(n) (2)

其中:s(n)、i(n)與n(n)分別是期望信號(hào)、干擾信號(hào)與噪聲。假設(shè)期望信號(hào)、干擾信號(hào)與噪聲互不相關(guān)，且期望信號(hào)與干擾信號(hào)為平穩(wěn)信號(hào)，噪聲為高斯白噪聲。

最大化信干噪比(SINR)為

SINR=wHRswwHRi+nw (3)

其中:

Rs=E{s(n)sH(n)} (4)

Ri+n=E{(i(n)+n(n))(i(n)+n(n))H} (5)

分別是M×M維的信號(hào)協(xié)方差陣與干擾加噪聲協(xié)方差陣，E{#8226;}是統(tǒng)計(jì)期望。通常，矩陣Rs有任意的秩，即

1≤rank{Rs}≤M

在點(diǎn)信源情況下(Rs秩為1)，有

s(n)=s(n)a(θs)

Rs=σ2sa(θs)a(θs)H

其中:a(θs)是期望信號(hào)導(dǎo)向向量;θs是期望信號(hào)相對(duì)于陣列法線方向的入射角度;s(n)是零均值信號(hào)，σ2s=E{|s(n)|2}是s(n)的方差。式(3)可以簡(jiǎn)化為

SINR=σ2s|wHa(θs)|2wHRi+nw(6)

在滿秩信源不相干散射情況下，Rs有如下形式

Rs=σ2s∫π/2-π/2ρ(θ)a(θ)aH(θ)dθ (7)

其中:ρ(θ)是歸一化角度功率譜密度[5]。

1.2 傳統(tǒng)自適應(yīng)波束形成器

通用秩信號(hào)模型下，期望信號(hào)方向無(wú)失真的最小方差波束形成器

minw wH Ri+nws.t.wHRsw=1 (8)

用拉格朗日乘子法解出

wopt=P{R-1i+nRs}(9)

其中:P{#8226;}是生成矩陣主特征向量(最大特征值對(duì)應(yīng)的特征向量)的算子。在點(diǎn)信源的情況下，即為MVDR波束形成器

wopt=R-1i+na(θs)a(θs)HR-1i+na(θs)(10)

當(dāng)沒有系統(tǒng)誤差時(shí)，用數(shù)據(jù)協(xié)方差陣R=E{x(n)xH(n)}=Ri+n+Rs替代式(10)的Ri+n，結(jié)果不變[1]。

實(shí)際應(yīng)用中，不能獲得真實(shí)數(shù)據(jù)協(xié)方差陣R，只能得到其估計(jì)值，常用的無(wú)偏估計(jì)為

R^=1K∑kn=1x(n)xH(n)(11)

其中:K為快拍數(shù)，用估計(jì)值R^代替R或Ri+n，即為采樣矩陣求逆波束形成器(SMI)

wSMI=P{R^-1Rs}(12)

當(dāng)存在期望信號(hào)導(dǎo)向向量與角度功率譜密度誤差時(shí)，上述自適應(yīng)波束形成器性能嚴(yán)重下降。廣泛使用的穩(wěn)健技術(shù)是對(duì)角加載波束形成算法(LSMI)

wLSMI=P{(R^+γI)-1Rs}(13)

其中，最優(yōu)加載量通常通過(guò)數(shù)值模擬的方式選擇，沒有更優(yōu)的選擇方法。

1.3 通用秩信號(hào)模型的穩(wěn)健自適應(yīng)波束形成算法[5]

假定信號(hào)協(xié)方差陣與數(shù)據(jù)協(xié)方差陣都存在有限誤差，誤差矩陣都為Hermitian矩陣，即

Rs～=Rs+Δ1 其中，‖Δ1‖≤ξ (14)

R～=R+Δ2 其中，‖Δ2‖≤γ (15)

其中:Rs～與R～分別是實(shí)際信號(hào)協(xié)方差陣與數(shù)據(jù)協(xié)方差陣，Rs與R是對(duì)應(yīng)的假定協(xié)方差陣。實(shí)際情況下，輸出信干噪比可表示為

SINR=wHR～swwHR～w(16)

最大化SINR優(yōu)化模型中考慮誤差

minw max‖Δ2‖≤γ wH(R+Δ2)w

s.t. min‖Δ1‖≤ξ wH(Rs+Δ1)w=1 (17)

得到最差情況下存在解析解的最優(yōu)穩(wěn)健自適應(yīng)波束形成器

wSGLW03=P{(R+γI)-1(Rs-ξI)}(18)

其中:wSGLW03記做文獻(xiàn)[5]中算法得到的最優(yōu)權(quán)向量，γ與ξ由環(huán)境參數(shù)估計(jì)。

1.4 信號(hào)子空間

假定x是協(xié)方差矩陣為C=E{xxH}的復(fù)值向量隨機(jī)過(guò)程，標(biāo)量函數(shù)

J(W)=E{‖x-WWHx‖2}=

tr(C)-2tr(wHCW)+tr(WHCWWHW)(19)

其中:變量W為M×r復(fù)值矩陣，不失一般性，假設(shè)其秩為r。由文獻(xiàn)[15]中的定理可知，當(dāng)且僅當(dāng)W=UrQ時(shí)，W是標(biāo)量函數(shù)J(W)的平穩(wěn)點(diǎn)。其中:Ur是C的任意r個(gè)獨(dú)立的特征向量，Q是任意的r×r維酉矩陣，而且在每個(gè)平穩(wěn)點(diǎn)處，標(biāo)量函數(shù)J(W)的值都等于其余特征向量對(duì)應(yīng)的特征值之和。因此，當(dāng)且僅當(dāng)Ur包含C的r個(gè)主特征向量(矩陣特征按從大到小排序后，前r個(gè)特征值對(duì)應(yīng)的特征向量)時(shí)，標(biāo)量函數(shù)達(dá)到最小值，即當(dāng)標(biāo)量函數(shù)達(dá)到最小值時(shí)，矩陣變量W收斂為C的r個(gè)主特征向量構(gòu)成子空間的某個(gè)正交基上。若已知數(shù)據(jù)向量是由r個(gè)信號(hào)加噪聲構(gòu)成的，最小化標(biāo)量函數(shù)將會(huì)得到信號(hào)子空間的一個(gè)正交基。需要注意的是，由于目標(biāo)函數(shù)的酉不變性，這個(gè)正交基不是惟一的，但信號(hào)子空間的投影矩陣WWH是惟一的。當(dāng)然，在r=1時(shí)，最小化目標(biāo)函數(shù)將會(huì)得到C的歸一化主特征向量。因此，可以通過(guò)迭代的方式最小化目標(biāo)函數(shù)，得到矩陣C的歸一化的主特征向量。

2 遞推穩(wěn)健自適應(yīng)波束形成算法

2.1 算法推導(dǎo)

存在解析解的穩(wěn)健自適應(yīng)波束形成算法 ，計(jì)算量集中在計(jì)算復(fù)雜度都為O(M3)的矩陣求逆與特征值分解，實(shí)際應(yīng)用中并不高效。本文通過(guò)矩陣求逆引理與子空間跟蹤理論推導(dǎo)出一種高效算法——遞推穩(wěn)健自適應(yīng)波束形成算法。

推導(dǎo)時(shí)借助文獻(xiàn)[5]的引理:對(duì)于任意的M×M Hermitian矩陣X與任意M×M滿秩Hermitian矩陣Y有如下關(guān)系成立

P{XY}=Y-1/2P{Y1/2XY1/2} (20)

顯然，+γI是Hermitian矩陣，雖然在點(diǎn)信源等情況下Hermitian矩陣Rs非滿秩，但對(duì)任意ξ>0，矩陣Rs-ξI卻總是滿秩的。需要注意的是，要使矩陣Rs-ξI非負(fù)定，ξ須小于Rs的最大特征值，因?yàn)楫?dāng)矩陣負(fù)定時(shí)，式(17)的約束條件不成立。在這種情況下，穩(wěn)健的自適應(yīng)波束形成器式(18)(用代替R)為

wSGLW03=(Rs-ξI)-1/2P{(Rs-ξI)1/2(+γI)-1(Rs-ξI)1/2}=

(Rs-ξI)-1/2P{G-1}

(21)

其中，矩陣G定義為

G=(Rs-ξI)-1/2(+γI)(Rs-ξI)-1/2(22)

為了推導(dǎo)遞推算法，用長(zhǎng)度為N的固定窗，將第n步對(duì)角加載的數(shù)據(jù)協(xié)方差陣(n)展開為

DL(n)=DL(n-1)+1Nx(n)xH(n)-1Nx(n-N)xH(n-N)(23)

其中:DL(n)=(n)+γI。即為秩2更新，γI可以選做DL(n)的初始值，可將G(n)同樣展開表示為

G(n)=G(n-1)+(n)H(n)-(n-N)H(n-N)(24)

其中:變換了的數(shù)據(jù)快拍表示為

(i)=1N(Rs-ξI)-1/2x(i)(25)

(i)=1N[(Rs-ξI)-1/2]Hx(i)(26)

之所以用(i)與(i)兩個(gè)不同的符號(hào)表示，是因?yàn)橐话闱闆r下，矩陣Rs-ξI是不定的，即對(duì)于任意的ξ>0且小于Rs的最大特征值，矩陣Rs-ξI同時(shí)存在正負(fù)特征值，使得其平方根矩陣(Rs-ξI)1/2一定不是Hermitian陣，因此用兩個(gè)符號(hào)以示區(qū)別。

用矩陣求逆引理[16]得到

H-1(n-1)=G-1(n-1)+G-1(n-1)(n-N)H(n-N)G-1(n-1)1-H(n-N)G-1(n-1)(n-N)(27)

G-1(n)=H-1(n-1)-H-1(n-1)(n)H(n)H-1(n-1)1+H(n)H-1(n-1)(n)(28)

其中，中間變量H(n-1)=G(n-1)-(n-N)H(n-N)，G-1(n)的初始值為1γ(Rs-ξI)。

在1.4節(jié)中，在當(dāng)變量W為矩陣C的主特征向量時(shí)，無(wú)約束目標(biāo)函數(shù)J(W)取得最小值，最直接的迭代求解法是梯度下降類算法。求標(biāo)量函數(shù)J(W)關(guān)于W的梯度

J(W)=(-2C+CWWH+WWHC)W(29)

因此，依據(jù)梯度下降法的主特征向量更新有

W(n)=W(n-1)-μ[-2C+CW(n-1)WH(n-1)+W(n-1)WH(n-1)C]W(n-1) (30)

用wpro(n)、wpro(n-1)與G-1(n)分別代替上式中的W(n)、W(n-1)與C，得到更新波束形成器權(quán)向量的遞推式為

wpro(n)=wpro(n-1)-μ[-2G-1(n)+G-1(n)wpro(n-1)wHpro(n-1)+

wpro(n-1)wHpro(n-1)G

-1(n)]wpro(n-1)

(31)

其中:wpro(n)為本文算法的權(quán)向量，μ為步長(zhǎng)參數(shù)。

2.2 算法描述及計(jì)算復(fù)雜度分析

將遞推穩(wěn)健自適應(yīng)波束形成算法總結(jié)如下:

initialize

G-1(0)=1γ(Rs-ξI)

wpro(0)=(Rs-ξI)-1/2P{G-1(0)}

while(new snapshots x(n))

Compute transformed snapshots (2M2-2M)

計(jì)算式(25)(26)

if n

G-1(n)=G-1(n-1)-G-1(n-1)(n)H(n)G-1(n-1)1+H(n)G-1(n-1)(n)

else(10M2)

計(jì)算式(27)(28)

Compute new weighted vector (7M2-4M)

計(jì)算式(31)

end

其中:2M2-2M、5M2、10M2、7M2-4M表示對(duì)應(yīng)計(jì)算需要的浮點(diǎn)數(shù)操作的次數(shù)[16]，每次浮點(diǎn)數(shù)操作為一次浮點(diǎn)數(shù)加法或浮點(diǎn)數(shù)乘法。一般情況下，總的浮點(diǎn)數(shù)操作為19M2-6M，計(jì)算復(fù)雜度為O(M2)，大大降低了文獻(xiàn)[5]算法的計(jì)算量。

3 仿真實(shí)驗(yàn)結(jié)果

算法推導(dǎo)過(guò)程中，并未限定陣列結(jié)構(gòu)，因此遞推穩(wěn)健自適應(yīng)波束形成算法適用于任何陣列，但為分析簡(jiǎn)便起見，實(shí)驗(yàn)考慮陣元數(shù)M=20，陣元間隔為半波長(zhǎng)的標(biāo)準(zhǔn)均勻線性陣列[1]， 每一點(diǎn)都是100次的均值。假定只有一個(gè)期望信號(hào)與一個(gè)干擾信號(hào)，訓(xùn)練序列中包含期望信號(hào)，干擾噪聲比INR=20 dB。固定窗長(zhǎng)度N=20，數(shù)據(jù)快拍總數(shù)K=500。

實(shí)驗(yàn)比較以下算法的輸出信干噪比:采樣矩陣求逆(SMI)，即式(12);對(duì)角加載的采樣矩陣求逆(LSMI)，即式(13);文獻(xiàn)[5]的算法，即式(18)，記做SGLW03;本文的遞推算法。

同時(shí)在每個(gè)圖中顯示最優(yōu)信干噪比，即式(9)(Ri+n與Rs精確已知的情況)。

LSMI、SGLW03與本文算法的對(duì)角加載量為γ=30，按照加載量噪聲比為10 dB的原則選取[1]。

實(shí)驗(yàn)1 考慮點(diǎn)信源且系統(tǒng)沒有誤差。

設(shè)參數(shù)ξ=16，期望信號(hào)與干擾的入射角度為相對(duì)于陣列法線方向30°與-30°，圖1顯示了在固定信噪比SNR=0 dB時(shí)，本文算法與SGLW03算法的輸出信干噪比。可以看出，本文算法與SGLW03算法有相似的性能。

實(shí)驗(yàn)2 考慮點(diǎn)信源且期望信號(hào)存在DoA估計(jì)誤差。

設(shè)參數(shù)ξ=16，實(shí)際的期望信號(hào)與干擾的入射角度為30°與-30°，假定的期望信號(hào)入射角度為33°，即估計(jì)誤差為3°。本文算法使用窗長(zhǎng)N=20。圖2給出了SNR=0時(shí)平均輸出信干噪比SINR隨快拍數(shù)的變化。圖3表示快拍數(shù)為200時(shí)平均信干噪比SINR隨信噪比SNR的變化。從兩圖中看出，即使DoA失配時(shí)，本文的算法和SGLW03算法性能相近，不同快拍和信噪比下都有較高的輸出信干噪比，而SMI隨著快拍增加或信噪比提高性能嚴(yán)重下降，高信噪LSMI算法的性能也有所下降。

實(shí)驗(yàn)3 考慮滿秩信源下信號(hào)協(xié)方差陣下有偏差。

設(shè)參數(shù)ξ=9，假設(shè)實(shí)際的期望信號(hào)與干擾都存在局部不相干散射[5]，其實(shí)際的角度功率譜密度分別服從中心角度為30°與-30°的高斯與均勻分布，角度擴(kuò)散為4°，相當(dāng)于高斯分布的標(biāo)準(zhǔn)差為2°，均勻分布的范圍為[-32°，-28°]。而假定的信號(hào)協(xié)方差陣存在偏差，假設(shè)其估計(jì)的角度功率譜密度也是服從高斯分布，中心角度為32°，標(biāo)準(zhǔn)差為3°，即估計(jì)的中心角度與標(biāo)準(zhǔn)差都有偏差。圖4顯示當(dāng)SNR=0時(shí)，各算法輸出信干噪比隨快拍數(shù)的變化。圖5顯示了快拍數(shù)為200時(shí)輸出信干噪比隨信噪比的變化，其中算法假設(shè)窗口長(zhǎng)度為N=500，即相當(dāng)于數(shù)據(jù)協(xié)方差陣采用秩1更新。同實(shí)驗(yàn)2有相似的結(jié)果，即使是在信號(hào)協(xié)方差陣滿秩的情況下，本文的遞推算法和SGLW03算法有同樣優(yōu)越的穩(wěn)健性能。

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

利用信號(hào)子空間跟蹤理論，針對(duì)通用秩信號(hào)模型的穩(wěn)健自適應(yīng)波束形成算法計(jì)算復(fù)雜度高的問(wèn)題，本文提出了一種遞推穩(wěn)健自適應(yīng)波束形成算法，并分析了該算法每次更新權(quán)向量只需19M2-6M次浮點(diǎn)數(shù)操作，即計(jì)算復(fù)雜度為O(M2)。從數(shù)值仿真中可以看到，本文算法在沒有系統(tǒng)誤差時(shí)，能快速收斂于文獻(xiàn)[5]的算法，而且無(wú)論是點(diǎn)信源存在DoA估計(jì)誤差，還是滿秩信源存在角度功率譜密度參數(shù)估計(jì)誤差，本文算法都有優(yōu)越的穩(wěn)健性能。

參考文獻(xiàn):

[1]VAN TREES H L. Optimum array processing， detection， estimation， and modulation theory， vol 4[M]. New York: Wiley， 2002.

[2]LI J， STOICA P. Robust adaptive beamforming[M]. Hoboken， NJ: Wiley， 2006.

[3]GERSHMAN A B. Robust adaptive beamforming in sensor arrays[J]. AEUInnt Journal of Electronics and Communications， 1999， 53(6): 305-314.

[4]VOROBYOV S A， GERSHMAN A B， LUO Z Q. Robust adaptive beamforming using worstcase performance optimization: a solutionto the signal mismatch problem[J]. IEEE Trans on Signal Processing， 2003， 51 (2): 313-324.

[5]SHAHBAZPANAHI S， GERSHMAN A B， LUO Z Q. Robust adaptive beamforming for generalrank signal models[J]. IEEE Trans on Signal Processing， 2003， 51(9): 2257-2269.

[6]KIM S J， MAGNANI A， MUTAPCIC A， et al. Robust beamforming via worstcase SINR maximization[J]. IEEE Trans. on Signal Processing， 2008， 56(4): 1539-1547.

[7]YU Z L， SER W， ER M H， et al. Robust response control for adaptive beamformers against arbitrary array imperfections[J]. Signal Processing， 2008， 88(9): 2234-2245.

[8]YU Z L， ER M H， SER W. A novel adaptive beamformer based on semidefinite programming (SDP) with magnitude response constraints[J]. IEEE Trans on Antennas and Propagation， 2008， 56(5): 1297-1307.

[9]ZHANG Y， NG B P， WAN Q. Sidelobe suppression for adaptive beamforming with sparse constraint on beam pattern[J]. IEE Electronics Letters， 2008，44(10): 615-616.

[10] LI J， XIE Y， STOICA P， et al. Beampattern synthesis via a matrix approach for signal power estimation[J]. IEEE Trans on Signal Processing， 2007， 55(12): 5643-5657.

[11]CHEN C， VAIDYANATHAN P P. Quadratically constrained beamforming robust against directionofarrival mismatch[J]. IEEE Trans on Signal Processing， 2007， 55(8): 4139-4150.

[12]SELN Y， ABRAHAMSSON R， STOICA P. Automatic robust adaptive beamforming via ridge regression[J]. Signal Processing， 2008， 88(1): 33-49.

[13]De LAMARE R C. Adaptive reducedrank LCMV beamforming algorithms based on joint iterative optimisation of filters[J]. IEE Electronics Letters， 2008， 44(9): 565-566.

[14]ElKEY A， KIRUBARAJAN T， GERSHMAN A B. Robust adaptive beamforming based on Kalman filter[J]. IEEE Trans on Signal Processing， 2005， 53(8):3032-3041.

[15]YANG B. Projection approximation subspace tracking[J]. IEEE Trans on Signal Processing， 1995， 43(1): 95-107.

[16]GOLUB G H， VAN LOAN C F. Matrix computations[M]. 3rd ed. Baltimore: Johns Hopkins University Press， 1996.