自然环境下的树木易受各种不利环境因素影响或因自身生理原因而出现衰退现象。树木受侵害产生衰退早期，一般从内部开始空腐，无法从树木外部获得准确的监控和测定。探地雷达是一种无损探测技术，对探测特体本身无任何破坏作用，相较于目前在林业领域主要应用的其他木质体无损探测技术，如应力波、超声波和X射线扫描法等，具有完全非损伤性、无环境污染和探测速度快等优点; 但是由于树木内部结构的复杂性(心材、边材和树皮及结疤、开裂、腐朽和中空等缺陷)以及含水率和环境导致的个体差异性，使得雷达波对于木质体内缺陷的检测还存在解析困难、探测精度不高、可视化程度低等不足，这些问题严重影响和制约雷达波在树木领域的应用。国内外学者对此展开了广泛研究，Devaru等(2005)阐述了探地雷达在树木无损检测方面的优点，并设计试验探究了探地雷达对树木躯干内部缺陷的自动检测算法；Butnor等(2009)研究发现探地雷达在松柏类、针叶类树木中对其树皮表面孔洞、内部有空气孔洞的检测情况效果很好，并采用阈值法对针叶树木内部缺陷进行了探测，面积误差在18.9%左右；Lorenzo等(2010)研究了雷达波对树干以及树根检测的可行性；Halabe等(2009)通过介电常数、含水率的精确测量，结合阈值法对木材厂的木材结疤以及内嵌铁钉深度进行检测，误差率在10%左右；Riggio等(2015)和Perez-Gracia等(2014)分别对建筑木房梁上的腐朽、损坏探测进行了定位研究；Ježová等(2016)通过分析规则型类圆柱体内部多次反射波的传输规律，研究了树木躯干成像的可能性；甘明旭等(2016)采用树木雷达对黄帝陵古侧柏(Platycladus orientalis)树干空腐进行探测识别，主要在雷达波仪器的使用以及腐朽的判定；邸向辉等(2013)探讨了探地雷达在木材无损检测领域的理论可行性，指出电磁波信号解释分析的困难性是影响树木无损探测的主要制约因素；陈勇平等(2017a；2017b)采用人工试验方法利用雷达对马尾松(Pinus massoniana)木材内部孔洞探测的定量检测进行了可行性分析。上述研究定性分析了各种因素对木材内部缺陷定位的影响，但并没有给出较为精准的定位方法。目前对木材内部缺陷的检测主要集中在应力波检测，岳小泉等(2016)采用应力波对二维检测腐朽缺陷面积进行估算，误差率在17.78%~52.61%之间，而对于树木躯干雷达波无损检测尚处于仪器应用和试验研究阶段，主要是对树木内部缺陷的定性分析，而对树木内部缺陷的定量分析如定位和分布成像的相关研究较少。为了实现雷达波对树木内部缺陷的定位和准确表征，本文研究了树木内部层位分层定位算法，采用雷达波探测层位分析技术与三维激光扫描(薛金林等，2014)相结合的方法，从而实现了对树木的无损检测、内部缺陷的定位以及二维图像可视化。

1 材料与方法 1.1 雷达波扫描图像获取

采用美国TreeWin公司树木雷达系统TRU，其检测工作主要包含2个独立步骤：一是野外数据采集；二是离线数据分析软件模块对数据进行分析。TRU系统主要包括2个组成部分：美国GSSI公司的SIR系列数据采集管理器和雷达介质耦合天线(900 MHz)。数据采集管理器(图 1a)用于雷达波数据数字化波形显示及测量数据的存储；雷达天线(图 1b)与主机互为一体，天线绕树外围每移动5 mm就产生1次信号对树木进行测量。

本研究以北京颐和园的柳树(Salix)为试验数据采集对象(2016年4月20日)。在试验探测过程中，首先标记所测截面的高度(图 2中横线表示所测截面的高度，竖线表示检测的起始点)，然后在所测截面进行匀速360°的圆周扫描，数据自动存储到野外数据采集器中。

1.2 外形轮廓获取

