一种用于松质骨变形测量的数字图像相关方法

熊慧芳 李姗栅 肖 潭

1(广东石油化工学院计算机与电子信息学院 广东 茂名 525000)2(中山大学工学院 广东 广州 510006)3(广东石油化工学院力学教学与实验中心 广东 茂名 525000)



一种用于松质骨变形测量的数字图像相关方法

熊慧芳1李姗栅2肖 潭3*

1(广东石油化工学院计算机与电子信息学院 广东 茂名 525000)2(中山大学工学院 广东 广州 510006)3(广东石油化工学院力学教学与实验中心 广东 茂名 525000)

针对松质骨这一类多孔材料，直接利用被测物体自身的微结构或纹理代替人工散斑标记，然后使用Matlab编写变形测量计算程序，尝试将数字图像相关方法DIC(digital image correlation)应用于松质骨变形测量。生成模拟散斑图像对该算法程序的可行性进行模拟验证，再利用猪股骨制作压缩试件对整个变形测量方案进行实验验证。模拟验证和实验验证结果表明：利用松质骨自身微结构纹理作为标记并编写算法程序进行DIC变形测量的方案是可行的，该方案也可用于和松质骨具有类似微结构的其他多孔材料的变形测量中。

变形测量 松质骨 数字图像相关方法 多孔材料 散斑标记

0 引 言

近十几年来，中国人口老龄化问题日渐凸显，老年人骨折、骨质疏松等问题得到了更广泛的关注。骨骼是组成人体的重要部分，按结构可分为密质骨和松质骨两类。骨的受载变形是骨力学研究的重要内容，变形测量方法有接触式和非接触式两种[1]。接触式是将应变片等测量装置直接固定在试件上进行测量，要求试件具有足够的刚度，只能测量某点的变形，但测量结果可靠、精度高。非接触式则是以光、声、电磁波等为媒介，间接对试件进行测量。该方式对试件的限制较少，为全场测量，但精度相对较低。其中，数字图像相关方法DIC具有测量设备易搭建和不易受周围环境干扰等优点，应用日益广泛。

早期,Zhang和Arola[2]曾将DIC应用于牛蹄角及牛动脉血管的应变测量。近期,Begonia等[3]应用DIC测量老鼠前臂在轴向压缩载荷作用下的应变；Srinivasan等[4]使用新开发的DIC装置测量骨水泥与骨小梁之间的微动与应变，研究微动所导致的组织液增多对骨小梁的吸收效应；Kaviani等[5]在应力松弛试验中利用共聚焦显微镜和DIC测量软骨生长板应变分布，以评估静、动载荷以及动载荷频率和幅度对软骨组织形态的影响。基于DIC的变形测量一般需要在试件表面人工喷漆制备散斑标记，只能测量表面变形。由于松质骨具有多孔结构，人们尝试将DIC拓展到三维来分析其力学性能。采用高分辨率CT进行扫描，Gillard[6]等建立了松质骨试件骨小梁三维结构模型，然后通过结构子集的匹配来计算试件变形，并分析CT图像噪声以及结构子集的划分对测量结果的影响。这种以试件内部结构为散斑标记计算试件变形的方法称为数字体图像相关法DVC(digital volume correlation)[7]，其算法的选择与计算量的控制很重要。除松质骨外，DIC还可应用于其他多孔材料的变形测量，例如金属纤维烧结板的面内和横向剪切变形测量[8]、多孔陶瓷[9]和泡沫金属[10]材料单轴拉伸时试件全场位移测量。

目前，像松质骨这一类多孔材料在受载情况下的内部变形，只能选取孔洞边缘特征作为位移标记进行测量。有别于常规的漆类散斑标记，本文直接利用松质骨自身的微结构特征作为标记，进行松质骨二维DIC变形测量方案研究，分析形函数、相关准则、图像子集、对比度等参数对松质骨变形测量的精度和计算时间的影响并优化参数，然后编写针对松质骨的DIC算法程序，进行模拟验证和实际试验验证。

1 DIC方法及改进

