一元及多元信号分解发展历程与展望

陈启明 文青松 郎 恂 谢 磊 苏宏业

在一个飞速发展的信息社会中,信号是无处不在的,例如人们随时可以听到的语音信号、随时可以看到的视频信号、伴随生命始终存在的生理信号、新冠肺炎感染患者的数量、工业生产中的控制监测信号、黑洞碰撞的引力波信号等[1].这些信号是变化的,变化的信号构成了五彩斑斓的世界,如何描述这些变化的信号以及揭示这些信号中潜在的信息,是信号处理研究的重要任务[1].一般而言,在传统的信号处理中,人们会从时域或频域两个角度去描述信号.但单纯从时域或频域角度出发描述信号较为片面,无法刻画非平稳和非线性信号的重要性质.在现实世界中,所获得的信号大部分都是非平稳和非线性的,为了满足现实信号处理需求,时频联合分析方法相关研究应运而生[2].

1998 年,Huang 等[3]提出一种自适应的非平稳非线性信号分解方法,成为信号处理发展的一个重要里程碑,随后很快发展为“后小波”时代时频分析领域的研究热点[4].信号分解假设复杂的非平稳非线性信号由多个简单的子信号组成,通过分析这些子信号的特征,可以间接或直接揭示原始复杂信号的时频信息,进而为各个领域的信号处理任务提供有力工具.近年来,多元/多变量/多通道信号分解理论方兴未艾,在诸多领域得到了成功应用,但目前尚未见到相关综述报道.为了填补这个空缺,本文从单变量和多变量两个方面对主流信号分解方法及其主要改进进行综述,分析比较这些方法的原理和优缺点,并在最后进行了研究展望.

1 传统信号分解与时频表示

信号处理的一个重要任务就是分析信号中的频谱成分,同时还希望知道不同频率成分所出现的时间.传统的信号处理方法主要集中于信号变换,例如通过傅里叶变换(Fourier transform,FT)及其反变换,人们构建起了信号时域与频域之间转换的桥梁.以傅里叶变换为基础的分析方法已经统治了线性时不变与平稳信号处理领域近200 年,尤其是Cooley 等[5]利用傅里叶算子的周期性和对称性,提出了快速傅里叶变换算法,将N点傅里叶变换的乘法计算量从N2次降为 (N/2×log2N) 次,这成为数字信号处理发展史上的转折点和里程碑.以此为契机,伴随着超大规模集成电路和计算机领域的迅猛发展,傅里叶变换不但已成为一个重要的数学分支,而且也成为信号分析和处理的重要工具,并在众多领域得到了广泛应用.但是,在科学研究与工程应用中,研究人员发现,傅里叶变换具有较多不足,主要体现在三个方面[1]: 1)傅里叶变换缺乏时间和频率的定位功能.傅里叶变换得到的结果是信号在整个积分区间的时间范围内所具有的频率特征的平均表示,无法通过傅里叶变换知道在某一个特定时刻或较短时间范围内的信号频率信息.2)傅里叶变换对于非平稳信号的局限性.只有时不变的信号才能够展开为无穷多个复正弦函数的和,而且这无穷多个复正弦信号的幅度、频率和相位都不随时间变化,即取某一特定常数.因此,傅里叶变换只适合于处理平稳的时不变信号,但从实际过程中采集到的数据往往具有时变特性[1],属于非平稳信号[5-6].3)傅里叶变换在分辨率上的局限性.傅里叶变换受到不定原理的制约[1],无法根据信号的特点来自适应调整时域分辨率和频域分辨率.但是,正是傅里叶变换的这些不足成为了几十年来推动人们寻找新的信号分析与处理方法的动力.

Gabor[7]提出短时傅里叶变换(Short-time Fourier transform,STFT)进行时域和频域的联合分析.这种方法简单易懂,但是不能自动调节时域窗口和频域窗口.Cohen[8]给出了各种时频分布的统一表示形式,称为Cohen 类时频分布,其中最核心的就是Wigner-Ville 分布(Wigner-Ville distribution,WVD).Wigner 分布是由Wigner[9]在1932年提出的概念,但直到1948 年,Ville 才将其应用于信号分析,因此命名为Wigner-Ville 分布.实际上,WVD 可理解为在一个特定区间上的傅里叶变换,因此它仍然受到不定原理的制约.WVD 的时频分析性能也易受交叉项的影响.此外,WVD 的瞬时频率是时间的单值函数,在处理多分量信号时,只能给出多个频率在一个时间点上的均值,这样就无法刻画多分量信号的频率分布.小波变换(Wavelet transform,WT)是过去20 多年信号处理领域最重要的进展之一[5],它最大的特点是在基本小波中引入了尺度因子,使得小波具有自动调节时域和频域分辨率的能力.然而,在对信号进行小波变换前如何选择小波基函数,仍是一个未解决的现实问题,一般需要使用者不断试用各类小波.

2 一元信号分解

上述传统信号处理方法都是基于基函数展开的思路,具有简单、唯一和对称等优点[10],但它们的基函数都是预先定义好且固定不变的,缺乏灵活性,而且受到不定原理的制约,其时频分析结果也比较模糊[11].在许多涉及非平稳信号的实际应用中,这些缺陷都是亟待解决的.因为基于数据驱动的信号分解与时频分析方法对输入数据很少有或几乎没有先验假设,所以人们对这些方法产生了极大的兴趣.

