基于L 型延遲陣列調(diào)制寬帶轉(zhuǎn)換器的信號載頻和二維到達角聯(lián)合估計*

姜思儀 付寧 喬立巖 彭喜元

(哈爾濱工業(yè)大學電子與信息工程學院, 哈爾濱 150001)

1 引 言

陣列信號處理已被廣泛應用于雷達[1]、聲納[2,3]和無線通信[4,5]等領域, 通過按一定樣式在空間中分布的傳感器陣列接收空間電磁信號, 并進行一系列處理得到空時頻各域的參數(shù).在實際應用中, 頻率 和 二 維 到 達 角(two-dimensional direction of arrival, 2D-DOA)是對電磁波進行識別的重要特征, 對載頻和2D-DOA 的聯(lián)合估計在陣列信號處理領域備受關注.

隨著信號頻率范圍越來越寬, 基于傳統(tǒng)奈奎斯特定理的信號采集方法給采樣、存儲和處理設備帶來了巨大的壓力.近來新興的壓縮感知(compressed sensing, CS)理論[6-9]能夠在少量樣本的情況下準確地恢復稀疏信號, 該理論大大降低了采樣率, 具有很高的應用價值.目前已有不少學者將CS 理論應用到模擬信號的采樣中[10-15].其中, 由Mishali 和Eldar[14,15]提出的MWC 系統(tǒng)是一種專門針對多頻帶模擬信號采樣的模擬信息轉(zhuǎn)換系統(tǒng).

為了降低陣列中多個通道同時采樣產(chǎn)生的龐大數(shù)據(jù)量, 近來許多學者提出了基于欠采樣的載頻和到達角聯(lián)合估計方法, 黃翔東等[16]利用互素稀疏陣列, 結合閉式中國余數(shù)定理實現(xiàn)了欠采樣下頻率和DOA 的聯(lián)合估計.沈志博等[17]提出一種CS奇異值分解算法, 但配對步驟較為繁瑣.Liu 等[18,19]提出了一種基于多陪集結構的聯(lián)合估計方法, 分別基于子空間分解和正則分解/平行因子分析(CANDECOMP/PARAFAC, CP)[20]方法, 但接收陣列結構較為復雜.接著, 他們又進行了擴展[21,22], 利用稀疏陣列減少了所需的陣元個數(shù), 但這種結構對ADC 仍有比較高的要求.Stein 等[23,24]將MWC技術與陣列信號處理相結合, 提出了L 型陣列MWC 結構來接收信號, 并提出了一種基于旋轉(zhuǎn)不變子空間(estimating signal parameter via rotational invariance techniques, ESPRIT)方法的聯(lián)合載波DOA 恢復算法.Cui 等[25]提出了一種基于均勻線性MWC 陣列由兩階段估計和參數(shù)配對構成的估計方法, 算法復雜且精度有限.Chen 等[26]提出了一種基于多個MWC 結構的DOA 和載頻聯(lián)合估計方法, 魯棒性較好, 但該結構所需的通道數(shù)較多, 硬件成本高.然而在實際中電磁波存在于三維空間中, 使用2D-DOA 來描述目標的方向更為準確, 但目前對基于欠采樣的載頻和2D-DOA聯(lián)合估計研究較少.陳玉龍等[27]提出了一種基于L 型陣列的2D-DOA 估計的欠采樣算法, 但所需的天線陣元個數(shù)較多.Esmaeil 等[28]提出了一種基于雙L 型陣列結構的網(wǎng)格CS 估計方法, 但是參數(shù)估計精度受網(wǎng)格的大小限制.他們還提出了一種貝葉斯壓縮感知算法[28], 方法不依賴網(wǎng)格但步驟復雜且計算量較大.

從以上分析可以發(fā)現(xiàn), 現(xiàn)有的方法大多都只考慮一維DOA 的估計, 需要假設目標同時位于一個平面上, 然而在實際應用中目標一般分散在三維空間中, 因此需要2D-DOA 估計.而基于欠采樣的載頻與2D-DOA 聯(lián)合估計方法研究較少, 現(xiàn)有的方法中存在結構冗余, 精度受網(wǎng)格限制和算法復雜等問題.因此在壓縮采樣條件下對信號的2D-DOA和頻率參數(shù)聯(lián)合估計還有待研究.

