基于数字图像相关的古建筑彩绘梁损伤识别方法

刘纲，顾嘉伟，李孟珠，贺成华

(重庆大学 a.山地城镇建设与新技术教育部重点实验室; b.土木工程学院，重庆 400045)

木结构是中国古建筑的主要结构形式，是传承中华灿烂文明的重要载体。受木材老化、环境及生物侵蚀的长期作用和突发灾害的过载作用，古建筑木结构普遍存在结构性能退化的问题，并伴有不同程度的残损[1]。但作为中国特有的历史文化遗产和宝贵的旅游资源，不少古建筑仍需对游客开放，长期“带病”工作可能加剧木结构的损伤。因此，利用先进技术手段实时掌握木结构的安全状态，是实现古建筑木结构预防性保护的根本前提[2-3]。

随着传感及信息技术的快速发展，结构监测及损伤识别技术已在航空航天、机械及土木工程结构中得到广泛应用，其通过传感设备实时获取结构响应和环境信息，利用损伤识别和安全评估及时掌握结构的安全状态，从而可有效预防结构安全事故[4-5]。近年来，该技术被逐步引入古建筑木结构领域，Jiang等[6]利用光纤光栅传感技术对古建筑木结构的梁偏转、柱倾斜、关节错位等进行监测，并开发了适用于穿斗式木结构的监测系统。杨娜等[7]基于某藏式古建木结构监测系统，采用统计过程控制理论对监测数据进行异常诊断，结果表明，游客数量对木结构的稳定存在较大影响。薛建阳等[8]通过殿堂式古建筑缩尺模型振动台试验与西安钟楼简化力学模型，提出了柱脚滑移状态下的结构等效刚度识别方法，利用速度、加速度传感器的响应识别出结构层间抗侧刚度值，结果表明，该方法可在地震倒塌前发出预警。Li等[9]在两层木框架模型的6个梁柱节点处分别布设加速度传感器，在随机激励下利用小波单谱熵指标识别出了木结构的损伤。王鑫等[10]建立了木结构框架模型，对环境激励下梁上各节点的加速度响应信号进行小波包分解，采用小波包能量曲率差指标实现了木梁的损伤定位。

但以上监测方法均需在古建筑木结构中安装传感器，对木结构自身有一定损伤，不利于古建筑木结构的保护。近十年来，随着摄像技术的快速发展，基于数字图像相关(digital image correlation，DIC)的监测技术得到长足发展。该技术采用数字相机非接触拍摄被测对象表面图像，通过图像的像素匹配获得被测对象的全场变形信息，进而提取应变、模态等测试参量，具有非接触、全场测量及便于实施等优势，并采用人工喷涂散斑等方式增密被测对象的表面特征，以大幅提高测量精度，目前已在现代结构变形、振动监测中展示了巨大潜力[11]。Yu等[12]在悬臂铝板上喷涂散斑点，利用高速相机进行全场振动测量，准确获取了被测物体的自振频率、振型等模态参数；Khadka等[13]在风力涡轮机叶片表面制作散斑点，利用DIC方法辨识了叶片在旋转状态下的模态参数；Molina-Viedma等[14]以一榀框架结构为研究对象，利用DIC方法提取了随机振动激励下结构的工作模态，采用局部模态滤波方法识别和定位了框架梁的损伤。

DIC测量方法可借助被测物体表面特定的灰度特征实现变形、模态参数等测试，而古建筑木结构表面往往绘制了彩色图案。为满足古建筑木结构保护要求，笔者直接利用木结构表面自身彩绘图案，探索采用DIC非接触式测量方法辨识木梁工作模态的可行性。在此基础上，选取转角模态和曲率模态作为损伤指标，利用DIC测试的工作模态沿梁长方向连续的特征，实现木结构早期损伤的准确定位。数值模拟和实验结果表明，DIC方法可辨识木梁工作模态并实现微小损伤定位，从而为古建筑木结构预防性保护提供新的监测手段。

1 基于DIC原理的模态参数辨识方法

1.1 DIC位移测试原理

数字图像相关法(DIC)是近年发展起来的基于视觉的新测试技术。以二维图像位移解算为例，DIC方法先建立参考图像(未变形)与目标图像(变形后)统一的x-y坐标系，再进行整、亚像素位移解算。

