CN111487318A - 一种时变结构瞬时频率提取方法 - Google Patents

Abstract

本发明涉及一种时变结构瞬时频率提取方法，对于桥梁结构，首先将广义解调定理与解析模态分解理论相结合，将具有密集模态甚至叠混的多分量信号分解为多个单分量信号，再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，从而实现调幅调频信号到纯调频信号的转变，最后通过小波系数模局部极大值方法精确提取瞬时频率。该方法有利于对桥梁时变结构的瞬时频率进行提取。

Description

一种时变结构瞬时频率提取方法

技术领域

本发明属于结构健康监测和安全评估领域，具体涉及一种时变结构瞬时频率提取方法。

背景技术

实际桥梁结构在环境侵蚀、材料老化、工作荷载、疲劳效应等多重因素作用下不可避免地出现损伤并逐渐累积，在极端情况下甚至有可能引发突发性的灾难事故。因此，为保障桥梁结构的安全性、完整性、适用性和耐久性，减少人员伤亡和财产损失，对在役桥梁结构采取有效的手段进行结构健康监测和安全评估显得十分必要和迫切，同时也具有较大的社会价值和经济效益。

作为桥梁结构健康监测和安全评估领域的一项核心技术，参数识别属于系统识别范畴，一般包括以质量、刚度和阻尼为特征参数的物理参数识别和以模态频率、模态振型和模态阻尼为特征参数的模态参数识别。时不变结构模态参数识别方法包括频域法和时域法，它们假定输入激励为白噪声、结构响应为平稳信号，因此无法反映模态参数的时变特征，同时也无法追踪结构的损伤演化趋势。截至目前，现有模态参数识别方法主要集中在时不变结构系统，针对环境激励下的时变结构开展密集模态参数和瞬时特征参数识别研究尚不多见。

具有刚度或质量突变的结构系统的动力特征参数必将随时间发生变化，其模态频率有可能呈现密集或者叠混特性。因此，环境激励下的在役桥梁结构，其响应信号不仅是非平稳的，而且极有可能是模态密集分布和叠混的。显然，从非平稳信号处理的角度识别环境激励下桥梁结构的时变密集模态参数更符合实际情况，对深入理解结构损伤诊断、有限元模型修正、振动控制和安全评估具有重要的理论意义和工程应用价值。

发明内容

本发明的目的在于提供一种时变结构瞬时频率提取方法，该方法有利于对桥梁时变结构的瞬时频率进行提取。

为实现上述目的，本发明采用的技术方案是：一种时变结构瞬时频率提取方法，对于时变钢桥或组合桥梁结构，首先将广义解调定理与解析模态分解理论相结合，将具有密集模态甚至叠混的多分量信号分解为多个单分量信号，再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，从而实现调幅调频信号到纯调频信号的转变，最后通过小波系数模局部极大值方法精确提取瞬时频率。

进一步地，该方法具体包括如下步骤：

步骤1)对于瞬时相位为φ(t)的渐进信号x(t)，其相应的解析信号表达式x_A(t)和瞬时频率f(t)分别如式(1)和(2)所示：

x_A(t)＝x(t)+jH[x(t)]＝x(t)e^j[φ(t)] (1)

式中，j为虚数单位，H[·]表示希尔伯特变换算子；通过解析信号x_A(t)来表示信号x(t)的广义傅里叶变换，即信号x(t)的调制，如式(3)所示：

式中，X_G(ω)为调制后的信号，φ(t)为与时间相关的实相位信号；为了将信号的时变瞬时频率调制到固定频率f₀，构造φ(t)，如式(4)所示：

式中，

为原多分量信号各阶瞬时频率估计值；

步骤2)对调制后的信号X_G(ω)进行解析模态分解，选取该信号的实部X_G(t)，假设X_G(t)是由n个分量x_i(t)组成的多分量信号：

X_G(t)＝x₁(t)+x₂(t)+…+x_i(t)+…+x_n(t)i＝1，2，…，n (5)

若分量信号的频率ω₁，ω₂，…，ω_n满足：|ω₁|＜ω_b1，ω_b1＜|ω₂|＜ω_b2，…，ω_b(n-2)＜|ω_n-1|＜ω_b(n-1)，和ω_b(n-1)＜|ω_n|，其中ω_bi∈(ω_i，ω_i+1)，i＝1，2，…，n-1为n-1个二分时不变截止频率，则分量信号x_i(t)通过式(6)和(7)解析得出：

