激光散斑的亞像素位移法計(jì)算及比較

王希軍

(中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林長(zhǎng)春130033)

1 引言

數(shù)字散斑相關(guān)方法(Digital Speckle Correlation Method，DSCM)是一種有效的光學(xué)測(cè)量手段。與傳統(tǒng)的干涉測(cè)量方法不同，它直接從物體表面隨機(jī)分布的人工或自然散斑場(chǎng)中提取位移和變形信息，不需要復(fù)雜的干涉條紋處理過(guò)程，具有全場(chǎng)測(cè)量、非接觸等優(yōu)點(diǎn)，且其光路簡(jiǎn)單，對(duì)測(cè)量環(huán)境要求不高，因而被廣泛應(yīng)用于實(shí)驗(yàn)力學(xué)及其他科學(xué)研究領(lǐng)域［1-4］。

在數(shù)字散斑相關(guān)測(cè)量中，獲取的數(shù)字散斑圖是離散的灰度信息，經(jīng)相關(guān)運(yùn)算后，得到整像素為單位的位移值，而物體的真實(shí)位移值不一定恰好是整像素，因此，亞像素相關(guān)測(cè)量方法已成為數(shù)字散斑相關(guān)測(cè)量中的關(guān)鍵技術(shù)和研究重點(diǎn)。

目前，國(guó)內(nèi)外常用的亞像素相關(guān)測(cè)量方法有插值法［5-6］、曲面擬合法［7］、牛頓-拉夫森法［8］、梯度算法［9-10］、頻域相關(guān)算法［10-11］以及遺傳算法［12］等。雖然對(duì)這些計(jì)算算法的精度分析研究開(kāi)展得較早，但很少有利用算法分析和計(jì)算激光散斑圖并比較各種算法精度的工作，因此，根據(jù)不同的測(cè)量要求選擇合適的亞像素算法對(duì)于提高激光散斑相關(guān)測(cè)量技術(shù)的精度和計(jì)算速度等工作顯得尤為重要。

本文逐一比較了插值算法、梯度算法和牛頓-拉夫森迭代算法的精度和效率。分析了處理激光散斑圖時(shí)這些算法的特點(diǎn)及適用性，為滿(mǎn)足激光散斑相關(guān)測(cè)量數(shù)據(jù)處理和計(jì)算精度的特定要求，提供選擇亞像素計(jì)算方法的依據(jù)。

2 亞像素計(jì)算散斑位移圖