整像素解算阶段，在参考图像中选取任意像素作为计算点P(x0,y0)，并以计算点为中心选取(2M+1)×(2M+1)像素大小的矩形区域为参考子集，其中M为参考子集的半径，如图1所示；然后，在目标图像中选取任意点为待测点P′(x0′,y0′)，以该点为中心形成与参考子集相同尺寸的目标子集，并计算相关系数CZNSSD。

图1 基于子集相关的整像素搜索原理Fig.1 Whole pixel search principle based on subset

(1)

在亚像素计算阶段，利用已获得的整像位移值，通过插值和亚像素搜索方法进行亚像素位移解算。对于二维图像，需要在x、y方向分别进行插值。其中常用的三次样条插值函数为[16]

∇f(xi,yi)(x′-xi)+∇f(xi,yi) (i=1,…,4)

(2)

∇′f(xi,yi)+∇′f(xi+1,yi+1)

2∇′f(xi,yi)-∇′f(xi+1,yi+1)

(3)

式(2)表示在x方向的插值函数，其中，i表示插值区域的第i行；(xi,yi)为插值区域整像素点的绝对坐标；∇f(xi,yi)为整像素点(xi,yi)的灰度梯度。插值过程：首先，以待插值目标点(x′,y′)为中心选取4×4像素大小的插值区域，如图2所示；然后利用第i行(i=1…4)4个整像素点，插值得到该行亚像素辅助点(x′,yi)的灰度值f(x′,yi)；最后，利用4个亚像素辅助点的灰度值，在y方向进行三次样条插值，得到待插值目标点的亚像素灰度值f(x′,y′)。

图2 双三次样条插值原理Fig.2 Principle of bicubic spline

最后，选择常用的逆向组合高斯牛顿法进行亚像素位移搜索解算，基本步骤为[17]：

1)首次迭代时，位移映射参数p=[u,ux,uy,v,vx,vy]T中，u、v分别为整像素位移u0、v0，一阶位移梯度ux=∂u/∂x、uy=∂u/∂y、vx=∂v/∂x、vy=∂v/∂y均为零。通过式(4)～式(5)计算子集映射函数W(x,y;p)与Hessian矩阵H。

(4)

(5)

式中：Δx、Δy为参考子集的中心点P与子集内部任意一点Q的坐标差，如图1所示；∇f为灰度梯度。

2)利用亚像素级相关函数计算位移增量映射参数Δp。

(6)

将目标子集映射函数乘以增量映射函数的逆，获得新映射函数W′(x,y;p)。

W′(x,y;p)=W(x,y;p)·W-1(x,y;Δp)

(7)

3)若位移增量映射参数Δp=[Δu,Δux,Δuy,Δv,Δvx,Δvy]T不满足收敛条件，则将新映射函数W′(x,y;p)代入步骤2)继续迭代。收敛条件为[18]

(8)

若满足收敛条件，则输出新映射函数W′(x,y;p)中的p。

4)利用整像素位移与位移映射参数p，通过式(9)获得目标图像匹配点P′(x0′,y0′)的坐标值。

(9)

1.2 工作模态辨识方法

当相机采集帧率满足那奎斯特采样定理时[19]，利用DIC方法解算出一系列图片中目标点的位移时程曲线；再通过傅里叶变换得到各像素点位移时程的幅频、相频曲线；然后根据各阶固有频率处幅频曲线的幅值大小、相频曲线的相位角判断各像素点相对位置，从而获得被测结构工作模态。

多自由度体系在模态空间下的运动微分方程可表示为[20]

(10)

式中：qi为第i阶正规坐标；mi、ci和ki分别为第i阶(i=1,2...n)正规坐标对应的广义质量、阻尼和刚度；f(t)为外荷载向量，t为时间；Φ为模态振型矩阵，上标T表示矩阵转置；n为体系的模态阶数。

令f(t)=Fejw t，则qi=Qiejw t，式(10)可转换为频域方程

(11)

由式(11)可得系统的频响函数H(w)=ΦiQi/F，表示系统位移响应与外部激励频率之间的关系，例如，在r点的y方向进行激励时，引起s点的y方向位移响应为

