一种基于相位的工程结构双向动态位移测量方法

1.本发明涉及土木工程结构动态监测技术领域，具体涉及一种基于相位的工程结构双向动态位移测量方法。


背景技术：

2.结构健康监测(简称shm)能够实现低成本条件的维护，执行结构状态预测，并消除意外灾难性故障。用于shm的标准有线加速度计是劳动密集型的仪器结构，非接触式测量方法避免了这些缺点。数码摄像机作为一种非接触式监测方法，因其成本较低，操作灵活，并能够提供非常高的空间分辨率的优势成为了近期研究的热点。结合图像处理方法、基于视频的测量已经成功地用于各种类型结构的振动测量。最近，shm中引入了数字图像相关(简称dic)和点跟踪(简称dpt)技术。
3.然而，这些方法通常需要在结构表面预制散斑模式图或高对比度的特征点标记，并使用基于图像灰度的相关性来估计结构运动场。一方面，这些要求提出了表面预处理的问题，特别是当测量结构面积大或环境复杂时这一要求很难达到。另一方面，基于图像灰度本身的相关性对环境噪声和干扰非常敏感，如场景光照的突然变化会导致结构运动场的估计产生较大的误差。
4.因此，相位运动估计法(简称pme)作为图像估计的一种，无需任何表面预处理，同时该方法是基于图像的相位来实现结构的运动估计的，而相位行为变化对光照变化相对不敏感(如阴影、高光等)。因此，相位轮廓确实提供了一个很好的近似运动场。
5.aral sarrafi等人在图像强度上应用一维gabor小波来提取结构运动。因为gabor小波是由高斯包络的复正弦函数组成，因此使用gabor小波获取的相位为局部相位，这可能导致结构的运动估计精度不足。为了提高计算精度，本发明人团队曾在中国专利申请cn113076517a中提出了一种基于希尔伯特变换的相位运动估计法，该方法利用hilbert变换代替gabor小波获取瞬时结构的相位以提高原始pme的精度。之后，aral sarrafi等人对gabor小波进行扩展，应用2d gabor波提取图像中结构的局部相位变化，并使用2.3m长风力机叶片的损伤来验证方法的损伤检测能力。虽然采用了二维gabor小波估计运动，但上述这些pme方法都集中在一维相位估计上。也就是说，一次只能考虑结构的一个方向运动，但图像通常是二维的，所以这些pme方法至少在理论上是不完善的，难以适应于动态结构监测。并且，在土木工程或机械工程中，结构的双向运动是常见的要求，例如高层建筑在地震作用下的强弱轴向位移就是双向运动，导致现有的这些pme方法都不适用。
6.解决以上问题成为当务之急。


技术实现要素：

