生雪莉，穆夢飛，殷敬偉，楊超然，劉婷

(1.哈爾濱工程大學(xué) 水聲技術(shù)重點實驗室，黑龍江 哈爾濱 150001；2.海洋信息獲取與安全工業(yè)和信息化部重點實驗室(哈爾濱工程大學(xué))，黑龍江 哈爾濱 150001；3. 哈爾濱工程大學(xué) 水聲工程學(xué)院，黑龍江 哈爾濱 150001)

隨著科技的發(fā)展，水聲目標識別面臨著越來越大的挑戰(zhàn)，為提高目標識別的準確度，需要綜合運用各種反應(yīng)目標信息的特征量。多普勒效應(yīng)是水下運動目標的重要特征之一，通過對主動聲吶回波進行多普勒頻移估計，可有效提升聲吶系統(tǒng)的檢測與識別能力。水聲環(huán)境的復(fù)雜和回波信號的微弱給目標回波多普勒頻移的估計帶來很大困難。因此，抗干擾能力強的頻移估計方法逐漸成為水聲目標識別領(lǐng)域的研究熱點。

在主動聲吶探測領(lǐng)域，匹配濾波算法是比較常用的檢測回波多普勒頻移的方法，但由于計算量的限制，該方法的估計精度有限，為提高頻移估計精度，對目標回波直接進行高精度的頻率參數(shù)估計并與發(fā)射信號的頻率參數(shù)作比較來得到頻移信息的方法也是多普勒頻移估計的重要手段之一。

信號的頻率估計被廣泛應(yīng)用在雷達、聲吶的動目標檢測與參數(shù)估計問題中，自20世紀50年代以來，各種有效的頻率估計方法不斷被研究和發(fā)展起來。最經(jīng)典的頻率估計方法就是傳統(tǒng)的離散傅里葉變換(discrete Fourier transform, DFT)[2]，但該方法受觀測樣本長度的影響較大，短樣本長度下，頻率分辨力較差。AR模型(auto-regressive，AR)估計是一種高分辨現(xiàn)代譜估計方法[3]，但在使用時需要選擇合適的階數(shù)，階數(shù)選取過小或過大會導(dǎo)致譜峰不明顯或者出現(xiàn)偽峰，且在分析短樣本數(shù)據(jù)時會產(chǎn)生譜偏移和譜分裂的問題。參數(shù)法中的最大似然估計法(maximum-likelihood，ML)理論上可以達到頻率估計的最優(yōu)性能，但計算量過于龐大，不適合實時應(yīng)用[4]。以MUSIC算法和ESPRIT算法為代表的子空間分解算法[5]，根據(jù)信號與噪聲子空間的正交性將信號子空間提取出來，可以實現(xiàn)高信噪比下的高精度頻率估計，但低信噪比下誤差較大，且計算復(fù)雜，計算量大。

線性調(diào)頻信號(LFM)作為一種瞬態(tài)信號，對這類信號的頻率參數(shù)進行估計，傳統(tǒng)的離散傅里葉分析方法已不再適用。一般使用聯(lián)合時頻分析方法進行分析。短時傅里葉變換是最常用的一種時頻分析方法，但分辨力受窗函數(shù)的約束，窗長取較短時，時間分辨力較高但頻率分辨力較低，窗長取較長時，頻率分辨力高但時間分辨力較低[6]。Wigner-Ville 分布(WVD)屬于二次型時頻分析方法[7]，比短時傅里葉變換擁有更加優(yōu)良的時頻分辨性能，但容易受到交叉項的干擾，后續(xù)又提出了一系列的改進算法來抑制交叉項的干擾，取得了良好的效果。近年來分數(shù)階傅里葉變換(fractional Fourier transform, FRFT)受到了人們的廣泛關(guān)注[8]，利用信號在某一分數(shù)階次域的聚焦特性，來得到LFM信號的中心頻率與調(diào)頻率的參數(shù)估計，具有較強的抗干擾能力。

