随着电子对抗技术的高速发展，为了应对战场环境的变化，不断涌现出各种性能完善的新体制辐射源，电子侦察设备在同一时刻所截获的信号可能包含一个或多个不同辐射源信号. 在上述复杂电磁环境中，如何从多分量的脉冲信号中提取出敌方有用信号信息成为无源电子侦察领域迫切需要解决的问题.

多分量混合信号的分离及其频谱分析一直是学者们关注的问题，目前，关于信号分解的方法有很多. 经验小波变换（empirical wavelet transform，EWT）是Gilles^[1]在小波变换的基础上提出的可以自适应地对信号的频谱进行划分的方法. Bhattacharyya等^[2]指出，当信号的频谱结构很复杂时，EWT不能很好地对频谱进行分割. 经验模态分解（empirical mode decomposition，EMD）是Huang等^[3]总结出来的经验性信号分解类算法，但利用EMD对输入信号进行模态分解时，会产生“模态混叠”现象. 变分模态分解（variational mode decomposition，VMD）是Dragomiretskiy等^[4]在EMD的研究基础上提出的新的自适应时频分析的方法，与EMD相比，该方法在处理非平稳信号和抑制噪声等方面具有较突出的优势. 在实际场景中，信号形式多种多样，且雷达信号多数为大带宽的频率调制信号. 变分非线性调频模态分解（variational nonlinear chirp mode decomposition，VNCMD）是Chen等^[5]在研究VMD处理变分问题的基础上提出的针对非线性调频信号的处理方法；该算法摆脱了VMD中对本征模态函数带宽的限制，可以从多分量信号中分解提取出大带宽的调幅-调频信号，更适合于雷达信号的分析.

在时频分析领域，对于多分量信号的频谱分析，Reddy等^[6]针对多分量线性调频信号，利用魏格纳-霍夫变换估计信号的相关参数和对信号进行分离. 施闻明等^[7]提出将自适应高斯分布核的魏格纳分布算法应用在多分量调频信号的时频特性分析中. Khan等^[8]提出自适应分数阶谱图的方法，对多分量非线性调频信号进行分析. Stankovic等^[9]提出利用时频分布的L-统计特征，对多分量信号中非平稳分量进行提取. Chen 等^[10]根据时频分布图中“脊”的分布信息，利用“脊”重组（ridge path regrouping，RPRG）算法对“脊”路提取的瞬时频率进行重组校正，估计多分量混合调频信号的频率信息，结合固有调频分解算法对多分量混合信号进行分解.

本文针对大带宽的多分量雷达信号，提出基于脊路跟踪的VNCMD多分量信号分离算法. 利用VNCMD分解过程中各个分量瞬时频率不断逼近真实值的特点，将文献[10]的脊路重组算法进行相应的改进，估计输入信号的瞬时频率；将其作为VNCMD的频率预设值，通过不断更新VNCMD的频率预设值提高VNCMD算法的分离性能，完成对多分量信号的分离.

从时频分布图中提取出多分量信号的瞬时频率信息是多分量信号分离处理中关键的一步. 脊路重组算法^[10]是从多分量信号的时频分布图中检测多分量信号瞬时频率的脊线（局部最大值），根据脊线在交叉点处的变化率对脊线进行重新分组并连接修复，提取所需的瞬时频率.

为了有效地得到多分量信号的时频内容，脊路重组算法使用短时傅里叶变换得到信号的时间-频率分布图，定义为

(1) ${\rm{TF}}(t,f) = \int\limits_R {x(\tau ){g_\sigma }(\tau - t){{\rm{exp}}\;({ - 2{\rm{j}} \text{π}f\tau }}}) {\rm{d}}\tau .$

式中： $x(t)$为多分量信号； ${\rm{TF}}(t,f)$为得到的时频分布； ${g_\sigma }(\tau )$为长度为 $\tau $的高斯窗，信号的瞬时频率可以从时频分布中提取出来.

使用“脊”检测算法^[8,11]检测 ${\rm{TF}}(t,f)$中的所有脊线，定义 $\Delta f$表示2个连续点之间的最大可允许频率变化，找到时频分布中的瞬时频率最大值点：

(2) $({t_m},{f_m}) = \arg \max\; \left| {{\rm{TF}}(t,f)} \right|.$

式中： ${f_m}$为提取到瞬时频率的最大值， ${t_m}$为 ${f_m}$处的时间.

令 $\overline{{\rm{IF}}}({t_m}) = {f_m},{\bar f_{\rm{R}}} = {f_m},{\bar f_{\rm{L}}} = {f_m}$， $\overline{\rm{{IF}}}(t)$为频率估计值序列， ${\bar f_{\rm{R}}}$、 ${\bar f_{\rm{L}}}$为定义的中间变量， $t = {t_0}, \cdots ,{t_{N - 1}}$，N为信号截止时间，令 $r = m + 1,m + 2, \cdots N - 1$， $l = m - 1, m - 2, \cdots ,0$. 当 $r$与 $l$每次变化时，找到 ${t_r}$与 ${t_l}$时刻 ${\bar f_{\rm{R}}}$、 ${\bar f_{\rm{L}}}$在2 $\Delta f$范围内的频率最大值：

