基于MODIS数据的青藏高原地表反射率重建方法研究

陈善静， 张文娟， 张 兵， 康 青， 徐 旭

（1.中国科学院 空天信息创新研究院，数字地球重点实验室，北京 100094；2.可持续发展大数据国际研究中心，北京 100094；3.中国科学院 空天信息创新研究院，北京 100094；4.中国科学院大学，北京 100049；5.陆军勤务学院，重庆 401311）

1 引言

青藏高原位于亚洲内陆，是中国最大、世界海拔最高的高原，被称为“世界屋脊”。青藏高原的地表反射率在高原植被指数、冰川变化、水土流失和环境保护等方面都有着重要影响，是高原科考和区域生态环境分析的基础性科学数据［1］。通过遥感手段可以大范围获得地表多/高光谱图像，是获取典型地物目标反射率的重要途径［2］。但是由于探测器噪声和气象条件的影响，部分影像数据难以真实反映地表反射率，需要对其进行修复和重建［3-4］。

在地表反射率数据生成和重建方面，Xiao等人针对AVHRR （Advancedvery High Resolution Radiometer） 的地表反射率产品受云和气溶胶影响，造成数据时相和空间不连续，影响后续应用的问题，提出一种基于无云区域时序NDVI（Normalized Difference Vegetation Index）上包络线的遥感图像地表反射率重建方法［5］。Wu等人以青藏高原湖泊区域的分类图为基础，采用内腐蚀掩膜法对MOD09A1反射率图像中受云影响造成数据缺失的区域进行了自动填充［6］。Tseng等人针对高分辨率SPOT卫星遥感影像受云影响问题，提出以多时相辅助数据为基础，采用传统的亮度增强、颜色增强以及小波变换等方法重建高质量无云图像［7］。Cheng等人提出利用多时相遥感影像中的相似像元对受云影响目标区域进行修复和重建［8］。Wang等人利用地表覆盖分类数据对单时相MODIS数据中受云影响的异常缺失像元进行了快速修复［9］。文献［10］和文献［11］以多景时序遥感影像作为辅助数据，结合传统的PCA方法和深度学习方法，对遥感影像中受云影响区域进行了反射率数据的修复与重建。申茜等人以GF-2数据为基础，运用全色和多光谱融合、几何精校正和相对辐射一致化算法，生成了0.8 m空间分辨率2015年—2019年北京市平原区范围地表反射率影像产品［12］。孙华生等人利用无人机获取的多光谱遥感数据对地表反射率进行了定量反演和重建［13］。何兴伟等人借助MODIS地表反射率产品分析中国西北沙漠地区地表反射率特征，再结合采用手持和机载地物光谱仪实测数据，构建了中国西北地区10个沙漠辐射定标场的反射率光谱模型［14］。

整体而言，各种卫星和无人机遥感图像是地表反射率制图的重要基础［15］。由于受云、雾和探测器噪声的影响，地表反射率数据可能存在异常。许多学者利用图像空间信息、时序信息和地面采集信息开展了地表反射率数据修复与重建工作。在众多地表反射率重建方法中，基于时序信息的地表反射率重建方法对纹理特征和细节特征修复效果较好，但是对时序遥感影像的完整性要求较高。而在大范围地表反射率修复重建过程中，部分区域的多时相辅助数据可能均为异常像元，表现为多时相数据残缺，无可参考的有效地表反射率数据。此时，该类方法对地表反射率的修复质量将受到较大影响。本文以高实时性、大宽幅的MOD09A1时序遥感影像和MCD12Q1地表覆盖分类数据为基础，考虑到时序遥感影像相关性和邻近同类地物光谱相似性，以及深度学习网络模型在数据驱动下的高性能非线性图像重建与数据拟合能力，基于深度学习中的残差网络（ResNet）［16］开展多时相数据残缺情况下的大范围高原地表反射率数据重建方法研究，为青藏高原光学环境和地理科学研究提供精确完整的基础数据。

2 研究区域与数据

青藏高原地处北半球中纬度，位于我国西南部，西起帕米尔高原，东至横断山脉，东西长约2 945 km，南抵喜马拉雅山，北至昆仑山，南北宽约1 532 km，其主体区域经纬度大致为：26°15′12"N～39°46′17"N，73°19′36"E～104°38′56"E。

