王瑞，馬艷

(西北工業(yè)大學(xué) 航海學(xué)院，陜西 西安710072)

0 引言

寬帶線性調(diào)頻(LFM)信號在聲納、魚雷自導(dǎo)、通信和海底勘測等探測設(shè)備中有著廣泛的應(yīng)用，針對此類信號的波達(dá)方向(DOA)估計問題也日益受到人們的重視［1-2］。但是由于寬帶LFM 類是一種典型的非平穩(wěn)信號，現(xiàn)有的非平穩(wěn)LFM 信號DOA估計算法都是建立在窄帶信號的假設(shè)基礎(chǔ)之上，因此它們不適合于寬帶信號。而分?jǐn)?shù)階傅里葉變換(FRFT)作為傳統(tǒng)傅里葉變換(FT)的推廣，可以借助于快速傅里葉變換(FFT)來實(shí)現(xiàn)，計算更方便。其次是對LFM 信號的能量在最佳旋轉(zhuǎn)角度具有很好的聚集性［3］。因此，將傳統(tǒng)的高分辨DOA 估計方法——ESPRIT 和MUSIC 算法與FRFT 相結(jié)合在LFM 的DOA 估計領(lǐng)域中引起了廣泛重視［4-9］。其中MUSIC 類算法需要依據(jù)信號的中心頻率構(gòu)造出FRFT 域的陣列流行向量和相關(guān)矩陣［7-8］。這些方法都是以觀測時間等于信號時寬為前提的，此時信號的中心頻率就是信號在t =0 時刻的頻率，基于FRFT 直接估計出的信號頻率直接對應(yīng)于信號的中心頻率［10］。

但在主動聲納和雷達(dá)中，常用發(fā)射脈沖信號來完成對目標(biāo)距離、速度和DOA 的檢測和估計。觀測時間內(nèi)是否包含目標(biāo)回波信號、包含回波信號的長度等都是未知的，對于這種脈沖LFM 不完全包含于觀測時間(包含的脈沖信號時寬不等于觀測時間)的情況，由于中心頻率不再對應(yīng)于t=0 時刻的中心頻率，因此基于FRFT 方法估計的中心頻率將有一定的誤差，從而給方位估計帶來較大的偏差。對于這個問題，文獻(xiàn)［11］從理論上分析了信號時寬不等于信號觀測時間時對分?jǐn)?shù)階譜位置的影響，并進(jìn)行了修正，但其中要求的信號確切位置在實(shí)際應(yīng)用中是不知道的，因此在實(shí)際中不太適用。

文中針對雷達(dá)、聲納和魚雷自導(dǎo)中這種觀測時間與信號時寬不相等的LFM 脈沖信號(忽略多普勒頻移)，提出了一種新的基于FRFT 的中心頻率估計方法，并在此基礎(chǔ)上對分?jǐn)?shù)階傅里葉域的陣列流行向量進(jìn)行了修正，進(jìn)而實(shí)現(xiàn)了對基于分?jǐn)?shù)階傅里葉域MUSIC 算法的寬帶LFM 信號DOA 估計性能的改善。并分析了信號寬度和觀測時間寬度比值對分?jǐn)?shù)階傅里葉域信噪比及估計性能的影響。該方法提高了基于FRFT 的寬帶LFM 脈沖信號的參數(shù)估計準(zhǔn)確性，并進(jìn)一步擴(kuò)大了用FRFT 處理寬帶LFM 的適用范圍。

1 分?jǐn)?shù)階傅里葉變換

FRFT 有很多種不同的定義方式，其中一種是從線性積分的角度給出了FRFT 的基本定義。信號x(t)的p 階FRFT 定義為

式中:α = pπ/2;Fα［·］為FRFT 的算子符號;Kα(t，u)為FRFT 的變換核，

式中:Aα=