(12)

由式(12)可得结构e个测点的位移频响函数

(13)

由式(13)可知，频响函数矩阵中任意一行或一列包含所有的模态参数，例如，式(13)中第1列表示在1号测点y方向激励，能够获得1～e号测点的y方向位移响应。因此，可通过两种方式进行工作模态的辨识：其一，在结构的某一固定点进行激振，在其他各点进行拾振，即可求得H的一列；其二，在结构的某一固定点拾振，在其他各点激振，即可求得H的一行。

在固定点进行激励，利用DIC方法识别出各测点的位移响应后，通过傅里叶变换获取各测点在同一固有频率下幅频曲线的幅值大小并进行归一化处理，再利用各测点相频曲线的相位角判断其相对位置，即可实现工作模态的辨识。

2 基于模态参数的损伤识别方法

采用模态参数进行损伤识别的原理在于，损伤将造成结构质量、刚度改变，进而导致结构模态参数发生变化，且在损伤处模态振型往往会出现奇异性峰值，从而判断结构损伤部位[21]。大量研究表明[22-24]，转角模态和曲率模态对结构早期损伤较为敏感，故选取这两个指标进行损伤识别和定位。

2.1 转角模态

转角模态是位移模态的一阶导数。因结构损伤处刚度减少，导致刚度变化截面处左右侧转角不同，但结构仍满足整体变形协调条件[25]，即

(14)

式中：xl、xr分别为截面的左、右侧；v(x)为y方向的位移模态；EI为梁的抗弯刚度；Kτ为刚度变化截面处的等效转动刚度，具体公式见文献[26]。式(14)表明，在刚度变化截面处位移模态的一阶导数不连续，即转角模态将出现畸变，从而可通过该特性识别梁截面发生的损伤。

2.2 曲率模态

曲率模态为位移模态的二阶导数，一般通过中心差分法近似计算[27]。

(15)

式中：ρ(x)为梁弯曲振动曲线的曲率半径；Δ为测点间距。当结构局部出现损伤或损伤程度增加时，梁损伤部位的抗弯刚度减少，故损伤位置处的曲率将增大，因而通过检测曲率模态在损伤前后的变化即可判断损伤程度与损伤位置。

3 数值模拟

3.1 数值模型

建立2 m单跨简支梁有限元模型，验证密集测点条件下损伤识别的有效性，矩形梁截面宽×高为5 cm×1 cm，如图3所示。在梁底部设置宽1 mm的横向通长裂纹模拟损伤，以裂纹深度h与梁高H的比值表示损伤程度，以左端支座为零点表示损伤位置。损伤工况如表1所示。

图3 简支梁有限元模型Fig.3 Finite element model of simply supported

表1 损伤工况Table 1 Damage cases of beam

为对比密集、稀疏测点下损伤识别效果，假设在简支梁上以0.25 m等间距布设加速度传感器为稀疏测点情况，相当于在被测结构1/8等分点布设传感器；而DIC方法能实现的测点最小间距由相机像素和拍摄物距决定，当前单台主流相机针对2 m长木梁能够识别出的物理距离(相当于测点间距)可高达0.34 mm，但考虑到计算工作量及必要性，以0.05 m作为DIC测试的最大间距(密集测点)。若0.05 m间距无法实现预期的损伤识别，可进一步选取更小的测点间距进行试算。梁的损伤模型及传感器测点布置如图3所示。在数值算例中，直接采用有限元计算的模态振型进行损伤识别，没有采用DIC进行模态振型识别。

3.2 损伤识别结果

考虑到古建筑木结构梁具有截面大、跨度小的特点，实际结构中往往无法有效对高阶振型进行激励，故选取稀疏、密集测点下简支梁的第一阶模态振型，其工况1～工况4下模态振型的变化如图4所示，其中工况0表示无损伤状态。稀疏测点下的转角、曲率模态在损伤前后无明显差别，而密集测点获得的转角模态与曲率模态在损伤位置处存在明显的奇异性峰值，可以判断损伤位置与损伤程度。

图4 基于第一阶振型的转角及曲率模态Fig.4 Rotation and curvature mode based on

