基于模糊C均值聚類算法的去云處理

許 晨,徐 洋,康 雪,呂達(dá)仁

(1.成都市氣象局，四川 成都 610000；2.四川省氣象局，四川 成都 610000;3.中國科學(xué)院大氣物理研究所，北京 100029)

在遙感數(shù)據(jù)中，薄云覆蓋不僅會影響遙感數(shù)據(jù)的解譯精度，嚴(yán)重時還會導(dǎo)致遙感數(shù)據(jù)不能得到有效利用。因此，去云處理是遙感圖像改善技術(shù)的一個重要方面。本文重點(diǎn)探討云區(qū)數(shù)據(jù)聚類的有效性，通過對同一軌不同通道的云區(qū)進(jìn)行聚類處理，將云區(qū)與植被、土壤、水域等進(jìn)行分離，生成云區(qū)和非云區(qū)數(shù)據(jù)，再運(yùn)用中值濾波和云區(qū)替換的方法實(shí)現(xiàn)去云。

1 資料來源與方法

1.1 資料來源

選取NOAA-17衛(wèi)星的1、4、5通道，分別對應(yīng)可見光和兩個遠(yuǎn)紅外通道，如表1所示。其中，通道1獲得白天云圖及地表數(shù)據(jù)；通道4、5獲得晝夜云圖、海表、地表溫度和土壤濕度[1]。

表1 NOAA-17衛(wèi)星1、4、5通道的波長范圍

1.2 模糊C均值聚類

1.2.1 算法介紹

模糊C均值(FCM，F(xiàn)uzzy C-Means)算法是基于目標(biāo)函數(shù)的聚類法，每個數(shù)據(jù)點(diǎn)屬于某個聚類的程度用隸屬度來確定，利用均方逼近理論構(gòu)造帶約束的非線性規(guī)劃函數(shù)，該算法設(shè)計簡單、容易實(shí)現(xiàn)，能借助非線性理論解決優(yōu)化算法，應(yīng)用面廣[2-3]。

(1)隸屬度函數(shù)

一個對象X隸屬于集合A的程度用隸屬度函數(shù)uA(X)表示，0≤uA(X)≤1，其自變量范圍是所有可能屬于集合的對象(即集合所在空間中的所有點(diǎn))[4]，uA(X)=1代表X∈A，說明向量完全屬于集合A；uA(X)=0代表XA，說明X完全不屬于集合A。

(2)模糊集合

設(shè)X={x}是一個集合，uA(X)定義在X上，并且uA(X)∈[0,1]，則X上的一個模糊集合wA是指：對于任意x∈X，都確定了一個數(shù)uA(X)，稱uA(X)為X對wA的隸屬度。模糊集合完全由其隸屬函數(shù)刻畫，使一個元素隸屬的集合類別非硬性[5]。

(3)模糊集合的C劃分

把n個模式中的觀測樣本集X={x1,x2,…,xn}分為C個模糊組，求每個模糊組的聚類中心，使非相似性的價值函數(shù)最小，根據(jù)每個數(shù)據(jù)點(diǎn)的隸屬函數(shù)得到樣本屬于各個類別的不確定性程度[6]。

根據(jù)隸屬度定義可知，一個數(shù)據(jù)集的隸屬度之和為1，即

(1)

因此，F(xiàn)CM的目標(biāo)函數(shù)(價值函數(shù))表示為：

(2)

其中uij∈[0,1]，ci表示模糊組的第i個聚類中心，dij=‖ci-xj‖是ci與第j個數(shù)據(jù)點(diǎn)間的歐幾里德距離，m是加權(quán)指數(shù),m∈[1,∞)[7]。

接下來求使式(2)達(dá)到最小值的必要條件，構(gòu)造如下新目標(biāo)函數(shù)：

(3)

其中，式(1)的n個約束式的拉格朗日乘子表示為λj(j=1,2,…,n)。對所有輸入?yún)⒘壳髮?dǎo)后得到使式(3)達(dá)到最小值的必要條件為[8]：

(4)

(4)加權(quán)指數(shù)m

在FCM算法中，加權(quán)指數(shù)m非常重要，它控制隸屬度的分配和聚類模糊程度，m∈[1,∞)。m的選取以合理聚類為最終目標(biāo)，要求FCM算法同時滿足兩個條件：實(shí)現(xiàn)模糊聚類及進(jìn)一步將數(shù)據(jù)集劃分清晰，這樣才能準(zhǔn)確分辨樣本集中各對象的類屬關(guān)系。關(guān)于參數(shù)m的優(yōu)選方法可表示為：

(5)

Jm是聚類目標(biāo)函數(shù)，U、P是不同m取值下FCM算法得到的最優(yōu)劃分，將最佳加權(quán)指數(shù)m的確定轉(zhuǎn)化為一個帶約束的非線性規(guī)劃問題[9]。

1.2.2 FCM實(shí)現(xiàn)

(1)初始化[10]：設(shè)置聚類類別數(shù)c，2≤c≤n，n表示數(shù)據(jù)個數(shù)。設(shè)置迭代停止閾值ε，將聚類原型模式P(0)初始化，將迭代計數(shù)器設(shè)為b=0。

(2)用值在[0,1]的隨機(jī)數(shù)初始化劃分矩陣U，使其滿足式(1)的約束條件。

(3)用式(4)更新劃分矩陣U(b)。

(4)根據(jù)式(4)更新聚類原型模式矩陣P(0)。

(5)重復(fù)以上過程進(jìn)行迭代，若‖c(b)-c(b+1)‖<ε，算法停止并輸出聚類原型P和劃分矩陣U。否則令b=b+1，轉(zhuǎn)向步驟(3)。其中,‖·‖為某種合適的矩陣范數(shù)[11]。

上述算法也可先初始化模糊劃分矩陣，然后再進(jìn)行迭代。由于不能確保FCM收斂于一個最優(yōu)解，因此算法的性能依賴于初始聚類中心[12]。

綜上，F(xiàn)CM算法的聚類結(jié)果受參數(shù)控制，其主要影響參數(shù)分別是：聚類數(shù)目c，加權(quán)系數(shù)m和迭代閾值ε。經(jīng)過FCM算法輸出兩個量：一個c*n的模糊劃分矩陣，以及c個聚類中心點(diǎn)向量。前者表示各樣本點(diǎn)屬于每個類的隸屬度，可根據(jù)模糊集合中的最大隸屬原則確定各樣本點(diǎn)的歸類。聚類中心則是這個類的代表點(diǎn)，表示每個類的平均特征[13]。

1.3 去云算法實(shí)現(xiàn)

鑒于模糊C均值算法的聚類有效性，本文將其與去云技術(shù)有機(jī)地聯(lián)系起來對遙感數(shù)據(jù)進(jìn)行聚類處理，并在此基礎(chǔ)上進(jìn)行類屬提取，將云與植被、土壤、水域等進(jìn)行分離，生成云區(qū)和非云區(qū)數(shù)據(jù)，再運(yùn)用中值濾波和云區(qū)替換的方法實(shí)現(xiàn)去云，如圖1所示。

圖1 算法整體實(shí)現(xiàn)框圖

1.3.1 聚類初始化

讀取HDF格式的NOAA-17遙感數(shù)據(jù)，顯示通道1、4、5。據(jù)相關(guān)資料分析[14-15]，可找出NOAA衛(wèi)星在可見光通道綠色植被的最小反射率和云區(qū)最高反射率之間的相對關(guān)系和特征，以此來確定云區(qū)反射率閾值，如表2所示。

表2 NOAA衛(wèi)星在可見光波段的反射率判云閾值表(%)