散斑場(chǎng)中各點(diǎn)周?chē)泥徑鼌^(qū)域通常稱(chēng)為子集或窗口。選取散斑場(chǎng)中某一點(diǎn)為中心的窗口作為該點(diǎn)位移和形變的信息，分析和搜索該窗口的移動(dòng)和變化，可獲得該點(diǎn)的位移和形變信息。在進(jìn)行數(shù)字散斑相關(guān)測(cè)量過(guò)程中，首先在物體變形前后拍攝兩幅散斑圖，如圖1所示。在變形前的散斑圖(參考圖)中以待測(cè)點(diǎn)P為中心取m×m的子區(qū)A，當(dāng)被測(cè)物體發(fā)生形變后，即在物體變形后所拍的散斑圖中(目標(biāo)圖)，子區(qū)A移至子區(qū)B的位置，相應(yīng)P點(diǎn)移動(dòng)到P'。由統(tǒng)計(jì)學(xué)可知，A與B兩個(gè)樣本空間的相關(guān)系數(shù)最大。因此可以根據(jù)相關(guān)系數(shù)的最大值來(lái)確定子區(qū)B的位置，從而得到P 點(diǎn)的位移［(x'－x)，(y'－y)］。

圖1 數(shù)字散斑相關(guān)方法測(cè)量原理示意圖Fig.1 Schematic diagram of digital speckle correlation method

但通過(guò)上述方法只能得到整像素位移，而物體實(shí)際的位移值往往不是整像素，并且整像素的位移定位精度是遠(yuǎn)遠(yuǎn)不夠的，需要進(jìn)行亞像素相關(guān)測(cè)量。常用的亞像素相關(guān)方法有插值法、牛頓-拉夫森法及梯度法。

2.1 相關(guān)系數(shù)插值法

該方法先用十字搜索法或其它整像素搜索方法算出子區(qū)的整像素位移，然后對(duì)整像素位移點(diǎn)相鄰9×9點(diǎn)的相關(guān)系數(shù)矩陣進(jìn)行二維插值。插值函數(shù)可以為雙線性插值、雙三次樣條插值函數(shù)等，一般雙三次樣條插值和雙線性插值具有較高的精度［13］。相關(guān)系數(shù)插值法的算法比較容易實(shí)現(xiàn)，但計(jì)算量較大。本文主要討論采用雙三次樣條插值計(jì)算激光散斑圖的結(jié)果。

2.2 梯度算法

基于梯度的亞像素位移算法最初由Davis和Freeman［9］提出，該方法假設(shè)物體只發(fā)生剛體運(yùn)動(dòng)。當(dāng)物體發(fā)生微小變形且選擇的計(jì)算子區(qū)足夠小時(shí)，則該子區(qū)即可近似剛體運(yùn)動(dòng)。令f(x，y)，g(x*，y*)分別表示物體變形前后子區(qū)的灰度，則:

式中:u，v分別為子區(qū)在x，y方向上的整像素位移，Δx，Δy分別為子區(qū)在x，y方向上的亞像素位移。選取最小平方距離相關(guān)函數(shù)作為子區(qū)匹配的算子，則:

圖2 鄰域相關(guān)系數(shù)矩陣差值示意圖Fig.2 Schematic diagram of neighborhood correlation coefficient matrix interpolation

真實(shí)的亞像素位移Δx，Δy應(yīng)使式(2)取得最小值，因此由極值條件:

式中，gx，gy為灰度的一階導(dǎo)數(shù)。

2.3 牛頓-拉夫森迭代算法

牛頓-拉夫森迭代法最初由 Bruck［8］提出，由于在算法中引入了位移的一次導(dǎo)數(shù)，因此其可以對(duì)具有形變的圖像子區(qū)進(jìn)行計(jì)算。

令f(x，y)為樣品變形前拍攝的散斑圖的灰度分布，g(x'，y')為變形后對(duì)拍攝的散斑圖的灰度分布。對(duì)于變形前散斑圖的某一子區(qū)灰度分布fs(x，y)，設(shè)其在變形后對(duì)應(yīng)于灰度分布 gs(x'，y')，則應(yīng)有 gs(x'，y')=fs(x，y)，且應(yīng)有以下對(duì)應(yīng)關(guān)系:

式中:u為變形后子區(qū)中心(x0，y0)在x方向上的位移，v為變形后子區(qū)中心(x0，y0)在y方向上的位移，Δx，Δy為點(diǎn)(x，y)到子區(qū)中心(x0，y0)的距離，?u/?x，?u/?y，?v/?x，?v/?y 為位移的一階偏導(dǎo)數(shù)。選取最小平方距離相關(guān)函數(shù):

由式(6)可見(jiàn)，子區(qū) f(x，y)與 g(x'，y')的相關(guān)系數(shù)C是關(guān)于子區(qū)位移和形變參數(shù)的函數(shù)，令P=(u，v，ux，uy，vx，vy)T，即 C=C(P)。顯而易見(jiàn)，當(dāng)P取得子區(qū)真實(shí)的位移和形變參量時(shí)，相關(guān)系數(shù)應(yīng)為最小值。因此可以通過(guò)式(6)的極值條件:

求出子區(qū)的位移和形變參數(shù)。式(7)可以用迭代法求解。令P0為迭代初值，P'為式(7)的解，對(duì)相關(guān)系數(shù)C(P')在P0點(diǎn)做泰勒級(jí)數(shù)展開(kāi)，得:

牛頓-拉夫森迭代法需要用其他算法獲得整像素位移值作為迭代初值，否則將嚴(yán)重影響收斂速度。

3 計(jì)算結(jié)果及討論

采用數(shù)值模擬的仿真散斑圖對(duì)上述3種算法的精度進(jìn)行了比較。激光散斑圖來(lái)自于磁流體干燥過(guò)程實(shí)驗(yàn)，如圖 3所示［14］。

圖3 磁流體干燥實(shí)驗(yàn)中的激光散斑圖Fig.3 Typical laser speckle diagram of the magnet fluid drying

在相關(guān)搜索算法中需要首先選擇模板窗口，選取大模板窗口可以降低噪聲的影響，提高計(jì)算精度，但大的窗口意味著要計(jì)算更多的點(diǎn)，因此計(jì)算效率會(huì)降低。

Zhang在文獻(xiàn)［15］中分析了窗口大小對(duì)計(jì)算精度和計(jì)算效率的影響，他指出窗口選取在31 pixel×31 pixel到51 pixel×51 pixel之間可以在精度和效率之間獲得較好的平衡。本文選取窗口為 41 pixel×41 pixel。

