金适宽，马盈盈，龚威，叶志伟，夏小鱼

1.武汉大学 测绘遥感信息工程国家重点实验室,武汉 430072;

2.武汉大学 电子信息学院,武汉 430072;

3.湖北工业大学 计算机学院,武汉 430068

1 引 言

云检测是卫星遥感影像预处理的关键步骤之一（周红妹等，1995；陈曦东等，2019；伟乐斯等，2021）。对光学遥感影像而言，由于云对短波辐射的显著影响，将云像元误认为地表或气溶胶将对许多传统遥感任务，如大气参数反演，土地利用变化、异常检测和分类，农作物监测和海洋生态调查等，产生负面影响（Zhu 和Woodcock，2012；张弛等，2021）。因此，可靠的云检测对陆地、海洋和冰冻层特性的卫星遥感研究至关重要。

利用云在可见光波段具有高反射的特点，在单一场景或少量的遥感影像中，云像元通常可以被人工识别和去除（刘成林和吴炳方，2004）。但对于日益增长的海量遥感资料，人工检测远远达不到应用需要。因此，针对这一问题，许多基于复杂阈值和逻辑测试的传统云识别方法被开发出来，例如用于陆地资源卫星Landsat ETM+传感器的自动云覆盖评估（Automated Cloud Cover Assessment）（Irish 等，2006），适用于MODIS 等中分辨率成像仪的MxD35 云检测算法（Ackerman 等，1998）和应用于AVHRR 辐射仪的扩展CLAVR−x 算法（Stowe 等，1999）等。这些传统算法通常依赖于特定波段反射或辐射信号的强度阈值和规则，通过辐射模拟、云场景统计和专家经验等手段，针对不同传感器精心设计。但这意味着这些方法在其他不同传感器中使用时具有较大的限制，并且由于强度阈值和规则是固定的，对气溶胶变化、季节循环和绝对地表反射率值也较为敏感（Segal−Rozenhaimer 等，2020）。在一些复杂的场景，例如地表阴影、高反射表面和浑浊水体中，传统方法会产生较大混淆（Hughes 和Hayes，2014；彭康龙 等，2021）。

为了改善卫星影像中复杂场景下云识别的效果，许多新颖的方法被开发出来。其中一些使用了动态阈值法，如Zhu 和Woodcock（2012）和Sun等（2016），前者根据不同场景自动设置阈值，并以云对象为单元识别Landsat和Sentinel卫星影像中的云；后者则在使用动态阈值的同时更多地考虑了地物反射特征，建立了一个先验的地表反射率库。此外，另一些方法充分利用了卫星重访周期之间获取的时间序列影像的特点，根据短时间内地物反射特性不会发生显著变化的假设，从反射信号中进一步分离出云像元（Li 等，2014；Lin等，2015）。近年来，随着机器学习和深度学习在分类问题上的优异表现越来越被重视，相关学者也开始尝试利用它们解决云识别问题，例如，Meng 等（2017）使用基于补丁修复的稀疏字典学习方法实现了厚云的检测和去除；Li 等（2019）设计了一种可用于多颗高分辨率传感器云识别的深度学习方法；Wei 等（2020）结合随机森林和超像素提取改善了Landsat 影像的云识别效果；Segal−Rozenhaimer 等（2020）使用卷积神经网络增强了可见光和短波红外影像晴空和云像素的区分程度。以上这些研究都给自动云识别算法和策略带来了新的解决思路。

