1 引言

目标或目标结构部件除过主体平动之外的旋转、振动等微小运动称为微动^[1]。微动反映了目标的精细运动信息与结构信息，因此微动特征通常被认为是区分不同目标独一无二的特征。通过从目标回波中提取精细的微动特征，可为目标的分类识别提供可靠依据^[2]。近年来，基于目标微动特征提取的空间、空中、地面目标、海面舰船等的分辨与识别技术得到了广泛的研究^[3,4]。旋转、振动、进动等微动引起回波的周期性附加多普勒调制，针对这一周期性特征，基于联合时频分布(Joint Time-Frequency Distribution, JTFD)的微动特征提取技术得到了广泛的应用^[5–8]。此类方法首先计算目标回波的时频分布，并通过其周期性参数来获得目标的微动特征。由于回波频率时变，因此JTFD的窗长设定应远小于微动周期，同时为了获得精确的微动频率，窗长又需尽可能地涵盖足够多的采样点，可见微动特征提取精度易受窗长影响；此外，JTFD方法对噪声较为敏感，在低信噪比条件下算法鲁棒性难以保证。归一化幅度差函数法与自相关法不受窗长影响，以其简单易行、运算速度快得到了广泛的应用^[9,10]。然而，这类方法存在多周期混淆的现象，且微动频率提取精度受到采样频率的制约，不适用于精度要求较高的场合。此外，由于常见周期性运动引起的回波附加多普勒调制，多表现为正弦调频或相位项为多个正弦信号叠加的形式，因此，基于微动回波为正弦调频形式的特征提取技术近年来也得到了发展。Hough变换是微多普勒特征提取较为有效的方法之一，它通过谱图域提取曲线参数信息，再通过换算至正弦域得到目标的微多普勒特征^[11,12]。

以上方法主要针对孤立目标回波进行微动特征提取，然而在实际场景中常常面临由多个目标组成的群目标，此时以上方法将不再适用。对于群目标回波，正弦调频傅里叶变换通过将回波分解在正交三角函数基上计算出回波的微动频谱^[13]，从而能够完成两个微动目标的回波频率提取，然而对于3个或3个以上的目标其频谱将出现干扰项，难以获得准确的微动频率^[14]。近年来，由于Bessel函数良好的正交性及其与信号频率成分的线性关系，使得以Bessel函数为基函数的傅里叶-贝塞尔变换(Fourier-Bessel Transform, FBT)在信号分离中得到了良好的应用^[15,16]。然而，由于与Bessel函数基对应的信号频率成分不可细分，因此FBT很少应用于特征提取中。为克服这一问题，k分辨参数被引入Bessel函数基中，使得一定时间积累下的信号频率分辨率得到了提升^[17]。

针对以上问题，本文在传统FBT及文献[17]的基础上，综合考虑调制指数、信号幅度等参数对分解结果的影响，提出了一种正弦调频Fourier-Bessel变换(Sinusoidal Frequency Modulation Fourier-Bessel Transform, SFMFBT)，分析了SFMFBT的相关性质，并将其引入群目标微动特征提取与回波分离重构中，具体工作内容如下：(1)将k分辨参数引入Bessel函数基，基于该函数基提出SFMFBT变换对，并对相关性质进行了推导；(2)将SFMFBT引入微动特征提取中，分析了误差的产生原因，并讨论了离散信号处理中的相关问题。与现有部分算法相比，基于SFMFBT的微动频率提取方法不需时域加窗，因此能够获得更高的积累增益与提取精度；此外，该方法同样适用于时频域交叠的微动群目标回波。

2 正弦调频Fourier-Bessel变换 2.1 变换定义

SFMFBT将信号的相位历程分解至各项 $\alpha $ 阶Bessel函数基上，信号 $s\left( t \right)$ 的SFMFBT定义为：

其中， ${J_\alpha }\left( {\nu t} \right)$ 为 $\alpha $ 阶Bessel函数， $\nu $ 为SFMFB域变量。相应地，正弦调频傅里叶-贝塞尔反变换(Inverse SFMFBT, ISFMFBT)定义为：

对于有限时域信号， $s\left( t \right),t \in \left( {0,T\;} \right)$ 可被分解为有限项加权的Bessel函数：

