基于網(wǎng)格點(diǎn)擴(kuò)散的二維DIC 連續(xù)變形測(cè)量方法研究

高山，陳泰錚，陳立偉，崔穎

哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的快速興起，圖像處理技術(shù)得到了蓬勃發(fā)展，憑借著其獨(dú)特的優(yōu)勢(shì)融入各個(gè)研究領(lǐng)域當(dāng)中[1-2]。在此背景下，數(shù)字圖像相關(guān)法（digital image correlation, DIC）逐漸進(jìn)入人們的視野。DIC 是一種非接觸式變形測(cè)量方法，通過(guò)比較形變前后物體表面圖像的灰度變化計(jì)算得到全場(chǎng)的位移與應(yīng)變。相比于傳統(tǒng)接觸式（如傳感器、應(yīng)變片等）應(yīng)變測(cè)量方法，DIC 具有操作簡(jiǎn)便、適應(yīng)性強(qiáng)、計(jì)算迅速、測(cè)量準(zhǔn)確等優(yōu)勢(shì)，目前已作為一種有效的光學(xué)測(cè)量方法廣泛應(yīng)用于固體材料的變形測(cè)量中[3-6]。

傳統(tǒng)的DIC 算法以遠(yuǎn)小于圖像尺寸的子集為單位進(jìn)行相關(guān)搜索，保證了測(cè)量精度卻限制了計(jì)算速度。為了改善搜索性能，粗細(xì)搜索法、爬山搜索法等多種方法被提出，然而在位移場(chǎng)測(cè)量中仍須多次使用搜索算法。因此，提高DIC 算法整體的運(yùn)算速度依舊是一個(gè)重要的問(wèn)題。

本文對(duì)連續(xù)變形區(qū)域位移場(chǎng)與應(yīng)變場(chǎng)的計(jì)算方法進(jìn)行研究，實(shí)現(xiàn)了基于網(wǎng)格點(diǎn)擴(kuò)散的二維DIC 變形測(cè)量；將該算法與傳統(tǒng)的逐網(wǎng)格點(diǎn)相關(guān)匹配的DIC 方法進(jìn)行對(duì)比，分析在測(cè)量二維連續(xù)變形時(shí)，網(wǎng)格點(diǎn)擴(kuò)散法在計(jì)算速度及位移場(chǎng)計(jì)算精度上的優(yōu)勢(shì)。

1 基于DIC 的應(yīng)變場(chǎng)測(cè)量原理

相關(guān)匹配是DIC 算法的核心，由整像素位移搜索以及亞像素位移計(jì)算2 部分組成；在相關(guān)匹配的基礎(chǔ)上，使用擴(kuò)散法快速遍歷所有網(wǎng)格點(diǎn)；對(duì)離散的位移值進(jìn)行插值得到全場(chǎng)的位移，使用多項(xiàng)式擬合進(jìn)行平滑，通過(guò)差分計(jì)算應(yīng)變場(chǎng)。

1.1 相關(guān)匹配

通過(guò)圖像子集的相關(guān)匹配計(jì)算，可得到樣件表面同一點(diǎn)在形變前后圖像中的位置，將二者的位置坐標(biāo)相減便可獲得該點(diǎn)以像素為單位的位移值。為了獲得更精確的位移結(jié)果，相關(guān)匹配分為整像素算法和亞像素算法。

1.1.1 整像素位移計(jì)算

數(shù)字圖像相關(guān)法在處理過(guò)程中針對(duì)有一定散斑特征的灰度圖像進(jìn)行處理，對(duì)圖像中某點(diǎn)的相關(guān)匹配是圍繞以該點(diǎn)為中心的圖像子集進(jìn)行的。普遍使用的相關(guān)函數(shù)有2 種[7]：零均值歸一化互相關(guān)（zero-normalized cross-correlation, ZNCC）及歸一化最小二乘互相關(guān)（zero-normalized sum of squared differences criterion, ZNSSD），均具有較高的可靠性與靈敏度。當(dāng)子集形狀是半徑為n、邊長(zhǎng)為(2n+ 1)的矩形窗口時(shí)，以子集中心為原點(diǎn)的2 種相關(guān)函數(shù)數(shù)學(xué)表示為

