CN112102180B - 一种基于Landsat光学遥感图像的云识别方法 - Google Patents

Abstract

本发明公开了一种基于Landsat光学遥感图像的云识别方法，属于图像处理领域，特别是遥感图像的清晰化处理。本发明中，提出了如何提高Landsat地表反射率(SR)云图像的除云性能的方法。从逻辑上讲，在执行去云操作前对QA进行进一步处理是解决这个问题的直接方法。因而开发了一种新的QA波段预处理方法称为Auto‑PCP，新方法旨在通过自动识别剩余的潜在云污染像素来修改QA波段。

Description

一种基于Landsat光学遥感图像的云识别方法

技术领域

本发明属于图像处理领域，特别是遥感图像的清晰化处理。

背景技术

Landsat提供了长时间的中等空间分辨率时间序列图像(Landsat 5/7/8:1984-)，已被广泛用于区域和全球陆地表面制图等。然而，在Landsat图像应用中一个严重障碍是频繁的云污染，平均而言，在地球表面每时每刻大约35％的面积都会被云层所覆盖，从而导致Landsat等光学遥感影像数据受到污染。云层覆盖是光学遥感图像缺失信息的主要来源，成为后续影像解译的难题。因而云的去除提高了光学遥感图像在不同应用和分析中的潜力，去云相关算法一直是光学遥感领域的研究热点，人们开发了大量的方法来去除光学图像中的云。

在现有的大多数厚云去除方法中，重建一个云区域的缺失值需要使用云区域附近的无云像素提供必要的辅助信息。例如，基于空间的去云方法主要是通过空间插值从邻近的无云像素中估计云像素的值。一般来说，基于空间的方法不适用于大云或者非均匀的云污染区域。最近，越来越多的方法使用在其他日期获取的一张或几张无云图像(称为“参考图像”)来辅助云图像中的云去除。这些基于多时间的方法原理大多是利用云图像中的无云像素来量化云图像和参考图像之间的反射率变化。从逻辑理论上讲，如果使用云图像中的无云像素，这些像素的质量将很大程度上影响云去除的性能。换句话说，即去云的前提是首先识别云污染像素，包括云和云阴影，才能满足后续去除厚云的要求。

目前，Landsat反射率产品中包含质量评估波段(QA)，其中标记了云和云阴影像素。美国地质调查局(USGS)基于Fmask(Function of mask)算法(3.3版)生成了Landsats4-8的QA波段。通过详细的比较和验证表明了Fmask在陆地卫星图像中云和云阴影检测方面的优越性。由于每景图像都提供了QA波段，利用现有的厚云去除方法恢复Landsat云图像中的云似乎是一项简单的任务。然而在的实践中，发现在模拟云图像上的去云性能通常很好(例如，假设无云图像中的某些区域被云覆盖)，但在一些真实云图像上发现去云效果较差。在这里给出一个示例来清楚地说明这个问题，对于2019年8月1日获取的华北平原Landsat云污染图像(图1A)，采用了改进邻域相似像素插值方法(MNSPI)进行恢复云区域，发现去云图像的视觉效果较差(见图1B中的“白色”和“黑色”块)。然后，将该云污染图像上的云覆盖在2019年8月17日获取的无云图像上，生成了一幅模拟云图像(图1C)。将MNSPI应用于模拟云图像去云中，取得了更好的去云性能(图1D)。由于在模拟云图像上的除云效果较好，因而可知在真实的Landsat 云图像上，云去除性能差的原因并不是云的形状或景观的异质性。最后通过对比这两幅图像，发现在真实的Landsat云图像中有一些云污染的像素没有被正确地标记在QA波段里。相反，在模拟云图像中可以确保所有的云污染的像素被正确地标记。由于云和云阴影掩膜的漏检误差，导致在真实云图像上的云去除性能较差。通过实验结果表明，Landsat产品中的QA波段还不能满足实际应用中去云的要求。

发明内容

本发明中，提出了如何提高Landsat地表反射率(SR)云图像的除云性能的方法。从逻辑上讲，在执行去云操作前对QA进行进一步处理是解决这个问题的直接方法。因而开发了一种新的QA波段预处理方法称为Auto-PCP，新方法旨在通过自动识别剩余的潜在云污染像素来修改QA波段。

