增量數字圖像相關方法的誤差分析及應用考量1)

池宇希 楊昊天 潘 兵

(北京航空航天大學航空科學與工程學院，強度與結構完整性全國重點實驗室，北京 100191)

數字圖像相關方法（digital image correlation,DIC）通過跟蹤變形前后被測物體表面數字圖像中的信息載體（如人工散斑、自然紋理等）的位置變化來獲得其表面的全場變形[1-2]。由于具有硬件設備和測量過程簡單、測量精度高、適用測量范圍廣等優點，這種基于數字圖像的非接觸式、全場光測力學方法已成為實驗力學領域中最受歡迎的主流變形測試技術。然而在實際工程應用中，當被測物體發生大變形[3]，或應用于高溫環境時，劇烈的光照及散斑變化[4-7]導致被測物表面信息變化較大，變形后圖像會出現所謂的“退相關效應（decorrelation effect）”。此時，采用固定參考圖像的傳統DIC（即序列DIC，下同）計算策略往往會失效。為了解決DIC 在上述場景中的應用問題，研究人員發展了增量計算策略[3,5,8-9]，即增量DIC 方法。通過更新參考圖像，增量DIC方法提高了當前參考圖像與變形圖像之間的相關性，使DIC 計算能夠順利進行。在更新參考圖像后，通過將各個階段的位移場累加以得到相對于初始參考圖像的位移場，避免了固定初始參考圖像時因“退相關效應”導致的計算失效的問題。

除了在上述易發生圖像退相關的應用場景外，增量計算策略在其他DIC 領域還有廣泛應用。例如，在3D–DIC 中，左右相機分別采集了一組序列圖像，為了將左右相機的計算點關聯起來，需要在當前幀中進行左右匹配[9]，并通過左相機進行時序匹配；或通過初始參考幀進行左右匹配，并以左右視圖各自的初始幀為參考圖像進行時序匹配[10]，這些過程均會發生參考圖的更新。再如，在多相機DIC 中，相鄰相機之間兩兩匹配時，其參考圖像也各不相同，為了在多個視角中確定同一個計算點，也需要對參考圖像進行更新。除此之外，在使用數字體圖像相關（digital volume correlation, DVC）測量內部大變形[11-13]時，增量算法也被用來匹配多個加載階段的體圖像。在這些需要更新參考圖像的場景中，每次更新參考圖像時，新的圖像子區中心的坐標往往位于亞像素位置，因此需要對參考圖像子區進行插值，會帶來一定的插值誤差。由于標準的反向組合高斯–牛頓算法中[14]，圖像梯度和Hessian 矩陣均在參考圖像中計算，這些微分操作會進一步放大插值誤差，在更新參考圖像后，計算結果會呈現出混疊現象[15-16]，故從測量精準度考慮，插值操作應當盡量避免。

為了減小增量算法中由于對參考圖像插值帶來的插值誤差，Zhou 等[17]提出了自適應參考圖像子區平移計算策略。在這一策略中，在更新參考圖像時，其使用距離亞像素位置計算點最近的整像素點作為新的圖像子區中心，避免了對更新參考圖像子區的灰度值的插值計算，這一增量策略也被用于DVC 中[13]。由于平移量較小，且由平移帶來的位移差可以通過形函數進行補償。這一策略不僅消除了更新參考圖像時由于參考圖像插值導致的系統誤差（通過自適應平移將其轉換成了隨機誤差），同時也減少了更新過程中由于插值帶來的額外計算量。這一策略雖然減輕了由于參考圖像插值導致的系統誤差，但由變形圖插值帶來的插值誤差仍不可避免。另外，現有參考圖更新的策略有：連續更新、固定步長更新和條件更新（即發生圖像退相關再進行更新）等策略。其中，連續更新由于頻繁更新參考圖會帶來額外的計算量，也會帶來額外的誤差累積，而固定步長更新的步長又往往難以選取。所以，在應用增量DIC 計算時使用者會面臨兩個令人困惑的問題：（1）更新參考圖像會不會引入額外的測量誤差，以及這一誤差的分布范圍及分布規律是怎樣的？（2）為盡量減少這一誤差和計算量，應該采用何種策略更新參考圖像？顯然，這兩個問題是增量DIC/DVC 算法實際使用時必需研究和討論的問題。

