一种基于激光雷达的林窗太阳辐射监测方法与流程

本发明涉及一种基于激光雷达的林窗太阳辐射监测方法。
背景技术：
：1947年林窗(canopygap)概念由英国生态学家Watt首次提出，它主要指森林群落中主林层受人为干扰(择伐)或自然干扰(大风、雪、水灾、泥石流、虫害等灾害)在林冠层形成的不连续的林中空隙。林窗是森林群落中经常发生的重要中小尺度干扰,是森林演替的一个重要阶段，它不仅是森林群落演替的驱动要素,也在森林结构、物种组成、动态和演替中起着重要作用，已成为当前森林生态学研究最活跃的领域之一。林地环境中，太阳辐射分为直射和散射辐射，直接光是太阳直射到冠层上方的光，散射光是光束在大气层中以不同方向散射到冠层上的光。太阳辐射是林窗最重要的生态环境因子之一，光照直接影响着动植物的生理生态过程，同时，光照也会对林窗地表温度、空气湿度、土壤水分等产生影响，从而导致林窗环境异质性。环境异质性决定了林窗更新，在森林动态及演替过程中起着重要作用。所以，林窗太阳辐射的测量是生态学家十分关注的问题，由于受到林窗大小、林窗形状、林冠高度以及林地坡度、坡向的影响，林窗光照表现出复杂的空间异质性，加之太阳位置时刻变化，使得太阳辐射的测量变得十分困难。传统林窗太阳辐射的监测分为直接监测法和间接相片法。直接监测法(如照度计)只能监测有限点的光照信息，费时费力，人力成本高。间接相片法(如鱼眼镜头影像)虽然较为简单快速，但在监测林窗光环境异质性时，仍然需要在林窗不同位置拍摄多个的半球面影像。上述两类方法在不同程度上存在工作量大、过程冗繁、破坏性强、耗时耗力等问题。机载激光雷达(AirborneLaserScanning，ALS)是一种新兴的主动式遥感技术，能在多时空尺度上获取森林生态系统高精度的植被结构信息、三维地形特征。ALS对林窗的高精度监测在认识森林干扰状况、更新规律以及反演生态环境参数等方面具有巨大潜力。国内研究多集中在利用ALS进行森林生物量监测，而对林窗太阳辐射的测量极少涉及。因此，找到一种快速、客观、准确测量林窗光环境的方法对林窗研究非常重要，有助于解释林窗物种组成、生物多样性与光照强度空间异质性的关系，有助于揭示控制林窗中植被演替和物种更新过程的机制，对认识开发我国森林事业具有重要的作用。因此，有必要设计一种基于激光雷达的林窗太阳辐射监测方法。技术实现要素：本发明所要解决的技术问题是提供一种基于激光雷达的林窗太阳辐射监测方法，该方法是基于机载激光雷达的快速、客观、准确的监测林窗光照强度的方法，能弥补现有技术中林窗太阳辐射监测方法的技术空白。发明的技术解决方案如下：一种基于激光雷达的林窗太阳辐射监测方法，包括以下步骤：步骤1：基于激光雷达ALS点云数据识别林窗，以确定林窗边缘；步骤2：计算复杂地形条件下冠层上方接收到的局部入射角方向太阳辐射总量Gi；步骤3：基于Gi计算得出林窗内不同位置的太阳辐射强度Globalgap。步骤1中，识别林窗的模型如下：式中，G(x，y)＝1，表示点(x，y)对应的位置为林窗；x，y分别为栅格(x，y)的横坐标和纵坐标，CHM(x，y)是栅格(x，y)上的冠层高度值；a为冠层边缘平均高度的判别阈值。判别阈值由经验值或实地抽样调查确定，如判别阈值为5米。步骤2中：Gi＝Bi+Di；式中，Gi为局部入射角方向太阳辐射总量，Bi,Di分别为局部入射角方向太阳直射辐射量和散射辐射量。局部入射角方向太阳直射辐射量的计算公式如下：Bi＝In·δi·hillshade；式中δi为光线局部入射角，hillshade为山体阴影影响系数；具体地，有：光线局部入射角：山体阴影影响系数：hillshade＝(cosα·cosSlop)+[sinα·sinSlop·cos(αs-Asp)]；式中Slop为地形的坡度，Asp为地形的坡向(坡度是指坡面的倾斜程度，坡向是指地形坡面的朝向)；δs为太阳赤纬角，有δs＝23.45sin(360°·(284+N)/365)，式中N为积日(一年中的第几天)；L为地理纬度；hs为太阳时角。局部入射角方向太阳散射辐射量Di的计算方法如下：考虑到周围山体可能对散射光的影响，计算山区地形中的Di时，分为入射角方向上无山体遮挡δi＞0和有山体遮挡δi≤0两种情况：式中变量说明：(1)skyview为可视域范围内天穹比例；该参数由ENVI软件计算得出；(2)hillshade为山体阴影系数；hillshade＝(cosα·cosSlope)+[sinα·sinSlope·cos(αs-Asp)]；(3)Io为大气层外的太阳辐射量，Io＝S0·(1+0.0344cos(360°·N/365))，式中，S0为太阳辐射常数1367(w/m2)，N为积日(一年中的第几天)。(4)τb为晴朗天气太阳直射辐射的穿透率，τb＝0.56(e-0.65M+e-0.095M)，式中，M为大气质量系数，具体有，M＝[1229+(614sinα)2]1/2-614sinα式中，α为太阳高度角；(5)So为太阳辐射常数1367(w/m2)(6)δi为光线局部入射角，有：δi＝sinδs·(sinL·cosSlope-cosL·sinSlope·cosAsp)+cosδscoshs·(cosL·cosSlope+sinL·sinSlope·cosAsp)+cosδs·sinSlope·sinAsp·sinhs式中Slope为地形的坡度；Asp为地形的坡向；δs为太阳赤纬角，L为待测区域的地理纬度，hs为太阳时角；(7)αs为太阳方位角；(8)Slope为地形坡度，计算公式详见公式11；(9)Dio为初始局部入射角方向太阳散射辐射强度：Dio＝Io·τd·(cosSlop)2/2sinα。林窗太阳辐射Globalgap的计算公式如下：Globalgap＝crownshade(1-LPI)·Gi；式中，crownshade为林窗边缘冠层产生的阴影，LPI为光束穿透系数，Gi为冠层上方太阳辐射总量；crownshade＝(cosα·cosSlopeDSM)+[sinα·sinSlopeDSM·cos(αs-AspDSM)]式中，SlopeDSM为冠层表面模型DSM的坡度，反映冠层表面的倾斜程度，AspDSM为冠层表面模型DSM的坡向，反映冠层表面的朝向。DSM的坡度和坡向计算参考公式11；α、αs为太阳高度角和方位角。光束透射系数number(地表点云)和number(地表点云+冠层点云)分别表示单位面积(如2m×2m)中地表点云数量与点云总量；LPI为光束透射系数，该系数与crownshade相结合计算林冠下方的太阳辐射，因为crownshade借鉴的是山体阴影计算公式，其假设条件是山体遮挡不透光，但实际是林冠层由于疏密度不同是有孔隙的，具有一定透光性，所以利用LPI系数模拟冠层的透光性。另一方面，理论上LPI考虑的是单位面积上垂直方向的点云比值，但由于实际操作中机载扫描仪发射的激光束与冠层顶部之间往往存在一定的夹角，植被点云的分布可能存在一定的变形，所以，本发明提出利用邻域分析法(neighborhoodanalysis)来优化LPI的计算，邻域范围由林窗边缘木高度和激光束扫描角度确定。目标对象(图5中黑色单元)的LPI为：式中，N为邻域内栅格单元的个数。目标对象的LPI即为邻域范围内每个栅格单元LPI值的均值。有益效果：本发明的基于激光雷达的林窗太阳辐射监测方法，其关键点如下：(1)在冠层上方太阳辐射总量的计算基础上，根据冠层表面模型DSM的起伏状态，顾及冠层遮蔽因素计算得出林窗内不同位置的太阳辐射强度；(2)构建基于激光点云的光束穿透系数模型(LPI)，提高太阳辐射监测精度。(1)通过激光点云生成能准确反映冠层起伏状况的表面模型DSM，借鉴山体阴影(已有技术)的计算方法，利用太阳高度角和方位角，计算冠层间的遮蔽关系，再根据冠层上方太阳辐射总量，可准确计算出林窗内不同位置的太阳辐射强度。采用数字高程模型DEM可计算出传统意义上山体遮蔽关系，在山体阴影的基础上，再利用DSM估算出小尺度范围内冠层间的遮蔽关系，该方法适合计算复杂地形环境中不同大小林窗、不同位置的太阳辐射。(2)冠层有孔隙，能让部分光线透射到地面，所以光线透射因素在估算由冠层所产生的阴影区光照时显得非常关键。激光束穿透冠层的过程中能产生高密度、高精度的点云，可模拟晴朗天气自然光线照射冠层情形。本发明提出通过计算单位面积上地表点云数量与点云总量(冠层点云与地表点云之和)的比值来求光束穿透系数LPI，并利用邻域分析法来优化LPI。该方法操作简单、计算准确，不仅可有效避免由点云分布畸变引起的系统误差，而且比传统算法(如光线追踪法)的计算复杂度要小得多。本发明的基于机载激光雷达的林窗太阳辐射监测方法，填补了现有技术中林窗太阳辐射监测方法的技术空白，是一种准确、可重复操作的测量方法，能快速测量晴朗天气情况下林窗内任意位置的太阳辐射强度。针对林窗太阳辐射强度的测定，现有技术如直接测量法和相片法，都只能在单个样方内进行实地操作，当获取大范围样方数据时，传统方法不同程度上存在工作量大、过程冗繁、耗时耗力等问题，且很难在短时间内重复操作，从而导致结果包含时空异质性误差。本发明利用机载激光雷达遥感技术，借助高密度的地表和植被冠层激光点云，计算林窗范围内不同位置的太阳辐射，突破实地测量的技术瓶颈，具有快速、准确、可重复操作等特点，且检测范围广，具有良好的应用前景。附图说明图1为冠层激光点云和地表激光点云示意图；图2为太阳高度角和方位角示意图；图3为2m×2m单位面积上LPI示意图；图4为激光束与林窗边缘木冠层遮蔽关系示意图；图5为邻域分析计算LPI示意图；图6为待测区域冠层表面模型DSM示意图；图7为待测区域数字高程模型DEM示意图；图8为待测区域林窗识别示意图；图9为2016年6月7日11时待测区域林窗太阳辐射计算结果示意图(太阳高度角70.23°，方位角102.092°)；图10为2016年6月7日16时待测区域林窗太阳辐射计算结果示意图(太阳高度角41.40°，方位角276.24°)；图11为本方法与野外调查结果的线性拟合图；具体实施方式以下将结合附图和具体实施例对本发明做进一步详细说明：实施例1：如图1～11，本发明利用机载激光雷达(ALS)点云数据监测晴朗天气状况下林窗内不同位置的太阳辐射，同时顾及复杂地形、冠层形态和透光率的影响，提高监测精度。第一步，利用ALS点云数据判别林窗，确定林窗边缘；；第二步，根据太阳高度角和方位角估测待测区域晴朗天气状况下太阳辐射强度，结合数字高程模型DEM，计算复杂地形条件下冠层上方接受到的太阳直射和散射辐射量；第三步，根据冠层表面模型DSM的起伏状态，顾及冠层遮蔽关系以及光束穿透系数计算得出林窗内不同位置的太阳辐射强度。1、基于ALS点云数据的林窗判别所述ALS数据是指由固定翼飞机或无人机搭载激光扫描仪获取的遥感数据，以点云格式存储。每个点云包括X，Y，Z坐标值以及回波强度信息。利用现有遥感数据处理软件TheEnvironmentforVisualizingImages(ENVI,v5.3)设定点云数据正确的投影坐标系统，计量单位为米。林窗识别模型：式中，G(x，y)＝1，表示点(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)分类。然后，对林窗识别结果进行后处理，以提高识别精度，利用ENVI软件在G(x，y)二值图基础上进行形态学滤波操作，确定最终的林窗范围。最后对林窗范围进行矢量化操作，得到林窗边缘的多边形矢量图。2、晴朗天气山区太阳辐射强度估算理论上来讲，达到地表的太阳辐射总量由局部入射角方向的太阳直射、散射以及地表反射构成。实际上林窗大都分布在森林覆盖率较高的林地中，地表对太阳辐射的反射很小，为了提高计算效率，在计算林窗冠层上方太阳辐射总量时，只考虑局部入射角方向的太阳直射和散射两部分，即：Gi＝Bi+Di(3)式中，Gi为局部入射角方向太阳辐射总量，Bi,Di分别为太阳直射和散射辐射量。(1)局部入射角方向太阳直射辐射计算：Bi＝Io·τb·δi·hillshade(4)式中，I0为大气层外的太阳辐射量，τb为晴朗天气中太阳直射辐射在大气层的穿透率，δi为光线局部入射角，hillshade为山体阴影系数。(上述参数由公式5-12计算得出)Io＝S0·(1+0.0344cos(360°·N/365))(5)式中，S0为太阳辐射常数1367(w/m2)，N为积日(一年中的第几天)。根据Kreith等研究人员在国际权威期刊发表的论文结果，晴朗天气太阳直射辐射的穿透率可简化为：τb＝0.56(e-0.65M+e-0.095M)(6)式中，M为大气质量系数，具体有，M＝[1229+(614sinα)2]1/2-614sinα(7)式中，α为太阳高度角。太阳高度角和方位角(图2)计算：α＝arcsin(sinL·sinδs+cosL·cosδscoshs)(8)式中δs为太阳赤纬角，L为待测区域的地理纬度，hs为太阳时角。光线局部入射角：δi＝sinδs·(sinL·cosSlope-cosL·sinSlope·cosAsp)+cosδscoshs·(cosL·cosSlope+sinL·sinSlope·cosAsp)+cosδs·sinSlope·sinAsp·sinhs(10)式中Slope为地形的坡度，Asp为地形的坡向，具体计算公式：Asp＝Slopesn/Slopewe(11)式中，slopewe为数字高程模型DEM上X轴方向的梯度，slopesn为DEM上Y轴方向的梯度；梯度取值范围0～90°。e5e2e6e1ee3e8e4e7针对上表的数据，有：e以及e1-e8表示在3×3大小的栅格中9个不同位置的高程值，cellsize为栅格大小，即栅格分辨率，如2米。山体阴影系数：hillshade＝(cosα·cosSlope)+[sinα·sinSlope·cos(αs-Asp)](12)(2)局部入射角方向太阳散射辐射计算：根据Gates等研究人员在国际权威期刊发表的论文结果，局部入射角方向太阳散射辐射与散射透射率τd，地形坡度Slope，太阳高度角α有关：Dio＝Io·τd·(cosSlop)2/2sinα(13)式中，Io为大气层外的太阳辐射量(参见公式5),τd＝0.271-0.294·τb(14)考虑到周围山体可能对散射光的影响，计算山区地形中散射辐射量时，分为入射角方向上无山体遮挡δi＞0和有山体遮挡δi≤0两种情况：式中skyview为可视域范围内天穹比例，可由ENVI软件计算得出。3、林窗太阳辐射计算林窗太阳辐射的计算应考虑冠层形态、冠层上方太阳辐射总量以及冠层的透光率等因素的影响，即：Globalgap＝crownshade(1-LPI)·Gi(16)式中，crownshade为林窗边缘冠层产生的阴影，LPI为光束穿透系数，Gi为冠层上方太阳辐射总量。crownshade＝(cosα·cosSlopeDSM)+[sinα·sinSlopeDSM·cos(αs-AspDSM)](17)式中，SlopeDSM为冠层表面模型DSM的坡度，反映冠层表面的倾斜程度，AspDSM为冠层表面模型DSM的坡向，反映冠层表面的朝向，具体计算参考公式11。α、αs为太阳高度角和方位角。crownshade借鉴的是山体阴影计算公式，其假设条件是遮挡部分不透光，但林冠层由于疏密度的不同，具有一定透光性，本发明利用LPI系数模拟冠层的透光性。LPI光束透射系数指太阳光线(包括直射和散射)穿透冠层达到地表的能力，在机载激光雷达应用中，LPI通过计算单位面积(2m×2m)中地表点云数量与点云总量的比值来得到LPI(图3)，即：LPI取值范围是(0,1),0表示光束完全被冠层遮挡，效果如同光束被不透明的物体挡住；1表示光束没有受冠层或植被构件的影响，直接能照射到地面。理论上，LPI计算考虑的是单位面积上垂直方向的点云比值，但由于实际操作中机载扫描仪发射的激光束与冠层顶部之间往往存在一定的夹角，植被点云的分布可能存在一定的变形。所以我们提出利用邻域分析法(neighborhoodanalysis)来优化LPI的计算。经统计，待测区域内林窗边缘木的平均高度为20米，激光束的扫描角度平均为15°，冠层遮蔽的水平距离为5.36米(图4)，由于我们设定的栅格分辨率为2米，所以单边遮蔽距离相当于3个栅格单元长度。由此，构建7×7栅格方形邻域(目标对象1个单元，其旁边各有3个单元)以计算目标对象的LPI，如图5所示，目标对象(图5中黑色单元)的LPI为：式中，N为邻域内栅格单元的个数。目标对象的LPI即为邻域范围内每个栅格单元LPI值的均值。以下通过具体实例说明本发明的应用：研究区概况：研究地位于湖南省东北部的福寿山林场(28°3′00″—28°32′30″N，113°41′15″—113°45′00″E)，地处罗霄山脉连云山支脉，地势南高北低，平均海拔1200余米，平均坡度为22—27度，呈现群山重叠的地貌。年平均气温12.1℃，年日照1500小时，无霜期217天。主要植被类型是典型的中亚热带常绿阔叶林，上层乔木树种主要有杉木、黄山松、青冈栎、苦槠、檫木、江南桤木、山核桃及壳斗科植物。实地林窗的抽样调查采用半球面影像法，根据鱼眼镜头投影原理确定林窗面积；采用角规法或伸缩式测高器测量林窗边界木高度；差分GPS或全站仪测量林窗中心海拔和位置。区域内林窗密度为12.77个/hm-2，占森林面积11.27％；林窗面积、边界木高均值分别为85.06m2和20.33m，区域内多为较小面积、边缘效应不太显著的林窗。野外调查方法：在固定样地内设置80个30m×30m的样方，用差分GPS确定每个样方的位置，在每个样方中心和对角线四分位处用数码相机外接鱼眼镜头(广角为183°，正投影)获取半球面林冠影像，影像方向与磁北方向重合。用GapLightAnalyzer(GLA，V2.0，图像处理软件)对冠层照片进行分析并计算得出林窗太阳辐射总量，计量单位是w/m2。调查工作选择当年6月天气晴朗状况下进行，调查时间段为北京时间6：00至18：00，每隔一小时测量一次，然后求和计算出总的太阳辐射强度。80个样方数据中40个作为训练样本数据，其余40个作为验证样本数据。机载激光雷达数据：机载激光雷达数据ALS获取时间为2016年夏季，选择晴朗少云天气进行飞行，根据试验地区范围和地形高度、梯度等情况，沿着地势设计多条航线。机载激光扫描系统为ALTM2050，激光束扫描角度平均15°，测量精度±20mm，平均光斑大小15cm，点云密度为2～10个/m2。利用ENVI软件设定激光点云数据的投影坐标为墨卡托投影(UTM)方式，参考椭球为WGS84，投影代号38，激光雷达数据以点云存储格式记录(表1)。表1激光点云存储格式，每个激光点对应一行记录基于本发明方法对研究区林窗太阳辐射监测过程：(1)基于ALS点云数据的林窗判别利用ENVI软件设定激光点云数据的投影坐标为墨卡托投影(UTM)方式，参考椭球为WGS84，由待测区域的经纬度确定投影代号为北38。利用自适应不规则三角网滤波技术根据高程值将原始点云划分为冠层点云和地表点云，采用栅格插值法分别对冠层点云和地表点云插值生成DSM(图6)和DEM(图7)栅格图像，并设定DSM和DEM的栅格分辨率同为2米。利用ENVI软件对DSM和DEM进行栅格差值运算，生成2米分辨率的冠层高度模型CHM。利用林窗识别模型(公式1)在CHM上进行林窗范围的判别，根据经验和实地调查结果，设定林窗边缘冠层高度阈值a为5米，生成林窗和非林窗范围的二值图。利用ENVI软件在G(x，y)二值图的基础上进行形态学滤波操作，采用“腐蚀”算法消除林木之间的小间隙(2米以内)或林冠中的小“空洞”，采用“膨胀”算法还原得到真实林冠层空隙，确定最终的林窗范围(图8)，在此基础上进行“栅格转化矢量”操作，勾勒并保存林窗边缘的多边形矢量图。(2)晴朗天气山区太阳辐射强度计算按照公式(3)计算复杂地形环境中林窗冠层上方太阳辐射总量，计量单位为w/m2。利用ENVI软件在DEM基础上计算地形坡度和坡向，结合待测区域的地理纬度、太阳位置、山体阴影等因素，按照公式(4)～公式(15)计算局部入射角方向的太阳直射辐射和散射辐射，生成整个待测区域冠层上方太阳辐射时空分布图。(3)林窗太阳辐射计算根据步骤(2)计算的冠层上方太阳辐射总量为基础，充分考虑冠层遮蔽影响和冠层光线透射率，按照公式(16)计算林窗范围内太阳辐射，生成整个待测区域内林窗太阳辐射强度的时空分布图(图9，10)。以冠层表面模型DSM为基础，按照公式(17)计算由冠层形态造成的阴影。按照公式(18)计算单位面积(2m×2m)内光束透射系数LPI，生成整个待测区域的LPI栅格图像。为了减少由于点云分布引起的系统误差，按照公式(19)及邻域分析法优化LPI的取值。(4)结果验证专业数据统计分析软件StatisticalProductandServiceSolutions(SPSS，V19)对2016年6月7日6：00～18：00时间段林窗内太阳辐射总量的估测与野外调查结果(验证样本)进行线性拟合精度分析(图11)。分别按照林窗大小和林窗内不同位置进行验证，利用配对T检验(pairedt-test)比较4种情况下激光雷达估测与野外调查结果的差异(表2)。根据野外调查样方位置，分别设定面积>200m2的为大林窗，面积≤200m2为小林窗，林窗内位置分为林窗中心区和边缘区。表2本方法与野外调查结果比较x为野外调查值，y为本方法估测值。林窗太阳辐射强度差值正态分布显著性Sig.≥0.05，说明利用配对T检验法具有统计学意义。配对T检验的原假设是光照强度差值的分布符合平均值为0的t分布，小林窗、大林窗和林窗中心的双侧显著性P都大于0.05(表2)，说明这3种情况下激光雷达监测的太阳辐射强度与野外调查结果不存在显著差异，而林窗边缘的光照强度监测与野外调查结果有明显差异。从相对均方根误差来看,大林窗和林窗中心的太阳辐射强度的监测结果优于小林窗和林窗边缘。从回归方程和相关系数来看，激光雷达监测的结果与野外调查的结果呈正相关性，但存在一定程度的低估。综合以上分析，本方法测得的林窗太阳辐射具有较高的精度，可用于快速测量不同大小的林窗、林窗内不同位置的光环境。通过对40个验证样本的数据分析表明，林窗边缘木树种及平均高度和激光点云密度对光束透射率LPI都有影响，进而可能影响到最终的监测结果。所以，实际操作中应增加飞行航带的重叠率以提高点云密度，同时注意多选择激光扫描角度小的数据以提高点云LPI的精度。当前第1页1 2 3