这一趋势始于90 年代末,当时Huang 等[3]提出一种递归算法,称为经验模态分解(Empirical mode decomposition,EMD).EMD 通过利用信号极值的递归筛选过程,将输入信号分解为固有的振荡模式,称为本征模态函数(Intrinsic mode function,IMF).一个振荡信号能够被称为IMF 需满足两个条件: 1)在其时间区间内,模态的极值点数目和过零点的数目应当相等或最多相差一个;2)在其时间区间内,分别由信号的局部极大值和极小值确定的上包络和下包络的均值为零.其中,第1 个条件保证了IMF 是一个窄带信号;第2 个条件则从信号局部特征的角度考虑,避免了由于信号波形不对称而引起的瞬时频率波动.图1 是一个典型的IMF 波形示意图,其极值点数目和过零点数目相同,且上包络和下包络的均值为零[1,5].

图1 一个IMF 的波形示意图Fig.1 Waveform diagram of an IMF

由上述讨论和图1 可以看出,IMF 类似于傅里叶变换得到的一个具有固定幅度和固定频率的分量,但是IMF 包含了幅度调制和频率调制的特性,因此更具一般性.

EMD 的筛分过程如图2 所示[1].首先,对信号极值进行插值,并对上包络和下包络进行平均,从而获得信号的局部均值,这些局部均值可视为信号中的低频成分估计;然后,将低频成分从输入信号中迭代分离出来,得到高频(快速振荡)成分.这样就完成了一次筛分[12].重复筛分过程,直到输入信号中的所有主要振荡模态都被提取出来.由于EMD完全是由数据驱动的,避免了线性时频变换方法的很多局限性,如受Heisenberg 不定原理限制而导致的有限时频分辨率和由使用固定基函数而获得的无意义的模态等.因此,EMD 自1998 年问世,便在科学界产生了重大影响,揭开了后小波时代时频分析领域的新篇章,并被应用到了众多工程领域,如语音增强[13]、图像处理[14]、设备诊断学[15]、生物医学[16]、气候学[17]、地球物理学[18]等.

图2 EMD 的筛分过程示意图Fig.2 Schematic diagram of sifting process of EMD

虽然EMD 在处理非线性非平稳信号方面具有较多优点,但是仍然有一些待解决的问题[18-20],主要包括迭代停止准则与样条函数的选择问题、端点效应与模态混叠问题、Hilbert 变换与分量正交性问题等.EMD 的一些改进版本,例如集成EMD (Ensemble empirical mode decomposition,EEMD)[21]、互补EMD (Complementary ensemble empirical mode decomposition,MCEEMD)[22]、中值EEMD(Median ensemble empirical mode decomposition,MEEMD)[23]和中值互补EEMD[24]等,在一定程度上修复了这些问题,但 EMD (包括分解算法、IMF的定义以及上述的EMD 算法)的各种改进版本都是建立在经验基础之上的,目前仍然缺乏坚实的数学理论分析做支撑.对此,Huang 等[18,20]指出,目前EMD 的理论水平犹如小波变换在20 世纪80 年代初的水平,他特别期望能有如Daubechies I.那样,为小波变换奠定坚实理论基础的学者出现,能将基于EMD 的时频分析方法置于坚实的数学理论之上[5].为了克服EMD 理论分析上的困难,部分学者尝试用基于鲁棒约束优化的方法取代EMD 中包络和局部均值估计[25-26],从理论上保证了算法能收敛到全局最优,但这些方法会遇到没有可行解的问题[27].Lin 等[28]提出迭代滤波分解(Iterative filtering decomposition,IFD),试图弥补EMD 数学理论上的不足.IFD 使用特定的移动均值替代EMD筛分过程中的包络均值,实现了自适应滤波,提升了分解稳定性,并严格证明了在一定条件下的算法收敛性[29].Cicone 等[30]进一步为IFD 设计了滤波器长度自适应更新策略,使得滤波器满足迭代滤波收敛的充分条件,为非线性和非平稳信号处理提供了一个完整的局部分析工具箱.

还有一些其他相关方法,例如匹配追踪[31]这种经典的原子分解算法.与常见的正交基函数相比,原子字典具有过完备性,因此可以更加灵活地表征复杂信号.匹配追踪算法将原子库中与当前信号最相关的原子作为当前最优原子,经过多次迭代,可将信号表示成多个最优原子的线性叠加模型.匹配追踪算法的时频分析效果与原子字典复杂程度密切相关,当分析复杂信号时,往往需要复杂原子模型,使得分解结果不稀疏,降低了算法的效率,也会导致时频特征间断问题[32].文献[33-36]受EMD 和压缩感知理论的启发,在本征模态函数组成的字典中,寻找多尺度数据的稀疏表示,将信号分解问题转化为非线性L1优化问题,并提出一种迭代算法递归求解该非线性稀疏优化问题,实现了稀疏时频分析.这些稀疏时频分析方法的分解结果与EMD 较类似,但其性能不受EMD 中停止准则的影响,抗噪声和抑制端点效应的能力也优于EMD.此外,Peng等[37]和Guo 等[38]提出基于算子的信号分解方法,该方法采用零空间追踪的方式,自适应估计算子和参数,将信号中的成分分离.这些方法的性能取决于所构造的算子,常见的算子有奇异局部线性算子[37]和复数微分算子[38]等.