风云3D（FY−3D）是中国研发的第二代极轨气象卫星之一（Yang 等，2019），其上搭载的先进中分辨率光谱成像仪MERSI Ⅱ（advanced Medium Resolution Spectral Imager）设置有与MODIS相似的光谱通道和空间分辨率，能完成对云，气溶胶，水汽，陆地表面特性和海洋水色的多种观测（Xu 等，2018）。起初，我们关注于MERSI Ⅱ对气溶胶的探测能力，然而由于仪器光谱特性的不同，适用于MODIS 气溶胶反演的云识别技术在MERSI Ⅱ上会产生大量晴空碎片，降低了卫星数据的使用率（Jin 等，2021）；对此，我们设计了一种基于机器学习和反射率模拟的方法——CRMC （Combine Reflectance simulation and Machine Learning for Cloud Detection）。该方法利用辐射传输模型模拟云在不同通道的反射率并作为样本训练机器学习算法。该方法从卫星影像中识别云像元，并利用CALIOP/CALIPSO 的同时观测数据进行验证。CRMC法最大的优点在于使用辐射传输模拟结果作为神经网络训练样本：在避免人为主观意识干预的同时能得到复杂下垫面条件下不同情况的云信号，并充分的挖掘了云在可见光和近红外波段的反射特性。

2 实验数据和方法

2.1 MERSI Ⅱ/FY-3D 数据

FY−3D 卫星于2017 年11 月发射升空，其上搭载的MERSI Ⅱ是主要的传感器之一。该传感器每天以45°角扫描地球表面并提供约2900 km覆盖范围的影像。在发射前，针对MERSI Ⅱ进行了详细的辐射校准以及空间和光谱表征能力测试，以确保数据产品的质量（Xu 等，2018）；反射率（reflectance）和亮温（brightness temperature）数据的不确定度分别约为3%和0.3 K。本研究使用2018年7—12月的MERSI Ⅱ数据，由风云赛（www.fysai.com/［2020−09−20］）提供。同时它们也可以从中国国家卫星气象中心的官方网站（http：//satellite.nsmc.org.cn/［2020−09−20］）获取。研究使用了MERSI Ⅱ传感器的前7个波段，具体的波段及其基础参数如表1所示。

表1 FY-3D MERSI Ⅱ传感器1—7波段基础参数Table 1 Spectral and spatial characteristics of FY-3D MERSI Ⅱ1—7 channels

2.2 MODIS产品

本研究主要使用MODIS 发布的3 种产品数据：MCD43C1、MCD12C1和MYD35。其中MCD43CI是由搭载在Terra和Aqua两颗卫星上MODIS传感器联合确定的地表二项反射函数BRDF（Bidirectional Reflectance Distribution Function）；MCD12C1 是基于国际地圈—生物圈计划（IGBP）的地表覆盖类型，用于模拟不同下垫面情况的反射率；而MYD35则是来自Aqua卫星上MODIS云掩膜的标准产品，用于对比验证。由于FY−3D 和Aqua 同为下午星，因此这里不再考虑来自上午星Terra 的云掩膜产品。

2.3 CALIOP/CALIPSO 产品

归功于主动遥感的优势，CALIOP 获得的廓线信息通常能更为准确的从不同大气环境中区分出云、气溶胶和洁净空气。这里使用4.0 版本、空间分辨率约为5 km 的垂直特性掩膜VFM（Vertical Feature Mask）产品作为云检测算法的验证数据。由于CALIPSO和Aqua同属于A−Train（The Afternoon Train）计划，因此它们的轨道相似，两颗卫星接近同步观测；而FY−3D和CALIPSO轨道区别较大，因此我们限制过境时间在30 min 之内的数据为同步观测数据。最终我们分别得到了约500万条样本用于对MODIS云掩膜产品和CRMC算法进行验证。

2.4 表面反射特性的确定

陆地地表反射特征由Ross−Li 核驱动的半经验BRDF模型确定（Roujean 等，1992），如式（1）：

式中，θi，θv，φ和λ分别表示太阳天顶角、卫星天顶角、相对方位角和波长；Kvol和Kgeo为体散射核和几何光学散射核；fiso，fvol和fgeo分别代表光谱散射、体散射和几何光学散射的权重，而最终地表反射率以ρs表示。为了确定不同下垫面条件下各种散射的权重参数，我们首先对2018 年全年的MODIS地表反射参数产品（MCD43C1）进行了系统抽样，以获取全球范围内时间和空间均匀分布的样本；其次，对这些权重参数进行了归一化处理，以突出光谱散射、体散射和几何光学散射的特征（焦子锑等，2011）；最后，我们对这些归一化之后的权重参数进行了聚类分析，获得了11 种不同类型下垫面的地表BRDF 反射权重参数，如图1 所示。其中，类型1—4 主要由稠密的植被构成，例如常绿阔叶林、常绿针叶林、热带雨林和热带稀树草原；类型5 和6 则是以开放的灌木和草地为主；而类型7 和8 的植被覆盖度则更低，主要是灌木、草地和裸土的混合；类型9 则全为裸土构成；此外，类型10 和11 分别表示陈雪和新雪，用于表示有雪覆盖时的下垫面反射特征。

