基于射線分割的林木應力波斷層成像算法*

劉 濤 李光輝

(江南大學人工智能與計算機學院 無錫 214122)

應力波斷層成像技術(shù)是一種新的無損檢測方法，在林木無損檢測領(lǐng)域得到了廣泛應用。與其他無損檢測技術(shù)(如X射線、超聲波和CT圖像)相比，應力波斷層成像技術(shù)有其特殊優(yōu)勢，如設備便攜、安全可靠，且適合各種環(huán)境條件等(王立海等， 2002； 楊學春等， 2002)。林木應力波無損檢測的基本思想是利用脈沖錘敲擊分布在林木周圍的傳感器，在林木內(nèi)部產(chǎn)生應力波傳播，根據(jù)傳感器記錄的應力波傳播時間，分析應力波在林木截面中的速度分布，實現(xiàn)對林木內(nèi)部缺陷的快速定位和檢測。通常，應力波在腐朽木質(zhì)中的傳播速度慢于健康木質(zhì)，當遇到裂縫、空洞等缺陷時，應力波會跨過該區(qū)域，導致傳播速度明顯降低，從而可有效區(qū)分出健康木質(zhì)和腐朽木質(zhì)(Dackermannetal.， 2014)。目前，國內(nèi)外關(guān)于應力波技術(shù)的研究已有大量成果(梁善慶等， 2008； 王立海等， 2011； 翁翔等， 2016)，Ross等(1992)最早對患有真菌感染的橡樹(Quercus)進行試驗，驗證了應力波技術(shù)檢測林木缺陷的可行性。楊學春等(2005)基于動力學理論對原木中的應力波傳播進行研究，建立了位移、速度、應力和應變方程。張厚江等(2011)與美國農(nóng)業(yè)部林產(chǎn)品實驗室合作，使用微鉆阻力儀和應力波實現(xiàn)了古建筑木構(gòu)件材料力學性能的快速檢測和評估。梁善慶等(2010)基于應力波斷層成像技術(shù)獲取木材斷層圖像，能夠?qū)δ静娜毕葸M行有效診斷。

由于木材內(nèi)部應力波傳播速度因材質(zhì)變化而發(fā)生改變，因此現(xiàn)有應力波斷層成像算法通常僅根據(jù)波速差異判斷材質(zhì)健康狀況，成像軟件的設計開發(fā)過程并未考慮具體的木材含水率、密度和缺陷類型等因素(安源， 2013； 杜曉晨， 2019)。一些學者專門研究了上述因素對應力波傳播速度的影響，如徐華東等(2011)分析含水率對紅松(Pinuskoraiensis)木材應力波傳播速度的影響，結(jié)果發(fā)現(xiàn)應力波傳播速度隨含水率增加而下降； 劉昊等(2014)研究密度對木材應力波傳播速度的影響，結(jié)果發(fā)現(xiàn)同一樹種內(nèi)，應力波傳播速度與木材密度存在正相關(guān)線性關(guān)系，不同樹種之間，應力波傳播速度與木材密度沒有明顯關(guān)系； Li等(2011)分析應力波傳播速度與木材硬度之間的映射關(guān)系，以評估應力波斷層成像技術(shù)對樹木腐爛的檢測能力，結(jié)果表明應力波斷層圖可以顯示出林木內(nèi)部缺陷，但與木材硬度圖不完全匹配。以往研究成果對加深應力波在木材中傳播規(guī)律的理解、設計更好的應力波斷層成像算法具有重要指導意義。

目前，關(guān)于應力波斷層成像算法的研究相對較少。Divos等(2005)指出，應力波傳感器數(shù)量和所應用的反演技術(shù)是影響斷層成像精度的重要因素。Feng等(2014)提出一種圖像重建算法，利用已知點速度估算未知網(wǎng)格的速度，并通過試驗證明了算法的可行性。Du等(2015)提出一種基于橢圓插值的斷層成像方法，根據(jù)應力波傳播射線相對應的橢圓區(qū)域估算網(wǎng)格單元的速度，提高了空間插值的精度。焦治等(2018)采用最小二乘法反演計算網(wǎng)格單元的速度分布，并使用誤差校正機制優(yōu)化斷層圖像，其檢測精度高于PICUS檢測儀。Liu等(2018)基于傳統(tǒng)的應力波反演原理，結(jié)合曲線傳播模型提出一種混合傳播模型(hybrid wave propagation model，HWPM)，能夠較好反映林木橫截面的缺陷情況。

