矢量序列信号的瞬时频率分析具有重要的应用需求，诸如雷达海杂波分析、海洋声信号分析等。瞬时频率计算对噪声污染非常敏感，因此矢量序列噪声抑制是不可回避的第一步环节。但是，从噪声污染信号中提取信号，面临的是一个“病态”投影信号复原的挑战性问题。对于海杂波尤其如此。噪声污染又加上海面结构所导致的信号非线性和非平稳性，更加剧了信号分析的困难。

现代信号处理越来越依赖于概率和统计以解决信号处理中的挑战问题^[1]。对于矢量形式的高斯噪声，人们认识到它的相位服从均匀分布(白色)而矢径服从高斯分布。借助信息熵的术语，均匀分布的信息熵最大，也意味着它的随机性最强，因此需要最为精细尺度的处理方法。而对于服从非均匀分布的矢径序列，有效的途径是处理算法的尺度能够依赖于数据的局部性尺度。

纵观近年来已有的分解处理方法，全变分TV(Total Variation)方法具有最为精细的时间变化分析性能。全变分降噪模型^[2]自从提出之后，由于物理意义的明显和精细处理特点，在图像分解、复原、分割等多方面得到广泛应用，目前已经形成了一种信号和图像处理领域的艺术级别的流派。TV方法的数字化计算过程中有2个参数，调整因子和迭代次数，需要进行选择，文献[2]给出了这参数调整因子的计算公式。而后，诸多学者研究了对处理效果的影响、诸如图像的边缘保持特性等。文献[3]在超声信号降噪的TV应用中认为该参数计算的耗费太大，经验性的给出了该参数的选择范围。

1998年，Nudun Huang等^[4]提出了适应于非平稳非线性信号的EMD(Empirical Mode Decomposition)分解方法。EMD将信号看作为“零均值快变震荡信号与慢变震荡信号的叠加”，其分解的主要特点是依据待分解数据的“局部尺度”，将原始信号分解为不同震荡频率的细节和內禀性模式IMFS。而后，通过Hilbert-Huang Spectrum计算信号的瞬时频率。文中也证明了利用分解量可以完全重构原始信号。但是，Nudun Huang的EMD方法仅仅适用于1D实数序列分解。

为了能够进行矢量序列分析，近年来逐步发展了RIEMD^[5]、BEMD^[6]等基于EMD的分解方法。BEMD认为当分析数据是两变量数据时，振动的表示是模糊的。他们以旋转表示两变量信号：“两变量信号＝快速旋转信号叠加于慢速旋转信号”之上。文献[6]同时给出了包络均值的求取算法。文献[7]注意到，对于规定两个方向的情况RIEMD与BEMD是等价的。文献[8]应用BEMD研究了在认知雷达场景分析的应用。

根据噪声统计特性及TV和EMD分析方法的特点，本文提出由TV相位降噪和矢径EMD分解相组合的矢量序列降噪、分解方法。在相位降噪中，以最大信杂比为准则，优化选择迭代计算中的调整因子和迭代次数，以获得最大信杂比的降噪效果。在矢径EMD分解中，为了防止重构矢量的相位跳变，对分解增加了非负判定环节。最后，将经过降噪的相位序列和EMD分解的矢径序列，以采样时间(样点序号)为参考，一一对应组合，重构矢量。考虑到历史的传承性，我们将该种矢量分解方法称之为VEMD(Vector Empirical Mode Decomposition)方法。

1 矢量时间序列的分解 1.1 相位序列TV降噪

任意一个矢量时间序列可以表示为

$ Z(t)=a(t)\exp\{jφ(t)\}。$

(1)

其中：a(t)是矢量的矢径；φ(t)是矢量在时刻t的相位。

将噪声污染的相位序列φ(t)表达为

$ φ=y+n。$

(2)

式中n为噪声。TV降噪模型为

$ \arg {\min _s}({\smallint _\mathit{\Omega }}[\nabla \varphi + \frac{\lambda }{2}{\left( {\varphi - y} \right)^2}]){\rm{d}}x。$

