融合多光谱影像的高光谱影像厚云去除方法

王蓝星，王群明,2，童小华

1.同济大学测绘与地理信息学院，上海 200092; 2.上海市数字光学前沿科学研究基地，上海 200092

遥感影像对于全球地表监测具有重要的作用。其中，高光谱影像能提供地物丰富的光谱信息，被广泛应用于环境监测、精细农业、异常目标检测、海洋监测及地表矿物制图等领域[1-5]，具有不可替代的地位。GF-5卫星于2018年5月发射，其作为全谱段的高光谱卫星，首次同时实现了对大气和陆地的综合观测，是中国高分专项中一颗重要的科研卫星，具有重大的应用价值。然而，云层(尤其是厚云)遮挡对高光谱影像的应用造成了较大的影响。研究高光谱影像的云去除问题具有重要意义。

现有的厚云去除(以下简称为“云去除”)方法大多针对于多光谱影像，主要分为3类：①基于空间信息的方法。此类方法仅利用云影像上非云区域中的有效信息(一般为云覆盖区域的邻近有效区域)对云覆盖区域进行重建[6-7]。文献[8]通过分析影像结构和纹理的空间模式，计算待重建影像上已知的最相似图块对之间的偏移量，在能量函数最小化的前提下，通过复制已知像元的方式填补缺失值。然而，由于此类方法利用的已知信息较为有限，预测结果的不确定性较大(尤其在云覆盖区域较大且纹理较为复杂时)。因此，此类方法更适用于小范围的云层去除。②基于时间信息的方法[9-11]。此类方法利用同一地区其他时间点获取的无云影像(即辅助影像)提供辅助信息，其关键在于利用辅助影像和云影像中共同无云区域的有效信息对两者之间的关系进行准确建模。然而，由于不同的地物往往随时间发生不同程度的变化，辅助影像与目标云影像之间的关系会随时间间隔增大变得愈加复杂，这给两者之间关系模型的构建造成了困难。因此，与目标云影像在时间上最接近的无云影像是辅助影像的最佳选择。能否获得时间上与目标云影像足够接近的辅助影像是此类方法能否获得令人满意的重建结果的关键因素。全局直方图匹配(global linear histogram match，GLHM)方法[12]基于辅助影像与云影像之间共同的无云区域数据构建两者之间的线性关系，从而对剩余缺失区域进行预测。GLHM方法最初为解决Landsat 7 ETM+条带问题而提出，但对于云去除问题同样适用，该方法简单易实现，能够快速对大范围缺失数据进行重建。③基于时空信息的方法。此类方法综合利用时间和空间信息，相较于前两种方法，具有更好的稳定性，能够应对不同情况下云下缺失信息的重建任务[13-16]。例如，增强的近邻相似像元插值法(modified neighborhood similar pixel interpolator,MNSPI)[17]在预测每个云像元时，利用云影像本身的无云信息获得一个空间预测项，同时利用辅助影像上对应位置的有效信息估计一个时间预测项，最后根据影像中的景观同质程度和像元信息随时间变化程度为空间项和时间项赋予不同的权重，得到最终结果。MNSPI方法计算效率高且具有较好的准确性，是一种常用的云去除方法。此外，基于机器学习(如深度学习)的方法也得到了越来越广泛的应用[18-19]，但其一般对训练样本的数量有较高的要求，如通常需要大量待预测区域之外的数据用于模型训练。文献[20]提出一种自适应回归方法，综合利用长时间序列数据(部分影像包含局部云污染)进行云去除。文献[21]基于光谱时间度量，使用K-近邻回归方法利用一年的无云影像预测缺失值。

与上述针对多光谱影像云去除的研究相比，高光谱影像的云去除面临着更多的挑战，其研究整体上远落后于多元谱影像云去除。分析上述云去除方法可以发现，辅助影像的利用是云去除的必要手段。然而，现有的高光谱数据时间分辨率往往较低。具体来说，GF-5卫星的重访周期为51 d[22]，EO-1高光谱数据的时间分辨率被设计为200 d[23]。因此，时域上最邻近的同源辅助影像与目标高光谱影像之间通常存在着较大的时间间隔，其间地物覆盖可能发生了较大的变化(例如植被的自然枯萎，农作物的轮作，城区的扩张等)，使得辅助影像的参考价值大大降低。因而，对于高光谱影像(如GF-5和EO-1高光谱影像)的云去除，寻求其他具有更高时间分辨率的辅助影像(如多光谱影像)十分必要。