為了進一步提高應力波斷層成像精度，本研究提出一種基于射線分割的林木應力波斷層成像算法(tomography imaging algorithm based on ray segmentation for stress wave of wood，簡記為RSIA)：首先，生成用于描述林木橫截面的網(wǎng)格圖； 然后，對成像區(qū)域內(nèi)每個網(wǎng)格單元進行應力波速度反演估算； 最后，通過網(wǎng)格圖反映林木橫截面應力波速度的空間分布。整個網(wǎng)格圖可以看作一個矩陣，由于應力波傳播射線數(shù)量有限，所以該速度矩陣通常是稀疏的，如何正確估算每個網(wǎng)格單元的速度是成像精度的關(guān)鍵。RSIA基于射線分割思想，將一條射線分割成多條線段，精細化利用每一條應力波傳播射線，能夠提高初始網(wǎng)格單元速度的準確性和相關(guān)性。利用HWPM和RSIA對多株樣本進行斷層成像對比發(fā)現(xiàn)，RSIA具有更高的成像精度。

1 應力波斷層成像原理

1.1 應力波在林木橫截面的傳播原理 假設林木橫截面為規(guī)則的圓形，在該截面均勻布置一定數(shù)量的傳感器獲取應力波時間數(shù)據(jù)。圖1中，S表示信號發(fā)射端，R1和R2表示信號接收端，應力波由信號發(fā)射端傳播至信號接收端。SR1為應力波沿弦線方向的傳播路徑，VC為路徑SR1的應力波傳播速度；SR2為應力波沿直徑方向的傳播路徑，VD為路徑SR2的應力波傳播速度；β為弦線方向與徑向方向的夾角(rad)。

圖1 林木橫截面應力波傳播示意Fig.1 Schematic diagram of stress wave propagation in tree cross sectionS： 信號發(fā)射端Signal transmitting end； R1、R2： 信號接收端Signal receiving end； SR1： 應力波沿弦線方向的傳播路徑The propagation path of the stress wave along the chord line direction； VC： 路徑SR1的應力波傳播速度The stress wave propagation velocity in the path SR1； SR2： 應力波沿直徑方向的傳播路徑The propagation path of the stress wave along the diameter direction； VD： 路徑SR2的應力波傳播速度The stress wave propagation velocity in the path SR2； β： 弦線方向與徑向方向的夾角The angle between the chord line direction and the radial direction(rad).

由于林木內(nèi)部為非均質(zhì)各向異性體，因此應力波在林木內(nèi)部的傳播速度并不均勻。以林木橫截面為例，應力波在林木邊緣的傳播速度小于沿徑向的傳播速度，而林木缺陷是根據(jù)應力波在林木內(nèi)部的傳播速度分布進行可視化成像，林木邊緣會因速度較小判定為缺陷區(qū)域，造成誤判。為了解決這一問題，Li等(2014)通過大量試驗分析推導出應力波在健康林木中的傳播模型，并基于該模型校正不同路徑上的應力波速度：

VC≈(1-0.2β2)×VD。

(1)

如圖2a為一株健康松木樣本的截面圖，圖2b為該樣本的原始應力波傳播射線圖，圖2c為經(jīng)過速度校正后的應力波傳播射線圖。

如圖3a為一株含有人工缺陷的松木樣本的截面圖，圖3b為該樣本的原始應力波傳播射線圖，圖3c為經(jīng)過速度校正后的應力波傳播射線圖。

在圖2和圖3中，綠色線段代表在該傳播路徑上不存在缺陷，紅色線段代表在該傳播路徑上存在缺陷，可見，應力波速度必須根據(jù)其傳播方向經(jīng)過校正后才能反映林木內(nèi)部的真實健康狀況。

圖2 健康樣本中應力波速度校正前后的射線圖對比Fig.2 The ray diagram comparison before and after stress wave velocity correction in healthy sample