針對以上問題，本文通過一組亞像素平移模擬實驗，分析并證明了圖像子區平移策略用于亞像素平移測量時由于亞像素插值導致的系統誤差的累積效應。并通過控制圖像子區平移的步長，總結了增量算法在不同步長條件下對亞像素插值誤差的累積規律。另外，通過實際實驗，研究了DIC 在實際實驗中亞像素插值誤差的分布規律，進一步證明了自適應圖像子區平移增量算法在實際實驗中將系統誤差轉化為隨機誤差的性質。最后，結合亞像素子區平移和經典的高斯低通預濾波方法，極大地降低了增量算法中的累積誤差，并給出了應用增量DIC 方法時的一些建議。

1 增量數字圖像相關方法

1.1 增量數字圖像相關計算

序列DIC 算法[14]可實現高效、精準而魯棒的位移場測量，但在一些容易出現退相關的場景，序列DIC 算法容易失效，需引入增量計算策略。在條件更新策略中，增量算法在當前圖像出現圖像退相關時會將當前圖像的前一幀圖像更新為當前參考圖像，由于DIC 匹配后，計算點往往處在亞像素位置，直接的做法是對當前亞像素位置的圖像子區像素重新采樣，以獲得新的參考圖像子區的灰度信息。

為避免參考像素在亞像素位置插值，文獻[17]提出將圖像子區中心平移到最鄰近的整像素位置。由于計算量小、實現簡單，且可避免插值誤差，這一自適應圖像子區平移策略在增量DIC 方法普遍采用。圖1 為自適應圖像子區平移策略示意圖。以計算點 (x,y) 為例，在更新參考圖像時，其在當前參考圖中的匹配點位于 (xref,yref) ，如圖1 中綠色箭頭所示。在進行圖像子區平移時，將其平移到最接近的整像素點，即 (round(xref), round(yref)) ，其中 round(x) 函數將x取整到最接近的整數。則當前計算點在平移后圖像子區的局部坐標系中的坐標 (δx,δy) 可以表示為

圖1 自適應圖像子區平移策略示意圖Fig.1 Schematic diagram of the adaptive subset offset strategy

若當前參考圖像的匹配點 (xref,yref) 相對初始參考子區的累積位移為 (uref,vref) 。更新參考圖像后，當前幀的計算結果為相對當前參考圖像的結果，其匹配點位于 (xcurr,ycurr) ，如圖1 中紅色箭頭所示。其相對當前參考子區的位移參數為{ucurr,vcurr,ux,vx,uy,vy}。另外，還需要結合形函數，得到平移子區的位移補償量，如圖1 中藍色箭頭部分所示。故在自適應圖像子區平移策略中，當前幀的計算點相對初始參考圖像的位移將由3 個部分構成。即當前參考圖像的累積位移，當前幀相對當前參考圖像的位移，以及圖像子區平移的補償位移。當前計算點相對初始參考圖像的位移可由下式給出

1.2 增量DIC 計算中的參考圖像更新策略

增量DIC 計算中更新參考圖像的策略可以歸納為兩類，即定步長更新策略和條件增量策略。其中，經典的連續更新策略可以視作步長為1 的定步長更新策略。如圖2 所示，定步長增量策略在經過給定次數的計算后，即將當前變形圖作為下一次相關計算的參考圖像。而條件增量策略一直使用當前參考圖進行計算，僅在后續計算中出現“退相關效應”時，將參考圖像更新為上一幀圖像。如果相鄰的參考圖之間相關系數較高，這些增量更新策略都能保證DIC 計算的順利進行。

圖2 增量DIC 方法中的定步長增量與條件增量策略示意圖Fig.2 Schematic diagram of the fixed-step and the conditional strategies in incremental DIC

