多层感知机结合辐射传输模型的复杂陆地表面云检测

邓梦娇，徐 新，2，马盈盈，龚 威，金适宽，胡瑞敏

（1. 武汉科技大学计算机科学与技术学院，湖北武汉 430065；2. 武汉科技大学智能信息处理与实时工业系统重点实验室，湖北武汉 430065；3. 武汉大学测绘遥感信息工程国家重点实验室，湖北武汉 430072；4. 武汉大学计算机学院，湖北武汉 430072）

1 引言

随着遥感技术的发展，卫星遥感数据在测绘、环境、国土、海洋、林业、农业等领域得到广泛的应用. 遥感卫星对地观测时会受到云的遮挡，导致遥感影像上的信息造成不同程度的缺失，降低其数据质量，这会直接影响到大气参数反演、异常检测和分类以及农作物检测等许多传统的遥感任务［1～3］. 此外，云层类型的多样性和地表的复杂性也阻碍了遥感影像中云的高精度检测［4，5］. 因此，准确地将云和晴空进行分离对遥感数据的处理和应用有着十分重要的意义.

目前，研究人员提出了大量基于阈值的传统云检测算法. 基于阈值的方法具有精度高、结果稳定等优点，其大多利用云具有较高的反射率以及较低的温度的特性，在不同的波段设定固定的阈值实现云检测. 例如用于AVHRR 的云掩膜算法应用5 个阈值测试，对每个像素而言，只有当所有的测试都证明其为无云像素时，才将其识别为无云像素［6，7］. 然而，对于复杂地表和云类型，很难从任何波长中确定合适的阈值来准确检测云. 针对这一问题，研究人员开发了动态阈值、自适应阈值以及多波段组合阈值等云检测方法. 例如，Jedlovec 等人［4］提出了一种用于GOES-12 成像仪的双谱复合阈值技术，利用20 天合成的差分图像来表征时空变化的晴空阈值. 随着机器学习在解决二分类问题上的良好发展，研究人员也开始将其应用到云检测上. 如Wei 等人［8］提出一种将随机森林与超像素提取结合的云检测方法；Li等人［9］提出了一种基于弱监督深度学习的云检测方法. 然而，由于地表的复杂性以及云层类型的多样性，如何用正确的数据训练机器学习模型变得尤为重要. 为了解决这个问题，本文通过SBDART 辐射传输模型模拟不同地表类型下的反射率，将其作为多层感知机模型的训练样本.

风云三号D 星（FY-3D）是我国的第二代极轨气象卫星，于2017年12月发射，2018年11月交付.FY-3D星上搭载的中分辨率光谱成像仪（advanced Medium Resolution Spectral Imager，MERSI II）是主要的传感器之一，可以获得地表、土壤、植被、水体、积雪等光谱范围内的反射率数据以及大气中粒子的吸收和散射辐射信息［10］. 本文以MERSI II 传感器为例，针对其覆盖面积广、扫描角度宽和像元空间分辨率较低等特点，设计了一套结合辐射传输模型和多层感知机模型的云检测方案. 该方案以辐射传输模型模拟的卫星表观反射率为训练样本，将训练好的多层感知机模型用来检测MERSI II影像中的云，最后使用CALIOP/CALIPSO 的同时观测数据对其进行验证，并将其云检测正确率与MYD35 产品进行对比. 相对于真实样本，本研究使用模拟样本有以下好处：（1）无人为选择的影响；（2）可应用于任何传感器；（3）可在传感器发射之前测试其性能. 本文方法所采用的思路为通用思路，可以为其他新的卫星或机载数据在试验阶段设定一套处理方案.

2 相关工作

云和地物在不同的波段具有不同的光谱特征，根据该特征，研究者们展开了一系列基于阈值的云检测工作. 例如，Rossow 等人［11］提出用于国际人造卫星云气候学计划（ISCCP）的云检测方法；Irish 等人［12，13］提出了一种用于Landsat 7 的自动云覆盖评估方法；Ackerman等人［14］提出了用于MODIS的云检测方法. 然而，此类方法依赖特定的传感器，对于复杂的地表和云层类型，很难找到合适的阈值来检测云［15］，并且当高分辨率卫星传感器波段较少且无法提供更多的特征来区分云和地表时，算法识别云的准确率会下降. 针对以上问题，Sun 等人［15］利用6S（Second Simulation of the Satellite Signal in the Solar Spectrum）［16，17］模型，模拟了不同观测和大气条件下的地表反射率和表观反射率变化的关系，提出了基于先验地表反射率数据库的通用动态阈值云检测算法；Zhu 等人［1］提出一种适合Landsat 影像的云检测方法，该方法根据不同的场景自动设置阈值并且结合云对象识别从影像中识别出云.