图1 利用MODIS地表参数产品聚类得到的11种不同下垫面类型中Ross−Li核BRDF权重Fig.1 Ross−Li core BRDF weight parameters in 11 different underlying surface types calculated from MODIS product by using clustering method

2.5 表观反射率模拟

SBDART （Santa Barbara DISORT Atmospheric Radiative Transfer）是由美国加州大学圣巴巴拉分校发开的一种辐射传输模拟软件，可以计算在晴朗和多云条件下地球大气层和地表的平行辐射，主要用于分析卫星遥感和大气能量收支研究中遇到的各种辐射传输问题（Ricchiazzi 等，1998）。在耦合大气和地表的系统中，辐射的传播可以通过求解辐射传输方程得到（Stamnes 等，2011）：

式中，I(τ，μ，ϕ)为散射辐射的分布，ω为单次散射反照率，p（cosΘ）为散射相函数，Θ 为散射角，μ和μ'分别为入射光和散射光极坐标角的余弦，ϕ和ϕ'分别为入射光和散射光的方位角，τ为光学厚度，而S*(τ，μ'，ϕ')则为单次散射源项。在本研究中，根据输入的不同内在光学特性IOPs（Inherent Optical Properties），利用SBDART 模型计算出各种复杂情况下卫星接收到的理论表观反射率，进而作为机器学习算法的训练样本。

大气IOPs 参数的设置主要考虑不同特性的痕量气体、气溶胶和云。其中，痕量气体的设置考虑了3 种6S（Second Simulation of the Satellite Signal in the Solar Spectrum）大气辐射代码种定义的标准大气模式：热带、中纬度冬季和中纬度夏季。这些大气模式已广泛应用于各种大气研究中（Kotchenova 等，2008），并提供了气压、温度、水蒸气和臭氧密度的标准垂直廓线（McClatchey 等，1971）；而气溶胶模型则选用了城市型、农村型、大陆型和海洋型（Vermote 等，1997）。每一种模型都定义了一种典型的气溶胶光学特性，用于模拟不同情况下气溶胶对卫星信号的影响，其中，海洋型气溶胶被认为只在海面上出现；对于云参数，利用SBDART 模型自带的球形云滴参数数据库，考虑水云和冰云的高度、光学厚度COD（Cloud Optical Depth）和云滴有效半径，其中，水云高度被设置为地面以上2 km 以上的高空，而冰云高度则被设置为恒定的8 km。

表面IOPs 考虑陆地地表和海面两种类型。陆地地表反射特性由Ross−Li 核驱动的BRDF 模型决定，如2.4节中所述，根据MODIS地表参数产品划分出了不同植被覆盖程度下的9种地表和不同冰雪覆盖情况下的2种地表；其光谱散射、体散射和几何光学散射的权重用于直接输入SBDART 模型计算地表反射率，海面的二项反射则通过3个参数确定：海洋色素浓度，海面风速和海洋盐度（Vermote 等，1997），其中，色素浓度可以反映海洋富营养化程度，海面风速影响海洋泡沫和反辉区对反射率的贡献，而海洋盐度则影响菲涅尔系数。7个关键IOPs参数的具体范围如表2所示。

表2 关键IOPs参数设置Table 2 IOPs and their setting range

最后，将上述IOPs 输入SBDART 模型，并随机进行外部混合，得到了1000 万个带标签的理论样本。这些样本将作为训练数据集在下一步中对浅层神经网络进行训练。