MODIS数据是美国国家航空航天局（National Aeronautics and Space Administration，NASA）在Terra和Aqua两颗卫星上搭载MODIS遥感器获取的对地观测数据［17-18］。根据MODIS的全球网格划分和青藏高原主体区域情况，本文所用到的MODIS数据主要包括：h24v05，h25v05，h25v06，h26v05和h26v06共计5个数据块。以青藏高原为目标区域对图像拼接和裁剪后得到研究区域如图1所示。

图1 研究区域整体情况（MO09A1 2019.7.12）Fig.1 Overview of the study area （MO09A1 2019.7.12）

本文主要利用MODIS数据中的反射率数据MOD09A1和地表覆盖分类数据MCD12Q1（http：//modis.gsfc.nasa.gov）。MOD09A1是覆盖可见光到短波红外范围的7个波段（B1-B7） 8天合成的反射率数据。MCD12Q1是根据Terra和Aqua观测所得的数据经过处理得到土地覆盖的年产品。本文采用MCD12Q1中的IGBP （International Geosphere Biosphere Programme） 分类标准，将地物分为：水体、草地、常绿针叶林、落叶阔叶林、永久冰雪区等17类。MOD09A1地表反射率数据和MCD12Q1地表覆盖分类数据基本情况如表1所示。本文中拟修复重建的地表反射率数据为2019年7月12日的MOD09A1数据，考虑到相同或相邻时间内，自然地物具有年际变化的规律，辅助数据为最邻近时间点的影像数据、前一年对应时间点影像数据和土地覆盖分类数据，详见表2。

表1 MOD09A1和MCD12Q1基本情况表Tab.1 Basic information of MOD09A1 and MCD12Q1

表2 拟重建的反射率数据及辅助数据情况Tab.2 Reflectance data to be reconstructed and auxiliary data

3 遥感数据预处理

3.1 反射率异常像元去除

MOD09A1遥感影像由于受探测器成像噪声和云的影响部分像元的地表反射率出现异常。此类异常像元影响了整个MOD09A1反射率产品的真实性和完整性。通过分析发现，MOD09A1遥感影像中异常像元主要有3类，分别为：无效像元、噪声像元和受云影响的像元。根据MOD09用户手册和其自带的波段质量控制层（surf_refl_qc_500m）以及状态层（surf_refl_state_500m）对异常像元进行定位和标记。在各数据层中异常像元标志位及处理规则如表3所示。对遥感数据所有异常像元去除以后将该位置标记为0，作为后续需要重建和修复的像元标志位。

表3 异常像元标志位及处理规则Tab. 3 Flag bit of abnormal pixels and processing rules

3.2 有效数据层提取与投影转换

MOD09A1数据文件除包括反射率数据层以外还通常包括数据质量、云覆盖、观测角度、过境时间等多个数据图层，同时MCD12Q1数据文件也包含了IGBP，UMD和LAI等多种地物分类图层。因此，为获得青藏高原地表反射率数据和对应的分类图层，对MOD09A1的B1-B7波段和MCD12Q1的IGBP图层进行提取。同时，MODIS数据产品一般为正弦投影，图像扭曲形变较为严重，需要统一转化常用的WGS-84地理坐标投影。因此采用MODIS Reprojection Tool（MRT）软件对于HDF-EOS格式、正弦投影下的MODIS产品进行重投影和采样，形成空间分辨率为0.004 231°的GeoTiff格式影像数据。在此基础上，利用ENVI按地理坐标对涉及青藏高原主体区域的5个GeoTiff影像数据进行拼接。去除反射率异常像元后的目标区域如图2所示，对应的辅助数据和地表覆盖分类信息如图3～图5所示。对比图2～图4可以看出，部分区域在目标图像和多时相辅助数据中均为异常像元（标记为黑色区域），无可参考的有效地表反射率数据，表现为多时相数据残缺。针对此类区域本文引入地表覆盖分类信息后可有效利用同类邻近地物光谱相似性，实现地表反射率修复以及重建效果的改善与优化。