随着机器学习的发展，研究人员开始尝试将机器学习应用到大气遥感领域，开始利用该技术进行云检测. 例如，Wang 等人［18］使用支持向量机进行云检测；Meng 等人［19］为恢复影像中因云遮挡而造成的信息缺失，以去除影像中的云，提出了一种基于稀疏字典学习的图像修复方法；Wei 等人［8］为识别复杂地表上的碎云和薄云，先使用随机森林方法得到初步的分类结果，考虑云的空间特征，采用SEEDS（Super-pixels Extracted via Energy-Driven Sampling）分割方法对初步分类结果进行再处理，得到最后的云识别结果. 深度学习是机器学习的一个分支，以卷积神经网络为代表的深度学习算法在计算机视觉、自然语言处理等领域得到十分广泛的应用［20～22］. 基于卷积神经网络强大的表征学习能力，研究人员也将其应用到云检测. 例如，Xie 等人［23］采用简单线性迭代聚类方法构建超像素，再使用卷积神经网络从而识别出云；Segal-Rozenhaimer等人［24］针对World-View-2（WV-2）和Sentinel-2（S-2）卫星提出了一种基于卷积神经网络的域自适应云检测方法. 全卷积网络（Fully Convolutional Networks，FCN）［25］将传统卷积神经网络中的全连接层转化为一个个卷积层，通过上采样恢复到与输入图像相同的尺寸，实现基于像素级的分类. 随着FCN的提出以及图像分割任务的发展，研究人员也将云检测当作语义分割问题来解决. 例如，Chai 等人［26］将云和云影检测作为一种语义分割问题，基于SegNet 提出了一种应用于Landsat 系列卫星的云和云影检测方法；Mateo-García 等人［27］基于FCN 模型并利用两种域适应变换策略分别实现从Landsat-8 卫星迁移到Proba-V卫星和从Proba-V迁移到Landsat-8的跨平台云检测；Li 等人［28］基于编码-解码结构提出了一种可用于多种传感器的多尺度卷积特征融合方法，该方法融合6 种不同尺度的特征，得到基于像素级的分类结果；Wielan 等人［29］提出了一种基于经典的Unet 结构的云和云影分割方法，并展示了不同的训练策略和光谱波段组合对分割性能的影响.Yang 等人［30］提出了一种CDnet（Cloud Detection neural network）云检测方法，该方法基于ResNet50 网络进行特征提取，运用一个特征金字塔模块（Feature Pyramid Module，FPM）提取不同尺度的上下文信息，再通过一个边界细化模块（Boundary Refinement，BR）来细化云的边界，从而实现云检测，然而，由于遥感影像尺寸较大，为降低计算一般将其裁剪为256×256，512×512 或1024×1024 大小的影像；为了支持更大尺度的空间特征，同时减少对大型卫星图像的计算需求，Luotamo等人［31］提出了一种两级联CNN模型组件结构依次处理欠采样和全分辨率图像，从而实现云的检测.

对于机器学习和深度学习方法，训练数据的数量以及训练样本的正确性是十分重要的，然而由于城镇用地、裸土等地表的复杂性，该地表类型上的云通常难以被检测和正确标注. 而辐射传输模型计算的样本能充分考虑到各种复杂的情况，可以避免人工选择样本的主观性. 因此，本研究采用另一种思路，通过SBDART辐射传输模型模拟卫星表观反射率，得到多层感知机模型的训练样本，将训练好的多层感知机检测MERSI II影像中的云，用CALIOP/CALIPSO 的同时观测数据对其进行验证，并从月份、纬度两个方面将其云检测正确率与MYD35 产品进行对比. 为了进一步验证本文方法的有效性，将模拟得到的训练样本分别用于k-最近邻、朴素贝叶斯以及随机森林算法，并从总正确率和F1 评分两个指标上与本文多层感知机算法进行对比.

