基于矩陣束的MIMO雷達(dá)低仰角快速估計(jì)方法

謝 榮 劉 崢 劉 俊

(西安電子科技大學(xué)雷達(dá)信號(hào)處理國(guó)家重點(diǎn)實(shí)驗(yàn)室 西安 710071)

1 引言

多輸入多輸出(MIMO)[1]雷達(dá)是目前國(guó)內(nèi)外雷達(dá)技術(shù)領(lǐng)域的一個(gè)新的研究熱點(diǎn)，與傳統(tǒng)雷達(dá)相比具有明顯的優(yōu)勢(shì)，已出現(xiàn)了很多針對(duì)MIMO雷達(dá)的參數(shù)估計(jì)方法[2,3]。

由于低空目標(biāo)的直達(dá)波與地面的多徑信號(hào)在同一波束寬度范圍內(nèi)疊加，產(chǎn)生多徑衰落現(xiàn)象，出現(xiàn)了低仰角目標(biāo)檢測(cè)和參數(shù)估計(jì)困難的問題。近年來很多學(xué)者研究了常規(guī)相控陣?yán)走_(dá)低角跟蹤時(shí)的角度估計(jì)問題[4?8]，包括基于數(shù)據(jù)協(xié)方差矩陣估計(jì)算法[4?7]和基于最大似然估計(jì)算法[8]。但是對(duì)多徑環(huán)境下的MIMO雷達(dá)參數(shù)估計(jì)問題，研究的學(xué)者很少[9]。由于MIMO雷達(dá)是多輸入多輸出雷達(dá)，其多徑信號(hào)更為復(fù)雜。基于數(shù)據(jù)協(xié)方差矩陣的算法不能用于單樣本數(shù)的情況，對(duì)信噪比要求也比較高，且基于傳統(tǒng)空間平滑技術(shù)的空間譜估計(jì)方法不能直接應(yīng)用。由于最大似然算法具有低信噪比門限和所需快拍數(shù)少的優(yōu)點(diǎn)，在相關(guān)信號(hào)的角度估計(jì)中能得到較好的結(jié)果，文獻(xiàn)[9]提出了一種基于最大似然的低仰角估計(jì)算法，算法考慮了收、發(fā)多徑信號(hào)，并且利用了直達(dá)波角與多徑反射角的關(guān)系將2維搜索簡(jiǎn)化為1維角度搜索，但當(dāng)陣元較多和搜索角度間隔較小時(shí)，其計(jì)算量仍然很可觀，特別是同時(shí)估計(jì)多目標(biāo)時(shí)，最大似然算法通常需要多維譜峰搜索，計(jì)算量將急劇增大。實(shí)際上為了便于工程應(yīng)用，提高測(cè)角精度與減小算法的運(yùn)算量同樣重要。Sarkar等人[10,11]提出了矩陣束(Matrix Pencil, MP)方法，該方法不僅可減少角度估計(jì)時(shí)的運(yùn)算量，而且適用于常規(guī)相控陣?yán)走_(dá)的相干目標(biāo)情況。

本文針對(duì)MIMO雷達(dá)低仰角估計(jì)問題，提出了一種基于矩陣束的MIMO雷達(dá)低仰角快速估計(jì)方法。該方法同時(shí)考慮了收、發(fā)多徑信號(hào)，根據(jù)接收數(shù)據(jù)矢量的特點(diǎn)，構(gòu)造一個(gè)前后向矩陣束，在單個(gè)樣本數(shù)下，直接對(duì)MIMO雷達(dá)多徑下的目標(biāo)角度進(jìn)行估計(jì)，避免了譜峰搜索。仿真結(jié)果表明，在低信噪比、單樣本數(shù)情況下，本文算法在保證低仰角估計(jì)性能的基礎(chǔ)上，大大減少了計(jì)算量，且在基本不增加計(jì)算量的情況下，可以同時(shí)估計(jì)多目標(biāo)角度值，利于工程實(shí)時(shí)處理。

2 信號(hào)模型

