一種全天相機云圖的自動處理方法*

魏詩雅，邱 波，曹子皇，范曉東，劉園園，穆永歡

(1.河北工業(yè)大學(xué)，天津 300401；2.中國科學(xué)院國家天文臺，北京 100101)

圖1 云圖類別區(qū)域劃分Fig.1 Division of cloud image region

對云的觀測與分析是研究大氣各種特性的重要手段之一，也是天文觀測中不可或缺的一步。比如在天文臺選址工作中，只有那些全年云圖中晴天多的地方才是合適的臺址。天文觀測中一般使用全天相機[1]對天空進(jìn)行拍照，從而獲得云圖；然后對云圖中的云彩占天空的面積比例(云量)進(jìn)行估算。在云量估算過程中，按照30 m口徑望遠(yuǎn)鏡判讀全天相機的方法，分別在天頂距44.7°和65°處畫圓，在圖1中用藍(lán)色和綠色標(biāo)注，將全天云量分為：

Clear：外圈以內(nèi)無云；

Outer：內(nèi)圈無云，外圈到內(nèi)圈之間有云；

Inner：內(nèi)圈有云；

Covered：外圈以內(nèi)厚云覆蓋超過50%。

在歷史文獻(xiàn)中也出現(xiàn)過一些自動處理云圖的方法，但沒有一種是根據(jù)30 m口徑望遠(yuǎn)鏡的原則進(jìn)行自動計算的。這些自動處理云圖的工作包括兩方面：(1)云狀的分類。通過提取云的顏色、紋理、位置等特征對卷云、層云、積云等云狀進(jìn)行識別，進(jìn)而可對天氣狀況做出快速準(zhǔn)確的判斷與預(yù)測[2]；(2)云量的觀測。目前，云量觀測主要通過遙感成像或者紅外成像儀器觀測完成。文[3]通過對全天空成像儀重新選定紅藍(lán)比閾值進(jìn)行云量計算，實現(xiàn)白天全天空云量的持續(xù)自動監(jiān)測；文[4]基于全天空紅外成像儀獲取圖像，利用統(tǒng)計晴空閾值對圖像進(jìn)行云像素識別和總云量計算；文[5]通過地基紅外云儀，獲取全天空紅外輻射亮溫圖像，并利用閾值分割方式得出全天空云分布及云量信息，有效減少地面環(huán)境參數(shù)及太陽光照對云圖的影響，能夠全天實時運行；文[6]利用一種基于優(yōu)化的神經(jīng)網(wǎng)絡(luò)分類程序和遺傳算法將云圖分為晴空和兩個云類：薄云和不透明的云，該方法對薄云檢測正確率只有61%；文[7]提出了基于MAP-MRF框架的薄云檢測算法，將綜合特征和空間信息考慮在內(nèi)，在云檢測方面有所提高，但是部分云圖受光照影響較大。因此本文在實現(xiàn)云量自動計算的過程中，主要針對云圖中存在月光影響的云檢測和不同厚度云層分類進(jìn)行了研究。

本文研究了如何進(jìn)行云量的自動計算與云圖分類。云圖的處理流程如圖2，首先針對多云和少云云圖分別使用時間分割法和差分法去除云圖中月亮影響的區(qū)域，然后對去除月亮影響區(qū)域后的多云云圖進(jìn)行二值化處理，將云與背景進(jìn)行分割，并使用基于灰度值的聚類算法對少云云圖的云的厚薄進(jìn)行量化分類，接著分別對多云和少云云圖計算總云量，最后依據(jù)30 m口徑望遠(yuǎn)鏡判讀云圖的方法對云圖進(jìn)行自動分類。

圖2 云圖處理流程圖Fig.2 Cloud images processing flow chart

1 云圖預(yù)處理

由于云圖數(shù)據(jù)類型多變，需先對樣本數(shù)據(jù)進(jìn)行預(yù)處理，提高算法的魯棒性和云量計算的精確度。云圖預(yù)處理主要包括云圖增強和云圖的預(yù)分類兩部分。

1.1 云圖增強

