鄧超兵，張程瀟，白玉磊，周延周

(廣東工業(yè)大學(xué)自動(dòng)化學(xué)院，廣東廣州510006)

數(shù)字體相關(guān)算法是在DIC技術(shù)基礎(chǔ)上發(fā)展起來(lái)，計(jì)算物體內(nèi)部三維位移場(chǎng)分布的一種算法［1-2］.20世紀(jì)90年代B.K.Bay等［3-4］提出了數(shù)字體相關(guān)方法，并得到了廣泛的關(guān)注，DVC需要的三維數(shù)據(jù)，大體上有兩種來(lái)源:(1)計(jì)算機(jī)層析三維成像;(2)光學(xué)共焦顯微三維成像.潘兵等［5-6］提出基于梯度的測(cè)量物體內(nèi)部變形的DVC算法，牛永強(qiáng)、樊雪松等［7-9］也做了大量有關(guān)方面的研究，J.Huang等［10］利用三維共焦熒光顯微鏡和DVC技術(shù)，測(cè)量加載前后隨機(jī)嵌入微米熒光珠凝膠內(nèi)部三維變形場(chǎng)分布［5］.

目前計(jì)算三維位移場(chǎng)的方法較多［11］.面臨的問(wèn)題就是如何在大量數(shù)據(jù)的前提下，做到最大程度提高位移計(jì)算精度和效率.在本文中，通過(guò)利用圖像一階灰度梯度來(lái)計(jì)算三維位移場(chǎng)，避免采用高階的灰度梯度帶來(lái)的麻煩，計(jì)算效率大幅提升，同時(shí)采用了在亞體素位置三元三次灰度插值法來(lái)提高位移計(jì)算精度的方法，做到精度和效率兼顧.

1 算法

1.1 數(shù)字體圖像相關(guān)(DVC)模型

數(shù)字體圖像是由空間體像素構(gòu)成的三維圖像.當(dāng)圖像內(nèi)部發(fā)生微小形變時(shí)，通過(guò)數(shù)字體圖像相關(guān)法可求出體圖像變形前后的三維位移場(chǎng)分布.其中變形前圖像稱為參考體圖像，變形后圖像稱為目標(biāo)體圖像.從變形前后圖像中各抽取一個(gè)小的體積塊，變形前體積塊稱為參考體積元f(x，y，z)，變形后的稱為目標(biāo)體積元g(x，y，z)，如圖1所示.

圖1 數(shù)字體圖像相關(guān)位移模型Fig.1 The model of digital image correlation displacement

測(cè)量體圖像中某一個(gè)點(diǎn)p(x，y，z)的三維位移，通常是選擇以該點(diǎn)為中心的一個(gè)包括(2N+1)×(2N+1)×(2N+1)個(gè)體素的體積塊的位移，稱為參考體積元，中心點(diǎn)坐標(biāo)為p(x，y，z)，目標(biāo)體圖像中的目標(biāo)體積元，其中心點(diǎn)坐標(biāo)為p＇(x，y，z)，f，g為圖像的灰度，變形前后產(chǎn)生的三維位移為Δx，Δy，Δz.u，v，w代表整體素位移，Δu，Δv，Δw代表亞體素位移.體積元做微小移動(dòng)近似剛體運(yùn)動(dòng)，變形前后同一點(diǎn)的灰度保持不變，f(p)=g(p＇).通過(guò)在目標(biāo)體圖像中逐體素的搜索與參考體積元最大相關(guān)的目標(biāo)體積元，求出變形前后的位移.最大相關(guān)是相關(guān)函數(shù)C的最大值:

式(1)中，fm和gm分別為參考和目標(biāo)體積元的平均灰度，當(dāng)C值越大越接近1，則表示相似程度越大，兩個(gè)體積元的灰度分布越接近，得到的整體素位移u，v，w值會(huì)越精確，但此時(shí)只能搜索到整體素位移，屬于粗測(cè)量，為達(dá)到更好效果，提高測(cè)量精度，所以引進(jìn)亞體素的位移測(cè)量方法.