在增量DIC 的實際使用中，當圖像變化較為平緩時，各圖之間的相關系數較高，應該選取較大的增量步長。而當圖像變化較快時，相關系數下降較快，則應選取較小的增量步長。故定步長的增量方法難以適應復雜多變的實際應用場景。而對于條件增量策略，為了確定更新參考圖的時機，需要使用一定的判據。為此，Pan 等[3]選擇使用計算區域中某個種子點的相關系數作為判據，當種子點的相關系數小于給定的閾值時，即可視為圖像發生“退相關效應”，需使用上一幀圖像作為新的參考圖像進行后續計算。這一條件增量策略可以自適應圖像的變化，并盡可能地減少更新次數，僅在需要的時候進行更新。由于更新參考圖后需要對新的參考圖的各個計算點重新預計算灰度梯度、Hessian 矩陣等，會帶來額外的計算量，而各階段的匹配誤差在更新中也會不斷累積，故在實際應用中，為了更好地使用增量DIC 方法，本文將通過實驗驗證來回答前文提出的兩個問題，即：（1）更新參考圖會不會帶來額外的誤差，以及這一誤差的分布范圍和分布規律是怎樣的；（2）為了減少這一誤差和計算量，應該采用何種策略更新參考圖像。

2 實驗與分析

2.1 模擬平移實驗

在模擬平移實驗中，本文選取了DIC challenge[18]中Sample 5 的參考圖像作為基準，并使用快速傅里葉變換方法以0.1 像素為步長平移了20 步，共平移了2 個像素。通過一組模擬平移實驗來評估增量算法在不同亞像素處進行更新參考圖像帶來的誤差的累積規律與分布特性。

對于定步長增量策略，本文設計了如下的對比實驗。其中使用傳統序列DIC 的對照組，始終以初始參考圖像為參考進行DIC 匹配。而使用增量DIC 的測試組，則以特定的步長進行更新參考圖像。同時，選取了不同的計算步長（1～9 像素）以研究更新步長的影響。其中，在每組實驗中均采用相同的計算點和計算參數。參考圖像大小為512×512 像素，計算區域為邊長460 像素的正方形區域，區域左上角像素坐標為 (30, 30)。以10像素為步長進行計算點采樣，共有2209 個計算點。由于不同的插值方法會導致不同水平的插值誤差，實驗統一采用雙三次B 樣條插值作為圖像灰度插值方式。

使用不同的計算方案完成計算后，統計了圖像中各計算點的平均位移，并與真實位移做差，得到了各組實驗的平均誤差。將各組平移實驗的平均誤差繪制如圖3 所示，其中橫軸為亞像素平移量，縱軸為各組實驗的平均誤差。Seq 表示使用序列DIC 方法，當前圖像相對初始參考圖像的亞像素位移分布，而Inc 表示使用增量DIC 方法時，當前圖像相對上一張參考圖像的亞像素位移分布。圖中黑色粗線為標準序列DIC 方法的誤差計算結果，可以看到，序列DIC 方法的平均誤差呈現出類似正弦函數的波動，這與已有DIC 插值誤差的認識一致。而其余各組使用不同增量步長計算方案的結果均表現出對這一插值誤差的累積現象，即每次更新參考圖后，該參考圖的插值誤差會被累積到后續計算中，而最后的總誤差即為各階段的誤差之和。這一結論是符合直覺的，增量算法不僅累積了位移，也同時累積了各個階段的插值誤差。這進一步說明了在使用增量DIC 算法時，應盡可能減少更新參考圖的數量，僅在必要時更新參考圖，以避免插值誤差的累積。值得注意的是，從式（2）可知，增量DIC 的誤差還包括使用形函數傳遞的補償位移導致的誤差，但這一誤差相對亞像素插值誤差的量級而言可以忽略不計，故不再贅述。

圖3 采用不同增量步長的增量DIC 方法與傳統序列DIC 方法的平均誤差對比Fig.3 Comparison of average error between incremental DIC method with different step sizes and traditional sequential DIC method

2.2 真實平移實驗

在實際實驗中，各計算點的亞像素位移并非如模擬平移實驗一樣集中分布在給定位移附近，而是會存在一定的分布，故插值誤差也會呈現出與模擬平移實驗不同的分布。由前述分析可知，插值誤差存在累積效應，且可以根據計算點的亞像素位移來估計其對應的插值誤差。為研究真實平移實驗的誤差分布特性，對一固定在平移臺上的平板進行平移實驗。相機型號為FLIR-GS3-PGE-91QS6M，分辨率為3376×2704 像素，正對試樣進行拍攝。實驗中，試件沿著水平方向以0.2 mm 的步長進行了5 次平移。對于各組圖像，分別采用序列DIC 方法和連續增量DIC 方法進行計算，并統計各計算點的位移分布（計算區域共有74×71 個計算點，計算步長為15 像素，圖像子區為51×51 像素）。

