基于朴素贝叶斯方法的FY-4A/AGRI云检测模型

鄢俊洁， 郭雪星， 瞿建华， 韩 旻

(北京华云星地通科技有限公司，北京 100081)

0 引言

云在地球大气系统的辐射收支平衡中的作用是由云的光学性质决定的[1]。云层的遮挡会导致卫星影像出现噪音，从而使遥感反演的参数出现误差，干扰后续的图像分类或目标识别。因此，对云的准确检测是遥感产品精准反演的前提。

传统的云检测方法主要利用云的高反射率和云顶的低温等特点[2]，结合物理阈值法[3]或同态滤波去云法[4]等进行检测。物理阈值法包括国际卫星云气候计划(international satellite cloud climatology project，ISCCP)云检测方法[5-8]、高级甚高分辨率辐射计云、陆、海处理方案(AVHRR processing scheme over cloud, land and ocean，APOLLO)检测算法[9-11]、美国国家海洋和大气管理局云高级甚高分辨率辐射计(the NOAA cloud advanced very high reso-lution radiometer，CLAVR)云检测算法[12-13]和二氧化碳(CO2)薄片法[14-15]。ISCCP法基于可见光和红外窗区波段的数据，把每一个像元的观测辐射值与晴空辐射值比较，若两者的差大于晴空辐射值本身的变化时，判定该像元点是云点; APOLLO法利用高级甚高分辨率辐射计(advanced very high resolution radiometer, AVHRR)5个全分辨率探测通道资料，包括5个阈值检测因子; CLAVR法关注的是极区，利用光谱及空间变化特征检测云的存在; CO2薄片法主要针对全球高云覆盖，包括薄卷云，利用对CO2敏感的红外光谱辐射，把穿透性云从不透明的云和晴空中区分出来。阈值法简单易行，但是对阈值的敏感性高，而阈值的选择具有一定的主观性，同时阈值法以像元为主要处理单元，没有考虑遥感影像的结构信息[16]。同态滤波去云法将频率过滤与灰度变化结合起来，分离云与背景地物，最终从影像中去除云的影响。这种方法由于涉及到滤波器以及截止频率的选择，在滤波过程中有时会导致一些有用信息的损失。随着模式识别理论的发展应用，正演模式被引入到云检测研究，如大气辐射传输模式[17]，考虑了大气温湿和云微物理状况，能够正演模拟云顶亮温[18]。传统的云检测算法需要大量的先验知识，经过大量的人工判读和调整，才能实现全球范围的云判识。针对这一问题，机器学习方法被逐渐引入到云检测研究中[19-20]。机器学习方法利用训练样本进行特征训练，将先验知识和人工判识的经验提供给计算机，由计算机来综合相关制约因素，完成回归学习过程，同时使算法的适应性增强[21]。利用机器学习进行云检测时，检测效果通常取决于选取的特征描述能力的强弱，因此需要寻找一种好的描述特征[20]。

以往气象卫星光学载荷的云检测主要还是采用与美国地球观测系统/中分辨率成像光谱仪(earth observing satellites/moderate resolution imaging spectroradiometer, EOS/MODIS)载荷类似的方法，通过多光谱在云上的反射和辐射特性进行云的识别，需要大量人为干预。据此，本文基于风云四号A星(FY-4A)多通道扫描成像辐射计(advanced geosynchronous radiation imager，AGRI)数据，提出了一种基于朴素贝叶斯方法的全自动云检测模型，使用朴素贝叶斯算法作为核心结构，基于光学载荷基本云检测原理选择合适的红外通道作为特性分类器参数，同时针对不同的地表类型和不同月份分别分类训练构建，最终得到基于朴素贝叶斯算法的全自动云检测模型，实现光学载荷影像云像元的高效、精准识别。

1 云检测方法流程

基于朴素贝叶斯的云检测方法从类别上分属于统计模型，因而整体上分为模型训练和产品推断2个阶段(图1)。

图1 朴素贝叶斯云检测方法流程Fig.1 Methodology of naive bayes cloud detection algorithm

