一种基于伞骨法与冠高比的树木冠层结构信息提取方法

摘要

本发明公开了一种基于伞骨法与冠高比的树木冠层结构信息提取方法，该度量方法包括如下步骤：步骤1)基于单次回波点云数据构建树冠高度模型，即在数字地面模型的基础上滤除非植被信息，获得包含仅植被分布及高度信息的基准图；步骤2)基于活动窗口的局部最大值搜索法，从所述基准图中提取单株立木的顶点；步骤3)在树顶检测结果基础之上，采用“伞骨算法”进行树冠投影边界识别，分离单体树木并确定各自树冠的边界位置，得到树冠边界识别结果；步骤4)结合树顶检测结果与树冠边界识别结果，基于伸缩圆与冠高比的树冠体积计算方法生成单体树冠体积分布图。步骤5)通过计算图像平均梯度反演冠层孔隙率，修正因树冠内部枝叶间存在的大量空隙而造成的计算误差，进一步提升步骤4)中的树冠体积反演精度。

说明书

技术领域

本发明涉及一种基于伞骨法与冠高比的树木冠层结构信息提取方法，特别 是一种应用机载激光雷达系统检测树冠边界和度量树冠体积的方法。

背景技术

绿地是城市中唯一有生命的基础设施，也是城市生态系统中的核心组成部 分。已有多项研究证明城市绿地在改善城市生态环境、促进人体健康等方面起 着不可或缺的积极作用。

遥感作为新型对地观测技术，具有大面积同步覆盖、实时连续、经济有效 等优点。高分辨率遥感与激光雷达等多源遥感数据为空间探测提供了新的视角 和研究方向。其中，LiDAR拥有的全自动、高精度立体扫描技术，使快速直接获 取地形表面模型成为可能，尤其适合植被等具备三维空间结构且立体形态不规 则的信息获取，正逐步成为推动城市绿地研究向立体拓展的重要载体。

周坚华等利用彩红外航片和计算机模拟技术，通过度量城市植物茎叶所占 空间体积估算城市绿地对城市的绿化效益。周延刚等利用彩红外航空影像建立 了植被高度模型、冠径-冠高关系模型和三维绿量计算方程。韦雪花等提出树 冠不能用一个或多个实心规则几何体来描述，提出了“体元模拟法”。

Hecht R等使用机载LiDAR数据估算城市植被体积。Liu和Shen等使用机 载LiDAR点云数据基于坡度过滤法和地面增长算法提取城市树木。Qi Chen等 基于机载LiDAR数据研发了一套实现单株树高、树冠面积、树冠率及叶面积等 提取的软件TiFFs。Wei Yao借助全波形机载LiDAR数据实现了树种分类与树冠 结构信息提取。

目前，国内外基于LiDAR的植被研究多以单体植被信息提取为主，融合多 源数据，在实现城市植被三维立体信息提取基础上研究其空间分布特征的较少。

考虑到以上问题，本发明针对针对城市绿化分布零散、结构多样，和以复 杂城市地物为背景的三维绿量提取技术研究不足、自动化弱、精度偏低等问题， 尝试建立一套基于机载LiDAR单次回波数据的单株立木识别与冠层结构信息提 取方法，用于解决不同冠层体积度量算法适应性差、城市大范围绿度空间提取 及度量效率差等难题，为更有效地提取诸如城市等大范围的植被冠层结构信息 提供技术参考。

发明内容

针对城市绿化分布零散、结构多样，和以复杂城市地物为背景的三维绿量 提取技术研究不足、自动化弱、精度偏低等问题，本发明提供一种无需伐倒立 木的可用于实现不同单体树种的树冠边缘检测、树冠体积估算、冠层结构信息 提取的技术流程。

本发明的目的通过以下技术步骤实现：

步骤1)基于单次回波点云数据构建树冠高度模型，即在数字地面模型的基 础上滤除非植被信息，获得包含仅植被分布及高度信息的基准图；

步骤2)基于活动窗口的局部最大值搜索法，从所述基准图中提取单株立木 的顶点(最高点)；

步骤3)在树顶检测结果基础之上，采用“伞骨算法”进行树冠投影边界识 别，分离单体树木并确定各自树冠的边界位置，得到树冠边界识别结果；

步骤4)结合树顶检测结果与树冠边界识别结果，基于伸缩圆与冠高比的树 冠体积计算方法生成单体树冠体积分布图。

步骤5)通过计算图像平均梯度反演冠层孔隙率，修正因树冠内部枝叶间存 在的大量空隙而造成的计算误差，进一步提升步骤4)中的树冠体积反演精度。

进一步，所述步骤1)的具体方法为：

