CN113887360A - 一种基于迭代扩展频散模态分解的频散波提取方法 - Google Patents

一种基于迭代扩展频散模态分解的频散波提取方法 Download PDF

Info

Publication number
CN113887360A
CN113887360A CN202111114136.1A CN202111114136A CN113887360A CN 113887360 A CN113887360 A CN 113887360A CN 202111114136 A CN202111114136 A CN 202111114136A CN 113887360 A CN113887360 A CN 113887360A
Authority
CN
China
Prior art keywords
frequency
dispersion
time
key
frequency dispersion
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111114136.1A
Other languages
English (en)
Other versions
CN113887360B (zh
Inventor
牛刚
姜源
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tongji University
Original Assignee
Tongji University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tongji University filed Critical Tongji University
Priority to CN202111114136.1A priority Critical patent/CN113887360B/zh
Publication of CN113887360A publication Critical patent/CN113887360A/zh
Application granted granted Critical
Publication of CN113887360B publication Critical patent/CN113887360B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Signal Processing (AREA)
  • Discrete Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Algebra (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种基于迭代扩展频散模态分解的频散波提取方法，包括：获取频散波时域信号；分别进行快速傅里叶变换和离散短时傅里叶变换，获得频域信号和时频表达式；通过双向快速单脊线提取法，初步估计能量最显著的关键时频脊线；通过迭代扩展频散模态分解，分离重构关键频散模态及其群延迟；重复获取残差信号以及分离重构关键频散模态，直至满足停止迭代准则；对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波时域波形。与现有技术相比，本发明在无需预知频散波模态数量的前提下，通过准确提取和分析信号中的频散模态，高效地对多分量非平稳频散波进行分解，在结构健康监测、水下声学、生物医学、地球物理学等领域有重要应用。

Description

一种基于迭代扩展频散模态分解的频散波提取方法
技术领域
本发明涉及信号处理领域，尤其是涉及一种基于迭代扩展频散模态分解的频散波提取方法。
背景技术
频散信号在结构健康监测、水下声学、生物医学、地球物理学等领域中都有广泛的研究。与时域调频信号不同，频散信号表现出的是与频率相关的特性，其中最重要的是信息是群延迟(Group Delay)。群延迟定义为相位函数对频率的导数，包含众多有关信号源与系统的重要信息，因此可以利用群延迟进行信号源的定位与系统辨识。然而，当宽带信号通过非线性系统时，每个频率都会产生不同的时间延迟，从而导致群延迟在频率上发生变化，从时频图上来看，信号的时频脊线不再是与频率轴平行的直线，而是与频率轴相关的曲线，这使得窄带滤波方法不再适用。因此，准确分离这些非平稳频散波并提取它们的群延迟曲线对于实际应用至关重要。
虽然众多学者和工程人员提出了有关频散信号的提取方法，但是现有的频散信号提取方法存在着算法计算消耗大、时频分辨率低、计算时间长等问题。目前，国内有关频散信号提取技术的发明专利申请有：
授权日为2021年8月6日，授权公告号为CN110045014B的中国专利中，公开了一种名称为“基于贝叶斯学习的Lamb波频散消除方法及其系统”。其通过构造多模态频散传播字典和非频散传播字典，将频散多模态Lamb波信号进行稀疏表示，并利用贝叶斯学习方法对字典系数进行系数求解，从而得到频散波的稀疏表达并进行频散消除。然而，该方法的计算结果与字典构建的匹配性有非常强的依赖，同时字典系数的稀疏化求解需要耗费大量的计算时间和计算资源消耗。
公开日为2021年5月28日，公开号为CN112858481A的中国专利中，公开了一种名称为“一种基于时间频率域评估板材中裂纹深度的方法”。其采用时频分析方法求解Lamb波频散曲线。然而传统的非参数化时频分析方法，在群延时快速波动时，会造成分辨率降低的现象；参数化时频分析方法虽然能够提高时频分辨率，但是计算复杂度大大提高，给实际工程应用带来巨大困难。
发明内容
本发明的目的就是为了克服上述现有技术存在计算消耗大、时频分辨率低、计算时间长的缺陷而提供一种基于迭代扩展频散模态分解的频散波提取方法。
本发明的目的可以通过以下技术方案来实现：
一种基于迭代扩展频散模态分解的频散波提取方法，包括以下步骤：
步骤1：获取频散波时域信号；
步骤2：对所述频散波时域信号进行快速傅里叶变换，获得对应的频域信号；
步骤3：对所述频散波时域信号进行离散短时傅里叶变换，获得时频表达式；
步骤4：初始化预设的停止迭代准则，设定频散波信号成分数量m＝0；
步骤5：通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线；
步骤6：通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线的初步估计值，分离重构关键频散模态及其群延迟；
步骤7：根据重构的关键频散模态，从频域信号中得到残差信号，对残差信号重复依次执行步骤5至步骤7，重复获取残差信号以及进行关键频散模态的分离重构，直至满足所述停止迭代准则；
步骤8：对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形。
进一步地，步骤2具体包括以下步骤：
S201：对频散波时域信号S(t)进行快速傅里叶变换，得到对应的频域信号S(f)；
S202：通过频散波时域信号的长度Nt，确定频域信号长度Nf，该频域信号长度Nf的计算表达式为：
Nf＝floor(Nt/2)+1
式中，floor(·)为向下取整符号；
S203：保留频域信号的前Nf项，第Nf+1项至第Nt项置零。
进一步地，步骤3具体包括以下步骤：
S301：设置窗长度Nwin，确定傅里叶点数为2Nf-1，构造相应窗函数W(t)；
S302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)。
进一步地，所述步骤5包括以下步骤：
步骤501：进行参数初始化；
步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(tem,fem)并沿频率方向双向搜索Δt范围内时间维度的能量极大值点；
步骤503：获得时频脊线作为关键模态群延时的初始估计值
Figure BDA0003274874080000031
进一步地，步骤501中，所述参数为最大时间允许波动范围Δt；
步骤502中，所述能量最大值点(tem,fem)的计算表达式为：
Figure BDA0003274874080000032
式中，t为时间，f为幅值；
所述沿频率方向双向搜索Δt范围内时间维度的能量极大值点的计算表达式为：
tR＝tem,tL＝tem
Figure BDA0003274874080000033
式中，Nf为频域信号长度。
进一步地，所述步骤6具体包括以下步骤：
步骤601：初始化参数，获取频域频散波信号S(f)和关键模态群延时初始值
Figure BDA0003274874080000034
步骤602：构造关键频散模态
Figure BDA0003274874080000035
并建立解析幅值
Figure BDA0003274874080000036
的冗余傅里叶展开式；
步骤603：将所述
Figure BDA0003274874080000037
的冗余傅里叶展开式带入关键频散模态
Figure BDA0003274874080000038
中构造矩阵表达式，通过求解最优化问题获得傅里叶参数矩阵
Figure BDA0003274874080000039
并重构得
Figure BDA00032748740800000310
步骤604：利用所述重构幅值
Figure BDA00032748740800000311
通过解卷绕求导方法估计群延时误差
Figure BDA00032748740800000312
对其进行滤波光滑然后更新群延时
Figure BDA00032748740800000313
步骤605：将
Figure BDA00032748740800000314
代入步骤603更新核矩阵Kk+1、参数矩阵
Figure BDA00032748740800000315
和解析幅值矩阵
Figure BDA00032748740800000316
重复计算达到预设的迭代次数后停止，然后求解重构关键频散模态skey(f)；
步骤606：更新频散波信号成分数量m＝m+1。
进一步地，步骤601中，所述参数包括正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L；
步骤602中，所述关键频散模态
Figure BDA0003274874080000041
的计算表达式为：
Figure BDA0003274874080000042
其中，k＝0,...,end为迭代次数，
Figure BDA0003274874080000043
为解析幅值：
Figure BDA0003274874080000044
其中，
Figure BDA0003274874080000045
为初始相位；
解析幅值
Figure BDA0003274874080000046
的冗余傅里叶展开式为：
Figure BDA0003274874080000047
其中，
Figure BDA0003274874080000048
fs为时域信号的采样频率；
步骤603中，构造的矩阵表达式为：
Figure BDA0003274874080000049
Figure BDA00032748740800000410
Kk＝Ek·T为Nf×(2L+1)的核矩阵，其中，
Figure BDA00032748740800000411
Figure BDA00032748740800000412
所述最优化问题为：
Figure BDA00032748740800000413
获得的所述傅里叶参数矩阵
Figure BDA00032748740800000414
的计算表达式为：
Figure BDA00032748740800000415
所述重构幅值
Figure BDA00032748740800000416
的计算表达式为：
Figure BDA00032748740800000417
式中，I为单位矩阵，(·)H为矩阵的共轭转置；
步骤604中，取
Figure BDA00032748740800000418
的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差
Figure BDA00032748740800000419
Figure BDA00032748740800000420
式中，∠·为求复数的相位运算，unwrap(·)为相位解卷绕运算，将相位值由[-π,π]区间叠加到[0,+∞)区间；
对
Figure BDA0003274874080000051
进行滤波光滑并更新瞬时频率
Figure BDA0003274874080000052
该步骤的计算表达式为：
Figure BDA0003274874080000053
Figure BDA0003274874080000054
其中
Figure BDA0003274874080000055
为改进二阶差分算子；
步骤605中，所述关键频散模态skey(f)的重构表达式为：
Figure BDA0003274874080000056
式中，end为预设的迭代次数，Kend为更新end次后的核矩阵，
Figure BDA0003274874080000057
为更新end次后的参数矩阵。
进一步地，所述步骤7包括以下步骤：
步骤701：将所述重构关键频散模态skey(f)作为第m个频散模态分量sm(f)，并将其从频散波频域信号S(f)中提取，得到残差信号Sr(f)：
Sr(f)＝S(f)-sm(f)
步骤702：计算残差信号能量与频散波频域信号能量之比：
Figure BDA0003274874080000058
步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复依次执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。
进一步地，所述步骤8包括以下步骤：
步骤801：扩展所有重构的频散模态分量sm(f)；
步骤802：对所有扩展后的频散模态分量进行快速傅里叶变换逆变换，得到频散波的时域波形。
进一步地，步骤801中，扩展后的频散模态分量
Figure BDA0003274874080000059
的计算表达式为：
Figure BDA00032748740800000510
式中，conj[·]表示求复数的共轭；
步骤802中，取快速傅里叶变换逆变换计算结果的实部作为对应的频散波时域波形
Figure BDA0003274874080000061
该频散波时域波形
Figure BDA0003274874080000062
的计算表达式为：
Figure BDA0003274874080000063
式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。
本发明新提出了一种迭代扩展频散模态分解法(IEDMD)，通过双向快速单脊线提取法精确提取频散模态时频脊线，实现了精确高分辨率的非稳态频散波提取。与现有技术相比，本发明具有以下优点：
(1)本发明相比同类方法，群延迟提取精度提高；
(2)本发明无需预知频散波模态的数量，通过能量比方法自适应地确定提取数量，更有利于实际工程的应用；
(3)本发明在强噪声背景下精确重构频散波的时域波形和群延时，有效地抑制了原信号中的噪声干扰和成分混叠，信息丰富完整，有利于对其进行深入分析；
(4)本发明计算速度相比同类方法提升2倍以上，瞬时频率提取精度提高20％以上，模态重构精度提高15％以上，特征阶次误差小于1％并且计算精度提高20％以上，无需较多人为干预，有利于嵌入式开发和实际工程应用。
附图说明
图1为本发明实施例中提供一种基于迭代扩展频散模态分解的频散波提取方法的流程示意图；
图2(a)为本发明实施例中频散波仿真信号时域图；
图2(b)为本发明实施例中频散波仿真模态1时域图；
图2(c)为本发明实施例中频散波仿真模态2时域图；
图3为本发明实施例中频散波的时频图；
图4(a)为本发明实施例中第一频散模态IEDMD群延迟估计图；
图4(b)为本发明实施例中第一频散模态IEDMD群延时计算误差曲线图；
图4(c)为本发明实施例中第一频散模态IEDMD时域信号及计算误差图；
图4(d)为本发明实施例中第一频散模态GDMD时域信号及计算误差图；
图5(a)为本发明实施例中第二频散模态IEDMD群延迟估计图；
图5(b)为本发明实施例中第二频散模态IEDMD群延时计算误差曲线图；
图5(c)为本发明实施例中第二频散模态IEDMD时域信号及计算误差图；
图5(d)为本发明实施例中第二频散模态GDMD时域信号及计算误差图；
图中，Time为时间，Amplitude为振幅，Freq为频率，Frequency为频率，Dispersivemode为频散模态，Residual为残差信号，True为真实信号。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
实施例1
如图1所示，本实施例提供一种基于迭代扩展频散模态分解的频散波提取方法，在无需预知频散波模态数量的前提下，通过准确提取和分析信号中的频散模态，高效地对多分量非平稳频散波进行分解，在结构健康监测、水下声学、生物医学、地球物理学等领域有重要应用。
假设采集到的频散波时域信号为S(t)，该方法具体包含以下步骤：
步骤1：根据实际应用场景，按需要安装传感器，获取相关领域的频散波时域信号S(t)；
步骤2，对所述时域信号进行快速傅里叶变换，获得频散波对应的频域信号，具体包括以下步骤：
步骤201：对时域频散波信号S(t)进行快速傅里叶变换(FFT)，得到对应的频域信号S(f)；
步骤202：通过时域信号长度Nt，确定频域信号长度Nf：
Nf＝floor(Nt/2)+1
式中，floor(·)为向下取整符号；
步骤203：保留频域信号的前Nf项，第Nf+1项至第Nt项置零。
步骤3，对所述时域信号进行离散短时傅里叶变换，获得时频表达式，具体包括以下步骤：
步骤301：设置窗长度Nwin，确定傅里叶点数为2Nf-1，构造相应窗函数W(t)；
步骤302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换(STFT)，获得时频表达式TFR(t,f)。
步骤4，设定迭代停止判别值ε，设定频散波信号成分数量m＝0。
步骤5，通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线，具体包括以下步骤：
步骤501：初始化参数最大时间允许波动范围Δt；
步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(tem,fem)：
Figure BDA0003274874080000081
沿频率方向双向搜索Δt范围内时间维度的能量极大值点：
tR＝tem,tL＝tem
Figure BDA0003274874080000082
步骤503：获得时频脊线
Figure BDA0003274874080000083
作为关键模态群延时的初始值。
步骤6，通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线初始值，分离重构关键频散模态及其群延迟，具体包括以下步骤：
步骤601：初始化正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L，输入频域频散波信号S(f)和关键模态群延时初始值
Figure BDA0003274874080000084
步骤602：构造关键频散模态
Figure BDA0003274874080000085
Figure BDA0003274874080000086
其中，k＝0,...,end为迭代次数，
Figure BDA0003274874080000087
为解析幅值：
Figure BDA0003274874080000088
其中，
Figure BDA0003274874080000089
为初始相位；然后建立解析幅值
Figure BDA00032748740800000810
的冗余傅里叶展开式：
Figure BDA00032748740800000811
其中，
Figure BDA00032748740800000812
fs为时域信号的采样频率；
步骤603：将所述
Figure BDA00032748740800000813
展开模型带入
Figure BDA00032748740800000814
中构造矩阵表达式：
Figure BDA0003274874080000091
Figure BDA0003274874080000092
Kk＝Ek·T为Nf×(2L+1)的核矩阵：
Figure BDA0003274874080000093
Figure BDA0003274874080000094
然后求解最优化问题：
Figure BDA0003274874080000095
获得傅里叶参数矩阵
Figure BDA0003274874080000096
并重构得
Figure BDA0003274874080000097
Figure BDA0003274874080000098
Figure BDA0003274874080000099
其中，I为单位矩阵，(·)H为矩阵的共轭转置；
步骤604：取
Figure BDA00032748740800000910
的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差
Figure BDA00032748740800000911
Figure BDA00032748740800000912
式中，∠·为求复数的相位运算，unwrap(·)为相位解卷绕运算，即将相位值由[-π,π]区间叠加到[0,+∞)区间；
对
Figure BDA00032748740800000913
进行滤波光滑并更新瞬时频率
Figure BDA00032748740800000914
Figure BDA00032748740800000915
Figure BDA00032748740800000916
其中
Figure BDA00032748740800000917
为改进二阶差分算子；
步骤605：将所述
Figure BDA00032748740800000918
跳转至步骤603更新核矩阵Kk+1、参数矩阵
Figure BDA00032748740800000919
和解析幅值矩阵
Figure BDA00032748740800000920
重复计算end次后停止，然后求解重构关键频散模态skey(f)：
Figure BDA00032748740800000921
步骤606：更新频散波信号成分数量m＝m+1；
步骤7，从所述频散波频域信号中提取出重构的关键频散模态分量，对残余信号重复进行步骤5和步骤6，直到达到停止迭代准则，具体包括以下步骤：
步骤701：将所述重构关键频散模态skey(f)作为第m个频散模态分量sm(f)，并将其从频散波频域信号S(f)中提取，得到残差信号Sr(f)：
Sr(f)＝S(f)-sm(f)
步骤702：计算残差信号能量与频散波频域信号能量之比：
Figure BDA0003274874080000101
步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。
步骤8，对所述所有重构关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形，具体包括以下步骤：
步骤801：将所有重构频散模态分量sm(f)在频域上进行扩展，得到扩展频散模态分量
Figure BDA0003274874080000102
Figure BDA0003274874080000103
式中，conj[·]表示求复数的共轭；
步骤802：对所述的扩展频散模态分量
Figure BDA0003274874080000104
逐一进行快速傅里叶变换逆变换(IFFT)，并取IFFT计算结果的实部作为对应的频散波时域波形
Figure BDA0003274874080000105
Figure BDA0003274874080000106
式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。
本实施例的具体实施结果如下：
如图2(a)所示，频散波仿真信号S(t)的信噪比为0dB，采样频率100Hz，采样时间15s；如图2(b)和图2(c)所示，频散波仿真信号S(t)中存在2个频散模态。如图3所示，两个频散模态在时间上和频率上都存在重叠。采用双向快速单脊线提取法，初步估计图3中能量最强的时频脊线，并利用IEDMD算法重构，可以得到精确的频散波模态分量及其群延迟。如图4(a)和图4(b)所示，IEDMD提取的第一频散模态群延迟与真实群延迟匹配度极高，与同类方法相比，群延迟的提取精度有显著提高，IEDMD方法的整体群延迟估计误差在1％以内，而GDMD方法在信号端点的群延时估计误差达到了15％；如图4(c)和图4(d)所示，IEDMD提取的第一频散模态时域信号与真实信号的误差极小，而GDMD提取的第一频散模态时域信号在信号端点产生较大误差。采用平均相对误差(RE_mean)指标，评价IEDMD对第一模态分量群延迟估计的准确性，如表1所示。与同类方法相比，第一模态分量群延迟估计的相对误差降低了36％。
表1IEDMD第一模态分量群延迟提取精度比较
Figure BDA0003274874080000111
将第一频散模态估计值从原始信号中提取，并对残差信号继续采用IEDMD方法提取第二模态分量。如图5(a)和图5(b)所示，IEDMD提取的第二频散模态群延迟与真实群延迟匹配度极高，与同类方法相比，群延迟的提取精度有显著提高，IEDMD方法的整体群延迟估计误差在0.2％以内，而GDMD方法在信号端点的群延时估计误差达到了50％；如图5(c)和图5(d)所示，IEDMD提取的第二频散模态时域信号与真实信号的误差极小，而GDMD提取的第二频散模态时域信号在信号端点产生较大误差。采用平均相对误差(RE_mean)指标，评价IEDMD对第二模态分量群延迟估计的准确性，如表2所示。与同类方法相比，第二模态分量群延迟估计的相对误差降低了95％。
表2IEDMD第二模态分量群延迟提取精度比较
Figure BDA0003274874080000112
采用均方根误差(RMSE)、信噪比(SNR)和相关系数(CC)对IEDMD方法的信号提取准确度进行评价，如表3所示。与同类方法相比，RMSE指标降低了32％，SNR指标提高了9.3％，CC指标提高了3.5％，表明IEDMD方法在大幅度降低频散波信号估计误差的同时，提取信号与原信号存在更高的相关度，并抑制了噪声。
表3IEDMD信号提取精度比较
Figure BDA0003274874080000113
综上所述，与同类方法相比，本方法的各项性能指标均有可观的提升。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，包括以下步骤：
步骤1：获取频散波时域信号；
步骤2：对所述频散波时域信号进行快速傅里叶变换，获得对应的频域信号；
步骤3：对所述频散波时域信号进行离散短时傅里叶变换，获得时频表达式；
步骤4：初始化预设的停止迭代准则，设定频散波信号成分数量m＝0；
步骤5：通过双向快速单脊线提取法，初步估计所述时频表达式中能量最显著的关键时频脊线；
步骤6：通过迭代扩展频散模态分解，根据所述频域信号和关键时频脊线的初步估计值，分离重构关键频散模态及其群延迟；
步骤7：根据重构的关键频散模态，从频域信号中得到残差信号，对残差信号重复依次执行步骤5至步骤7，重复获取残差信号以及进行关键频散模态的分离重构，直至满足所述停止迭代准则；
步骤8：对重构的所有关键频散模态进行快速傅里叶变换逆变换，得到频散波的时域波形。
2.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤2具体包括以下步骤：
S201：对频散波时域信号S(t)进行快速傅里叶变换，得到对应的频域信号S(f)；
S202：通过频散波时域信号的长度Nt，确定频域信号长度Nf，该频域信号长度Nf的计算表达式为：
Nf＝floor(Nt/2)+1
式中，floor(·)为向下取整符号；
S203：保留频域信号的前Nf项，第Nf+1项至第Nt项置零。
3.根据权利要求2所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤3具体包括以下步骤：
S301：设置窗长度Nwin，确定傅里叶点数为2Nf-1，构造相应窗函数W(t)；
S302：通过所述窗函数W(t)对时域频散波信号S(t)进行离散短时傅里叶变换，获得时频表达式TFR(t,f)。
4.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤5包括以下步骤：
步骤501：进行参数初始化；
步骤502：定位所述时频表达式TFR(t,f)中能量最大值点(tem,fem)并沿频率方向双向搜索Δt范围内时间维度的能量极大值点；
步骤503：获得时频脊线作为关键模态群延时的初始估计值
Figure FDA0003274874070000021
5.根据权利要求4所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤501中，所述参数为最大时间允许波动范围Δt；
步骤502中，所述能量最大值点(tem,fem)的计算表达式为：
Figure FDA0003274874070000022
式中，t为时间，f为幅值；
所述沿频率方向双向搜索Δt范围内时间维度的能量极大值点的计算表达式为：
tR＝tem,tL＝tem
Figure FDA0003274874070000023
式中，Nf为频域信号长度。
6.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤6具体包括以下步骤：
步骤601：初始化参数，获取频域频散波信号S(f)和关键模态群延时初始值
Figure FDA0003274874070000024
步骤602：构造关键频散模态
Figure FDA0003274874070000025
并建立解析幅值
Figure FDA0003274874070000026
的冗余傅里叶展开式；
步骤603：将所述
Figure FDA0003274874070000027
的冗余傅里叶展开式带入关键频散模态
Figure FDA0003274874070000028
中构造矩阵表达式，通过求解最优化问题获得傅里叶参数矩阵
Figure FDA0003274874070000029
并重构得
Figure FDA00032748740700000210
步骤604：利用重构幅值
Figure FDA00032748740700000211
通过解卷绕求导方法估计群延时误差
Figure FDA00032748740700000212
对其进行滤波光滑然后更新群延时
Figure FDA00032748740700000213
步骤605：将
Figure FDA00032748740700000214
代入步骤603更新核矩阵Kk+1、参数矩阵
Figure FDA00032748740700000215
和解析幅值矩阵
Figure FDA00032748740700000216
重复计算达到预设的迭代次数后停止，然后求解重构关键频散模态skey(f)；
步骤606：更新频散波信号成分数量m＝m+1。
7.根据权利要求6所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤601中，所述参数包括正则化参数λ，惩罚参数β,μ和傅里叶展开阶数L；
步骤602中，所述关键频散模态
Figure FDA0003274874070000031
的计算表达式为：
Figure FDA0003274874070000032
其中，k＝0,...,end为迭代次数，
Figure FDA0003274874070000033
为解析幅值：
Figure FDA0003274874070000034
其中，
Figure FDA0003274874070000035
为初始相位；
解析幅值
Figure FDA0003274874070000036
的冗余傅里叶展开式为：
Figure FDA0003274874070000037
其中，
Figure FDA0003274874070000038
fs为时域信号的采样频率；
步骤603中，构造的矩阵表达式为：
Figure FDA0003274874070000039
Figure FDA00032748740700000310
Kk＝Ek·T为Nf×(2L+1)的核矩阵，其中，
Figure FDA00032748740700000311
Figure FDA00032748740700000312
所述最优化问题为：
Figure FDA00032748740700000313
获得的所述傅里叶参数矩阵
Figure FDA00032748740700000314
的计算表达式为：
Figure FDA00032748740700000315
重构幅值
Figure FDA00032748740700000316
的计算表达式为：
Figure FDA00032748740700000317
式中，I为单位矩阵，(·)H为矩阵的共轭转置；
步骤604中，取
Figure FDA00032748740700000318
的相位值并进行解卷绕，通过对解卷绕相位进行求导得到群延时误差
Figure FDA00032748740700000319
Figure FDA0003274874070000041
式中，∠·为求复数的相位运算，unwrap(·)为相位解卷绕运算，将相位值由[-π,π]区间叠加到[0,+∞)区间；
对
Figure FDA0003274874070000042
进行滤波光滑并更新瞬时频率
Figure FDA0003274874070000043
该步骤的计算表达式为：
Figure FDA0003274874070000044
Figure FDA0003274874070000045
其中
Figure FDA0003274874070000046
为改进二阶差分算子；
步骤605中，所述关键频散模态skey(f)的重构表达式为：
Figure FDA0003274874070000047
式中，end为预设的迭代次数，Kend为更新end次后的核矩阵，
Figure FDA0003274874070000048
为更新end次后的参数矩阵。
8.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤7包括以下步骤：
步骤701：将所述重构关键频散模态skey(f)作为第m个频散模态分量sm(f)，并将其从频散波频域信号S(f)中提取，得到残差信号Sr(f)：
Sr(f)＝S(f)-sm(f)
步骤702：计算残差信号能量与频散波频域信号能量之比：
Figure FDA0003274874070000049
步骤703：当所述能量比R大于预先设定的迭代停止判别值ε时，重复依次执行步骤5、步骤6和步骤7；否则，得到频散模态分量总数M＝m，执行步骤8。
9.根据权利要求1所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，所述步骤8包括以下步骤：
步骤801：扩展所有重构的频散模态分量sm(f)；
步骤802：对所有扩展后的频散模态分量进行快速傅里叶变换逆变换，得到频散波的时域波形。
10.根据权利要求9所述的一种基于迭代扩展频散模态分解的频散波提取方法，其特征在于，步骤801中，扩展后的频散模态分量
Figure FDA0003274874070000051
的计算表达式为：
Figure FDA0003274874070000052
式中，conj[·]表示求复数的共轭；
步骤802中，取快速傅里叶变换逆变换计算结果的实部作为对应的频散波时域波形
Figure FDA0003274874070000053
该频散波时域波形
Figure FDA0003274874070000054
的计算表达式为：
Figure FDA0003274874070000055
式中，IFFT(·)为快速傅里叶变换逆变换操作，real[·]为取复数实部操作。
CN202111114136.1A 2021-09-23 2021-09-23 一种基于迭代扩展频散模态分解的频散波提取方法 Active CN113887360B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111114136.1A CN113887360B (zh) 2021-09-23 2021-09-23 一种基于迭代扩展频散模态分解的频散波提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111114136.1A CN113887360B (zh) 2021-09-23 2021-09-23 一种基于迭代扩展频散模态分解的频散波提取方法

Publications (2)

Publication Number Publication Date
CN113887360A true CN113887360A (zh) 2022-01-04
CN113887360B CN113887360B (zh) 2024-05-31

Family

ID=79010261

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111114136.1A Active CN113887360B (zh) 2021-09-23 2021-09-23 一种基于迭代扩展频散模态分解的频散波提取方法

Country Status (1)

Country Link
CN (1) CN113887360B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114674410A (zh) * 2022-03-24 2022-06-28 安徽大学 一种分量数时变的水声信号瞬时频率估计方法
CN114866159A (zh) * 2022-04-01 2022-08-05 华南理工大学 一种多分量线性调频信号时频分析方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100157731A1 (en) * 2008-12-22 2010-06-24 Shuchin Aeron Automatic dispersion extraction of multiple time overlapped acoustic signals
US20130272548A1 (en) * 2012-04-13 2013-10-17 Qualcomm Incorporated Object recognition using multi-modal matching scheme
CN103413134A (zh) * 2013-07-11 2013-11-27 四川大学 基于稀疏分解的地面移动目标微动信号特征提取
CN104568435A (zh) * 2015-01-23 2015-04-29 国电联合动力技术有限公司 一种基于时频能量谱提取风机状态信号调制成分的方法
CN111543946A (zh) * 2020-05-08 2020-08-18 南京邮电大学 基于改进变分模态分解算法的癫痫脑电信号自动检测方法
CN113125179A (zh) * 2021-03-11 2021-07-16 同济大学 一种旋转机械转速波动无键相阶次跟踪方法
CN113295413A (zh) * 2021-06-24 2021-08-24 北京交通大学 一种基于间接信号的牵引电机轴承故障诊断方法
CN113405823A (zh) * 2021-05-17 2021-09-17 同济大学 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100157731A1 (en) * 2008-12-22 2010-06-24 Shuchin Aeron Automatic dispersion extraction of multiple time overlapped acoustic signals
US20130272548A1 (en) * 2012-04-13 2013-10-17 Qualcomm Incorporated Object recognition using multi-modal matching scheme
CN103413134A (zh) * 2013-07-11 2013-11-27 四川大学 基于稀疏分解的地面移动目标微动信号特征提取
CN104568435A (zh) * 2015-01-23 2015-04-29 国电联合动力技术有限公司 一种基于时频能量谱提取风机状态信号调制成分的方法
CN111543946A (zh) * 2020-05-08 2020-08-18 南京邮电大学 基于改进变分模态分解算法的癫痫脑电信号自动检测方法
CN113125179A (zh) * 2021-03-11 2021-07-16 同济大学 一种旋转机械转速波动无键相阶次跟踪方法
CN113405823A (zh) * 2021-05-17 2021-09-17 同济大学 一种基于迭代扩展本征模态分解的旋转机械故障诊断方法
CN113295413A (zh) * 2021-06-24 2021-08-24 北京交通大学 一种基于间接信号的牵引电机轴承故障诊断方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
XUEJIE CAO等: "research on FDIR of Key sensors for single electromagnet suspension system", 2020 GLOBAL RELIABILITY AND PROGNOSTICS AND HEALTH MANAGEMENT, 18 October 2020 (2020-10-18), pages 1 - 10 *
赵立杰;邹世达;郭烁;黄明忠;: "基于正则化随机配置网络的球磨机工况识别", 控制工程, vol. 27, no. 01, 20 January 2020 (2020-01-20), pages 1 - 7 *
郝国成;李飞;白雨晓;王巍;: "基于NDSST的非平稳信号时频分析算法", 武汉大学学报(信息科学版), vol. 44, no. 06, 5 June 2019 (2019-06-05), pages 941 - 948 *
陈洪磊: "参数化兰姆波检测信号处理与缺陷定位算法研究", 中国博士学位论文全文数据库参数化兰姆波检测信号处理与缺陷定位算法研究, no. 6, 15 June 2021 (2021-06-15), pages 022 - 6 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114674410A (zh) * 2022-03-24 2022-06-28 安徽大学 一种分量数时变的水声信号瞬时频率估计方法
CN114866159A (zh) * 2022-04-01 2022-08-05 华南理工大学 一种多分量线性调频信号时频分析方法
CN114866159B (zh) * 2022-04-01 2023-04-21 华南理工大学 一种多分量线性调频信号时频分析方法

Also Published As

Publication number Publication date
CN113887360B (zh) 2024-05-31

Similar Documents

Publication Publication Date Title
CN107061996B (zh) 一种供水管道泄漏检测定位方法
CN100554917C (zh) 获取系统特征函数和信号特征值的方法
Li et al. Dictionary learning and shift-invariant sparse coding denoising for controlled-source electromagnetic data combined with complementary ensemble empirical mode decomposition
CN113887360A (zh) 一种基于迭代扩展频散模态分解的频散波提取方法
CN110967599A (zh) 一种电能质量扰动检测与定位算法
Chen et al. Applications of empirical mode decomposition iin random noise attenuation of seismic data
CN106226407A (zh) 一种基于奇异谱分析的超声回波信号在线预处理方法
CN103675617A (zh) 一种用于高频局部放电信号检测的抗干扰方法
CN101251446B (zh) 基于离散分数余弦变换的碰摩声发射信号降噪方法
Pan et al. A noise reduction method of symplectic singular mode decomposition based on Lagrange multiplier
Yang et al. Research on ultrasonic signal processing algorithm based on CEEMDAN joint wavelet packet thresholding
CN111222088B (zh) 一种改进的平顶自卷积窗加权电力谐波幅值估计方法
CN104808243B (zh) 一种叠前地震贝叶斯反演方法和装置
Battistel et al. On the physical definition of dynamic impedance: how to design an optimal strategy for data extraction
CN103675544A (zh) 基于优化算法的电力系统故障信号检测与波形识别方法
Cheng et al. Enhanced symplectic characteristics mode decomposition method and its application in fault diagnosis of rolling bearing
CN105044794A (zh) 一种核磁共振回波数据的压缩方法及装置
CN109630908B (zh) 一种多次降噪的管道泄漏定位方法
CN101334482B (zh) 一种预测地震波中的多次波和一次波信号的方法
CN111505709B (zh) 一种基于稀疏谱分解的衰减定性分析的方法
Chen et al. A signal-enhancement algorithm for the quantification of NMR data in the time domain
CN106980722B (zh) 一种脉冲响应中谐波成分的检测和去除方法
CN107367760B (zh) 基于加速线性Bregman算法的表面多次波和子波估计方法及系统
CN116401513A (zh) 一种基于深度残差网络的磁共振工频谐波噪声抑制方法
Faisal et al. Suppression of false-terms in wigner-ville distribution using time and frequency windowing

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant