DVC 中內(nèi)部散斑質(zhì)量評價及計算體素點的優(yōu)化選擇1)

鄒 翔 張軒豪 王延珺 潘 兵,2)

* (北京航空航天大學固體力學所,北京 100191)

? (中國工程物理研究院化工材料研究所,四川綿陽 621999)

引言

數(shù)字體圖像相關(guān)方法(digital volume correlation,DVC)是二維數(shù)字圖像相關(guān)方法(two-dimensional digital image correlation,2D DIC)在三維體圖像上的拓展.通過比較體成像設備獲取的被測試樣變形前后數(shù)字體圖像,該方法可測量物體內(nèi)部三維全場變形[1-2].由于可以提供被測試樣內(nèi)部(而不僅僅是表面)更為豐富的變形信息,DVC 在多孔材料(如骨[3]、金屬泡沫[4]、木材[5])、半透明生物材料(如細胞[6]等)、建筑材料(如混凝土[7]等)、復雜結(jié)構(gòu)(如3D 打印結(jié)構(gòu)[8]、生物植入物[9]等)的內(nèi)部變形場測量、力學行為表征和有限元模擬互證等方面都有著重要應用[10-12].

DVC 的基本原理是從參考體圖像中圍繞感興趣的計算點選擇一個立方體的圖像子體塊(subvolume),通過跟蹤該圖像子體塊在變形后圖像中的位置以獲得圖像子體塊中心點的三維位移矢量.在實際分析計算時,被測試樣體圖像需呈現(xiàn)隨機的灰度分布(通常稱作內(nèi)部散斑),該散斑場作為變形信息的載體隨試樣一起變形.因此被測試樣內(nèi)部散斑質(zhì)量和DVC 位移測量精度密切相關(guān).為對散斑圖質(zhì)量進行評價,在DIC 方法研究中已提出多個參數(shù)[13],這些評價參數(shù)可分為局部參數(shù)和全局參數(shù)[14-16].局部參數(shù)主要用于單個圖像子區(qū)內(nèi)散斑質(zhì)量的評價,如Pan 等[17]通過理論分析和模擬實驗證明了DIC 圖像子區(qū)的灰度梯度平方和(sum of square subset intensity gradient,SSSIG)參數(shù)與其位移測量結(jié)果的隨機誤差成反比,因此圖像子區(qū)SSSIG 參數(shù)越大,其內(nèi)部散斑質(zhì)量越好;Bomarito 等[18]將SSSIG 參數(shù)和二次自相關(guān)峰高參數(shù)結(jié)合提出了一種組合式的散斑評價參數(shù),該參數(shù)能夠解決原有SSSIG 參數(shù)不適用于規(guī)則散斑圖像(如棋盤格圖案)的情況.此外,如假設散斑質(zhì)量在整個圖像內(nèi)的分布是均勻的,則可在散斑質(zhì)量局部評價參數(shù)的基礎(chǔ)上提出可評價整幅散斑圖質(zhì)量的全局參數(shù),如Pan 等[19]在SSSIG 局部參數(shù)的基礎(chǔ)上提出了用于整幅散斑圖評價的平均灰度梯度(mean intensity gradient,MIG)全局參數(shù)等.