圖3 缺陷樣本中應力波速度校正前后的射線圖對比Fig.3 The ray diagram comparison before and after stress wave velocity correction in defective sample

1.2 基于射線分割的林木應力波斷層成像算法 假設應力波在林木橫截面中沿直線傳播，橫截面周圍均勻布置多個傳感器，依次敲擊傳感器獲取任意2個傳感器之間的傳播時間，由于傳感器之間的距離是已知的，因此可以得到每2個傳感器之間的應力波傳播速度。各傳感器之間的應力波傳播速度可用以下矩陣表示：

(2)

式中：n為傳感器數(shù)量； NaN為空值；Vij為第i號傳感器與第j號傳感器之間的應力波傳播速度；Dij為第i號傳感器與第j號傳感器之間的距離；Tij為第i號傳感器與第j號傳感器之間的應力波傳播時間，其中1≤i，j≤n，當i=j時，表示相同編號的傳感器。

根據(jù)應力波傳播速度，繪制應力波傳播射線圖，如圖4所示。綠色射線速度較快，表示該路徑上材質(zhì)健康；黃色射線速度較慢，表示該路徑上存在輕微腐朽；紅色射線速度很慢，表示該路徑上存在嚴重腐朽或空洞。黃色和紅色射線僅能表示該路徑上存在腐朽，但并不能確定腐朽區(qū)域的具體位置，腐朽可能位于射線中央附近，也可能位于射線兩端附近。在計算成像區(qū)域內(nèi)網(wǎng)格單元的速度時，不能使用整條射線的平均速度表示單元內(nèi)的速度。為了解決這一問題，本研究提出一種基于射線分割的林木應力波斷層成像算法，該算法的基本思想是利用每條射線與其他射線的交點，將其分割成多條線段，根據(jù)已知射線傳播速度估算待分割射線上多條線段的速度，從而得到更為精細的速度分布，提高初始網(wǎng)格單元速度的準確性。

圖4 基于射線分割的林木應力波斷層成像算法示意Fig.4 Schematic diagram of stress wave tomography imaging algorithm based on ray segmentation for nondestructive testing of wood

首先，校正采集的應力波速度，并對校正后的速度進行歸一化和離散化，選取合適閾值劃分林木中的健康區(qū)域、輕微腐朽區(qū)域和嚴重腐朽區(qū)域。由于不同樹種木材特性(如密度和含水率等)不同，應力波在不同樹種內(nèi)的傳播速度必然存在差異，因此選取的閾值也不相同。設閾值為θ1、θ2，區(qū)間[0,θ1]表示嚴重腐朽區(qū)域，區(qū)間[θ1,θ2]表示輕微腐朽區(qū)域，區(qū)間[θ2,1]表示健康區(qū)域。如圖5中有空洞缺陷的紅松樣本，含水率為24.3%，經(jīng)速度校正后，應力波傳播速度(m·s-1)區(qū)間為[1 174,1 980]，將其歸一化至[0,1]，經(jīng)多次試驗得到θ1為0.21、θ2為0.5。

根據(jù)相鄰兩交點的速度重新估算該線段的速度，線段的速度可表示為：

Sv=(start+end)/2。

(3)

式中： start和end分別表示兩交點的速度；Sv表示線段的速度。

循環(huán)執(zhí)行該操作直至所有射線均經(jīng)過上述分割處理，得到改進后的應力波傳播射線圖，如圖5所示。

接著，估算成像區(qū)域內(nèi)每個網(wǎng)格單元的速度，判斷某一點是否在成像區(qū)域內(nèi)按下式進行：

(4)

式中：r為被檢測林木橫截面半徑；x和y分別代表網(wǎng)格單元的橫縱坐標，若f(x,y)=1，則表示該網(wǎng)格在成像區(qū)域內(nèi)； 若f(x,y)=0，則表示該網(wǎng)格在成像區(qū)域外。

當網(wǎng)格單元有應力波射線穿過時，該網(wǎng)格單元的速度估計值為：

圖5 射線分割結(jié)果對比Fig.5 Comparison of ray segmentation results

(5)

