一种全天相机云图的自动处理方法*

魏诗雅，邱 波，曹子皇，范晓东，刘园园，穆永欢

(1.河北工业大学，天津 300401；2.中国科学院国家天文台，北京 100101)

图1 云图类别区域划分Fig.1 Division of cloud image region

对云的观测与分析是研究大气各种特性的重要手段之一，也是天文观测中不可或缺的一步。比如在天文台选址工作中，只有那些全年云图中晴天多的地方才是合适的台址。天文观测中一般使用全天相机[1]对天空进行拍照，从而获得云图；然后对云图中的云彩占天空的面积比例(云量)进行估算。在云量估算过程中，按照30 m口径望远镜判读全天相机的方法，分别在天顶距44.7°和65°处画圆，在图1中用蓝色和绿色标注，将全天云量分为：

Clear：外圈以内无云；

Outer：内圈无云，外圈到内圈之间有云；

Inner：内圈有云；

Covered：外圈以内厚云覆盖超过50%。

在历史文献中也出现过一些自动处理云图的方法，但没有一种是根据30 m口径望远镜的原则进行自动计算的。这些自动处理云图的工作包括两方面：(1)云状的分类。通过提取云的颜色、纹理、位置等特征对卷云、层云、积云等云状进行识别，进而可对天气状况做出快速准确的判断与预测[2]；(2)云量的观测。目前，云量观测主要通过遥感成像或者红外成像仪器观测完成。文[3]通过对全天空成像仪重新选定红蓝比阈值进行云量计算，实现白天全天空云量的持续自动监测；文[4]基于全天空红外成像仪获取图像，利用统计晴空阈值对图像进行云像素识别和总云量计算；文[5]通过地基红外云仪，获取全天空红外辐射亮温图像，并利用阈值分割方式得出全天空云分布及云量信息，有效减少地面环境参数及太阳光照对云图的影响，能够全天实时运行；文[6]利用一种基于优化的神经网络分类程序和遗传算法将云图分为晴空和两个云类：薄云和不透明的云，该方法对薄云检测正确率只有61%；文[7]提出了基于MAP-MRF框架的薄云检测算法，将综合特征和空间信息考虑在内，在云检测方面有所提高，但是部分云图受光照影响较大。因此本文在实现云量自动计算的过程中，主要针对云图中存在月光影响的云检测和不同厚度云层分类进行了研究。

本文研究了如何进行云量的自动计算与云图分类。云图的处理流程如图2，首先针对多云和少云云图分别使用时间分割法和差分法去除云图中月亮影响的区域，然后对去除月亮影响区域后的多云云图进行二值化处理，将云与背景进行分割，并使用基于灰度值的聚类算法对少云云图的云的厚薄进行量化分类，接着分别对多云和少云云图计算总云量，最后依据30 m口径望远镜判读云图的方法对云图进行自动分类。

图2 云图处理流程图Fig.2 Cloud images processing flow chart

1 云图预处理

由于云图数据类型多变，需先对样本数据进行预处理，提高算法的鲁棒性和云量计算的精确度。云图预处理主要包括云图增强和云图的预分类两部分。

1.1 云图增强

为了提高处理的正确率，图像增强很重要。本文进行图像增强的方法为形态学中的高低帽法[8]。其中，高帽变换的特性主要是高通滤波，由此可突出云的灰度峰值，从而增强云图中云的边界信息；低帽变化的主要功能是检测云图中的低谷，强调了距离较近或范围较小的目标云的界限，使之不易被忽略。二者的结合可以使图像中的前景与背景在灰度方面被进一步拉开，既突出细节信息，又强调边缘信息，提高图像的对比度。

其中，f为原图；°为开运算；·为闭运算；b为预先确定的形态学算子，在本算法中b是半径为5像素点的圆盘。图3(a)为原图，(b)为增强后的云图，可以看出图3(b)中云的边界和小块薄云比图3(a)有明显增强。

图3 云图增强Fig.3 Cloud images enhancement

1.2 云图预分类

本文中的180张云图数据是由佳能550D相机以及佳能EF 8 mm f 2.8鱼眼镜头构成的全天相机连续两天内拍摄采集的。曝光时间为60 s，时间间隔为5 min保持恒定。由于不同云图中月亮对云的影响各不相同，少云云图中的云与月亮边界能被清晰地分割识别，而多云云图中的月亮部分或者完全被云覆盖，影响了月亮与云的区分。因此需将180张云图分为多云和少云两类，提高云图处理的精度。由图4中多云状态(a)和少云状态(b)对比可以看出，在100到150的灰度范围内图4(a)明显比图4(b)的数值高，在0到50的灰度范围内图4(a)的数值低于图4(b)。利用这一特点对所有云图进行分类。

图4 多云和少云云图直方图对比。(a)多云云图的原图和直方图；(b)少云云图的原图和直方图Fig.4 Histogram comparison of cloudy and less cloud images(a)shows the original cloudy image and its histogram；(b)shows the original less cloud images and its histogram