將信號在一組過完備原子庫上進行稀疏分解來實現(xiàn)信號的參數(shù)估計是一類新的方法。用來表示信號的基原子可以根據(jù)信號的內(nèi)在特性靈活選取，分解的結(jié)果可以得到一個關(guān)于信號的稀疏表示，該過程被稱為信號的稀疏分解。目前稀疏分解算法得到了快速發(fā)展，其中最常用的方法是Pati等基于匹配追蹤(matching pursuit,MP)算法[9]提出的正交匹配追蹤(orthogonal matching pursuit,OMP)算法[10]。本文引入了一種基于稀疏分解的頻率估計方法，并使用OMP算法進行信號的稀疏重構(gòu)，可以在低信噪比下對目標回波的頻率參數(shù)進行高精度的估計，進而估計出回波的多普勒頻移。針對精度越高，原子庫越龐大，計算量越大的問題，使用了一種快速算法，在保證精度的同時有效降低了稀疏分解的計算量。通過數(shù)值仿真，驗證了算法的有效性。

1 運動目標回波信號模型

在主動探測中，當聲吶系統(tǒng)與目標之間存在徑向速度時，會產(chǎn)生多普勒效應(yīng)，表現(xiàn)為回波信號在時域上發(fā)生了伸縮，在頻域上發(fā)生了頻率偏移。現(xiàn)對發(fā)生多普勒頻移后的回波信號模型進行推導(dǎo)。

假設(shè)一臺收發(fā)合置的主動聲吶靜止在某處，此時一個目標正以徑向速度v向聲吶運動過來。假設(shè)聲吶的發(fā)射信號為：

(1)

初始時刻，目標與聲吶之間的距離為L，發(fā)射脈沖遇到目標返回的時間記為t1，在這段時間內(nèi)目標與聲吶之間的距離縮短了vt1/2，由圖1(a)中的幾何關(guān)系可知：

(2)

圖1 聲吶與目標存在徑向運動時脈沖前沿和后沿接收示意Fig.1 Schematic diagram of receiving pulse front and back edges when sonar and target have radial motion

進一步得到：

(3)

當脈沖后沿剛好離開聲吶時，聲吶與目標之間的距離縮短了vT，設(shè)脈沖后沿從離開聲吶到回到聲吶經(jīng)過的時間為t2，由圖1(b)的幾何關(guān)系可以得到：

(4)

進一步得到：

(5)

由式(3)、(5)可以看出，脈沖前沿和脈沖后沿的往返時間并不相同，設(shè)二者差值為t′，得到：

(6)

已知發(fā)射脈寬為T，則接收信號脈寬為：

(7)

式中s為多普勒伸縮因子，表示為：

(8)

在水聲環(huán)境下，目標運動速度一般滿足v?c，則式(8)可以進一步等效為：

(9)

式中Δ=2v/c為多普勒頻移因子。接收信號r(t)可以表示為：

(10)

多普勒伸縮因子s=1-Δ表征目標發(fā)生多普勒效應(yīng)后時域發(fā)生伸縮的尺度。而多普勒頻移因子Δ=2v/c表征目標發(fā)生多普勒效應(yīng)后頻域頻移的尺度。

對于窄帶信號多普勒效應(yīng)可視為簡單的載頻偏移，滿足公式：

(11)

式中：fd為頻移；fecho為回波信號的頻率；f0為發(fā)射信號頻率。

對于LFM這種寬帶信號，不同頻率成分多普勒頻移的大小不同，設(shè)起始頻率為fL，截止頻率為fH，調(diào)頻斜率為：

(12)

發(fā)生多普勒效應(yīng)之后，回波信號起始頻率為sfL，截止頻率為sfH，調(diào)頻斜率變?yōu)椋?/p>

(13)

式中：kecho為回波信號的調(diào)頻率；k為發(fā)射信號的調(diào)頻率。

進一步得到回波的多普勒頻移因子的估計表達式為：

(14)

2 基于稀疏分解的頻移估計

2.1 利用稀疏分解估計CW回波頻移