a)基于改进的分水岭算法对航空影像进行边缘特征提取；b)基于航空影 像的红色和近红外波段计算得到NDVI；c)利用最大类间差算法(OTSU)确定 NDVI最佳分割阈值；d)根据阈值，划定植被区域；e)借助投票法修正二维植 被信息提取结果。

进一步，所述步骤2)的具体方法为：

a)基于一定数量的树冠高度模型采样点，构建树高-冠径关系模型；b) 根据已建立的关系模型确定逐像素遍历过程中的搜索范围；c)依据搜索范围， 判断像素是否满足树顶条件；d)依次遍历图像中像素点，直至找出全部树顶点。

进一步，所述步骤3)的具体方法为：

a)生成一系列以树顶点为对称中心的等旋转角的辐射状线段；b)为中心 点为起始点，沿方向的辐射状线段向外围遍历，依次判断位于同一条辐射状线 段上相邻单元的高度差，若高差值满足判定规则，即判断该外围像素点为树冠 边界点；c)重复b)过程直至获得全部截断点，通过逆时针依次连接各截断点， 最终形成封闭多边形，形成单体树冠轮廓线。

进一步，所述步骤4)的具体方法为：

a)在树冠边界识别结果的基础上，利用伸缩圆算法获得椭球体树冠结构的 长半轴与短半轴参数；b)结合树顶检测结果的高度信息，通过冠高比概念计算 得到极半径参数；c)经由椭球体积计算公式最终获得树冠体积的估算量。

进一步，所述步骤5)的具体方法为：

a)依次遍历每个树冠单元，计算其平均梯度值；b)获取上一步计算结果 中的最大值与最小值，采用线性函数转换法，将平均梯度值归一化并作二次转 化计算，得到每株立木的冠层孔隙率修正系数；c)将步骤4)中的树冠体积计 算结果乘以该系数，获得修正后的单体树冠体积值。

附图说明

图1为影像分割结果；

图2为树冠高度模型；

图3为树顶检测结果；

图4为树冠边界识别结果；

图5为单体树冠体积分布图；

具体实施方式

下面结合附图对本发明“一种基于伞骨法与冠高比的树木冠层结构信息提 取方法”作进一步阐述说明。

本发明所涉及的“基于伞骨法的树冠投影边界识别算法”和“基于伸缩圆 与冠高比的树冠体积计算方法”，是城市植被空间结构参量提取的重要技术创 新。其中，伞骨法具有“分离邻接树冠的相交界限”与“识别结果准确且符合 树冠形态”的优点，同时提升了算法的运行效率，为区域尺度的树冠投影边界 识别提供一种可行的方法；基于伸缩圆与冠高比的树冠体积计算方法则有效解 决了不同冠层体积度量算法适应性差、城市大范围绿度空间提取及度量效率差 等难题，可为同领域的相关研究者借鉴使用；冠层孔隙率修正因子可用于修正 因采用几何实心体模拟树冠形态而造成的计算误差，使树冠体积反演结果具备 一定的可靠性。

(一)树冠高度模型构建

采取“先二维，后三维”的操作步骤获得树冠高度模型。首先，基于高分 辨率遥感图像的面向对象图像分类方法，同时充分考虑高分辨率航空影像提取 边界信息的优势，采用改进的分水岭分割算法进行影像边缘特征提取，得到标 记后图像分割结果(图1)。然后基于航空影像的红色和近红外波段计算得到 NDVI，利用最大类间差算法(OTSU)确定NDVI最佳分割阈值，进行植被提取。 在完成图像分割和阈值提取的基础上，借助投票法修正二维植被信息提取结果。 具体步骤为：遍历图像分割结果中的每个对象；统计每个对象中植被像元累计 个数占整个对象像元总数的百分比；若该对象内植被像元所占百分比超过50％， 即判断该对象为植被。经由随机样本检验，基于以上步骤得到的提取结果，精 度优于95％。再将生成的植被二值掩膜图像与数据预处理步骤中得到的DHM图层 作叠加，最终计算得到树冠高度模型(图2)。

(二)树顶检测

树顶通常是指在一定聚簇范围内植被的最高点，同时认定树冠直径与树高 存在一定的相关性。本发明基于活动窗口的局部最大值搜索法，通过构建树高 -冠径的关系假设实时推算出窗口变量大小。首先，将图像等分为N×N的格网 大小，依次从每个格网中随机选择任意一棵树木，量算其横纵向树冠直径长度 的平均值和最大高度值。随着树高的增加，树冠直径也随之扩大，二者基本呈 现线性关系，由此建立树高-冠径关系模型。该模型将用于逐像素遍历过程中 确定搜索范围。具体步骤为：