(3) ${\bar f_{\rm{R}}} = \mathop {\arg \max }\limits_{f \in [{{\bar f}_{\rm{R}}} - \Delta f,{{\bar f}_{\rm{R}}} + \Delta f]} \left| {{\rm{TF}}({t_r},f)} \right|,$

(4) ${\bar f_{\rm{L}}} = \mathop {\arg \max }\limits_{f \in [{{\bar f}_{\rm{L}}} - \Delta f,{{\bar f}_{\rm{L}}} + \Delta f]} \left| {{\rm{TF}}({t_l},f)} \right|,$

使 $\overline{{\rm{IF}}}({t_r}) = {\bar f_{\rm{R}}}$， $\overline{{\rm{IF}}}({t_l}) = {\bar f_{\rm{L}}}$， $r$与 $l$变化结束后得到一条脊线 $\overline{{\rm{IF}}}(t),t = {t_0}, \cdots ,{t_{N - 1}}$. 每得到一条脊线时，将这条脊线的所有瞬时频率置零，重复以上方法得到所有脊线 ${\overline{{\rm{IF}}}_1}(t), \cdots ,{\overline{{\rm{IF}}}_M}(t)$. 其中，M为脊线个数.

求任意2条重叠“脊”的所有交点区间，设 ${\bar t_s} < {\bar t_e}$为2个采样时刻， $s,e = 0, \cdots ,N - 1$，若满足