1.2 三元三次插值

如果在上述體積元內(nèi)有一點(diǎn)，它不在整體素點(diǎn)(x，y，z)位置，而是處于亞體素位置.于是每次逐步搜索的不再是整體素，而是亞像素.而體圖像在亞像素位置并不存在灰度值，所以每次移動(dòng)搜索前，要對(duì)亞體素所在的位置和該點(diǎn)灰度大小進(jìn)行重建.本文采用高精度的三元三次灰度插值法，原理如圖2所示.圖2(a)表示所要插值點(diǎn)周圍的一個(gè)單元組成，大小為4×4×4個(gè)體素.圖2(b)是從z軸上方向下看得到的xy平面圖，每一個(gè)xy平面層都通過(guò)雙三次插值得到一個(gè)亞體素構(gòu)造點(diǎn).圖2(c)是經(jīng)過(guò)雙三次插值之后，每一層都有亞體素重建，通過(guò)一次一元三次插值，就能得到最終預(yù)插值點(diǎn)的亞體素重建.

圖2 亞體素重建中三元三次插值原理圖Fig.2 Schematic principle diagram of tricubic interpolation of subvoxel reconstruction

雙立方插值［12］，稱立方卷積插值法，是一種高精度插值方法，它考慮了待求點(diǎn)周圍16個(gè)整像素點(diǎn)的灰度值，還考慮了這16個(gè)點(diǎn)在決定待求點(diǎn)灰度值時(shí)的權(quán)重，權(quán)重與這些整像素點(diǎn)離該點(diǎn)的距離密切相關(guān)，滿足某一插值內(nèi)核函數(shù)S(w)，表達(dá)式以及雙三次插值公式如下:

其中，A，B，C均為矩陣，形式如下:

其中f(i，j)表示源圖像在該像素點(diǎn)的灰度值.u，v均為在［0，1)區(qū)間的浮點(diǎn)數(shù)，分別表示待插值點(diǎn)與最臨近像素點(diǎn)在水平和豎直方向的距離.通過(guò)圖2(b)中的雙三次灰度插值，再構(gòu)造一元三次插值函數(shù)把前面得到在z軸方向上的4個(gè)亞像素插值點(diǎn)灰度值代入公式，解出插值函數(shù)各個(gè)系數(shù)，從而得到欲插值點(diǎn)的灰度值.

1.3 基于梯度的亞體素位移算法

在三維位移場(chǎng)測(cè)量中，體積元做微小位移u，v，w，視為剛體運(yùn)動(dòng)，于是變形前后同一點(diǎn)的灰度值近似不變，此時(shí)有

將式(5)在g(x+u，y+v，z+w)處進(jìn)行一階泰勒展開，得到

其中，gx，gy，gz分別為g(x+u，y+v，z+w)沿x，y，z方向的一階灰度梯度，參考文獻(xiàn)綜合比較各種梯度算子，選用Barron算子計(jì)算灰度梯度，能得到比較高的精度.通過(guò)對(duì)式(6)采用最小二乘法原理并取得駐值，最終求得亞體素位移為Δ(Δu，Δv，Δw)，被測(cè)體圖像的三維位移場(chǎng)為整體素位移與亞像素位移之和，所以被測(cè)的整體素位移為Δx，Δy，Δz.

2 算法仿真

2.1 3D數(shù)字散斑體圖像