在計算完成后，統計了兩種方法得到的亞像素位移的分布，如圖4 所示?？梢钥吹叫蛄蠨IC方法的亞像素位移隨著平移量的增大，逐漸趨于均勻分布U(0,1)。這是因為在真實實驗中，相機光軸與試樣平面往往不是完全垂直，且存在鏡頭的放大率、畸變等多種因素的影響，經過一定的位移后，圖像的亞像素位移將趨于隨機的均勻分布。然而，在增量DIC 方法中，由于每次均更新參考圖像，每一幀相對上一幀的位移較小，其增量位移的亞像素位移分布相對集中。圖中可以看出，各次增量的位移主要分布在大于0.5 像素的區間，經過多次累積后會導致所得位移偏小。接下來，將通過前述的分析來估計其累積誤差的分布。

圖4 使用傳統序列及增量DIC 方法在各個平移階段的亞像素位移分布Fig.4 Subpixel displacement distributions obtained by traditional sequential and incremental DIC methods in real translation experiment

根據前述的模擬平移實驗結果可知，亞像素位移的插值誤差與亞像素位移之間的關系可用正弦函數來近似，而插值誤差的量級與圖像的噪聲水平和所使用的插值函數有關。根據前述的實驗結果，插值誤差可以近似表示為

其中 du為亞像素位移部分，e為估計的插值誤差。

對于序列DIC 方法，其插值誤差可以直接通過式（3）進行估計。而對于增量DIC 方法，則需要將各增量步的插值誤差進行累加，使用式（3）估計的插值誤差的分布如圖5 所示。由于正弦函數的非線性效應，序列DIC 方法中均勻分布的亞像素位移的插值誤差呈現出兩邊高中間低的特性。而對于增量DIC 方法，由于每次更新參考圖像均會帶來一定的增量誤差，故在累加后會得到一個范圍更廣的累積誤差分布。另外，由圖4 的結果可知，增量DIC 方法的亞像素位移分布較為集中且偏向于大于0.5 像素的區間。由式（3）計算得出其對應的插值誤差會有明顯的偏移，這導致累積后的插值誤差相對序列DIC 方法呈現出較為明顯的偏離。在經過4 次參考圖像更新后，累積插值誤差的最大值已經達到 0.02 像素，且有50%以上的計算點的插值誤差大于0.01 像素。這一累積誤差的量級已經足夠體現出明顯的差異，會進一步影響后續DIC 應變計算的精度。所以，在實際使用增量DIC 時，應盡可能地減少增量次數，以減少連續幀之間亞像素位移分布集中導致的累積誤差。

圖5 使用傳統序列和增量DIC 方法的累積插值誤差分布Fig.5 Cumulative error distributions obtained by traditional sequential and incremental DIC methods

為了降低這一累計誤差的影響，Pan[19]提出使用高斯低通預濾波對圖像進行預處理。這一預處理手段可以將插值的系統誤差轉換為隨機誤差，體現為位移結果的平均偏差降低，但標準差變大。文中指出，通過選取合適的高斯核大小，可以在不顯著增加標準差的同時大幅降低偏差。為此，本研究對前述實驗圖片使用5×5 像素高斯核函數進行了高斯低通預濾波處理，并統計了增量DIC與序列DIC 方法之間的相對誤差。為了直觀展示二者之間的誤差分布，圖6 展示了誤差的累積概率分布曲線，其中橫軸為相對誤差的分布范圍，縱軸為對應的累積概率。其中實線為圖像未經高斯濾波的結果，其相對誤差較大且分散。另外，隨著增量次數的增加，增量算法的最大相對誤差也逐漸增大，經過4 次增量計算后最大相對誤差約為0.04 像素。虛線為圖像經過高斯預濾波之后的結果，經過高斯預濾波后，絕大多數點的相對誤差均落在±0.01 像素之間，且增量算法并未導致明顯的相對誤差增加。結果證明，在使用增量DIC 算法時配合高斯低通預濾波算法，除了顯著降低由于亞像素插值帶來的誤差之外，還可以降低由于增量算法帶來的累積誤差。

