CN114564685A - 一种基于相位的工程结构双向动态位移测量方法 - Google Patents
一种基于相位的工程结构双向动态位移测量方法 Download PDFInfo
- Publication number
- CN114564685A CN114564685A CN202210214636.0A CN202210214636A CN114564685A CN 114564685 A CN114564685 A CN 114564685A CN 202210214636 A CN202210214636 A CN 202210214636A CN 114564685 A CN114564685 A CN 114564685A
- Authority
- CN
- China
- Prior art keywords
- image
- phase
- frequency
- formula
- monocycle
- 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.)
- Pending
Links
- 238000006073 displacement reaction Methods 0.000 title claims abstract description 17
- 230000002457 bidirectional effect Effects 0.000 title claims abstract description 16
- 238000000691 measurement method Methods 0.000 title abstract description 7
- 238000000034 method Methods 0.000 claims abstract description 74
- 230000009466 transformation Effects 0.000 claims abstract description 25
- 239000011159 matrix material Substances 0.000 claims abstract description 23
- 238000012545 processing Methods 0.000 claims abstract description 12
- 238000004364 calculation method Methods 0.000 claims abstract description 9
- 238000011156 evaluation Methods 0.000 claims abstract description 5
- 238000000844 transformation Methods 0.000 claims abstract description 4
- 230000008859 change Effects 0.000 claims description 10
- 238000001914 filtration Methods 0.000 claims description 7
- 238000012546 transfer Methods 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 3
- 230000021615 conjugation Effects 0.000 claims description 2
- 239000000203 mixture Substances 0.000 claims description 2
- 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 1
- 238000012544 monitoring process Methods 0.000 abstract description 9
- 230000000052 comparative effect Effects 0.000 description 9
- 230000000694 effects Effects 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000005286 illumination Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 101100242814 Caenorhabditis elegans parg-1 gene Proteins 0.000 description 1
- 101100407084 Caenorhabditis elegans parp-2 gene Proteins 0.000 description 1
- 102100037834 Protein phosphatase methylesterase 1 Human genes 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000011089 mechanical engineering Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000003534 oscillatory effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 108010086028 protein phosphatase methylesterase-1 Proteins 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10016—Video; Image sequence
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20048—Transform domain processing
- G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Discrete Mathematics (AREA)
- Multimedia (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于相位的工程结构双向动态位移测量方法，按照以下步骤进行：S1、计算得到图像序列In(x,y,t)的谐振频率fs；S2得到频率矩阵F(u,v,t)；S3、得到滤波后的频率矩阵FB(u,v,t)；S4、得到t时刻的单周期图像序列Is(x,y,t)；S5、基于Riesz变换处理单周期图像Is(x,y,t)，同时计算评估图像中结构的双向运动。采用以上技术方案的一种基于相位的工程结构双向动态位移测量方法，将Riesz变换引入到相位运动估计中，提出了一种可以同时获取两个方向运动的视频运动估计方法，解决了现有PME方法中的局限，使计算效率得到大幅提升，更适合于动态结构监测，为民用基础设施的动态监测提供了一种简单实用的新方法。
Description
技术领域
本发明涉及土木工程结构动态监测技术领域，具体涉及一种基于相位的工程结构双向动态位移测量方法。
背景技术
结构健康监测(简称SHM)能够实现低成本条件的维护，执行结构状态预测，并消除意外灾难性故障。用于SHM的标准有线加速度计是劳动密集型的仪器结构，非接触式测量方法避免了这些缺点。数码摄像机作为一种非接触式监测方法，因其成本较低，操作灵活，并能够提供非常高的空间分辨率的优势成为了近期研究的热点。结合图像处理方法、基于视频的测量已经成功地用于各种类型结构的振动测量。最近，SHM中引入了数字图像相关(简称DIC)和点跟踪(简称DPT)技术。
然而，这些方法通常需要在结构表面预制散斑模式图或高对比度的特征点标记，并使用基于图像灰度的相关性来估计结构运动场。一方面，这些要求提出了表面预处理的问题，特别是当测量结构面积大或环境复杂时这一要求很难达到。另一方面，基于图像灰度本身的相关性对环境噪声和干扰非常敏感，如场景光照的突然变化会导致结构运动场的估计产生较大的误差。
因此，相位运动估计法(简称PME)作为图像估计的一种，无需任何表面预处理，同时该方法是基于图像的相位来实现结构的运动估计的，而相位行为变化对光照变化相对不敏感(如阴影、高光等)。因此，相位轮廓确实提供了一个很好的近似运动场。
Aral Sarrafi等人在图像强度上应用一维Gabor小波来提取结构运动。因为Gabor小波是由高斯包络的复正弦函数组成，因此使用Gabor小波获取的相位为局部相位，这可能导致结构的运动估计精度不足。为了提高计算精度，本发明人团队曾在中国专利申请CN113076517A中提出了一种基于希尔伯特变换的相位运动估计法，该方法利用Hilbert变换代替Gabor小波获取瞬时结构的相位以提高原始PME的精度。之后，Aral Sarrafi等人对Gabor小波进行扩展，应用2D Gabor波提取图像中结构的局部相位变化，并使用2.3m长风力机叶片的损伤来验证方法的损伤检测能力。虽然采用了二维Gabor小波估计运动，但上述这些PME方法都集中在一维相位估计上。也就是说，一次只能考虑结构的一个方向运动，但图像通常是二维的，所以这些PME方法至少在理论上是不完善的，难以适应于动态结构监测。并且，在土木工程或机械工程中，结构的双向运动是常见的要求，例如高层建筑在地震作用下的强弱轴向位移就是双向运动，导致现有的这些PME方法都不适用。
解决以上问题成为当务之急。
发明内容
为解决以上的技术问题，本发明提供了一种基于相位的工程结构双向动态位移测量方法。
其技术方案如下：
一种基于相位的工程结构双向动态位移测量方法，其要点在于，按照以下步骤进行：
S1、任选图像序列In(x,y,t)中的一张灰度图像Ir(x,y)，基于该灰度图像Ir(x,y)计算得到图像序列In(x,y,t)的谐振频率fs；
S2、利用二维快速傅里叶变换处理图像序列In(x,y,t)，使其由时域转移到频域，得到频率矩阵F(u,v,t)；
S3、使用中心频率为fs、带宽为1Hz的巴特沃思理想带通滤波器对频率矩阵F(u,v,t)在频率区间[fs-0.5,fs+0.5]进行滤波，得到滤波后的频率矩阵FB(u,v,t)；
S4、利用二维傅里叶逆变换处理滤波后的频率矩阵FB(u,v,t)，得到t时刻的单周期图像序列Is(x,y,t)；
S5、基于Riesz变换处理单周期图像Is(x,y,t)，同时计算评估图像中结构的双向运动；
所述步骤S5按照以下步骤进行：
S51、将t时刻的单周期图像Is(x,y,t)视为一个矩阵；
S52、对单周期图像Is(x,y,t)进行Riesz变换；
S53、根据步骤S52的Riesz变换结果，评估图像中结构同时在x方向和y方向运动的相位差；
S54、根据步骤S53得到的结构同时在x方向和y方向运动的相位差，建立结构双向运动与相位差之间的关系，从而获取结构在x方向和y方向的运动。
与现有技术相比，本发明的有益效果：
采用以上技术方案的一种基于相位的工程结构双向动态位移测量方法，将Riesz变换引入到相位运动估计中，提出了一种可以同时获取两个方向运动的视频运动估计方法，解决了现有PME方法中的局限，使计算效率得到大幅提升，更适合于动态结构监测，为民用基础设施的动态监测提供了一种简单实用的新方法。
附图说明
图1为频域中Riesz变换算子的示意图；
图2为输入信号I与其Riesz变换(R1,R2)的组合示意图；
图3为采用巴特沃思理想带通滤波器进行滤波的示意图；
图4在对比例一的其中一个振荡位置处的抛物面示意图；
图5为对比例一的使用不同参数的PME方法和本发明方法(RPME)识别的x方向运动与真实运动对比图；
图6为对比例一的使用不同参数的PME方法和本发明方法(RPME)识别的y方向运动与真实运动对比图；
图7为对比例一的采用PME方法和本发明方法(RPME)在x方向和y方向运动识别效果的平均绝对误差和标准偏差的对比图；
图8为对比例二的高层建筑缩尺模型实验图；
图9为对比例二的使用不同参数的PME方法和本发明方法(RPME)的运动识别曲线对比图((a)为x方向，(b)为y方向)；
图10为对比例二的采用PME方法和本发明方法(RPME)在x方向和y方向运动识别效果的平均绝对误差和标准偏差的对比图。
具体实施方式
以下结合实施例和附图对本发明作进一步说明。
一种基于相位的工程结构双向动态位移测量方法，按照以下步骤进行：
S1、任选图像序列In(x,y,t)中的一张灰度图像Ir(x,y)，基于该灰度图像Ir(x,y)计算得到图像序列In(x,y,t)的谐振频率fs；
具体地说，步骤S1按照以下步骤进行：
S11、任选图像序列In(x,y,t)中的一张灰度图像Ir(x,y)；
S12、利用快速傅里叶变换处理灰度图像Ir(x,y)，使其由时域转移到频域，得到频率域Xr(u,v)，Xr(u,v)表达式为：
式(1)中，x和y分别表示选取位置在空域坐标系上的坐标，u和v分别表示空域中x和y在频域中对应的自变量，j表示虚部，M和N分别表示灰度图像Ir(x,y)的长度和宽度，Xr(u,v)表示灰度图像Ir(x,y)的频率域；
S13、利用功率谱密度法处理频率序列Xr(u,v)，得到以下表达式：
计算PSD(u,v)最大值对应的频率u和v，从而计算得到谐振频率fs＝(u2+v2)1/2，即：选择峰值点为谐振频率fs＝(u2+v2)1/2。
S2、利用二维快速傅里叶变换处理图像序列In(x,y,t)，使其由时域转移到频域，得到频率矩阵F(u,v,t)，F(u,v,t)的表达式为：
式(3)中，x和y表示空域中的自变量，u和v表示空域中x和y在频域中对应的自变量，j表示虚部；
S3、请参见图3，使用中心频率为fs、带宽为1Hz的巴特沃思理想带通滤波器对频率矩阵F(u,v,t)在频率区间[fs-0.5,fs+0.5]进行滤波，得到滤波后的频率矩阵FB(u,v,t)。
具体地说，滤波后的频率矩阵FB(u,v,t)与频率矩阵F(u,v,t)的关系式为：
FB(u,v,t)＝F(u,v,t)*BH(u,v) (4)
式(4)中，BH(u,v)表示频域中的巴斯沃特思理想带通滤波器，其中BH(u,v)的表达式为：
式(5)中，D(u,v)＝[(u-M/2)2+(v-N/2)2]/2，D(u,v)表示在频域中的点(u,v)到圆心的距离，M和N分别表示图像的长度和宽度，D1＝fs+0.5,D2＝fs-0.5。
S4、利用二维傅里叶逆变换处理滤波后的频率矩阵FB(u,v,t)，得到t时刻的单周期图像序列Is(x,y,t)，Is(x,y,t)的表达式为：
S5、基于Riesz变换处理单周期图像Is(x,y,t)，同时计算评估图像中结构的双向运动；
具体地说，步骤S5按照以下步骤进行：
S51、将t时刻的单周期图像Is(x,y,t)视为一个矩阵；
S52、对单周期图像Is(x,y,t)进行Riesz变换；
对单周期图像Is(x,y,t)进行Riesz变换的公式为：
式(7)中，R[Is(x,y,t)]表示对单周期图像Is(x,y,t)进行Riesz变换，其中，Riesz算子R[.]＝(R1,R2)，R1表示x方向的Riesz变换结果，R2表示y方向的Riesz变换结果，i表示虚数，u＝[u1,u2]T表示频域信号维数，u1和u2分别表示单周期图像Is(x,y,t)在频域中与空域中的自变量x和y对应的变量，上标T表示转置运算，Is(u)表示t时刻的单周期图像Is(x,y,t)的傅里叶变换结果。
如图1所示的一个Riesz算子，其中在求关于u1轴的Riesz变换时，u2＝0，通过符号sign(u1)＝u1/|u1 2+0|1/2得到，这可以通过沿图1中的u2＝0等高线来可视化。
S53、根据步骤S52的Riesz变换结果，评估图像中结构同时在x方向和y方向运动的相位差；
步骤S53中，将任意t时刻的单周期图像Is(x,y，t)及其Riesz变换结果(R1,R2)共同构成三元，并构建球坐标系，如图2所示，同时将Is(x,y，t)记作I，可以写成如下方程：
I＝Acos(φ),R1＝Asin(φ)cos(θ),R2＝Asin(φ)sin(θ) (8)
式(8)中，A为振幅，θ为方向，φ为空间相位。
从图2可以看出，基于上述球坐标系，实轴I与空间向量(A,φ,θ)之间的夹角定义为空间相位φ，实轴I在R1轴和R2轴上的投影分别能够作为在x方向和y方向的瞬时相位；因此，为了计算空间相位φ，可以将三元(I,R1,R2)写成一个四元数q，其中q是由一个实部I和两个虚部(R1,R2)组成，公式如下：
q＝q0+iq1+jq2 (9)
式(9)中，其中q0＝I，q1＝R1，q2＝R2；故根据四元数表示法，i和j分别表示两个虚数单位，振幅A＝|q|；并且将相邻两帧图像分别表示为q＝q0+iq1+jq2和r＝r0+ir1+jr2，则相邻两帧图像的相位变化等于r/q，即:
Phase(φ(r/q))＝Phase(φ(r×q*/|q|2))＝Phase(φ(r×q*)) (10)
式(10)中，Phase(φ(r/q))表示以四元数形式表示的相邻两帧图像q和r的相邻图像相位差，q*表示q的共轭，其中对于r×q*表示为：
式(11)中，r0、r1和r2分别表示与相邻帧图像对应的I、R1和R2；
将r×q*写成一个四元数p，即：r×q*＝p＝p0+ip1+jp2，故基于式(11)可知：p0＝r0q0+r1q1+r2q2，p1＝-r0q1+r1q0，p2＝-r0q2+r2q0，四元数p的实部是Re(p)＝p0，虚部是Im(p)＝pv＝ip1+jp2，因此有：
式(12)中，pv表示p的虚部；
结合公式(10)和公式(12)可知，相邻两帧的相位变化Phase(φ(r/q))＝Phase(φ(r×q*))，也就是等于公式(12)的相位，由于只有虚部对计算相位有用，故：
因此，相邻两帧图像的相位差Phase(φ(r/q))为：
方向θ为：
在x方向和y方向的相位差Δφx和Δφy为：
S54、根据步骤S53得到的结构同时在x方向和y方向运动的相位差，建立结构双向运动与相位差之间的关系，从而获取结构在x方向和y方向的运动。
步骤S54中，任意时刻t的单周期图像序列Is(x,y，t)可以表示为：
I(x,t)＝L1cos(ωx) (17)
式(17)中，用I(x,t)表示Is(x,y，t)，其中向量x＝[x,y]T表示图像像素灰度所在位置，L1表示t时刻图像的灰度幅度，ω＝[ωx,ωy]T表示图像的频率，ωx和ωy分别表示x方向和y方向的频率；
由式(7)和(17)可知，正弦函数L1 cos(ωx)的Riesz变换对为：
正弦函数L1 cos(ωx)的正交对Q(x,t)为：
Q(x,t)＝L1sin(ωx) (19)
与一维情况一致，故瞬时相位φ1(x,t)为:
φ1(x,t)＝ωx (20)
相同地，在x方向和y方向运动为δ＝[δx,δy]T的下一帧图像I(x+δ,t+Δt)可以表示为：
I(x+δ,t+Δt)＝L2cos(ωx) (21)
式(21)中，L2表示下一帧图像的灰度幅度，δx和δy分别表示为x方向和y方向的运动；
在x方向和y方向运动为δ＝[δx,δy]T的下一帧图像I(x+δ,t+Δt)的瞬时相位为：
φ2(x,t+Δt)＝ω(x+δ) (22)
故相位差Δφ为：
Δφ＝φ2(x,t+Δt)-φ1(x,t)＝ωδ (23)
式(23)中，Δφ＝[Δφx,Δφy]T,ω＝[ωx,ωy]T，δ＝[δx,δy]T，结合公式(16)和公式(23)可得：
Δφx＝ωxδx,Δφy＝ωyδy (24)
由式(24)可知，能够根据空间相位变化计算得到x方向运动δx和y方向运动δy。
对比例一：有限元实验
采用一个同时沿x方向和y方向移动的二维抛物曲面来模拟一组图像/或视频：
Z(x,y,t)＝-((x-δx(t))2+(y-δy(t))2)+120 (25)
式(25)中，Z(x,y,t)表示二维抛物曲面，δx(t)和δy(t)分别是任意的运动位移,在这里请参见图4,采用运动函数δx(t)＝e-ξωntsin(ωnt)来模拟粘性阻尼振荡的抛物曲面表面沿x方向和y方向的轴运动，其中沿x轴的运动设置为ωn＝15和ξ＝0.09,沿y轴的运动设置为ωn＝4和ξ＝0.3，其中，ωn和ξ表示运动函数中的自变量。
采用本方法(简称RPME)和传统PME方法分别对视频进行处理，提取运动信息。在RPME中，使用上述提到的巴斯沃特理想带通滤波器将原始视频分解成单周期图像，对应的fs＝2。在PME中，当计算x方向运动时，设置Gabor参数中的θ＝0，波长λ随机设置为80、50和30，其中，θ表示Gabor波中的高斯方向；当计算y方向运动时，设置Gabor参数中的θ＝π/2，波长λ设置与x方向一样。识别结果对比图如图5和6所示。这里的PME-1,PME-2和PME-3分别指波长λ设置为80、50和30的识别结果。为了定量地比较这两种方法，使用平均绝对误差(Meanabsolute error,简称MAE)和标准偏差(Standard deviation,简称STD)来评价两种识别方法的效果，其定义如下:
式(26)中，obvervedt表示第t个计算点的预估运动值，truet表示第t个计算点的真实运动值，n表示计算点的总数。
使用MAE和STD来评价RPME和PME方法的识别效果如附图7所示。
从图5和图6可以看出，本发明方法提出的RPME方法在提取x方向和y方向的运动信息时，识别效果最好，最大误差在1pixel以内。这是因为在Riesz变换中采用了瞬时相位，Riesz变换是希尔伯特的二维扩展。因此，该Riesz变换提高了所提方法的识别精度。而PME方法则是基于局部高斯包络的正弦信号来提取图像的局部相位，因此PME方法识别的运动信息的精度不如RPME方法高。此外，波长的选择也降低了PME方法的识别精度。当选择合适的波长时，识别效果可以达到合理的程度，如图5(a)和6(a)所示，否则可能会产生较大的误差，如图5(b)和6(b)所示。换句话说，使用PME的一个限制是它需要预先知道要测量的对象。另一方面，与每次只考虑图像单一方向相位的PME相比，本发明方法提出的RPME方法提出的空间相位同时识别两个方向的相位更加有效。
由图7可以看出，RPME方法在y方向上的MAE和STD均小于x方向上的MAE和STD。这是因为x方向的移动幅度(18.5像素)大于y方向的移动幅度(13.5像素)。实际上，x和y方向的变异系数分别为0.194和0.198。因此，该方法在两个方向上的识别能力是相同的。但可以看出，RPME方法在x方向和y方向上的MAE至少比PME方法低50.83％和32.17％。RPME方法得到的x方向和y方向的STD比PME方法至少低28.57％和11.81％。因此，该方法的识别精度优于传统的PME方法。
对比例二：高层建筑缩尺模型实验
如图8所示，实验室验证采用一个长宽高分别为1000mm×200mm×55mm的高层建筑缩尺模型来实现。将智能手机安装在模型的正上方，平行于xy平面，完成对模型顶部双向运动信息的录像。手机距离模型顶部平面67mm。手机的分辨率为3840×2160，帧率为60帧/秒。本实验中智能手机的分辨率为5.454pixel/mm。作为对比，采用两个KEYENCE LK-G3001 CCD(Charge Coupled Devices)激光位移传感器监测建筑模型顶部的双向运动(在x方向和y方向上)，其采样频率为500Hz。使用LED照明，保证拍摄地点的照明均匀，最大限度地减少照明对测试结果的影响。利用冲击锤作用于缩尺模型顶部使其在x和y方向同时振动并录制10s视频。然后使用PME和RPME对录像视频进行识别和比较。图像处理中的监测点为建筑模型顶部的中心点。采用i5-8500 CPU@3.00GHz的个人计算机执行PME和RPME方法。由于相机与激光位移传感器数据集之间不可能进行时间同步，所以在数据中，时间序列是手工对齐的。
在RPME中使用本发明中提到的巴斯沃特理想带通滤波器将原始视频分解成单周期图像，对应的fs＝7。在PME中，当计算x方向运动时，设置Gabor参数中的θ＝0，当计算y方向运动时，设置Gabor参数中的θ＝π/2，波长λ都设置为80。识别结果的对比图如图9所示。使用MAE和STD评价运动识别的效果如图10所示。表1给出了两种方法的计算时间。
表1.不同算法运动识别的计算时间
从图9中可以看出，在0s附近，RPME和PME识别出的x、y方向结果与激光传感器的结果不同。这可能是由于测试环境引起的智能手机的振动，造成了图像识别的运动与参考之间的显著误差。在0.5s～10s范围内，RPME方法与参考点之间的时间历史曲线变化比传统PME方法的时间历史曲线变化小，从x和y两个方向的曲线峰值可以明显反映出这一点。这是因为传统的PME采用高斯包络正弦提取局部相位，而提出的RPME采用Riesz变换获得瞬时相位。因此，该方法有效地提高了运动识别的准确性
本试验的RPME和PME方法的MAE和STD见图10。在x、y方向上，RPME法的MAE分别比PME法低45％和40％。而两个方向的STD也比PME低。
RPME和PME方法的计算时间如表1所示。结果表明，PME方法的计算时间是RPME方法的6倍，即效率提高了83.41％。这是因为RPME方法可以同时获得结构在两个方向上的运动，而传统的PME方法需要两次计算才能获得结构在两个方向上的运动。此外，Riesz变换的实现是基于在频域内的点积关系，而PME的实现是实现图像和2D Gabor在空间域卷积，这也会导致显著增加计算成本。
最后需要说明的是，上述描述仅仅为本发明的优选实施例，本领域的普通技术人员在本发明的启示下，在不违背本发明宗旨及权利要求的前提下，可以做出多种类似的表示，这样的变换均落入本发明的保护范围之内。
Claims (7)
1.一种基于相位的工程结构双向动态位移测量方法，其特征在于，按照以下步骤进行：
S1、任选图像序列In(x,y,t)中的一张灰度图像Ir(x,y)，基于该灰度图像Ir(x,y)计算得到图像序列In(x,y,t)的谐振频率fs；
S2、利用二维快速傅里叶变换处理图像序列In(x,y,t)，使其由时域转移到频域，得到频率矩阵F(u,v,t)；
S3、使用中心频率为fs、带宽为1Hz的巴特沃思理想带通滤波器对频率矩阵F(u,v,t)在频率区间[fs-0.5,fs+0.5]进行滤波，得到滤波后的频率矩阵FB(u,v,t)；
S4、利用二维傅里叶逆变换处理滤波后的频率矩阵FB(u,v,t)，得到t时刻的单周期图像序列Is(x,y,t)；
S5、基于Riesz变换处理单周期图像Is(x,y,t)，同时计算评估图像中结构的双向运动；
所述步骤S5按照以下步骤进行：
S51、将t时刻的单周期图像Is(x,y,t)视为一个矩阵；
S52、对单周期图像Is(x,y,t)进行Riesz变换；
S53、根据步骤S52的Riesz变换结果，评估图像中结构同时在x方向和y方向运动的相位差；
S54、根据步骤S53得到的结构同时在x方向和y方向运动的相位差，建立结构双向运动与相位差之间的关系，从而获取结构在x方向和y方向的运动。
2.根据权利要求1所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S52中，对单周期图像Is(x,y,t)进行Riesz变换的公式为：
式(7)中，R[Is(x,y,t)]表示对单周期图像Is(x,y,t)进行Riesz变换，其中，Riesz算子R[.]＝(R1,R2)，R1表示x方向的Riesz变换结果，R2表示y方向的Riesz变换结果，i表示虚数，u＝[u1,u2]T表示频域信号维数，u1和u2分别表示单周期图像Is(x,y,t)在频域中与空域中的自变量x和y对应的变量，上标T表示转置运算，Is(u)表示t时刻的单周期图像Is(x,y,t)的傅里叶变换结果。
3.根据权利要求2所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S53中，将任意t时刻的单周期图像Is(x,y，t)及其Riesz变换结果(R1,R2)共同构成三元，并构建球坐标系，其中，将Is(x,y，t)记作I，得到如下方程：
I＝Acos(φ),R1＝Asin(φ)cos(θ),R2＝Asin(φ)sin(θ) (8)
式(8)中，A为振幅，θ为方向，φ为空间相位；
基于上述球坐标系，实轴I与空间向量(A,φ,θ)之间的夹角定义为空间相位φ，实轴I在R1轴和R2轴上的投影分别能够作为在x方向和y方向的瞬时相位；因此，为了计算空间相位φ，可以将三元(I,R1,R2)写成一个四元数q，其中q是由一个实部I和两个虚部(R1,R2)组成，公式如下：
q＝q0+iq1+jq2 (9)
式(9)中，其中q0＝I，q1＝R1，q2＝R2；故根据四元数表示法，i和j分别表示两个虚数单位，振幅A＝|q|；并且将相邻两帧图像分别表示为q＝q0+iq1+jq2和r＝r0+ir1+jr2，则相邻两帧图像的相位变化等于r/q，即:
Phase(φ(r/q))＝Phase(φ(r×q*/|q|2))＝Phase(φ(r×q*)) (10)
式(10)中，Phase(φ(r/q))表示以四元数形式表示的相邻两帧图像q和r的相邻图像相位差，q*表示q的共轭，其中对于r×q*表示为：
式(11)中，r0、r1和r2分别表示与相邻帧图像对应的I、R1和R2；
将r×q*写成一个四元数p，即：r×q*＝p＝p0+ip1+jp2，故基于式(11)可知：p0＝r0q0+r1q1+r2q2，p1＝-r0q1+r1q0，p2＝-r0q2+r2q0，四元数p的实部是Re(p)＝p0，虚部是Im(p)＝pv＝ip1+jp2，因此有：
式(12)中，pv表示p的虚部；
结合公式(10)和公式(12)可知，相邻两帧的相位变化Phase(φ(r/q))＝Phase(φ(r×q*))，也就是等于公式(12)的相位，由于只有虚部对计算相位有用，故：
因此，相邻两帧图像的相位差Phase(φ(r/q))为：
方向θ为：
在x方向和y方向的相位差Δφx和Δφy为：
4.根据权利要求3所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S54中，任意时刻t的单周期图像序列Is(x,y，t)可以表示为：
I(x,t)＝L1cos(ωx) (17)
式(17)中，用I(x,t)表示Is(x,y,t)，其中向量x＝[x,y]T表示图像像素灰度所在位置，L1表示t时刻图像的灰度幅度，ω＝[ωx,ωy]T表示图像的频率，ωx和ωy分别表示x方向和y方向的频率；
由式(7)和(17)可知，正弦函数L1cos(ωx)的Riesz变换对为：
正弦函数L1cos(ωx)的正交对Q(x,t)为：
Q(x,t)＝L1sin(ωx) (19)
与一维情况一致，故瞬时相位φ1(x,t)为:
φ1(x,t)＝ωx (20)
相同地，在x方向和y方向运动为δ＝[δx,δy]T的下一帧图像I(x+δ,t+Δt)可以表示为：
I(x+δ,t+Δt)＝L2cos(ωx) (21)
式(21)中，L2表示下一帧图像的灰度幅度，δx和δy分别表示为x方向和y方向的运动；
在x方向和y方向运动为δ＝[δx,δy]T的下一帧图像I(x+δ,t+Δt)的瞬时相位为：
φ2(x,t+Δt)＝ω(x+δ) (22)
故相位差Δφ为：
Δφ＝φ2(x,t+Δt)-φ1(x,t)＝ωδ (23)
式(23)中，Δφ＝[Δφx,Δφy]T,ω＝[ωx,ωy]T，δ＝[δx,δy]T，结合公式(16)和公式(23)可得：
Δφx＝ωxδx,Δφy＝ωyδy (24)
由式(24)可知，能够根据空间相位变化同时计算获取x方向运动δx和y方向运动δy。
5.根据权利要求1所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S1按照以下步骤进行：
S11、任选图像序列In(x,y,t)中的一张灰度图像Ir(x,y)；
S12、利用快速傅里叶变换处理灰度图像Ir(x,y)，使其由时域转移到频域，得到频率域Xr(u,v)；
S13、利用功率谱密度法处理频率序列Xr(u,v)，选择峰值点为谐振频率，fs＝(u2+v2)1/2。
6.根据权利要求5所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S12中，对灰度图像Ir(x,y)进行快速傅里叶变换得到的频率序列Xr(u,v)表达式为：
式(1)中，x和y分别表示选取位置在时域坐标系上的坐标，u和v分别表示空域中x和y在频域中对应的自变量，j表示虚部，M和N分别表示灰度图像Ir(x,y)的长度和宽度，Xr(u,v)表示灰度图像Ir(x,y)的频率域；
所述步骤S13中，利用功率谱密度法处理频率域Xr(u,v)得到以下表达式：
计算PSD(u,v)最大值对应的频率u和v，从而计算得到谐振频率fs＝(u2+v2)1/2。
7.根据权利要求1所述的一种基于相位的工程结构双向动态位移测量方法，其特征在于，所述步骤S2中，对图像序列In(x,y,t)进行二维快速傅里叶变换得到的频率矩阵F(u,v,t)表达式为：
式(3)中，x和y表示空域中的自变量，u和v表示空域中x和y在频域中对应的自变量，j表示虚部；
所述步骤S3中，滤波后的频率矩阵FB(u,v,t)与频率矩阵F(u,v,t)的关系式为：
FB(u,v,t)＝F(u,v,t)*BH(u,v) (4)
式(4)中，BH(u,v)表示频域中的巴斯沃特思理想带通滤波器，其中BH(u,v)的表达式为：
式(5)中，D(u,v)＝[(u-M/2)2+(v-N/2)2]/2，D(u,v)表示在频域中的点(u,v)到圆心的距离，M和N分别表示图像的长度和宽度，D1＝fs+0.5,D2＝fs-0.5；
所述步骤S4中，对滤波后的频率矩阵FB(u,v,t)进行二维傅里叶逆变换得到的单周期图像Is(x,y,t)的表达式为：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210214636.0A CN114564685A (zh) | 2022-03-07 | 2022-03-07 | 一种基于相位的工程结构双向动态位移测量方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210214636.0A CN114564685A (zh) | 2022-03-07 | 2022-03-07 | 一种基于相位的工程结构双向动态位移测量方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114564685A true CN114564685A (zh) | 2022-05-31 |
Family
ID=81717118
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210214636.0A Pending CN114564685A (zh) | 2022-03-07 | 2022-03-07 | 一种基于相位的工程结构双向动态位移测量方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114564685A (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2002010788A1 (en) * | 2000-07-31 | 2002-02-07 | Koninklijke Philips Electronics N.V. | Magnetic resonance method for forming a fast dynamic image |
| US20140153692A1 (en) * | 2012-11-30 | 2014-06-05 | Canon Kabushiki Kaisha | Combining Differential Images by Inverse Riesz Transformation |
| EP2770322A1 (en) * | 2013-02-26 | 2014-08-27 | C.R.F. Società Consortile per Azioni | Method and system for detecting defects in painting of components, in particular of motor-vehicle bodies |
| CN109949251A (zh) * | 2019-04-03 | 2019-06-28 | 天津中航亿达科技有限公司 | 一种基于黎兹金字塔的大气湍流退化视频快速稳定方法 |
| CN110068388A (zh) * | 2019-03-29 | 2019-07-30 | 南京航空航天大学 | 一种基于视觉和盲源分离的振动检测方法 |
| GB202103298D0 (en) * | 2021-03-10 | 2021-04-21 | Univ College Cardiff Consultants Ltd | Method and Apparatus for Non-Destructive Testing |
| CN113076517A (zh) * | 2021-04-01 | 2021-07-06 | 重庆大学 | 基于希尔伯特变换的土木工程结构动态监测相位评估方法 |
-
2022
- 2022-03-07 CN CN202210214636.0A patent/CN114564685A/zh active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2002010788A1 (en) * | 2000-07-31 | 2002-02-07 | Koninklijke Philips Electronics N.V. | Magnetic resonance method for forming a fast dynamic image |
| US20140153692A1 (en) * | 2012-11-30 | 2014-06-05 | Canon Kabushiki Kaisha | Combining Differential Images by Inverse Riesz Transformation |
| EP2770322A1 (en) * | 2013-02-26 | 2014-08-27 | C.R.F. Società Consortile per Azioni | Method and system for detecting defects in painting of components, in particular of motor-vehicle bodies |
| CN110068388A (zh) * | 2019-03-29 | 2019-07-30 | 南京航空航天大学 | 一种基于视觉和盲源分离的振动检测方法 |
| CN109949251A (zh) * | 2019-04-03 | 2019-06-28 | 天津中航亿达科技有限公司 | 一种基于黎兹金字塔的大气湍流退化视频快速稳定方法 |
| GB202103298D0 (en) * | 2021-03-10 | 2021-04-21 | Univ College Cardiff Consultants Ltd | Method and Apparatus for Non-Destructive Testing |
| CN113076517A (zh) * | 2021-04-01 | 2021-07-06 | 重庆大学 | 基于希尔伯特变换的土木工程结构动态监测相位评估方法 |
Non-Patent Citations (4)
| Title |
|---|
| ARAL SARRAFI 等: "Uncertainty Quantification of Phase-Based Motion Estimation on Noisy Sequence of Images", 《HEALTH MONITORING OF STRUCTURAL AND BIOLOGICAL SYSTEMS》, 31 December 2017 (2017-12-31), pages 1 - 8 * |
| J LI等: "Whole-field thickness strain measurement using multiple camera digital image correlation system", 《OPTICS AND LASERS IN ENGINEERING》, vol. 90, 31 March 2017 (2017-03-31), pages 19 - 25, XP055839406, DOI: 10.1016/j.optlaseng.2016.09.012 * |
| M.Z. LI 等: "Two-dimensional motion estimation using phase-based image processing with Riesz transform", 《MECHANICAL SYSTEMS AND SIGNAL PROCESSING》, vol. 188, 1 April 2023 (2023-04-01), pages 1 - 20 * |
| 周根: "基于视觉放大光流跟踪的微幅机械振动测量方法研究", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》, no. 06, 15 June 2021 (2021-06-15), pages 029 - 48 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chen et al. | Modal identification of simple structures with high-speed video using motion magnification | |
| Liu et al. | Structural motion estimation via Hilbert transform enhanced phase-based video processing | |
| Miao et al. | Phase-based displacement measurement on a straight edge using an optimal complex Gabor filter | |
| CN111784647B (zh) | 基于视频振动放大的高精度结构模态测试方法 | |
| Chen et al. | Cable force determination using phase-based video motion magnification and digital image correlation | |
| Li et al. | Two-dimensional motion estimation using phase-based image processing with Riesz transform | |
| Fu et al. | Non-contact optical dynamic measurements at different ranges: a review | |
| Wang et al. | Measurement of sinusoidal vibration from motion blurred images | |
| Zang et al. | Phase-based vibration frequency measurement from videos recorded by unstable cameras | |
| CN114862809A (zh) | 一种基于移动终端与图像处理的振动监测方法和装置 | |
| CN104897065A (zh) | 一种板壳结构表面位移场的测量系统 | |
| Cosco et al. | Towards phase-based defect detection: A feasibility study in vibrating panels | |
| Li et al. | Dynamic characteristics identification of an arch dam model via the phase-based video processing | |
| Shang et al. | Multi-point vibration measurement for mode identification of bridge structures using video-based motion magnification | |
| CN110532725B (zh) | 基于数字图像的工程结构力学参数识别方法及系统 | |
| CN114564685A (zh) | 一种基于相位的工程结构双向动态位移测量方法 | |
| CN113076517B (zh) | 基于希尔伯特变换的土木工程结构动态监测相位评估方法 | |
| Thiyagarajan et al. | Implementation of video motion magnification technique for non-contact operational modal analysis of light poles | |
| Wang et al. | Real-time vibration visualization using GPU-based high-speed vision | |
| Shimo et al. | Development of dynamic shape and strain measurement system by sampling moire method | |
| Zhao et al. | Precise position detection of linear motor movers based on extended joint transformation correlation | |
| Peng et al. | Full-field visual vibration measurement of rotating machine under complex conditions via unsupervised retinex model | |
| Zhong et al. | Three-dimensional translation vibration measurement system based on linear array sensor and composite fringe pattern | |
| Shen et al. | Video-based vibration measurement for large structure: A spatiotemporal disturbance-adaptive morphological component analysis | |
| Holak | A motion magnification application in video-based vibration measurement |
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 |