大多数树木为不规则柱体结构，为了实现树木内部状况的准确定位，可以用近似圆柱形和不规则柱体的方法获取树木外部轮廓。圆柱形外形在软件分析中，其截面可采用极坐标转换的方式近似为圆形，而不规则柱体的树木外形，可通过三维激光扫描仪获取其空间三维坐标来对其外形轮廓进行追踪。三维激光扫描仪内部有激光发射和激光探测装置，2个步进电机控制激光光束在水平和垂直方向的移动。扫描仪发射一束激光，该光束遇到障碍物发生反射，扫描仪中的激光探测器探测到反射光，则完成一个扫描点的测量。通过测量激光发射与接收之间的时间差T，计算出激光传播距离，即可获得空间物体表面的三维坐标。树木通过激光扫描测量方法，可以快速获取被测目标表面的三维坐标数据，从而得到树木外形轮廓数据。本研究采用自行开发的三维激光扫描系统进行树木轮廓定位。

1.3 雷达波原理

探地雷达探测树木的基本原理是高频电磁波以脉冲形式通过发射天线被定向地送入树木体内。木材本身的介电常数与较多因素有关，其中含水率与介电性质直接相关。木材内部各部位的含水率是不同的，导致木材各个层位的介电常数也有差异，当木材内部出现腐朽、空洞等缺陷时，缺陷部分含水率和介电常数与周围介质有较大差异。当电磁波在含有缺陷(腐烂、空洞)的树木内传播时，其路径、强度将发生变化。探地雷达波在树木内不同介质中传播形成的反射和折射不同，遇到存在介电性质差异的腐朽或空洞时，电磁波便发生强反射使雷达反射波振幅变化较大。因此在对雷达波回波数据进行处理和分析的基础上，根据接收到的雷达波波形幅度、时间等并采用阈值法推断树木内缺陷的空间位置、结构、介电性质和几何形态，可实现对树木内部缺陷的分层探测识别。电磁波在介质中传播时，其路径和波形将随所通过介质的电性质和几何形态发生变化，故可通过对电磁波回波信号的处理，分析回波时间和波形幅度，来推断树木内缺陷位置和层介质的结构(Solla et al., 2013)。

本研究建立了雷达波在结构层中传播的正演模型(Lambot et al., 2004)，通过得到的反射波模拟验证雷达波在结构层中的传播规律(图 3)。第1层的相对介电常数和厚度分别为5和20 cm，第2层的相对介电常数和厚度分别为13和15 cm，第3层的相对介电常数和厚度分别为3和10 cm。

正演模型选择的初始波为rikerz波，900 MHz(图 4a)。电磁波垂直入射时，根据反射系数的菲涅耳(Fresnel)公式(Halabe et al., 2009)可知：1)界面两侧介质的电磁学性质差异与反射波强度呈正相关，且从反射振幅上可以判定两侧介质的性质和属性；2)波从介电常数小的介质进入介电常数大的介质，即从高速介质进入低速介质、从光疏介质进入光密介质时，反射系数为负，即反射波振幅反向；反之，从低速介质进入高速介质时，反射波振幅与入射波同向。因此，雷达反射波的振幅和方向是雷达波判别的重要依据。已知正演模型的各个参数、每个层位的厚度(d)以及对应的介电常数(ε)，根据式(2)可以推算电磁波在每层中传播的速度，根据式(3)可以推知每个层位的厚度。通过正演模型推算出来的厚度与实际设定的层位厚度对比及菲涅尔公式的结论，可以推知层位的大致位置。如图 3，第1/2层之间的波形与第2/3层之间的波形相反(图 4b)，进而验证了菲涅尔公式的判断依据。

对于脉冲型GPR系统而言，接收到的回波信号可以近似认为是各层界面反射波的叠加(Baili et al., 2009)：

$ {y_{\rm{r}}}\left( t \right) = \sum\limits_{i = 0}^N {{A_i}x\left( {t - {\tau _i}} \right)} + n\left( t \right)。$

(1)

式中：x(t)为入射GPR脉冲；N为层界面数；A_i为每层界面的反射波幅值；n(t)为所加的高斯白噪声；τ_i为第i层的双程回波时间。

反射波幅值(A_i)可以运用振幅比的方法求解介电常数，回波时延可用于确定雷达波在介质层内传播的时间，进而可以根据这些参数估计出每层的厚度。电磁波在介质中的传播速度(v_p)计算公式如下：