Auto-PCP使用与去云除操作相同的数据，从而可避免数据收集时的额外负担。Auto-PCP 只需要一张无云的参考图像，这也是基于多时间的去云方法中所需的最小图像数量。此外， Auto-PCP不使用热红外图像和卷云图像，在光学图像进行去云操作时不收集这些图像。 Auto-PCP是基于QA波段的基础上发展的，因此需要预定义更少的阈值，这使得在QA波段上进行预处理成为一个简单的操作。Auto-PCP在Landsat地表反射率云污染图像中进行开展，为了更清楚的说明这个方法，使用符号P_C和P_CS表示云图像QA波段中标记为云和云阴影的像素，其余像素表示为符号P_other，Auto-PCP的目标是标记P_other中漏检的云和云影像素以提高云去除性能。因而本发明一种基于Landsat光学遥感图像的云识别方法，该方法包括：

步骤1：识别初始云像素

关于云图像和参考图像之间，云像素在蓝光波段有较大的反射率变化；云像素在蓝光波段的反射率变化值ρ_blue(change)为：

ρ_blue(change)＝max(ρ_blue(cloud)-ρ_blue(ref′))-(ρ_blue(cloud)-ρ_blue(ref′)) (1)其中，

和

分别表示像素i在云图像和参考图像中蓝光波段的反射率值；使用公式(2)来减小由于地表变化引起的云污染图像和参考图像之间的反射率差异，ρ_blue(cloud) 表示在云图像中像素蓝光波段的反射率值，ρ_blue(ref)表示参考图像中像素蓝光波段的反射率值，P_other表示云图像QA波段中标记为非云和云阴影的其他像素像素；在公式(1)中 max(ρ_blue(cloud)-ρ_blue(ref′))表示图像像素间反射率变化的最大值，这一项的引入是为了使像素间反射率变化值ρ_blue(change)大于0，与此同时，云像素的反射率变化值ρ_blue(change)更小；

云像素的另一个特征是云图像中蓝光波段的反射率值大；对于云像素，结合ρ_blue(cloud)值和ρ_blue(change)值，计算云指数：

结合ρ_blue(cloud)和ρ_blue(change)的归一化形式，CI的取值范围在-1到1之间，对于云像素来说CI值比非云像素大；基于以下两个标准在P_other像素中确定了初始云像素：

1)：初始云像素i的CI值

与P_C像素集的CI中值median(CI(P_C))的加权做比较，表示如下：

其中，median()表示取中值；

2)：使用相对阈值来分别确定每类的初始云像素；利用k-means方法对参考图像进行聚类；分类的数目设置为4-5类；以其中的一类为例，在类A中初始云像素i的CI值表示为

必须满足以下条件：

其中mean(·)和std(·)分别是均值函数和标准差函数，CI(P_other)表示P_other的CI值；参数a为标准差倍数，可知a的值越小将会识别更多的初始云像素，但同时也增加了错检误差；经过测试实验，将a参数值设为2；

步骤2：识别初始的云阴影像素；提出云阴影指数CSI识别初始云阴影，CSI的核心理念与 CI相似，CSI考虑云图像中云阴影像素近红外波段的反射率值较低(ρ_nir(cloud))；并且同时，CSI 描述了云阴影像素在云图像和参考图像之间近红外反射率下降的特征，CSI表示如下：

其中：

ρ_nir(change)＝(ρ_nir(ref′)-ρ_nir(cloud))-min(ρ_nir(ref′)-ρ_nir(cloud)) (7)

ρ_nir(cloud)表示云图像中像素在近红外波段的反射率值，ρ_nir(ref′)表示参考图像中像素在近红外波段的反射率值，

表示云图像中，属于P_other集合的像素i在近红外波段的反射率值，

表示参考图像中，属于P_other集合的像素i在近红外波段的反射率值， min(·)是最小值函数，CSI取值范围为-1至1，对于云阴影像素来说，它的值更小，越靠近-1；如果在P_other中的像素i