式中：i、j為表示子集點(diǎn)在子集窗口中沿X與Y坐標(biāo)軸方向的上坐標(biāo)序列；xi和yj為形變前子集點(diǎn)的坐標(biāo)；和為形變后子集點(diǎn)的坐標(biāo)；函數(shù)f()和g()分別返回的是形變前后某一像素點(diǎn)的灰度值；fn與gn分別為形變前后子集中圖像點(diǎn)的平均灰度值，二者的計(jì)算結(jié)果對(duì)子集之間的相關(guān)性有不同的描述，ZNCC 的計(jì)算結(jié)果最大為1，且越接近于1 相關(guān)性越強(qiáng)，小于1 則相關(guān)性減弱，直到-1 時(shí)相關(guān)性最低；而ZNSSD 的結(jié)果接近0 時(shí)相關(guān)性最強(qiáng)。

考慮形變前的散斑圖像（參考圖像）和形變后的散斑圖像（形變圖像），當(dāng)要計(jì)算參考圖像上某一點(diǎn)的位移時(shí)，以該點(diǎn)為中心建立一定大小的子集窗口（參考子集），該窗口可以是矩形也可以是圓形；接下來(lái)，在形變圖像中逐像素點(diǎn)移動(dòng)該子集并計(jì)算相關(guān)函數(shù)，根據(jù)相關(guān)性大小得到參考子集在形變圖像中的位置（形變子集）；最后，將形變子集與參考子集中心點(diǎn)坐標(biāo)相減得到該點(diǎn)的整像素位移。

相關(guān)函數(shù)所識(shí)別的散斑特征可以是物體表面自然存在的隨機(jī)紋理，也可以通過(guò)對(duì)試件表面噴涂油漆或有色金屬等方法來(lái)人工獲得，且應(yīng)盡量避免散斑局部重復(fù)的情況[8]。無(wú)論如何，數(shù)字圖像相關(guān)法的一條基本前提就是散斑特性會(huì)隨著試件的變形而一同變形。

1.1.2 亞像素位移計(jì)算

為了進(jìn)一步提高DIC 測(cè)量結(jié)果的精確度，在整像素位移的基礎(chǔ)上通過(guò)非線性迭代優(yōu)化實(shí)現(xiàn)亞像素位移的計(jì)算。引入描述子集變形情況的有限元形函數(shù)[9]為

式中：一階形函數(shù)p包含6 個(gè)變形參數(shù)，其中，u、v為子集中心點(diǎn)沿圖像坐標(biāo)系X、Y方向的位移，其余為對(duì)應(yīng)的一階偏導(dǎo)。子集變形示意如圖1所示。

圖1 子集變形示意

圖1 中P點(diǎn)為參考子集中心，Q點(diǎn)為參考子集中任意一點(diǎn)，u0、v0分別為子集中心點(diǎn)的位移，P′、Q′分別為P、Q2 點(diǎn)變形后的位置。可見，確定了形函數(shù)p，便確定了子集變形的情況。在得到參考子集中心點(diǎn)的整像素的位移后，u、v被確定，將其余未定的變形參數(shù)設(shè)置為0，此時(shí)的形函數(shù)被作為迭代初值，記為p0。本文使用反向組合高斯牛頓法[10]（ inverse compositional gauss-newton method, IC-GN）對(duì)p0進(jìn)行非線性迭代。

通過(guò)引入子集變形函數(shù)W，將歸一化最小二乘互相關(guān)準(zhǔn)則ZNSSD 改寫成輸入為 Δp的函數(shù)形式；對(duì)ZNSSD 求導(dǎo)得到取極值時(shí)的 Δp，由ZNSSD的相關(guān)性描述可知，在該 Δp下形變子集的相關(guān)性達(dá)到最大；使用 Δp更新p0，完成一次迭代優(yōu)化；將更新過(guò)后的形函數(shù)重復(fù)上述過(guò)程，直到滿足迭代停止的條件為止；更新得到的最終的變形參數(shù)u、v便是亞像素位移測(cè)量結(jié)果[11-12]。子集的非線性優(yōu)化迭代過(guò)程如圖2 所示。

圖2 非線性優(yōu)化流程

1.2 位移場(chǎng)應(yīng)變場(chǎng)測(cè)量