(5) $\begin{array}{l} \left| {{{\rm{{\overline{{IF}}}}}_k}(t) - {{\rm{{\overline{{IF}}}}}_m}(t)} \right| \left\{\begin{aligned} & \leqslant {D_{\rm{f}}},\;t \in [{{\bar t}_s},{{\bar t}_e}], \\ & > {D_{\rm{f}}},\;t = {{\bar t}_s} - 1/f_{\rm{{s}}},t = {{\bar t}_e} + 1/f_{\rm{{s}}}, \\ \end{aligned}\right. \end{array} $

则 $[{\bar t_s},{\bar t_e}]$称为2条脊线 ${\overline{{\rm{IF}}}_k}(t)$和 ${\overline{{\rm{IF}}}_m}(t)$的交点区间. 式中： ${D_{\rm{f}}}$为预定义阈值，阈值 ${D_{\rm{f}}}$控制交点区间的长度； $f_{\rm{s}}$为采样频率.

对这些交点区间进行编号，将第i个区间表示为

(6) ${{{{\it\Lambda }}}^{(i)}} = [{{{R}}^{(i)}},{{{T}}^{(i)}}] = [\{ R_1^{(i)},R_2^{(i)}\} ,\{ \bar t_s^{(i)},\bar t_e^{(i)}\} ].$

式中： ${{{R}}^{(i)}} = \{ R_1^{(i)},R_2^{(i)}\} \subseteq \{ 1, \cdots ,M\}$表示重叠“脊”的集合， ${{{T}}^{(i)}}{\rm{ = \{ }}\bar t_s^{(i)},\bar t_e^{(i)}\}$为间隔的起始和结束时刻的集合.

由于多条脊线相交于一点会出现多余的交点区间，合并在时频分布中彼此接近的交点区间，检测时频分布中任意2个交点区间的距离 ${\rm{dist(}}{\it\Lambda ^{(i)}},{\it\Lambda ^{(q)}})$. 若 ${\rm{dist(}}{\it\Lambda ^{(i)}},{\it\Lambda ^{(q)}}) \leqslant {{D}_{t - f}}$，其中 ${{D}_{t - f}}$是阈值，则2个区间可以合并为

(7) $\begin{array}{l} {\it\Lambda ^{(i)}} \oplus {\it\Lambda ^{(q)}} = [{{{R}}^{(i)}} \cup {{{R}}^{(q)}},\{ \min \;({{{T}}^{(i)}} \cup {{{T}}^{(q)}}), \\ \begin{array}{*{20}{c}} {}&{}&{}&{} \end{array}\max \;({{T}^{(i)}} \cup {{T}^{(q)}})\} ]. \\ \end{array} $

式中： $ \oplus $表示合并操作. 定义I为合并后 ${{{R}}^{(i)}}$中元素的个数，Q为合并后 ${{{R}}^{(q)}}$中元素的个数，则合并后2个交点区间的距离为

(8) $\begin{array}{l} {\rm{dist\;(}}{\it\Lambda ^{(i)}},{\it\Lambda ^{(q)}}) = \\ \sqrt {{{\left( {\bar t_e^{(i)} - \bar t_e^{(q)}} \right)}^2} + \Big[\sum\limits_{k = 1}^I {{{\rm{{\overline {IF} }}}_{R_k^{(i)}}}(\bar t_e^{(i)})/I - \sum\limits_{k = 1}^Q {{{\rm{{\overline {IF} }}}_{R_k^{(q)}}}(\bar t_e^{(q)})/{{Q}}{\Big]^2}} } } . \end{array}$

在式（8）中重复合并操作，直到任意2个交点区间的距离大于 ${{D}_{t - f}}$.

在每个交点区间重新组合“脊”，将检测到的脊线断开，根据区间中脊线的变化率重新连接.

在 $t_s^{ - (i)}$和 $t_e^{ - (i)}$时刻的“脊”的斜率（或变化率）可以近似为

(9) ${\rm{slp}}_{R_k^{(i)}}^ - = \Big[{\overline {\rm{{IF}}} _{R_K^{(i)}}}(\bar t_s^{(i)}) - {{\rm{\overline {IF}}} _{R_K^{(i)}}}(\bar t_s^{(i)} - \Delta t)\Big]\Big/\Delta t,$

(10) ${\rm{slp}}_{R_k^{(i)}}^ + = \Big[{{\rm{\overline {IF}}} _{R_K^{(i)}}}(\bar t_e^{(i)} + \Delta t) - {{\rm{\overline {IF}}} _{R_K^{(i)}}}(\bar t_e^{(i)})\Big]\Big/\Delta t.$

式中： $\Delta t > 0$是一个小的时间增量.

测量 $( - \infty ,\bar t_s^{\,(i)}]$和 $[\bar t_e^{\,(i)}, + \infty )$中不同“脊”段的匹配程度的连接矩阵元素定义为

(11) ${\rm{C}}{{\rm{M}}^{(i)}}(m,n) = \left| {{\rm{slp}}_{R_m^{(i)}}^ - - {\rm{slp}}_{R_n^{(i)}}^ - } \right|;\;m,n = 1, \cdots ,I.$

式中： $I$为区间 ${\it\Lambda ^{(i)}}$中交叉脊线的数量. 通过找到矩阵的最小元素来确定正确的连接对^[10]：

(12) $({m_0},{n_0}) = \mathop {\arg \min }\limits_{m,n \in \{ 1, \cdots ,I{\rm{\} }}} {\rm{C}}{{\rm{M}}^{(i)}}(m,{n}).$

$( - \infty ,\bar t_s^{\,(i)}]$中的“脊” ${{\rm{\overline {IF}}}_{R_{{m_0}}^{\,(i)}}}$和 $[\bar t_e^{\,(i)}, + \infty )$中的“脊” ${{\rm{\overline {IF}}}_{R_{{n_0}}^{(i)}}}$具有相似的变化率，应在这个交点区间内连接. 在 $(\bar t_s^{\,(i)},{{\rm{\overline {IF}}}_{R_{{m_0}}^{(i)}}}(\bar t_s^{\,(i)}))$和 $(\bar t_e^{(i)},{\overline {\rm{IF}}_{R_{{n_0}}^{(i)}}}(\bar t_e^{(i)}))$两时频点之间对脊线进行线性插值，在找到第一个连接对后，将矩阵 ${\bf{C}}{{\bf{M}}^{(i)}}$的 ${m_0}$行和 ${n_0}$列的所有元素设置为足够大的值，重复式（12）及式（12）后的操作，连接该区间内的其他“脊”. 对所有交点区间重复上述脊路重组算法，可以提取出各分量正确的瞬时频率脊线.

由脊路重组算法^[10]原理可知，在对多个分量进行瞬时频率脊线提取的过程中，脊路信息是逐条检测的. 当提取出第1条脊线时，利用调制和解调原理，将该条脊线包含的瞬时频率信息去除，剩余信号的时频信息会在原脊线交叉处出现能量缺口，可能会导致RPRG无法完整地提取出第2条脊线，进而导致提取的脊线轨迹出现严重的偏差.

针对以上问题，提出采用逐步检测法对脊路重组算法进行相应的改进. 当脊线出现能量缺口时，以脊线交叉处为中心 $O$，设置一个长为 $H$、宽为 $W$的矩形区域. 将该矩形沿着中心点 $O$，在 $0\sim 2\text{π}$区域内进行扫描，计算每一个角度下矩形区域所包含的能量总和，其中最大的能量所对应的角度即为脊线所对应的角度. 找到最大能量所对应的角度后，提出按下式所示的能量分布关系对矩形区域内的能量进行修补：

(13) $E(t,f) = \alpha {E_{\max}}\exp\;\Bigg[ - {\rm{j}}\left(\frac{{{{(t - {t_0})}^2}}}{{2{\sigma _t}^2}} + \frac{{{{(f - {f_0})}^2}}}{{2{\sigma _f}^2}}\right) \Bigg].$

式中： ${E_{\rm{{max}}}}$为时频分布图中的能量最大值， ${t_0}$和 ${f_0}$为中心点坐标， ${\sigma _t}$和 ${\sigma _f}$为标准差， $\alpha $为可调参数.将“脊”线连接后可以根据时频分布图中的“脊”信息提取分量的瞬时频率.

与VMD^[4]算法的分解结果相比，VNCMD^[5]算法的目的是将多分量信号分解为非线性调频模态的和的形式. 非线性调频模态为调幅调频信号，表达式为

(14) $c(t) = a(t)\cos\;\Big[2\text{π} \int_0^t {f(s){\rm{d}}s + \phi } \Big].$

式中： $a(t)$为 $c(t)$的瞬时幅度， $a(t) > 0$； $f(t)$为 $c(t)$的瞬时频率， $f(t) > 0$； $\phi $为 $c(t)$的初始相位. 与相位函数相比， $f(t)$和 $a(t)$的变化均是缓慢的.

在实际场景中，接收信号往往包含着多个非线性调频模态，且往往伴随着高斯白噪声，数学模型为

(15) $x(t) = \sum\nolimits_{i = 1}^K{a_i}(t)\cos \;\Big[2\text{π} \int_0^t {{f_i}(s){\rm{d}}s} + {\phi _i}\Big] + n(t).$

式中： $K$为信号中所包含的非线性调频模态的个数； $n(t)$表示高斯白噪声，均值为0，标准差为 $\sigma $，即 $n\sim {\rm{N}}(0,{\sigma ^2})$.

对于式（14）所示的非线性调频模态，解析形式^[12]可以表示为

(16) ${c_{\rm{A}}}(t) = a(t)\exp\;\Big[{\rm{j}}\Big(2\text{π} \int_0^t {f(s){\rm{d}}s} + {\phi _0}\Big)\Big].$

定义解调算子^[13]（demodulation operator，DO）为

(17) ${\cal{D}}(t) = \exp\;\Big[ - {\rm{j}}2\text{π} \Big(\int_0^t {{f_{\rm{d}}}(s){\rm{d}}s} - {f_{\rm{c}}}t\Big)\Big].$

调制算子^[13]（modulation operator，MO）为

(18) ${\cal{M}}(t) = \exp\;\Big[{\rm{j}}2\text{π} \Big(\int_0^t {{f_{\rm{d}}}(s){\rm{d}}s} - {f_{\rm{c}}}t\Big)\Big].$

当调制算子的瞬时频率 ${f_{\rm{d}}}(t) = f(t)$时，根据调制和解调的过程可知，式（15）可以表述为

(19) $\begin{split} x(t) =& \sum\nolimits_{i = 1}^K\Big[({u_i}(t)\cos \;\Big(2\text{π} \int_0^t {{{\tilde f}_i}(s){\rm{d}}s} \Big) + \\ & \begin{array}{*{20}{c}} {}&{} \end{array}{v_i}(t)\sin \;\Big(2\text{π} \int_0^t {{{\tilde f}_i}(s){\rm{d}}s} \Big)\Big] + n(t). \end{split} $

式中： ${u_i}(t)$、 ${v_i}(t)$为解调信号， $\{ {\tilde f_i}(t):i = 1,\cdots,K\} $为各个解调算子的瞬时频率.

根据 ${u_i}(t)$和 ${v_i}(t)$二阶导数的L2范数的平方估计原始信号的带宽^[14]，VNCMD算法中的变分问题可以表述为

(20) $\left.\begin{split}& \mathop {\min }\limits_{{{u}}_i , {{{v}}_i} , {{{\tilde f}}_i} } \sum\nolimits_{i = 1}^K(||{{{u}}_i}^{\prime \prime }||_2^2 + ||{{{v}}_i}^{\prime \prime }||_2^2) ;\\ & {\rm{s}}.{\rm{t}}.\;\Bigg|\Bigg|c(t) - \sum\nolimits_{i = 1}^K\Bigg[{u_i}(t)\cos \;\left(2\text{π} \int_0^t {{{\tilde f}_i}(s){\rm{d}}s} \right) + \\ &\quad {v_i}(t)\sin \;\left(2\text{π} \int_0^t {{{\tilde f}_i}(s){\rm{d}}s} \right)\Bigg]{\Bigg|\Bigg|_2} \leqslant \varepsilon . \\ \end{split}\right\} $

式中：上界 $\varepsilon $由噪声的强度决定， $\varepsilon = \sqrt {N{\sigma ^2}} $.

引入修正二阶差分算子 $\varOmega $，式（20）中约束性变分问题的离散形式可以等价于

(21) $\left.\begin{split}& \mathop {\min }\limits_{ {{{u}}_i} , {{{v}}_i} , {{\tilde {f}}_i} ,{{m}}} \{ {{\cal{I}}_{{{\cal{Q}}_\varepsilon }}}({{m}}) + {\sum\nolimits_i}(||\varOmega {{{u}}_i}||_2^2 + ||\varOmega {{{v}}_i}||_2^2)\}; \\ &\quad {\rm{s}}.{\rm{t}}.\;{{m}} = {{c}} - {\sum\nolimits_i}({{{P}}_i}{{{u}}_i} + {{{Q}}_i}{{{v}}_i}). \end{split} \right\}$