其中，k为SFMFB核函数分辨率，C_m为第m项SFMFB系数， ${\mu _m}$ 为 $\alpha $ 阶Bessel函数的第m项升序正根。将以k为分辨率的SFMFB级数称为k-SFMFB级数(简称SFMFB级数)。式(3)中第m项SFMFB级数C_m计算为：

2.2 性质分析

下面以正弦调频信号 $s\left( t \right)$ 为例，进行SFMFBT以分析变换结果与信号频率之间的关系。设有限时域信号

其中， $\sigma $ 为信号散射系数， $\omega $ 为信号调制频率，a和b分别为正弦项与余弦项的调制指数。将信号 $s\left( t \right)$ 代入式(4)中得

其中，C₀为虚数，C₁, C₂分别为：

考虑Bessel函数的不定积分性质^[18]

及以下积分^[19]

将k分辨参数引入式(11)，并与式(9)、式(10)代入式(6)，可得C₁和C₂分别为：

联立式(6)、式(12)和式(13)可得第m项SFMFB级数为

其中，C₀为虚数项， $\tan {\varphi _m} = \left( {\mu _m^2 - {k^2}{T\;^2}{\omega ^2}} \right)\big/kT\;\omega $ 且 $\tan \phi = b/a$ 。由于T 和 $\omega $ 均为常数，因此SFMFB级数实值部分的幅值随项数m而变化。当 ${\mu _m} \to \omega kT\;$ 时， $\left| {{\mathop{\rm Re}\nolimits} \left( {{C_m}} \right)} \right|$ 取得最大值，因此可通过级数实部最大值项数进行频率的估计，即

其中， ${m_{\max }} = \mathop {\arg \max }_m \left\{ {\left| {{\mathop{\rm Re}\nolimits} \left( {{C_m}} \right)} \right|} \right\}$ 。

由Bessel函数根的性质^[20]，当 $m \to \infty $ 时，相邻两项Bessel函数正根之差 $\Delta {\mu _m}$ 趋近于 ${{π}}$ ：

因此，相邻两项SFMFB级数对应的频率分辨率为：

令信号实际频率f与提取频率 ${f\!_{{m_{\max }}}}$ 之差的绝对值为绝对误差 ${f\!_e} = \left| {{f\!_{{m_{\max }}}} - f} \right|$ ，则任一信号实际频率的最大绝对误差满足 $\max {f\!_e} < 1/2kT\;$ ，则当信号足够长时，k分辨越大分辨率 $\Delta {f\!_k}$ 越为精细，频率提取结果越精确。

3 基于SFMFBT的微动特征提取 3.1 信号模型

设空间存在某自旋群目标，雷达与群目标质心 $o'$ 的初始距离为R₀，该群目标包含P个子目标，第p个子目标质心为 $ o\!_p$ ，其上含有 ${q_p}\left( {p = 1,2, ·\!·\!· ,P} \right)$ 个散射点。以雷达O为原点建立全局坐标系 $\left( {X,Y,Z} \right)$ ，对于每个子目标分别以目标质心o_p为原点建立本地坐标系 $\left( {x,y,z} \right)$ ，坐标 $\left( {X,Y,Z} \right)$ 的3个坐标轴分别与 $\left( {x,y,z} \right)$ 的坐标轴平行。目标p以角速度 ${{{ω}}_p} = {\Big( {{\omega _{px}},{\omega _{py}},{\omega _{pz}}} \Big)^{\rm{T}}}$ 进行自旋，其自旋轴与雷达视线(Line-Of-Sight, LOS)方向单位向量n的夹角为 ${\varepsilon _p}$ 。雷达与群目标的坐标几何关系示意图如图1所示。

雷达发射载频为f_c的单频连续波，则在t时刻散射点q_p的回波为：

其中， ${\sigma\! _{{q_p}}}$ 为散射点q_p的散射系数， ${r\!_{{q_p}}}\left( t \right)$ 为t时刻点q_p与雷达之间的距离

其中， ${\varOmega _p} = {\left\| {{{{ω}}_p}} \right\|_2}$ 为目标p的旋转角速度， ${\theta _{{q_p}}}$ 为q_p点初相，对 ${s_{{q_p}}}\left( t \right)$ 作基带变换后得到

