資新運(yùn)， 耿 帥， 趙姝帆， 夏軍劍， 舒 放

（1.軍事交通學(xué)院工程實(shí)驗(yàn)中心，天津300161； 2.天津財(cái)經(jīng)大學(xué)理工學(xué)院，天津300222）

3種亞像素位移測量算法的比較研究

資新運(yùn)1， 耿 帥1， 趙姝帆1， 夏軍劍1， 舒 放2

（1.軍事交通學(xué)院工程實(shí)驗(yàn)中心，天津300161； 2.天津財(cái)經(jīng)大學(xué)理工學(xué)院，天津300222）

梯度法、曲面擬合法和N-R法是數(shù)字散斑相關(guān)方法中提高位移測量精度的3種主要亞像素位移算法。研究了3種算法的原理并進(jìn)行了數(shù)值模擬仿真。首先使用計(jì)算機(jī)模擬生成一系列位移為0.01 pixel的散斑圖，并在生成的圖像中添加方差不同的高斯噪聲，對3種算法在噪聲條件下的精度進(jìn)行比較，結(jié)果顯示在精度要求為0.01 pixel時(shí)，梯度法的噪聲方差上限為0.006，N-R法的噪聲方差上限為0.000 5，而曲面擬合法的計(jì)算精度具有較大不確定性，無法保證精度，梯度法更適用于工程應(yīng)用。

計(jì)量學(xué)；數(shù)字散斑相關(guān)；曲面擬合；梯度法；N-R法；亞像素

1 引 言

數(shù)字散斑相關(guān)方法（Digital Speckle Correlation Method，DSCM）又稱數(shù)字圖像相關(guān)方法（Digital Image Correlation Method，DICM），是上世紀(jì)80年代提出的一種可以測量物體變形和應(yīng)變的光學(xué)非接觸測量方法［1］，它通過分析變形前后物體表面散斑圖像的灰度來獲得被測物體的力學(xué)性能。該方法具有非接觸、全場測量、實(shí)驗(yàn)設(shè)備簡單，對環(huán)境要求低等優(yōu)點(diǎn)［2］，經(jīng)過多年的研究，其在理論和應(yīng)用方面取得了較大的發(fā)展。近年來，基于數(shù)字圖像相關(guān)方法的測量軟件已經(jīng)面世，應(yīng)用領(lǐng)域包括材料力學(xué)、振動(dòng)力學(xué)、沖擊力學(xué)、破壞力學(xué)、微納米力學(xué)、疲勞及可靠性等［3～5］，并可實(shí)時(shí)監(jiān)測零部件應(yīng)力分布情況［6，7］，預(yù)防故障發(fā)生。

數(shù)字散斑相關(guān)技術(shù)通過記錄物體變形前后的圖像并進(jìn)行相關(guān)計(jì)算得到物體的位移和變形。由于該方法處理的是數(shù)字化的圖像，因此獲得的位移只能是像素的整數(shù)倍，但實(shí)際的位移值通常不恰好為整像素，而且整像素的位移精度在固體材料或構(gòu)件表面的變形測量中是遠(yuǎn)遠(yuǎn)不夠的。為獲得更高測量精度，在不提高光學(xué)系統(tǒng)放大倍數(shù)和CCD分辨率的前提下，需要對數(shù)字散斑進(jìn)行亞像素處理。目前，相關(guān)研究人員提出了亞像素灰度插值法、曲面擬合法、梯度法、N-R法（Newton-rapshon method）［8～10］等亞像素搜索算法，其中應(yīng)用最為廣泛的是曲面擬合法、梯度法和N-R法。在無噪聲條件下3種算法的測量精度已有定論［11］，然而在工程測量中，圖像采集和傳輸過程中，電器系統(tǒng)和外界影響都可產(chǎn)生噪聲。由系統(tǒng)外部環(huán)境干擾引起的噪聲可以通過改善外部環(huán)境來消除或減小，由光、電、元器件本身及系統(tǒng)內(nèi)部設(shè)備電路引起的噪聲，有些內(nèi)部噪聲是可以減小的，但有些是不可避免的。噪聲會(huì)使散斑圖像表面產(chǎn)生不同的灰度特征，從而對數(shù)字圖像相關(guān)計(jì)算結(jié)果產(chǎn)生一定影響，因此有必要對噪聲情況下的3種亞像素位移算法進(jìn)行比較研究。