L 型陣列結構有設計簡單, 陣列方向圖主瓣窄, 估計性能較好等優(yōu)點, 因此本文提出了一種基于MWC 技術的L 型延遲陣列結構, 在L 型陣列結構的基礎上增加延遲通道來增加可估計到達角的維度, 利用延遲通道與未延遲通道的采樣值之間的已知的相位差可直接估計載頻, 進而計算2DDOA 參數(shù), 無需額外的載頻和2D-DOA 的參數(shù)配對操作, 避免了配對所引入的誤差和復雜度的提升.并結合L 型延遲陣列結構的特點構造相關矩陣和三線性模型, 提出了兩種聯(lián)合載波頻率和2DDOA 估計方法.其中一種基于ESPRIT 算法, 另一種基于CP 分解技術.基于ESPRIT 的方法有計算量小, 無需譜峰搜索的優(yōu)點, 適用于實時性要求較高的應用場合; 基于CP 分解技術的方法利用三維數(shù)據(jù)的結構特性, 有魯棒性較好的優(yōu)點, 適用于信噪比較低的應用場合, 本文主要基于這兩種方法進行研究.本文第2 節(jié)介紹了信號模型; 第3 節(jié)首先介紹了L 型延遲陣列接收結構, 接著給出了接收信號模型; 第4 節(jié)首先對所提的兩種方法的原理和步驟進行了詳細的闡述.接著給出了參數(shù)估計條件, 最后對復雜度進行了詳細的分析.第5 節(jié)給出了本文所提算法的驗證結果以驗證該算法的正確性和有效性.

2 問題描述

考慮M個遠場目標發(fā)出的窄帶信號si(t)ej2πfi,其中i=1,2,··· ,M.信號的頻率范圍為F=[-fNyq/2,fNyq/2] , 其中fNyq為奈奎斯特頻率.si(t)為互不相關的基帶信號, 帶寬不超過B且遠小于調(diào)制載頻fi.假設信源發(fā)出的信號以方位角θi ∈(-90°,90°) 和俯仰角φi ∈(0°,90°) 的入射方向被天線陣列接收, 為了避免空間模糊, 假設載頻、方位角與俯仰角滿足[24]:

本文的目標是設計一種基于MWC 技術的陣列結構, 使其能夠從欠采樣樣本中完全估計上述信號的載頻, 俯仰角和方位角參

數(shù).

3 L 型延遲陣列與接收信號模型

3.1 L 型延遲陣列結構

提出圖1 所示的L 型延遲天線陣列結構, 利用延遲通道的采樣值可以先估計載頻, 進而計算2D-DOA 參數(shù), 避免了額外的參數(shù)配對操作.該結構由兩個分別沿x軸和y軸正方向均勻分布的N個線性天線陣列{x1,x2,··· ,xN}和{y1,y2,··· ,yN}組成, 兩個軸在原點共用同一個天線陣元.相鄰陣元之間的間距d滿足d≤c/fNyq, 其中c為光速.

圖1 L 型陣列MWC 結構圖Fig.1.L shaped array MWC.

每個天線陣元后依次連接混頻器、低通濾波器和采樣模塊, 類似MWC 通道結構.在x軸的每個陣元后增加一個延遲通道, 在混頻之前對接收信號進行固定延時τ, 如圖2 所示.

圖2 x 軸延遲通道結構圖Fig.2.x-axis delay channel.

3.2 接收信號模型

定義xn(t)和yn(t) 分別為x軸和y軸的第n個天線陣元接收的信號,(t) 為xn(t) 經(jīng)過延遲τ后的信號, 由于目標信號為窄帶信號, 有si(t+τ)≈si(t)[24], 則接收信號可以表示為

圖3 接收信號頻譜示意圖Fig.3.Spectrum of received signal.

陣元接收信號在混頻器與周期為Tp=1/fp的± 1 偽隨機序列p(t) 相乘后, 再經(jīng)過截止頻率為fp/2 的低通濾波器, 只留下Fp=[-fp/2,fp/2] 內(nèi)的頻譜.濾波后的信號經(jīng)過傅里葉變換可得:

濾波后的信號以fs=fp的速率進行低速采樣得到觀測值, 將其傅里葉變換后寫為矩陣形式可得:

其 中X(f) 為N×1 的矩陣, 第n個元素為Xn(ej2πfTs).W(f) 是M×1 的 矩 陣, 第i個 元 素為矩陣為陣列流型矩陣, 其中

定義x軸通道的觀測值為x[k]=[x1[k],··· ,xN[k]]T, 對(4)式取離散時間傅里葉逆變換, 在時域有:

類似地, 對于y軸和x軸延遲通道的采樣值y[k]=[y1[k], y2[k]··· , yN[k]]T和有:

本文的目的是從欠奈奎斯特采樣值x[k],y[k]和x′[k] 中 估計信號的載頻俯仰角和方位角參數(shù).

4 載頻與二維DOA 聯(lián)合估計方法

本節(jié)給出了基于圖1 接收結構所提的兩種載頻與二維DOA 聯(lián)合估計方法, 一種方法基于ESPRIT 技術, 另一種方法基于CP 分解.對兩種方法的原理和步驟進行了詳細的闡述.接著給出了參數(shù)估計條件, 最后對復雜度進行了詳細的分析.

4.1 基于ESPRIT 的估計方法

首先給出基于ESPRIT 的參數(shù)聯(lián)合估計方法,將采樣值分為前N—1 項和后N—1 項組成的兩個子矩陣, 并構造如下的相關矩陣:

其中Ax1和Ay1分別為Ax和Ay的前N—1 行,

構造如(9)式的矩陣R,

對R進行奇異值分解, 得到前M個非零奇異值對應的左奇異向量U=[U1;U2;U3;U4;U5;U6].U是一個 6 (N -1)×M的矩陣,Ui為 (N -1)×M的矩陣,i=1,2,··· ,6.通過推導可得矩陣Ui和對角陣Φx,Φy和Ψ之間的關系, 定義矩陣V1,V2和V3,

對矩陣 (V1+V2+V3) 進行特征值分解, 這樣可以保證矩陣Φx,Φy和Ψ中元素順序?qū)? 可得到特征向量矩陣, 進而通過(13)式—(15)式計算可以得到

目標信號的載頻fi, 方位角θi和俯仰角φi就可以計算,

其中 ∠ (·) 為求復數(shù)的相位角, 綜上所述, 現(xiàn)將基于ESPRIT 的估計方法的步驟總結為如下.

輸入各通道Q快拍采樣值x[k] ,y[k]和x′[k] ,其中k=1,··· ,Q;

輸出載頻俯仰角和方位角

步驟1根據(jù)(8)式和(9)式計算協(xié)方差矩陣R;

步驟2對R進行奇異值分解, 得到左奇異向量U, 根據(jù)(10)式—(12)式計算矩陣V1,V2和V3;

步驟3對矩陣 (V1+V2+V3) 進行特征值分解, 根據(jù)(13)式—(15)式計算得到

步驟4根據(jù)(16)式—(18)式計算載頻俯仰角和方位角

4.2 基于CP 分解的估計方法

CP 分解[17]可以將一個高維的張量分解成多個秩為1 的張量的和, 從而降低待處理參數(shù)的維度.為了使(9)式中的矩陣R滿足CP 分解的模型,定義一個大小為 6×M的矩陣R=[R1R2···RM] , 其第i列向量Ri為

由于信號滿足互不相關, 因此自相關矩陣Rw中除對角線元素外均為0.定義一個三階張量χ(N-1)×(N-1)×6, 其正向切片Xk=Rk, k=1,2,...,6,Rk可以表示為

通過觀察可以發(fā)現(xiàn)(20)式符合典型的CP 分解模型, 陣列流型矩陣Ax1和Ay1均為范德蒙德矩陣, 當N≥M+1 時,Ax1和Ay1均 為滿k-秩 等于M.由于信源滿足(1)式, 由(19)式的定義可以看出,矩陣R一定有任意兩列獨立, 即kR≥2.因此有

