Preliminary analysis on the summer habitat selection of Mongolia subspecies of goitred gazelle in Wulate, Inner Mongolia of northern China
-
摘要:目的研究鹅喉羚蒙古亚种与夏季栖息地的关系，为探索鹅喉羚蒙古亚种的生境选择和生境恢复提供科学依据。方法于2016年7—8月，采用样线法和随机样方法，在内蒙古乌拉特后旗的蒙古野驴−梭梭林国家级自然保护区展开了鹅喉羚蒙古亚种夏季生境选择研究。共获得71个鹅喉羚数据样方，测定了15种生境因子，计算Vanderploeg和Scavia选择指数同时进行主成分分析。结果鹅喉羚喜好高平原和残丘，坡向半阴，灌木类型以白刺灌木林为主。倾向选择灌木数量为0 ~ 10株、20 ~ 30株，植被盖度为0 ~ 20%、60% ~ 80%，距道路的距离为10 ~ 15 km，距居民点的距离为10 ~ 15 km的生境。主成分分析结果显示，前5个主成分的累计贡献率为76.504%，其中第1主成分的贡献率为22.058%，地形因子（坡向、坡位和坡度）是第1主成分的影响因子，其他几个主成分的生境选择因子主要是干扰、水源和食物及植被因子。结论因此，我们认为夏季鹅喉羚对生境的选择呈现一定的偏好和回避，是由于夏季保护区内气候、地形、食物以及干扰等因子特有的季节性特征导致，也体现出其对保护区实际生境情况和资源分布进行的择优和取舍行为。Abstract:ObjectiveDue to the lack of research on goitred gazelle (Gazella subgutturosa hillieriana) at home and abroad, it has some adverse effects on its protection. Therefore, this paper studies the habitat selection of Mongolia subspecies of goitred gazelle in summer, and will provide a basis for exploring habitat selection and habitat restoration of the species.MethodWe investigated the habitat selection of Mongolia subspecies of goitred gazelle in the Ass-Saxoul National Nature Reserve of Wulate, Inner Mongolia of northern China from July to August 2016, using methods of transects and random samples. Totally 71 quadrats with gazelle and its tracks were obtained and 15 habitat selection factors were measured. The Vanderploeg and Scavia selection index and principal component analysis were used to determine the main habitat factors.ResultThe results showed that goitred gazelle chose high plain and torso mountain, semi-shady slope, shrub with nitraria sibirica (Nitraria tangutorum); shrub numbers of no more than 10 or 20–30; the vegetation coverage of less than 20% or 60%–80%; and with the distance from the road and resident sites 10–15 km. Principal component analysis (PCA) indicated that the total cumulative contribution rate for 5 components was 76.504%, the contribution rate of the first principal component was 22.058%, in which the highest was terrain factors (slope orientation, slope position and slope degree); high contribution rate factors in other principal components were disturbance, water distance and vegetation factors.ConclusionTherefore, we suppose that gazelles ’ preference and avoidance for habitat selection is due to the unique seasonal characteristics of climate, topograghy, food and disturbance in summer and reflects its preference and trade-off for the actual habitat and resource distribution in the protected area.
-
无损检测对于发现林木内部缺陷，防止腐朽扩散，提高木材利用率具有重要意义。应力波无损检测技术具有成本低，操作方便，无辐射，不受被测木材的大小形状影响等优点，已经在林木无损检测中得到了广泛的应用。应力波无损检测技术通过测量应力波在树木横截面中多个传播方向的传播时间，然后反演计算应力波在内部单元的传播速度，并生成木材的断层图像，从而准确、直观地检测到木材内部缺陷大小与位置[1]。国内外关于木材应力波无损检测技术的研究已有大量成果。Dikrallah等[2]通过实验分析了湿材的声学各向异性，研究了应力波速度与方向角之间的数学关系。梁善庆等[3]与美国农业部林产品实验室合作，将弹性波断层成像技术运用到活树的检测，结果表明该技术能够直观显示活树内部腐朽位置、大小、形状和严重程度。杨学春等[4]对林木内部的应力波传播规律进行了研究，建立了相关的波动方程及其数值求解模型。Li等[5]研究了健康树木横截面内部的应力波传播速度变化规律，建立了数学模型，利用美国黑樱桃树(Prunus serotina)的无损检测实验结果验证了该模型的有效性。Socco等[6]提出了一种较为复杂的数据处理方法用于对活树的健康监测，并在地震波图像算法的基础上研究了一种反演算法用来生成树干的应力波传播时间图像，该算法对一些小的缺陷(如小空洞)不敏感，而且没有很好地考虑木材的各向异性问题。Du等[7]在考虑木材各向异性的基础上，提出了一种改进型的椭圆空间插值方法，并利用速度补偿获得更精确的空间插值数据，该方法能够有效抵抗由于缺陷区域密度变化引起的信号干扰。目前，国际市场上已有不少林木应力波断层成像检测设备(例如，匈牙利产的FAKOPP，德国产的PiCUS Sonic Tomograph和Arbotom)，均能够生成木材横截面的二维或三维断层图像，但成像效果和检测精度仍然有待改善[8]。
为了进一步提高林木应力波断层图像的准确性，提出了一种应力波速度反演的误差校正机制(ECM)。在此基础上，利用最小二乘QR分解法(LSQR)设计了一种新的林木应力波无损检测断层成像算法(ECIA)，并且利用大量的真实样本实验验证了该算法的有效性。
1. 应力波速度模型及反演过程
介绍应力波在健康树木中的传播速度模型及无损检测基本思想，并描述林木无损检测断层成像的速度反演过程。
1.1 应力波在健康树木中的传播速度模型
为简单起见，我们将树干视为一个圆柱体，其横截面为正圆形，如图 1所示，其中S代表应力波发射端传感器的位置，R1和R2代表应力波接收端传感器的位置。在检测过程中，通过脉冲锤敲击发射端，产生应力波，并从发射端传播到接收端。检测仪测量发射端与接收端传感器之间的应力波的传播时间。SR2代表应力波沿直径的传播路径，其传播速度称为径向传播速度，用VR表示。SR1代表应力波沿与直径成β角的弦的传播路径，其传播速度称为弦向传播速度，用VT表示[9]。通常在健康树木中，应力波的径向传播速度最快，并且随着弦向角β的增大，传播速度会变慢。Li等[5]通过研究建立了应力波在健康树木中的速度模型：
VT/VR≈1−0.2β2 (1) ![]()
图 1 树木横截面内应力波传播示意图S代表应力波发射端传感器的位置，R1和R2代表应力波接收端传感器的位置；VR、VT分别表示径向传播速度和弦向传播速度；β表示弦向角。Figure 1. Stress wave propagation in the cross-section of tree trunkS is a source sensor, and R1 and R2 are receiver sensors. VT is the tangential velocity and VR is the radial velocity. β is the angle between radial direction and tangential direction.式中：VT/VR表示弦向传播速度和径向传播速度的比值，β表示弦向角。
树木产生腐朽或空洞往往从内部开始，树干外层材质相对而言更健康，所以实际检测过程中，经常使用各相邻传感器之间的应力波传播速度平均值作为健康参考值，并根据公式(1)计算横截面内任意方向的参考速度Vref[10]。对任意的传播路径，将实际测量得到的传播速度Vtest和参考传播速度Vref进行比较，可以初步判断该路径是否存在缺陷。根据国内外学者的研究经验[1]，若Vtest和Vref的比值小于0.9，可以认为该路径存在缺陷。
1.2 速度反演原理及LSQR迭代算法
速度反演是指根据检测仪测量得到的不同传播路径上的应力波平均速度，反演计算应力波在横截面内不同单元的速度分布的过程[11]。目前大多数的木材断层成像反演算法都是基于地震波传播的模拟算法[12]。
如图 2所示，反演算法需将树干横截面划分为若干网格单元，且每个网格单元内需要至少有一条路径穿过。从S到R为一条应力波传播路径，该路径的传播时间t可以用慢度p(速度的倒数)从起点S到终点R的线积分表示：
t=∫RSpdl (2) ![]()
图 2 树木横截面网格单元划分示意图S、R分别代表应力波发射端传感器和应力波接收端传感器的位置；1, 2, …, M代表网格单元编号；P1, P2, P3, P4代表传播路径与网格单元的交点；dij代表第i条路径被第j个网格单元分割的距离。Figure 2. Grid cells of a cross-section of a tree trunkS is a source sensor, and R is a receiver sensor. 1, 2, …, M represent the No. of grid cells. P1, P2, P3, P4 mean the intersections of stress wave paths and grid cells. dij is segment intercepted by the ith wave path on the jth cell.式中：dl表示路径长度的增量，p是应力波的传播慢度(速度的倒数)。
为了计算慢度分布，把图 2所示的横截面划分为M个网格单元，从左上角第一个网格单元开始，自左上至右下依次编号，将连续函数p(x, y)近似成每个网格单元内离散的慢度分布。假定每个网格单元内的慢度不变，将式(2)离散表示为[13]：
ti=M∑j=1pjdij(i=1,2,⋯,N) (3) 式中：dij表示第i条路径被第j个网格单元分割的长度，N表示总的应力波传播路径数。式(3)可用矩阵的形式表示：
\mathit{\boldsymbol{T}}{\rm{ = }}\mathit{\boldsymbol{DP}} (4) 式中：T和P分别代表N×1的传播时间列向量和M×1的传播慢度列向量，D代表传播路径被网格单元分割的长度距离矩阵，其大小为N×M。
对于给定的T和D，只需要计算出慢度列向量P后就可以生成初步的断层图像。由于一条传播路径只会被少数几个网格单元分割，所以对于任意路径i在大多数网格单元中，被分割的路径长度dij的值是0，因此式(4)中的系数矩阵D通常是一个大型稀疏矩阵。LSQR方法是一种利用Lanczos迭代法求解最小二乘问题的方法，由于其计算量小，运用矩阵的稀疏性简化计算，因此十分适合求解大型稀疏矩阵问题[14]。
由于应力波检测仪的传感器数量有限，应力波传播路径总数N通常小于网格单元总数M，故矩阵D的行数小于列数，式(4)为欠定方程组[15]。欠定方程组有多个解，为了求欠定方程组的特解，LSQR需要先将系数矩阵D对角化[14, 16]。
\left\{ \begin{array}{l} \begin{array}{*{20}{l}} {\mathit{\boldsymbol{D}} = \mathit{\boldsymbol{UB}}{\mathit{\boldsymbol{V}}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{U}}_k} = \left[ {{u_1}, {u_2}, \cdots , {u_k}} \right]}\\ {{\mathit{\boldsymbol{V}}_k} = \left[ {{v_1}, {v_2}, \cdots , {v_k}} \right]} \end{array}\\ {\mathit{\boldsymbol{B}}_k} = \left[ {\begin{array}{*{20}{c}} {{\alpha _1}}&{}&{}&{}\\ {{\beta _2}}&{{\alpha _2}}&{}&{}\\ {}&{{\beta _3}}& \ddots &{}\\ {}&{}& \ddots &{{\alpha _k}}\\ {}&{}&{}&{{\beta _{k + 1}}} \end{array}} \right] \end{array} \right. (5) 通过逐列对比，得到迭代关系式，利用此迭代方法实现矩阵D的下双对角分解。
\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{U}}_{k + 1}}\left( {{\beta _1}{e_1}} \right) = b}\\ {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{V}}_k} = {\mathit{\boldsymbol{U}}_{k + 1}}{\mathit{\boldsymbol{B}}_k}}\\ {{\mathit{\boldsymbol{D}}^{\rm{T}}}{\mathit{\boldsymbol{U}}_{k + 1}} = {\mathit{\boldsymbol{V}}_k}\mathit{\boldsymbol{B}}_k^{\rm{T}} + {\alpha _{k + 1}}{v_{k + 1}}\mathit{\boldsymbol{e}}_{k + 1}^{\rm{T}}} \end{array}} \right. (6) 式中：ek+1T表示n阶单位矩阵的第k+1行。同理，矩阵D可以分解为上双对角矩阵。系数矩阵对角化完成后，设：
\left\{ {\begin{array}{*{20}{l}} {{r_k} = b - \mathit{\boldsymbol{D}}{x_k}}\\ {{x_k} = {\mathit{\boldsymbol{V}}_k}{y_k}}\\ {{t_{k + 1}} = {\beta _1}{e_1} - {\mathit{\boldsymbol{B}}_k}{y_k}} \end{array}} \right. (7) 联立式(6)和式(7)可得:
\begin{array}{*{20}{c}} {{r_k} = b - \mathit{\boldsymbol{D}}{x_k} = {\mathit{\boldsymbol{U}}_{k + 1}}\left( {{\beta _1}{e_1}} \right) - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{V}}_k}{y_k} = }\\ {{\mathit{\boldsymbol{U}}_{k + 1}}\left( {{\beta _1}{e_1}} \right) - {\mathit{\boldsymbol{U}}_{k + 1}}{\mathit{\boldsymbol{B}}_k}{y_k} = {\mathit{\boldsymbol{U}}_{k + 1}}{t_{k + 1}}} \end{array} (8) 由于Uk+1是正交矩阵，所以通过QR方法解yk使‖tk+1‖取极小值，求解最小二乘问题min‖β1e1- Bkyk‖，通过反复迭代，最终会收敛于极小值，从而得到式(4)的最终解。
2. 基于误差校正机制的断层成像算法
2.1 误差校正机制
经过速度反演得到每个网格单元的速度分布后，树木内部的断层成像就可以初步生成了，但是Vtest和Vref的比值的阈值0.9仅是一个经验值，可能会造成实际情况和断层图像的偏差。提出一种速度误差校正机制，用来提高成像算法的精度。
以往的断层成像算法中[13]，计算每个网格单元的传播速度Vtest和参考速度Vref后，通过计算两者的比值Vtest/Vref直接决定该网格单元是否存在缺陷区域：若比值小于0.9，则该网格单元标记为缺陷区域，否则标记为健康区域。这种评判标准有时会因为速度反演的不准确导致成像误差。事实上，树木内部的缺陷通常会从一个区域向其邻近的区域扩散，所以邻居网格单元之间的健康状况是互相关联的。提出了一种基于邻居网格单元的误差校正机制ECM，对反演得到的δ分布进行校正，用来减少上述的潜在误差。ECM分析由LSQR方法求出的每个网格单元的速度分布及其相应的Vtest和Vref比值δ分布(式9)，对当前网格单元的δ值进行调整，从而减少误差。
{\delta _i} = {V_{{\rm{ test }}}}/{V_{{\rm{ ref }}}} (9) 式中：δi代表第i个网格单元的Vtest和Vref比值，全部网格单元的δ值构成δ分布(用向量R表示)。
由于ECM是基于邻居网格单元的δ值进行校正的，在横截面边缘的网格单元会产生位于横截面外部的邻居网格单元(如图 3中编号为1、2、M-1、和M的网格单元)，称这种网格单元为虚拟单元。因此，不同于速度反演过程，在误差校正过程中，需要将包括虚拟单元在内的所有网格单元重新编号。由于虚拟单元内没有路径穿过，故将虚拟单元内的传播速度及参考速度均设置为0，其δ值为无穷大，但为了防止计算时产生数据异常，算法将其δ值校正为0。
![]()
图 3 树干横截面内网格单元的分类网格单元内的变量或数字代表该单元的编号。Figure 3. Grid cell classification of a cross-section of a tree trunkThe variable or number represent the No. of grid cells.完全处于横截面内部的网格单元以及与横截面边缘相交的网格单元统称为横截面单元，ECM对两种横截面单元应用了不同的校正策略。如图 3所示，对于完全处于横截面内部的第i个网格单元，只需要考虑其4个紧邻的网格单元(m，n，i1，i2)的δ值大小。对于第j个网格单元，其4个紧邻的网格单元中有虚拟单元l，此时需要考察包括其对角邻居单元在内的所有8个邻居单元(k1，k，k2，j1，j2，l1，l，l2)。这8个邻居单元中，第l1、l和l2个网格单元为虚拟单元，第k1、k、k2、j1和j2个网格单元为横截面单元，故将上述5个网格单元作为第j个网格单元的新邻居网格单元。
ECM是根据当前网格单元及其邻居网格单元的δ值，对其δ值进行校正，校正方案如式(10)：
{\delta _i} = \left\{ \begin{array}{l} 1, \;\;\;\;s > l/2\\ 1, \;\;\;\;s = l/2且{\delta _i} > 0.9\\ 0, \;\;\;s = l/2且{\delta _i} \le 0.9\\ 0, \;\;\;\;s < l/2 \end{array} \right. (10) 式中：δi代表第i个网格单元的δ值，s代表其邻居网格单元中δ>0.9的个数，l代表邻居网格单元的数量。该校正机制从第一个网格单元开始，逐个扫描所有网络单元，并根据公式(10)进行校正。这个过程反复进行，直到所有网格单元的δ值保持稳定为止。ECM算法的误差校正过程如下所示:
算法名：ECM误差校正机制
输入：反演生成的δ分布(用向量R表示);
输出：误差校正后的δ分布(用向量R′表示);
开始:
1：do
2：for i=1，2，…，M do
3： if第i个网格单元的四邻域均为横截面单元
4： then l=4;
5： else
6： l=第i个网格单元的八邻域中横截面单元的个数;
7： end if
8： 根据式(10)对R进行误差校正，得到R′;
9：end for
10：if R′≠R
11： then flag=1;
12：else
13： flag=0;
14： end if
15： R=R′;
16：while(flag==1)
结束
2.2 断层成像算法
基于速度反演和误差校正机制设计了断层成像算法(ECIA)。算法将树干横截面模拟成一个正圆，通过检测仪器测量应力波在树木内部的传播速度，利用LSQR方法反演计算应力波在树干横截面各个网格单元内的速度分布，并使用ECM校正速度误差，最终生成树木横截面的断层图像。
ECIA算法步骤如下：
输入.传感器个数、树木横截面周长、应力波传播时间、传感器两两间距离。
输出.树木横截面断层图像
Step1.根据输入数据建立坐标系，模拟树木横截面，划分网格单元，计算每个网格单元内的参考速度。
Step2.计算应力波传播路径被网格单元分割的距离，得到距离矩阵D。
Step3.利用LSQR反演求解速度分布，并求相应的传播速度与参考速度的比值δ分布向量R。
Step4.调用ECM算法对δ分布向量R进行误差校正。
Step5.根据各网格单元的δ值进行颜色填充。
Step6.利用区域生长算法分割缺陷区域与健康区域，对缺陷区域边缘进行曲线拟合以优化图像。
3. 结果与分析
3.1 实验材料
分别采用实验室原木和活树完成了两组无损检测成像实验。原木试样包括含空洞的柳杉(Cryptomeria fortunei，含水率11.5%)、含空洞的乌桕(Sapium sebiferum，含水率12.0%)和健康的柳杉(含水率12.5%)作为原木组实验对象，并分别记为试样1、2、3。此外，选取了江苏扬州市区6棵有内部缺陷的古树(香樟(Cinnamomum camphora)、枫香(Liquidambar formosana)、桑树(Morus alba)、圆柏(Sabina chinensis)、垂柳(Salix babylonica)、圆柏)作为实验对象，分别编号为试样4, 5, …, 9。实地检测的场景如图 4所示。实验过程中使用德国生产的PICUS应力波断层成像仪采集应力波传播数据，使用国产泰克曼TM410木材水分仪测量原木试样的含水率。除此之外，实验中还使用了Resistograph微钻阻力测试仪测量树木内部的阻力值，作为评估活树内部缺陷大小及位置的参考标准，从而验证成像算法的准确性。
3.2 原木的成像结果
利用PICUS对3个原木试样进行无损检测，采集应力波在试样内部的传播速度，并生成断层图像，3个原木试样所测截面周长分别为101、99和100 cm。基于相同的应力波传播速度数据，利用前文设计的ECIA算法产生断层图像，并与原木试样实际缺陷以及PICUS的断层成像进行对比。本组实验中，对于两种断层图像和实际的原木横截面照片，利用AutoCAD中不规则图形面积的计算方法得到原木横截面面积及缺陷面积[17]，从而计算缺陷区域的面积与整个横截面面积之比(RS)，比较断层图像与试样实际情况的RS，从而定量评估算法的准确性。
图 5~7分别展示了试样1、2、3的实际情况及两种断层成像算法所生成的图像。其中，ECIA生成的断层图像中，仅用棕色区域表示缺陷区域，PICUS生成的断层图像中，绿色、紫色和蓝色都代表不同程度的缺陷区域。因此，计算PICUS的缺陷面积时，需要计算这几部分的面积之和。
定量比较结果如表 1所示，其中，ΔRS表示断层图像的RS与真实RS的相对误差，计算方法为用它俩的差值除以真实的RS。
表 1 原木断层成像结果的定量评价Table 1. Quantitative evaluation to tomography images of logs试样编号
No. of
sample试样横截面积
Cross-sectional area
of the sample/cm2真实S
Actual S/cm2真实RS
Actual RS/%ECIA PICUS S/cm2 RS/% ΔRS/% S/cm2 RS/% ΔRS/% 1 812.593 5 297.284 1 36.584 6 324.558 8 39.941 1 -9.174 6 232.514 7 28.613 9 21.787 6 2 794.226 0 78.019 2 9.823 3 102.108 9 12.856 4 -30.876 6 285.921 4 36.000 0 -266.475 6 3 804.247 7 0 0 0 0 0 88.467 2 11.000 0 >0 注：S代表试样或断层图像中的缺陷面积；RS代表试样或断层图像中缺陷面积与横截面积的比值；ΔRS代表断层图像的RS与真实RS的相对误差。Notes：S is the defective area, RS is the ratio of the defective area to the cross-section area, and ΔRS is the relative error between the RS of ECIA or PICUS and the actual RS. 表 1中，S代表试样或断层图像中的缺陷面积，ΔRS≠0代表断层图像与试样的真实情况存在误差。从表 1可知，ECIA的断层图像比PICUS更接近于真实情况，但是其断层图像仍存在误检区域(如图 6)。较大的误检区域是由原木的低含水量造成的。国内外大量的研究与实验表明：应力波技术更适用于活树的健康检测[6]。
3.3 活树的成像结果
与原木实验类似，使用PICUS对6棵古树进行无损检测，采集应力波传播速度数据。由于无法砍伐树木获得真实的横截面图像，我们在同一高度的横截面上使用Resistograph进行微钻阻力测试。对于选定的应力波传播路径，Resistograph检测仪能够准确地测量该路径上缺陷区域的长度，经过多路径测量，可以确定缺陷区域的大体位置。图 8~19分别给出了6棵古树的实物图、ECIA和PICUS的断层图像，以及在同一横截面中若干路径的阻力曲线。断层图像和阻力曲线中Sij到Eij的部分代表试样在第i条阻力曲线中第j段缺陷区域，S和E分别代表该段缺陷区域的起点和终点。阻力曲线图中横坐标代表路径长度，蓝色部分为去掉树皮部分的真实长度，纵坐标代表相对阻力。在每条阻力曲线中，路径m-n表示与断层图像中第m个传感器到第n个传感器的路径相对应。
从试样4的外观可以观察出在1号和2号传感器之间有明显的缺陷区域，PICUS的检测结果也是如此。但通过Resistograph的阻力检测，发现1号和2号传感器之间仅有树皮缺陷，缺陷区域向树木内部扩散不深，反而在树木内部其他位置(如图 8中S12-E12、S22-E22、S31-E31)，有若干较小的缺陷区域，ECIA通过速度反演算法检测出了这些缺陷区域，并通过ECM保证了对轻微缺陷区域的检测精度。
如图 10a所示，试样5有明显的空洞，由于这个空洞过大，使得应力波只能沿接近树皮的边缘部分传播，造成了应力波数据异常，导致PICUS生成的断层图像中缺陷区域比实际的空洞小，如路径3-8中S11-E11部分在阻力曲线图中为缺陷区域，而在图 10c中不全为缺陷。ECIA使用ECM校正了由于数据异常造成的误差，使得断层图像更精确。
由于试样6的树干胸径较大，Resistograph无法穿透树木，所以只能取得路径上的部分阻力值数据，图 12b和12c中从路径起点到L点为取得阻力值数据的部分。如图 12a，试样6外部有一处明显的缺陷区域，并且被养护人员填充了特殊的材料，这些填料造成了Resistograph测量的阻力值较高。另外，如图 13a，S11-E11部分阻力值相对较低，可视为缺陷区域。ECIA生成的断层图像准确地检测出了S11-E11段的缺陷区域，但PICUS生成的断层图像中S11-E11部分仅存在长度较小的轻微腐朽，与阻力曲线图相差较大。
图 14和图 15分别给出了试样7的断层图像比较及3条应力波路径的阻力曲线。从图 14可以看出，ECIA和PICUS产生的断层图像对缺陷区域的位置及形状判断大致相同，但ECIA生成的断层图像中缺陷区域较大。路径1-7中两种算法都有不同程度增大了缺陷区域的长度，但ECIA生成的断层图像中路径1-7中的误检长度相对较小。路径2-9中，ECIA成功的检测到了S23-E23部分的缺陷而PICUS在该部分只检测到了较小长度的缺陷。同样在路径4-10中，PICUS生成的断层图像中S31-E31部分的缺陷长度比阻力曲线图 15中给出的该部分缺陷长度小，而ECIA生成的断层图像在该部分的缺陷长度与阻力曲线图中的缺陷长度较为符合。
对于试样8，两种算法生成的断层图像(如图 16b、c所示)十分相近，相比阻力曲线反映的实际情况，断层图像的路径3-9中，ECIA和PICUS对缺陷区域的位置长度判断都较为准确，在路径5-9中由于树木节子的存在(如图 17a，60~70 mm处阻力值突然增大的部分)，导致应力波传播数据异常，使得两种算法生的断层图像在该条路径上的缺陷区域的长度都增大了。
如图 18a，试样9的外部存在一个空洞，虽然空洞较小但一直扩散至树干内部，两种成像算法都成功地检测出了空洞。但在路径10-6中，S12-E12部分存在较小的缺陷区域，两种算法生成的断层图像中该部分都标记为健康区域。ECIA生成的断层图像在路径10-5中的缺陷区域与阻力曲线图 19基本一致，达到了较高的精确度，但PICUS无论在路径10-5中还是路径7-12中都存在一定程度的漏检。
活树无法像原木一样测量缺陷区域的面积，为了评估成像算法的精确度，计算了部分路径上缺陷区域的长度与路径长度的比(RL)，利用Resistograph测量真实的缺陷长度，断层图像中对应路径上的缺陷长度可以通过几何计算得到。如同原木实验中的ΔRS，定义了断层图像的RL与真实RL的相对误差ΔRL。对于6棵活树试样，ECIA和PICUS生成的断层图像的实验结果如表 2所示。实验结果证明ECIA能够更准确地检测出轻微腐朽的区域。对于试样4、5、7、9，ECIA成像结果的相对误差绝对值小于8%。由于试样6的树干存在较大空洞，试样8树干中存在节子，这些特殊情况造成检测到的应力波数据异常，使得两种算法的成像结果与实际情况相比都存在较大的误差，但ECIA的相对误差值仍小于PICUS。总体来说，ECIA的精确度比PICUS高，能够较为准确地反映实际情况。但是ECIA产生的断层图像中缺陷区域的边缘真实感略显不足，需要更完善的校正机制和算法来改进断层图像的真实感。
表 2 活树断层成像结果的定量评价Table 2. Quantitative evaluation to tomography images of standing trees试样编号
Sample No.路径
Path真实RL
Actual RLECIA PICUS RL ΔRL RL ΔRL 4 2-10 17.910 4 17.254 1 3.664 4 12.537 3 29.999 9 2-9 13.432 8 12.455 2 7.277 7 0 100.000 0 5-1 5.970 1 6.275 8 -5.120 5 18.209 0 -205.003 2 5 3-8 88.549 6 87.042 3 1.702 2 67.605 6 23.652 3 4-7 83.333 3 77.500 0 7.000 0 53.928 6 35.285 7 5-10 74.626 8 76.417 9 -2.400 1 68.656 7 7.999 9 6 9-1 75.555 5 93.043 5 -23.145 9 43.478 3 42.455 1 2-11 84.444 4 88.888 9 -5.263 2 88.888 9 -5.263 2 7 1-7 89.937 1 93.548 4 -4.015 4 95.913 9 -6.645 5 2-9 75.496 7 81.073 5 -7.386 8 56.989 2 24.514 3 4-10 91.463 4 93.566 7 -2.299 6 80.444 4 12.047 4 8 5-9 3.125 0 13.437 5 -330.000 0 14.375 0 -360.000 0 3-9 3.703 7 4.521 1 -22.069 8 12.592 6 -240.000 0 9 10-6 32.558 1 31.395 3 3.571 5 27.665 6 15.027 0 10-5 27.907 0 27.702 3 0.733 5 36.511 6 -30.833 1 7-12 46.774 2 48.387 1 -3.448 3 43.037 7 7.988 4 注：RL代表试样或断层图像中路径上缺陷区域的长度与整个路径长度的比值；ΔRL代表断层图像的RL与真实RL的相对误差。Notes: RL is the ratio of the defective length to the whole length of the path under test, and ΔRL is the relative error between the RL of ECIA or PICUS and the actual RL. 4. 结论与讨论
提出了一种适用于林木应力波无损检测的断层成像算法ECIA，算法运用LSQR方法求解速度反演方程，反演计算出横截面内每个网格单元的应力波传播速度。为了减少速度反演产生的误差，提出了一种基于邻居网格单元的误差校正机制ECM，使用ECM减少了速度反演产生的误差，从而提高了断层图像的准确度。通过对9个样本的两组实验，证实了本算法的准确度高于PICUS的准确度。本算法已在扬州市区古树名木健康监测项目中得到了实际应用。今后，我们将会完善ECM及成像算法，进一步提高断层图像的真实感。
-
![]()
图 1 内蒙古乌拉特蒙古野驴−梭梭林国家级自然保护区地理位置、功能区划及样线分布图
Figure 1. Location, functional zoning and transects of Wulate Haloxylon-Ammodendron Forest National Nature Reserve
表 1 乌拉特自然保护区鹅喉羚夏季生境选择特征值
Table 1 Eigenvalue of Gazella subgutturosa hillieriana summer habitat selection in Wulate Nature Reserve
主成分序号
Mean
component No.特征值
Eigenvalue贡献率
Ratio of contribution/%累计贡献率
Accumulative ratio of contribution/%1 2.647 22.058 22.058 2 2.176 18.132 40.190 3 1.850 15.416 55.606 4 1.384 11.529 67.135 5 1.124 9.369 76.504 6 0.723 6.026 82.530 7 0.538 4.483 87.013 8 0.504 4.201 91.214 9 0.362 3.016 94.229 10 0.294 2.450 96.679 11 0.231 1.922 98.601 12 0.168 1.399 100.000 
下载: 导出CSV
表 2 喉羚生境选择转置矩阵
Table 2 Transpose matrix of Gazella subgutturosa hillieriana habitat selection
变量 Variable 第1主成分 First principal component 第2主成分 Second principal component 第3主成分 Third principal component 第4主成分 Fourth principal component 第5主成分 Fifth principal component 坡向 Slope orientation 0.901 − 0.002 − 0.027 − 0.036 − 0.014 坡位 Slope position 0.835 0.004 0.064 − 0.055 − 0.111 坡度 Slope degree 0.811 − 0.078 − 0.014 0.012 0.103 距居民点距离
Distance from residents− 0.113 0.846 0.186 − 0.133 0.050 植被总盖度
Total vegetation coverage0.063 0.838 0.074 0.070 − 0.137 海拔 Altitude − 0.115 0.680 − 0.594 0.108 0.170 距道路距离
Distance from road− 0.070 − 0.022 0.833 0.126 0.056 距水源距离
Distance from water0.068 0.367 0.809 − 0.006 0.145 灌木高度 Shrub height 0.018 − 0.202 − 0.010 0.842 − 0.098 灌木类型 Shrub type 0.084 − 0.197 − 0.109 − 0.812 − 0.099 灌木距离 Shrub distance 0.162 − 0.071 0.350 0.018 0.820 隐蔽级 Hiding cover degree 0.429 − 0.061 0.395 0.037 − 0.631
下载: 导出CSV
表 3 乌拉特自然保护区鹅喉羚夏季对15种生态因子的生境选择
Table 3 Summer habitat selection of Gazella subgutturosa hillieriana on 15 ecological factors in Wulate Nature Reserve
生态因子
Ecological factor生态因子项目划分
Item classification of
ecological factor (I)选择/总样方数量
Selected/total quadrat
number (ri/Pi)选择系数
Selection
coefficient (Wi)选择指数
Selection
index (Ei)生境选择
Habitat selection地貌特征
Geomorphic feature高平原 High plain 8/9 0.255 0.120 P 沙地−高平原 Desert-high plain 3/5 0.172 − 0.076 AR 宽浅凹地 Wide shallow concave ground 2/4 0.143 − 0.165 NP 残丘 Torso mountain 8/9 0.255 0.120 P 沙地 Desert 27/44 0.176 − 0.065 AR 海拔
Altitude800 ~ 900 m 25/35 0.349 0.023 AR 900 ~ 1 000 m 14/24 0.285 − 0.078 AR 1 000 ~ 1 100 m 9/12 0.366 0.047 AR 坡度
Slope degree0° 31/48 0.196 − 0.121 NP 1° ~ 5° 11/17 0.197 − 0.120 NP 5° ~ 10° 5/5 0.304 0.097 AR > 10° 1/1 0.304 0.097 AR 坡位
Slope position无坡位 No slope 31/48 0.230 − 0.042 AR 中坡位 Middle slope 12/16 0.267 0.032 AR 下坡位 Down slope 3/4 0.267 0.032 AR 中下坡位 Middle-down slope 2/3 0.237 − 0.027 AR 坡向
Slope aspect无坡向 No slope 31/48 0.209 0.022 AR 阳坡 Sunny slope 3/4 0.243 0.097 AR 半阳坡 Semi-sunny slope 9/13 0.224 0.057 AR 阴坡 Shady slope 0/1 0.000 − 1.000 N 半阴坡 Semi-shady slope 5/5 0.324 0.236 P 灌木类型
Shrub type梭梭 Haloxylon ammodendron 33/50 0.103 0.016 AR 霸王 Sarcozygium xanthoxylon 14/24 0.091 − 0.046 AR 白刺 Nitraria tangutorum 12/15 0.125 0.112 P 泡泡刺 N. sphaerocarpa 27/36 0.117 0.080 AR 驼绒藜 Krascheninnikovia ceratiodes 1/3 0.052 − 0.315 NP 沙拐枣 Calligonum mongolicum 1/1 0.156 0.220 P 鹰爪柴 Convolvulus gortschakovii 3/5 0.094 − 0.032 AR 中间锦鸡儿 Caragana intermedia 2/3 0.104 0.021 AR 裸果木 Gymnocarpos przewalskii 0/1 0.000 − 1.000 N 无灌木 No shrub 4/4 0.156 0.220 P 灌木数量
Shrub number0 ~ 10 26/37 0.328 0.134 P 10 ~ 20 18/28 0.300 0.090 AR 20 ~ 30 4//4 0.373 0.197 P > 30 0/1 0.000 − 1.000 N 灌木高度
Shrub height0 ~ 0.5 m 50/79 0.172 − 0.075 AR 0.5 ~ 1.0 m 28/36 0.211 0.028 AR 1.0 ~ 1.5 m 10/14 0.194 − 0.015 AR 1.5 ~ 2.0 m 5/9 0.151 − 0.140 NP 灌木距离
Shrub distance0 ~ 2 m 20/26 0.260 0.019 AR 2 ~ 4 m 10/22 0.154 − 0.239 NP > 4 m 14/19 0.249 − 0.002 AR 隐蔽级
Hiding cover degree0 ~ 20% 6/8 0.212 0.029 AR 20% ~ 40% 15/26 0.163 − 0.102 NP 40% ~ 60% 12/18 0.188 − 0.030 AR 60% ~ 80% 10/12 0.235 0.081 AR > 80% 5/7 0.202 0.004 AR 植被盖度
Vegetation coverage0 ~ 20% 7/8 0.247 0.106 P 20% ~ 40% 8/18 0.126 − 0.228 NP 40% ~ 60% 7/15 0.132 − 0.205 NP 60% ~ 80% 11/12 0.259 0.129 P > 80% 15/18 0.236 0.082 AR 草本植物类型
Herbaceous plant type藜科 Chenopodiaceae 48/62 0.099 − 0.057 AR 菊科 Asteraceae 20/21 0.122 0.047 AR 蒺藜科 Zygophyllaceae 25/26 0.123 0.051 AR 禾本科 Gramineae 16/20 0.102 − 0.040 AR 豆科 Leguminosae 8/8 0.128 0.071 AR 百合科 Liliaceae 9/11 0.105 − 0.029 AR 锁阳科 Cynomoriaceae 2/2 0.128 0.071 AR 列当科 Orobanchaceae 1/2 0.064 − 0.269 NP 无草本 No herb 1/1 0.128 0.071 AR 距水源距离
Distance from water0 ~ 10 km 9/15 0.227 − 0.048 AR 10 ~ 20 km 6/9 0.252 0.005 AR 20 ~ 30 km 15/24 0.237 − 0.028 AR > 30 km 18/24 0.284 0.064 AR 距道路距离
Distance from road0 ~ 5 km 27/43 0.289 − 0.071 AR 5 ~ 10 km 10/16 0.288 − 0.073 AR 10 ~ 15 km 11/12 0.423 0.118 P 距居民点距离
Distance from residents0 ~ 5 km 11/17 0.223 − 0.057 AR 5 ~ 10 km 20/34 0.203 − 0.104 NP 10 ~ 15 km 11/12 0.316 0.116 P 15 ~ 20 km 6/8 0.258 0.017 AR 注：P表示偏好或选择，AR表示几乎随机选择，NP表示没有偏好，N表示回避选择。Notes: P indicates preference or choice, AR means almost random selection, NP indicates no preference, and N indicates no choice at all.
下载: 导出CSV
-
[1] 李晓军, 龚大洁, 郑智, 等. 鹅喉羚[J]. 生物学通报, 2015, 50(4):13−14. Li X J, Gong D J, Zheng Z, et al. Goitred gazelle[J]. Bulletin of Biology, 2015, 50(4): 13−14.
[2] 汪松. 中国物种红色名录(脊椎动物)[M]. 北京: 高等教育出版社, 2009: 549−550. Wang S. China species red list (vertebrata)[M]. Beijing: Higher Education Press, 2009: 549−550.
[3] Ostrowski S, Williams J B. Heterothermy of free-living Arabian sand gazelles (Gazella subgutturosa marica) in a desert environment[J]. The Journal of Experimental Biology, 2006, 209: 1421−1428. doi: 10.1242/jeb.02151
[4] Cunningham P L, Wronski T. Seasonal changes in group size and composition of Arabian sand gazelle Gazella subgutturosa marica Thomas, 1897 during a period of drought in central western Saudi Arabia[J]. Current Zoology, 2011, 57(1): 36−41. doi: 10.1093/czoolo/57.1.36
[5] Dittmann M T, Hebel C, Arif A, et al. Metabolic rates of three gazelle species (Nanger soemmerringii, Gazella gazella, Gazella spekei) adapted to arid habitats[J]. Mammalian Biology, 2015, 80: 390−394. doi: 10.1016/j.mambio.2015.05.008
[6] Cunningham P L. Observations of the seasonal dietary preference of male Gazella subgutturosa marica Thomas, 1897 (Cetartiodactyla: Bovidae) along foraging trails of central Saudi Arabia[J]. Journal of Threatened Taxa, 2009, 1(9): 445−448. doi: 10.11609/JoTT
[7] Lerp H, Klaus S, Allgöwer S. Data on phylogenetic analyses of gazelles (genus Gazella) based on mitochondrial and nuclear intron markers[J]. Data in Brief, 2016, 7: 551−557. doi: 10.1016/j.dib.2016.02.062
[8] BlankD A, Ruckstuhl K, Yang W. Antipredator strategy of female goitered gazelles (Gazella subgutturosa Guld. 1780) with hiding fawn[J]. Behavioural Processes, 2015, 119: 44−49. doi: 10.1016/j.beproc.2015.07.013
[9] Fadakar D, Rezaei H R, Naseri M, et al. Phylogenetic analysis of Persian Gazella, Gazella subgutturosa (Artiodactyla: Bovidae) based on cytochrome b in central Iran[J]. Molecular Biology Research Communications, 2013, 2(4): 151−159.
[10] 高行宜, 姚军. 新疆哈密盆地初冬鹅喉羚的地理分布与种群数量[J]. 干旱区地理, 2006, 29(2):213−218. doi: 10.3321/j.issn:1000-6060.2006.02.008 Gao X Y, Yao J. Study on the geography distribution and population of Gazella subgutturosa in the Hami Basin, Xinjiang in early winter[J]. Arid Land Geography, 2006, 29(2): 213−218. doi: 10.3321/j.issn:1000-6060.2006.02.008
[11] 初红军, 蒋志刚, 蒋峰, 等. 鹅喉羚夏季和冬季卧息地选择[J]. 动物学研究, 2009, 30(3):312−317. Chu H J, Jiang Z G, Jiang F, et al. Summer and winter bed-site selection by goitred gazelle (Gazella subgutturosa sairensis)[J]. Zoological Research, 2009, 30(3): 312−317.
[12] 杨维康, 徐文轩, 刘伟, 等. 卡拉麦里山有蹄类保护区鹅喉羚的采食地选择[J]. 干旱区研究, 2010, 27(2):236−240. Yang W K, Xu W X, Liu W, et al. Feeding habitat characteristics of Gazella subgutturosa in the Karamaili Mountain Nature Reserve, Xinjiang[J]. Arid Zone Research, 2010, 27(2): 236−240.
[13] 李莹, 徐文轩, 乔建芳, 等. 卡拉麦里山有蹄类自然保护区鹅喉羚的时空分布与生境选择[J]. 干旱区地理, 2009, 32(2):261−265. Li Y, Xu W X, Qiao J F, et al. Spatio-temporal distribution and habitat selection of Gazella subgutturosa in Kalamaili Mountain Nature Reserve in four seasons[J]. Arid Land Geography, 2009, 32(2): 261−265.
[14] 艾孜孜江•乃比, 夏米西丁•阿不都热依木, 热木图拉•阿卜杜克热木, 等. 艾比湖国家级湿地自然保护区鹅喉羚种群数量与分布[J]. 四川动物, 2015, 34(3):447−452. Ezizjan N, Shamshidin A, Rahmutulla A, et al. Populatin size and distribution of Gazella subgutturosa in Ebinur Lake National Wetland Nature Reserve[J]. Sichuan Journal of Zoology, 2015, 34(3): 447−452.
[15] 塔吉古丽•吐热甫, 艾孜孜江•乃比, 宁礼群, 等. 新疆艾比湖鹅喉羚mtDNA D-Loop区序列多态性研究[J]. 中国畜牧兽医, 2016, 43(5):1316−1322. Tajigul T, Ezizjan N, Ning L Q, et al. Genetic diversity of Gazella subgutturosa based on mitochondrial DNA D-Loop sequence in Xinjiang Ebinur Lake[J]. Chinese Animal Husbandry and Veterinary Medicine, 2016, 43(5): 1316−1322.
[16] 赵宠南, 苏云, 刘振生, 等. 贺兰山牦牛冬春季的生境选择[J]. 生态学报, 2012, 32(6):1763−1764. Zhao C N, SU Y, Liu Z S, et al. Habitat selection of feral yak in winter and spring in the Helan Mountains, China[J]. Acta Ecological Sinica, 2012, 32(6): 1763−1764.
[17] 宋朝枢, 贾昆峰. 乌拉特梭梭林自然保护区科学考察集[M]. 北京: 中国林业出版社, 2000: 1−111. Song C S, Jia K F. Scientific survey of Wulate Haloxylon Ammodendron Forest Nature Reserve[M]. Beijing: China Forestry Publishing House, 2000: 1−111.
[18] 鲁庆彬, 胡锦矗. 岷山黑熊生境选择的初步分析[J]. 兽类学报, 2003, 23(2):98−103. doi: 10.3969/j.issn.1000-1050.2003.02.002 Lu Q B, Hu J C. Preliminary analysis on the habitat selection of black bears in the Minshan Mountains[J]. Acta Theriologica Sinica, 2003, 23(2): 98−103. doi: 10.3969/j.issn.1000-1050.2003.02.002
[19] 徐文轩, 乔建芳, 夏参军, 等. 卡拉麦里山保护区鹅喉羚卧息地特征的季节变化[J]. 生态学杂志, 2010, 29(4):688−691. Xu W X, Qiao J F, Xia C J, et al. Seasonal variations of bedding site characteristics of Gazella subgutturosa in Kalamaili Mountain Nature Reserve[J]. Chinese Journal of Ecology, 2010, 29(4): 688−691.
[20] 朱云云, 王孝安, 王贤, 等. 坡向因子对黄土高原草地群落功能多样性的影响[J]. 生态学报, 2016, 36(21):6823−6831. Zhu Y Y, Wang X A, Wang X, et al. Effect of slope aspect on the functional diversity of grass communities in the Loess Plateau[J]. Acta Ecological Sinica, 2016, 36(21): 6823−6831.
-
期刊类型引用(13)
1. 夏甜甜，杜秋硕，刘鹏. 基于数字赋能的泰山古柏生态修复技术应用探究. 山东林业科技. 2024(02): 77-83 . 
百度学术
2. 叶程浩，杜晓晨，鲍宇. 一种基于YOLOv5和轮廓约束的木材空洞缺陷应力波层析成像算法. 木材科学与技术. 2023(05): 61-68 . 百度学术
3. 侯艳芳，王鑫，侯婷婷，郭靖，张艳. 古建筑木结构损伤检测方法研究进展. 粘接. 2022(01): 123-126+130 . 百度学术
4. 宋辉，姬晓文，高琦，缑鹏超，王志强. 基于数据驱动的不规则用电检测误差校正系统设计. 电子设计工程. 2022(06): 109-112+117 . 百度学术
5. 杨露露，董喜斌，徐华东. 电阻法和应力波法在活立木内部腐朽缺陷检测中的对比. 森林工程. 2022(04): 82-88 . 百度学术
6. 戴禄君，吴瑞瑞. 数据驱动的非线性光学显微成像误差校正研究. 激光杂志. 2022(07): 159-163 . 百度学术
7. 路程琳，何宇航，谢雨宏，刘依萍，刘嘉仪，王正. 应力波断层扫描法在立木检测方面的研究现状. 林业机械与木工设备. 2022(10): 37-40 . 百度学术
8. 刘涛，李光辉. 基于射线分割的林木应力波断层成像算法. 林业科学. 2021(09): 181-192 . 百度学术
9. 赵松，王勐. 移动AR+VR支持下全景图像拼接均匀性校正方法. 计算机仿真. 2021(11): 189-192+398 . 百度学术
10. 何宇航，陈偲，王弘历，黄席阳，王正. 断层扫描法在木材无损检测的研究进展. 木工机床. 2020(01): 13-16 . 百度学术
11. 张慧娟. 木材应力波检测精度的影响因素研究. 中国林副特产. 2020(01): 27-30 . 百度学术
12. 魏喜雯，孙丽萍，许述正，杨扬，杜春晓. 基于应力波传播速度模型的原木缺陷定量检测. 北京林业大学学报. 2020(05): 143-154 . 本站查看
13. 刘嘉新，高景泉，李超. 应用兰德韦伯算法的木材缺陷图像重建. 东北林业大学学报. 2019(12): 125-128 . 百度学术
其他类型引用(8)
计量
- 文章访问数: 3469
- HTML全文浏览量: 940
- PDF下载量: 94
- 被引次数: 21
