基于多源卫星遥感数据的森林高度反演方法

本发明属于卫星遥感图像处理与应用的
技术领域：
。
背景技术：
：森林是陆地系统的主体，在全球水文、生态、碳循环及气候变化中起着重要作用。树高是林层划分的依据，是森林调查中重要因子之一，不仅用于反映林地生产力，而且用于确定立木材积和材积生长率。传统的森林资源调查以地面调查为主，工作量大、周期长，且难以实现大区域全面调查，而随着遥感技术应用日益深入和成熟，其覆盖范围大、重复周期短、成本低等独特的优势，为在大尺度上进行森林结构信息探测提供了可能。传统的光学遥感数据源由于受到云、雨、雾等影响，无法实现区域及全球范围内连续无缝的森林参数提取；另外，普通的光学传感器难以提供森林垂直分布信息。合成孔径雷达(sar)可实现对观测对象全天时、全天候的连续观测，其多频率、多极化与多入射角的特征为大面积的森林结构参数反演带来了应用前景。riom与hoffer等利用机载与星载l波段hh极化sar,研究表明树高、树龄等与雷达后向散射具有正相关关系(参见therelationbetweentypesofmaritimepineforestandl-bandhhpolarizationradarbackscatter)。基于sar后向散射系数估算森林参数的主要局限在于一定生物量水平上的信号饱和。干涉sar(insar)具有对森林垂直结构敏感的特性，能够有效地补充目前遥感技术森林高度估测中的不足。极化干涉sar(polinsar)在干涉sar的基础上增加了极化信息，进一步拓展了干涉sar在森林高度估测中的应用范围。最先用于描述森林植被的干涉sar模型是“水云模型”。treuhaft等提出了随机方向性散射(rv)模型、包含地面散射成分的随机方向性(rvog)模型和方向性散射(ovog)模型(参见estimatingforestvericialstructurefrommultialtitude,fixed-baselineradarinterferometricandpolarimetricdata)。cloude等基于rvog和ovog建立了极化干涉sar森林树高反演的模型和方法，并进行了理论分析和应用。其提出的三阶段植被参数反演方法，是一种从物理角度出发进行求解的新的几何方法，其利用衰减系数补偿森林密度和结构的变化以获取较高的估测精度(参见robustparameterestimationusingdualbase-linepolarimetricsarinterferometry)。另外，cloude提出了极化相干层析(pct)方法，进一步拓展了polinsar用于植被结构参数信息提取的理论和方法(参见three-stageinversionprocessforpolarimetricsarinterferometry).激光雷达(lidar)是近十年来发展比较迅速的主动遥感技术，并在森林结构参数的测量与反演方面得到了成功应用。sun等利用icesatglas对森林冠层高度数据进行了评估(参见forestverticalstructruefromglas:anevaluationusinglvisandsrtmdata)。到目前为止，国内外研究学者提出了很多森林高度估测算法，但仍存在一些明显的缺陷：(1)算法复杂度高；(2)现有算法中由于冠层和地表相位中心的不确定性、森林微波散射模型的局限性使得其在区域和全球尺度大范围开展还需要进一步研究；(3)验证区域森林类型单一，地形平坦，对于复杂地形的影响多未考虑；(4)使用数据源单一，未联合多源遥感数据进行森林高度反演；(5)lidar数据虽可获得垂直结构信息，进而估测森林高度，但由于其获取费用高而使得其大面积应用受到很大局限。技术实现要素：为了解决polinsar技术用于森林高度提取中地相位估计不准确、算法复杂度高及激光雷达数据高成本的问题，本发明采用一种基于星载sar和光学影像联合反演森林高度方法，首先对sar影像进行极化处理得到整幅影像的后向散射系数、入射角，利用insar测量技术获取主副影像相位、相干性、数字高程模型(dem)，并从dem中提取坡度信息，然后基于多光谱影像计算生物物理变量和光谱植被指数，使用机器学习算法将野外测量、极化信息、干涉测量、生物物理变量、光谱植被指数和地形变量联系起来，有效提取了森林高度信息。本发明采用的技术方案具体步骤如下：,步骤一、影像预处理：1.sar影像极化处理：获取到的高分辨率合成孔径雷达遥感影像是单视复数据(slc)的sentinel-1原始数据，依次经过多视、滤波、辐射定标和地理编码得到vv后向散射系数、vh后向散射系数和入射角；2.sar影像干涉处理：同时对主副slc影像依次进行干涉图生成和干涉去平、自适应滤波和相干系数生成、相位解缠、轨道精炼和重去平，及相位转高程，提取出坡度数据；对经过自适应滤波和相干系数生成处理的主副slc影像进行地理编码得到相位和相干系数；3.多光谱数据预处理：对sentinel-2原始多光谱数据，进行辐射定标和大气校正获得level-2a级产品数据；4.数据配准和剪裁：对预处理后sar和多光谱影像进行重采、配准和剪裁，得到实验区域的多源遥感数据。步骤二、提取特征向量：基于步骤一中得到的多光谱影像计算生物物理变量和光谱植被指数，生物物理变量包括包括叶面积指数(lai)，植被覆盖率(fvc)，叶绿素叶片中的糖含量(cab)，冠层水含量(cwc)和吸收的光合有效辐射的比例(fapar)；光谱植被指数包括：比值植被指数(rvi)、增强型植被指数(evi)、归一化植被指数(ndvi1、ndvi2、ndvi3)。步骤三、基于梯度提升决策树的xgboost回归学习方法估测树高：将步骤一与步骤二中提取的后向散射系数、入射角、干涉相位、相干系数、坡度、生物物理变量及光谱植被指数构成多维特征向量作为预测变量，结合野外实测数据，送入xgboost回归器，训练得到森林高度反演模型；步骤四、将提取的后向散射系数、入射角、干涉相位、相干系数、坡度、生物物理变量及光谱植被指数输入模型得到最终的森林高度预测结果。步骤一:影像预处理中1.sar影像极化处理：利用完整的遥感图像处理平台(envi)中的雷达图像基本处理工具(sarscape)进行多视处理、极化滤波、辐射定标和地理编码。(a)多视处理：目的是为了抑制sar图像的斑点噪声。具体方法是，使用envi中的多视工具(multilooking)处理，处理得到的多视强度图像是距离向和/或方位向像元分辨率的平均值，提高了多视图像的辐射分辨率，降低了空间分辨率。(b)极化滤波：目的是去除雷达图像的斑点噪声抑制其对影像上地物信息的干扰。具体方法是利用滤波工具(filtering)，选取refined-lee滤波方式进行极化滤波处理，以抑制相干斑噪声对影像上地物信息的干扰。(c)辐射定标和地理编码的方法：利用地理编码和辐射定标工具对雷达图像进行辐射定标，并为图像中的每个像元附上相应的经纬度地理信息，同时生成vv后向散射系数图、vh后向散射系数图和入射角图；2.sar影像干涉处理：利用完整的遥感图像处理平台(envi)中的雷达图像基本处理工具(sarscape)。(a)获取参考数字高程模型(dem)，根据原始sar数据的矢量文件下载数字高程模型。(b)干涉图生成和干涉去平：利用interferogramgeneration工具，输入两幅影像，一个作为主影像，一个作为副影像，同时输入参考dem文件，生成干涉图和去平后的干涉图。(c)自适应滤波和相干系数生成：利用adaptivefilterandcoherencegeneration工具，生成滤波后的干涉图和相干系数图。(d)相位解缠：利用phaseunwrapping工具，生成解缠后的相位图。优选地，选择最小费用流算法“minimumcostflow”，设置相干性阈值为0.2，分解等级为1。(e)轨道精炼和重去平：利用refinementandre-flattening工具，在步骤c)生成的相干系数图的平滑区域上选择控制点，避免选择地形残差条纹区域，进一步消除可能存在的平地效应，纠正相位偏移。(f)相位转高程以及地理编码：利用phasetoheightconversionandgeocoding工具，将步骤(e)得到的相位图转换为高程数据，根据地理编码到制图坐标系统，完成相位图转换为高程图的过程，生成高精度的dem。(g)地理编码：使用geocoding工具完成对步骤(c)获得的干涉图和相干系数图进行地理编码。(h)提取坡度：利用遥感图像处理平台(arcmap10.2)中的surfaceslope工具，完成从dem数据提取坡度信息的过程。3.多光谱数据预处理：在欧空局(esa)网站上下载sentinel-2level-1c的原始多光谱数据，并利用遥感图像处理平台(snap)中的sen2cor插件完成对原始多光谱数据的辐射定标和大气校正，并重采样到同一分辨率。4.数据配准和剪裁：将sar数据、重采样后的多光谱数据以及高精度的dem数据进行配准和剪裁，得到实验区域的多源遥感数据。步骤二：提取特征向量步骤中(a)提取生物物理变量：生物物理变量的独特信息，有助于增强森林高度的预测能力。利用遥感图像处理平台(snap6.0)中的biophysicalprocessor工具，得到多光谱数据的生物物理变量，包括：叶面积指数(lai)、植被覆盖度(fvc)、叶片叶绿素含量(cab)、林冠层含水量(cwc)和吸收的光合有效辐射的比例(fapar)。(b)提取光谱植被指数：根据sentine-2数据中的多个波段的反射值计算出包括比值植被指数(rvi)、增强型植被指数(evi)、归一化植被指数(ndvi1、ndvi2和ndvi3)在内的5个植被指数，具体计算公式如下：rvi＝nir/r(1)evi＝2.5*((nir-r)/(nir+6*r-7.5*b+1))(2)ndvi1＝(nir-r)/(nir+r)(3)ndvi2＝(nir2-re3)/(nir+re3)(4)ndvi3＝(re2-re1)/(re2+re1)(5)其中，r为红光波段(b4)的反射值，g为绿光波段(b3)的反射值，b为蓝光波段(b2)的反射值，nir为近红外波段(b5)的反射值，nir2为近红波段(b8)的反射值，re1为红边波段(b5)的反射值，re2为红边波段(b6)的反射值，re3为红边波段(b7)的反射值。步骤三：基于梯度提升决策树的xgboost回归学习方法估测树高步骤中将上述步骤获得的两个后向散射系数(vv、vh)、干涉相位、相干系数、入射角、坡度、五个生物物理变量(lai、fvc、cab、cwc、fapar)以及五个光谱植被指数(rvi、evi、ndvi1、ndvi2、ndvi3)，一共16个特征作为模型的输入，送入python机器学习库中的xgboost回归器，并设定参数为max_depth＝5，learning_rate＝0.1，n_estimators＝100，silent＝false，objective＝'reg:gamma'。将实测数据中随机选取4/5用于模型参数的训练，剩下的数据用决定系数、均方根误差来对比分析算法对精度进行验证，得到森林高度反演模型。本发明的有益效果：本发明根据多源遥感信息可以快速、准确的提取东北地区森林高度，解决了现有森林高度估测算法中地相位估计不准确、算法复杂度高及激光雷达数据高成本的问题。此外，本发明使用的数据源易于获取，利用polsar和insar技术，提取对森林高度敏感的极化信息和干涉信息，结合坡度和多光谱信息，为森林高度的计算提供大量有效的数据。本发明为森林高度估测提供有效的算法，为森林生物量、森林管理及碳循环提供一定的技术支持。附图说明图1是本发明基于多源遥感数据的东北地区森林高度反演模型流程图。图2是本发明实施例1使用的研究区。图3是本发明实施例1的后向散射系数图和入射角图。图4是本发明实施例1的干涉相位图和相干系数图。图5是本发明实施例1的坡度图。图6是本发明实施例1的生物物理变量图。图7是本发明实施例1的光谱植被指数图。图8是本发明实施例1的经过xgboost回归之后真实值和预测值的关系图。具体实施方式实施例1：本发明联合多源遥感数据，其中sar数据采用哨兵1号(sentinel-1)卫星数据，其空间分辨率为20m，如表1所示。光学数据采用2019年2月25日和3月17日哨兵二号(sentinel-2)卫星数据,其波段信息如表2所示。实验区域为吉林省长春市净月潭国家级风景区(图2)，森林类型主要包括针叶林、阔叶林、针阔混交林等，利用sar和光学影像提取的特征构建多维特征变量，结合野外测量数据，采用xgboost算法完成森林高度反演。表1卫星获取日期产品类型极化方式产品idsentinel-1b2019.02.18level-1slc-sdvvv、vh272fsentinel-1b2019.03.02level-1slc-sdvvv、vh7a83sentinel-1b2019.03.14level-1slc-sdvvv、vh2d24sentinel-1b2019.03.26level-1slc-sdvvv、vh239b表2sentinel-2波段中心波长(um)分辨率(m)b1-海岸/气溶胶波段0.44360b2-蓝波段0.4910b3-绿波段0.5610b4-红波段0.66510b5-红边波段0.70520b6-红边波段0.7420b7-红边波段0.78320b8-近红外波段0.84210b8a-近红波段0.86520b9-水蒸气波段0.94560b10-短波红外波段1.37560b11-短波红外波段1.6120b12-短波红外波段2.1920步骤一:影像预处理1.sar影像极化处理：利用完整的遥感图像处理平台(envi)中的雷达图像基本处理工具(sarscape)，依次进行如下处理；(a)多视处理：利用多视处理工具(multilooking)从sentinel-1sar影像数据中读取距离向和方位向分辨率，并计算出距离向视数和方位向视数。(b)滤波：利用滤波工具(filtering)，选取refinedlee滤波方式，窗口大小选为5x5，可以较好的抑制相干斑噪声对影像上地物信息的干扰。(c)辐射定标和地理编码：利用滤波工具(geocodingandradiometriccalibration)能自动从滤波后的数据文件中读取参数，从而完成辐射定标。地理编码选取由envi软件自动下载的30m分辨率的srtm的数字高程模型(dem)作为参考标准，为每个像元附上相应的经纬度地理信息，同时生成vv后向散射系数图、vh后向散射系数图和入射角图，如图3(a)、3(b)、3(c)。2.sar影像干涉处理：利用完整的遥感图像处理平台(envi)中的雷达图像基本处理工具(sarscape)，依次进行如下处理；(a)下载参考dem文件：打开工具/sarscape/generaltools/digitalelevationmodelextraction/srtm-3version2，在inputfile输入两个时相sentinel-1的矢量文件(slc_list文件)，在dem/cartographicsystem中设置geogloblewgs84，其他参数按照默认，下载srtm-3version2dem。(b)干涉图生成和干涉去平：打开/sarscape/interferometry/phaseprocessing/interferogramgeneration工具，在input面板，输入两个时相sentinel-1的矢量文件(slc_list文件)和srtm-3version2dem文件，距离向视数和方位向视数采用自动添加，制图分辨率按照默认的20。生成干涉图和去平后的干涉图。(c)自适应滤波和相干系数生成：打开/sarscape/interferometry/phaseprocessing/adaptivefilterandcoherencegeneration工具，采用goldstein滤波，这种滤波方法的滤波器是可变的，提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声，这种方法是最常用的方法，生成滤波后的干涉图和相干系数图。(d)相位解缠：干涉相位只能以2π为模，所以只要相位变化超过了2π，就会重新开始和循环。相位解缠是对去平和滤波后的位相进行相位解缠，解决2π模糊的问题。打开/sarscape/interferometry/phaseprocessing/phaseunwrapping工具，选择最小费用流算法“minimumcostflow”，设置相干性阈值为0.2，分解等级为1，生成解缠后的相位图。(e)轨道精炼和重去平：打开/sarscape/interferometry/phaseprocessing/refinementandre-flattening工具，对滤波后的相干系数图，在控制点生成面板上，点击next，打开控制点选择工具，鼠标变为选点状态，在图像上平滑的地方单击鼠标左键，选择控制点，进行轨道精炼和相位偏移的计算，消除可能的斜坡相位，对卫星轨道和相位偏移进行纠正，生成轨道精炼和重去平后的相位图。(f)相位转高程：打开/sarscape/interferometry/phaseprocessing/phasetoheightconversionandgeocoding工具，将步骤(e)得到的相位图转换为高程数据，根据地理编码到制图坐标系统，完成相位图转换为高程图的过程，生成高精度的dem。(g)地理编码：使用geocoding工具对步骤(c)完成对干涉图和相干系数图的地理编码，如图4(a)、(b)所示。(h)提取坡度：利用遥感图像处理平台(arcmap10.2)中的arctoolbox/3danalyst/surfacetriangulation/surfaceslope工具，输入insar流程得到的dem图像，完成从dem数据提取坡度信息的过程(图5)。3.多光谱数据预处理：利用欧空局(esa)发布的sen2cor插件完成对sentinel-2level-1c原始多光谱数据的辐射定标和大气校正，得到level-2a产品，并重采样到同一分辨率20m。4.数据配准和剪裁：将sar数据、重采后的多光谱数据以及dem数据进行配准和剪裁，得到实验区域的多源遥感数据。步骤二：提取特征向量(a)提取生物物理变量：利用遥感图像处理平台(snap6.0)中的optical/thematiclandprocessing/biophysicalprocessor工具，对重采样到20m的level-2asentine-2数据进行辐射传输模型和神经网络算法计算，得到多光谱数据的生物物理变量(lai、fvc、cab、cwc、fapar)，如图6(a)、(b)、(c)、(d)、(e)所示。(b)提取光谱植被指数：植被指数是对地表植被状况的简单、有效和经验的度量,目前已经定义了40多种植被指数,广泛地应用在全球与区域土地覆盖、植被分类和环境变化。一共提取了包括rvi、evi、ndvi1、ndvi2、ndvi3在内的5个植被指数，如图7(a)、(b)、(c)、(d)、(e)所示。步骤三：基于梯度提升决策树的xgboost回归学习方法估测树高xgboost算法是基于梯度提升决策树的一种集成学习模型，该算法中的决策树具有先后关联，当前预测以上一轮的预测误差为基础，利用各轮预测误差迭代构建模型，提升预测的准确性。2019年3月上旬获得的野外测量数据包含一些林分特征，例如平均高度，胸径(dhb)，经度，纬度，树冠郁闭度，树木种类，树木数量和树冠宽度。实地测量了30个样地，样地大小为20×20m^2，由于在同一个实测点有四景sar影像，因此实测点扩充到120个。此外，森林高度范围：9–27m。在本实验中没有获得涵盖整个研究区域的完整参考数据。但是，现场数据源的收集提供了森林状况的信息。从步骤一和步骤二中提取了两个后向散射系数(vv、vh)、干涉相位、相干系数、入射角、坡度、五个生物物理变量(lai、fvc、cab、cwc、fapar)以及五个光谱植被指数(rvi、evi、ndvi1、ndvi2、ndvi3)，一共16个特征作为模型的输入，送入python机器学习库中的xgboost回归器，其中参数选取为(max_depth＝5，learning_rate＝0.1，n_estimators＝100，silent＝false，objective＝'reg:gamma')。从120组实测数据中随机选取4/5用于模型参数的训练，剩下的24组数据用于精度验证。实验结果如图8所示，均方根误差(rmes)＝2.0352,决定系数(r2)＝0.6840。光学遥感数据可获取大范围的光谱信息，但是其可见光及红外波段仅与叶水平生物量产生作用；而微波具有穿透树冠的能力，不仅能与树叶发生作用，还能与树干、树枝发生作用，获取森林内部结构信息。本发明联合多源遥感数据，分别从水平结构和垂直结构获取森林信息，提升了森林高度预估的准确性，可以实现对森林高度较为稳定精确的反演效果。当前第1页12