一種超像素上Parzen窗密度估計的遙感圖像分割方法

張大明， 張學勇， 李璐， 劉華勇

(1.安徽建筑大學數(shù)理學院，合肥 230022; 2.安徽省建筑聲環(huán)境重點實驗室，合肥 230601)

0 引言

高分辨率多光譜影像能提供幾何形狀、光譜和紋理等豐富的地物信息，對地物目標識別精度的提高有很大幫助?，F(xiàn)在，已經(jīng)有許多衛(wèi)星能夠提供米級和亞米級的高分辨率遙感圖像，它們在農業(yè)、采礦、防災減災和環(huán)境調查等各個領域中有廣泛應用[1-3]。圖像分割是指將圖像分割成若干區(qū)域，每個區(qū)域內部像素之間具有相似性，同時不同區(qū)域之間具有較大的互異性。分割質量的好壞對后續(xù)的目標識別和解譯等工作具有重要作用。

一般的圖像分割方法，多是基于像素的，即根據(jù)像素本身的灰度、顏色、紋理等特征進行分割； 根據(jù)具體方法不同，又有基于閾值、基于邊緣、基于聚類，以及基于圖論等許多方法[4-5]。對于遙感圖像，不可避免地容易受到各種噪聲干擾。而基于像素的方法，除了計算量大、實時性差以外，抗噪性也較差。而基于分裂-合并的分割策略則能有效地克服這些問題?；诜至?合并的分割算法是從初始分割出發(fā)，基于區(qū)域特征的相似度，迭代地合并相似的2個相鄰區(qū)域，直到滿足某個代價函數(shù)而停止合并。這里涉及到3個方面： ①初始分割，產(chǎn)生的過分割結果直接影響到最終分割結果的質量； ②區(qū)域特征提取，所提取特征是后續(xù)區(qū)域合并的依據(jù); ③區(qū)域合并，包括區(qū)域的相似度獨立和合并準則的確定。學者們從這3方面對合并-分裂算法進行了廣泛而深入的研究[6-8]。

與一些傳統(tǒng)方法(如分水嶺算法)相比，超像素將圖像中局部特征(如紋理、灰度和顏色等)相似作為一個視覺意義上的整體，替代大量原始像素來表達影像特征，更符合人類對圖像的感知方式，減少了冗余信息，大大降低了計算的復雜度。因此，采用超像素方法進行初始分割成為主流[9-11]。Achanta等[12]提出簡單線性迭代聚類(simple linear iterative clustering, SLIC)算法生成超像素。該算法利用顏色和空間特征，生成形狀規(guī)則、尺寸均勻的超像素塊，其算法簡單，計算速度快，而且超像素塊數(shù)目調節(jié)方便，是超像素分割的經(jīng)典算法。

目前使用最廣泛的遙感圖像分割方法是eCognition軟件中的分形網(wǎng)絡演化算法(fractal net evolution approach, FNEA )。該方法是一種多尺度分割方法，通過計算和比較圖像光譜和形狀的相似性，自下而上的逐層合并對象。但該方法需要手動反復調節(jié)尺度、緊致度和異質度等多個參數(shù)，才可能得到較好的分割效果[13]。概率密度估計是概率統(tǒng)計和模式識別中基本問題之一，主要分為參數(shù)估計和非參數(shù)估計。密度估計在圖像分割的應用也很廣泛，最為著名的是基于非參數(shù)估計的均值漂移算法[14-15]； 此外非參數(shù)估計還應用于基于閾值的圖像分割中[16]； 而基于參數(shù)估計的有限混合模型也經(jīng)常用于圖像分割的建模中[17-18]。

本文基于分裂-合并的分割策略，在對圖像進行超像素初始分割的基礎上，對每個超像素塊進行非參數(shù)密度估計，并基于估計的密度，重新將超像素塊合并聚類，提出一種新的遙感圖像分割算法。

1 基于超像素及其最小生成樹上Parzen窗密度估計的圖像分割算法

本文方法的主要思路是： 先將圖像進行初始超像素分割，把每個超像素塊看作圖上一個頂點，然后在圖的最小生成樹上測量頂點間的距離，采用Parzen窗估計方法，在圖上估計每個頂點的概率密度，進而尋找不同密度函數(shù)的模態(tài)進行聚類。其基本原理示意圖如圖1所示。