为了检测松质骨试件受外载荷作用后的变形情况，需要自行编写以DIC原理为基础的测量分析计算程序并在MATLAB环境下实现。

1.1 算法原理

图1 DIC算法原理示意图

开发DIC算法程序之前需要确定以下几个问题：

1) 形函数

(1)

(2M+1)×(2M+1)为参考子集的大小，α(xi,yi)、β(xi,yi)即为形函数。当且仅当测试试件为刚体时，方可使用零阶形函数，表达式如下：

(2)

如果试验过程中试件还发生了转动、伸缩变形，则在计算过程中需选用一阶或更高阶的形函数[13]。一阶形函数的一般表示式如下：

(3)

其中，u、v分别为参考子集中心点P(x0,y0)在x、y方向上的整像素位移，Δu、Δv分别为中心点P(x0,y0)在亚像素级别的位移，ux、uy、vx、vy分别为u、v在x、y方向上的一阶梯度，Δxi,Δyi分别为Q(xi,yi)点在x、y方向上和中心点P(x0,y0)之间的距离。

松质骨实际受压时的变形不仅包含刚体位移，同时还包含伸缩等变形，需选取一阶或更高阶的形函数。相对高阶的形函数会相应地延长计算时间，综合考虑算法准确度和计算效率，本系统选用一阶形函数。

2) 相关准则

在变形图像中对参考子区进行特征匹配时，需要事先选定一个相关准则作为判定条件。相关准则主要有两类，第一类为互相关准则，包括：

(1) 互相关(CC)

(4)

(2) 标准化互相关(NCC)

(5)

(6)

第二类为平方和准则，包括：

(1) 平方和(SSD)

(7)

(2) 标准化平方和(NSSD)：

(8)

(9)

(10)

考虑到图像获取过程中光强的稳定性、相机曝光与畸变、周围环境噪声的影响，需考虑各准则的稳定性及准确性。CC/SSD准则对光强线性放大及补偿均比较敏感，NCC/NSSD准则仅对光强补偿敏感，而且这四种准则均易受环境噪声影响。然而ZNCC/ZNSSD准则对光强线性放大及补偿均不敏感，且具有良好的抗噪性[14]。Pan等的研究表明引入梯度变换模型后，相较于上述四种相关准则，ZNCC/ZNSSD具有更强的抗噪性、稳定性及更高的精确度[15]。由于本DIC测量装置是通过手动方式来调整光源，会导致图像光强分布的不稳定性及不平均性，所以优先考虑ZNCC或ZNSSD准则。采用模拟散斑图进行对比计算，上述两种准则精确度相似，但ZNCC准则计算耗时较短，故最终选用ZNCC准则作为特征匹配时的判定标准。

3) 搜索方法

利用DIC进行位移测量，具体过程可分为两部分：整像素搜索和亚像素搜索。目前比较常用的亚像素搜索方法有三种：曲面拟合法、牛顿迭代法、梯度法。梯度法是一种基于光流法的新型计算方法，它将空间梯度与迭代相结合进行计算。相较于曲面拟合法，梯度法考虑了形函数，具有更高的计算准确度；相较于牛顿迭代法，梯度法只需计算一阶梯度，大大减少了计算量。2009年Pan[16]提出一种基于迭代的梯度法ILS(iterative least squares algorithm)，经证明与改进的牛顿迭代法是等价的。综合分析目前主流的亚像素搜索方法的优缺点，考虑松质骨试件微结构特征的特点，本文算法程序选用了梯度法进行亚像素级别的位移搜索。

1.2 纹理代替散斑的可行性

为了验证利用松质骨自身多孔微结构特征替代人工散斑的可行性，我们制备了一块松质骨试件然后进行拍摄，如图2(a)所示；图2(b)是图2(a)某局部放大灰度图。对图2(b)进行灰度转换，灰度直方图如图2(c)所示，可看出灰度分布比较均匀，具有良好的对比度，图片可以直接用于DIC变形分析。

图2 试件的原始拍摄图像、局部灰度图像及灰度直方图