图2 去除反射率异常像元后的目标区域图像（2019.7.12，黑色为标记为0的无值区域，占全图的64.98%）Fig.2 Image of the target area after removing pixels with abnormal reflectance （2019.7.12， Black is the valueless area marked with 0， accounting for 64.98% of the whole image）

图3 去除反射率异常像元后的多时相辅助数据（2019.7.4，黑色为标记为0的无值区域，占全图的54.82%）Fig.3 Multi-temporal auxiliary data after removing pixels with abnormal reflectance （2019.7.4， Black is the valueless area marked with 0， accounting for 54.82% of the whole image）

图4 去除反射率异常像元后的多时相辅助数据（2018.7.12，黑色为标记为0的无值区域，占全图的50.13%）Fig.4 Multi-temporal auxiliary data after removing pixels with abnormal reflectance （2018.7.12， Black is the valueless area marked with 0， accounting for 50.13% of the whole image）

图5 目标区域地表覆盖分类图（2019 MCD12Q1- IGBP）Fig.5 Land cover classification map of target area（2019 MCD12Q1-IGBP）

4 基于深度学习的缺失数据重建

4.1 残差网络构建

目标区域地表反射率数据去除异常像元后出现大面积数据缺失，因此，本文利用深度残差网络模型对数据缺失区域进行修复重建，提出的残差网络基本模型如图6所示。该模型将模拟数据缺失的训练样本、多时相辅助数据、地表覆盖分类信息输入残差网络，利用多层卷积和跳过连接进行特征提取、融合和重建，采用Smooth L1损失函数计算输出结果与对应无数据缺失训练样本之间的误差，通过误差反向传播，实现模型训练和参数优化。本文方法以深度学习网络模型为基础，利用多时相辅助数据和地表覆盖分类信息，对反射率大面积缺失的遥感影像进行修复和重建。引入多时相辅助数据可以有效重建待修复区域的纹理特征和细节特征，但此类方法通常对多时相辅助数据的完整性有较高要求。本文方法在多时相辅助数据残缺情况下，基于空间邻近同类地物光谱具备高相似性的普适地学规律，同步引入地表覆盖分类信息，实现了将同类地物时-空-谱多维特征信息融合用于地表反射率的修复与重建。

图6 本文提出的残差网络基本结构Fig.6 Basic architecture of the proposed residual network in this paper

本文提出的残差网络由1个前置卷积（Preconv）、16个残差模块（Residual Module）、1个后置卷积（Post-conv）组成。每个残差模块由2个3×3卷积核、1个ReLU函数和1个跳过连接构成。残差网络中卷积通道数为64，步长为1，边缘填充为1。对于输入特征图层x，经过残差模块输出结果可表示为：

其中：w为权重系数，b为偏离量，f(·)为激活函数。输出结果H(x)和输入x具有相同的特征维度。

本文中用于计算误差的Smooth L1损失函数具有收敛平稳、快速，在误差较小时避免出现持续震荡等优势，其计算公式如式（2）所示：

其中：T0为64×64像元无数据缺失的样本图像，Tout为深度残差网络重建后输出的64×64像元预测图像，γ为分界系数。

4.2 基于云掩膜的训练样本生成

训练样本是深度学习网络发挥作用的基础。本文基于拟重建的MOD09A1遥感影像中数据完整区域，随机裁剪尺度为64×64像元图像T0作为无数据缺失的真实训练样本。在MOD09A1遥感影像中选取真实有云区域进行数据抠除形成云掩膜，如图7（f）所示。该云掩膜中有云区域数值为0，无云区域数值为1。在T0上按云掩膜在随机位置进行数据抠除，模拟去除异常像元后存在数据缺失的掩膜样本TM。辅助数据为对应区域去除异常像元后的多时相MOD09A1图像T1（2019.7.4）和T2（2018.7.12），以及地表覆盖分类图TLC。TM，T1，T2和TLC同步输入到深度残差网络，得到重建后输出结果图像Tout。利用Smooth L1函数计算Tout和T0的误差。训练样本及辅助数据样例如图7所示。

图7 训练样本及辅助数据样例图Fig.7 Example of training sample and auxiliary data

4.3 基于地表覆盖分类和K-means聚类的样本数据增广

