基于解析模态分解的时变与弱非线性结构密集模态参数识别

王佐才，任伟新，邢云斐

（合肥工业大学 土木与水利工程学院，合肥 230009）

结构健康监测和安全评估已成为国内外众多学者致力研究的一个重要课题。而其中关键问题之一是对结构的参数，特别是时变与非线性系统的时变参数进行有效的识别。因此，近10年来，基于希尔伯特变换在时变和非线性系统识别和损伤检测领域中得到了广泛的关注［1－8］。如：Feldman［5］提出的 Hilbert Vibration Decomposition（HVD）方法，并利用该方法对非线性结构的时变参数进行了较为深入的研究。然而，HVD分解对噪声非常敏感，微小的噪声对分解的结果也会产生较大的影响，另一方面，如果一个信号具有多个同样能量量级或者密集频率的信号分量，HVD也很难将其分解。

为了实现对信号瞬时频率的识别，黄锷等［6－8］提出了经验模态分解（EmpiricalModeDecomposition，EMD）和希尔伯特 －黄变换（Hilbert-HuangTransform，HHT）。EMD能够将复杂信号分解为有限个本征模函数（IntrinsicModeFunction，IMF），对每一个IMF进行希尔伯特变换和谱分析，可以实现复杂信号的瞬时特征提取，是一种有效的非平稳信号分析工具［9－11］。

虽然HHT在非平稳和非线性信号的特征提取方面取得了广泛的应用，但在许多实际应用中也遇到了各种困难［12］。例如：HHT无法分离密集的模态响应，特别是具有模态频率叠混（overlapping）的信号；对于时变的具有模态叠混的结构响应，结构的模态响应也往往存在于多个邻近的分解信号中，需要进行进一步的重组。为了解决这些难题，国内外的学者也做了相关的研究工作。例如，Chen［13］研究了HHT方法在密集模态结构中模态参数的识别方法，在EMD分解过程中采用间歇检查（Intermittencycheck）来分离密集模态。Yang［14－15］在进行 EMD分解前，利用带通滤波器对信号进行滤波预处理，从而把结构的密集模态分离出来。Wang等［16］提出了波组分解的方法，通过把信号的频率从高频转化到低频，增大相邻频率的比值，从而使得EMD有可能把低频的信号分解出来。程远胜等［17］提出将HHT与数学规划方法相结合的方法，用于时变多自由度系统的参数识别。黄天立等［18］提出采用波组信号前处理和正交化经验模式分解的方法予以改进，并将此方法称为改进的HHT。然而，对于具有密集模态和时变频率的信号分解，到目前为止，几乎还没有一种有效的方法能够很好的解决这些问题。另一方面，目前的方法往往简单地把识别出的响应信号的瞬时频率等同于结构本身的瞬时频率，但是，对于时变或者非线性结构系统在环境激励或者地震荷载作用下的响应信号的瞬时频率并不等同于结构本身的瞬时频率，如何利用信号的瞬时频率来推导时变与非线性结构本身的瞬时频率也并没有很好的解决。

Chen等［12］进一步结合希尔伯特变换提出了解析模态分解方法（AnalyticalModeDecompositionAMD）。其本质是利用Hilbert变换把每一具有特定频率成分的信号解析地分解出来。对于具有多个密集时不变频率信号分量的复杂信号，AMD通过构造一对具有相同特定频率的正交函数，并利用这对正交函数与原复杂信号的乘积的Hilbert变换，把任意在频率小于正交函数频率的信号解析地分解出来，并最终实现密集模态的参数识别。Chen等［19－20］进一步提到，通过选取一个时变的二分截止频率代替时不变的截止频率，AMD可以进一步应用于非平稳信号的分解。