(a) 圖(b) 最小生成樹(c) Parzen窗密度估計(d) K-means聚類

以下先概述SLIC算法，然后給出超像素塊的特征提取以及度量方法，再給出基于Parzen窗的密度估計方法，最后詳細討論提出的算法，并給出算法流程圖。

1.1 SLIC算法

SLIC超像素分割算法的步驟為[12]：

1)初始化種子點。根據(jù)輸入圖像的像素點個數(shù)L，均勻地在圖像中設置k個種子點，則超像素大小為L/k。

2)種子點處理。在種子點的3×3鄰域內對所有像素進行搜索，將種子點移動到該鄰域內梯度最小的位置，以避免種子點落在邊緣位置，從而確定更好的聚類效果。

3)距離測量。將圖像轉化到Lab顏色空間，分別計算每一個像素點i與該種子點j的顏色距離和空間距離，公式分別為：

(1)

(2)

(3)

4)迭代聚合與連通性處理。重復上述步驟，直到聚類結果穩(wěn)定，一般迭代十余次即可終止。但可能會存在一些孤立像素或是多連通像素，使用分量連通方法以產(chǎn)生緊致連通的超像素塊。

1.2 區(qū)域特征的提取與度量

與傳統(tǒng)傅里葉方法不同的是，Gabor濾波器可以同時在頻率域和空間域都有較好的分辨力，可以方便地在方向和尺度2個參數(shù)上進行調整，以符合人類的視覺特性。Gabor濾波器的小波特性使其可以詳盡地描述圖像紋理特征[19]。

設m,n分別為尺度參數(shù)和方向參數(shù)的Gabor濾波器組表示為gmn(x,y)， 給定圖像I(x,y)，將該圖像與Gabor濾波器組的卷積得到該圖像的Gabor變換Wmn(x,y)，即

Wmn(x,y)=I(x,y)?gmn(x,y)。

(4)

對于濾波后的圖像Wmn，可得其均值μmn和標準差σmn，則圖像I(x,y)的紋理特征向量為[μ11,σ11,μ12,σ12,…,μmn,σmn]，共計2mn維。實驗中，取3個尺度(m=3)，8個方向(n=8)，因此紋理特征為48維向量。

1.3 Parzen窗密度估計

Parzen窗密度估計(Parzen windows density estimation)又稱核密度窗估計，是一種重要的非參數(shù)密度估計方法，最早由Parzen[20]提出，它不需要向參數(shù)估計方法那樣，假定數(shù)據(jù)分布的概率密度函數(shù)參數(shù)形式已知，它可以對任意分布的數(shù)據(jù)作密度估計?，F(xiàn)有N個d維空間樣本x1,x2,…,xN，Parzen窗密度估計即以這N個已知樣本來對空間中任一點x的概率密度進行估計，即

(5)

式中： 參數(shù)h為核函數(shù)的帶寬(bandwidth)；K(·)為核函數(shù)(kernel function)，一般為對稱和單峰的函數(shù)，且滿足如下條件，即

(6)

1.4 提出的算法

經(jīng)過超像素初始分割后，2個超像素塊直接的紋理相似度定義為：

(7)

根據(jù)超像素塊之間的紋理特征的相似度，每個頂點(超像素塊)vi的概率密度由與它相連的頂點按照Parzen窗密度估計方法可以估計為：

(8)

式中G(·)為Gaussian核函數(shù)。在式(8)中僅僅包含了紋理信息，而沒有考慮空間信息，而實際上鄰接的超像素塊極有可能屬于同類，即具有近似的概率密度。因此，本文提出同時包含紋理信息、光譜信息和空間信息的Parzen窗密度估計方法。

最小生成樹是最小權重生成樹的簡稱，它是無向圖中包含所有頂點且權值和最小的子圖，任意2個頂點之間距離同時蘊含了特征的相似度(權值)和空間的信息[25]。因此，用最小生成樹上距離帶入Parzen窗密度估計方法中描述頂點(超像素塊)的概率密度。

