CN114596234B - 一种复杂山地的ndvi地形阴影效应校正方法 - Google Patents

Abstract

本发明公开了一种NDVI地形阴影效应校正方法，以解决受太阳高度角和方位角影响而导致的复杂山地阴影区域NDVI较光照区域偏低的问题。该方法以光照区域和阴影区域的太阳辐射差异为基础，引入反映阴影强度的可变系数(α)来模拟阴影区域所接受不到的太阳直射光反射率，阴影区域的模拟太阳直射光反射率作为补偿量与原始反射率一起用于计算校正后的NDVI。本发明公开的方法能够避免本影区域的过度校正，同时，改善落影区域的校正不足，校正后的NDVI不会存在异常值，以一种简单而稳健的方式提高复杂山地NDVI的校正精度。

Description

一种复杂山地的NDVI地形阴影效应校正方法

技术领域

本发明属于遥感影像地形校正技术领域，具体地说，涉及一种复杂山地的NDVI地形阴影效应校正方法。

背景技术

NDVI是多种植被指数中应用最多最广泛的一种，相对于其他植被指数，NDVI的优点在于增强了对植被的响应能力，同时可以消除一部分地形及阴影的影响。然而，这些影响的消除是不完全的，特别是地形起伏较大的复杂山地区域。

山地区域的遥感影像因地形引起目标辐射特性差异，面向太阳的表面通常比背离太阳的表面接收到更多的辐射，即地形阴影效应。反映在遥感影像上则形成光照区和阴影区，其中，光照区由面向太阳的表面构成，阴影区则由背离太阳的表面以及山体的影子投射在邻近区域的表面构成。阴影的存在导致同一地物的光谱响应存在差异，进而导致同一区域内相同地物的NDVI值出现明显差异，影响地物解译、分类与反演的结果。因此，如何消除阴影对NDVI的影响这一问题显得格外重要。

目前，对于NDVI地形阴影效应的校正方法通常是基于地形校正模型的校正方法，此类方法需要先对各个波段的光谱信息进行改正，再计算NDVI。其缺点在于，影像获取时，如果太阳高度角较低，并且，研究区位于地势较为陡峭的山地中，则在本影区域极易出现过度校正的情况，同时极易产生异常值，而落影区域极易出现校正不足的情况。

综上，提供一种简单高效的NDVI地形阴影效应校正方法以反演更精确的NDVI是急需解决的一个问题。

发明内容

本发明要解决的技术问题在于克服现有技术的不足，针对现有方法存在的复杂山地阴影区域NDVI校正后失真和异常的问题，提供一种复杂山地的NDVI地形阴影效应校正方法，以光照区域和阴影区域的太阳辐射差异为基础，通过阴影区域地物固有的反射率特性来模拟阴影区域所接受不到的太阳直射光反射率，并将其作为补偿量与原始反射率一起用于计算校正后的NDVI，从而实现NDVI地形阴影效应的校正。

为解决上述技术问题，本发明采用技术方案的基本构思是：

一种复杂山地的NDVI地形阴影效应校正方法，包括以下步骤：

步骤1，获取复杂山地的多光谱影像数据和同区域的数字高程模型DEM数据。

步骤2，根据DEM数据计算天空观测因子和地形观测因子。

步骤3，通过影像成像时间获取太阳天顶角和方位角，并由6S模型计算太阳直射辐射、天空散射辐射、地形反射辐射。

步骤4，构建阴影指数，并通过阴影指数计算反映阴影强度的可变系数α。

步骤5，利用影像的红波段和近红外波段反射率并结合天空观测因子、地形观测因子、太阳直射辐射、天空散射辐射、地形反射辐射和可变系数α计算红波段和近红外波段的模拟太阳直射光反射率补偿量。

步骤6，构建NDVI校正公式，实现NDVI地形阴影效应校正。

步骤7，评估校正结果。

