時?空欠采樣下多個信號的頻率和DOA 聯合估計

郜麗鵬，李佳林

哈爾濱工程大學信息與通信工程學院，黑龍江哈爾濱150001

隨著電子戰技術的日益發展，電磁環境也變得越發復雜，針對多信號的多參數估計是雷達、聲吶、電子戰等領域的研究熱點，能夠快速、準確地偵察檢測多個雷達信號的頻率、方向信息，對于電子對抗、雷達干擾技術而言至關重要[1?2]。在經典的空間譜估計算法MUSIC和ESPRIT的基礎上提出的改進算法可以完成對多個信號的頻率和DOA聯合估計及配對[3?4]，但是這些算法的運算量都較大，且均需要滿足信號在時域和空域上的理論前提。時域方面，香農采樣定理表明，奈奎斯特采樣率是被采樣信號的最大頻率的2倍，比如為了適應電子戰2~18GHz的工作頻段，需要的采樣頻率會很高，必須相應地增加模數轉換器（analog-to-digital converter，ADC）的硬件成本和整體結構的功耗，甚至在某些特定情況下是不可實現的。同樣，在空域方面若要適應電子戰的工作頻段，為滿足相鄰陣元的間距不得大于信號波長的二分之一的采樣定理，需要密集的安裝陣元天線，這樣會帶來相鄰陣元間的互耦和硬件安裝困難等問題，進而影響測向精度。因此，亟需尋求一種在時?空域欠采樣條件下，能夠對多信號多參數（頻率、DOA）聯合估計的新方法。

文獻[5?6]提出了一種基于閉式中國余數定理（CRT）的頻率、DOA聯合估計方法，建立了基于閉式CRT的頻率和DOA估計的模型，利用不同速率ADC采樣樣本（離散傅里葉變換，Discrete Fourier Transform）變換后搜索的譜峰位置構建余數完成頻率估計，通過譜峰位置求得信號到達每個陣元的初相位及其相位差構建余數完成DOA估計，實現了高精度的聯合估計。但該算法是針對單目標信號的參數估計，限制了其應用場合。文獻[7]在閉式CRT的頻率和DOA估計模型的基礎上，提出了針對多目標信號的參數估計，利用信號的幅度不一致的信息，在多個信號經不同速率ADC采樣樣本DFT變換后得到的多個譜峰中，尋找屬于同一信號的一組譜峰以構建頻率余數和相位余數完成多個信號的頻率和DOA聯合估計，但該方法在有2個及2個以上信號的幅度一致的條件下會失效。文獻[8]提出了一種能夠重構多個整數的中國余數定理，并且將其應用在窄帶信號的頻率估計上，但其未考慮當信噪比較低時或樣本數較少時譜峰位置估計不準確的問題，導致其測頻精度的降低。

本文建立了一種基于重構多個整數的CRT算法和AM估計器的多信號的頻率、DOA參數聯合估計模型，利用多整數估計的CRT算法對稀疏線陣欠采樣后的少量樣本做DFT、譜峰校正，余數構建等，完成多個信號的頻率估計，并根據估計完成的頻率值和已經構建的余數信息在多組相位差中確定每個信號對應的一組相位差，完成頻率和相位差的參數配對，進一步使用閉式CRT算法魯棒地估計出多個信號的DOA值。本文算法只需對各個陣元接收到的多個信號做一次并行欠采樣，并且譜校正后得到的譜峰信息可同時用于頻率估計和DOA估計，滿足目前復雜電磁環境下需要快速參數估計的需求。仿真結果證明了本文算法的高精度和低信噪比下的頑健性。

1 重構多個整數的頑健CRT算法

模數 {m1=ΓM1,m2=ΓM2,···,ms=ΓMs}由 s個整數構成，其中 M1,M2,···,Ms互素。假設 n個 整數X1,X2,···,Xn均勻分布在( L,G]內 ，其中G