式中：

(22) $\left.\begin{array}{l}{{{P}}_i} = {\rm{diag}}\;[\cos\;({\varphi _i}({t_0})) ,\cdots, \cos\;({\varphi _i}({t_{N - 1}}))],\\ {{Q}_i} = {\rm{diag}}\;[{\rm{sin}}\;({\varphi _i}({t_0})), \cdots, \sin \;({\varphi _i}({t_{N - 1}}))]; \end{array}\right\}$

(23) $ \left.\begin{aligned} {\varOmega} = \begin{bmatrix}-1&1&0&\ldots&0\\ 1&-2&1&\ldots&0\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\ldots&1&-2&1\\ 0&\ldots&0&1&-1\\ \end{bmatrix};\\ {{\cal{I}}_{{{\cal{Q}}_\varepsilon }}}({{y}}) \triangleq \left\{ \begin{array}{l} 0,\quad {{y}} \in {{\cal{Q}}_\varepsilon }; \\ + \infty ,\quad{{y}} \notin {{\cal{Q}}_\varepsilon }; \\ \end{array} \right. \end{aligned}\right\} $

其中， ${{{\cal{Q}}_\varepsilon }} $表示以原点为中心、半径为 $ {\varepsilon}$的球体， ${{{{\cal{Q}}_\varepsilon }}} \triangleq \left\{ \begin{array}{l}\!\!\!\! {{q}}\in {{\bf{R}}^{{N}\times{1}} },\;||{{q}}||_2 \leqslant \varepsilon \!\!\!\!\end{array} \right\};$ ${{m}}$为辅助变量（或噪声变量），表示噪声对分解结果的影响.