7.为解决以上的技术问题，本发明提供了一种基于相位的工程结构双向动态位移测量方法。
8.其技术方案如下：
9.一种基于相位的工程结构双向动态位移测量方法，其要点在于，按照以下步骤进行：
10.s1、任选图像序列in(x,y,t)中的一张灰度图像ir(x,y)，基于该灰度图像ir(x,y)计算得到图像序列in(x,y,t)的谐振频率fs；
11.s2、利用二维快速傅里叶变换处理图像序列in(x,y,t)，使其由时域转移到频域，得到频率矩阵f(u,v,t)；
12.s3、使用中心频率为fs、带宽为1hz的巴特沃思理想带通滤波器对频率矩阵f(u,v,t)在频率区间[f
s-0.5,fs 0.5]进行滤波，得到滤波后的频率矩阵fb(u,v,t)；
[0013]
s4、利用二维傅里叶逆变换处理滤波后的频率矩阵fb(u,v,t)，得到t时刻的单周期图像序列is(x,y,t)；
[0014]
s5、基于riesz变换处理单周期图像is(x,y,t)，同时计算评估图像中结构的双向运动；
[0015]
所述步骤s5按照以下步骤进行：
[0016]
s51、将t时刻的单周期图像is(x,y,t)视为一个矩阵；
[0017]
s52、对单周期图像is(x,y,t)进行riesz变换；
[0018]
s53、根据步骤s52的riesz变换结果，评估图像中结构同时在x方向和y方向运动的相位差；
[0019]
s54、根据步骤s53得到的结构同时在x方向和y方向运动的相位差，建立结构双向运动与相位差之间的关系，从而获取结构在x方向和y方向的运动。
[0020]
与现有技术相比，本发明的有益效果：
[0021]
采用以上技术方案的一种基于相位的工程结构双向动态位移测量方法，将riesz变换引入到相位运动估计中，提出了一种可以同时获取两个方向运动的视频运动估计方法，解决了现有pme方法中的局限，使计算效率得到大幅提升，更适合于动态结构监测，为民用基础设施的动态监测提供了一种简单实用的新方法。
附图说明
[0022]
图1为频域中riesz变换算子的示意图；
[0023]
图2为输入信号i与其riesz变换(r1,r2)的组合示意图；
[0024]
图3为采用巴特沃思理想带通滤波器进行滤波的示意图；
[0025]
图4在对比例一的其中一个振荡位置处的抛物面示意图；
[0026]
图5为对比例一的使用不同参数的pme方法和本发明方法(rpme)识别的x方向运动与真实运动对比图；
[0027]
图6为对比例一的使用不同参数的pme方法和本发明方法(rpme)识别的y方向运动与真实运动对比图；
[0028]
图7为对比例一的采用pme方法和本发明方法(rpme)在x方向和y方向运动识别效果的平均绝对误差和标准偏差的对比图；
[0029]
图8为对比例二的高层建筑缩尺模型实验图；
[0030]
图9为对比例二的使用不同参数的pme方法和本发明方法(rpme)的运动识别曲线对比图((a)为x方向，(b)为y方向)；
[0031]
图10为对比例二的采用pme方法和本发明方法(rpme)在x方向和y方向运动识别效果的平均绝对误差和标准偏差的对比图。
具体实施方式
[0032]
以下结合实施例和附图对本发明作进一步说明。
[0033]
一种基于相位的工程结构双向动态位移测量方法，按照以下步骤进行：
[0034]
s1、任选图像序列in(x,y,t)中的一张灰度图像ir(x,y)，基于该灰度图像ir(x,y)计算得到图像序列in(x,y,t)的谐振频率fs；
[0035]
具体地说，步骤s1按照以下步骤进行：
[0036]
s11、任选图像序列in(x,y,t)中的一张灰度图像ir(x,y)；
[0037]
s12、利用快速傅里叶变换处理灰度图像ir(x,y)，使其由时域转移到频域，得到频率域xr(u,v)，xr(u,v)表达式为：
[0038][0039]
式(1)中，x和y分别表示选取位置在空域坐标系上的坐标，u和v分别表示空域中x和y在频域中对应的自变量，j表示虚部，m和n分别表示灰度图像ir(x,y)的长度和宽度，xr(u,v)表示灰度图像ir(x,y)的频率域；
[0040]
s13、利用功率谱密度法处理频率序列xr(u,v)，得到以下表达式：
[0041][0042]
计算psd(u,v)最大值对应的频率u和v，从而计算得到谐振频率fs＝(u2 v2)
1/2
，即：选择峰值点为谐振频率fs＝(u2 v2)
1/2
。
[0043]
s2、利用二维快速傅里叶变换处理图像序列in(x,y,t)，使其由时域转移到频域，得到频率矩阵f(u,v,t)，f(u,v,t)的表达式为：
[0044][0045]
式(3)中，x和y表示空域中的自变量，u和v表示空域中x和y在频域中对应的自变量，j表示虚部；
[0046]
s3、请参见图3，使用中心频率为fs、带宽为1hz的巴特沃思理想带通滤波器对频率矩阵f(u,v,t)在频率区间[f
s-0.5,fs 0.5]进行滤波，得到滤波后的频率矩阵fb(u,v,t)。
[0047]
具体地说，滤波后的频率矩阵fb(u,v,t)与频率矩阵f(u,v,t)的关系式为：
[0048]
fb(u,v,t)＝f(u,v,t)*bh(u,v)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0049]
式(4)中，bh(u,v)表示频域中的巴斯沃特思理想带通滤波器，其中bh(u,v)的表达式为：
[0050][0051]
式(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＝f
s-0.5。
[0052]
s4、利用二维傅里叶逆变换处理滤波后的频率矩阵fb(u,v,t)，得到t时刻的单周期图像序列is(x,y,t)，is(x,y,t)的表达式为：
[0053][0054]
s5、基于riesz变换处理单周期图像is(x,y,t)，同时计算评估图像中结构的双向运动；
[0055]
具体地说，步骤s5按照以下步骤进行：
[0056]
s51、将t时刻的单周期图像is(x,y,t)视为一个矩阵；
[0057]
s52、对单周期图像is(x,y,t)进行riesz变换；
[0058]
对单周期图像is(x,y,t)进行riesz变换的公式为：
[0059][0060]
式(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)的傅里叶变换结果。
[0061]
如图1所示的一个riesz算子，其中在求关于u1轴的riesz变换时，u2＝0，通过符号sign(u1)＝u1/|u
12
0|
1/2
得到，这可以通过沿图1中的u2＝0等高线来可视化。
[0062]
s53、根据步骤s52的riesz变换结果，评估图像中结构同时在x方向和y方向运动的相位差；
[0063]
步骤s53中，将任意t时刻的单周期图像is(x,y，t)及其riesz变换结果(r1,r2)共同构成三元，并构建球坐标系，如图2所示，同时将is(x,y，t)记作i，可以写成如下方程：
[0064]
i＝acos(φ),r1＝asin(φ)cos(θ),r2＝asin(φ)sin(θ)
ꢀꢀꢀꢀ
(8)
[0065]
式(8)中，a为振幅，θ为方向，φ为空间相位。
[0066]
从图2可以看出，基于上述球坐标系，实轴i与空间向量(a,φ,θ)之间的夹角定义为空间相位φ，实轴i在r1轴和r2轴上的投影分别能够作为在x方向和y方向的瞬时相位；因此，为了计算空间相位φ，可以将三元(i,r1,r2)写成一个四元数q，其中q是由一个实部i和两个虚部(r1,r2)组成，公式如下：
[0067]
q＝q0 iq1 jq2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0068]
式(9)中，其中q0＝i，q1＝r1，q2＝r2；故根据四元数表示法，i和j分别表示两个虚数单位，振幅a＝|q|；并且将相邻两帧图像分别表示为q＝q0 iq1 jq2和r＝r0 ir1 jr2，则相邻两帧图像的相位变化等于r/q，即:
[0069]
phase(φ(r/q))＝phase(φ(r
×q*
/|q|2))＝phase(φ(r
×q*
))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0070]
式(10)中，phase(φ(r/q))表示以四元数形式表示的相邻两帧图像q和r的相邻图像相位差，q
*
表示q的共轭，其中对于r
×q*
表示为：
[0071][0072]
式(11)中，r0、r1和r2分别表示与相邻帧图像对应的i、r1和r2；
[0073]
将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，因此有：
[0074][0075]
式(12)中，pv表示p的虚部；
[0076]
结合公式(10)和公式(12)可知，相邻两帧的相位变化phase(φ(r/q))＝phase(φ(r
×q*
))，也就是等于公式(12)的相位，由于只有虚部对计算相位有用，故：
[0077][0078]
因此，相邻两帧图像的相位差phase(φ(r/q))为：
[0079][0080]
方向θ为：
[0081][0082]
在x方向和y方向的相位差δφ
x
和δφy为：
[0083][0084]
s54、根据步骤s53得到的结构同时在x方向和y方向运动的相位差，建立结构双向运动与相位差之间的关系，从而获取结构在x方向和y方向的运动。
[0085]
步骤s54中，任意时刻t的单周期图像序列is(x,y，t)可以表示为：
[0086]
i(x,t)＝l1cos(ωx)
ꢀꢀꢀꢀ
(17)
[0087]
式(17)中，用i(x,t)表示is(x,y，t)，其中向量x＝[x,y]
t
表示图像像素灰度所在位置，l1表示t时刻图像的灰度幅度，ω＝[ω
x
,ωy]
t
表示图像的频率，ω
x
和ωy分别表示x方向和y方向的频率；
[0088]
由式(7)和(17)可知，正弦函数l
1 cos(ωx)的riesz变换对为：
[0089][0090]
正弦函数l
1 cos(ωx)的正交对q(x,t)为：
[0091]
q(x,t)＝l1sin(ωx)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0092]
与一维情况一致，故瞬时相位φ1(x,t)为:
[0093]
φ1(x,t)＝ωx
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0094]
相同地，在x方向和y方向运动为δ＝[δ
x
,δy]
t
的下一帧图像i(x δ,t δt)可以表示为：
[0095]
i(x δ,t δt)＝l2cos(ωx)
ꢀꢀꢀꢀꢀꢀ
(21)
[0096]
式(21)中，l2表示下一帧图像的灰度幅度，δ
x
和δy分别表示为x方向和y方向的运
动；
[0097]
在x方向和y方向运动为δ＝[δ
x
,δy]
t
的下一帧图像i(x δ,t δt)的瞬时相位为：
[0098]
φ2(x,t δt)＝ω(x δ)
ꢀꢀꢀꢀꢀ
(22)
[0099]
故相位差δφ为：
[0100]
δφ＝φ2(x,t δt)-φ1(x,t)＝ωδ
ꢀꢀꢀꢀ
(23)
[0101]
式(23)中，δφ＝[δφ
x
,δφy]
t
,ω＝[ω
x
,ωy]
t
，δ＝[δ
x
,δy]
t
，结合公式(16)和公式(23)可得：
[0102]
δφ
x
＝ω
x
δ
x
,δφy＝ωyδyꢀꢀꢀꢀꢀ
(24)
[0103]
由式(24)可知，能够根据空间相位变化计算得到x方向运动δ
x
和y方向运动δy。
[0104]
对比例一：有限元实验
[0105]
采用一个同时沿x方向和y方向移动的二维抛物曲面来模拟一组图像/或视频：
[0106]
z(x,y,t)＝-((x-δ
x
(t))2 (y-δy(t))2) 120
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0107]
式(25)中，z(x,y,t)表示二维抛物曲面，δ
x
(t)和δy(t)分别是任意的运动位移,在这里请参见图4,采用运动函数δ
x
(t)＝e-ξωnt
sin(ωnt)来模拟粘性阻尼振荡的抛物曲面表面沿x方向和y方向的轴运动，其中沿x轴的运动设置为ωn＝15和ξ＝0.09,沿y轴的运动设置为ωn＝4和ξ＝0.3，其中，ωn和ξ表示运动函数中的自变量。
[0108]
采用本方法(简称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的识别结果。为了定量地比较这两种方法，使用平均绝对误差(mean absolute error,简称mae)和标准偏差(standard deviation,简称std)来评价两种识别方法的效果，其定义如下:
[0109][0110]
式(26)中，obverved
t
表示第t个计算点的预估运动值，true
t
表示第t个计算点的真实运动值，n表示计算点的总数。
[0111]
使用mae和std来评价rpme和pme方法的识别效果如附图7所示。
[0112]
从图5和图6可以看出，本发明方法提出的rpme方法在提取x方向和y方向的运动信息时，识别效果最好，最大误差在1pixel以内。这是因为在riesz变换中采用了瞬时相位，riesz变换是希尔伯特的二维扩展。因此，该riesz变换提高了所提方法的识别精度。而pme方法则是基于局部高斯包络的正弦信号来提取图像的局部相位，因此pme方法识别的运动信息的精度不如rpme方法高。此外，波长的选择也降低了pme方法的识别精度。当选择合适的波长时，识别效果可以达到合理的程度，如图5(a)和6(a)所示，否则可能会产生较大的误差，如图5(b)和6(b)所示。换句话说，使用pme的一个限制是它需要预先知道要测量的对象。另一方面，与每次只考虑图像单一方向相位的pme相比，本发明方法提出的rpme方法提出的
空间相位同时识别两个方向的相位更加有效。
[0113]
由图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方法。
[0114]
对比例二：高层建筑缩尺模型实验
[0115]
如图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方法。由于相机与激光位移传感器数据集之间不可能进行时间同步，所以在数据中，时间序列是手工对齐的。
[0116]
在rpme中使用本发明中提到的巴斯沃特理想带通滤波器将原始视频分解成单周期图像，对应的fs＝7。在pme中，当计算x方向运动时，设置gabor参数中的θ＝0，当计算y方向运动时，设置gabor参数中的θ＝π/2，波长λ都设置为80。识别结果的对比图如图9所示。使用mae和std评价运动识别的效果如图10所示。表1给出了两种方法的计算时间。
[0117]
表1.不同算法运动识别的计算时间
[0118][0119]
从图9中可以看出，在0s附近，rpme和pme识别出的x、y方向结果与激光传感器的结果不同。这可能是由于测试环境引起的智能手机的振动，造成了图像识别的运动与参考之间的显著误差。在0.5s～10s范围内，rpme方法与参考点之间的时间历史曲线变化比传统pme方法的时间历史曲线变化小，从x和y两个方向的曲线峰值可以明显反映出这一点。这是因为传统的pme采用高斯包络正弦提取局部相位，而提出的rpme采用riesz变换获得瞬时相位。因此，该方法有效地提高了运动识别的准确性
[0120]
本试验的rpme和pme方法的mae和std见图10。在x、y方向上，rpme法的mae分别比pme法低45％和40％。而两个方向的std也比pme低。
[0121]
rpme和pme方法的计算时间如表1所示。结果表明，pme方法的计算时间是rpme方法的6倍，即效率提高了83.41％。这是因为rpme方法可以同时获得结构在两个方向上的运动，而传统的pme方法需要两次计算才能获得结构在两个方向上的运动。此外，riesz变换的实现是基于在频域内的点积关系，而pme的实现是实现图像和2d gabor在空间域卷积，这也会导致显著增加计算成本。
[0122]
最后需要说明的是，上述描述仅仅为本发明的优选实施例，本领域的普通技术人员在本发明的启示下，在不违背本发明宗旨及权利要求的前提下，可以做出多种类似的表示，这样的变换均落入本发明的保护范围之内。