一种从激光雷达高斯回波数据中提取树高的方法

申请号 CN201310598367.3 申请日 2013-11-22 公开(公告)号 CN104062644A 公开(公告)日 2014-09-24
申请人 董立新; 发明人 董立新;
摘要 本 发明 涉及遥感测绘领域,提供了一种从 激光雷达 高斯回 波数 据中提取树高的方法。该方法包括:选取与树高待测地点同属于一个地区类型的若干校准参考地点,根据与该若干地点一一对应的实测树高值H、数字高程模型下的地形指数g、从激光雷达高斯回波数据中提取得到的 波形 质心到回波结束 位置 的距离wc和边缘长度l,对H=a0×(wc-a1×g+a2×l)进行函数拟合,得到一组最佳拟合参数a0、a1和a2;对于该树高待测地点,根据波形质心到回波结束位置的距离w’c、边缘长度l’、该地点的地形指数g’和所述最佳拟合参数a0、a1和a2由公式H’=a0×(w’c-a1×g’+a2×l’)计算得出该地点的测量树高值H’。本发明可以实现坡度较大地区激光雷达高斯回波数据中树高的提取。
权利要求

1.一种从激光雷达高斯回波数据中提取树高的方法,其特征在于,该方法包括:
步骤101:选取与树高待测地点同属于一个地区类型的若干校准参考地点,并获取该若干校准参考地点处的实测树高值;
步骤102:获取该若干地点的激光雷达高斯回波数据和数字高程模型下与该若干地点一一对应的地形指数g;
步骤103:从所述激光雷达高斯回波数据中提取得到波形质心到回波结束位置的距离wc和边缘长度l;
步骤104:以H=a0×(wc-a1×g+a2×l)为拟合函数对该对应若干地点的若干组H、g、wc和l进行函数拟合,得到一组最佳拟合参数a0、a1和a2;
步骤105:对于该树高待测地点,先获取其激光雷达高斯回波数据,并从中提取波形质心到回波结束位置的距离w’c和边缘长度l’;
步骤106:结合数字高程模型下的该地点的地形指数g’和所述最佳拟合参数a0、a1和a2由公式H’=a0×(w’c-a1×g’+a2×l’)计算得出该地点的树高测量值H’。
2.根据权利要求1所述的方法,其特征在于,所述提取波形质心到回波结束位置的距离包括:
波形质心的位置为植被回波波形面积一半处所在的横坐标位置;
设定高斯回波结束的阈值为结束噪声的平均值加上其标准偏差的和,从最后一个回波位置开始向后搜索,如果波形数据的值小于高斯回波结束的阈值,则确定此处的横坐标为回波结束位置;
所述提取波形质心到回波结束位置的距离为所述波形质心的位置与所述回波结束位置之差。
3.根据权利要求1所述的方法,其特征在于,所述提取边缘长度包括:
设定高斯回波结束的阈值为结束噪声的平均值加上其标准偏差的和,从最后一个回波位置开始向后搜索,如果波形数据的值小于高斯回波结束的阈值,则确定此处的横坐标为回波结束位置;
从所述回波结束位置开始向前查找一个回波位置,使其与所述回波结束位置之差大于激光脉冲的半宽,则所述边缘长度大小等于该回波位置与所述回波结束位置之差。
4.根据权利要求1所述的方法,其特征在于,对于某地点进行的所述获取地形指数g步骤包括:
对应该地点,获取在数字高程模型中5×5大小的采样窗口下地表高程的平均大小作为地形指数g。
5.根据权利要求1所述的方法,其特征在于,所述进行函数拟合包括使用Lenvenberg-Marquardt非线性拟合方法进行函数拟合。
6.根据权利要求1所述的方法,其特征在于,所述同属于一个地区类型包括同属于一个坡度类型。
7.根据权利要求1所述的方法,其特征在于,所述同属于一个地区类型包括同属于一个植被类型。

说明书全文

