基于TIDW的木材內部缺陷三維應力波成像方法*

陳方翔，馮海林，杜曉晨，方益明，翁 翔

（1.浙江農林大學信息工程學院，杭州311300；2.浙江省林業智能監測與信息技術研究重點實驗室，杭州311300）

應力波斷層成像技術在木材無損檢測中得到了廣泛應用。Ross等［1］最早利用應力波檢測技術對紅橡樹腐朽區域進行檢測。通過應力波成像軟件觀察所獲得的圖像能得到木材內部腐朽位置、面積和腐朽程度［2］。如王立海等［3］利用應力波與X射線二維CT圖像的結合診斷能高效準確地確定原木內部腐朽區域；戚大偉等［4］對X射線圖進行濾波、銳化，從而使得木材內部缺陷的細節更加明顯；Wang等［5］利用應力波成像技術采用8個傳感器對原木進行初期腐朽檢測，實驗結果表明該技術檢測木材內部缺陷準確率達到62%，其中有8.5%的健康木材被誤檢為腐朽材。梁善慶等［6］通過應力波斷層成像技術獲取木材二維斷層圖像，通過斷層圖像獲取木材內部缺陷的大小、形狀等信息。楊學春等［7］利用應力波法對原木內部腐朽檢測的相關內容進行研究，結果表明應力波測試儀能得到原木內部腐朽基本形狀的二維圖像，并判斷不同樹種內部的嚴重腐朽情況。

這些研究大部分都是利用現有的商業化成像軟件進行二維成像應用研究，對應力波成像算法本身的研究還比較少見。馮海林等［8］提出一種圖像重建算法，利用周圍點值估計未知網格點的速度值，測試結果證明了該方法的可行性。Choi等［9］利用基于模型的損傷檢測算法來定位缺陷位置及腐朽程度。木材三維應力波成像對算法本身的研究也較少，研究人員常常利用無損檢測技術獲取木材二維斷層圖進行三維研究。Brancheriau等［10］利用超聲波成像技術，獲取木材內部結構。王在山等［11］利用線性插值法對木材CT圖像進行三維重建，從而提高木材利用效率。

反距離加權算法（Inverse Distance Weighted）最早用于氣象和地質成像研究。劉廣明等［12］利用光譜指數法與IDW法對三維土壤鹽分空間變異特征進行了解析與評價，提升了區域土壤鹽分的三維預測精度。Wang等［13］利用IDW方法獲取長江及黃河流域氣象數據，得到人類活動與氣象的關系。陳冬花等［14］分別利用IDW法與Kringing法分析氣溫與海拔的關系，實驗結果表明在低溫環境下IDW模型模擬效果比Kringing模型略好。該算法根據空間分布規律對未知區域進行推演，具有較好的實驗測試效果，并且其運算速度快、使用范圍廣、計算時所需存儲空間小等特征，在氣象、地理成像研究中較為常見。

將IDW算法應用于木材內部缺陷應力波三維成像會存在如下問題：①預估點的計算結果與鄰域內所有已知點均有關，增加了算法復雜性；②預估點的計算結果與鄰域內所有已知點之間的距離有關，忽視了插值點的空間結構性。因此，本文提出了一種TIDW（Top-kInverse Distance Weighted）算法，該算法將預估點鄰域關系擴展到三維空間，增加預估點的搜索半徑，找出鄰域內與其距離小于搜索半徑的所有已知點，引入Top-k查詢，查找搜索半徑內與預估點影響最大的k個已知點，從而計算得到預估點的值并進行應力波三維成像。利用TIDW算法與IDW算法對樣本1至樣本5進行三維應力波成像實驗對比，結果表明TIDW成像算法具有更高的成像精度和檢測準確率。

1 反距離加權算法

空間插值是一種利用已知點屬性求得未知點屬性的方法，通過空間中已知點與未知點的空間結構關系，從而計算得出未知區域中點的屬性。