x₁(t)＝s₁(t)，…，x_i(t)＝s_i(t)-s_i-1(t)，…，x_n(t)＝X_G(t)-s_n-1(t)i＝1，2，…，n(6)

s_i(t)＝sin[ω_bi(t)]H{X_G(t)cos[ω_bi(t)]}-cos[ω_bi(t)]H{X_G(t)sin[ω_bi(t)]}

i＝1，2，…，n-1 (7)

步骤3)对调制后的信号X_G(ω)进行逆广义傅里叶变换，即重构，可得：

步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，以得到纯调频信号；对于一个分解得到的单分量信号x_i(t)，其希尔伯特变换通过式(9)表示：

式中，P为柯西积分主值；

步骤5)对于不满足Bedrosian定理的单分量信号，希尔伯特解调存在较大误差，通过迭代希尔伯特解决这一问题；通过所述单分量信号x₁(t)构造出如式(10)所示的解析信号z(t)，z(t)的实部为原始信号，虚部为原始信号的希尔伯特变换：

式中，

和φ₁＝tan^-1{H[x_i(t)]/x_i(t)}分别为幅值函数和调频函数；解析信号的实部表示为幅值函数A₁和调频函数cosφ₁的乘积，如式(11)所示：

表示第一次迭代希尔伯特变换后的原始分量信号。将生成的调频信号

作为新的信号并继续对其进行希尔伯特变换，产生新的幅值函数A₂和调频函数cosφ₂，即：

式中，

表示第二次迭代希尔伯特变换后的原始分量信号，

步骤6)不断重复步骤5)，得到如下迭代公式：

表示第m次迭代希尔伯特变换后的原始分量信号。由于每一次迭代过程均会产生新的调频函数及幅值函数，因此只有当新的幅值函数A_m趋近于1时，迭代才会停止；此时得到的调频函数

即为单分量信号x_i(t)的纯调频信号，也近似为渐进单分量信号；

步骤7)对经过迭代希尔伯特变换处理后的纯调频信号

进行连续小波变换得到如式(14)所示的小波系数；

式中，W_x(a,b)表示小波系数；a表示尺度因子，与频率成反比关系；b表示平移因子，与时间有关；

为

的共轭复数，

表示小波母函数；小波系数表征小波母函数与信号相似的程度，小波系数越大表示两者越相似；

步骤8)由于小波系数模值在小波脊线上取得局部极大值，对任意时刻b_n，通过式(15)寻找不同尺度下小波系数模的局部极大值，从而得到对应的小波脊点(a_m,b_n)；按照同样的方法继续提取下一时刻b_n+1所对应的小波脊点(a_m+1,b_n+1)，直到信号的末端终点，最后，将提取到的局部极大值点连接起来作为信号的小波脊线，即求得渐进单分量信号的瞬时频率；

|W_x(a_r,b)|＝max|W_x(a_m,b_n)| (15)

其中，W_x(a_r,b)表示局部极大值点。

相较于现有技术，本发明具有以下有益效果：提出了一种基于广义解调解析模态分解、迭代希尔伯特变换和小波系数模局部极大值法相结合的时变结构瞬时频率提取方法，该方法能针对时变和非线性结构系统仅利用其非平稳响应信号即可实现瞬时频率的提取。在参数识别的基础上，能进一步构建损伤位置及时变损伤指标，可用于设备故障检测、结构健康监测和安全评估领域。该方法仅需测得响应信号即可实现参数提取，无需测量输入激励，具有很强的实用性和广阔的应用前景。

附图说明

图1为本发明方法的实现流程图。

图2为本发明实施例中时变悬臂梁试验装置示意图。

图3为本发明实施例中原多分量信号x(t)调制前后时频图。

图4为本发明实施例中分量信号x₁(t)，x₂(t)识别值与理论值对比图。

图5为本发明实施例中分量信号x₁(t)，x₂(t)解调前后对比图。

图6为本发明实施例中分量信号x₁(t)，x₂(t)瞬时频率识别值与理论值对比图。

图7为本发明实施例中信号原多分量信号y(t)调制前后时频图。