$ {v_{\rm{p}}} = c/\sqrt {\varepsilon '} ; $

(2)

厚度(Z)计算公式如下：

$ Z = \frac{{c \times t}}{{2 \times \sqrt {\varepsilon '} }}。$

(3)

式中：ε′为相对介电常数；t为电磁波双程的时间。

由于树木内部结构的复杂性及个体化差异性，无法采用先验知识获取树木内部的电磁特性，需要研究在线无损参数反演算法。本研究采用层波反演方法估计各层介质的介电常数。以图 3所示木材典型的3层介质为例，A₁为第1层与第2层界面间的反射波幅度，A₂为第2层与第3层界面间的反射波幅度。则第1层介质介电常数的估计公式为：

$ {\varepsilon _1} = {\left( {\frac{{1 + {A_1}/{A_{\rm{m}}}}}{{1 - {A_1}/{A_{\rm{m}}}}}} \right)^2}。$

(4)

式中：A_m为发射波的幅值，实际测试中发射波难以获得，可以采用雷达波对金属板的全反射进行反向获得的回波幅值。

第2层介质介电常数的估计公式为：

$ {\varepsilon _2} = {\varepsilon _1}{\left[ {\frac{{1 - {{\left( {{A_1}/{A_{\rm{m}}}} \right)}^2} + \left( {{A_2}/{A_{\rm{m}}}} \right)}}{{1 - {{\left( {{A_1}/{A_{\rm{m}}}} \right)}^2} - \left( {{A_2}/{A_{\rm{m}}}} \right)}}} \right]^2}。$

(5)

在式(5)中，第1层介质假定为均匀的，该层电介质在某种意义上是一个常数，且在其厚度内不发生变化，因此该公式也即假定雷达信号在表层内的传播是不衰减的(Loizos et al., 2007)。根据层剥反演获得的相对介电常数可以对树木内部的传输速度进行在线估测，分析雷达信号在不同层位反射信号的时间间隔即分层处理，从而实现树木层位的识别。

1.4 雷达波对树木分层处理方法

在雷达波反射信号中提取不同层位的反射波形，需要研究雷达波信号分层处理方法。时域雷达波信号分层算法主要包括阈值法、匹配滤波器法和希尔伯特积算法。

1.4.1 阈值法

阈值法是最常用的层介面检测方法(Lahouar，2003)，其通过比较所分析信号和固定先验阈值之间的大小来判断分目标信号是否存在。假定最小可检测的反射信号(通常可以选择作为一个合适的阈值)高于噪声水平，则最佳阈值是指阈值的设置保证单道信号采样点中信号误判率最低。阈值V_t(Lahouar，2003)选取根据公式：

$ {V_{\rm{t}}} = \delta \sqrt { - 2\lg {P_{\rm{f}}}} 。$

(6)

式中：δ为噪声偏差；P_f为虚警率，本文使用的雷达仪器，单道波为512个采样点。由于实际采集中的噪声是未知的，因此一般使用单道波的尾部作为噪声信号。

GPR信号中的回波脉冲通常是多峰的，单个反射脉冲有多个局部峰值点，对缺陷点定位分层产生了干扰。信号包络可以解决多重峰的问题，包络信号通常是在通信系统中使用以提取调制包络和拒绝高频载波信号的。本文使用的GPR信号虽然属于无载波脉冲体制，但是依然可以提取出信号幅度变化趋势的包络：

$ {x_{\rm{e}}}\left( t \right) = \left| {{x_{\rm{a}}}\left( t \right)} \right| = \left| {x\left( t \right) + j\mathop {x\left( t \right)}\limits^ \wedge } \right|。$

(7)

式中：$ \mathop {x\left( t \right)}\limits^ \wedge $为回波信号x(t)的希尔伯特变换；x_a(t)为分析信号；x_e(t)为信号的包络。

根据以上分析，阈值法检测步骤大致为：1)提取探地雷达回波信号的包络信号；2)搜索包络信号的极值点；3)确定信号的阈值；4)将包络峰值与阈值点相比较，若峰值大于阈值点则判定为检测到反射层，并计算相对应的时间延迟；反之，则没有检测到反射层。

采用阈值法对正演模型中的回波单道信号进行定位分层(图 5)，平行于时间轴的实线为阈值线。图 5中超过阈值线的有3处，通过计算3处时间延迟并与正演模型中层位时间延迟进行对比，可知阈值法可以定位出层位。

1.4.2 匹配滤波器法

匹配滤波器(Lahouar，2003; Mucciardi，2002)可以使输出信噪比最大。匹配滤波器的脉冲响应可以表示为：

$ h\left( t \right) = x\left( {T - t} \right)。$

(8)

式中：T为信号x(t)持续的时间。

匹配滤波器法中，在时间T内，当信号x(t)超过某个阈值S_t时，匹配滤波器将会输出最大值或最小值(取决于被检测信号与原始信号的极性)。依据N-P准则，匹配滤波器的阈值S_t取决于容许的最大虚警率P_f(Lahouar，2003)：

