用于监测山地的遥感影像近红外波段地形阴影校正方法

技术领域

本发明涉及地形阴影补偿技术领域，特别是一种用于监测山地的遥感影像近红外波段地形阴影校正方法。

背景技术

地形阴影严重干扰山区地物监测、分类和提取的精度，给山区土地覆被解译造成巨大困难，因此研究地形阴影光谱信息补偿方法具有极其重要的意义。目前的地形阴影校正方法主要有基于DEM数据的经验校正模型、物理校正模型以及波段组合优化计算模型(植被指数)等。基于DEM数据的经验校正方法主要考虑太阳直射辐射，以反射率与cos i的统计关系为基础进行经验校正，包括统计-经验模型、归一化模型、C校正、SCS校正等，但该类方法阳坡校正效果好，而地形阴影校正效果有限，尤其落影校正效果不佳。物理校正模型同时考虑太阳直射辐射、大气散射辐射以及周围地形反射辐射，效果较好，但模型参数较多且计算过程复杂，推广难度大。基于波段组合优化计算模型的方法通过构建特殊的植被指数来获取消除地形影响后的植被信息，如波段比模型、VBSI指数、地形调节植被指数(TAVI)、归一化差值山地植被指数、改进的增强植被指数、阴影消除植被指数(SEVI)等，然而该类模型只能获取单一波段的指数信息，缺乏多光谱信息。近年提出的“一种遥感影像可见光波段地形落影校正方法(专利号2021101469628)”可以校正红-绿-蓝可见光地形阴影，而不能校正近红外波段地形阴影。近红外波段光谱信息是植被遥感的重要数据，对山区地物监测、植被分类和提取至关重要。

发明内容

有鉴于此，本发明的目的在于提供一种用于山区地物监测的遥感影像近红外波段地形阴影校正方法，获取能有效消除地形阴影干扰的遥感影像近红外波段反射率信息。

为实现上述目的，本发明采用如下技术方案：用于监测山地的遥感影像近红外波段地形阴影校正方法，包括以下步骤：

步骤S1：获取研究区的卫星遥感数据和数字高程模型DEM数据，并做数据预处理；

步骤S2：采用机器学习和水体指数即NDWI提取地形阴影；

步骤S3：采用大气校正后的地表反射率进行SCS+C地形校正；

步骤S4：采用地表反射率计算地形阴影消除植被指数即SEVI；

步骤S5：进行红光波段地形阴影校正；

步骤S6：利用阴影区红光波段反射率采用以下公式计算阴影区近红外波段反射率：

ρ

式中，ρ

步骤S7：校正效果评价，采用目视分析、光谱统计分析、相对误差分析等进行地形校正效果评价，光谱统计分析具体包括反射率与cosi散点分析、玫瑰图。

在一较佳的实施例中，所述步骤S1具体为：

步骤S11：收集研究区的卫星遥感数据和DEM数据；

步骤S12：对卫星遥感数据进行辐射定标和大气校正，得到研究区地表反射率数据；

步骤S13：对DEM数据进行投影变换、重采样及裁剪操作，得到研究区DEM数据，最后将获取数据规范为同一数据格式。

在一较佳的实施例中，所述步骤S2具体为：

步骤S21：在研究区分别选取阴影区、光照区及水体三类样本点；

步骤S22：采用随机森林即RF分类方法结合三类样本点对研究区进行分类，提取总体阴影，并利用水体指数NDWI进一步分离水体的干扰。

在一较佳的实施例中，所述步骤S3具体为：

步骤S31：结合研究区DEM数据，遥感影像太阳高度角、太阳方位角数据，计算太阳入射角余弦值cos i，计算公式如下：

cos i＝cosσcosθ+sinσsinθcos(β-ω) (2)

式中，i表示太阳入射角；σ表示地形坡度角；θ表示太阳天顶角；β表示地形坡向角；ω表示太阳方位角；

步骤S32：结合cos i计算矫正参数c，计算公式如下：

L

c＝a/b (4)

式中，L

步骤S33：对红光波段地表反射率数据进行SCS+C校正，计算公式如下：

式中，L

在一较佳的实施例中，所述步骤S4具体为：

