熊慧芳 李姍柵 肖 潭

1(廣東石油化工學院計算機與電子信息學院 廣東 茂名 525000)2(中山大學工學院 廣東 廣州 510006)3(廣東石油化工學院力學教學與實驗中心 廣東 茂名 525000)

?

一種用于松質(zhì)骨變形測量的數(shù)字圖像相關方法

熊慧芳1李姍柵2肖 潭3*

1(廣東石油化工學院計算機與電子信息學院 廣東 茂名 525000)2(中山大學工學院 廣東 廣州 510006)3(廣東石油化工學院力學教學與實驗中心 廣東 茂名 525000)

針對松質(zhì)骨這一類多孔材料，直接利用被測物體自身的微結(jié)構(gòu)或紋理代替人工散斑標記，然后使用Matlab編寫變形測量計算程序，嘗試將數(shù)字圖像相關方法DIC(digital image correlation)應用于松質(zhì)骨變形測量。生成模擬散斑圖像對該算法程序的可行性進行模擬驗證，再利用豬股骨制作壓縮試件對整個變形測量方案進行實驗驗證。模擬驗證和實驗驗證結(jié)果表明：利用松質(zhì)骨自身微結(jié)構(gòu)紋理作為標記并編寫算法程序進行DIC變形測量的方案是可行的，該方案也可用于和松質(zhì)骨具有類似微結(jié)構(gòu)的其他多孔材料的變形測量中。

變形測量 松質(zhì)骨 數(shù)字圖像相關方法 多孔材料 散斑標記

0 引 言

近十幾年來，中國人口老齡化問題日漸凸顯，老年人骨折、骨質(zhì)疏松等問題得到了更廣泛的關注。骨骼是組成人體的重要部分，按結(jié)構(gòu)可分為密質(zhì)骨和松質(zhì)骨兩類。骨的受載變形是骨力學研究的重要內(nèi)容，變形測量方法有接觸式和非接觸式兩種[1]。接觸式是將應變片等測量裝置直接固定在試件上進行測量，要求試件具有足夠的剛度，只能測量某點的變形，但測量結(jié)果可靠、精度高。非接觸式則是以光、聲、電磁波等為媒介，間接對試件進行測量。該方式對試件的限制較少，為全場測量，但精度相對較低。其中，數(shù)字圖像相關方法DIC具有測量設備易搭建和不易受周圍環(huán)境干擾等優(yōu)點，應用日益廣泛。

早期,Zhang和Arola[2]曾將DIC應用于牛蹄角及牛動脈血管的應變測量。近期,Begonia等[3]應用DIC測量老鼠前臂在軸向壓縮載荷作用下的應變；Srinivasan等[4]使用新開發(fā)的DIC裝置測量骨水泥與骨小梁之間的微動與應變，研究微動所導致的組織液增多對骨小梁的吸收效應；Kaviani等[5]在應力松弛試驗中利用共聚焦顯微鏡和DIC測量軟骨生長板應變分布，以評估靜、動載荷以及動載荷頻率和幅度對軟骨組織形態(tài)的影響?；贒IC的變形測量一般需要在試件表面人工噴漆制備散斑標記，只能測量表面變形。由于松質(zhì)骨具有多孔結(jié)構(gòu)，人們嘗試將DIC拓展到三維來分析其力學性能。采用高分辨率CT進行掃描，Gillard[6]等建立了松質(zhì)骨試件骨小梁三維結(jié)構(gòu)模型，然后通過結(jié)構(gòu)子集的匹配來計算試件變形，并分析CT圖像噪聲以及結(jié)構(gòu)子集的劃分對測量結(jié)果的影響。這種以試件內(nèi)部結(jié)構(gòu)為散斑標記計算試件變形的方法稱為數(shù)字體圖像相關法DVC(digital volume correlation)[7]，其算法的選擇與計算量的控制很重要。除松質(zhì)骨外，DIC還可應用于其他多孔材料的變形測量，例如金屬纖維燒結(jié)板的面內(nèi)和橫向剪切變形測量[8]、多孔陶瓷[9]和泡沫金屬[10]材料單軸拉伸時試件全場位移測量。