3 方法

3.1 表观反射率的模拟

SBDART（Santa Barbara DIS ORT Atmospheric Radiative Transfer）是一个可以用于计算在晴朗和有云条件下地球大气和地面间平面平行辐射传输的软件工具，可用于分析卫星遥感和大气能量收支研究中遇到的各种辐射传输问题［32］. 在耦合大气和地表的系统中，辐射传播可以通过求解下列辐射传输方程得到［33］：

其中，I(τ，μ，Φ)为散射辐射的分布，τ为光学厚度，ω和Φ分别为单次散射反照率以及散射角，p(cosθ)为散射相函数，μ和μ'分别为入射光和散射光的方位角，S*(τ，μ'，Φ')为单次散射源项. 因此，通过SBDART 来模拟各种复杂的大气条件下卫星所接收到的表观反射率.

输入不同的内在光学特性（Inherent Optical Properties，IOPs）到SBDART模型，关键IOPs参数的设置如表1所示. 其中，大气IOPs参数的设置主要考虑不同特性的气溶胶和云. 对于云参数，使用SBDART 模型自带的球形云滴参数数据库，主要考虑冰云和水云的高度、云光学厚度（Cloud Optical Depth，COD）和云滴有效半径.而表面IOPs 主要考虑陆地类型，陆地地表的反射特性由Ross-Li核驱动的BRDF模型确定. 其中，地表下垫面各向异性反射由聚类方法从MODIS 产品中得到，每一种地物类型聚类得到两套反射参数，并在辐射传输模型中将其以随机的线性比例混合，以模仿不同季节可能的反射特征.

表1 关键IOPs参数设置

通过SBDART 模型的模拟，得到约560 万个带标签的样本，每条样本除有关键IOPs 参数特征外，还含有MERSI II 传感器1～7 波段以及19 波段的反射率特征.将该数据作为多层感知机模型的训练数据.

3.2 多层感知机

多层感知机（Multilayer Perceptron，MLP）是一种将一组输入数据映射到一组输出数据的前馈神经网络［34］，其一重要特点是多层，因此也被称为深度神经网络（Deep Neural Networks，DNN）. 云检测可看作一种二分类问题，多层感知机能够用来解决复杂的分类问题，并且具有准确率高的特点.

多层感知机由3 层组成，分别为输入层、隐藏层以及输出层. 输入层接收数据的输入，隐藏层对输入数据特征进行处理，而输出层输出预测结果. 多层感知机的输出描述如下［35］：

其中，xi为给定样本的i个特征的输入，wij和wjp分别为输入层与隐藏层的连接权重以及隐藏层与输出层的连接权重，φH为隐藏层的激活函数，而φo为输出层的激活函数.

通过上一步辐射传输模型的模拟，得到约560万条带标签的样本，每条样本含有关键IOPs 参数特征、MERSI II 传感器1～7 波段以及19 波段的反射率特征.样本标签0 表示无云，1 表示有云. 本研究选取其中的13 个特征作为多层感知机模型的输入参数，分别为太阳天顶角、卫星天顶角、相对方位角、8 个波段的反射率、地表高程以及地表类型. 因此，多层感知机模型的输入神经元个数设置为13. 将多层感知机的隐藏层个数设置为2，其每层神经元的个数分别设置为10 和8，各层的激活函数选择ReLU函数. 将输出层的神经元个数设置为2，输出层使用sigmoid 函数. 由于在线性回归问题中，常选用均方误差（Mean Squared Error，MSE）函数作为损失函数，而在分类问题中，常选用交叉熵函数作为损失函数［36］，因此，本研究选用交叉熵损失函数来计算传播误差.

多层感知机模型最后一层的输出结果按照以下方式来确定：由于输出层通过sigmoid 函数将输出值映射到0 和1 之间，输出层两个神经元的值分别代表着无云和有云的概率，而每一层神经元从0 开始编号，即输出层两个神经元的编号分别为0 和1，与样本标签数值一致，那么可选取概率大的神经元编号作为整个模型的输出结果，即0表示无云，1表示有云.

将多层感知机模型训练的学习率设置为0.001，将batch_size 的大小设置为256，训练次数设置为200. 为了防止过拟合，在隐藏层的激活函数后面加上dropout层，dropout层随机丢弃神经元的概率分别设置为0.2 和0.5. 整个方法的流程图如图1所示.

