一种快速修复Landsat影像条带色差的方法

闫星光, 李 晶, 闫萧萧, 马天跃, 苏怡婷, 邵嘉豪, 张 瑞

中国矿业大学(北京)地球科学与测绘工程学院, 北京 100083

引 言

近几十年, 卫星遥感影像成为自然资源调查、 植被生态反演、 生态环境调查的主要数据来源[1-3]。 高质量的遥感影像应具备镶嵌精度高、 信息丰富和色调和谐3个条件; 但由于影像在拍摄过程中的受到气候变化、 季节更替等外部因素以及传感器自身和卫星重返周期的内部因素影响, 不可避免地导致影像产生光照分布不均匀的情况, 使得各景影像间的灰度值存在不同程度的差异[4-5]。 特别是在长时间序列的中、 大尺度遥感监测中, 通常需要对比多期遥感影像来判别地物的变化情况。 以月、 季和年等不同时间尺度镶嵌而成的影像数据, 多会产生影像灰度值不同的问题和相邻影像条带的色调不均匀的现象。

遥感影像匀光处理技术方法归结起来可以分为三类: 基于照度与反射模型算法、 基于加性噪声模型算法和基于统计型的方法[6]。 基于照度与反射模型算法和加性噪声模型算法方法都是按照数学模型来模拟影像中的亮度分布, 对单景影像内的色调进行不同程度的补偿。 例如Mask匀光算法[7]、 Retinex算法[8]和Wallis算法[9]等, 虽然使用Mask匀光算法和Retinex算法能得到较好的效果, 但影像内部较暗的区域容易出现色差和影像边缘偏色的现象, 从而导致影像灰度值失真的情况; Wallis算法依靠线性关系进行传递, 能使原有影像的整体亮度和色调不发生较大变化, 但因线性累积误差的传递, 使得离基准影像较远的影像容易出现偏色现象。 基于统计型的方法是统计影像的特征和灰度值, 按照一定的规则使影像间不同位置的色调和亮度达到一致, 多以线性变换法、 直方图匹配法和基于均值和标准偏统计的方法为主, 其中, 线性变换法处理后的影像虽能同时考虑区域范围内多景影像的色调一致性, 但无法保持影像非线性的特点, 会导致影像产生局部色调偏差; 基于直方图方法的匀色处理, 能很好的解决多景影像灰度值分布不连续的现象, 适用于影像内部特征简单的影像; 基于均值和标准偏统计的影像修复方法, 仅能通过单一的统计值来进行影像修复, 缺乏影像不同特征区域的精确表达。 有多景影像灰度值的遥感影像的余光处理过程中, 主要存在以下三个方面问题: 其一, 中、 大尺度多景影像拼接所造成的条带色差研究中, 现有的匀光算法在处理后很难实现影像条带间亮度和色调的基本一致; 其二, 在长时序遥感监测中, 基于多源遥感的影像修复规则不一, 无法有效保证当期影像质量的客观性和长时序影像反演的精度; 其三, 现有的研究方法多基于本地ERDAS和GeoDodging等软件, 在长时间序列的中、 大尺度遥感影像匀光处理中, 难以高效和快速的实现影像修复的批量化处理[10-12]。 在长时间序列和中、 大尺度范围的遥感影像匀光处理中, 如何高效快速地实现基于当期同源遥感影像的修复问题亟需解决。 因此, 提出了一种基于随机森林方法的直方图匹配方法, 减少影像直方图形态特征差异, 从而降低影像失真的现象。

随着云计算的兴起, 大量的遥感监测和光谱指数分析都在云平台上实现。 Google Earth Engine(GEE)是卫星遥感影像可视化计算和分析处理的云平台, 包含了Landsat、 MODIS和Sentinel系列等众多卫星影像数据, 现有超600个公共数据集[13]。 特别是在长时间序列、 大范围的遥感监测研究中, GEE平台凭借其强大的云计算能力, 极大的提升影像处理时间和工作效率。 已有诸多学者在GEE平台完成众多的遥感监测和生态反演[14-15], 但现有多数研究中, 都只是将影像进行辐射修复和大气修复后直接进行植被反演和分析, 忽视由影像拼接造成的条带色差的问题。 针对上述问题, 我们以山西省为研究区, 基于GEE云计算平台采用随机森林的直方图匹配方法, 消除镶嵌后NDVI影像产生的色调不均衡现象, 并将1986年—2020年修复前后的NDVI影像进行对比分析, 旨在为长时间序列的遥感影像自动化匀色处理寻求一种简单、 高效、 适应性强的方法。

1 研究区域及数据预处理

1.1 研究区域

如图1所示, 山西省位于中国华北地区(图1), 介于北纬34°34′ N—40°44′ N, 东经110°14′ E—114°33′ E之间, 总面积15.67 km2。 山西省呈现“两山夹一川”的态势, 内部起伏不平, 山区面积占总面积的80.1%, 是典型的黄土覆盖的山地高原, 地势东北高西南低。 山西省地跨黄河、 海河两大水系, 河流属于自产外流型水系。 黄河流域生态保护和高质量发展成为国家重要战略, 山西省是黄河流域的重要组成部分, 2021年10月08日中共中央、 国务院印发《黄河流域生态保护和高质量发展规划纲要》中新的战略布局当中, 在构建黄河流域“一轴两区五极”的发展动力格局中, 山西省处于“两区”即粮食主产区和能源富集区的重要地位。