满足以下两个条件则为初始云阴影像元：

1)：

与P_C像素集的CI中值median(CSI(P_C))的加权做比较，表示为：

2)：分别确定每类的初始云阴影像素；以其中一类为例，在类A中初始阴影像素i的CSI值

必须满足下列条件：

其中参数b是标准差倍数；其值为2。

步骤3：匹配初始云和云阴影像素；

步骤3.1：云总是伴随着云阴影:利用云、云阴影和太阳之间的几何位置关系来进一步细化初始云和云阴影像素；设初始云像素i位于(x,y)，根据下面的方程对应的云阴影像素为(x′,y′)：

x′＝x-H_{cloud_i}×tanθ×sinφ (10)

y′＝y+H_{cloud_i}×tanθ×cosφ (11)

其中θ是太阳天顶角，φ是太阳方位角，这两者都由Landsat图像元数据文件提供，H_{cloud_i}表示位于图像(x,y)处此像素的云高；以前估算云高的方法大多是根据大气温度从地面到云的递减率来估计每个云块的高度；然而，本发明没有使用热红外图像来估计H_{cloud_i}，这是因为QA 波段中漏检的云通常都是薄云，薄云的亮度温度容易受到地表的影响；作为一个解决方案，

步骤3.2：基于QA波段估算H_{cloud_i}范围，假设H_{cloud_i}在QA波段标记云像素的云高度范围内；为此，根据8邻域的连通性在QA波段生成云斑块和云阴影斑块；对于云斑块j的云高度H_{cloud_pj}可以从云斑块和其对应的阴影斑块之间的水平距离进行近似估计denoted asDIS_pj, 表示为：

步骤3.3：计算DIS_pj(图2)；

步骤a：将阴影斑块沿着云投影方向移动，当云与阴影斑块的重叠区域达到阴影斑块面积的三分之一则停止移动；此时移动的水平距离是假定为x₁，阴影斑块与云斑块重叠的边界记为曲线ab；

步骤b：将曲线ab与阴影斑块对应的外边界a′b′进行匹配，然后基于两条曲线相似的形状的特性，通过线性插值使两条曲线具有相同的长度，并计算它们之间的相关系数；如果相关系数大于等于0.9，就认定匹配成功，此时两条曲线之间的水平距离表示为x₂；最终，DIS_pj是x₁与 x₂的水平距离之和(图2)；如果相关系数小于0.9，则计算QA波段里n个云斑的高度 {H_{cloud_p1},..,H_{cloud_pn}}，则可能的云高范围H_{cloud_i}为满足公式13条件的任意值：

min{H_cloudp1,..,H_cloudpn}≤H_{cloud_i}≤max{H_{cloud_p1},..,H_{cloud_pn}} (13)

为了避免H_{cloud_i}受到异常值的影响，从{H_{cloud_p1},..,H_{cloud_pn}}中删除了最大1％的云高值； Auto-PCP允许初始云像素在可能的云高度范围内变化，从而通过使用不同的高度匹配所有初始云和云阴影像素；

若初始云阴影像素找到相匹配的云像素，则保留初始云阴影像素；如果初始云像素的投影像素被其它云覆盖；则从被覆盖的云像素开始投影并确定下一个投影像素，该过程一直持续到投影像素为云阴影像素或无云像素；如果最后投影的像素为云阴影像素，则保留该初始云像素，否则，初始云像素将被移除；此外，在图像边缘的一些初始云和阴影像素的投影位置可能在图像之外，保留这些初始的云和阴影像素；

步骤3.4：删除8邻域连接的像素集小于7个像素的小云斑块和阴影斑块。

进一步的，所述步骤1公式5中a参数值设为2，步骤2公式10中b参数值设为2。

为了说明该方法的技术效果，在四个实验区域对该方法分别进行了真实实验测试与模拟实验测试，通过使用三个不同的QA波段来比较云去除的性能:(1)原始的QA波段(简称“QA_original”)；(2)云和云阴影边缘周围扩充两个像素处理的QA波段(简称“QA_dilation”)；(3)Auto-PCP处理的QA波段(简称“QA_Auto-PCP”)，以此来说明Auto-PCP 方法的有效性。这四个实验区域分别代表两季农作物，常绿林，城市，高寒草地四种植被类型。