根据VMD算法求解约束性变分问题的原理，VNCMD引入增广Lagrange函数，将上述约束性变分问题转化为非约束性变分问题，式（21）的增广Lagrange表达式^[15]如下：

(24) $\begin{split}& {L_\alpha }( {{{u}}_i} , {{{v}}_i} ,{{{f}}_i} ,{{m}},{{\lambda}} ) = \\ & \;\;\;\;{{\cal{I}}_{{{\cal{Q}}}_\varepsilon }}({{m}}) + {\sum\nolimits_i}(||\varOmega {{{u}}_i}||_2^2 + ||\varOmega {{{v}}_i}||_2^2) + \\ &\;\;\;\; ({\alpha }/{2})\Big|\Big|m + {\sum\nolimits_i}({{{P}}_i}{{{u}}_i} + {{{Q}}_i}{{{v}}_i}) - {{c}}^i + {{\lambda }}/{\alpha }\Big|\Big|_2^2 - ||{{\lambda}} ||_2^2/({{2\alpha }}). \\ \end{split} $

式中： $\alpha $为惩罚因子， ${{\lambda}} $为Lagrange乘子， ${{{{c}}}^i}$为信号分量.为了解决上述变分问题，VNCMD算法使用ADMM算法^[16]，求得 ${L_\alpha }$的极小值点.

对于噪声变量 ${{m}}$， ${{{m}}^{k + 1}}$可以表示为

(25) $\begin{split} {{{m}}^{k + 1}} =& \mathop {\arg \min\; }\limits_{{m}} \{ {L_\alpha }( {{{u}}_i} ,{{{v}}_i} , {{{f}}_i} ,{{m}},{{\lambda}} )\} =\\& \mathop {\arg \min }\limits_{{m}}\; \Big\{ {{\cal{{I}}}_{{{\cal{Q}}}_\varepsilon }}({{m}}) + ({\alpha }/{2})\times \\ &\Big|\Big|{{m}} + {\sum\nolimits_i}({{{P}}_i}{{u}_i} + {{Q}_i}{{v}_i}) - {{{c}}^i} + {{\lambda}} /{\bf{\alpha }}\Big|\Big|_2^2\Big\} . \\ \end{split} $

根据邻近算子^[17]的定义可知，式（25）的解为

(26) ${{{m}}^{k + 1}} = {{\cal{P}}_{{{\cal{Q}}_\varepsilon }}}(\hat{{c}} - {\sum\nolimits_i}({{{P}}_i}{{{u}}_i} + {{{Q}}_i}{{{v}}_i}) - {{\lambda} }/{\bf{\alpha} }).$

(27) ${{\cal{P}}_{{{\cal{Q}}_\varepsilon }}}{\rm{(}}{{x}}{\rm{)}} \triangleq \left\{ {\begin{array}{*{20}{c}} {\dfrac{\varepsilon }{{{{\left\| {{x}} \right\|}_2}}} \cdot {{x}}},&{{{\left\| {{x}} \right\|}_2} > \varepsilon }; \\ {{x}},&{{{\left\| {{x}} \right\|}_2} \leqslant \varepsilon }. \end{array}} \right.$

式中： $\hat{{c}} $为剩余的信号分量.

对于2个正交的解调信号 ${{{u}}_i}$和 ${{{v}}_i}$的更新问题，可以求得 ${{u}}_i^{k + 1}$和 ${{v}}_i^{k + 1}$的更新表达式为

(28) ${{u}}_i^{k + 1} = {\left(\frac{2}{\alpha }{\varOmega ^{\rm{T}}}\varOmega + {({{P}}_i^{k + 1})^{\rm T}}{{P}}_i^{k + 1}\right)^{ - 1}}{({{P}}_i^{k + 1})^{\rm{T}}}\hat{{c}},$

(29) ${{v}}_i^{k + 1} = {\left(\frac{2}{\alpha }{{\varOmega} ^{\rm{T}}}\varOmega + {({{Q}}_i^{k + 1})^{\rm{T}}}{{Q}}_i^{k + 1}\right)^{ - 1}}{({{Q}}_i^{k + 1})^{\rm{T}}}\hat{{c}}.$

式中：

(30) ${{P}}_i^{k + 1} = {\rm{diag}}\;\Big[\cos \;\Big(2\text{π} \int {f_i^{k + 1}(t){\rm{d}}t} \Big)\Big],$

(31) ${{Q}}_i^{k + 1} = {\rm{diag}}\;\Big[\sin \;\Big(2\text{π} \int {f_i^{k + 1}(t){\rm{d}}t} \Big)\Big].$

根据式（28）、（29），可得

(32) ${{c}}_i^{k + 1} = {{P}}_i^{k + 1}{{u}}_i^{k + 1} + {{Q}}_i^{k + 1}{{v}}_i^{k + 1},$

(33) ${{f}}_i^{k + 1} = {{f}}_i^k + \gamma \Delta {{f}}_i^{k + 1}.$

式中： $\gamma $为比例因子， $0 < \gamma < 1.0$.

对于Lagrange乘子 ${{\lambda}} $的更新问题， ${{{\lambda}} ^{k + 1}}$可以表述为

(34) ${{{\lambda}} ^{k + 1}} = {{{\lambda}} ^k} + \alpha ({{{m}}^{k + 1}} + \sum {{c}}_i^{k + 1} - \hat{{c}}).$

VNCMD算法^[5]的具体实施过程如下.

算法1：VNCMD算法^[5]

输入：多分量信号 $x(t)$，瞬时频率估计值 ${ f}_i^1$

输出：信号分量 ${ c}^{i}$

1) 初始化 $\alpha $， $\mu $， ${ f}_i^1$， ${ P}_i^1$， ${ Q}_i^1$， ${ u}_i^1$， ${ v}_i^1$， ${{\lambda} ^1} ={\bf 0}$， $k = 0$；