獲取網(wǎng)格點(diǎn)的亞像素位移之后，對(duì)相關(guān)離散點(diǎn)的位移值進(jìn)行三次樣條插值，最終獲得感興趣區(qū)域（region of interest, ROI）連續(xù)的位移場(chǎng)。對(duì)網(wǎng)格點(diǎn)的步長(zhǎng)進(jìn)行合理的設(shè)置，并依托相關(guān)匹配算法準(zhǔn)確且穩(wěn)定的性能，該方法可以得到十分準(zhǔn)確的計(jì)算結(jié)果。

計(jì)算中通過(guò)對(duì)位移場(chǎng)進(jìn)行差分計(jì)算來(lái)獲得應(yīng)變場(chǎng)，而插值得到的位移場(chǎng)具有噪聲[13]，且差分計(jì)算會(huì)放大噪聲產(chǎn)生的影響，從而增大應(yīng)變場(chǎng)的測(cè)量誤差。因此，在計(jì)算應(yīng)變場(chǎng)之前要對(duì)連續(xù)位移場(chǎng)進(jìn)行中值濾波、曲面擬合等操作平滑去噪，減小對(duì)后續(xù)應(yīng)變場(chǎng)測(cè)量的影響。

本文使用格林-拉格朗日應(yīng)變形式（greenlagrangian strain）計(jì)算二維應(yīng)變場(chǎng)為

式中：Exx、Eyy分別為沿圖像坐標(biāo)系X方向與Y方向的線應(yīng)變大小；Exy為剪應(yīng)變，它們通過(guò)4 個(gè)位移梯度分量計(jì)算得到。

2 基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 方法實(shí)現(xiàn)

在實(shí)現(xiàn)基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 應(yīng)變場(chǎng)測(cè)量方法過(guò)程中，本文使用逐點(diǎn)搜索法實(shí)現(xiàn)整像素位移測(cè)量，并將整像素位移作為IC-GN 法的迭代初值進(jìn)行優(yōu)化，從而實(shí)現(xiàn)相關(guān)匹配。在使用一階形函數(shù)的IC-GN 算法中，變形函數(shù)W以及改寫過(guò)后的ZNSSD 分別為

在W中， ξ為表示子集內(nèi)各像素點(diǎn)相對(duì)于子集中心點(diǎn)偏移量 Δx、 Δy的增廣向量；p為一階形函數(shù)，其余參量為變形參數(shù)。將CLS進(jìn)行一階泰勒展開，對(duì)得到的式子求導(dǎo)并置零計(jì)算得到 Δp，結(jié)果為

式中： ?f為參考子集的圖像梯度；fn、gn分別為參考子集中所有子集點(diǎn)的平均灰度和形變子集中所有子集點(diǎn)的平均灰度；函數(shù)f與g分別返回參考圖像與形變圖像某點(diǎn)處的灰度值；自變量X為圖像坐標(biāo)的增廣向量；H為黑森矩陣（Hessian）；W的梯度矩陣、黑森矩陣H以及X、 Δf與 Δg分別為

計(jì)算得到 Δp之后，IC-GN 法通過(guò)變形函數(shù)W來(lái)更新變形參數(shù)為

直到迭代終止，變形函數(shù)W中子集中心點(diǎn)處的變形參數(shù)u、v即為該點(diǎn)的亞像素位移。

DIC 的測(cè)量結(jié)果是連續(xù)的位移場(chǎng)與應(yīng)變場(chǎng)，而整像素搜索算法相對(duì)比較耗時(shí)，因此不能采用對(duì)參考圖像中每一點(diǎn)都進(jìn)行相關(guān)匹配的操作來(lái)獲得位移場(chǎng)。常用方法是在參考圖像中選取一定大小的ROI，在其中均勻生成具有一定步長(zhǎng)的網(wǎng)格點(diǎn)，對(duì)這些點(diǎn)進(jìn)行相關(guān)匹配得到亞像素位移，再通過(guò)插值得到整個(gè)ROI 區(qū)域的位移場(chǎng)[14]。

通過(guò)減少相關(guān)匹配的次數(shù)，雖然在一定程度上提高了DIC 算法的運(yùn)行速度，但仍有較長(zhǎng)的時(shí)間花費(fèi)在網(wǎng)格點(diǎn)的整像素搜索上。