反距離加權算法是空間插值方法之一，常用于地理空間插值［15］。該方法利用預估點與已知點之間的距離為權重，離預估點越近的已知點權重越大，其權重大小與距離成反比。設空間預估點為pi（xi，yi）（i=1，2，3，…，n），其鄰域內有已知點qi（xi，yi）（i=1，2，3，…，n），其中Zqi為鄰域內已知點的屬性值。待預估點Zpi屬性值是通過待預估點鄰域內已知點qi（xi，yi）的屬性值Zqi加權平均值計算得出，其權的大小與待插點和鄰域內點之間的距離di有關，是一種與距離的倒數成反比關系的插值方法。其基本公式為［16］：

其中Zpi為待預估點的估算值，Zi為鄰域內第i個已知點的屬性值，di為第i個已知點與待預估點之間的距離，n為參與插值計算鄰域內已知點的個數。為權重系數，m越大則說明越靠近預估點的已知點對插值結果的影響越大。

2 TIDW成像算法

本文改進了IDW算法，提出了TIDW（Top-kIn?verse Distance Weighted）的應力波成像算法，對預估點與其鄰域內已知點的空間結構性進行了優化，增加了預估點的搜索半徑，引入Top-k查詢技術，提高了預估點的計算精度。

2.1 鄰域空間擴展方法

算法根據木材周圍傳感器采集到的應力波傳播速度數據，得到鄰域內已知點的速度集，并通過已知點數據集，計算得到預估點的屬性。在計算預估點屬性時需對其鄰域內已知點進行篩選，本文中的篩選方法參考了空間定位［17］和方位搜索法，其中方位搜索法一般包括四方搜索、六方搜索、八方搜索等，方位的劃分越多，則選擇的參估點越多。四方搜索法，是根據預估點的橫、縱坐標把平面分成四個象限，在每一個象限中查找與預估點距離最近的已知樣本點。文中將二維空間中的四方搜索法擴展至三維空間中，當預估點位置確定后，利用其搜索半徑查找相關已知點。若空間區域為N，鄰域內的已知點X屬于N，預估點根據搜索半徑篩選X中的已知點。預估點計算方法如下：

①設預估點pi（xi，yi）周圍8個鄰域空間內存在個已知點qi（xi，yi）（i=1，2，…，M，M＜X），且qi的屬性為Zqi。

②每個鄰域空間中qi的選取應滿足搜索半徑r，r滿足r＜R（R為木材的半徑），r的取值應保證預估點pi（xi，yi）在空間鄰域內存在已知點qi與之相關，并設閾值δ，當搜索到已知點的數目小于該閾值δ時，可擴大搜索半徑，直到搜索到的已知點數達到該閾值為止。

③令λpi表示預估點pi（xi，yi）到其鄰域內有關點的權重，則λpi可以表示為：

其中，表示點pi（xi，yi）和qi（xi，yi）之間的距離倒數，m為適合的常數，通常m=1。可將式（1）轉化為：

2.2 基于Top-k查詢方法

Top-k查詢技術主要用于查詢大量數據中相關性最大的k個結果［18-19］，TIDW算法利用Top-k查詢技術，找出在預估點鄰域范圍內與其影響最大的k個已知點及其屬性。Top-k查詢定義如下：給定M個元組的集合T，每個元組具有m′=（u1，u2，…，um）個屬性，將集合T，存儲為列文件的集合S={S1，…，Sm}，每個列文件為二元組合Si（rid，ui），其中rid表示對象的標示符，ui表示對象在屬性處的屬性值。每個列文件的存儲方式為各元組的屬性值的單調非增序列。定義F為m′個屬性的評分函數，其式如下［20］：

λi是評分函數定義在屬性值ui上的權重。通常，F是單調函數，即 ?a1,?a2∈T，如果對所有1≤i≤m′，a1u1≤a2u2，那么。利用Top-k查詢技術查詢集合T，中各元組的k個子集，通過讀取m列已經降序排列的列文件集合S，順序讀取序列中，當元組rid出現時，隨機讀的方式在另外一個m-1個列文件獲取其他屬性值，然后計算他們的評分值，如果這個值是目前見過最大的k個，用優先隊列維護k個元組及其相關信息。對每一列序列，設其當前讀取位置ui，設閾值τ=F（u1，u2，…，um），當優先隊列里k個元組分數值的最小值不小于τ時，查詢結束。將式（3）變形如下：