因此，本文为了解决具有密集模态的时变与非线性多自由度体系的时变参数识别问题，首先推导了单自由度与多自由度体系在自由振动和受迫振动下模态响应信号的瞬时频率与结构本身瞬时频率之间的关系。本文提出了把解析模式分解方法扩展到时变与非线性结构的模态分解。该方法通过小波变换选取时变的二分截止频率，对结构的时变模态响应进行分离，从而实现多自由结构时变参数识别。

1 时变与非线性结构系统瞬时频率识别

1.1 单自由度系统（SDOF）

时变或非线性单自由度系统结构的自由振动运动方程［21］可以简单的表述为：

对于非线性结构，非线性恢复力随时间变化的函数可以转换成（t）x（t）的形式［21］，（t）是结构瞬时频率，x（t）是结构位移响应。同样，非线性阻尼力也可以转换成2h0（t）x·（t）的形式，h0（t）为快时变瞬时阻尼，x·（t）为速度。

对式（1）进行希尔伯特变换，可得

式中：h0L和h0H分别为时变非线性阻尼h0（t）的低通（慢变部分）和高通（快变部分）滤波值；和分别为时变固有频率（t）的低通（慢变部分）和高通（快变部分）滤波值。

对于单自由度系统的自由振动，位移响应x（t）可以测得，位移的解析信号Z（t）可以表述为

Z（t）＝x（t）＋jH［x（t）］＝A（t）ej∫ω（t）dt（3）式中，A（t）和ω（t）分别是解析信号的瞬时振幅和瞬时频率。对解析信号求一次和二次导数，可以得到

因此，速度和加速度以及它们的希尔伯特变换均可以表述为关于位移的函数

将等式（5）和式（6）代入式（1）和式（2）中，可得

由于任何信号x和它的希尔伯特变换H［x］的相位差为90°，所以它们不可能同时为零。因此，式（7）和式（8）中系数矩阵的行列式必须等于零，从而导出：

对于单自由度时变或非线性系统的受迫振动，其运动方程可以表述为：

与自由振动类似，如果瞬时阻尼h0（t）和瞬时固有频率（t）是关于时间的慢时变函数，那么类似的可以推导出瞬时阻尼 h0（t）和固有频率（t）：

工程实际中，由于响应信号包络线的导数值都比固有频率小得多，其影响可以忽略不计，即0。因此，响应信号的瞬时频率可以写成：

环境激励或者地震作用f可以假设为具有零均值的时程函数，所以式（13）中的后一部分也是零均值快变的时变函数。如果（t）是弱非线性系统的慢变时变函数，通过AMD方法［12］过滤掉信号瞬时频率 ω2（t）快变时变的部分，就可以获得结构系统的固有频率（t）。

1.2 多自由度系统（MDOF）

对于有n个自由度的时变或弱非线性系统，其运动方程可以表述为

式中：M（t），C（t）和 K（t）分别是时变质量、阻尼和刚度矩阵：f（t）是外荷载。在模态坐标中，式（14）可以转换为

式中：是第i阶模态质量；i是第i阶振型向量；（t）是第i阶模态响应的固有频率。模态响应的信号瞬时频率（t）可以写成

此外，环境激励或者地震作用f可假设为具有零均值的时程函数，同样的，式（16）的后一部分是零均值快变的时变函数。通过AMD方法［12］过滤掉模态响应qi（t）的瞬时频率（t）快变时变的部分，就可以获得结构系统的第i阶固有频率（t）。

对于多自由度体系，所测得的第I个自由度的位移xl（t）为模态响应 qi（t）函数

式中：li为第I个自由度的第i阶振型向量。因此，从xl（t）中分解而得的第i阶模态响应x（i）l（t）可以写成：

对于有n个自由度的时变或弱非线性系统，li为慢变时变函数且未知，因此，无法获得第i阶模态响应qi（t），而第 i阶响应（t）＝liqi则可以通过下文将要介绍的AMD分解方法获得，其相应的解析信号可表示为