2 數(shù)字散斑相關(guān)方法的基本原理

在外力作用下，物體表面發(fā)生應(yīng)變位移，散斑的灰度值保持不變，因此可以計(jì)算散斑位置的變化，從而求得物體表面的位移。利用數(shù)字圖像相關(guān)方法時(shí)，首先采集試件表面變形前后的兩幅散斑圖，在變形前的圖像中，取以待求點(diǎn)P（x1，y1）為中心的一塊矩形圖像子區(qū)f（x，y），這塊區(qū)域的灰度特征保持不變，因此，根據(jù)灰度信息在變形后的目標(biāo)圖像中尋找相同的矩形區(qū)域g（x′，y′）和它的中心點(diǎn)P（x2，y2），計(jì)算得到兩中心點(diǎn)的位移u，v，如圖1所示。在尋找過程中是通過一定的搜索方法，并按某一相關(guān)函數(shù)來進(jìn)行相關(guān)計(jì)算，求得相關(guān)系數(shù)C（u，v）取最大值的矩形區(qū)域的［4］。

式中：x′＝x＋u，y′＝y(tǒng)＋v，u和v分別為中心點(diǎn)在x方向和y方向上的位移；Π是刻畫f（x，y）和g（x′，y′）在某種程度上相似的函數(shù)。目前相關(guān)函數(shù)有近10種，可以根據(jù)要求選擇合適的相關(guān)函數(shù)。

求相關(guān)函數(shù)最值的過程就是求子區(qū)中心點(diǎn)位移和變形的過程，選擇不同的相關(guān)函數(shù)，兩區(qū)域完全相關(guān)時(shí)相關(guān)系數(shù)取最大或最小。根據(jù)峰值相關(guān)系數(shù)計(jì)算該子區(qū)域在變形后的位置，由此得到該區(qū)塊的變形量。對變形前散斑中的所有子區(qū)域進(jìn)行類似運(yùn)算，就可以得到整個(gè)位移場。

圖1 數(shù)字圖像相關(guān)原理示意圖

3 數(shù)字散斑相關(guān)中的亞像素搜索算法

3.1 曲面擬合法

整像素運(yùn)算后得到整像素位移結(jié)果（u，v），同時(shí)也得到了參考窗口與變形后圖像中目標(biāo)窗口間的相關(guān)系數(shù)值，C（u，v）為全部相關(guān)系數(shù)值中的最大值。曲面擬合法是在整像素目標(biāo)點(diǎn)及其鄰近點(diǎn)的相關(guān)系數(shù)陣的基礎(chǔ)上，對其進(jìn)行曲面擬合，擬合曲面的極值點(diǎn)位置即為所求圖像區(qū)域的中心位置。常用的曲面擬合法有高斯曲面擬合、二次曲面擬合、二維拉格朗日曲面擬合等方法。本文選取二次曲面擬合法來實(shí)現(xiàn)相關(guān)系數(shù)場的模擬，再通過求導(dǎo)的方法獲得模擬函數(shù)的極大值，得到亞像素位移。

以整像素搜索算法得到的相關(guān)系數(shù)峰值點(diǎn)為中心，取3×3 pixels擬合窗口，通過這9點(diǎn)相關(guān)系數(shù)值及坐標(biāo)位置進(jìn)行曲面擬合。二次曲面的擬合公式一般采用二元二次多項(xiàng)式：

式中有6個(gè)待定系數(shù)，A＝［a0，a1，a2，a3，a4，a5］T。以相關(guān)系數(shù)峰值點(diǎn)為中心，取3×3 pixels擬合窗口進(jìn)行擬合，可得到包含9個(gè)式（2）的方程組，因此可利用最小二乘法來求解。