類似地,在利用DVC 測量內(nèi)部變形時,由于被測對象、體成像設備類型和參數(shù)以及實驗條件等的差異[20],獲取的被測物體內(nèi)部散斑圖像千差萬別.因此,如何定量評價圖像內(nèi)部散斑的優(yōu)劣也是DVC 中一個重要的基本問題.將DIC 中已有的散斑評價參數(shù)直接推廣到DVC 中是最有效的方法之一.在DIC 中的諸多散斑評價參數(shù)中,SSSIG 參數(shù)是從DIC 位移測量精度的理論模型中推導出來,具有嚴格的理論基礎(chǔ),具有普適性.并且該參數(shù)是一種典型的局部評價參數(shù),能夠克服全局參數(shù)在遇到材料內(nèi)部出現(xiàn)空洞等非連續(xù)和非均勻分布時失效的問題[21],因此更適合于材料內(nèi)部散斑質(zhì)量的評價.然而,在將SSSIG 參數(shù)在DVC 內(nèi)部散斑評價中推廣應用需要解決以下幾個問題:(1) SSSIG 參數(shù)在DVC 中最新的三維反向組合高斯?牛頓亞體素配準算法[22](3D inverse compositional Gauss-Newton,3D IC-GN)中應用的可行性問題.SSSIG 最初被提出來時是基于DIC 經(jīng)典亞像素配準算法牛頓?拉普森算法[17],因具有更優(yōu)的抗噪性能和更高的計算效率,IC-GN 算法已成為DIC/DVC 方法中的標準算法,因此系統(tǒng)地研究SSSIG 參數(shù)與該算法測量精度的關(guān)系是有必要的;(2)由于DVC 內(nèi)部散斑類型的多樣性,這就要求評價指標應具有廣泛的適用性,因此需要研究SSSIG參數(shù)在不同類型內(nèi)部散斑中的適用性;(3)與DIC 通過人工制作散斑不同,DVC 主要通過被測物體自身內(nèi)部紋理的灰度對比度來形成散斑場.由于被測物體內(nèi)部組成和結(jié)構(gòu)難以改變,無法通過人工優(yōu)化DVC 測量對象的內(nèi)部散斑質(zhì)量,因此從計算方法層面研究減弱低質(zhì)量散斑對DVC 位移測量精度的影響無疑是有價值的.

基于以上考慮,本文研究了DVC 中體圖像的內(nèi)部散斑質(zhì)量評價參數(shù)及其適用性,并通過優(yōu)化選取子體塊內(nèi)的體素點來實現(xiàn)增大子體塊尺寸以減小測量誤差的同時使計算量沒有顯著增加.論文首先從DVC 方法的位移測量精度理論分析出發(fā),系統(tǒng)介紹了SSSIG 評價參數(shù)的理論基礎(chǔ);然后,通過數(shù)值模擬平移實驗,檢驗了該參數(shù)在不同類型內(nèi)部散斑中的有效性;最后,為了在算法層面上改善DVC 子體塊的SSSIG 參數(shù),提出了基于灰度梯度的體素點選擇DVC 方法,并通過模擬實驗和實際實驗驗證了該方法的有效性.

1 DVC 位移測量誤差的理論分析及模擬實驗驗證

本節(jié)首先對最新的3D IC-GN 亞體素配準算法[22]進行簡單介紹,然后給出基于3D IC-GN 算法的DVC 位移測量誤差的理論公式推導[23-25].該公式顯示位移測量誤差與SSSIG 參數(shù)、體圖像噪聲、亞體素灰度插值算法直接相關(guān).最后通過模擬數(shù)值平移實驗驗證了SSSIG 參數(shù)用于不同類型內(nèi)部散斑質(zhì)量評價的有效性.

1.1 3D IC-GN 算法

DVC 方法在具體實現(xiàn)過程中需執(zhí)行兩個步驟:初值估計和亞體素位移測量.對任一計算點,首先需要通過簡單的整體素搜索(種子點)或可靠性導向的位移初值傳遞策略(非種子點)來確定其位移迭代初值.為進一步提高位移測量精度,需進行亞體素配準計算[26-27].在眾多DVC 亞體素配準算法中,3D ICGN 算法因具有高精度、高效率的優(yōu)點而獲得廣泛使用,該算法利用迭代運算使非線性相關(guān)函數(shù)由初始估計值逐漸收斂至局部最佳極值(即使參考和變形子體塊獲得最佳匹配的變形參數(shù)矢量)以獲得具有亞體素精度的位移計算結(jié)果.

3D IC-GN 算法計算流程圖如圖1 所示,該算法通過優(yōu)化如式(1)所示的零均值歸一化最小平方距離相關(guān)函數(shù)(zero-mean normalized sum-of-square difference,ZNSSD)獲得亞體素位移

圖1 3D IC-GN 算法計算流程圖[22]Fig.1 Flow chart of the 3D IC-GN algorithm[22]