Landsat系列数据是至今应用最广泛的多光谱数据之一[24-30]。目前，Landsat 8卫星仍在稳定运行，其全球周期性覆盖，较高的时间分辨率(16 d)[31]使得其更有可能提供时间上更接近于目标高光谱云影像的辅助影像。同时，其30 m空间分辨率和GF-5及EO-1高光谱影像一致。此外，Landsat 8卫星数据可供用户免费下载，是一种易获取的辅助影像。然而，现有云去除方法通常基于辅助影像与云影像波段一一区间对应的假设，即对任一含云波段，选取辅助影像中同区间的波段进行重建。但是，由于光谱分辨率不一致，多光谱影像与高光谱影像各波段的光谱区间(即光谱分辨率)设置存在较大的差异。现有方法无法用于基于多光谱辅助影像的高光谱云去除。为此，必须寻求更有效的波段映射模型以充分利用多光谱数据进行高光谱数据的云去除。

文献[32]提出一种基于空谱信息的随机森林(spatial-spectral-based random forest,SSRF)方法，利用同源多光谱辅助影像对多光谱(即Landsat和Sentinel-2)遥感影像进行云去除。SSRF方法继承了经典的随机森林(RF)方法[33]在描述因变量与自变量之间复杂非线性关系时的优势，其在对多光谱云覆盖影像的各个波段分别进行云去除时，可同时利用时域近邻影像的多个波段构建拟合模型，无须假设辅助影像与云影像波段一一对应。这些特点为异源数据的利用提供了可能。本文基于已有的SSRF方法，研究将其扩展至异源影像间的云去除情形，提出一种利用Landsat 8 OLI影像作为辅助影像对高光谱影像(如GF-5和EO-1高光谱影像)进行云去除的方法(记为SSRF_M方法)，以解决同源高光谱辅助影像因时间间隔较长而存在较大地物变化的难题。

1 本文方法

SSRF_M方法同时利用Landsat 8 OLI数据的6个波段(blue、green、red、NIR、SWIR 1和SWIR 2波段)对GF-5或EO-1高光谱影像的每个云覆盖波段单独进行重建。具体来说，SSRF_M方法在构建样本时，对于辅助数据中的有效像元，利用以其为中心的3×3像元大小的图像块在6个波段的Landsat反射率值作为样本的输入(即自变量)，对应输出(即因变量)为待预测数据中同一位置下的有效高光谱像元的反射率值。本文对每个波段缺失像元的预测如下

(1)

图1 SSRF_M流程(以高光谱单波段为例)

SSRF_M的实施主要包括两个阶段：SSRF模型训练和云下信息重建(即预测阶段)。其对高光谱影像中每个云覆盖波段单独处理，具体如下。

1.1 训练阶段

1.1.1 构建训练样本

(1)获取待重建高光谱影像的云掩膜，用于区分云覆盖区域和无云区域。原始云掩膜中通常赋予云及其阴影不同取值以便区分，在本文中将云和阴影均视为缺失区域生成新的掩膜进行预测。

(2)获取与待重建高光谱影像在时间上最接近的覆盖同一地区的无云Landsat 8 OLI影像，作为辅助影像。

(3)利用云掩膜提取辅助影像与云影像上共同的无云区域数据。同时，利用云掩膜在辅助影像上提取云影像中云覆盖区域对应位置的有效数据。

(4)利用上一步中的无云区域数据构建训练数据。具体来说，自变量(即输入)为辅助的Landsat数据中的图像块，因变量(即输出)为目标云覆盖高光谱波段中对应的无云数据。根据常用经验值，本文随机选择30%的样本进行训练，剩下的70%的样本用于评价训练模型的准确性。同时，将树的数量设置为100。

1.1.2 SSRF模型训练

将得到的训练数据输入SSRF，获得自变量(即6个Landsat波段下的3×3像元图像块)和因变量(即高光谱单个波段中的缺失像元)之间的非线性关系模型(即式(1)中的f)。本文基于R软件运行SSRF程序。

1.2 预测阶段

1.2.1 构建预测样本

利用构建训练样本阶段中步骤(3)中得到的云覆盖区对应Landsat辅助影像上的有效数据，构建测试数据(即图块自变量)，作为预测数据的输入。

1.2.2 SSRF模型预测

将上一步中的测试数据输入通过步骤1.2训练好的SSRF模型，得到云像元的预测数据。以上步骤依次用于高光谱影像各波段进行云去除。