根据红光波段和近红外波段地表反射率计算研究区SEVI，其公式如下：

式中，ρ

其中，f(Δ)计算步骤如下：利用DEM计算坡度图，将坡度图重采样为6×6平方公里空间分辨率的块，将坡度最大的前1％的块选为陡坡块集；按照每个陡坡块的空间范围采用遥感影像反射率数据逐块计算SEVI及其信息熵H，公式为：

式中，p

令调节因子f(Δ)从0开始，以预设步长0.001依次递增，计算陡坡块SEVI信息熵，当H取值最大时，取与最大H对应的调节因子f(Δ)作为该陡坡块的调节因子优化解f

f

比较不同陡坡块的信息熵，选择信息熵最大的纯植被陡坡块对应的调节因子作为整景影像的最优调节因子f

f

式中，m为陡坡块的数量。

在一较佳的实施例中，所述步骤S5具体为：

步骤S51：以整个研究区为处理范围，采用ArcGIS随机生成点工具生成一定数目的随机点，以光照区图层为掩膜提取随机点；

步骤S52：采用ArcGIS中的多值提取至点工具将光照区随机样本点对应的SEVI及SCS+C校正后红光波段反射率记录到该点的属性表中，结合光照区随机样本点属性表制作样本集，并将样本集按预设比例分为训练样本集和测试样本集；

步骤S53：采用训练样本集的SEVI作为自变量，以训练样本集SCS+C校正后红光波段地表反射率为因变量构建RF回归模型；

步骤S54：调节超参数使RF回归模型最优；

步骤S55：以地形阴影区的SEVI为自变量，采用调参优化后的RF回归模型预测地形阴影区红光波段反射率。

与现有技术相比，本发明具有以下有益效果：

本发明能有效消除地形阴影，尤其是落影对遥感影像近红外波段的干扰，弥补了SEVI只能获取消除地形阴影干扰的单波段信息的不足，也弥补了常规基于DEM的地形校正方法在阴影区效果不佳尤其是落影失效的不足。

附图说明

图1为本发明优选实施例的方法流程示意图；

图2为本发明优选实施例的近红外波段地形阴影校正前后效果对比图；

图3为本发明优选实施例的近红外波段反射率光谱特征分析图；

图4为本发明优选实施例的近红外波段地形阴影(本影、落影)反射率相对阳坡反射率误差柱状图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

应该指出，以下详细说明都是示例性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式；如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

请参照图1-4，本发明提供一种用于监测山地的遥感影像近红外波段地形阴影校正方法，包括以下步骤：

步骤S1：获取数据并预处理，所需数据包括：卫星遥感数据如Landsat 8OLI、DEM数据如ASTER GDEM V2。对卫星遥感数据进行辐射定标和大气校正，得到地表反射率数据；对GDEM V2数据进行投影变换、重采样及裁剪等操作，得到DEM数据，最后将获取数据规范为同一数据格式；

步骤S2：地形阴影提取，首先采用RF分类提取总体阴影，然后计算NDWI去除阴影中的水体，得到地形阴影提取结果；

步骤S3：对地表反射率进行SCS+C校正；

步骤S4：采用地表反射率进行SEVI计算；

步骤S5：进行红光波段地形阴影校正；

步骤S6：根据得到的阴影区红光波段反射率采用以下公式计算阴影区近红外波段反射率：

ρ

式中，ρ

步骤S7：校正效果评价，采用目视分析、光谱统计分析(反射率与cosi散点分析、玫瑰图)、相对误差分析等进行地形校正效果评价。

在本实施例中，所述步骤S1具体为：

步骤S11：获取研究区的卫星遥感数据和DEM数据；

步骤S12：对卫星遥感数据进行辐射定标和大气校正，得到研究区地表反射率数据；

步骤S13：对DEM数据进行投影变换、重采样及裁剪操作，得到研究区DEM数据，最后将获取数据规范为同一数据格式。

在本实施例中，所述步骤S2具体为：

步骤S21：在研究区分别选取阴影区、光照区及水体三类样本点；

步骤S22：采用RF分类方法结合三类样本点对研究区进行分类，提取总体阴影，并利用NDWI进一步分离水体的干扰；

在本实施例中，所述步骤S3具体为：