首先，可以用經(jīng)典的算法，如Krusal算法、Prim算法[26]，在以超像素塊為頂點的區(qū)域鄰接圖上建立最小生成樹; 然后在最小生成樹上尋找任意2頂點vi和vj的路徑{vi,vi+1,…,vj}，其中vi和vi+1是路徑中相鄰的頂點，這樣，頂點vi和vj的距離DMST為:

(9)

式中wp,p+1為最小生成樹上連接相鄰頂點vp和vp+1的權值。然后，將式(8)中的wij替換為DMST(i,j)進行Parzen窗估計，即

(10)

可以看到，在Parzen窗密度估計中的距離DMST(i,j)中，既含有SLIC算法的光譜信息，也融入了超像素塊的紋理信息，同時在最小生成樹上度量也自然引入了空間近鄰信息。在Parzen窗密度估計中是以每個頂點為中心的窗函數(shù)(這里是Gaussian核函數(shù))到待估計點的疊加。由Gaussian核函數(shù)的性質，對概率密度有所貢獻的點只能是與待估計點空間距離較近(累計權重較小)的點，以及與待估計點特征(光譜特征、紋理特征)相似(邊上權重較小)的點。所以，本文提出的這種概率密度估計方法是合理的。

對于萬能軋機的重軌生產(chǎn)，軋機對型鋼斷面的軋制要高于普通軋面20 m左右。在離鋼軌尾部10 m的范圍內，會存在一個高于正常軌道0.5 mm以上的“高點”，該“高點”會在鋼軌有0.5 m左右的持續(xù)長度。對于這些影響列車運行的“高點”，傳統(tǒng)“高點”處理方式為器具打磨，但器具打磨會嚴重影響鋼體的質量與壽命。而使用全軋程熱力耦合數(shù)值模擬系統(tǒng)，能夠完成E孔型和UF孔型的連軋工作。其中UF孔型屬于半封閉的軌道孔型，軋件在離開E軋機進入UF軋機的過程中，會由于半封閉孔型而發(fā)生“甩尾”現(xiàn)象，軋件尾部會明顯高于脫離軋機時的軌高。

最后，將具有相似的概率密度的頂點被認為屬于同一模態(tài)(類)，采用K-means算法完成聚類，得到最終分割結果。整個分割算法流程如圖2所示。

圖2 提出的多光譜遙感圖像分割算法流程

2 實驗結果及討論

本文多光譜遙感圖像分割的實驗平臺為2.8 GHz的Pentium4處理器、8 G RAM的計算機，采用軟件為MATLAB R2004a。實驗1是討論參數(shù)的選擇對提出算法分割效果的影響； 實驗2和實驗3是將提出的算法與eCognition軟件中不同參數(shù)FNEA算法作對比，驗證提出算法的有效性。

為了更有效地評價分割結果，從定性和定量2個方面進行評價。定性即為主觀的目視對比判別； 定量則基于手動標注的實際地物邊界采用P-R方法[27]計算，公式分別為：

Precesion=TP/(TP+FP)，

(11)

Recall=TP/(TP+FN)，

(12)

式中：Precesion為準確率，其值越大，表示分割精度越高；Recall為召回率，其值越大，則表示邊界的附著度越好；TP為被正確分割的地物目標像素數(shù)目；FP為被誤分為地物目標的背景像素數(shù)目；FN為被誤分為背景的地物目標像素數(shù)目。

2.1 實驗分析一

該實驗遙感圖像選自NWPU VHR-10數(shù)據(jù)集[28]。圖像采集自Google Earth，包括10個類別目標的高分辨率RGB彩色圖像和紅外圖像。實驗圖像為數(shù)據(jù)集中178號圖像的局部，其中包括船只、水體、植被、建筑物、橋梁等地物目標，圖像大小為266像素×385像素，如圖3所示。該實驗目的是討論提出算法中的參數(shù)設置對分割效果的影響。在提出的算法中主要需要設置3個參數(shù)： ①SLIC算法中的S，用于控制超像素塊大小； ②Parzen窗帶寬h； ③K-means的聚類數(shù)。聚類數(shù)的選擇無疑會對聚類效果有較大影響，許多文獻討論如何自動選擇聚類數(shù)的問題[29-30]，這不是本文主要關注點(實驗中發(fā)現(xiàn)，只要聚類數(shù)在一個大致合理的范圍內，分割結果相對穩(wěn)定，限于篇幅，沒有給出實驗圖像)，實驗中還是依據(jù)經(jīng)驗采用手動給定聚類數(shù)。該實驗中，所有的聚類數(shù)手動選擇為5。