在實(shí)際工程中，信號都是經(jīng)過采樣后的離散時間信號，通常采用離散FRFT(DFRFT)，Pei 和Ozaktas 等在這方面做出了卓越的貢獻(xiàn)［12-16］。Ozaktas等在文獻(xiàn)［16］中提出了一種計算DFRFT 的快速算法，這種方法是對連續(xù)FRFT 的核函數(shù)在FRFT 和時域直接采樣得到的，該算法將FRFT 分解為先乘以LFM 信號，然后和LFM 信號卷積，最后再乘以一個LFM 信號。在利用該算法對信號進(jìn)行離DFRFT 之前必須對信號要進(jìn)行量綱歸一化處理。對離散信號直接用該算法，暗含著取信號的時域?qū)挾葹椋?Td/2，Td/2］，其中Td為觀測時間，而頻率范圍為［-fs/2，fs/2］，其中fs為采樣頻率，則時寬帶寬積為N=Tdfs，得到的尺度因子S 和歸一化寬度U［17］分別為

則u 域的范圍為(-U/2，U/2)，采樣間隔為Δu =1/U=(Tdfs)-1/2.

2 基于FRFT 的MUSIC 算法

鑒于FRFT 對LFM 信號良好的能量聚焦性能，對于多個寬帶LFM 的方位估計提出了用分?jǐn)?shù)階傅里葉域的MUSIC 算法［7-8］，在分?jǐn)?shù)階傅里葉域構(gòu)造相關(guān)矩陣和陣列流行向量，實(shí)現(xiàn)對多個寬帶LFM 的DOA 估計。

2.1 均勻陣列信號模型

在二維平面內(nèi)M 個陣元以等間隔d 排列在一條直線上，假設(shè)信號為遠(yuǎn)場源，可看成是平行入射，入射波與陣列的法線夾角θ 為DOA，第一個陣元作為參考陣元。假設(shè)有L 個寬帶LFM 入射到均勻線列陣，第j 個寬帶LFM 為sj(t)= Cjexp［jπ(2fjt +Kjt2)］，j =1，2，…，L，式中:Cj是第j 個信號的幅度;fj為第j 個信號的中心頻率;Kj為第j 個信號的調(diào)頻率。

則第i 個陣元上收到信號為

式中:τij=(i-1)dsin(θj)/c 為第i 個陣元相對參考陣元的時延，θj為第j 個LFM 的入射角，c 為海洋中的聲速，為1 500 m/s.假設(shè)信號相互之間、信號與噪聲中間、各個陣元噪聲之間都不相關(guān)。由(4)式可以看出，LFM 信號的時延僅僅改變了信號的中心頻率和初始相位，但調(diào)頻率保持不變。

2.2 基于FRFT 的MUSIC 算法

對參考陣元的信號為第j 個LFM 信號sj(t)，由FRFT 的定義式(1)式，得sj(t)關(guān)于α 角的FRFT 為

式中:T 為信號時間。由(2)式可以看出，在αj=-arccot Kj時具有最好的能量聚集特性。此時sj(αj，顯然，sj(αj，u)在uj=fj/csc αj具有最大值:sj(αj，uj)=最大值點(diǎn)的幅值［18］為

信號分?jǐn)?shù)階傅里葉域幅度的最大值與信號的時間寬度、信號幅度呈正比，而與|sin αj|呈反比。由FRFT 的時移性質(zhì)［3］可得第i 個陣元上信號的最大值［8］為

式中:

式中:Xαj= ［Xαj1(u)，Xαj2(u)，…，Xαji(u)，…，為M×N 維的陣元接收信號的αj階FRFT;Sαj=［Sαj1(u)，Sαj2(u)，…，Sαjj(u)，…，SαjL(u)］T為L× N 維信號的αj階FRFT;Aαj=［Aαj1，Aαj2，…，為M×L 維αj階分?jǐn)?shù)階傅里葉域的陣列流行向量;Nαj=［Nαj1(u)，Nαj2(u)，…，NαjM(u)］T為M×N 維噪聲信號的αj階FRFT.

則Xαj的自協(xié)方差矩陣為

式中:Us是由大特征值對應(yīng)的特征向量構(gòu)成的信號子空間;UN是由小特征值對應(yīng)的特征向量構(gòu)成的噪聲子空間;Es是由大特征值構(gòu)成的對角矩陣;EN是由小特征值構(gòu)成的對角矩陣。

根據(jù)MUSIC 算法得到第j 個信號基于FRFT 的MUSIC 空間譜為