式中：Vxy為成像區(qū)域內(nèi)坐標(x，y)的網(wǎng)格單元速度估計值；Vi為穿過該網(wǎng)格單元的射線速度；N為穿過該網(wǎng)格單元的射線數(shù)量。

實際上，由于應力波傳播射線較少，網(wǎng)格單元較多，會出現(xiàn)部分網(wǎng)格單元未被射線穿過的情況，導致無法計算應力波速度，此時選用最近鄰插值法進行數(shù)據(jù)補全，該網(wǎng)格單元的速度估計值為：

Gxy=w1v1+w2v2+w3v3+w4v4。

(6)

式中：Gxy為無應力波傳播射線的網(wǎng)格速度；v1、v2、v3、v4為距該網(wǎng)格歐氏距離最近的4個非空網(wǎng)格單元的速度；w1、w2、w3、w4為對應速度的權(quán)重。

最后，通過分析每個網(wǎng)格單元的速度判定該網(wǎng)格狀態(tài)，確定截面中的健康區(qū)域、輕微腐朽區(qū)域和嚴重腐朽區(qū)域，并進行可視化成像。

綜上所述，基于射線分割的林木應力波斷層成像算法步驟如下：

步驟1，利用FAKOPP應力波檢測儀采集樣本橫截面應力波數(shù)據(jù)，根據(jù)式(1)校正原始應力波速度，并對校正后的速度進行歸一化和離散化。

步驟2，劃分網(wǎng)格單元，根據(jù)校正后的應力波速度繪制傳播射線圖。

步驟3，分割任意一條射線，根據(jù)端點速度估算該射線上多條線段的速度。

步驟4，重復步驟3，直到所有射線均經(jīng)過分割處理，將分割后的射線可視化，得到改進的應力波傳播射線圖。

步驟5，計算成像區(qū)域內(nèi)網(wǎng)格單元的速度，并對無數(shù)據(jù)的網(wǎng)格單元通過最近鄰插值法進行數(shù)據(jù)補全，確保所有網(wǎng)格單元均被處理。

步驟6，根據(jù)網(wǎng)格單元的速度判定網(wǎng)格狀態(tài)，結(jié)合圖像處理方法，實現(xiàn)林木內(nèi)部缺陷的圖像重建。

2 樣本與數(shù)據(jù)采集

為評估基于射線分割的林木應力波斷層成像算法的有效性，本研究選取4株實驗室原木樣本和4株江南大學校內(nèi)活樹樣本進行試驗。

2.1 原木樣本 原木樣本為1株健康雪松(Cedrusdeodara)、1株含人工缺陷的雪松、1株含人工缺陷的烏桕(Sapiumsebiferum)和1株含人工缺陷的紅松。原木樣本基本信息如表1所示，其中樣本高度為原木圓盤的厚度，測量高度為被檢測截面在原木圓盤的高度，如圖6所示。

2.2 活樹樣本 在江南大學校園內(nèi)選取1株柳樹(Salixspp.)、1株香樟(Cinnamomumcamphora)和2株圓柏(Sabinachinensis)共4株活樹樣本，因活樹截面情況無法直觀觀察到，故使用Resistograph微鉆阻力儀評估活樹內(nèi)部健康狀況。活樹基本信息如表2所示。

表1 原木樣本基本信息Tab.1 General information of log sample

表2 活樹基本信息Tab.2 General information of live trees

3 結(jié)果與分析

本研究選取4株原木樣本和4株活樹樣本進行應力波斷層成像試驗，利用FAKOPP應力波檢測儀采集樣本橫截面應力波數(shù)據(jù)，并與HWPM進行對比評估RSIA的效果?；趯嶒炇以緲颖具M行試驗，能夠直觀觀察到樣本內(nèi)部缺陷，便于現(xiàn)場檢驗RSIA的檢測效果。通過實驗室原木樣本試驗證明算法有效性后，再將其應用于活立木內(nèi)部缺陷檢測，由于活立木內(nèi)部缺陷無法直觀觀察到，故使用Resistograph微鉆阻力儀評估RSIA的準確性。

