本发明涉及一种基于射线分割的检测树木内部缺陷的断层成像方法，属于林木的无损检测领域技术领域。
背景技术：
应力波断层成像技术是常见的林木无损检测技术之一。利用分布在树干周围的传感器记录应力波在树木内部的传播时间，分析应力波在树木截面中的速度分布，实现树木内部缺陷的快速检测。
通常情况下，应力波在健康木质中的传播速度快于腐朽木质，当遇到空洞、裂缝等情况时，应力波会跨过该区域，导致传播速度明显降低，从而区分出木质的健康程度。
目前，国内外对于树木内部缺陷的断层成像检测已经进行了较为广泛的研究，可以对林木缺陷进行有效诊断，但是由于应力波传播射线的数量有限，导致成像区域内网格单元的参考速度值过少，甚至部分网格单元没有参考速度值，由此导致待插值网格单元的估计值与实际情况相差较大，从而导致成像精度都不高，包括缺陷面积、形状和位置等信息。
技术实现要素：
为了解决现有树木无损检测方法对树木进行检测时，对于树木内部缺陷、面积、形状和位置成像精度低的问题，本发明提供了一种基于射线分割的检测树木内部缺陷的断层成像方法，所述方法采用射线分割的方法，得到更为精细的应力波传播速度分布，提高了网格单元速度参考值的相关性和精度，从而得到树木内部更为精确的成像效果，以便检测出树木内部缺陷、面积、形状和位置的准确信息。
一种基于射线分割的检测树木内部缺陷的断层成像方法，所述方法包括：
s1获取应力波在待检测树木内部的传播速度并进行修正，对修正后的应力波速度值进行归一化和离散化处理；
s2在成像区域内划分网格单元，并根据归一化后的应力波传播速度绘制应力波传播射线图；
s3分割根据归一化后的应力波传播速度绘制应力波传播射线图中的所有射线，并根据已知传播射线速度值估算待分割射线上多条线段的速度值；将分割后的射线进行可视化，得到改进的应力波传播射线图；
s4计算成像区域内每个网格单元的速度值；对无数据的网格单元采用最近邻插值法获取对应的速度值；
s5根据成像区域内每个网格单元的速度值判定网格状态，结合图像处理方法实现树木内部缺陷的重建，得到树木内部缺陷的图像。
可选的，所述s3包括：
s31选择s2中根据归一化后的应力波传播速度绘制应力波传播射线图中的任意一条传播射线作为待分割射线，利用已知应力波传播速度的多条射线对待分割射线进行分割，根据多个交点将待分割射线分割成多条线段；
s32通过相邻两交点的速度值重新估算待分割射线上各条线段的速度值，线段估算值可表示为：
sv＝(start end)/2公式(2)
其中，start和end分别表示每条线段两端的交点的速度值，sv表示该线段的速度值；
s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图。
可选的，所述s4包括：
根据下述公式(3)判断各个网格是否在成像区域内：
其中，r为该树木样本的半径，(x，y)为网格左下角的坐标，若f(x，y)＝1，则表示该网格在成像区域内，若f(x，y)＝0，则表示该网格在成像区域外；
计算成像区域内的网格的速度估计值。
可选的，若网格单元有应力波射线穿过，该网格的速度估计值为所有穿过该网格的射线速度值的均值，根据下述公式(4)计算该网格的速度估计值。
其中，vxy为某一个成像区域内网格单元的速度估计值，vi为穿过该网格单元的第i条射线的速度值，n为穿过该网格单元的射线总数，1≤i≤n；
若网格单元没有应力波射线穿过，该网格的速度估计值根据下述公式(5)计算得到：
gxy＝w1v1 w2v2 w3v3 w4v4公式(5)
其中gxy为无应力波射线穿过的网格速度值，v1，v2，v3，v4表示距该网格单元欧氏距离最近的4个非空网格单元速度值，w1，w2，w3，w4为对应速度值的权重，空网格单元为无应力波射线穿过的网格。
可选的，所述s1中获取应力波在待检测树木内部的传播速度并进行修正，包括：
在待检测树木的横截面周围均匀地布置若干传感器，通过敲击传感器产生应力波来获取任意两传感器间的应力波传播时间，并根据传感器间的距离估算出每条应力波传播射线的速度值；
根据下述公式(1)进行速度修正：
vt≈(1-0.2β2)*vr公式(1)
其中，vt为应力波沿弦向方向的传播速度，即修正后的速度；vr为应力波沿直径方向的传播速度，β为弦线方向与径向方向间的夹角。
可选的，所述s1中对修正后的应力波速度值进行归一化和离散化处理，包括：
将修正后的应力波传播速度值归一化至区间[0，1]内；
设定阈值θ1、θ2，将归一化后应力波传播速度值进行离散化；区间[0，θ1]表示严重腐朽区域，区间[θ1，θ2]表示轻微腐朽区域，区间[θ2，1]表示健康区域。
可选的，所述s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图，包括：
判断各射线每条线段的速度值处于区间[0，θ1]、区间[θ1，θ2]还是区间[θ2，1]；
根据各射线每条线段所属区间绘制出改进后的应力波传播射线图。
可选的，所述设定阈值θ1、θ2为：
根据先验知识确定不同树木对应的阈值θ1、θ2。
可选的，所述w1，w2，w3，w4为根据对应的欧氏距离的大小确定。
可选的，所述方法选用12个传感器获取应力波在待检测树木内部的传播速度。
本发明还提供一种基于射线分割的树木内部缺陷检测方法，所述方法采用上述断层成像方法获取树木内部缺陷的图像，进而根据树木内部缺陷的图像测量得到树木内部缺陷面积、形状和位置信息。
本发明有益效果是：
本申请方法在获得根据归一化后的应力波传播速度绘制的应力波传播射线图后，对其中的每条射线进行分割，分割时根据已知传播射线速度值估算该待分割射线上多条线段的速度值，将每条射线分割为不同速度的多条线段，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图，进而判断出各个网格的网格状态，从而使得各个网格所反应的树木内部情况更贴近于真实的情况，得到树木内部更为精确的成像效果，为后续检测出树木内部缺陷面积、形状和位置的准确信息提供更精确的信息。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明一个实施例中的无损检测的实验平台的示意图。
图2是树木横截面应力波传播示意图。
图3是本发明一个实施例中的速度修正示意图。
图4是本发明一个实施例中的射线分割算法示意图。
图5是本发明一个实施例中的射线分割示意图。
图6是本发明一个实施例中的断层成像结果仿真图。
图7是本发明一个实施例中分别采用hwpm方法和本申请方法得到的健康雪松样本的成像结果。
图8是本发明一个实施例中分别采用hwpm方法和本申请方法得到的缺陷红松样本的成像结果。
图9是本发明一个实施例中分别采用hwpm方法和本申请方法得到的缺陷乌桕样本的成像结果。
图10是本发明一个实施例中分别采用hwpm方法和本申请方法得到的缺陷雪松样本的成像结果。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
实施例一：
本实施例提供一种基于射线分割的检测树木内部缺陷的断层成像方法，所述方法包括：
s1获取应力波在待检测树木内部的传播速度并进行修正，对修正后的应力波速度值进行归一化和离散化处理；
s2在成像区域内划分网格单元，并根据归一化后的应力波传播速度绘制应力波传播射线图；
s3分割根据归一化后的应力波传播速度绘制应力波传播射线图中的所有射线，并根据已知传播射线速度值估算待分割射线上多条线段的速度值；将分割后的射线进行可视化，得到改进的应力波传播射线图；
s31选择s2中根据归一化后的应力波传播速度绘制应力波传播射线图中的任意一条传播射线作为待分割射线，利用已知应力波传播速度的多条射线对待分割射线进行分割，根据多个交点将待分割射线分割成多条线段；
s32通过相邻两交点的速度值重新估算待分割射线上各条线段的速度值，线段估算值可表示为：
sv＝(start end)/2公式(2)
其中，start和end分别表示每条线段两端的交点的速度值，sv表示该线段的速度值；
s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图。
s4计算成像区域内每个网格单元的速度值；对无数据的网格单元采用最近邻插值法获取对应的速度值；
s5根据成像区域内每个网格单元的速度值判定网格状态，结合图像处理方法实现树木内部缺陷的重建，得到树木内部缺陷的图像。
实施例二：
本实施例提供一种基于射线分割的检测树木内部缺陷的断层成像方法，以对原木进行无损检测为例进行说明。
所述方法包括：
步骤一：获取应力波在待检测树木内部的传播速度并进行修正，对修正后的应力波速度值进行归一化和离散化处理；具体的，包括以下步骤：
s11在进行无损检测时，搭建图1所示的检测平台，在待检测原木的横截面周围均匀地布置一定数量的传感器，通过敲击传感器产生来获取任意两传感器间的应力波传播时间，并根据传感器间的距离估算出每条应力波传播射线的速度值，实际应用中可采用fakopp应力波检测仪来采集应力波在树木内部的传播时间。
为进行区分，可以对所使用的传感器进行编号，本实施例中使用了12个传感器，编号为1至12。
s12由于树木的各向异性体特征导致应力波在树干内的传播速度不均匀。通常情况下，应力波在树木边缘的传播速度小于沿径向传播的速度，而应力波传播速度是林木缺陷成像的主要依据，在进行缺陷可视化时，树干边缘会因为速度较慢而导致误判。所以需要通过使用公式(1)进行速度修正。
vt≈(1-0.2β2)*vr公式(1)
请参考图2，其中，vt为应力波沿弦向方向的传播速度，vr为应力波沿直径方向的传播速度，β为弦线方向与径向方向间的夹角；图2中，sr1表示弦向方向，sr2表示直径方向。
s13将修正后的应力波传播速度值归一化至区间[0，1]内。
s14设定阈值θ1、θ2，将应力波传播速度值进行离散化；区间[0，θ1]表示严重腐朽区域，区间[θ1，θ2]表示轻微腐朽区域，区间[θ2，1]表示健康区域。
步骤二：在成像区域内划分网格单元，并根据归一化后的应力波传播速度绘制应力波传播射线图；
根据归一化后的应力波传播速度绘制应力波传播射线图，如图3所示，图3(a)为待检测原木的横截面，图3(b)为根据原始速度绘制得到的应力波传播射线图，原始速度即为直接根据传感器间的距离估算出每条应力波传播射线的速度值进行归一化处理后的数值，图3(c)为根据修正后的速度进行归一化处理后的数值绘制得到的应力波传播射线图。
步骤三：分割根据修正后的速度进行归一化处理后的数值绘制得到的应力波传播射线图中的所有射线，根据已知传播射线速度值估算待分割射线上多条线段的速度值。将分割后的射线进行可视化，得到改进的应力波传播射线图。
具体的，包括以下步骤：
s31选择任意一条待分割射线，利用已知应力波传播速度的多条射线对待分割射线进行分割，根据多个交点将待分割射线分割成多条线段，如图4所示，图4中，对应7号传感器和1号传感器之间的射线被其他7条射线分割，加上端点共计有7个交点。
s32通过相邻两交点的速度值重新估算待分割射线上各条线段的速度值，线段估算值可表示为：
sv＝(start end)/2公式(2)
其中，start和end分别表示每条线段两端的交点的速度值，sv表示该线段的速度值。
s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图。
判断各射线每条线段的速度值处于区间[0，θ1]、区间[θ1，θ2]还是区间[θ2，1]，根据所属区间绘制出改进后的应力波传播射线图。
如图5所示，图5(a)与图3(a)所示的待检测原木的横截面；图5(b)为图3(c)所示的根据修正后的速度进行归一化处理后的数值绘制得到的应力波传播射线图；图5(c)为采用上述分割步骤进行射线分割后的射线图。
步骤四：计算成像区域内每个网格单元的速度值，并对无数据的网格单元采用最近邻插值，确保所有网格单元都被处理。最后根据网格单元的速度值判定网格状态，结合图像处理方法实现树木内部缺陷的重建，得到树木内部缺陷的图像。
具体的，包括：
s41计算成像区域内每个网格单元的速度值。
首先，根据下述公式(3)判断各个网格是否在成像区域内：
其中，r为该树木样本的半径，(x，y)为网格左下角的坐标，若f(x，y)＝1，则表示该网格在成像区域内，若f(x，y)＝0，则表示该网格在成像区域外。
若网格单元有应力波射线穿过时，该网格的速度估计值为所有穿过该网格的射线速度值的均值，根据下述公式(4)计算该网格的速度估计值。
其中，vxy为某一个成像区域内网格单元的速度估计值，vi为穿过该网格单元的第i条射线的速度值，n为穿过该网格单元的射线总数，1≤i≤n。
但在实际情况中，由于应力波传播射线过少，网格单元较多，会出现部分网格单元未被射线穿过，即存在部分网格单元没有参考速度值，对于没有参考速度值的网格单元，则选用最近邻插值法进行数据补全，具体的，根据下述公式(5)计算该网格单元的速度估计值：
gxy＝w1v1 w2v2 w3v3 w4v4公式(5)
其中gxy为无参考射线的网格速度值，v1，v2，v3，v4表示距该网格单元欧氏距离最近的4个非空网格单元速度值，w1，w2，w3，w4为对应速度值的权重。
s42当成像区域内所有网格单元都被处理后，根据网格单元的速度值来判定该网格单元为健康区域或腐朽区域。最后结合图像处理的方法，生成最终的树木横截面断层图像，可以直观地看到缺陷的位置、大小和形状。
为验证本申请方法能够获得较为精确的树木横截面断层图像，以便能够更精确的确定树木内部缺陷、面积、形状和位置的准确信息，特在多个原木样本上进行实验，并与hwpm方法得到的成像结果进行对比，通过统计缺陷面积来评估成像的准确性；多个原木样本分别为健康雪松样本(简称雪松1)、缺陷红松样本(简称红松)、缺陷乌桕样本(简称乌桕)和缺陷雪松样本(简称雪松2)。
hwpm方法参见“leil,guanghuil.acoustictomographybasedonhybridwavepropagationmodelfortreedecaydetection[j].computers&electronicsinagriculture,2018,151:276-285.”一文中的介绍。
实验过程中，根据先验知识得到的各个原木样本对应的阈值θ1、θ2的取值分别如下表1所示，其中雪松和红松对应的θ1、θ2的取值相同：
表1：四种原木样品分别对应的阈值θ1、θ2的取值
具体的，实验结果如图7-10和表2所示，图7-10分别对应雪松1、红松、乌桕和雪松2的原木情况和对应的成像结果，其中(a)为对应的原木样本的横截面图像，(b)为采用hwpm方法得到的成像结果，(c)为采用本申请方法得到的成像结果，所统计的各个参数缺陷面积如下表2所示：
表2：四种原木样品的实际缺陷面积以及分别采用hwpm方法和本发明方法得到的缺陷面积数据
由上表2可知，本申请方法根据所得到的成像图像测得的缺陷面积值，比根据hwpm方法的成像结果测得的缺陷面积更贴近真实值。
根据图7-10可知，本申请方法所得到的成像图像中，缺陷的位置以及形状都更贴近于真实的缺陷位置和形状，由此可知，本申请方法能够得到树木内部更为精确的成像效果，为后续检测出树木内部缺陷面积、形状和位置的准确信息提供更精确的信息。
本发明实施例中的部分步骤，可以利用软件实现，相应的软件程序可以存储在可读取的存储介质中，如光盘或硬盘等。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
1.一种基于射线分割的检测树木内部缺陷的断层成像方法，其特征在于，所述方法包括：
s1获取应力波在待检测树木内部的传播速度并进行修正，对修正后的应力波速度值进行归一化和离散化处理；
s2在成像区域内划分网格单元，并根据归一化后的应力波传播速度绘制应力波传播射线图；
s3分割根据归一化后的应力波传播速度绘制应力波传播射线图中的所有射线，并根据已知传播射线速度值估算待分割射线上多条线段的速度值；将分割后的射线进行可视化，得到改进的应力波传播射线图；
s4计算成像区域内每个网格单元的速度值；对无数据的网格单元采用最近邻插值法获取对应的速度值；
s5根据成像区域内每个网格单元的速度值判定网格状态，结合图像处理方法实现树木内部缺陷的重建，得到树木内部缺陷的图像。
2.根据权利要求1所述的方法，其特征在于，所述s3包括：
s31选择s2中根据归一化后的应力波传播速度绘制应力波传播射线图中的任意一条传播射线作为待分割射线，利用已知应力波传播速度的多条射线对待分割射线进行分割，根据多个交点将待分割射线分割成多条线段；
s32通过相邻两交点的速度值重新估算待分割射线上各条线段的速度值，线段估算值可表示为：
sv＝(start end)/2公式(2)
其中，start和end分别表示每条线段两端的交点的速度值，sv表示该线段的速度值；
s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图。
3.根据权利要求2所述的方法，其特征在于，所述s4包括：
根据下述公式(3)判断各个网格是否在成像区域内：
其中，r为该树木样本的半径，(x，y)为网格左下角的坐标，若f(x,y)＝1，则表示该网格在成像区域内，若f(x,y)＝0，则表示该网格在成像区域外；
计算成像区域内的网格的速度估计值。
4.根据权利要求3所述的方法，其特征在于，
若网格单元有应力波射线穿过，该网格的速度估计值为所有穿过该网格的射线速度值的均值，根据下述公式(4)计算该网格的速度估计值。
其中，vxy为某一个成像区域内网格单元的速度估计值，vi为穿过该网格单元的第i条射线的速度值，n为穿过该网格单元的射线总数，1≤i≤n；
若网格单元没有应力波射线穿过，该网格的速度估计值根据下述公式(5)计算得到：
gxy＝w1v1 w2v2 w3v3 w4v4公式(5)
其中gxy为无应力波射线穿过的网格速度值，v1，v2，v3，v4表示距该网格单元欧氏距离最近的4个非空网格单元速度值，w1，w2，w3，w4为对应速度值的权重，空网格单元为无应力波射线穿过的网格。
5.根据权利要求4所述的方法，其特征在于，所述s1中获取应力波在待检测树木内部的传播速度并进行修正，包括：
在待检测树木的横截面周围均匀地布置若干传感器，通过敲击传感器产生应力波来获取任意两传感器间的应力波传播时间，并根据传感器间的距离估算出每条应力波传播射线的速度值；
根据下述公式(1)进行速度修正：
vt≈(1-0.2β2)*vr公式(1)
其中，vt为应力波沿弦向方向的传播速度，即修正后的速度；vr为应力波沿直径方向的传播速度，β为弦线方向与径向方向间的夹角。
6.根据权利要求5所述的方法，其特征在于，所述s1中对修正后的应力波速度值进行归一化和离散化处理，包括：
将修正后的应力波传播速度值归一化至区间[0,1]内；
设定阈值θ1、θ2，将归一化后应力波传播速度值进行离散化；区间[0,θ1]表示严重腐朽区域，区间[θ1,θ2]表示轻微腐朽区域，区间[θ2,1]表示健康区域。
7.根据权利要求6所述的方法，其特征在于，所述s33重复上述步骤，至所有射线都经过分割处理，并根据各射线每条线段的速度值绘制出改进后的应力波传播射线图，包括：
判断各射线每条线段的速度值处于区间[0,θ1]、区间[θ1,θ2]还是区间[θ2,1]；
根据各射线每条线段所属区间绘制出改进后的应力波传播射线图。
8.根据权利要求7所述的方法，其特征在于，所述w1，w2，w3，w4为根据对应的欧氏距离的大小确定。
9.根据权利要求8所述的方法，其特征在于，所述方法选用12个传感器获取应力波在待检测树木内部的传播速度。
10.一种基于射线分割的树木内部缺陷检测方法，其特征在于，所述方法采用权利要求1-9任一所述的断层成像方法获取树木内部缺陷的图像，进而根据树木内部缺陷的图像测量得到树木内部缺陷面积、形状和位置信息。
技术总结