其中， ${\varPhi _{{q_p}}}\left( t \right)$ 为相位项。对相位项关于时间求导，可得目标回波多普勒频率为：

由式(21)可知，自旋目标的回波频率为正弦形式，自旋群目标回波为多分量正弦调频信号的叠加。

3.2 微动频率提取算法

在实际应用中，雷达回波经数字采样后为离散形式，已知离散回波的相位测量范围为 $\left( { - {{π}},{{π}}} \right)$ ，当回波总相移超过 $2{{π}}$ 时产生相位模糊，因此在进行微动频率提取之前应首先对回波进行相位解模糊。当回波相邻两采样点的最大相移小于 $2{{π}}$ 时可通过相邻采样点相位之差进行解模糊处理^[17]。此外，如式(3)所述，离散有限时域信号可展开为有限项SFMFB级数，称SFMFB级数最大项数为最大有效项数M，则最大有效项数M可通过以下计算得出。设回波信号采样频率为 $ f\!_s $ ，信号序列长为N，f_M为第M项级数对应频率。根据Nyquist采样定律，可表征的信号最大频率为 ${f\!_M} = {f\!_s}/2$ ，则 $f\!_M $ 可由首项系数对应频率f₁与 $\Delta {f\!_k}$ 表示为：

则有

由式(15)可知，在对信号进行SFMFBT的计算过程中，微动频率可根据信号的SFMFB系数实部绝对值的最大项得到。接下来将对SFMFBT算法的频率误差进行分析，并在此基础上给出一种误差减小的方法进行微动频率修正。

如我们所知，若对信号进行傅里叶变换(Fourier Transform, FT)，其FT频谱谱线与信号实际频率的差值不超过相邻谱线频率差值(即频率分辨率)的二分之一，即最大误差与FT基函数的频率分辨率相同，当且仅当信号实际频率落在相邻两函数基等分处时取得最大误差。然而，这一情况仅适用于单频率成分投影在单项级数或单根谱线上。与FT频谱不同，在SFMFB级数中单频率成分投影在若干项SFMFB系数上，即第m项SFMFB系数对应的理论频率f_m在第m项系数上的投影值虽为最大，但其在与第m项邻近的若干项系数上的投影不为零。对信号 $s\left( t \right)$ 进行SFMFBT，若第m项SFMFB系数的实部绝对值与第m+1项SFMFB系数的实部绝对值相等，即 $\left| {{\mathop{\rm Re}\nolimits} \left( {{C_m}} \right)} \right| = \left| {{\mathop{\rm Re}\nolimits} \left( {{C_{m + 1}}} \right)} \right|$ ，则有(推导过程见附录)

由基本不等式性质，则信号 $s\left( t \right)$ 实际频率满足

其中， ${\mu _m} \ne {\mu _{m + 1}}$ 。因此，当 $\left| {{\mathop{\rm Re}\nolimits} \left( {{C_m}} \right)} \right| = \left| {{\mathop{\rm Re}\nolimits} \left( {{C_{m + 1}}} \right)} \right|$ 时，信号实际频率大于第m项与第m+1项SFMFB系数对应理论频率的平均值，其最大绝对误差也大于频率分辨率的二分之一：

下面给出一种将最大绝对误差减小至 $\Delta {f_k}/2$ 的方法。首先计算中心误差 ${f_{{\rm{cen}}}}$ ，并将初始提取频率 ${f_{{m_{\max }}}}$ 与 ${f_{{\rm{cen}}}}$ 之差作为修正的频率。具体计算方法为：(1)根据 ${C_{{m_{\max }}}}$ 计算频率 ${f^{\left( 1 \right)}} = {f_{{m_{\max }}}} - \Delta {f_k}/2$ 和 ${f^{\left( L \right)}} = {f_{{m_{\max }}}} + \Delta {f_k}/2$ ，并将频率离散化为L项： ${f^{\left( l \right)}} \in \left\{ {{f^{\left( 1 \right)}}, ·\!·\!· ,{f^{\left( L \right)}}} \right\}$ ; (2)依次以 ${f^{\left( l \right)}}$ 为调制频率代入 $s\left( t \right)$ ，计算得到最大项系数对应频率 $f_{{m_{\max }}}^{\left( l \right)} \in \left\{ {{f_{{m_{\max }} - 1}},{f_{{m_{\max }}}},{f_{{m_{\max }} + 1}}} \right\}$ ; (3) ${f_{{\rm{cen}}}}$ 即为各项频率 ${f^{\left( l \right)}}$ 与提取频率 $f_{{m_{\max }}}^{\left( l \right)}$ 的最大误差与最小误差的平均值^[16]。