$ {S_{\rm{t}}} = {\rm{erf}}{{\rm{c}}^{ - 1}}\left( {{P_{\rm{f}}}} \right)\sqrt {{\sigma ^2}E} 。$

(9)

式中：E为已知信号x(t)的能量；erfc^-1(P_f)为虚警率P_f的逆互补误差函数。

匹配滤波器法可归结为以下步骤：1)利用式(8)计算出匹配滤波器对应的脉冲响应h(t)；2)利用h(t)对检测到的回波信号y_r(t)进行匹配滤波，获得滤波信号y_MF(t)；3)找到y_MF(t)绝对值的最大值y_{max MF}以及最大值所对应的时间延迟t_max；4)比较y_{max MF}与S_t之间的大小；如果y_{max MF} > S_t，记录相对应的时间延迟t_max，然后减去在带宽内的初始信号x(t)，从最大的时间延迟t_max可以推出每个层位的时间延迟；5)将被检测到的探测波段在y_r(t)内置为0；6)重复步骤2。

采用匹配滤波器法时，正演模型中单道信号的定位分层见图 6。图 6中4幅子图从上到下依次为模型的单道回波信号图以及从回波信号中提取出来的第2层、第2层、第3层的反射信号。

1.4.3 希尔伯特积算法

Huang等(2008)提出了一种新的时频分析方法——希尔伯特-黄变换(Hilbert-Huang transfrom，HHT)，该方法对于非平稳、非线性信号的分析比较直观，且自适应强。在这个理论中，通过经验模态分解(empirical mode decomposition，EMD)将信号自适应分解成有限多个内在模分量(intrinsic mode function，IMF)和1个表征信号趋势变化的残余信号，并且提出对得到的各个IMF运用希尔伯特变换进行时频分析。

根据Rosefeld的子带乘积理论，如果噪声信号可以分解成多频子带分量，则每个分量的逐点乘积运算可以凸显有用信号而减弱噪声。基于这一思想，Kim等(2003)提出了一种低信噪比下检测生物医学信号的小波检测器。受到小波检测器的启发，本文使用IMF分量来构造检测算法。希尔伯特积算法主要包括以下4个步骤：

1) 对回波信号进行EMD分解，得到相应的N个IMF分量c₁(t)~c_N(t)。

2) 计算所有IMF分量绝对值的点积，记为：

$ P\left( t \right),P\left( t \right) = \left| {\prod\limits_{i = 1}^N {{c_i}\left( t \right)} } \right|。$

(10)

3) 因为回波信号包含正峰和负峰，使得接收信号由多峰脉冲组成，易造成虚假检测，因此需要对信号进行窗函数的平滑处理。用P(t)和窗函数W(t)进行卷积来完成平滑处理：

$ T\left( t \right) = W\left( t \right) * P\left( t \right)。$

(11)

判断T(t)和阈值V_t的大小，T(t)信号中大于阈值的判定为检测到反射层，并计算相对应的时间延迟；反之，则没有检测到反射层。

运用希尔伯特积算法，对模型正演单道波进行分析处理，见图 7。图 7中呈现出的3个波峰，根据时间延迟计算，就是正演模型中的第1层、第2层和第3层。该方法验证通过希尔伯特积算法，可以对不同的介质进行定位分层。

2 结果与分析 2.1 分层处理算法比较

为了验证分层算法的准确性，本研究采用楔形正演模型对3种算法进行对比。运用蒙特卡洛法随机选取36个不同的点，对每个点的距离进行层位误差分析。如图 8，通过FDTD算法正演模拟rikerz波进入楔形正演模型来比较阈值法、匹配滤波器法和希尔伯特积算法的误差率，雷达波从右到左、从上到下对正演模型进行探测扫描。共获得36组试验数据，对每组数据进行误差分析。通过运用3种算法，预估出每个正演模型第1层和第2层之间的距离dx，用其与实际距离dy之间的误差来计算3种层位识别算法的精确度。实际距离dy为正演模型设定，算法的误差计算公式为：

$ {\rm{error\% = }}\frac{{{\rm{d}}y - {\rm{d}}x}}{{{\rm{d}}y}} \times 100\% 。$

(12)

经过误差比较(图 9)，可知希尔伯特积算法相较于其他2种算法误差率最低。

2.2 试验数据分析

