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

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

Info

Publication number
CN111487318A
CN111487318A CN202010478847.6A CN202010478847A CN111487318A CN 111487318 A CN111487318 A CN 111487318A CN 202010478847 A CN202010478847 A CN 202010478847A CN 111487318 A CN111487318 A CN 111487318A
Authority
CN
China
Prior art keywords
signal
frequency
wavelet
component
time
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
CN202010478847.6A
Other languages
English (en)
Other versions
CN111487318B (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.)
Fujian Agriculture and Forestry University
Original Assignee
Fujian Agriculture and Forestry 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 Fujian Agriculture and Forestry University filed Critical Fujian Agriculture and Forestry University
Priority to CN202010478847.6A priority Critical patent/CN111487318B/zh
Publication of CN111487318A publication Critical patent/CN111487318A/zh
Application granted granted Critical
Publication of CN111487318B publication Critical patent/CN111487318B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04Analysing solids
    • G01N29/045Analysing solids by imparting shocks to the workpiece and detecting the vibrations or the acoustic waves caused by the shocks
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04Analysing solids
    • G01N29/12Analysing solids by measuring frequency or resonance of acoustic waves
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4454Signal recognition, e.g. specific values or portions, signal events, signatures
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4472Mathematical theories or simulation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/46Processing the detected response signal, e.g. electronic circuits specially adapted therefor by spectral analysis, e.g. Fourier analysis or wavelet analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/02Indexing codes associated with the analysed material
    • G01N2291/023Solids
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/02Indexing codes associated with the analysed material
    • G01N2291/023Solids
    • G01N2291/0234Metals, e.g. steel
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/02Indexing codes associated with the analysed material
    • G01N2291/028Material parameters
    • G01N2291/0289Internal structure, e.g. defects, grain size, texture

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Pathology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biochemistry (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Mathematical Physics (AREA)
  • Acoustics & Sound (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

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

Description

一种时变结构瞬时频率提取方法
技术领域
本发明属于结构健康监测和安全评估领域,具体涉及一种时变结构瞬时频率提取方法。
背景技术
实际桥梁结构在环境侵蚀、材料老化、工作荷载、疲劳效应等多重因素作用下不可避免地出现损伤并逐渐累积,在极端情况下甚至有可能引发突发性的灾难事故。因此,为保障桥梁结构的安全性、完整性、适用性和耐久性,减少人员伤亡和财产损失,对在役桥梁结构采取有效的手段进行结构健康监测和安全评估显得十分必要和迫切,同时也具有较大的社会价值和经济效益。
作为桥梁结构健康监测和安全评估领域的一项核心技术,参数识别属于系统识别范畴,一般包括以质量、刚度和阻尼为特征参数的物理参数识别和以模态频率、模态振型和模态阻尼为特征参数的模态参数识别。时不变结构模态参数识别方法包括频域法和时域法,它们假定输入激励为白噪声、结构响应为平稳信号,因此无法反映模态参数的时变特征,同时也无法追踪结构的损伤演化趋势。截至目前,现有模态参数识别方法主要集中在时不变结构系统,针对环境激励下的时变结构开展密集模态参数和瞬时特征参数识别研究尚不多见。
具有刚度或质量突变的结构系统的动力特征参数必将随时间发生变化,其模态频率有可能呈现密集或者叠混特性。因此,环境激励下的在役桥梁结构,其响应信号不仅是非平稳的,而且极有可能是模态密集分布和叠混的。显然,从非平稳信号处理的角度识别环境激励下桥梁结构的时变密集模态参数更符合实际情况,对深入理解结构损伤诊断、有限元模型修正、振动控制和安全评估具有重要的理论意义和工程应用价值。
发明内容
本发明的目的在于提供一种时变结构瞬时频率提取方法,该方法有利于对桥梁时变结构的瞬时频率进行提取。
为实现上述目的,本发明采用的技术方案是:一种时变结构瞬时频率提取方法,对于时变钢桥或组合桥梁结构,首先将广义解调定理与解析模态分解理论相结合,将具有密集模态甚至叠混的多分量信号分解为多个单分量信号,再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调,从而实现调幅调频信号到纯调频信号的转变,最后通过小波系数模局部极大值方法精确提取瞬时频率。
进一步地,该方法具体包括如下步骤:
步骤1)对于瞬时相位为φ(t)的渐进信号x(t),其相应的解析信号表达式xA(t)和瞬时频率f(t)分别如式(1)和(2)所示:
xA(t)=x(t)+jH[x(t)]=x(t)ej[φ(t)] (1)
Figure BDA0002516658080000021
式中,j为虚数单位,H[·]表示希尔伯特变换算子;通过解析信号xA(t)来表示信号x(t)的广义傅里叶变换,即信号x(t)的调制,如式(3)所示:
Figure BDA0002516658080000022
式中,XG(ω)为调制后的信号,φ(t)为与时间相关的实相位信号;为了将信号的时变瞬时频率调制到固定频率f0,构造φ(t),如式(4)所示:
Figure BDA0002516658080000023
式中,
Figure BDA0002516658080000024
为原多分量信号各阶瞬时频率估计值;
步骤2)对调制后的信号XG(ω)进行解析模态分解,选取该信号的实部XG(t),假设XG(t)是由n个分量xi(t)组成的多分量信号:
XG(t)=x1(t)+x2(t)+…+xi(t)+…+xn(t)i=1,2,…,n (5)
若分量信号的频率ω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个二分时不变截止频率,则分量信号xi(t)通过式(6)和(7)解析得出:
x1(t)=s1(t),…,xi(t)=si(t)-si-1(t),…,xn(t)=XG(t)-sn-1(t)i=1,2,…,n(6)
si(t)=sin[ωbi(t)]H{XG(t)cos[ωbi(t)]}-cos[ωbi(t)]H{XG(t)sin[ωbi(t)]}
i=1,2,…,n-1 (7)
步骤3)对调制后的信号XG(ω)进行逆广义傅里叶变换,即重构,可得:
Figure BDA0002516658080000025
步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调,以得到纯调频信号;对于一个分解得到的单分量信号xi(t),其希尔伯特变换通过式(9)表示:
Figure BDA0002516658080000031
式中,P为柯西积分主值;
步骤5)对于不满足Bedrosian定理的单分量信号,希尔伯特解调存在较大误差,通过迭代希尔伯特解决这一问题;通过所述单分量信号x1(t)构造出如式(10)所示的解析信号z(t),z(t)的实部为原始信号,虚部为原始信号的希尔伯特变换:
Figure BDA0002516658080000032
式中,
Figure BDA0002516658080000033
和φ1=tan-1{H[xi(t)]/xi(t)}分别为幅值函数和调频函数;解析信号的实部表示为幅值函数A1和调频函数cosφ1的乘积,如式(11)所示:
Figure BDA0002516658080000034
Figure BDA0002516658080000035
表示第一次迭代希尔伯特变换后的原始分量信号。将生成的调频信号
Figure BDA0002516658080000036
Figure BDA0002516658080000037
作为新的信号并继续对其进行希尔伯特变换,产生新的幅值函数A2和调频函数cosφ2,即:
Figure BDA0002516658080000038
式中,
Figure BDA0002516658080000039
表示第二次迭代希尔伯特变换后的原始分量信号,
Figure BDA00025166580800000310
Figure BDA00025166580800000311
步骤6)不断重复步骤5),得到如下迭代公式:
Figure BDA00025166580800000312
Figure BDA00025166580800000313
表示第m次迭代希尔伯特变换后的原始分量信号。由于每一次迭代过程均会产生新的调频函数及幅值函数,因此只有当新的幅值函数Am趋近于1时,迭代才会停止;此时得到的调频函数
Figure BDA00025166580800000314
即为单分量信号xi(t)的纯调频信号,也近似为渐进单分量信号;
步骤7)对经过迭代希尔伯特变换处理后的纯调频信号
Figure BDA00025166580800000315
进行连续小波变换得到如式(14)所示的小波系数;
Figure BDA00025166580800000316
式中,Wx(a,b)表示小波系数;a表示尺度因子,与频率成反比关系;b表示平移因子,与时间有关;
Figure BDA0002516658080000041
Figure BDA0002516658080000042
的共轭复数,
Figure BDA0002516658080000043
表示小波母函数;小波系数表征小波母函数与信号相似的程度,小波系数越大表示两者越相似;
步骤8)由于小波系数模值在小波脊线上取得局部极大值,对任意时刻bn,通过式(15)寻找不同尺度下小波系数模的局部极大值,从而得到对应的小波脊点(am,bn);按照同样的方法继续提取下一时刻bn+1所对应的小波脊点(am+1,bn+1),直到信号的末端终点,最后,将提取到的局部极大值点连接起来作为信号的小波脊线,即求得渐进单分量信号的瞬时频率;
|Wx(ar,b)|=max|Wx(am,bn)| (15)
其中,Wx(ar,b)表示局部极大值点。
相较于现有技术,本发明具有以下有益效果:提出了一种基于广义解调解析模态分解、迭代希尔伯特变换和小波系数模局部极大值法相结合的时变结构瞬时频率提取方法,该方法能针对时变和非线性结构系统仅利用其非平稳响应信号即可实现瞬时频率的提取。在参数识别的基础上,能进一步构建损伤位置及时变损伤指标,可用于设备故障检测、结构健康监测和安全评估领域。该方法仅需测得响应信号即可实现参数提取,无需测量输入激励,具有很强的实用性和广阔的应用前景。
附图说明
图1为本发明方法的实现流程图。
图2为本发明实施例中时变悬臂梁试验装置示意图。
图3为本发明实施例中原多分量信号x(t)调制前后时频图。
图4为本发明实施例中分量信号x1(t),x2(t)识别值与理论值对比图。
图5为本发明实施例中分量信号x1(t),x2(t)解调前后对比图。
图6为本发明实施例中分量信号x1(t),x2(t)瞬时频率识别值与理论值对比图。
图7为本发明实施例中信号原多分量信号y(t)调制前后时频图。
图8为本发明实施例中分量信号y1(t),y2(t)识别值与理论值对比图。
图9为本发明实施例中分量信号y1(t),y2(t)解调前后对比图。
图10为本发明实施例中分量信号y1(t),y2(t)瞬时频率识别值与理论值对比图。
图11为本发明实施例中信号原多分量信号z(t)调制前后时频图。
图12为本发明实施例中分量信号z1(t),z2(t)识别值与理论值对比图。
图13为本发明实施例中分量信号z1(t),z2(t)解调前后对比图。
图14为本发明实施例中分量信号z1(t),z2(t)瞬时频率识别值与理论值对比图。
图15为本发明实施例中信号原多分量信号s(t)调制前后时频图。
图16为本发明实施例中分量信号s1(t),s2(t)识别值与理论值对比图。
图17为本发明实施例中分量信号s1(t),s2(t)解调前后对比图。
图18为本发明实施例中分量信号s1(t),s2(t)瞬时频率识别值与理论值对比图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步的详细说明。
本发明提供了一种时变结构瞬时频率提取方法,对于时变钢桥或组合桥梁结构,首先将广义解调定理与解析模态分解理论相结合,将具有密集模态甚至叠混的多分量信号分解为多个单分量信号,由于分解所得单分量信号通常不是纯调频信号,再通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调,从而实现调幅调频信号到纯调频信号的转变,最后通过小波系数模局部极大值方法精确提取瞬时频率。具体包括如下步骤:
步骤1)对于瞬时相位为φ(t)的渐进信号x(t),其相应的解析信号表达式xA(t)和瞬时频率f(t)分别如式(1)和(2)所示:
xA(t)=x(t)+jH[x(t)]=x(t)ej[φ(t)] (I)
Figure BDA0002516658080000051
式中,j为虚数单位,H[·]表示希尔伯特变换算子;通过解析信号xA(t)来表示信号x(t)的广义傅里叶变换,即信号x(t)的调制,如式(3)所示:
Figure BDA0002516658080000052
式中,XG(ω)为调制后的信号,φ(t)为与时间相关的实相位信号;为了将信号的时变瞬时频率调制到固定频率f0,构造φ(t),如式(4)所示:
Figure BDA0002516658080000053
式中,
Figure BDA0002516658080000054
为原多分量信号各阶瞬时频率估计值,即通过其他算法得到瞬时频率大概的估计值,可通过同步挤压小波变换或其他瞬时频率识别方法得到;与理论值相比,虽然识别值
Figure BDA0002516658080000055
存在一定的误差,但经过调制后的信号的瞬时频率总体上仍然会集中在固定频率f0附近;
步骤2)对调制后的信号XG(ω)进行解析模态分解,选取该信号的实部XG(t),假设XG(t)是由n个分量xi(t)组成的多分量信号:
XG(t)=x1(t)+x2(t)+…+xi(t)+…+xn(t) i=1,2,…,n (5)
若分量信号的频率ω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个二分时不变截止频率,则分量信号xi(t)通过式(6)和(7)解析得出:
x1(t)=s1(t),…,xi(t)=si(t)-si-1(t),…,xn(t)=XG(t)-sn-1(t) i=1,2,…,n(6)
si(t)=sin[ωbi(t)]H{XG(t)cos[ωbi(t)]}-cos[ωbi(t)]H{XG(t)sin[ωbi(t)]}
i=1,2,…,n-1 (7)
步骤3)对调制后的信号XG(ω)进行逆广义傅里叶变换,即重构,可得:
Figure BDA0002516658080000061
步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调,以得到纯调频信号;对于一个分解得到的单分量信号xi(t),其希尔伯特变换通过式(9)表示:
Figure BDA0002516658080000062
式中,P为柯西积分主值;
步骤5)对于不满足Bedrosian定理的单分量信号,希尔伯特解调存在较大误差,通过迭代希尔伯特解决这一问题;通过所述单分量信号xi(t)构造出如式(10)所示的解析信号z(t),z(t)的实部为原始信号,虚部为原始信号的希尔伯特变换:
Figure BDA0002516658080000063
式中,
Figure BDA0002516658080000064
和φ1=tan-1{H[xi(t)]/xi(t)}分别为幅值函数和调频函数;解析信号的实部表示为幅值函数A1和调频函数cosφ1的乘积,如式(11)所示:
Figure BDA0002516658080000065
Figure BDA0002516658080000066
表示第一次迭代希尔伯特变换后的原始分量信号。将生成的调频信号
Figure BDA0002516658080000067
Figure BDA0002516658080000068
作为新的信号并继续对其进行希尔伯特变换,产生新的幅值函数A2和调频函数cosφ2,即:
Figure BDA0002516658080000069
式中,
Figure BDA00025166580800000610
表示第二次迭代希尔伯特变换后的原始分量信号,
Figure BDA00025166580800000611
Figure BDA0002516658080000071
步骤6)不断重复步骤5),得到如下迭代公式:
Figure BDA0002516658080000072
Figure BDA0002516658080000073
表示第m次迭代希尔伯特变换后的原始分量信号。由于每一次迭代过程均会产生新的调频函数及幅值函数,因此只有当新的幅值函数Am趋近于1时,迭代才会停止;此时得到的调频函数
Figure BDA0002516658080000074
即为单分量信号xi(t)的纯调频信号,也近似为渐进单分量信号;
步骤6)不断重复步骤5),得到如下迭代公式:
Figure BDA0002516658080000075
由于每一次迭代过程均会产生新的调频函数及幅值函数,因此只有当新的幅值函数An趋近于1时,迭代才会停止;此时得到的调频函数xn(t)=cosφn即为单分量信号x1(t)的纯调频信号,也近似为渐进单分量信号;
步骤7)对经过迭代希尔伯特变换处理后的纯调频信号x(t)进行连续小波变换得到如式(14)所示的小波系数;
Figure BDA0002516658080000076
式中,Wx(a,b)表示小波系数;a表示尺度因子,与频率成反比关系;b表示平移因子,与时间有关;
Figure BDA0002516658080000077
Figure BDA0002516658080000078
的共轭复数,
Figure BDA0002516658080000079
表示小波母函数;小波系数表征小波母函数与信号相似的程度,小波系数越大表示两者越相似;
步骤8)由于小波系数模值在小波脊线上取得局部极大值,对任意时刻bn,通过式(15)寻找不同尺度下小波系数模的局部极大值,从而得到对应的小波脊点(am,bn);按照同样的方法继续提取下一时刻bn+1所对应的小波脊点(am+1,bn+1),直到信号的末端终点,最后,将提取到的局部极大值点连接起来作为信号的小波脊线,即求得渐进单分量信号的瞬时频率;
|Wx(ar,b)|=max|Wx(am,bn)| (15)
其中,Wx(ar,b)表示局部极大值点。
下面举具体实施例对本发明进行说明。
1)现有多分量调幅调频信号x(t)=x1(t)+x2(t)=(t+6)cos[2π(t2+12t)]+(3t+18)cos[2π(0.6t2+7t)]。设定采样频率fs=100Hz,采样时间t=10s,采样间隔dt=0.01s,考虑到实际工程中有噪声干扰,因此添加20%高斯白噪声(7dB),噪声强度由信噪比定义:
Figure BDA0002516658080000081
式中,Asignal和Anoise分别代表信号和噪声的均方根值,噪声水平是指
Figure BDA0002516658080000082
Figure BDA0002516658080000083
之间的比值。
首先选用复Morlet小波对添加了20%高斯白噪声原信号进行同步挤压小波变换变换得到原信号的时频图,如图3(a)所示。从图中可以看出:该多分量信号包括了两个单分量信号,且两个分量信号的频率随时间呈线性变化,不仅是时变而且也是模态叠混的,因此AMD无法分解该信号。
采用广义解调解析模态分解定理对该信号进行处理。首先对原信号进行同步挤压小波变换得到原多分量信号的各阶瞬时频率的估计值,然后根据估计值选取合适的固定频率f0=4Hz对原多分量信号进行信号调制,调制后的时频图如图3(b)所示。从图3(b)中可以看出:分量信号x2(t)经调制后的瞬时频率在4Hz左右浮动,而分量信号x1(t)的频率大概介于(10Hz,20Hz)之间。此时两分量信号之间未出现模态叠混,因此选择截止频率4.1Hz对调制后的多分量信号进行AMD分解,即可成功对信号进行分离。然后进行逆广义傅里叶变换便后可得到原分量信号x1(t)和x2(t),两分量信号的分解结果分别如图4(a)和(b)所示。
紧接着对得到的单分量调幅调频信号进行迭代希尔伯特变换解调,解调结果如图5(a)和(b)所示,从图中可以看出分量信号在解调后的幅值已接近1,可近似认为渐进单分量信号。
最后通过小波系数局部模极大值法对得到的近似渐进单分量信号进行瞬时频率提取,结果如图6(a)和(b)所示。由图6(a)和(b)可以看出:瞬时频率识别值与理论值吻合得很好该算法的正确性得到了验证。
2)考虑多分量调幅调频信号y(t)=y1(t)+y2(t)=(t+9)cos(10πt+2πt2)+(2t+1)cos[20πt+3πt2+sin(2πt)],已添加20%高斯白噪声(7dB),由信号表达式可知该多分量信号由两个调频调幅信号组成。设定采样频率fs=100Hz,采样时间t=10s,采样间隔dt=0.01s。
首先选用复Morlet小波对该信号进行同步挤压小波变换得到其时频图,如图7(a)所示。从图7(a)中可获知:该多分量信号由两个时变信号组成,分别为时变正弦信号和一个时变线性信号。由于两分量信号出现模态叠混,传统AMD无法分离模态叠混信号。
采用广义解调解析模态分解定理对该信号进行处理。首先对原信号进行同步挤压小波变换得到原多分量信号的各阶瞬时频率的估计值,然后根据估计值选取合适的固定频率f0=4Hz对原多分量信号进行信号调制,调制以后时频图如图7(b)所示:分量信号y1(t)经调制以后,瞬时频率在4Hz左右浮动,而分量信号y2(t)经调制以后,瞬时频率变化范围大致在(10Hz,20Hz)之间呈正弦规律变化,而且两分量信号未模态叠混。选取二分截止频率5Hz,通过AMD对信号进行分解。随后对得到的分量信号进行逆广义傅里叶变换便可得到原调幅调频分量信号y1(t)和y2(t)。两分量信号的分解结果如图8(a)和(b)所示。由于添加了白噪声,信号在初始阶段吻合度不高,但是在其他阶段吻合较好。
由于分解得到的分量信号并不是纯调频信号,因此通过迭代希尔伯特变换对两分量信号进行解调,结果如图9(a)和(b)所示,分量信号在解调后的幅值已接近1,可近似认为渐进单分量信号。
最后采用小波系数局部模极大值对解调后所得的近似纯调频信号进行瞬时频率提取,结果如图10(a)和(b)所示。由于信号中添加了20%白噪声,初始阶段识别值不可避免出现了端点效应,但整体和理论值还是比较吻合。
3)继续考虑多分量调幅调频信号z(t)=z1(t)+z2(t)=(3t+1)cos(18πt+4πt2)+(2t+2)cos[2πt2+sin(2πt)+8πt],同样添加20%高斯白噪声(7dB)。设定采样频率fs=100Hz,采样时间t=10,采样间隔dt=0.01s。
首先选用复Morlet小波对原添加白噪声的多分量信号进行同步挤压小波变换得到其小波量图,如图11(a)所示。从图中可以看出:原多分量信号包含了两个时变信号,其中一个是呈正弦变化,一个呈线性变化,且两分量信号出现叠混。
通过同步挤压小波变换得到分量信号的瞬时频率估计值以后,根据瞬时频率估计值选择固定频率f0=4Hz对原信号进行调制,调制后的效果如图11(b)所示:单分量信号z2(t)调制后,其瞬时频率在4Hz左右浮动,而单分量信号z1(t)的频率大致介于(9Hz,30Hz)之间。选定截止频率为4.5Hz,采用AMD对信号进行分解,并对分解的信号进行逆广义傅里叶变换得到原分量信号z1(t)和z2(t),如图12(a)和(b)所示。
在得到两个调幅调频信号以后,引入迭代希尔伯特变换对单分量信号进行解调,解调过程结束后最终得到两个近似纯调频分量信号,结果如图13(a)和(b)所示。
最后通过小波系数局部模极大值对分量信号的瞬时频率进行提取,结果如图14(a)和(b)所示,同样地除了不可避免的端点效应以外,整体上和理论值也是很吻合的。
4)继续考虑多分量调幅调频信号s(t)=s1(t)+s2(t)=(3t+1)cos[πt2+sin(2πt)+12πt]+(2t)cos[2πt2+sin(2πt)+30πt],已添加20%高斯白噪声。设定采样频率fs=100Hz,采样时间t=10s,采样间隔dt=0.01s。
首先选择复Morlet小波对该信号做同步挤压小波变换,得到小波量图,如图15(a)所示。由图15(a)可知:原多分量信号包含了两个正弦规律变化的时变分量信号,且两个信号出现了模态叠混,因此AMD无法直接处理该信号。
首先通过同步挤压小波变换得到信号的瞬时频率估计值,然后通过估计值确定固定频率f0=3Hz,紧接着对原多分量信号进行调制,调制以后的信号时频图如图15(b)所示。由图15(b)可知:原来两条正弦规律变化的信号经过调制以后均成了近似线性信号,其中信号s1(t)的瞬时频率在3Hz左右浮动,而信号s2(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),其相应的解析信号表达式xA(t)和瞬时频率f(t)分别如式(1)和(2)所示:
xA(t)=x(t)+jH[x(t)]=x(t)ej[φ(t)] (1)
Figure FDA0002516658070000011
式中,j为虚数单位,H[·]表示希尔伯特变换算子;通过解析信号xA(t)来表示信号x(t)的广义傅里叶变换,即信号x(t)的调制,如式(3)所示:
Figure FDA0002516658070000012
式中,XG(ω)为调制后的信号,φ(t)为与时间相关的实相位信号;为了将信号的时变瞬时频率调制到固定频率f0,构造φ(t),如式(4)所示:
Figure FDA0002516658070000013
式中,
Figure FDA0002516658070000014
为原多分量信号各阶瞬时频率估计值;
步骤2)对调制后的信号XG(ω)进行解析模态分解,选取该信号的实部XG(t),假设XG(t)是由n个分量xi(t)组成的多分量信号:
XG(t)=x1(t)+x2(t)+…+xi(t)+…+xn(t)i=1,2,…,n (5)
若分量信号的频率ω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个二分时不变截止频率,则分量信号xi(t)通过式(6)和(7)解析得出:
x1(t)=s1(t),…,xi(t)=si(t)-si-1(t),…,xn(t)=XG(t)-sn-1(t)i=1,2,…,n (6)
si(t)=sin[ωbi(t)]H{XG(t)cos[ωbi(t)]}-cos[ωbi(t)]H{XG(t)sin[ωbi(t)]}i=1,2,…,n-1 (7)
步骤3)对调制后的信号XG(ω)进行逆广义傅里叶变换,即重构,可得:
Figure FDA0002516658070000021
步骤4)通过迭代希尔伯特变换对分解得到的调幅调频信号进行解调,以得到纯调频信号;对于一个分解得到的单分量信号xi(t),其希尔伯特变换通过式(9)表示:
Figure FDA0002516658070000022
式中,P为柯西积分主值;
步骤5)对于不满足Bedrosian定理的单分量信号,希尔伯特解调存在较大误差,通过迭代希尔伯特解决这一问题;通过所述单分量信号xi(t)构造出如式(10)所示的解析信号z(t),z(t)的实部为原始信号,虚部为原始信号的希尔伯特变换:
Figure FDA0002516658070000023
式中,
Figure FDA0002516658070000024
和φ1=tan-1{H[xi(t)]/xi(t)}分别为幅值函数和调频函数;解析信号的实部表示为幅值函数A1和调频函数cosφ1的乘积,如式(11)所示:
Figure FDA0002516658070000025
Figure FDA0002516658070000026
表示第一次迭代希尔伯特变换后的原始分量信号;将生成的调频信号
Figure FDA0002516658070000027
Figure FDA0002516658070000028
作为新的信号并继续对其进行希尔伯特变换,产生新的幅值函数A2和调频函数cosφ2,即:
Figure FDA0002516658070000029
式中,
Figure FDA00025166580700000210
表示第二次迭代希尔伯特变换后的原始分量信号,
Figure FDA00025166580700000211
Figure FDA00025166580700000212
步骤6)不断重复步骤5),得到如下迭代公式:
Figure FDA00025166580700000213
Figure FDA00025166580700000214
表示第m次迭代希尔伯特变换后的原始分量信号;由于每一次迭代过程均会产生新的调频函数及幅值函数,因此只有当新的幅值函数Am趋近于1时,迭代才会停止;此时得到的调频函数
Figure FDA00025166580700000215
即为单分量信号xi(t)的纯调频信号,也近似为渐进单分量信号;
步骤7)对经过迭代希尔伯特变换处理后的纯调频信号
Figure FDA00025166580700000216
进行连续小波变换得到如式(14)所示的小波系数;
Figure FDA0002516658070000031
式中,Wx(a,b)表示小波系数;a表示尺度因子,与频率成反比关系;b表示平移因子,与时间有关;
Figure FDA0002516658070000032
Figure FDA0002516658070000033
的共轭复数,
Figure FDA0002516658070000034
表示小波母函数;小波系数表征小波母函数与信号相似的程度,小波系数越大表示两者越相似;
步骤8)由于小波系数模值在小波脊线上取得局部极大值,对任意时刻bn,通过式(15)寻找不同尺度下小波系数模的局部极大值,从而得到对应的小波脊点(am,bn);按照同样的方法继续提取下一时刻bn+1所对应的小波脊点(am+1,bn+1),直到信号的末端终点,最后,将提取到的局部极大值点连接起来作为信号的小波脊线,即求得渐进单分量信号的瞬时频率;
|Wx(ar,b)|=max|Wx(am,bn)| (15)
其中,Wx(ar,b)表示局部极大值点。
CN202010478847.6A 2020-05-29 2020-05-29 一种时变结构瞬时频率提取方法 Active CN111487318B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010478847.6A CN111487318B (zh) 2020-05-29 2020-05-29 一种时变结构瞬时频率提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010478847.6A CN111487318B (zh) 2020-05-29 2020-05-29 一种时变结构瞬时频率提取方法