为了进一步提高计算效率，还可以在整像素搜索中进行一些改进。传统的整像素搜索以参考子区为模板，在变形图像中逐点计算相关系数，系数最高/最低的点即为变形后所在整像素位置。逐点搜索法结果准确但比较费时，考虑算法整体的计算效率，我们在计算过程中添加了一个人机交互功能，可根据试验过程中不同的变形情况人为给定初始位移，从而省略冗余的搜索过程，缩短计算时间。

2 模拟图像验证

实际拍摄所获得的图片含有光线、镜头等各种外界因素引起的误差，而通过数值模拟生成的散斑图片则能避免大多数外界因素的影响。所以通过人工模拟生成移动前后的散斑图，用以验证作者所开发的DIC变形分析程序计算结果的准确性。散斑图的生成原理为利用高斯函数模拟生成若干随机分布的散斑，散斑的光强由对应点的函数值大小决定[17]：

(11)

(12)

其中Ir、Id为移动前后的图像光强，S为散斑颗粒数，a为散斑颗粒的平均大小，rk为每一个随机分布的散斑的位置，I0为高斯光斑的中心光强，U(r)为给定的位移场。

第一次试验所生成的模拟散斑图参数如下：图像大小480×480 pixel，散斑数量3 000，散斑平均大小3 pixel，给定位移为竖直方向(Y方向)向下移动1个像素。移动前后的模拟散斑图像如图3(a)、(b)所示。

图3 模拟散斑图

利用作者所编写的程序对上述两幅模拟图像进行位移分析，结果如图4(a)所示。然后在水平方向(X方向)分别给定0.5、0.1、0.05 pixel的位移，生成相应的散斑图，再利用本程序对其进行位移分析，所得结果如图4(b)、(c)、(d)所示，结果的均值误差及标准差列在表3中。

(a) Y方向上给定1 pixel位移

(b) X方向给定0.5 pixel位移

(c) X方向给定0.1 pixel位移

(d) X方向给定0.05 pixel位移图4 程序分析计算结果

表3 测量结果均值误差、标准差

图4和表3的结果表明作者所编写的DIC变形分析程序在1 pixel级别的准确性和稳定性最好，在亚像素级别的准确度一般，且测量数据的波动也较大。按普通数码相机的分辨率和标准试件尺寸计算，散斑图片中1 pixel的位移所对应的试件真实位移已达微米级别，能够满足松质骨压缩试验变形测量的需要。

3 松质骨压缩试验

3.1 试件准备

考虑到人骨样本的获取有难度，本试验选用易获得且性质比较接近的猪股骨制作压缩试件，试件准备过程如下：

(1) 预处理

将购入的猪股骨进行初步处理，去除表面残余软组织，用纯水反复冲洗至水清；

(2) 切割

预处理后通过人工切割，制成压缩试件所需方块，尺寸约为1×1×1 cm3；

(3) 脱脂

使用氯仿/甲醇溶液对试件进行脱脂，然后试件自然风干12 h；

(4) 保存

风干后的试件置于-16 ℃的冰箱中封存，试验前取出静置，至试件温度升至室温后开始试验。

3.2 试验过程

采用型号LLOYD LR10K Plus万能材料试验机进行松质骨压缩试验，选择10 kN传感器；图像拍摄采用佳能EOS 5D Mark II相机，镜头型号佳能EF 100 mm f/2.8 USM。

试验时材料试验机参数设定加载速度3 mm/min，加载时间1 min。将试验样品放置于特制压盘上，启动万能材料试验机，同时使用相机拍摄样品整个加载过程，视频拍摄速度30帧/秒。

3.3 试验结果分析

从试验视频的第5 s开始，每隔5 s抽取一幅图像，整个加载过程共抽取七幅图像，第一幅作为试件变形前的参考图像，余下六幅作为试件逐级变形后的图像。使用作者所开发的程序进行分析计算，试件逐级变形的位移场显示在图5中，相对应的应变场显示在图6中，各图像中1 pixel所对应的实际长度约为18 μm。

