随着科技的发展，水声目标识别面临着越来越大的挑战，为提高目标识别的准确度，需要综合运用各种反应目标信息的特征量。多普勒效应是水下运动目标的重要特征之一，通过对主动声呐回波进行多普勒频移估计，可有效提升声呐系统的检测与识别能力。水声环境的复杂和回波信号的微弱给目标回波多普勒频移的估计带来很大困难。因此，抗干扰能力强的频移估计方法逐渐成为水声目标识别领域的研究热点。

在主动声呐探测领域，匹配滤波算法是比较常用的检测回波多普勒频移的方法，但由于计算量的限制，该方法的估计精度有限，为提高频移估计精度，对目标回波直接进行高精度的频率参数估计并与发射信号的频率参数作比较来得到频移信息的方法也是多普勒频移估计的重要手段之一。

信号的频率估计被广泛应用在雷达、声呐的动目标检测与参数估计问题中，自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 运动目标回波信号模型

在主动探测中，当声呐系统与目标之间存在径向速度时，会产生多普勒效应，表现为回波信号在时域上发生了伸缩，在频域上发生了频率偏移。现对发生多普勒频移后的回波信号模型进行推导。

假设一台收发合置的主动声呐静止在某处，此时一个目标正以径向速度v向声呐运动过来。假设声呐的发射信号为：

$ {g_T}(t) = \left\{ {\begin{array}{*{20}{c}} {g(t),}&{t \in [0,T]}\\ {0,}&{{\rm{其他}}} \end{array}} \right. $

(1)

初始时刻，目标与声呐之间的距离为L，发射脉冲遇到目标返回的时间记为t₁，在这段时间内目标与声呐之间的距离缩短了vt₁/2，由图 1(a)中的几何关系可知：

$ L = \frac{{v{t_1}}}{2} + \frac{{c{t_1}}}{2} $

(2)

进一步得到：

$ {t_1} = \frac{{2L/c}}{{1 + v/c}} $

(3)

当脉冲后沿刚好离开声呐时，声呐与目标之间的距离缩短了vT，设脉冲后沿从离开声呐到回到声呐经过的时间为t₂，由图 1(b)的几何关系可以得到：

$ L = vT + \frac{{c{t_2}}}{2} + \frac{{v{t_2}}}{2} $

(4)

进一步得到：

$ {t_2} = \frac{{\frac{{2L}}{c} - \frac{{2vT}}{c}}}{{1 + \frac{v}{c}}} $

(5)

由式(3)、(5)可以看出，脉冲前沿和脉冲后沿的往返时间并不相同，设二者差值为t′，得到：

$ {t^\prime } = {t_1} - {t_2} = \frac{{2vT/c}}{{1 + v/c}} $

(6)

已知发射脉宽为T，则接收信号脉宽为：

$ {T_1} = T - {t^\prime } = \frac{{1 - \frac{v}{c}}}{{1 + \frac{v}{c}}}T = sT $

(7)

式中s为多普勒伸缩因子，表示为：

$ s = \frac{{1 - \frac{v}{c}}}{{1 + \frac{v}{c}}} = \frac{{c - v}}{{c + v}} $

(8)

在水声环境下，目标运动速度一般满足v≪c，则式(8)可以进一步等效为：

$ s \approx 1 - \frac{{{v^2}}}{c} \approx 1 - \frac{{2v}}{c} = 1 - \Delta $

(9)

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

$ r(t) = \sqrt s * g(s * (t - \tau )) $

(10)

多普勒伸缩因子s=1－Δ表征目标发生多普勒效应后时域发生伸缩的尺度。而多普勒频移因子Δ=2v/c表征目标发生多普勒效应后频域频移的尺度。

对于窄带信号多普勒效应可视为简单的载频偏移，满足公式：

$ {f_{\rm{d}}} = \frac{{2v}}{c}{f_0} = {f_{{\rm{echo}}}} - {f_0} $

(11)

式中：f_d为频移；f_echo为回波信号的频率；f₀为发射信号频率。

对于LFM这种宽带信号，不同频率成分多普勒频移的大小不同，设起始频率为f_L，截止频率为f_H，调频斜率为：

$ k = \frac{{{f_H} - {f_L}}}{T} $

(12)

发生多普勒效应之后，回波信号起始频率为sf_L，截止频率为sf_H，调频斜率变为：

$ {k_{{\rm{echo}}}} = {s^2}\left( {\frac{{{f_H} - {f_L}}}{T}} \right) = {s^2}k $

(13)

式中：k_echo为回波信号的调频率；k为发射信号的调频率。

进一步得到回波的多普勒频移因子的估计表达式为：

$ \hat \varDelta = s - 1 = \sqrt {\frac{{{k_{{\rm{echo}}}}}}{k}} - 1 $

(14)

2 基于稀疏分解的频移估计 2.1 利用稀疏分解估计CW回波频移

CW回波信号的模型为：

$ x(t) = A\cos (2{\rm{ \mathsf{ π} }}f + {\varphi _0}) + n(t) $

(15)

式中：A、f、φ₀分别为CW回波信号的幅值、频率和初相；n(t)为背景噪声。本文的目的是从背景噪声中估计出回波信号的频率和初相，进一步得到回波的多普勒频移。

传统的信号分解是将信号分解到一组完备正交基上，但这种方法并不是最优的，近年来“冗余字典”不断被提出并代替传统的正交基，通过将信号分解到一组过完备的非正交基上来得到信号的稀疏表示^[11]。冗余字典中的基被称为原子，字典原子可以根据待分解信号的内在特性进行构造。关于信号稀疏分解的方法有很多，本文采用应用比较广泛的正交匹配追踪(OMP)算法。

首先通过信号模型(15)来构造过完备原子库：

$ {{g_{{\gamma _i}}}(t) = \cos (2{\rm{ \mathsf{ π} }}{f_m}t + {\varphi _j}),n = 1,2, \cdots ,N} $

(16)

$ {{\gamma _i} = (m,j),i = 1,2, \cdots ,M \times J} $

(17)

式中：N为观测信号的长度，也是字典的行数；f_m为频率参数，按照期望的搜索精度均匀取值；M为频率搜索个数，在固定的频率搜索范围内，M越大，频率估计精度越高；φ_j为初相参数，J为初相搜索个数，在φ_j∈[0, 2π]内，J越大，初相估计精度越高，M×J为字典的列数，即字典的长度。对这由M×J个原子组成的过完备原子库进行稀疏分解，可以求得待估计信号的频率参数和初相参数。

设回波信号采样序列为y，OMP算法的实现步骤为：

1) 初始化残差余量R₀=y，匹配原子集赋空值activ－set=Ø，匹配原子记录矩阵赋空值Aug－t=Ø，迭代次数计数器置1，即time=1，最大迭代次数为k，系数数组hat－y赋零值；