(3)

式中Ω信号定义域，调整参数λ>0控制平滑量。(3)式中的第1项称为正则项，第2项称为保真项。(3)式最速下降法的Euler-Lagrange方程是

$ \frac{{\partial \varphi }}{{\partial t}} = {\rm{div}}\left( {\frac{{\nabla \varphi }}{{\left| {\nabla \varphi } \right|}}} \right) + \lambda \left( {\varphi - y} \right)。$

(4)

但对于单变量信号，由$ \frac{{\partial \varphi }}{{\partial t}} = 0$，可导出离散化的迭代计算公式为：

$ {y^{n + 1}} = \varphi - \frac{1}{\lambda }y_{tt}^n。$

(5)

式中y_ttⁿ为第n次迭代yⁿ的二阶差分，初始化y⁰=φ，迭代过程结束后，以yⁿ⁺¹作为目标函数φ_Denoise。文献[9]中给出了yⁿ迭代计算的MATLAB程序网址。

TV降噪的计算过程涉及调整因子和迭代次数的选择。以自适应方法选择将会导致大计算量耗费。最大信杂比是信号处理所追求的一般性准则。本文采用最大信杂比为准则的方法指导和迭代次数的选择。信杂比估计采用了经验性估计方法。假设信号频率低于噪声频率，将低频功率谱峰面积与背景杂波的功率谱面积之比作为估计的信杂比。

最大信杂比TV相位降噪及参数过程如下：首先设定λ和的初始数值及其变化步长；在每一次迭代计算之后，估计信杂比。最后，以最大信杂比下的调整因子和迭代次数作为参数，应用于TV降噪。一个实际过程的信杂比估计的图形如1所示。

1.2 修改的EMD矢径分解

考虑到矢量重构，关于矢径分解的一个重要约束是不能够改变原数值的符号(正/负)。因为改变一个矢量的正负，等价于矢量相位的变化。本文在原始EMD分解框架的基础上，增加一个IMF的判别环节，∀lmf(t_i) < 0，则分解停止。

EMD分解框架：

记矢径序列为{x(i)}_i=1^N。

(1) 寻{x(i)}_i=1^N的局部极大值/极小值：使得局部极大和极小值的数量相等或最多差1。

(2) 求局部包络：在相邻的极大值之间(极小值之间同样处理)，利用立方样条插值，形成极大值包络e_max(和极小值包络e_min)。

(3) 局部均值曲线：对每一个样点，均值曲线：

$ Imf\left( i \right) = \frac{1}{2}[{e_{\max }}\left( i \right) + {e_{\min }}\left( i \right)]。$

(6)

(4) Imf非负判断：如果∀Imf(t_i) < 0, 停止。

(5) 从信号中减均值，得到分解的细节Detail(i)=x(i)-Imf(i)。

经过上述过程，从原始信号得到第一层分解的Imf₁和Detail₁。

将Imf₁作为新的待分解信号，重复(1)~(5)，直到L层的Imf_L。其中，(4)条为本文加入的约束条件。

1.3 矢量重构

将VEMD分解的最后一层内禀性模式与降噪相位序列φ_Denoise一一对应，构成矢径、相位对，则称为重构的分解矢量序列，

$ {z_{re}}\left( i \right) = \{ im{f_L}\left( i \right){{\rm{e}}^{j{\varphi _{Denoise}}\left( i \right)}}\} _{i = 1}^N。$

(7)

它保留了原始信号的最大能量并且满足TV相位保真性约束条件，因此将其看作为纯净信号的近似。

1.4 瞬时相位梯度与瞬时功率谱计算