从式（19）可以看出，分离出来的模态响应信号（t）的瞬时频率等于模态响应qi（t）的瞬时频率。因此，过滤掉AMD分解得到模态响应（t）的瞬时频率的快变时变的部分，就可以获得结构系统的第i阶固有频率（t）。

2 基于AMD方法的时变非线性结构响应信号分解

Chen等［12］首先提出了具有密集模态的时不变信号的解析模式分解定理，并给出了相应的证明。解析模式分解可以描述如下：

对于任意由n个信号分量（t）（i＝1，2，…，n）组成的原信号x（t），如果它的每一分量的频率 ω1，ω2，…，ωn（ωi＞0；i＝1，2，…，n）满足：（｜ω1｜＜ωb1），（ωb1＜｜ω2｜＜ωb2），…，（ωb（n－2）＜｜ωn－1｜＜ωb（n－1）和（ωb（n－1）＜｜ωn｜），其中，ωbi∈（ωi，ωi＋1）（i＝1，2，…，n－1）为n－1个二分截止频率，那么它的每一信号分量可以解析地给出：

式中，H［·］表示希尔伯特变换运算。

对于具有时变频率的非平稳信号，本文将解析模式分解定理拓展成如下：

对于任意由n个信号分量（t）（i＝1，2，…，n）组成的原信号x（t），如果它的每一分量的时变频率ω1（t），ω2（t），…，ωn（t）满足：｜ω1（t）｜＜ωb1（t），ωb1（t）＜｜ω2（t）｜＜ωb2（t），…，ωb（n－2）（t）＜｜ωn－1（t）｜＜ωb（n－1）（t），和 ωb（n－1）（t）＜｜ωn（t）｜。其中，ωbi（t）∈（ωi（t），ωi＋1（t））（i＝1，2，…，n－1）是选取的时变截止频率。那么它的每一信号分量可以解析地给出：

式中：dτ为时变截止频率的积分，即相位角。

根据上述AMD分解的表述，可以设计如图1所示的自适应低通滤波器。AMD的本质是利用Hilbert变换把每一具有特定频率成分的信号解析的分解出来。对于多个时变密集频率信号叠加的复杂信号，AMD定理通过构造一对具有相同特定时变频率的正交函数，并利用这对时变正交函数与原复杂信号的乘积的Hilbert变换，把任意在频率时间平面内低于正交函数时变频率的信号解析地分解出来。

图1 基于AMD分解的时频低通滤波器框图Fig．1 Block diagram of a time-frequency lowpass filter with AMD method

3 自由振动下时变与非线性系统模态参数识别

3.1 自由振动下杜芬系统的瞬时频率识别

本文首先用硬化弹簧和线性阻尼组成的杜芬系统进行数值模拟，其杜芬方程［21］为

系统初始位移为10 cm，初始速度为0的自由振动信号可由四阶龙格库塔法计算获得。本算例取时间间隔为0．1 s，其计算得到的自由振动响应和傅里叶谱如图2所示。

图2 非线性杜芬系统的自由振动响应和傅里叶谱Fig．2 Free response and its Fourier spectrum of nonlinear Duffing system

对于单自由的体系，通过AMD方法过滤掉信号瞬时频率的快变时变的部分，就可以获得结构系统的固有频率。响应信号的瞬时频率，阻尼系数以及利用本文方法识别出来的结构瞬时频率、阻尼如图3所示。从图3中可以看出，在振幅较大的时间范围内，结构非线性硬化程度比较大，结构的瞬时频率也比较大。自由振动开始时，瞬时频率迅速下降，最后逐渐趋于相应的线性系统的频率值（0．16 Hz）。因此，响应信号的瞬时频率的慢变部分可以准确地反映弱非线性系统的频率变化。响应信号的瞬时阻尼系数的慢变部分为水平直线（0．025），说明识别出来的阻尼系数为常数0．025，这与杜芬系统的阻尼系数也是一致的。

图3 识别出的瞬时频率和阻尼系数Fig．3 Identified instantaneous frequency and damping coefficient