Publications (2)

Publication Number Publication Date
CN111487318A true CN111487318A (zh) 2020-08-04
CN111487318B CN111487318B (zh) 2023-03-24

Family

ID=71811405

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010478847.6A Active CN111487318B (zh) 2020-05-29 2020-05-29 一种时变结构瞬时频率提取方法

Country Status (1)

Country Link
CN (1) CN111487318B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112069945A (zh) * 2020-08-25 2020-12-11 大连理工大学 工程结构时变频率和阻尼比的一种识别方法
CN114880627A (zh) * 2022-07-08 2022-08-09 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN115687862A (zh) * 2022-10-19 2023-02-03 北京科技大学 一种基于时变滤波的旋转机械信号时频分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1819572A (zh) * 2006-03-23 2006-08-16 上海交通大学 基于希尔伯特-黄变换的二进制频移键控系统解调方法
CN103455470A (zh) * 2013-09-03 2013-12-18 上海交通大学 一种瞬时频率含交叉点的信号时频分解方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1819572A (zh) * 2006-03-23 2006-08-16 上海交通大学 基于希尔伯特-黄变换的二进制频移键控系统解调方法
CN103455470A (zh) * 2013-09-03 2013-12-18 上海交通大学 一种瞬时频率含交叉点的信号时频分解方法