进一步地，步骤1中，多光谱影像数据为经过消除大气干扰的反射率数据，数字高程模型DEM数据替换为坡度和坡向数据。

进一步地，步骤2包括：通过步骤1中的DEM数据按照以下公式计算天空观测因子V_d和地形观测因子C_t：

V_d≈(1+cosθ_S)/2

C_t＝1-V_d

式中，θ_s是每个像元的坡度数据。

进一步地，步骤3为：根据影像拍摄区域的大气模式、气溶胶含量参数(或能见度)、目标高度参数、传感器高度参数、光谱参数和地表反射率类型通过6S大气校正模型计算红波段和近红外波段的太阳直射辐射E_d、天空散射辐射E_f、地形反射辐射E_a。

进一步地，步骤4中：根据步骤1得到的反射率数据，所述阴影指数SI按照公式SI＝(ρ_Coastal-ρ_Green)/(ρ_Coastal+ρ_Green+2×ρ_NIR)计算获得，求取SI的阈值c；

式中，ρ_Coastal、ρ_Green和ρ_NIR分别为海岸波段、绿波段和近红外波段的反射率；

当SI小于c时，所述α＝0；

当SI大于c时，α＝(SI-c)/(SI_max-c)

式中SI_max为SI的最大值。

进一步地，步骤5中，红波段和近红外波段的模拟太阳直射光反射率补偿量按照以下公式计算：

ρ_{dc_Red}＝α×ρ_Red×E_{d_Red}/(E_{f_Red}×V_d+E_{a_Red}×C_t)

ρ_{dc_NIR}＝α×ρ_NIR×E_{d_NIR}/(E_{f_NIR}×V_d+E_{a_NIR}×C_t)

式中，ρ_{dc_Red}和ρ_{dc_NIR}分别为红波段的模拟太阳直射光反射率补偿量和近红外波段的模拟太阳直射光反射率补偿量，ρ_Red和ρ_NIR分别为红波段和近红外波段反射率，E_{d_Red}和E_{d_NIR}分别为红波段和近红外波段的太阳直射辐射，E_{f_Red}和E_{f_NIR}分别为红波段和近红外波段的天空散射辐射，E_{a_Red}和E_{a_NIR}分别为红波段和近红外波段的地形反射辐射。

进一步地，步骤6中：第五步中模拟的红波段和近红外波段太阳直射光反射率补偿量与原始反射率一起用于计算校正后的NDVI，从而获得NDVI的校正公式：

NDVI_corr＝(ρ_NIR+ρ_{dc_NIR}-ρ_Red-ρ_{dc_Red})/(ρ_NIR+ρ_{dc_NIR}+ρ_Red+ρ_{dc_Red})

式中，NDVI_corr为地形阴影效应校正后的NDVI值。

进一步地，步骤7为：

校正结果的评估方式包括：校正前后的NDVI与验证NDVI的密度散点图对比；根据阴影类型(本影和落影)对比阴坡与阳坡的NDVI；校正前后的NDVI与验证NDVI的剖面线对比。评价标准为：密度散点图的动态范围关于1：1线高度对称，且不发生“外溢”现象；校正后本影和落影的NDVI接近阳坡的NDVI；校正后的NDVI剖面线接近验证NDVI剖面线。

采用上述技术方案后，本发明与现有技术相比具有以下有益效果。

本发明提供的复杂山地的NDVI地形阴影效应校正方法，以光照区域和阴影区域的太阳辐射差异为基础，引入反映阴影强度的可变系数(α)来模拟阴影区域所接受不到的太阳直射光反射率，阴影区域的模拟太阳直射光反射率作为补偿量与原始反射率一起用于计算校正后的NDVI。校正后，避免了本影区域NDVI过度校正的现象，校正后的NDVI不会存在异常值，改善了落影区域NDVI校正不足的现象。

下面结合附图对本发明的具体实施方式作进一步详细的描述。

附图说明