除了EMD 的改进版本以外,学术界还提出一些类似于EMD 迭代筛分过程的新分解算法,来尝试解决这些问题.例如 Smith[39]通过分离调频和调幅成分,提出局部均值分解(Local mean decomposition,LMD)算法.Frei 等[40]通过分段计算旋转分量在每一个局部时间区间上的瞬时频率和瞬时幅值,提出计算复杂度较低的本征时间尺度分解(Intrnsic time-scale decomposition,ITD)算法,可以实现数据的在线分解与实时处理.ITD 的分解过程[1]如图3 所示,其中Xk表示输入信号x(t) 的极值点;虚线L(t) 和虚线H(t) 分别表示基线和旋转分量.需要指出的是,上述分解算法仍然都是基于经验的分解算法,分解结果也会受到端点效应和模态混叠的影响.另外,这些方法都需要准确估计信号的局部极值或均值,但实际数据往往会被噪声干扰,不易准确估计这些极值或均值,这使得这些方法的抗噪声能力较差.

图3 ITD 分解过程示意图Fig.3 Schematic diagram of ITD decomposition process

另一类类似于EMD 的数据驱动信号处理方法包括同步压缩变换(Synchro squeezed transform,SST)[41]和经验小波变换(Empirical wavelet transform,EWT)[12].与EMD 不同,SST 是通过一套合理且方便的数学框架,来提取信号中的模态,它既可以在STFT 域,也可以在小波域中工作[42].SST首先通过频率重整算子[43]锐化信号的STFT 谱图或小波谱图;然后,在模态总数已知条件下,采用脊线提取技术,来估计瞬时频率;最后,通过在相应的脊线附近对重整STFT 或重整小波变换积分,实现模态重构.需要注意的是,SST 要求模态在时频面上各自独立,无法处理时频信息交叉的情况.

EWT 是加州大学洛杉矶分校的Gilles[12]提出的,这种方法融合了经验模态分解和小波变换的优点,通过峰值检测机制对频谱进行分割,并在每个分割区间构建小波滤波器,实现将一个复杂信号分解为一系列具有紧支撑频谱的调频调幅信号.图4是EWT 的模态频谱分割示意图[32],其中竖虚线表示检测到的模态边界频率.EWT 的效果与基于频谱分割所构造自适应的小波滤波器是否精确有很大关系,但对实际信号进行准确的频谱分割是一项颇有挑战性的任务[44].

图4 EWT 的模态频谱分割示意图Fig.4 Schematic diagram of modal spectrum division of EWT

同样来自加州大学洛杉矶分校的Dragomiretskiy等[45]在2013 年提出变分模态分解算法(Variational mode decomposition,VMD).VMD 的目标是将输入信号分解为一系列具有稀疏特性的模态,这里的稀疏特性指的是所有模态都是集中在各自中心频率附近的窄带信号.为了实现这个目标,VMD 构建了一个带约束的变分优化问题,其中目标函数是最小化所有模态的带宽,约束条件是分解得到的模态能够完全重构输入信号.VMD 目标函数的构造分为3 个步骤: 1)对每个模态进行Hilbert 变换,得到其解析信号;2)将解析信号的频谱平移到零中频,得到基带信号;3)利用H1高斯平滑度,估算每个模态的基带信号带宽,并将这些带宽的和最小化作为目标函数.VMD 所建立的优化问题可以通过交替方向乘子法在频域求解,最终得到输入信号中包含的模态和对应的中心频率.

与之前基于经验的信号分解方法不同,VMD的目标函数具有较完备的数学理论基础,它的求解公式表明,VMD 在本质上是一个自适应最优Wiener滤波器组,其中心频率如图5 中竖虚线所示[32].由于其优良的特性,VMD 一经提出,就受到了极大的关注,目前已成功在机械故障诊断[46]、风速预测[47]、疾病诊断[48]、金融数据分析[49]、新型冠状病毒肺炎预测[50]等领域得到应用.

图5 VMD 原理示意图Fig.5 Schematic diagram of VMD principle

VMD 的出现是信号分解领域研究的一个分水岭.在此之前,以EMD 为代表的基于经验的分解方法虽然得到了广泛的关注与应用,但在数学基础理论方面的进展艰难而缓慢,以至于EMD 最开始被提出的时候颇有争议,甚至受到了质疑[51].VMD的出现改变了信号分解方法研究全凭经验的现状,提供了具有坚实数学优化理论的研究视角.但是VMD 也有一些不足,例如VMD 的分解性能与模态数量K和惩罚系数α这两个参数密切相关,并且这两个参数无法自适应获得,需要用户事先指定.

目前,对于VMD 的主要改进工作也集中在如何确定这两个参数的问题上.这些改进VMD 参数依赖性的方法可分为两类: 1)固定惩罚系数α,通过穷举法迭代优化模态数量K.例如Li 等[52]使用近似完全重构的判据,来确定合适的模态数量;Lian等[53]通过判断提取出来的模态特性,来选取合适的模态数量;Cai 等[54]利用傅里叶变换得到的频谱,作为判断模态数量是否合适的标准.上述这类改进方法虽然简单,但由于VMD 的性能受到模态数量K和惩罚系数α的共同调节[55],只考虑模态数量的作用存在一定的风险.2)同时对模态数量K和惩罚系数α进行优化.这类方法的思路是,首先,针对特定的应用场景,构造适应度函数;然后,使用一些启发式智能算法,对参数进行寻优.如人工鱼群算法[56]和蝗虫算法[57]等.尽管这类基于智能寻优的改进方法考虑了模态数量K和惩罚系数α这两个参数,但这些改进方法没有触及VMD 的本质,且有效性受限于特定场景下的适应度函数和优化算法,难以推广到其他场景.为了解决VMD 的参数选取问题,Chen等[58]提出一种自整定变分模态分解算法(Self-tuning variational mode decomposition,SVMD).SVMD 将VMD 的目标函数由同时提取K个模态改进为每次只提取一个模态,且在提取模态的过程中,惩罚系数α不再采用固定值,而是会根据信号特点自适应地对α进行更新.