對于擬合函數(shù)C，擬合曲面極值點(diǎn)滿足對x和y的偏導(dǎo)數(shù)為0，即：

由式（3）可知曲面極值點(diǎn)位置為：

3.2 梯度法

基于梯度的亞像素位移算法是假設(shè)圖像子區(qū)做微小位移時(shí)，若子區(qū)足夠小，則圖像子區(qū)運(yùn)動(dòng)可近似為面內(nèi)剛體運(yùn)動(dòng)［11］，此時(shí)有：

式中：f（x，y）和g（x′，y′）分別表示位移前、位移后的子區(qū)圖像；u和v分別為與x、y方向上的整像素位移；Δu，Δv為相應(yīng)的亞像素位移。

式（5）可通過一階的泰勒級(jí)數(shù)展開，同時(shí)舍棄高階小量，可得：

gx和gy為子區(qū)灰度一階梯度，這里選用抗噪聲能力強(qiáng)的標(biāo)準(zhǔn)化協(xié)方差互相關(guān)函數(shù)的平方：

3.3 N-R法

N-R法考慮子區(qū)變形的情況，假設(shè)變形前子區(qū)灰度分布為f（x，y），變形后相應(yīng)分布為g（x′，y′），則存在以下對應(yīng)關(guān)系：

式中：u和v為圖像子區(qū)中心在x和y方向上的位移；Δx和Δy為點(diǎn)（x，y）到子區(qū)中心（x0，y0）的距離；?u／?x，?u／?y，?v／?x，?v／?y為位移一階偏導(dǎo)數(shù)。

為了提高精度，選擇零均值歸一化的最小平方距離相關(guān)函數(shù)：

兩子區(qū)最相似即為相關(guān)系數(shù)C取最小值，令P＝（u，v，?u／?x，?u／?y，?v／?x，?v／?y）T，即公式（21）可以表示為C＝C（P），參數(shù)可以通過對該式求極值獲得：

6個(gè)未知參數(shù)構(gòu)成一個(gè)非線性方程，取迭代的估計(jì)初值為：

其中u0，v0為整像素搜索得到的x和y方向上的位移。

通過N-R法求解：

其中▽▽C（P）為相關(guān)系數(shù)的二次偏導(dǎo)數(shù)，也被稱為Hessian矩陣。Hessian矩陣可近似為［12］：

因此有：

經(jīng)推導(dǎo)，可得迭代關(guān)系為：

其分量形式為：

4 亞像素位移測量算法的數(shù)值模擬驗(yàn)證

4.1 模擬散斑生成方法

為了避免在實(shí)際圖像采集過程中引入誤差，更好地驗(yàn)證算法，可采用計(jì)算機(jī)模擬生成數(shù)字散斑，進(jìn)行仿真實(shí)驗(yàn)驗(yàn)證。文獻(xiàn)［13］對模擬散斑的顆粒大小、顆粒數(shù)目等做了研究，并得到以下結(jié)論：最佳散斑顆粒尺寸在3～6 pixels，最佳計(jì)算窗口為31×31～51×51 pixels。借助模擬數(shù)字散斑，通過精確的控制位移量、自由添加噪聲，可以仿真不同位移、噪聲環(huán)境下的變形實(shí)驗(yàn)。

根據(jù)文獻(xiàn)［13］提供的算法，生成一對位移為（u，v）的模擬散斑，位移分量的表達(dá)式為：

式中：u0，v0為剛體位移量；ux，uy，vx，vy為位移梯度。

假設(shè)物體變形前后的數(shù)字散斑呈現(xiàn)高斯分布，以高斯函數(shù)模擬散斑顆粒光強(qiáng)分布，則變形前后圖像灰度分布函數(shù)分別為I1（x，y）和I2（x，y）。