滿足CP 分解的唯一性條件[17], 在不考慮列置換模糊和尺度模糊的條件下, 因子矩陣Ax1,Ay1和R可以從張量χ中唯一恢復.通過固定除一個因子矩陣外的所有矩陣的方式, 將CP 轉(zhuǎn)化為以下的線性最小二乘的問題:

其中X(n)為張量X的n-模式矩陣化, 因子矩陣和分別為Ax1,Ay1和R的估計值.利用交替最小二乘法進行循環(huán)求解, 可解得三個因子矩陣.對中的元素計算可得(23)式—(25)式, 為了方便將其記為

目標信號的載頻fi, 方位角θi和俯仰角φi參數(shù)就可以計算:

綜上所述, 現(xiàn)將基于CP 分解的估計方法的步驟總結為如下.

輸入各通道Q快拍采樣值x[k],y[k]和x′[k] ,其中k=1,··· ,Q;

輸出載頻俯仰角和方位角

步驟1根據(jù)(8)式和(20)式計算三線性模型χ;

步驟2利用交替最小二乘法對三線性模型χ進行求解, 可求得因子矩陣

步驟3利用中元素根據(jù)(23)式—(25)式計算得和

步驟4根據(jù)(26)式—(28)式計算載頻俯仰角和方位角

4.3 參數(shù)估計條件

由于L 型延遲陣列MWC 結構是基于MWC理論[14,15]展開的, 因此估計條件有一定的一致性,結合文獻[14]中MWC 理論的重構條件給出如下的參數(shù)成功估計的條件.

對于滿足(1)式的M個互不相關的窄帶遠場信號使用圖1 所示的共有2N—1 個陣元的L 型延遲陣列MWC 結構接收, 當滿足以下條件:

1)陣元間距d≤c/fNyq, 陣元個數(shù)N≥M+1 ,通道延遲τ≤1/fNyq;

2)fs≥fp≥B;

4.4 復雜度分析

接下來分析基于ESPRIT 的參數(shù)估計算法的復雜度.計算矩陣R大約需要N2×Q次復乘, 其中Q為每通道快拍數(shù),N為一個軸上的陣元個數(shù).對R進行奇異值分解大約需要 6 (N -1)3次復乘,對矩陣 (V1+V2+V3) 進 行特征值分解, 需要M3次復乘,M為信源個數(shù).因此, 該算法的時間復雜度大約為O(N3+N2Q+M3) , 所需的空間復雜度約O(N2+M2).

同樣, 對基于CP 分解的參數(shù)估計算法的復雜度進行分析.交替最小二乘算法的每次循環(huán)中都要對3 個矩陣分別進行求解, 其中計算R⊙Ay1大約需 要 (N -1)2×M次 復 乘, 計 算R⊙Ay1的 偽 逆則需要2(N-1)×M2+M3次復乘,計算X(1)與的乘積則需要(N-1)3×M次復乘.假設交替最小二乘算法循環(huán)I次, 則CP 分解的時間復雜度大約為O(IM ·N3+INM2+N2Q) , 其大小和迭代循環(huán)次數(shù)有直接關系, 迭代的次數(shù)由設定的閾值和最大循環(huán)次數(shù)有關.在計算過程中需要存儲矩陣X(k)和因子矩陣等, 所需的空間復雜度約為O(N3).可以看出, CP 分解方法的單次循環(huán)的計算復雜度與ESPRIT 方法相當, 但由于CP分解方法是求解優(yōu)化問題, 需要通過不斷重復迭代直到達到收斂條件, 而ESPRIT 方法通過一次特征值分解即可求得解析解.因此無論是計算的時間復雜度還是空間復雜度, 基于ESPRIT 的方法均小于CP 分解方法.

表1 不同方法的復雜度對比Table 1.Complexity comparison of different methods.