VMD 实际上只是对窄带信号的分解效果较好,这意味着基于VMD 的方法处理带宽较大时变信号的能力非常有限.为了解决这个问题,2017 年,Chen等[59]提出非线性调频模态分解算法(Nonlinear chirp mode decomposition,NCMD).NCMD 通过解调技术[60],将带宽较宽的时变信号进行频率解调,可以使之转换为窄带信号.NCMD 频率解调过程如图6 所示[59],图中曲线、上方横线和下方横线分别表示原信号、解调信号和基带信号的时频曲线,BW和 B Wmin分别表示原始信号和解调信号的带宽.利用这一特性,NCMD 将信号分解问题转化为解调问题,并建立了带约束的非线性调频模态变分优化模型.

图6 NCMD 频率解调过程示意图Fig.6 Schematic diagram of NCMD frequency demodulation process

与VMD 类似,NCMD 亦采用了交替方向乘子法求解所提出的变分优化模型.求解所得更新公式表明,NCMD 可以被视为一个时频滤波器组,该滤波器组的中心频率就是估计得到的瞬时频率.NCMD 的计算过程可以概括为: 首先,用当前估计得到的瞬时频率对解调的正交信号进行迭代更新;然后,通过反正切解调技术[61],利用正交信号的相位信息进一步更新瞬时频率;重复上述两个步骤,直到解调后的信号具有最窄的频带.与VMD 相比,NCMD 的最大优势在于能够处理非平稳信号,并提供准确的时频信息.虽然NCMD 提出的时间很短,但是由于其在非平稳信号时频分析上的出色性能,目前已被应用于碰撞摩擦故障检测[62]、生理信号监测[63]、高速铁路故障诊断[64]等领域.值得指出的是,虽然NCMD 相比于VMD 有较大优势,但NCMD也具有和VMD 类似的不足,即分解性能依赖先验知识,要求使用者预先指定分解模态数量和惩罚系数.

为了解决这个问题,Chen 等[65]提出NCMD 的改进版本,称为自适应调频模态分解算法(Adaptive chirp mode decomposition,ACMD),并将其成功应用到旋转机械的碰撞摩擦故障检测中[62].ACMD采用一种类似于匹配追踪[31,34]的贪婪算法,对信号模态逐一进行递归估计,直到提取出信号中的所有的模态.虽然ACMD 在自适应时频分析上显示出重大优势,但Chen 等[66-68]发现,ACMD 算法的停止条件对旋转机械故障信号比较有效,但在分解工业过程振荡信号时,频繁出现难以收敛导致过分解问题,因此改进了ACMD 算法的收敛判据,提出快速自适应调频模态分解算法[67],缓解了ACMD 过分解问题,并减少了分解用时.

与NCMD 同年,Chen 等[69-70]提出本征调频模态分解算法(Intrinsic chirp mode decomposition,ICMD).考虑到很多应用场景采集的信号模态的瞬时频率和瞬时幅值具有连续光滑的特性,ICMD 采用傅里叶级数或多项式函数逼近模态的瞬时频率和瞬时幅值,从而构造了一种通用的非线性调频模态参数化模型,能够有效刻画不同调制程度的信号分量.在此模型基础上,ICMD 利用广义参数化时频变换方法[71]估计模态的瞬时频率,进而通过正则化最小二乘法估计模态的瞬时幅值,最终实现信号模态的重构.ICMD 在非平稳信号的分解效果上,与NCMD 难分伯仲,但由于ICMD 并未采用NCMD所使用的交替方向乘子法这种迭代式的优化方法,而是采用简单的最小二乘法,即可得到结果,因此在算法复杂度上,ICMD 具有显著优势.目前,ICMD已被成功应用于信号消噪[72]、雷达信号处理[73]、机械故障诊断[74]等领域.

本节对众多单变量信号分解方法进行了介绍和分析,选择EMD、VMD 和NCMD 作为代表方法进行案例验证,其中EMD 是最经典的信号分解方法,已得到广泛关注和应用,因此选之作为经验性信号分解方法的代表;VMD 是信号分解研究由经验性转向数学优化建模的标志性方法,因此选之作为基于优化的信号分解方法的代表;NCMD 是最近提出的能处理时变信号的分解方法,推动了基于优化的信号分解方法向纵深发展,因此选之作为近年提出的有较大影响力的信号分解方法的代表.本文以单变量信号(1)为例,它包含3 个模态,其中2 个模态是固定频率的正弦函数,频率分别是10 Hz 和30 Hz,最后一个模态具有时变瞬时频率(50+20t)Hz,信号中加入了一定程度的噪声:

EMD、VMD 和NCMD 的分解结果如图7、图8和图9 所示,图中x为原始信号,IMF 为分解所得模态.可以看出,EMD 的第1 个模态提取得较好,但第2 个模态和第3 个模态之间出现了明显的模态混叠效应,这是EMD 的典型问题;VMD 由于其原理所限,难以处理带宽较大的时变信号,因此其分解结果受时变模态的影响,出现了很大误差;而NCMD 由于其算法中采用了解调算子处理时变瞬时频率,能够做到不受时变瞬时频率的影响,正确提取出了3 个模态.这些实验结果与本文中对各个方法的分析是一致的.

图7 EMD 分解结果Fig.7 The decomposition results of EMD

图8 VMD 分解结果Fig.8 The decomposition results of VMD

图9 NCMD 分解结果Fig.9 The decomposition results of NCMD