Zqi為已知點的屬性值，λpi為定義在屬性值上的權重，為最終屬性值。根據式（4），將式（5）變形如下：

Zqi為已知點的屬性值，為定義在屬性值上的權重，F(pi)為評分函數。

2.3 TIDW算法成像步驟

步驟1：輸入傳感器采集到的應力波速度數據；

步驟2：輸入傳感器坐標X（x，y，z），通過采集的時間矩陣Tm×n，計算出線速度Vm×n矩陣；

步驟3：設預估點X為2 000，計算出鄰域內已知點集M′，計算權重 λpi；

步驟4：設定r與閾值δ根據式（2）計算權重，保存為序列形式；

步驟5：利用式（5）與已知點集M′求Zpi，令i=1，轉步驟③；

步驟6：令 λi=λpi；帶入式（6）；

步驟8：根據式（6）中的屬性值F（pi），將預估點根據不同屬性值進行顏色賦值，并三維可視化。

3 實驗分析與結果

本實驗采用項目組自主研發的便攜式木材斷層成像設備，在對木材進行檢測時，需要在被測木材的不同高度處隨機安裝傳感器，每次使用12個傳感器進行數據采集，實驗時依次敲擊各個傳感器，對同一樣本需進行2次數據采集，共采集24組數據，將數據用于三維成像分析。圖1為敲擊12個傳感器采集樣本實驗數據圖。

圖1 成像設備及檢測方式

3.1 實驗材料與方法

本實驗樣本編號、樣本種類、樣本周長、樣本測量高度、樣本含水率等信息如表1所示。

表1 實驗樣本基本信息

本文實驗采用意大利Klortner公司生產的KT-R重錘式木材測濕儀來測量含水率。在樣本缺陷內注入適量的水，直到水量達到測量高度時停止，取出缺陷內水并將其倒入量杯測量缺陷內水的體積V′得到樣本缺陷體積，樣本體積V是由樣本橫截面的面積乘以測量高度計算得到，公式如下：

圖2為樣本1至樣本6的實物圖，本實驗利用人工挖鑿空洞與自然腐朽的孔洞進行測試效果對比。其中樣本1、3、6中的孔洞為人工挖鑿，用于模擬真實腐朽情況；樣本2、4、5中的孔洞為自然腐朽。

圖2 6個實驗樣本木樁

3.2 IDW算法與TIDW算法成像結果對比

將兩種算法分別對不同缺陷類型樣本進行成像實驗對比，實驗結果如表2至表6所示。成像結果中以深色點組成的區域來模擬樣本中實際腐朽的區域，淺色點組成的區域來模擬樣本中健康的區域。

表2為對缺陷相對規則并且缺陷為工人挖鑿的樣本進行二維成像測試，由測試結果可知，IDW算法的成像結果并不能顯示缺陷位置，而TIDW算法的成像效果較好，并能反映出缺陷的基本位置與缺陷情況。

表3為對樣本1-2進行三維成像測試。樣本缺陷體積相對均勻，并由表中實驗結果可知，利用IDW算法的成像結果并不能顯示缺陷的腐朽程度，而TIDW算法三維成像效果較好，能反映出缺陷的位置與缺陷的腐朽程度。

表4為對缺陷不規則并且為自然腐朽的樣本進行二維成像測試。由測試結果可知，IDW算法成像結果并不能顯示缺陷位置，而TIDW算具有較好的成像效果，并能反映出缺陷的基本位置與缺陷情況。

表5為對樣本3至樣本5進行三維成像實驗，樣本3和樣本4的缺陷體積相對均勻，而樣本5的缺陷不均勻。由表中實驗結果可知，利用IDW算法成像效果并不能顯示缺陷的腐朽程度，而TIDW算法的三維成像效好，并能反映出缺陷的位置與腐朽程度。

