一种基于机载激光雷达的林窗时空变化监测与量化方法与流程

本发明属于林业遥感
技术领域：
，具体涉及一种基于机载激光雷达的林窗时空变化监测与量化方法。
背景技术：
：林窗是近年来国内外森林经营管理关注的热点问题之一。“林窗(canopygap)”一词最初由英国学者watt1947年提出，用以表示森林中一株以上冠层乔木死亡后所产生的林中空地或小地段，是新个体入侵、占据和更新的空间。森林循环理论将森林看作空间上异质、时间上变动的“流动镶嵌体”，林窗空间分布(大小、形状以及位置)的动态变化在某种程度上增强了干扰生境的异质性，影响了群落内斑块的性质及其镶嵌状况，从而影响群落结构与动态过程，最终影响景观的结构和功能，因此对林窗空间分布监测与变化状况的量化，是全面系统地认识森林生态系统长期动态变化过程的关键，是科学进行森林经营管理的基础。在林窗发育的不同阶段，空间结构的异质性和功能的异质性是交织变化或发展的。林窗的生成、稳定(维持)、扩展、萎缩、位移或闭合等时空变化增强了林窗功能异质化，势必对植物侵入、种子萌发、幼苗定居、幼树生长起选择作用。受到空间数据采集技术的限制，学者对林窗的研究主要集中在单个观测时期内“静态”地分析林窗空间结构的异质性对植被的物种组成、形态与生理、更新与生长等方面的影响，而“动态”地监测林窗时空变化的研究较少。林窗空间结构的监测主要依靠的还是野外调查，包括直接测量法和相片估测法，该类方法费时费力，人工成本高且受外界因素影响较大；又由于缺乏变化模式的量化指标，导致林窗空间分布变化状态的表征不客观、不系统，影响了林窗动态监测的规范化处理。因此，有必要研究一种方法或技术，实现对林窗分布的时空动态监测，具体解决林窗生成、稳定、扩展、萎缩、位移或闭合等变化状态的科学表征问题，最终实现林窗时空变化监测与量化。机载激光雷达(airbornelaserscanning，als)是一种新兴的主动式遥感技术，能在多时空尺度上获取森林生态系统高精度的植被结构信息、三维地形特征。als对林窗的高精度监测在认识森林干扰状况、更新规律以及反演生态环境参数等方面具有巨大潜力，但我国在这方面尚属起步阶段。因此，加强林窗空间分布动态监测方面的研究，尤其是利用先进的遥感技术对林窗时空动态模式进行量化，有助于提高林窗研究的理论水平，进一步揭示森林循环的动态规律与生物多样性维持的机制，对提高森林经营管理水平具有重要的现实意义。技术实现要素：本发明目的是提供了一种基于机载激光雷达的快速、客观、准确监测林窗空间分布及其变化状态的量化方法，为了弥补现有技术中林窗时空变化监测与量化的技术空白。为实现以上目的，本发明采用如下技术方案：一种基于机载激光雷达的林窗时空变化监测与量化方法，包括以下步骤：步骤1：基于多时相的激光雷达点云数据，识别林窗、绘制林窗边缘矢量图并保存林窗矢量数据集；所述的激光雷达点云数据是指由固定翼飞机或无人机搭载激光雷达扫描仪获取遥感数据，每一个激光点的数据包括x，y，z坐标数据以及回波强度数据；利用现有遥感图像处理软件theenvironmentforvisualizingimages(envi,v5.3)设置点云数据的正确投影坐标系统，计量单位为米。步骤2：比较两个时期的林窗矢量数据集，构建林窗新生成、稳定、扩展、萎缩、位移或闭合状态的逻辑数学模型；步骤3：将两个时期的林窗矢量数据集叠合，进行林窗矢量多边形叠置分析，根据所述步骤2构建的逻辑数学模型，构建林窗矢量多边形变化状态的判别指标。优选的，所述步骤1中，设林窗矢量数据集为gn，设林窗中某一个点(x，y)的取值为g(x，y)，识别林窗的模型如下：式中，g(x，y)＝1，表示点(x，y)对应的位置为林窗；g(x，y)＝0，表示点(x，y)对应的位置为非林窗；x，y分别为栅格(x，y)的横坐标和纵坐标，chm(x，y)是栅格(x，y)上的冠层高度值；a为冠层边缘平均高度的判别阈值。判别阈值由经验值或实地抽样调查确定，如判别阈值为5米。优选的，所述步骤1中，chm为冠层表面模型与地表高程模型的差值，即：chm(x，y)＝dsm(x，y)—dem(x，y)；式中dsm为冠层表面模型，dem为数字高程模型。首先，利用自适应不规则三角网滤波方法将原始点云划分为冠层点云和地表点云。利用克里格插值法将冠层点云内插成dsm，反距离权值内插法将地表点云插值生成dem。考虑到点云的密度和冠层表征的特点，设定dsm和dem具有相同栅格分辨率。其次，利用envi软件对dsm和dem进行差值运算，得到冠层高度栅格模型chm，并在此基础上根据林窗边缘冠层的高度阈值a进行林窗和非林窗的二值图g(x，y)分类。然后，对林窗识别结果进行后处理，以提高识别精度，利用arcgis软件在g(x，y)二值图基础上进行形态学滤波操作，最终确定林窗边缘并对其进行矢量化操作，绘制林窗边缘矢量图并保存数据集gn。优选的，所述步骤2的具体步骤包括：步骤21：设定先后两个时期ti,tj获取的林窗矢量数据集为式中i＜j；i,j＞0；步骤22：比较ti时期的林窗矢量数据集与tj时期的林窗矢量数据集，ti时期的林窗数据集与tj时期对应的林窗数据集进行比较，可能出现三种情况：两期的林窗多边形没有重叠，完全重叠与部分重叠。没有重叠表示ti林窗存在，而tj林窗消失(闭合)了；或ti林窗不存在，而tj林窗新生成了。完全重叠表示两个时期的林窗个体都存在，林窗维持稳定状态没有发生几何变形，或者在原位置上扩展或萎缩。部分重叠的现象是由于林窗边缘植被更新或死亡的空间分布的异质性，导致tj林窗可能发生位移(没变形)，萎缩并位移或扩展并位移等动态。分别定义林窗新生成、稳定、扩展、萎缩、位移或闭合状态的逻辑数学模型为：新生成：闭合：稳定：扩展无位移：萎缩无位移：扩展有位移：萎缩有位移：优选的，所述步骤3中构建林窗时空动态变化的判别指标的具体步骤包括：步骤31：分别对ti,tj两个时期的林窗矢量多边形进行拓扑结构编码，式中下标i或j表示第i或j时期获取的数据，而且定义i＜j；i,j＞0；步骤32：将ti,tj两个时期的林窗矢量多边形进行叠合操作，输出图层综合ti,tj两个时期林窗矢量数据集的属性，保留林窗在ti,tj两个时期内所有多边形的特征；图层叠合过程中将原来的多边形要素分割成新要素，新要素综合了原来两个数据集的属性，其结果通常就是把一个多边形按另一个多边形的空间分布状态进行交集运算，从而划分成多个多边形，同时在属性分配过程中将输入对象的属性拷贝到新对象属性表中，这样做的目的是完全保留林窗空间分布中几何特征(林窗边缘)和属性特征(林窗、非林窗)的变化信息；步骤33：对叠合后的林窗矢量多边形的共享弧段进行编码，然后计算单个林窗多边形属性特征ti,tj叠加面积与原有林窗面积的比值；步骤34：利用所述步骤32和步骤33中的叠合处理、共享弧段以及面积比值构建林窗变化的判别指标，监测在不同时期林窗新生成、稳定、扩展、萎缩、位移或闭合的变化状态以及量化其变化程度。优选的，所述步骤31中分别对ti,tj两个时期的林窗矢量多边形进行拓扑结构编码时，采用e表示林窗范围内或表示已经存在的林窗，o表示林窗范围外或表示非林窗。优选的，所述步骤33中，对叠合后的林窗矢量多边形的共享弧段编码为bb、bi、ib、bn、nb，其中b表示多边形的边界，i表示多边形的内部，n表示非边界或内部。优选的，所述步骤33中，面积比值的具体计算公式为：p(eeij)＝area(eeij)/area(ei)，p(eoij)＝area(eoij)/area(ei)，p(oeij)＝area(oeij)/area(ei)；式中，p表示面积的比值，area表示林窗范围或面积，ei：第ti时期，单个林窗的多边形范围，eeij：假设某个林窗ti时期存在，tj时期也存在，ee表示ti～tj时期内林窗的重叠部分，oeij：假设ti时期某个林窗未出现，在tj时期已出现，oe表示ti～tj期间内林窗的变化范围，eoij，假设ti时期某个林窗存在，在tj时期已消失，eo表示ti～tj期间内林窗的消失范围。优选的，所述步骤31中，利用arcgis软件对不同时期内林窗矢量多边形进行拓扑结构编码。本发明的技术方案具有以下有益效果：本发明利用机载激光雷达技术与地理空间数据拓扑分析方法，实现林窗分布(矢量多边形)快速、准确的监测和绘制，对单个林窗在不同观测期内变化特征(诸如新生成、稳定、扩展、萎缩、位移或闭合等状态)进行客观、系统地量化，以期提高森林空间结构精准监测的水平，促进林窗时空变化状态的科学化表征。本发明加强了林窗空间分布动态监测方面的研究，尤其是利用先进的遥感技术对林窗时空动态模式进行量化，有助于提高林窗研究的标准化水平，进一步揭示森林循环的动态规律与生物多样性维持的机制，提高森林经营管理水平。(1)林窗空间分布状况是林窗研究的基础，本发明利用机载激光雷达技术与地理空间数据拓扑分析方法，实现林窗多边形快速、准确的监测和绘制，对单个林窗在不同观测期内变化特征(诸如新生成、稳定、扩展、萎缩、位移或闭合等状态)进行客观、系统地量化，以期提高森林空间结构精准监测的水平，促进林窗空间变化状况表征的标准化程度；解决了现有技术中存在的景观尺度上林窗空间分布的动态监测难题以及单个林窗具体变化程度的准确表述等问题。(2)针对林窗空间分布的监测，本发明采用的机载激光雷达技术是一种主动式遥感技术，通过激光雷达穿透能力获得准确反映森林冠层空间结构的三维点云数据，在此基础上进行林窗识别与空间分布的绘制，该方法能快速、准确绘制多尺度林窗分布多边形，利用多时相的激光点云数据相比较，能准确反映林窗空间分布的变化状况；解决了现有技术主要集中在实地调查方法中存在的费时费力、过程冗繁、精度不高，受自然环境的约束较大等难题。(3)本发明的方案能将复杂的地理对象和现象简化和抽象到计算机中进行表示、处理和分析，利用地理空间数据拓扑分析方法对林窗空间分布变化状况进行量化，能系统、准确表述林窗新生成、稳定、扩展、萎缩、位移或闭合等变化特征；解决了目前针对林窗空间分布变化状况，缺乏一套完整、科学的量化指标，实际操作过程主观判断较多，人为因素干扰了变化状况的准确表述等技术难题。附图说明图1为冠层激光点云和地表激光点云示意图；图2为多边形叠合分析示意图；图3为待测区域冠层表面模型chm示意图；图4为chm剖面图以5米高度为阈值判别林窗；图5为待测区域林窗识别示意图；图6为“扩展”变化状态的林窗边缘示意图，a)、b)分别为同一林窗t1和t2时期数据；图7为“萎缩”变化状态的林窗边缘示意图，a)、b)分别为同一林窗t1和t2时期数据。具体实施方式以下提供本发明的优选实施例，以助于进一步理解本发明。本领域技术人员应了解到，本发明实施例的说明仅是示例性的，并不是为了限制本发明的方案。实施例1：本方案基于机载激光雷达的林窗时空变化监测与量化方法，具体实施分为以下3个步骤：步骤1：基于多时相的激光雷达点云数据，识别林窗、绘制林窗边缘矢量图并保存林窗矢量数据集；激光雷达点云数据是指由固定翼飞机或无人机搭载激光雷达扫描仪获取遥感数据，每一个激光点的数据包括x，y，z坐标数据以及回波强度数据；利用现有遥感图像处理软件theenvironmentforvisualizingimages(envi,v5.3)设置点云数据的正确投影坐标系统，计量单位为米。其中，设林窗矢量数据集为gn，设林窗中某一个点(x，y)的取值为g(x，y)，识别林窗的模型如下：式中，g(x，y)＝1，表示点(x，y)对应的位置为林窗；g(x，y)＝0，表示点(x，y)对应的位置为非林窗；x，y分别为栅格(x，y)的横坐标和纵坐标，chm(x，y)是栅格(x，y)上的冠层高度值；a为冠层边缘平均高度的判别阈值。判别阈值由经验值或实地抽样调查确定，如判别阈值为5米。chm为冠层表面模型与地表高程模型的差值，即：chm(x，y)＝dsm(x，y)—dem(x，y)；(2)式中dsm为冠层表面模型，dem为数字高程模型。首先，利用自适应不规则三角网滤波方法将原始点云划分为冠层点云和地表点云，参见附图1所示。利用克里格插值法将冠层点云内插成dsm，反距离权值内插法将地表点云插值生成dem。考虑到点云的密度和冠层表征的特点，设定dsm和dem具有相同栅格分辨率。其次，利用envi软件对dsm和dem进行差值运算，得到冠层高度栅格模型chm，并在此基础上根据林窗边缘冠层的高度阈值a进行林窗和非林窗的二值图g(x，y)分类。然后，对林窗识别结果进行后处理，以提高识别精度，利用arcgis软件在g(x，y)二值图基础上进行形态学滤波操作，最终确定林窗边缘并对其进行矢量化操作，绘制林窗边缘矢量图并保存数据集gn。步骤2：比较两个时期的林窗矢量数据集，构建林窗新生成、稳定、扩展、萎缩、位移或闭合状态的逻辑数学模型；具体步骤包括：步骤21：设定先后两个时期ti,tj获取的林窗矢量数据集为式中i＜j；i,j＞0；步骤22：比较ti时期的林窗矢量数据集与tj时期的林窗矢量数据集，ti时期的林窗数据集与tj时期对应的林窗数据集进行比较，可能出现三种情况：两期的林窗多边形没有重叠，完全重叠与部分重叠。没有重叠表示ti林窗存在，而tj林窗消失(闭合)了；或ti林窗不存在，而tj林窗新生成了。完全重叠表示两个时期的林窗个体都存在，林窗维持稳定状态没有发生几何变形，或者在原位置上扩展或萎缩。部分重叠的现象是由于林窗边缘植被更新或死亡的空间分布的异质性，导致tj林窗可能发生位移(没变形)，萎缩并位移或扩展并位移等动态。分别定义林窗新生成、稳定、扩展、萎缩、位移或闭合状态的逻辑数学模型为：新生成：闭合：稳定：扩展无位移：萎缩无位移：扩展有位移：萎缩有位移：步骤3：将两个时期的林窗矢量数据集叠合，进行林窗矢量多边形叠置分析，根据所述步骤2构建的逻辑数学模型，构建林窗矢量多边形变化状态的判别指标。其中，构建林窗时空动态变化的判别指标的具体步骤包括：步骤31：利用arcgis软件分别对ti,tj两个时期的林窗矢量多边形进行拓扑结构编码，式中下标i或j表示第i或j时期获取的数据，而且定义i＜j；i,j＞0，采用e表示林窗范围内或表示已经存在的林窗，o表示林窗范围外或表示非林窗。步骤32：将ti,tj两个时期的林窗矢量多边形进行叠合操作，输出图层综合ti,tj两个时期林窗矢量数据集的属性，保留林窗在ti,tj两个时期内所有多边形的特征。参见附图2所示，图层叠合过程中将原来的多边形要素分割成新要素，新要素综合了原来两个数据集的属性，其结果通常就是把一个多边形按另一个多边形的空间分布状态进行交集运算，从而划分成多个多边形，同时在属性分配过程中将输入对象的属性拷贝到新对象属性表中，这样做的目的是完全保留林窗空间分布中几何特征(林窗边缘)和属性特征(林窗、非林窗)的变化信息；步骤33：对叠合后的林窗矢量多边形的共享弧段进行编码，对叠合后的林窗矢量多边形的共享弧段编码为bb、bi、ib、bn、nb，其中b表示多边形的边界，i表示多边形的内部，n表示非边界或内部；然后计算单个林窗多边形属性特征ti,tj叠加面积与原有林窗面积的比值，面积比值的具体计算公式为：p(eeij)＝area(eeij)/area(ei)，p(eoij)＝area(eoij)/area(ei)，p(oeij)＝area(oeij)/area(ei)；式中，p表示面积的比值，area表示林窗范围或面积，ei：第ti时期，单个林窗的多边形范围，eeij：假设某个林窗ti时期存在，tj时期也存在，ee表示ti～tj期间内林窗的重叠部分，oeij：假设ti时期某个林窗未出现，在tj时期已出现，oe表示ti～tj期间内林窗的变化范围，eoij，假设ti时期某个林窗存在，在tj时期已消失，eo表示ti～tj期间内林窗的消失范围。步骤34：利用所述步骤32和步骤33中的叠合处理、共享弧段以及面积比值构建林窗变化的判别指标，参见下表1所示，监测林窗的新生成、稳定、扩展、萎缩、位移或闭合等变化状态及量化其变化程度。表1林窗动态变化拓扑分析表以下通过具体实例说明本发明的应用：研究区概况：研究地位于湖南省东北部的福寿山林场(28°3′00″—28°32′30″n，113°41′15″—113°45′00″e)，地处罗霄山脉连云山支脉，地势南高北低，平均海拔1200余米，平均坡度为22—27度，呈现群山重叠的地貌。年平均气温12.1℃，年日照1500小时，无霜期217天。主要植被类型是典型的中亚热带常绿阔叶林，上层乔木树种主要有杉木、黄山松、青冈栎、苦槠、檫木、江南桤木、山核桃及壳斗科植物。实地林窗调查采用半球面影像法，根据鱼眼镜头投影原理确定林窗大小；采用角规法或伸缩式测高器测量林窗边缘木高度；差分gps或全站仪测量林窗中心海拔和边缘位置。野外调查方法：在固定样地内设置80个30m×30m的样方，用差分gps确定每个样方的位置，在每个样方中心和对角线四分位处用数码相机外接鱼眼镜头(广角为183°，正投影)获取半球面林冠影像，影像方向与磁北方向重合。用gaplightanalyzer(gla，v2.0，图像处理软件)对冠层照片进行分析并绘制出林窗边缘。调查工作分别在2016和2009年夏季，选择天气晴朗状况下进行。80个样方数据中40个作为训练样本数据，其余40个作为验证样本数据。机载激光雷达数据：t1,t2激光雷达数据获取时间分别是2009年6月和2016年9月，t1机载激光扫描系统为lms-q560，激光束扫描角度平均22.5°，平均光斑大小50cm，点云密度为2～6个/m2。t2机载激光扫描系统为altm2050，激光束扫描角度平均15°，平均光斑大小25cm，点云密度为2～10个/m2。每束激光都包含第一回波和最后回波的坐标值、高度值、强度值等信息。激光雷达点云数据都采用las格式，利用envi软件设定激光点云数据的地理坐标系统为cgcs2000，参考椭球为wgs84，投影代号为38，激光雷达数据以点云存储格式记录(参见下表2所示)。表2激光点云存储格式示意表(每个激光点对应一行记录)航带号x坐标/my坐标/m海拔高度/m脉冲回波强度138476525.6413152896.212864.682.4138476816.6833153028.5041023.530.8138476790.2253152909.4411012.421.6138477068.0383152737.4611018.281.2138477147.4133152790.3781045.592.9…………………………基于本发明方法对研究区林窗进行时空变化监测与量化过程：(1)基于激光雷达点云数据的林窗监测利用自适应不规则三角网滤波方法根据高程值将t1机载激光雷达点云数据划分为冠层点云和地表点云，并分别进行栅格插值生成dsm和dem，设定dsm和dem的栅格分辨率同为2米。利用envi软件对dsm和dem进行栅格差值运算，生成2米分辨率的冠层高度模型chm，参见附图3所示。利用林窗识别模型(公式1)在chm上进行林窗判别，设定林窗边缘冠层高度阈值a为5米，参见附图4所示，在chm剖面图上，植被高度小于5米的为林窗，反之为林冠(非林窗)。利用envi软件生成g(x，y)林窗和非林窗范围的二值图并进行形态学滤波操作，采用“腐蚀”和“膨胀”算法消除林木之间的小间隙(2米以内)或林冠中的小“空洞”，或栅格图像噪音，确定最终的林窗范围，参见附图5所示，在此基础上进行“栅格转化矢量”操作，绘制并保存林窗边缘多边形矢量数据g2009。t2机载激光雷达数据操作流程如同上述t1操作，绘制并保存林窗的边缘多边形矢量数据g2016。(2)基于空间分析的林窗动态模式的量化利用arcgis软件中的arccatalog模块创建地理数据库，新建多边形要素集，选择cgcs2000大地坐标系统，分别将g2009和g2016林窗多边形矢量数据导入到地理数据库中，创建矢量数据的拓扑结构，设置xy容差值为0.1m。对g2009和g2016多边形要素进行叠合分析，具体操作为“联合”，根据林窗变化的逻辑数学模型，也即公式(3)～(9)，以及表1内容示意添加拓扑规则，输出要素类将包含代表所有输入的几何并集的多边形以及所有输入要素类的所有字段，然后对叠合结果进行拓扑关系验证，根据空间分析的结果，依次选择符合林窗变化条件的结果进行统计并绘制结果。(3)结果分析专业数据统计分析软件statisticalproductandservicesolutions(spss，v19)对t1和t2激光雷达数据监测的林窗边缘与野外调查结果(验证样本)进行线性拟合精度分析，利用配对t检验(pairedt-test)比较本文方法与野外调查的差异，参见下表3所示。表3本方案的方法与野外调查结果的比较表x为野外调查值，y为本方法的方法估测值。林窗边缘位置差值正态分布显著性sig.≥0.05，说明利用配对t检验法具有统计学意义。配对t检验的原假设是位置差值的分布符合平均值为0的t分布，t1和t2数据的双侧显著性p都大于0.05，说明激光雷达监测的林窗位置与野外调查结果不存在显著差异。从相对均方根误差来看，林窗位置的监测结果误差较小，且与野外调查的结果相关性较高。本方法测得的林窗位置及分布状况具有较高的精度，可用于快速测量不同大小、不同形状的林窗。t1机载激光雷达数据监测到的林窗密度为11.67个/hm-2，林窗面积均值为81.02m2，最小4.06m2，最大727.43m2，林窗面积集中分布在四分位距间内，特征分布的偏度为2.542，说明其具有正偏离性，即多为较小面积的林窗，参见下表4所示。t2林窗也有相似情况，其密度为12.81个/hm-2，也以小面积林窗为主。表4林窗面积(m2)基本特征统计表根据变化状态的逻辑数学模型，对t1-t2林窗变化进行量化，参见表5所示，新生成林窗以5～50m2占绝大多数，干扰方式主要是折枝；由于林窗的侧向更新和垂直更新，小级别的林窗更容易完全闭合；在“稳定”状态下，大于300m2的大林窗比例最大，“扩展无位移”状态下，同样是大林窗为主，其余级别的林窗比例相差不大，其主要原因可能是大林窗内生境恶化，部分边缘木产生新的干扰面积，使得林窗面积变大，参见附图6所示；“萎缩无位移”状态150m2以下林窗为主，原因可能是林窗侧向更新占主导作用，参见附图7所示；“扩展有位移”和“萎缩有位移”两种变化状态，各级林窗都有所比例，主要原因是林窗边缘木的更新和干扰同时作用，最终导致了林窗边缘的空间分布变化。上述变化过程不仅验证了林窗干扰的生态过程，同时完成了各种变化状态的量化工作，为森林生态学研究提供了准确的监测数据，为进一步揭示森林循环的动态规律与生物多样性维持的机制提供信息化支撑和技术保障。表5林窗动态变化比例统计表最后应当说明的是，以上实施例仅用于说明本申请的技术方案而非对其保护范围的限制，尽管参照上述实施例对本申请进行了详细的说明，所述领域的普通技术人员应当理解：本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或等同替换，但以上变更、修改或等同替换，均在本申请的待授权或待批准之权利要求保护范围之内。当前第1页12