為了提高處理的正確率，圖像增強很重要。本文進(jìn)行圖像增強的方法為形態(tài)學(xué)中的高低帽法[8]。其中，高帽變換的特性主要是高通濾波，由此可突出云的灰度峰值，從而增強云圖中云的邊界信息；低帽變化的主要功能是檢測云圖中的低谷，強調(diào)了距離較近或范圍較小的目標(biāo)云的界限，使之不易被忽略。二者的結(jié)合可以使圖像中的前景與背景在灰度方面被進(jìn)一步拉開，既突出細(xì)節(jié)信息，又強調(diào)邊緣信息，提高圖像的對比度。

其中，f為原圖；°為開運算；·為閉運算；b為預(yù)先確定的形態(tài)學(xué)算子，在本算法中b是半徑為5像素點的圓盤。圖3(a)為原圖，(b)為增強后的云圖，可以看出圖3(b)中云的邊界和小塊薄云比圖3(a)有明顯增強。

圖3 云圖增強Fig.3 Cloud images enhancement

1.2 云圖預(yù)分類

本文中的180張云圖數(shù)據(jù)是由佳能550D相機以及佳能EF 8 mm f 2.8魚眼鏡頭構(gòu)成的全天相機連續(xù)兩天內(nèi)拍攝采集的。曝光時間為60 s，時間間隔為5 min保持恒定。由于不同云圖中月亮對云的影響各不相同，少云云圖中的云與月亮邊界能被清晰地分割識別，而多云云圖中的月亮部分或者完全被云覆蓋，影響了月亮與云的區(qū)分。因此需將180張云圖分為多云和少云兩類，提高云圖處理的精度。由圖4中多云狀態(tài)(a)和少云狀態(tài)(b)對比可以看出，在100到150的灰度范圍內(nèi)圖4(a)明顯比圖4(b)的數(shù)值高，在0到50的灰度范圍內(nèi)圖4(a)的數(shù)值低于圖4(b)。利用這一特點對所有云圖進(jìn)行分類。

圖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 云圖中月亮影響區(qū)域的去除

通過云圖預(yù)處理，云的邊界得到明顯的增強。然而受到月光影響，月亮影響區(qū)域與云的灰度值極為相近，容易將月亮的影響區(qū)域誤判為云，因此本文分別采用時間分割法和差分法對月亮影響區(qū)域進(jìn)行去除處理。

2.1 時間分割法

時間分割法是利用某些運動的物體在相鄰幀里灰度值變化不大的原理，在連續(xù)拍攝的圖片中識別物體與其運動軌跡。尤其在云圖中，月亮的移動速度遠(yuǎn)低于云的移動速度。研究中利用云圖時間上的連續(xù)性，在云圖中判斷目標(biāo)云圖與其鄰近30張云圖中某點的灰度值大小。步驟如下：

(1)月亮區(qū)域點：判斷目標(biāo)云圖中某點的灰度值與其鄰近的30張云圖中對應(yīng)點的灰度值差值的絕對值是否小于某一固定閾值。當(dāng)目標(biāo)云圖中該點灰度值較大且與鄰近至少30張云圖中對應(yīng)點的灰度值差值的絕對值均小于某一分割閾值時，則該點可能為月亮中的一點，如(3)式：

其中，I(t)和I(t-1)為連續(xù)拍攝的云圖；t為拍攝時間；D(x，y)為連續(xù)拍攝的云圖之間的差分云圖；T為經(jīng)過觀察后選取的分割閾值。當(dāng)I(t)與I(t-1)相減的絕對值大于所選的T時，D(x，y)＝1，即該差分云圖為月亮，反之D(x，y)＝0，該差分云圖為背景天空。

(2)月亮邊界點：檢查該點的灰度值是否在目標(biāo)云圖該點位置局部范圍內(nèi)取得極大值。當(dāng)該點灰度值取得極大值時，則該點是月亮中的一點，并將其灰度值置0；

(3)對目標(biāo)云圖中所有點依次進(jìn)行以上處理，直到將月亮識別出來。

通過時間分割法處理后的云圖中，月亮能得到很好的去除，但是由于月光影響使得月亮周圍的灰度值在較大的區(qū)域仍然存在。云圖二值化后，在少云云圖中該高亮區(qū)域易被誤判為云；而在多云云圖中由于月亮周圍被云覆蓋，相應(yīng)的月光影響被云抵消，云與背景天空的分割效果較好，因此該方法適合處理多云云圖。

2.2 差分法