待聚類的有云數(shù)據(jù)是灰度數(shù)據(jù)，將控制聚類結(jié)果的加權(quán)系數(shù)m轉(zhuǎn)化為灰度值，并利用云、植被、水域、土壤等在三通道的灰度值的固有特性，將其聚成不同的類別，如表3所示。

表3 衛(wèi)星數(shù)據(jù)反照率值

先將初始聚類中心灰度值設(shè)定在反射率灰度閾值和亮溫閾值范圍內(nèi)，此時要針對不同季節(jié)和不同地區(qū)的地形特征進(jìn)行具體分析，這樣能使聚類效果更加準(zhǔn)確，聚類速度更快。并將設(shè)定的初始聚類中心灰度值作為初始聚類的起始位置，計算聚類迭代次數(shù)，生成N個類別。

讀入數(shù)據(jù)流程如圖2所示，單通道歸一化主要對圖像像素做處理，以保證選擇的待聚類數(shù)據(jù)有云且云區(qū)較清楚。

圖2 讀數(shù)據(jù)流程圖

1.3.2 云區(qū)聚類

先讀取待聚類數(shù)據(jù)，設(shè)置參數(shù)進(jìn)行云區(qū)聚類。并分析聚類后數(shù)據(jù)，若不佳則改變參數(shù)，重新迭代聚類。然后提取類屬，為去云做準(zhǔn)備(圖3、圖4)。

圖3 FCM算法初始化

圖4 聚類提取云區(qū)

1.3.3 云區(qū)處理

鑒于云區(qū)的覆蓋面積小，且它和椒鹽噪聲較類似，故選取對椒鹽噪聲有很好消除作用的中值濾波法處理云區(qū)。中值濾波的優(yōu)點(diǎn)在于既能過濾噪聲，又能保證細(xì)節(jié)不模糊、很好的保護(hù)邊緣輪廓信息。它消除噪聲的效果與模板尺寸和參與運(yùn)算的像素數(shù)有關(guān)。其實(shí)現(xiàn)方法為[16]：首先用一個窗口在圖像上掃描，并將窗口中心與待處理像素位置相重合，讀取窗口內(nèi)包含的圖像像素的灰度值，將其按升序排列起來，取灰度值居中的像素灰度為窗口中心像素的灰度。

待中值濾波處理云區(qū)后，將有云區(qū)域進(jìn)行替換，得到去云數(shù)據(jù)。需注意的是，在中值濾波處理云區(qū)時，使用的模板越大，去云會呈現(xiàn)一定的模糊。流程如圖5所示。

圖5 云區(qū)處理流程圖

2 結(jié)果與分析

2.1 讀取數(shù)據(jù)

通過ENVI軟件讀取HDF格式的NOAA-17遙感數(shù)據(jù)，選取7月25日在四川盆地的部分區(qū)域數(shù)據(jù)，并合成1、4、5通道數(shù)據(jù)，如圖6所示。

圖6 讀取NOAA-17數(shù)據(jù)

分析圖6，四幅圖中的亮色區(qū)域都較小，故認(rèn)為此亮色區(qū)域全為云區(qū)，忽略雪山等非云因素的影響。其中，(a)、(b)、(c)單通道數(shù)據(jù)很好地體現(xiàn)了云的特征：在可見光通道1中，白色區(qū)域?yàn)樵茀^(qū)；在紅外通道4、5中，亮色的小突起部分為云區(qū)；在合成通道數(shù)據(jù)中，亮色區(qū)域?yàn)樵茀^(qū)。將上述四組數(shù)據(jù)作為待聚類數(shù)據(jù)。

2.2 數(shù)據(jù)聚類

采用FCM算法對圖6的數(shù)據(jù)進(jìn)行聚類。初始化聚類中心時，依據(jù)云在可見光和紅外通道表現(xiàn)出的反射率差異大小，即數(shù)據(jù)中的灰度值差異，對數(shù)據(jù)進(jìn)行聚類，效果如圖7。

圖7 FCM算法進(jìn)行數(shù)據(jù)聚類