在真实实验中，对于四个研究区域，比较了真实云图像去云的效果(图3-6)。每个图中展示了分别使用三种不同的QA波段，利用MNSPI方法去云得到的结果图，在原始的QA波段中，每个区域都存在着不同程度的云漏检现象(见每张图的放大图)。然而在青藏高原区域中， QA原始波段错误地将许多雪像元识别为云像元，但同时也存在云像元的漏检情况(图6)。最后通过比较三种方法恢复云的视觉效果，发现Auto-PCP去除云的性能明显更好，这要归功于对云和阴影的有效识别。在模拟实验中，基于模拟云图像进行定量评价，使用Fmask3.3生成每幅模拟云图像的原始QA波段。采用两项指标进行定量评价：一是均方根误差(RMSE)，它主要量化预测的反射率值与真实值之间的差异，二是结构相似度指数(SSIM)，主要量化预测图像与真实图像之间的整体图像结构相似度。图8显示了六个波段的平均RMSE和SSIM值，发现，QA_Auto-PCP在所有测试区域中的RMSE值最小，SSIM值最高(图8)。

附图说明

图1中(A)为2019年8月1日华北平原实测陆地卫星云图，(B)为真实云图像上的厚云去除，(C)为将(A)图像中的云叠加在2019年8月17日获取的无云图像上生成的云模拟图像， (D)为模拟云图像上的厚云去除。

图2中显示云斑块与对应阴影斑块之间几何匹配的示例图，DIS_pj是云斑块j与匹配阴影斑块之间的水平距离。

图3为中国华北平原测试区真实实验结果，上面三张图为真实云图像，由Auto-PCP识别的云污染像素(QA_Auto-PCP)，以及它们的放大视图；下面三张图为通过MNSPI使用不同的 QA波段生成的云去除图像；注:QA_Auto-PCP中的白色多边形代表了在原始QA波段中识别出的云污染的像素，黄色多边形代表了Auto-PCP识别出的剩余云像素。

图4为中国东南部测试区真实实验结果图。

图5为中国杭州测试区真实实验结果图。

图6为青藏高原测试区真实实验结果，每个图中的蓝色矩形突出明显的差异。

图7为四个测试区域的模拟云图像(标准假彩色合成)。

图8为实验二中模拟云图像去云结果的RMSE和SSIM值，其值为六个波段的平均值。

具体实施方式

在四个区域测试了Auto-PCP。首先是在华北平原的两季农作物耕地，在冬小麦收获后，夏季玉米或大豆通常在六月中旬种植，地表反射率值变化迅速；第二个站点是中国东南部的常绿林区，属于亚热带季风气候，云污染频繁；第三个站点是中国杭州，这个城市正经历着快速的城市化，景观构成复杂并有一条河流穿过城市；第四个站点在青藏高原，处于高海拔地区并且具有积雪和较复杂的地形。

在四个区域收集了landsat8地表反射率图像。设计了两组实验，使用图像信息如表1所示。

表1.不同的实验中使用的Landsat-8地表反射率数据

(1)实验I

(2)实验II:方括号中为参考图像和云模拟图像之间的SSIM

在实验I中，对真实云图像进行了云去除。在每个测试区域，分别使用了一幅陆地卫星云图。通过使用三个不同的QA波段来比较云去除的性能:(1)原始的QA波段(简称“QA_original”)；(2)云和云阴影边缘周围扩充两个像素处理的QA波段(简称“QA_dilation”)；(3)Auto-PCP处理的QA波段(简称“QA_Auto-PCP”)，对最终的恢复结果进行了视觉评估。

在四个测试区域的真实云图像上比较了云去除的性能(图3-6)。在每个图中，分别展示了 Auto-PCP处理的QA波段和结合改进的邻域相似像素插值方法(MNSPI)生成的去云图像。结果表明，在QA_original波段中，一些云和云的阴影没有被识别出来(见每张图的放大图)。在青藏高原站点中，QA原始波段错误地将许多雪像元识别为云像元，但同时也存在云像元的漏检情况(图6)。因此，基于QA_original的云去除图像在空间上具有明显的不连续性，在图像上存在显而易见的“白色”和“黑色”块。通过膨胀对QA波段进行处理，虽然在一定程度上提高了云去除的性能(参见各图中的QA_dilation)，但QA_dilation仍然不能完全覆盖原始QA 波段中被忽略的云和阴影斑块。然而通过应用Auto-PCP，发现云去除的性能有了明显提高，这要归功于Auto-PCP对剩余云和阴影像元的有效识别。

在实验二中，基于模拟云图像进行定量评估。使用Fmask 3.3生成每幅模拟云图像的原始 QA波段，将采用两项指标对云恢复结果进行定量评价：一是均方根误差(RMSE)，它主要量化预测的反射率值与真实反射率值之间的差异。由于QA_original、QA_dilation和QA_Auto-PCP 重构的像素个数不同，因此对三个QA波段并集(即QA_original∪QA_dilation∪QA_Auto-PCP) 中的像素计算RMSE。第二个评价指标是结构相似度指数(SSIM)，主要量化预测图像与真实图像之间的整体图像结构相似度。

图7显示了四个测试区域的模拟云图像。对六个波段的恢复结果进行了定量评估,包括蓝光、绿光、红光、近红外、和两个短波红外波段(即地球资源观测卫星8中的波段2-7)。图8 显示了所有六个波段的平均RMSE和SSIM值，从中可至QA_Auto-PCP在所有站点中RMSE值最小，SSIM值最高。例如，在华北平原通过MNSPI去云方法的恢复结果中，QA_original、QA_dilation和QA_Auto-PCP的RMSE值分别为0.0368、0.0321和0.0152。

与QA_original相比，QA_dilation和QA_Auto-PCP中标记了更多的云污染像素，这可能增加了云错检误差的概率。为了帮助理解使用不同QA波段去除云的性能，进一步分析了用于检测云污染像素的漏检误差和错检误差，定义为：

结果表明，QA_Auto-PCP显著降低了云和阴影的漏检误差(如在华北平原漏检误差由8.7％降至0.6％)。与此同时，错检误差没有增加，甚至略有下降(表2)。相反，虽然QA_dilation也减少了漏检误差，但其错检误差有明显的增加。

表2.不同QA波段的云污染像素错检误差和漏检误差

。

Claims (2)

1.一种基于Landsat光学遥感图像的云识别方法，该方法包括：

步骤1：识别初始云像素

关于云图像和参考图像之间，云像素在蓝光波段的反射率变化值ρ_blue(change)为：

ρ_blue(change)＝max(ρ_blue(cloud)-ρ_blue(ref′))-(ρ_blue(cloud)-ρ_blue(ref′)) (1)

其中，

和

分别表示像素i在云图像和参考图像中蓝光波段的反射率值；ρ_blue(cloud)表示在云图像中像素蓝光波段的反射率值，ρ_blue(ref)表示参考图像中像素蓝光波段的反射率值，P_other表示云图像QA波段中除开标记为云和云阴影的其他像素；在公式(1)中max(ρ_blue(cloud)-ρ_blue(ref′))表示图像像素间反射率变化的最大值；

对于云像素，结合ρ_blue(cloud)值和ρ_blue(change)值，计算云指数：

结合ρ_blue(cloud)和ρ_blue(change)的归一化形式，CI的取值范围在-1到1之间，对于云像素来说CI值比非云像素大；基于以下两个标准在P_other像素中确定了初始云像素：

1)：初始云像素i的CI值

与P_C像素集的CI中值median(CI(P_C))的加权做比较，P_C表示云图像QA波段中标记为云的像素，表示如下：

其中，median()表示取中值；

2)：使用相对阈值来分别确定每类的初始云像素；利用k-means方法对参考图像进行聚类；分类的数目设置为4-5类；以其中的一类为例，在类A中初始云像素i的CI值表示为

必须满足以下条件：