2.6 浅层神经网络

浅层神经网络SNN（Shallow Neural Network）是指层数较少的人工神经网络，它具有计算速度较快的特点，能有效解决大部分复杂程度较低的分类问题（骆剑承等，2001）。本研究以一个具有一个隐含层和一个输出层的两层神经网络为例，训练辐射传输模型模拟得到的样本。该网络隐含层包含3个节点，激励函数使用正切函数；输出层包含一个节点，激励函数使用线性函数。对于需要学习大量数据的神经网络而言，简单的结构通常能在保证可靠性的同时增加运算速度。CRMC方法的详细流程图如图2所示。

图2 CRMC算法流程图Fig.2 Flow chart of CRMC algorithm

3 实验结果与讨论

3.1 云光谱特征

图3 显示了4 种典型下垫面条件下辐射传输模型模拟得到水云和冰云的光谱特征，其中，太阳天顶角和卫星天顶角为30°，相对方位角为60°。从图3可见，水云在可见光和近红外波段都会产生一个较高的反射，使得影像看起来更“白”；而冰云仅在可见光有强反射，且随着波长增加反射强度呈下降趋势。对于裸土或沙漠地表而言（图3（b）），冰云在近红外的吸收可能会降低光谱反射强度，但除此之外，冰云信号的强度和水云一样，都是远远大于地表的。相比利用经验判断，更希望机器学习算法能学习这些云特征。除此之外，在不同波段的光谱反射中，1.38 μm 处的反射对高层冰云具有较好的识别效，有助于将其与地表积雪区分开来。除了冰云外，水云和地表的信号在这个波段几乎十分微弱。这是因为该波段处于一个极强水汽吸收带（Gao 等，1993），导致来自地表和低空水云的反射信号被水汽吸收而无法到达卫星；相反，冰云通常处于较高的海拔高度，高海拔水汽含量低，水汽吸收作用低，在卫星影响上会产生一个强烈的反射。因此，我们利用辐射传输模型模拟的结果生成训练样本，训练构建的SNN，并尝试进行云识别。

图3 不同下垫面条件下的水云和冰云的反射光谱特征Fig.3 Spectral reflectance of water and ice clouds under different underlying surface type

3.2 MERSI II光学影像云识别

利用训练完的SNN 识别MERSI Ⅱ光学影像的云的结果如图4所示。其中，SNN 的输入参数分别为MERSI Ⅱ传感器前7 个波段表观反射率、太阳和观测几何、地表高程以及海陆掩膜。我们选择了4个案例来表明CRMC 方法在不同下垫面条件下的云识别能力，见图4。图4 中，从上至下第一个例子（20180428_0655）为蒙古国西部地区，海拔高度在1—2 km 之间，地表植被覆盖较为稀疏，裸土为主；第二个例子（20180610_0815）为印度半岛南部及印度洋，海拔高度在500 m 以下，植被覆盖茂密；最后两个例子（20180829_0455 和20181031_0455）为中国东部沿海地区和朝鲜半岛，陆地主要平原地区海拔在200 m 以下，植被覆盖较为茂密，海洋包括渤海湾，黄海和东海。受人类活动影响，渤海海水富营养化日益严重（于春艳等，2013），在光谱上显得更绿。

总体而言，CRMC 法展现了较好的对FY−3D MERSI Ⅱ影像的云识别能力。能识别包括高亮裸露地表在内的大部分云；然而，相对于在陆地的云识别能力，CRMC 在海面上效果较差。如图4 第3 个例子（20180829_0455） 所示，由于较高的反射率，图像左侧的海面反辉区（Sunglint）被CRMC 算法误分为云像元；另外，在长江入海口（32°N，122°E），由于浑浊的河水汇入了近海，海水混合大量泥沙导致光谱特征发生了较大的变化；而海陆掩膜数据判断该区域依旧是海面，因此，这种泥沙导致的光谱变化也容易使海水像元被误分为云像元。此外，我们发现CRMC 方法具有一定抵挡灰霾干扰的能力，如图4 中最后一个例子（20181031_0455），黄色虚线框内显示出现了明显的区域灰霾，但它被认为是云的概率仅接近0.1。