1.1 余數誤差情況下的頑健重構算法

首先考慮只有一個帶重構整數的情況，先前許 多 CRT算 法[9?11]的 理 論 前 提 為|r?il?ril|=然而在實際應用中，有可能會發生的 情 況 ， 這 會 導 致進一步使得CRT算法失效。文獻[12]給出了這種問題的解決方法，其主要思想是將每個模數中的誤差 ?il調整為一致。首先根據 ml=ΓMl，可以得到：

式中 β為整數，因此

可以發現，只要重構出 Q就可以求得待重構的整數 X。接下來將重點放在如何魯棒、頑健地重 構 整 數 Q 上 。分 別用 a和 b表 示集 合最小和最大的元素。接下來分以下幾種情況進行討論：

1.2 重構多個整數的算法

對于多個整數的情況，為了方便討論，假設待重構整數的集合 {Xi}(i=1,2,···,n)和模數的集合{ml|l=1,2,···,s}都 是以升序排列的，會得到 n×s個余數 {ri,l}，如文獻[13]所提到的，2個及2個以上的整數重構的難點是整數和余之間的對應關系是未知的，即由每個模數得到的 n個 余數與 n個整數的對應關系是未知的，所以無法使用CRT算法進行重構。

對于這種問題，首先概述主要的算法思想。由于最大的模值 ms>d+2δ，所以可以在區間[X1? δ,Xn+δ]找 到 X ，并且通過余數集 r?is和 ?X?ms正確地估計 {Xi?X}的值。因∑此，選擇估計待重構整數集的平均值 ˉ，即X同樣地，在整數集與余數集對應關系未知的前提下，可以得到每個整 數 模 除 以 同 一 模 值 的 余 數 的 和并通過 r?和上述的CRT算法頑健的重l構整數集的和為了頑健地重構需要假設 M并且屬于同一模值的余數誤差的總和需要滿足

接下來就是如何通過 X?ˉ 和 r?is頑健地重構各個整數 Xi。由于考慮到若有3個整數A、 B 和 m， 如果 |A?B|

引理 對于任意的 f∈R ， 在區間 [f,f+ms)內至多 有 n個 元 素 κp,κp+1,···,κp+n?1， 其 中p∈{1,2,···,h?n+1}。

關于引理的證明見文獻[8]。

2 欠采樣的頻率估計方法

首先給出系統模型，如圖1所示，由 s(s>2)個陣元天線組成非均勻線陣，假設有 m個同時到達的遠場信號，不失一般性，將其表示為單頻復指數形式：

式中： Ai和 fi分別為各個信號的幅度和頻率；φi為信號的初相位； i=1,2,···,m。假設各個信號的頻率fi都不相同，由于噪聲和陣列陣元的間距，第l個陣元接收到的信號為：

式中： φil為 第i個 信號到達第l個陣元的初相位；ζl(t)為 高 斯 白 噪 聲 ； l=1,2,···,s。 從任 一 時刻 開始，s個陣元分別以 fs1～ fss的采樣速率對同時到達的 m個 信號進行欠采樣，并且得到了 s路樣本，假設每路的樣本數為 N ，則第l個陣元得到的樣本序列為：

圖1 稀疏陣列模型

對 s路欠采樣樣本序列做DFT，并進行譜峰搜索，得到譜峰位置 kil，則各個信號的頻率 fi可表示為：

式中： nil為 模糊倍數； δil為頻偏小數，且 |δil|?0.5，也是影響測頻精度的重要因素，δil可由AM估計器進行校正，進而得到更精確的余數；i=1,2,···,M;l=1,2,···,s。可將式（3）轉變為中國余數定理模型

式中：待測的信號頻率 fi對應CRT中的被除數 Xi；采樣頻率 fsl對 應CRT中的模數 ml=ΓMl，Γ為最大公約數， Ml為互素的整數；i =1,2,···,M ；l =1,2,···,s。而重構所需的余數為

