代延輝，許 蔚，李 坤

（昆明理工大學 建筑工程學院，云南 昆明 650500）

0 引 言

如何精準獲得材料的力學性能始終是當今研究的熱點問題之一，根本原因在于無法準確測量其內(nèi)部變形，傳統(tǒng)的測量方法（如應變片法、散斑干涉法、光彈性法等）具有一定的局限性，導致測量結果與實際變形相差甚大。近年來，隨著三維層析設備與計算機圖像處理技術的成熟，許多新型測量技術噴涌而現(xiàn)，最具代表性的是數(shù)字體圖像相關方法（Digital Volume Correlation，DVC）。數(shù)字體圖像相關方法是一種有效的光學測量技術，具有無損性、抗干擾能力強、測量范圍廣等多種優(yōu)勢，廣泛地應用于生物醫(yī)學、實驗力學和材料科學等多學科交叉領域。

數(shù)字體圖像相關方法通過對比變形前后物體圖像的位置變化，以獲取其內(nèi)部的變形場信息。該方法在計算過程中通常分為如下兩步，即整體素位移搜索和亞體素位移測量。整體素位移搜索可采用空域或頻域搜索以獲得物體內(nèi)部的位移，然而該方式測得的位移精度不足以達到預期理想，因此需要采用亞體素位移算法進一步提高測量精度。亞體素位移測量算法是將整體素搜索的位移初值不斷進行優(yōu)化直至達到亞體素精度，是提高測量精度的關鍵性技術，對此研究人員先后提出擬合算法、灰度梯度算法、前向累加牛頓-拉弗森（FA-NR）算法以及三維反向高斯-牛頓（IC-GN）算法等多種算法用于實現(xiàn)亞體素位移的測量?；诓煌愋偷乃惴?，其測量精度和計算效率各不相同，如何選取最佳算法尤為重要，因此定量對比多種算法的計算性能具有十分重要的指導意義。

本文介紹了在數(shù)字體圖像相關方法中最為常見的三種亞體素位移測量算法，通過計算機模擬三維散斑體圖像的剛體平移實驗，對比不同算法的測量精度和計算效率，隨后利用椎骨零應變原位掃描實驗進一步對比三種算法的測量誤差，實驗均表明IC-GN算法的性能最優(yōu)。

1 數(shù)字體圖像相關方法的基本原理

采用數(shù)字體圖像相關方法測量物體內(nèi)部變形時，需要結合三維層析設備記錄測試對象在不同狀態(tài)下的體圖像，通過匹配技術尋找這些圖像中灰度相同或相似的體素點，對比該點在變形前后位置的變化以獲取物體內(nèi)部變形場信息。具體步驟如下：首先在變形前的體圖像中選取一個尺寸為(2+1)voxel×(2+1)voxel×(2+1)voxel的正方體為參考子體塊，通過搜索尋找變形圖像中灰度最為相似的目標子體塊，并且采用相關系數(shù)（如最小平方距離函數(shù)、歸一化互相關函數(shù)等）評價兩者的相似程度。由于搜索時只能以體素為單位，該方式測得位移無法達到亞體素級別的精度，因此需要借用亞體素位移算法提高測量精度。

1.1 三元二次多項式擬合算法

數(shù)字體圖像相關法的擬合算法與二維算法類似，如圖1所示。首先通過整體素搜索獲得位移初值，選取一個以整體素點為中心、大小為3×3×3的子體塊為研究對象，采用三元二次多項式對該點周圍的幾個離散點進行插值擬合，該多項式極值點對應的坐標便是所求點的位置，進而計算出亞體素位移。

圖1 擬合點示意圖

三元二次多項式擬合函數(shù)為：

式中共有10個未知參數(shù)（～），為了求取其對應的數(shù)值，通常選取等距的12個離散點依次代入三元二次多項式中，可計算出所有未知量及多項式表達式。擬合獲取多項式的極值點應滿足下列關系：

最終計算得出的亞體素位移為：

1.2 前向累加牛頓-拉弗森（FA-NR）算法

FA-NR算法是通過逐步非線性迭代優(yōu)化整體素位移的初值直至收斂到最佳值，其實質(zhì)是將相關系數(shù)求導問題轉化為參數(shù)求最優(yōu)值問題。FA-NR算法的基本原理與擬合算法截然不同，需要用一階形函數(shù)描述子體塊的形狀變化，其定義如下：

式中：(,,)和(′,′,′)分別表示子體塊變形前后的中心點；，，為中心點在，，方向上的位移變化量；Δ，Δ，Δ為子體塊內(nèi)一點到中心點的距離；u，u，u，v，v，v，w，w，w分別表示子體塊在，，方向上的一階位移梯度。

一階形函數(shù)=[,,,u,u,u,v,v,v,w,w,w]共有12個參數(shù)，將整體素位移作為迭代的初始值，采用抗噪性能強的零均值歸一化互相關函數(shù)（ZNCC）作為相關系數(shù)對離散數(shù)據(jù)進行處理，當變形前后的圖像子體塊灰度分布最為相似時，相關系數(shù)達到極值。

零均值歸一化互相關函數(shù)（ZNCC）如下：