在对群目标回波进行微动频率提取与回波分离重构的过程中，以能量大小作为判别SFMFB系数中是否含有目标回波的依据。由于SFMFB系数与信号频率成分一一对应，则信号能量与SFMFB系数的平方成正比。令 ${\varsigma _{{m_{\max }}}} = {\Large\sum} {{{\left( {{C_{{{m'}_{\max }}}}} \right)}^2}} \Big/ $ $ {\Large\sum}_{m = 1}^M {{{\left( {{C_m}} \right)}^2}}$ ，其中 ${C_{{{m'}_{\max }}}} \!=\! \Big[ {{C_{{m_{\max }}}} \!-\! w,·\!·\!·,{C_{{m_{\max }}}} \!+\! w} \Big]$ 为以 ${C_{{m_{\max }}}}$ 为邻域的若干项系数，存在阈值 $\eta $ ，若 ${\varsigma _{{m_{\max }}}} > \eta $ ，则判定SFMFB域上仍存在目标回波，若 ${\varsigma _{{m_{\max }}}} \le \eta $ ，则不再包含目标。

综上所述，利用SFMFBT算法进行微动群目标频率提取及信号分离重构的步骤如下：

Step 1  初始化：目标个数i=0，幅值最大项系数对应项数 ${m_{\max ,i}}{\rm{ = }}1$ , SFMFBT分辨率 $k \ge 1$ ，信号能量阈值 $\eta > 0$ ；

Step 2  计算最大有效项数 $M = kN - $ $ \displaystyle\frac{{{\mu _1}}}{{{π}}} + 1$ ，频率分辨率 $\;\Delta {f_k} \approx {f_s}/2kN$ ；

Step 3  计算 ${C_m} = {\mathop{\rm SFMFB}\nolimits} \left[ {s\left( n \right)} \right]$ ，寻找 ${m_{\max ,i}} = $ $ \mathop {\arg \max }_m \!\left\{ {\left| {{\mathop{\rm Re}\nolimits} \left( {{C_m}} \right)} \right|} \right\}$ ，计算对应频率 $ {f_{{m_{\max ,i}}}} \!\!=\!\! {\mu _{{m_{\max ,i}}}}\;{f_s}\big/ $ $2{{π}}kN$ ；

Step 4  离散化 ${f^{\left( l \right)}} \in \left[ {{f_{{m_{\max ,i}}{\rm{ - }}1}},{f_{{m_{\max ,i}} + 1}}} \right]$ ，计算中心误差频率 ${f_{{\rm{cen}}}}$ 、修正频率 ${\hat f_i} = {f_{{m_{\max ,i}}}} - {f_{{\rm{cen}}}}$ ；

Step 5  计算 ${\varsigma _{{m_{\max ,i}}}} \!=\!\! {\Large\sum} \!\!{{{\left(\!\! {{C_{{{m'}\!_{\max ,i}}}}} \right)}^2}} \!\! \bigg/\!\!{\Large\sum}_{m = 1}^M {{{\left(\!\! {{C_m}} \right)}^2}} $ ；若 ${\varsigma _{{m_{\max ,i}}}} > \eta $ ，令i=i+1，转至Step 6；若 ${\varsigma _{{m_{\max ,i}}}} \le \eta $ ，停止迭代；

Step 6  令 ${C_m}{\rm{ = }}{C_m}\backslash {C_{{{m'}\!_{\max ,i}}}}$ ；计算 $ {s_i}\left( n \right)= $ $ {\mathop{\rm ISFMFB}\nolimits} \!\left[\! {{C_{{{m'}_{\max ,i}}}}} \right]$ , ${s_i}\left( n \right)$ 即为第i目标回波信号，重复Step 4至Step 6。

4 仿真实验 4.1 算法有效性验证

设空间存在3个自旋锥形目标，每个目标各含有1个锥顶散射点和2个锥底边缘散射点，其中锥顶散射点位于自旋轴上对回波不存在微多普勒调制。目标1以3.0 Hz进行自旋，自旋轴与LOS方向夹角为 ${{π}}/6$ ，底面半径为0.5 m；目标2以6.5 Hz进行自旋，自旋轴与LOS方向夹角为 ${{π}}/9$ ，底面半径为0.3 m；目标3以10.2 Hz进行自旋，自旋轴与LOS方向夹角为 ${{π}}/12$ ，底面半径为0.4 m。雷达照射时间为2 s，脉冲重复频率(Pulse Repetition Frequency, PRF)为256 Hz。分别在k=1, k=4和k=10的条件下，回波信号的SFMFBT结果如图2所示。

从图2可以看出，信号经SFMFBT后表现为SFMFB级数，在不同k分辨率下SFMFB级数均分别在3处取得峰值。按照级数峰值项数计算得到3个目标的初始频率、修正频率和绝对误差如表1所示。在k=10条件下，目标回波及采用ISFMFBT得到的各目标分离重构回波的时频分析如图3所示。

从图3可以看出，从原始目标回波的时频图像难以直接读取群目标的有效频率信息，而各子目标重构信号的时频分布则较好地反映了3个目标各自的频率信息。

4.2 性能分析 4.2.1 算法鲁棒性分析

本节在不同信噪比(Signal-to-Noise Ratio, SNR)条件下对算法的鲁棒性进行验证，并与联合时频分布结合扩展Hough变换(TF-Hough)的微动参数提取方法^[11]进行对比分析。TF-Hough方法在对回波进行时频分析的基础上，通过时频分析结果与扩展Hough空间的映射关系，完成对正弦曲线进行检测及参数提取。以4.1节旋转目标1回波信号为例，分别在SNR=–30 dB至SNR=20 dB的高斯白噪声条件下进行300次Monte Carol实验，频率提取正确概率p_e随SNR及k分辨率的变化曲线如图4所示。

从图4可以看出，取k=1且当SNR>0 dB时，所提算法的频率提取正确概率可达90%以上。即使当信号序列足够长时，频率提取精确程度与k分辨率成正比，但随着k的增大，级数的划分也越为精细，算法也越易受噪声干扰。可见，虽然k分辨率越小频率提取精度越低，但算法在高斯白噪声下鲁棒性越强。对于TF-Hough方法，由于时频分析结果对噪声较为敏感且精度易受窗长影响，当SNR>15 dB时，正弦检测及频率提取的正确概率可达90%。由此可见，基于SFMFBT的方法作为一种参数化微动频率提取方法，与TF-Hough方法相比具有更高的检测精度，在低SNR条件下算法的鲁棒性更强。

下面对阈值 $\eta $ 进行讨论。设噪声分解得到的系数幅值在全部级数项上均匀分布，则有

取w=2k，以4.1节信号模型为例：①在接收信号全部为高斯白噪声(不包含目标回波)，有 ${\varsigma _{{m_{\max }}}} = 0.0081$ ，则阈值下限可设为 ${\eta _L} = 0.0081$ ；②在无噪声条件下，仿真得到三分量信号的比值分别为 ${\varsigma _{{m_{\max 1}}}} \!=\!\! 0.6824$ , ${\varsigma _{{m_{\max 2}}}} \!=\! 0.6596$ 和 ${\varsigma _{{m_{\max 3}}}} \!=\! 0.5895$ ；③在SNR=0 dB条件下，仿真得到三分量信号的比值分别为 ${\varsigma _{{m_{\max 1}}}} = 0.1680$ , ${\varsigma _{{m_{\max 2}}}} = 0.1530$ 和 ${\varsigma _{{m_{\max 3}}}} =$ 0.1179。又由图4的仿真结果可知，在SNR>0 dB且k=1时检测正确概率p_e>90%，则阈值上限可设为 ${\eta _{_U}} = 0.1$ 。此外，门限 $\eta $ 还受到子目标个数、各子目标的强散射点个数、各散射点散射系数、散射点旋转半径、旋转半径所在直线与LOS所呈角度等一些因素的影响，因此阈值可取区间约为 $\eta \in \left[ {0.01, \; 0.1} \right]$ 。

在实际应用中，对于信号 $s\left( t \right)$ ，所取k值越大频率分辨率越小，微动频率估计精度也得到一定程度的提升；另一方面，k值越小算法在相同SNR条件下的鲁棒性越强。因此k的取值应由所需微动频率估计精度、接收信号序列长及SNR共同决定。对于所需精度高、噪声干扰较小的环境，k取值应相对较大；对于噪声干扰较大的环境，k取值应相对较小。

4.2.2 旋转半径对算法精度的影响

同样以自旋目标为例，讨论目标旋转半径对微动频率提取精度的影响。以归一化均方误差(Normalized Root-Mean-Square Error, NRMSE)进行算法提取频率 $\hat f$ 与信号实际频率f 接近程度的度量，即 ${\rm{NRMSE}} \!=\!\! \left| {\hat f - f} \right| \! \Big/ f$ 。由式(21)可知，SFMFB系数的最大值与信号正弦调频形式的调制指数有关。同样以4.1节旋转目标1回波信号为例，在不同旋转半径及SNR=–4 dB至SNR=8 dB下分别进行300次Monte Carol实验，算法提取的旋转频率与实际频率的NRMSE曲线如图5所示。

从图5可以看出，在相同旋转半径下，SNR越大提取频率与实际频率越为接近，对应的NRMSE越小；在相同SNR条件下，旋转半径越大则最大SFMFB级数幅值越为显著，提取频率与实际频率的NRMSE越小，这与式(14)结论一致。

4.2.3 运算量分析

对所提算法进行运算量分析。由式(24)可知，信号的SFMFB级数最大有效项数M与k分辨和信号序列长N的乘积成正比，则计算序列长为N的信号的k分辨SFMFB级数，其算法复杂度为 $O(kN)$ 。下面分别在 $k = {2^0},{2^1}, ·\!·\!· ,{2^5}$ 的条件下，对序列长为N=20, 40, 80, 160的信号进行了500次仿真实验，算法运行时间如图6所示。

由图6可知，SFMFB级数的运算量受到信号序列长和k分辨参数的共同影响。对于某一信号，算法采用的函数基的k分辨值越大，运算量越大；相应地，对于相同k分辨函数基，信号序列越长则运算量越大。这一结论与算法复杂度为 $O(kN)$ 的结论相符。

5 结论

本文提出了一种正弦调频Fourier-Bessel变换，对变换的相关性质进行了推导，并将SFMFBT引入群目标微动频率提取中，并对离散信号处理中的相关问题进行了讨论。基于SFMFBT的微动频率提取方法利用全部信号序列长，相比于加窗算法其具有更高积累增益与提取精度，此外，该方法同样适用于群目标的微动频率提取与信号分离重构。下一步可将该方法与其他微动特征提取算法进行结合以提取更多的微动特征。

附录

由式(24)所述，当第m项与第m+1项SFMFB系数的实部绝对值相等，即

其中 $\tan {\varphi _m} = \left( {\mu _m^2 - {k^2}{T\;^2}{\omega ^2}} \right)/kT\;\omega $ 。由于第m项系数实部绝对值与第m+1项实部绝对值相等，将 $\omega = {\mu _m}/kT\;$ , $\omega = {\mu _{m + 1}}/kT\;$ 分别代入式(A-1)，则 ${\varphi _m}$ 计算为：

于是有

令 ${\gamma _m} = \displaystyle\frac{{{\mu _{m + 1}}\sin \left( {{\lambda _{m + 1}} - \phi } \right)/{J_{\alpha + 1}}\left( {{\mu _{m + 1}}} \right)}}{{{\mu _m}\sin \left( {{\mu _m} - \phi } \right)/{J_{\alpha + 1}}\left( {{\mu _m}} \right)}}$ ，经数值仿真可知 $\mathop {\lim }\limits_{m \to \infty } \left| {{\gamma _m}} \right| \to 1$ ，联立式(A-1)–式(A-4)，则有

即为式(24)所述。