圖3 實驗1的原始圖像

表1中給出了參數(shù)S和h取不同數(shù)值時的分割結果。其中，從左至右，第1—3列依次是S設置為5，20和50時的分割結果； 從上到下，第1行為超像素預分割結果，第2—5行分別是帶寬h取0.5，5，20和50時的分割結果。參數(shù)S越小，預分割的超像素

表1 實驗1的分割結果

塊數(shù)目就越多; 反之，S越大，超像素塊越大，其數(shù)目就越少。因為采用的是分裂-合并的分割策略，若S設置過大，會使得后續(xù)合并的分割結果過于粗糙，容易產(chǎn)生欠分割問題(如表1第3列所示)。由于合并過程是基于Parzen窗密度估計，這也要求樣本點不可以過少。由于樣本點過少，導致密度估計不夠準確，也是導致表1第3列中分割效果較差的原因之一。但另一方面，若S設置太小，超像素數(shù)目過大，會使得計算量(包括超像素分割、紋理特征測量和Parzen窗密度估計)顯著增加。

表2中是參數(shù)S取不同值時的平均計算時間的比較，容易見到S取5時比取20時計算時間增加了近2倍，同時分割效果并不更好，且更容易出現(xiàn)過分割問題(見表1第1列與第2列的比較)。

表2 實驗1中參數(shù)S的選取對計算時間的影響

另一個參數(shù)帶寬h，涉及到密度估計的準確性，自然會影響到最后聚類結果(如表1第5行h=50時，帶寬取得過大，嚴重降低了分割效果)。但同時，實際的分割結果對帶寬的選擇相對而言并不敏感，如第2—4行，帶寬差異并不小，可分割結果卻遠沒有預計的大。實際上，實驗中選擇了多個帶寬，結果表明，只要帶寬不是過大(如第5行)或過小，分割結果差異不會太大。而實驗中采用文獻[24]的方法都能夠選擇到理想的帶寬。所以，后面的實驗都是采用該方法自動選擇帶寬。

從表3的P-R定量評價中也可以看出，這2個參數(shù)的選擇在Precision指標和Recall指標上基本也可以反映出于目視判別解讀相同的結論。

表3 實驗1的定量評價

2.2 實驗分析二

實驗2圖像是2010年4月15日青海玉樹7.1級地震發(fā)生后由QuickBird衛(wèi)星拍攝的多光譜圖像，分辨率為0.6 m[31]。實驗中采用了其局部圖，大小為650像素×838像素，包含植被、建筑物、道路、山體等地物目標，如圖4(a)所示。圖4(b)和(c)分別是FNEA算法尺度參數(shù)分別設置為50和150的分割結果，其中形狀異質度和緊致度參數(shù)分別設置為默認值shape=0.1和compactness=0.5。圖4(d)—(f)都是本文方法的分割結果，其中帶寬h由文獻[24]的算法自動選擇，類別數(shù)設置為6，圖4(d)—(f)參數(shù)S分別設置為5，20和30。

(a) 原始圖像(b) FNEA(50)(c) FNEA(150)

(d) 本文方法(S=5)(e) 本文方法(S=20)(f) 本文方法(S=30)

實驗3圖像是NWPU VHR-10數(shù)據(jù)集[28]中286號圖像的局部圖，圖像大小為751像素×979像素，包括飛機、水體、艦船、植被、機場地物標志等地物目標，如圖5(a)所示。與圖4中實驗2類似，圖5(b)和(c)分別是FNEA算法尺度參數(shù)分別設置為150和200的分割結果，其中形狀異質度和緊致度參數(shù)也都設置為默認值shape=0.1和compactness=0.5。而圖5(d)—(f)也是參數(shù)S分別設置為5，20，30的分割結果，帶寬h由文獻[24]中算法自動選擇，類別數(shù)設置為5。

