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

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

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

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

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

1 基于DIC 的應變場測量原理

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

1.1 相關匹配

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

1.1.1 整像素位移計算

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

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

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

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

1.1.2 亞像素位移計算

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

式中：一階形函數(shù)p包含6 個變形參數(shù)，其中，u、v為子集中心點沿圖像坐標系X、Y方向的位移，其余為對應的一階偏導。子集變形示意如圖1所示。

圖1 子集變形示意

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

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

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

1.2 位移場應變場測量

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

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

本文使用格林-拉格朗日應變形式（greenlagrangian strain）計算二維應變場為

式中：Exx、Eyy分別為沿圖像坐標系X方向與Y方向的線應變大小；Exy為剪應變，它們通過4 個位移梯度分量計算得到。

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

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

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

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

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

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

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

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

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

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

圖3 網(wǎng)格點擴散法流程

3 實驗結果分析

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

3.1 實驗過程

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

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

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

3.2 實驗結果

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

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

表1 網(wǎng)格點匹配時間對比

接下來，設置矩形ROI 尺寸為220×220，對2 種DIC 算法的u、v方向位移場測量結果進行對比，如圖5 和圖6 所示。

圖5 傳統(tǒng)逐點相關匹配所得位移場

圖6 網(wǎng)格點擴散法所得位移場

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

對圖6 所示位移場進行全場的二階多項式擬合，使用格林-拉格朗日應變形式計算二維應變場，結果如圖7 所示。

圖7 應變場測量結果

在參考圖像中選取8 個位置點，計算這8 個位置處位移測量值的絕對誤差。如表2 所示。

表2 傳統(tǒng)逐點相關匹配與網(wǎng)格點擴散方法測量結果pixel

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

如圖8 所示測量結果相對誤差曲線，實驗點處u、v方向位移測量值的相對誤差均處在1%以下。

圖8 相對誤差曲線

接下來，為了從總體上把握2 種方法的測量偏差，計算得到220×220 的矩形ROI 內共48 400個位移計算點的平均偏差，如表3 所示。

表3 位移平均偏差計算結果pixel

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

4 結論

本文針對傳統(tǒng)的逐點相關匹配法計算效率較低的問題，提出了一種高效測量連續(xù)變形的基于網(wǎng)格點擴散的二維DIC 方法。通過實驗分析，得出結論：

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

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

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