對(12)式進(jìn)行一維搜索，即可得到第j 個脈沖信號的入射角。從(8)式可以看出陣列流行向量不僅和聲程差τij有關(guān)，還和信號的中心頻率fj及最佳變換角度αj= - arccot Kj有關(guān)。在陣列結(jié)構(gòu)確定后，聲程差僅與信號的入射角有關(guān)。

由(8)式可以看出，LFM 的FRFT 方向向量與信號的中心頻率fj有關(guān)，對于觀測時間不等于信號時寬的脈沖信號，由于信號成分不完全包含于觀測時間內(nèi)，因此信號經(jīng)過FRFT 后估計得到的頻率不再是脈沖信號的中心頻率，而是對應(yīng)于0 時刻的頻率，這就給傳統(tǒng)的基于FRFT 的MUSIC 算法帶來估計上的誤差，為此文中提出了一種針對脈沖LFM 的FRFT 的基于MUSIC 算法，對估計的中心頻率進(jìn)行修正，重新構(gòu)造分?jǐn)?shù)階域的方向向量。

3 基于FRFT 的脈沖信號MUSIC 算法

3.1 脈沖信號的中心頻率估計

LFM 的時頻分布線如圖1所示，信號的時寬為［-T/2，T/2］，中心頻率為f0，調(diào)頻率為K.觀測時間為［-Td/2，Td/2］，當(dāng)T =Td，并且信號的起點(diǎn)和觀測時間的起點(diǎn)相同時，信號的中心頻率正好是t=0 時刻點(diǎn)的信號頻率，可以用信號αj= -arccot Kj階FRFT 幅度的最大值點(diǎn)的uj估計:fj=ujcsc αj.但是當(dāng)LFM 脈沖信號沿t 軸平移τ，即信號不完全包含于觀測時間［-Td/2，Td/2］內(nèi)的情況，如圖1所示，在觀測時間內(nèi)LFM 時寬為T -τ，此時信號的中心頻率變?yōu)閒0c，而經(jīng)過FRFT 后估計出的信號頻率為t=0 時刻的頻率f0b(0 ＜τ ＜T/2)或f0a(τ ＞T/2)，二者已不再相等。另一方面，延遲后LFM 的調(diào)頻率保持不變，也就是FRFT 幅度出現(xiàn)最大值的αj不變。

由圖1中的幾何關(guān)系得出:

當(dāng)τ ＞0 時，

式中:當(dāng)τ ＞T/2 時，fH為在觀測時間內(nèi)時延之后信號的最高頻率，fL為在觀測時間內(nèi)時延之后信號的最低頻率。如圖1所示，雖然t=0 時刻沒有信號成分，但是由于p 階FRFT 相當(dāng)于信號的Wigner 分布在時間-頻率平面上順時針旋轉(zhuǎn)α=pπ/2，所以其在u 軸對應(yīng)u0a，根據(jù)(13)式能得出t =0 時刻的頻率f0a，因此(15)式仍適用。當(dāng)τ ＜0 時，同理可得

3.2 修正后基于FRFT 的脈沖信號MUSIC 算法

在實(shí)際應(yīng)用中，通常采用DFRFT.文中選用Ozaktas 提出的計算DFRFT 的快速算法［16］。對參考陣元信號sj(t)做DFRFT 得到sj(α，u)，用多分辨方法［19］在(α，u)平面對sj(α，u)譜峰搜索，得到峰值點(diǎn)對應(yīng)的位置(αj，uj).由于DFRFT 算法在對信號進(jìn)行DFRFT 之前必須對信號要進(jìn)行量綱歸一化處理，因此uj將變?yōu)檫M(jìn)行量綱歸一化后的值，結(jié)合(3)式可得做N 點(diǎn)DFRFT 后LFMsj(t)的中心頻率為

考慮到量綱歸一化的影響，(8)式變?yōu)?/p>

由圖1可以看出，脈沖信號平移后，信號的中心頻率不再對應(yīng)于t =0 時刻。由(16)式估計出來的頻率j不再是第j 個脈沖信號的中心頻率fj，變?yōu)閳D1中所示的f0bj.而第j 個脈沖信號中心頻率則對應(yīng)于圖1中所示的f0cj.因此需對于(17)式中的中心頻率進(jìn)行修正。

