基于Landsat8 OLI影像的山区植被地形落影校正方法研究

吴勇锋，江 洪

（福州大学数字中国研究院，空间数据挖掘与信息共享教育部重点实验室，卫星空间信息技术综合应用国家地方联合工程研究中心，福建 福州 350108）

地形阴影是地形影响的一种，它的存在严重干扰山区植被信息提取精度，给山区土地覆被解译以及生态参量的遥感反演带来巨大挑战［1］.目前，地形阴影校正方法主要可以概括为3类：1）基于山地辐射传输模型方法同时考虑太阳直射辐射、大气散射辐射以及周围地形反射辐射，效果较好，如Proy校正模型［2］、Sandmeier模型［3］和光学遥感反射率计算模型［4］，但模型参数较多且计算过程复杂，推广难度大；2）基于DEM数据的经验校正方法以反射率与cosi的统计关系为基础进行经验校正，包括统计-经验模型（如Teil⁃let-回归校正［5］、b校正［6］、VECA模型［7］）、归一化模型（如二阶校正模型［8］、坡度匹配模型［9］）、朗伯体反射率模型（如余弦校正模型、C校正模型［3］、SCS模型［10］、SCS+C模型［11］等）以及非朗伯体反射率模型（如Minnaert模型［12］），但该类方法只对本影有效，对落影校正效果不佳；3）基于波段组合优化计算模型的方法通过构建特殊的植被指数来获取消除地形影响后的植被信息，如波段比模型［13］、地形调节植被指数（Topography Adjusted Vegetation Index，TAVI）［14-15］、归一化差值山地植被指数（Normalized Difference Mountain Vegeta⁃tion Index，NDMVI）［16-17］、阴影消除植被指数（Shadow Eliminated Vegetation Index，SEVI）［18-21］、植被区分阴影消除植被指数（Vegetation Distinguished and Shadow Eliminated Vegetation Index，VDSEVI）［22］等，但该类模型只能获取落影的指数信息，缺乏多光谱信息.多光谱信息是遥感影像分类的基础，对山区土地利用分类至关重要.

为了获取消除落影影响的多波段反射率信息，提升山区土地利用分类精度，笔者以福建省连江县内某山地丘陵区域的Landsat8 OLI影像为例，提出一种地形落影校正方法.首先采用地表反射率计算的SE⁃VI作为阴影校正信息源，通过构建光照区SCS+C校正后红绿蓝波段反射率与SEVI间的随机森林回归模型，进而利用落影区SEVI校正地形落影，从而获取消除地形落影影响的红绿蓝波段地表反射率的光谱信息.

1 研究区与数据源

1.1研究区概况研究区位于福建省福州市东部（图1），地理坐标范围介于119°25′46″E～119°30′55″E，26°7′21″N～26°12′10″N之间，总面积为72 km2，最小坡度为0°，最大坡度大于60°，平均坡度21.7°，标准方差10.37°.该地区以山地为主，地形崎岖复杂，区内植被以常绿阔叶林和竹林为主.

图1 研究区地理位置

1.2数据源及预处理本研究采用Landsat8 OLI卫星影像、ASTERGDEM V2高程数据2种数据源，如表1所示.其中Landsat8 OLI卫星影像过境日期为2019年12月11日，行列号为119/042，太阳高度角36.98°，太阳方位角155.56°.

表1 数据源

为消除大气和光照等因素对地物反射的影响，利用Landsat8 OLI卫星影像数据元文件提供的辐射定标参数将遥感影像的DN值转换为辐射亮度值.采用ENVI5.3中的FLAASH大气校正模块对影像的辐射亮度进行大气校正，从而获得更接近真实地物的地表反射率数据，并对DEM数据进行投影变换，使其与Landsat8 OLI卫星影像具有相同的投影坐标系统.

2 研究方法

为了实现山区遥感影像地形落影校正，首先利用随机森林分类方法将研究区分为阴影区、光照区和水体，利用DEM数据和影像角度信息相结合提取本影，落影则由阴影区和本影的差运算得到；其次，利用研究区地表反射率数据进行SCS+C校正与SEVI计算，以光照区为掩膜对SCS+C校正结果和SEVI进行随机采样，生成样本集，进而利用训练样本集生成随机森林回归模型，并结合测试样本集完成随机森林回归模型精度评价.最后，结合落影区域的SEVI和随机森林回归模型，实现地形落影校正，从而获取消除地形落影影响的地表反射率信息，之后对地形落影校正效果进行评价.研究的技术流程图如图2所示.