圖像的減法運算也稱為差分法[9]，常用于檢測相同場景下圖像之間的差異。將差分法應(yīng)用于云圖的處理，利用兩幅云圖相鄰時間段內(nèi)月亮位置的不變性，將月亮及月光影響區(qū)域去除。

通過對全天相機云圖的觀察，月亮在連續(xù)40 min內(nèi)位置的變化可以忽略不計。因此首先利用云圖采集的時間信息，從可用云圖中篩選某天的0點到上午6點30分之間，兩兩間隔約為40 min的10幅有月無云的云圖，稱之為背景圖，其它的180張云圖為連續(xù)兩天拍攝，拍攝時間為0點到上午6點30分，稱之為原圖。對原圖和背景圖進(jìn)行預(yù)處理后，將相同時間段內(nèi)的原圖與背景圖進(jìn)行減運算，可以很好地去除月亮影響區(qū)域，并將云從背景天空中分割出來，如圖5。從圖5可以看出，當(dāng)云圖中月亮未被云覆蓋且清晰可見時，利用差分法處理效果較好，因此該方法適合處理少云云圖。

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

2.3 不同云圖處理效果對比

針對少云和多云云圖分別使用差分法和時間分割法，差分法不僅去除了月亮影響區(qū)域，也將云從背景天空中分割出來(如圖6)；時間分割法處理多云云圖，能很好地去除未被云遮蓋的月亮影響區(qū)域(如圖7)。

3 云量的計算

經(jīng)過云圖預(yù)處理和去除月亮影響區(qū)域后，有效地將云從背景天空中分割出來。因此接下來分別對多云和少云云圖采用不同的方法進(jìn)行云量計算。

3.1 少云云圖云量的計算

對于少云云圖，由于云圖中不同區(qū)域云的厚度差別較大，因此研究基于灰度值利用K-Means[10]對云圖中云的厚薄進(jìn)行聚類分析。K-Means聚類方法的基本思想是先隨機選定k個聚類中心，計算每個樣本聚類中心到新的聚類中心的距離，迭代計算直到群集中心的位移距離小于給定值為止。然后根據(jù)云圖像素間的灰度值等特征方差，對云圖像素樣本進(jìn)行聚類。相比于其它的聚類算法，K-Means算法收斂速度快，聚類效果較優(yōu)。本文根據(jù)聚類中心灰度值的大小將云分為4類，分別為天空、薄云、中層云和厚云，并通過調(diào)節(jié)迭代次數(shù)，觀察云分類的效果。

圖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左邊是原圖，右邊是聚類結(jié)果圖，其中(a)，(b)，(c)，(d)4幅圖的聚類中心的灰度值分別為36、71、0、104，對應(yīng)薄云、中層云、天空、厚云。其中聚類中心為0的結(jié)果圖中，白色代表背景天空，黑色代表有云的區(qū)域；其余的3幅云圖白色為云，黑色為背景天空。通過分析云圖中不同區(qū)域厚云、中層云和薄云的分布及云量的大小判斷云圖的分類。

圖8 云層聚類結(jié)果Fig.8 Cloud Clustering Results

目前國際上主要針對衛(wèi)星云圖根據(jù)云檢測得到云像元，將云像元分為不同的種類，計算每種像元的概率，從而得到總云量。云量計算的方法主要有ISCCP，CLAVR-1，CLAVR-X，MODIS等。其中CLAVR-X[11]是通過云檢測的方法，將像元分為厚云、混合云、混合晴空及晴空4種情況。由于該方法分類比較細(xì)，因此研究中通過改進(jìn)CLAVR-X方法中的參數(shù)，并利用厚云、中層云、薄云、天空4種情況云量的加權(quán)和，得到少云云圖的總云量。其中天空的云量為0，如下式：

其中， Nthick＿cloud， Nmiddle＿cloud， Nthin＿cloud分別為厚云、 中層云、 薄云的像素個數(shù)； Ntotal為區(qū)域內(nèi)總的像素個數(shù)；f(c)為區(qū)域內(nèi)的總云量。經(jīng)過多次試驗，并且與人工觀測的結(jié)果對比得到了中層云和薄云的權(quán)重系數(shù)分別是0.9和0.6。

3.2 多云云圖的云量計算