在图5和图6中，图像I表示水平方向(X方向)测量结果，图像II表示竖直方向(Y方向)的测量结果。在整个试验过程中，试件大部分区域竖直方向上的位移和应变值都比水平方向上的值要大。因为试件受载方向为竖直方向，通常材料的泊松比小于1，所以测量值总体上表现为竖直方向上较大。

从图像(a)到图像(f)，随着压缩试验的进行，图5中的位移和图6中的应变在水平和竖直方向呈增大的趋势，同时位移和应变在水平和竖直方向都呈不均匀分布，试件中部及上部的位移和应变数值较大，符合实际试验时所观察到的试件变形情况。在试验后期，当松质骨试件开始出现结构破裂时，我们观测到这种破裂大多数出现在试件上半部分。这种现象应该与松质骨试件本身力学性能的不均匀分布有关。松质骨的外围是密质骨，在松质骨内靠近密质骨处和距离密质骨有一定距离处，骨小梁的粗细、疏密分布以及力学性能是有差异的。

在图5II中，竖直位移向下为正。其中图5II(f)矩形框内的区域取正值，是整个试件主要发生向下位移的区域；如试件上端圆框所示，分散存在着若干局部位移为负的区域。由于松质骨并不是一个符合连续化假设的材料，而是一个形态和力学性质分布都不均匀的结构。当不均匀结构受压变形时，局部区域可能发生和其他区域位移不同甚至相反的不稳定变形，即力学上的局部失稳现象。进一步对图像进行放大观察，相关区域的骨小梁结构存确实在失稳或断裂现象，该局部的位移方向和周边结构的位移方向不一致。

Ⅰ为X方向

Ⅱ为Y方向(a)-(f)分别为t=10、15、20、25、30、35 s时的位移测量结果(原图为彩色)图5 试件位移测量结果

Ⅰ为X方向

Ⅱ为Y方向(a)-(f)分别为t= 10、15、20、25、30、35 s时的应变测量结果(原图为彩色)图6 试件应变测量结果

4 结 语

尽管DIC方法的提出才短短三十多年，但是因其特有的非接触性、全场性以及相对简单的测量要求，该方法已在多种材料、尤其是生物材料的位移应变测量中有广泛的应用。骨骼的机械性能一直是生物医学工程研究者们关注的重点。作者设计并搭建了一套针对松质骨等多孔材料微小变形的测量系统，以DIC算法为基础开发了相关的变形分析程序并进行模拟和试验验证，结果表明：可借助松质骨自身微结构特征替代人工制备散斑进行位移测量，位移测量精度可达到微米级别，应变精度可达到με级别。本测量系统不仅能应用于松质骨等多孔材料的变形测量，通过对算法进行相应的改变与调试，还可拓展至其它具有天然纹理材料如木材、石材的变形测量中。利用材料自身所具有的微结构或纹理特征作为标记进行变形测量是实现材料内部三维变形测量的技术途径之一。作为初步探索，本文实现了对松质骨试件的二维变形测量，并为相关算法拓展至三维预留了空间。

[1] Grassi L, Isaksson H. Extracting accurate strain measurements in bone mechanics: A critical review of current methods[J]. Journal of the Mechanical Behavior of Biomedical Materials, 2015, 50:43-54.

[2] Zhang D, Arola D D. Applications of digital image correlation to biological tissues.[J]. Journal of Biomedical Optics, 2004, 9(4):691-699.

[3] Begonia M T, Dallas M, Vizcarra B, et al. Non-Contact Strain Measurement in the Mouse Forearm Loading Model Using Digital Image Correlation (DIC)[J]. Bone, 2015, 81:593-601.

[4] Srinivasan P, Miller M A, Verdonschot N, et al. Experimental and computational micromechanics at the tibial cement-trabeculae interface.[J]. Journal of Biomechanics, 2016, 49(9):1641-1648.

[5] Kaviani R, Londono I, Parent S, et al. Growth plate cartilage shows different strain patterns in response to static versus dynamic mechanical modulation[J]. Biomechanics and Modeling in Mechanobiology, 2016, 15(4):1-14.