相關系數(shù)極值點應滿足下列關系：

采用正向迭代的表達式為：

式中：p為第次迭代值；?(p)為相關函數(shù)的一階導數(shù) 向 量，即Jacobian矩陣；??(p)為二階導數(shù)，即Hessian矩陣。控制迭代的收斂條件，當滿足 |p-p|≤0.001或達到規(guī)定迭代次數(shù)時，迭代終止。由于迭代過程中需要亞體素位置的灰度及灰度梯度信息，通常采用高精度的三元三次插值法對子體塊的灰度進行插值處理以獲得所需信息。

1.3 三維反向高斯-牛頓（IC-GN）算法

IC-GN算法的原理與FA-NR算法類似，但是迭代方向相反，IC-GN每次迭代過程中以目標子體塊為不變量，參考子體塊按照扭曲函數(shù)發(fā)生改變，其迭代參數(shù)以相乘的形式進行更新，如圖2所示。IC-GN算法常采用零均值歸一化最小平方距離相關函數(shù)（ZNSSD）用于評價子體塊變形前后相似程度的衡量標準。

圖2 IC-GN算法

零均值歸一化最小平方距離相關函數(shù)（ZNSSD）為：

與正向算法類似，IC-GN算法達到收斂條件時停止迭代。除此之外，IC-GN算法同樣也需要進行三元三次插值處理獲得灰度梯度。

2 數(shù)值模擬平移實驗

在實際應用數(shù)字體圖像相關法測量位移時會受到三維層析圖像采集系統(tǒng)、加載條件、室外環(huán)境等眾多因素的影響，因而采用數(shù)值模擬實驗的效果會比較理想。文獻[12]將二維計算機模擬散斑圖制作的方法擴展至三維可有效地模擬出三維散斑體圖包含的灰度信息，該方法分為以下三步：

1）在一定范圍區(qū)域內(nèi)隨機分布給定數(shù)量的點作為散斑的中心。

2）利用下式計算出該區(qū)域任意一點灰度：

式中：為散斑顆粒的總數(shù)；是散斑顆粒中心的光強強度且隨機分布；為散斑顆粒的大??；(x,y,z)表示散斑顆粒中心的位置。

3）將以上所有信息保存為數(shù)字圖像，便可以得到相應的散斑體圖像。

通過設置以上參數(shù)可獲得散斑數(shù)量為12 000個、半徑為3 voxel、尺寸為150 voxel×150 voxel×150 voxel的三維散斑體圖像，將該圖像連續(xù)平移10次，每次移動0.1 voxel，直到最大位移為1 voxel時停止，得到10組剛性平移的體圖像，并分別給平移前后的圖像添加均值為0、方差為2的高斯噪聲，如圖3所示。將平移后的圖像作為變形體圖像，利用三種不同的亞體素算法計算圖像的變形。

圖3 計算機制作的散斑體圖像

2.1 三種亞體素算法的精度比較

為了精準評價數(shù)字體圖像相關方法在平移實驗的測量誤差，通常將計算結果與真實施加的位移對比。

均值誤差（系統(tǒng)誤差）為：

標準差（隨機誤差）為：

式中：u為數(shù)字體圖像相關方法計算出第點的位移；為預加載的真實位移值；為計算點的個數(shù)。

使用三種方法對變形前后的體圖像進行運算時，計算子體塊的尺寸設置為21 voxel×21 voxel×21 voxel和41 voxel×41 voxel×41 voxel，計算步長均為10 voxel。不同算法對應的測量誤差如圖4所示，三種算法的均值誤差都呈現(xiàn)出周期為1 voxel、近似于正弦曲線的分布規(guī)律，這與文獻[13]中的誤差規(guī)律基本一致。位移位于0、0.5 voxel和1 voxel時，測得的位移平均值最接近于真實值，均值誤差接近于0；當施加的位移為0.2 voxel和0.8 voxel時，F(xiàn)A-NR算法和IC-GN算法測量位移的均值與真實位移相差最大，即均值誤差為最大，而擬合算法在平移0.3 voxel，0.7 voxel對應的均值誤差為最大。其次對比兩個不同的計算子體塊，位移的均值誤差變化甚小，而標準差卻成倍的減少，即計算子體塊大小的選擇對數(shù)字體圖像相關法測量位移的均值誤差影響不大，對標準差的影響卻十分顯著。

圖4 不同計算子體塊對應的測量誤差

比較三種算法的測量誤差，均值誤差方面，IC-GN算法在不同尺寸的計算子體塊中測量誤差最大為0.006 0 voxel，約為擬合算法的1 5，F(xiàn)A-NR算法測量的最大值約為其2倍，而擬合算法的測量誤差最大達到0.028 2 voxel，這說明IC-GN算法的測量精度要優(yōu)于其他兩種算法。標準差方面，IC-GN算法和FA-NR算法相近，且均小于擬合算法。在子體塊為21 voxel×21 voxel×21 voxel時，IC-GN算法和FA-NR算法測得標準差均小于0.002 voxel，而擬合算法最大為0.013 8 voxel，由此可見擬合算法的穩(wěn)定性明顯更低。當子體塊增大1倍時，三種算法的標準差減小了約4倍。