图8为本发明实施例中分量信号y₁(t)，y₂(t)识别值与理论值对比图。

图9为本发明实施例中分量信号y₁(t)，y₂(t)解调前后对比图。

图10为本发明实施例中分量信号y₁(t)，y₂(t)瞬时频率识别值与理论值对比图。

图11为本发明实施例中信号原多分量信号z(t)调制前后时频图。

图12为本发明实施例中分量信号z₁(t)，z₂(t)识别值与理论值对比图。

图13为本发明实施例中分量信号z₁(t)，z₂(t)解调前后对比图。

图14为本发明实施例中分量信号z₁(t)，z₂(t)瞬时频率识别值与理论值对比图。

图15为本发明实施例中信号原多分量信号s(t)调制前后时频图。

图16为本发明实施例中分量信号s₁(t)，s₂(t)识别值与理论值对比图。

图17为本发明实施例中分量信号s₁(t)，s₂(t)解调前后对比图。

图18为本发明实施例中分量信号s₁(t)，s₂(t)瞬时频率识别值与理论值对比图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步的详细说明。

本发明提供了一种时变结构瞬时频率提取方法，对于时变钢桥或组合桥梁结构，首先将广义解调定理与解析模态分解理论相结合，将具有密集模态甚至叠混的多分量信号分解为多个单分量信号，由于分解所得单分量信号通常不是纯调频信号，再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，从而实现调幅调频信号到纯调频信号的转变，最后通过小波系数模局部极大值方法精确提取瞬时频率。具体包括如下步骤：

步骤1)对于瞬时相位为φ(t)的渐进信号x(t)，其相应的解析信号表达式x_A(t)和瞬时频率f(t)分别如式(1)和(2)所示：

x_A(t)＝x(t)+jH[x(t)]＝x(t)e^j[φ(t)] (I)

式中，j为虚数单位，H[·]表示希尔伯特变换算子；通过解析信号x_A(t)来表示信号x(t)的广义傅里叶变换，即信号x(t)的调制，如式(3)所示：

式中，X_G(ω)为调制后的信号，φ(t)为与时间相关的实相位信号；为了将信号的时变瞬时频率调制到固定频率f₀，构造φ(t)，如式(4)所示：

式中，

为原多分量信号各阶瞬时频率估计值，即通过其他算法得到瞬时频率大概的估计值，可通过同步挤压小波变换或其他瞬时频率识别方法得到；与理论值相比，虽然识别值

存在一定的误差，但经过调制后的信号的瞬时频率总体上仍然会集中在固定频率f₀附近；

步骤2)对调制后的信号X_G(ω)进行解析模态分解，选取该信号的实部X_G(t)，假设X_G(t)是由n个分量x_i(t)组成的多分量信号：

X_G(t)＝x₁(t)+x₂(t)+…+x_i(t)+…+x_n(t) i＝1，2，…，n (5)

若分量信号的频率ω₁，ω₂，…，ω_n满足：|ω₁|＜ω_b1，ω_b1＜|ω₂|＜ω_b2，…，ω_b(n-2)＜|ω_n-1|＜ω_b(n-1)，和ω_b(n-1)＜|ω_n|，其中ω_bi∈(ω_i，ω_i+1)，i＝1，2，…，n-1为n-1个二分时不变截止频率，则分量信号x_i(t)通过式(6)和(7)解析得出：

x₁(t)＝s₁(t)，…，x_i(t)＝s_i(t)-s_i-1(t)，…，x_n(t)＝X_G(t)-s_n-1(t) i＝1，2，…，n(6)

s_i(t)＝sin[ω_bi(t)]H{X_G(t)cos[ω_bi(t)]}-cos[ω_bi(t)]H{X_G(t)sin[ω_bi(t)]}

i＝1，2，…，n-1 (7)

步骤3)对调制后的信号X_G(ω)进行逆广义傅里叶变换，即重构，可得：

步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，以得到纯调频信号；对于一个分解得到的单分量信号x_i(t)，其希尔伯特变换通过式(9)表示：

式中，P为柯西积分主值；

步骤5)对于不满足Bedrosian定理的单分量信号，希尔伯特解调存在较大误差，通过迭代希尔伯特解决这一问题；通过所述单分量信号x_i(t)构造出如式(10)所示的解析信号z(t)，z(t)的实部为原始信号，虚部为原始信号的希尔伯特变换：