图2 地形落影校正技术流程图

2.1落影提取山区阴影依据其形成原因可分为本影和落影，如图3所示.

图3 地形阴影示意图

本影可根据公式计算，落影由总体阴影和本影的差值求得.

其中，Sself表示阴影检测结果，若为1则为本影，若为0则为非本影；σ表示地形坡度角，π表示180°，ω表示太阳方位角，β表示地形坡向角，γ表示太阳高度角.

2.2落影校正

2.2.1SCS+C校正SCS+C校正［11］计算公式

其中，LSCS+C表示SCS+C校正之后的像元值，LT表示SCS+C校正之前的像元值，θ表示太阳天顶角，σ表示地形坡度角，c表示地形校正参数，i表示太阳入射角.

cosi的计算公式

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

2.2.2SEVI计算SEVI［18］计算公式

其中，RVI表示比值植被指数，f(Δ)表示地形调节因子，SVI表示阴影植被指数，Bnir为近红外波段反射率，Br为红光波段反射率.

f(Δ)的计算步骤：

步骤1选择具有明显地形阴影效应的样区，确保阴阳坡阴影对等；

步骤2使f(Δ)从0开始，以0.001为间隔进行循环叠加，同时计算SEVI与RVI的相关系数r1以及SEVI与SVI的相关系数r2；

步骤3当r1和r2的差值无限接近时，得到最优的f(Δ).

其中，r1为SEVI和RVI的相关系数，r2为SEVI和SVI的相关系数，n为参与f(Δ)计算的影像像元数，x为SE⁃VI，y1为RVI，y2为SVI.

2.2.3随机森林回归校正

1）在光照区随机选取样本点提取SEVI以及SCS+C校正后红绿蓝波段地表反射率数据，将SEVI数据作为自变量，将与其对应的SCS+C校正后红绿蓝波段地表反射率数据作为因变量，制作原始样本集.

2）使用Python3.7的Sklearn模块将原始样本集随机分为训练样本集和测试样本集，其中训练样本和测试样本的比例为7∶3.在训练样本集中通过Bootstrap重抽样得到K个与训练样本集相等的训练样本并生成K棵回归树，组成随机森林回归模型［23］（见式9）.在Bootstrap重抽样时，未被抽取的概率为当N→∞时，表明每次有36.8%的数据未被抽取，这些数据被称为袋外数据（Out of Bag，OOB）［24］，可用于估计预测误差.

其中，x为输入向量，θk为独立同分布随机变量，K为回归树的个数，N为训练样本集样本数.

3）回归树生长过程中，每个分裂节点随机抽取所有变量中的M个变量作为当前节点的特征子集（一般选取总变量的1/3，本文中M为1）进行分裂，并使每棵回归树得到最大限度生长，分裂过程中不需要剪枝.

4）将所有回归树的预测平均值作为最终预测结果［25］.其最终分类决策如下

其中，H(x)表示K棵回归树得到的预测值的平均值，hk(x)表示在训练集上训练出的第k棵回归树得到的估计值.

5）使用Python3.7的bayes_opt模块的贝叶斯优化（Bayesian Optimization）调参方法寻找使随机森林回归模型表现最优的超参数，并用最优超参数进行模型训练，分别得到红绿蓝波段SCS+C校正后地表反射率与SEVI的随机森林回归模型，结合落影区域的SEVI数据即可预测出落影区红绿蓝波段的地表反射率数据.

2.3落影校正效果评价模型评价是建模过程中的核心工作之一，通过计算随机森林回归模型的决定系数（公式（11））来评价模型的预测能力.

其中，yi为实际观测值为模型预估值为样本平均数，n为样本数.

红绿蓝波段随机森林回归模型的预测精度（决定系数）如表2所示.

表2 Landsat8 OLI影像各波段随机森林回归模型的预测精度

采用目视分析、统计特征分析、光谱特征分析、相对误差分析和反射率与cosi相关性分析等方法进行地形校正效果评价.相对误差绝对值（Absolute Relative Error，ARE）计算公式

其中，ARE表示相对误差绝对值，Bshadow表示阴影区域内光谱波段的平均值，Bsunlit表示相邻非阴影阳坡光谱波段的平均值.

3 结果分析

3.1目视分析原始影像、SCS+C校正后影像地以及本文方法校正后真彩色影像如图4.对比影像校正效果可知，原始影像中地形阴影效应十分明显，阴影区地表反射率偏低，光谱信息缺失，如图4a所示.经SCS+C校正后，地形阴影效应削弱，本影区域光谱信息得以恢复，但仍有少部分落影区域光谱信息难以恢复，如图4b所示.经本文方法校正后，地形阴影效应基本消除，落影区光谱信息进一步恢复，光照区和阴影区的光谱差异减小，如图4c所示.