對于多云云圖，由于云圖中大部分區(qū)域被厚云覆蓋，云與背景的灰度值差別較大，中層云和薄云云量較少，因此將多云云圖云層分為厚云和天空兩類。利用自適應(yīng)閾值法中的Otsu算法對去除月亮影響區(qū)域的多云云圖進(jìn)行二值化處理，將云與背景進(jìn)行分割。

Otsu算法利用聚類的思想，把云圖的灰度值按灰度級分成兩類，使得兩類間的灰度值差異最大，每類中的灰度差異最小，通過方差的計算尋找一個合適的閾值分割云圖。該算法計算簡單，不受云圖亮度和對比度的影響，魯棒性較高。

最后通過統(tǒng)計二值處理后的云圖中所有的灰度值為1的像素點的個數(shù)，并根據(jù)(5)式計算總云量：

其中，Ncloudy為灰度值為1的像素點的個數(shù)；Ntotal為區(qū)域內(nèi)總的像素個數(shù)；f(c)為區(qū)域內(nèi)的總云量。

3.3 結(jié)果討論

作為方法研究，針對180張云圖首先進(jìn)行人工分類以得到可對比的真實判別結(jié)果，其次使用本文提出的自動分類方法進(jìn)行分類。最終以人工分類為準(zhǔn)，得到每種分類結(jié)果的識別準(zhǔn)確率如表1。

從表1可以看出，人工觀測結(jié)果與自動識別的結(jié)果在Inner，Clear和Outer 3種類別存在一定的誤差。由于受云圖中星點和全天相機亞克力罩反光的影響，Outer和 Clear類容易被誤判為 Inner類。其中有2張應(yīng)為Outer類，被誤判為Inner，如圖9(a)；有25張應(yīng)為Clear類，被誤判為Inner，如圖9(b)。

圖9(a)，(b)左圖為原圖，右圖為自動檢測結(jié)果。圖9(a)中內(nèi)圈的白點為星點，圖9(b)內(nèi)圈的白色區(qū)域為全天相機亞克力罩上的反光。

表1 分類結(jié)果對比Table 1 Comparison of classification results

圖9 云圖結(jié)果分析Fig.9 Analysis of Cloud image results

綜上可知，本文提出的自動處理云圖的方法對于Inner，Covered類的識別率較高，分別為93.9%和88.46%。然而Inner，Outer和Clear類別的識別精度較低，容易受噪聲的影響。但是總體的自動識別率為76.67%，基本能夠達(dá)到自動處理分類的效果。在判讀速度方面，人工判讀180張云圖大約需要30 min，而自動處理分類所需時間為2.1 min，較人工判讀速度大大提高，隨著云圖張數(shù)的增多，自動處理的時間優(yōu)勢也更加明顯。

4 結(jié) 論

本文提出的自動處理分類方法判讀速度快，平均識別準(zhǔn)確率較高。由于天文觀測者在人工觀測云量的過程中，對于Inner，Outer和Clear 3類沒有給定一個量化的分類標(biāo)準(zhǔn)，因此自動處理分類結(jié)果相對于人工觀測結(jié)果存在一定的誤差。但是根據(jù)30 m口徑望遠(yuǎn)鏡判讀全天相機的分類結(jié)果中，各類別的平均識別準(zhǔn)確率為76.67%，基本能夠滿足自動處理分類的要求。

本文實現(xiàn)了云圖中月亮影響區(qū)域去除和云量的量化計算分類，有效緩解了人工判讀云圖類別的低效性。但是，目前云圖數(shù)據(jù)量有限，本文僅對兩天內(nèi)的180張云圖進(jìn)行處理研究。由于是方法研究，雖然一些參數(shù)是固定的，但是在數(shù)據(jù)足夠的情況下，這些參數(shù)(如月亮檢測閾值T)可采用自適應(yīng)最大方差閾值法確定。云量計算式的權(quán)重系數(shù)可通過統(tǒng)計不同厚度的云的灰度值，選用灰度百分比作為一個基本參數(shù)，其它參數(shù)都與此參數(shù)進(jìn)行關(guān)聯(lián)確定。因此在接下來需要增加樣本量且對云圖中的不同噪聲建立不同的模型，進(jìn)行有效的識別去除。同時需要對強月光情況下Clear類被錯分為Inner類的問題進(jìn)行專門研究，以提高整體識別率。