圖1 信號的瞬時頻率Fig.1 Instantaneous frequency of signal

結(jié)合(14)式得修正后的中心頻率

將(18)式代入(17)式得

故修正后的方向向量:

算法具體計算步驟如下:

1)對參考陣元上的接收信號進(jìn)行DFRFT 并進(jìn)行二維搜索，得到第j 個脈沖信號出現(xiàn)，峰值坐標(biāo)

2)由(18)式對中心頻率fj進(jìn)行修正，得到修正后的f0cj.由(19)式構(gòu)造第j 個脈沖信號的方向向量。

3)由(10)式構(gòu)造第j 個信號的分?jǐn)?shù)階相關(guān)矩陣RαjXX.

4)對RαjXX進(jìn)行特征分解，得到噪聲子空間EN.由(20)式求得基于FRFT 的MUSIC 算法的空間譜P(θ)cor.

5)找出極大值對應(yīng)的角度θcor就是信號的入射方向。

6)存在多個脈沖信號時，重復(fù)步驟1 ～5，分別得到各個信號的DOA 估計。

3.3 脈寬對估計結(jié)果的影響

幅度為C 的LFM 延遲τ 后，在觀測時間T 中包含的信號寬度為T - τ，那么它的“最佳”階αj=-acrcotKj階FRFT 的幅度為一個Sinc 函數(shù)，最大值的幅度為

從(21)式可以看出信號的幅度與延遲τ 成反比，τ 越大，觀測時間中包含的信號成分越少，幅度就越低，也就是信號的能量越小。而在觀測時間內(nèi)，噪聲的寬度始終是T，對于均值為0，方差為σ2的白噪聲，它的αj= - acrcotKj階FRFT 后依然是白噪聲，它的均值為0，方差為所以FRFT 域的信噪比將隨著延遲的增大而降低，因而將會對估計精度產(chǎn)生影響。

3.4 運(yùn)算量討論

與傳統(tǒng)的時域MUSIC 算法相比，基于FRFT 的MUSIC 算法增加了第1 步對參考陣元進(jìn)行DFRFT并搜索其峰值點(diǎn)的位置(αj，uj)，并估計信號成分的中心頻率，和對其他陣元完成αj的DFRFT 并構(gòu)造分?jǐn)?shù)階相關(guān)矩陣RαjXX，這增加了一定的運(yùn)算量。但是對LFM 來說，可以用多分辨搜索方法來實(shí)現(xiàn)第1 步的搜索［19］，并利用基于FFT 的快速DFRFT來完成對各個陣元上信號的分析［16］。長度為N 的信號，每進(jìn)行一個角度的DFRFT，所需的計算量為O(Nlog2N).搜索出參考陣元的峰值點(diǎn)位置后，其他陣元的信號只需要做這個角度的DFRFT 就可以了。

4 仿真結(jié)果

考慮陣元個數(shù)為6 的均勻線陣。信號中心頻率為10 kHz，帶寬為5 kHz，時寬為0.017 s，采樣頻率為60 kHz，快拍數(shù)N 為1 023，入射角為20°，信噪比SNR=20 dB，信號時延為τ =0.008 s.圖2給出了參考陣元上接收信號的α 角的DFRFT 仿真結(jié)果。得延遲后其中心頻率理論值為f0=8.827 0 kHz，由其峰值位置可得未修正的信號中心頻率估計值f0b=7.651 2 kHz，修正后的中心頻率估計值為f0c=8.825 6 kHz，時延估計值為τ=0.008 s.圖3是在-90°～90°范圍內(nèi)以0.01°為間隔得到的未修正中心頻率的傳統(tǒng)算法和本文算法的DOA 空間譜pθ.由其峰值可得傳統(tǒng)算法入射角估計值^θ =23.35°和修正后的入射角估計值^θcor=20.04°.從圖3中可以看出，時延后由于信號不完全包含于觀測時間使得中心頻率偏移，傳統(tǒng)基于FRFT 的MUSIC 算法很難準(zhǔn)確地估計出入射角。但本文算法大大提高了此類脈沖信號入射角估計的準(zhǔn)確性。