2 试 验

2.1 试验数据

本文利用GF-5和EO-1高光谱影像，选取3个不同的区域进行模拟的厚云去除试验。每个区域均包含一景用于生成模拟云影像的无云高光谱影像，还包括用于提供时间辅助信息的一景同源高光谱影像以及一景异源(Landsat 8 OLI)多光谱影像。试验中使用的3类数据的空间分辨率均为30 m(图2)。使用数据的具体信息见表1。如图2所示，对于每个区域，本文基于无云参考影像生成模拟云影像。特别地，区域3的云掩膜基于其他区域真实的云生成，最大限度地接近于真实的云覆盖情况，而其余两个区域中，云均为人工随机生成。各区域中高光谱影像和多光谱影像显示时RGB波段均一一对应。区域1和区域2均选自中国北京，为典型的城市区域，异质性较强，纹理复杂。其中，区域1中同源辅助影像和异源辅助影像均与模拟的云覆盖影像较为相似，其原因可能在于该地区种植的植被四季常青。然而，在区域2中，同源的辅助影像与云覆盖影像之间存在明显的色调差异，其原因在于两者之间时间间隔较大，地物发生了一定程度的季节性变化。相比之下，与云覆盖影像之间时间间隔更小的异源辅助影像与前者在视觉上十分接近。区域3选自美国法戈，其作为典型的农业区域，匀质性较好，然而对比同源辅助影像和云覆盖影像看可以看出，该区域地物发生了显著的变化，整体色彩差异较大，给同源高光谱辅助信息的利用造成了一定的困难。反之，Landsat 8 OLI影像的获取时间更接近于云覆盖影像，二者在视觉上更为接近。试验中使用的高光谱影像均剔除了质量较差的波段。此后，本文分别利用与参考影像同源的无云影像以及Landsat 8 OLI无云影像作为辅助影像，重建模拟云影像的云下缺失信息。

注：3个区域高光谱影像分别以波段150、10、3，163、111、65和161、104、58为RGB组合；区域1的Landsat影像以red、blue、green，区域2和区域3的Landsat影像均以SWIR1、NIR、blue为RGB组合。

表1 厚云去除试验数据

2.2 试验内容

为论证SSRF_M的有效性，本文将其与经典的MNSPI方法进行了对比。MNSPI方法作为一种典型的主流方法，利用同源数据提供辅助信息，要求辅助影像与云影像波段之间存在一一对应的关系，不适用于利用异源辅助影像的情况。考虑到MNSPI方法仅能利用单个波段的辅助数据参与对应波段的云去除，为进一步在方法层面验证SSRF_M的有效性，本文测试了同样计算简单，且能同时利用多个波段进行预测的GLHM方法。为使GLHM方法在数据利用层面上与SSRF_M更接近，同样基于3×3局部图像块进行运算，故将其称为利用多光谱影像的空谱GLHM(spatial-spectral GLHM，SSGLHM_M)方法。总之，对于每个区域的模拟云影像，同时测试了3种云去除的方法：①SSRF_M：基于SSRF方法，利用Landsat 8 OLI作为辅助影像；②MNSPI_H：基于MNSPI方法，利用与云影像同源的其他时间点获取的无云高光谱影像作为辅助影像。③SSGLHM_M：基于GLHM方法，利用Landsat 8 OLI辅助影像所提供的空间和光谱信息。需要说明的是，本文提出的SSRF_M仅需利用单景辅助影像，具有简单易实现的优势，故并未将其与需要大量辅助影像或训练数据(例如需多景或一年内的所有可用辅助影像)的方法[18-21]进行对比。

对3种方法的预测结果进行了定性和定量评价，前者包括两种不同波段组合的视觉效果展示，后者包括4种评价指标，即均方根误差(RMSE)，相关系数(correlation coefficient,CC)，通用图像质量指数(universal image quality index,UIQI)[34]和光谱角距离(spectral angle mapper,SAM)。评价指标均基于云区域的所有像元计算。值得注意的是，前3项指标均基于单个波段进行计算，而SAM先基于单个像元运算然后取所有像元结果的均值。4种评价指标的计算公式如下。

(1)RMSE。RMSE能够衡量预测影像与参考影像之间的差异性，其值越小表明预测结果越理想，理想值为0。RMSE定义如下

(2)

式中，Pb(x,y)和Rb(x,y)分别表示大小为m×n像元的影像P和R在位置(x,y)处波段b的像元值。