考慮收發(fā)共置的單基地MIMO雷達(dá)。如圖1所示，假設(shè)MIMO雷達(dá)架設(shè)在菲涅耳區(qū)較為平坦的表面上。由于是考慮俯仰上的測(cè)角問題，因此假設(shè)天線收發(fā)陣列均是垂直放置等距線陣，發(fā)射陣元數(shù)目為M，接收陣元數(shù)目為N，陣元間隔均為d,ha為天線高度，ht為遠(yuǎn)場(chǎng)目標(biāo)T的高度，Rd為天線和目標(biāo)間的直線距離，Rs為目標(biāo)經(jīng)過地面反射后到達(dá)雷達(dá)的多徑距離，θd和θs分別為進(jìn)入天線的直達(dá)波和多徑信號(hào)的入射角。假設(shè)水平方向以上的方向?yàn)榻嵌鹊恼较颉D1中，實(shí)線代表發(fā)射多徑信號(hào)，虛線表示接收多徑信號(hào)。

圖1 低仰角目標(biāo)多徑幾何模型

由于 MIMO雷達(dá)各發(fā)射信號(hào)se(t)相互正交，即滿足

式中sem(t),sel(t)分別表示第m個(gè)和第l個(gè)發(fā)射陣元的發(fā)射信號(hào)，Te表示發(fā)射信號(hào)的脈沖寬度，上標(biāo)H為矩陣共軛轉(zhuǎn)置運(yùn)算。假設(shè)發(fā)射信號(hào)采用變頻窄帶信號(hào)，且各發(fā)射信號(hào)功率相同，滿足式(1)的第m個(gè)天線發(fā)射的信號(hào)為

式中g(shù)(t)為門控信號(hào)用于對(duì)變頻信號(hào)進(jìn)行截?cái)啵籄為發(fā)射信號(hào)幅度；fm=f0+cmΔf為第m個(gè)天線的發(fā)射頻率，cm為發(fā)射信號(hào)頻率編碼，Δf為發(fā)射信號(hào)頻率間隔，滿足 Δf=1/Te，且f0?MΔf; φm為初始相位，不失一般性可設(shè)其為零。

考慮單目標(biāo)的多徑情況，發(fā)射信號(hào)sem(t)到達(dá)目標(biāo)和經(jīng)過目標(biāo)后反射到接收陣列的信號(hào)均包含直達(dá)和反射兩條路徑，此時(shí)第n個(gè)接收陣元接收的第m個(gè)發(fā)射陣元的回波信號(hào)為

式中F=exp (?j2 πfmτ), τ為陣列中心到目標(biāo)的距離 產(chǎn)生 的時(shí) 延 ，an(θ)=exp(j2 πd(n?1) sinθ/ λ),(θ=θd, θs) , τdm=(m?1)dsin θd/c為第m個(gè)發(fā)射陣元到陣列中心的直達(dá)波波程時(shí)延差，τsm=(m?1)dsin θs/c為第m個(gè)發(fā)射陣元到陣列中心的反射 波 波 程 時(shí) 延 差 ， τsn=τdn+Δτ,Δτ=Rd/c·(cos θd/cos θs? 1)為直達(dá)波和反射波的波程差所對(duì)應(yīng)的時(shí)延，c為光速，λ為發(fā)射信號(hào)載波波長(zhǎng)，ε=ρ ·ex p (?j2πf0Δτ), ρ 為地面反射系數(shù)，nmn(t)為此時(shí)回波信號(hào)的背景噪聲，假設(shè)它為零均值高斯白噪聲且與信號(hào)不相關(guān)。

假設(shè)在遠(yuǎn)場(chǎng)平面波情況下，忽略包絡(luò)時(shí)延，并假設(shè)在觀測(cè)時(shí)間內(nèi)反射系數(shù)不變。利用各發(fā)射信號(hào)之間的正交性，將發(fā)射信號(hào)se(t)與第n個(gè)接收陣元接收的第m個(gè)發(fā)射陣元的回波信號(hào)rmn(t)進(jìn)行相乘，可以分離出一個(gè)信號(hào)矢量Γnm：