式中，

和φ₁＝tan^-1{H[x_i(t)]/x_i(t)}分别为幅值函数和调频函数；解析信号的实部表示为幅值函数A1和调频函数cosφ₁的乘积，如式(11)所示：

表示第一次迭代希尔伯特变换后的原始分量信号。将生成的调频信号

作为新的信号并继续对其进行希尔伯特变换，产生新的幅值函数A₂和调频函数cosφ₂，即：

式中，

表示第二次迭代希尔伯特变换后的原始分量信号，

步骤6)不断重复步骤5)，得到如下迭代公式：

表示第m次迭代希尔伯特变换后的原始分量信号。由于每一次迭代过程均会产生新的调频函数及幅值函数，因此只有当新的幅值函数A_m趋近于1时，迭代才会停止；此时得到的调频函数

即为单分量信号x_i(t)的纯调频信号，也近似为渐进单分量信号；

步骤6)不断重复步骤5)，得到如下迭代公式：

由于每一次迭代过程均会产生新的调频函数及幅值函数，因此只有当新的幅值函数A_n趋近于1时，迭代才会停止；此时得到的调频函数x_n(t)＝cosφ_n即为单分量信号x₁(t)的纯调频信号，也近似为渐进单分量信号；

步骤7)对经过迭代希尔伯特变换处理后的纯调频信号x(t)进行连续小波变换得到如式(14)所示的小波系数；

式中，W_x(a,b)表示小波系数；a表示尺度因子，与频率成反比关系；b表示平移因子，与时间有关；

为

的共轭复数，

表示小波母函数；小波系数表征小波母函数与信号相似的程度，小波系数越大表示两者越相似；

步骤8)由于小波系数模值在小波脊线上取得局部极大值，对任意时刻b_n，通过式(15)寻找不同尺度下小波系数模的局部极大值，从而得到对应的小波脊点(a_m,b_n)；按照同样的方法继续提取下一时刻b_n+1所对应的小波脊点(a_m+1,b_n+1)，直到信号的末端终点，最后，将提取到的局部极大值点连接起来作为信号的小波脊线，即求得渐进单分量信号的瞬时频率；

|W_x(a_r,b)|＝max|W_x(a_m,b_n)| (15)

其中，W_x(a_r,b)表示局部极大值点。

下面举具体实施例对本发明进行说明。

1)现有多分量调幅调频信号x(t)＝x₁(t)+x₂(t)＝(t+6)cos[2π(t²+12t)]+(3t+18)cos[2π(0.6t²+7t)]。设定采样频率f_s＝100Hz，采样时间t＝10s，采样间隔dt＝0.01s，考虑到实际工程中有噪声干扰，因此添加20％高斯白噪声(7dB)，噪声强度由信噪比定义：

式中，A_signal和A_noise分别代表信号和噪声的均方根值，噪声水平是指

与

之间的比值。

首先选用复Morlet小波对添加了20％高斯白噪声原信号进行同步挤压小波变换变换得到原信号的时频图，如图3(a)所示。从图中可以看出：该多分量信号包括了两个单分量信号，且两个分量信号的频率随时间呈线性变化，不仅是时变而且也是模态叠混的，因此AMD无法分解该信号。

采用广义解调解析模态分解定理对该信号进行处理。首先对原信号进行同步挤压小波变换得到原多分量信号的各阶瞬时频率的估计值，然后根据估计值选取合适的固定频率f₀＝4Hz对原多分量信号进行信号调制，调制后的时频图如图3(b)所示。从图3(b)中可以看出：分量信号x₂(t)经调制后的瞬时频率在4Hz左右浮动，而分量信号x₁(t)的频率大概介于(10Hz，20Hz)之间。此时两分量信号之间未出现模态叠混，因此选择截止频率4.1Hz对调制后的多分量信号进行AMD分解，即可成功对信号进行分离。然后进行逆广义傅里叶变换便后可得到原分量信号x₁(t)和x₂(t)，两分量信号的分解结果分别如图4(a)和(b)所示。

紧接着对得到的单分量调幅调频信号进行迭代希尔伯特变换解调，解调结果如图5(a)和(b)所示，从图中可以看出分量信号在解调后的幅值已接近1，可近似认为渐进单分量信号。