式中f(x)和g(x) 分別表示參考和變形子體塊中點x=(x,y,z)T處的灰度值; ξ=(Δx,Δy,Δz)T表示參考子體塊中各體素點的局部坐標;fm和gm分別表示參考和變形子體塊的灰度平均值;Δf和Δg分別是參考和變形子體塊中所有體素點的灰度標準差;p=(u,ux,uy,uz,v,vx,vy,vz,w,wx,wy,wz)是描述目標子體塊的位置和形狀變化的變形矢量;Δp=(Δu,Δux,Δuy,Δuz,Δv,Δvx,Δvy,Δvz,Δw,Δwx,Δwy,Δwz)表示參考子體塊的增量變形矢量;W(ξ;p) 為目標子體塊的翹曲函數(shù);W(ξ;Δp)為參考子體塊的增量翹曲函數(shù).

采用高斯牛頓算法對ZNSSD 相關(guān)函數(shù)進行迭代優(yōu)化,則變形矢量增量可寫為

其中 ?f是參考體圖像中各點的灰度梯度,H為Hessian 矩陣.在確定變形矢量增量Δp之后,就可以確定子體塊的增量翹曲函數(shù)W(ξ;Δp).然后將該增量翹曲函數(shù)與現(xiàn)有的翹曲函數(shù)W(ξ;p) 反向組合,以更新目標子體塊的翹曲函數(shù).對于每個參考子體塊,重復進行式(1)的迭代運算,直至滿足預設的收斂條件,從而獲取最終優(yōu)化的變形矢量.此外,為表征目標子體塊中更為復雜的局部變形模式,也可采用包含二階變形分量的翹曲函數(shù)[28].

1.2 位移測量誤差的理論分析

實際實驗中獲得的體圖像不可避免含有隨機噪聲[3,28-29],且各體素點的圖像噪聲相互獨立,設參考體圖像和變形體圖像實際的灰度值分別為(x) 和(x),真實的圖像灰度值為f(x) 和g(x),則

式中 δf(x) 和 δg(x) 分別表示在參考體圖像和變形體圖像中的圖像噪聲(假設標準差分別為 σf(x) 和σg(x))引起的灰度偏差.因此,在實際實驗中,式(2)的實際形式為

根據(jù)高斯白噪聲的統(tǒng)計特性,變形矢量增量Δp的期望和方差為

式中h(x1)=f(x1)?g(x2) 表示灰度插值偏差.

為了更加直觀地展示式(6),將其應用到沿x方向的平移工況中,由上節(jié)分析可知,如果u為子體塊真實變形量,當采用3D IC-GN 算法迭代1 次時,則位移結(jié)果為u′≈u+Δu+uxΔu+uyΔv+uzΔw≈u+Δu,此時移u測量誤差ue=u′?u≈Δu,因此由式(6)可知

由式(7)可知,當不考慮形函數(shù)欠匹配問題時,在3D IC-GN 算法中位移u的系統(tǒng)誤差E(Δu) 和隨機誤差Var(Δu) 均與子體塊灰度梯度和(即SSSIG 參數(shù))大小有關(guān).但系統(tǒng)誤差主要受灰度插值誤差h(x1) 影響,當采用精度較高的插值算法,并對圖像進行預濾波,3D IC-GN 算法的系統(tǒng)誤差可控制在較低水平,而隨機誤差Var(ue) 與子體塊SSSIG 參數(shù)呈反相關(guān)關(guān)系.也就是說,在DVC 計算中,某方向的SSSIG 參數(shù)越大,其位移測量結(jié)果越精確.

1.3 內(nèi)部散斑質(zhì)量評價參數(shù)的模擬實驗驗證

為驗證上節(jié)推導的理論公式的正確性,本文對3 種不同類型內(nèi)部散斑進行了數(shù)值模擬平移實驗.與實際實驗相比,數(shù)值模擬實驗不但可以精確控制散斑圖的位移和應變信息,還能夠排除由體圖像采集系統(tǒng)(如設備自熱)、非理想加載條件等非理想實驗因素對位移測量的影響[15],最大程度地保證了計算結(jié)果的可靠性.驗證實驗選取的參考體圖像是由不同X-ray CT 獲取的3 種具有代表性的材料(銅顆粒填充環(huán)氧樹脂基復合材料、柚子皮多孔材料和多聚物粘結(jié)糖)的體圖像,圖2 給出了3 種材料體圖像的切片圖、以及選取的感興趣區(qū)域(volume of interest,VOI)及其灰度分布直方圖.