瞬时频率的计算在历史上是一个有争议的问题。概念上，瞬时频率定义为$\frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}} $。瞬时频率是频率域最为精细的表示，但是其精细性同时也会导致一些物理意义上不可解释的现象。例如，可能会出现负的“瞬时频率”，以及某一时刻的瞬时频率对此时刻前后的影响等问题。一般来说，分析结果的“鲁棒性”与精细性是一对互相矛盾的两个侧面，追求精细性必然会损失“鲁棒性”。时频联合分析的测不准原理是很好的例证。为了避免不必要的争议，我们这里将传统的“瞬时频率”理解为瞬时相位梯度。在规定了方向之后，负值的瞬时相位梯度表示与规定的方向相反，而其物理意义又可以近似为瞬时频率。对于具有非平稳性和非线性的信号，瞬时相位梯度也对应于瞬时频率意义的精细刻画。

对重构的，通过(8)式计算瞬时相位梯度。为了与瞬时频率对应，乘以因子$\frac{1}{{2{\rm{ \mathsf{ π} }}}} $。

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;f({t_i}) = \frac{1}{{2{\rm{ \mathsf{ π} }}}}\frac{{\Delta \varphi ({t_i})}}{{\Delta t}} = \\ \frac{1}{{2{\rm{ \mathsf{ π} }}}}[Angle({Z_{re}}({t_i})) - Angle({Z_{re}}({t_{i + 1}}))]。\end{array} $

(8)

为了计算瞬时功率谱，在(8)式的基础上，将

$ \begin{array}{l} P({t_i}) = \frac{1}{2}[P({t_i}) + P({t_{i + 1}})] = \\ \frac{1}{2}[im{f^2}\left( i \right) + im{f^2}\left( {i + 1} \right)]。\end{array} $

(9)

赋给联合坐标，作为瞬时功率谱。

2 分解实验

为了验证本文所提出方法的有效性，我们以人工产生的线性频率调制信号，线性+正弦频率调制的信号进行了分解噪声实验。实验过程如下：首先生成纯净的矢量时间序列信号Z_pure(t)，然后利用MATLAB功能函数，在信号上叠加信噪比为1dB的高斯白噪声(请注意这是信杂比较低的情况)，然后利用本文的VEMD分解方法进行分解。为了表示的统一性，记矢径分解的最后一层为Imf(t)。将Imf(t)与经过TV降噪的相位序列一一对应组合，构成近似矢量Z_re(t)=Imf(t)exp{jφ_Denoise(t)}，计算瞬时频率。

实验用例1，线性频率调制信号Z_pure(t)=3.0exp{j2π300(t²)/2}，实验结果如图 2所示。

实验用例2，线性频率调制+正弦频率调制信号，过程同实验用例1。

Z(t)=3exp{j2π300(t²)/2+j5sin(2π10t)}实验结果如图 3所示。

由上述2个实验结果的可见，本文的方法可以有效抑制矢量高斯噪声。虽然在瞬时频率的初始和最后阶段有所形变，但是基本上保持了纯净信号的变化规律。

实验用例3，加拿大IPIX雷达数据这里采用的试验数据来自于加拿大McMaster大学的IPIX雷达数据库网站，其雷达数据包括了HH、HV、VH，VV四种极化方式，每一个数据文件含有14个距离门的回波信号，其中目标为一个直径1 m的球形密封器件，表面包裹了铝箔，以增强反射回波信号。

图 4是IPIX实测数据的降噪实验结果，使用了IPIX数据文件31#HV极化方式的1 500个样点，主要用以实验本文方法的有效性。由图 4可见，本文的分解降噪方法的确较为有效的抑制了噪声干扰。以此为基础，可进行降噪数据的分析，建模等研究。由于时间限制，这些研究工作有待后续进行。

3 结语

信号的统计特性是处理方法选择的重要依据。高斯噪声的相位服从均匀分布，意味着其具有最强的随机性。因此，本文利用细致的TV方法对矢量序列的相位进行降噪；又根据矢量的矢径序列一般服从非均匀分布特点，利用适用于非线性非平稳序列的EMD方法进行矢径分解；据此形成一种新的矢量序列分解、降噪方法。这种针对不同统计特性“分而治之”的信号处理方法，希望能够为读者提供一定的参考。以此为基础，后续将进行降噪之后的数据分析，建模等研究。