式中Q為一個(gè)脈沖寬度內(nèi)的采樣點(diǎn)數(shù)，nnm=為噪聲分量。對(duì)N個(gè)接收陣元接收的所有發(fā)射陣元的回波信號(hào)進(jìn)行上述處理，可獲得一個(gè)MN×1維的綜合信號(hào)矢量[9]：

3 算法實(shí)現(xiàn)

根據(jù)式(5)可知，此時(shí)每個(gè)接收陣元只獲得一個(gè)樣本數(shù)，大量的基于協(xié)方差矩陣的方法不能適用，文獻(xiàn)[9]采用最大似然算法，運(yùn)算量很大，且估計(jì)精度與搜索步長(zhǎng)有關(guān)。本文通過構(gòu)造矩陣束進(jìn)行低仰角估計(jì)。

式中

其中Yl,l=1,2,… ,L+ 1為一個(gè)大小為M(N?L)×1的信號(hào)子向量，即加窗后所得的信號(hào)矢量，L為自由參數(shù)，取值范圍為P≤L≤N?P，其中P為目標(biāo)個(gè)數(shù)(包括所有直達(dá)波和反射波的角度)。由于這里分析多徑情況下單目標(biāo)的情況(即P=2)，其取值范圍為2 ≤L≤N? 2。在這個(gè)取值范圍內(nèi)，不同的L值對(duì)應(yīng)著不同的角度估計(jì)性能，根據(jù)文獻(xiàn)[10]的分析，L取為N/ 3 ～N/ 2時(shí)估計(jì)性能最優(yōu)。

根據(jù)前向Hankel矩陣Zf，定義如下兩個(gè)Hankel子矩陣：

將式(5)的信號(hào)形式代入，可以將子矩陣和進(jìn)一步展開，寫成如下形式：

式中

通過分析，根據(jù)式(9a)和式(9b)可以得到如下的矩陣束：

式中I2為2×2的單位陣。由于2 ≤L≤N? 2，一般情況下矩陣束的秩為2，但當(dāng)γ=exp(j2πdsin θd/λ)或 γ=exp(j2 πdsin θs/λ)時(shí) ，中將會(huì)有一行為零，此時(shí)矩陣束的秩將為1，可利用基于廣義特征值分解的方法來估計(jì)目標(biāo)的角度[10,11]。先由式(12)可以得到如下矩陣：

再對(duì)矩陣Σ進(jìn)行特征值分解，便可以得到低仰角目標(biāo)的角度值：

其中γ1,2為矩陣Σ經(jīng)特征值分解后所得最接近單位圓的兩個(gè)特征值，其中正角度值即是目標(biāo)的低仰角θd。以上分析中均沒有考慮噪聲的影響，為了進(jìn)一步減小單樣本數(shù)下觀測(cè)噪聲對(duì)角度估計(jì)精度的影響，利用綜合信號(hào)矢量Γ，重新構(gòu)造一個(gè)后向Hankel矩陣：

式中上標(biāo)*為矩陣共軛運(yùn)算。然后結(jié)合式(6)構(gòu)造一個(gè)前后向矩陣為

同理，根據(jù)前后向矩陣Zfb，定義兩個(gè)前后向子矩陣為

在噪聲環(huán)境中，文獻(xiàn)[11]利用了Hankel矩陣的性質(zhì)來計(jì)算式(12)的廣義特征值，而式(18)不滿足文獻(xiàn)[11]中的計(jì)算條件，這里利用廣義特征值分解的總體最小二乘法來計(jì)算。對(duì)進(jìn)行奇異值分解，表示為

式中E2s由非零的大奇異值組成。在不改變廣義特征值的條件下，可以通過對(duì)前后向矩陣束進(jìn)行左乘和右乘V2s運(yùn)算進(jìn)行降維，得到新的矩陣束：

降維處理后，原維數(shù)為 2M(N?L)×L的矩陣束的廣義特征值問題轉(zhuǎn)化為維數(shù)為2×2的矩陣束的廣義特征值問題，且變換中只需要一次奇異值分解。

綜上所述，本文提出的基于矩陣束的MIMO雷達(dá)低仰角估計(jì)方法的主要步驟如下：