Non-Patent Citations (16)

* Cited by examiner, † Cited by third party
Title
HENG LIU ET AL.: "Moire´ interferogram phase extraction: a ridge detection algorithm for continuous wavelet transforms", 《APPLIED OPTICS》 *
JINGLIANG LIU ET AL.: "A combined method for instantaneous frequency identification in low frequency structures", 《ENGINEERING STRUCTURES》 *
JING-LIANG LIU ET AL.: "A new instantaneous frequency extraction method for nonstationary response signals in civil engineering structures", 《JOURNAL OF LOW FREQUENCY NOISE,VIBRATION AND ACTIVE CONTROL 》 *
ZUO-CAI WANG ET AL.: "Time-Varying Linear and Nonlinear Structural Identification with Analytical Mode Decomposition and Hilbert Transform", 《J. STRUCT. ENG》 *
刘景良等: "变分模态分解和同步挤压小波变换识别时变结构瞬时频率", 《振动与冲击》 *
刘景良等: "基于改进同步挤压小波变换识别信号瞬时频率", 《振动.测试与诊断》 *
刘景良等: "基于最大坡度法提取非平稳信号小波脊线和瞬时频率", 《工程力学》 *
戴豪民等: "瞬时频率计算方法的比较研究和改进", 《四川大学学报(自然科学版)》 *
杨仁树等: "基于CEEMD与TQWT组合方法的爆破振动信号精细化特征提取", 《振动与冲击》 *
汪赵华等: "基于改进小波脊提取算法的数字信号瞬时频率估计方法", 《中国科学院研究生院学报》 *
王超等: "基于动态规划提取信号小波脊和瞬时频率", 《中南大学学报(自然科学版)》 *
王超等: "基于复小波变换的结构瞬时频率识别", 《振动工程学报》 *
程军圣: "基于广义解调时频分析的多分量信号分解方法", 《振动工程学报》 *
胡志祥等: "基于递归希尔伯特变换的振动信号解调和瞬时频率计算方法", 《振动与冲击》 *
郑近德: "非平稳信号分析的广义解析模态分解方法", 《电子学报》 *
黎海妹等: "一种改进的小波脊提取算法", 《通信技术》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112069945A (zh) * 2020-08-25 2020-12-11 大连理工大学 工程结构时变频率和阻尼比的一种识别方法
CN112069945B (zh) * 2020-08-25 2024-01-05 大连理工大学 工程结构时变频率和阻尼比的一种识别方法
CN114880627A (zh) * 2022-07-08 2022-08-09 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN114880627B (zh) * 2022-07-08 2022-09-23 西南交通大学 一种基于迭代解调时变滤波的自适应瞬时频率估计方法
CN115687862A (zh) * 2022-10-19 2023-02-03 北京科技大学 一种基于时变滤波的旋转机械信号时频分析方法
CN115687862B (zh) * 2022-10-19 2023-08-01 北京科技大学 一种基于时变滤波的旋转机械信号时频分析方法

