基于陣列劃分的近場DOA估計(jì)算法

陳光輝，曾孝平，焦 爽

（1.重慶大學(xué)微電子與通信工程學(xué)院，重慶 400044；2.中國大唐集團(tuán)科學(xué)技術(shù)研究院有限公司中南電力試驗(yàn)研究院，河南鄭州 450000）

1 引言

信源定位是陣列信號(hào)處理領(lǐng)域中的一個(gè)關(guān)鍵技術(shù)，其在雷達(dá)、聲納、電子對(duì)抗、醫(yī)療電子等領(lǐng)域中有著廣泛的應(yīng)用［1～3］.根據(jù)信源與接收陣列之間的距離，信源定位可以分為遠(yuǎn)場和近場信源定位［4，5］.對(duì)于遠(yuǎn)場信源定位，信源的傳輸模型被假設(shè)為平面波，僅用到達(dá)方向（Direction Of Arrival，DOA）就可以完成信源定位的任務(wù).在過去的幾十年中，學(xué)者們提出了很多針對(duì)遠(yuǎn)場DOA估計(jì)的經(jīng)典算法.例如最大似然法［6］，多重信號(hào)分類（Multiple Signal Classification，MUSIC）算法［7，8］，旋轉(zhuǎn)不變技術(shù)（Estimation of Signal Parameters via Rotational Invariance Technique，ESPRIT）方法［9，10］.對(duì)于近場信源定位，信源的傳輸模型被假設(shè)為球面波，需要估計(jì)信源的DOA和距離才可以完成信源定位的任務(wù).針對(duì)近場DOA估計(jì)，學(xué)者們提出了二維 近 場MUSIC算 法［11，12］，ESPRIT高 階 近 場 估 計(jì) 算法［13，14］，以及一些在MUSIC和ESPRIT算法基礎(chǔ)上改進(jìn)的算法［15～19］.但這些算法在低信噪比（Signal-to-Noise Ratio，SNR）下，因測量數(shù)據(jù)的微小變化，將出現(xiàn)DOA估計(jì)的精度和分辨率嚴(yán)重下降的問題.針對(duì)這一問題，本文提出了一種基于陣列劃分的近場DOA估計(jì)算法.

該算法首先利用四階累積量和陣列劃分構(gòu)造只包含DOA信息的托普利茲矩陣，并利用一維MUSIC算法對(duì)托普利茲矩陣進(jìn)行特征值分解和構(gòu)造空間譜.然后，由于MUSIC算法的空間譜在譜峰處是一個(gè)斷點(diǎn)、且趨近較大值.因此，為提高DOA估計(jì)分辨率，該算法求解空間譜的一階導(dǎo)數(shù)來構(gòu)造一個(gè)新“空間譜”.最后，利用協(xié)方差矩陣建立二維近場MUSIC空間譜，并將第二步得到的DOA估計(jì)值代入二維近場MUSIC空間譜，然后利用本文構(gòu)造的新“空間譜”逐個(gè)估計(jì)近場信源的距離.此外，本文理論分析了算法的計(jì)算復(fù)雜度，并利用仿真實(shí)驗(yàn)驗(yàn)證了該算法在低信噪比下有效的提高了DOA的估計(jì)精度.

2 近場模型

近場模型如圖1所示.其中，2M+1表示接收陣列的總天線個(gè)數(shù)，N表示接收陣列劃分的第一個(gè)子陣列的天線個(gè)數(shù)，信源的數(shù)量為K，噪聲為高斯噪聲.第m個(gè)天線的輸出表示為：

圖1 近場模型的接收陣列

其中，sk(t)表示第k個(gè)近場窄帶信號(hào)源，nm(t)表示加性高斯噪聲，τmk表示第m個(gè)接收天線與第0個(gè)接收天線之間的傳輸延遲.定義rm表示信源到第m個(gè)接收天線的距離.由余弦定理可知：

由圖1的近場模型可知，信源傳輸模型是球面波，因此，τmk被表示如下：

rm被帶入到τmk，可得：

其中，λ、d、θk、rk分別表示信源的波長、天線間距、DOA和距離.為簡化式（4），f(d/rk)被定義為：

式（5）在變量為0的點(diǎn)被麥克勞林展開，f(d/rk)被重寫為：

其中，f1(0)、f2(0)、f3(0)分別表示f(d/rk)在零點(diǎn)處的一階、二階、三階導(dǎo)數(shù).通常d/rk的值比較小，故可以被忽略［20］.因此，τmk被重寫為：

其中，wk和φk分別表示為：

由式（8）和式（9）可知，wk僅與近場信源的DOA有關(guān)，而φk既與近場信源的DOA有關(guān)，也與近場信源的距離有關(guān).

3 本文算法的原理

3.1 DOA估計(jì)

對(duì)于近場信源的DOA估計(jì)，本文提出的算法主要包括：陣元選擇、構(gòu)造托普利茲矩陣和新“空間譜”.陣元選擇消除四階累積量中與距離有關(guān)的信息，避免近場信源估計(jì)中的二維譜峰搜索，降低算法的計(jì)算量；構(gòu)造托普利茲矩陣可以使矩陣經(jīng)過特征值分解之后其特征向量相互正交.在這種情況下，特征向量組成的矩陣僅包含一個(gè)方向向量相關(guān)的信息；構(gòu)造新“空間譜”的目的是利用MUSIC算法的空間譜在譜峰處是一個(gè)斷點(diǎn)、且趨近較大值的性質(zhì)來提高信源估計(jì)的分辨率.

3.1.1 陣元選擇

其次，近場模型的接收陣列被劃分為三個(gè)子陣列：[-N，N]、[N+1，M]、[-M，-N-1]，且三個(gè)子陣列互不重疊.此外，m，n，p，q的取值是為了使四階累積量中φk的系數(shù)為0，這樣四階累積量只與近場信源的DOA有關(guān)，可以有效的避免近場信源估計(jì)中的二維譜峰搜索，降低算法的計(jì)算量.具體的陣元選擇方法如下：

令m∈[-M，-N-1]，n∈[N+1，M]，p∈[-N，0]，q∈[0，N]，且m=-n，p=-q，可以得到一個(gè)具有(MN)(N+1)個(gè)元素的向量c1.其中，c1的第(M-m)(N+1)+N+1-p個(gè)元素表示為：

令m∈[N+1，M]，n∈[-M，-N-1]，p∈[0，N]，q∈[-N，0]，且m=-n，p=-q，可以得到一個(gè)具有(MN)(N+1)個(gè)元素的向量c2.其中，c2的第(m-N-1)(N+1)+p+1個(gè)元素表示為：

令p∈[-N，N]，q∈[N，-N]，m=n=0，p=-q，可以得到一個(gè)具有2N+1個(gè)元素的向量c3，其中c3的第p個(gè)元素表示為：

令p=-N-1，q=N+1，m=-M，n=-N，可以得到一個(gè)標(biāo)量c4，其被表示為：

令p=N+1，q=-N-1，m=M，n=N，可以得到一個(gè)標(biāo)量c5，其被表示為：

由式（12）和（13）可知，c1與c2的倒序共軛對(duì)稱.由式（14）可知，c3的2N+1個(gè)元素共軛對(duì)稱.由式（15）和（16）可知，c4和c5也共軛.c1，c2，c3，c4，c5的共軛對(duì)稱關(guān)系是構(gòu)造托普利茲矩陣的基礎(chǔ).

3.1.2 托普利茲矩陣構(gòu)造

構(gòu)造c1，c2，c3，c4，c5之后，需要構(gòu)造一個(gè)托普利茲矩陣.因?yàn)橥衅绽澗仃嚱?jīng)過特征值分解之后其特征向量相互正交.在這種情況下，特征向量組成的矩陣僅包含一個(gè)與方向向量相關(guān)的信息.具體的托普利茲矩陣構(gòu)造過程見下文.

首先，一個(gè)(2(M-N+1)(N+1)+1)維的向量c被定義為：

其中，T表示矩陣轉(zhuǎn)置.

其次，定義一個(gè)矩陣CH，把向量c中第((M-N+1)(N+1)+2-m)到(2(M-N+1)(N+1)+2-m)的元素組成CH的第m列.其中，CH的具體表達(dá)式如下所示：

由c1，c2，c3，c4，c5的共軛對(duì)稱關(guān)系知，托普利茲矩陣CH的 維 度 是((M-N+1)(N+1)+1)×((M-N+1)(N+1)+1).為簡化式（18），定義Q=(M-N)(N+1).此外，式（18）可以被表示為矩陣形式，如下所示：