(1)用se(t)與整個(gè)回波信號(hào)相乘進(jìn)行分離，獲得綜合信號(hào)矢量Γ；

(5)根據(jù)式(14)求解低仰角目標(biāo)的角度值。

4 計(jì)算量分析

本文算法和文獻(xiàn)[9]的最大似然算法都是對(duì)一個(gè)MN× 1維的綜合信號(hào)矢量Γ進(jìn)行處理。根據(jù)前面的分析可知，本文算法的計(jì)算量主要集中在2M(N?L) ×L維矩陣的奇異值分解和矩陣束{E2s的構(gòu)造上，不需要額外的搜索過程，其 計(jì) 算 量(次 復(fù) 乘)約 為O(L3+ 2M(N?L)L2+4M(N?L)L+ 4L)，當(dāng)收發(fā)陣元共置M=N，自由參 數(shù)L=N/2時(shí) ， 計(jì) 算 量 約 為O(N4/4 + 17N3/8+ 2N)。而文獻(xiàn)[9]算法的計(jì)算量主要集中在投影矩陣PA和最大似然函數(shù)TN(θ)的計(jì)算上，由于仍需要1維搜索，搜索步長(zhǎng)Δθ與精度需求有關(guān)，且在每個(gè)搜索角度下就需要計(jì)算一次PA和TN(θ)，如果搜索角域被劃分為Kd等分，則計(jì)算量約為Kd·O(5 (MN)2+ 33MN)，當(dāng)收發(fā)陣元共置M=N時(shí)，計(jì)算量為Kd·O(5N4+ 33N2)，而一般Kd＞＞N，所以本文方法的計(jì)算量要明顯低于文獻(xiàn)[9]方法。

5 算法仿真和性能分析

下面對(duì)本文算法的有效性和估計(jì)性能進(jìn)行仿真和分析。仿真中采用圖1所示的低仰角目標(biāo)多徑幾何模型，假設(shè)天線中心頻率f0=150 MHz，目標(biāo)距離Rd=120 km ，天線中心高度ha=110 m ，本文算法無需已知反射系數(shù)，不失一般性假設(shè)反射系數(shù)ρ=0.9exp(jπ)。

這里定義角度估計(jì)值的均方根誤差為

圖2給出了信噪比SNR=10 dB時(shí)在不同目標(biāo)高度ht下本文算法和文獻(xiàn)[9]方法對(duì)多徑目標(biāo)直達(dá)波入射角度θd估計(jì)的均方根誤差的變化曲線。仿真中，收發(fā)共置陣元個(gè)數(shù)M=N=16，陣元間隔d=λ /2，采樣點(diǎn)數(shù)Q=100，取自由參數(shù)L=N/2=8，其中文獻(xiàn)[9]算法的搜索步長(zhǎng)為 Δθ=0.01o。仿真采用100次Monte-Carlo試驗(yàn)，其中目標(biāo)高度ht從1000 m變化到8000 m。

從圖2中可以看出本文算法可以有效地估計(jì)多徑下的目標(biāo)角度，且隨著目標(biāo)高度的增大，本文算法的估計(jì)精度越來越接近文獻(xiàn)[9]最大似然算法的精度。當(dāng)目標(biāo)高度ht為2000 m時(shí)，直達(dá)波入射角約為θd≈ 0.9o，此時(shí)本文算法的角度估計(jì)誤差為 0.046o，對(duì)應(yīng)的高度誤差約為207 m。此時(shí)與最大似然算法相比，角度估計(jì)誤差僅相差 0.027o，對(duì)應(yīng)的高度誤差僅相差約57 m，但本文方法在計(jì)算量上大大減少，實(shí)時(shí)性更好。

圖2 θd估計(jì)均方根誤差與目標(biāo)高度的變化關(guān)系

圖3為本文算法和文獻(xiàn)[9]算法對(duì)多徑下目標(biāo)直達(dá)波入射角度θd估計(jì)的均方根誤差(RMSE)隨信噪比SNR變化的關(guān)系。仿真中，M=N=16,d=λ/2,Q=100,L=8,ht=5000 m，直達(dá)波入射角 θd≈2.33o，反射波入射角度 θ≈? 2.44o，其中文獻(xiàn)[9]s算法的搜索步長(zhǎng)為 0.01o。仿真采用100次Monte-Carlo試驗(yàn)。從圖中可以看出本文算法可以有效地估計(jì)多徑下的目標(biāo)角度，且隨著信噪比的增大，本文算法的估計(jì)精度越來越接近文獻(xiàn)[9]算法的精度。需要注意的是，在16 dB以上時(shí)，圖中本文算法的精度超過了最大似然算法，這是因?yàn)榉抡嬷凶畲笏迫凰惴ú捎玫慕嵌人阉鞑介L(zhǎng)為 0.01o，而本文算法是直接求解角度，因此計(jì)算出來的精度更高。

