基于遺傳算法的數字圖像相關變形初值估計

劉 禹，肖世德*，張 睿，張若凌，張 磊

(1.西南交通大學 機械工程學院 機電測控系，成都 610031；2.中國空氣動力研究與發展中心，綿陽 621000)

引 言

數字圖像相關法(digital image correlation,DIC)是由YAMAGUCHI等人[1-2]提出的一種非接觸圖像測量方法,該測量方法光路簡單、抗干擾、適合全場測量、使用范圍廣[3]。隨著數字圖像處理技術的發展和圖像采集分辨率的提高，DIC方法在變形測量中的應用愈加廣泛[4-8]。DIC方法按其測量精度可分為整像素DIC測量方法和亞像素DIC測量方法，其中整像素DIC方法有粗細搜索法、十字搜索法、人工魚群算法、粒子群算法等相關搜索方法，亞像素DIC方法有曲面擬合、灰度梯度等。牛頓-拉夫森迭代法(Newton-Raphson,N-R)是一種常用的亞像素搜索方法，其收斂范圍一般在幾個像素[4]，但其初值估計的準確度對其收斂速度和計算精度影響較大。測量小變形物體表面時，目標點變形前后只有幾個像素的位移量，即使有測量誤差，N-R法通過一定的迭代次數仍能找到最優解[9]，但測量大變形物體表面時，通過數字圖像相關法匹配的目標點變形前后像素位移量會驟增，測量誤差會被放大，使得N-R迭代次數增加且增大計算量，甚至可能陷入局部峰值[10-13]。

遺傳算法(genetic algorithm,GA)是一種全局搜索方法，該方法搜索速度快，測量結果準確，適用于大變形對象的整像素位移測量[14-17]。該方法的計算時間受搜索區域增大的影響小，是一種有效的大范圍搜索方法，但在末端定位時存在局部震蕩現象。本文中提出一種基于遺傳算法的數字圖像相關變形初值估計法。首先，選取待測數據點鄰域內若干估測點，通過一種遺傳算法和粗細搜索法的混合方法進行整像素搜索，匹配出變形前后估測點的坐標值變化。等概率選取3組或以上不共線的估測點對進行仿射變換并擬合出具有6個參量的仿射變換模型，依據仿射變換結果得到待測數據點的變形初值估計。最后，將變形初值代入N-R迭代算法中，結合雙三次樣條插值方法獲取最終更精確的亞像素位移值。仿真和橡膠圓柱壓縮實驗結果表明，本文中的方法有效減少匹配時間，在大變形情況下測量精度仍然保持穩定，與傳統方法相比有明顯優勢。

1 DIC的基本原理

物體表面常存在自然地斑點紋理，以此作為特征可以估算兩幅圖片的相似度。DIC方法是以人工制作的細小顆粒狀散斑圖像為特征，匹配變形前后兩幅圖像上對應子區相似程度的方法。其基本原理如圖1所示。

Fig.1 Subset region of image before and after deformation

采集變形前后兩幅圖像進行對比，變形前圖像記為參考圖像，變形后圖像記為變形圖像。在參考圖像中選取的一定大小的矩形參考子區，其中心為P(x0,y0)，依據變形前后像素點灰度值不變的原則，變形圖像必然存在一與參考子區相似度最高的子區圖像，記為變形子區，其中心點坐標為P′(x0′,y0′)，參考子區和變形子區的相似度通常采用相關函數來度量，依據TONG等人[18-19]對相關系數的研究可知，零均值歸一化互相關函數對圖像灰度值的線性變化不敏感，適合在光照環境有微弱變化的環境下使用，一定程度降低外界光照干擾對匹配結果的影響，本文中采用該方法，表達式為:

C=

(1)

式中，R為圖像子區半徑，f(x,y)為參考圖像中點(x,y)處的灰度值；g(x′,y′)為變形圖像中點(x′,y′)處的灰度值；〈f〉和〈g〉為對應子區的灰度平均值。變形子區和參考子區的匹配度與相關系數C值正相關，C取值在[0,1]之間，C=1代表二者完全匹配。

2 基于遺傳算法的數字圖像相關法

為了解決傳統DIC方法迭代運算受初值影響較大且容易陷入局部最優的問題，提出基于遺傳算法的數字圖像相關法。首先，根據待測點坐標選取鄰域內3組不共線的估值點，通過遺傳算法在大區域匹配各個估值點在變形后的空間域信息，得到待測數據點周圍變形前后的若干估值點對，以此為依據根據仿射變換原理求取待測數據點位移變化估計初值和待測數據點變形后坐標值。最后，將數據點位移變化估計初值代入N-R迭代算法求取最終亞像素位移值。算法流程如圖2所示。

2.1 遺傳算法描述

遺傳算法是模擬自然生物進化過程的算法，本文中選取100組初始解S=[u,v]作為初始種群，其中u,v分別代表x,y方向上的位移值，其取值為[-25,25]間隨機數。通過一定方式對其編碼，對個體進行選擇、交叉、變異等操作，依據適應度函數計算數值篩選出性能優良的個體并保留其基因至下一代，代代之間不斷優化，直至滿足最優條件。