其中，H表示矩陣共軛轉(zhuǎn)置.A(θ)和CS公式如下所示：

基于式（19）～（21），A(θ)可以看作是K個(gè)遠(yuǎn)場信源的虛擬方向矩陣，且其對(duì)應(yīng)(M-N+1)(N+1)+1天線（虛擬天線）.另外，因?yàn)閟k(t)統(tǒng)計(jì)獨(dú)立，且CS和A(θ)均為滿秩矩陣，因此CH也是滿秩矩陣.所以，可以采用一維遠(yuǎn)場MUSIC算法估計(jì)近場信源的DOA.其中，基于一維遠(yuǎn)場MUSIC算法，a(θ)與Un之間的關(guān)系可以描述為［8］：

其中，Un表示由CH矩陣特征值分解之后得到的噪聲空間.因此，MUSIC算法的空間譜函數(shù)定義為：

3.1.3 構(gòu)造新“空間譜”

MUSIC算法空間譜的物理意義表示為：P(θ)不是真實(shí)的譜，可以理解為掃描向量a(θ)到噪聲子空間Un的歐幾里德距離的倒數(shù).由MUSIC算法空間譜的物理意義和aH(θ)Un=0可知，P(θ)函數(shù)在信源的DOA處存在尖峰（也可以認(rèn)為是斷點(diǎn)）.利用此性質(zhì)，對(duì)P(θ)求一階導(dǎo)數(shù)，銳化譜峰，可以提高DOA估計(jì)的分辨率.因此，P(θ)的一階導(dǎo)數(shù)d(θ)定義為：

為方便計(jì)算，選擇合適的步長h可以用一階差分函數(shù)y(θ)代替一階導(dǎo)數(shù)d(θ)，其中一階差分函數(shù)y(θ)被定義為：

其中，式（25）即是本文算法定義的新“空間譜”的表達(dá)式.

3.2 距離估計(jì)

為了估計(jì)近場信源的距離，本文算法又構(gòu)造一個(gè)協(xié)方差矩陣R=E{y(t)y*(t)}，對(duì)其特征值分解可以構(gòu)造一個(gè)二維空間譜，具體的表達(dá)式如下：

其中，U0表示協(xié)方差矩陣R=E{y(t)y*(t)}經(jīng)過特征值分解之后的噪聲子空間.在3.1節(jié)中獲得的DOA估計(jì)值{θk，k=1，…，K}被代入式（26）中，近場信源的距離可以被一一估計(jì).因此，式（26）可以被重寫為：

在估計(jì)近場信源的距離時(shí)，同樣利用本文構(gòu)造的新“空間譜”估計(jì)近場信源的距離，公式如下：

3.3 本文算法的步驟

本文算法的步驟如表1所示：

表1 本文算法的步驟

4 仿真分析

4.1 DOA估計(jì)分辨率理論分析和仿真

本節(jié)首先分析本文算法定義的新“空間譜”如何確定近場信源的DOA，然后分析本文算法的DOA估計(jì)分辨率.假設(shè)θ0是近場信源的DOA，由式（22），（23）和MUSIC算法的空間譜的物理意義可知，P(θ0)是一個(gè)較 大 的 正 值，P(θ0-h)，P(θ0+h)，P(θ0+2h)遠(yuǎn) 小 于P(θ0).因此，可以計(jì)算出y(θ0-h)是一個(gè)較大的正值，y(θ0+h)是一個(gè)較大的負(fù)值.所以，根據(jù)式（25）中y(θ)函數(shù)是否先有一個(gè)較大的正值，后有一個(gè)較大的負(fù)值來確定近場信源的DOA.此外，本文算法通過對(duì)空間譜函數(shù)求一階導(dǎo)數(shù)，一定會(huì)先產(chǎn)生一個(gè)極大值，然后產(chǎn)生一個(gè)極小值.這反映到譜峰上是先產(chǎn)生一個(gè)正的譜峰，后又產(chǎn)生一個(gè)負(fù)的譜峰.由于本文算法在求空間譜一階導(dǎo)數(shù)的過程中，采用從-90°到90°的遞增順序，所以在譜峰搜索時(shí)會(huì)舍棄搜索到的負(fù)譜峰，這樣即不會(huì)對(duì)DOA估計(jì)結(jié)果造成影響，也不會(huì)造成DOA估計(jì)分辨率下降.當(dāng)θ1和θ2是兩個(gè)近場信源的DOA時(shí)，MUSIC算法的DOA估計(jì)分辨率定義為［7］：