文獻[28]中的基于CS-OMP 的2D-DOA 與載頻聯(lián)合估計方法, 利用雙L 型陣列MWC 陣列結構進行信號接收, 用OMP 算法進行參數(shù)估計,其所需的時間復雜度大約為O(NM2P2+N3+N2(M+Q)), 其中P為分割的網(wǎng)格個數(shù).空間復雜度大約為O(NP2+N2+NM) , 可見CS-OMP算法復雜度很大程度上由網(wǎng)格的大小決定, 網(wǎng)格越小, 參數(shù)估計精度雖然會變高, 但時間復雜度和空間復雜度會急速增長.將3 種方法的復雜度總結如表1 所示.圖4 為M= 3,Q= 160, 迭代次數(shù)I=50 時3 種算法的時間復雜度, 結果表明ESPRIT方法的復雜度最小, 但是從第5 節(jié)仿真實驗可以看出CP 分解方法在魯棒性方面有更好的表現(xiàn).

圖4 不同方法在M = 3, Q = 160, 迭代次數(shù)I = 50 時的復雜度對比圖Fig.4.Multiplications comparison vs.N with M = 3, Q =160, and I = 50.

5 仿真實驗

在本節(jié)中, 對本文提出的基于L 型延遲陣列MWC 結構的ESPRIT 算法、CP 分解算法以及CS-OMP 方法[28]進行對比分析, 探究3 種估計方法在不同陣元數(shù)、快拍數(shù)以及不同信噪比情況下的參數(shù)估計表現(xiàn).實驗中使用的信號均根據(jù)2 節(jié)中定義的窄帶信號模型由仿真生成.定義

分別作為載頻和DOA 參數(shù)估計準確度的評價指標[24], 其中和φ?i分別為載頻、方位角和俯仰角的估計值.

設置信號個數(shù)M=3 , 奈奎斯特頻率fNyq=10 GHz, 帶寬B= 150 MHz.偽隨機序列p(t) 每周期65 個點, 周期頻率fp= 1.1B= 154 MHz.調(diào)制載頻fi、方位角θi和俯仰角φi均在定義范圍內(nèi)隨機選取.設置每個軸的陣元數(shù)N=6 , 間距d=0.03 m.低通濾波器的截止頻率為fp/2=77 MHz,每通道采樣率為fs=fp.

首先對傳感器陣元個數(shù)對估計效果的影響進行探究.設置信噪比SNR=10 dB, 快拍數(shù)Q= 100,每個軸的天線個數(shù)從臨界最小個數(shù)N=M+1=4到N=10 遞增, 步進值為1.仿真結果如圖5 所示,可以看出3 種方法隨著天線個數(shù)的增加, 參數(shù)估計效果逐漸變好, 這主要因為傳感器數(shù)量的增加提高了系統(tǒng)對噪聲的魯棒性, 并使其能夠處理更多的源信號.同時可以發(fā)現(xiàn)在信噪比為10 dB 的情況下,CP 分解方法的估計效果要比ESPRIT 算法好, 因為CP 分解方法可以充分利用多維信息, 直接處理接收數(shù)據(jù)形成的張量模型, 具有更好的性能.結合復雜度分析的結論可知, 相比另外兩種方法,ESPRIT 方法的計算復雜度更低, 運算量較小.而基于雙L 型陣列的CS-OMP 方法估計精度受網(wǎng)格大小限制, 因此效果不如提出的兩種方法.