類(lèi)似文獻(xiàn)［9］方法，截圖以512×512為參考，不失一般性地模擬出平移0.01～0.09 pixel(行方向上)，間隔0.01 pixel的 9幅散斑圖和平移0.1～0.9 pixel(列方向)，間隔為0.1 pixel的9幅位移散斑圖，共模擬出1幅基準(zhǔn)散斑圖和18幅位移散斑圖。

然后，分別用3種算法計(jì)算18幅相對(duì)于基準(zhǔn)散斑圖的位移散斑圖，得到了3種亞像素算法偏差均值圖如圖4所示和標(biāo)準(zhǔn)差圖如圖5所示。

從圖4中可以看出，在這3種算法中梯度算法的偏差均值在整個(gè)預(yù)設(shè)位移范圍0.01～0.9 pixels內(nèi)絕對(duì)值都是最小的，在±0.001 pixels之內(nèi);牛頓-拉夫森迭代算法在位移為0.01～0.1 pixels時(shí)較插值算法的均值誤差小，在位移為0.1～0.9 pixels時(shí)，牛頓-拉夫森迭代算法和插值算法的平均誤差相當(dāng)，達(dá)到±0.003 pixels。但就處理激光散斑圖的算法而言，梯度算法處理精度最好，插值法的處理精度最差，牛頓-拉夫森迭代算法的處理精度居中。標(biāo)準(zhǔn)差表示數(shù)據(jù)的誤差分布，標(biāo)準(zhǔn)值越小，誤差分布越小，故可用它來(lái)評(píng)價(jià)算法的穩(wěn)定性。從圖5標(biāo)準(zhǔn)差曲線可以看出，梯度算法在3種算法中標(biāo)準(zhǔn)差值最小，但在0.5 pixel位移圖計(jì)算時(shí)是個(gè)例處，插值算法的標(biāo)準(zhǔn)差最大。因此，在3種算法中，梯度算法最穩(wěn)定，而插值算法最不穩(wěn)定。至于梯度算法在位移0.5 pixels處產(chǎn)生了較大的標(biāo)準(zhǔn)差，潘兵在文獻(xiàn)［16］給出的解釋:式(3)在做泰勒展開(kāi)時(shí)的截?cái)嗾`差為o(h)2，這是位移為0.5 pixel時(shí)達(dá)到最大值所致。

圖4 3種算法的均值誤差Fig.4 Means bias error of three methods

圖5 3種算法的標(biāo)準(zhǔn)差Fig.5 Standard deviation of three methods

表1 3種算法的計(jì)算時(shí)間(以插值算法的用時(shí)為歸一單位)Tab.1 Computing consuption of three methods

表1給出了3種算法的相對(duì)時(shí)間比較，3種算法是用Matlab編程實(shí)現(xiàn)的，計(jì)算機(jī)的CPU為Intel Pentium D，主頻2.8 GHz，內(nèi)存2 G。通過(guò)比較中可以看出，梯度算法和插值算法的用時(shí)相當(dāng)，略高于插值算法;而牛頓-拉夫森迭代算法用時(shí)最長(zhǎng)，約是插值算法的9倍。這是由于牛頓-拉夫森算法中每次迭代都要對(duì)相關(guān)系數(shù)矩陣進(jìn)行插值運(yùn)算，而在插值算法和基于梯度的算法中只進(jìn)行了一次插值運(yùn)算。綜合比較3種算法的精度、穩(wěn)定性和用時(shí)，可以看出，在處理激光散斑圖時(shí)，梯度算法優(yōu)于其它兩種算法。

4 結(jié)論

本文介紹了數(shù)字散斑圖像相關(guān)法的位移測(cè)量原理以及較常用的插值算法、梯度算法和牛頓-拉夫森迭代算法3種亞像素位移算法及特點(diǎn)。同時(shí)利用算法得到了激光散斑圖計(jì)算的偏差均值曲線和標(biāo)準(zhǔn)差曲線。利用偏差均值的物理意義，比較出梯度算法在亞像素計(jì)算散斑圖時(shí)的精度較高。通過(guò)標(biāo)準(zhǔn)差曲線得出梯度算法在亞像素計(jì)算散斑圖時(shí)的穩(wěn)定性較好。盡管梯度算法用時(shí)略高于插值法，但對(duì)于激光散斑測(cè)量小位移實(shí)驗(yàn)或連續(xù)剛性位移實(shí)驗(yàn)，梯度算法比較有效，應(yīng)作為優(yōu)先考慮的算法。

當(dāng)像素位移量高達(dá)0.5時(shí)，梯度算法產(chǎn)生了較大誤差。此時(shí)牛頓-拉夫森迭代法和插值法在0.5～0.9 pixels區(qū)間的偏差均值較小，對(duì)于激光散斑測(cè)量大位移或振動(dòng)實(shí)驗(yàn)，它們的標(biāo)準(zhǔn)差在0.003 pixels之內(nèi)。

［1］ PETERS W H，RANSON W F.Digital imaging techniques in experimental stress analysis［J］.Opt Eng.，1981，21:427-431.

［2］ 邢冬梅，李鴻琦，李林安，等.梯度功能材料穩(wěn)態(tài)溫度場(chǎng)下的熱應(yīng)力分析［J］.功能材料，2002，33(5):505-507.XING D M，LI H Q，LI L A，et al..Study on thermal strain field of functional gradient material in steady temperature［J］.Functional Material，2002，33(5):505-507.(in Chinese)

［3］ 云禮寧，馬少鵬，李霞鎮(zhèn)，等.用數(shù)字散斑相關(guān)方法研究竹材在拉伸載荷下的斷裂行為［J］.北京理工大學(xué)學(xué)報(bào)，2011，31(3):258-261.YUN L N，MA SH P，LI X ZH，et al..Experimental study of the fraxturc behavior of bamboo under tension based on digital speckle correlation method［J］.Transactions of Beijing Institute of Technology，2011，31(3):258-261.(in Chinese)

［4］ CHU T C，RANSON W F，SUTTON M A，et al..Applications of digital image-correlation techniques to experimental mechanics［J］.Exp.Mech.，1985，25:232-244.

［5］ PEARSON J J，HINES D C，GOLOSMAN S.Video-rate image correlation processor［J］.SPIE，1977，119:197-205.

［6］ DVORCHENKO V N.Bounds on(deterministic)correlation functions with applications to registration［J］.IEEE Trans Pattern Anal Mach Intell，1983，5(2):206-213.

［7］ 潘兵，續(xù)伯欽.數(shù)字圖像相關(guān)中亞像素位移測(cè)量的曲面擬合法［J］.計(jì)量學(xué)報(bào)，2005，26(2):128-134.PAN B，XU B Q.Subpixel registration using qradratic surface fitting in digital image correlation［J］.Acta Metrologica Sinica，2005，26(2):128-134.(in Chinese)

［8］ BRUCK H A，MCNEILL S R，SUTTON M A，et al..Digital image correlation using Newton Raphson method of partial di6erential correction［J］.Exp Mech，1989，29:261-267.

［9］ AVIS C Q，F(xiàn)REEMAN D M.Statistics of subpixel registration algorithms based on spatiotemporal gradients or block matching［J］.Opt Eng.，1998，37(4):1290-1298.

［10］ ZHOU P，GOODSON K E.Subpixel displacement and deformation gradient measurement using digital image/speckle correlation(DISC)［J］.Opt.Eng.，2001，40(8):1613-1620.

［11］ CHEN D J，CHIANG F P，TAN Y S，et al..Digital speckle-displacement measurement using a complex spectrum method［J］.Appl.Opt.，1993，32:1839-1849.

［12］ ORIAT L，LANTZ E.Subpixel detection of the center of an object using a spectral phase algorithm on the image［J］.Pattern Recognition，1998，31:761-771.

［13］ 李慶揚(yáng)，王能超.易大義編著，數(shù)值分析［M］.北京:清華大學(xué)出版社，2008.LI Q Y，WANG N CH，YI D Y.Numerical Analysis［M］.Beijing:Tsinghua University Press，2008.(in Chinese)

［14］ 李新忠，岱欽，王希軍，等.用動(dòng)態(tài)散斑法研究磁流體薄膜的干燥過(guò)程［J］.光電子·激光，2007，18(2):959-962.LI X ZH，DAI Q，WANG X J，et al..Dynamic speckle method for studying the drying process of magnet fluid film［J］.J.Optoelectronics·Laser，2007，18(2):959-962.(in Chinese)

［15］ ZHANG J，JIN G CH，MA SH P，et al..Application of an improved subpixel registration algorithm on digital speckle correlation measurement［J］.Optics Laser Technology，2003，35:533-542.

［16］ 潘兵，謝惠民，戴福隆.數(shù)字圖像相關(guān)中亞像素位移測(cè)量算法的研究［J］.力學(xué)學(xué)報(bào)，2007，39(2):245-252.PAN B，XIE H M，DAI F L.An investigation of sub-pixel displacements registration algorithms in digital image correlation［J］.Chinese J.Theoretical and Appl.Mechanics，2007，39(2):245-252.(in Chinese)