其中，

基于式（29），本文算法的DOA估計(jì)分辨率為：

由式（31）可知，本文算法DOA估計(jì)分辨率近似等于MUSIC算法的一階導(dǎo)數(shù).由MUSIC算法的空間譜函數(shù)的物理意義可知，空間譜函數(shù)在DOA處是斷點(diǎn)，且其在DOA處的值遠(yuǎn)大于相鄰角度的值.因此，本文算法DOA估 計(jì) 分 辨 率φ(θ)高 于MUSIC的 估 計(jì) 分 辨 率φmusic(θ).此外，本文又通過實(shí)驗(yàn)驗(yàn)證了上述分析，實(shí)驗(yàn)結(jié)果如圖2和圖3所示.

由圖2（a）可知，在DOA=（10°，13°），SNR=0 dB條件下，傳統(tǒng)的MUSIC方法未估計(jì)出兩個(gè)DOA，僅僅在角度等于10°處估計(jì)出一個(gè)譜峰，且不明顯.通過對(duì)5°～20°范圍的空間譜放大處理，顯示在圖2（a）的子圖中，發(fā)現(xiàn)在角度等于13°處的空間譜的譜峰也不明顯.但由圖2（b）可知，在DOA=（10°，13°），SNR=0 dB條件下，本文算法可以在角度等于10°和13°處估計(jì)出兩個(gè)DOA，且空間譜的譜峰明顯.因此，本文算法在SNR=0 dB條件下，可以有效的提高DOA的估計(jì)分辨率.

由圖3（a）可知，在DOA=（10°，13°），SNR=20 dB條件下，傳統(tǒng)的MUSIC方法估計(jì)出兩個(gè)DOA，且空間譜的譜峰較明顯.由圖3（b）可知，在DOA=（10°，13°），SNR=20 dB條件下，本文算法也可估計(jì)出兩個(gè)DOA，且空間譜的譜峰更明顯.但由于MUSIC算法在高SNR條件下，DOA估計(jì)分辨率已經(jīng)較高，所以本文算法在SNR=20 dB條件下，DOA估計(jì)分辨率和MUSIC算法是相似的.最終，由圖2和圖3可知，本文算法提高了低信噪比下的DOA估計(jì)分辨率.

圖2 DOA=(10°,13°),信噪比SNR=0 dB的對(duì)比圖

圖3 DOA=(10°,13°),信噪比SNR=20 dB的對(duì)比圖

4.2 DOA和距離估計(jì)分辨率仿真

仿真參數(shù)如表2所示，噪聲為高斯白噪聲，分別利用3個(gè)不同的近場信源對(duì)比了本文算法和文獻(xiàn)［19］中的算法.其中，3個(gè)近場信源分別由K1、K2、K3表示，K1是（［0°，3λ］，［1°，4λ］），K2是（［10°，3λ］，［13°，4λ］），K3是（［20°，3λ］，［25°，4λ］），DOA角度差Δθ分別為1°、3°、5°.DOA和距離估計(jì)結(jié)果如圖4～6所示.

②Stufflebeam，D．L．，“The CIPP Model for Evaluation”，International Handbook of Educational Evaluation，Springer Netherlands，2003．

表2 仿真參數(shù)表

在圖4中，當(dāng)Δθ為1°和3°時(shí)，文獻(xiàn)［19］的算法只估計(jì)出一個(gè)DOA.當(dāng)Δθ為5°時(shí)，文獻(xiàn)［19］算法估計(jì)兩個(gè)DOA.在圖5中，當(dāng)Δθ為1°時(shí)，本文算法估計(jì)出一個(gè)DOA，當(dāng)Δθ為3°和5°時(shí)，本文算法估計(jì)兩個(gè)DOA，且空間譜的譜峰較明顯.原因是當(dāng)信源的DOA相近時(shí)，信源之間的空間相關(guān)性較強(qiáng).但本文算法利用空間譜函數(shù)在DOA處是斷點(diǎn)并具有較大值的物理特性和空間譜函數(shù)的一階導(dǎo)數(shù)構(gòu)造一個(gè)新“空間譜”來減少DOA處的空間相關(guān)性，因此提高了DOA估計(jì)分辨率.