附图作为本申请的一部分，用来提供对本发明的进一步的理解，本发明的示意性实施例及其说明用于解释本发明，但不构成对本发明的不当限定。显然，下面描述中的附图仅仅是一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。在附图中：

图1为本发明具体流程示意图；

图2为本发明实施例多光谱影像数据和DEM数据示意图；

图3为本发明实施例阴影指数SI和可变系数α示意图；

图4为本发明实施例NDVI校正结果示意图；

图5为本发明实施例校正前后的NDVI与验证NDVI的密度散点图；

图6为本发明实施例子区域NDVI校正结果示意图；

图7为本发明实施例子区域NDVI剖面线示意图。

需要说明的是，这些附图和文字描述并不旨在以任何方式限制本发明的构思范围，而是通过参考特定实施例为本领域技术人员说明本发明的概念。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对实施例中的技术方案进行清楚、完整地描述，以下实施例用于说明本发明，但不用来限制本发明的范围。

在本发明的描述中，需要说明的是，术语“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

实施例一

本发明提供了一种复杂山地的NDVI地形阴影效应校正方法，如图1所示，本发明包括如下步骤：

步骤1，获取复杂山地的多光谱影像数据和同区域的数字高程模型DEM数据，如图2所示。

具体地，获取云南省高黎贡山地区的两期Landsat8 OLI多光谱影像，成像时间为2014年1月22日和2014年4月12日，两幅影像的太阳高度角分别为39.26°和62.60°。图2(a)为1月份的影像，用于NDVI地形阴影效应校正实验；图2(b)为4月份的影像，用于验证校正结果。两幅影像均做了辐射定标和黑暗像元法大气校正处理。图2(c)为DEM数据。

步骤2，根据DEM数据计算天空观测因子和地形观测因子。

通过DEM数据按照以下公式计算天空观测因子V_d和地形观测因子C_t：

V_d≈(1+cosθ_S)/2

C_t＝1-V_d

式中，θ_s是每个像元的坡度数据。

步骤3，通过影像成像时间获取太阳天顶角和方位角，并由6S模型计算太阳直射辐射、天空散射辐射、地形反射辐射。

具体地，根据影像拍摄区域的大气模式、气溶胶含量参数(或能见度)、目标高度参数、传感器高度参数、光谱参数和地表反射率类型通过6S大气校正模型计算红波段和近红外波段的太阳直射辐射E_d、天空散射辐射E_f、地形反射辐射E_a。

步骤4，构建阴影指数，并通过阴影指数计算反映阴影强度的可变系数α，如图3所示。

根据步骤1得到的反射率数据，图3(a)为所述阴影指数SI，按照公式SI＝(ρ_Coastal-ρ_Green)/(ρ_Coastal+ρ_Green+2×ρ_NIR)计算获得。

式中，ρ_Coastal、ρ_Green和ρ_NIR分别为海岸波段、绿波段和近红外波段的反射率。

图3(b)为SI结合其阈值c所计算的可变系数α，具体计算方式如下：

当SI小于c时，所述α＝0；

当SI大于c时，α＝(SI-c)/(SI_max-c)，

式中，SI_max为SI的最大值。

步骤5，利用影像的红波段和近红外波段反射率并结合天空观测因子、地形观测因子、太阳直射辐射、天空散射辐射、地形反射辐射和可变系数α计算红波段和近红外波段的模拟太阳直射光反射率补偿量。

具体地，红波段和近红外波段的模拟太阳直射光反射率补偿量按照以下公式计算：

ρ_{dc_Red}＝α×ρ_Red×E_{d_Red}/(E_{f_Red}×V_d+E_{a_Red}×C_t)

ρ_{dc_NIR}＝α×ρ_NIR×E_{d_NIR}/(E_{f_NIR}×V_d+E_{a_NIR}×C_t)