(a) 原始圖像(b) FNEA(150)(c) FNEA(200)

(d) 本文方法(S=5)(e) 本文方法(S=20)(f) 本文方法(S=30)

從圖4(b)和(c)以及圖5(b)和(c)中可以看到，在FNEA方法中，當尺度參數(shù)變大，則過分割的問題會得到一定程度克服，但很難完全避免。如圖4(b)和(c)中的植被部分，由于植被的光譜異質度較大，產(chǎn)生了明顯的過分割，而圖5(b)中的水體部分，由于光譜相似度較大，也產(chǎn)生了比較明顯的誤分割問題。而尺度參數(shù)大，又很容易產(chǎn)生欠分割，如圖5(c)。同時，圖4中建筑物的光譜信息“同譜異物”的情況比較明顯，大多數(shù)尺度情況下，都容易出現(xiàn)過分割情況。所以，很難找到合理的尺度參數(shù)，同時避免以上各種問題的出現(xiàn)。從表4和表5中的定量評價來看，實驗2和實驗3中分別取比較典型的大小2個不同尺度參數(shù)FNEA方法分割結果，其大尺度下的Precision和Recall指標均較小尺度下的有所提高，但非常有限。這也從另一個方面印證了上面的分析。

表4 實驗2定量評價

表5 實驗3定量評價

在本文方法中，初始分割階段，SLIC算法考慮了像素的光譜特征； 然后又融入了超像素塊的紋理信息； 距離測量時，利用最小生成樹結構，整合了超像素塊的空間近鄰信息。在此融合了各種信息的距離測度基礎上，用Parzen窗估計頂點密度，距離較近的頂點(超像素塊)，其密度應該更接近，后續(xù)的K-means算法更傾向于將其分為一類。因而，從圖4和圖5的(d)—(f)中也可以看到，本文方法能夠在一定程度上處理好過分割、欠分割和誤分割等問題，即使地物目標復雜的遙感影像，也可以產(chǎn)生較為合理的分割結果。

但是，正如實驗1中所指出的，提出的算法中超像素預分割中參數(shù)S還是對分割結果有一定影響，如圖4和圖5中(d)比(e)和(f)存在更明顯的過分割問題。所以參數(shù)S取得過小，超像素塊數(shù)目過多，合并的效果并不好； 當然，參數(shù)S過大也容易出現(xiàn)分割效果粗糙(實驗2和3中沒有再給出分割效果圖)。好在分割結果對參數(shù)S并不十分敏感，只要位于合理范圍(實驗中一般在10～30之間)，其分割結果相對比較接近(如圖4和圖5(e)和(f))。在表4和表5的定量分析中也可以看到，提出的算法中參數(shù)S取5，20和30時，Precision和Recall指標均比FNEA方法好，而參數(shù)S取20和30的結果則提高更明顯，且比較接近。

3 結論

針對傳統(tǒng)基于像素的遙感圖像分割算法易受噪聲干擾、計算效率低等缺點，在分裂-合并的框架下，本文提出了一種將超像素和Parzen窗密度估計相結合的遙感圖像分割方法，該方法預分割階段是采用SLIC算法將遙感圖像粗分割為超像素，每個超像素塊視作圖結構中的一個頂點，在圖的最小生成樹上測量頂點間距離，并采用Parzen窗估計每個頂點密度，然后以這些概率密度為指標對超像素塊進行K-means聚類，以達到“融合”目的，得到最終分割結果。提出的方法中有2個方面的主要工作： 一是構建了基于最小生成樹的整合了光譜信息、紋理特征和空間信息的距離測量方式； 二是提出了基于圖上頂點的Parzen窗密度估計的聚類方法。

為驗證提出算法的有效性，在多種多光譜遙感圖像上進行分割實驗。通過算法中不同參數(shù)下的分割結果，以及與不同參數(shù)的FNEA算法的分割結果進行比較研究，可以發(fā)現(xiàn)提出的算法結果穩(wěn)定，魯棒性較好，同時能夠克服過分割、誤分割和欠分割等問題，得到較好的分割效果。