DVC 中內部散斑質量評價及計算體素點的優化選擇1)

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

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

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

引言

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

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

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

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

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

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

1.1 3D IC-GN 算法

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

3D IC-GN 算法計算流程圖如圖1 所示,該算法通過優化如式(1)所示的零均值歸一化最小平方距離相關函數(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) 為目標子體塊的翹曲函數;W(ξ;Δp)為參考子體塊的增量翹曲函數.

采用高斯牛頓算法對ZNSSD 相關函數進行迭代優化,則變形矢量增量可寫為

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

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

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

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

根據高斯白噪聲的統計特性,變形矢量增量Δp的期望和方差為

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

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

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

1.3 內部散斑質量評價參數的模擬實驗驗證

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

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

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

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

圖3 3 種不同類型內部散斑測量結果Fig.3 Measurement results of three speckle patterns

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

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

2 DVC 方法中計算體素點的優化選擇

2.1 SSSIG 參數的優化

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

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

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

2.2 體素點優化選擇方法

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

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

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

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

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

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

2.3 模擬實驗驗證與討論

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

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

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

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

2.4 泡沫銅材料重掃描實驗

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

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

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

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

圖8 重掃描實驗結果Fig.8 Rescan results

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

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

3 結論

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