(2)CC。CC反映预测影像与参考影像之间的相关性，其值越接近于1表明预测影像与参考影像越接近。CC定义如下

(3)

(3)UIQI。UIQI用于评价预测影像与参考影像之间空间结构的相似程度，其值越接近于1表明预测影像的重建质量越好。UIQI定义为

(4)

式中，σPbRb为影像P和R在波段b的协方差；σPb和σRb分别为影像P和R在b波段的标准差。

(4)SAM。SAM用于衡量预测影像与参考影像在光谱维度上的接近程度

SAM=arccos

(5)

2.3 试验结果及分析

图3中以假彩色(同一区域用同种波段组合)显示出3种云层去除方法的结果。由图3可以看出，3种方法在城市区(区域1和区域2)都以较好的视觉连续性重建了云下区域的信息。对于地物发生更为复杂变化的农业区(区域3)，SSGLHM_M的结果中非云区域与云区域之间存在较为明显的接缝，而MNSPI_H的结果中存在着明显的色调异常。例如，在图3中，MNSPI_H将区域3中部紫色的地物错误地预测为粉色和蓝色。相反地，本文提出的SSRF_M在区域3中同样表现出令人满意的性能，其重建结果在视觉上具备较好的连续性，且无明显的色调异常。为更清楚地对比3种方法的重建结果，图3对每个区域放大展示了一个子区域。可以看出，MNSPI_H方法的重建结果均存在着明显的噪声，且在区域1和区域3中存在着明显的色调异常。相较之下，两种利用多光谱影像作为辅助影像的方法(即SSRF_M和SSGLHM_M)能得到更接近于参考影像的重建结果。这表示与获取时间距目标影像更远的同源高光谱影像相比，在时间上更接近的多光谱影像能提供更准确的辅助信息。进一步地，与SSGLHM_M相比，SSRF_M的重建结果在色调上与参考影像更为接近。

注：3个区域分别以波段150、10、3，163、111、65和161、104、58为RGB组合；白色线框内为模拟的云区域。

为定量评估3种方法，将4种评价指标对应的值绘制成条形图，如图4所示，其中RMSE、CC和UIQI为其在所有波段的平均值。由图4可知，在3个试验区域中，SSRF_M重建结果的RMSE最小，CC和UIQI最大，而MNSPI_H方法的RMSE最大，CC与UIQI最小。对于SAM，SSRF_M在区域1的结果与GLHM_M几乎相同，而在区域2和区域3中，前者的结果具有较为明显的优势。此外，图5展示了3个区域数据中每个波段的RMSE值。由图5可知，在3种方法的预测结果中，SSRF_M在各个波段下的RMSE几乎始终处于最下方(即RMSE值最小)。由此可以说明，本文提出的SSRF_M能以更高的精度恢复云下缺失信息。

注：RMSE、CC和UIQI为各波段的平均值，SAM先基于单个像元计算然后取所有像元结果的均值。

图5 区域1—区域3模拟云层去除结果各波段RMSE折线图

此外，为进一步检验各方法预测光谱的准确性，在每个区域中随机选取了两个像元，并展示其在参考影像及不同方法预测结果中的光谱曲线，如图6所示。可以发现，SSRF_M预测的光谱曲线普遍与参考影像更为接近，表明其能够更大程度地恢复缺失像元的光谱信息。

图6 区域1—区域3中3种方法预测的像元光谱曲线

3 结 论

本文提出一种利用时间分辨率较高的多光谱影像(Landsat 8 OLI影像)对时间分辨率较低的高光谱影像(GF-5和EO-1高光谱影像)进行云去除的方法(SSRF_M)，一定程度上解决了同源高光谱辅助影像因时间间隔过大而可用价值偏低的问题。SSRF_M在描述因变量与自变量之间复杂非线性关系时的优势使其能够同时利用辅助多光谱影像的多个波段数据对高光谱各个云覆盖波段进行重建，无须辅助波段和云覆盖波段区间一一对应。试验表明，与利用同源高光谱数据作为辅助影像的方法(MNSPI_H)及利用Landsat 8 OLI数据作为辅助影像的线性拟合方法(SSGLHM_M)相比，SSRF_M能得到更准确的结果。值得注意的是，SSRF方法对空间信息的挖掘能力仍有待进一步提高，在更大的空间范围内筛选一定数量的相似像元是一种可能的有效方法。此外，也可考虑发展其他有效的学习模型，在预测时充分利用空间邻域信息。