對于 m個 信號的情況，在 s路樣本序列DFT變換后，都會搜索得到 m個 譜峰，而在 m個信號的幅度有2個及2個以上是相同的情況下，不能在每路 m個 譜峰中尋找對應同一信號的 s個譜峰，即不能構造出 m組 余數 {ri1,ri2,···ris}， 進而分別進行 m次閉式CRT算法重構出 m個信號的頻率。

考慮到第1章提出的重構多個整數的CRT算法，其前提條件為 ms>d+2δ，即模數組中的最大模值需要大于多個整數中最大值與最小值的差值，結合本節的頻率估計模型，即最大的欠采樣速率 fss需 要大于信號帶寬（ fm?f1），所以該CRT算法適用于窄帶信號中的多個頻率估計，算法流程如下：

1）對 s路欠采樣樣本做 N點的DFT變換，譜峰搜索后得到譜峰位置 kil， i=1,2,···,m ， l=1,2,···,s，即每路樣本得到 m個譜峰位置。

2）校正每路樣本的譜峰位置對應的頻偏小數，對于第l 路欠采樣樣本，將 m頻偏小數的值進行初始化=0， i=1,2,···,m，按式（5）進行計算得到迭代值

式中： q=±0.5； Re(·)表示取其實部。將其迭代2次以上的收斂值做為頻偏小數的估計值，i=1,2,···,m ，l=1,2,···,s。

3）將譜峰位置 kil、頻偏小數代入式（4）得到校正后帶有誤差頻率余數估計值 r?il， i=1,2,···,m，l=1,2,···,s。

4）計算對應同一模數 ml的余數估計值的和

3 欠采樣的DOA估計方法

3.1 基于CRT算法的DOA估計模型

陣列結構如圖1所示，由相位干涉儀原理可知，由信號到達陣元的波程差會產生相位差，在任一時刻，陣元 j 和陣元 j+1的 接收到的信號 Si的相位差為：

式中：i=1,2,···,m；j=1,2,···,s?1；λi為 信號Si的波長； θi為 信號Si的 入射角度；dj為陣元間距，當則會產生相位模糊問題，即：

式中：nij為 模糊倍數； ?φij為真實情況下測得的相位差。由式（6）、（7）可得：

式中：C為非負常數；Mθ為任意正整數； Γθ1～ Γθ(s?1)為互素的整數。令待重構的整數Q別定義模為相位差構造的對應CRT中的余數為

所以，只要我們得到每個信號在陣元間的相位差，就可以通過CRT算法重構出Qθi，在頻率值和相位差余數配對后，通過式（9）求得DOA的估計值，其中 λ0的值可由第2章中測得頻率結果得到。

3.2 相位差余數與信號頻率的配對

考慮到DFT變換后譜峰位置對應的相位為信號Si到 達陣元l的 初相位，記為 φpil。由第2章分析可知，對于m個信號同時到達s個陣元的情況，可以通過s路樣本序列得到m×s個譜峰位置。同樣我們可以得到m×s個初相位，即每個陣元可以得到m個初相位。只有在相鄰陣元的初相位中，選擇屬于同一信號頻率的2個初相位，才能得到正確的相位差，進而得到正確的DOA估計值。

不同整數Xi模 除同一模值m得 到的余數為ri，定 義P0為 余 數ra和rb不 相 等 的 概 率 ， 其 中i=1,

表1中s為模值的個數。由式（10）可知，當模值m相 對較大時，余數ra和rb不相等的概率極大，特別是對應第2章中的欠采樣采樣頻率，一般都在100MHz以上，此時P0幾乎為1。所以，可以認為，對于任一陣元的樣本序列sl(n)，通過DFT變換后的譜峰位置得到的頻率余數值

式中i=1,2,···,m。

表1 s個陣元采集m個信號得到的余數及對應初相位

如表1所示，每列的余數都是不相同的，同時，我們由譜峰位置得也可以到對應的初相位 φpil。