圖3 θd估計(jì)均方根誤差與SNR的變化關(guān)系

最后分析多徑環(huán)境下同時(shí)存在多個(gè)目標(biāo)的情況。根據(jù)前面的計(jì)算量分析可知，此時(shí)文獻(xiàn)[9]最大似然算法中的聯(lián)合導(dǎo)向矢量A的未知量將增多，需要多維的搜索來進(jìn)行角度估計(jì)，計(jì)算量大大增加；而本文方法由于采用了降維處理，且一般N＞P，可以認(rèn)為計(jì)算量基本不變。下面通過仿真分析本文方法同時(shí)估計(jì)兩個(gè)目標(biāo)角度時(shí)的情況。

圖4為采用本文算法同時(shí)估計(jì)兩個(gè)低仰角目標(biāo)角度的情況，橫坐標(biāo)為Monte-Carlo試驗(yàn)次數(shù)，縱坐標(biāo)為采用本文算法估計(jì)的目標(biāo)角度。仿真中，M=N=24,d=λ /2,Q=100,L=12,Monte-Carlo試驗(yàn)次數(shù)為200。假設(shè)兩個(gè)目標(biāo)的信噪比均為15 dB，目標(biāo)1的高度為4000 m，直達(dá)波入射角θd≈ 1.86o；目標(biāo)2的高度為7000 m，直達(dá)波入射角θd≈ 3.29o。

圖4 同時(shí)兩個(gè)目標(biāo)時(shí)的估計(jì)結(jié)果

從圖4中可以看出，本文算法可以同時(shí)有效地實(shí)現(xiàn)兩個(gè)目標(biāo)的角度估計(jì)，此時(shí)目標(biāo)1和目標(biāo)2的角度估計(jì)均方根誤差分別為 0.1029o和 0.0584o，其中由于目標(biāo)2的高度稍高，估計(jì)的角度值誤差較小。根據(jù)估計(jì)的仰角值換算獲得目標(biāo)對(duì)應(yīng)的估計(jì)高度，其中目標(biāo)1的估計(jì)高度與真值的最大偏移約為637.6 m，目標(biāo)2的估計(jì)高度與真值的最大偏移約為372.1 m，能滿足一般測(cè)高誤差為距離的2%的設(shè)計(jì)要求。

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

本文針對(duì)MIMO雷達(dá)低仰角估計(jì)問題，提出了一種基于矩陣束的快速低仰角估計(jì)方法。該方法同時(shí)考慮了收、發(fā)多徑信號(hào)，并對(duì)前后向矩陣束進(jìn)行降維，最后直接利用廣義特征值分解的總體最小二乘法對(duì)目標(biāo)角度進(jìn)行估計(jì)。理論分析和仿真結(jié)果表明，本文算法可以有效地克服多徑效應(yīng)的影響，不需要估計(jì)協(xié)方差矩陣，在低信噪比、單樣本數(shù)時(shí)，算法在保證低仰角估計(jì)性能和不增加計(jì)算量的基礎(chǔ)上，可以同時(shí)估計(jì)多目標(biāo)角度值，利于工程實(shí)時(shí)處理。

[1] Fishler E, Alexander M, Blum R,et al.. MIMO radar: an idea whose time has come[C]. Proc. of the IEEE Radar Conference,Philadelphia, PA, United States, 2004: 71-78.