模型训练阶段的目标是生成云检测特征概率的查找表，具体步骤如下： ① 先执行预处理，主要内容包括初始化准备工作，按照文件时间确定模型的月度标志，读取L1b数据，按照下垫面归集云和晴空像元的位置； ② 基于第一步归集到各类下垫面的云和非云的像元，计算用于与云检测相关的特征，在本研究中使用了6类特征，分别为T11，Tstd，Emiss4，Tmax-T，TD_11-85和GeoColor，3.3节对各类特征进行了详细介绍。

产品推断阶段的目标是基于时次观测的L1b产品文件推断其云检测产品。与模型训练阶段相同，推断阶段同样执行了预处理操作，与之不同的是推断阶段依赖训练阶段的云检测特征查找表来判识当前观测的云像元，对每个像元的6种云属特征，查找与特征值对应的机率，按照条件概率公式推断云属概率值，并按照下垫面不同云属阈值得到最终的云产品。

2 创建数据集

FY-4气象卫星是新一代静止轨道定量遥感气象卫星，其载荷AGRI通道数由风云二号气象卫星的5个增加为14个(表1)，覆盖了可见光、近红外、短波红外、中波红外和长波红外等波段。星上辐射定标精度为0.5 K、灵敏度为0.2 K、可见光空间分辨率为0.5 km，能实现云、气溶胶、水汽、陆地表面特性、海洋水色等大气、陆地、海洋参量的高精度定量反演[22]。

表1 FY-4A AGRI载荷通道 [23]Tab.1 Channel setting of FY-4A AGRI

3 朴素贝叶斯云检测算法

3.1 算法原理

朴素贝叶斯是基于贝叶斯定理与特征条件独立假设的分类方法。对于给定的训练数据集，首先基于特征条件独立假设学习输入输出的联合概率分布。然后基于此模型对给定的输入利用贝叶斯定理求出后验概率最大的输出。

具体地，设云特征为n维向量集合X={x1,x2,x3,...,xn}，输出为云类别标记Y={0,1}。则云检测问题描述为：

(1)

基于特征条件独立假设，上式转化为：

(2)

式中：X(i)为向量集合X中的第i维；xi为X中第i维的取值。

本文利用极大似然估计计算学习联合概率分布P(X,Y)。假定训练样本集合容量为J，云属先验概率的极大似然估计为：

(3)

特征条件概率分布的极大似然估计为：

(4)

式中ajk为第j个样本的第k个可能取值。

由此，对给定像元的云属特征x，通过学习到的模型计算其后验概率分布，根据后验概率最大化准则判别其云属类型。另外应云检测产品的要求，本文对检测结果按照概率阈值划分了可能云和可能晴空2种中间类别。

3.2 地表覆盖分类

由于云在不同下垫面类型上的特征差异显著，本文将地球表面分为7类型(深海、浅水、陆地、积雪、北极、南极、沙漠)。对不同地表类型进行分类的目的是为了得到对晴空条件认识存在的系统性偏差，这种偏差因地表类型的不同而有很大差异。分类的输入数据来自FY-4气象卫星导航文件，图2显示了2019年1月和6月FY-4A全圆盘中下垫面类型的分布。

图2 全球地表覆盖分类

从图2中可以看出，地表类型变化中最大的驱动因素是积雪的变化，本文中的积雪覆盖使用近实时积雪覆盖产品(near-real-time ice and snow extent, NISE)，积雪覆盖产品包括由微波成像仪(SSM/I)生成的北半球和南半球海冰、积雪覆盖范围。积雪覆盖产品以南、北半球2个25 km方位角、等面积投影(EASE grid)方式提供。首先，将EASE grid投影转换为全球等经纬度投影，然后根据经纬度信息与FY-4A全圆盘4 km网格数据进行匹配，得到全圆盘范围的积雪覆盖信息。图3显示了2019年具体各月份积雪和海冰的空间覆盖变化，从图中可以看出，积雪和海冰的空间覆盖度随季节变化差异显著。

图3-1 2019年全球各月积雪覆盖

图3-2 2019年全球各月积雪覆盖

3.3 云检测分类器

提取良好的云属特征是构建云检测分类模型的关键步骤。本文在调研现有云检测算法的基础上，针对FY-4A卫星AGRI载荷的观测特性选取了6种云相关特征：

