基于位移場局部最小二乘擬合的綺光軸運動應變求解方法

基于位移場局部最小二乘擬合的綺光軸運動應變求解方法

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