本发明涉及一种联合低密度机载lidar点云和同步遥感影像的森林单木树高估算方法。
背景技术：
森林树高信息作为森林结构参数中最为重要的参数之一，是反演估测其他森林参数的基础。目前，国内外学者对利用激光雷达数据来提取森林树高信息作了大量研究，但由于受到点云采样密度和森林生境的影响，多集中于高密度点云数据对树高信息的提取，而缺乏对低密度点云数据的研究。为此，本发明提出将lidar数据和同步获取的高分辨率航空影像相结合，充分挖掘和利用低密度激光雷达数据提供的林冠高程信息和航空影像提供的单木林冠结构信息。
树冠形状通常具有一定的形状特征，如针叶林外观上呈现出圆锥形的特点，而阔叶林林冠则表现为顶部比较平坦的圆形伞状形态。基于林冠这一几何形状特征，pollock(1996)等提出一种从航空摄影所获取的单目影像中重构林冠的成像模型。但这种成像方法需要假设地形是水平的，且所用的单目影像并不能够很好的重建林冠的三维表面。sheng等人(2002)为了实现针叶林树冠表面的三维重建，提出了5参数的三维树模型(见附图2)。图中(xt,yt,zt)为林冠顶点的地面坐标；bh为树的基部高度；ch为林冠深度；cr为林冠半径；cc为林冠曲率的调节系数。一旦这五个参数确立，便可以利用下面的方程来模拟出林冠表面上任一点处对应的坐标值。
本发明中所采用的三维树模型正是这一模型在机载lidar点云数据上的拓展应技术。目前已有相关研究将这一模型成功的运用到森林参数的估测中来。如persson和morsdorf等(2004)利用该模型来从高密度lidar数据中获取树高、冠幅和冠基高度等参数，其树高估测的平方根误差(rmse)达到0.6米左右。而paris等人(2013)则首次将该模型用于高采样密度lidar数据和航拍影像联立进行单木水平树高的反演，并取得了不错的估测效果。但目前该研究主要局限于地形平坦且树木相对独立地区，对于郁闭度较高且坡度较陡的山地地区以及针对低采样密度的机载lidar数据则缺乏相应的解决方案。
技术实现要素：
本发明的目的在于提供一种联合lidar点云和同步遥感影像的森林单木树高估算方法，将低密度机载lidar数据与同步获取的航空影像相结合，通过采用三维参数树模型的方法来实现对林分单木水平树高信息的估算。
为实现上述目的，本发明采用以下技术方案：一种联合lidar点云和同步遥感影像的森林单木树高估算方法，其特征在于，包括以下步骤：步骤s1：获取森林区域低密度机载lidar点云数据和同步高分辨率航空遥感影像步骤s2：通过对目标区同步航空遥感影像进行镶嵌和正射校正，获得目标区数字正射影像图dom，该影像与低密度机载lidar点云数据的匹配误差要求小于a；步骤s3：从数字正射影像图中获取单木林冠位置、林冠边界信息以及单木冠幅信息；步骤s4：对研究区的机载lidar点云数据进行处理，结合获取的单木林冠信息得到对应单木林冠内的点云数据；步骤s5：基于s3和s4步骤结果，构建三维树高模型；步骤s6：基于步骤s5构建的三维树高模型，根据单木林冠内激光点的个数大于1、等于1和无激光点三种情况进行三维树模型的优化和重建，从而计算得到对应林冠顶点的高度值ztop。
在本发明一实施例中，所述步骤s2中a为0.5m。
在本发明一实施例中，所述步骤s3的具体方法如下：步骤s31：基于dom影像的绿光波段，通过一个固定窗口探测样地内潜在的林冠顶点位置，然后采用自适应的动态窗口对获取的潜在顶点进行判断，如当前顶点为对应窗口区域的最大值则保存，否则删除；步骤s32：采用基于标记控制的分水岭分割方法来从dom影像中精确的提取单木林冠边界信息；步骤s33：利用寻谷法从影像中估测每个顶点的冠幅信息。
进一步的，所述动态窗口大小通过计算潜在顶点八个剖面方向半方差值的变程值来自适应的确定。
在本发明一实施例中，所述步骤s4的具体方法如下：步骤s41：采用基于渐进三角网tin算法进行点云滤波处理；步骤s42：对植被点进行高程归一化处理以获取到相对准确的林冠表面点云数据，即通过对获取的样地区域地面点进行tin插值得到dem高程，将植被点的高程减去dem高程值，从而获取相对准确的地物形态和高度信息。
进一步的，步骤s41包括以下具体步骤：提取格网内最低点作为地面种子点构建初始tin模型，进而对格网内点云按高程进行排序，然后通过不断向上加密三角网的方法提取地面点，最终完全分离出植被点云和地面点云。
在本发明一实施例中，所述步骤s5的具体方法如下：步骤s51：基于s3和s4参数提取结果，构建一个椭球体来描述针叶林的林冠几何外形，其中(xtop,ytop,ztop)为对应的林冠顶点坐标，cc为林冠表面曲率的调节系数，ch为林冠深度，bh为林冠的冠基高度，cr为林冠冠幅半径；所述林冠椭球体主要由位于林冠表面的激光点云数据来进行拟合得到，即用每个林冠内对应的lidar数据的三维坐标来唯一确定一个林冠三维树模型，每个林冠包络用如下的公式(1)来表示：
式中z的高程区间满足于ztop-ch<z<ztop，(x,y,z)为对应的林冠表面激光点的坐标值；步骤s52：确定林冠顶点的平面坐标值(xtop,ytop)和对应林冠的冠幅半径cr，设定一组固定的cc和ch参数值。
在本发明一实施例中，当单木林冠内激光点的个数大于1时三维树模型的优化和重建，林冠为d，包括以下步骤：步骤s611：林冠内每一个激光点都可以构建一个方程来计算得到对应的高度值ztop；步骤s612：通过计算每个林冠对应的最优树模型来估测林冠高度值，即利用林冠内激光点云数据来拟合得到一个与林冠点云最接近的三维树模型；利用拟合得到的三维树模型与激光点云之间距离残差最小时的参数组合作为最优三维树模型，其每个林冠激光点对应的残差值由下面方程式表示：
式中rj(z'top)为当前林冠内第j个lidar数据的地面坐标(xj,yj,zj)代入模型方程所计算得到的残差值；步骤s613：对当前林冠内每个激光点都计算得到这样一个残差值，最后将林冠内所有lidar数据对应残差值的平方和达到最小值时的cc和ch参数值作为当前林冠的最佳三维树模型；步骤s614：采用一个最小二乘最优化方法计算林冠内所有lidar数据对应残差平方和取得最小值时的z′top值，将其作为当前最佳树模型下的树顶高度值；其计算公式如下：
针对所有的ch和cc组合，计算出一组残差矩阵和对应的z′top值，最终选择残差矩阵中最小值所对应的ch和cc组合构建当前树冠的最优三维参数模型，而对应的z′top值即为最优的树高值。
进一步的，当单木林冠内激光点的个数等于1时三维树模型的优化和重建，林冠为h，具体包括以下步骤：步骤s621：在d对应的林冠点云中查找与当前林冠内激光点高程相近的林冠点云；步骤s622：计算查找到的林冠点云到对应林冠顶点位置的距离d和当前林冠内激光点到林冠顶点的距离dref，距离计算公式如下：
步骤s623：在计算得到的所有林冠点云到相应林冠顶点的距离d中查找与dref最相似的前n个林冠的树模型；步骤s624：计算n个林冠对应的树高值的均值作为当前林冠h的树高值，并将n个林冠中高度与树高均值最相近的林冠三维树模型作为当前林冠h的三维树模型。
进一步的，当不存在单木林冠内激光点时三维树模型的优化和重建，林冠为l，具体包括以下步骤：步骤s631:计算当前林冠l的顶点与d情况下林冠中所有顶点之间的距离；步骤s632:对步骤s631获取的距离矩阵按升序进行排序，并选取其中前k个林冠顶点；步骤s633:将当前林冠l与选取的k个林冠的冠幅面积做差，然后对面积差值做升序排列，再选取其中前n个林冠作为与林冠l相似的林冠；
步骤s634：计算前n个林冠对应树高值的均值作为当前林冠l的树高值，并将n个林冠中高度与树高均值相近的林冠三维树模型作为当前林冠l的三维树模型。
本发明与现有技术相比具有以下有益效果：本发明有效解决了低采样密度机载lidar数据树高估算结果偏低而只能用于林分平均树高估算的问题，可将树高估算等级提升到单木水平，并有效提高林分单木树高估算的精度。
附图说明
图1是本发明技术流程。
图2是本发明所采用的三维树模型参数。
图3是本发明实施案例福建省将乐县针叶林树高估算各阶段结果图；其中图3a实施案例福建将乐县试验区机载lidar数据植被点云；图3b实施案例将乐县试验区同步dom影像；图3c实施案例树冠分割结果；图3d实施案例冠幅估测结果；图3e实施案例将乐县实验区单木林冠内植被点云分布图；图3f实施案例将乐林场实验区单木树高重建结果。
图4为本发明实施案例三维树模型重建树高值与实测树高值的对比图。
具体实施方式
为了让一般技术人员更好的理解本发明的技术方案，下面结合附图及具体实施例对本发明做进一步说明。
本发明采用如下技术方案：一种联合低密度机载lidar点云和同步遥感影像的森林单木树高估算方法方法，其包括以下步骤：
步骤s1：获取森林区域低密度机载lidar点云数据和同步高空间分辨率航空遥感影像
步骤s2：通过对目标区同步航空遥感影像进行镶嵌和正射校正，获得目标区数字正射影像图(dom)，该影像与低密度机载lidar点云数据的匹配误差要求小于a；
步骤s3：从数字正射影像图(dom)中获取单木林冠位置、林冠边界信息以及单木冠幅信息；
步骤s4：对研究区的机载lidar点云数据进行处理，结合获取的单木林冠信息得到对应单木林冠内的点云数据；
步骤s5：基于s3和s4步骤结果，构建三维树高模型；
步骤s6：在步骤s5构建的三维树高模型基础上，根据单木林冠内激光点的个数大于1、等于1和无激光点三种情况进行三维树模型的优化和重建，从而计算得到对应林冠顶点的高度值ztop。
本发明主要流向示意图参见图1。
所述低密度机载lidar数据点云密度通常小于1个点/m2，同步高分辨率航空遥感影像分辨率一般为0.2-0.5m。本发明将低密度机载lidar数据与同步获取的航空影像相结合，通过采用一个三维参数树模型的方法来实现对森林林分单木水平树高信息的估算。
将同步航空数字正射影像(dom)与低密度机载lidar点云数据精确匹配，匹配误差要求小于0.5m
进一步的，所述步骤s3的具体方法如下：
步骤s31：基于dom影像的绿光波段，通过一个较小的固定窗口探测样地内潜在的林冠顶点位置，然后采用自适应的动态窗口对获取的潜在顶点进行判断，如当前顶点为对应窗口区域的最大值则保存，否则删除。较佳的，动态窗口大小主要通过计算潜在顶点八个剖面方向半方差值的变程值来确定；
步骤s32：采用基于标记控制的分水岭分割方法来从dom影像中精确的提取单木林冠边界信息；
步骤s33：利用寻谷法从影像中估测每个顶点的冠幅信息。
进一步的，所述步骤s4的具体方法如下：
步骤s41：采用基于渐进三角网(tin)算法进行点云滤波处理。在本发明一实施例中，该算法提取格网内最低点作为地面种子点构建初始tin模型，进而对格网内点云按高程进行排序，然后通过不断向上加密三角网的方法提取地面点，这样最终完全分离出植被点云和地面点云。
步骤s42：对植被点进行高程归一化处理以获取到相对准确的林冠表面点云数据。即通过对获取的样地区域地面点进行tin插值得到dem高程，将植被点的高程减去dem高程值，从而获取相对准确的地物形态和高度信息。
进一步的，所述步骤s5的具体方法如下：
步骤s51：基于s3和s4参数提取结果，构建一个广义椭球体来描述针叶林的林冠几何外形。其中(xtop,ytop,ztop)为对应的林冠顶点坐标，cc为林冠表面曲率的调节系数，ch为林冠深度，bh为林冠的冠基高度，cr为林冠冠幅半径。本发明所采用的三维树模型如图2所示。与原模型不一样的是，该林冠椭球体主要由位于林冠表面的激光点云数据来进行拟合得到，即用每个林冠内对应的lidar数据的三维坐标来唯一确定一个林冠三维树模型。每个林冠包络可以用式1的数学表达式来描述：
式中z的高程区间满足于ztop-ch<z<ztop，这样可以保证参与计算的点云数据都位于林冠表面。(x,y,z)为对应的林冠表面激光点的坐标值。步骤s52：确定林冠顶点的平面坐标值(xtop,ytop)和对应林冠的冠幅半径cr，设定固定的一组cc和ch的参数值。
进一步的，所述步骤s6的具体方法如下：
步骤s61：林冠内大于一个激光点的树高重建，这里用d表示：
在此种情况下，林冠内每一个激光点都可以构建一个方程来计算得到对应的高度值ztop。为了避免用一组cc和ch的参数值来代替所有激光点拟合的林冠三维模型，本发明通过计算每个林冠对应的最优树模型来估测林冠高度值，即利用林冠内激光点云数据来拟合得到一个与林冠点云最接近的三维树模型。这里主要利用拟合得到的三维树模型与激光点云之间距离残差最小时的参数组合作为最优三维树模型，其每个林冠激光点对应的残差值由下面方程式表示：
式中rj(z'top)为当前林冠内第j个lidar数据的地面坐标(xj,yj,zj)代入模型方程所计算得到的残差值。然后对当前林冠内每个激光点都计算得到这样一个残差值，最后将林冠内所有lidar数据对应残差值的平方和达到最小值时的cc和ch参数值作为当前林冠的最佳三维树模型。本发明主要采用了一个最小二乘最优化方法(nls优化算法)计算林冠内所有lidar数据对应残差平方和取得最小值时的z′top值，将其作为当前最佳树模型下的树顶高度值。其计算公式如下：
针对所有的ch和cc组合，计算出一组残差矩阵和对应的z′top值，最终选择残差矩阵中最小值所对应的ch和cc组合为当前树冠模型的最优三维参数，而对应的z′top值即为最优的树高值。
步骤s62：林冠内仅有一个激光点的树高重建，这里用h表示：
当pt等于1，即林冠内只有一个激光点时，其对应的三维树模型方程的残差值为0，因而不能用计算最小残差平方和的方式来获取最优树模型。本发明中主要根据林冠内激光点与树顶点距离，在已获取树高林冠中查找相似的树冠来进行估算，其计算原理如图1所示，主要包括以下几个步骤：
(1)在d对应的林冠点云中查找与当前林冠内激光点高程相近的林冠点云；
(2)计算查找到的林冠点云到对应林冠顶点位置的距离d和当前林冠内激光点到林冠顶点的距离dref，距离计算公式如下：
(4)在计算得到的所有林冠点云到相应林冠顶点的距离d中查找与dref最相似的前n个林冠的树模型；
(5)计算n个林冠对应的树高值的均值作为当前林冠h的树高值，并将n个林冠中高度与树高均值最相近的林冠三维树模型作为当前林冠h的三维树模型。
步骤s63：林冠内不存在激光点云的树高重建，这里用l表示：
当识别的林冠内没有点云数据分布时，本发明对该类型树高的重建主要采用一个类似k-nn的方法来获取最佳三维树模型。即假设在一般的森林场景中，同一生境下的树木与其周边树木在树高和冠幅上都存在一定的相似性。因此，通过查找与对应林冠距离上比较接近的k个林冠，再根据林冠之间冠幅大小的相似性来估测当前林冠的树顶高度值。其估测过程主要包括以下几个步骤：
(1)计算当前林冠l的顶点与d情况下林冠中所有顶点之间的距离；
(2)对上一步获取的距离矩阵按升序进行排序，并选取其中前k个林冠顶点；
(3)将当前林冠l与选取的k个林冠的冠幅面积做差，然后对面积差值做升序排列，再选取其中前n个林冠作为与林冠l相似的林冠；
(4)计算前n个林冠对应树高值的均值作为当前林冠l的树高值，并将n个林冠中高度与树高均值相近的林冠三维树模型作为当前林冠l的三维树模型。
以福建省将乐县局部林区2014年12月机载lidar点云数据及同步获取的dom影像为例进行试验分析。lidar平均扫描点距小于2米，整个激光点云密度约为0.7个/m2(图3a)，同步获取的dom正射影像分辨率为0.5m(图3b)。
首先从数字正射影像图(dom)中获取单木林冠位置、林冠边界信息(图3c)以及单木冠幅信息(图3d)
在对将乐实验区进行单木树高重建时，先根据影像获取的林冠边界等信息将研究区内的单木林冠分为三种情况，然后分别进行单木树高的重建。如图3e所示为试验区所获取的单木林冠内激光点云的分布图。根据单木林冠内激光点的数量分为三种情况分别进行统计，其中，林冠内大于一个激光点的单木共计268棵；林冠内只有一个激光点的单木共计133棵；林冠点云缺失的单木共计184棵。
将乐实验区内单木树高重建的过程也分为三种情况分别进行。
(1)在对林冠内大于一个激光点的单木最优参数进行拟合时，基于研究区内树种主要为长势比较均一的人工杉木幼林，其高度普遍位于4～6米之间，且林冠几何外形基本一致，故直接设置ch和cc参数组合分别为：ch设置为1～3米，迭代计算的步长为0.5米；cc设置为1.1～1.9，迭代计算步长为0.1。
(2)对于林冠内仅有一个激光点和无激光点的情况，采用一个类似k-nn的方法来获取单木的最佳三维树模型。在林冠内仅有一个激光点时，设置点云高差为±0.5米，在林冠内有多个激光点的林冠内查找与当前林冠内激光点高程相似的林冠点云，然后计算林冠点云到对应顶点的距离，并选取距离差值在0.5米以内的林冠作为与当前林冠相似的林冠，并取相似林冠对应高度值的均值作为当前林冠高度值；
(3)对林冠内无激光点时，首先计算当前林冠顶点与有多个激光点林冠顶点的距离，然后将距离进行排序，并选择前10个顶点对应的冠幅值与当前林冠冠幅面积做差，对面积差值进行排序后选择前三个林冠最为相似林冠，并取相似林冠对应高度值的均值作为当前林冠高度值。
通过以上步骤，便完成对研究区内所有单木的树高重建工作，通过计算，最终得到将乐试验区内的单木树高值(图3f)。
为验证三维树高重建方法在将乐实验区的有效性与估测精度，根据野外实测的单木树高数据对重建的单木树高和直接利用机载lidar数据提取的树高进行对比精度分析。附表1所示为两种方法所获取树高与实测树高的精度评价。结果表明：直接从点云中提取树高的最低估测精度57.4％，最高为87.8％，均小于采用三维树高重建方法估算树高的精度，其最高达到了96.4％。在整体估测精度上，三维重建树高的平均估测精度达到83.6％，也大于直接点云提取的精度(75.24％)。此外，两种方法估测的结果与实测树高相比，树高值均偏低，但直接点云提取树高的均方根误差(rmse)为1.31米，而三维重建方法则相应减少为0.82米。因此，采用三维树高重建方法，根据林冠内的点云估测树高可以达到较好的估测精度，可以有效的解决机载激光雷达低点云密度时森林树高过低估测的问题。
表1将乐林场实验区树高重建结果与机载点云数据估测树高的精度
在前面的结果验证中，采用的三维重建树高结果均为林冠内大于一个激光点的单木，其在一定程度上能达到较高的估测精度。但对于林冠内只有一个激光点和林冠点云缺失的情况，采用一个类似k-nn的方法，从林冠内有多个激光点时获得的最优三维树模型中计算出相似的k棵树模型，进行树高估测。为此，为了更好的验证本发明所采用三维树高重建方法的估测性能与普适性，将树高重建过程中三种情况下的单木树高估测精度分别进行对比分析与验证。如图4所示为将乐实验区三种情况下单木树高与实测树高的散点图。从图中可以看出，三种情况下的树高估测值与对应实测树高的相关性都很明显，其决定系数(r2)均达到0.61以上，且随着林冠内激光点数量的不同，对应的树高估测值与实测值的相关关系也逐渐增强，呈现出从林冠点云缺失到有一个激光点再到大于一个激光点逐渐递进的发展趋势。而这一结论是与实际情况相符的，即林冠内存在的激光点个数越多，其获取树的几何三维模型与实际树木的三维外形越接近，从而极大地提高树高的估测精度。此外，虽然后两种情况下树高估测值与实测值的相关性相对较低，但整体树高估测水平还是达到与d情况下决定系数相近的80.8％，这也间接的表现出d情况下树高估测精度对树高重建模型的整体估测水平具有一定的影响。其次，在估测精度上h与l情况下的单木估测精度相似，前者为80.2％，后者为76.3％。这主要是由于该区域内森林长势基本一致，均为同一时期栽种的人工杉木，其林冠在三维几何外形上基本相似，因此，当获取的冠幅基本相似时，林冠内只有单一激光点对于重建三维树高的影响较小。
以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。