CN107479065A - 一种基于激光雷达的林窗立体结构量测方法 - Google Patents
一种基于激光雷达的林窗立体结构量测方法 Download PDFInfo
- Publication number
- CN107479065A CN107479065A CN201710574787.6A CN201710574787A CN107479065A CN 107479065 A CN107479065 A CN 107479065A CN 201710574787 A CN201710574787 A CN 201710574787A CN 107479065 A CN107479065 A CN 107479065A
- Authority
- CN
- China
- Prior art keywords
- data
- canopy
- laser radar
- woods window
- point cloud
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 93
- 238000005259 measurement Methods 0.000 title claims abstract description 48
- 238000012545 processing Methods 0.000 claims abstract description 13
- 230000008569 process Effects 0.000 claims abstract description 6
- 238000002203 pretreatment Methods 0.000 claims abstract description 4
- 238000001914 filtration Methods 0.000 claims description 7
- 241000630841 Anthracoceros coronatus Species 0.000 claims description 5
- 230000003044 adaptive effect Effects 0.000 claims description 5
- GSJBKPNSLRKRNR-UHFFFAOYSA-N $l^{2}-stannanylidenetin Chemical compound [Sn].[Sn] GSJBKPNSLRKRNR-UHFFFAOYSA-N 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 2
- 238000011835 investigation Methods 0.000 abstract description 3
- 238000010276 construction Methods 0.000 abstract description 2
- 239000013256 coordination polymer Substances 0.000 description 11
- 238000009826 distribution Methods 0.000 description 10
- 239000002023 wood Substances 0.000 description 9
- 238000010586 diagram Methods 0.000 description 8
- 238000004422 calculation algorithm Methods 0.000 description 7
- 238000012544 monitoring process Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 5
- 240000007594 Oryza sativa Species 0.000 description 4
- 235000007164 Oryza sativa Nutrition 0.000 description 4
- 235000011609 Pinus massoniana Nutrition 0.000 description 4
- 241000018650 Pinus massoniana Species 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 235000009566 rice Nutrition 0.000 description 4
- 241000251468 Actinopterygii Species 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000009499 grossing Methods 0.000 description 3
- 239000000047 product Substances 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 238000003860 storage Methods 0.000 description 3
- 238000010998 test method Methods 0.000 description 3
- 240000001548 Camellia japonica Species 0.000 description 2
- 241001365762 Castanopsis eyrei Species 0.000 description 2
- 244000050510 Cunninghamia lanceolata Species 0.000 description 2
- 240000000604 Lespedeza bicolor Species 0.000 description 2
- 235000016677 Lespedeza bicolor Nutrition 0.000 description 2
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 230000009897 systematic effect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 239000011800 void material Substances 0.000 description 2
- 235000006467 Camellia japonica Nutrition 0.000 description 1
- 241000282693 Cercopithecidae Species 0.000 description 1
- 241000633850 Cinnamomum bodinieri Species 0.000 description 1
- 241000723346 Cinnamomum camphora Species 0.000 description 1
- 241000930581 Dicranopteris pedata Species 0.000 description 1
- 241000610361 Eurya Species 0.000 description 1
- 241001189733 Eurya hebeclados Species 0.000 description 1
- 238000006424 Flood reaction Methods 0.000 description 1
- 241000238631 Hexapoda Species 0.000 description 1
- 240000005893 Pteridium aquilinum Species 0.000 description 1
- 235000009936 Pteridium aquilinum Nutrition 0.000 description 1
- 241000219492 Quercus Species 0.000 description 1
- 241001480058 Quercus glauca Species 0.000 description 1
- 241000208422 Rhododendron Species 0.000 description 1
- 244000083724 Rhododendron simsii Species 0.000 description 1
- 241001335286 Viburnum fordiae Species 0.000 description 1
- 241000607479 Yersinia pestis Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002146 bilateral effect Effects 0.000 description 1
- 239000006227 byproduct Substances 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000005253 cladding Methods 0.000 description 1
- 235000018597 common camellia Nutrition 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 238000013211 curve analysis Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 239000011521 glass Substances 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000007427 paired t-test Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 210000003462 vein Anatomy 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
- G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Electromagnetism (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于激光雷达的林窗立体结构量测方法，步骤1：激光雷达数据的获取步骤以及预处理步骤；所述的激光雷达数据是指由固定翼飞机或无人机搭载激光雷达仪获取遥感数据；每一个激光点的数据包括X，Y，Z坐标数据以及回波强度数据；步骤2：林冠层数据的处理步骤；步骤3：林窗识别步骤；步骤4：林窗立体结构测量步骤；该基于激光雷达的林窗立体结构量测方法易于实施，不受调查地林分类型和地面状况的影响，能快速、准确、大范围、多维度测量林窗结构。
Description
技术领域
本发明涉及一种基于激光雷达的林窗立体结构量测方法。
背景技术
林窗(canopy gap)指森林群落中主林层受人为干扰(择伐)或自然干扰(大风、雪、水灾、泥石流、虫害等灾害)在林冠层形成的不连续的林中空隙，是促进森林植被更新的重要空间。林窗研究作为森林生态学重要内容之一，具有重要的理论和实践意义。
林窗立体结构的测量是林窗研究的基础，林窗立体结构主要包括林窗面积、林窗形状和林窗边缘木高度。现阶段林窗测量方法可分为两类：一是基于实地测量的直接法，二是基于摄影技术或遥感的间接法。
基于实地测量的直接法，主要包括椭圆法、等角多边形、等角椭圆扇形法和三角形法等，基本思路是将单个林窗近似分割成多个同类型的几何形状，然后将多个几何形状面积求和以实现林窗总面积的估测。直接法费时费力，受外界因素影响大，精度也不高，只能估测林窗的面积，不能测量林窗形状和边缘木高度。基于摄影技术或遥感的间接法，主要包括基于摄影技术的平面相片法和半球面影像法，以及基于遥感技术的航片法。基于摄影技术的间接法采用装配普通镜头或鱼眼镜头在林窗内同一位置不同高度处垂直向上拍摄的2张影像(2次拍摄的方位角相同)，根据镜头投影原理和2次拍摄点的高差测量林窗立体结构。该类方法具有客观、精度较高以及可重复的优点，但对测量人员的摄影技术要求较高，特别是涉及鱼眼镜头的影响因素较多。
航片法采用航天或航空光学遥感影像，根据数字图像处理技术勾勒或识别林窗边缘，达到测量林窗的目的。航片法最大的优点是能在景观尺度上快速识别和测量森林中大量林窗的面积，但不能获取林窗立体结构特征，由于受部分冠层阴影的影响，航片法通常对小林窗的测量精度不高。
因此，有必要设计一种新的基于激光雷达的林窗立体结构量测方法。
发明内容
本发明所要解决的技术问题是提供一种基于激光雷达的林窗立体结构量测方法，该基于激光雷达的林窗立体结构量测方法易于实施。
发明的技术解决方案如下：
一种基于激光雷达的林窗立体结构量测方法,其特征在于，包括以下步骤：
步骤1：激光雷达数据的获取步骤以及预处理步骤；
所述的激光雷达数据是指由固定翼飞机或无人机搭载激光雷达扫描仪获取遥感数据；
每一个激光点的数据包括X，Y，Z坐标数据，以及回波强度数据；
预处理步骤如下：(1)读取激光点云原始数据(LAS格式)，获取每个点云的x轴数值、y轴数值、z轴数值、回波强度。(2)设定投影坐标系为墨卡托投影(UTM)方式，参考椭球为WGS84，投影带号T由调查地区的中央经度换算得到。(3)采用基于数据本身的航带间重叠区误差分析，提高点云数据定位精度，并修正航带的变形，其次，对各航带数据进行系统误差消除，使航带间平均高程差值在-2～2厘米之间。
步骤2：林冠层数据的处理步骤；
步骤3：林窗识别步骤；
步骤4：林窗立体结构测量步骤。
步骤2包括获取CHM、CSM和CPI数据的步骤；
所述的CHM是指冠层高度模型，CSM是指冠层梯度模型，CPI是指林冠空隙模型。
步骤3中，构建基于CHM、CSM和CPI数据的林窗识别模型，并确定林窗识别模型参数，并得到最终识别结果；最终识别结果为林窗识别结果栅格图；
步骤4中，基于步骤3中的最终识别结果，对林窗面积、形状和边缘木高度特征进行三维空间测量，得到测量结果。
所述的步骤2的林冠层数据的处理步骤中：
(1)基于自适应不规则三角网(TIN)滤波方法将激光点云分类成地表点云和冠层点云，(a)采用第一栅格数据插值法将地表点云插值生成分辨率为R1的数字高程模型(DEM)数据，记数字高程模型数据为DEM数据；(b)采用第二数据插值法将冠层点云插值生成分辨率为R2的冠层表面模型(DSM)数据，记冠层表面模型(DSM)数据为DSM数据；
(2)基于DEM数据和DSM数据，采用差值法生成冠层高度模型(CHM)数据；
即有CHM(x，y)＝DSM(x，y)—DEM(x，y)；
x，y分别为栅格(x，y)的横坐标和纵坐标；
CHM(x，y)是栅格(x，y)上的高程差值；
(3)在CHM数据的基础上采用拟合曲面法计算冠层梯度，生成分辨率为R3的冠层梯度模型(CSM)数据；
(4)利用遥感图像处理软件，根据冠层点云占总点云数量的比例生成林冠空隙模型(CPI)数据，林冠空隙模型(CPI)数据记为CPI数据；CPI数据的分辨率为R4；CPI的每一个元素代表预定面积(如每平方米)内的冠层点云占总点云数量的比例，取值范围0～100％，取值越大，越是代表林间空地；CPI的每一个元素代表预定面积(如每平方米)内高度在3米以下的冠层点云占总点云数量的比例，取值范围0～100％，取值越大，越是代表林间空地。
步骤(3)中，采用拟合曲面法计算冠层梯度，是指基于CHM数据，以N1×N1(N1是核元素的大小，一般可以是3、5、7等奇数)核元素进行栅格平滑处理：计算公式为：
式中Slope为梯度，slopewe为X方向上的梯度，slopesn为Y方向上的梯度；梯度取值范围0～90度。
| e5 | e2 | e6 |
| e1 | e | e3 |
| e8 | e4 | e7 |
针对上表的数据，有：
e以及e1-e8表示在3×3大小的栅格中9个不同位置的高程值(具体为高程差值)。
所述的步骤3中的林窗识别步骤为采用多条件组合法进行林窗识别：
记林窗模型为W，设林窗中某一个点(x，y)的取值为W(x，y)，有
式中，otherwise表示其余，&表示逻辑与。
W(x，y)＝1，表示点(x，y)对应的位置为林窗。
精确起见，CHM，CSM和CPI应该都采用相同的分辨率或换算程相同的分辨率，后文中有参数设定的表格，交代了分辨率都是设定为1米。
具体地：
(1)先借助专业激光雷达处理软件Terrasolid(v017)(现有技术)中自适应不规则三角网(TIN)滤波方法(现有主要激光点云滤波方法之一，其在森林区域的滤波效果最佳)将激光点云分类成地表点云和冠层点云，采用栅格数据插值法1(现有的地理信息科学中有大量数据插值算法，比较有代表性有克里金插值、反向距离加权插值法、样条插值法等，插值算法是通过一组具有Z值的分散点生成估计表面的高级地统计过程)将地表点云插值生成分辨率为R1的数字高程模型(DEM)，采用数据插值法2；数据插值法2为克里金插值、反向距离加权插值法、样条插值法中的一种)；1)为了优化栅格插值效果，优选采用不同插值方法对地表点云和冠层点云这两类数据进行操作，2)分辨率R1和R2也是根据处理对象有所变化的；后文也给出了相应的具体参数。将冠层点云插值生成分辨率为R2的冠层表面模型(DSM)；
(2)DSM与DEM差值法生成冠层高度模型(CHM)，即CHM＝DSM—DEM；两个栅格数据相减，即是同一栅格位置上两个高程值的差值。实际操作中DSM和DEM的分辨率是一致的，当DSM和DEM的分辨率不同时，首先是将不同分辨率转化成同一分辨率，即转化成较高的那一个分辨率。
(3)在CHM基础上采用拟合曲面法计算冠层梯度，生成分辨率为R3的冠层梯度模型(CSM)，拟合曲面法计算梯度即以N1×N1核元素进行栅格平滑处理：例如式中Slope为梯度，slopewe为X方向上的梯度，slopesn为Y方向上的梯度。梯度取值范围0～90度，林窗边缘的冠层梯度一般比较大。
| e5 | e2 | e6 |
| e1 | e | e3 |
| e8 | e4 | e7 |
e表示高程值，e以及e1-e8表示在3×3大小的栅格中9个不同位置的高程值。
(4)借助专业遥感图像处理软件The Environment for Visualizing Images(ENVI,v5.3)(现有技术)，根据冠层点云占总点云数量的比例生成林冠空隙模型(CPI)，即计算每平米高度在3米以下的点云占总点云数量的比例，取值范围0～100％，取值越大，越是代表林间空地。CPI分辨率为R4。
理论上说，点云总数量就是冠层点云数量加上地表点云数量。利用ENVI遥感软件可计算点云密度，实际中先计算冠层点云密度，然后是点云总密度，将这两个密度相除就得到比例，其中冠层点云、地表点云的分类采用前文表述的自适应不规则三角网(TIN)滤波方法。CPI模型中选择3米高度阈值，是因为3米以上点云大都是成树的树冠点云，林窗边缘的冠层大都小于3米，林窗中心的灌木或草本冠层点云高度都在1米左右，所以实践中根据3米以下的点云占总点云数量的比例来判断是否是林窗，CPI取值越大，越是代表林窗。
步骤3的林窗识别步骤中：
(1)林窗识别采用上述模型的多条件组合法。
多条件组合法数学形式：林窗＝CHM(冠层高度参数CH)+CSM(林窗边缘梯度参数CS)+CPI(冠层空隙参数CP)。
式中冠层高度参数CH：在冠层高度模型(CHM)中小于高度参数CH的视为林窗，反之亦然；
林窗边缘梯度参数CS：在冠层梯度模型(CSM)中小于梯度参数CS的视为林窗，反之亦然；
冠层空隙参数CP：在林冠空隙模型(CPI)中小于空隙参数CP的视为非林窗，反之亦然。
在冠层高度模型(CHM)上设定冠层高度参数CH阈值(米)，如在CHM中设定CH阈值为5米，即表示在CHM中高度小于5米的范围是林窗，高度大于5米的是林冠层；
在冠层梯度模型(CSM)上设定林窗边缘梯度参数CS阈值(度)，如CSM中的冠层梯度范围是0～90°，设定CS阈值为50°，即表示在CSM中小于50°的范围为林窗，大于50°的为林冠层；
在林冠空隙模型(CPI)上设定空隙参数CP阈值(％)，如CPI中林冠空隙范围是0～100％，设定CP阈值为60％，即表示在CPI中大于60％的范围为林窗，小于60％的为林冠层；
步骤4中，(1)借助专业地理信息系统处理软件ArcGIS(v10.2)(现有技术)将林窗识别结果栅格图与CHM叠加，并在CHM基础上对林窗立体结构进行测量。。林窗边缘木高度即为CHM上东西南北四个方向上冠层高度的平均值，单位为米。林窗形状指数为即为林窗面积(A)和周长(p)的函数，该指数为无量纲，取值范围≥1，取值越大，代表林窗形状越复杂，取值接近1，林窗近似圆形。(2)将测量结果以图形格式或文本格式存储。
本发明采用如下技术方案：
[1]遥感数据预处理模块。获取测量地区原始激光雷达遥感数据，按照标准格式读入计算机软件；利用地面控制点和几何校正数学模型矫正非系统因素产生的误差；同时将点云数据坐标系统转换成指定地图投影坐标系统，并根据植被分布情况，裁剪或拼接遥感数据。该模块主要包括激光雷达遥感数据的读取，几何校正，数据拼接等。
[2]林冠层处理模块。应用遥感图像处理方法生成能准确表征冠层水平和垂直方向特征的冠层高度模型(Canopy height model，CHM)；采用拟合曲面法在CHM上计算冠层边缘梯度以反映林窗边缘梯度变化特征；采用激光点云密度比较法计算冠层空隙率从冠层点云空间分布角度检验林窗特征。该模块主要包括冠层高度模型(CHM)，冠层梯度模型(Canopy slope model，CSM)和林冠空隙模型(Canopy porosity index，CPI)等。
[3]林窗识别模块。在上述模块[2]基础上，构建由CHM、CSM和CPI组成的林窗识别模型，根据样地检测数据和模型参数优化算法确定林窗识别模型参数，在得到林窗对象(待定)的识别结果上进行识别后处理，排除误差得到最终识别结果；应用栅格混淆矩阵进行识别精度评价。
[4]林窗立体结构测量模块。根据林窗识别模块[3]的结果，应用计算机图形学方法在相关软件中对林窗面积、形状和边缘木高度等特征进行三维空间测量，测量结果以图形文件或文本文件形式存储。
本发明利用激光雷达采样精度高(点云垂直精度约0.6米，平面精度约0.5米，平均密度为2-8个/平方米)及其能对林冠层进行分层监测的特点，应用相关的遥感和计算机技术，实现不同类型森林(如南方针阔混交林、阔叶林或北方针叶林等)中林窗立体结构的高精度、多尺度、多维度的测量。
有益效果：
本发明提出一种新的基于遥感技术的林窗立体结构测量方法，利用激光雷达高精度采样及其分层监测的特点，实现林窗立体结构的测量。本发明不受调查地林分类型和地面状况的影响，能快速、准确、大范围、多维度测量林窗结构，能够应用在林业、生态监测、以及灾害监测等遥感应用部门。
本发明的使用机载激光雷达获取数据，机载激光雷达(light detection andranging，LiDAR)是一种新兴的主动式遥感技术，能在多时空尺度上获取森林生态系统高精度的植被结构信息、三维地形特征。激光雷达在林冠层快速、准确监测和模拟方面具有巨大潜力，虽然国际上关于激光雷达在林窗研究中应用已有介绍，但是国内还鲜有利用激光雷达数据进行林窗立体结构测量方面的报道。
[1]与传统的基于实地测量的直接法相比较，本发明具有客观、快速、精度高以及可重复的优点。
[2]与基于摄影技术的间接法相比较，本发明不受测量场地环境限制和个人摄影水平高低的影响，能对单个或景观尺度内多个林窗的立体结构进行快速、准确、流程化的测量。
[3]与基于光学航片的间接法相比较，本发明利用激光雷达技术的特点，不仅能准确测量不同大小、不同复杂度的林窗面积，还能测量林窗边缘木平均高度和林窗形状。
附图说明
图1为基于激光雷达的林窗立体结构量测方法对应的流程图；
图2为原始点云示意图；
图3为数字高程模型(DEM)示意图；
图4为数字表面模型(DSM)示意图；
图5为冠层高度模型(CHM)示意图；
图6为冠层梯度模型(CSM)示意图；
图7为林冠空隙模型(CPI)示意图；
图8为CHM高度二值图；
图9为CSM梯度二值图；
图10为CPI百分比二值图；
图11为识别后的林窗多边形示意图；
图12为图11中A区域对应的林窗几何结构量测示意图。
具体实施方式
以下将结合附图和具体实施例对本发明做进一步详细说明：
实施例1：如图1-12，具体实施分为4步骤：
(一)激光雷达数据预处理模块
(1)读取激光点云原始数据(LAS格式)，获取每个点云的x轴数值、y轴数值、z轴数值、回波强度。(2)设定投影坐标系为墨卡托投影(UTM)方式，参考椭球为WGS84，投影带号T由调查地区的中央经度换算得到。(3)采用基于数据本身的航带间重叠区误差分析，提高点云数据定位精度，并修正航带的变形，其次，对各航带数据进行系统误差消除，使航带间平均高程差值在-2～2厘米之间。
由固定翼飞机或无人机搭载激光雷达扫描仪获取遥感数据，选择晴朗少云天气进行飞行，根据试验地区范围和地形高度、梯度等情况，沿着地势设计多条航线，航带间重叠率为15～25％，飞行高度800～1100m。激光雷达系统通过发射单束窄带激光脉冲并接收反射的回波信号记录地物空间信息，测量精度±20mm，光斑大小约25cm，波形采样间隔1ns，脉冲重复频率最高100KHz，激光束散射角0.5mrad，扫描角范围±22.5deg，扫描速度5-160行/s，扫描角精度0.0025deg，多对象分辨率0.6m，平均点云密度2-10个/m2。
激光点云存储格式，每个激光点对应一行记录
(二)林冠层处理模块
(1)借助专业激光雷达处理软件Terrasolid(v017)(现有技术)中自适应不规则三角网(TIN)滤波方法(现有主要激光点云滤波方法之一，其在森林区域的滤波效果最佳)将激光点云分类成地表点云和冠层点云，采用栅格数据插值法1(现有的地理信息科学中有大量数据插值算法，比较有代表性有克里金插值、反向距离加权插值法、样条插值法等，插值算法是通过一组具有Z值的分散点生成估计表面的高级地统计过程)将地表点云插值生成分辨率为R1的数字高程模型(DEM)，采用数据插值法2将冠层点云插值生成分辨率为R2的冠层表面模型(DSM)；(2)DSM与DEM差值法生成冠层高度模型(CHM)，即CHM＝DSM—DEM；(3)在CHM基础上采用拟合曲面法计算冠层梯度，生成分辨率为R3的冠层梯度模型(CSM)，拟合曲面法计算梯度即以N1×N1核元素进行栅格平滑处理：例如式中Slope为梯度(又名坡度)，slopewe为X方向上的梯度，slopesn为Y方向上的梯度。梯度取值范围0～90度，林窗边缘的冠层梯度一般比较大。
(4)借助专业遥感图像处理软件The Environment for Visualizing Images(ENVI,v5.3)(现有技术)，根据林冠点云占总点云数量的比例生成林冠空隙模型(CPI)，即计算每平米高度在3米以下的植被点云占总点云数量的比例，取值范围0～100％，取值越大，越是代表林间空地。CPI分辨率为R4。
上述模型(DEM、DSM、CHM、CSM以及CPI)的生成，可自行编程序实现也可通过相关商业软件完成。典型商业软件如Terrasolid、ENVI等。林冠层处理模块中已介绍了各种模型的生成方法。
(三)林窗识别模块
(1)林窗识别采用上述模型的多条件组合法。
该方法可以表示为：林窗＝CHM(冠层高度参数CH)+CSM(林窗边缘梯度参数CS)+CPI
(冠层空隙参数CP)。
即：记林窗模型为W，设林窗中某一个点(x，y)的取值为W(x，y)，有
otherwise表示其余，&表示逻辑与。
W(x，y)＝1，表示点(x，y)对应的位置为林窗。
式中冠层高度参数CH：在冠层高度模型(CHM)中小于高度参数CH的视为林窗，反之亦然；
林窗边缘梯度参数CS：在冠层梯度模型(CSM)中小于梯度参数CS的视为林窗，反之亦然；
冠层空隙参数CP：在林冠空隙模型(CPI)中小于空隙参数CP的视为非林窗，反之亦然。
在冠层高度模型(CHM)上设定冠层高度参数CH阈值(米)，如在CHM中设定CH阈值为5米，即表示在CHM中高度小于5米的范围是林窗，高度大于5米的是林冠层；
在冠层梯度模型(CSM)上设定林窗边缘梯度参数CS阈值(度)，如CSM中的冠层梯度范围是0～90°，设定CS阈值为50°，即表示在CSM中小于50°的范围为林窗，大于50°的为林冠层；
在林冠空隙模型(CPI)上设定空隙参数CP阈值(％)，如CPI中林冠空隙范围是0～100％，设定CP阈值为60％，即表示在CPI中大于60％的范围为林窗，小于60％的为林冠层；
(2)林窗识别模型参数的确定，模型参数由受试者工作特征曲线(receiveroperating characteristic curve，ROC曲线)阈值决定，以实测训练样本为基础，分别对冠层高度参数CH、林窗边缘梯度参数CS和冠层空隙参数CP进行ROC曲线分析。通常将模型参数取值分成10等级(阈值)，每个阈值对应ROC曲线上的真正类率(纵轴值)和假正类率(横轴值)，选取真正类率高同时假正类率低的临界值作为模型参数的取值。例如：林窗(待定)＝CHM(5米)+CSM(50度)+CPI(60％)，即同时满足三个条件(在CHM上高度小于5米，在CSM上梯度小于50度，在CPI上空隙比例小于60％)的为林窗(待定)。
根据冠层高度参数CH，在CHM上进行二值化处理，即满足条件的范围为林窗，设定布尔值1；
同理，将CSM和CPI进行二值化处理，最后将3个二值化图像进行“交集”运算，交集为布尔值1的范围视为林窗(待定)结果。
(3)林窗识别后处理。采用形态学方法对林窗(待定)进行检测，借助专业遥感图像处理软件ENVI(v5.3)(现有技术)在栅格分辨率设定N2×N2圆形形态学算子的基础上，利用“腐蚀”算法消除林木之间的小间隙(2米以内)或林冠中的小“空洞”，利用“膨胀”算法还原得到真实林冠层空隙，最后将大于25平方米的冠层空隙确定为林窗。
(4)借助专业遥感图像处理软件ENVI(v5.3)(现有技术)，应用栅格混淆矩阵进行林窗识别精度评价与分析，具体评价指标是Kappa系数、总体精度、错判比例和漏判比例。
具体实现可以现有商业软件或自行编程实现。
Kappa系数：
它是通过把所有地表真实分类中的像元总数(N)乘以混淆矩阵对角线(xkk)的和，再减去某一类地表真实像元总数与被误分成该类像元总数之积对所有类别求和的结果，再除以总像元数的平方减去某一类中地表真实像元总数与该类中被误分成该类像元总数之积对所有类别求和的结果所得到的。代表着分类与完全随机的分类产生错误减少的比例。
Kappa系数计算公式：
总体分类精度(OA)：指被正确分类的类别像元数与总的类别个数的比值；OA值虽然能很好的表征分类精度，但是对类别像元个数极度不平衡的多类地物来说，其值收到像元数据较多类别的影响较大，不能很好的表征每个类别地物。
错判比例：被分为地物类别实际属于另一类的像元比例。
漏判比例：指本身属于地表真实分类，但没有分到相应类别中的像元比例。
(四)林窗立体结构测量模块
(1)借助专业地理信息系统处理软件ArcGIS(v10.2)(现有技术)将林窗识别结果栅格图与CHM叠加，并在CHM基础上对林窗立体结构进行测量。林窗面积即为林窗识别结果栅格图的范围，单位为平方米。林窗边缘木高度即为CHM上东西南北四个方向上冠层高度的平均值，单位为米。林窗形状指数为即为林窗面积(A)和周长(p)的函数，该指数为无量纲，取值范围≥1，取值越大，代表林窗形状越复杂，取值接近1，林窗近似圆形。(2)将测量结果以图形格式或文本格式存储。
实施例(2个)
样地情况：
研究地区位于湖南省雪峰山东南部的武冈林场(26°25′—27°00′N,110°22′—113°3′E),地处祁邵丘陵区,三面环山,南高北低,海拔多在550～850m.该区属中亚热带季风湿润气候,雨量充沛,年均气温16℃,年均降水量1400mm,年无霜期250d,土壤以山地红壤和黄棕壤为主。主要植被类型有暖性针叶林、常绿阔叶林、针阔混交林和落叶阔叶林。实地林窗调查采用半球面影像法，根据鱼眼镜头投影原理确定林窗面积；采用角规法或伸缩式测高器测量林窗边界木高度；DGPS或全站仪测量林窗中心海拔和位置。
样地1，林分群落上层乔木有甜槠(Castanopsis eyrei)、青冈栎(Cyclobalanopsis glauca)、猴樟(Cinnamomum bodinieri)、马尾松(Pinus massoniana)等。主要灌木有南方荚蒾(Viburnum fordiae)、山茶(Camellia japonica L.)、胡枝子(Lespedeza bicolor)、映山红(Rhododendron simsii)、微毛柃(Eurya hebeclados)等。
样地2，林分群落上层乔木有杉木(Cunninghamia lanceolata)、马尾松(Pinusmassoniana)。草本主要有蕨(Pteridium aquilinum)和芒萁(Dicranopteris dichotoma)等。
遥感数据情况：遥感数据来自湖南省自然科学基金项目(1、南方生态公益林森林空间结构多源遥感量化分析(No:2015JJ2201)；2、顾及复相干性分布机制的极化干涉SAR森林树高反演研究(No:2017JJ3515))
关键参数设定：
林窗识别率：
本方案针对不同林分(样地1：常绿阔叶林，样地2：人工针叶林)中林窗识别率都很高，在栅格混淆矩阵中Kappa系数都大于0.8，说明林窗栅格分类效果好。林窗面积与实地鱼眼镜头测量结果比较：
Sig.:配对差值正态分布显著性；P:配对T检验双侧显著性；y:LiDAR测量；x:实地测量
借助专业数据统计分析软件Statistical Product and Service Solutions(SPSS，v19)，利用配对T检验(paired t-test)比较LiDAR与实测林窗面积的差异，两者配对的面积差值正态分布显著性Sig.都大于0.05，说明利用配对T检验法具有统计学意义。配对T检验的原假设是面积差值的分布符合平均值为0的t分布，双侧显著性P都大于0.05，说明LiDAR与实测林窗面积没有显著差异。分别以LiDAR和实测林窗面积为因变量和自变量做线性回归分析，得到的相关系数都很高，说明两者具有强线性相关性。回归方程的斜率大于1，说明LiDAR测量的林窗面积平均比实地测量的大，就增加的程度而言，阔叶林林窗稍微大于针叶林林窗。
林窗边缘木高度与实地伸缩式测高器测量结果比较：
用配对T检验法具体比较LiDAR估测与实地测量林窗边界木高的差异。两者配对的高度差值正态分布显著性Sig.大于0.05，说明利用配对T检验法具有统计学意义，双侧显著性为P大于0.05，维持原检验假设，即两者测量的林窗边界木高没有显著差异。回归模型的相关系数R2有较高值，说明两者具有较强线性相关性，总体而言，LiDAR估测比实地测量的略低，阔叶林林窗边缘木高度误差程度小于针叶林。
林窗形状系数与实地测量结果比较：
根据T检验的差值正态分布显著性和双侧显著性结果，说明LiDAR计算的林窗形状指数与实地测量值没有显著差异，但LiDAR计算的通常小于实地测量值。
总之，利用本发明的方法测量林窗面积、边缘木高度和形状指数与实地测量值之间没有显著差异。具体而言，①遥感测量的林窗面积通常偏大，阔叶林偏大的程度比针叶林的稍微明显些；②遥感测量的林窗边缘木高度通常偏小，针叶林偏小的程度比阔叶林的稍微明显些；③遥感测量计算的林窗形状指数通常略小与实地测量数据计算值。
主要原因可能是激光点云的分布密度还不够多以及地形起伏、激光扫描角度等因素，造成冠层边缘部分点云数量偏少，使得遥感监测的林窗边缘比实际的稍微大些，另外遥感监测的角度是近似垂直向下，而实地测量的半球面影像法角度是垂直向上的，这种角度差异以及冠层生态特征，使得两种方法测量结果稍有偏差，这种情况在复层型阔叶林中可能更加明显些。遥感测量的林窗边缘木高度偏小的原因也有激光点云密度有关，树冠顶部有可能没有点云，另外栅格化处理过程高度信息有可能缺失。遥感监测的林窗形状指数偏小，说明遥感监测的林窗形状复杂度比实际测量的要低，这主要是遥感数据的分辨率与实地测量的半球面影像法的分辨率上有差距。解决上述问题，就有效的途径是提高激光雷达点云的密度，并结合高分辨率的遥感影像进行林窗监测。
本发明方案中是直接在CHM中识别林窗；现有技术中，林窗判别条件为单一的，即CHM高度在特定阈值(例如5米)以下的部分判定为林窗，而发明方案中判定条件是多个组合，即从冠层高度差、林窗边缘梯度以及冠层点云比例等3方面综合考虑，考虑的影响因素更加全面。现有技术中，比较适合人工针叶林林窗特征提取，而发明方案适用更多类型的林分，如常绿阔叶林或针阔混交林等，而且林窗判别和特征提取效果要比替代方案稳定、准确。
Claims (7)
1.一种基于激光雷达的林窗立体结构量测方法，其特征在于，包括以下步骤：
步骤1：激光雷达数据的获取步骤以及预处理步骤；
所述的激光雷达数据是指由固定翼飞机或无人机搭载激光雷达扫描仪获取遥感数据；
每一个激光点的数据包括X，Y，Z坐标数据，以及回波强度数据；
步骤2：林冠层数据的处理步骤；
步骤3：林窗识别步骤；
步骤4：林窗立体结构测量步骤。
2.根据权利要求1所述的基于激光雷达的林窗立体结构量测方法,其特征在于，步骤2包括获取CHM、CSM和CPI数据的步骤；
所述的CHM是指冠层高度模型，CSM是指冠层梯度模型，CPI是指林冠空隙模型。
3.根据权利要求2所述的基于激光雷达的林窗立体结构量测方法,其特征在于，步骤3中，构建基于CHM、CSM和CPI数据的林窗识别模型，并确定林窗识别模型参数，并得到最终识别结果；最终识别结果为林窗识别结果栅格图。
4.根据权利要求3所述的基于激光雷达的林窗立体结构量测方法,其特征在于，步骤4中，基于步骤3中的最终识别结果，对林窗面积、形状和边缘木高度特征进行三维空间测量，得到测量结果。
5.根据权利要求2所述的基于激光雷达的林窗立体结构量测方法,其特征在于，所述的步骤2的林冠层数据的处理步骤中：
(1)基于自适应不规则三角网(TIN)滤波方法将激光点云分类成地表点云和冠层点云，(a)采用第一栅格数据插值法将地表点云插值生成分辨率为R1的数字高程模型(DEM)数据，记数字高程模型数据为DEM数据；(b)采用第二数据插值法将冠层点云插值生成分辨率为R2的冠层表面模型(DSM)数据，记冠层表面模型(DSM)数据为DSM数据；
(2)基于DEM数据和DSM数据，采用差值法生成冠层高度模型(CHM)数据；
即有CHM(x，y)＝DSM(x，y)—DEM(x，y)；
x，y分别为栅格(x，y)的横坐标和纵坐标；
CHM(x，y)是栅格(x，y)上的高程差值；
(3)在CHM数据的基础上采用拟合曲面法计算冠层梯度，生成分辨率为R3的冠层梯度模型(CSM)数据；
(4)利用遥感图像处理软件，根据冠层点云占总点云数量的比例生成林冠空隙模型(CPI)数据，林冠空隙模型(CPI)数据记为CPI数据；CPI数据的分辨率为R4；CPI的每一个元素代表预定面积(如每平方米)内的冠层点云占总点云数量的比例，取值范围0～100％，取值越大，越是代表林间空地。
6.根据权利要求5所述的基于激光雷达的林窗立体结构量测方法,其特征在于，步骤(3)中，采用拟合曲面法计算冠层梯度，是指基于CHM数据，以N1×N1核元素进行栅格平滑处理：计算公式为：
式中Slope为梯度，slopewe为X方向上的梯度，slopesn为Y方向上的梯度；梯度取值范围0～90度。
针对上表的数据，有：
e以及e1-e8表示在3×3大小的栅格中9个不同位置的高程值。
7.根据权利要求6所述的基于激光雷达的林窗立体结构量测方法,其特征在于，所述的步骤3中的林窗识别步骤为采用多条件组合法进行林窗识别：
记林窗模型为W，设林窗中某一个点(x，y)的取值为W(x，y)，有
otherwise表示其余，&表示逻辑与。
W(x，y)＝1，表示点(x，y)对应的位置为林窗。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710574787.6A CN107479065B (zh) | 2017-07-14 | 2017-07-14 | 一种基于激光雷达的林窗立体结构量测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710574787.6A CN107479065B (zh) | 2017-07-14 | 2017-07-14 | 一种基于激光雷达的林窗立体结构量测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107479065A true CN107479065A (zh) | 2017-12-15 |
| CN107479065B CN107479065B (zh) | 2020-09-11 |
Family
ID=60595647
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710574787.6A Expired - Fee Related CN107479065B (zh) | 2017-07-14 | 2017-07-14 | 一种基于激光雷达的林窗立体结构量测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107479065B (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108132096A (zh) * | 2017-12-20 | 2018-06-08 | 中南林业科技大学 | 一种基于激光雷达的林窗太阳辐射监测方法 |
| CN108375985A (zh) * | 2018-02-06 | 2018-08-07 | 南京和图地理信息工程有限公司 | 一种土地三维规划设计平台及其设计方法 |
| CN109407113A (zh) * | 2018-11-19 | 2019-03-01 | 中南林业科技大学 | 一种基于机载激光雷达的林窗时空变化监测与量化方法 |
| CN110110641A (zh) * | 2019-04-29 | 2019-08-09 | 中国水利水电科学研究院 | 一种流域性洪涝场景的无人机监测方法及系统 |
| CN110569786A (zh) * | 2019-09-06 | 2019-12-13 | 中国农业科学院农业资源与农业区划研究所 | 一种基于无人机数据采集的果树识别和数量监测方法及系统 |
| CN111913169A (zh) * | 2019-05-10 | 2020-11-10 | 北京四维图新科技股份有限公司 | 激光雷达内参、点云数据的修正方法、设备及存储介质 |
| CN113009481A (zh) * | 2021-01-15 | 2021-06-22 | 扬州哈工科创机器人研究院有限公司 | 一种基于干涉sar雷达的森林地物成像反演方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2040093A2 (en) * | 2003-07-17 | 2009-03-25 | NovAtel Inc. | A seismic measuring system including GPS receivers |
| CN101413875A (zh) * | 2008-11-25 | 2009-04-22 | 中山大学 | 一种树木冠层分析仪的数据采集装置 |
| CN104502919A (zh) * | 2015-01-13 | 2015-04-08 | 南京大学 | 利用机载激光雷达点云提取城市植被三维覆盖的方法 |
| CN105866792A (zh) * | 2016-05-31 | 2016-08-17 | 中国科学院遥感与数字地球研究所 | 一种新的星载激光雷达树高提取方法 |
-
2017
- 2017-07-14 CN CN201710574787.6A patent/CN107479065B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2040093A2 (en) * | 2003-07-17 | 2009-03-25 | NovAtel Inc. | A seismic measuring system including GPS receivers |
| CN101413875A (zh) * | 2008-11-25 | 2009-04-22 | 中山大学 | 一种树木冠层分析仪的数据采集装置 |
| CN104502919A (zh) * | 2015-01-13 | 2015-04-08 | 南京大学 | 利用机载激光雷达点云提取城市植被三维覆盖的方法 |
| CN105866792A (zh) * | 2016-05-31 | 2016-08-17 | 中国科学院遥感与数字地球研究所 | 一种新的星载激光雷达树高提取方法 |
Non-Patent Citations (2)
| Title |
|---|
| 刘峰等: "基于机载激光雷达的中亚热带常绿阔叶林林窗特征", 《应用生态学报》 * |
| 周梦维等: "基于机载小光斑全波形LIDAR的作物高度反演", 《农业工程学报》 * |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108132096A (zh) * | 2017-12-20 | 2018-06-08 | 中南林业科技大学 | 一种基于激光雷达的林窗太阳辐射监测方法 |
| CN108375985A (zh) * | 2018-02-06 | 2018-08-07 | 南京和图地理信息工程有限公司 | 一种土地三维规划设计平台及其设计方法 |
| CN109407113A (zh) * | 2018-11-19 | 2019-03-01 | 中南林业科技大学 | 一种基于机载激光雷达的林窗时空变化监测与量化方法 |
| CN110110641A (zh) * | 2019-04-29 | 2019-08-09 | 中国水利水电科学研究院 | 一种流域性洪涝场景的无人机监测方法及系统 |
| CN111913169A (zh) * | 2019-05-10 | 2020-11-10 | 北京四维图新科技股份有限公司 | 激光雷达内参、点云数据的修正方法、设备及存储介质 |
| CN111913169B (zh) * | 2019-05-10 | 2023-08-22 | 北京四维图新科技股份有限公司 | 激光雷达内参、点云数据的修正方法、设备及存储介质 |
| CN110569786A (zh) * | 2019-09-06 | 2019-12-13 | 中国农业科学院农业资源与农业区划研究所 | 一种基于无人机数据采集的果树识别和数量监测方法及系统 |
| CN110569786B (zh) * | 2019-09-06 | 2022-03-29 | 中国农业科学院农业资源与农业区划研究所 | 一种基于无人机数据采集的果树识别和数量监测方法及系统 |
| CN113009481A (zh) * | 2021-01-15 | 2021-06-22 | 扬州哈工科创机器人研究院有限公司 | 一种基于干涉sar雷达的森林地物成像反演方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107479065B (zh) | 2020-09-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107479065A (zh) | 一种基于激光雷达的林窗立体结构量测方法 | |
| Akay et al. | Using LiDAR technology in forestry activities | |
| CN113034689A (zh) | 基于激光点云的地形三维模型及地形图构建方法和系统、存储介质 | |
| CN113673181A (zh) | 一种基于多源风场数据的机场区域风切变智能识别方法 | |
| CN110427857A (zh) | 一种基于遥感数据融合的输电线路地质灾害分析方法 | |
| Farid et al. | Using airborne lidar to predict Leaf Area Index in cottonwood trees and refine riparian water-use estimates | |
| CN103363962A (zh) | 一种基于多光谱影像的湖泊水储量遥感估算方法 | |
| CN110988909B (zh) | 基于tls进行高寒脆弱区沙地植被的植被盖度测定方法 | |
| An et al. | Three-dimensional point cloud based sky view factor analysis in complex urban settings. | |
| Kent et al. | Urban morphology parameters from global digital elevation models: Implications for aerodynamic roughness and for wind-speed estimation | |
| Packalen et al. | Edge-tree correction for predicting forest inventory attributes using area-based approach with airborne laser scanning | |
| CN111091079A (zh) | 基于tls的高寒脆弱区植被优势单株结构参数测定方法 | |
| CN106780586B (zh) | 一种基于地面激光点云的太阳能潜力评估方法 | |
| Fu et al. | Assessment of approaches for monitoring forest structure dynamics using bi-temporal digital aerial photogrammetry point clouds | |
| CN108959705A (zh) | 一种预测大面积亚热带森林生物量的混合效应模型 | |
| CN108981616A (zh) | 一种由无人机激光雷达反演人工林有效叶面积指数的方法 | |
| Ferraz et al. | Canopy density model: A new ALS-derived product to generate multilayer crown cover maps | |
| CN109146951A (zh) | 一种基于无人机激光雷达孔隙度模型估测银杏人工林叶面积指数的方法 | |
| Fu et al. | Prediction of subtropical forest parameters using airborne laser scanner | |
| CN108132096B (zh) | 一种基于激光雷达的林窗太阳辐射监测方法 | |
| Sun et al. | Forest canopy closure estimation in greater khingan forest based on gf-2 data | |
| Jia et al. | Effects of point density on DEM accuracy of airborne LiDAR | |
| Guo et al. | Construction of 3D landscape indexes based on oblique photogrammetry and its application for islands | |
| Straatsma | Quantitative mapping of hydrodynamic vegetation density of floodplain forests under leaf-off conditions using airborne laser scanning | |
| Jin et al. | UAV-RGB-image-based aboveground biomass equation for planted forest in semi-arid Inner Mongolia, China |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200911 |