3.2 自由振动下索的时变频率识别试验

为了验证本文方法在时变参数识别上的有效性，本文利用文献［22］的一时变拉索试验数据进行了验证。试验索为一根75的钢绞线，一端用反力架锚固，另一端固定在MTS加载系统上。两锚固点间的索长为4．55 m。将加速度传感器竖向安装在索的中部，其采样频率600 Hz，试验装置如图4所示。拉索由锤子敲击产生自由振动，本文分别对拉力线性变化和正弦变化两种情况的时变频率进行了识别。

图4 拉索试验装置Fig．4 Cable test setup

拉索拉力与索的频率关系可以通过固定索的拉力，测得索的自由振动响应，并利用峰值法识别出结构在固定拉力下的固有频率，并把固定拉力下测得的索的固有频率作为准确值。索的拉力与固有频率的关系如图5所示。

图5 同定值拉力作用下识别出的一阶固有频率Fig．5 Identified first natural frequency with fixed tension forces

3．2．1 拉力线性变化时索的瞬时频率识别

试验时索的拉力从20 kN开始以1．67kN／s的速率线性增加，索力变化过程中同时采集索的冲击加速度响应，共采集6 s，响应信号如图6所示。由于用锤子产生自由振动的方式无避免噪声效应的产生，如果这个噪声激励是零均值的，那么所测模态响应中慢变时变部分的瞬时频率可以看作是索的时变频率。因此，本文选用二分截止频率为20 Hz，利用AMD把第一阶模态响应分解出来，识别出的响应信号的瞬时频率和结构的瞬时频率如图7所示。从图中可以看出，测试获得模态响应信号中的瞬时频率等于系统的慢变时变频率加上环境振动的零均值快变的时变频率。因而模态响应信号中的瞬时频率的慢变时变频率部分可以有效地识别出索的时变频率。

图6 拉力线性变化时测得索中点处的自由响应Fig．6 Measured free response at middle of cable with linear varying tension force

图7 拉力线性变化时识别出的瞬时频率Fig．7 Identified instantaneous frequency with linear varying tension force

图8 拉力正弦变化时测得索中点处的自由响应Fig．8 Measured free response at middle of cable with sinusoidal tension force

3．2．2 拉力正弦变化时索的瞬时频率识别

此试验中，拉力按的方式随时间正弦变化。同样，索力变化过程中同时采集索的冲击加速度响应，共采集6 s。响应信号如图8所示。通过本文方法识别出的瞬时频率如图9所示，从图中可以看出，模态响应瞬时频率慢时变的部分与系统瞬时频率是较为吻合的，可以有效的识别出索的时变频率。

图9 拉力正弦变化时识别出的瞬时频率Fig．9 Identified instantaneous frequency with sinusoidal tension force

4 受迫振动下具有密集模态的时变系统参数识别

对一具有密集模态的时变两层框架系统（图10）进行数值模拟。下层和上层的集中质量分别为m1＝2．63×105kg和 m2＝1．75×103kg，结构的第一和第二阶阻尼比为2%。结构的第一阶和第二阶固有时变频率为。对应的结构下层和上层刚度系数为：

分别对结构在1940 El Centro地震激励下和具有零均值和0．1 g（g为重力加速度）方差的高斯白噪声激励下的响应进行数值模拟，并把第一层模拟的位移响应x1作为实测的时程响应。为了利用AMD对模态响应进行分离，首先选用中心频率为5 Hz，带宽为8的复Morlet小波函数，对模拟的信号进行小波变换，获得小波量图，并粗略利用模最大的方法提取瞬时频率脊线（如图11），从而通过取小波脊线的平均值作为AMD分解的时变截止频率。利用AMD分解获得每一阶的模态响应如图12所示。利用本文方法识别的模态响应信号的瞬时频率和结构瞬时频率如图13所示。