CW回波信號的模型為：

x(t)=Acos(2πf+φ0)+n(t)

(15)

式中：A、f、φ0分別為CW回波信號的幅值、頻率和初相；n(t)為背景噪聲。本文的目的是從背景噪聲中估計出回波信號的頻率和初相，進一步得到回波的多普勒頻移。

傳統(tǒng)的信號分解是將信號分解到一組完備正交基上，但這種方法并不是最優(yōu)的，近年來“冗余字典”不斷被提出并代替?zhèn)鹘y(tǒng)的正交基，通過將信號分解到一組過完備的非正交基上來得到信號的稀疏表示[11]。冗余字典中的基被稱為原子，字典原子可以根據(jù)待分解信號的內(nèi)在特性進行構(gòu)造。關(guān)于信號稀疏分解的方法有很多，本文采用應(yīng)用比較廣泛的正交匹配追蹤(OMP)算法。

首先通過信號模型(15)來構(gòu)造過完備原子庫：

gγi(t)=cos(2πfmt+φj)，n=1,2,…，N

(16)

γi=(m,j)，i=1,2,…，M×J

(17)

式中：N為觀測信號的長度，也是字典的行數(shù)；fm為頻率參數(shù)，按照期望的搜索精度均勻取值；M為頻率搜索個數(shù)，在固定的頻率搜索范圍內(nèi)，M越大，頻率估計精度越高；φj為初相參數(shù)，J為初相搜索個數(shù)，在φj∈[0,2π]內(nèi)，J越大，初相估計精度越高，M×J為字典的列數(shù)，即字典的長度。對這由M×J個原子組成的過完備原子庫進行稀疏分解，可以求得待估計信號的頻率參數(shù)和初相參數(shù)。

設(shè)回波信號采樣序列為y，OMP算法的實現(xiàn)步驟為：

1)初始化殘差余量R0=y，匹配原子集賦空值activ-set=?，匹配原子記錄矩陣賦空值A(chǔ)ug-t=?，迭代次數(shù)計數(shù)器置1，即time=1，最大迭代次數(shù)為k，系數(shù)數(shù)組hat-y賦零值；

2)R0與字典gγ的所有列向量gγi求內(nèi)積，將內(nèi)積最大值對應(yīng)的列向量位置記為pos；

3)記錄字典gγ中pos對應(yīng)的列，將其擴展至Aug-t矩陣，同時在字典gγ中將pos列去除；

6)判斷time>K是否成立，若成立，則回到步驟2)繼續(xù)執(zhí)行，若不成立，執(zhí)行步驟7)；

OMP算法流程圖如圖2所示。

圖2 OMP算法執(zhí)行流程Fig.2 OMP algorithm execution flow chart

2.2 利用稀疏分解估計LFM回波頻移

考慮到稀疏分解的特性，建立的原子要盡可能的逼近LFM信號的結(jié)構(gòu)，根據(jù)LFM的信號形式，建立原子[12]：

(18)

式中：γ=(fu,Kv)為原子參數(shù)組。fu、Kv分別對應(yīng)LFM信號的起始頻率和調(diào)斜率，且fu∈(fmin,fmax)，u=1，2，…，U，U為起始頻率搜索個數(shù)；Kv∈(Kmin,Kmax)，v=1，2，…，V，V為調(diào)頻斜率搜索個數(shù)。fu和Kv根據(jù)范圍和搜索精度均勻取值，可以構(gòu)造出過完備原子庫Gf，原子庫的長度為U×V。由于信號OMP分解的特性，將淹沒在噪聲中的LFM回波信號 進行OMP分解，將會從原子庫中挑選出最逼近真實回波信號的原子，該原子對應(yīng)的參數(shù)組即為回波信號起始頻率和調(diào)頻率的估計值，算法實現(xiàn)步驟與2.1節(jié)中所述一致，由上述描述可以看出該算法會假設(shè)信號均勻存在整個觀測時間段內(nèi)，所以在不知道回波時延的情況下無法準確估計起始頻率，但調(diào)頻率的估計仍然準確，所以可以通過回波信號的調(diào)頻率k與發(fā)射信號的調(diào)頻率kecho的關(guān)系來得到多普勒頻移因子的估計，如式(14)所示。

2.3 快速算法

3 實驗仿真

3.1 CW回波信號的頻移估計仿真

3.1.1 仿真條件

CW發(fā)射信號頻率為7 kHz，脈寬為100 ms，采樣頻率為32 kHz，目標的徑向運動速度為5 kn且靠近聲吶運動，根據(jù)公式計算回波的多普勒頻移為24.007 2 Hz，信噪比為10 dB，取2 048點目標回波數(shù)據(jù)進行處理，一次回波頻率估計結(jié)果如圖3所示。

圖3 稀疏分解與FFT估計結(jié)果對比Fig.3 Comparison of sparse decomposition and FFT estimation results

由圖3可以看出，F(xiàn)FT的回波頻率估計結(jié)果為7 031.25 Hz，頻移估計值為31.25 Hz，估計誤差為7.007 2 Hz；OMP算法的回波頻率估計結(jié)果為7 023.963 9 Hz，頻移估計值為23.963 9 Hz，誤差為-0.043 3 Hz；OMP算法的估計精度要遠高于FFT，且有著極低的旁瓣和極尖銳的譜峰。

3.1.2 性能分析

為了驗證稀疏分解算法的性能，進行蒙特卡洛實驗。作為對比，本文還計算應(yīng)用了工程上常用的一種頻率估計方法，自適應(yīng)Notch濾波器法。該方法由Griffiths最早提出[12]，自適應(yīng)Notch濾波器分為2個部分：窄帶濾波和瞬時頻率估計。本文使用最小均方算法(LMS)作為濾波器的自適應(yīng)學(xué)習(xí)算法。

CW發(fā)射信號頻率為7 kHz，采樣率為32 kHz。噪聲是均值為零，方差為σ2的高斯白噪聲。

1)不同信噪比下的性能分析：信號脈寬為100 ms，目標徑向運動速度設(shè)置為5 kN且目標做靠近聲吶的運動，計算得到回波多普勒頻移為24.007 2 Hz，目標回波采樣點數(shù)為3 200點，取其中的2 048點數(shù)據(jù)進行處理，比較2種算法在不同信噪比下的頻移估計性能，信噪比設(shè)置范圍[-20 dB,20 dB]，步進為2 dB，做1 000次蒙特卡洛實驗，得到的2種算法頻移估計均方根誤差對比結(jié)果如圖4所示。

圖4 不同信噪比下頻移估計均方根誤差對比Fig.4 Comparison of root mean square error of frequency shift estimation under different signal-to-noise ratios

由圖4可知，從估計精度來看，在高信噪比條件下，這2種方法均具有比FFT的譜線分辨力高的多的分辨力和估計精度。從抗噪性能來看，隨著信噪比的降低，自適應(yīng)Notch濾波器法在信噪比為-8 dB以下時出現(xiàn)性能惡化；而稀疏分解法在信噪比為-20 dB時仍能較為準確的估計出回波的頻移，噪聲抑制能力較為穩(wěn)健。

2)不同徑向速度下的算法性能分析：信號脈寬為100 ms，信噪比設(shè)置為10 dB，目標回波采樣點數(shù)為3 200點，取其中的2 048點數(shù)據(jù)進行處理，比較2種算法在多普勒頻移下的頻移估計性能，目標徑向速度設(shè)置范圍[-40 kN,40 kN]，步進為4 kN，做1 000次蒙特卡洛仿真，結(jié)果如圖5所示。

從圖5可以看出，自適應(yīng)Notch濾波器的估計誤差隨多普勒的絕對值增加而增大，比較適合低多普勒下的頻移估計；而稀疏分解法在不同的多普勒下，估計性能都比較穩(wěn)定，且始終高于自適應(yīng)Notch濾波器的估計精度。

