CN111353400B - 一种基于视觉测振的全场景振动强度图谱分析方法 - Google Patents
一种基于视觉测振的全场景振动强度图谱分析方法 Download PDFInfo
- Publication number
- CN111353400B CN111353400B CN202010111389.2A CN202010111389A CN111353400B CN 111353400 B CN111353400 B CN 111353400B CN 202010111389 A CN202010111389 A CN 202010111389A CN 111353400 B CN111353400 B CN 111353400B
- Authority
- CN
- China
- Prior art keywords
- image
- expression
- filter
- phase
- 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.)
- Active
Links
- 238000005259 measurement Methods 0.000 title claims abstract description 22
- 238000000034 method Methods 0.000 title claims abstract description 21
- 230000000007 visual effect Effects 0.000 title claims abstract description 15
- 238000010183 spectrum analysis Methods 0.000 title claims abstract description 5
- 230000033001 locomotion Effects 0.000 claims abstract description 89
- 230000008859 change Effects 0.000 claims abstract description 85
- 238000012545 processing Methods 0.000 claims abstract description 21
- 238000001514 detection method Methods 0.000 claims abstract description 19
- 238000010586 diagram Methods 0.000 claims abstract description 12
- 230000014509 gene expression Effects 0.000 claims description 72
- 238000006073 displacement reaction Methods 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 21
- 230000009466 transformation Effects 0.000 claims description 18
- 238000013178 mathematical model Methods 0.000 claims description 13
- 238000000354 decomposition reaction Methods 0.000 claims description 10
- 238000009826 distribution Methods 0.000 claims description 9
- 238000001914 filtration Methods 0.000 claims description 9
- 238000000605 extraction Methods 0.000 claims description 8
- 230000010354 integration Effects 0.000 claims description 6
- 238000007794 visualization technique Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 239000003086 colorant Substances 0.000 claims description 3
- 239000000284 extract Substances 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 230000010363 phase shift Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 abstract description 7
- 230000036541 health Effects 0.000 abstract description 5
- 238000012544 monitoring process Methods 0.000 abstract description 5
- 238000009434 installation Methods 0.000 abstract description 2
- 238000005516 engineering process Methods 0.000 description 7
- 238000007689 inspection Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000002159 abnormal effect Effects 0.000 description 3
- 238000000691 measurement method Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012800 visualization Methods 0.000 description 2
- 206010010356 Congenital anomaly Diseases 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000005339 levitation Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000003362 replicative effect Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000011179 visual inspection Methods 0.000 description 1
- 230000009012 visual motion Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H9/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M13/00—Testing of machine parts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
- G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Abstract
一种基于视觉测振的全场景振动强度图谱分析方法。属于机械振动的结构健康检测的领域。本发明针对以上问题，提出了一种具备亚像素级振动测量精度且能反应机械全局振动强度的基于视觉测振的全场景振动强度图谱分析方法。按以下步骤进行分析：A、利用局部相位变化提取运动信息；B、可视化全场景振动强度图谱的图像处理。本发明选取合适的滤波器通带频率可以得到不同频率下的振动幅值强度图，该可视化方法不仅解决了传统接触式检测中传感器安装位置不确定的问题，而且不同频率下的振动幅值强度图也为基于视觉的结构件的模态振型分析以及健康监测提供了条件。
Description
技术领域
本发明属于机械振动的结构健康检测的领域，这种方法涉及一种基于相位变化提取物体运动信息的视觉测振方法。
背景技术
旋转机械在所有机械设备中占比最大，约占80％，是国民生产生活中的基础性设施。机械的异常振动对机械的正常运转状态和正常工作寿命都会造成极大的负面影响。对于旋转机械而言，异常振动始终是影响其可靠运行和工作寿命的一道难关。有文献表明，约60％的旋转机械故障是由异常振动造成，此类故障不仅会带来经济损失还有可能酿成重大事故。由此可见，振动监测和相应的故障检测对保证各类大型机械尤其是旋转动力机械的安全稳定运行至关重要。针对旋转机械的振动监测和故障检测，不管在大型动力机械，如航空发动机，汽轮机，压缩机，还是在精密机械，如卫星飞轮，动量轮，高速分子泵，磁悬浮高速电机，振动信号始终是反映其运行状态的重要参数，是其故障检测的重要指征，所以可以基于振动信号对指定机械进行故障检测。
基于振动信号的故障检测方法主要有三种，分别是接触式传感器测量、非接触式传感器测量和视觉测量。接触式加速度传感器测量方法是一种传统的基于振动信号的故障检测方法。此类方法虽有较高响应速度。这种方法易于理解和实现，但只能提供单一方向的有限输入，并且无法量化传感器布置位置和数量，在特殊情况下，如存在高温、高压、油气、腐蚀条件下，接触式传感器可能无法进行有效安装，或是与待测机械距离过远导致信号衰减大，造成大的测量误差，此外对于轻型薄壁结构，接触式振动测量会产生质量负载效应。考虑到接触式传感器测量的局限性，非接触式传感器测量以及多种非接触式传感器应运而生，这些传感器包括激光多普勒测振仪，声学测振仪等，虽然实现了非接触式传感器测量，但还是面临全局多点检测和环境噪声的挑战。与此同时，计算机视觉领域和工业视觉设备得到长足的发展，搭载CCD和CMOS成像元件的工业相机能够获得待测设备图像，视觉测量方法的雏形由此显现。视觉检测方法在诸多如缺陷检测，运动检测，目标检测的场合中有广泛应用，国外视觉方案提供商如美国的康耐视，日本的基恩士都有较为成熟的解决方案，但都仅限于像素级的大运动，鲜有在振动检测方面应用，因此针对振动检测的视觉测量方法是有研究的必要和价值的。针对亚像素级的振动，传统的视觉运动算法如光流法，数字图像相关法等并无法有效检测。针对故障检测，传统视觉方法多数仅适用于静态对象。近年，随着测量需求的提高，对微小运动的非接触视觉测量发展迅速，国内外高校也进行了大量研究工作，如国内的华中科技大学、香港理工大学，美国的麻省理工学院都有课题组在进行相关课题研究。
针对旋转机械故障检测，基于视觉的非接触测量技术具有能够多点、多方向测量等先天优势，图像能够全面反映设备运行状态。因此，研究基于机器视觉的旋转机械故障检测对实现旋转机械安全稳定运行的目标具有重要意义。
发明内容
本发明针对以上问题，提出了一种具备亚像素级振动测量精度且能反应机械全局振动强度的基于视觉测振的全场景振动强度图谱分析方法。本发明中采用局部相位运动提取技术提取亚像素级的机械振动的运动信息，并采用可视化技术将全视场的机械振动以色谱图的形式进行呈现出来。
本发明的技术方案为：按以下步骤进行分析：
A、利用局部相位变化提取运动信息；
A1、建立图像强度的数学模型，得到图像序列在时域和频域内的数学模型，结合傅里叶变换的相关性质，从而寻求图像强度随时间变化的影响因素；
A2、分析傅里叶变换处理图像的局限性，转而用复值可控金字塔进行图像分割；
A3、根据复值可控金字塔的性质，通过构造Gabor图像滤波器，将图像转换到时域；
A4、确定Gabor滤波器的表达式，并建立Gabor滤波器处理图像序列的数学模型，从而最终确定A2中图像强度与相位之间存在的关系；
A5、根据步骤A4中确定的图像强度与相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化；
B、可视化全场景振动强度图谱的图像处理；
B1、运算图像中每一个像素点在一定时间段内的相位变化，将相位变化构成全场景的局部相位变化矩阵；
B2、利用归一化和可视化方法，将步骤B1中得到的全场景局部相位变化矩阵转化成利用颜色表征振动强度的色谱图；
B3、根据步骤B2得到全场景的振动强度色谱图后，可选择进行后续的二值化处理和时域滤波处理。
步骤A1中，时刻t的图像强度可假设为I(x,t)，该式为f(x)其随时间t变换的函数，即I(x,t)＝f(x-δ(t))；f(x)与I(x,t)经傅里叶变换后得到以下表达式：
对比(1)(2)两个表达式可以发现，φω→φω+ωδ(t)，因此能够得到变换前后的相位差ωδ(t)，结合傅里叶变换的相依理论，表明空间中物体的任何运动都会导致频域中相位的变化，得到随时间变化的图像强度本质上与相位变化有关。
步骤A4中，Gabor滤波器表达式的一般形式如下：
(3)表达式可分解为实数部分和虚数部分并简化为g(x,y,λ,θ,ψ,σ,γ)＝Gθ+iHθ，其中
其中x′＝xcosθ+ysinθ，y′＝-xsinθ+ycosθ，(3)式中参数λ是以像素为单位的波长，决定了滤波器的空间分辨率；θ是滤波器方向因子，范围为[0,2π]，决定了滤波器空间分解得到相位的方向；γ为滤波器的轮廓因子，决定了滤波器的形状，当γ＝1时，滤波器为圆形，γ<1时，为椭圆；σ为滤波器的高斯因子标准差，该值随滤波器带宽B而变化，一般不直接设置；滤波器的带宽B的表达式如下：
不同参数对滤波器存在不同的影响，通过调整不同的参数，可以得到不同尺寸和方向下的Gabor滤波器；
确定Gabor滤波器表达式后，建立Gabor滤波器处理图像序列的数学模型；假设在t时刻，输入二维图像为I(x,y,t)，利用(4)式中的二维Gabor滤波器对其进行处理，将图像与复值可控金字塔的每一层ψr,θ进行卷积操作，即能得到该层的复值可控金字塔的系数如下：
其中r,θ分别为当前层的尺度和方向；由于复值可控金字塔的基函数为Gabor变换，因此(6)式可转化为：
将图像I(x,y,t)与滤波器Gθ+iHθ进行卷积操作后，图像I在复域内的幅值为A(γ,θ,x,y,t)，相位为φ(γ,θ,x,y,t)；由于Gabor变换为加窗后的傅里叶变换，故此处的相位φ(γ,θ,x,y,t)为局部相位；随时间变化，在时刻t+△t，产生运动(δx,δy)，此时图像I(x,y,t)随时间变化为I(x+δx,y+δy,t+△t)，假设△t极小，即可得到假设特征提取方向为水平方向，也即θ＝0，(7)式中x′,y′即可简化为x′＝x,y′＝y，由卷积的基本定义可将(7)式写成如下积分的连续形式：
由卷积的基本定义将连续形式的表达式转化为离散形式，离散表达式如下：
因为x,y在(8)式和(9)式中是积分变量，可以将积分后的(8)式和(9)式简化为如下表达式：
C(u,v,t)＝T{I(x,y,t),g(x,y,λ,θ,ψ,σ,γ)} (10)
将(3)式代入(8)式，运算(10)式的卷积积分，得到表达式如下：
(11)式和(12)式中的C(u,v,t)，C(u,v,t+△t)分别为图像I经g(x,y,λ,θ,ψ,σ,γ)变换后在复数域内的系数；利用换元积分简化(12)式，即令α,β分别为(12)式中x+δx,y+δy，可得积分表达式如下：
化简(11)式提出与积分无关的复指数项可得到积分表达式如下：
同理化简在时刻t+△t的积分表达式(13)，得到结果如下：
分别计算(14)式和(15)式的相位，分别得到结果如下：
对(16)式和(17)式所得的相位角做差，得到频域内在时间间隔△t下的相位变化的表达式如下：
根据傅里叶变换的相移理论，时域内的运动变化导致频域内的相位变化；因此，图像经过复值可控金字塔(Gabor滤波器)分解后可以得到局部运动信息；由(18)式可知，频域内的相位变化与时域内的位移变化构成正比例关系；这表明，频域内的相位信息可以保留时域内的运动信息；因此，利用复值可控金字塔空间分解后得到不同尺度和方向的局部相位信息，可以由此测量出时域内的运动信息。
步骤A5中，根据步骤A4确定的时域内运动与频域内相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化等其他运动信息；因此，利用复值可控金字塔对图像序列I(x,y,t)进行不同尺度和方向的分解，得到包含局部运动的幅值和相位信息的表达式如下：
其中θ表示金字塔当前层检测方向，γ表示尺度；可以假设时间变化△t→0时局部相位近似为常数，局部相位的表达式如下：
φθ(x,y,t)＝c (13)
对于指定常数c，对(13)式关于时间求偏导，得到偏导表达式如下：
在(14)式中，令u＝dx/dt，v＝dy/dt，u,v分别为x和y方向的运动速度；以水平和竖直方向为例，当检测水平方向运动时，即θ＝0时，当检测竖直方向运动时，即θ＝π/2，/>此时，图像中的像素运动速度表达式如下：
在(15)式的基础上，在水平和竖直方向上分别对时间进行积分，从而得到图像中的像素运动位移表达式如下：
在(16)式中，dx(t0)，dy(t0)分别为图像中的像素在x,y方向上的位移；对于(15)式和(16)式得到的图像中像素的速度和位移，可以通过相机标定的方法将像素的速度和位移转换为实际的速度和位移。
步骤B1中，对二维图像上的任意一个像素点计算当前帧I(x,y,t)和图像序列的第一帧I(x,y,0)之间的局部相位差，就可以得到在像素点Dt(x,y)处的相位变化△φr,θ(x,y,t)，表示如下：
△φr,θ(x,y,t)＝φr,θ(x,y,t)-φr,θ(x,y,0) (17)
在(17)式中，可求解出像素点Dt(x,y)频域内的相位变化，从而得到该点在时域内的运动信息；因为Dt(x,y)可以是图像内的任意一个像素点，所以可以通过运算图像中每个像素点的相位变化，从而获取全场景的振动信息；假设图像分辨率为H×V，通过计算，可以由(17)式得到在0<x<H,0<y<V范围内所有像素点相对于0时刻的像素运动信息，将相位差填入像素坐标Dt(x,y)处即可得到全场景局部相位变化矩阵表达式如下：
在(18)式中，θ＝θ0表示复值可控金字塔的运动提取方向为θ0。
步骤B2中，对于步骤B1得到的全场景局部相位变化矩阵对该矩阵进行归一化和可视化处理，从而将全场景局部相位变化矩阵/>转化为全场景振动强度色谱图，从而通过颜色深浅直观地得到全场景的振动强度。
步骤B3中，对于步骤B2得到的全场景振动强度色谱图，利用二值化和时域滤波进行进一步处理；
首先，对于需要不同振动强度分布的情况，通过改变二值化的门限值，从而得到反映不同振动强度分布的全场景振动强度分布图；
其次，时域滤波可以提取感兴趣的运动频率，因此将经过时域带通滤波器的全场景局部相位变化信息构成局部运动相位差矩阵其中fpa,fpb分别为带通过滤波器的通带上限和通带下限频率。
本发明选取合适的滤波器通带频率可以得到不同频率下的振动幅值强度图，该可视化方法不仅解决了传统接触式检测中传感器安装位置不确定的问题，而且不同频率下的振动幅值强度图也为基于视觉的结构件的模态振型分析以及健康监测提供了条件。
本发明采用局部相位运动提取技术提取亚像素级的机械振动的运动信息，并采用可视化技术将全视场的机械振动以色谱图的形式进行呈现出来；有效克服了现有视觉检测技术存在的种种缺陷，相较于传统视觉检测技术，首先，本案无需在待测物体表面预先布置散斑图或者高对比标记点，节省时间和人力成本，且适用于更多测量环境；其次，本案解决传统视觉方法易受噪声和干扰的问题；最后，本案可通过高效的运算得到全场景的振动强度并将其可视化为振动强度图谱。
附图说明
图1为本发明中复值可控金字塔的Gabor滤波器受不同参数影响的示意图；
图2为本发明中复制可控金字塔的Gabor变换的示意图；
图3为本发明中复值可控金字塔的多尺度、多方向分解示意图。
具体实施方式
为能清楚说明本专利的技术特点，下面通过具体实施方式，并结合其附图，对本专利进行详细阐述。
本发明提出一种基于视觉测振的全场景振动强度图谱分析方法，该方法基于利用局部相位变化提取运动信息的视觉测量技术得出时域内运动信息与频域内相位变化之间的线性关系，这表明频域内的相位变化保留了时域内的运动信息，该方法又基于这一点，运算图像中每个像素的相位变化，从而构成全场景局部相位变化矩阵，再通过可视化技术将全场景局部相位变化矩阵转化为全场景振动强度色谱图，从而实现全场景振动强度的直观呈现。
A、利用局部相位变化提取运动信息；
A1、建立图像强度的数学模型，得到图像序列在时域和频域内的数学模型，结合傅里叶变换的相关性质，从而寻求图像强度随时间变化的影响因素；得到随时间变化的图像强度本质上与相位变化有关。
步骤A1中，时刻t的图像强度可假设为I(x,t),I(x,t)是一维图像的图像强度数学模型，I是强度intensity的缩写，x是像素点位置，t是时间，该表达式表示一维图像中x位置的像素点在t时刻的图像强度；该式为f(x)其随时间t变换的函数，即I(x,t)＝f(x-δ(t)),f(x)表示整个图像序列中第一帧(初始帧)图像的图像强度，等同于I(x,0)；f(x)与I(x,t)经傅里叶变换后得到以下表达式：
对比(1)(2)两个表达式可以发现，φω→φω+ωδ(t)，因此能够得到变换前后的相位差ωδ(t)，结合傅里叶变换的相依理论，表明空间中物体的任何运动都会导致频域中相位的变化，得到随时间变化的图像强度本质上与相位变化有关。
A2、分析傅里叶变换处理图像的局限性，转而考虑用复值可控金字塔进行图像分割；复值可控金字塔可将图像在不同的方向、尺度和位置上进行分解，从而得到空间内局部运动的幅值相位变化。
结合步骤A1经过傅里叶变换得到的两个表达式，注意到用傅里叶变换处理图像存在局限性；由于傅里叶变换的核函数是正弦函数ejωx，因此只能得到全局相位变化；此外，傅里叶变换后的图像在空间不同位置的频率特性相互混叠，无法得到特定方向的相位信息。因此，考虑到傅里叶变换的局限性，采用复系数可控金字塔从而获取局部的相位变化信息。
A3、根据复值可控金字塔的性质，利用一组具有不同方向和尺度参数的Gabor滤波器对图像进行处理，对图像在各个方向、空间尺度和位置上进行分解，从而得到空间内局部运动的幅值相位变化。
由于复值可控金字塔的基函数是一系列经过高斯函数调制后的正弦信号的和函数，即为Gabor变换，所以可通过构造Gabor图像滤波器，将图像转换到时域，并且得到不同方向、尺度和位置的局部相位信息。
A4、确定Gabor滤波器的表达式，并建立Gabor滤波器处理图像序列的数学模型，从而最终确定A2中图像强度与相位之间存在的关系；
步骤A4中，Gabor滤波器表达式的一般形式如下：
(3)表达式可分解为实数部分和虚数部分并简化为g(x,y,λ,θ,ψ,σ,γ)＝Gθ+iHθ，其中
其中x′＝xcosθ+ysinθ，y′＝-xsinθ+ycosθ，(3)式中参数λ是以像素为单位的波长，决定了滤波器的空间分辨率。θ是滤波器方向因子，范围为[0,2π]，决定了滤波器空间分解得到相位的方向。γ为滤波器的轮廓因子，决定了滤波器的形状，当γ＝1时，滤波器为圆形，γ<1时，为椭圆。σ为滤波器的高斯因子标准差，该值随滤波器带宽B而变化，一般不直接设置。滤波器的带宽B的表达式如下：
不同参数(即Gabor滤波器参数，指上文中的θ、γ、σ、B参数，分别表示滤波器方向、轮廓、标准差、带宽)对滤波器存在不同的影响，通过调整不同的参数，可以得到不同尺寸和方向下的Gabor滤波器。图1为调整不同参数得到的Gabor滤波器示意图。
图1中第一行为方向因子θ对滤波器的影响，可见不同θ使得滤波器在空间内能够提取不同方向的运动。第二行为波长λ对滤波器的影响，波长决定了滤波器能够分辨最小的像素区域。可见波长λ越长，能够分辨出的最小像素区域越小。第三行为带宽对滤波器的影响，可见，随着带宽的减小，滤波器的范围变大，相互平行的增强和衰减条纹也越多。由(5)式可知，带宽大小收波长和高斯标准差影响，因此Gabor滤波器的性能最终由多参数共同决定。图2为Gabor变换对于简单一维运动的作用示意图。
确定Gabor滤波器表达式后，建立Gabor滤波器处理图像序列的数学模型。假设在t时刻，输入二维图像为I(x,y,t)，利用(4)式中的二维Gabor滤波器对其进行处理，将图像与复值可控金字塔的每一层ψr,θ进行卷积操作，即能得到该层的复值可控金字塔的系数如下：
其中r,θ分别为当前层的尺度和方向。由于复值可控金字塔的基函数为Gabor变换，因此(6)式可转化为：
将图像I(x,y,t)与滤波器Gθ+iHθ进行卷积操作后，图像I在复域内的幅值为A(γ,θ,x,y,t)，相位为φ(γ,θ,x,y,t)。由于Gabor变换为加窗后的傅里叶变换，故此处的相位φ(γ,θ,x,y,t)为局部相位。随时间变化，在时刻t+△t，产生运动(δx,δy)，此时图像I(x,y,t)随时间变化为I(x+δx,y+δy,t+△t)，假设△t极小，即可得到假设特征提取方向为水平方向，也即θ＝0，(7)式中x′,y′即可简化为x′＝x,y′＝y，由卷积的基本定义可将(7)式写成如下积分的连续形式：
由卷积的基本定义将连续形式的表达式转化为离散形式，离散表达式如下：
因为x,y在(8)式和(9)式中是积分变量，可以将积分后的(8)式和(9)式简化为如下表达式：
C(u,v,t)＝T{I(x,y,t),g(x,y,λ,θ,ψ,σ,γ)} (10)
将(3)式代入(8)式，运算(10)式的卷积积分，得到表达式如下：
(11)式和(12)式中的C(u,v,t)，C(u,v,t+△t)分别为图像I经g(x,y,λ,θ,ψ,σ,γ)变换后在复数域内的系数。利用换元积分简化(12)式，即令α,β分别为(12)式中x+δx,y+δy，可得积分表达式如下：
化简(11)式提出与积分无关的复指数项可得到积分表达式如下：
同理化简在时刻t+△t的积分表达式(13)，得到结果如下：
分别计算(14)式和(15)式的相位，分别得到结果如下：
对(16)式和(17)式所得的相位角做差，得到频域内在时间间隔△t下的相位变化的表达式如下：
根据傅里叶变换的相移理论，时域内的运动变化导致频域内的相位变化。因此，图像经过复值可控金字塔(Gabor滤波器)分解后可以得到局部运动信息。由(18)式可知，频域内的相位变化与时域内的位移变化构成正比例关系。这表明，频域内的相位信息可以保留时域内的运动信息。因此，利用复值可控金字塔空间分解后得到不同尺度和方向的局部相位信息，可以由此测量出时域内的运动信息。
图3为图像经过复值可控金字塔在不同尺度和方向分解后的结果示意图，其中左下图为输入图像，经过8个方向空间分解后得到前两行结果，右下图为经过2倍降采样后分别在八个方向的空间分解结果。可见，经过复系数可控金字塔空间分解后可以得到不同尺度下各个空间方向的运动信息。
A5、根据步骤A4中确定的图像强度与相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化。
步骤A5中，根据步骤A4确定的时域内运动与频域内相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化等其他运动信息。因此，利用复值可控金字塔(Gabor滤波器)对图像序列I(x,y,t)进行不同尺度和方向的分解，得到包含局部运动的幅值和相位信息的表达式如下：
其中θ表示金字塔当前层检测方向，γ表示尺度。可以假设时间变化△t→0时局部相位近似为常数，局部相位的表达式如下：
φθ(x,y,t)＝c (13)
对于指定常数c，对(13)式关于时间求偏导，得到偏导表达式如下：
在(14)式中，令u＝dx/dt，v＝dy/dt，u,v分别为x和y方向的运动速度。以水平和竖直方向为例，当检测水平方向运动时，即θ＝0时，当检测竖直方向运动时，即θ＝π/2，/>此时，图像中的像素运动速度表达式如下：
在(15)式的基础上，在水平和竖直方向上分别对时间进行积分，从而得到图像中的像素运动位移表达式如下：
在(16)式中，dx(t0)，dy(t0)分别为图像中的像素在x,y方向上的位移。对于(15)式和(16)式得到的图像中像素的速度和位移，可以通过相机标定的方法将像素的速度和位移转换为实际的速度和位移。
B、可视化全场景振动强度图谱的图像处理；
B1、基于上述利用局部相位变化提取运动信息的视觉测量技术，运算图像中每一个像素点在一定时间段内的相位变化，将相位变化构成全场景的局部相位变化矩阵。
步骤B1中，由于上述的利用局部相位变化提取运动信息的视觉测量技术中已经证明时域内的运动信息与频域内的相位变化存在线性关系，因此，可以通过像素点在频域内的相位变化表征该点在时域内的运动信息。基于此观点，对二维图像上的任意一个像素点计算当前帧I(x,y,t)和图像序列的第一帧I(x,y,0)之间的局部相位差，就可以得到在像素点Dt(x,y)处的相位变化△φr,θ(x,y,t)，表示如下
△φr,θ(x,y,t)＝φr,θ(x,y,t)-φr,θ(x,y,0) (17)
在(17)式中，可求解出像素点Dt(x,y)频域内的相位变化，从而得到该点在时域内的运动信息。因为Dt(x,y)可以是图像内的任意一个像素点，所以可以通过运算图像中每个像素点的相位变化，从而获取全场景的振动信息。假设图像分辨率为H×V，通过计算，可以由(17)式得到在0<x<H,0<y<V范围内所有像素点相对于0时刻的像素运动信息，将相位差填入像素坐标Dt(x,y)处即可得到全场景局部相位变化矩阵表达式如下：
在(18)式中，θ＝θ0表示复值可控金字塔的运动提取方向为θ0。
B2、利用归一化和可视化方法，将步骤B1中得到的全场景局部相位变化矩阵转化成利用颜色表征振动强度的色谱图。
步骤B2中，对于步骤B1得到的全场景局部相位变化矩阵可以对该矩阵进行归一化和可视化处理，因为图像中各像素点的相位变化可以反映该点的运动信息，所以可以将全场景局部相位变化矩阵/>转化为全场景振动强度色谱图，该色谱图的颜色可以表征相位点的振动强度，颜色越鲜明表示该区域的振动强度越大，颜色越暗淡表示该区域的振动强度越小，从而可以通过颜色深浅直观地得到全场景的振动强度。
B3、根据步骤B2得到全场景的振动强度色谱图后，可选择进行后续的二值化处理和时域滤波处理。
步骤B3中，对于步骤B2得到的全场景振动强度色谱图，可以利用二值化和时域滤波进行进一步处理。首先，对于需要不同振动强度分布的情况，可以通过改变二值化的门限值，从而得到反映不同振动强度分布的全场景振动强度分布图。其次，时域滤波可以提取感兴趣的运动频率，因此将经过时域带通滤波器的全场景局部相位变化信息构成局部运动相位差矩阵其中fpa,fpb分别为带通过滤波器的通带上限和通带下限频率。选取合适的滤波器通带频率可以得到不同频率下的振动幅值强度图，该可视化方法不仅解决了传统接触式检测中传感器安装位置不确定的问题，而且不同频率下的振动幅值强度图也为基于视觉的结构件的模态振型分析以及健康监测提供了条件。
本发明具体实施途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。
Claims (1)
1.一种基于视觉测振的全场景振动强度图谱分析方法，其特征在于，按以下步骤进行分析：
A、利用局部相位变化提取运动信息；
A1、建立图像强度的数学模型，得到图像序列在时域和频域内的数学模型，结合傅里叶变换的相关性质，从而寻求图像强度随时间变化的影响因素；
A2、分析傅里叶变换处理图像的局限性，转而用复值可控金字塔进行图像分割；
A3、根据复值可控金字塔的性质，通过构造Gabor图像滤波器，将图像转换到时域；
A4、确定Gabor滤波器的表达式，并建立Gabor滤波器处理图像序列的数学模型，从而最终确定A2中图像强度与相位之间存在的关系；
A5、根据步骤A4中确定的图像强度与相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化；
B、可视化全场景振动强度图谱的图像处理；
B1、运算图像中每一个像素点在一定时间段内的相位变化，将相位变化构成全场景的局部相位变化矩阵；
B2、利用归一化和可视化方法，将步骤B1中得到的全场景局部相位变化矩阵转化成利用颜色表征振动强度的色谱图；
B3、根据步骤B2得到全场景的振动强度色谱图后，选择进行后续的二值化处理和时域滤波处理；
步骤A1中，时刻t的图像强度假设为I(x,t)，该式为f(x)其随时间t变换的函数，即I(x,t)＝f(x-δ(t))；f(x)与I(x,t)经傅里叶变换后得到以下表达式：
对比(1)(2)两个表达式发现，φω→φω+ωδ(t)，因此能够得到变换前后的相位差ωδ(t)，结合傅里叶变换的相依理论，表明空间中物体的任何运动都会导致频域中相位的变化，得到随时间变化的图像强度本质上与相位变化有关；
步骤A4中，Gabor滤波器表达式的一般形式如下：
(3)表达式分解为实数部分和虚数部分并简化为g(x,y,λ,θ,ψ,σ,γ)＝Gθ+iHθ，其中
其中x′＝xcosθ+ysinθ，y′＝-xsinθ+ycosθ，(3)式中参数λ是以像素为单位的波长，决定了滤波器的空间分辨率；θ是滤波器方向因子，范围为[0,2π]，决定了滤波器空间分解得到相位的方向；γ为滤波器的轮廓因子，决定了滤波器的形状，当γ＝1时，滤波器为圆形，γ<1时，为椭圆；σ为滤波器的高斯因子标准差，该值随滤波器带宽B而变化，一般不直接设置；滤波器的带宽B的表达式如下：
不同参数对滤波器存在不同的影响，通过调整不同的参数，得到不同尺寸和方向下的Gabor滤波器；
确定Gabor滤波器表达式后，建立Gabor滤波器处理图像序列的数学模型；假设在t时刻，输入二维图像为I(x,y,t)，利用(4)式中的二维Gabor滤波器对其进行处理，将图像与复值可控金字塔的每一层ψr,θ进行卷积操作，即能得到该层的复值可控金字塔的系数如下：
其中r,θ分别为当前层的尺度和方向；由于复值可控金字塔的基函数为Gabor变换，因此(6)式转化为：
将图像I(x,y,t)与滤波器Gθ+iHθ进行卷积操作后，图像I在复域内的幅值为A(γ,θ,x,y,t)，相位为φ(γ,θ,x,y,t)；由于Gabor变换为加窗后的傅里叶变换，故此处的相位φ(γ,θ,x,y,t)为局部相位；随时间变化，在时刻t+△t，产生运动(δx,δy)，此时图像I(x,y,t)随时间变化为I(x+δx,y+δy,t+△t)，假设△t极小，得到假设特征提取方向为水平方向，也即θ＝0，(7)式中x′,y′简化为x′＝x,y′＝y，由卷积的基本定义将(7)式写成如下积分的连续形式：
由卷积的基本定义将连续形式的表达式转化为离散形式，离散表达式如下：
因为x,y在(8)式和(9)式中是积分变量，将积分后的(8)式和(9)式简化为如下表达式：
C(u,v,t)＝T{I(x,y,t),g(x,y,λ,θ,ψ,σ,γ)}(10)
将(3)式代入(8)式，运算(10)式的卷积积分，得到表达式如下：
(11)式和(12)式中的C(u,v,t)，C(u,v,t+△t)分别为图像I经g(x,y,λ,θ,ψ,σ,γ)变换后在复数域内的系数；利用换元积分简化(12)式，即令α,β分别为(12)式中x+δx,y+δy，得积分表达式如下：
化简(11)式提出与积分无关的复指数项得到积分表达式如下：
同理化简在时刻t+△t的积分表达式(13)，得到结果如下：
分别计算(14)式和(15)式的相位，分别得到结果如下：
对(16)式和(17)式所得的相位角做差，得到频域内在时间间隔△t下的相位变化的表达式如下：
根据傅里叶变换的相移理论，时域内的运动变化导致频域内的相位变化；因此，图像经过复值可控金字塔(Gabor滤波器)分解后得到局部运动信息；由(18)式知，频域内的相位变化与时域内的位移变化构成正比例关系；这表明，频域内的相位信息保留时域内的运动信息；因此，利用复值可控金字塔空间分解后得到不同尺度和方向的局部相位信息，由此测量出时域内的运动信息；
步骤A5中，根据步骤A4确定的时域内运动与频域内相位之间的关系，由空间内局部运动的相位变化推导出速度变化和位移变化等其他运动信息；因此，利用复值可控金字塔对图像序列I(x,y,t)进行不同尺度和方向的分解，得到包含局部运动的幅值和相位信息的表达式如下：
其中θ表示金字塔当前层检测方向，γ表示尺度；假设时间变化△t→0时局部相位近似为常数，局部相位的表达式如下：
φθ(x,y,t)＝c(13)
对于指定常数c，对(13)式关于时间求偏导，得到偏导表达式如下：
在(14)式中，令u＝dx/dt，v＝dy/dt，u,v分别为x和y方向的运动速度；以水平和竖直方向为例，当检测水平方向运动时，即θ＝0时，当检测竖直方向运动时，即θ＝π/2，/>此时，图像中的像素运动速度表达式如下：
在(15)式的基础上，在水平和竖直方向上分别对时间进行积分，从而得到图像中的像素运动位移表达式如下：
在(16)式中，dx(t0)，dy(t0)分别为图像中的像素在x,y方向上的位移；对于(15)式和(16)式得到的图像中像素的速度和位移，通过相机标定的方法将像素的速度和位移转换为实际的速度和位移；
步骤B1中，对二维图像上的任意一个像素点计算当前帧I(x,y,t)和图像序列的第一帧I(x,y,0)之间的局部相位差，就得到在像素点Dt(x,y)处的相位变化△φr,θ(x,y,t)，表示如下：
△φr,θ(x,y,t)＝φr,θ(x,y,t)-φr,θ(x,y,0)(17)
在(17)式中，求解出像素点Dt(x,y)频域内的相位变化，从而得到该点在时域内的运动信息；因为Dt(x,y)是图像内的任意一个像素点，所以通过运算图像中每个像素点的相位变化，从而获取全场景的振动信息；假设图像分辨率为H×V，通过计算，由(17)式得到在0<x<H,0<y<V范围内所有像素点相对于0时刻的像素运动信息，将相位差填入像素坐标Dt(x,y)处得到全场景局部相位变化矩阵表达式如下：
在(18)式中，θ＝θ0表示复值可控金字塔的运动提取方向为θ0；
步骤B2中，对于步骤B1得到的全场景局部相位变化矩阵对该矩阵进行归一化和可视化处理，从而将全场景局部相位变化矩阵/>转化为全场景振动强度色谱图，从而通过颜色深浅直观地得到全场景的振动强度；
步骤B3中，对于步骤B2得到的全场景振动强度色谱图，利用二值化和时域滤波进行进一步处理；
首先，对于需要不同振动强度分布的情况，通过改变二值化的门限值，从而得到反映不同振动强度分布的全场景振动强度分布图；
其次，时域滤波提取感兴趣的运动频率，因此将经过时域带通滤波器的全场景局部相位变化信息构成局部运动相位差矩阵其中fpa,fpb分别为带通过滤波器的通带上限和通带下限频率。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010111389.2A CN111353400B (zh) | 2020-02-24 | 2020-02-24 | 一种基于视觉测振的全场景振动强度图谱分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010111389.2A CN111353400B (zh) | 2020-02-24 | 2020-02-24 | 一种基于视觉测振的全场景振动强度图谱分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111353400A CN111353400A (zh) | 2020-06-30 |
| CN111353400B true CN111353400B (zh) | 2024-04-02 |
Family
ID=71195706
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010111389.2A Active CN111353400B (zh) | 2020-02-24 | 2020-02-24 | 一种基于视觉测振的全场景振动强度图谱分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111353400B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113076517B (zh) * | 2021-04-01 | 2022-09-30 | 重庆大学 | 基于希尔伯特变换的土木工程结构动态监测相位评估方法 |
| CN113724334B (zh) * | 2021-07-19 | 2023-09-15 | 广东工业大学 | 一种基于机器视觉的电梯振动检测方法和系统 |
| CN113781522B (zh) * | 2021-08-25 | 2023-10-24 | 西安交通大学 | 一种基于计算机视觉的射击工况下火炮身管振动测量方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106251344A (zh) * | 2016-07-26 | 2016-12-21 | 北京理工大学 | 一种基于视觉感受野的多尺度红外目标自适应检测方法 |
| CN110084127A (zh) * | 2019-03-29 | 2019-08-02 | 南京航空航天大学 | 一种基于视觉的磁悬浮转子振动测量方法 |
| CN110261083A (zh) * | 2019-06-13 | 2019-09-20 | 南京航空航天大学 | 一种基于视觉的磁悬浮转子振动力抑制效果测量方法 |
-
2020
- 2020-02-24 CN CN202010111389.2A patent/CN111353400B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106251344A (zh) * | 2016-07-26 | 2016-12-21 | 北京理工大学 | 一种基于视觉感受野的多尺度红外目标自适应检测方法 |
| CN110084127A (zh) * | 2019-03-29 | 2019-08-02 | 南京航空航天大学 | 一种基于视觉的磁悬浮转子振动测量方法 |
| CN110261083A (zh) * | 2019-06-13 | 2019-09-20 | 南京航空航天大学 | 一种基于视觉的磁悬浮转子振动力抑制效果测量方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111353400A (zh) | 2020-06-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111353400B (zh) | 一种基于视觉测振的全场景振动强度图谱分析方法 | |
| Chen et al. | Modal identification of simple structures with high-speed video using motion magnification | |
| Diamond et al. | Accuracy evaluation of sub-pixel structural vibration measurements through optical flow analysis of a video sequence | |
| CN112254801B (zh) | 一种微小振动视觉测量方法及系统 | |
| CN110084127B (zh) | 一种基于视觉的磁悬浮转子振动测量方法 | |
| Liu et al. | Structural motion estimation via Hilbert transform enhanced phase-based video processing | |
| Sarrafi et al. | Mode extraction on wind turbine blades via phase-based video motion estimation | |
| Peng et al. | Camera-based micro-vibration measurement for lightweight structure using an improved phase-based motion extraction | |
| CN106441896A (zh) | 滚动轴承故障模式识别及状态监测的特征向量提取方法 | |
| Javed et al. | Vibration measurement of a rotating cylindrical structure using subpixel-based edge detection and edge tracking | |
| CN112381860B (zh) | 一种旋转叶片动频测量的无标记计算机视觉方法 | |
| CN105388012A (zh) | 基于非线性调频小波变换的阶次跟踪方法 | |
| Zhou et al. | Vibration measurement with video processing based on alternating optimization of frequency and phase shifts | |
| Havaran et al. | Markers tracking and extracting structural vibration utilizing Randomized Hough transform | |
| US20220026311A1 (en) | Tire testing apparatus | |
| CN111060315A (zh) | 一种基于视觉的机械故障诊断方法 | |
| Zhu et al. | A visual measurement method of structural body vibration displacement combined with image deblurring | |
| Yang et al. | Sparse representation of complex steerable pyramid for machine fault diagnosis by using non-contact video motion to replace conventional accelerometers | |
| Deng et al. | Synchronous monitoring of axial vibration and rotation speed of rotating cylinder by linear array scanning | |
| Lu et al. | Observation of tower vibration based on subtle motion magnification | |
| Sarrafi et al. | Using 2D phase-based motion estimation and video magnification for binary damage identification on a wind turbine blade | |
| Yang et al. | Image analyses for video-based remote structure vibration monitoring system | |
| US11977096B2 (en) | Motion, vibration and aberrant condition detection and analysis | |
| Ragulskis et al. | Contrast enhancement of time-averaged fringes based on moving average mapping functions | |
| Gwashavanhu et al. | Shape principal component analysis as a targetless photogrammetric technique for condition monitoring of rotating machines |
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 |