2 云图中月亮影响区域的去除

通过云图预处理，云的边界得到明显的增强。然而受到月光影响，月亮影响区域与云的灰度值极为相近，容易将月亮的影响区域误判为云，因此本文分别采用时间分割法和差分法对月亮影响区域进行去除处理。

2.1 时间分割法

时间分割法是利用某些运动的物体在相邻帧里灰度值变化不大的原理，在连续拍摄的图片中识别物体与其运动轨迹。尤其在云图中，月亮的移动速度远低于云的移动速度。研究中利用云图时间上的连续性，在云图中判断目标云图与其邻近30张云图中某点的灰度值大小。步骤如下：

(1)月亮区域点：判断目标云图中某点的灰度值与其邻近的30张云图中对应点的灰度值差值的绝对值是否小于某一固定阈值。当目标云图中该点灰度值较大且与邻近至少30张云图中对应点的灰度值差值的绝对值均小于某一分割阈值时，则该点可能为月亮中的一点，如(3)式：

其中，I(t)和I(t-1)为连续拍摄的云图；t为拍摄时间；D(x，y)为连续拍摄的云图之间的差分云图；T为经过观察后选取的分割阈值。当I(t)与I(t-1)相减的绝对值大于所选的T时，D(x，y)＝1，即该差分云图为月亮，反之D(x，y)＝0，该差分云图为背景天空。

(2)月亮边界点：检查该点的灰度值是否在目标云图该点位置局部范围内取得极大值。当该点灰度值取得极大值时，则该点是月亮中的一点，并将其灰度值置0；

(3)对目标云图中所有点依次进行以上处理，直到将月亮识别出来。

通过时间分割法处理后的云图中，月亮能得到很好的去除，但是由于月光影响使得月亮周围的灰度值在较大的区域仍然存在。云图二值化后，在少云云图中该高亮区域易被误判为云；而在多云云图中由于月亮周围被云覆盖，相应的月光影响被云抵消，云与背景天空的分割效果较好，因此该方法适合处理多云云图。

2.2 差分法

图像的减法运算也称为差分法[9]，常用于检测相同场景下图像之间的差异。将差分法应用于云图的处理，利用两幅云图相邻时间段内月亮位置的不变性，将月亮及月光影响区域去除。

通过对全天相机云图的观察，月亮在连续40 min内位置的变化可以忽略不计。因此首先利用云图采集的时间信息，从可用云图中筛选某天的0点到上午6点30分之间，两两间隔约为40 min的10幅有月无云的云图，称之为背景图，其它的180张云图为连续两天拍摄，拍摄时间为0点到上午6点30分，称之为原图。对原图和背景图进行预处理后，将相同时间段内的原图与背景图进行减运算，可以很好地去除月亮影响区域，并将云从背景天空中分割出来，如图5。从图5可以看出，当云图中月亮未被云覆盖且清晰可见时，利用差分法处理效果较好，因此该方法适合处理少云云图。

图5 差分法处理结果。(a)原图；(b)背景图；(c)去除月亮后的结果图Fig.5 Difference method for processing results(a)shows original image； (b)shows background image； (c)shows the result after removing moon

2.3 不同云图处理效果对比

针对少云和多云云图分别使用差分法和时间分割法，差分法不仅去除了月亮影响区域，也将云从背景天空中分割出来(如图6)；时间分割法处理多云云图，能很好地去除未被云遮盖的月亮影响区域(如图7)。

3 云量的计算

经过云图预处理和去除月亮影响区域后，有效地将云从背景天空中分割出来。因此接下来分别对多云和少云云图采用不同的方法进行云量计算。

3.1 少云云图云量的计算

对于少云云图，由于云图中不同区域云的厚度差别较大，因此研究基于灰度值利用K-Means[10]对云图中云的厚薄进行聚类分析。K-Means聚类方法的基本思想是先随机选定k个聚类中心，计算每个样本聚类中心到新的聚类中心的距离，迭代计算直到群集中心的位移距离小于给定值为止。然后根据云图像素间的灰度值等特征方差，对云图像素样本进行聚类。相比于其它的聚类算法，K-Means算法收敛速度快，聚类效果较优。本文根据聚类中心灰度值的大小将云分为4类，分别为天空、薄云、中层云和厚云，并通过调节迭代次数，观察云分类的效果。

图6 差分法去除月亮。(a)，(b)，(c)原图；(d)，(e)，(f)差分法处理少云云图Fig.6 Difference method to remove the moon(a)，(b)，(c)shows the original image；(d)，(e)，(f)shows the result of the difference method for processing less clouds

图7 时间分割法去除月亮。(a)，(b)，(c)原图；(d)，(e)，(f)时间分割法处理多云云图Fig.7 Time Division method to remove the moon(a)，(b)，(c)shows the original image；(d)，(e)，(f)shows the result of the difference method for processing less clouds