最后通过小波系数局部模极大值法对得到的近似渐进单分量信号进行瞬时频率提取，结果如图6(a)和(b)所示。由图6(a)和(b)可以看出：瞬时频率识别值与理论值吻合得很好该算法的正确性得到了验证。

2)考虑多分量调幅调频信号y(t)＝y₁(t)+y₂(t)＝(t+9)cos(10πt+2πt²)+(2t+1)cos[20πt+3πt²+sin(2πt)]，已添加20％高斯白噪声(7dB)，由信号表达式可知该多分量信号由两个调频调幅信号组成。设定采样频率f_s＝100Hz，采样时间t＝10s，采样间隔dt＝0.01s。

首先选用复Morlet小波对该信号进行同步挤压小波变换得到其时频图，如图7(a)所示。从图7(a)中可获知：该多分量信号由两个时变信号组成，分别为时变正弦信号和一个时变线性信号。由于两分量信号出现模态叠混，传统AMD无法分离模态叠混信号。

采用广义解调解析模态分解定理对该信号进行处理。首先对原信号进行同步挤压小波变换得到原多分量信号的各阶瞬时频率的估计值，然后根据估计值选取合适的固定频率f₀＝4Hz对原多分量信号进行信号调制，调制以后时频图如图7(b)所示：分量信号y₁(t)经调制以后，瞬时频率在4Hz左右浮动，而分量信号y₂(t)经调制以后，瞬时频率变化范围大致在(10Hz，20Hz)之间呈正弦规律变化，而且两分量信号未模态叠混。选取二分截止频率5Hz，通过AMD对信号进行分解。随后对得到的分量信号进行逆广义傅里叶变换便可得到原调幅调频分量信号y₁(t)和y₂(t)。两分量信号的分解结果如图8(a)和(b)所示。由于添加了白噪声，信号在初始阶段吻合度不高，但是在其他阶段吻合较好。

由于分解得到的分量信号并不是纯调频信号，因此通过迭代希尔伯特变换对两分量信号进行解调，结果如图9(a)和(b)所示，分量信号在解调后的幅值已接近1，可近似认为渐进单分量信号。

最后采用小波系数局部模极大值对解调后所得的近似纯调频信号进行瞬时频率提取，结果如图10(a)和(b)所示。由于信号中添加了20％白噪声，初始阶段识别值不可避免出现了端点效应，但整体和理论值还是比较吻合。

3)继续考虑多分量调幅调频信号z(t)＝z₁(t)+z₂(t)＝(3t+1)cos(18πt+4πt²)+(2t+2)cos[2πt²+sin(2πt)+8πt]，同样添加20％高斯白噪声(7dB)。设定采样频率f_s＝100Hz，采样时间t＝10，采样间隔dt＝0.01s。

首先选用复Morlet小波对原添加白噪声的多分量信号进行同步挤压小波变换得到其小波量图，如图11(a)所示。从图中可以看出：原多分量信号包含了两个时变信号，其中一个是呈正弦变化，一个呈线性变化，且两分量信号出现叠混。

通过同步挤压小波变换得到分量信号的瞬时频率估计值以后，根据瞬时频率估计值选择固定频率f₀＝4Hz对原信号进行调制，调制后的效果如图11(b)所示：单分量信号z₂(t)调制后，其瞬时频率在4Hz左右浮动，而单分量信号z₁(t)的频率大致介于(9Hz，30Hz)之间。选定截止频率为4.5Hz，采用AMD对信号进行分解，并对分解的信号进行逆广义傅里叶变换得到原分量信号z₁(t)和z₂(t)，如图12(a)和(b)所示。

在得到两个调幅调频信号以后，引入迭代希尔伯特变换对单分量信号进行解调，解调过程结束后最终得到两个近似纯调频分量信号，结果如图13(a)和(b)所示。

最后通过小波系数局部模极大值对分量信号的瞬时频率进行提取，结果如图14(a)和(b)所示，同样地除了不可避免的端点效应以外，整体上和理论值也是很吻合的。