在頻率估計后，m個信號的頻率估計值f?i和s個陣元的采樣頻率fsl均為已知的，故可以通過f?imodfsl，l=1,2,· · ·,s， 得到信號Si的s個 余數ril，并依次用余數ri1～ris，在表中每列找到r?i1～r?is對應的的位置，進而得到信號Si到 達s個陣元的s個初相位 φpi1～ φpis， 并 求 得 信 號Si的s?1個 相 位 差 ?φij，j=1,2,···s?1。 同理，依次可以得到m個信號的m×(s?1)個相位差。所以此時第1章的重構多個整數 的C RT算法適用，而考慮到該算法的理論前即要求通過譜峰位置得到所有信號到達同一陣元的初相位的誤差的和小于DOA估計模型中的這個條件對于相位測量是比較苛刻的。但與頻率估計模型不同的一點是，對于DOA估計來說，經過上述方法的參數配對，屬于同一信號帶有DOA信息的待重構整數Qθi的

s個 相 位 差 余 數個待重構整數與余數集的對應關系是已知的。所是已知的，也就是每以我們可以分別使用m次傳統的閉式CRT算法求得m個 整數Qθi的值，這樣相較于直接重構多個整數的算法，我們就把相位測量誤差的限度降低了m倍。故DOA估計算法流程如下：

1）對于第l路樣本，根據譜峰位置對應的幅值求得個信號到達陣元l的 初相位 φpil,i=1,2,···,m，并由頻率估計過程中得到的頻偏小數通過式（11）對信號初相位進行校正：

式中 i=1,2,···,m。

式中 j=1,2,···,l?1。

3）通過相位差 {?φij,j=1,2,···,l?1}構成的相位差余數 {rθij,j=1,2,···,l?1}和閉式CRT算法求得信號 Si的 DOA估計值θ?i， 分別對 m組相位差余數使用CRT算法即可求得 m個信號的DOA估計值。

4 仿真實驗

設有3個遠場信號同時到達陣列，陣列由4個陣元組成，假設信號的頻率在1 500～2150MHz，根據本文提出的欠采樣下的頻率和DOA聯合估計方法，最大模值需大于多個整數中最大值與最小值的差的條件。陣元中最大的采樣頻率 fs4應設置為 650MHz，選取對應CRT算法參數中的最大公約數 Γ=50×106，則一組互素模數中最大模數M4=13。為滿足時域欠采樣的條件，選取3個小于 M4的模數構成互素的模數組[M1,M2,M3,M4]=[3,5,11,13]，所以對應的4個陣元ADC的采樣頻率分別為 fs1=150MHz、 fs2=250MHz、 fs3=550MHz、fs4=650MHz。顯然陣元ADC的采樣速率遠小于信號的頻率，屬[于時域欠]采樣。非負常數 C取為0.1，互質參數 Γθ1,Γθ2,Γθ3=[3,5,7]， 即陣元間距為[d1,d2,d3]=[0.35m,0,21m,0.15m]，小于入射信號的半波長，屬于空域欠采樣。

本節實驗采用檢測概率 pd和均方根誤差（root mean square error,RMSE）對算法的頑健性、抗噪性以及參數測量精度進行考量。假設每個陣列的接收特性相同并且是相互獨立的，欠采樣樣本序列的長度 N=2000，噪聲是均值為零的高斯白噪聲，在每個信噪比（SNR）條件下進行1000次Monte Carlo實驗。

4.1 仿真實驗1

在時域欠采樣的條件下，單獨對多信號的頻率估計的頑健性、抗噪聲性能進行考察。在該仿真中，固定3個信號的入射角均為38°，設信號的頻率 fi,i=1,2,3，在15002~2150MHz內隨機選取。

圖2 不同信噪比下的頻率檢測概率

4.2 仿真實驗2