3.1 原木樣本成像結(jié)果分析 圖7a為一株健康雪松樣本。2種算法的斷層成像結(jié)果如圖7b和圖7c所示，均檢測出該雪松樣本是健康的，沒有誤判。

圖7d為一株含人工缺陷的紅松樣本，缺陷位置位于截面右側(cè)。2種算法的斷層成像結(jié)果如圖7e和圖7f所示，雖然均能檢測到樣本內(nèi)部的空洞和腐朽區(qū)域，但是HWPM斷層圖像中的空洞和腐朽區(qū)域明顯大于實際情況，誤判較多； 而RSIA斷層圖像中缺陷位置和大小更加符合實際情況。

圖7 原木樣本斷層成像結(jié)果Fig.7 The tomographic results of log samples

圖7g為一株含人工缺陷的烏桕樣本，缺陷位置位于截面底部。2種算法的斷層成像結(jié)果如圖7h和圖7i所示，HWPM未檢測到樣本內(nèi)部的空洞區(qū)域，僅能判斷出輕微腐朽區(qū)域，且腐朽區(qū)域的形狀和大小也有較大偏差； 而RSIA能夠準確判斷出空洞和腐朽區(qū)域的位置和形狀，不過缺陷區(qū)域略小于實際情況，這是由于樹干不規(guī)則導致采集的應力波數(shù)據(jù)不準確造成的。

圖7j為一株含人工缺陷的雪松樣本，缺陷位置位于截面右側(cè)。2種算法的斷層成像結(jié)果如圖7k和圖7l所示，均能檢測到樣本內(nèi)部的空洞和腐朽區(qū)域，HWPM斷層圖像中缺陷區(qū)域略大于實際情況； 而RSIA斷層圖像中缺陷的位置、形狀和大小更加符合實際情況。

對于原木樣本，通過計算網(wǎng)格單元數(shù)量統(tǒng)計缺陷面積，面積超過網(wǎng)格單元的50%記1，小于網(wǎng)格單元的50%記0，原木缺陷面積如表3所示。

表3 原木缺陷面積Tab.3 Defect area of log sample

由表3可知，RSIA較HWPM在缺陷面積上與實際缺陷面積大小更加接近。

綜上，HWPM檢測含缺陷的原木樣本時，斷層圖像中缺陷面積往往大于實際情況，且形狀與實際情況存在偏差； 而RSIA通過分割原始射線，提高初始網(wǎng)格單元速度的準確性，得到了更為精細的速度分布，可更準確地還原樣本內(nèi)部情況。對比試驗結(jié)果可知，RSIA的平均成像精度有了較大提高，在缺陷位置、大小和形狀方面更加符合實際情況。

3.2 活樹樣本成像結(jié)果分析 Resistograph阻力曲線能夠準確反映指定路徑上的內(nèi)部情況，通過多路徑鉆探，可以大致推斷出缺陷的大小和位置。

圖8a為一株內(nèi)部具有缺陷的圓柏樣本，2種算法的斷層成像結(jié)果如圖8b和圖8c所示。圖9a為使用微鉆阻力儀沿路徑1-7鉆探得到的阻力曲線，其中，區(qū)域Ⅰ為缺陷部分，區(qū)域Ⅱ為健康部分，2種算法均未檢測到缺陷區(qū)域Ⅰ所示健康區(qū)域。圖9b為使用微鉆阻力儀沿路徑4-10鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分，區(qū)域Ⅱ中出現(xiàn)部分阻力急劇上升的情況，這是樹干內(nèi)部存在節(jié)子導致的。在2種算法的斷層圖像中，均顯示出大面積紅色和黃色區(qū)域，表明該圓柏樣本中存在大面積芯材腐朽。

圖8d為一株表面具有明顯裂縫的圓柏樣本，2種算法的斷層成像結(jié)果如圖8e和圖8f所示。圖9c為使用微鉆阻力儀沿路徑6-12鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分，2種算法均未檢測到該路徑中的缺陷，這是因為該缺陷部分面積小且腐朽不嚴重，導致應力波數(shù)據(jù)表現(xiàn)不明顯。圖9d為使用微鉆阻力儀沿路徑7-12鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分。圖9e為使用微鉆阻力儀沿路徑10-5鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分。由于林木表面具有一條較大裂縫，應力波傳播過程中需繞過該部分，導致應力波速度下降。在2種算法的斷層圖像中，均含有從樹干表面延伸至內(nèi)部的黃色和紅色區(qū)域，與實際情況相符。在路徑7-12和路徑10-5的缺陷長度上，RSIA的斷層圖像與阻力曲線更加吻合。