圖2 3 種不同類型內(nèi)部散斑的切片圖、VOI 及其灰度分布直方圖Fig.2 Slice figures,VOI and histograms of three speckle patterns

在平移實驗中,利用傅里葉變換對3 種材料的體圖像做精確的位置平移.具體為:將原體圖像作為參考體圖像,參考體圖像沿x方向平移0.2 體素而得到平移變形體圖像.為模擬真實情況下的圖像噪聲,在參考和變形體圖像中分別加入標準差為3 灰度值的高斯白噪聲.DVC 分析中均選取大小為200×200×200 體素的體圖像中心區(qū)域VOI 進行研究,計算采用3D IC-GN 亞體素配準算法和一階形函數(shù),子體塊尺寸為41×41×41 體素,計算步長為10 體素.

圖3(a)給出了3 種材料位移u測量結(jié)果的標準差(standard deviation (SD) error)與子體塊圖像x方向SSSIG (簡記為SSSIGx)參數(shù)平均值的關(guān)系的理論結(jié)果和測量結(jié)果,由表可知,在不同材料中,位移測量結(jié)果的標準差與理論預測值吻合良好;不同材料間,散斑圖像SSSIGx參數(shù)越大的材料,位移測量結(jié)果標準差也較小,這和上面的理論公式相一致.需要說明的是,由式(7) 可知,均值誤差不僅與圖像SSSIG 有關(guān)還與子體塊內(nèi)各體素點灰度插值偏差有關(guān),而灰度插值偏差在分析中難以統(tǒng)計,因此本文不再對均值偏差進行重點分析.但是為了降低位移測量結(jié)果的均值誤差,在DVC 計算前需要對體圖像進行低通預濾波[31].

圖3 3 種不同類型內(nèi)部散斑測量結(jié)果Fig.3 Measurement results of three speckle patterns

圖3(b)還給出了3 種材料中子體塊測量值u與其子體塊散斑圖像參數(shù)SSSIGx之間的關(guān)系,從圖3中可以看出:盡管各材料之間測量結(jié)果差別較大,然而對同一散斑材料中,測量結(jié)果的變化幅度均整體上與散斑圖像SSSIGx參數(shù)呈現(xiàn)負相關(guān)關(guān)系,即隨子體塊散斑圖像參數(shù)SSSIGx的增大,其測量結(jié)果變化幅度越小;此外,由圖3 可知,3 種材料中,柚子皮多孔材料的位移標準差最大,其次是銅顆粒填充環(huán)氧樹脂基復合材料,而多聚物粘結(jié)糖最小.

以上的研究結(jié)果表明,散斑圖像SSSIGx參數(shù)能夠反映出材料內(nèi)部散斑質(zhì)量,對同一種材料而言,其散斑圖像SSSIGx值越大,則DVC 位移測量精度越高.

2 DVC 方法中計算體素點的優(yōu)化選擇

2.1 SSSIG 參數(shù)的優(yōu)化

為提高位移測量精度,在DIC 中通常采用增大子區(qū)尺寸來優(yōu)化子區(qū)SSSIG 參數(shù)[17,21],該方法也可以推廣到DVC 中.但由于DVC 處理的是三維數(shù)據(jù),簡單地通過增大子體塊尺寸來優(yōu)化子體塊SSSIG 參數(shù)將大幅度增加計算量,進而影響計算效率.通常,在DVC 計算時子體塊內(nèi)所有的體素點均參與相關(guān)性計算,這將耗費大量的計算時間.另外一方面,通過式(7)可以看出子體塊的SSSIG 參數(shù)實質(zhì)上其內(nèi)部所有體素點灰度梯度平方的總和,子體塊內(nèi)灰度梯度小的體素點對提高SSSIG 值影響較小,即這些點對子體塊計算精度影響較小,因此可在增大子體塊尺寸的同時將子體塊內(nèi)灰度梯度平方偏小的體素點剔除出DVC 計算,以實現(xiàn)在增大子體塊SSSIG的同時盡可能的降低計算量.在實際實驗需要對3 個方向均進行優(yōu)化,因此選擇了3 個方向的SSSIG參數(shù)的均值(簡記為SSSIGmean)作為目標函數(shù),見式(8).