图1 研究区示意(a), 山西省位置 (b)

1.2 数据处理

Landsat和MODIS系列影像数据获取均来自于GEE云平台公共数据集。 图2为技术流程, 主要分为数据预处理、 影像匹配和验证分析三个部分。 数据预处理过程主要包括: 上传山西省矢量边界; 逐年影像时间筛选(每年1月1日至12月31日)和波段选择(“Red”、 “NIR”和“pixel_qa”); “pixel_qa”波段去云和NDVI波段计算; 最后, 按照qualityMosaic函数进行影像拼接和镶嵌。 影像匹配部分, 首先, 通过目视解译的方法, 判断逐年NDVI影像是否需要进行影像匀光处理, 将需影像修复的部分作为目标影像, 参考影像为目标影像相邻的区域, 目的是让目标影像获取和相邻影像的一致的色调; 其次, 分别统计参考影像和目标影像NDVI的DN (digital number)值, 进而计算概率密度函数和累积分布函数, 按照随机森林方法进行直方图匹配, 从而获得匀光处理目标影像。 验证分析分为三个部分: 第一, 山西省逐年NDVI影像修复前后对比; 第二, 1986年—2020年匀光处理前后NDVI影像的时序分析; 第三, 匀光处理后的NDVI影像和2000年后MODIS影像对比分析。

图2 卫星遥感影像数据处理流程和分析图

山西省边界远超出单景影像的覆盖范围, 因此在影像拼接中将由多幅不同轨迹的影像组成, 目标影像选取的原则是按照小于研究区总面积的50%进行修复, 以最大程度的减小匀光处理后影像对原始数据的影响。 经统计, 所有年份的目标影像面积占比均小于总研究区的30%, 参考影像色差所选取的区域均为目标影像的衔接条带。

在影像数据的选取上, 根据现有长时间序列的研究, 所使用的影像数据多为Landsat TOA和SR影像数据, 最终选取Landsat 5 TM、 Landsat 7 ETM+和Landsat 8OLI作为研究数据, MODIS(MOD13Q1、 MOD13A1和 MOD13A2)数据作为2000年修复后Landsat影像的验证对比, 具体影像数据信息如表1所示。

表1 Landsat和MODIS影像数据列表

2 研究方法和原理

2.1 直方图匹配算法

直方图匹配算法, 又称直方图规定化匹配方法, 是一种通过直方图变化对图像进行视觉增强的方法, 多用于图像增强。 直方图匹配用于卫星遥感影像之间的均衡匀色处理时, 是以一幅影像为参考影像, 将其通过直方图映射的方式修复目标影像, 从而达到多幅影像的色彩一致性[16]。

直方图匹配主要有3个步骤: (1)统计参考影像和目标影像逐像素DN值个数; (2)建立两者的灰度累积直方图; (3)通过一定的映射规则, 从而使目标影像获取和参考影像相近的累积分布概率[17]。 设参考影像和目标影响的概率密度函数与对应的累积分布函数分别为Pr(gr),Pt(gt),Hr(gr),Ht(gt), 首先, 要使目标影像的灰度值等于参考影像的灰度值首先满足式(1)—式(3)

Hr(gr)=Ht(gt)

(1)

(2)

(3)

2.2 随机森林方法

随机森林是通过建立多个独立的决策树, 建立分类树的集合, 其目的就是通过大量的基础树模型找到最稳定可靠的结果[18]。 无论所需处理的是离散型数据, 还是连续性数据, 随机森林方法都可以有效的处理。 本文所需匀光处理的逐年参考影像位置都是不定的, 影像波段DN值数据皆无固定规律。 因此, 调用GEE云平台中随机森林函数ee.Classifier.smileRandomForest()参与直方图匹配, 分别将参考影像和目标影像DN值及其概率分布作为训练样本, 将分类后的目标影像的值按照参考影像DN值的概率分布作为直方图匹配的映射规则进行直方图匹配。

3 结果与讨论

3.1 逐年影像分析

通过对1986年—2020年逐年NDVI影像进行筛选, 最终共找出1986、 1987、 1988、 1989、 1991、 1993、 1994、 1995、 1996、 1997、 1999、 2000、 2003、 2005、 2006、 2008、 2009、 2010、 2011、 2017年共计20年的影像存在条带色差问题。 分别对逐年的Landsat TOA和SR影像数据进行直方图匹配修复, 并对比修复前后的NDVI均值和标准差进行了分析, 结果如表2所示。

表2 山西省陆地卫星图像修复前后的NDVI值和标准偏差

1986年—2020年, 经修复后的影像精度有不同程度的提高, 其中修复后的Landsat TOA影像NDVI平均值比修复前提高了0.008 1, 影像的平均标准差下降了0.005 8; 修复后的Landsat SR影像NDVI平均值比修复前提高了0.009 2, 影像的平均标准差下降了0.004 4。