圖8g為一株表面具有較小裂縫的香樟樣本，2種算法的斷層成像結(jié)果如圖8h和圖8i所示。圖9f為使用微鉆阻力儀沿路徑1-7鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分。圖9g為使用微鉆阻力儀沿路徑3-10鉆探得到的阻力曲線，該路徑均健康，無缺陷部分。圖9h為使用微鉆阻力儀沿路徑6-2鉆探得到的阻力曲線，其中，區(qū)域Ⅰ和Ⅲ為健康部分，區(qū)域Ⅱ為缺陷部分。圖9f和圖9h均存在極小的缺陷部分，一般是因蟲蛀導致的。在2種算法的斷層圖像中，均得到1號傳感器和2號傳感器之間存在黃色和紅色區(qū)域，與實際存在的裂縫位置相吻合。

圖8j為一株表面存在明顯真菌感染的柳樹樣本，2種算法的斷層成像結(jié)果如圖8k和圖8l所示。圖9i為使用微鉆阻力儀沿路徑3-7鉆探得到的阻力曲線，其中，區(qū)域Ⅰ為缺陷部分，區(qū)域Ⅱ為健康部分。圖9j為使用微鉆阻力儀沿路徑5-9鉆探得到的阻力曲線，其中，區(qū)域Ⅰ為缺陷部分，區(qū)域Ⅱ為健康部分。由阻力曲線可以看出，該樣本僅在樹干表面存在極小缺陷。在2種算法的斷層圖像中，均顯示樹干邊緣存在部分黃色區(qū)域，表明樹干表面存在輕微腐朽。

在斷層圖像中，虛線Pi-j表示微鉆阻力儀的鉆探路徑。為了定量評估RSIA的有效性，比較斷層圖像與實際阻力曲線在同一路徑下的缺陷長度，通過對多路徑下的數(shù)據(jù)進行比較，反映斷層成像算法在缺陷位置、大小等方面的準確性。

表4中，用l來表示使用Resistograph微鉆阻力儀沿鉆探路徑出現(xiàn)的缺陷長度，通過多路徑測量l可以反映出該截面中的大概缺陷情況。l0為阻力曲線中的缺陷長度，l1、l2分別為HWPM和RSIA算法圖中計算得到的l，Δl表示l1或l2與l0之間的相對誤差，Δl的定義如下：

(7)

式中：li代表l1或l2，Δl越小，算法越準確。

綜上，在活樹應力波斷層成像試驗中，RSIA基于射線分割思想，精細化利用每一條傳播射線，提高初始網(wǎng)格單元速度的準確性，得到了更為精細的速度分布。在圖像重建時，較HWPM有著更高的精確度，同時RSIA在斷層圖像中可以更好地還原缺陷情況。

圖8 活樹斷層成像結(jié)果Fig.8 The tomographic results of live trees

圖9 活樹阻力曲線Fig.9 Resistance curves of live trees

4 結(jié)論

本研究提出一種基于射線分割的林木應力波斷層成像算法(RSIA)，該算法通過校正原始應力波傳播速度，繪制應力波傳播射線圖并劃分網(wǎng)格單元，將每一條應力波傳播射線進行分割，重新估算該射線上多條線段的速度，得到改進的應力波傳播射線圖，提高初始網(wǎng)格單元速度的準確性和相關(guān)性，從而得到更加精細的速度分布，結(jié)合圖像處理方法，實現(xiàn)了對樣本截面的高精度成像。選取4株原木樣本和4株活樹樣本進行應力波斷層成像試驗，證實了該算法的有效性，且在檢測缺陷的位置、大小等方面準確度高于HWPM； 但RSIA對微小的節(jié)子、蟲蛀等缺陷不能有效檢測，今后還需開展大量試驗并完善算法，進一步提高斷層成像精度。