发明名称 基于激光雷达综合波形模型反演森林生物物理特性的方法
摘要 一种基于激光雷达综合波形模型反演森林生物物理特性的方法,通过进行树冠分割和单木范围内的波形信息汇总,然后再汇总林分尺度上,使得反演尺度变换更为灵活合理;借助数字地形和高度误差阈值来判断地面返回脉冲,从而间接提升冠层返回脉冲的判别效果;从多个维度上提取综合波形特征信息,从而更加深入挖掘了波形数据中丰富的森林冠层信息。借助以上自主开发技术及创新内容,本发明在提升森林生物物理特性估算精度的同时也为其在机理上解释提供了理论及方法支撑。有效地避免了由于舍弃了原始波形数据中记录丰富的森林连续垂直分布信息从而大大降低了波形数据的应有价值的缺陷。
申请公布号 CN104180754A 申请公布日期 2014.12.03
申请号 CN201410362870.3 申请日期 2014.07.28
申请人 南京林业大学 发明人 曹林;代劲松
分类号 G01B11/02(2006.01)I;G01B11/28(2006.01)I;G01S17/89(2006.01)I 主分类号 G01B11/02(2006.01)I
代理机构 南京钟山专利代理有限公司 32252 代理人 戴朝荣
主权项 一种基于激光雷达综合波形模型反演森林生物物理特性的方法,其特征在于,步骤如下:步骤1:借助机载小光斑全波形LiDAR的传感器进行数据采集,数据采集中获得的LiDAR波形数据参数包括有脉冲发射频率、扫描频率和激光脉冲的光斑半径,由此所述的机载小光斑全波形LiDAR的传感器记录了每束激光脉冲返回的完整的LiDAR波形信息,所述的LiDAR带有高斯滤波器;步骤2:进行LiDAR波形数据的预处理,进行LiDAR波形数据的预处理的具体方式为首先把每束激光脉冲返回的完整的LiDAR波形信息导出到LiDAR的信息处理系统,LiDAR的信息处理系统就对该LiDAR波形信息进行噪声水平估计和数据平滑的处理,所述的噪声水平估计和数据平滑的处理的具体方式为把LiDAR波形信息转换到频率域,再将频率比设定的基准频率更高的低值部分作为噪声水平的判断标准,然后对转换到频率域的LiDAR波形信息进行还原,再用高斯滤波器进行平滑处理,然后进入高斯拟合(分解)及波形数据点云化的处理,所述的高斯拟合(分解)及波形数据点云化的处理的具体方式为对进行平滑处理后所得的波形数据采用非线性最小二乘法进行拟合,然后通过基于局部最大峰值的检测滤波算法从进行拟合处理后的波形数据中提取离散点云,每个离散点中记录了每束激光脉冲返回的完整的LiDAR波形信息的高度和能量强度信息;然后进入生成数字地形的处理,所述的生成数字地形的处理的具体方式为首先对从进行拟合处理后的波形数据中提取出的离散点云进行分类,然后对末次回波进行迭代选权滤波法处理用以去除非地面点,最后使用迭代选权滤波法后的末次回波数据并借助自然邻近法插值生成数字地形模型;步骤3:进行LiDAR波形数据校正,所述的进行LiDAR波形数据校正的具体方式为到LiDAR的信息处理系统根据公式(1)和公式(2)导出校正后返回波内每个高斯波的波宽和校正后返回波内每个高斯波内每个高斯波的能量强度<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>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000021.GIF" wi="1259" he="77" /></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>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000022.GIF" wi="1347" he="79" /></maths>上式中,<img file="FDA0000544430090000023.GIF" wi="68" he="71" />为校正后返回波内每个高斯波的波宽,<img file="FDA0000544430090000024.GIF" wi="68" he="71" />为校正后返回波内每个高斯波内每个高斯波的能量强度,W<sub>i</sub>为原始返回波内每个高斯波的波宽,W<sup>e</sup>为发射波的波宽,I<sub>i</sub>为原始返回波的能量强度,I<sup>e</sup>为发射波的能量强度,<img file="FDA0000544430090000025.GIF" wi="69" he="79" />传感器到反射物体的距离,<img file="FDA0000544430090000026.GIF" wi="68" he="79" />为标称距离,k为变化系数,i为正整数;进入单木容积内的原始波形信息汇总的处理,所述的单木容积内的原始波形信息汇总的处理的第一步为单木分割的处理,所述的单木分割的处理的具体方式为LiDAR的信息处理系统将地面以上点云插值生成数字表面模型DSM,然后数字表面模型DSM减去数字地形模型DTM从而得到归一化植被高度CHM,再通过高斯滤波器结合公式(3)对归一化植被高度CHM进行平滑处理后得到植被的高斯二维分布函数G(x,y),公式(3)如下所示:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mi>G</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><msqrt><mn>2</mn><mi>&pi;</mi></msqrt></mfrac><msup><mi>e</mi><mrow><mo>-</mo><mfrac><mrow><msup><mi>x</mi><mn>2</mn></msup><mo>+</mo><msup><mi>y</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>&sigma;</mi><mn>2</mn></msup></mrow></mfrac></mrow></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000031.GIF" wi="831" he="169" /></maths>σ为标准差,x为植被的横坐标,y为植被的纵坐标,然后借助局部最大值法从高斯二维分布函数G(x,y)中提取单木树冠,即在设定高度以上的像元内布设种子并允许其爬向最大坡度的方向,当种子到达的某一像元内高度都高于周边像元时,则将此像元作为树顶,树顶所对应的种子的爬行区域为冠幅,所述的单木容积内的原始波形信息汇总的处理的第二步为LiDAR的信息处理系统借助在单木分割的处理中提取的冠幅和树高参数,在椭圆柱体范围内汇总其中的所有原始波形脉冲能量,构成单木综合脉冲,再借助数字地形模型DTM对单木综合脉冲进行归一化处理;步骤4:进入计算单木综合波形的特征参数并汇总至样地尺度的处理,所述的计算单木综合波形的特征参数并汇总至样地尺度的处理的第一步为LiDAR的信息处理系统通过设置差异阈值来对比数字地形模型DTM提取高度和对应位置上最后一个波峰所在高度,依此判断地形对应的波峰位置从而推测冠层对应的波峰,提取包括地面到脉冲能量质心的距离HOME、波形起点到地面距离所表示的波形距离WD、地面到脉冲能量质心的距离HOME同波形距离WD的比值所表示的高度中位数比、波形起点到第一波峰的距离所表示的冠层外层粗糙度ROUGH、波形起点到第一波峰连线的波形垂直角度所表示的前坡度角FS以及冠层高度和地面到脉冲能量质心的距离HOME的差值除以冠层高的值VDR所构成的特征变量,特征变量是用来作为后续建模的备选因子;所述的计算单木综合波形的特征参数并汇总至样地尺度的处理的第二步为LiDAR的信息处理系统把特征变量汇总为样地尺度上的均值和标准差;步骤5:进行对各生物物理特性的特征变量敏感性分析、反演模型构建和精度验证,LiDAR的信息处理系统首先借助Pearson’s相关系数分析2组LiDAR特征变量对于生物物理特性的敏感性r,具体如公式(4)所示:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><mi>r</mi><mo>=</mo><mfrac><mrow><mrow><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mo>(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub><mo>)</mo></mrow></mrow><mrow><mo>(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>y</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub><mo>)</mo></mrow></mrow><msqrt><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>&CenterDot;</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>y</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000041.GIF" wi="1491" he="289" /></maths>其中,x<sub>i</sub>为地面实测的某生物物理特性,y<sub>i</sub>为某LiDAR特征变量,<img file="FDA0000544430090000042.GIF" wi="63" he="68" />为x<sub>i</sub>的平均值,<img file="FDA0000544430090000043.GIF" wi="61" he="68" />为y<sub>i</sub>的平均值,n为样地数量,然后应用逐步回归法构建预测模型,逐步回归的实施过程是每一步都要对已引入回归方程的变量计算其偏回归平方和,然后选一个偏回归平方和最小的变量,在预先给定的F水平下进行P&lt;0.05的显著性检验,P值为用于检验显著性的参数,如果显著便保留该变量;相反,则要剔除该变量,为了保证模型因变量和自变量之间的线性相关关系,自然对数变换被应用于所有的因变量和自变量,所述的预测模型因变量位地面调查的生物物理特性,所述的自变量为LiDAR提取的特征变量,接着预测模型采用决定系数R<sup>2</sup>、均方根误差RMSE和相对均方根误差rRMSE并结合公式(5)、公式(6)和公式(7)进行评价:<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><msup><mi>R</mi><mn>2</mn></msup><mo>=</mo><mn>1</mn><mo>-</mo><mfrac><mrow><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow><mrow><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000051.GIF" wi="1398" he="284" /></maths>其中,x<sub>i</sub>为地面实测的某生物物理特征,<img file="FDA0000544430090000052.GIF" wi="65" he="68" />为x<sub>i</sub>的平均值,<img file="FDA0000544430090000053.GIF" wi="66" he="76" />为模型估算的某生物物理特征,n为样地数量。<maths num="0006" id="cmaths0006"><math><![CDATA[<mrow><mi>RMSE</mi><mo>=</mo><msqrt><mfrac><mn>1</mn><mi>n</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mi>i</mi></msub><mo>)</mo></mrow><mn>2</mn></msup></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000544430090000054.GIF" wi="1415" he="162" /></maths><maths num="0007" id="cmaths0007"><math><![CDATA[<mrow><mi>rRMSE</mi><mo>=</mo><mfrac><mi>RMSE</mi><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>i</mi></msub></mfrac><mo>&times;</mo><mn>100</mn><mo>%</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow><mo>.</mo></mrow>]]></math><img file="FDA0000544430090000055.GIF" wi="1436" he="142" /></maths>
地址 210037 江苏省南京市龙蟠路159号