Also Published As

Publication number Publication date
CN111487318B (zh) 2023-03-24

Similar Documents

Publication Publication Date Title
CN111487318B (zh) 一种时变结构瞬时频率提取方法
Yan et al. Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis
Zhang et al. Damage detection method based on operating deflection shape curvature extracted from dynamic response of a passing vehicle
Wang et al. Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis
Li et al. Time-varying characteristics of bridges under the passage of vehicles using synchroextracting transform
Huang et al. A method for tachometer-free and resampling-free bearing fault diagnostics under time-varying speed conditions
CN101586997A (zh) 一种拉索振动基频的计算方法
CN105787655B (zh) 超高层结构模态参数识别方法
Zhang et al. Variational mode decomposition based modal parameter identification in civil engineering
Dion et al. Harmonic component detection: Optimized Spectral Kurtosis for operational modal analysis
CN105928702B (zh) 基于形态分量分析的变工况齿轮箱轴承故障诊断方法
Rodopoulos et al. A parametric approach for the estimation of the instantaneous speed of rotating machinery
Le et al. Distinction between harmonic and structural components in ambient excitation tests using the time–frequency domain decomposition technique
CN104165742A (zh) 一种基于互谱函数的运行模态分析实验方法及装置
CN104050147A (zh) 将时域信号转换成频域信号的方法与系统
Wang et al. Weak fault diagnosis of rolling bearing under variable speed condition using IEWT-based enhanced envelope order spectrum
Liu et al. Damage-sensitive and domain-invariant feature extraction for vehicle-vibration-based bridge health monitoring
Zhang et al. Improved local cepstrum and its applications for gearbox and rolling bearing fault detection
US8995230B2 (en) Method of extracting zero crossing data from full spectrum signals
CN106980722B (zh) 一种脉冲响应中谐波成分的检测和去除方法
CN117470484A (zh) 刹车片模态频率测量方法及测量系统
Li et al. Modal parameter identification for closely spaced modes using an Empirical Fourier decomposition-based method
Al-hababi et al. Time-frequency domain methods for the identification of breathing cracks in beam-like structures
Kandula et al. Field testing of indirect displacement estimation using accelerometers
Kumar et al. Damage Identification of Beam Structure Using Discrete wavelet transform

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