基于稀疏分解的水下運動目標多普勒頻移估計方法

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

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

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

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

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

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

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

1 運動目標回波信號模型

在主動探測中，當聲吶系統與目標之間存在徑向速度時，會產生多普勒效應，表現為回波信號在時域上發生了伸縮，在頻域上發生了頻率偏移?，F對發生多普勒頻移后的回波信號模型進行推導。

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

(1)

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

(2)

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

進一步得到：

(3)

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

(4)

進一步得到：

(5)

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

(6)

已知發射脈寬為T，則接收信號脈寬為：

(7)

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

(8)

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

(9)

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

(10)

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

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

(11)

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

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

(12)

發生多普勒效應之后，回波信號起始頻率為sfL，截止頻率為sfH，調頻斜率變為：

(13)

式中：kecho為回波信號的調頻率；k為發射信號的調頻率。

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

(14)

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

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

CW回波信號的模型為：

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

(15)

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

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

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

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

(16)

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

(17)

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

設回波信號采樣序列為y，OMP算法的實現步驟為：

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

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

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

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

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

圖2 OMP算法執行流程Fig.2 OMP algorithm execution flow chart

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

考慮到稀疏分解的特性，建立的原子要盡可能的逼近LFM信號的結構，根據LFM的信號形式，建立原子[12]：

(18)

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

2.3 快速算法

3 實驗仿真

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

3.1.1 仿真條件

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

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

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

3.1.2 性能分析

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

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

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

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

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

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

從圖5可以看出，自適應Notch濾波器的估計誤差隨多普勒的絕對值增加而增大，比較適合低多普勒下的頻移估計；而稀疏分解法在不同的多普勒下，估計性能都比較穩定，且始終高于自適應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信號發射頻帶為6～8 kHz，脈寬為250 ms，采樣頻率為32 kHz，目標徑向速度為6 kN且向聲吶靠近，對應的回波多普勒頻移因子Δ=0.004 115 5，信噪比為3 dB。結果如下：

圖6 稀疏分解結果Fig.6 Sparse decomposition results

3.2.2 性能分析

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

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

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

1)不同信噪比下的性能分析：目標徑向運動速度設置為6 kN且目標做靠近聲吶的運動，計算得到回波的多普勒頻移因子為Δ=0.004 115 5，目標回波采樣點數為8 000點，比較各算法在不同信噪比下的頻移因子估計性能，信噪比設置范圍[-20 dB,20 dB]，步進為5 dB，做100次蒙特卡洛實驗，結果如圖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以下時出現性能惡化，而稀疏分解法的噪聲抑制能力較為穩健，在信噪比為-10 dB以下才時出現性能惡化。

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

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

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

4 結論

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

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

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