发明名称 一种基于全波形LiDAR单木冠层容积分解的树种分类方法
摘要 本发明公开了一种基于全波形LiDAR单木冠层容积分解的树种分类方法,包括:借助机载小光斑全波形LiDAR传感器进行数据采集;传感器记录每束激光脉冲返回的完整波形信息;LiDAR波形数据预处理;单木定位和冠幅提取;基于发射能量及传感器与地物的距离信息对LiDAR波形数据进行校正;构建体元框架并进行LiDAR波形的结构化分解;对单木进行冠层容积分解并提取特征变量;使用随机森林分类器进行树种分类。本发明验证结果表明,总体精度提升了11%左右;Kappa系数提升了0.1左右。
申请公布号 CN104808191A 申请公布日期 2015.07.29
申请号 CN201510234714.3 申请日期 2015.05.08
申请人 南京林业大学 发明人 曹林;许子乾;代劲松
分类号 G01S7/48(2006.01)I 主分类号 G01S7/48(2006.01)I
代理机构 江苏圣典律师事务所 32237 代理人 贺翔
主权项 一种基于全波形LiDAR单木冠层容积分解的树种分类方法,其特征在于,包括以下步骤:1)借助机载小光斑全波形LiDAR传感器进行数据采集;传感器记录每束激光脉冲返回的完整波形信息;2)LiDAR波形数据预处理A)噪声水平估计和数据平滑:首先把原始数据转换到频率域,再将频率较高的低值部分作为噪声水平的判断标准;然后选用高斯滤波器进行平滑。B)高斯拟合分解及波形数据点云化:基于回波数据是多个高斯函数的累加这一假设,对波形数据采用非线性最小二乘法进行拟合;然后通过局部最大峰值检测滤波算法从处理后的波形数据上提取离散点云,每个离散点中记录了返回信号的能量和振幅信息;C)生成数字地形:首先对从波形数据中提取出离散点云进行分类,然后对末次回波进行Kraus滤波处理用以去除非地面点,最后使用滤波后的末次回波数据并借助自然邻近法插值生成数字地形模型DTM;3)单木定位和冠幅提取A)对地面以上点云进行中值滤波,然后将点云中的高度信息栅格化生成数字表面模型DSM。将DSM减去数字地形模型,得到归一化植被高度CHM;B)通过局部最大值法确定单木树顶所在位置;C)单木冠幅的确定:单木冠幅的确定借助冠幅半径来进行描述,方法为首先以树顶为中心拟合16个半径方向上的冠幅剖面,然后计算到局部最小值的水平距离,最后将这些距离值进行水平方向上的平均从而得到冠幅半径。4)基于发射能量及传感器与地物的距离信息对LiDAR波形数据进行校正:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><msubsup><mi>W</mi><mi>i</mi><mi>c</mi></msubsup><mo>=</mo><msub><mi>W</mi><mi>i</mi></msub><mo>/</mo><msup><mi>W</mi><mi>e</mi></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000713943480000011.GIF" wi="1204" he="85" /></maths><maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><msubsup><mi>I</mi><mi>i</mi><mi>c</mi></msubsup><mo>=</mo><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mo>&CenterDot;</mo><msubsup><mi>s</mi><mi>i</mi><mi>k</mi></msubsup><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><msup><mi>I</mi><mi>e</mi></msup><mo>&CenterDot;</mo><msubsup><mi>s</mi><mn>0</mn><mi>k</mi></msubsup><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000713943480000012.GIF" wi="1328" he="85" /></maths>式中,<img file="FDA0000713943480000013.GIF" wi="65" he="78" />为校正后返回波内每个高斯波的波宽,<img file="FDA0000713943480000014.GIF" wi="60" he="75" />为校正后返回波内每个高斯波内每个高斯波的能量强度,W<sub>i</sub>为原始返回波内每个高斯波的波宽,W<sup>e</sup>为发射波的波宽,I<sub>i</sub>为原始返回波的能量强度,I<sup>e</sup>为发射波的能量强度,<img file="FDA0000713943480000015.GIF" wi="67" he="79" />传感器到反射物体的距离,<img file="FDA0000713943480000016.GIF" wi="68" he="81" />为标称距离,k为变化系数;5)构建体元框架并进行LiDAR波形的结构化分解首先利用三维体积单元划分地表以上的三维空间;然后在每个单元内,汇总穿过其中的最大脉冲能量并基于数字地形模型进行高度归一化;6)对单木进行冠层容积分解并提取特征变量首先借助体元框架,对样地内森林冠层内的脉冲信息进行三维分解;根据每个体元单元内的能量信息,分类为“填充”和“空缺”;“填充”体元单元又分为“透光层”和“微明层”;“透光层”的体元数量为“特征变量1”,占整个单木冠幅空间内体元总数的百分比为“特征变量2”,“透光层”的体元内的能量强度总和占整个冠幅空间内能量强度总和的比例为“特征变量3”;“微明层”的体元数量为“特征变量4”,占整个单木冠幅空间内体元总数的百分比为“特征变量5”,“微明层”的体元内的能量强度总和占整个冠幅空间内能量强度总和的比例为“特征变量6”;“填充”即透光层+微明层的体元数量为“特征变量7”,占整个单木冠幅空间内体元总数的百分比为“特征变量8”,“填充”内的能量强度总和占整个冠幅空间内能量强度总和的比例为“特征变量9”;“空缺”的体元数量为“特征变量10”,占整个单木冠幅空间内体元总数的百分比为“特征变量11”,“空缺”体元内的能量强度总和占整个冠幅空间内能量强度总和的比例为“特征变量12”;7)使用随机森林分类器进行树种分类A)随机森林分类是由很多决策树分类模型{h(X,Θ<sub>k</sub>),k=1,2,...}组成的组合分类模型,且参数集{Θ<sub>k</sub>}是独立同分布的随机向量,在给定自变量X下,每个决策树分类模型都由一票投票权来选择最优的分类结果;通过k轮训练,得到一个分类模型序列{h<sub>1</sub>(X),h<sub>2</sub>(X),…h<sub>k</sub>(X)},再用它们构成一个多分类模型系统,最终的分类决策函数为:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mi>H</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mi>arg</mi><munder><mi>max</mi><mi>Y</mi></munder><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>k</mi></munderover><mi>I</mi><mrow><mo>(</mo><msub><mi>h</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mi>Y</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000713943480000021.GIF" wi="1392" he="154" /></maths>式中,H(x)表示组合分类模型,h<sub>i</sub>是单个决策树分类模型,Y表示输出变量,I为示性函数;B)使用特征变量1‑12和随机森林分类器进行分类,并借助总体分类精度和卡帕系数这两个定量化指标对分类精度进行评价;卡帕系数计算公式:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><msub><mi>K</mi><mi>hat</mi></msub><mo>=</mo><mfrac><mrow><mi>N</mi><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>r</mi></munderover><msub><mi>x</mi><mi>ii</mi></msub><mo>-</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>r</mi></munderover><mrow><mo>(</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo></mrow></msub><msub><mi>x</mi><mrow><mo>+</mo><mi>i</mi></mrow></msub><mo>)</mo></mrow></mrow><mrow><msup><mi>N</mi><mn>2</mn></msup><mo>-</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>r</mi></munderover><mrow><mo>(</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo></mrow></msub><msub><mi>x</mi><mrow><mo>+</mo><mi>i</mi></mrow></msub><mo>)</mo></mrow></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000713943480000031.GIF" wi="1452" he="309" /></maths>式中,r为总的类别数,x<sub>ii</sub>为对角线上的像元数,x<sub>i+</sub>和x<sub>+i</sub>是列和行的总像元,N是总像元数。
地址 210037 江苏省南京市玄武区龙蟠路159号