為了更加直觀地展示材料中灰度梯度的分布情況,圖4 給出了多孔材料中一子體塊內(nèi)部灰度梯度和的分布情況,由圖可知,在子體塊內(nèi)部有較大一部分體素點的灰度梯度平方和處于較低的水平,因此將這些體素點排除出DVC 計算將能夠有效減輕計算量.

圖4 多孔材料灰度梯度分布圖Fig.4 Intensity gradient distribution of porous materials

2.2 體素點優(yōu)化選擇方法

基于上述的SSSIG 參數(shù)優(yōu)化思想,提出了一種基于灰度梯度的計算體素點優(yōu)化選擇方法,圖5 給出了該方法在計算一子體塊位移時采取的計算流程圖.該方法的具體計算步驟為:

圖5 體素點優(yōu)化選擇方法計算流程圖Fig.5 Calculation flow chart of voxel selection DVC method

(1) 在每一個子體塊開始DVC 計算前首先計算其內(nèi)部所有體素點3 個方向的灰度梯度平方的總和?f,只有總和大于所設定灰度梯度閾值 ?f0的體素點(有效點)才會參與到該子體塊SSSIGmean參數(shù)的計算和之后的DVC 計算;

(2) 計算該子體塊SSSIGmean參數(shù),將該子體塊SSSIGmean參數(shù)與所設定的SSSIGmean參數(shù)閾值T0對比,如果低于該閾值,則將該子體塊半邊長M增大1 體素;

(3) 重復以上步驟,直到子體塊SSSIGmean參數(shù)大于所設定的閾值T0,然后基于選定的體素點,完成該子體塊的DVC 計算.

在該方法中,體素點灰度梯度閾值 ?f0和SSSIGmean參數(shù)閾值T0的設定是重點,在子體塊中,3 個方向的灰度梯度平方的總和小于灰度梯度閾值?f0的體素點將被刪除,而只有高于灰度梯度閾值?f0的體素點才會參與到DVC 計算.灰度梯度閾值?f0過小,將起不到提高測量精度的作用,灰度梯度閾值 ?f0過大,將會使子體塊有效體素點大大減小,進而影響測量結(jié)果.同時,SSSIGmean參數(shù)閾值T0設置過大,則使計算量增加,太小又起不到優(yōu)化作用,因此實際應用中要根據(jù)不同的情況選擇合適的灰度梯度閾值 ?f0和SSSIGmean參數(shù)閾值T0.

2.3 模擬實驗驗證與討論

為驗證體素點優(yōu)化選擇方法的有效性,將該方法應用于1.3 節(jié)中模擬平移實驗,通過3 種材料原有的?f的分布情況確定了其灰度梯度閾值 ?f0和SSSIGmean參數(shù)閾值T0,其中子體塊灰度梯度閾值 ?f0的取法為:將子體塊所有體素點(記為N)按灰度梯度總和?f大小遞增排列,選用編號N1=int(N/2) 體素點對應的灰度梯度總和為閾值 ?f0,SSSIGmean參數(shù)閾值T0設為該子體塊原有SSSIGmean的1.5 倍.

DVC 計算中,子體塊尺寸設為41×41×41 體素,計算步長10 體素.圖6 給出了常規(guī)DVC 和體素點選擇DVC 位移測量結(jié)果的對比情況.

圖6 兩種方法測量結(jié)果對比Fig.6 Comparison of the results of two methods