一种从激光雷达高斯回波数据中提取树高的方法

技术领域

[0001] 本发明涉及遥感测绘领域,具体涉及一种从激光雷达高斯回波数据中提取树高的方法。

背景技术

[0002] 激光雷达作为采用光电探测技术手段的主动遥感设备,以其高分辨率的特点广泛应用于林业、测绘学等领域。例如GLAS(Geoscience Laser Altimeter System,地球科学激光测高仪系统),搭载于美国2003年发射的科学试验卫星ICEsat上的第一台星载激光雷达传感器,能连续获取大气和地面的回波数据。其中的数字化仪记录从卫星到地表765km范围内的回波信号(时间-强度曲线),经过滤波与分析可提取发射脉冲与表面回波。其1064nm波段回波数据包含了包括高度信息在内的垂直结构信息。研究者将激光雷达回波用数学形式表达为数个高斯回波函数的叠加加上一个估计偏差,每个高斯回波数据分量的波峰位置对应一个高度信息,从而可从激光雷达高斯回波数据中提取如树高等植被垂直结构信息,透过该类数据可以很好地评价全球植被生物量和循环。
[0003] 一般情况下,现有的树高提取方案中多数取高斯回波的初始位置与地面回波位置(均在高斯回波数据中判断找出)之差为树高提取结果,这一结果被普遍认为是与森林冠层高度最接近的计算结果。
[0004] 但是有研究表明,坡度达到一定程度会对所述高斯回波数据产生造成很大影响(Harding,1994;Heyder,2005)。一般的,当地面存在坡度时,地面回波波峰位置不明显,此时直接使用该方法从波形数据中提取树高会造成很大的误差。其原因主要在于,当地形坡度超过一定范围时,地面回波与植被回波相互混合,加之激光雷达传感器的衰减造成回波获取与分解都很困难。尤其是目前很多这一方面的实验都是在北美洲比较平坦的地势下获取数据并完成的,而相比较而言,国内局部地区的地势起伏更大,现有的数据提取算法上可能存在一些不适用的情形。因此,根据激光雷达高斯回波数据对地势起伏大的地区,比如山地地区的森林调查与研究是目前的研究难点之一。

发明内容