首先由計(jì)算機(jī)模擬生成一幅3D散斑圖作為變形前的散斑圖［13］，然后對(duì)于體圖像中每一點(diǎn)，根據(jù)預(yù)設(shè)位移找到變形后其對(duì)應(yīng)的位置，將該點(diǎn)的灰度值作為該位置處的灰度值.以高斯函數(shù)模擬散斑顆粒的光強(qiáng)分布，如圖像背景光強(qiáng)均勻，則變形前后3D散斑體圖像灰度分布函數(shù)如式(8)，其中，s為空間散斑顆粒數(shù)目，R為空間散斑顆粒大小，是變參數(shù)為第k個(gè)散斑顆粒中心的隨機(jī)分布光強(qiáng).(xk，yk，zk)和(x′k，y′k，z′k)分別是變形前后第k個(gè)空間散斑顆粒的中心位置.

對(duì)參考體圖像施加任意的平移公式為

式中u，v，w，Δx，Δy，Δz為預(yù)先設(shè)定的x，y和z方向整體像素位移分量和亞像素位移分量，據(jù)式(9)對(duì)參考圖像施加任意平移，平移前圖像如圖3所示.

圖3 參考體圖像的3D圖Fig.3 3 D reference volume image

2.2 模擬剛體平移

通過(guò)計(jì)算機(jī)模擬仿真，模擬參數(shù):模擬散斑顆粒數(shù)目s=1000，散斑圖像大小為100voxel×100voxel×100voxel，剛體平移實(shí)驗(yàn)中，設(shè)定x，y，z方向上的平移量分別為u=3.28voxel，v=2.32voxel，w=2.35voxel.生成參考體圖像的3D散斑圖，如圖3所示.剛體平移5組模擬參數(shù)分別為:(1)散斑半徑R=2pixel，子體塊邊長(zhǎng)ranr為5，10，…50voxel.(2)散斑半徑R=11pixel，子體塊邊長(zhǎng)ranr為5，10，15…50voxel.(3)顆粒半徑R=2pixel，ranr為5，8，11…32voxel.(4)顆粒半徑R=11 pixel，ranr為5，8，11…32voxel.采用不同散斑顆粒大小和不同邊長(zhǎng)的子體快來(lái)做剛體平移模擬實(shí)驗(yàn)，并且計(jì)算預(yù)設(shè)值與實(shí)驗(yàn)值之間的絕對(duì)誤差、相對(duì)誤差來(lái)比較散斑半徑和子體塊大小對(duì)精度的影響.(5)顆粒R的大小為2到11之間的隨機(jī)數(shù)，子體塊邊長(zhǎng)ranr為3，5，7，…，21voxel，在每一個(gè)子體塊邊長(zhǎng)中，采用10組隨機(jī)位置的子體塊，來(lái)得到預(yù)設(shè)值與實(shí)驗(yàn)值之差的標(biāo)準(zhǔn)誤差以及總搜索的計(jì)算效率等，結(jié)果如圖4～7所示.

3 結(jié)論

本文介紹了測(cè)量物體內(nèi)部三維矢量場(chǎng)的方法基本原理簡(jiǎn)單，在梯度法的基礎(chǔ)上引入3次圖像灰度插值的方法，計(jì)算精度高.從圖4～圖7可得出，基本達(dá)到預(yù)定效果，空間整體素與亞體素之和與預(yù)設(shè)值之間的絕對(duì)誤差、相對(duì)誤差都很小，計(jì)算精度較高.通過(guò)圖4和圖5可知，在同樣大小子體塊邊長(zhǎng)下，散斑顆粒半徑R=2的性能比R=11的稍好，基本驗(yàn)證了算法的正確性;在圖6中，對(duì)于同一子體塊邊長(zhǎng)下，散斑顆粒半徑隨機(jī)，從參考體圖像中隨機(jī)抽取10個(gè)位置的子體塊來(lái)計(jì)算位移場(chǎng)與預(yù)設(shè)值之間誤差的標(biāo)準(zhǔn)誤差，標(biāo)準(zhǔn)誤差比較小;同時(shí)，對(duì)于子體塊體積大小的遞增，每次搜索和進(jìn)行相關(guān)計(jì)算的時(shí)間都會(huì)變得增加，如圖7所示，對(duì)于越大的子體塊，搜索時(shí)所要進(jìn)行灰度值相關(guān)性計(jì)算的就越多，插值點(diǎn)也相應(yīng)增加，故而時(shí)間增加.

