基于OTSU 在MODIS 数据中的云检测应用∗

解本巨孙 岩于龙振

（1.青岛科技大学信息科学技术学院 青岛 266100）（2.青岛科技大学经济与管理学院 青岛 266100）

1 引言

知识的积累，推动人类对新领域的探索，不断扩展数据获取的渠道，人们通过发射卫星来得到更多的地球信息，以便更好地了解自然规律。MODIS是搭载在Terra 和Aqua 两颗卫星上的中分辨率成像光谱仪［1］，通过观察地球表面获取相关信息，传输到地面进而对数据反演利用［2］。然而地球表面超过60%的部分都被云层覆盖［3］，因此通过卫星取得的信息，就不得不面对云层带来的干扰［4］，所以云检测是充分利用卫星遥感信息的前提［5］，为地质观察、气象预报以及自然灾害预警［6~9］等都提供极大的帮助。

我国地处中低纬度地区，受季风洋流影响较大，尤其在沿海区域，雨量充沛，上空更是常年笼罩着大量的云，对遥感信息的利用和转化产生了极大阻力。目前，云检测的方法主要分为两大类［10］，一是比较常见的物理方法［11~12］，也称之为阈值法，如何全军等［13］提出了基于MODIS 数据的阈值的云检测识别，该方法简单易操作，但在不同季节不同地点，检测效果有很大出入，并且在薄云和低云效果很差。又如马浩等［14］在2018 年提出基于某一地区的多波段固定阈值云检测，在对冬季的黄渤海地区进行云检测能达到一个较好的效果，但是适用的时相性较差。另一类为模式识别方法，常见的有SVM统计学方法，人工神经网络等，如瞿建华等［15］在2019 年提出基于深度学习的神经网络云检测方法，该方法准确率高，普适性强，能够达到较好的效果，但该方法模型训练周期长，使用的数据量多，并且对于数据集的选择比较困难。本文在多波段阈值法研究的基础上，融合改进后的OTSU 模型，实现阈值的自动获取，同时也提高了检测精度。

2 资料处理

MODIS 共有36 个中等分辨率的光谱波段［16］，每一到两天对地球表面观测一次，获取陆地、海洋、汽溶胶、以及云层［17］等众多信息。MODIS完全免费公开，可以从美国国家地质调查局地球探索者数据库网站上获取，根据自己的实际需求，选择相应的产品类型。为了方便研究，本文采用MODIS 1B 级别产品，选择在35.0°~42.0°纬度和117.0°~124.6°经度之间的区域数据，地理位置处于山东半岛以及黄渤海区域，该地区四季分明，且下垫面种类丰富，具有典型的代表性，其中用到的波段信息如表1 所示。

表1 选用波段的相关信息

前四个波段，只需要计算其反射值，而第29 波段则需要计算每个像元的亮温。本文选取2019 年11 月7 日、2020 年1 月5 日、2020 年2 月2 日以及2020 年4 月6 日 的 数 据 为 例，分 别 以20201107、20200105、20200202 和20200406 作为图像在文中的编号，经过预处理后得到的影像如图1所示。

图1 预处理后的影像

3 云检测方法

云检测主要原理是根据云和下垫面如水、土地、雪、植被等在各个的波段中所产生的辐射率以及辐射亮度的不同，来进行区别检测。云通常在可见光的范围反射率较其他下垫面物体高很多，而在红热波段亮温却很低，充分利用云这一特性，就是云检测的关键。以往的多波段阈值云检测方法，通过人眼观察，确定相应段数据的阈值，再综合分析每一幅影像检测效果得到统一数值，虽然方法实现简单，但是在不同时间不同高度，云的辐射率和亮温是有很大差异的，阈值的固定会使云检测鲁棒性差，具有极大缺陷。

云检测的本质同样是从背景信息中提取目标，不同的是在遥感数据中背景物体很多，并且各波段数据分布情况不同，因此需要进行一些转换和处理。通过融合OTSU 思想，将遥感里的数据里的云看成前景，把非云看成背景，选择合适的波段，满足算法模型的需求，在多波段固定阈值法的基础上改进，从而达到阈值的自动确定，减少主观性误差，适应不同季节不同类型的数据信息，以此实现影像分割。