2) 令 $k = k + 1$；

3） 更新 ${{ m}^{k + 1}}$，初始化 $i = 0$;

4） 令 $i = i + 1$；

5） 更新 ${ u}_i^{k + 1}$、 ${ v}_i^{k + 1}$、 ${ c}_i^{k + 1}$和 ${ f}_i^{k + 1}$；

6） 重复步骤3）、4），直到 $i$等于信号分量总数 $K$；

7） 更新 ${{\lambda} ^{k + 1}}$；

8） 当条件 $||{{{m}}^{k + 1}} + {\sum\nolimits_i}{{c}}_i^{k + 1} - {{c}}||_2^2 > ||{{c}}||_2^2$成立时，执行步骤9）~12）；

9） 初始化 ${{\lambda} ^{k + 1}} = {\bf 0}$， $i = 0$；

10） 令 $i = i + 1$；

11） 更新 ${ u}_i^{k + 1}$和 ${ v}_i^{k + 1}$；

12） 重复步骤10）~12），直到 $i$等于信号尺度数 $K$；

13） 重复步骤2）~12），直到满足约束条件 ${\sum\nolimits_i}\Big|\Big|{{c}}_i^{k + 1} - {{c}}_i^k\Big|\Big|_2^2 \Big/ {\Big|\Big|{{c}}_i^k\Big|\Big|_2^2} \leqslant \delta$.

原信号经过VNCMD分解，得到对应于原始信号中每个信号分量的非线性调频模态.

VNCMD算法^[5]在分解过程中，通过迭代更新各个信号分量的瞬时幅度和瞬时频率，达到对多分量混合信号进行分离的目的. 实际上，该算法在每一次的迭代更新过程中得到的结果均是对原始信号中各个信号分量瞬时频率和瞬时幅度的一次逼近. 在该过程中，瞬时频率的误差是无法预知的，因此在VNCMD算法^[5]的基础上，提出基于脊路跟踪的VNCMD（ridge track VNCMD，RT-VNCMD）算法. 该算法在每次分解过程中，在一次VNCMD分解完成之后，对分解得到的各个信号分量求和；利用改进的脊路重组算法估计重构信号的瞬时频率，更新VNCMD的频率预设值；重新对原始信号进行分解，达到“脊”路跟踪的目的，减小分解过程中的误差.

引入相关误差（reference error，RE）和均方根误差（root mean square error，RMSE），衡量VNCMD分解得到的各个模态的瞬时频率和瞬时幅度的误差. 瞬时频率的相关误差与瞬时幅度的均方根误差表达式如下：

(35) $\left. \begin{array}{l} {\rm{R}}{{\rm{E}}_i} = {{{{\left\| {{{\tilde {{f}}}_i} - {{{f}}_i}} \right\|}_2}}}\Big/{{{{\left\| {{{{f}}_i}} \right\|}_2}}}, \\ {\rm{RMS}}{{\rm{E}}_i} = 10\lg \sqrt {\displaystyle\sum\limits_{t = 1}^N {{{\left| {{c_i}(t) - {{\tilde c}_i}(t)} \right|}^2\Bigg/N}} } . \\ \end{array} \right\}$

式中： ${{{f}}_i}$为多分量信号中某信号分量的瞬时频率， ${c_i}(t)$为瞬时幅度， ${\tilde { f}_i}$为分解得到的相应分量的瞬时频率， ${\tilde c_i}(t)$为分解得到的相应分量的瞬时幅度， $N$为信号长度. 对于 $K$分量的多分量信号 $x(t)$，设在一次RT-VNCMD分解过程中共经历J次VNCMD分解，则RT-VNCMD分解算法的流程如图1所示.

具体的计算过程如下.

算法2：RT-VNCMD算法

输入：多分量信号 $x(t)$

输出：信号分量 ${{c}} = [c_J^1(t),c_J^2(t), \cdots ,c_J^K(t)]^{\rm{T}}$

1） 利用STFT变换获得信号 $x(t)$的时频分布图，初始化 $j = 0$.