圖6 使用傳統序列和增量DIC 方法的相對誤差的累積概率分布Fig.6 Cumulative probability distributions of the relative displacement error between traditional sequential and incremental DIC methods

2.3 橡膠拉伸大變形的增量DIC 測量

在材料測試實驗中，增量DIC 常用于存在大變形的實驗場景。為了展示增量DIC 算法在實際實驗中的適用性，使用增量DIC 方法計算了一組橡膠材料的大變形拉伸實驗圖像。橡膠試樣按ISO 37:2005 標準切割為標準啞鈴型拉伸試件，其尺寸和外形如圖7(a)所示，厚度為2 mm，并在表面使用噴漆制作了隨機的散斑。如圖7(b)所示，實驗中使用一臺SUNS-UTM5 萬能實驗機以2.5 mm/min 的速度勻速拉伸試件，最大拉伸位移為68 mm。實驗中采用的相機型號為FLIRGS3-PGE-91QS6M（分辨率為3376×2704 像素，裝有Schneider 23mm F/1.4 鏡頭），其固定在試件前方約1.2 m 處，正對試件表面進行采集，按每分鐘15 幀的采集速率進行采集，實驗全過程一共采集了407 幀圖像。

圖7 (a) 橡膠試件的尺寸與實際圖像和(b) 實驗現場裝置Fig.7 (a) The dimensions and physical image of the rubber specimen and (b) the experimental setup

在匹配計算中，圖像子區大小為41×41 像素，步長為5 像素，所有圖像均經過了3×3 像素的高斯低通預濾波處理。由于橡膠試件變形過大，常規DIC 方法在初值估計階段即會失效。即使采用時序傳遞策略（使用上一幀的結果作為初值），在變形過大時仍會由于形函數欠匹配、采樣失真等原因導致相關系數下降，從而發生“退相關”。而采用增量算法，則可以順利完成全部計算。在實驗中采用了條件增量策略，即當計算區域中相關系數最低的計算點的相關系數低于0.95 時，則采用上一幀的圖像作為參考圖。這一條件判據相對使用特定種子點的判據[3]，更能體現全場的匹配質量，避免由于部分計算區域已“退相關”而種子點仍保持較高相關系數而帶來的誤判。所得部分計算結果如圖8 所示，可見增量DIC 算法在經過多次增量后能得到實驗全過程的變形，克服了傳統DIC 方法在大變形時難以繼續計算的問題。

圖8 實驗過程中橡膠試件沿拉伸方向的位移(單位：mm)Fig.8 The displacement of the rubber specimen in the stretching direction during the tensile test (unit: mm)

除此之外，我們還選取了不同的增量步長來計算同一組數據，以展示不同的增量步長在實際實驗中的累積效應。在這里，我們設計了7 組定步長的測試組，其步長分別為1, 2, 4, 8, 16, 32,64 幀。圖9(a) 展示了對實驗的最后一幀，計算區域橫向中心線上的數據點（如圖8 中第0 幀所示），使用不同增量策略得到的位移結果。其中，中間的黑色線條為采用條件增量策略的增量方法，各計算點位移沿著拉伸方向表現出非常好的線性關系。在這一均勻拉伸實驗中，各點位移應與其坐標成正比，故這一結果符合預期。而對于其他步長的結果，則表現出不同水平的累積誤差，并沿著拉伸方向有周期性的變化。其中，在序號為50 的計算點處，各種步長的計算結果均相同。將這一點使用不同固定增量步長及使用條件增量步長得到的相關系數繪制如圖9(b)所示，這里采用了零均值歸一化互相關（zero-mean normalized cross correlation, ZNCC）相關系數?？梢?，在不更新參考圖時，隨著變形逐漸增大，相關系數會不斷下降，而更新參考圖后，相關系數又會回到1 附近。使用條件增量策略時，僅更新了4 次參考圖，相比使用不同增量步長的其他實驗組，條件增量策略可以盡可能少地進行增量計算，同時也能獲得盡可能少的累積誤差。除此之外，還統計了使用不同增量步長計算全部406 幀圖像的計算時間，列于表1 中?？梢?，使用條件增量策略由于減少了增量的次數，可以盡可能地減少計算量，提高計算效率。