3.1 多波段阈值方法

云检测算法对阈值判断分析如下：

1）在可见光波段波段1 中，植被反射率最低，然后是水、陆地、冰雪，它们的反射率在0.2上下，且数值之间相接近，分布状况相似，可视为同类型数据，其与云的反射率有明显差距，并且随着云层的高度厚度，差距将更明显，因此将该波段反射率的阈值设为T1，当B1大于此阈值时判定为云。

2）在MOODIS 数据的26 波段中，由于陆地，植被、水等下垫面的辐射将会被水汽吸收，卫星能接收的辐射数据只有冰雪和云，而高云和冰雪的辐射率值差距很大，因此可以用来检测高云，将该波段反射率的阈值设为T2，当B26大于此阈值时判定为云。

3）云和雪在可见光范围内的反射性比较相似，而在第6 波段处，云对太阳辐射的吸收能力较弱，反射率较高，而雪却恰好相反。计算波段1 和波段6 的归一化指数V1=（B1-B6）/（B1+B6），雪、云以及水这一指数是正数，而植被土壤是负数，利用V1>0 这个特性将它们分离出来，同时归一化也可以消除一定的大气辐射以及仪器影响。在波段29中，其他物体相比冰雪的亮温较低，利用V2>C29将冰雪分离出来。在波段1 中，由于云的反射率明显高于其他下垫面很多，且水的辐射率不会超过0.2，所以可以设定B1>0.2将水去除，经过以上操作后样本中只剩云像元即可判断为云。由于归一化的这个设定可能会将海岸线、湖岸判定为云，因此再利用在波段8 中对云的出现很敏感的特性，设置B8>T3来进行对海岸线的误判进行修正。

3.2 动态阈值实现

OTSU 算法也被称为最大类间方差法，被誉为是图像分割中确定阈值的最佳算法［19~20］。其根据图片里像素点的灰度值不同，将图片分成两个部分，分别为前景（目标）和背景，当这两部分的类间方差越大时，则可说明这两部分的差别越大，以此进行图像分割。

对于一个图像I（a，b），其大小记为W*L，其灰度范围为0～M-1，记i 为某点的灰度值，把分离前景和背景的阈值记为T，将i小于阈值T的所有像素占图像的比例记为p0，平均灰度记为u0，大于阈值T的像素占图像全部像素点的比例为p1，记其平均灰度为u1，图像的总平均灰度记为u，类间方差记为V则有：

对于上述公式，当V 取得最大值时，能够使得分割出的目标和背景都与图像中心距离最远，对应的T便是得到的阈值。但由于MODIS数据量大，云的边界不明显，传统OTSU 算法对于这种情况分割后的结果不理想。除了均值外，还有灰度平均方差也能够反映灰度分布，边界附近的灰度值往往波动较大，而平均方差值更能反映其离散程度，因此这里采用平均方差来代替传统OTSU 中的均值［21~22］，改进后的阈值T应满足：

其中g2、g02、g12如下，Pi代表图像中灰度值i 出现的概率。

由于图像的灰度值是在0 到255 区间内的整数，在计算图像的类间方差时只需直接带入公式计算就可以，对于遥感数据里的反射率和亮温并不在这个整数值区间内，所以需要进行标准化处理，从而映射到相应区间，适应改进后的OTSU 模型，对应公式如下，其中min是样本中的最小值，max则为样本中最大值。

4 实现过程

云检测实现过程借助VS2016 开发工具，通过C#语言搭建云检测平台，使用IDL语言对遥感数据的进行处理，基于组件技术模型技术调用COM_IDL_connect 进行读取IDL 编译后的SAV 文件，引入ArcGIS 的插件MapControl 将处理后的数据转换为为图片显示在页面中，无需再通过GIS 查看地图，具体实现步骤如下。

1）首先将遥感数据文件通过C#程序读取，进行预处理，分别为辐射校正、几何校正和区域裁剪。