[6] Gillard F, Boardman R, Mavrogordato M, et al. The application of digital volume correlation (DVC) to study the microstructural behaviour of trabecular bone during compression[J]. Journal of the Mechanical Behavior of Biomedical Materials, 2013, 29C(1):480-499.

[7] Roberts B C, Perilli E, Reynolds K J. Application of the digital volume correlation technique for the measurement of displacement and strain fields in bone: a literature review.[J]. Journal of Biomechanics, 2014, 47(5):923-934.

[8] Zhao T F, Chen C Q. The shear properties and deformation mechanisms of porous metal fiber sintered sheets[J]. Mechanics of Materials, 2014, 70(70):33-40.

[9] Pandey A, Shyam A, Watkins T R, et al. The uniaxial tensile response of porous and microcracked ceramic materials [J]. Journal of the American Ceramic Society, 2014, 97(3):899-906.

[10] Jung A, Wocker M, Chen Z, et al. Microtensile testing of open-cell metal foams — Experimental setup, micromechanical properties[J]. Materials & Design, 2015, 88:1021-1030.

[11] Yamaguchi I. A laser-speckle strain gauge[J]. Journal of Physics E Scientific Instruments, 2000, 14(11):1270.

[12] Peters W H, Ranson W F. Digital Imaging Techniques In Experimental Stress Analysis[M]// Environmental jurisprudence in India. Kluwer Law International, 1999:427-431.

[13] Lu H, Cary P D. Deformation measurements by digital image correlation: Implementation of a second-order displacement gradient[J]. Experimental Mechanics, 2000, 40(4):393-400.

[14] Pan B, Xie H, Wang Z. Equivalence of digital image correlation criteria for pattern matching[J]. Applied Optics, 2010, 49(28):5501-5509.

[15] Pan B, Qian K, Xie H, et al. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review[J]. Measurement Science & Technology, 2009, 20(6):152-154.

[16] Pan B, Asundi A, Xie H, et al. Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements[J]. Optics & Lasers in Engineering, 2009, 47(7-8):865-874.

[17] Zhou P, Goodson K E. Subpixel displacement and deformation gradient measurement using digital image/speckle correlation (DISC)[J]. Optical Engineering, 2001, 40(8):1613-1620.

A DIGITAL IMAGE CORRELATION METHOD FOR DEFORMATION MEASUREMENT OF CANCELLOUS BONE

Xiong Huifang1Li Shanshan2Xiao Tan3*

1(SchoolofComputerandElectronicInformation,GuangdongUniversityofPetrochemicalTechnology,Maoming525000,Guangdong,China)2(SchoolofEngineering,SunYat-senUniversity,Guangzhou510006,Guangdong,China)3(CenterforMechanicalTeachingandTesting,GuangdongUniversityofPetrochemicalTechnology,Maoming525000,Guangdong,China)

Aimed at the porous materials such as cancellous bone, directly using the measured object’s own micro-structure or texture instead of artificial speckle markers, and then use Matlab to prepare the deformation measurement calculation program, and try to apply the digital image correlation(DIC) method to the measurement of cancellous bone deformation. Generating the simulated speckle image to simulate the feasibility of the algorithm, and the whole deformation measurement scheme is validated by using the compression test piece made of porcine femur. The results of simulation and experiment show that it is feasible to use the microstructure texture of the cancellous bone as the mark and use the algorithm program to carry out DIC deformation measurement. The protocol can also be used for deformation measurement of other porous materials having similar microstructures to cancellous bone.

Deformation measurement Cancellous bone Digital image correlation method Porous material Speckle mark

2016-06-20。国家自然科学基金项目(11272360);广东省自然科学基金项目(2014A030313793);茂名市科技计划项目(2016001)。熊慧芳，工程师，主研领域：计算机软件与应用。李姗栅，硕士生。肖潭，教授。

TP39

A

10.3969/j.issn.1000-386x.2017.07.030