圖4 預(yù)設(shè)值與實(shí)驗(yàn)值的絕對(duì)誤差Fig.4 Absolute error between setting values and experinental values

圖5 預(yù)設(shè)值與實(shí)驗(yàn)值的相對(duì)誤差Fig.5 Relative error between setting values and experinental values

圖6 標(biāo)準(zhǔn)誤差Fig.6 Standard error

圖7 計(jì)算效率Fig.7 Computational efficiency

［1］Franck C，Hong S，Maskarinec S A，et al.Three-dimensionalfull-field measurements of large deformations in soft materials using confocal microscopy and digital volume correlation［J］.Experimental Mechanics，2007，47(3):427-438.

［2］Gates M，Lambros J，Heath M.Towards high performance digital volume correlation［J］.Experimental Mechanics，2011，51(4):491-507.

［3］Bay B K，Smith T S，F(xiàn)yhrie D P，et al.Digital volume correlation:Three-dimensional strain mapping using X-ray tomography［J］.Experimental Mechanics，1999，39(3):217-226.

［4］Bay B K.Methods and applications of digital volume correlation［J］.Journal of Strain Analysis for Engineering Design，2008，43(8):745-760.

［5］潘兵，吳大方，謝惠民，等.基于梯度的數(shù)字體圖像相關(guān)方法測(cè)量物體內(nèi)部變形［J］.光學(xué)學(xué)報(bào)，2011，31(6):120-126.Pan B，Wu D F，Xie H M，et al.Spatial-gradient-based digital volume correlation technique for internal deformation measurement［J］.Acta Optica Sinica，2011，31(6):120-126.

［6］潘兵，吳大方，郭保橋.數(shù)字體圖像相關(guān)方 法中基于迭代最小二乘法的物體內(nèi)部變形測(cè)量［J］.實(shí)驗(yàn)力學(xué)，2011，26(6):665-673.Pan B，Wu D F，Guo B Q.Digital volume image correlation by using iterative Least-Squares for Internal deformation measurement of an object［J］.Journal of Experimental Mechanics，2011，26(6):665-673.

［7］樊雪松.數(shù)字散斑相關(guān)方法的研究［D］.天津:天津大學(xué)機(jī)械工程學(xué)院，2004.

［8］牛永強(qiáng).物體內(nèi)部三維位移場(chǎng)測(cè)量算法研究［D］.合肥:中國(guó)科學(xué)技術(shù)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，2010.

［9］夏桂鎖.數(shù)字圖像相關(guān)測(cè)量方法的理論及應(yīng)用研究［D］.天津:天津大學(xué)機(jī)械工程學(xué)院，2003.

［10］Huang J，Pan X，Li S，et al.A digital volume correlation technique for 3D deformation measurements of soft gels［J］.Inter-national Journal of Applied Mechanics，2011，3(2):335-354.

［11］汪敏，胡小方，伍小平.物體內(nèi)部三維位移場(chǎng)分析的數(shù)字圖像相關(guān)方法［J］.物理學(xué)報(bào)，2006，55(10):5135-5139.Wang M，Hu X F，Wu X P.Digital image correlation method for the analysis of 3D internal displacement field in object［J］.Acta Physica Sinica，2006，55(10):5135-5139.

［12］王會(huì)鵬，周利莉，張杰.一種基于區(qū)域的雙三次圖像插值算法［J］.計(jì)算機(jī)工程，2010，36(19):216-218.Wang H P，Zhou L L，Zhang J.Region-based bicubic image interpolation algorithm［J］.Computer Engineering，2010，36(19):216-218.

［13］Zhou P，Goodson K E.Subpixel displacement and deformation gradient measurement using digital image/speckle correlation［J］.Optical Engineering，2001，40(8):1613-1620.