2） 利用改进的脊路重组算法估计信号 $x(t)$中各个分量的瞬时频率 ${{{F}}_j}(t) = [f_j^1(t),f_j^2(t),\cdots,f_j^K(t)]^{\rm{T}}$，将 ${{{F}}_j}(t)$设置为VNCMD的频率预设值.

3） 对信号 $x(t)$进行VNCMD分解，获得各个信号分量的瞬时频率 ${{{F}}_{j + 1}}(t) = [f_{j + 1}^1(t),f_{j + 1}^2(t),\cdot\cdot\cdot ,f_{j + 1}^K(t)]^{\rm{T}}$.

4） 计算 ${{{F}}_{j + 1}}(t)$中所有分量与 ${{{F}}_j}(t)$中相应分量的相关误差，记为 ${\bf{R}}{{\bf{E}}_j}$，且 ${\bf{R}}{{\bf{E}}_j} = [e_j^1,e_j^2,\cdots,e_j^K]^{\rm{T}}$.

5） 对步骤3）中VNCMD获得的各个信号分量求和，得重构信号 $\tilde x(t)$.

6） 估计重构信号 $\tilde x(t)$中每个信号分量的瞬时频率，并重新将其设置为VNCMD的频率预设值.

7）令 $j = j + 1$，重复步骤3）~6），直到 ${\bf{R}}{{\bf{E}}_{j}}$中相关误差均小于0.001为止，此时迭代次数记为J.

为了验证提出算法的分离效果，采用下式所示的由线性调频信号和正弦调频信号组成的混合信号进行仿真分析：

(36) $\left. \begin{array}{l} {s_1} = {A_1}{{\rm{exp}}\;\left\{{{\rm{j}}[2\text{π} ({{{f_{01}}}}t + {\rm{0.5}}{k_1}{t^2})]}\right\}}, \\ {s_2}\!\! = \!\!{A_2}{{\rm{exp}}\;\left\{{{\rm{j}}[2\text{π} ({f_{02}}t \!\!+ 0.5{k_2}{t^2} \!\!+ {m_{\rm{f}}}{\rm{sin}}\;(2\text{π} {f_{\rm{m}}}t))]}\right\}}, \\ {s_3} = {A_3}{{\rm{exp}}\;\left\{{{\rm{j}}[2\text{π} ({f_{03}}t + 0.5{k_3}{t^2})]}\right\}}. \end{array} \!\!\right\}$

式中：A₁、A₂、A₃为信号幅值， ${s_1}$表示初始频率 ${f_{01}}$=15 MHz，带宽 ${B_1}$=5 MHz的线性调频信号，脉冲宽度 ${T_1}$=10 μs，调频斜率 ${k_1} = {B_1}/{T_1}$； ${s_2}$表示线性频率调制和正弦频率调制的混合调制信号，线性调频部分的初始频率 ${f_{02}}$=1 MHz，带宽 ${B_2}$=40 MHz，脉冲宽度 ${T_2}$=10 μs，调频斜率 ${k_2} = {B_2}/{T_2}$，正弦调频部分的调制指数 ${m_{\rm{f}}}$=5，调制频率 ${f_{\rm{m}}}$=0.2 MHz； ${s_3}$表示初始频率 ${f_{03}}$=25 MHz，带宽 ${B_3}$=25 MHz的线性调频信号，脉冲宽度 ${T_3}$=10 μs，调频斜率 ${k_3} = {B_3}/{T_3}$； $\eta $表示高斯白噪声，均值为0，标准差为 $\sigma $，采样率为100 MHz.

在多分量信号进行分解之前，需要在多分量信号的时频分布图中估计出各分量信号的瞬时频率，作为VNCMD的预设频率. 为了验证提出的改进脊路重组算法的性能，采用式（36）所示 ${s_1}$和 ${s_2}$组成的多分量信号，在不同信噪比下分别使用未改进和改进的脊路重组算法，从多分量时频分布图中估计信号的瞬时频率，在每个噪声强度下开展1 000次蒙特卡洛实验，结果如图2所示.图中，RE_f为瞬时频率相关误差.

从图2可以看出，对于信号分量 ${s_1}$，当信噪比大于5 dB时，改进的脊路重组算法的估计误差略低于原始算法；当信噪比小于5 dB时，改进的脊路重组算法的估计误差略高于原始算法. 这是由于在实验过程中，不同噪声强度下对噪声的能量进行判定时采用统一的门限标准，因此在低信噪比的条件下，改进算法的性能略低于原始算法，但是两者之间相差不多. 对于信号分量 ${s_2}$，与原始算法相比，在不同的噪声强度下，改进的脊路重组算法的估计误差远远低于原始算法；当信噪比大于5 dB时，改进的脊路重组算法的估计误差接近于0. 利用改进的脊路重组算法可以在整体上提升算法的频率估计性能.

采用由 ${s_1}$和 ${s_2}$组成的二分量混合信号，保持信号的其他参数不变，在不同噪声强度下进行仿真实验. 分别计算RT-VNCMD和VNCMD分解得到的各个分量瞬时频率的相关误差和各个分量与原始信号中相应分量瞬时幅度的均方根误差RMSE，在每个噪声强度下均进行1 000次蒙特卡洛实验，结果如图3、4所示. 经过RT-VNCMD和VNCMD分解后得到的各分量时域图如图5所示，为了便于分析，将时域图截取到3 μs处.