圖4 文獻(xiàn)[19]中DOA估計(jì)結(jié)果圖

圖5 本文算法的DOA估計(jì)結(jié)果圖

圖6（a）～（c）是文獻(xiàn)［19］的距離仿真結(jié)果，圖6（d）～（f）是本文算法的距離仿真結(jié)果圖.圖6（a）、圖6（b）、圖6（d）中僅有一個(gè)空間譜，而圖6（c）、圖6（e）、圖6（f）中有兩個(gè)空間譜.原因是文獻(xiàn)［19］對(duì)K1和K2近場信源只估計(jì)出一個(gè)DOA，對(duì)K3近場信源估計(jì)出兩個(gè)DOA，而本文算法對(duì)K1近場信源也只估計(jì)出一個(gè)DOA，而對(duì)K2和K3近場信源估計(jì)出兩個(gè)DOA.

圖6 文獻(xiàn)[19]和本文算法的距離估計(jì)結(jié)果圖

4.3 DOA和距離估計(jì)精度仿真

統(tǒng)計(jì)意義的均方根誤差被定義為DOA的估計(jì)精度.其中，近場信源的均方根誤差被定義為：

其中，T是蒙特卡洛實(shí)驗(yàn)的次數(shù)（T=500），K是近場信源的個(gè)數(shù)，βi，j是第j次蒙特卡洛實(shí)驗(yàn)中αi，j的估計(jì)值.本文的算法和文獻(xiàn)［13，16，19］中的方法進(jìn)行了比較.同時(shí)為了更好的驗(yàn)證本文算法的性能，文獻(xiàn)［21］中CRB方法也被比較，結(jié)果如圖7～10所示.

在圖7和圖8中，隨著信噪比的增加，均方根誤差逐漸減少，近場信源的DOA和距離估計(jì)精度被逐漸提高.對(duì)于DOA估計(jì)，文獻(xiàn)［16，19］和本文算法的均方根誤差低于文獻(xiàn)［13］中的ESPRIT-Like算法，但都高于文獻(xiàn)［21］中CRB方法.當(dāng)信噪比低于10 dB時(shí)，本文算法的均方根誤差低于文獻(xiàn)［16，19］和ESPRIT-Like算法.當(dāng)信噪比高于10 dB時(shí)，本文算法和文獻(xiàn)［16，19］的均方根誤差是相似的.因?yàn)樵诘托旁氡葧r(shí)，接收陣列的統(tǒng)計(jì)矩陣的托普利茲結(jié)構(gòu)被破壞.但本文算法首先利用陣列劃分和四階累積量構(gòu)造一個(gè)新的托普利茲矩陣，同時(shí)又利用空間譜函數(shù)的一階導(dǎo)數(shù)構(gòu)造一個(gè)新“空間譜”來減小DOA處的空間相關(guān)性.因此，低信噪比下，本文算法提高了近場信源的DOA的估計(jì)精度.

圖7 DOA的均方根誤差與信噪比關(guān)系圖

對(duì)于距離估計(jì)，當(dāng)信噪比低于10 dB時(shí)，本文算法的均方根誤差低于文獻(xiàn)［13，19］，相似于文獻(xiàn)［16］，當(dāng)信噪比高于10 dB時(shí)，文獻(xiàn)［16］的均方根誤差低于本文算法和文獻(xiàn)［13，19］.因?yàn)閷?duì)于距離估計(jì)，文獻(xiàn)［19］和本文算法的接收陣列的統(tǒng)計(jì)矩陣由協(xié)方差矩陣構(gòu)造的，而文獻(xiàn)［16］的接收陣列的統(tǒng)計(jì)矩陣由四階累積量構(gòu)造的.眾所周知，四階累積量矩陣可以有效的抵抗高斯噪聲，而協(xié)方差矩陣不具有此特性.所以，文獻(xiàn)［16］在高信噪比下的均方根誤差低于本文算法和文獻(xiàn)［13，19］.此外，本文算法同樣采用新“空間譜”估計(jì)近場信源的距離，這可以提高低信噪比下距離的估計(jì)精度.因此，在低信噪比下，本文算法的均方根誤差相似于文獻(xiàn)［16］，低于文獻(xiàn)［13，19］.