图4 利用CRMC识别FY−3D卫星MERSI Ⅱ光学影像中云的4个案例（图（b）真彩色影像中直线代表CALIPSO过境线，红色为CALIPSO识别的云区域，绿色为非云区域；最后一个例子中，黄色虚线框表示可能的灰霾区域）Fig.4 Four cases of cloud recognition in FY−3D satellite MERSI Ⅱoptical image by using CRMC method（（b）The straight line in the true color image represents the CALIPSO transit line，the red is the cloud area identified by CALIPSO，and the green is the cloud−free area.In the last example，the yellow dashed box indicates the possible haze area）

3.3 结合MODIS和CALIOP数据的对比验证

为了定量的研究CRMC 方法的准确性，我们将其和CALIPSO 的VFM 产品进行了对比。以VFM产品划分的结果作为真值，并且只要VFM 产品中标记了任意一层为云，本研究即认为该像元为云像元。然而，由于CALIPSO 分辨率较粗，并且过境时间和FY−3D 略有不同，导致数据匹配上会出现一定的系统误差。为了最小化误差带来的影响，我们还利用VFM 产品验证了MODIS 传感器MYD35的云掩膜产品，作为对比。实验范围包括亚洲中东部及临海区域。

为了找到一个合适的阈值去划分不同云概率条件下CRMC 方法的准确性，我们绘制了云识别正确率和阈值选择的关系，并认为某像元云概率小于阈值时为无云像元，而大于时为云像元。正确率由正确识别的像元数量占总像元数量的比值计算。从图5中可以看出，随着阈值不断增大，无云像元的正确率上升，而有云像元的正确率下降。利用选择不同的阈值，可以控制云像元筛选的严格程度，以适用不同的研究之中。

图5 云识别正确率随云概率阈值划分的变化关系Fig.5 Relationship between hit rate of CRMC method and threshold selection

基于CALIPSO VFM 数据集，我们分别计算了CRMC 方法在总体样本、陆地样本和海面样本中，云识别正确率和阈值变化的关系（图5）。随着阈值的变化，最大正确率在总体、陆地和海面样本中分别为79.6%、78.5%和81.2%；而当对云和非云像元的区分效果相同时，CRMC方法对云和非云像元的识别精度分别为78.3%和78.5%。由于地表覆盖和云的变化通常和气候有关，为了进一步验证实验效果，我们计算了不同月份和纬度条件下，利用CRMC 法对FY−3D MERSI Ⅱ影像云识别的正确率（阈值为0.2）和Aqua MODIS 云掩膜产品的正确率。从图6 中可以看到，CRMC 算法对陆地表面的云识别较为准确，在7 月达到最高85.9%；并且，除在1 月、2 月和4 月正确率略低于MODIS 产品外，其余月份正确率均较高。从纬度上来看，CRMC算法在中低纬度云识别精度较高，最高在赤道附近达到约91.5%；相反，在高纬度地区它的精度则较MODIS 产品低。这是因为夏季和低纬度地区植被生长通常更为茂密，而茂密的植被具有鲜明的光谱反射特性，更容易从中区分出云。然而，在海面上，CRMC 算法云识别精度在70%—80%之间，小于MODIS 产品，这可能是因为海面反射更加规则，固定的阈值法有利于识别反射特征不明显的薄云，并且MODIS 产品使用了多个亮温波段的劈窗测试（Split Window）。海面的亮温同样相对于地表而言更加有规律，例如冰冻阈值测试以270 K 为阈值，根据11 μm 波段的亮温识别热带海面的高层厚云（Ackerman 等，1998）。MERSI Ⅱ作为与MODIS 设计相似的传感器也具有6 个亮温波段，因此，在未来的研究中我们可以同样增加亮温通道的劈窗测试以提高云识别的正确率。

