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

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

引 言

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

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

1 DIC的基本原理

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

Fig.1 Subset region of image before and after deformation

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

C=

(1)

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

2 基于遺傳算法的數(shù)字圖像相關(guān)法

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

2.1 遺傳算法描述

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

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

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

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

(2)

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

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

(3)

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

經(jīng)過上述選擇、交叉、變異產(chǎn)生的下一代群體，若最優(yōu)個體的適應(yīng)度值不滿足要求，則重新進(jìn)行個體適應(yīng)度評估并進(jìn)入下一輪循環(huán)，直到最優(yōu)個體滿足要求。遺傳過程的停止條件為兩代種群中的最優(yōu)個體和當(dāng)代種群中的最優(yōu)個體相同或者進(jìn)行到最大迭代次數(shù)。

2.2 遺傳算法與粗細(xì)搜索法的混合算法

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

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

Fig.3 Schematic of hybrid algorithm

2.3 仿射變換初值估計

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

(4)

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

2.4 亞像素位移計算

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

(5)

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

3 實(shí)驗(yàn)驗(yàn)證

3.1 模擬散斑圖仿真實(shí)驗(yàn)

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

I1(x,y)=

(6)

(7)

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

分別利用傳統(tǒng)N-R迭代法、曲面擬合法與本文中算法對目標(biāo)圖像中隨機(jī)數(shù)據(jù)點(diǎn)進(jìn)行匹配。比較不同散斑顆粒尺寸、不同散斑顆粒數(shù)目和不同應(yīng)變值對測量平均誤差和標(biāo)準(zhǔn)差的影響，結(jié)果如圖5、圖6所示。最后對迭代次數(shù)和匹配時間進(jìn)行對比，如表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

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

Fig.6 Standard deviation of measurement error

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

3.2 壓縮變形實(shí)驗(yàn)

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

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

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 結(jié) 論

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