裂纹深度的增加会导致稀疏测点下的转角模态和曲率模态幅值均产生微小变化，但在损伤位置并未存在奇异性峰值，故采用损伤前后模态差作为判断指标，以提高识别灵敏度，稀疏测点与密集测点的模态差如图5所示，工况1/0.25表示工况1/测点间距为0.25 m。无噪声情况下，稀疏测点的模态差在损伤位置处均未出现明显的奇异性峰值，而密集测点的模态差可以准确判断出损伤位置；峰值的相对大小与损伤程度呈正比关系，由此可区分损伤程度，且曲率差指标更为准确。

图5 基于第一阶振型的转角及曲率差Fig.5 Rotation and curvature difference based on the first mode

图6为稀疏、密集测点对简支梁多处损伤识别的对比结果。对于图6(a)中的转角模态，稀疏测点无法实现损伤识别；当间距为0.01 m时，密集测点出现奇异性峰值，但对于靠近跨中处的损伤，转角模态幅值出现减小的情况，并无奇异性峰值，对损伤不够敏感。对于图6(b)中的曲率模态，当损伤发生在稀疏测点位置时，不能识别损伤，而当损伤发生在测点之间时，会导致损伤的错误定位以及损伤程度的辨识错误；当间距为0.05 m时，根据密集测点曲率模态的峰值位置和大小，能够准确区分不同位置和不同程度的损伤。

图6 基于第一阶振型的两处损伤识别Fig.6 Two damage identifications based on the first mode

4 实验验证

选取长×宽×高为1.5 m×0.05 m×0.01 m的简支木梁为研究对象，在梁侧面从左至右分别粘贴旋子、苏式、和玺3种常见的古建筑彩绘图，共同作为DIC拍摄区域，3种彩绘图的长度均为0.25 m。采用1 920×1 080像素的SONY HDR-CX670摄像机进行拍摄，帧率取为50 Hz。先以0.05 m等间距布置的测点作为密集测点，同时在梁上以0.25 m等间距布置5个加速度传感器，设置其采样频率为64 Hz。在梁跨中施加竖向初位移，以木梁自由振动方式激振。实验布置如图7所示。

图7 实验布置图

以数值模拟算例定义的损伤方式在木梁底部设置1 mm深度的横向通长裂缝，损伤工况如表2所示。

表2 试验损伤工况Table 2 Damage cases of test

摄像机光轴垂直于木梁表面进行拍摄时，表示物理位移与图像位移之间关系的尺度因子λ=dknown/Iknown[28]，其中，dknown为物体表面已知长度，mm；Iknown为在图像上对应的像素长度，pixel。实验木梁的厚度为10 mm，在图像上对应像素长度为16 pixel，故尺度因子λ=0.625 mm/pixel，像素间距为0.625 mm，因此，DIC方法能够精确实现等间距密集测量。

在木梁跨中布置激光位移传感器(LDS)以验证DIC方法测量位移的准确性。LDS的型号为LR-TB2000，采样频率为64 Hz。根据计算的尺度因子，将图像位移转换为物理位移，跨中处自由衰减的位移时程曲线如图8(a)所示。由图8(a)可知，激光位移传感器与DIC方法获取的跨中位置处位移时程曲线基本一致，从而说明DIC测量结构位移具有较好的准确性。

采用DIC方法识别出的位移时程曲线如图8(b)所示，图例中数值代表距左端支座的距离。从图8(b)可知，0.25与1.25 m处、0.5与1 m处的测点振动幅值基本吻合，即对称位置处的位移时程曲线一致，与理论相符，表明不同类型的古建筑彩绘图对DIC的识别精度影响不大。

图8 位移时程曲线Fig.8 Displacement time-history

通过识别梁上密集测点的位移时程曲线，进行FFT变换,获取各像素点的幅频曲线与相频曲线。根据相位角正负与幅值大小进行归一化处理，获得梁的一阶工作模态，无损状态下的工作模态如图9所示。由该图可知，通过DIC密集测点获得的工作模态曲线更光滑，包含更多的局部信息。

图9 一阶工作模态Fig.9 First order operation