表2 缺陷規則的樣本二維測試結果

表3 缺陷規則的樣本三維測試結果

表4 缺陷不規則的樣本二維測試結果

表6為對缺陷相對規則并且存在2個工挖鑿的空洞樣本進行二維與三維成像測試。由測試結果可知，IDW算法成像結果并不能顯示缺陷位置與缺陷腐朽程度，而TIDW算法成像效果相對較好，能反映出缺陷的基本位置，但不能很好的反映出缺陷的腐朽程度。

3.3 成像結果分析

表7中顯示缺陷的實測量體積、算法計算得到的缺陷體積以及相對誤差，用于驗證檢測結果的正確性。通過計算圖中紅色缺陷部分點的比率乘以實測樣本體積，可得到缺陷部分體積，如下：

其中Vi為通過算法計算得到缺陷部分的體積，V為對應樣本體積，n′為缺陷點個數，X為插值點總數。

V1為利用IDW算法計算所得體積，V2為利用TIDW算法計算所得體積。Dt1為IDW算法計算得到的缺陷體積與實測缺陷體積的相對誤差，Dt2為TIDW算法計算得到的缺陷體積與實測缺陷體積的相對誤差，其公式如下：

表5 缺陷不規則的樣本三維測試結果

表6 樣本中存在多個缺陷的測試結果

表7 樣本缺陷檢測結果

本文中1至6號實驗樣本的實測缺陷體積、計算得到缺陷體積及相對誤差信息如表7所示。其中相對誤差Dt1的范圍在14.07%～422.2%之間，相對誤差Dt2的范圍在0.18%～52.56%之間，樣本1至樣本5的相對誤差Dt2小于相對誤差Dt1，并且相對誤差Dt2在16.1%以內，而樣本6中的相對誤差Dt2大于相對誤差Dt1。實驗結果顯示樣本1，2，3，4的相對誤差Dt2在10%以內，表明利用TIDW算法成像的精度較高；相對誤差，說明利用IDW算法成像效果并不理想，與樣本實際缺陷體積有較大誤差。樣本5的相對誤差Dt2為16.1%，其誤差偏大的原因可能是由于樣本5的上下端腐朽程度相差較大，腐朽不均勻，使得TIDW算法成像結果與實際情況存在一定誤差，即利用算法計算所得的缺陷體積與實際缺陷體積有較大誤差；樣本6的相對誤差Dt1為14.07%，相對誤差Dt1小于相對誤差Dt2，但由表6中實驗圖像結果可知，IDW算法成像結果中的缺陷位置與實際缺陷位置不同，缺陷的腐朽程度與實際樣本的缺陷腐朽程度也不同；樣本6的相對誤差Dt2為52.56%，其原因可能是樣本中存在2個人工開鑿的孔洞，利用TIDW算法對此類樣本進行成像時會存在較大誤差，并結合由表6中成像結果可知，TIDW算法成像中的缺陷位置基本符合實際缺陷位置，但其腐朽程度與實際缺陷腐朽程度不同，因此兩種算法對此類樣本成像的準確率都不高。由樣本1至樣本6的檢測結果可知，利用IDW算法檢測樣本中缺陷體積，其平均檢測準確率僅為38.26%，而利用TIDW算法檢測樣本中缺陷體積，其平均檢測準確率為85.24%。

4 結論

本文提出了一種基于TIDW的木材內部缺陷三維應力波成像方法，應用于應力波木材無損檢測，并對木材內部缺陷進行三維成像研究。通過與IDW算法成像結果比較，該算法無論從定性還是定量上都能更好反映木材內部腐朽情況，具有較高的成像精度。利用傳感器采集木材周圍應力波傳播速度，通過對不同樣本缺陷情況進行實驗，驗證了該算法對木材內部缺陷三維成像的可行性。檢測單個空洞的缺陷樣本，算法成像準確率可達到85.24%以上，而針對多孔洞的樣本實驗檢測準確率不高，能基本反映樣本缺陷位置，但缺陷體積計算有一定誤差，需要進一步研究。