2) R₀与字典g_γ的所有列向量g_γi求内积，将内积最大值对应的列向量位置记为pos；

3) 记录字典g_γ中pos对应的列，将其扩展至Aug－t矩阵，同时在字典g_γ中将pos列去除；

4) 根据最小二乘法估计得到$\hat{\alpha}_{\text {time }}$=(Aug－t^T*Aug－t)^－1*Aug－t^T*y，且time=time+1；

5) 更新残差R=y－Aug－t*$\hat{\alpha}_{\text {time }}$，记录最大投影系数的位置，稀疏表示向量表示为hat－y(pos)=$\hat{\alpha}_{\text {time }}$；

6) 判断time>K是否成立，若成立，则回到步骤2)继续执行，若不成立，执行步骤7)；

7) 记录hat－y中最大值对应的位置p，原子库中p位置处原子对应的参数集($\hat{f}, \hat{\varphi}$)即为所求信号的参数估计值(频率和初相)。

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

2.2 利用稀疏分解估计LFM回波频移

考虑到稀疏分解的特性，建立的原子要尽可能的逼近LFM信号的结构，根据LFM的信号形式，建立原子^[12]：

$ {g_\gamma } = \exp \left\{ {{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{f_u}t + \frac{{{K_v}}}{2}{t^2}} \right)} \right\} $

(18)

式中：γ=(f_u, K_v)为原子参数组。f_u、K_v分别对应LFM信号的起始频率和调斜率，且f_u∈(f_min, f_max)，u=1，2，…，U，U为起始频率搜索个数；K_v∈(K_min, K_max)，v=1，2，…，V，V为调频斜率搜索个数。f_u和K_v根据范围和搜索精度均匀取值，可以构造出过完备原子库G_f，原子库的长度为U×V。由于信号OMP分解的特性，将淹没在噪声中的LFM回波信号进行OMP分解，将会从原子库中挑选出最逼近真实回波信号的原子，该原子对应的参数组即为回波信号起始频率和调频率的估计值，算法实现步骤与2.1节中所述一致，由上述描述可以看出该算法会假设信号均匀存在整个观测时间段内，所以在不知道回波时延的情况下无法准确估计起始频率，但调频率的估计仍然准确，所以可以通过回波信号的调频率k与发射信号的调频率k_echo的关系来得到多普勒频移因子的估计，如式(14)所示。

2.3 快速算法