圖3 信號DOA 估計的空間譜Fig.3 Estimated spatial spectrum of DOA

其他條件同上，在不同信噪比下做100 次Monte Carlo 獨(dú)立實(shí)驗(yàn)，結(jié)果如表1所示。從表1可以看出，在延遲較小時(τ =0.008 s)，盡管未修正的傳統(tǒng)FRFT 估計DOA 算法，其估計結(jié)果的方差與修正后估計的方差相差無幾，但是估計的均值和實(shí)際的入射角有較大差異;當(dāng)延遲較大時(τ =0.015 s)，未修正的算法不僅估計的均值偏差較大，而且方差也隨之增大許多，但是修正后的估計均值與理論值相差無幾，方差增大了。這正如3.3 節(jié)中所述，隨著觀測時間中包含的信號脈寬的減少，F(xiàn)RFT 域的信噪比下降有關(guān)。

表1 不同SNR 和時延情況下的估計誤差比較Tabl.1 Estimated errors

5 結(jié)論

針對雷達(dá)、聲納和魚雷自導(dǎo)中觀測時間與信號時寬不相等的LFM 脈沖信號導(dǎo)致對寬帶LFM 信號DOA 估計誤差偏大的問題，提出了一種新的基于FRFT 的中心頻率估計方法，并在此基礎(chǔ)上對分?jǐn)?shù)階傅里葉域的陣列流行向量進(jìn)行了修正，改善了基于FRFT 的MUSIC 算法的寬帶LFM 信號DOA 估計性能，分析了信號寬度和觀測時間寬度比值對分?jǐn)?shù)階傅里葉域信噪比及估計性能的影響。仿真結(jié)果表明，修正后的算法較傳統(tǒng)的基于FRFT 的MUSIC 算法:1)對脈沖信號中心頻率估計比傳統(tǒng)算法準(zhǔn)確;2)DOA 估計的準(zhǔn)確性較高;3)在脈沖信號寬度特別小時，也能夠準(zhǔn)確估計DOA，而且方差較小。該方法擴(kuò)大了用FRFT 處理寬帶LFM 信號的適用范圍。

References)

［1］ Xiang L，Ye Z，Xu X，et al.Direction of arrival estimation for uniform circular array based on fourth-order cumulants in the presence of unknown mutual coupling［J］.IET Microwaves，Antennas＆ Propagation，2008，2(3):281 -287.

［2］ Liu Z M，Huang Z T，Zhou Y Y.Direction-of-arrival estimation of wideband signals via covariance matrix sparse representation［J］.Signal Processing，2011，59(9):4256 -4270.

［3］ Ozaktas H M，Kutay M A，Zalevsky Z.The fractional Fourier transform with applications in optics and signal processing［M］.New York:John Wiley ＆ Sons，2001:117 -183.

［4］ Cui Y，Liu K，Wang J.Direction of arrival estimation of coherent wideband LFM signals in multipath environment［C］∥2010 IEEE 10th International Conference on Signal Processing (ICSP).Beijing:IEEE，2010:58 -61.

［5］ Chong H，Xiaomin Z.DOA estimation of multi-component LFM in complex environment using ESPRIT based on FRFT［C］∥2011 IEEE International Conference on Signal Processing，Communications and Computing (ICSPCC).Xi’an:IEEE，2011:1 -4.

［6］ Jin X，Zhang T，Bai J，et al.DOA estimation of coherent wideband LFM signals based on fractionalFourier transform and virtual array［C］∥2010 3rd International Congress on Image and Signal Processing (ICISP).Yantai:IEEE，2010:4380 -4384.

［7］ 陶然，周云松.基于分?jǐn)?shù)階傅里葉變換的寬帶LFM 信號波達(dá)方向估計新算法［J］.北京理工大學(xué)學(xué)報，2005，25(10):895-899.TAO Ran，ZHOU Yun-song.A novel method for the direction of arrival estimation of wideband linear frequency modulated sources based on fractional Fourier transform［J］.Transactions of Beijing Institute of Technology，2005，25(10):895 -899.(in Chinese)