根据数值模拟实验，考虑转角模态与曲率模态对损伤的相对敏感程度，以及实验中环境噪声与光照条件的影响，采用较为敏感的曲率差损伤指标判断损伤识别效果，基于第一阶振型不同损伤程度的单处损伤识别结果如图10所示。图10(a)中，根据不同工况下、不同测点位置处的曲率差峰值并不能确定损伤位置在0.5 m处，相反，在其他测点出现了错误的损伤定位。原因在于传感器与木梁的接触面较大，安装传感器时不能准确布置在目标采集点位置上，导致利用传感器获取的位移模态进行中心差分时，测点间距的较大误差会造成曲率模态峰值发生明显改变，从而导致损伤位置的定位错误，并且根据峰值的相对大小也难以判断损伤程度，因此，利用传统加速度传感器难以识别木梁的损伤。

采用DIC方法识别结构模态时，虽然图像像素的最小间距为0.625 mm，但考虑到计算效率，先选取0.05 m作为密集测点，识别结果如图10(b)所示。从图10可知，通过DIC方法计算的曲率差在木梁距左端0.5 m处存在明显的奇异性峰值，可实现损伤的准确定位，且不同损伤情况下的峰值也存在差别，其大小与损伤程度呈正比关系。因此，可通过奇异性峰值大小判断木梁损伤的相对大小。

图10 基于第一阶振型的单处损伤曲率差Fig.10 Curvature difference of single damage based on the first mode

在木梁存在两处损伤的情况下，稀疏、密集测点的识别结果如图11所示。图11(a)中存在3处明显峰值，其中0.75 m处的峰值最大，与实际损伤位置存在较大的误差，原因在于根据曲率模态计算式(15)，损伤造成的位移模态突变会导致其相邻测点的曲率模态值发生改变，因此，0.5 m处的损伤导致位移模态在该处发生突变，0.75 m处测点的曲率模态由于该突变产生了峰值，同时0.8 m处的损伤也会对0.75 m处测点造成相同的影响，因此，稀疏测点易受多处损伤的影响而导致错误的损伤定位，并且会造成损伤位置处的奇异性峰值叠加，从而导致损伤程度的误判。

由图11(b)可知，密集测点下，在损伤位置0.5、0.8 m处存在明显奇异性峰值，可以判断在这两处位置存在损伤，并且对应位置处的峰值大小与预设的裂纹深度成正比。因此，当木梁存在多处损伤时，该方法不但可实现小尺度损伤的准确识别，还能判断不同损伤的相对大小。将测点的间距进一步减小，能够实现两处损伤引起的位移模态突变不会影响同一测点的曲率模态值，从而保证密集测点下的损伤准确定位、损伤位置处奇异性峰值保持稳定。

图11 基于第一阶振型的两处损伤曲率差Fig.11 Curvature difference of two damages based on the first mode

综上，实验表明，使用传统加速度传感器进行损伤识别受限于测点间距过大、布设位置偏离等因素，造成损伤的错误定位；而DIC方法的密集测点能够实现多处损伤的定位以及损伤相对程度的辨别。

5 结论

针对接触式传感器进行古建筑木结构彩绘梁动态损伤识别易造成木结构损伤的问题，利用木结构梁表面自身彩绘图案，选取不同模态参数作为损伤指标，提出基于数字图像相关(DIC)原理的非接触式损伤识别方法，数值模拟和彩绘简支梁实验表明：

1)以古建筑彩绘图案为对象，DIC方法能实现彩绘梁的无损全场动态位移精准测量，通过密集测点获取损伤信息更多的工作模态，有效克服了传统接触式测量造成木梁表面损坏、测点不足的缺陷。

2)以规则排列的图像像素为测点，能够实现绝对等间距测量，避免接触式传感器测点布置偏离导致的信息采集错误。

3)DIC方法测点间距可达毫米级，有效避免了两处损伤对同一测点的累积影响，能够实现木梁小尺度直裂缝的准确定位和损伤相对大小的辨识。

4)受限于古建筑内部场地大小，单个相机可能存在无法拍摄大跨度木梁的缺点，后期需进一步研究使用多个相机进行图像拼接，从而实现大型古建筑结构整体的模态测试。