由圖可知,常規(guī)DVC 中每個子體塊中體素點的數(shù)目是一定的,具體為41×41×41=68 921 個體素點,而在體素點選擇DVC 中,盡管子體塊尺寸增加了,但3 種材料子體塊中有效體素點的平均數(shù)量依次為58 619,61 016 和64 521,比常規(guī)DVC 略低,這說明優(yōu)化后的DVC 方法計算量不僅沒有增加,還存在一定程度的降低.而在測量誤差方面,優(yōu)化后3 種材料散斑圖像的測量標準差分別下降了15.63%,25.00%和18.75%.通過以上的分析結(jié)果可知,體素點優(yōu)化選擇方法能夠在保證計算精度的前提下有效降低由于子體塊尺寸增大所引起的計算量的急劇增加.

2.4 泡沫銅材料重掃描實驗

為進一步檢驗體素點選擇DVC 方法的有效性,將該方法應用到泡沫銅材料重掃描實驗.

重掃描實驗是DVC 中一種常用的測量精度評價方法,在重掃描實驗中,被測試件靜止,且不對其施加任何載荷,利用體成像設備對試件進行兩次連續(xù)掃描,通過DVC 分析兩次掃描獲得的體圖像可以獲得試件在靜止狀態(tài)下的變形量測量值,通過此實驗可以檢驗DVC 方法的測試精準度.實驗所用試件是孔隙率為97.5%泡沫銅材料,試件尺寸為11 mm×11 mm×8 mm.體成像設備選用的是 Skyscan1172 臺式 X 射線 CT 系統(tǒng)(Bruker Corporation,美國),空間采樣率選取為 2000×1336 像素,體素尺寸設定為9 μm,旋轉(zhuǎn)步長為 0.4°.每次CT 掃描持續(xù) 45 min,共獲取試件2 組灰度體圖像.

圖7 給出了CT 實驗裝置圖和試件切片圖.在DVC 分析中,采用 3D IC-GN 亞體素配準算法和一階形函數(shù).感興趣計算區(qū)域選在圖像中心,大小為300×300×500 體素,子體塊尺寸預設為41×41×41 體素,計算步長為20 體素,共有15×15×25=5625 個規(guī)則分布的計算點.體素點選擇采用了和3.3 節(jié)一致的灰度梯度閾值 ?f0和SSSIGmean參數(shù)閾值T0.

圖8 對比給了采用常規(guī)和體素點選擇DVC 計算獲得的位移場,由圖可以很直觀地看出,采用體素點選擇DVC 的測量誤差要小于常規(guī)DVC.表1 還給出了兩種方法計算量和位移測量結(jié)果的詳細對比情況.由表1 可知,優(yōu)化后位移u,v,w的標準差分別下降了22.80%,27.73%和28.55%;與此同時,優(yōu)化前后各位移分量的均值誤差變化較小.由此可見,在泡沫銅重掃描實驗中,本文提出的體素點選擇DVC可有效降低DVC 的位移測量誤差,但沒有明顯增加其計算量.

圖8 重掃描實驗結(jié)果Fig.8 Rescan results

圖8 重掃描實驗結(jié)果 (續(xù))Fig.8 Rescan results (continued)

表1 重掃描實驗位移測量結(jié)果對比Table 1 Comparison of displacement measurement results in rescan experiment

3 結(jié)論

提出了適用于DVC 的內(nèi)部散斑質(zhì)量評價參數(shù)?SSSIG,并通過模擬平移實驗驗證了SSSIG 參數(shù)在3 種(銅顆粒填充環(huán)氧樹脂基復合材料、柚子皮多孔材料和多聚物粘結(jié)糖)不同類型內(nèi)部散斑圖像中的有效性.結(jié)果表明,在DVC 中,子體塊散斑圖像SSSIG 參數(shù)越大,測量誤差越小;更進一步,為了在增大子體塊SSSIG 參數(shù)的同時盡可能的降低計算量,提出了一種體素點優(yōu)化選擇方法,該方法通過將子體塊中灰度梯度較小的體素點剔除出DVC 計算來實現(xiàn)在增大子體塊尺寸的同時而不增加計算量,從而實現(xiàn)高精度、高效率的DVC 分析.3 種不同類型材料散斑圖像的模擬平移實驗和泡沫銅的重掃描實驗均驗證了體素點優(yōu)化選擇方法在有效降低DVC方法的測量誤差同時提高計算效率.