為了進(jìn)一步提高運(yùn)行速度，Pan[15]提出了一種通過(guò)種子點(diǎn)（seed）擴(kuò)散的方法來(lái)計(jì)算位移場(chǎng)：首先，在ROI 內(nèi)任意選取一個(gè)網(wǎng)格點(diǎn)作為種子，通過(guò)相關(guān)匹配得到該種子點(diǎn)的亞像素位移；接著，將亞像素位移值作為該點(diǎn)周圍4 個(gè)鄰近點(diǎn)的形函數(shù)初值進(jìn)行非線性優(yōu)化，從而得到這些點(diǎn)的位移結(jié)果；之后，將4 個(gè)鄰近點(diǎn)以同樣的方法進(jìn)行擴(kuò)散直到獲得所有網(wǎng)格點(diǎn)的亞像素位移；同時(shí)，在上述過(guò)程中每一次擴(kuò)散都通過(guò)相關(guān)函數(shù)的值來(lái)決定擴(kuò)散順序，將相關(guān)性差的點(diǎn)留到最后進(jìn)行計(jì)算。種子點(diǎn)擴(kuò)散法只需要對(duì)一個(gè)網(wǎng)格點(diǎn)進(jìn)行完整的相關(guān)匹配，如此便能省去對(duì)所有網(wǎng)格點(diǎn)進(jìn)行整像素搜索的時(shí)間，在保證了測(cè)量準(zhǔn)確性的同時(shí)，顯著提高了計(jì)算效率。同時(shí)，在擴(kuò)散中使用ZNCC 進(jìn)行相關(guān)性排序，使得該算法一并適用于陰影區(qū)域、非連續(xù)變形區(qū)域的位移、應(yīng)變測(cè)量，極大地改善了DIC 算法的性能。盡管如此，在連續(xù)變形的測(cè)量當(dāng)中，種子點(diǎn)擴(kuò)散法在計(jì)算速度上仍然具有提升空間。

本文采用網(wǎng)格點(diǎn)擴(kuò)散法完成對(duì)所有網(wǎng)格點(diǎn)的相關(guān)匹配。相比于種子點(diǎn)擴(kuò)散，該方法使用隊(duì)列的方式實(shí)現(xiàn)網(wǎng)格點(diǎn)的無(wú)順序擴(kuò)散，無(wú)需ZNCC 相關(guān)性大小對(duì)隊(duì)列中的點(diǎn)進(jìn)行排序。這種改變雖然不能用于非連續(xù)變形的測(cè)量，但卻仍然適用于連續(xù)區(qū)域，且因減少了額外的相關(guān)函數(shù)計(jì)算以及隊(duì)列排序功能，降低程序計(jì)算量及運(yùn)行時(shí)間。基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 方法程序流程如圖3 所示。

圖3 網(wǎng)格點(diǎn)擴(kuò)散法流程

3 實(shí)驗(yàn)結(jié)果分析

本節(jié)通過(guò)對(duì)數(shù)字散斑圖像進(jìn)行仿真實(shí)驗(yàn)，將基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 算法與傳統(tǒng)的逐網(wǎng)格點(diǎn)相關(guān)匹配的DIC 方法進(jìn)行比較，驗(yàn)證擴(kuò)散法在計(jì)算速度以及測(cè)量準(zhǔn)確度上的優(yōu)勢(shì)。

3.1 實(shí)驗(yàn)過(guò)程

仿真圖像設(shè)置參數(shù)為:分辨率256×256，u、v方向分別具有-1.1 pixel 以及0.7 pixel 的位移。生成的數(shù)字散斑圖及其變形情況如圖4 所示。

圖4 數(shù)字散斑圖及變形情況

將圖4 所示的參考圖像與形變圖像作為DIC程序的輸入，設(shè)置子集半徑為15 pixel、網(wǎng)格點(diǎn)間隔10 pixel。繪制不同大小的矩形ROI，對(duì)比2 種DIC 算法完成所有網(wǎng)格點(diǎn)相關(guān)匹配的用時(shí)。

3.2 實(shí)驗(yàn)結(jié)果

繪制4 種不同大小的矩形ROI，在不同網(wǎng)格點(diǎn)數(shù)量下對(duì)比2 種算法的相關(guān)匹配用時(shí)。

如表1 所示，網(wǎng)格點(diǎn)擴(kuò)散法的計(jì)算時(shí)間要遠(yuǎn)低于傳統(tǒng)的逐網(wǎng)格點(diǎn)相關(guān)匹配算法。在不同數(shù)量網(wǎng)格點(diǎn)的計(jì)算情況下，圖像分辨率為256×256時(shí)，擴(kuò)散法的計(jì)算速度幾乎為傳統(tǒng)法的10 倍。而且，整像素算法的搜索時(shí)間與圖像大小呈正相關(guān)，所以，在實(shí)際應(yīng)用中這種速度差距會(huì)更加明顯。