在圖9和圖10中，隨著采樣快拍數(shù)的增加，均方根誤差逐漸減少，近場信源的DOA和距離估計(jì)精度被逐漸提高.對(duì)于DOA估計(jì)，本文算法的均方根誤差低于文獻(xiàn)［13，16，19］，但均高于文獻(xiàn)［21］CRB方法.對(duì)于距離估計(jì)，本文算法的均方根誤差低于文獻(xiàn)［13，19］，相似于文獻(xiàn)［16］，但均高于文獻(xiàn)［21］中CRB方法.原因是圖9和10的仿真結(jié)果是在信噪比SNR=0 dB條件下得到的，由圖7可知，低信噪比下本文算法的DOA估計(jì)的均方根誤差低于文獻(xiàn)［13，16，19］，由圖8可知，低信噪比下本文算法的距離估計(jì)的均方根誤差相似于文獻(xiàn)［16］，低于文獻(xiàn)［13，19］.因此，隨著采樣快拍數(shù)的增加，在低信噪比下，本文算法的DOA估計(jì)精度高于文獻(xiàn)［13，16，19］，本文算法的距離估計(jì)精度相似于文獻(xiàn)［16］，高于文獻(xiàn)［13，19］.

圖9 DOA的均方根誤差與采樣快拍數(shù)關(guān)系圖

圖10 距離的均方根誤差與采樣快拍數(shù)關(guān)系圖

4.4 計(jì)算復(fù)雜度分析

在計(jì)算復(fù)雜度分析中，考慮了幾個(gè)主要步驟：統(tǒng)計(jì)矩陣的構(gòu)建，特征值分解和譜峰搜索.為了更好的比較本文算法的計(jì)算復(fù)雜度，對(duì)比了本文算法和文獻(xiàn)［11，13，15，16，19］中算法，具體的計(jì)算復(fù)雜度結(jié)果如表3所示.

在表3中，M、K、L分別表示接收陣列的天線個(gè)數(shù)，信源個(gè)數(shù)和采樣快拍數(shù).ng表示在DOA空間搜索的次數(shù)，ng=180/Δ，其中Δ是搜索步長，Δ=0.01，nl表示在距離空間搜索的次數(shù)其中D是陣列中天線的孔徑，λ為波長.

表3 計(jì)算復(fù)雜度

本文算法的計(jì)算復(fù)雜度高于文獻(xiàn)［13］的ESPRITLike方法，相近于文獻(xiàn)［16］中的方法，低于文獻(xiàn)［11］中的2D-MUSIC算法和文獻(xiàn)［15，19］中的算法.原因是文獻(xiàn)［13］中的ESPRIT-Like算法不需要譜峰搜索，所以計(jì)算量是低于其他算法的，但從圖7和圖8可知，ESPRITLike算法的DOA和距離估計(jì)精度是遠(yuǎn)低于其它算法的；文獻(xiàn)［16］中算法只構(gòu)造了一個(gè)四階累積量，但增加了K次施密特正交化，所以計(jì)算量和本文相似，但從圖7可知，文獻(xiàn)［16］的DOA估計(jì)精度在低信噪比下低于本文算法；文獻(xiàn)［11］需要二維譜峰搜索，文獻(xiàn)［15，19］構(gòu)造的四階累計(jì)量的維度高于本文算法，而本文算法不需要二維譜峰搜索，只需要K+1次一維譜峰搜索（其中，一次譜峰搜索是對(duì)DOA估計(jì)，K次譜峰搜索是對(duì)距離估計(jì)）.

5 結(jié)論

本文基于四階累積量、協(xié)方差矩陣以及空間譜的物理性質(zhì)提出一種基于陣列劃分的近場DOA估計(jì)算法.通過理論分析驗(yàn)證了該算法可以有效的提高DOA估計(jì)分辨率.仿真結(jié)果表明：該算法在低信噪比下提高了DOA的估計(jì)精度，且不需要二維譜峰搜索，只需K+1次一維譜峰搜索.