技术领域
[0001] 本
发明涉及
电池测试领域，特别是一种基于两期无人机影像的针叶林采伐蓄积量估算方法。
背景技术
[0002] 森林蓄积量是指一定森林面积上存在着的林木树干部分的总材积(孟宪宇，2006)。它是反映森林资源的丰富程度、衡量森林生态环境优劣的重要依据，也是森林固
碳能
力的重要标志和森林地上
生物量部分的重要组成。因此，蓄积量是重要的林分调查因子，它为森林的采伐利用和经营提供重要的依据(谢士琴等，2017)。由于劳动强度大，调查周期长，费时费力等缺点，传统的一、二类森林资源调查方法获取的森林蓄积量数据不能及时准确地反映大面积森林蓄积量的分布状况及动态变化(高广磊等，2013；刘清旺等，2017；史洁清等，2017)。
[0003] 随着遥感技术在林业调查中的应用，国内外学者利用RS、GIS、 GPS技术进行了森林蓄积量估测的相关研究。按照数据源与研究方法的不同，这些研究大致可以分为两类：数据源上，主要有遥感影像数据，LiDAR点
云数据和多源数据的相结合。然而机载LiDAR点云获取的成本较高，无人机遥感具有操作简单、高效快捷、作业成本低等优势，无人机遥感所获取的影像
分辨率可达厘米级，是传统卫星遥感无法比拟的，有利于林分参数的目视判读，可以大幅度替代繁重的人工野外调查工作，提高调查效率。
[0004] 综合国内外学者研究内容发现，蓄积量估算的研究区多数为地势平坦的林业样地，而对于如福建省内山地海拔高差大，地势陡峭的林业小班区域关注较少，特别是对于森林采伐蓄积量无人机遥感估算目前鲜有研究报道。如何消除地形对树高提取的影响从而提高蓄积量估算的
精度仍是研究重点。
发明内容
[0005] 有鉴于此，本发明的目的是提出一种基于两期无人机影像的针叶林采伐蓄积量估算方法，能够有效地消除地形对树高提取的影响，从而提高蓄积量估算的精度。
[0006] 本发明采用以下方案实现：一种基于两期无人机(Unmaned Aerial Vehicle：UAV)影像的针叶林采伐蓄积量估算方法，具体包括以下步骤：
[0007] 步骤S1：利用无人机获取伐区采伐前后两期的航拍影像数据，经处理得到两期镶嵌影像和影像点云数据；
[0008] 步骤S2：对两期影像点云数据进行匹配；
[0009] 步骤S3：对步骤S2匹配后的数据进行布料模拟
算法滤波处理，分离出林冠点云和地面点云；
[0010] 步骤S4：基于步骤S3获取的点云数据，利用自然邻域插值法生成数字表面模型和数字高程模型，进而获取冠层高度模型；
[0011] 步骤S5：基于步骤S4中构建的冠层高度模型，采用动态窗口局域最大值法进行林冠
顶点探测，从而得到株数信息与对应树高信息；
[0012] 步骤S6：基于实测数据研建胸径-树高回归模型，利用步骤S5 的结果获取对应的单木胸径信息；并利用影像抽样或监督分类获得针叶林树种组成比例信息；
[0013] 步骤S7：利用二元立木材积公式计算相应树种的单木材积，结合株数信息得到采伐蓄积量。
[0014] 进一步地，步骤S1中，所述无人机飞行设置的航向重叠度优于 80％，旁向重叠度优于60％；所获的无人机影像分辨率优于10cm；点云
密度高于10/m3。
[0015] 进一步地，步骤S2具体为：对两期影像点云数据选取多个同位点进行配准，数据的匹配误差要求在两个像元之内，保证研究区域中每棵树木点云的
位置与相应地面点云位置完好对应。
[0016] 进一步地，步骤S2进一步具体为：将获得的两期影像点云数据看做从P,Q两个不同视
角下获得的具有重叠区域的数据点集，若在数据点集P中提取一个子集{mi|mi∈P,i＝1,2,...,N}，在数据点集Q中有一子集{mi′|mi′∈Q,i＝1,2,...,N}与之各点对应，配准就是求解在不同视角 (
坐标系)下，三维数据点之间的坐标转换关系R和T，其中R为旋转矩阵，T为平移矩阵，使得目标函数达到最小化；即多视角对齐的实质是计算满足如下目标函数的旋转和平移变换矩阵：
[0017] f(R,T)＝min∑[R·pi+T-qi]2；
[0018] 式中，pi与qi为需要对齐的点。
[0019] 进一步地，步骤S3具体包括以下步骤：
[0020] 步骤S31：转换点云几何坐标，使配准之后的点云数据翻转；
[0021] 步骤S32：初始化布料格网，通过格网分辨率确定格网
节点数量；
[0022] 步骤S33：将点云点和格网点投影到同一
水平面，寻找每个格网点对应的点云点，并记录对应点云点的高程值；
[0023] 步骤S34：计算格网点受重力作用位移的位置，并比较该节点位置高程与其对应点云点高程，如果该节点位置高程小于或者等于点云点高程，则将该节点位置替换到对应点云点位置，并标记为不可动点；其中格网点受重力作用位移的位置计算公式为：
[0024]
[0025] 式中，m是布料的
质量，被设为常数1；X是某个时刻的节点位置；Δt是时间步长，G是重力常数；给出时间步长和节点初始位置，就能够得到当前点的位置；
[0026] 步骤S35：计算每个格网点受邻近节点影响而移动的位置，每个布料点的修正位移计算公式为：
[0027]
[0028] 式中，为节点的移动向量； 为待移动节点的当前位置； 为 的邻近节点位置；为垂直向上的标准向量， b是判断节点是否移动的参量，当节点可移动时，设b＝1，否则b＝0；
[0029] 步骤S36：重复步骤S34和步骤S35，当所有节点的最大高程变化足够小或者超出最大
迭代次数时，模拟过程终止，并进入步骤S37；
[0030] 步骤S37：分离出林冠点云(非地面点云)和地面点云；计算格网点和相应点云点之间的距离，对于点云点，如果该距离小于
阈值h，被分为地面点，否则被分为非地面点。
[0031] 进一步地，步骤S4具体包括以下步骤：
[0032] 步骤S41：分别对林冠点云和地面点云采用自然邻域插值法，生成数字表面模型DSM和数字高程模型DEM；自然邻域插值特征在于利用最邻近样本点，并基于区域大小按比例对这些样本应用权重来进行插值，与常规方法相比，能够较好的还原地形表面的局部特征；
[0033] 步骤S42：利用数字表面模型DSM和数字高程模型DEM做差值运算获得采伐区的冠层高度模型CHM。
[0034] 进一步地，步骤S5具体包括以下步骤：
[0035] 步骤S51：通过一个平均树冠尺寸的固定窗口探测样地内潜在的林冠顶点位置；
[0036] 步骤S52：采用自适应的动态窗口对获取的顶点进行判断，如当前顶点为对应窗口区域的最大值则保存，否则删除；其中动态窗口的大小通过计算潜在顶点八个剖面方向半方差值的变化值来自适应的确定(Wulder et al.,2000)，其公式计算如下：
[0037]
[0038] 式中，γ(h)为经验半方差值，xi为影像的像元位置，h为两个像元的空间分割距离，Z(xi)为对应影像xi处的像元值，N为在一定分隔距离下像元对的对数。
[0039] 进一步地，步骤S6具体包括以下步骤：
[0040] 步骤S61：利用已有的5种树高—胸径模型，反推得到5种胸径—树高模型；
[0041] 步骤S62：将已有的实测树高、胸径数据带入所述胸径—树高模型，利用决定系数R2、均方根误差RMES指标对胸径—树高模型的精度进行评价，从而选择出最优的胸径—树2 2
高模型；其中，R的正常取值范围为[0,1],R 越接近1，表明胸径—树高模型的拟合效果越好；其中RMES值越小越好，用拟合值与原始值之间的偏差来评价模型的回归效果；其中，R2与RMES的计算采用下式：
[0042]
[0043]
[0044] 其中，yi为原始值， 为模型拟合值，为原始值的平均值，n 为样本数；
[0045] 步骤S63：利用步骤S5的结果和选择的最优胸径—树高模型能够获得采伐区的单木胸径信息。
[0046] 进一步地，步骤S61中，5种胸径—树高模型分别为：
[0047] 线性回归模型(Linear regression model)：H＝a+bD；
[0048] 异速生长模型(Allometric growth model)：H＝aDb；
[0049] Meyer model：H＝a(1-e-bD)；
[0050] WeiBull model：H＝a(1-e-bDc)；
[0051] Bates and Watts model：
[0052] 较佳的，所述针叶林树种分类提取可采用
随机森林，
支持向量机等分类器实现，对于高度差异明显的树种可以直接用高度阈值的方法进行分类，例如
马尾松和杉木。
[0053] 进一步地，步骤S7中，所述二元立木材积公式为：
[0054]
[0055] 式中，Di为木材的胸径，Hi为木材的树高，Vi为木材的单木材积， k、hi、di均为系数。
[0056] 较佳的，所述二元立木材积公式以马尾松和杉木为例：
[0057] Vi＝0.00008721Di1.785388607Hi0.9313923697；
[0058] Vj＝0.0000942941Dj1.832223553Hj0.8197255549；
[0059] 式中，Di、Hi、Vi分别为杉木的胸径、树高和单木材积，Dj、 Hj、Vj分别为马尾松的胸径、树高和单木材积。
[0060] 与
现有技术相比，本发明有以下有益效果：本发明解决了地形起伏较大的针叶林区域单期影像不能提取高精度树高以及蓄积量信息的问题。通过采集伐区两期无人机影像，解决了树高参数精确估算的问题，建立胸径—树高模型估算胸径数据，可有效提高森林蓄积量的估算精度；有利于伐区森林资源调查和森林经营管理的信息化和智慧化。
附图说明
[0061] 图1为本发明
实施例的方法流程示意图。
[0062] 图2为本发明实施例的采伐前后无人机影像数据，图2中，(a)为采伐前UAV影像，(b)为采伐后UAV影像。
[0063] 图3为本发明实施例的两期点云匹配结果，图3中，左边为俯视图，右边为侧视图。
[0064] 图4为本发明实施例的布料模拟滤波处理后林冠点云和地面点云。图 4中，(a)为林冠点云，(b)为地面点云。
[0065] 图5为本发明实施例的基于林冠点云和地面点云自然邻域法插值生成数字表面模型(DSM)与数字高程模型(DEM)。图5中，(a)为DSM 模型，(b)为DEM模型。
[0066] 图6为本发明实施例的插值生成的DEM与实测数据对比。
[0067] 图7为本发明实施例的采用DSM和DEM差值生成的冠层高度模型 (CHM)。
[0068] 图8为本发明实施例的CHM自适应滤波结果，其中，(a)是固定窗口提取林冠顶点结果，(b)是动态窗口提取林冠顶点结果。
[0069] 图9为本发明实施例的采伐区马尾松杉木分类与顶点提取结果。
具体实施方式
[0070] 下面结合附图及实施例对本发明做进一步说明。
[0071] 应该指出，以下详细说明都是示例性的，旨在对本
申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
[0072] 需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本
说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
[0073] 如图1所示，本实施例提供了一种基于两期无人机(Unmaned Aerial Vehicle：UAV)影像的针叶林采伐蓄积量估算方法，具体包括以下步骤：
[0074] 步骤S1：利用无人机获取伐区采伐前后两期的航拍影像数据，经处理得到两期镶嵌影像和影像点云数据；
[0075] 步骤S2：对两期影像点云数据进行匹配；
[0076] 步骤S3：对步骤S2匹配后的数据进行布料模拟算法滤波处理，分离出林冠点云和地面点云；
[0077] 步骤S4：基于步骤S3获取的点云数据，利用自然邻域插值法生成数字表面模型和数字高程模型，进而获取冠层高度模型；
[0078] 步骤S5：基于步骤S4中构建的冠层高度模型，采用动态窗口局域最大值法进行林冠顶点探测，从而得到株数信息与对应树高信息；
[0079] 步骤S6：基于实测数据研建胸径-树高回归模型，利用步骤S5 的结果获取对应的单木胸径信息；并利用影像抽样或监督分类获得针叶林树种组成比例信息；
[0080] 步骤S7：利用二元立木材积公式计算相应树种的单木材积，结合株数信息得到采伐蓄积量。
[0081] 在本实施例中，步骤S1中，所述无人机飞行设置的航向重叠度优于80％，旁向重叠度优于60％；所获的无人机影像分辨率优于10cm；点云密度高于10/m3。
[0082] 在本实施例中，步骤S2具体为：对两期影像点云数据选取多个同位点进行配准，数据的匹配误差要求在两个像元之内，保证研究区域中每棵树木点云的位置与相应地面点云位置完好对应。
[0083] 在本实施例中，步骤S2进一步具体为：将获得的两期影像点云数据看做从P,Q两个不同视角下获得的具有重叠区域的数据点集，若在数据点集P中提取一个子集{mi|mi∈P,i＝1,2,...,N}，在数据点集Q中有一子集{mi′|mi′∈Q,i＝1,2,...,N}与之各点对应，配准就是求解在不同视角(坐标系)下，三维数据点之间的坐标转换关系R和T，其中R为旋转矩阵，T为平移矩阵，使得目标函数达到最小化；即多视角对齐的实质是计算满足如下目标函数的旋转和平移变换矩阵：
[0084] f(R,T)＝min∑[R·pi+T-qi]2；
[0085] 式中，pi与qi为需要对齐的点。
[0086] 在本实施例中，步骤S3具体包括以下步骤：
[0087] 步骤S31：转换点云几何坐标，使配准之后的点云数据翻转；
[0088] 步骤S32：初始化布料格网，通过格网分辨率确定格网节点数量；
[0089] 步骤S33：将点云点和格网点投影到同一水平面，寻找每个格网点对应的点云点，并记录对应点云点的高程值；
[0090] 步骤S34：计算格网点受重力作用位移的位置，并比较该节点位置高程与其对应点云点高程，如果该节点位置高程小于或者等于点云点高程，则将该节点位置替换到对应点云点位置，并标记为不可动点；其中格网点受重力作用位移的位置计算公式为：
[0091]
[0092] 式中，m是布料的质量，被设为常数1；X是某个时刻的节点位置；Δt是时间步长，G是重力常数；给出时间步长和节点初始位置，就能够得到当前点的位置；
[0093] 步骤S35：计算每个格网点受邻近节点影响而移动的位置，每个布料点的修正位移计算公式为：
[0094]
[0095] 式中， 为节点的移动向量； 为待移动节点的当前位置； 为 的邻近节点位置；为垂直向上的标准向量， b是判断节点是否移动的参量，当节点可移动时，设b＝1，否则b＝0；
[0096] 步骤S36：重复步骤S34和步骤S35，当所有节点的最大高程变化足够小或者超出最大迭代次数时，模拟过程终止，并进入步骤S37；
[0097] 步骤S37：分离出林冠点云(非地面点云)和地面点云；计算格网点和相应点云点之间的距离，对于点云点，如果该距离小于阈值h，被分为地面点，否则被分为非地面点。
[0098] 在本实施例中，步骤S4具体包括以下步骤：
[0099] 步骤S41：分别对林冠点云和地面点云采用自然邻域插值法，生成数字表面模型DSM和数字高程模型DEM；自然邻域插值特征在于利用最邻近样本点，并基于区域大小按比例对这些样本应用权重来进行插值，与常规方法相比，能够较好的还原地形表面的局部特征；
[0100] 步骤S42：利用数字表面模型DSM和数字高程模型DEM做差值运算获得采伐区的冠层高度模型CHM。
[0101] 在本实施例中，步骤S5具体包括以下步骤：
[0102] 步骤S51：通过一个平均树冠尺寸的固定窗口探测样地内潜在的林冠顶点位置；
[0103] 步骤S52：采用自适应的动态窗口对获取的顶点进行判断，如当前顶点为对应窗口区域的最大值则保存，否则删除；其中动态窗口的大小通过计算潜在顶点八个剖面方向半方差值的变化值来自适应的确定(Wulder et al.,2000)，其公式计算如下：
[0104]
[0105] 式中，γ(h)为经验半方差值，xi为影像的像元位置，h为两个像元的空间分割距离，Z(xi)为对应影像xi处的像元值，N为在一定分隔距离下像元对的对数。
[0106] 在本实施例中，步骤S6具体包括以下步骤：
[0107] 步骤S61：利用已有的5种树高—胸径模型，反推得到5种胸径—树高模型；
[0108] 步骤S62：将已有的实测树高、胸径数据带入所述胸径—树高模型，利用决定系数R2、均方根误差RMES指标对胸径—树高模型的精度进行评价，从而选择出最优的胸径-树高模型；其中，R2的正常取值范围为[0,1],R2越接近1，表明胸径—树高模型的拟合效果越好；其中RMES值越小越好，用拟合值与原始值之间的偏差来评价模型的回归效果；其中，R2与RMES的计算采用下式：
[0109]
[0110]
[0111] 其中，yi为原始值， 为模型拟合值， 为原始值的平均值，n 为样本数；
[0112] 步骤S63：利用步骤S5的结果和选择的最优胸径—树高模型能够获得采伐区的单木胸径信息。
[0113] 在本实施例中，步骤S61中，5种胸径—树高模型分别为：
[0114] 线性回归模型(Linear regression model)：H＝a+bD；
[0115] 异速生长模型(Allometric growth model)：H＝aDb；
[0116] Meyer model：H＝a(1-e-bD)；
[0117] WeiBull model：H＝a(1-e-bDc)；
[0118] Bates and Watts model：
[0119] 较佳的，在本实施例中，所述针叶林树种分类提取可采用随机森林，支持向量机等分类器实现，对于高度差异明显的树种可以直接用高度阈值的方法进行分类，例如马尾松和杉木。
[0120] 在本实施例中，步骤S7中，所述二元立木材积公式为：
[0121]
[0122] 式中，Di为木材的胸径，Hi为木材的树高，Vi为木材的单木材积， k、hi、di均为系数。
[0123] 在本实施例中，所述二元立木材积公式以马尾松和杉木为例：
[0124] Vi＝0.00008721Di1.785388607Hi0.9313923697；
[0125] Vj＝0.0000942941Dj1.832223553Hj0.8197255549；
[0126] 式中，Di、Hi、Vi分别为杉木的胸径、树高和单木材积，Dj、 Hj、Vj分别为马尾松的胸径、树高和单木材积。
[0127] 特别的，本实施例以下结合一个实例对本发明进行详细介绍。本实施例选择的针叶林采伐区位于福建省将乐县邓坊
林场管理区，海拔 215米，坡度28°，伐区面积为22亩,主要树种为马尾松和杉木，采伐株数为366株，实测蓄积量为199/m3。
[0128] 采伐前使用瑞士SenseFly公司的小型固定翼eBee无人机获取待伐区无人机影像。无人机配备SONYWX220相机，1800万
像素，等效焦距为24mm。采伐前无人机影像采集时间是
2016年4月14日下午，飞行设置的航向重叠率为80％，旁向重叠率为60％，共获得313张有效照片，镶嵌和正射处理后获得影像
覆盖面积约为0.5平方千米，影像正射镶嵌后分辨率为
0.065米(图2中(a)所示)。小班采伐后使用大疆精灵3无人机获取采伐后小班影像，拍摄时间为2017年5月 16日上午，影像重叠度90％，共获得162张有效照片，处理后的正射影像覆盖面积约为0.19平方千米，影像正射镶嵌后分辨率为0.046 米(图2中(b)所示)。通过两期的无人机影像(图2)，可以清晰的看到森林小班区域采伐前的林冠和采伐后的地形，其中边界为采伐区矢量边界。
[0129] 为了消除地形对后续树高提取的影响，对两期点云数据选取6个同位点进行匹配(图3)，配准精度小于0.1m，保证研究区域中每棵树木点云的位置偏差控制在0.1米范围内。
[0130] 图4是进行布料模拟滤波(CSF)算法处理之后分离出的林冠点云数据和地面点云数据。实施例中设置的布料模拟滤波参数为：格网分辨率为0.1m；布料松紧程度分为三个级别，平坦地形硬度为1，略微起伏地形硬度为2，陡峭地形硬度为3。根据实验地区的实际情况，选为3；距离阈值设置为0.5m,用于地面点与和非地面点分类；迭代次数，通常设置为100-150以内即可满足滤波要求，本次实验设置为 100。
[0131] 图5是利用自然邻域插值法对地面和林冠点云处理得到的数字高程模型(DEM)与数字表面模型(DSM)；图6是插值生成的高程与实测的高程数据对比，均方根误差为0.49m，满足冠层高度模型对DEM 数据的要求。图7是利用DSM与DEM的差值运算得到的冠层高度模型 (CHM)。
[0132] 图8为CHM自适应滤波结果。其中，图8中(a)是应用固定窗口局域最大值法进行林冠顶点检测的结果，存在部分林冠检测到多个顶点的问题；图8中(b)是在固定窗口检测结果
基础上，应用自适应窗口最大值法进行林冠顶点检测的结果，可以发现消除了部分林冠出现多个顶点的现象。实验提取采伐株数为339，实际的采伐株数为 366，株数提取精度为93％。
[0133] 基于获取的高精度的冠层高度模型(CHM)，以及单木的树高(H) 信息，利用已有的5个树高—胸径模型(线性回归模型、异速生长模型、Meyer模型、WeBll模型、Bates and Watts模型)反推获得5 个胸径—树高模型(D-Hmodel)，结合实测400组树高胸径数据进行实验，其中利用150组马尾松和150组杉木数据进行实验发现马尾松的树高和胸径表现出了较强的相关性R2基本在0.82以上，而杉木的树高和胸径的相关性较弱，R2基本在0.60左右。
根据R2最大，RMSE 最小的标准，试验中选取异速生长模型为马尾松的D-H回归模型， WeiBull模型为杉木的D-H回归模型。
[0134] 图9为马尾松与杉木分类结果与顶点检测结果。经实地勘测发现，该地的马尾松树高均在16.5m以上，杉木树高均在16.5m以下，以该数值为阈值可将马尾松与杉木区分，结合D-H回归模型，得到相应的胸径信息，代入福建省人工林二元立木材积公式中：
[0135] Vi＝0.00008721Di1.785388607Hi0.9313923697
[0136] Vj＝0.0000942941Dj1.832223553Hj0.8197255549
[0137] 将所有单木材积相加得到该采伐区的总蓄积量为182/m3。而伐区调查显示实测蓄积量为199m3。比较结果显示本次试验的估算的采伐蓄积量的精度为91％，满足林业资源调查的精度要求。
[0138] 以上所述仅为本发明的较佳实施例，凡依本发明申请
专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。