［8］ 劉小河，王建英，楊美英.基于分?jǐn)?shù)階傅里葉變換的寬帶LFM相干信號的DOA 估計［J］.數(shù)據(jù)采樣與處理，2008，23(5):547 -550.LIU Xiao-he，WANG Jian-ying，YANG Mei-ying.DOA estimation of coherent wideband LFM sources based on fractional Fourier transform［J］.Journal of Data Acquisition ＆ Processing，2008，23(5):547 -550.(in Chinese)

［9］ 張艷紅，齊林，穆曉敏，等.基于分?jǐn)?shù)階傅立葉變換的WLFM信號DOA 估計［J］.信號處理，2005，21(4):57 -60.ZHANG Yan-hong，QI Lin，MU Xiao-min，et al.DOA estimation of WLFM signal based on fractional Fourier transforms［J］.Signal Processing，2005，21(4):57 -60.(in Chinese)

［10］ 馬艷，羅美玲.基于分?jǐn)?shù)階傅里葉變換水下目標(biāo)距離及速度的聯(lián)合估計［J］.兵工學(xué)報，2011，32(8):1030 -1053.MA Yan，LUO Mei-ling.FRFT-based joint range and radial velocity estimation of underwater target［J］.Acta Armamentarii，2011，32(8):1030 -1053.(in Chinese)

［11］ 吳曉濤.基于分?jǐn)?shù)階傅里葉變換的信道估計算法研究［D］，哈爾濱:哈爾濱工業(yè)大學(xué)，2010.WU Xiao-tao.Research on channel estimation algorithm based on fractional Fourier transform［D］.Harbin:Harbin Institute of Technology，2010.(in Chinese)

［12］ Hsue W L，Pei S C.Rational-ordered discrete fractional Fourier transform［C］∥2012 Proceedings of the 20th European Signal Processing Conference (EUSIPCO).Bucharest:IEEE，2012:2124 -2127.

［13］ Pei S C，Liu C L，Lai Y C.The generalized fractional Fourier transform［C］∥2012 IEEE International Conference on Acoustics，Speech and Signal Processing (ICASSP).Kyoto:IEEE，2012:3705 -3708.

［14］ Ko? A，Ozaktas H M，Candan C，et al.Digital computation of linear canonical transforms［J］.Signal Processing，2008，54(6):2383 -2394.

［15］ Oktem F S，Ozaktas H M.Linear algebraic analysis of fractional Fourier domain interpolation［C］∥17th Signal Processing and Communications Applications Conference (SIU)，Antalya:IEEE，2009:873 -875.

［16］ Ozaktas H M，Arikan O，Kutay M A，et al.Digital computation of the fractional Fourier transform［J］.IEEE Transactions on Signal Processing，1996，44(9):2141 -2150.

［17］ 趙興浩，鄧兵，陶然.分?jǐn)?shù)階傅立葉變換數(shù)值計算中的量綱歸一化［J］.北京理工大學(xué)學(xué)報，2005，25(4):360 -364.ZHAO Xing-hao，DENG Bing，TAO Ran.Dimensional normalization in the digital computation of the fractional Fourier transform［J］.Transactions of Beijing Institute of Technology，2005，25(4):360 -364.(in Chinese)

［18］ 劉峰，徐會法，陶然，等.分?jǐn)?shù)階Fourier 域多分量LFM 信號間的分辨研究，中國科學(xué):信息科學(xué)，2012，42(2):136 -148.LIU Feng，XU Hui-fa，TAO Ran，et al.Research on resolution among multi-component LFM signals in the fractional Fourier domain［J］.Scientia Sinica Informationis，2012，42(2):136 -148.(in Chinese)

［19］ 馬艷，王瑞，陳航.一種FRFT 變換幅度峰值的快速搜索［C］∥2012 魚雷自導(dǎo)引信技術(shù)研討會論文集.北京:中國造船工程學(xué)會水中兵器委員會，2012，9:95 -98.MA Yan，WANG Rui，CHEN hang.A method of transform amplitude quick search in FRFT［C］∥2012 Torpedo Homing Letter Technology Conference Proceedings.Beijing:China Shipbuilding Engineering Society Underwater Weapons Committee，2012，9:95-98.(in Chinese)