一種用于松質骨變形測量的數字圖像相關方法

熊慧芳 李姍柵 肖 潭

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