综上所述,历经20 余年的发展,单变量信号分解研究取得了丰硕成果,已成为现代信号处理领域的研究热点.本文将常见方法根据其作用域分为时域、频域和时频域等多个类别,并在表1 总结了它们的优点和局限性.

表1 常见单变量信号分解方法归类总结Table 1 Classification and summary of common univariate signal decomposition methods

3 多元信号分解

3.1 多元信号分解方法

虽然基于数据驱动的非平稳信号分解和时频分析工作在蓬勃发展,但拓展现有的数据驱动方法来处理非平稳多变量/多通道/多元信号,也引起了学术界和工程界的极大兴趣[75-76].由于传感器和计算机技术的进步,在现代科学和工程应用中,广泛存在对多变量/多通道/多元类型的数据处理方法的需求[77],例如基于多通道脑电[78]或心电信号[79]的分类、多变量信号消噪[80-81]、图像融合[14]等.针对非平稳数据的多变量信号分解与时频分析技术的主要挑战有以下两点: 1)具有模态齐整特性.即不同通道中具有相同或相似的频率的模态出现在同一尺度[82].2)提取出多通道之间的相关信息[83].在信号分解领域的研究中,多元、多通道和多变量表示相同的意思.图10 给出了多元/多通道/多变量信号分解领域的一些术语解释,以便理解本文后续内容.

图10 多元/多通道/多变量信号分解领域术语的图形化解释Fig.10 Graphical interpretation of terms in multivariate signal decomposition

最直接的多变量信号处理方法是使用单变量的方法,逐一、单独分析多变量信号的每一个通道.但文献[75,83-84]指出,对于多变量信号,如果采用单变量信号分解技术逐一处理每个变量,会导致信息泄露和变量之间相关信息缺失等问题,且不满足模态齐整[75]的要求.因此,需要为单变量信号分解方法应用至多变量情形开发特定的拓展方法,以便直接在多变量信号所在的多维空间里处理多变量信号.这也是研究者们开展多元信号分解研究的目的.参照文献[85]的案例,本文以单变量ICMD 和多变量ICMD 处理多变量信号(2)为例:

单变量和多变量信号分解方法处理多变量信号(2)结果见图11 和图12.可以看出,用单变量方法分解多变量信号时,不仅无法满足模态齐整要求,而且分解误差也较大;反之,多变量信号分解方法的性能令人满意,分解所得模态的误差较小,不同通道中具有相似频率的成分也位于同一个尺度,满足模态齐整性质要求.因此,研究多变量信号分解方法处理多变量信号是十分必要和有意义的.

图11 单变量ICMD 分解多变量信号的结果Fig.11 The decomposition results of multivariate signals by the univariate ICMD

图12 多变量ICMD 分解多变量信号的结果Fig.12 The decomposition results of multivariate signals by the multivariate ICMD

多变量信号分解工作是从EMD 应用于复数信号处理开始的.2007 年,Tanaka 等[86]基于复数域的性质,巧妙地利用了原始的EMD,来分解二元时间序列,并把这种方法命名为CEMD (Complex empirical mode decomposition),揭开了数据驱动信号分解与时频分析方法向多变量拓展的序幕.但Tanaka 等[86]的CEMD 仅局限于处理复数信号,不适用于三变量以上情况.此外,CEMD 无法保证复数信号的实部和虚部分解所得的IMF 数量一致,因此这种思路的适用性非常有限.随后,Altaf 等[87]将EMD 中的单变量的“振荡”与双变量的“旋转”概念相对应,认为双变量信号是快速旋转和慢速旋转的叠加,并给出了复数空间中极值的定义,在此基础上,提出旋转复数EMD (Rotation complex empirical mode decomposition,RCEMD).

遗憾的是,Altaf 等[87]的工作只将这个思路应用到了复数信号上,并未充分挖掘其潜能.很快,Rilling等[84]受RCEMD 的启发,将双变量输入信号通过投影映射到单位圆上,并在此基础上,清晰明确地给出了在双变量空间中信号极值、均值和包络的定义,正式提出了双变量EMD (Bivariate empirical mode decomposition,BEMD).BEMD 的分解原理见图13,两种二维包络的均值计算示意图见图14[84],图中曲线上蓝点表示极值点,中心红点表示均值点.图15 进一步对二维局部极值点进行了说明,根据微积分中局部极值点的相关定义,当选择Y轴方向为投影方向时,A点将被视为局部极大值点;同样,B点可视为相应的局部极小值点.因此,以投影方向为基准,可将多变量数据的极值点与投影方向的极值点一一对应,构成了后来基于高维空间投影的多变量信号分解技术的基础.

图15 二维局部极值点示例Fig.15 Example of two-dimensional local extreme points

Ur Rehman 等[88]意识到Rilling 等[84]的双变量拓展思路可以被推广到更多变量的情形.Ur Rehman等[88]通过在高维空间中,建立均匀分布的投影向量集,分别计算输入信号在各个方向上的投影包络线,然后通过计算包络线的均值,定义多变量信号的局部均值函数,并在此基础上,提出三变量EMD (Trivariate empirical mode decomposition,TEMD)[88]、四变量EMD (Quadrivariate empirical mode decomposition,QEMD)[89]和多变量EMD (Multivariate empirical mode decomposition,MEMD)算法[75].Ur Rehman 等[75]提出的MEMD 算法正式开启了多变量信号分解与时频分析的新时代,迅速受到了来自各个领域研究人员的关注.目前,已被应用于脑机接口[90]、设备诊断[91]、因果分析[92]、地球物理[93]、生物医学[78]等诸多领域.然而,MEMD 也继承了原始EMD 的所有局限性,例如对采样频率敏感、噪声鲁棒性差以及EMD 算法的经验性.此外,MEMD随着输入信号数目的增加,计算量会出现指数型增长.为了提高MEMD 的计算效率,Lang 等[94]提出快速MEMD 算法(Fast multivariate empirical mode decomposition,FMEMD),重新定义了多变量IMF 的基本概念,并在此基础上,架构了FMEMD与EMD 方法之间的映射关系,使得FMEMD 计算量与输入信号的数量无关.FMEMD 中一个双变量信号及其投影信号、局部均值和多变量IMF 的示意图见图16、图17、图18 和图19.