式中：s為散斑顆粒數(shù)量；d為散斑顆粒尺寸；I0為散斑圖像背景光強(qiáng)，計(jì)算中取值為1；（xk，yk）為散斑顆粒的位置，在程序中可由隨機(jī)分布函數(shù)生成。

根據(jù)上述方法，可生成模擬數(shù)字散斑圖如圖2所示，其中散斑圖像大小為256×256 pixels，散斑顆粒尺寸為3 pixels，顆粒數(shù)為1 000個(gè)。

圖2 數(shù)字散斑圖像及其灰度分布圖

圖3 理想情況的數(shù)字散斑圖平移0～0.1 pixel時(shí)3種亞像素算法測量誤差

4.2 理想條件下亞像素位移模擬驗(yàn)證

為了驗(yàn)證以上3種亞像素測量算法的性能，用計(jì)算機(jī)模擬生成一幅沒有噪聲的數(shù)字散斑圖像作為基準(zhǔn)，然后依次平移0.01 pixel，連續(xù)生成10幅數(shù)字散斑序列圖作為目標(biāo)圖像。計(jì)算點(diǎn)陣為11×11共121個(gè)點(diǎn)，計(jì)算子區(qū)大小均采用41×41。使用本文的3種亞像素處理算法計(jì)算位移，得到位移的均值誤差和標(biāo)準(zhǔn)差如圖3所示。

由圖3可見，在理想情況下，N-R法計(jì)算精度最高，梯度法計(jì)算精度稍劣于N-R方法，曲面擬合法最差。這與文獻(xiàn)［11］中得到的結(jié)論相符。

4.3 噪聲影響下亞像素位移模擬驗(yàn)證

在實(shí)際工程應(yīng)用中，測量過程難免會(huì)受到噪聲的影響，根據(jù)文獻(xiàn)［14］，噪聲主要的來源有：CCD暗電流影響、CCD復(fù)位噪聲、圖像卡噪聲、外部環(huán)境振動(dòng)等。為模擬更加真實(shí)的情況，討論噪聲對3種亞像素算法的影響，分別為基準(zhǔn)圖像添加均值為0，方差σ為0.000 1、0.001、0.01、0.1的高斯噪聲，添加噪聲后的數(shù)字散斑圖像如圖4所示。

以相同的方法，分別以每幅含噪聲的數(shù)字散斑圖為新的基準(zhǔn)圖像，依次平移0.01 pixel，連續(xù)生成10幅數(shù)字散斑序列圖作為目標(biāo)圖像。計(jì)算其均值誤差與標(biāo)準(zhǔn)差結(jié)果如圖5所示。圖5（a）、（b）分別為σ＝0.000 1時(shí)3種算法均值誤差和標(biāo)準(zhǔn)差，圖5（c）、（d）分別為σ＝0.001時(shí)3種算法均值誤差和標(biāo)準(zhǔn)差，圖5（e）、（f）分別為σ＝0.01時(shí)3種算法均值誤差和標(biāo)準(zhǔn)差。

由圖5（a）中可以看到，當(dāng)σ＝0.000 1時(shí)梯度法和N-R法的測量精度均能達(dá)到0.01 pixel，曲面擬合法已達(dá)不到精度要求；由圖5（c）中可以看到，當(dāng)σ＝0.001時(shí)梯度法的測量精度能達(dá)到0.01 pixel，N-R法也達(dá)不到精度要求；由圖5（e）中可以看到，當(dāng)σ＝0.01時(shí)梯度法也達(dá)不到精度要求。3種亞像素算法的均值誤差和標(biāo)準(zhǔn)差均呈隨著噪聲增大而增加的趨勢，其中N-R對噪聲最敏感，而梯度法受噪聲影響最小，曲面擬合法介于兩者之間。當(dāng)噪聲方差在0.01時(shí)，3種算法計(jì)算誤差均大于了0.01 pixel。對于精度為0.01 pixel的測量而言，3種算法都無法實(shí)現(xiàn)，繼續(xù)增加噪聲方差已無意義，本文不再計(jì)算在噪聲方差為0.1時(shí)3種算法的測量精度。