目前，像松質(zhì)骨這一類多孔材料在受載情況下的內(nèi)部變形，只能選取孔洞邊緣特征作為位移標記進行測量。有別于常規(guī)的漆類散斑標記，本文直接利用松質(zhì)骨自身的微結(jié)構(gòu)特征作為標記，進行松質(zhì)骨二維DIC變形測量方案研究，分析形函數(shù)、相關準則、圖像子集、對比度等參數(shù)對松質(zhì)骨變形測量的精度和計算時間的影響并優(yōu)化參數(shù)，然后編寫針對松質(zhì)骨的DIC算法程序，進行模擬驗證和實際試驗驗證。

1 DIC方法及改進

為了檢測松質(zhì)骨試件受外載荷作用后的變形情況，需要自行編寫以DIC原理為基礎的測量分析計算程序并在MATLAB環(huán)境下實現(xiàn)。

1.1 算法原理

圖1 DIC算法原理示意圖

開發(fā)DIC算法程序之前需要確定以下幾個問題：

1) 形函數(shù)

(1)

(2M+1)×(2M+1)為參考子集的大小，α(xi,yi)、β(xi,yi)即為形函數(shù)。當且僅當測試試件為剛體時，方可使用零階形函數(shù)，表達式如下：

(2)

如果試驗過程中試件還發(fā)生了轉(zhuǎn)動、伸縮變形，則在計算過程中需選用一階或更高階的形函數(shù)[13]。一階形函數(shù)的一般表示式如下：

(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)之間的距離。

松質(zhì)骨實際受壓時的變形不僅包含剛體位移，同時還包含伸縮等變形，需選取一階或更高階的形函數(shù)。相對高階的形函數(shù)會相應地延長計算時間，綜合考慮算法準確度和計算效率，本系統(tǒng)選用一階形函數(shù)。

2) 相關準則

在變形圖像中對參考子區(qū)進行特征匹配時，需要事先選定一個相關準則作為判定條件。相關準則主要有兩類，第一類為互相關準則，包括：

(1) 互相關(CC)

(4)

(2) 標準化互相關(NCC)

(5)

(6)

第二類為平方和準則，包括：

(1) 平方和(SSD)

(7)

(2) 標準化平方和(NSSD)：

(8)

(9)

(10)

考慮到圖像獲取過程中光強的穩(wěn)定性、相機曝光與畸變、周圍環(huán)境噪聲的影響，需考慮各準則的穩(wěn)定性及準確性。CC/SSD準則對光強線性放大及補償均比較敏感，NCC/NSSD準則僅對光強補償敏感，而且這四種準則均易受環(huán)境噪聲影響。然而ZNCC/ZNSSD準則對光強線性放大及補償均不敏感，且具有良好的抗噪性[14]。Pan等的研究表明引入梯度變換模型后，相較于上述四種相關準則，ZNCC/ZNSSD具有更強的抗噪性、穩(wěn)定性及更高的精確度[15]。由于本DIC測量裝置是通過手動方式來調(diào)整光源，會導致圖像光強分布的不穩(wěn)定性及不平均性，所以優(yōu)先考慮ZNCC或ZNSSD準則。采用模擬散斑圖進行對比計算，上述兩種準則精確度相似，但ZNCC準則計算耗時較短，故最終選用ZNCC準則作為特征匹配時的判定標準。

3) 搜索方法

利用DIC進行位移測量，具體過程可分為兩部分：整像素搜索和亞像素搜索。目前比較常用的亞像素搜索方法有三種：曲面擬合法、牛頓迭代法、梯度法。梯度法是一種基于光流法的新型計算方法，它將空間梯度與迭代相結(jié)合進行計算。相較于曲面擬合法，梯度法考慮了形函數(shù)，具有更高的計算準確度；相較于牛頓迭代法，梯度法只需計算一階梯度，大大減少了計算量。2009年Pan[16]提出一種基于迭代的梯度法ILS(iterative least squares algorithm)，經(jīng)證明與改進的牛頓迭代法是等價的。綜合分析目前主流的亞像素搜索方法的優(yōu)缺點，考慮松質(zhì)骨試件微結(jié)構(gòu)特征的特點，本文算法程序選用了梯度法進行亞像素級別的位移搜索。

1.2 紋理代替散斑的可行性

為了驗證利用松質(zhì)骨自身多孔微結(jié)構(gòu)特征替代人工散斑的可行性，我們制備了一塊松質(zhì)骨試件然后進行拍攝，如圖2(a)所示；圖2(b)是圖2(a)某局部放大灰度圖。對圖2(b)進行灰度轉(zhuǎn)換，灰度直方圖如圖2(c)所示，可看出灰度分布比較均勻，具有良好的對比度，圖片可以直接用于DIC變形分析。