4)继续考虑多分量调幅调频信号s(t)＝s₁(t)+s₂(t)＝(3t+1)cos[πt²+sin(2πt)+12πt]+(2t)cos[2πt²+sin(2πt)+30πt]，已添加20％高斯白噪声。设定采样频率f_s＝100Hz，采样时间t＝10s，采样间隔dt＝0.01s。

首先选择复Morlet小波对该信号做同步挤压小波变换，得到小波量图，如图15(a)所示。由图15(a)可知：原多分量信号包含了两个正弦规律变化的时变分量信号，且两个信号出现了模态叠混，因此AMD无法直接处理该信号。

首先通过同步挤压小波变换得到信号的瞬时频率估计值，然后通过估计值确定固定频率f₀＝3Hz，紧接着对原多分量信号进行调制，调制以后的信号时频图如图15(b)所示。由图15(b)可知：原来两条正弦规律变化的信号经过调制以后均成了近似线性信号，其中信号s₁(t)的瞬时频率在3Hz左右浮动，而信号s₂(t)的频率则随着时间呈现线性增大的趋势。接下来设定截止频率为8Hz对调制后的多分量信号进行分解，再对分解的结果进行逆傅里叶变换得到原单分量调幅调频单分量信号，结果如图16(a)和(b)所示。

然后引入迭代希尔伯特变换对单分量调幅调频分量信号进行解调，结果如图17(a)和(b)所示，解调完成后，两分量信号近似接近渐进单分量信号。

最后通过小波局部模极大值方法提取各分量信号的瞬时频率，结果如图18(a)和(b)所示。除了最开始的端点效应以外，识别值与理论值整体比较吻合，验证了该算法的准确性和适用性。

为了验证上述方法，本发明还提供了基于上述时变结构瞬时频率提取方法的时变悬臂梁试验装置，该装置

包括一根悬臂梁(本实施例为铝合金悬臂梁)、一块悬挂的磁铁、质量块、加速度传感器、力锤及采集仪；所述悬臂梁每隔一段距离设置一个锤击点，并通过锤击提供激励；悬挂的磁铁用于在设定时刻将放置于悬臂梁上的质量块吸起以实现结构质量的改变；通过放置于悬臂梁跨中的加速度传感器采集加速度响应信号，然后采用广义解调解析模态分解定理、迭代希尔伯特变换和小波系数模局部极大值相结合的方法来提取时变悬臂梁结构的瞬时频率，最后与有限元分析结果或理论解进行比较，以验证该方法的正确性和有效性。该试验装置的工作步骤如下：

步骤1)试验装置参数设定如下：所用铝合金悬臂梁的重量为0.81kg，截面尺寸为40mm×15mm，长度为500mm；悬臂端的质量块重1kg，在质量块上方放置一块悬挂的磁铁；

步骤2)梁上每隔100mm设置一锤击点，然后通过力锤敲击悬臂梁的自由端并在某一个时刻放下细绳使永磁铁垂直靠近质量块并将其吸起以实现结构质量的改变；设定采样频率为2000Hz，在梁的跨中位置安装加速度传感器以采集加速度冲击响应；试验前，基于“冻结法”提前测得带质量块与不带质量块时悬臂梁的基频作为理论值；

步骤3)试验时，先在悬臂梁末端用力锤施加激振力，然后采集响应信号，试验开始2s时通过磁铁吸引质量块以改变悬臂梁质量，最后在试验开始2.3s后再次用力锤锤击悬臂梁以避免响应信号衰减过快；

步骤4)通过所提方法对传感器采集到的响应信号进行分析处理，最终识别时变铝梁结构的瞬时频率，并与理论值相比较来验证所提方法的有效性。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。

Claims (2)

1.一种时变结构瞬时频率提取方法，其特征在于，对于时变钢桥或组合桥梁结构，首先将广义解调定理与解析模态分解理论相结合，将具有密集模态甚至叠混的多分量信号分解为多个单分量信号，再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，从而实现调幅调频信号到纯调频信号的转变，最后通过小波系数模局部极大值方法精确提取瞬时频率。

2.根据权利要求1所述的一种时变结构瞬时频率提取方法，其特征在于，具体包括如下步骤：

步骤1)对于瞬时相位为φ(t)的渐进信号x(t)，其相应的解析信号表达式x_A(t)和瞬时频率f(t)分别如式(1)和(2)所示：