该方法首先利用小波变换近似将密集模态的时变频率脊线提取出来，再根据近似的时变脊线选择每一信号分量的时变截止频率，利用AMD分解实现每一信号的分离，而不再需要对信号进行重建。这样的优点是，根据Heisenberg-Gabor不确定性定理，小波分析不能同时在时域与频域内达到较高的精度［23］，那么对具有密集模态的时变响应，不仅需要较高的频率精度，实现模态的分离，同时也需要较高的时间精度，从而可以准确地得到信号的瞬时特征。结合AMD分解的优势，先通过选择合适的具有较高频率精度的小波参数，先选择失去时间精度的瞬时频率，利用AMD实现每一信号的分离和本文提取的时变参数识别方法，就可以较精确地得到同时具有较高时间与频率精度的瞬时频率。

从图11中也可以看出，简单的通过模拟最大提取的小波脊线并不能有效的识别出结构的瞬时频率。而如图13所示，结合AMD与本文提出的方法识别出来的结构瞬时频率与理论值十分吻合，从而有效地识别出了具有密集模态的多自由的瞬时频率。

图10 两层框架结构Fig．10 Two-story shear building

图11 基于小波变换的结构瞬时频率参数识别Fig．11 Instantaneous frequencies by wavelet transform method

图12 AMD分解的模态响应Fig．12 Decomposed modal responses by AMD

图13 基于AMD分解方法的瞬时频率识别Fig．13 Identified instantaneous frequencies by AMD

5 结 论

本文推导了单自由度与多自由度体系在自由振动和受迫振动下模态响应信号的瞬时频率与结构本身的瞬时频率的关系。对于具有密集模态的时变与非线性的多自由度体系，通过小波变换选取二分时变截止频率，并利用本文提出的扩展解析模式分解方法对时变与非线性结构的模态响应进行分解，从而实现对多自由度结构时变参数的识别。通过数值模拟和试验验证，可以得到如下结论：

（1）对于时变的线性结构和弱非线性结构，模态响应信号的瞬时频率缓慢变化的部分与结构系统的瞬时频率近似相等。因此可以利用模态响应信号中的瞬时频率的慢变部分有效地识别出结构瞬时频率。

（2）对于具有密集模态的多自由度结构体系，利用小波变换近似选择二分时变截止频率，并利用扩展的解析模式分解方法可以有效地对结构的模态响应信号进行分离，从而实现对具有密集模态的多自由度体系瞬时频率的识别。

（3）数值模拟和试验结果表明，本文方法对时变线性和非线性系统结构的瞬时模态参数识别具有较高的时间和频率精度。

［1］Liu B，Riemenschneider S，Xu Y．Gearbox fault diagnosis using empirical mode decomposition and Hilbert spectrum［J］．Mechanical Systems and Signal Processing，2006，20（3）：718－734．

［2］Chen H G，Yan Y J，Jiang J S．Vibration-based damaged detection in composite wingbox structures by HHT［J］．Mechanical Systems and Signal Processing，2007，21（1）：307－321．

［3］Feldman M． Non-linear system vibration analysis using Hilbert transform-I：free vibration analysis method［J］．Mechanical Systems and Signal Processing，1994，8（2）：119－127．

［4］Feldman M．Non-linear free-vibration identification via the Hilbert transform［J］．Journal of Sound and Vibration，1997，208（3）：475－489．

［5］Feldman M． Time-varying vibration decomposition and analysis based on the Hilbert transform［J］．Journal of Sound and Vibration，2006，295（3－5）：518－530．

［6］Huang N E，Shen Z，Long S R．A new view of nonlinear water waves：the Hilbert spectrum［J］．Annual Review of Fluid Mechanics，1999，31：417－457．

［7］Huang N E，Shen Z，Long S R，et al．The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis［C］．Proceedings of the Royal Society of London-Series A，1998，454：903－995．

［8］Huang N E，Wu M C，Long S，et al．Confidence limit for the empirical mode decomposition and Hilbert spectral analysis［C］．Proceeding of the Royal Society of London，Series A，2003，459：2317－2345．