圖5 不同陣元個數(shù)下參數(shù)估計效果 (a) 載頻估計效果; (b) 方位角估計效果; (c) 俯仰角估計效果Fig.5.Performance of estimated parameters under different N: (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

圖6 不同快拍數(shù)Q 下參數(shù)估計效果 (a) 載頻估計效果; (b) 方位角估計效果; (c) 俯仰角估計效果Fig.6.Performance of estimated parameters under different Q : (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

圖7 不同信噪比下參數(shù)估計效果 (a) 載頻估計效果; (b) 方位角估計效果; (c) 俯仰角估計效果Fig.7.Performance of estimated parameters under different SNR: (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

接下來探究快拍個數(shù)Q的大小對參數(shù)估計精確度的影響.設置信噪比為10 dB, 陣元個數(shù)N= 10,Q的值從40 到110, 步進值為10.參數(shù)估計結果如圖6 所示, 可以看出隨著快拍數(shù)Q的增大, 各方法的信號參數(shù)估計效果變好.所提的陣列結構的總陣元個數(shù)為2N— 1=19 個, 文獻[28]中的雙L 型陣列的總陣元個數(shù)為3N— 1=29 個.由于所提方法使用陣元個數(shù)較少, 因此在快拍數(shù)較少的時候參數(shù)估計的表現(xiàn)不如CS-OMP 方法, 但是在快拍數(shù)大于80 時, 所提的兩種方法的載頻和方位角參數(shù)估計效果均比CS-OMP 方法好, 這是由于CS-OMP方法的參數(shù)估計精度受網(wǎng)格大小的限制.

最后探究系統(tǒng)信噪比大小對各個方法估計精度的影響.設置信噪比從—20 dB 到30 dB 遞增, 步進值為5 dB, 快拍數(shù)Q= 100, 陣元數(shù)N= 10.實驗結果如圖7 所示, 隨著信噪比的增大, 3 個方法的參數(shù)估計效果均越來越好.通過對比可以發(fā)現(xiàn),基于ESPRIT 的參數(shù)估計方法在信噪比小于20 dB時, 魯棒性不如CP 分解方法, 這是因為ESPRIT方法需要先計算通道采樣值間的協(xié)方差函數(shù), 而CP 分解方法可以直接處理多維的信息, 具有更好的性能.由于CS-OMP 方法陣元個數(shù)比所提的陣列結構多, 因此在信噪比較低時效果較好.但在信噪大于5 dB 時, CS-OMP 方法的精度受網(wǎng)格限制, 參數(shù)估計的誤差較大.同時CS-OMP 算法只能重構θi ∈(0°,90°) 內(nèi)的角度, 提出的兩種方法均可以恢復在 (-90°,90°) 之內(nèi)的角度, 范圍較大, 更加適合于實際應用.

從仿真結果和復雜度分析結果可以看出, 提出的基于ESPRIT 的方法相對計算量更小, 實時性更高, 但估計效果比CP 分解方法稍差.提出的CP 分解方法的相對魯棒性更好, 但計算量大約是ESPRIT 方法的IM倍.CS-OMP 算法在信噪比較低時對2 D-DOA 的估計效果較好, 但信噪比較高時精度受限于網(wǎng)格大小, 計算復雜度隨網(wǎng)格的增加會呈指數(shù)增長.同時提出的方法可以估計的方位角范圍為 (-90°,90°) , 而CS-OMP 算法可估計范圍僅為 ( 0°,90°) , 因此提出的算法更加具有實用性.

6 結 論

本文提出了一種L 型延遲陣列MWC 的接收結構, 通過增加延遲結構直接估計信號載頻, 進而估計2D-DOA 參數(shù), 無需額外的參數(shù)配對, 避免了配對引入的誤差和復雜度的提升.接著在該陣列結構特點的基礎上分別提出了兩種參數(shù)估計方法, 一種基于ESPRIT 方法, 另一種基于CP 分解技術,針對兩種方法進行了復雜度分析和仿真, 探究了陣元間距、快拍數(shù)以及信噪比對方法估計精度的影響, 并對提出的兩種方法以及文獻[28]中的基于網(wǎng)格的CS-OMP 方法進行對比分析, 從實驗結果中可以總結出所提方法均可從欠奈奎斯特樣本中估計信號的載頻和二維DOA 參數(shù).基于ESPRIT方法的計算量較小, 但不足之處是低信噪比條件下表現(xiàn)欠佳, 適用于主動雷達等要求實時性的場合;而基于CP 分解技術的方法魯棒性更好, 適用于低信噪比的場合, 例如被動雷達、遙感等被動感知的目標探測領域, 但所需計算量較大, 實時處理困難.本文研究了空間電磁波信號在欠采樣條件下的載頻和二維到達角估計, 如果在已知信號的某些特征的情況下, 還可以進一步改進來減小算法復雜度或提升抗噪性等.未來在減小算法復雜度方面可以結合信號已知信息來合理的選取初始值, 或者利用梯度下降等方法加快收斂速度, 減少循環(huán)次數(shù); 在增強參數(shù)估計精度方面可以利用先驗信息設計合理的去噪算法; 也可以利用互質(zhì)陣列、嵌套陣列等稀疏陣列來減少接收陣列的陣元個數(shù), 降低硬件復雜度和成本.