图1 算法流程图

4 实验数据与结果

4.1 数据介绍

4.1.1 MERSI II/FY-3D 数据和MODIS 产品

MERSI II共配备25个通道，包括16个可见光-近红外通道、3个短波红外通道和6个中长波红外通道［37］，其扫描宽幅约为2900 km，并能每天完成一次对全球的覆盖. 本研究使用了MERSI II 传感器的前7个通道以及第19号通道，具体的波段及其基础参数如表2所示.

表2 MERSI II 传感器1-7及19波段基础参数

中分辨率成像光谱仪（Moderate Resolution Imaging Spectroradiometer，MODIS）是搭载在TERRA 和AQUA卫星上的主要传感器之一，两颗卫星相互配合每1～2天可以重复观测整个地球表面［38］，其观测到的数据可帮助我们理解全球陆地、海洋等的动态变化过程. 由于MODIS 和MERSI II 两种传感器的设计相似，因此本文主要使用MODIS 发布的MYD35，MCD43C1 以及MCD12C1 三种产品数据，分别作为云掩膜对比以及辐射传输模型的输入数据. 由于FY-3D 为下午星，因此选用同为下午星的AQUA 卫星上的MODIS 云掩膜产品MYD35，用于与本文方法所得结果进行对比. 其中，MYD35 为MODIS 二级产品，该产品算法采用一系列阈值测试来检测1 km 像素内是否存在云.MCD43C1 产品提供了光谱散射、体散射和几何光学散射三个权重参数，用于BRDF 模型以确定陆地地表的反射特性，本研究根据每一个地表类型从该产品中聚类出了两种地物反射参数，在模拟时随机线性混合这两种地表参数. 而MCD12C1 产品是Terra 和Aqua 卫星遥感获得的地表覆盖类型数据，用于确定不同经纬度下的土地覆盖类型，该参数作为辐射传输模型的输入数据之一.

4.1.2 CALIOP/CALIPSO 产品

CALIPSO（The Cloud-Aerosol LIDAR Infrared Pathfinder Satellite Observations）卫 星 是 由 美 国NASA 于2006年4月成功发射的太阳同步轨道卫星，分别于当地时间01∶30 和13∶30 过境，16 天可覆盖全球. CALIPSO由3 个部分构成，分别为双波长偏振激光雷达（Cloud-Aerosol Lidar with Orthogonal Polarization，CALIOP）、红外成像辐射计（Imaging Infrared Radiometer，IIR）和宽幅照相机（Wide Field Camera，WFC）. CALIOP 每24 h发射170 万次激光. 垂直特性掩膜（Vertical Feature Mask，VFM）产品数据由CALIOP 获得，采用一系列算法确定云和气溶胶的空间位置信息. 本研究使用4.0 版本、空间分辨率约为5 km 的VFM产品，最终得到近400万个样本，用作本文云检测算法和MODIS 云掩膜产品的验证数据.

4.1.3 实验数据

在辐射传输模型模拟阶段，利用MODIS 发布的MCD43C1和MCD12C1两种产品数据，得到辐射传输模型的输入数据. 在多层感知机模型训练阶段，使用辐射传输模型模拟得到的约560 万条带标签的样本数据.而在多层感知机模型的测试阶段，使用MERSI II 传感器获取的真实数据. 该数据的真实标签来源于CALIOP 的VFM 产品数据.CALIOP 是一种主动探测设备，发射激光脉冲来进行大气探测，主动信号对云滴颗粒很敏感，因此它的VFM 数据可以看作真实值，故可作为测试集的标签. 最终得到近400 万个样本作为测试数据.

本研究将辐射传输模型模拟得到的约560 万条数据作为训练集，将MERSI II传感器获取的近400万条数据作为测试集. 数据下载链接如表3所示.

表3 所用数据

4.2 云光谱特征