表1 網(wǎng)格點(diǎn)匹配時(shí)間對(duì)比

接下來(lái)，設(shè)置矩形ROI 尺寸為220×220，對(duì)2 種DIC 算法的u、v方向位移場(chǎng)測(cè)量結(jié)果進(jìn)行對(duì)比，如圖5 和圖6 所示。

圖5 傳統(tǒng)逐點(diǎn)相關(guān)匹配所得位移場(chǎng)

圖6 網(wǎng)格點(diǎn)擴(kuò)散法所得位移場(chǎng)

通過(guò)觀察，2 種DIC 算法測(cè)量所得位移場(chǎng)與散斑圖的理論位移值基本一致，u方向位移在-1.12~-1.09 pixel 起伏，總體接近-1.1 pixel；v方向位移在0.68~0.71 pixel 起伏，總體接近0.7 pixel。

對(duì)圖6 所示位移場(chǎng)進(jìn)行全場(chǎng)的二階多項(xiàng)式擬合，使用格林-拉格朗日應(yīng)變形式計(jì)算二維應(yīng)變場(chǎng)，結(jié)果如圖7 所示。

圖7 應(yīng)變場(chǎng)測(cè)量結(jié)果

在參考圖像中選取8 個(gè)位置點(diǎn)，計(jì)算這8 個(gè)位置處位移測(cè)量值的絕對(duì)誤差。如表2 所示。

表2 傳統(tǒng)逐點(diǎn)相關(guān)匹配與網(wǎng)格點(diǎn)擴(kuò)散方法測(cè)量結(jié)果pixel

對(duì)比表2 中的計(jì)算結(jié)果，2 種DIC 算法所得u、v方向的位移結(jié)果近乎一致；在圖5 和圖6 的位移場(chǎng)熱力圖也可以得到相同的結(jié)果。

如圖8 所示測(cè)量結(jié)果相對(duì)誤差曲線，實(shí)驗(yàn)點(diǎn)處u、v方向位移測(cè)量值的相對(duì)誤差均處在1%以下。

圖8 相對(duì)誤差曲線

接下來(lái)，為了從總體上把握2 種方法的測(cè)量偏差，計(jì)算得到220×220 的矩形ROI 內(nèi)共48 400個(gè)位移計(jì)算點(diǎn)的平均偏差，如表3 所示。

表3 位移平均偏差計(jì)算結(jié)果pixel

由表3 可見，2 種方法的位移平均值均接近理論值，且平均偏差在0.003 像素左右，表明各計(jì)算點(diǎn)位移值與平均值的差異程度偏小。因此，基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 方法的計(jì)算結(jié)果，與傳統(tǒng)的逐網(wǎng)格點(diǎn)相關(guān)匹配相比幾乎一致，保證了測(cè)量結(jié)果的準(zhǔn)確性。

4 結(jié)論

本文針對(duì)傳統(tǒng)的逐點(diǎn)相關(guān)匹配法計(jì)算效率較低的問(wèn)題，提出了一種高效測(cè)量連續(xù)變形的基于網(wǎng)格點(diǎn)擴(kuò)散的二維DIC 方法。通過(guò)實(shí)驗(yàn)分析，得出結(jié)論：

1）使用ZNSSD 一階泰勒展開實(shí)現(xiàn)的ICGN 算法在連續(xù)位移測(cè)量中具有較高的準(zhǔn)確性，實(shí)驗(yàn)結(jié)果相對(duì)誤差在1%以下；

2）相比于傳統(tǒng)的逐網(wǎng)格點(diǎn)相關(guān)匹配，基于網(wǎng)格點(diǎn)擴(kuò)散的DIC 算法計(jì)算時(shí)間更短，極大地提高了程序總體的運(yùn)行速度；

3）本文提出的基于網(wǎng)格點(diǎn)擴(kuò)散的二維DIC 方法只適用于連續(xù)變形的高效測(cè)量，對(duì)于更加復(fù)雜的非連續(xù)變形問(wèn)題，該方法還需進(jìn)一步完善，這也是今后研究的重點(diǎn)。