1)选择树冠高度模型图层的左上角第一个像素点，其值为树高TH；

2)经由树高-冠径关系公式计算得到冠径CD；

3)以被选点为中心，CD/2/2个像素点为半径，确定搜索范围；

4)如果被选点是搜索范围内的最高点，则标记为树顶点，并将其位置信息 和高度信息存入列表；

5)重复以上过程，依次遍历图像中的全部像素点，直至找出全部树顶点。

基于上述步骤所得树顶检测结果如图3所示。

(三)树冠投影边界识别

鉴于树冠高度模型上的树冠高度值在一定区域范围内具有稳定均一的特 质，若其变化幅度超出既定阈值，便可认定其为边界所在处。伞骨算法通过构 建一系列以树顶点p为对称中心且旋转角相等的辐射状线段，每条辐射状线段 长度为M个像素，则以不同像素点为中心的一系列辐射状线段可定义为 W(p，n，m)。由于树冠边界形状多样且无法预测，故设树顶点为中心p，并作为起 始点，顺序沿各方向的辐射状线段向外围遍历，通过判断位于同一条辐射状线 段上前一像素h_i与后一像素h_i+1的高度值关系，决定是否停止搜索。停止条件如 下：①h_i+1＝0；②h_i-h_i+1≥h_t/3；③h_i-h_i+1＜0(i＝1，2，…，n)。为防止因检测半 径过长等不合理结果，还可加入绝对限制条件：n＝h_t/3，即n的大小取决于研究 区内冠径与树高关系的实际情况。通过逆时针依次连接各截断点，最终形成封 闭多边形，实现单体树冠轮廓识别，进而获得的树冠投影边界识别结果(图4)。

(四)树冠体积计算

伸缩圆算法是在树冠边界识别结果基础上，以树冠投影多边形的重心为圆 心，由小及大，生成不同半径长r的检测圆C_r。当检测圆半径r由1开始，以 单个像素为步长累加的过程中，当首次发生S_i时，对应的r₁值即为 椭球体的短半轴b；当首次发生时，对应的r₂值即为椭球体的长 半轴a。依次遍历完全部树冠识别结果，将生成一系列与树冠对应的椭圆长、短 半轴值。

冠高比是冠高与树高的比值。假设树冠的几何形态符合椭球特征，即椭球 割面积最大值出现在冠高的二分之一处，该椭圆面的边界即为树冠边界，而落 在该边界点上的像素高度平均值h_c已知，结合边界点高度与树顶高度h_t的几何 关系，经由公式c＝h_t-h_c，便可计算求出椭球体的极半径c。最后代入到椭球体 积计算公式V＝4πabc/3，便可得到树冠所占据的空间几何体积。

(五)冠层孔隙率修正系数

截至上一步骤，推算所得树冠体积计算值仅仅衡量了其外缘廓线围起来的 几何实心体所占据的空间范围，与“绿化三维量”的定义仍存在细微差别，这 项误差来自树冠内部枝叶间的大量空隙。基于“树冠枝叶数量越多，相邻叶片 间距越小，树冠整体就越紧实，反映在图像上便表现为冠层上表面纹理越平滑” 的假设，本发明通过计算图像平均梯度反演冠层孔隙率，平均梯度定义如下：

$\bar{G}=\frac{1}{M\times N}\underset{i=0}{\overset{M-1}{\Sigma }}\underset{j=0}{\overset{N-1}{\Sigma }}\sqrt{\frac{{\left(\frac{\partial f(i,j)}{\partial x}\right)}^2+{\left(\frac{\partial f(i,j)}{\partial y}\right)}^2}{2}}$

式中，f(i，j)代表由M行N列组成的二维矩阵中位于第i行、第j列的像 素值；与分别表示在像素点(i，j)沿x和y方向的一阶导数，用于体现 邻域数值间的突变程度。用差分来近似代替导数，则在点(i，j)处沿x和y方向 的差分可表示为：

$\left(\begin{array}{cc}\frac{\partial f(i,j)}{\partial x}\approx f(i+1,j)-f(i,j)& \frac{\partial f(i,j)}{\partial y}\end{array},\approx f(i,j+1)-f(i,j)\right)$

进一步采用线性函数转换法，将平均梯度值归一化，使之成为0-1之间的 小数。获取上一步计算结果中的最大值与最小值并依据下方公式对 结果作二次计算，即得到冠层孔隙率修正系数ε。

$ϵ=1-\frac{{\bar{G}}_k-{\bar{G}}_{\min }}{{\bar{G}}_{\max }-{\bar{G}}_{\min }},(k=\mathrm{1,2},...n)$

对研究区内全部树冠体积计算值乘以该系数，获得修正后的单体树冠体积 值(图5)。