图2 是辐射传输模型在不同太阳天顶角范围内模拟得到的水云和冰云的反射光谱特征，每个图均为模拟的在不同的云光学厚度下表1 所列出的波段处的反射特征. 其中，图2（a）是太阳天顶角取值为0°～30°，图2（b）是太阳天顶角取值为30°～60°，图2（c）是太阳天顶角取值为60°～90°. 卫星天顶角和相对方位角分别为30°和60°. 由图2可知：第一，当云光学厚度和波长一定时，在不同的太阳天顶角范围下，云的反射率有所不同，并且随着太阳天顶角的增大，反射率先增大后减小，这可能是因为太阳天顶角过大，倾斜路径的散射会更强；第二，对于不同波长而言，有云和无云的反射率有明显的差异，在1.38 μm 波段处，由于该波段处于一个极强水汽吸收带［39］，水汽的强烈吸收导致地面和低空水云的辐射难以到达传感器，进而导致水云的反射信号非常微弱，而无云的反射率几乎为0；第三，随着云光学厚度的增大，反射率也随着增大，这是由于云越来越密实，云体反射率逐渐增大而造成的.

图2 不同角度下水云和冰云的反射光谱特征.由于1.38微米波段处主要受水汽吸收的影响，故将该波段单独列出放置最右侧

基于此，通过多层感知机来学习这些特征差异，进而检测出MERSI II影像中的云.

4.3 MERSI II光学影像云检测

本研究选取了不同的地区进行实验，其中包含了植被、冰雪、裸土等多种下垫面. 本研究选取了4 个案例，用训练好的多层感知机检测FY-3D 的MERSI II 影像中的云，识别结果如图3 所示. 其中，第1 列到第3 列分别为真彩色影像、云检测结果以及CALIOP 观测的VFM 产品. 真彩色影像中红色线为CALIPSO 过境轨道线. 图3（a）为我国西南部及印度北部地区的影像，该区域地貌复杂，横跨青藏高原；图3（b）是包含我国西部及土库曼斯坦邻近几个国家的影像，其植被较少，以温带荒漠和温带草原为主. 从这两幅影像的云检测结果可以看出，影像中的云几乎均被多层感知机检测出来. 图3（c）为蒙古东部地区及黑龙江省的影像，植被类型以森林为主，从其云检测结果可以看出，有部分地物被误判成了云. 图3（d）为我国中西部地区，其主要植被类型为草原和荒漠，从其云检测结果可以发现，部分云没有被检测出来. 从图3 可以看到，本研究方法能够较好地从影像中识别出云，尤其对于复杂的陆地表面，这可能是因为该地表类型特征易于被多层感知机捕捉.

图3 云检测结果图

4.4 基于CALIOP数据的对比实验

本研究利用CALIPSO 的VFM 产品对其进行验证，将VFM 产品得出的结果作为云检测结果，也就是说，如果VFM 产品划分其为云像元，则认为其为云像元，进而来计算本研究云检测方法的正确率.

4.4.1 多层感知机方法与MYD35产品的对比

利用训练好的多层感知机对MERSI II 影像进行云检测，计算了在不同月份和纬度下的正确率，并与MYD35 产品的正确率进行对比. 图4（a）为云检测正确率随月份变化的关系图，图4（b）为云检测正确率随纬度变化的关系图. 从图4（a）中可以看到，本文算法云检测结果有很明显的季节特征，除1 月、2 月、4 月和12 月之外，本文方法在其他月份下的云检测正确率均高于MODIS 产品，尤其是夏季.7 月份的云检测正确率能达到86.8%，而在冬季时的云检测效果最差. 这可能是因为冬季地表被冰雪覆盖，而冰雪与云的光谱特征非常相似，导致容易出现误判和漏判的情况. 增加亮温测试可能可以缓解冬季识别效果差的状况，但亮温信息很大程度上受下垫面而不是太阳光影响，复杂的地表类型使得亮温模拟更加困难. 此外，从图4（b）中可以看到，本文云检测方法在低纬度的云检测正确率较MODIS 产品高，在赤道附近的精度可达到91.7%. 这可能是因为低纬度地区，植被较为复杂，多层感知机能够较好地检测出复杂地表的云.

图4 正确率随月份和纬度变化的关系图

为了更好地验证本文的方法，计算了在不同地表类型下云检测的正确率. 同样地，将其与MYD35 产品正确率进行对比，结果如表4 所示. 表中加粗字体表示本文算法云检测正确率高于MYD35产品的地方.