由于对信号进行稀疏分解时，需要在冗余的原子库中逐个搜索与观测信号最匹配的原子，这就导致了OMP算法在使用时产生巨大的运算量，但是同时又因为OMP不同于其他信号分解方法，其原子库可以根据回波信号的先验知识及分解过程自适应的构造。使得OMP算法可以更灵活更快速地使用。由于水下运动目标速度一般不超过30 m/s，对应的多普勒频移因子的取值范围为Δ∈[－0.04, 0.04]，以CW信号为例，若发射信号的频率为f_CW，回波信号的频率为$\hat{f}$，则回波信号的频率范围为$\hat{f}$∈[0.96f_CW, 1.04f_CW]，本文所用的快速算法是，先在此频率范围内以较大的频率间隔为步长进行搜索，得到回波频率的初估计值，然后在初估计值附近做更加精细的频率搜索，精细搜索的步长依据期望得到的精度来自设置，这个过程可以多次进行，直到精度满足要求。以回波频率变化范围为100 Hz为例，期望得到的频率估计精度为0.01 Hz，若直接在该范围内进行搜索，则需要搜索10 000次，若先以1 Hz进行粗搜索再以0.01 Hz进行精细搜索，搜索次数只需要300次，与直接在大范围频率内进行精细频率搜索相比，该方法的字典长度大大降低，对应的运算量也大大减小。对初相的搜索也可使用类似的快速算法。

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的回波频率估计结果为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。噪声是均值为零，方差为σ²的高斯白噪声。

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可知，从估计精度来看，在高信噪比条件下，这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滤波器的估计精度。

3.2 LFM回波信号的频移估计仿真 3.2.1 仿真条件

LFM信号发射频带为6~8 kHz，脉宽为250 ms，采样频率为32 kHz，目标径向速度为6 kN且向声呐靠近，对应的回波多普勒频移因子Δ=0.004 115 5，信噪比为3 dB。结果如下：

图 6是稀疏分解的结果，根据被挑选出的原子所在的位置求出对应的参数组估计结果，回波起始频率估计值为6 024.720 Hz，调频率估计值为8 065.640 Hz/s。图 7为回波的瞬时频率真实值与估计值的对比，可以看到二者几乎重合，由于在不知道回波时延的情况下无法准确估计起始频率，但调频率的估计仍然准确，所以我们通过回波信号的调频率k_echo与发射信号的调频率k的关系来得到多普勒频移因子的估计，估计值为$\mathit{\hat \Delta }$=0.004 054 2，与理论值基本相符。

3.2.2 性能分析

为了验证稀疏分解算法的性能，进行蒙特卡洛实验。这里将文献[14]中的WVD峰值检测算法，作为一种对比算法，WVD是LFM信号的最佳估计器，可以对信号的瞬时频率进行无偏的估计，在无噪声的情况下，LFM信号的WVD分布表现为时频平面的一条“脊线”，信号的瞬时频率值就对应着脊线的峰值，因此对信号的瞬时频率进行估计即是对WVD时频分布的最大值进行提取。

仿真信号参数：LFM发射信号频带为6~8 kHz，脉宽为250 ms，采样率为32 kHz，噪声是均值为零，方差为σ²的高斯白噪声。

1) 不同信噪比下的性能分析：目标径向运动速度设置为6 kN且目标做靠近声呐的运动，计算得到回波的多普勒频移因子为Δ=0.004 115 5，目标回波采样点数为8 000点，比较各算法在不同信噪比下的频移因子估计性能，信噪比设置范围[-20 dB, 20 dB]，步进为5 dB，做100次蒙特卡洛实验，结果如图 8所示。

如图 8所示，在高信噪比下WVD峰值检测法和稀疏分解法的估计精度较高且基本相当，随着信噪比的降低，WVD峰值检测法在信噪比为-5 dB以下时出现性能恶化，而稀疏分解法的噪声抑制能力较为稳健，在信噪比为-10 dB以下才时出现性能恶化。

2) 不同径向速度下的性能分析：信噪比设置为3 dB，目标回波采样点数为8 000点，比较各算法在不同目标径向速度下的频移因子估计性能，速度设置范围[－20 kN, 20 kN]，步进为4 kN，做100次蒙特卡洛实验，结果如图 9所示。

从图 9可以看到稀疏分解法和WVD峰值检测法在不同的多普勒下，频移因子的估计精度都较为稳定，且稀疏分解的估计精度要高于WVD峰值检测法。

4 结论

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

2) 与常规方法相比，稀疏分解方法体现出了优良的棒性，在极低的信噪比条件下仍能对回波的频移进行有效估计，抗噪声性能稳健。在不同的目标运动速度下，估计性能稳定。

该方法对解决在强干扰环境下水声目标探测与识别的问题提供了思路和依据，具有一定的理论意义与工程应用价值。