從圖5（a）和（c）中可以大致判斷，當(dāng)測量精度要求為0.01時(shí)，N-R法對應(yīng)的噪聲方差上限σ應(yīng)在0.000 1～0.001之間。為了確定N-R法的噪聲方差上限值，現(xiàn)將一張模擬散斑圖依次加入均值為0，方差大小為0.000 1～0.001的高斯噪聲，每次增量為0.000 1，得到10加入噪聲后的圖像作為目標(biāo)圖像，原始圖像為參考圖像。按照同樣的方法計(jì)算均值誤差和標(biāo)準(zhǔn)差。

由圖6（a）中可以看到，N-R法在x方向上要滿足0.01 pixel的測量精度，噪聲方差應(yīng)控制在0.000 6以下，由圖6（c）中可以看到，在y方向上要滿足0.01 pixel的測量精度，噪聲方差應(yīng)控制在0.000 5以下。因此x、y方向要同時(shí)滿足測量精度，噪聲方差應(yīng)控制在0.000 5以下。

從圖5（c）和（e）中可以大致判斷，梯度法對應(yīng)的噪聲方差上限σ應(yīng)在0.001～0.01之間，為了確定梯度法的噪聲方差上限值，運(yùn)用同樣的方法在一張模擬散斑圖依次加入均值為0，方差大小為0.001～0.01的高斯噪聲，每次增量為0.001。

由圖7（a）中可以看到，梯度法在x方向上要滿足0.01 pixel的測量精度，噪聲方差應(yīng)控制在0.007以下，由圖7（c）中可以看到，在y方向上要滿足0.01 pixel的測量精度，噪聲方差應(yīng)控制在0.006以下。因此x、y方向要同時(shí)滿足測量精度，噪聲方差應(yīng)控制在0.006以下。

圖5 噪聲影響的數(shù)字散斑圖平移0～0.1 pixel時(shí)3種亞像素算法測量誤差

圖6 噪聲對N-R法計(jì)算精度的影響

圖7 噪聲對梯度法計(jì)算精度的影響

5 結(jié) 論

在無噪聲條件下，N-R法計(jì)算精度最高，梯度法次之，曲面擬合法最差。在精度要求為0.01 pixel時(shí)，梯度法的噪聲方差上限為0.006，N-R法的噪聲方差上限為0.000 5，而曲面擬合法的計(jì)算精度具有較大不確定性，無法保證精度。3種方法的均值偏差都隨噪聲增大而增大，其中梯度法的均值偏差最小，曲面擬合法次之，N-R的均值誤差最大；3種方法的標(biāo)準(zhǔn)差隨著噪聲的增大而增加，N-R法的標(biāo)準(zhǔn)差增加的幅度最大。由此可見，梯度法的抗噪能力最強(qiáng)，曲面擬合法次之，N-R法抗噪能力最弱。由此可見，N-R法在無噪聲情況下精度很高，但是計(jì)算速度慢并且對噪聲非常敏感，所以N-R法適用于實(shí)驗(yàn)條件下采集的圖像質(zhì)量較高的情況；曲面擬合法的計(jì)算速度快但精度低，因此適用于計(jì)算速度要求高、精度要求一般的情況；梯度法的抗噪聲能力較好，精度較高，計(jì)算速度在兩者之間，可廣泛運(yùn)用于工程實(shí)踐。

］

［1］ PetersW H，Ranson W F.Digital imaging techniques in experimental stress analysis［J］.Optical Engineering，1982，21（3）：427-431.

［2］ 潘兵，謝惠民，續(xù)伯欽，等.數(shù)字圖像相關(guān)中的亞像素位移定位算法進(jìn)展［J］.力學(xué)進(jìn)展，2005，35（3）：345-352.

［3］ 宋義敏，馬少鵬，楊小彬，等.巖石變形破壞的數(shù)字散斑相關(guān)方法研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2011，30（1）：170-175.