圖5 不同徑向速度下頻移估計均方根誤差對比Fig.5 Comparison of root mean square error of frequency shift estimation at different radial velocities

3.2 LFM回波信號的頻移估計仿真

3.2.1 仿真條件

LFM信號發(fā)射頻帶為6～8 kHz，脈寬為250 ms，采樣頻率為32 kHz，目標徑向速度為6 kN且向聲吶靠近，對應(yīng)的回波多普勒頻移因子Δ=0.004 115 5，信噪比為3 dB。結(jié)果如下：

圖6 稀疏分解結(jié)果Fig.6 Sparse decomposition results

3.2.2 性能分析

為了驗證稀疏分解算法的性能，進行蒙特卡洛實驗。這里將文獻[14]中的WVD峰值檢測算法，作為一種對比算法，WVD是LFM信號的最佳估計器，可以對信號的瞬時頻率進行無偏的估計，在無噪聲的情況下，LFM信號的WVD分布表現(xiàn)為時頻平面的一條“脊線”，信號的瞬時頻率值就對應(yīng)著脊線的峰值，因此對信號的瞬時頻率進行估計即是對WVD時頻分布的最大值進行提取。

圖7 重構(gòu)信號與真實信號的瞬時頻率對比Fig.7 Comparison of instantaneous frequency of reconstructed signal and real signal

仿真信號參數(shù)：LFM發(fā)射信號頻帶為6～8 kHz，脈寬為250 ms，采樣率為32 kHz，噪聲是均值為零，方差為σ2的高斯白噪聲。

1)不同信噪比下的性能分析：目標徑向運動速度設(shè)置為6 kN且目標做靠近聲吶的運動，計算得到回波的多普勒頻移因子為Δ=0.004 115 5，目標回波采樣點數(shù)為8 000點，比較各算法在不同信噪比下的頻移因子估計性能，信噪比設(shè)置范圍[-20 dB,20 dB]，步進為5 dB，做100次蒙特卡洛實驗，結(jié)果如圖8所示。

圖8 不同信噪比下的頻移因子估計均方根誤差對比Fig.8 Comparison of root mean square error of frequency shift factor estimation under different signal-to-noise ratios

如圖8所示，在高信噪比下WVD峰值檢測法和稀疏分解法的估計精度較高且基本相當，隨著信噪比的降低，WVD峰值檢測法在信噪比為-5 dB以下時出現(xiàn)性能惡化，而稀疏分解法的噪聲抑制能力較為穩(wěn)健，在信噪比為-10 dB以下才時出現(xiàn)性能惡化。

2)不同徑向速度下的性能分析：信噪比設(shè)置為3 dB，目標回波采樣點數(shù)為8 000點，比較各算法在不同目標徑向速度下的頻移因子估計性能，速度設(shè)置范圍[-20 kN,20 kN]，步進為4 kN，做100次蒙特卡洛實驗，結(jié)果如圖9所示。

圖9 不同徑向速度下頻移因子估計均方根誤差對比Fig.9 Comparison of root mean square error of frequency shift factor estimation at different radial velocities

從圖9可以看到稀疏分解法和WVD峰值檢測法在不同的多普勒下，頻移因子的估計精度都較為穩(wěn)定，且稀疏分解的估計精度要高于WVD峰值檢測法。

4 結(jié)論

1)本文利用稀疏分解方法估計CW回波和LFM回波的多普勒頻移，具有較高的估計精度。同時使用快速算法，在保證估計精度的同時，有效降低了稀疏分解算法的計算量。

2)與常規(guī)方法相比，稀疏分解方法體現(xiàn)出了優(yōu)良的棒性，在極低的信噪比條件下仍能對回波的頻移進行有效估計，抗噪聲性能穩(wěn)健。在不同的目標運動速度下，估計性能穩(wěn)定。

該方法對解決在強干擾環(huán)境下水聲目標探測與識別的問題提供了思路和依據(jù)，具有一定的理論意義與工程應(yīng)用價值。