圖2 試件的原始拍攝圖像、局部灰度圖像及灰度直方圖

為了進一步提高計算效率，還可以在整像素搜索中進行一些改進。傳統(tǒng)的整像素搜索以參考子區(qū)為模板，在變形圖像中逐點計算相關系數(shù)，系數(shù)最高/最低的點即為變形后所在整像素位置。逐點搜索法結(jié)果準確但比較費時，考慮算法整體的計算效率，我們在計算過程中添加了一個人機交互功能，可根據(jù)試驗過程中不同的變形情況人為給定初始位移，從而省略冗余的搜索過程，縮短計算時間。

2 模擬圖像驗證

實際拍攝所獲得的圖片含有光線、鏡頭等各種外界因素引起的誤差，而通過數(shù)值模擬生成的散斑圖片則能避免大多數(shù)外界因素的影響。所以通過人工模擬生成移動前后的散斑圖，用以驗證作者所開發(fā)的DIC變形分析程序計算結(jié)果的準確性。散斑圖的生成原理為利用高斯函數(shù)模擬生成若干隨機分布的散斑，散斑的光強由對應點的函數(shù)值大小決定[17]：

(11)

(12)

其中Ir、Id為移動前后的圖像光強，S為散斑顆粒數(shù)，a為散斑顆粒的平均大小，rk為每一個隨機分布的散斑的位置，I0為高斯光斑的中心光強，U(r)為給定的位移場。

第一次試驗所生成的模擬散斑圖參數(shù)如下：圖像大小480×480 pixel，散斑數(shù)量3 000，散斑平均大小3 pixel，給定位移為豎直方向(Y方向)向下移動1個像素。移動前后的模擬散斑圖像如圖3(a)、(b)所示。

圖3 模擬散斑圖

利用作者所編寫的程序?qū)ι鲜鰞煞M圖像進行位移分析，結(jié)果如圖4(a)所示。然后在水平方向(X方向)分別給定0.5、0.1、0.05 pixel的位移，生成相應的散斑圖，再利用本程序?qū)ζ溥M行位移分析，所得結(jié)果如圖4(b)、(c)、(d)所示，結(jié)果的均值誤差及標準差列在表3中。

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

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

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

(d) X方向給定0.05 pixel位移圖4 程序分析計算結(jié)果

表3 測量結(jié)果均值誤差、標準差

圖4和表3的結(jié)果表明作者所編寫的DIC變形分析程序在1 pixel級別的準確性和穩(wěn)定性最好，在亞像素級別的準確度一般，且測量數(shù)據(jù)的波動也較大。按普通數(shù)碼相機的分辨率和標準試件尺寸計算，散斑圖片中1 pixel的位移所對應的試件真實位移已達微米級別，能夠滿足松質(zhì)骨壓縮試驗變形測量的需要。

3 松質(zhì)骨壓縮試驗

3.1 試件準備

考慮到人骨樣本的獲取有難度，本試驗選用易獲得且性質(zhì)比較接近的豬股骨制作壓縮試件，試件準備過程如下：

(1) 預處理

將購入的豬股骨進行初步處理，去除表面殘余軟組織，用純水反復沖洗至水清；

(2) 切割

預處理后通過人工切割，制成壓縮試件所需方塊，尺寸約為1×1×1 cm3；

(3) 脫脂

使用氯仿/甲醇溶液對試件進行脫脂，然后試件自然風干12 h；

(4) 保存

風干后的試件置于-16 ℃的冰箱中封存，試驗前取出靜置，至試件溫度升至室溫后開始試驗。

3.2 試驗過程

采用型號LLOYD LR10K Plus萬能材料試驗機進行松質(zhì)骨壓縮試驗，選擇10 kN傳感器；圖像拍攝采用佳能EOS 5D Mark II相機，鏡頭型號佳能EF 100 mm f/2.8 USM。

試驗時材料試驗機參數(shù)設定加載速度3 mm/min，加載時間1 min。將試驗樣品放置于特制壓盤上，啟動萬能材料試驗機，同時使用相機拍攝樣品整個加載過程，視頻拍攝速度30幀/秒。

3.3 試驗結(jié)果分析