图16 双变量信号Fig.16 Bivariate signal

图17 双变量信号的投影信号Fig.17 Projection signal of bivariate signal

图19 多变量IMFFig.19 Multivariate IMF

Lang 等[95-96]创造性地提出了两种多变量形式的ITD 算法,分别称为间接多变量ITD (Indirect multivariate intrinsic time-scale decomposition,IMITD)和直接多变量ITD (Direct multi-variate intrinsic time-scale decomposition,DMITD).IMITD 与FMEMD 的思路相似,即通过定义多变量与单变量运算空间映射关系,架构超定线性方程组,以求解多变量分解结果.需要指出的是,IMITD 使用Halton-Hammersley 采样技术获得高维空间的采样点.与使用等角度采样技术获得的采样点相比,Halton-Hammersley 采样可以使得采样点对高维空间的覆盖更加均匀.

等角度采样和Halton-Hammersley 序列采样在三维球体上的采样点见图20 和图21.可以看出,等角度采样结果在球体两极更为集中,不利于后续的均匀投影操作.IMITD 在局部特征处理上,比MEMD 的效果要好,计算效率也较低,但如果投影方向选择不恰当,会导致IMITD 不能正确提取出期望的基线,从而引起IMITD 分解产物没有物理意义.

图20 等角度采样Fig.20 Uniform angle sampling

图21 Halton-Hammersley 序列采样Fig.21 Halton-Hammersley sequences based sampling

DMITD 是通过恰当定义多变量极值点、多变量基线节点和多变量基线算子等概念,直接拓展ITD 算法至多变量情形,实现ITD 算法在多变量空间的运算操作.与IMITD 相比,DMITD 对投影方向的鲁棒性要好一些,但运算效率低.上述多变量信号分解方法都继承了与之相对应的单变量信号分解的缺陷,例如模态混叠与端点效应.此外,由于这些方法的思路是通过将多变量输入信号投影映射到高维的空间中再分解,因此效果均与投影向量的数目和方向有关.目前还没有明确的依据来确定投影方案,相关研究尚待完善.

多变量SST (Multivariate synchrosqueezed transform,MSST)[97]和多变量EWT (Multivariate empirical wavelet transform,MEWT)[98]放弃了将输入多变量信号投影映射到高维空间的做法,开辟了基于小波变换的多变量信号处理新思路.具体地,MSST 首先对每个信号通道分别应用标准SST 算法;然后,对时频域进行自适应划分,以分离输入数据中的单模态多变量振荡;最后,估计出多变量瞬时频率和幅值,并在此基础上,计算出一个多变量同步压缩变换算子.MSST 可以得到多变量信号清晰的时频谱,对于探索性的数据分析非常有用.但MSST 只给出了时频谱表达,不能重构模态,因此MSST 的应用范围是有限的[97].MEWT 首先采用模态估计过程来获取多变量数据中的最优信号;然后,对其相应的频谱进行分割,以恢复出所有输入信号通道中的所有模态.这种方法存在EWT 固有的缺陷,即需要基于有效的频谱分割来显式构造自适应小波滤波器组[98],这个缺陷对实际物理系统中的信号,难以做到频谱的有效分割.

2019 年,受单变量调制振荡信号推广至多变量形式的启发[99-100],Ur Rehman 等[83]提出多变量VMD算法(Multivariate variational mode decomposition,MVMD).虽然2017 年Wang 等[101]就提出了复数VMD (Complex variational mode decomposition,CVMD),但CVMD 遭遇了与CEMD 一样的问题,即仅局限于处理复数信号,不适用于三变量以上的情况,且无法保证复数信号的实部和虚部分解所得的IMF 数量一致,因此这种思路的适用性非常有限.MVMD 从输入多变量信号中,寻找一组共同的多变量调制振荡,这些振荡在完全重构输入数据所有通道的同时,具有最小的带宽和.MVMD是极具潜力的多变量信号分解方法,这种拓展形式不仅继承了标准VMD 的许多理想性质,还展现出优良的模态齐整特性.MVMD 现已在风机故障诊断[102]、脑电信号检测[103]、基因工程[104]和信号消噪[105]等领域得到成功应用.但 MVMD 直接将模态频谱的重心的估计作为中心频率,因此不适宜处理带宽较大的时变信号.此外,MVMD 也不能直观地提供时频信息.

目前,大多数多变量信号处理方法局限于处理窄带信号,对时变信号的分解能力有限,因此Chen 等[77]受MVMD 的启发,提出多变量非线性调频模态分解算法(Multivariate nonlinear chirp mode decomposition,MNCMD),解决了连续时变条件下多变量信号分解与时频分析问题.MNCMD 现已被成功应用于过程控制系统中的复杂多重厂级振荡根因分析[106-107].MNCMD 虽然在时变多元信号分解与时频表示任务上性能突出,但是它的计算复杂度较高,达到随后,Chen等[85,108]进一步结合傅里叶级数建模,提出多变量本征调频模态分解算法(Multivariate intrinsic chirp mode decomposition,MICMD),能够以 O (N) 复杂度达到与MNCMD相似的时变多元信号分解和时频分析效果.