2.2 三種亞體素算法的計算效率比較

在數(shù)字體圖像相關方法中，除了要考慮計算的精度外，計算效率也是值得考慮的因素之一，隨著計算子體塊的增大，計算效率均有不同程度的降低。三種算法均在計算子體塊為21 voxel×21 voxel×21 voxel時計算效率最高，在子體塊為41 voxel×41 voxel×41 voxel時計算速度最低。其次擬合算法的計算效率最高，IC-GN算法次之，F(xiàn)A-NR算法效率最低，如圖5所示，以FA-NR的計算效率為標準。

圖5 計算效率比

在子體塊為41 voxel×41 voxel×41 voxel時，擬合算法的計算效率是IC-GN的6.41倍，是FA-NR算法的20.34倍。探究其原因可能在于：擬合算法在計算過程中，不需要進行Hessian矩陣和子體塊灰度插值的計算，且每次多項式只需計算一次即可，無需迭代，因此擬合算法具有較高的計算效率。FA-NR算法在每次迭代過程中會反復計算Hessian矩陣，且亞體素位置的灰度和灰度梯度信息需要反復利用三元三次插值獲得，這將極大消耗計算時間和內(nèi)存。而在IC-GN算法的迭代過程中，由于采用反向迭代的策略，Hessian矩陣只需要計算一次即可，迭代次數(shù)顯著減少，提高了計算效率。綜上分析，IC-GN既保證了計算精度的同時也提高了計算效率，F(xiàn)A-NR算法則是以犧牲時間而獲得較高的測量精度，而擬合算法的計算效率最高但是計算精度最低。

3 零應變原位掃描實驗

零應變原位掃描的實驗對象為正常成年男性腰椎標本，由昆明醫(yī)科大學捐贈，如圖6所示。測試時，將實驗樣品放置在Micro CT-225上進行零應變原位掃描，掃描期間樣品沒有任何移動，掃描時長約為15 min，并將掃描圖像以39μm的各向同性體素尺寸重建。CT設備掃描試件時并沒有發(fā)生移動，理論上物體沒有發(fā)生變形，然而由于多種因素的影響，如CT射線熱誤差、算法誤差、體圖像質(zhì)量等，必然導致測量結果有少許偏差，且這些誤差并不能完全消除。由于運算時采用相同的體圖像，測量結果可直接用于對比三種算法的測量誤差大小。

圖6 椎骨體圖像及其灰度分布圖

利用數(shù)字體圖像相關方法測量椎骨變形時，計算子體塊大小選擇為21 voxel×21 voxel×21 voxel，計算步長設置為5 voxel。三種算法的測量誤差如圖7所示，其中擬合算法的誤差最大，F(xiàn)A-NR算法和IC-GN算法的誤差相對較小且基本相同。

圖7 三種算法的測量均值絕對誤差和標準差

雖然算法不同，但是誤差最值出現(xiàn)的方向一致，即均值絕對誤差最大值在軸，標準差則出現(xiàn)在軸方向。在計算精度方面，IC-GN算法精度最高，其均值絕對誤差在0.057 5～0.146 9 voxel范圍內(nèi)，F(xiàn)A-NR算法的誤差為0.058 9～0.148 5 voxel，而擬合算法精度最低，其誤差范圍為0.063 1～0.159 7 voxel。在計算穩(wěn)定性方面，IC-GN算法的表現(xiàn)最優(yōu)，其標準差為0.057 5～0.146 9 voxel，F(xiàn)A-NR算法的標準差為0.059 6～0.014 85 voxel，擬合算法穩(wěn)定性最差，其標準差范圍為0.062 1～0.155 2 voxel。需要特別注意的是，三種算法在零應變實驗中的誤差最值要明顯大于平移實驗中的最值，由此可見零應變實驗的誤差主要是由外部因素引起，而算法誤差所占比例相對較小。綜合分析，F(xiàn)A-NR算法和IC-GN算法在椎骨原位掃描下測量的精度及穩(wěn)定性也均優(yōu)于擬合算法。

4 結 論

本文簡要地介紹了三種常見的亞體素位移測量算法：擬合算法、FA-NR算法和IC-GN算法，通過兩組不同的實驗對比三種算法的計算精度及其穩(wěn)定性，結果顯示IC-GN算法的性能最優(yōu)。在散斑體圖像的平移實驗中，三種算法的計算精度均隨著子體塊的增大而提高，但是計算效率卻成倍的遞減。其中IC-GN算法的測量精度最高，而擬合算法的測量誤差和穩(wěn)定性均為最差，而其計算效率最高。在椎骨零應變測試中三種算法表現(xiàn)的性能也不盡相同，IC-GN算法和FA-NR算法的測量精度均優(yōu)于擬合算法。其次對比兩組實驗可發(fā)現(xiàn)平移實驗的測量精度要高于零應變測試的精度，主要原因是外部因素引起的誤差要大于算法誤差。綜合對比計算精度和計算效率，可確定IC-GN算法為數(shù)字體圖像相關方法中的最佳亞體素位移測量算法。

注：本文通訊作者為許蔚。