從試驗視頻的第5 s開始，每隔5 s抽取一幅圖像，整個加載過程共抽取七幅圖像，第一幅作為試件變形前的參考圖像，余下六幅作為試件逐級變形后的圖像。使用作者所開發(fā)的程序進行分析計算，試件逐級變形的位移場顯示在圖5中，相對應的應變場顯示在圖6中，各圖像中1 pixel所對應的實際長度約為18 μm。

在圖5和圖6中，圖像I表示水平方向(X方向)測量結(jié)果，圖像II表示豎直方向(Y方向)的測量結(jié)果。在整個試驗過程中，試件大部分區(qū)域豎直方向上的位移和應變值都比水平方向上的值要大。因為試件受載方向為豎直方向，通常材料的泊松比小于1，所以測量值總體上表現(xiàn)為豎直方向上較大。

從圖像(a)到圖像(f)，隨著壓縮試驗的進行，圖5中的位移和圖6中的應變在水平和豎直方向呈增大的趨勢，同時位移和應變在水平和豎直方向都呈不均勻分布，試件中部及上部的位移和應變數(shù)值較大，符合實際試驗時所觀察到的試件變形情況。在試驗后期，當松質(zhì)骨試件開始出現(xiàn)結(jié)構(gòu)破裂時，我們觀測到這種破裂大多數(shù)出現(xiàn)在試件上半部分。這種現(xiàn)象應該與松質(zhì)骨試件本身力學性能的不均勻分布有關。松質(zhì)骨的外圍是密質(zhì)骨，在松質(zhì)骨內(nèi)靠近密質(zhì)骨處和距離密質(zhì)骨有一定距離處，骨小梁的粗細、疏密分布以及力學性能是有差異的。

在圖5II中，豎直位移向下為正。其中圖5II(f)矩形框內(nèi)的區(qū)域取正值，是整個試件主要發(fā)生向下位移的區(qū)域；如試件上端圓框所示，分散存在著若干局部位移為負的區(qū)域。由于松質(zhì)骨并不是一個符合連續(xù)化假設的材料，而是一個形態(tài)和力學性質(zhì)分布都不均勻的結(jié)構(gòu)。當不均勻結(jié)構(gòu)受壓變形時，局部區(qū)域可能發(fā)生和其他區(qū)域位移不同甚至相反的不穩(wěn)定變形，即力學上的局部失穩(wěn)現(xiàn)象。進一步對圖像進行放大觀察，相關區(qū)域的骨小梁結(jié)構(gòu)存確實在失穩(wěn)或斷裂現(xiàn)象，該局部的位移方向和周邊結(jié)構(gòu)的位移方向不一致。

Ⅰ為X方向

Ⅱ為Y方向(a)-(f)分別為t=10、15、20、25、30、35 s時的位移測量結(jié)果(原圖為彩色)圖5 試件位移測量結(jié)果

Ⅰ為X方向

Ⅱ為Y方向(a)-(f)分別為t= 10、15、20、25、30、35 s時的應變測量結(jié)果(原圖為彩色)圖6 試件應變測量結(jié)果

4 結(jié) 語

盡管DIC方法的提出才短短三十多年，但是因其特有的非接觸性、全場性以及相對簡單的測量要求，該方法已在多種材料、尤其是生物材料的位移應變測量中有廣泛的應用。骨骼的機械性能一直是生物醫(yī)學工程研究者們關注的重點。作者設計并搭建了一套針對松質(zhì)骨等多孔材料微小變形的測量系統(tǒng)，以DIC算法為基礎開發(fā)了相關的變形分析程序并進行模擬和試驗驗證，結(jié)果表明：可借助松質(zhì)骨自身微結(jié)構(gòu)特征替代人工制備散斑進行位移測量，位移測量精度可達到微米級別，應變精度可達到με級別。本測量系統(tǒng)不僅能應用于松質(zhì)骨等多孔材料的變形測量，通過對算法進行相應的改變與調(diào)試，還可拓展至其它具有天然紋理材料如木材、石材的變形測量中。利用材料自身所具有的微結(jié)構(gòu)或紋理特征作為標記進行變形測量是實現(xiàn)材料內(nèi)部三維變形測量的技術途徑之一。作為初步探索，本文實現(xiàn)了對松質(zhì)骨試件的二維變形測量，并為相關算法拓展至三維預留了空間。

[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