1)T11： 11 μm亮温是经由AGRI载荷观测到的11 μm谱段辐亮度计算得到的亮度温度xbt_11。因中、高层云在11 μm谱段的亮温相较于地面低很多，因此可以直接使用该特征检测云[24]。

2)Tstd： 11 μm亮温局地标准差是为以目标像元为中心，局部窗口范围内的所有像元亮温的标准差。使用该特征能够较为准确地区分中、高层云的边缘，薄云和晴空背景。

3)Tmax-T： 11 μm亮温局地最大值差被定义为以目标像元为中心局部窗口范围内亮温最大值与目标像元亮温的差。与Tstd类似，中、高层云边缘处的该特征值显著高于晴空。

4)TD_11-85： 11 μm与8.5 μm亮温差定义为目标像元11 μm亮温与8.5 μm亮温之差。由于冰和水的吸收峰值位于在窗区通道不同的波长[25]，有云像元处该特征值一般为正值，而晴空像元处该特征值则普遍小于0。

5)Emiss4： 伪4 μm发射率定义为：

(5)

式中： e为自然对数； c2为常数；xbt_11为11 μm亮温；xbt_4为4 μm亮温。图4表现了FY-4A卫星AGRI载荷2020年1月1日昼夜时刻云与晴空发射率随太阳天顶角的关系。分析图4可以发现，云的伪4 μm发射率随昼夜变化十分明显。当太阳天顶角小于85°时(图4(a))，云的伪4 μm发射率多数大于1.8，晴空则集中于1附近，可能云、可能晴空存在于两者的过渡区间。而当太阳天顶角大于85°时(图4(b))，云的伪4 μm发射率的下界降低到1以下，与晴空分布存在较大的重叠。因此，本研究选择计算伪4 um发射率2种特征Emiss4_Day和Emiss4_Night，分别在太阳天顶角小于85°和大于85°时参与后验概率计算。

图4 伪4 μm发射率随太阳天顶角的变化

6)GeoColor： GeoColor是一种由AGRI载荷多通道观测值合成的3通道彩色图像，具有良好的视觉表现力[26]。图5是一幅彩色合成示例图，其制作方法为： 在太阳天顶角小于90°时，采用0.46 μm、0.64 μm和0.86 μm反射率经过拉伸、变换得到白昼图像； 当太阳天顶角大于90°时，采用3.9 μm和10.4 μm亮温经由拉伸、组合、变换得到夜间图像；按照太阳天顶角计算晨昏线位置处的渐变掩模，合并昼夜图像得到最终的彩色合成图。云和晴空在GeoColor合成图中具有较为明显的视觉差别，具体地，云表现为白色，晴空则表现为下垫面的色彩。由于云在GeoColor中的色彩是RGB通道共同决定的结果，且昼夜较为均匀，本研究决定在最大估计时学习通道值特征的联合概率分布，并将其引入到朴素贝叶斯的框架中作为一种新的独立特征。>

图5 2020年1月1日10时 UTC FY-4A/AGRI GeoColor彩色合成图Fig.5 Color composite image of FY-4A/ AGRIGeoColor at 10: 00 UTC January 1, 2020

3.4 模型构建

研究使用2019年12个月的FY-4A卫星AGRI载荷L1产品，定位(Geolocation， GEO)产品和云检测(cloud mask，CLM)产品。其中，L1记录谱段观测值，GEO提供时间和地理相关信息，CLM作为云检测的真实标记。每个月选取20%的时次构成训练集，80%的时次构成测试集合，图6是带有集合时次数量的数据划分统计情况。

图6 数据划分统计Fig.6 Data partition statistical chart

图7为特征分布与云属概率曲线修正图。为了降低载荷观测特征随时间变化对云属判断产生的影响，更加准确地建立其云属关系，本文以月度作为模型求取的最小时间尺度。研究求取云检测单一特征月度模型的步骤描述如下：

1)在收集数据集中随机选取20%的时次作为样本集，按云属和非云属统计每个时次各下垫面上特征的频数直方图，如图7(a)。