本节选用MEMD、MVMD 和MNCMD 作为代表性方法,进行案例验证实验.其中,MEMD 使用了高维空间投影的方式进行多变量拓展,也是第1 个多变量信号分解方法,因此选之作为以高维空间投影进行多变量拓展的代表方法;MVMD 使用了多元调制振荡的方式进行多变量拓展,因此选之作为这种多变量拓展方式的代表;MNCMD 是最近提出的能处理时变多元信号的分解方法,因此选之以展示对时变多元信号的处理能力.

本文以多变量信号(3)为例,该信号由四个通道组成:x1通道包含一个时变模态和两个时不变模态,x2通道只包含一个时变模态,x3通道包含两个时不变模态,x4通道包含一个时变模态和一个时不变模态.信号中加入了噪声:

MEMD、MVMD 和MNCMD 的分解结果见图22、图23 和图24,图中第1 行为原始信号,IMF1～IMF5为分解所得模态.可以看出,在MEMD的分解结果中,产生了很多冗余模态,虽然在该结果中,时不变模态提取的效果不错,但时变模态出现了严重的分裂现象;MVMD 虽然避免了MEMD产生冗余模态过多问题,但由于其原理所限,分解结果中的时变模态误差很大,尤其是在时变模态的高频部分,有一部分直接泄露到了第2 个模态中;而MNCMD 能够很好地提取出这些时变和时不变模态,并且表现出模态齐整性.这些实验结果与本文中对各方法优缺点的分析是一致的.

图22 MEMD 的分解结果Fig.22 The decomposition results of MEMD

图23 MVMD 的分解结果Fig.23 The decomposition results of MVMD

图24 MNCMD 的分解结果Fig.24 The decomposition results of MNCMD

综上所述,多元信号分解历经10 余年的发展,涌现了诸多多元信号分解拓展方式,例如高维空间投影和多元调制振荡等,提出了很多各具特点的多元信号分解方法.随着多元信号分解研究的兴起与蓬勃发展,研究者们在此过程中发现了新的问题,给出了新的定义,发展了新的概念,揭示了新的性质,极大丰富了信号分解的研究和应用范围.表2总结了主要多元信号分解方法的优点和局限性.

表2 多元信号分解方法归类总结Table 2 Classification and summary of multivariate signal decomposition methods

3.2 多元信号分解应用

多元信号分解方法已经在很多领域得到了广泛应用.本文以生物医学工程、工业控制系统、机械故障监测、时间序列预测等领域的研究为例,进行说明.

1)在生物医学工程领域,多元信号分解方法的典型应用是心电信号ECG (Electroen cephalo gram)和脑电信号EEG 的处理.ECG 信号是反映心脏基本功能和病理信息的重要参考值.ECG 信号的采集需要在胸部的右上、左上、左下和右下侧连接多个导联,以观察心电图的变化,这就构成了多个信号通道,形成了多变量信号.例如文献[109-110]使用MEMD 处理ECG 信号,移除了基线漂移的影响,提高了心血管疾病诊断的准确性.在进行脑电数据采集时,通常会根据需要选取不同数量的导联或电极点,在不同的脑区采集EEG 信号,这就形成了多变量EEG 信号.在脑机接口研究中,各种复杂脑电信号的非稳定性和多通道性一直阻碍着常规基于基函数信号处理方法的应用及多通道同步分析.文献[78,111]利用MEMD 提取多通道脑电信号中的特征,提升了癫痫疾病诊断效果.文献[112]利用MVMD 分析多通道脑电信号的时频特征,并在人类情感识别任务中取得了优异的性能.

2)在工业控制系统中,多变量信号分解方法在性能评估上也得到了广泛应用.由于控制系统的大规模、高集成、强耦合的特点[113-114],不同设备或单元之间会相互影响,极易在系统的多个部位表现出相似的故障特征.典型的例子是过程控制系统中的厂级振荡.控制系统规模庞大、机理复杂,使得厂级振荡表现出非线性、非平稳、多模态、强噪声等特征,这影响了厂级振荡的检测和诊断效果,严重时甚至威胁整个系统的稳定性和安全性[115].文献[95,116-117]分别使用了MEMD、MITD 和MNCMD,对不同的工业控制系统厂级振荡数据进行分析,一致认为,多元信号分解方法是目前厂级振荡检测最主流和最有效的手段之一,有利于提升控制系统性能水平.

3)近年来,多变量信号分解方法在机械故障监测领域得到了较多关注.例如文献[118]将MEMD与Teager 能量谱结合,获取信号的故障特征信息,实现微小故障特征的提取,并在核主泵轴承外圈早期故障检测实际任务中的表现优于基于EMD 的方法.文献[119]在改进MVMD 自适应性的基础上,实现了一种具有突出优势的轴承故障诊断方法.文献[120]利用MNCMD 能处理时变多元信号的特点,实现了对转子轴承系统在非平稳过程(如启动和关闭)中,瞬时振动状态的实时分析与监测.

4)多变量信号分解可以用于提升多变量时间序列预测性能.因为现实时间序列数据往往具有高度非线性和非平稳特征,多变量信号分解方法不仅能有效处理非线性和非平稳因素,还能充分保留变量之间的相关信息,以提升多变量时间序列预测性能.文献[121]提出一种基于MEMD 和支持向量回归的混合预测模型,它的新颖之处主要在于MEMD的应用,使得多元数据分解能够有效提取不同相关变量之间的固有信息.在多个数据集上的实验结果表明,基于MEMD 的混合模型是一种很有前景的电力峰值负荷预测方法.除了电力负荷预测之外,多元信号分解还广泛应用于各种行业的时间序列预测,例如经济领域的股价指数预测[122]和原油价格预测[123]、物理学中的太阳辐射预测[124-125]、生态环境中的空气质量预测[126]等.

