快速DIC形變測試技術與工程應用
1李得睿 2陳衛紅 2彭波 3晏班夫
（1上海交通大學船舶海洋與建筑工程學院，上海，200240；2北京智博聯科技股份有限公司，北京 100088；3湖南大學土木工程學院，長沙，410082）
[摘 要] 數字圖像相關（Digital image correlation, DIC）技術以其亞像素級別測試精度、可遠距多點非接觸測試、設備簡單、操作簡便等特點，在結構形變測試中受到較多關注。本研究采用基于傅里葉變換的互相關(Fourier transform-based cross correlation, FTCC)算法與反向組合高斯牛頓 (Inverse Compositional Gauss–Newton, IC-GN)算法的DIC運算框架，對DIC在工程應用中的可行性進行分析。通過精度分析以及實際工程應用案例，對基于FTCC和IC-GN算法的DIC計算架構進行測試，結果表明DIC技術在工程應用中的測試精度有理論保證，且得到了較好的實測結果。
[關鍵詞] 數字圖像相關；形變測試；亞像素匹配
Rapid DIC Deformation Measurement Technology and Engineering Application
1Li Derui 2Chen Weihong 2Peng Bo 3Yan Banfu
（1 School of Naval and Architectural , Ocean & Civil Engineering, Shanghai Jiaotong University ,Shanghai，200240； 2 Beijing ZBL Science & Technology Co. Ltd.，Beijing 100088； 3 College of Civil Engineering, Hunan University, Changsha,410082）
Abstract: Digital image correlation technology has attracted more attention in structural deformation test because of its sub-pixel level test accuracy, remote multi-point non-contact test, simple equipment and easy operation. In this study, DIC operation framework based on cross-correlation algorithm of Fourier transform and inverse Compositional Gauss-Newton algorithm is used to analyze the feasibility of DIC in engineering application. Through the accuracy analysis and practical engineering application cases, the DIC computing architecture based on FTCC and IC-GN algorithm is tested. The results show that DIC technology has theoretical guarantee for the test accuracy in engineering application, and good measurement results are obtained.
Keywords: digital image correlation; deformation measurement; subpixel matching
前 言
形變測試是研究各類土木工程材料與結構特性的重要手段。傳統形變測試方法普遍存在現場操作困難、多點測試不便、測試成本高、需人造靶標、測試效率低下等問題。在結構測試領域，近年來數字圖像相關（Digital image correlation, DIC）方法以其亞像素級別測試精度、可遠距多點非接觸測試、設備簡單、操作簡便等特點，在結構形變測試中受到較多關注
[1-3]。
依據算法原理，可將DIC技術分為局部DIC技術
[4]與全局DIC技術
[5], [6]；依據拍攝原理，又可將局部DIC技術分為二維局部DIC技術
[7]、立體局部DIC技術
[8], [9]與數字體積相關(Digital Volume Correlation, DVC)技術
[10]。其中，全局DIC技術并非主流DIC技術，目前大部分學術論文內容以及幾乎所有商業DIC軟件均基于局部DIC技術實現；而立體局部DIC技術與二維局部DIC技術并無本質區別，立體局部DIC技術需要進行多相機組的外參數標定，來實現三維位移的測量，其本身仍運行二維局部DIC的算法；DVC技術實質上是DIC技術的三維拓展，其基于類似X射線對物體進行三維掃描，得到物體的三維實驗數據，運行二維DIC算法的三維拓展版本，實現對物體內部的三維運動測試，DVC難以在土木工程測試領域推廣，其所需設備非常規設備，并且大部分土木工程測試情景不需要運用DVC進行三維體積位移追蹤。基于此，本文主要基于二維局部DIC技術，進行土木工程結構位移測試的研究與可行性分析。為了方便表述，本文一律將二維局部DIC技術稱之為DIC技術。
晏班夫等
[3]為減少基于快速歸一化互相關（Fast normalized cross correlation, FNCC）
[11]的整像素并行搜索算法的運算量，基于CUDA 并行計算平臺，使用消費級NVIDIA GPU 顯卡，將圖像配準領域廣泛應用的基于快速傅里葉變換的互相關（Fourier transform-based cross correlation，FTCC）
[12]整像素初值搜索算法引入GPU-DIC并行運算框架中，大大提升了初值搜索速度，結合反向組合高斯-牛頓（Inverse compositional Gauss–Newton，IC-GN）
[13]亞像素算法，位移測試達到了與傳統FNCC并行算法相當的精度與穩定性主要，通過實驗室實驗證明了DIC在土木工程形變測試領域的實用性。
本文主要采用理論精度分析與現場實驗對集成FTCC和IC-GN算法的DIC計算架構進行測試，以驗證該測試體系以及測試系統在工程應用場景下，尤其是土木工程室外場景下結構位移測試的可行性。實驗結果表明，該DIC計算架構具備在土木工程室外場景下的測試可行性。
1 原理簡介
1.1 硬件系統
圖 1為典型的數字圖像檢測系統，主要由筆記本電腦（計算設備）、數字圖像采集儀、被測物以及照明設備（光源）等硬件設備及數字圖像測量軟件組成，其中數字圖像測量軟件是核心。圖2為基于GPU并行運算的快速DIC形變測試研究框架及流程，其中DIC技術的核心算法包括整像素初值搜索算法（FNCC及FTCC）與亞像素匹配算法（IC-GN與FA-NR）。亞像素算法中應用最廣泛的是基于形函數的迭代型算法，例如IC-GN、FA-NR等。得到的位移場經過差分及平滑可以得到應變場，這部分內容在本文并未涉及。
圖 1 數字圖像檢測系統基本構成
圖 2 ZNCC遍歷運算示意圖
1.2 DIC整像素算法
1.2.1 傳統FNCC算法
如圖 2所示，
是尺寸為
像素的搜索區域，也即匹配圖像模板。從搜索區域中截取出其中藍色虛線區域，作為引用子集
，引用子集尺寸為
像素，
中黃色實線區域尺寸與引用子集
尺寸相同。整像素搜索的核心目標，即定量地找出
的中心點，在
中的整像素坐標
該坐標也就是整像素初值搜索的結果。容易想到，將
中與引用子集
尺寸相同的黃色實線區域，作為目標區域，從左至右、從上到下地逐個像素點遍歷全部的
，在每個遍歷路徑上的像素點處，將目標區域與
做式的ZNCC運算，即
其中，
，
，
為
的灰度均值，
為
中黃色實線目標區域的灰度均值，
隨
與
的變化而發生改變。
遍歷過后，即得到一個
維的
系數方陣，該方陣中，假設
最大值處的坐標為
，則整像素初值搜索結果
為
可以看到，由式直接表述的ZNCC遍歷算法，至少需要
次加法運算與
次乘法運算，且計算機乘法運算耗時遠大于加法運算。在實際DIC形變測試中，一般有，
，
，式的運算量隨著
的增大，呈平方倍增長，隨著
的增大，呈4次方趨勢增長。由此可知，式的直接實現方式并非高效。
Lewis
[11]于1995年提出了快速歸一化互相關(Fast normalized cross correlation, FNCC)整像素初值搜索算法，FNCC算法即式的快速算法。FNCC將式的分子基于快速傅里葉變換進行計算，采用圖像積分表快速計算式的分母部分。
實際運用FNCC算法進行整像素初值搜索時，
的值代表了搜索范圍，
越大，則FNCC的搜索能力越強，但相對應的運算量也就越大，且運算量隨著
與
的增大呈非線性快速增長。
越大則意味著做相關運算的數據越多，也即魯棒性越強，
越大則意味著搜索范圍越廣，增大了受干擾的可能性，也即魯棒性變差。實際運算時，
的取值范圍一般為15-71，普遍取
。
的取值可表示為
，
表示搜索范圍，單位為像素，
的大小則隨測試情況而定。
若兩幀圖像間，被測物位移變化較大，如混凝土加載試驗中，混凝土梁臨近破壞前的最后幾幀圖像，由于混凝土梁此時剛度很小，所以兩幀圖像間混凝土梁的變位較大，
取50也有可能達不到相應的搜索范圍要求，且此時FNCC計算量非常大，魯棒性也受到較大影響。由于土木工程測試情形，普遍具備小位移、慢位移、小旋轉的特點，FNCC算法還是能夠滿足絕大多數土木工程形變測試場景的。在進行高速沖擊試驗時，則需要運用高速攝像機進行圖像的采集，來實現兩幀間小位移的條件。
1.2.2 FTCC算法
針對土木工程試驗中的大剛體位移情形，晏班夫等
[3]引入基于傅里葉變換的互相關(Fourier transform-based cross correlation, FTCC)算法，作為DIC整像素初值搜索算法，并證實了FTCC算法可實現DIC技術在大剛體位移情形下的快速計算。
FTCC被Reddy等
[12]提出后，被廣泛用于計算視覺中的圖像配對與校準工作，潘兵等
[13]將FMT-CC引入RG-DIC算法中，實現了目標物存在較大平移、旋轉情形下的初值搜索。FTCC算法利用引用子集
與當前待匹配子集
的重疊區域，將其頻域內的互功率譜進行Fourier逆變換，得到脈沖函數，通過檢測該函數峰值的位置，實現空域內圖像二維位移的檢測。
設兩子集尺寸相同，且
為
在
和
方向分別平移
和
后的位置，即:
（2）
（3）
一般平移不會改變頻譜幅度，僅改變相位，進一步計算兩者的互功率譜，即
（4)
其中
代表
的復共軛。可以看出，互功率譜在頻域中的相位相當于引用子集與當前子集在空域中的相位差，于是對該互功率譜進行二維逆傅里葉變換，得到一脈沖函數，其峰值即對應子集在空域中的相對平移量
和
。
理論上FTCC假定兩個待匹配子集之間所有數據相同，有效數據周圍區域拓延補零，但具體應用到DIC中的FTCC程序，所處理的均是具有部分重疊區域的圖像，其余部分為不相關數據。計算仿真表明，FTCC算法要求相關子集之間至少要有大于25%的重疊區域，即
與
均小于相對應維度的子集尺寸的一半，且重疊區域越高，脈沖峰值越大。這意味著FTCC的搜索范圍為子集尺寸的一半。
1.3 DIC亞像素算法
Bruck等
[14]于上世紀八十年代提出了正向累加牛頓拉夫森(Forward additive Newton–Raphson, FA-NR)算法，用于將模板圖像快速匹配至目標圖像，該算法的核心思想在于迭代優化求局部最優解。潘兵
[15]將IC-GN算法
[16]做了相關改進并將其引入DIC理論，作為一種相對于FA-NR算法，效率更高、魯棒性更好的DIC亞像素匹配算法。本文基于ZNSSD相關度函數，采用一階形函數，對DIC中的IC-GN算法進行詳細介紹。
采用ZNSSD相關度函數，基于式，定義如下相關度函數
（5）
式中，
與
分別為引用子集與目標子集的灰度值，即圖 2中引用子集與黃色實線目標區域，
為各子集的中心點坐標的齊次形式。
為各子集內像素點到中心點的原始距離。
為
的變化量。
、
、
、
與式中分子、分母項一一對應。
對式基于
做一階泰勒展開并舍棄高階小量，可得
（6）
其中，
代表全子集區域求和，
為引用子集
在點
處的灰度值，
為引用子集
于點
處分別在x與y方向的灰度梯度值，其余參數含義同式。
與傳統FA-NR算法
[14]思路相同，若求得式局部極值，即求解
此時有
（7）
其中，Hessian矩陣具體形式為
經過式計算后，得到
，帶入warp函數形成增量warp函數
，基于
對原
進行更新，即
（8）
將
中的
帶入式開始新的迭代，迭代終止條件為迭代次數大于10次或歐式距離誤差小于0.001。
經過上述迭代過程，IC-GN算法可收斂至局部最優解，該最優解使得式達到局部最小值。與FA-NR相似，一般來講，迭代初值與正確結果相差幾個像素以內，IC-GN算法均能成功收斂至正確的亞像素結果。
通過上述IC-GN計算過程，可以看到，不同于FA-NR，IC-GN算法的每一次迭代運算中，與
有關的項僅為
，Hessian矩陣和灰度值梯度項
均在迭代開始前計算完成即可，所以IC-GN避免了大部分冗余計算，運算效率高于經典FA-NR算法。基于零階與二階形函數的IC-GN算法，原理上沒有改變，僅
形式以及相關運算矩陣維度發生了變化。基于零階形函數的IC-GN算法詳細推導過程請參考文獻[17]，基于二階形函數的IC-GN算法請參考文獻[18]，其中詳細記錄了二階形函數的
矩陣形式。
1.4 DIC算法分析
亞像素算法迭代收斂與否取決于整像素搜索算法提供的迭代初值是否合理。理想狀態下，整像素初值搜索算法所提供的整像素初值結果與正確結果之間的誤差不大于0.5個像素，收斂半徑一般在幾個像素以內。FTCC整像素算法利用引用子集與當前子集的重疊區域，將其頻域內的互功率譜進行Fourier逆變換，得到一個脈沖函數，通過檢測該函數峰值的位置，實現空域內圖像二維位移的檢測。
用整像素匹配算法所能達到的檢測精度，一般應通過迭代優化型亞像素匹配算法滿足測量需求，主要的亞像素算法包括正向算法FA-NR （Forward Additive Newton–Raphson）與反向算法IC-GN兩種，由潘兵等人于2013年提出的反向算法IC-GN是在正向算法FA-NR的基礎上發展起來，其主要思想就是采用迭代優化的方式求得某個判據的極值，目前普遍采用的判據為零均值歸一化最小平方差（zero-mean normalized sum of squared difference, ZNSSD），ZNSSD與ZNCC可互推導，除了具有ZNCC的優勢外，其應用于IC-GN中，更易于求導計算。
不同于FA-NR，IC-GN基于參考子集進行迭代計算，省去了FA-NR算法中每次迭代都要進行的Hessian矩陣的計算，從而大大降低了計算量，同時，IC-GN比FA-NR具有更好的魯棒性。
2 工程應用可行性分析
實際工程中，采用DIC技術進行測試時需要面對各種不同的測試場景，這些測試場景的主要區別在測試距離以及視場寬度，本節內容主要基于民用級數字相機級鏡頭的具體參數，結合不同的測試距離，對DIC技術進行室外場景測試時的精度進行具體分析，以此來驗證DIC技術在土木工程室外測試情景下的可行性。
DIC技術的精度驗證主要由硬件參數、理論算法以及測試距離決定，硬件參數方面，主要有分辨率、感光芯片像元尺寸以及鏡頭焦距等；理論算法方面，通過IC-GN算法可將整像素級別測試精度提高100倍，這也使得DIC技術能夠進行高精度測試；測試距離方面，一般在0.5-1000m不等。
基于小孔成像原理，結合相關硬件及拍攝參數，可計算得到測試精度如下：
（9）
式中，
為測試精度，單位為(距離/像素)，
為測試距離，
為感光芯片像元尺寸，
為相機焦距，
為亞像素精度放大系數。實際測試中，
一般可取為2.6
，
取為0.01，則可得不同焦距與拍攝距離下，DIC技術的理論測試精度，計算結果示于表 1。
同樣基于小孔成像原理，結合相關硬件及拍攝參數，可計算得到拍攝視場寬度如下：
（10）
式中，
為拍攝視場寬度，單位為公米制實際距離，
為測試距離，
為感光芯片像元尺寸，
為相機焦距，
為相機感光芯片橫向像素個數。實際測試中，
一般可取為2.6
，
取為1920像素，則可得不同焦距與拍攝距離下，進行DIC測試時的理論視場寬度，計算結果示于表 2。
表 1 不同焦距(mm)與測試距離(m)下的DIC理論測試精度(mm/pixel)
表 2 不同焦距(mm)與測試距離(m)下的DIC測試視場寬度(m)
需要強調的是，表 1與表 2的理論計算結果均基于一般情況下的硬件參數計算得到，實際應用中可通過改變
與
得到更高測試精度與不同的視場范圍，盡管如此，表 1與表 2仍具有較強的參考價值。通過表 1數據可以看出，DIC測試精度水平極高，在小焦距、長距離的極端情況下，也可達到毫米級的測試精度。基于表 2數據，結合表 1數據，可以看出，在滿足測試精度的前提下，DIC測試的視場范圍變動較大，可滿足絕大多數的實際室外測試情景。同時，表 1及表 2數據為正對拍攝情形下的理論計算結果，如果是斜對拍攝，則精度會受小孔成像投影影響而有所降低，拍攝視場寬度則會有所增大。
3 工程應用
3.1 蘇通大橋電渦流阻尼器運動姿態測試
圖 3 蘇通大橋電渦流阻尼器運動姿態測試
運用多點數字圖像檢測系統對蘇通大橋電渦流阻尼器的變位視頻進行了位移測試。由于一般情況下電渦流阻尼器的變位十分緩慢，肉眼難以察覺，所以該阻尼器變位視頻以30fps的幀率，1080p的分辨率對電渦流阻尼器進行長達十分鐘的拍攝，拍攝畫面以及處理結果示于圖 4 (a)，現場拍攝的電渦流阻尼器及其支架示于圖 4 (b)。
因為電渦流阻尼器的三角結構，所以僅需拍攝其兩斜邊的剛體位移，就可以把電渦流阻尼器在焦平面內最重要的二維運動姿態通過攝影的方式記錄下來，再利用DIC進行跟蹤即可。圖 4 (a)中阻尼器兩斜邊上方的藍點與綠點為追蹤目標點，又由于橋面上風速較高，會使得相機產生小幅高頻振動，為了消除相機自身振動的影響，在底邊橫軸上選取了紅點作為基準點。最終將兩目標點的絕對位移結果減去基準點絕對位移結果，即得到了圖 4 (a)中的兩條相對位移時程曲線，即阻尼器兩斜邊相對于底邊的相對位移結果。圖 4 (b)中的兩條位移時程曲線是阻尼器兩斜邊的豎向位移時程結果，均值化后的曲線極值分別為[-4.49,4.52]mm與[-3.16,4.19]mm，可以看出，該阻尼器在環境激勵的位移時程以及運動姿態均可較精確測得。
由于阻尼器變位緩慢，相機以30fps幀率拍攝得到的畫面在局部時間內基本處于靜止狀態，畫面穩定，質量較高，此時可將ZNCC閾值適當調高(如
，且因為兩幀之間變化緩慢，可采取小參數的FNCC算法結合IC-GN的單點DIC算法，進行阻尼器變位的多點DIC跟蹤測試。
3.2 虎門大橋渦振視頻跟蹤
2020年5月5日下午14時左右，廣東省虎門大橋的加勁梁發生明顯振動，據網絡信息報道，專家初步推測的原因為水馬改變了梁體的氣動外形。本文基于網絡視頻(視頻素材來源：B767-400ER bilibili)，運用多點數字圖像檢測系統對虎門大橋渦振振幅與頻率進行了分析，如圖8所示。在圖中選擇多個密集點進行DIC位移時程跟蹤，進而將每一時刻的空間位置點連線，可以得到大橋的振動模態。可以看出，虎門大橋渦振頻率較為單一，經過分析，在該視頻拍攝時刻，目標點處最大振幅為37cm(簡單標定結果)，振動頻率為0.365Hz，對應虎門大橋第三階加勁梁豎彎對稱振型。
由于該視頻較為模糊不清，所以ZNCC閾值設為0.7，運算過程中，采用零階形函數進行亞像素匹配運算，由于圖像畫面中的散斑數據穩定但較為模糊，可采取中等尺寸窗口(30×30)參數的FTCC算法結合IC-GN的單點DIC算法。
圖 4 DIC密集目標點跟蹤結果
4 結論
1）基于理論公式對DIC技術在實際工程應用場景下的測試可行性進行了分析，通過定量分析數據得到，DIC技術在室外場景的拍攝情形下具備較大的測試可行性，DIC測試視場和測試精度均可在一般拍攝條件下達到工程測試需求。
2）基于多點數字圖像檢測系統進行了蘇通大橋電渦流阻尼器運動姿態測試、虎門大橋渦振視頻測試等實際土木工程形變測試，實際試驗結果證明，DIC形變測試系統具備較高可行性與實用性。
參考文獻
[1]. 邵新星, 陳振寧, 戴云彤, et al. 數字圖像相關方法若干關鍵問題研究進展[J]. 實驗力學, 2017(3)：305-325.
[2]. 葉肖偉, 董傳智. 基于計算機視覺的結構位移監測綜述[J]. 中國公路學報, 2019, 32(11).
[3]. 晏班夫，李得睿，徐觀亞，朱平，邵旭東. 基于快速DIC與正則化平滑技術的結構形變測試；DOI:10.19721/j.cnki.1001-7372.2020.09.019；中國公路學報；2020年第9期；
[4]. Bruck H A, McNeill S R, Sutton M A, et al. Digital image correlation using Newton-Raphson method of partial differential correction[J]. Experimental mechanics, 1989, 29(3): 261-267.
[5]. Sun Y, Pang J H L, Wong C K, et al. Finite element formulation for a digital image correlation method[J]. Applied optics, 2005, 44(34): 7357-7363.
[6]. Réthoré J, Hild F, Roux S. Extended digital image correlation with crack shape optimization[J]. International Journal for Numerical Methods in Engineering, 2008, 73(2): 248-272.
[7]. Pan B, Qian K, Xie H, et al. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review[J]. Measurement science and technology, 2009, 20(6): 062001.
[8]. Orteu J J. 3-D computer vision in experimental mechanics[J]. Optics and lasers in engineering, 2009, 47(3-4): 282-291.
[9]. Sutton M A, McNeill S R, Helm J D, et al. Advances in two-dimensional and three-dimensional computer vision[M]//Photomechanics. Springer, Berlin, Heidelberg, 2000: 323-372.
[10]. Bay B K, Smith T S, Fyhrie D P, et al. Digital volume correlation: three-dimensional strain mapping using X-ray tomography[J]. Experimental mechanics, 1999, 39(3): 217-226.
[11]. Lewis J P. Fast Normalized Cross-Correlation. Industrial Light & Magic (1995).
[12]. Reddy B．S．, Chatterji B．N． An FFT-based technique for translation，rotation，and scale-invariant image registration． IEEE Transactions on Image Processing，1996，5(8)：1266～1271
[13]. Pan B, Wang Y, Tian L. Automated initial guess in digital image correlation aided by Fourier–Mellin transform[J]. Optical Engineering, 2017, 56(1):014103.
[14]. Bruck H A, McNeill S R, Sutton M A, et al. Digital image correlation using Newton-Raphson method of partial differential correction[J]. Experimental mechanics, 1989, 29(3): 261-267.
[15]. Pan B, Li K, Tong W. Fast, Robust and Accurate Digital Image Correlation Calculation Without Redundant Computations[J]. Experimental Mechanics, 2013, 53(7):1277-1289.
[16]. Baker S, Mattews I (2004) Lucas-Kanade 20 years on: a unifying framework. Int J Comput Vision 56:221–255.
[17]. Pan B., Tian L, Song X. Real-time, non-contact and targetless measurement of vertical deflection of bridges using off-axis digital image correlation[J]. Ndt & E International, 2016, 79:73-80.
[18]. Gao Y, Cheng T, Su Y, et al. High-efficiency and high-accuracy digital image correlation for three-dimensional measurement[J]. Optics and Lasers in Engineering, 2015, 65: 73-80.
相對應的二維傅里葉變換
有如下關系：