表1 使用不同固定增量步長及條件增量策略的計算耗時Table 1 Computation time of fixed-step strategy of different steps and conditional strategy

圖9 (a) 不同位置的計算點使用不同更新策略下的位移和(b) 70 號計算點使用不同更新策略下全過程的ZNCC 相關系數Fig.9 (a) Displacements of different calculation points with different updating strategy and (b) ZNCC of #70 calculation point with different updating strategy during the experiment

值得注意的是，圖9(a)中不同步長的累積誤差在空間上呈現出周期性的波動，且其周期與增量步長成反比例關系。這是由于在這一實驗中采用了穩定的加載速度和固定的采集幀率，各計算點的增量位移保持相對穩定。圖10(a)和圖10(b)分別展示了使用條件增量和連續增量策略時不同位置全過程的增量位移，即每一幀圖像位移與上一幀圖像的位移差。從圖10(a)可以觀察到，增量位移的大小也與其橫向坐標成正比，這一結果符合均勻拉伸的假設。且序號為20 和序號為80的計算點的增量位移分別保持在0.5 像素和1.5像素附近，這表示間距為60 個計算點的兩點的平均增量位移差約為1 像素，故這兩點應具有相同的誤差分布。這一現象與圖9(a)中不同步長的增量結果出現周期性變化相符。由于勻速加載和拉伸試驗的加持方式，增量位移與加載時間和空間坐標均呈線性關系，導致不同步長的累積誤差呈現出周期性變化。其中采用連續增量策略的結果變化周期約為60 個計算點，而對于其他步長的增量算法，由于其間隔時間變長，兩點的增量位移差會同比增長，而變化周期會同比縮短。

圖10 沿拉伸方向不同位置計算點各幀的增量位移：(a) 使用條件增量策略, (b) 使用連續增量策略Fig.10 Incremental displacement of different calculation points along stretch direction with (a) conditional strategy and (b)sequential strategy

另外，對于序號為50 的計算點而言，其各幀的增量位移恰好為1 個像素。按前述的理論估計，該點的累積誤差應該接近于0，這也與圖9(a)的結果相符。另外，在圖10(b)中也可看出，該點在連續增量中也保持穩定，但對于其他位置的計算點而言，均存在一定的計算誤差。圖中可以看出對于各位移均有向最近的整像素偏移的系統誤差，而經過約400 次連續增量后，累積誤差最大將近80 像素，如圖9(a)中所示。這一結論再次證明了增量算法中插值誤差的累積效應，也證明了條件更新策略的優越性。

3 結論

增量DIC 計算是解決大變形、環境光變化劇烈、散斑退相關等非常規場景中變形測量的有效手段。本文通過一系列的模擬及真實實驗，分析了增量DIC 計算誤差的理論累積規律，并通過真實實驗展示了增量DIC 方法在實際實驗中由亞像素位移分布所引起的累積誤差的分布和形式，證明了高斯低通預濾波對這一累積誤差的抑制效果。另外，還展示了增量DIC 算法在真實大變形實驗中的應用，并通過對比實驗證明了條件更新策略相比固定步長更新策略的優勢。

根據上述研究，本文給出使用增量DIC 的如下建議。

（1）增量DIC 方法由于包含亞像素位移插值過程，存在穩定的系統誤差，并會隨著增量的過程不斷累積。

（2）更新參考圖像會不可避免地增加測量誤差和計算量，應盡可能少地更新參考圖像。在實際使用增量DIC 時，建議采用條件更新策略，如使用種子點或計算區域中相關系數最小的計算點的相關系數作為更新判據。

（3）更新參考圖時，應采用隨機圖像子區平移策略將系統誤差轉化為隨機誤差，結合圖像高斯低通預濾波處理，能顯著降低由于插值帶來的系統誤差，以獲得更好的誤差抑制效果。

另外，增量DIC 算法不僅可用在容易發生退相關的測試環境，在3D–DIC、多相機DIC 等需要統一不同參考圖中的同一計算點坐標系的場景也適用。在這些場景中，由于存在多個視角，故需要選定一主相機，其他相機的參考圖均與主相機之間進行一次增量更新，再各自進行時序匹配。由于各相機時序匹配與2D 匹配過程相同，在多個相機之間的增量更新及后續的時序更新中，本文給出的建議仍然適用。