综上所述,多元信号分解方法并不局限应用于某一类场景,可以与多个学科交叉,为各行各业赋能,从而产生研究价值和应用价值.目前,多元信号分解方法的应用仍然在蓬勃发展中,表3 归纳总结了几种常见多元信号分解方法的适用场景.

表3 常见多元信号分解方法的适用场景Table 3 Applicable scenarios of common multivariate signal decomposition methods

4 总结与展望

非线性和非平稳是现实信号的普遍规律,时频联合分析技术正是应现实科学研究和工程应用需求而产生和发展的.瞬时频率是研究非平稳和非线性信号最有力的工具[127],但对于实际的非线性和非平稳信号,由于Bedrosian 定理[128]的限制,难以直接求解瞬时频率[129].一个自然的做法是将所研究的信号分解为一个个单分量的信号,每个单分量信号只包含一种振荡模态,这就是基于信号分解的时频分析研究最初的想法.本文从单变量和多变量两个方面,梳理了信号分解领域的发展历程与研究现状,比较分析了主流信号分解方法的优缺点,在此基础上,本文认为未来可以从以下五个方向进一步探索:

1)间歇信号分解理论与技术.虽然Xie 等[130]结合K均值聚类方法,提升了ITD 处理间歇信号的能力,但对于大多数方法,特别是VMD、NCMD、ICMD 等基于优化建模的信号分解方法,均要求瞬时频率和瞬时幅值都是连续光滑函数,这意味着间歇信号并不满足这些信号分解方法的假设条件,因此现有的信号分解方法处理间歇信号的能力极为有限.进一步地,目前多变量间歇信号分解理论与技术是一个开放性问题.对于MEMD、MITD、MVMD、MNCMD 和MICMD 等多元信号分解方法中出现的模态齐整特性,在多变量间歇信号中该如何定义和体现,仍然是一个尚未讨论的问题.因此,研究复杂间歇信号分解理论和开发复杂间歇信号时频分析技术,对丰富信号分解与时频分析领域研究具有重要意义,也将极大拓展信号分解理论方法应用范围[131].

2)自适应或无参化信号分解理论与技术.虽然SVMD 等方法在一定程度上为特定的信号分解技术提供了参数整定或自适应更新方案,但是这些方案依然会有一些超参数或阈值需要提前指定.需要注意的是,很多改进方案(如EEMD、CEEMD (Complementary ensemble empirical mode decomposition)等),虽然提升了EMD 抗模态混叠和抗端点效应的能力,但也引入了新的超参数;另外,还有一些多变量拓展方式也引入了超参数,例如基于高维空间投影思路的多元信号分解方法就会引入投影向量方向和数量这两个超参数.这些超参数会影响到信号分解方法在处理实际信号时的实用性和便捷性.本文认为至少可以从以下两个角度来缓解这个问题: a)将已有分解方法与参数寻优技术结合,为特定场景下的应用需求提供定制化的参数整定方案[53,132-134];b)研究新的参数自适应调整的或无参化的信号分解方法,这具有相当的挑战性.

3)新型多变量拓展技术.多变量/多元/多通道信号分解与时频分析研究方兴未艾.对于现有的多变量信号分解与时频分析技术,目前主要是基于高维空间投影[75,94-95]和多变量调制振荡[97,100]两种多元拓展思路.基于高维空间投影方法会带来投影向量方向和投影向量数目敏感性的问题,基于多变量调制振荡方法则有时会出现零频分量波动的问题.如何解决现有多变量拓展技术的不足和开发新的多变量拓展方案,是一件富有挑战而又令人期待的任务.

4)信号分解方法的理论完备性.对于经典的基于经验性的EMD、LMD、ITD 等方法,虽然应用领域广泛、使用效果较好,但是其缺乏严格的数学理论支持,亟待能有如Daubechies I.那样能为小波变换奠定坚实理论基础的学者出现,能将基于EMD的时频分析方法置于坚实的数学理论之上[1].对于新兴的基于优化的VMD、NCMD 等方法,虽然具备优化模型的理论支持,但其算法收敛性和结果唯一性的证明仍然是一个悬而未决的问题,这使得在使用这类方法时仍然把握不足.因此,进一步完备信号分解理论基础,明确算法收敛条件,对该领域的长期化与实用化发展具有重大意义.

5)新领域和交叉领域的探索融合.信号分解属于基础的信号处理技术,是工程学科的基础学科,是一门“使能”学科,它既可以探索新领域,发现新需求,又可以与传统学科领域结合起来,为其他学科的传统需求赋能.例如引力波的时频分析研究[135],机械、工业控制等领域的故障检测与诊断研究[136-137];复杂系统因果分析研究[106,117]、脑机接口与生理疾病特征识别[110]、新型冠状病毒肺炎相关问题分析与预测[50,138]、云计算集群中的复杂时序信号周期检测与分解应用[139-141].信号分解与Transformer[142]结合提出的新型神经网络,使单变量时序预测和多变量时序预测误差分别降低了14.8%和22.6%[143-144].诸如此类的旧问题与新需求广泛存在,蓬勃发展的信号分解理论方法为这些问题和需求的解决提供了有效的新途径,这些新领域和交叉领域的发展,也启发着信号分解技术新的突破方向.