从表4中可以看到，多层感知机在常绿阔叶林地表类型下的云检测正确率最高，而常绿阔叶林主要分布于低纬度地区，这正好与上文结论相符. 而地表类型为草地时，云检测正确率最低，这可能是因为该地表类型下的样本数较多，涵盖了可出现的多种情况，而多层感知机在训练阶段并没有学习到这些特征. 多层感知机在城镇用地、农田等复杂地表类型的正确率相差不大，并且均高于MYD35 产品的正确率，这可能归功于我们利用MODIS 产品MCD43C1 较为准确地计算了这些复杂地表类型下的地表各向异性反射. 然而，地表类型为浓密灌木丛时，本文方法云检测的正确率明显低于MYD35 产品，这可能是因为本文得到的该类型样本数较少，多层感知机没有完全学习或者过度学习了该类型下的特征.

表4 不同地表类型下的云检测正确率比较

4.4.2 多层感知机算法与k-最近邻、朴素贝叶斯、随机森林算法的对比

k-最近邻（K-Nearest Neighbours，KNN）［40］、朴素贝叶斯（Naive Bayesian，NB）［41］、随机森林（Random Forest，RF）［42～44］均为机器学习中高效的分类方法. 本文用该3种算法与本文多层感知机方法进行对比. 这3种算法所用的数据与多层感知机所用数据一致，即训练数据均为辐射传输模型模拟得到的模拟数据.

为了定量的评价这几种方法，本文采用正确率（Accuracy）和F1 评分这两种评价指标. 正确率表示的是正确识别的云和晴空像素的数量占总像素数量的比例.F1 评分为查准率（Precision）和查全率（Recall）的综合评价，用于综合反映整体的指标. 其中，查准率表示的是识别为云的像素中实际为云像素的比例，查全率表示的是识别为云的像素占真实云像素总数的比例.这两种评价指标的值越大，则算法的性能越好. 其计算公式如下：

其中，TP 为正确识别的云像素数量；FP 为误判为云的晴空像素数量；TN 为正确识别的晴空像素的数量；FN为漏判的云像素数量.

表5为多层感知机与KNN，NB，RF的对比结果. 表中在各评价指标上的最高值用加粗字体表示. 由表5可见，本文的多层感知机方法在正确率和F1评分这两个评价指标上是比较的方法中最高的，分别达到了76.25%和82.32%. 多层感知机的F1 评分值略高于KNN，NB，RF 3种算法，而这3种算法的F1评分值相差不大，这表明多层感知机算法对云像素具有更好的识别能力. RF算法云识别正确率仅69.51%，多层感知机正确率远高于随机森林算法，而随机森林算法F1评分却与其他对比算法相差不大，这表明随机森林算法不能较好地识别云像素，间接表明多层感知机更适合我们的方法.

表5 4种方法的对比

5 结语

本文结合多层感知机和辐射传输模型提出了一种云检测方法，该方法通过辐射传输模型模拟，得到了各种复杂陆地表面的反射率数据，以此作为多层感知机的训练样本，待训练完成后用以检测MERSI II 光学影像中的云像元，并用CALIPSO 的VFM 产品以及MODIS云掩膜产品对其进行对比验证. 本方法具有一定的通用性，并且云检测的结果较为理想，其总的正确率达到76.25%. 在不同的季节以及不同纬度下云检测正确率也有所不同. 通过与MYD35 产品对比发现，从季节上看，本文所提出的云检测方法在6月、7月和8月的效果最好，其中，7 月份的云检测正确率能达到86.8%. 从纬度上看，本文方法在低纬度地区的效果最好，其中，在赤道附近地表的准确率能达到91.7%. 本文计算了不同地表类型下的云检测正确率，与MYD35产品对比发现，本文方法在城市、农田和裸地等复杂下垫面条件下的云检测效果更好. 为了进一步验证本文方法的有效性，将同样的样本分别用于k-最近邻、朴素贝叶斯以及随机森林算法的训练，待训练完成后，将此3 种算法用于检测MERSI II影像中的云像元，并从总正确率和F1评分两个指标上与本文多层感知机算法进行对比. 结果表明，多层感知机结合辐射传输模型模拟具有更好的效果.

利用多层感知机可以捕捉到不同特征之间更深层次的联系，更为准确地检测出云. 在今后的研究中，可以考虑先将地表进行分类，确定不同地表类型下足够的样本数量，分别对不同地表类型进行训练，以进一步提高云检测的正确率.