通过对4.2节构建训练样本分析发现，在青藏高原的水体、草地和永久冰雪覆盖区域经常受到云的影响造成数据缺失，因而对应的地表反射率训练样本非常少，造成了各个类别的样本数量不均衡。针对该问题，本文提出基于地表覆盖分类和K-means聚类的样本数据增广方法，对高原水体、草地和永久冰雪覆盖区域的地表反射率样本进行扩增。该方法流程如下：

（1）在有数据缺失目标图像MOD09A1（2019.7.12）中随机裁剪的尺度为64×64像元图像，并在对应的MCD12Q1中裁剪的地表覆盖分类图TLC。此时可能存在无值区域。对TLC中水体（Label=17）、草地（Label=10）和永久冰雪覆盖区（Label=15）对应类标签数量进行统计，当水体、草地、冰雪大于全图50%时保留裁剪结果和对应的TLC；

（2）结合TLC对裁剪结果中所有同类地物有值像元进行光谱维K-means聚类。中对应像元=1，2，…，n；l=1，2，…，17}，其中，为类别标签为l的第i个像元，其光谱维度d=7；

结合地表覆盖分类和K-means聚类算法对高原水体、草地和永久冰雪覆盖区域训练样本增广结果样例如图8所示。

图8 训练样本增广样例图（各行分别对应水体、草地和冰雪覆盖区域样本）Fig.8 Examples of training sample augmentation （Each row corresponds to the samples of water bodies， grasslands， permanent snow and ice area respectively）

5 实验结果与分析

为了验证青藏高原地表反射率重建方法的有效性和精度，本文分别利用4.2节和4.3节方法生成的掩膜样本，以及真实MOD09A1反射率图像开展了两组定量化和可视化的对比验证实验。掩膜样本实验中在MOD09A1青藏高原地区共裁剪生成1 340景样本，其中真实的无数据缺失区域裁剪样本581景，利用地表覆盖分类和K-means聚类生成的增广样本759景。所有样本按云掩膜进行数据抠除，形成对应的缺失数据样本。将裁剪的1 340景样本随机打乱后按3∶1∶1的比例划分为训练集、验证集和测试集。模型学习率为10-5，梯度下降优化方式为Adam，训练总次数为300。模型训练时间为5.13 h。本文方法训练过程中损失函数收敛情况图9所示。模型运行的服务器配置为双核Intel Xeon（R） Gold 5128，内存128 G，4块NVIDIA GeForce RTX 3090，Windows 10操作系统。整个模型以Pytorch库为基础上用Python（version： 3.8.12）编程实现。

图9 本文方法训练过程中损失函数收敛情况Fig.9 Convergence of loss function in the training procedure of the proposed method

在地表反射率重建结果评价指标为峰值信噪比（Peak-Signal-to-Noise Ratio ， PSNR）［19］、结构相似度（Structural Similarity Index， SSIM）、相关系数（Correlation Coefficient，CC）和光谱夹角（Spectral Angle Mapper，SAM）。对比方法分别为：（1）基于多时相数据的反射率均值重建法（Temporal_average）［7-8］，该方法将多时相遥感影像中的反射率平均值用于对应异常像元反射率重建；（2）基于地表覆盖分类信息和残差网络的地表反射率重建方法（LC_ResNet）［9，16］，该方法将地表覆盖分类信息和残差网络相结合用于单时相MOD09A1反射率数据的修复与重建；（3）基于多时相数据和残差网络的地表反射率重建方法（Temporal_ResNet）［10-11，16］，该方法利用多景时序遥感影像和残差网络对MOD09A1反射率数据进行修复与重建。本文方法在以残差网络为基础，同步引入了地表覆盖分类信息和多时相辅助影像数据，通过多源信息深度融合实现MOD09 A1反射率数据的修复与重建。

5.1 基于掩膜样本的地表反射率重建实验