2）数据经过预处理后，将生成TIF 文件，使用MapControl 插件载入该文件即可看到原始影像图。使用IDL 读取该文件，其数据主要是各个波段中的像元的辐射率或亮温，通过HISTOGRAM 函数对这些数据进行密度的分析计算，均等分成256 个小区间，将标准化后的数据注入。

3）观察统计后的数据，在分成的256 个小区间中，会发现在数据的首端或者末尾存在极少量像元“独占”多个小区间情况，由于整个数据量很大，几万甚至几十万个，而区间大小是由数据的最值决定，这些点将会扩大区间的范围，从而影响结果的精度，因此可以将这些点视为噪点进行过滤，从而减少误差。经过多次实验分析，当前端或者尾部小区间里落入的像元个数连续少于或者等于3，记录这些像元，将其去除，会使得去云效果更好。采用设计的DIVIDI 函数进行噪点去除工作，具体内容为，第一步对数据首端遍历，直到遇到的小区间内元素的个数值大于3，则终止，并且统计噪点个数为N1；第二步对数据的尾部进行倒遍历，统计噪点个数为N2；第三步去掉数据前N1个，尾端的后N2个，再次将新数据通过HISTOGRAM 函数统计，重复上述步骤，经过多次迭代，直至N1和N2值均为0，此时为云检测使用数据。

4）通过上述操作就可以得到不同波段里每一个像元对应反射率或者亮温的分布情况，将数据输入到融合OTSU 算法模型的OTSU_Threshold 函数中，动态计算得出每个波段对应的阈值，结果如表2。

表2 各个图像阈值

5）将阈值输入到云检测模型中，利用其自带的掩膜工具分离云和其他下垫面物体数据，将生成的掩膜图像中云像元的值设为0，其他地方像元值不变，从而达到去云效果。

5 结果分析

5.1 检测结果

为了更好地对比和检验效果，将固定阈值［14］与改进后OTSU 模型产生的动态阈值带入云检测模型处理，影像如图2所示。

图2 云检测结果对比

共有四组不同日期的分图，其中每个分图包含的3 个小图从从左到右分别为经过预处理后的影像、固定阈值和OTSU 自动阈值的云识别结果。观察图像可以看出，固定阈值法在不同日期检测效果差异较大，其中在20200406 影像中漏云情况较大。固定阈值法中阈值确定是由人目视观察多幅影像的去云效果得出，无法覆盖全部类型云，导致其的桎梏性，对于边界区域的薄云、碎云漏判现象严重。

5.2 精度验证

本文从定量的角度进行评估检测结果，分别对两种方法的总体精度、精确率和召回率进行评价，其定义如表3所示。

表3 混淆矩阵

其中，正样本代表云像素点，负样本则代表非云点。总体精度计算方法如式（8），用来表示检测结果中判断正确云像元的比例，精确率计算如式（9），表示预测正样本占所有预测样本的比例，召回率计算如式（10），表示正样本中被预测出的比例。

对于检测后的图像进行计算，得出的结果如表4所示。

通过表4 可以看出，经过改良后OTSU 动态阈值法较传统的固定阈值法总体精度提升较大，对于漏判现象明显改善，结合之前影像效果图，可以看出同日期的低云薄云检测效果明显加强，尤其在云边界地区。其根据不同的影像中所包含云类型状况的不同计算出对应的阈值，更好地将云与非云区分出，在环渤海地区具有良好的适应性。

6 结语

该云检测方法对每一个波段阈值的确定都能根据不同的遥感数据自动计算出其最佳阈值，不需人为观察调整，能够快捷准确地得出阈值，且经过改进后的OTSU 算法，能够更好地将薄云区分开来，极大地提高了在黄渤海地区适用性，同时也满足了遥感数据反演需求，为后续的海冰识别，天气预警，环境监测等提供重要的基础。同时本文也存在一些不足，当遥感信息数据极大时，可能计算时间过长，并且模型只使用到五个波段的遥感信息，没有更大限度的利用遥感信息，且当出现极端的情况，即云像元占比极少的时，误差比较大，这些在日后的科研工作中，都将继续深入分析，不断完善。