輪盤賭法是遺傳算法中最常用的選擇方法之一，但該方法選擇過程中有一定概率跳過最優解。本文中將輪盤賭法和最佳保留選擇法結合使用，經過輪盤賭法選擇后如果當代中最優個體沒有滿足要求，則進行最佳保留策略，選出n個個體中最差的5個，用當代最優的5個個體將其替換。

遺傳算法中通過兩個體之間染色體交叉產生新個體。交叉算子有簡單交叉、啟發式交叉、算術交叉等，此處采用算數交叉，交叉產生的兩個新的個體表示為:

[Si1Si2]=β[S1S2]+(1-β)[S2S1]

(2)

式中,問題的規模i=1,2，…,n；Si為種群中編號為i的個體，β為[0,1]之間隨機數。此處交叉發生概率設為P1=0.8。

變異操作發生概率較小，此處設其概率為P2=0.01。變異算子常有均勻變異、非均勻變異、多維正態變異、邊界變異等，此處采用非均勻變異操作。取第l次迭代時種群中的個體Sl進行變異,變異生成的新個體表達式為:

(3)

式中,L為最大迭代次數；ΝM×M為對角矩陣且其對角線元素為[0,1]間隨機數；M為個體染色體向量的維數；Du為染色體分量最大值，Dl為染色體分量最小值；b代表對迭代數依賴程度。

經過上述選擇、交叉、變異產生的下一代群體，若最優個體的適應度值不滿足要求，則重新進行個體適應度評估并進入下一輪循環，直到最優個體滿足要求。遺傳過程的停止條件為兩代種群中的最優個體和當代種群中的最優個體相同或者進行到最大迭代次數。

2.2 遺傳算法與粗細搜索法的混合算法

傳統遺傳算法進行整像素搜索時，實現了全局尋優搜索，但其收斂位置受到交叉和變異操作的影響，有一定概率偏離真實位置1～2個像素，造成N-R迭代收斂效率的降低。本文中將遺傳算法和粗細搜索法結合，如圖3所示。首先通過遺傳算法在大范圍內進行全局搜索，借助其快速收斂的優勢，在短時間內找到真實位置的初步估算位置A。遺傳算法存在末端收斂局部震蕩問題，直接使用遺傳算法的搜索結果會存在一定像素偏差，故引入粗細搜索思想，在最終末端定位時以A為中心，選取周圍5×5像素的子區逐點搜索使得相關函數達到最大值的點作為最終結果，有效消除末端收斂的局部震蕩現象，提升整像素搜索的穩定性。

針對遺傳算法有一定概率陷入局部最優解的現象，本文中對比前后兩個待測點檢測位移值之差，若超出設定閾值，則重新計算，從而有效避免陷入局部最優。

Fig.3 Schematic of hybrid algorithm

2.3 仿射變換初值估計

在參考圖像中選取感興趣區域，以51×51[11]像素大小對該區域進行網格劃分，得到若干窗口，每個窗口的中心即為一個待測點。以其中一個待測點為例，取其鄰域內的若干待測點作為估值點，通過遺傳算法計算得到第i個估值點變形前后的位置坐標分別為hi=(xi,yi)T和hi′=(xi′,yi′)T。隨機選取3組估值點對代入仿射變換，其表達式為:

(4)

表達式中通過3組估值點對的坐標信息既可求解包含6個參量a1,a2,a3,a4,a5,a6的仿射變換。通過仿射變換即可求得變形后待測點的精確坐標位置，從而得到相對變形前待測點的像素位移估值，并以此作為N-R迭代的初始值。該方法相較于傳統方法，計算出的位移初值更接近待測點的實際位移，有效降低了后續N-R迭代運算次數，縮短收斂時間。

2.4 亞像素位移計算

為提高測量精度，常需要采用曲面擬合法、灰度梯度法、N-R迭代等方法計算亞像素位移值。其中N-R迭代法是一種常用的獲取亞像素位移值的方法，其具有精度高、計算結果穩定可靠等優點，因此，本文中選用該算法計算亞像素位移值。進行N-R迭代之前需要通過灰度差值方法獲取亞像素位置的灰度值和灰度梯度等信息，本文中采用精度較高的雙三次樣條插值法，其表達式為:

(5)

式中,f(x,y)為插值區域中位于(x,y)處的灰度值，aij為插值系數，對f(x,y)求偏導即得到各方向灰度梯度值。通過N-R迭代運算后求取的亞像素位移值即為最終結果。

3 實驗驗證

3.1 模擬散斑圖仿真實驗

為了驗證本文中提出算法在計算精度和計算性能上面的提升，在計算機上使用本文中算法和其它傳統算法對模擬散斑圖像進行特征點匹配，計算亞像素位移值。參考ZHOU[20]提出的模擬散斑圖模擬物體表面的變形過程，變形前后的模擬散斑圖灰度表示為:

I1(x,y)=

(6)

(7)

式中，s為散斑顆粒數量，r為散斑顆粒半徑，(x0,y0)為散斑圖中心位置，I0為光強分布;u,ux,uy,v,vx,vy為散斑的變形參量，決定圖像的位移量和應變量。為驗證不同的紋理對測量的影響，將散斑顆粒大小r、顆粒數目s分別設置為3pixels、2500和2pixels、4500兩組，如圖4所示。每組vx應變值按照0.005的步距由0.001過渡到0.05。

分別利用傳統N-R迭代法、曲面擬合法與本文中算法對目標圖像中隨機數據點進行匹配。比較不同散斑顆粒尺寸、不同散斑顆粒數目和不同應變值對測量平均誤差和標準差的影響，結果如圖5、圖6所示。最后對迭代次數和匹配時間進行對比，如表1所示。

Fig.4 The simulated speckle image generated by MATLAB

a—speckle parameterr=3pixels,s=2500 b—speckle parameterr=2pixels,s=4500

Fig.5 Average measurement error

a—speckle parameterr=3pixels,s=2500 b—speckle parameterr=2pixels,s=4500

Table 1 Comparison of computational performance at different data points

模擬散斑實驗結果表明，隨著應變的增加，3種算法的誤差平均值均呈現上升趨勢，且曲面擬合方法的誤差較大，本文中算法生成較為準確的初值估計，一定程度上降低了平均誤差值，但在應變值小于0.01時誤差值稍大于傳統N-R迭代法，這是仿射變換過程引入了計算誤差的緣故。本文中算法和傳統N-R迭代法的標準差較為穩定，曲面擬合法標準差數值較大且有震蕩現象。應變進一步增大后，本文中算法的標準差值相對傳統N-R迭代法的標準差值降低，說明在大變形情況下本文中算法能夠保持穩定性。同時，對比結果可看出散斑顆粒大小和顆粒數量的變化對變形測量的影響不大。從表1可看出，本文中算法一定程度上也減少了N-R迭代運算次數，對比兩種方法的平均計算時間可發現，本文中算法的匹配時間相對傳統方法平均降低37.54%,在計算效率上有一定的提高。

Fig.6 Standard deviation of measurement error

a—speckle parameterr=3pixels,s=2500 b—speckle parameterr=2pixels,s=4500

3.2 壓縮變形實驗

為驗證本文中算法在實際應用中的可靠性，選取直徑60mm、高度150mm的橡膠圓柱棒材進行壓縮變形測量實驗。在試件表面制作散斑圖案，將試件放置在位移分辨率為0.001mm的精密伺服壓力機上，以0.1mm為步距對試件壓縮，試件壓縮量分別為0.5mm,0.6mm,…,3.5mm。采用Baumer的LXG120M型號CCD相機采集圖像，幀頻為10frame/s。具體實驗裝置如圖7a所示。

為了獲得材料在壓縮后的變形場信息，在一系列壓縮圖像中選取兩幅圖像作為圖像處理的樣本，如圖7b、圖7c所示。采用本文中算法和傳統N-R迭代法對橡膠壓縮變形表面進行區域位移場分布計算，計算過程中數據點之間的間距設置為51個像素，計算結果如圖8所示。結果顯示位移場中矢量箭頭朝向試件的右下方，這是由于所選測量區域位于橡膠中部偏右，軸向的壓縮位移和徑向的膨脹位移同時作用。將位移場沿x和y方向分解，結果顯示y方向為主要變形方向，x方向由左向右位移值逐漸增大，符合橡膠圓柱壓縮過程中軸向變形為主，膨脹變形為輔的規律，等值線分布圖如圖9所示。通過位移場可以計算出應變分布，本文中算法感興趣區域內測量應變均值為0.0061，與采用ANSYS 17.0仿真分析中該區域0.0069應變值接近，說明測量算法在實際應用中可靠。

Fig.7 Experimental device and sample images

a—experimental device b—pre-deformation image c—post deformation image

Fig.8 Displacement field distribution

Fig.9 Displacement field distribution of x and y direction

a—xdirection displacement field distribution b—ydirection displacement field distribution

4 結 論

提出一種基于遺傳算法的數字圖像相關法，結合遺傳算法與粗細搜索混合算法的全局搜索性能準確匹配出整像素位移,后續通過仿射變換原理根據待測點周圍若干變形前后的估值點對推算出待測點更加準確的位移估值，以此代入N-R迭代算法，加快其收斂速度和精度。模擬實驗表明，本文中采取的混合算法使得各算法優勢優勢互補，對不同的散斑圖像均有適應性，且相較于傳統DIC方法，在大應變測量中能夠有效提高測量精度且保證測量結果穩定性。橡膠圓柱壓縮實驗結果驗證了本文中算法在實際大變形測量中的可行性。后期將使用該方法進行金屬材料的大變形測量和振動抗干擾分析。