［4］ Dong Y，Kakisawa H，Kagawa Y.Optical system for microscopic observation and strain measurement at high temperature［J］.Measurement Science and Technology，2014，25（2）：025002.

［5］ 俞海，郭榮鑫，夏海廷，等.數(shù)字圖像相關(guān)法在WC／Cu復(fù)合材料線膨脹系數(shù)測量中的應(yīng)用［J］.光學(xué)精密工程，2013，21（10）：2696-2703.

［6］ Cofaru C，PhilipsW，Van Paepegem W.A novel speckle pattern—adaptive digital image correlation approach with robust strain calculation［J］.Optics and Lasers in Engineering，2012，50（2）：187-198.

［7］ Sicard J，Sirohi J.Measurement of the deformation of an extremely flexible rotor blade using digital image correlation［J］.Measurement Science and Technology，2013，24（6）：065203.

［8］ 潘兵，續(xù)伯欽，陳丁，等.數(shù)字圖像相關(guān)中亞像素位移測量的曲面擬合法［J］.計(jì)量學(xué)報(bào)，2005，26（2）：128-134.

［9］ 潘兵，吳大方，謝惠民，等.基于梯度的數(shù)字體圖像相關(guān)方法測量物體內(nèi)部變形［J］光學(xué)學(xué)報(bào)，2011，31（6）：0612005

［10］ Cofaru C，Philips W，Van Paepegem W.Improved Newton-Raphson digital image correlation method for full-field displacement and strain calculation［J］.Applied Optics，2010，49（33）：6472-6484.

［11］ 潘兵，謝惠民，戴福隆.數(shù)字圖像相關(guān)中亞像素位移測量算法的研究［J］.力學(xué)學(xué)報(bào)，2007，39（2）：245-252.

［12］ Vendroux G，KnaussW G.Submicron deformation field measurements：Part 2.Improved digital image correlation［J］.ExperimentalMechanics.1998，38（2）：86-92.

［13］ Zhou P，Goodson K E.Sub-pixel displacement and deformation gradient measurement using digital image／speckle correlation［J］.Optical Engineering，2001，40（8）：1613-1620.

［14］ 楊勇，王琰蕾，李明.高精度數(shù)字圖像相關(guān)測量系統(tǒng)及其技術(shù)研究［J］.光學(xué)學(xué)報(bào)，2006，26（2）：197-201.

Com parison of Three Sub-pixel Displacement Algorithm

ZIXin-yun1， GENG Shuai1， ZHAO Shu-fan1， XIA Jun-jian1， SHU Fang2
（1.Engineering Experiment Center，Military Transportation University，Tianjin 300161，China；2.Institute of Technology，Tianjin University of Finance，Tianjin 300222，China）

Gradientmethod，surface fitting and N-R method are 3 main sub-pixel displacement algorithms of digital speckle correlationmethod to improve the precision of themeasuring system.The principles of three kinds of algorithm were studied and numerical simulation was carried out.To compare the different precision from 3 methods under noise，the computermodels generate a series of speckle gram whose displacement is 0.01 pixels and different Gaussian noise was added in generated graphics.The results showed that in the precision requirement in 0.01 pixels，the upper bounds on noise variance of gradientmethod and N-R method are respectively 0.006 and 0.000 5，and the calculation precision of surface fittingmethod has huge uncertainty，gradientmethod ismore suitable for engineering application.

Metrology；Digital speckle correlation；Quadratic surface fitting；Gradientmethod；N-R method；Subpixel

TB93

：A

：1000-1158（2015）03-0260-08

10.3969／j.issn.1000-1158.2015.03.09

2014-08-20；

：2014-11-27

國家863計(jì)劃項(xiàng)目（2013AA065303）；國家自然科學(xué)基金（91120306）；天津市自然科學(xué)基金（14JCQNJC01600）

資新運(yùn)（1971-），男，湖南衡陽人，軍事交通學(xué)院教授，博士生導(dǎo)師，主要研究方向?yàn)閭鞲信c檢測技術(shù)、數(shù)字圖像處理。lookttlook＠163.com