以1994年修复前后影像的进行分析为例, 1994年Landsat SR和TOA影像修复前后的NDVI平均值分析结果见表3。 Landsat TOA和Landsat SR图像修复前相比, 影像修复区的NDVI平均值分别提高了0.111 2和0.121 9, 提高比例的分别高达32.6%和29.03%。 因目标影像面积只占总面积的6.76%, 山西省整个研究区的NDVI值变化并不明显, Landsat TOA和Landsat SR图像NDVI值仅分别增加了0.007 7和0.008 5。 总体来看, 修复后的目标影像更加接近于山西省整体的NDVI值。 图3为修复前后的影像, 从图像的平滑度来看修复后的图像整体过度更加平滑, 能更好地反映山西省归一化植被指数的真实情况。

表3 1994年影像修复前后NDVI值

图3 1994年修复前后的NDVI图像

1994年影像修复前后NDVI剖面的对比分析如图4所示, 修复后的NDVI影像能平稳的反映该区域的NDVI值。 修复前Landsat TOA影像的拟合度为0.008 4, 修复后为0.170 7; 修复前的Landsat SR影像的拟合仅有0.002 4, 修复后为0.120 4。 结果表明, 修复后的影像拟合度更高, 表明修复后的影像比修复前更加符合影像的整体过度。

图4 NDVI指数剖面前后的1994年图像比较

3.2 影像时序分析

为了验证影像修复前后对于山西省34年间NDVI的影响, 分别对影像修复前后的NDVI值进行时序分析, 如图5所示。 一元线性回归的趋势性分析结果表明: Landsat TOA影像修复后的NDVI值斜率为0.006 2小于影像修复前, 拟合度R2为0.810 6高于修复前; Landsat SR影像修复前的NDVI斜率为0.007 1, 而修复后的斜率为0.006 7, 且拟合优度R2值0.836 3大于修复前的0.829 3, 表明影像修复后的拟合度更高。 整体上看Landsat影像修复后的结果在长时间序列的变化波动性更小, 趋势更加平滑。 修复后的Landsat SR影像比Laodsat TOA影像提升幅度更明显。

图5 1986年—2020年NDVI的时间序列分析

3.3 多源影像对比分析

为了验证影像修复后的准确性, 将2000年后修复的Landsat图像分别与MODIS系列数据(250 m、 500 m和1 km)进行相关性分析。 结果显示, 2000年、 2003年、 2005年、 2011年和2017年的影像都有明显改善。 图6和图7分别显示了2003年影像修复前后的Landsat SR和TOA影像与MODIS(250 m、 500 m和1 km)影像。 直方图方法修复后的Landsat影像能更好地反映色彩平衡, 整体视觉效果更加符合MODIS影像NDVI的分布特征。

图6 2003年修复前后的NDVI图像

图7 2003年不同分辨率的MODIS验证图像

2003年修复后的Landsat影像和MODIS影像的Pearson相关性分析如表4所示, 修复后的Landsat SR和TOA与MODIS影像的相关系数分别提高了0.049和0.061(p<0.05); 修复后的Landsat SR与MOD13Q1、 MOD13A1、 MOD13A2影像相关系数分别提高了0.050、 0.047和0.049; 修复后的Landsat TOA与MOD13Q1、 MOD13A1、 MOD13A2影像的相关系数分别提高了0.066、 0.060和0.059。 结果表明, 经图像修复后的Landsat数据与MODIS影像更加一致。

表4 Landsat图像和MODIS图像的Pearson相关分析

修复前后Landsat和MODIS影像NDVI平均时间序列变化的对比结果见图8。 与修复前相比, 修复后的Landsat影像与MODIS影像的整体趋势更加一致, 其中修复后的Landsat TOA影像的拟合度由0.812 5提高到0.871 1; 修复后的Landsat SR图像的拟合度由0.821 8提高到0.853 7。

图8 2000年至2020年的NDVI时间序列分析

4 结 论

目前, 在影像修复过程中, 现有的研究多是在研究区的影像镶嵌和NDVI计算之前完成, 这大大限制了影像的处理速度。 本工作针对归一化植被指数NDVI影像拼接后存在的影像带状斑块效应和色彩不均匀问题, 以山西省作为研究区, 利用GEE平台调用随机森林函数提出一种基于云端快速进行直方图影像匀光处理的方法, 极大的提高了影像修复的效率。

基于同源影像的直方图匹配能最大程度地保留当期原始影像的DN值和色彩亮度, 且在长时序研究中使当期影像的NDVI反演结果更加准确, 经过该方法处理后影像的色彩一致性较好, 同时无需考虑研究区的地理差异和空间异质性。 此外, 通过对比1986年—2020年逐年影像修复前后的结果, 经过本方法修复后的影像在长时间序列的植被监测过程中能更精确、 可靠的得出影像的修复结果, 有效减少NDVI值在长时间序列的突变, 提高长时间序列分析的准确性和稳定性。 本方法能有效改善影像条带色差较大的区域, 但对于影像条带边界不明显的区域识别仍需提升, 后续研究的重点将围绕影像色差边界的自动识别和修复展开。