2)根据步骤1)得到频数直方图计算特征的云属和非云属条件概率分布，随后按公式(4)求取特征的贝叶斯模型参数。

3)在统计样本特征的分布时，区间边缘的样本总数较小，晴空与云的比例易受噪声的干扰，云属概率值有时会表现异常。因此需要对步骤2)得到的模型参数进行手动的修正，图7(b)展示了修正后11 μm亮温特征对应的云属概率曲线。

3.5 云检测类别划分

在根据求解得到的模型检测产品时，首先需要判断像元的下垫面类型，然后计算算法定义的6类特征值，最后照式(2)融合多个特征的条件概率分布系数，由此可以得到该像元的云属概率值。

现有主流的云检测算法提供的是4类检测结果，晴空、可能晴空、可能云和云，如EOS/MODIS，FY3D/MERSI和NOAA/AVHRR云检测产品。为了生成类别产品，本研究选定0.5作为云与非云的概率阈值。为了得到较为准确的云和晴空类别，设定0.1的阈值划分晴空和可能晴空，0.9的阈值划分为可能云和云。

4 结果验证分析

4.1 FY-4A/AGRI数据云检测结果

以FY-4A业务CLM产品作为参照，本文采用分别针对云和晴空的召回率(probability of detection,POD)、误判率(false alarm ratio,FAR)并综合Kuipers评分(Kuiper’s skill score,KSS)来从多个角度反映算法的精度。假设参照为云的像元且实际被判别为云的像元数量为a，参照为云的像元但实际被判别为晴空的像元数量为b，参照为晴空的像元但实际被判别为云的像元数量为c，参照为晴空的像元且实际被判别为晴空的像元数量为d，算法对云判识的POD反映的是本文算法判别出的云占实际的云的比例，云的POD公式为：

(6)

云的FAR反映的是本文算法判别出的错误的云占判出的云的比例，公式为：

(7)

晴空的POD反映的是本文算法判别出的晴空占实际的云的比例，公式为：

(8)

晴空的FAR反映的是本文算法判别出的错误的晴空占实际的晴空的比例，公式为：

(9)

KSS是一种对错误分类敏感的补充性综合指标，常用于算法传统指标相近时更严格的评价，计算公式为：

(10)

2019年12个月测试集验证统计结果如表2。从表中可以看出，本文方法云检测的POD最高为6月份的98.9%，最低为9月份的90.4%。总体上，在2019年全年验证集中，云的平均POD为97%，平均FAR为6.3%，晴空的平均POD为89.0%，平均FAR为2.9%，方法的平均KSS为87.4%。

表2 2019年12个月朴素贝叶斯以业务CLM为真值的交叉比对结果

图8展示了测试集1月1日10时的云检测结果。整体而言，贝叶斯云检测算法与业务算法的结果大体一致。对比同时次的彩色合成图(图5)，细节上，贝叶斯算法在白天部分判断出的晴空多于业务算法，可能晴空少于业务算法，与合成图观察到的结果更一致。在耀斑区，贝叶斯算法出现了部分中间状态，业务算法表现更好。在澳大利亚中部，业务算法出现了错判，而贝叶斯算法表现更好。在夜晚部分，业务算法在中国南海海域判云较多，而贝叶斯算法可能云较多，更符合实际云图分布。

图8 2020年1月1日10时UTC云检测结果对比

4.2 模型分析

利用统计手段计算模型连续的参数值是本文算法区别于经典云检测方法的最大不同之处。根据不同的特征量和下垫面信息从模型中获得对应的有云概率值，并将不同特征量的概率值归一化到0～1之间获得云属概率并根据分类阈值获得云检测结果。下面针对T11(11 μm亮温)特征参数模型进行分析。