其中mean(·)和std(·)分别是均值函数和标准差函数，CI(P_other)表示P_other的CI值；参数a为标准差倍数；

步骤2：识别初始的云阴影像素；提出云阴影指数CSI识别初始云阴影，CSI表示如下：

其中：

ρ_nir(change)＝(ρ_nir(ref′)-ρ_nir(cloud))-min(ρ_nir(ref′)-ρ_nir(cloud)) (7)

ρ_nir(cloud)表示云图像中像素在近红外波段的反射率值，ρ_nir(ref′)表示参考图像中像素在近红外波段的反射率值，

表示云图像中，属于P_other集合的像素i在近红外波段的反射率值，

表示参考图像中，属于P_other集合的像素i在近红外波段的反射率值，min(·)是最小值函数，CSI取值范围为-1至1，对于云阴影像素来说，它的值更小，越靠近-1；如果在P_other中的像素

满足以下两个条件则为初始云阴影像元：

1)：

与P_C像素集的CI中值median(CSI(P_C))的加权做比较，表示为：

2)：分别确定每类的初始云阴影像素；以其中一类为例，在类A中初始阴影像素i的CSI值

必须满足下列条件：

其中参数b是标准差倍数；

步骤3：匹配初始云和云阴影像素；

步骤3.1：云总是伴随着云阴影:利用云、云阴影和太阳之间的几何位置关系来进一步细化初始云和云阴影像素；设初始云像素i位于(x,y)，根据下面的方程对应的云阴影像素为(x′,y′)：

x′＝x-H_{cloud_i}×tanθ×sinφ (11)

y′＝y+H_{cloud_i}×tanθ×cosφ (12)

其中θ是太阳天顶角，φ是太阳方位角，这两者都由Landsat图像元数据文件提供，H_{cloud_i}表示位于图像(x,y)处此像素的云高；

步骤3.2：基于QA波段估算H_{cloud_i}范围，假设H_{cloud_i}在QA波段标记云像素的云高度范围内；为此，根据8邻域的连通性在QA波段生成云斑块和云阴影斑块；从云斑块和其对应的阴影斑块之间的水平距离计算云斑块j的云高度H_{cloud_pj},表示为：

步骤3.3：计算DIS_pj；

步骤a：将阴影斑块沿着云投影方向移动，当云与阴影斑块的重叠区域达到阴影斑块面积的三分之一则停止移动；此时移动的水平距离是假定为x₁，阴影斑块与云斑块重叠的边界记为曲线ab；

步骤b：将曲线ab与阴影斑块对应的外边界a′b′进行匹配，然后基于两条曲线相似的形状的特性，通过线性插值使两条曲线具有相同的长度，并计算它们之间的相关系数；如果相关系数大于等于0.9，就认定匹配成功，此时两条曲线之间的水平距离表示为x₂；最终，DIS_pj是x₁与x₂的水平距离之和；如果相关系数小于0.9，则计算QA波段里n个云斑的高度{H_{cloud_p1},..,H_{cloud_pn}}，则可能的云高范围H_{cloud_i}为满足公式13条件的任意值：

为了避免H_{cloud_i}受到异常值的影响，从{H_{cloud_p1},..,H_{cloud_pn}}中删除了最大1％的云高值；Auto-PCP允许初始云像素在可能的云高度范围内变化，从而通过使用不同的高度匹配所有初始云和云阴影像素；

若初始云阴影像素找到相匹配的云像素，则保留初始云阴影像素；如果初始云像素的投影像素被其它云覆盖；则从被覆盖的云像素开始投影并确定下一个投影像素，该过程一直持续到投影像素为云阴影像素或无云像素；如果最后投影的像素为云阴影像素，则保留该初始云像素，否则，初始云像素将被移除；此外，在图像边缘的一些初始云和阴影像素的投影位置可能在图像之外，保留这些初始的云和阴影像素；

步骤3.4：删除8邻域连接的像素集小于7个像素的小云斑块和阴影斑块。

2.如权利要求1所述的一种基于Landsat光学遥感图像的云识别方法，其特征在于，所述步骤1公式5中a参数值设为2，步骤2公式10中b参数值设为2。

Images