电磁波在介质中传播时，其路径和波形将随所通过介质的介电特性和几何位置而变化。正常树质(树皮、边材和心材)与缺陷(腐朽、空洞和虫洞)的介电常数差异较大，反射波幅值与相邻两界面间的介电常数差值呈正比。雷达反射波的时延主要取决于木质体内部介电特性以及传播距离，因此可以通过分析回波信号的时延和幅度来推测介质层的介电常数和传播距离。

探地雷达多采用雷达波反射法，一般依据步进轮控制同一截面等间距(本研究采用5 mm间距)，根据获得的单道A扫描数据使用5 mm间距排列获得笛卡尔坐标系下的B扫描数据。由于树木多为不规则柱状结构，这种扫描结果导致解读上的困惑性，易造成树木内部图像成像效果较差。针对这一问题，本研究采用三维激光扫描树木三维图形数据，通过获取的相同高度界面的外形轮廓点云信息，采用图像追踪方法获取树木外形轮廓曲线与5 mm等间距的测试位置坐标，将B扫描数据和A扫描数据与外轮廓曲线的中测试点法线进行坐标映射，进而获得树木内部雷达波成像的绝对定位信息。

为了判断测试方法的准确性，选择颐和园柳树为试验对象。选取3个试件样本T1、T2、T3作为具体试验分析对像(图 10b)，3个样本树龄都在30年左右，T1含水率40.31%~89.87%, T2含水率35.83%~53.67%, T3:含水率42.51%~85.42%。

雷达波数据经过常规预处理后，使用希尔伯特积检测算法对树木内部缺陷(腐朽或空洞)进行定位分层。图 10a为3棵树木的B扫描图，图中绿色线为树干/空气的位置，红色点组成的线为树干内缺陷部分，可能为空洞或腐朽的位置。A′、B′、C′、D′、E′离树皮的距离由式(3)计算(图 11)，在所测截面使用阻力钻对腐朽的实际位置进行测量，试验数据与实际测量数据(取中间值)进行对比。试验中对A′、B′、C′、D′、E′多次测量，计算其深度误差率。

B扫描图通过雷达天线在树木某一高度进行360°环绕探测获取，实际树木缺陷面积通过网格法计算，试验结果与Treewin分析结果进行对比分析。Treewin采用单一介电常数结合阈值法对雷达回波数据进行定位，结合极坐标转换生成树木内部缺陷的二维图像(图 10c)，通常与实际树木异相差较大，缺陷区域面积误差率在20%左右。本研究采用希尔伯特积算法对树木进行分层，层剥反演介电常数实现树木内缺陷位置定位，结合三维激光对树木外形轮廓的数据，利用轮廓追踪和坐标变换方法生成不规则树木的二维缺陷分层图(图 10d)。令S_j为检测出的腐朽缺陷面积，S_z为实际缺陷的面积，通过对比试件图像或检测图像，用网格法来确定S_j和S_z。可得面积误差率S%为：

$ S{\rm{\% = }}\frac{{\left| {{S_{\rm{j}}} - {S_{\rm{z}}}} \right|}}{{{S_{\rm{z}}}}} \times 100\% 。$

(13)

从表 1中可以看出，采用希尔伯特积算法以及三维扫描数据和轮廓追踪法，可以将对树木内缺陷定位的深度误差控制在10%左右，二维缺陷分层图显示的缺陷面积与对所测量截面缺陷的面积误差控制在5%左右，树木内缺陷区域计算的精确度有很大提高。

3 结论

探地雷达在树木无损检测领域的应用受限于树木内部结构复杂、个体性差异化大以及不规则树木外形轮廓导致的绝对定位困难；同时雷达波在木质体媒介中的传输规律和传输特性受树木内部复杂结构影响，雷达波解析解译困难，精度有待提高。本研究在分析树木内部分层界面识别算法的基础上，通过比较阈值法、匹配滤波器法和希尔伯特积算法在树木内部层位探测的准确性，根据雷达波在树木不同介质以及缺陷分界面反射波的振幅和时延，实现了树木内部结构和缺陷的分层定位，得出希尔伯特积算法对树木内缺陷区域的检测更加准确，采用层剥反演在缺陷定位的基础上结合三维激光扫描和外形轮廓追踪技术实现了不规则树木轮廓的绝对定位和分布成像。对3种算法分别使用FTDT正演方法进行对比验证，并应用于颐和园柳木试件试验测试中，结果表明对于树木内部缺陷区域可以实现定位和表征，单点径向距离和缺陷面积误差可分别控制在10%和5%左右。本研究提出的方法可实现树木内部缺陷的准确定位和分布成像，同时在线反演方法对不同木质体具有良好的适应性。