[0005] (一)解决的技术问题
[0006] 针对现有技术的不足,本发明提供一种从激光雷达高斯回波数据中提取树高的方法,实现了坡度较大地区激光雷达高斯回波数据中树高的提取。
[0007] (二)技术方案
[0008] 一种从激光雷达高斯回波数据中提取树高的方法,其特征在于,该方法包括:
[0009] 步骤101:选取与树高待测地点同属于一个地区类型的若干校准参考地点,并获取该若干校准参考地点处的实测树高值;
[0010] 步骤102:获取该若干地点的激光雷达高斯回波数据和数字高程模型下与该若干地点一一对应的地形指数g;
[0011] 步骤103:从所述激光雷达高斯回波数据中提取得到波形质心到回波结束位置的距离wc和边缘长度l;
[0012] 步骤104:以H=a0×(wc-a1×g+a2×l)为拟合函数对该对应若干地点的若干组H、g、wc和l进行函数拟合,得到一组最佳拟合参数a0、a1和a2;
[0013] 步骤105:对于该树高待测地点,先获取其激光雷达高斯回波数据,并从中提取波形质心到回波结束位置的距离w’c和边缘长度l’;
[0014] 步骤106:结合数字高程模型下的该地点的地形指数g’和所述最佳拟合参数a0、a1和a2由公式H’=a0×(w’c-a1×g’+a2×l’)计算得出该地点的树高测量值H’。
[0015] 优选地,所述提取波形质心到回波结束位置的距离包括:波形质心的位置为植被回波波形面积一半处所在的横坐标位置;设定高斯回波结束的阈值为结束噪声的平均值加上其标准偏差的和,从最后一个回波位置开始向后搜索,如果波形数据的值小于高斯回波结束的阈值,则确定此处的横坐标为回波结束位置;所述提取波形质心到回波结束位置的距离为所述波形质心的位置与所述回波结束位置之差。
[0016] 优选地,所述提取边缘长度包括:设定高斯回波结束的阈值为结束噪声的平均值加上其标准偏差的和,从最后一个回波位置开始向后搜索,如果波形数据的值小于高斯回波结束的阈值,则确定此处的横坐标为回波结束位置;从所述回波结束位置开始向前查找一个回波位置,使其与所述回波结束位置之差大于激光脉冲的半宽,则所述边缘长度大小等于该回波位置与所述回波结束位置之差。
[0017] 优选地,对于某地点进行的所述获取地形指数g步骤包括:对应该地点,获取在数字高程模型中5×5大小的采样窗口下地表高程的平均大小作为地形指数g。
[0018] 优选地,所述进行函数拟合包括使用Lenvenberg-Marquardt非线性拟合方法进行函数拟合。
[0019] 优选地,所述同属于一个地区类型包括同属于一个坡度类型。
[0020] 优选地,所述同属于一个地区类型包括同属于一个植被类型。
[0021] (三)有益效果
[0022] 本发明至少具有如下的有益效果:
[0023] 本发明采用的公式中加入了树高待测地点处的地形指数这一参量,也就是对于树高提取问题中的地形效果也给予了一定考虑,使不同坡度下的树高提取问题都有对应的不同程度的该地点地形指数带来的修正,使树高提取精度对于坡度的敏感性降低,提高了坡度较大地区激光雷达高斯回波数据中树高的提取算法的精度;
[0024] 而且本发明在植被回波的选择中使用波形质心位置,有研究结果表明,这一处理使提取的树高值对坡度的敏感性较低,对坡度较大地区结果相对比较稳定(Heyder,2005);
[0025] 本发明还引入边缘长度l对较大坡度地区地面回波半宽增大的问题进行修正,改进后取得了比较好的效果。以三峡库区处于GLAS脚点的55个采样样点中树高提取的结果为例,最终在地形复杂地区各个坡度下区分植被类型时的树高值的均方根误差(RMSE)达到0.696,不区分植被类型时的树高值的均方根误差达到了0.737相较于背景方法一般能达到的0.59-0.69而言,结果比较令人满意。
[0026] 当然,实施本发明的任一产品或方法必不一定需要同时达到以上所述的所有优点。附图说明
[0027] 为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
[0028] 图1是本发明一个实施例的方法流程图
[0029] 图2是本发明一个实施例中的高斯回波数据处理结果示意图;
[0030] 图3是本发明一个实施例中三峡库区处于GLAS脚点的55个采样样点中树高提取的结果,及与实际测量值的比较。

具体实施方式