對比圖7分析發(fā)現(xiàn)：對單通道的聚類結(jié)果(a)、(b)、(c)能真實(shí)反映聚類信息，而對合成通道的聚類結(jié)果(d)則在某些位置會產(chǎn)生信息重疊，造成信息量的部分丟失。

2.3 云區(qū)檢測與提取

對聚類后數(shù)據(jù)進(jìn)行云區(qū)檢測和提取，云區(qū)檢測的目的是保證云區(qū)信息量不丟失，保證其完整性。如圖8所示。

圖8 云區(qū)檢測與提取

圖8對比分析，可得三點(diǎn)結(jié)論：

(1)在單通道云區(qū)數(shù)據(jù)(a)、(b)、(c)中，(b)中通道4在圖像左下方出現(xiàn)與(a)、(c)不相似的情況，這時參考云在三個通道的綜合特征，忽略通道4左下部分的異常情況。

(2)合成通道云區(qū)數(shù)據(jù)圖8(d)是從圖7(d)中直接提取的云區(qū)，對比圖8(e)發(fā)現(xiàn)，它丟失了少量云區(qū)信息。

(3)三通道云區(qū)交集數(shù)據(jù)(e)是三通道各自云區(qū)信息取其交集而得，真實(shí)全面地反映了云在合成通道的特征。

2.4 云區(qū)去除

將通道1作為輔助數(shù)據(jù)，通道4、5作為基準(zhǔn)數(shù)據(jù)。然后采用中值濾波和替換融合處理通道4、5的云區(qū)數(shù)據(jù)，如圖9所示。

圖9 去云結(jié)果

如圖9所示，(a)、(b)是對通道4、5的云區(qū)數(shù)據(jù)進(jìn)行中值濾波后的結(jié)果，濾波模板為3×3大小。對比圖6(b)、(c)分析，云區(qū)得到一定程度上的去除，但呈現(xiàn)一定模糊，丟失了細(xì)節(jié)信息；(c)、(d)是先對三通道進(jìn)行中值濾波、再進(jìn)行通道合成后的數(shù)據(jù)，分別采用3×3大小和7×7大小的濾波模板?？煽闯觯跒V波過程中，采用的模板越大，圖像越模糊；(e)、(f)是先對三通道進(jìn)行中值濾波、然后參考三通道交集云區(qū)數(shù)據(jù)圖8(e)，將通道1作為輔助數(shù)據(jù)、采用替換融合后的去云結(jié)果，分別對應(yīng)3×3大小和7×7大小的濾波模板?？梢钥吹?，通過此方法，有效地對原數(shù)據(jù)進(jìn)行了去云處理。

3 結(jié)論

提出了一種簡單有效的多時相去云方法，通過模糊C均值聚類算法對云區(qū)數(shù)據(jù)聚類，再采用中值濾波、云區(qū)替換達(dá)到去云效果。(1)在對數(shù)據(jù)進(jìn)行聚類時，對單通道的聚類結(jié)果能真實(shí)反映實(shí)際信息，而對合成通道的聚類結(jié)果則在某些位置會產(chǎn)生信息重疊，造成信息量的部分丟失。(2)進(jìn)行云區(qū)檢測和提取時，單通道云區(qū)信息可能出現(xiàn)不相似情況，這時參考云在三通道的綜合特征，忽略異常情況。三通道云區(qū)交集數(shù)據(jù)較合成通道數(shù)據(jù)更能真實(shí)全面反映云在合成通道的特征。(3)對去云后的單通道數(shù)據(jù)直接采用中值濾波，云區(qū)雖得到較好去除，但呈現(xiàn)一定模糊，丟失了細(xì)節(jié)信息；在濾波過程中，采用的模板越大，圖像越模糊；將中值濾波后的結(jié)果進(jìn)行云區(qū)替換，使細(xì)節(jié)更加清晰。此外，本文的算法只針對薄云區(qū)域。