[2] 陳金立, 顧紅, 蘇衛(wèi)民. 一種雙基地MIMO雷達(dá)快速多目標(biāo)定位方法[J]. 電子與信息學(xué)報(bào), 2009, 31(7): 1664-1668.Chen Jin-li, Gu Hong, and Su Wei-min. A method for fast multi-target localization in bistatic MIMO radar system[J].Journal of Electronics&Information Technology, 2009, 31(7):1664-1668.

[3] 謝榮, 劉崢. 基于多項(xiàng)式求根的雙基地 MIMO 雷達(dá)多目標(biāo)定位方法[J]. 電子與信息學(xué)報(bào), 2010, 32(9): 2197-2200.Xie Rong and Liu Zheng. Multi-target localization based on polynomial rooting for bistatic MIMO radar[J].Journal of Electronics&Information Technology, 2010, 32(9):2197-2200.

[4] 胡曉琴, 陳建文, 王永良, 等. 一種適用于米波雷達(dá)低角測(cè)高環(huán)境的快速算法[J]. 信號(hào)處理, 2009, 25(5): 761-765.Hu Xiao-qin, Chen Jian-wen, Wang Yong-liang,et al.. A fast algorithm for the meter-wave radar low-angle height-finding[J].Signal Processing, 2009, 25(5): 761-765.

[5] 董玫, 趙永波, 張守宏, 等. 傾斜面陣下的米波雷達(dá)測(cè)高方法[J]. 西安電子科技大學(xué)學(xué)報(bào), 2009, 36(1): 43-47.Dong Mei, Zhao Yong-bo, Zhang Shou-hong,et al.. Method for altitude measurement in the VHF radar under the oblique planar array[J].Journal of Xidian University, 2009, 36(1):43-47.

[6] 趙光輝, 陳伯孝, 吳向東, 等. 基于差分預(yù)處理的米波雷達(dá)低仰角處理算法[J]. 電子與信息學(xué)報(bào), 2009, 31(2): 363-365.Zhao Guang-hui, Chen Bai-xiao, Wu Xiang-dong,et al.. An algorithm based on differential preprocessing of low elevation estimation in VHF radar[J].Journal of Electronics&Information Technology, 2009, 31(2): 363-365.

[7] 吳向東, 趙永波, 張守宏等. 利用改進(jìn)的Toeplitz化技術(shù)實(shí)現(xiàn)米波雷達(dá)低仰角測(cè)高[J]. 電子與信息學(xué)報(bào), 2010, 32(3):570-574.Wu Xiang-dong, Zhao Yong-bo, Zhang Shou-hong,et al..Height finding of meter-wave radar using improved Toeplitz technique at low-angle environment[J].Journal of Electronics&Information Technology, 2010, 32(3): 570-574.

[8] 趙永波, 張守宏. 雷達(dá)低角跟蹤環(huán)境下的最大似然波達(dá)方向估計(jì)方法[J]. 電子學(xué)報(bào), 2004, 32(9): 1520-1523.Zhao Yong-bo and Zhang Shou-hong. Maximum likelihood DOA estimation in radar low-angle tracking environment[J].Acta Electronica Sinica, 2004, 32(9): 1520-1523.

[9] 吳向東, 趙永波, 張守宏, 等. 一種MIMO雷達(dá)低角跟蹤環(huán)境下的波達(dá)方向估計(jì)新方法[J]. 西安電子科技大學(xué)學(xué)報(bào), 2008,35(5): 793-798.Wu Xiang-dong, Zhao Yong-bo, Zhang Shou-hong,et al.. New method for DOA estimation for the MIMO radar in low-angle tracking environment[J].Journal of Xidian University, 2008,35(5): 793-798.

[10] Sarkar T K and Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[J].IEEE Antennas Propagation Magazine, 1995, 37(1): 48-54.

[11] Yilmazer N, Koh J, and Sarkar T K. Utilization of a unitary transform for efficient computation in the matrix pencil method to find the direction of arrival[J].IEEE Transactions on Antennas and Propagation, 2006, 54(1): 175-181.