式中，ρ_{dc_Red}和ρ_{dc_NIR}分别为红波段的模拟太阳直射光反射率补偿量和近红外波段的模拟太阳直射光反射率补偿量，ρ_Red和ρ_NIR分别为红波段和近红外波段反射率，E_{d_Red}和E_{d_NIR}分别为红波段和近红外波段的太阳直射辐射，E_{f_Red}和E_{f_NIR}分别为红波段和近红外波段的天空散射辐射，E_{a_Red}和E_{a_NIR}分别为红波段和近红外波段的地形反射辐射。

步骤6，构建NDVI校正公式，实现NDVI地形阴影效应校正，如图4所示。

具体地，第五步中模拟的红波段和近红外波段太阳直射光反射率补偿量与原始反射率一起用于计算校正后的NDVI，从而获得NDVI的校正公式：

NDVI_corr＝(ρ_NIR+ρ_{dc_NIR}-ρ_Red-ρ_{dc_Red})/(ρ_NIR+ρ_{dc_NIR}+ρ_Red+ρ_{dc_Red})

式中，NDVI_corr为地形阴影效应校正后的NDVI值。

图4(a)、(b)、(c)分别为验证NDVI，校正前的NDVI和地形阴影效应校正后的NDVI。可以看到，校正前后的NDVI值范围相同，表明所述方法校正后不会产生异常值，保证了校正后没有额外误差的引入。

步骤7，评估校正结果。

具体地，校正结果的评估方式包括：校正前后的NDVI与验证NDVI的密度散点图对比；根据阴影类型(本影和落影)对比阴坡与阳坡的NDVI；校正前后的NDVI与验证NDVI的剖面线对比。评价标准为：密度散点图的动态范围关于1：1线高度对称，且不发生“外溢”现象；校正后本影和落影的NDVI接近阳坡的NDVI；校正后的NDVI剖面线接近验证NDVI剖面线。

图5(a)和(b)分别为校正前和校正后的NDVI与验证NDVI之间的密度散点图。可以看到，原始NDVI与验证NDVI的散点图中，点分布集中，聚集在1:1线的两侧，具有一定对称性，但在验证NDVI＝0.8处，原始NDVI出现了降低，即动态范围“外溢”，表明阴影区域的NDVI较低。然而，校正后NDVI与验证NDVI的散点分布具有高度对称性，并且动态范围没有“外溢”，表明阴影区域的NDVI被成功校正。

为了进一步评估校正效果，选择一块子区域进行详细分析(图4中的红框部分)，如图6所示，图中的白色虚线为剖面线。图6(a)为实验影像及阴坡中本影和落影的分布情况，图6(b)为验证影像，图6(c)和(d)分别为校正前的NDVI和验证NDVI，图6(e)为校正后的NDVI。可以看到，校正前阴坡的NDVI明显低于阳坡，校正后阴坡的NDVI基本与阳坡相同，且校正后的NDVI更接近验证NDVI。另外，所述校正方法在光照非常弱的区域不会造成过度校正，且对本影和落影有着相同的校正效果，传统的地形校正模型很难达到这一效果。

剖面线验证结果如图7所示，图中的虚线为验证NDVI，实线为校正前后的NDVI，其中图7(a)为校正前，图7(b)为校正后。可以看到，校正后的NDVI剖面线与验证NDVI基本重合，表明所述方法在本影和落影区域都有良好的校正效果。

综上，本发明在此实施例中能够简单高效地校正NDVI地形阴影效应。

以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专利的技术人员在不脱离本发明技术方案范围内，当可利用上述提示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明方案的范围内。

Claims (5)

1.一种复杂山地的NDVI地形阴影效应校正方法，其特征在于，该方法包括以下骤：

步骤1，获取复杂山地的多光谱影像数据和同区域的数字高程模型DEM数据；

步骤2，根据DEM数据计算天空观测因子和地形观测因子；

步骤3，通过影像成像时间获取太阳天顶角和方位角，并由6S模型计算太阳直射辐射、天空散射辐射、地形反射辐射；