图4 地形校正前后效果对比

3.2统计特征分析为了定量分析地形校正方法的校正效果，分别统计地形校正前后影像均值、标准差以及变异系数（Coefficient of Variation，CV）［26］，并进行对比分析.影像标准差反映的是影像各波段的离散程度，地形校正后，各波段影像像元间的差异缩小.若校正后影像标准差小于原始影像，表明该校正方法具有一定的校正效果.影像变异系数是指其标准差和平均值之比，地形校正后变异系数越小，表明地形校正效果越好.由地形校正前后研究区影像各波段均值、标准差及变异系数可知，SCS+C校正后，影像各波段标准差及变异系数相对于原始影像均有所下降，说明SCS+C校正都有一定的校正效果.本文方法校正后，影像各波段标准差及变异系数均小于SCS+C校正，说明本文方法弥补了SCS+C校正的不足，使影像像元间的差异进一步缩小，如表3所示.

表3 地形校正前后影像统计参数对比

图5 局部地形校正前后效果对比

3.3光谱特征分析依据阴影检测结果，选取30组本影、落影及相邻非阴影阳坡样本，并分别统计30组样本在原始影像、SCS+C校正以及本文方法校正后红绿蓝三波段的光谱均值.统计结果显示，原始影像的红绿蓝波段光谱值在本影和落影处均低于相邻非阴影阳坡.SCS+C校正后本影处的红绿蓝波段光谱值恢复到阳坡水平，而落影依旧低于阳坡，证明SCS+C校正对本影的校正效果良好，但对落影的校正效果欠佳.经本文方法校正后，落影处的红绿蓝波段光谱值校正到阳坡水平，证明本文方法能够对落影进行校正.

3.4本影、落影相对误差分析地形校正前后，红绿蓝波段对应的本影、落影相对误差绝对值统计结果表明，如表4所示：

表4 红绿蓝波段本影落影区域相对误差绝对值

1）未经地形校正的各光谱波段本影、落影与其相邻阳坡的相对误差较大，相对误差绝对值在30%～50%左右.

2）SCS+C校正后，红绿蓝波段本影与其相邻阳坡的相对误差大幅降低，分别从40.40%、43.43%、29.28%降至8.75%、13.35%、6.28%，而落影与其相邻阳坡的相对误差降幅相对较小，分别从48.76%、51.30%、38.50%降至38.25%、41.23%、29.84%，表明SCS+C校正对本影校正效果较好，而对落影校正效果欠佳.

3）本文方法校正后，红绿蓝波段落影与其相邻阳坡的相对误差大幅降低，分别从48.76%、51.30%、38.50%降至0.43%、1.82%、1.91%，表明本文方法在落影区域能够取得较好的校正效果.

图6 30组样本红绿蓝波段折线图

3.5反射率与cosi相关性分析由于地形效应影响，原始遥感影像地表反射率与cosi之间存在一定相关性.地形校正后，此相关性将被削弱甚至消除，故地形校正前后两者决定系数的变化被广泛用于检验地形校正的效果［21］.地形校正前后各波段反射率与cosi的斜率及决定系数如表5所示.影像校正前，各波段像元值与其cosi都存在一定相关性.SCS+C校正后，研究区影像各波段斜率和决定系数均显著降低，表明SCS+C校正方法能够削弱地形阴影效应.本文方法校正后，斜率和决定系数的下降幅度都高于SCS+C校正，表明本文方法能够弥补SCS+C校正的不足，使整体校正效果更优.

表5 地表反射率与光照系数线性相关性分析

4 结束语

结合SCS+C模型和SEVI指数完成山区植被地形落影可见光波段校正，并对本文方法及传统地形校正方法进行比较.未经地形校正的原始影像地形阴影效应明显，本影和落影与相邻阳坡的相对误差在30%～50%左右.因此，利用遥感影像研究崎岖山区植被时，要考虑地形阴影效应.SCS+C等传统地形校正对本影的校正效果较好，但对落影校正效果欠佳.本文方法对山区落影区域植被可见光多光谱信息校正效果明显，弥补了传统地形校正方法对落影校正失效的缺陷和SEVI只能获取消除地形阴影干扰的单一指数信息的不足.本文地形落影校正方法目前对可见光波段效果较好，对其他波段的适用性还有待进一步研究.