關崢嶸 譚毅華 田金文
(華中科技大學自動化學院 武漢 430074)
隨著遙感技術的不斷發展,我們能獲得的遙感衛星影像數據量也越來越多,遙感衛星影像的分辨率也越來越高,這些衛星影像攜帶的地面目標信息也越來越多。但同時,地球表面是包裹著大氣層,遙感傳感器常常受到大氣或者云層的干擾[1],使得我們對地球表面的觀測受到遮擋,因此針對遙感衛星影像的云檢測也就成為了遙感衛星圖像處理的關鍵性技術之一,云檢測也被認為是遙感影像后續處理(識別和分類等)的預處理技術。此外在氣象學方面,針對云的檢測也被認為是進行氣象預測的重要手段之一[2]。
目前來說,針對云檢測是方法大致分為三種:閾值法,紋理分析法和基于機器學習的方法?;陂撝档姆椒ɡ昧瞬煌镔|反應在圖像上不同的明暗程度,又可以分為簡單的單閾值[3~4],動態閾值法[5];紋理分析法從本質上來說就是利用同類物質之間的相似性和不同物質之間的不連續性,往往需要結合其他信息一起進行判別,比如多時相[6~7]和物體反射率[8],這些方法雖然都取得了一定的進步,但是大都適應性不強,利用了一些先驗知識,導致了在場景復雜的情況下的檢測效果較差。而基于機器學習的方法往往能獲得比較好的效果,利用訓練樣本來進行特征訓練,使得算法的適應性大大加強,支持向量機(Support Vector Machine,SVM)是最常用的分類器[9~12],還有利用神經網絡模型來進行訓練[13]。此外一些比較新穎的圖像處理技術也被應用到云檢測當中,比如超像素分割[14~15],模糊聚類[16~17]。利用機器學習來進行云檢測時,檢測效果往往取決于選取的特征描述能力的強弱,因此如何尋找出一種好的描述特征,是我們目前面臨的挑戰之一。論文針對該問題,提出了一種基于改進MR 的描述特征來進行云檢測,并在高分一號衛星圖像數據集上獲得了令人滿意的效果。
最大響應濾波器(Maximum Response Filter)[18]采用了38個濾波器,但最終只選取8組最大的響應值作為輸出,因此又被稱為MR8。MR8 的38 個濾波器包含了18 個高斯一階導數濾波器(三個尺度六個方向),18 個高斯二階倒數濾波器(三個尺度六個方向),尺度為(σx,σy)={(1 ,3),( 2,6),(4,12)},方向角取值為θ={0°,30°,60°,90°,120°,150°};1個高斯濾波器和1個高斯型的拉普拉斯濾波器(Laplacian of Gaussian),標準差σ 取10。38 個濾波器的時域圖如圖1所示。
由于MR8 集合了高斯一階和二階導數濾波器,因此能有效地提取圖像的邊緣信息,此外由于采用了不同尺度的濾波器使得提取的特征能適應不同的圖像分辨率,不同的方向濾波器使得其在非均質的紋理下獲得更強的分類能力,而最終只選取最大的8個響應值,使得MR8特征具有旋轉不變的性質。綜上所述,MR8特征具有較強的特征描述能力。
在計算MR8 特征的過程中,由于需要使用高斯核與圖像進行大量的卷積操作,而通常遙感衛星圖像分辨率較高,幅面較大,導致整個提取特征的過程中非常耗時,因此如何快速地實現MR8 是急需解決的問題。受之前的工作啟發[19~20],本論文使用遞歸(Infinite Impulse Response,IIR)數字濾波器來快速近似,以設計一種快速實現MR8的方法。
一維信號的高斯核卷積過程可用兩次遞歸的形式來實現,第一個是前向遞歸(Forward,n遞增):

第二個是后向遞歸(Backward,n遞減):

可以表示為如圖2所示。

圖2 一維信號的卷積通過遞歸來實現
式(1)和式(2)中的參數B,b0,b1,b2,b3,由下面的式子決定:

其中m0,m1,m2為常數,取值分別為1.16680,1.10783,1.40586。而q 是與標準差σ 有關的數值,一般按照如下方法取值:

根據上面的方法把一維信號的高斯卷積操作變成了兩次遞歸來實現,這很容易推廣到二維圖像的高斯卷積操作上去,最后只需要做四次遞歸操作即可:從上到下,從左到右,從下到上和從右到左。
為了盡可能好地描述出云的圖像特征,我們這里并沒有直接采用原始的MR8 特征進行操作,而是利用了K-D 樹來構建MR8 直方圖特征用于判別。具體的步驟如下所示。
1)計算訓練數據集上的MR8 特征,并利用K均值聚類算法得到紋理特征的數據字典;
2)利用得到的數據字典構建K-D樹;
3)得到圖像塊的MR8 特征,將得到的紋理特征進行K-D樹查詢;
4)生成查詢結果的歸一化直方圖,并利用該直方圖特征作為云判別的特征。
本文提出的算法步驟可以如圖3所示:

圖3 算法流程圖
在訓練階段,分為兩步:第一步是構建數據字典,首先是提取圖像塊(大小為128*128)的MR8特征,進行K均值聚類以后為了避免特征維數過大導致的維數災難,我們這里使用了主成分分析(Principal Component Analysis,PCA)進行了降維處理,最后得到的這些聚類中心就構成了紋理特征的數據字典,然后利用數據字典構建K-D 樹,如果按照傳統的方法在得到數據字典后需要與聚類中心進行歐式距離遍歷,再去統計得到結果的直方圖特征,這樣的操作不僅耗時而且特征的描述效果也不甚理想;第二步是利用數據字典得到圖像塊的紋理特征,隨后把這些紋理特征利用支持向量機(Support Vector Machine,SVM)來訓練出一個分類器。
在測試階段時與上面類似,只需要得到圖像塊的直方圖特征,然后利用SVM去判別即可。
在實驗過程中,訓練所采用的數據來自與谷歌地圖(Google Map),測試所用的數據來自于高分一號衛星數據(拍攝于云南某地,分辨率為2m)。訓練數據集一共有6000 個樣本,包括2000 個正樣本(云區圖像),4000 個負樣本(包括森林,湖泊,城市,草地等區域圖像)。
為了確定實驗中的數據塊大小、詞典的單詞數量和是否采取PCA等參數,我們做了一系列的實驗來進行對比。
根據實驗結果,最終我們選定使用128*128 的圖像塊大小,數據字典數量為1024,采用51*51 的高斯卷積核,因此我們在圖像上得到的MR8 特征其實是78*78*8 維的特征,再進行K-D 樹查詢后,得到了78*78 維的查詢結果。接下來對其進行量化并統計其直方圖分布,作為該圖像塊的描述特征。從實驗中我們可以看到,進行PCA操作有效地縮短了算法所需的時間,而且去除了冗余數據后也使得檢測精度有小范圍的提升。

圖4 高分一號衛星的測試遙感圖像示例

表1 不同的單詞個數對準確率的影響(圖像塊大小為128*128)

表2 不同的圖像塊對準確率的影響(單詞數量為1024)

表3 PCA對時間和準確率的影響(2048*2048大小圖像)
測試圖像我們選用了35 幅2048*2048 大小的衛星圖像,并在不同的算法上進行了對比,其中我們用F_alarm 表示虛警率,準確率precision 通過TP/(TP+FP)得到,召回率通過TP/(TP+FN)得到,精度accuracy 通過(TP+TN)/(TP+FN+TN+FP)計算得到。

表4 不同方法的云檢測結果對比
從表中可以看到,我們的算法以微小的時間代價在各方面都得到了最好的結果,CI 和DT 的方法雖然耗費時間比較短,但是損失了較大的準確率,因此可以認為本論文提出的算法是優于這些方法的。
本文提出了一種通過快速MR8 和K-D 樹構建直方圖特征來進行云檢測的算法,與一般的檢測方法不同的是,該方法通過四個方向上的遞歸操作來快速實現MR8 特征的計算,并利用數據字典和K-D 樹構建直方圖特征,基于MR8 的直方圖特征能有效地描述云塊區域,具有較強的適應性而且執行速度較快,能滿足不同條件下快速進行遙感衛星圖像的自動云判別的要求。