接下來驗證多個信號的頻率、DOA聯合估計算法的精度性能。該仿真中，固定3個信號頻率和 入 射角 度 Si[fi,θi]分 別 為且同時到達陣列，基于AM估計器的頻率估計的RMSE仿真曲線如圖3所示。從圖3可以看出，本文提出的基于多整數估計CRT的頻率估計算法在各個信噪比下均有著較高的精度，當信噪比大于10dB后，3個信號的頻率估計誤差均小于1kHz，即小于信號頻率的

圖3 不同信噪比下3個信號頻率估計的均方根誤差

且由于信號頻率估計的精確度較高，我們通過頻率估計值和ADC采樣率可以得到準確的余數，并通過頻率估計中譜峰位置得到的余數正確配對同一信號到達4個陣元的初相位。圖4為在信號和相位參數配對后，基于閉式CRT算法分別對3個信號的DOA估計的RMSE曲線?？梢钥闯觯诟鱾€信噪比條件下均表現出較高的DOA估計精度，并且在(0,90°]測角范圍內，信號入射角度越大，DOA估計誤差越大，與理論DOA估計誤差一致。

圖4 不同信噪比下3個信號DOA估計的均方根誤差

4.3 仿真實驗3

考察頻率和DOA聯合估計時，DOA估計的頑健性和抗噪聲性能。在該仿真中，固定3個信號的頻率分別為1.573、1.900、2.016GHz，信號角度在(0,90°]內隨機選取。若DOA估計值則認為檢測成功；否則為檢測失敗。在DOA估計過程中，采用重構多個整數的CRT和獨立3次使用閉式CRT2種算法做比較，得到的檢測概率隨信號比SNR變化的仿真曲線如圖5所示。正如3.2節的討論，重構多個整數的CRT算法在DOA估計時對允許的相位測量誤差的限度較低，其在低信噪比條件下的性能不佳，而在參數配對后單獨3次使用閉式CRT算法在信噪比大于?9dB后，檢測概率均能達到90%以上，可以看出其較好的頑健性。

圖5 不同信噪比下DOA的檢測概率

4.4 仿真實驗4

最后對比基于MUSIC的空時二維譜頻率和DOA聯合估計算法的性能和算法耗時。與仿真2相同，固定3個信號頻率和入射角度 Si[fi,θi]分別為33°]，且同時到達陣列。對于空時二維譜算法，陣列陣元間距保持不變，4個陣元的采樣速率均為650MHz，滿足時域和空域為欠采樣的條件。角度 在(0,90°]的 范 圍 內 以0.1°為 步 長，頻 率 在(1500MHz,2150MHz]以 0.1MHz為步長搜索得到二維譜如圖6所示。

圖6 基于MUSIC的空時二維譜

可以看出由于時域和空域的欠采樣，空間二維譜會出現偽峰，無法對3個入射信號進行正確的頻率和DOA估計。且基于MUSIC的空時二維譜的估計算法的仿真時間為71.985s，并且會隨著搜索步長的減小而增大。而本文提出的時?空欠采樣下的頻率和DOA聯合估計算法的仿真時間為0.037s，極大地減少了算法運算時間，適合快速時變的工作場景。

5 結論

本文提出了一種在時?空欠采樣條件下，能夠對多個入射信號進行頻率和DOA聯合估計的算法，并且給出了對應的陣列模型及算法過程。

1）解決了頻率估計時，多組頻率余數與多個待估計信號頻率對應關系未知情況下，CRT算法失效的問題；

2）解決了多個信號頻率值與信號到達陣列初相位的參數配對問題，使得在多信號DOA估計過程中，可以獨立多次地使用CRT算法進行單個信號DOA的估計，保證了其頑健性與精確度。

本文提出的算法只需對各個陣元做一次并行欠采樣，即可完成頻率和DOA的聯合估計，有著低運算量和高估計精度的優點，適合雷達、電子戰、無線通信等對測量實時性，測量精度有要求等工作場合，具有較廣泛的應用場景。