x_A(t)＝x(t)+jH[x(t)]＝x(t)e^j[φ(t)] (1)

式中，j为虚数单位，H[·]表示希尔伯特变换算子；通过解析信号x_A(t)来表示信号x(t)的广义傅里叶变换，即信号x(t)的调制，如式(3)所示：

式中，X_G(ω)为调制后的信号，φ(t)为与时间相关的实相位信号；为了将信号的时变瞬时频率调制到固定频率f₀，构造φ(t)，如式(4)所示：

式中，

为原多分量信号各阶瞬时频率估计值；

步骤2)对调制后的信号X_G(ω)进行解析模态分解，选取该信号的实部X_G(t)，假设X_G(t)是由n个分量x_i(t)组成的多分量信号：

X_G(t)＝x₁(t)+x₂(t)+…+x_i(t)+…+x_n(t)i＝1，2，…，n (5)

若分量信号的频率ω₁，ω₂，…，ω_n满足：|ω₁|＜ω_b1，ω_b1＜|ω₂|＜ω_b2，…，ω_b(n-2)＜|ω_n-1|＜ω_b(n-1)，和ω_b(n-1)＜|ω_n|，其中ω_bi∈(ω_i，ω_i+1)，i＝1，2，…，n-1为n-1个二分时不变截止频率，则分量信号x_i(t)通过式(6)和(7)解析得出：

x₁(t)＝s₁(t)，…，x_i(t)＝s_i(t)-s_i-1(t)，…，x_n(t)＝X_G(t)-s_n-1(t)i＝1，2，…，n (6)

s_i(t)＝sin[ω_bi(t)]H{X_G(t)cos[ω_bi(t)]}-cos[ω_bi(t)]H{X_G(t)sin[ω_bi(t)]}i＝1，2，…，n-1 (7)

步骤3)对调制后的信号X_G(ω)进行逆广义傅里叶变换，即重构，可得：

步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调，以得到纯调频信号；对于一个分解得到的单分量信号x_i(t)，其希尔伯特变换通过式(9)表示：

式中，P为柯西积分主值；

步骤5)对于不满足Bedrosian定理的单分量信号，希尔伯特解调存在较大误差，通过迭代希尔伯特解决这一问题；通过所述单分量信号x_i(t)构造出如式(10)所示的解析信号z(t)，z(t)的实部为原始信号，虚部为原始信号的希尔伯特变换：

式中，

和φ₁＝tan^-1{H[x_i(t)]/x_i(t)}分别为幅值函数和调频函数；解析信号的实部表示为幅值函数A₁和调频函数cosφ₁的乘积，如式(11)所示：

表示第一次迭代希尔伯特变换后的原始分量信号；将生成的调频信号

作为新的信号并继续对其进行希尔伯特变换，产生新的幅值函数A₂和调频函数cosφ₂，即：

式中，

表示第二次迭代希尔伯特变换后的原始分量信号，

步骤6)不断重复步骤5)，得到如下迭代公式：

表示第m次迭代希尔伯特变换后的原始分量信号；由于每一次迭代过程均会产生新的调频函数及幅值函数，因此只有当新的幅值函数A_m趋近于1时，迭代才会停止；此时得到的调频函数

即为单分量信号x_i(t)的纯调频信号，也近似为渐进单分量信号；

步骤7)对经过迭代希尔伯特变换处理后的纯调频信号

进行连续小波变换得到如式(14)所示的小波系数；

式中，W_x(a,b)表示小波系数；a表示尺度因子，与频率成反比关系；b表示平移因子，与时间有关；

为

的共轭复数，

表示小波母函数；小波系数表征小波母函数与信号相似的程度，小波系数越大表示两者越相似；

步骤8)由于小波系数模值在小波脊线上取得局部极大值，对任意时刻b_n，通过式(15)寻找不同尺度下小波系数模的局部极大值，从而得到对应的小波脊点(a_m,b_n)；按照同样的方法继续提取下一时刻b_n+1所对应的小波脊点(a_m+1,b_n+1)，直到信号的末端终点，最后，将提取到的局部极大值点连接起来作为信号的小波脊线，即求得渐进单分量信号的瞬时频率；

|W_x(a_r,b)|＝max|W_x(a_m,b_n)| (15)

其中，W_x(a_r,b)表示局部极大值点。

Images