［9］罗奇峰，石春香．Hilbert-Huang变换理论及其计算中的问题 ［J］．同济大学学报（自然科学版），2003，31（6）：637－640．LUO Qi-feng，SHI Chun-xiang．Hilbert transform and several problems in its calculation method［J］．Journal of Tongji University，2003，31（6）：637－640．

［10］陈隽，徐幼麟．HHT方法在结构模态参数识别中的应用［J］．振动工程学报，2003，16（3）：383－388．CHEN Jun，XU You-lin．Application of HHT for modal parameter identification to civil structures［J］．Journal of Vibration Engineering，2003，16（3）：383－388．

［11］黄天立，楼梦麟．基于HHT的非线性结构系统识别研究［J］．地震工程与工程振动，2006，26（3）：80－83．HUANG Tian-li，LOU Meng-lin．System Identification of nonlinear structures based on HHT［J］． Earthquake Engineering and Engieering Vibration，26（3）：80－83．

［12］Chen G D，Wang Z C．A signal decomposition theorem with Hilbert transform and its application to narrowband time series with closely spaced frequency components［J］．Mechanical Systems and Signal Processing，2012，28：258－279．

［13］ Chen J，Xu Y L．Identification of modal damping ratios of structures with closely spaced modal frequencies［J］．Structural Engineering and Mechanics，2002，14（4）：417－434．

［14］Yang J N，Lei Y，Pan S，et al．System identification of linear structures based on Hilbert-Huang spectral analysis，part I：normal modes［J］．Earthquake Engineering and Structural Dynamics，2003，32（9）：1443－1467．

［15］Yang JN，Lei Y，Pan S，et al．System identification of linear structures based on Hilbert-Huang spectral analysis，part II：complex modes［J］．Earthquake Engineering and Structural Dynamics，2003，32（10）：1533－1554．

［16］Wang W．Decomposition of wave groups with EMD method in The Hilbert Transform in Engineering［M］．Taylor＆Francis Group，LLC，FL，U．S．A．2005，267－280．

［17］程远胜，熊飞，刘均．基于HHT方法的时变多自由度系统的参数识别［J］．华中科技大学学报（自然科学版），2007，35（5）：41－43．CHENG Yuan-sheng，XIONG Fei，LIU Jun．Identifying the parameters of multi-freedom degree systems with time varying using Hilbert-Huang transform ［J］．Journal of Huazhong University of Science and Technology，2007，35（5）：41－43．

［18］黄天立，邱发强，楼梦麟．基于改进HHT方法的密集模态结构参数识别［J］．中南大学学报（自然科学版），2011，42（7）：2054－2062．HUANG Tian-li，QIU Fa-qiang，LOU Meng-lin．Application of an improved HHT method for modal parameters identification of structures with closely spaced modes［J］．Journal of Central South University，2011，42（7）：2054－2062．

［19］Chen GD，Wang Z C．Response to the letter to editor by Dr．M．Feldman entitled a signal decomposition or lowpass filtering with Hilbert transform［J］．Mechanical Systems and Signal Processing，2011，25（8）：3204．

［20］Feldman M．A signal decomposition or lowpass filtering with Hilbert transform［J］． Mechanical Systems and Signal Processing，2011，25（8）：3205－3208．

［21］Feldman M．Hilbert transform applications in mechanical vibration［M］．Wiley，2011：292．

［22］王超，任伟新，黄天立．基于复小波变换的结构瞬时频率识别［J］．振动工程学报，2009，22（5）：492－496．WANG Chao，REN Wei-xin，HUANG Tian-li．Instantaneous frequency identification of a structure based on complex wavelet transform ［J］．Journal of Vibration Engineering，2009，22（5）：492－496．

［23］Flandrin P． Time-frequency／time-scale analysis［M］．Academic press，San Diego，CA，1999．