[0031] 下面将结合本发明的附图,对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0032] 实施例1:本发明实施例提出了一种从激光雷达高斯回波数据中提取树高的方法,演示了以三峡地区的GLAS高斯回波数据的树高提取过程。参见图1,该方法包括:
[0033] 步骤101:选取与树高待测地点同属于一个地区类型的若干校准参考地点,并获取该若干校准参考地点处的实测树高值;
[0034] 实际操作中,作为参考数据而言的所述实测树高值H可以通过很多手段获取,例如本发明实施例中采取了实地测量的方法,将三峡库区处于GLAS脚点的若干个采样样点作为校准参考地点,并在这些地点进行了树高测量。且本发明实施例中,所述“同一地区类型”具体指的是同一坡度类型,或者同一植被类型,或者同一坡度类型且同一植被类型。
[0035] 步骤102:获取该若干地点的激光雷达高斯回波数据和数字高程模型下与该若干地点一一对应的地形指数g;
[0036] 对GLAS数据进行处理得到高斯回波数据的方法有很多种,本实施例采取如下的方法:先通过高斯滤波器低通滤波,再基于直方图估计背景噪声,查找数据中拐点并利用拐点筛选后找到若干高斯波峰,然后剔除无效信号后按面积合并为六个高斯波峰,最后用滤波后数据对其进行Lenvenberg-Marquardt非线性拟合,得到最终的高斯回波数据。处理结果如图2所示(三峡地区的不同地形与植被类型的四组激光雷达高斯回波数据(实线)和原始平滑滤波函数(虚线),在说明书附图中分别表示为(a)植被:平缓地、(b)无植被:台地、(c)无植被:平地、(d)植被:陡坡)。
[0037] 地形指数g指的是对应某地点一定大小采样窗口(3×3,5×5与7×7)下数字地面模型(DEM)中地表高程的平均大小。其值代表了该位置处地面的垂直高度。在本发明实施例中采用5×5采样窗口下的地形指数数据。其中,数字地面模型(DEM)为公知的辅助数据库,对应地面的经纬度坐标可以读出该地点的地面高程。
[0038] 步骤103:从所述激光雷达高斯回波数据中提取得到波形质心到回波结束位置的距离wc和边缘长度l;
[0039] 与背景技术描述的相同,高斯回波数据用数学形式表达为数个高斯回波函数的叠加加上一个估计偏差,每个高斯回波数据分量的波峰位置(横坐标)对应一个高度信息。从而高斯回波数据如图2所示的一样为以时间为横坐标,信号强度为纵坐标的,按时间顺序排列的一系列高斯回波。其中最靠左的(对应时间最小)回波为垂直高度最高的,最靠右的则为最低的。从而按照常理而言,一般情况下最左边的高斯回波即为植被回波,最右边的高斯回波即为地面回波,回波位置指的就是回波波峰所在的横坐标位置,回波初始位置与结束位置就是回波数据中最左和最右边的点。但实际操作中为保证精度使用数学方法精确判断这些数据。
[0040] 波形质心的位置为植被回波波形面积一半处所在的横坐标位置,而回波结束位置由以下方法确定:设定高斯回波结束的阈值为结束噪声的平均值加上其标准偏差的和。从最后一个回波位置开始向后搜索,如果波形数据的值小于高斯回波结束的阈值,则确定此处的横坐标为回波结束位置。
[0041] 边缘长度l为地面回波的半宽,计算方法为从回波结束位置开始向前查找一个回波位置,使其与回波结束位置之差大于激光脉冲的半宽,则该回波位置即为地面回波的位置,边缘长度大小就等于该回波位置与回波结束位置之差,即为地面回波半宽。
[0042] H、g、wc和l四个参数的单位原本均为m,但其中wc和l需要参照GLAS产品说明进行由时间到高度的同比转换。但本方法中由于存在参数a0、a1和a2,故所有的同比转换均无必要进行,只需在采集数据时保持标度一致,其余的所有同比转换过程都可以隐含在参数a0、a1和a2中。
[0043] 步骤104:以H=a0×(wc-a1×g+a2×l)为拟合函数对该对应若干地点的若干组H、g、wc和l进行函数拟合,得到一组最佳拟合参数a0、a1和a2;
[0044] 本方法实施例采用的拟合方法为Lenvenberg-Marquardt非线性拟合方法。每组拟合参数都对应于一种地区类型,不同地区类型间的拟合参数会有一定差异,为保证树高提取精度须尽可能地采用同类型地区中的地点作校准参考地点。
[0045] 步骤105:对于该树高待测地点,先获取其激光雷达高斯回波数据,并从中提取波形质心到回波结束位置的距离w’c和边缘长度l’;
[0046] 该步骤中的提取方法与上面的步骤中的完全相同。
[0047] 步骤106:结合数字高程模型下的该地点的地形指数g’和所述最佳拟合参数a0、a1和a2由公式H’=a0×(w’c-a1×g’+a2×l’)计算得出该地点的测量树高值H’。
[0048] 与上面的步骤相同,地形指数g’仍然是数字地面模型(DEM)中地表高程的平均大小,对应于树高待测地点的经纬度坐标,采样窗口大小同样为5×5。公式中的拟合参数已经在函数拟合过程中给出,结合树高待测地点的激光雷达高斯回波数据和对应的地形指数就可以得到该地点树高的测量值。
[0049] 总的来看,本发明实施例提出的方法大致上可以分为参数拟合(步骤101至104)和数据代入(步骤105至106)两个步骤。而参数拟合仅仅是根据校准参考地形的采样结果对拟合函数中的拟合参数进行获取或者校准,也就是说对于树高提取而言只需要一组较为准确的参数即可,而并不是每一次树高的提取都需要进行一次校准。对本发明实施例而言,树高的提取过程实际上只包括步骤105至106,只是提取过程中用到的参数由步骤101至104进行获取或校准。
[0050] 在本发明实施例步骤103中所说的地形指数采用5×5采样窗口下的地形指数数据,因为在本发明实施例的实验数据表明森林地区中5×5采样窗口下的地形指数与实测地面高度差异的相关性相比3×3与7×7的更高,参见表1。表1中:0°、45°、90°、135°和全窗口代表同一大小采样窗口下的五种采样模式。采样点的数量为针叶林14个、阔叶林12个、混交林4个。总体来看,除阔叶林的采样结果中7×7采样窗口下的相关度较高(但仍未达到0.6,故本发明实施例中阔叶林部分数据不作重点),5×5采样窗口下的相关度都比较高,故选取该大小的采样窗口。
[0051] 表1不同采样窗口下的地形指数与实测地面高度差异的相关性
[0052]
[0053] 在步骤101至104的处理之后,在本发明实施例的数据中示意的拟合结果参见表2
2,其中R 为拟合度,代表采样点对拟合结果的贴合程度。
[0054] 表2示意的参数拟合结果
[0055]
[0056] 方法中采用植被回波的波形质心而不采用回波初始位置的原因在于,有研究指出,这一处理使提取的树高值对坡度的敏感性较低,对坡度较大地区结果相对比较稳定。而且在缓坡或粗糙度较大的情况下,地表回波半宽增大,回波结束位置以一定比例向后移动,会对树高的提取造成影响。从而引入边缘长度l,在本发明实施例的数据中引入前后的拟合效果如表3所示。可见,改进后对于同一植被类型的拟合度明显更好一些。
[0057] 表3引入边缘长度l的拟合效果示意
[0058]
[0059] 对应的树高提取结果参见图3,图3是三峡库区处于GLAS脚点的55个采样样点中树高提取的结果,同时将GLAS激光雷达回波数据中的树高提取结果与地面实测数据进行对比。可见,在三峡库区GLAS对森林冠层高度的测量结果除个别有偏差外,大部分比较令人满意。
[0060] 对于三峡库区处于GLAS脚点的55个采样样点中树高提取的结果,按坡度等级(坡度)划分的误差统计如表4所示。本发明实施例对小于20度坡度等级下的利用30个GLAS激光雷达高斯回波数据理直接测量所得树高提取值与地面实测数据进行对比,平均误差结果如表5-8。可见坡度小于5度时的精度最高;坡度大于5度时,由于地形的影响,造成林分冠层高度测量精度下降。同时,随着坡度增大,其测量误差也逐渐增大。
[0061] 表4不同植被类型各坡度等级树高提取值误差统计(单位:m)
[0062]
[0063] 总的来说,本发明实施例实现了较大坡度地区激光雷达高斯回波数据中树高的提取,且均方根误差达到了0.696,如前所述,原因主要是质心对坡度的敏感性较低。而通过对不区分森林类型,该方法的均方根误差达到了0.737,相较于背景方法一般能达到的0.59-0.69而言,结果比较令人满意。
[0064] 需要说明的是,在本文中,诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来,而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下,由语句“包括一个……”限定的要素,并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
[0065] 以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
QQ群二维码
意见反馈