基于掩膜样本进行地表反射率重建实验可以利用掩膜前反射率完整样本对重建结果进行定量化的精度评估。本文方法（Ours）和3种对比方法对地表反射率的重建精度如表4所示。针对典型裸地、湖水、草地和冰雪覆盖区域地表反射率的重建效果如图10所示。结合图表可以看出，传统的Temporal_average对地表纹理信息重建效果较好，但图像边缘效应明显，并且对时相数据完整性要求较高。当作为辅助数据的时相数据也存在数据缺失时，即多时相数据残缺时，这部分区域的反射率无法重建，相关结果如图10（h）、图10（n）和 图10（t）所 示。LC_ResNet，Temporal_ResNet和本文方法均是基于深度学习模型的反射率重建方法。在反射率重建过程中，3种方法逐步引入了地表覆盖分类信息和多时相辅助信息，重建结果的精度逐渐提高，误差也逐步缩小。相比单独依靠多时相数据的反射率重建方法，本文方法由于引入包含地理信息的地表覆盖分类数据，实现了利用残缺多时相数据对地表反射率进行精确恢复与重建，进而降低了对多时相数据的完整性的要求。

表4 地表反射率重建精度对比表Tab.4 Comparison of surface reflectance reconstruction accuracy

图10 典型地表反射率重建效果图（从上到下各行分别对应裸土、湖水、草地和冰雪覆盖区）Fig.10 Result of reflectance reconstruction of classic surface （Each row from top to bottom corresponds to barren， water bodies， grasslands， permanent snow and ice area respectively）

5.2 青藏高原地表反射率全图重建实验

本实验利用MOD09A1数据对青藏高原的真实地表反射率进行全图重建。由于受数据缺失的影响，没有大尺度地表反射率的真实参考值，因此，本部分实验主要进行定性化和可视化的重建效果评估。3种对方算法和本文方法对青藏高原地表反射率重建效果如图11～图14所示。结合修复前遥感影像图2～图4可以看出，青藏高原南部和东南部的部分区域由于长期受云的影响在多时相数据中均存在数据缺失，因此，Temporal_average重建后仍有部分区域数据不完整。在残差网络的基础上引入地表覆盖分类信息后，地表反射率重建效果得到改善，但在地表类型比较单一的区域重建反射率图像的纹理特征不够理想。利用多时相数据加残差网络的修复效果如图13所示。该方法对多时相辅助数据中均无值的区域重建效果较差，从图13右下角的放大区域可以看出，部分无值区域几乎无纹理特征，也无地物光谱变化信息。本文方法在残差网络的基础上融合了时相信息和地表覆盖信息，实现了将变化的遥感信息和较为固定的地理信息融合应用于地表反射率重建。从图14的放大区域可以看出，即使在时序数据残缺区域，结合地表覆盖信息后也仍能较好地实现对地物反射率信息的修复和重建。

图11 Temporal_average重建效果Fig.11 Result of Temporal_average

图12 LC_ResNet重建效果Fig.12 Result of LC_ResNet

图13 Temporal_ResNet重建效果Fig.13 Result of Temporal_ResNet

图14 本文方法重建效果Fig.14 Result of proposed method

6 结论

青藏高原地表反射率是自然灾害和生态环境监测的重要基础性参考数据。MODIS遥感器重访周期短，具备高效、快速、大范围获取高原地表反射率的能力。但是，MODIS遥感器获取的反射率数据由于受探测器噪声和云的影响出现大量的异常像元，影响到数据的完整性和可用性。针对该问题，本文基于邻近时序遥感影像和年际时序遥感影像具有高相关性，以及邻近同类地物光谱具备高相似性的普适的地学规律，提出了一种基于MODIS数据的青藏高原地表反射率重建方法。以MOD09A1多时相数据和MCD12Q1地表覆盖分类信息为基础，经过反射率异常像元去除、有效数据层提取和投影转换等预处理后，得到目标区域反射率图像及辅助数据。根据残差网络基本原理，构建了基于多时相数据与地表覆盖分类信息的深度残差网络。利用完整数据区域的云掩膜样本、基于地表覆盖分类和K-means聚类算法的增广样本对模型进行训练和优化，最终将训练好的模型用于地表反射率缺失数据区域的修复与重建。通过两组对比试验表明，本文方法在多时相数据残缺情况下，结合地表覆盖分类信息后重建的地表反射率的精度分别为35.199 6（PSNR），0.897 8（SSIM），0.915 6（CC），1.107 8°（SAM），均优于3种对比方法。本文方法降低了多时相影像数据量和完整性的要求，适合于青藏高原大范围地表反射率修复与重建。