图9为1月和6月7类下垫面类型所构建模型的T11特征概率曲线。其中，横坐标AGRI载荷在11 μm谱段的观测亮温，纵坐标为经由式(2)计算得出的后验云属概率。11 μm通道为长波红外通道，对于此通道，大气组成成分对辐射的影响可以忽略，水汽只有微弱的吸收和再发射作用，卫星接收到的辐射主要是云和下垫面的发射辐射。根据普朗克定律，通常情况下高云比低云温度低，低云比下垫面温度低。但是由图9可知，不同下垫面和季节，云的T11特征概率差别很大，在使用单一阈值进行区别时就需要特别注意。当下垫面为南极、北极和积雪时，云和下垫面差别较小，温度很低的高厚云才容易被识别； 在不同季节，当下垫面为深海和浅海时，云与晴空的曲线变化十分陡峭，说明云的区分明显，T11在此下垫面比较适用； 由于沙漠和积雪的下垫面主要分布在北半球，在不同季节时，沙漠和积雪的云属概率差异很明显，在1月时出现明显的逆温现象，这是因为此时沙漠和积雪的温度有可能比云低。另外，综合不同季节的模型和检验结果，在南极下垫面的有效情况样本点过少，此下垫面的模型正确性待确认。

图9 1月和6月各下垫面T11特征概率曲线

观察图10左侧上部的耀斑区，可以发现本文算法与FY-4A CLM产品表现接近一致，都检测到了成片的云区。但本文算法在亮度反应较为强烈的海洋上会产生一定程度地误判，其原因是GeoColor合成图中耀斑处。观察左侧中部的薄云区，可以发现本文算法在薄云边缘处的判识精度仍有待提升。观察右侧中下部的海陆边界处，可以发现FY-4A CLM产品沿着边界线出现了较为明显的误判，而本文算法表现良好。同时，由于FY-4A CLM产品使用了12 μm亮温，所以当谱段观测出现缺失时会形成明显的条带。本文算法使用的谱段较少，且基于概率的条件独立假设在信息缺失条件下仍能够得出后验云属概率，因而可以有效避免条纹现象。另外，本文对12个月份分别构建模型，因而需要评估不同月度模型之间的差异对云检测精度的影响。

图10 2019年9月26日5时 UTC云检测局地结果对比

图11直观展示了不同月度模型对云检测精度的影响。在图像的上部，6月份模型检测结果出现了大量的错判晴空。这主要是不同月度模型的下垫面划分不同，出现晴空错判的下垫面是积雪，而6月份积雪下垫面上云的特征与1月份有较大不同。2个模型在其他区域上的表现接近一致。

图11 不同月度模型检测差异

为了研究模型的月度差异，研究在1月的测试数据上定量对比基于1月和6月标记数据得到的2组模型，得到如表3的百分比混淆矩阵。分析表格可以发现，1月模型的云判识结果要高出6月模型3.84百分点，晴空甚至高出7百分点； 而明显错判的比例也明显较低，例如把云判为晴空的比例低6.91百分点，把晴空判为云的比例低4.02百分点。说明使用同期月度模型推断的云属类别会更加准确。

表3 不同月度模型云检测混淆矩阵统计Tab.3 Model cloud detection confusion matrix statistics for different months (%)

5 结论

本文针对FY-4A/AGRI云检测问题，提出了一种基于朴素贝叶斯算法的全自动云检测模型，实现了FY-4A/AGRI影像云像元的高效、精准识别。与已有的云检测方法相比，该模型引入了朴素贝叶斯理论，将云检测的阈值判断转化为概率评价，即不局限于是否有云，而是以云出现的概率做评价，评价结果更加合理。同时，新算法中涉及的数据只有红外通道，较好地解决了传统算法同时使用红外通道和可见光数据造成的晨昏区云检测分割线的问题。此外，该模型还引入了新的合成图特征，学习联合概率，取得了较好的云检测效果。

总体而言，基于朴素贝叶斯算法的全自动云检测模型展现出良好的识别光学载荷影像云像元的潜力，具有一定的理论意义和应用价值。不足之处在于使用的训练集真值来自FY-4业务产品，模型精度依赖于产品数据的精度，因此未来可以引入测云卫星CALIPSO的数据作为真值进行模型学习训练。另外，现有方法没有考虑到特征随太阳高度角和耀斑角的变化，未来可以加入M估计的方法，以改善分类评价效果。我们还计划将该模型应用于风云三号E星光学载荷的云检测，为解决气象卫星晨昏轨道云检测提供新的方法。