［1］Ross R J，Ward J C，Tenwolde A.Identifying bacterially infected oak by stress wave nondestructive evaluation［J］.USDA Forest Service，Forest Products Laboratory，1962（30）：277-281.

［2］Gilbert A E，Smiley E T.Picus Sonic Tomography for the quantifi?cation of decay in white oak（Quercus alba）and Hickory（Carya spp.）［J］.Journal of Arboriculture，2004，30（5）：277-280.

［3］王立海，孫天用.基于應力波與X射線二維CT圖像原木內部腐朽無損檢測［J］.森林工程，2011，27（6）：26-29.

［4］戚大偉.木材無損檢測圖像處理系統的研究［J］.林業科學，2001，37（6）：92-96.

［5］Wang X P，Wiedenbeck J，Ross R J，et al.Nondestructive Evalua?tion of Incipient Decay in Hardwood Logs［D］.Gen.Tech.Rep.FPL-GTR-162.Madison，WI：U.S.Department of Agriculture，For?est Service，Forest Products Laboratory，2005.

［6］梁善慶，傅峰，胡娜娜，等.三種應力波斷層成像診斷木材內部缺陷［A］.第二屆中國林業學術大會，S11木材及生物質資源高效增值利用與木材安全論文集［C］.北京：中國林學會，2009，457-462.

［7］楊學春，王立海.紅松木材結構缺陷對應力波傳播參數的影響［J］.東北林業大學學報，2005，33（1）：30-31.

［8］Feng H L，Li G H，Fu S，et al.Tomographic image reconstruction using an interpolation method for tree decay detection［J］.Biore?sources，2014，9（2）：3248-3263.

［9］Choi F C，Li J，Samali B，et al.Application of modal-based dam?age-detection method to locate and evaluate damage in timber beams［J］.Journal of Wood Science，2007，53（5）：394-400.

［10］Brancheriau L，Saadat-Nia M A，Gallet P，et al.Ultrasonic Imag?ing of Reaction Wood in Standing Trees［J］.Acoustical Imaging，2012，31（3）：399-411.

［11］王在山，戚大偉.木材CT圖像的三維重構［J］.森林工程，2014，30（6）：38-40.

［12］劉廣明，吳亞坤，楊勁松等.基于光譜指數的區域土壤鹽分三維空間變異研究［J］.光譜學與光譜分析，2013，33（10）：2758-2761.

［13］Wang Y，Ding Y J，Ye B s，et al.Contributions of climate and hu?man activities to changes in runoff of the Yellow and Yangtze riv?ers from 1950 to 2008［J］.Science China Earth Science，2013，56（8）：1398-1412.

［14］SONG J J，KWON S，LEE G W.Incorporation of parameter uncer?tainty into spatial interpolation using Bayesian trans-Gaussian kriging［J］.Advances in Atmospheric Sciences，2015，32（3）：413-423.

［15］Song J J，KWON S，LEE G W.Incorporation of parameter uncer?tainty into spatial interpolation using Bayesian trans-Gaussian kriging［J］.Advances in Atmospheric Sciences，2015，32（3）：413-423.

［16］范玉潔，余新曉，張紅霞，等.降雨資料Kriging與IDW插值對比分析-以漓江流域為例［J］.水文，2014，34（6）：61-65.

［17］張晉升，李一博，王偉魁，等.基于TFDA的有限空間多目標聲定位方法［J］.傳感技術學報，2013，26（11）：1508-1512.

［18］聶云峰，王長勝，陳崇毅.一種空間查詢高效的無線傳感網絡路由協議［J］.傳感技術學報，2015，28（5）：745-751.

［19］鄭瑾，賈維嘉，王國軍.無線傳感器網絡中基于數據分布表的Top-k查詢協議［J］.傳感技術學報，2010，23（9）：1340-1346.

［20］Shaikh S A，Kitagawa H.Top-koutlier detection from uncertain data［J］.International Journal of Automation and Computing，2014，11（2）：128-142.