步骤4，构建阴影指数，并通过阴影指数计算反映阴影强度的可变系数α；根据步骤1得到的反射率数据，所述阴影指数SI按照公式SI＝(ρ_Coastal-ρ_Green)/(ρ_Coastal+ρ_Green+2×ρ_NIR)计算获得，求取SI的阈值c；

式中，ρ_Coastal、ρ_Green和ρ_NIR分别为海岸波段、绿波段和近红外波段的反射率；

当SI小于c时，所述α＝0；

当SI大于c时，α＝(SI-c)/(SI_max-c)

式中SI_max为SI的最大值；

步骤5，利用影像的红波段和近红外波段反射率并结合天空观测因子、地形观测因子、太阳直射辐射、天空散射辐射、地形反射辐射和可变系数α计算红波段和近红外波段的模拟太阳直射光反射率补偿量；其中，红波段和近红外波段的模拟太阳直射光反射率补偿量按照以下公式计算：

ρ_{dc_Red}＝α×ρ_Red×E_{d_Red}/(E_{f_Red}×V_d+E_{a_Red}×C_t)

ρ_{dc_NIR}＝α×ρ_NIR×E_{d_NIR}/(E_{f_NIR}×V_d+E_{a_NIR}×C_t)

式中，ρ_{dc_Red}和ρ_{dc_NIR}分别为红波段的模拟太阳直射光反射率补偿量和近红外波段的模拟太阳直射光反射率补偿量，ρ_Red和ρ_NIR分别为红波段和近红外波段反射率，E_{d_Red}和E_{d_NIR}分别为红波段和近红外波段的太阳直射辐射，E_{f_Red}和E_{f_NIR}分别为红波段和近红外波段的天空散射辐射，E_{a_Red}和E_{a_NIR}分别为红波段和近红外波段的地形反射辐射；

步骤6，构建NDVI校正公式，实现NDVI地形阴影效应校正；其中，构建NDVI的校正公式：步骤5中模拟的红波段和近红外波段太阳直射光反射率补偿量与原始反射率一起用于计算校正后的NDVI，从而获得NDVI的校正公式：

NDVI_corr＝(ρ_NIR+ρ_{dc_NIR}-ρ_Red-ρ_{dc_Red})/(ρ_NIR+ρ_{dc_NIR}+ρ_Red+ρ_{dc_Red})

式中，NDVI_corr为地形阴影效应校正后的NDVI值；

步骤7，评估校正结果。

2.根据权利要求1所述的复杂山地的NDVI地形阴影效应校正方法，其特征在于，所述步骤1中，多光谱影像数据为经过消除大气干扰的反射率数据，数字高程模型DEM数据替换为坡度和坡向数据。

3.根据权利要求1所述的复杂山地的NDVI地形阴影效应校正方法，其特征在于，所述步骤2中，通过步骤1中的DEM数据按照以下公式计算天空观测因子V_d和地形观测因子C_t：

V_d≈(1+cosθ_S)/2

C_t＝1-V_d

式中，θ_s是每个像元的坡度数据。

4.根据权利要求1所述的复杂山地的NDVI地形阴影效应校正方法，其特征在于，所述步骤3为：根据影像拍摄区域的大气模式、气溶胶含量参数、目标高度参数、传感器高度参数、光谱参数和地表反射率类型通过6S大气校正模型计算红波段和近红外波段的太阳直射辐射E_d、天空散射辐射E_f、地形反射辐射E_a。

5.根据权利要求1所述的复杂山地的NDVI地形阴影效应校正方法，其特征在于，所述步骤7为：

校正结果的评估方式为：校正前后的NDVI与验证NDVI的密度散点图对比；根据阴影类型对比阴坡与阳坡的NDVI；校正前后的NDVI与验证NDVI的剖面线对比；评价标准为：密度散点图的动态范围关于1：1线高度对称，且不发生“外溢”现象；校正后本影和落影的NDVI接近阳坡的NDVI；校正后的NDVI剖面线接近验证NDVI剖面线。

Images