从图3、4可以看出，在不同信噪比下，RT-VNCMD分解得到的信号分量 ${s_1}$和 ${s_2}$的瞬时频率的相关误差均小于VNCMD的分解结果，RT-VNCMD分解得到的2个分量瞬时幅度的均方误差小于VNCMD的分解结果，在低信噪比下RT-VNCMD的分离效果更好；从图5可以看出，当SNR=15 dB时，由RT-VNCMD与VNCMD分解得到的各分量时域波形图与原信号时域波形图几乎完全重合，可以很好地将2个信号分量恢复出来；总体来说，在相同的噪声强度下，RT-VNCMD对多分量信号的分离性能优于VNCMD对多分量信号的分离性能.

为了进一步验证RT-VNCMD对多分量信号的分离效果，采用式（36）中由 ${s_1}$、 ${s_2}$和 ${s_3}$组成的三分量混合信号，保持信号的其他参数不变，在不同噪声强度下进行仿真实验. 分别计算RT-VNCMD和VNCMD分解得到的各个分量瞬时频率的相关误差和各个分量与原始信号中相应分量瞬时幅度的均方根误差，在每个噪声强度下均进行1 000次蒙特卡洛实验，结果如图6、7所示. 经过RT-VNCMD和VNCMD分解后得到的各分量时域图如图8所示.

从图6、7可以看出，当信噪比较高时，由RT-VNCMD分解得到的 ${s_3}$分量的瞬时频率的相关误差大于由VNCMD分解的结果；当信噪比较低时，RT-VNCMD的分离效果更好. 从图8可以看出，当SNR=15 dB时，由RT-VNCMD与VNCMD分解得到的各分量时域波形图与原信号时域波形图几乎完全重合，可以很好地将信号恢复出来. 总体来说，由RT-VNCMD和VNCMD分解得到的信号分量 ${s_1}$、 ${s_2}$和 ${s_3}$的瞬时频率的相关误差与瞬时幅度的均方根误差较对二分量混合信号的分离结果大一些，但从实验结果可以看出，不论是对二分量混合信号或三分量混合信号进行分离，RT-VNCMD对多分量信号的分离性能都优于VNCMD对多分量信号的分离性能.

采用某实测雷达信号数据，加入强度相等的大带宽频率调制信号对该实测雷达信号数据进行干扰. 混合信号的时域波形和STFT获得的二维时频信息如图9、10所示. 图中，A为幅度.干扰信号 ${s_j}$为线性频率调制和正弦频率调制的混合调制信号，具体的数学表达式为

(37) ${s_j} = {{\rm{exp}}\;\left\{{{\rm{j}}\left[2\text{π} ({{\hat f}_0}t + 0.5\hat k{t^2} + {{\hat m}_{\rm{f}}}\sin\;(2\text{π} {{\hat f}_{\rm{m}}}t))\right]}\right\}}.$

${s_j}$中线性调频部分的初始频率 ${\hat f_0}$为20 MHz，带宽 $\hat B$为60 MHz，脉冲宽度为6.24 μs，调频斜率 $\hat k = \hat B/\hat T$；正弦调频部分的调制指数 ${\hat m_{\rm{f}}}$=3，调制频率 ${\hat f_{\rm{m}}}$=0.5 MHz. 干扰信号中携带有信噪比为14 dB的高斯白噪声. 信号的采样频率为200 MHz.

在进行信号分解之前，根据时频分布图中的能量“脊”的信息，估计多分量信号的瞬时频率，并将此瞬时频率作为VNCMD和RT-VNCMD的频率预设值，经过RT-VNCMD和VNCMD分解后得到的各分量时域图如图11所示. 可以看出，对实际信号的分离效果较理想.

为了说明该算法的性能指标，分别计算VNCMD和RT-VNCMD分解得到的干扰信号瞬时频率的相关误差和瞬时幅度的均方根误差，结果如表1所示.

从表1可以看出，与VNCMD的结果相比，RT-VNCMD分解得到的干扰信号瞬时频率的相关误差大约下降了63%，瞬时幅度的均方根误差大约下降了45%. 与VNCMD相比，RT-VNCMD在实际信号的应用方面可以取得更好的分离效果.

本文针对大带宽多分量雷达信号的分离问题，提出基于脊路跟踪的VNCMD信号分离算法，即使用改进后的脊路重组算法提取出多分量信号的瞬时频率并将其作为VNCMD的预设频率. 利用VNCMD更新迭代时可以分解得到混合信号中各信号分量的原理，完成对多分量信号的分离. 通过仿真实验及使用实际雷达信号的实验表明，在对二分量混合信号进行分离时，使用提出的基于脊路跟踪的VNCMD多分量信号分离算法的分离效果优于VNCMD算法的分离效果，在对三分量混合信号进行分离时，虽然整体分离效果不如对二分量混合信号的效果好，但使用提出的基于脊路跟踪的VNCMD多分量信号分离算法的分离效果优于VNCMD算法. 由于实际情况中出现时域频域同时重叠且信号分量超过3个的混合信号概率较低，对信号分量超过3个时的混合信号暂不作讨论. 对实际信号的实验结果验证了基于脊路跟踪的VNCMD多分量信号分离算法的有效性，且分离效果优于使用VNCMD算法的分离效果.