图8左边是原图，右边是聚类结果图，其中(a)，(b)，(c)，(d)4幅图的聚类中心的灰度值分别为36、71、0、104，对应薄云、中层云、天空、厚云。其中聚类中心为0的结果图中，白色代表背景天空，黑色代表有云的区域；其余的3幅云图白色为云，黑色为背景天空。通过分析云图中不同区域厚云、中层云和薄云的分布及云量的大小判断云图的分类。

图8 云层聚类结果Fig.8 Cloud Clustering Results

目前国际上主要针对卫星云图根据云检测得到云像元，将云像元分为不同的种类，计算每种像元的概率，从而得到总云量。云量计算的方法主要有ISCCP，CLAVR-1，CLAVR-X，MODIS等。其中CLAVR-X[11]是通过云检测的方法，将像元分为厚云、混合云、混合晴空及晴空4种情况。由于该方法分类比较细，因此研究中通过改进CLAVR-X方法中的参数，并利用厚云、中层云、薄云、天空4种情况云量的加权和，得到少云云图的总云量。其中天空的云量为0，如下式：

其中， Nthick＿cloud， Nmiddle＿cloud， Nthin＿cloud分别为厚云、 中层云、 薄云的像素个数； Ntotal为区域内总的像素个数；f(c)为区域内的总云量。经过多次试验，并且与人工观测的结果对比得到了中层云和薄云的权重系数分别是0.9和0.6。

3.2 多云云图的云量计算

对于多云云图，由于云图中大部分区域被厚云覆盖，云与背景的灰度值差别较大，中层云和薄云云量较少，因此将多云云图云层分为厚云和天空两类。利用自适应阈值法中的Otsu算法对去除月亮影响区域的多云云图进行二值化处理，将云与背景进行分割。

Otsu算法利用聚类的思想，把云图的灰度值按灰度级分成两类，使得两类间的灰度值差异最大，每类中的灰度差异最小，通过方差的计算寻找一个合适的阈值分割云图。该算法计算简单，不受云图亮度和对比度的影响，鲁棒性较高。

最后通过统计二值处理后的云图中所有的灰度值为1的像素点的个数，并根据(5)式计算总云量：

其中，Ncloudy为灰度值为1的像素点的个数；Ntotal为区域内总的像素个数；f(c)为区域内的总云量。

3.3 结果讨论

作为方法研究，针对180张云图首先进行人工分类以得到可对比的真实判别结果，其次使用本文提出的自动分类方法进行分类。最终以人工分类为准，得到每种分类结果的识别准确率如表1。

从表1可以看出，人工观测结果与自动识别的结果在Inner，Clear和Outer 3种类别存在一定的误差。由于受云图中星点和全天相机亚克力罩反光的影响，Outer和 Clear类容易被误判为 Inner类。其中有2张应为Outer类，被误判为Inner，如图9(a)；有25张应为Clear类，被误判为Inner，如图9(b)。

图9(a)，(b)左图为原图，右图为自动检测结果。图9(a)中内圈的白点为星点，图9(b)内圈的白色区域为全天相机亚克力罩上的反光。

表1 分类结果对比Table 1 Comparison of classification results

图9 云图结果分析Fig.9 Analysis of Cloud image results

综上可知，本文提出的自动处理云图的方法对于Inner，Covered类的识别率较高，分别为93.9%和88.46%。然而Inner，Outer和Clear类别的识别精度较低，容易受噪声的影响。但是总体的自动识别率为76.67%，基本能够达到自动处理分类的效果。在判读速度方面，人工判读180张云图大约需要30 min，而自动处理分类所需时间为2.1 min，较人工判读速度大大提高，随着云图张数的增多，自动处理的时间优势也更加明显。

4 结 论

本文提出的自动处理分类方法判读速度快，平均识别准确率较高。由于天文观测者在人工观测云量的过程中，对于Inner，Outer和Clear 3类没有给定一个量化的分类标准，因此自动处理分类结果相对于人工观测结果存在一定的误差。但是根据30 m口径望远镜判读全天相机的分类结果中，各类别的平均识别准确率为76.67%，基本能够满足自动处理分类的要求。

本文实现了云图中月亮影响区域去除和云量的量化计算分类，有效缓解了人工判读云图类别的低效性。但是，目前云图数据量有限，本文仅对两天内的180张云图进行处理研究。由于是方法研究，虽然一些参数是固定的，但是在数据足够的情况下，这些参数(如月亮检测阈值T)可采用自适应最大方差阈值法确定。云量计算式的权重系数可通过统计不同厚度的云的灰度值，选用灰度百分比作为一个基本参数，其它参数都与此参数进行关联确定。因此在接下来需要增加样本量且对云图中的不同噪声建立不同的模型，进行有效的识别去除。同时需要对强月光情况下Clear类被错分为Inner类的问题进行专门研究，以提高整体识别率。