基于位移場(chǎng)局部最小二乘擬合的綺光軸運(yùn)動(dòng)應(yīng)變求解方法

基于位移場(chǎng)局部最小二乘擬合的綺光軸運(yùn)動(dòng)應(yīng)變求解方法

1數(shù)字圖像相關(guān)法材料或結(jié)構(gòu)的表面變形測(cè)量是實(shí)驗(yàn)固體力學(xué)的一項(xiàng)重要任務(wù)之一。由于具有非接觸、全場(chǎng)、高空間分辨力和高精度測(cè)量的優(yōu)點(diǎn),包括干涉方法和非干涉方法在內(nèi)的各種光學(xué)測(cè)量方法在實(shí)驗(yàn)固體力學(xué)中扮演著不可替代的重要角色。I.Yamaguchi、W.H.Peters等提出數(shù)字圖像相關(guān)(Digitalimagecorrelation,DIC)方法已經(jīng)成為實(shí)驗(yàn)固體力學(xué)中一種有眾多應(yīng)用的非接觸全場(chǎng)光學(xué)測(cè)量方法。在過去的20多年中對(duì)該方法的研究多集中在如何提高其位移測(cè)量的精確度上,并提出多種行之有效的亞像素位移測(cè)量算法。如同其他任何一種實(shí)驗(yàn)測(cè)量方法一樣,數(shù)字圖像相關(guān)方法在實(shí)際應(yīng)用中由于各種因素(如照明光強(qiáng)的波動(dòng),圖像采集過程中的各種噪聲及程序計(jì)算中的舍入誤差等)的影響不可能完全精確地恢復(fù)真實(shí)的位移場(chǎng)。各種自編或商業(yè)軟件報(bào)道可達(dá)到的位移測(cè)量精度在±0.02～±0.1pixel之間。如果直接對(duì)位移場(chǎng)進(jìn)行差分來計(jì)算應(yīng)變,則位移場(chǎng)中包含的噪聲將會(huì)被惡性放大,從而得到不可信的應(yīng)變計(jì)算結(jié)果。因此如何從含噪聲的離散位移數(shù)據(jù)中準(zhǔn)確地提取人們更關(guān)注的應(yīng)變信息成為該方法的一個(gè)關(guān)鍵問題。需要說明的是,在數(shù)字圖像相關(guān)方法中利用經(jīng)典的牛頓-拉夫森(Newton-Raphson,N-R)算法可直接給出各計(jì)算點(diǎn)的位移梯度(應(yīng)變)信息,但該應(yīng)變計(jì)算結(jié)果的誤差較大,文獻(xiàn)表明,僅當(dāng)被測(cè)物體的應(yīng)變大于1%時(shí)該方法的應(yīng)變計(jì)算結(jié)果才可靠。因此通過對(duì)原始的位移場(chǎng)信息平滑去噪后再通過數(shù)值差分計(jì)算會(huì)獲得比牛頓-拉夫森算法直接給出應(yīng)變信息更可靠的應(yīng)變計(jì)算結(jié)果。但利用有限元平滑或薄板樣條對(duì)含噪聲離散位移數(shù)據(jù)平滑的方法,其數(shù)學(xué)形式和編程執(zhí)行都顯得異常復(fù)雜和繁瑣,不適合于一般應(yīng)用和推廣。因此本文提出對(duì)數(shù)字圖像相關(guān)方法直接計(jì)算得到的含噪聲的離散位移場(chǎng)進(jìn)行逐點(diǎn)局部最小二乘擬合來計(jì)算應(yīng)變場(chǎng)。首先介紹了數(shù)字圖像相關(guān)方法的原理,闡述了基于局部最小二乘擬合的全場(chǎng)應(yīng)變求解方法,并討論了計(jì)算區(qū)域邊界、孔及裂紋附近區(qū)域等情況下的應(yīng)變計(jì)算方法。最后對(duì)模擬均勻應(yīng)變和中心帶圓孔的鋁板單向拉伸實(shí)驗(yàn)圖片進(jìn)行了計(jì)算。2圖像數(shù)字圖像的相關(guān)方法2.1基于局部位移場(chǎng)的數(shù)字圖像相關(guān)法數(shù)字圖像相關(guān)方法通過處理變形前后被測(cè)對(duì)象表面的數(shù)字圖像直接獲得位移和應(yīng)變信息。通常將變形前的數(shù)字圖像稱為“參考圖像”(Referenceimage),變形后的數(shù)字圖像稱為“變形后圖像”(Deformedimage)。在參考圖像中取以某待求點(diǎn)(x,y)為中心的(2M+1)pixel×(2M+1)pixel(M為局部位移場(chǎng)中各數(shù)據(jù)點(diǎn)的局部坐標(biāo))大小的矩形參考圖像子區(qū),在變形后圖像中通過一定的搜索方法按預(yù)先定義的互相關(guān)函數(shù)進(jìn)行相關(guān)計(jì)算,尋找與參考圖像子區(qū)的相關(guān)系數(shù)為最大值或最小值(取決于所選擇的相關(guān)函數(shù))的以(x′,y′)為中心的目標(biāo)圖像子區(qū),以確定參考圖像子區(qū)的位移u,v。在利用數(shù)字圖像相關(guān)方法進(jìn)行實(shí)際計(jì)算時(shí)通常將參考圖像中間的待計(jì)算區(qū)域劃分成虛擬網(wǎng)格形式,通過計(jì)算每個(gè)網(wǎng)格節(jié)點(diǎn)的位移得到全場(chǎng)位移信息。為提高位移測(cè)量的空間分辨率,網(wǎng)格節(jié)點(diǎn)間距離(計(jì)算步長(zhǎng))通常取在2～10pixel之間。2.2圖像正演分析如圖1所示,假設(shè)參考圖像子區(qū)中的內(nèi)各點(diǎn)(x,y)與變形后的目標(biāo)圖像子區(qū)中各點(diǎn)(x′,y′)按(1)式的函數(shù)關(guān)系一一對(duì)應(yīng):x′=x0+Δx+u+?u?xΔx+?u?yΔy,y′=y0+Δy+v+?v?xΔx+?v?yΔy,(1)式中u,v分別為參考圖像子區(qū)中心在x,y方向的位移值,(Δx,Δy)為點(diǎn)(x,y)到計(jì)算窗口中心(x0,y0)的距離,ux,uy,vx,vy為圖像子區(qū)的位移梯度。2.3圖像平均生成比例相關(guān)函數(shù)是評(píng)價(jià)參考圖像子區(qū)f(x,y)和目標(biāo)圖像子區(qū)g(x′,y′)相似程度的函數(shù),有近十種相關(guān)函數(shù)可供選擇。本文使用歸一化最小平方距離相關(guān)函數(shù):Cf,g(p)=Μ∑x=-ΜΜ∑y=-Μ[f(x,y)-fm√Μ∑x=-ΜΜ∑y=-Μ[f(x,y)-fm]2-g(x′,y′)-gm√Μ∑x=-ΜΜ∑y=-Μ[g(x′,y′)-gm]2]2,(2)式中p為迭代后的近似值,fm=1(2Μ+1)2Μ∑x=-ΜΜ∑y=-Μ[f(x,y)]2,gm=1(2Μ+1)2Μ∑x=-ΜΜ∑y=-Μ[g(x′,y′)]2,分別為參考圖像子區(qū)和目標(biāo)圖像子區(qū)的灰度平均值。該相關(guān)函數(shù)實(shí)際上可從最常用的標(biāo)準(zhǔn)協(xié)方差互相關(guān)函數(shù)推導(dǎo)得到。與其他互相關(guān)函數(shù)相比,該互相關(guān)函數(shù)的相關(guān)峰是全場(chǎng)唯一且尖銳,因而能更準(zhǔn)確地尋找到整個(gè)搜索區(qū)域的相關(guān)系數(shù)極值。并且該相關(guān)函數(shù)的抗干擾性強(qiáng)于其他相關(guān)函數(shù),它對(duì)目標(biāo)圖像子區(qū)灰度的線性變換不敏感,如對(duì)目標(biāo)子區(qū)的灰度做線性變換g′(x′,y′)=a×g(x′,y′)+b,a,b均為常量系數(shù),則按(2)式計(jì)算的相關(guān)系數(shù)值維持不變。2.4cp的初始估計(jì)從(2)式可知,相關(guān)函數(shù)Cf,g(p)是關(guān)于p=(u,ux,uy,v,vx,vy)T6個(gè)參量的函數(shù),其取值范圍是[0,4]。當(dāng)變形前后的圖像子區(qū)最為相似時(shí),相關(guān)系數(shù)Cf,g(p)應(yīng)取最小值,亦即相關(guān)系數(shù)Cf,g(p)的梯度趨近于0,因此有?Cf,g(p)=(?C?pi)i=1,2,?,6=-2Μ∑x=-ΜΜ∑y=-Μ{[f(x,y)-fm√Μ∑x=-ΜΜ∑y=-Μ[f(x,y)-fm]2-g(x′,y′)-gm√Μ∑x=-ΜΜ∑y=-Μ[g(x′,y′)-gm]2]×1√Μ∑x=-ΜΜ∑y=-Μ[g(x′,y′)-gm]2??g(x′,y′)?pi}i=1,2,?,6=0.(3)(3)式可用牛頓-拉夫森迭代法或其他優(yōu)化方法求解,并以整像素位移搜索結(jié)果為迭代的初始值。經(jīng)整理可得?C(p)=?C(p0)+??C(p0)(p-p0)=0,(4)因此p=p0-?C(p0)??C(p0).(5)式中p0為初值估計(jì),?C(p0)為相關(guān)函數(shù)的梯度,??C(p0)為相關(guān)函數(shù)的二次偏導(dǎo),通常被稱為赫斯(Hessian)矩陣。根據(jù)G.Vendroux等的研究,可以對(duì)赫斯矩陣做近似處理,在不影響精度的情況下簡(jiǎn)化計(jì)算,簡(jiǎn)化后的赫斯矩陣為??Cf,g(p)=(?2C?pi?pj)i=1,2,?,6j=1,2,?,6≈2Μ∑x=-ΜΜ∑y=-Μ[g(x′,y′)-gm]2Μ∑x=-ΜΜ∑y=-Μ[?2g(x′,y′)?pi?pj]i=1,2,?,6j=1,2,?,6.(6)2.5自然邊界條件下亞像素位置的計(jì)算牛頓-拉夫森方法在實(shí)現(xiàn)過程中需要亞像素位置的灰度和灰度梯度值,所以插值方法的選擇是該方法的關(guān)鍵,因其直接影響計(jì)算的精度和收斂特性。本文采用在自然邊界條件下雙三次樣條插值來計(jì)算亞像素位置的灰度和灰度偏導(dǎo),具體是通過兩個(gè)方向的一維樣條插值得到。程序運(yùn)行表明,該插值算法具有很高的精度和很好的收斂性,與文獻(xiàn)中的結(jié)論相符合。2.6組合條件牛頓-拉夫森方法的收斂條件為相鄰兩次迭代結(jié)果的位移波動(dòng)小于10-4pixel,應(yīng)變波動(dòng)小于5×10-6。3維離散位移擬合在數(shù)學(xué)理論中應(yīng)變和位移的關(guān)系通常被描述為一種數(shù)值差分操作過程,因此為得到應(yīng)變信息,首先對(duì)離散位移數(shù)據(jù)進(jìn)行數(shù)值差分。但數(shù)值差分是一個(gè)不穩(wěn)定和危險(xiǎn)的操作,在應(yīng)用時(shí)需要特別小心,因?yàn)樗鼘?duì)信號(hào)的噪聲(或誤差)具有明顯的放大作用(特別是高頻噪聲)。如果直接對(duì)含噪聲的位移場(chǎng)進(jìn)行差分計(jì)算,那么微小的位移測(cè)量誤差將會(huì)被惡性放大,從而導(dǎo)致更大的應(yīng)變計(jì)算誤差,得到的應(yīng)變信息變得極不可信。例如,假設(shè)位移的測(cè)量誤差為±0.02pixel,相鄰兩個(gè)計(jì)算節(jié)點(diǎn)的間距為5pixel,由直接差分計(jì)算的應(yīng)變誤差最大為Δε=15(|±0.02||±0.02|)=8000×10-6,而由中間差分計(jì)算的應(yīng)變誤差最大為Δε=110(|±0.02||±0.02|)=4000×10-6,如此之大的計(jì)算誤差將完全掩蓋真實(shí)的應(yīng)變信息,在大多數(shù)情況下難以接受。當(dāng)然,可以通過增加標(biāo)距(差分節(jié)點(diǎn)間的距離)來減小應(yīng)變的計(jì)算誤差,但是用直接差分的方法仍然無法計(jì)算區(qū)域邊界、孔洞及裂紋附近區(qū)域的應(yīng)變信息?？紤]到離散位移數(shù)據(jù)中的噪聲影響,本文對(duì)原始的離散位移數(shù)據(jù)用逐點(diǎn)最小二乘擬合的方法來求解位移的導(dǎo)數(shù)(應(yīng)變)。由于擬合的過程既包含了噪聲的去除,因此所得到的應(yīng)變精度較直接差分會(huì)有顯著提高。該方法的基本思想就是用完全二維多項(xiàng)式對(duì)離散位移數(shù)據(jù)的一個(gè)局部子域[稱之為應(yīng)變計(jì)算窗口,假設(shè)其中包含(2M+1)×(2M+1)個(gè)均勻分布的數(shù)據(jù)點(diǎn),如圖2所示]進(jìn)行分片逐點(diǎn)擬合,利用最小二乘法得到擬合多項(xiàng)式的系數(shù),就很容易得到擬合區(qū)域中間點(diǎn)的對(duì)應(yīng)值和各階導(dǎo)數(shù),作為該點(diǎn)的平滑和差分結(jié)果。接著,將局部子域移動(dòng)到下一個(gè)數(shù)據(jù)點(diǎn),計(jì)算新的數(shù)據(jù)子集的中心點(diǎn)平滑和差分結(jié)果。為了更好地濾除離散位移數(shù)據(jù)中的噪聲和得到較好的平滑效果,本文僅考慮擬合函數(shù)為二維一次多項(xiàng)式。具體的做法是,對(duì)位移場(chǎng)局部子域中的離散數(shù)據(jù)用二維一次多項(xiàng)式擬合u、v場(chǎng)位移,則有u(x,y)=a0+a1x+a2y,v(x,y)=b0+b1x+b2y,(7)式中x,y=M,a0,a1,a2,b0,b1,b2為待求的擬合多項(xiàng)式系數(shù),u(x,y),v(x,y)為離散位移數(shù)據(jù)點(diǎn)。假設(shè)局部數(shù)據(jù)子集u(x,y)為由(2M+1)×(2M+1)數(shù)據(jù)點(diǎn)組成,因此可將(2M+1)×(2M+1)數(shù)據(jù)點(diǎn)組成的二維離散位移數(shù)據(jù)矩陣u(x,y)(圖2)排列成一維列向量n=(2M+1)×(2M+1)-1。因此,(7)式的u(x,y)可以重新寫成矩陣形式:u=Xa?[1-m-m1-m+1-m1001m-1m1mm][a0a1a2]=[～u(-m,m)～u(-m+1,m)～u(0,0)～u(m-1,m)～u(m,m)],(8)于是可利用最小二乘法來求解待定系數(shù)向量[a0a1a2]=(XΤX)-1XΤu,(9)式中(XTX)-1XT為X的偽逆矩陣,對(duì)位移場(chǎng)v(x,y)可按同樣方式得到[b0b1b2]=(XΤX)-1XΤv。得到擬合多項(xiàng)式的系數(shù)后,小變形情況下柯西(Cauchy)應(yīng)變分量εx=?u/?x=a1,εy=?v/?y=b2,γxy=?u/?y+?v/?x=a2+b1,(10)有限變形情況下的格林(Green)應(yīng)變分量為Exx=?u?x+12[(?u?x)2+(?v?x)2]=a1+12(a21+b21),Eyy=?v?y+12[(?u?y)2+(?v?y)2]=b2+12(a22+b22),Exy=12(?u?y+?v?x)+12(?u?x?u?y+?v?x?v?y)=12(a2+b1)+12(a1a2+b1b2).(11)對(duì)于邊界點(diǎn)、孔洞及裂紋附近區(qū)域,由于局部位移數(shù)據(jù)所包含的有效數(shù)據(jù)點(diǎn)要小于(2M+1)×(2M+1)個(gè),在計(jì)算系數(shù)矩陣X和位移向量u(x,y)時(shí)可以避免(忽略)這些無效的數(shù)據(jù)點(diǎn)。實(shí)際上由于待求系數(shù)向量只有三個(gè)未知數(shù),因此只要應(yīng)變計(jì)算窗口中的有效數(shù)據(jù)點(diǎn)大于3個(gè)即可。由于通常的應(yīng)變計(jì)算窗口都取在11×11～21×21點(diǎn)之間,因此以上處理方式合理可靠,完全能夠保證(9)式有解。4應(yīng)變分布計(jì)算用具有均勻變形和非均勻變形的兩個(gè)典型實(shí)驗(yàn)驗(yàn)證本文方法的有效性,首先用數(shù)字圖像相關(guān)方法計(jì)算得到的離散位移數(shù)據(jù),隨后用局部最小二乘擬合計(jì)算全場(chǎng)應(yīng)變分布。所用的計(jì)算參量為:圖像子區(qū)大小為41pixel×41pixel,網(wǎng)格節(jié)點(diǎn)步長(zhǎng)為5pixel。4.1x方向的均勻應(yīng)變首先用PengZhou等提出的散斑模擬方法生成具有均勻應(yīng)變的兩幅模擬散斑圖用以檢驗(yàn)本文提出的應(yīng)變計(jì)算方法。模擬散斑圖如圖3所示,具體參量為:圖像大小為576pixel×576pixel,散斑顆粒大小為4pixel,散斑顆粒數(shù)目為4000,信噪比為20。其中變形后圖像相對(duì)參考圖像在x方向有2000×10-6的均勻應(yīng)變。圖4為由牛頓-拉夫森方法計(jì)算得到x方向的位移場(chǎng)u(x),該圖顯示了數(shù)字圖像相關(guān)方法的計(jì)算結(jié)果與預(yù)先施加的隨x坐標(biāo)線性變化的位移場(chǎng)分布十分吻合。對(duì)圖4所示的位移場(chǎng)用局部最小二乘擬合(應(yīng)變計(jì)算窗口為21×21點(diǎn))可得到x方向的應(yīng)變場(chǎng),如圖5所示。由于在均勻應(yīng)變的情況下,應(yīng)變分布在整個(gè)計(jì)算區(qū)域范圍內(nèi)應(yīng)該為一常量(2000×10-6),從圖5可以看出本文方法的計(jì)算結(jié)果是可靠和精確的。圖6為對(duì)圖5的離散位移場(chǎng)利用不同大小的應(yīng)變計(jì)算窗口(從3×3點(diǎn)到21×21點(diǎn))所計(jì)算得到的全場(chǎng)應(yīng)變值的統(tǒng)計(jì)結(jié)果。從圖6可以看出對(duì)于均勻應(yīng)變,盡管不同大小應(yīng)變計(jì)算窗口計(jì)算結(jié)果的均值都接近預(yù)先給定的真實(shí)應(yīng)變值,但是較大應(yīng)變計(jì)算窗口所計(jì)算得到的應(yīng)變的標(biāo)準(zhǔn)差明顯減小,意味著在均勻變形情況下應(yīng)該選擇大的應(yīng)變計(jì)算窗口,其計(jì)算結(jié)果將更逼近真值。4.2應(yīng)變計(jì)算方法的確定對(duì)中心帶孔的薄鋁板(材料為L(zhǎng)Y10)在單向拉伸情況下的全場(chǎng)變形情況進(jìn)行了計(jì)算,試件形狀及幾何尺寸如圖7所示,其中試樣中心有一直徑為1.5mm的通孔。試樣表面的人工散斑和計(jì)算區(qū)域如圖8(a)所示,淺色方框內(nèi)為要計(jì)算的區(qū)域,而淺色方框內(nèi)的橢圓形深色區(qū)域?yàn)橐苊庥?jì)算的無效區(qū)域(圓孔)。試樣加載前后表面的參考圖像和變形后圖像由CCD攝像機(jī)采集(