步骤S31：结合研究区DEM数据，遥感影像太阳高度角、太阳方位角数据，计算太阳入射角余弦值cos i，计算公式如下：

cos i＝cosσcosθ+sinσsinθcos(β-ω) (2)

式中，i表示太阳入射角；σ表示地形坡度角；θ表示太阳天顶角；β表示地形坡向角；ω表示太阳方位角；

步骤S32：结合cos i计算调节参数c，计算公式如下：

L

c＝a/b (4)

式中，L

步骤S33：对红光波段地表反射率数据进行SCS+C校正，计算公式如下：

式中，L

在本实施例中，步骤S4具体包括以下步骤：

根据红光波段和近红外波段地表反射率计算研究区SEVI，其公式如下：

式中，ρ

其中，f(Δ)计算步骤如下：利用DEM计算坡度图，将坡度图重采样为6×6平方公里空间分辨率的块，将坡度最大的前1％的块选为陡坡块集；按照每个陡坡块的空间范围采用遥感影像反射率数据逐块计算SEVI及其信息熵H，公式为：

式中，p

令调节因子f(Δ)从0开始，以预设步长0.001依次递增，计算陡坡块SEVI信息熵，当H取值最大时，取与最大H对应的调节因子f(Δ)作为该陡坡块的调节因子优化解f

f

比较不同陡坡块的信息熵，选择信息熵最大的纯植被陡坡块对应的调节因子作为整景影像的最优调节因子f

f

式中，m为陡坡块的数量。

在本实施例中，步骤S5具体包括以下步骤：

步骤S51：以整个研究区为处理范围，采用ArcGIS随机生成点工具生成一定数目的随机点，以光照区图层为掩膜提取随机点；

步骤S52：采用ArcGIS中的多值提取至点工具将与光照区随机样本点对应的SEVI及SCS+C校正后红光波段反射率记录到该点的属性表中，结合光照区随机样本点属性表制作样本集，并将样本集按预设比例分为训练样本集和测试样本集。

步骤S53：调用python中集成的Random Forest Regressor模块，以训练样本集的SEVI作为自变量，以训练样本集SCS+C校正后红光波段地表反射率为因变量构建RF回归模型；

步骤S54：调用python中集成的Bayesian Optimization模块调节RF回归模型中的超参数，使RF回归模型最优。

步骤S55：以地形阴影区的SEVI为自变量，采用调参优化后的RF回归模型预测得到地形阴影区红光波段反射率；

在本实施例中，步骤S7具体包括以下步骤：

步骤S71：图2目视分析可见，未做近红外等波段校正前地形影响显著，地形阴影遍布；采用本发明方法校正的近红外波段，结合“一种遥感影像可见光波段地形落影校正方法(专利号2021101469628)”校正的红、绿波段生成的近红-红-绿假彩色合成影像上地形阴影消失，地形影响消除，影像呈平面分布。

步骤S72：采用近红外波段反射率与cosi散点图分析表明：地形校正前，相关性强，r

步骤S73：玫瑰图分析表明未经地形校正的近红外波段本影、落影的数值小，位于图中内核；而阳坡的数值大，分布成半壳在外。SCS+C校正后，近红外波段反射率本影、落影的数值变大，形成从内核向外的放射状；阳坡数值分布形成半壳状，表明SCS+C校正对地形阴影校正有效果。本发明方法校正后，近红外波段反射率本影、落影与阳坡的数值接近，形成一个闭圏，均匀分布于0～360度坡向，表明本发明方法对近红外波段地形阴影校正效果突出。

步骤S74：采用本影、落影相对阳坡的近红外波段反射率误差评价近红外波段地形阴影校正效果，其公式如下：

校正结果(图4)表明未经地形校正的近红外波段地表反射率的本影、落影与阳坡的相对误差较大，分别为84.01％和86.34％。SCS+C校正后，近红外波段反射率的本影、落影与阳坡的相对误差分别降为16.38％和43.26％，表明SCS+C校正对本影校正效果较好，而对落影校正效果欠佳。本发明方法校正后，近红外波段反射率的本影、落影与阳坡的相对误差分别降为0.90％和0.83％，表明本发明方法对近红外波段地形阴影校正效果突出。