最后根据MODIS 地表类型气候学数据（MCD12C1），我们研究了当阈值设定为0.2 时CRMC 算法在不同下垫面条件下对云识别的正确率，并与MODIS 云掩膜（MYD35）进行了对比，如表3所示。其中，CRMC 算法优于MODIS云掩膜产品的地方用加粗字体表示。在当前阈值下，CRMC 在陆地上的云识别正确率为78.07%，略高于MODIS 官方云掩膜产品的77.88%；而在海面效果较差，正确率仅为78.75%，小于MODIS 官方产品82.93%的正确率。这里海面和陆地的划分使用的是MCD12C1 产品中定义的标准，和CALIPSO 的定义不完全一致，因此在结果上和图6显示出了一定的差距。此外，我们发现在农田、城镇用地和裸土等地表复杂的地区，CRMC算法的正确率高于MODIS 官方云产品。这一结果可能归功于我们利用地表二项分布产品较为准确的估算了这些复杂下垫面条件下的地表各向异性反射。根据以上研究，CRMC方法显示出了较好的云识别效果，并在未来有如下改进方向。

图6 CRMC法云识别正确率（阈值为0.2）和MYD35正确率随月份和纬度的变化Fig.6 Hit rates of CRMC（threshold of 0.2）and MYD35 change with variations of month and latitude

表3 不同下垫面下CRMC（阈值为0.2）与MYD35云识别正确率比较Table 3 Comparison of the hit rate of cloud recognition between CRMC（threshold of 0.2）and MYD35 under different underlying surfaces based on MCD12C1

（1）对于类似于FY−3D MERSI Ⅱ等具有亮温通道的传感器而言，尝试利用亮温信息并结合亮温通道劈窗法。

（2）在利用辐射传输模型模拟训练样本时，选用更精确的IOPs参数。

（3）选用深度学习方法代替本文使用的浅层神经网络，增加网络对云像元反射特征的学习并提高泛化能力。并且，由于云识别的准确性与下垫面类型相关，在训练样本和云识别时可以将下垫面类型作为一个额外的特征加入。

4 结 论

本文提出了一种从遥感卫星的光学影像中识别云的新方法CRMC。该方法结合了辐射传输模型和机器学习算法的优势，通过辐射传输模型模拟不同条件下的云和气溶胶理论表观反射率，并将其作为标记好的样本对浅层神经网络进行训练，以解决遥感影像中云识别这种分类问题。辐射传输模型计算的样本能充分考虑不同下垫面条件下、不同云类型和光学厚度的情况，避免人工选择样本的主观性。而利用机器学习算法能充分利用可见光和近红外波段的反射率特征，并根据后续实验选择合适的阈值，从而达到有效地从陆地地表识别云像元的目的。

通过与CALIPSO 卫星VFM 数据集的对比验证显示，CRMC方法的最大正确率在总体、陆地和海面样本中分别为79.6%、78.5%和81.2%。通过调整云阈值，它也能控制对云剔除的严格程度，生成适合于后续实验的云掩膜产品。当阈值选择为0.2 时，对于云和无云像元的总识别精度都达到约78.4%。通过进一步与MODIS云掩膜产品的横向对比发现，CRMC 方法对于海面的平均正确率比MODIS 云掩膜产品略低。这可能是因为我们仅考虑了可见光和近红外波段的反射率信息，而忽略了亮温的变化。对于表面特性较为均匀的海面，亮温信息对于云的识别尤为重要。而对于陆地云的识别正确率在春、夏、秋季和中低纬度地区较好，均高于MODIS 云掩膜产品。最高的正确率出现在赤道附近，为91.5%。此外在阔叶林、农田、城市地标和裸土等下垫面条件下，CRMC算法正确率均优于MODIS云产品，显示出了的较大潜力。

志 谢感谢MODIS 和CALIPSO 团队免费提供的用于本研究的地表特性、云掩膜和VFM 产品；感谢国家卫星气象中心和风云赛（http：//www.fysai.com/［2020-09-20］） 提供的本研究使用的FY-3D MERSI Ⅱ传感器数据；最后，感谢所有匿名评论人和编辑对本文提出的建设性意见。