一種改進的ICM遙感影像分割算法

楊 軍， 裴劍杰

(1.蘭州交通大學電子與信息工程學院,蘭州 730070； 2.蘭州交通大學測繪與地理信息學院/甘肅省地理國情監(jiān)測工程實驗室,蘭州 730070)

0 引言

隨著遙感技術的發(fā)展，不同傳感器獲取的遙感影像開始廣泛應用于農業(yè)、林業(yè)、環(huán)境監(jiān)測、城鄉(xiāng)規(guī)劃以及軍事偵察等領域[1-4]。為了充分利用遙感影像數(shù)據(jù)中包含的有用信息，研究人員進行了大量影像分析處理研究，其中影像分割成為該領域的研究熱點之一。目前，已有許多理論和技術被用于影像分割，其中馬爾科夫隨機場(Markov random field, MRF)模型方法融合了局部先驗模型和像元的鄰域結構信息，在影像分割中得到廣泛應用[5-7]。由貝葉斯理論可知，影像分割問題可以把計算已知觀測影像標記場的最大后驗概率問題轉換為等價的吉布斯隨機場(Gibbs random field, GRF)能量函數(shù)最小化問題[8]?，F(xiàn)有的能量函數(shù)最小化算法包括全局優(yōu)化和局部優(yōu)化2大類，全局優(yōu)化算法[9]有Metropolis算法、模擬退火算法(simulated annealing, SA)、圖割算法(graph cuts, GC)和遺傳算法(genetic algorithms, GA)等； 局部優(yōu)化算法[10]有迭代條件模式(iterated conditional model, ICM)算法、馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo，MCMC) 算法和神經網絡算法等。在能量函數(shù)最小化優(yōu)化的眾多算法中，ICM算法是一個經典的局部優(yōu)化算法，通過在求解后驗概率最大化的迭代過程中采取“貪婪”的策略，每次對標記的賦值采取條件概率最大的方式，因而能非常快速地達到能量函數(shù)在鄰域系統(tǒng)內的極值。與全局優(yōu)化算法(如SA算法和GA算法)相比，ICM算法所需計算時間較少，收斂速度較快； 但由于后驗能量函數(shù)具有非凸性，ICM算法很容易收斂到一個局部極值上，無法達到全局最優(yōu)。此外，傳統(tǒng)的ICM算法的優(yōu)化效果受初始標記點的影響較大，一方面，該算法使用K-means聚類法獲取初始標記時只考慮了像元的光譜信息，而忽略了像元鄰域的空間信息，當影像存在較大的噪聲時，分割后的影像往往會出現(xiàn)許多離散斑塊和孤立點； 另一方面，如果沒有一個合適的初始標記，優(yōu)化過程往往會陷于局部最優(yōu)，而對于維數(shù)很高的遙感影像數(shù)據(jù)來說，全局性的更新優(yōu)化很難實現(xiàn)。對于這個問題，一種解決方法是先對原始影像進行低通濾波，再用ICM算法將能量函數(shù)最小化優(yōu)化，但這種方式會導致分割影像邊緣模糊和紋理細節(jié)丟失。

近年來，國內外研究者針對影像分割噪聲抑制和分割結果準確性之間的平衡進行了廣泛研究，對傳統(tǒng)的ICM算法提出了許多改進方法。侯一民等[11]在傳統(tǒng)ICM算法的鄰域基團勢函數(shù)基礎上，引入了鄰域中像元的強度差值以及像元之間的距離因子用于分割。楊紅磊等[12]提出了用MRF描述像元的空間相關性，引入鄰域分類信息作為像元分類的空間約束。這2種方法雖然都很好地抑制了影像噪聲的影響，但是每一次迭代都要重新計算鄰域信息，算法復雜度高。謝昭等[13]利用K-近鄰思想構建圖像鄰域系統(tǒng)，通過設置合適的初始類別數(shù)，用ICM算法進行聚類，提高了聚類準確性； 該方法雖然弱化了初始標記的影響，但是用于分割時結果仍然具有不確定性。Baumgartner 等[14]提出了一種連續(xù)波段融合的影像分割方法，將雙邊濾波器(bilateral filter, BF)結合遙感影像的上下文信息用于分割，類似于運行一個ICM算法； 該算法幾乎可以不用迭代即可達到最終分割狀態(tài)，但分割精度一般。上述方法中沒有一個既能有效抑制分割噪聲，又能使分割結果具有較高精度、適合大多數(shù)遙感影像分割的算法。針對這些問題，本文將保邊去噪效果良好的BF引入到傳統(tǒng)的ICM算法中，用于保留遙感影像的地物邊緣并去除噪聲； 然后采用多閾值最大類間方差法(由日本學者大津(Nobuyuki Otsu)提出的一種自適應閾值確定方法，亦稱大津法(Otsu))遍歷影像直方圖主峰的方式獲取分割閾值，進而確定初始分割類別數(shù)，獲取初始標記； 引入MRF度量空間相關性，并采用航攝影像和多光譜影像數(shù)據(jù)分別進行實驗，以驗證本文算法的有效性。

1 經典ICM算法

基于MRF的影像分割中，分割影像常用2個隨機場來描述，一個是特征場(或稱觀測場)，另一個是標記場。特征場以標記場為條件，用標記類別的似然分布函數(shù)描述特征數(shù)據(jù)或特征向量的分布； 標記場則用先驗知識描述標記場的局部相關性。影像分割需要將特征場中的數(shù)據(jù)與標記場中的先驗知識概率建立有機的聯(lián)系，這就需要將先驗信息所期望的分割結果與真實觀測數(shù)據(jù)的細節(jié)信息進行平衡，而貝葉斯決策理論能夠很好地做到這一點[15]?；贛RF的影像分割通常用貝葉斯決策理論中的最大后驗概率(maximum a posterior, MAP)準則作為MRF建模的優(yōu)化準則，這就是經典的MAP-MRF框架[16]，該框架下的影像分割問題的實質即計算影像的最大后驗概率。由Markov和Gibbs分布的等價性可知，最大后驗概率的求解又可以轉化為對最小能量函數(shù)的計算，而ICM算法是計算最小能量函數(shù)的經典算法。

1.1 MRF影像分割問題

設Y為影像數(shù)據(jù)所表示的特征場或觀測場，y是Y的一個樣本； 標記場X為影像分割結果所描述的隨機場，x是X的一個樣本(即一種分割標記結果)。根據(jù)貝葉斯決策理論，可將影像分割問題表示為

(1)

式中：P(X=x|Y=y)為在影像數(shù)據(jù)Y=y給定的條件下，標記場X=x的后驗概率；P(Y=y|X=x)為在給定標記場X=x的條件下，觀測場Y=y的聯(lián)合分布或者似然函數(shù)；P(Y=y)為觀測場Y=y的聯(lián)合分布。因為觀測場Y=y為影像數(shù)據(jù)，是一個已知量，所以P(Y=y)不隨標記場中任意一個X=x的變化而變化。

如果直接根據(jù)式(1)進行影像分割，幾乎得不到正確的結果，并且計算過程非常復雜。因此，為了進行貝葉斯框架下的影像分割，還需要給出2個基本假設：

1)假設一(條件獨立性)[17]： 假設在給定標記場的一個樣本的條件下觀測場的每一成分相互獨立。

如果要從遙感影像中分割出K種地物，那么標記場中會有K個標記，即標記場中的每一位置可能有K種取值。換言之，觀測場的所有觀測值將被分為K種成分，于是式(1)可表示為

(2)

式中：P(yk|X=x)為在給定標記場X=x的條件下，觀測場成分yk的概率分布。

觀測場P(Y=y)是已知的影像數(shù)據(jù)，而標記場P(X=x)可用MRF模型[18]表示。要實現(xiàn)影像分割，只有P(yk|X=x)是一個未知量，所以還需要做進一步假設。

2)假設二(同質區(qū)域相同分布)[19]： 假設影像中同一類型的像元服從同一種分布(比如高斯分布、指數(shù)分布等)。

(3)

式中：P(Xi=k|XNi)為標記場的局部先驗概率；Ni為位置i的鄰域。

1.2 MAP計算

(4)

由于Markov-Gibbs的等價性，可以使用Gibbs分布函數(shù)表示MRF鄰域系統(tǒng)中的先驗概率P(xi|xNi)，從而通過優(yōu)化能量函數(shù)的方法計算式(4)的解。

由Gibbs分布表達式可得，標記場先驗概率P(xi|xNi)在Gibbs隨機場中的表達式[21]為

(5)

同樣，后驗概率P(xi|yi)也可以用后驗能量函數(shù)來表示，即

(6)

在式(5)和(6)中：U(xi|xNi)為影像分割標記問題的先驗能量函數(shù)；Z為一個稱為配分函數(shù)的歸一化常量；T為溫度常量，通常設置為1[18]。

式(5)和(6)式可以進一步轉化為

P(xi|xNi)∝e-U(xi|xNi),

(7)

P(xi|yi)∝e-U(xi|yi)。

(8)

將式(7)和(8)代入式(4)，并在等式兩邊同時取自然對數(shù)，便將乘積形式轉化為求和形式，這樣更利于計算最優(yōu)解。轉化過程為

(9)

(10)

于是，式(10)表示最大后驗概率P(xi|yi)等價于最小化后驗能量函數(shù)U(xi|yi)： U(yi|xi)為觀測場中地物類別的第i個像素點的似然函數(shù)能量；U(xi|xNi)為標記場中第i個像素點對應的先驗概率能量。因此，最終的能量關系可簡寫為

U(xFi|yi)=U(yi|xi)+U(xi|xNi)，

(11)

式中：xi為第i個像素點的當前分割標記；xFi為最終的分割標記。

1.3 能量函數(shù)計算

ICM算法執(zhí)行過程如圖1所示。

圖1 ICM算法執(zhí)行過程Fig.1 Execution process of ICM algorithm

2 算法改進

大部分遙感影像由于傳感器的姿態(tài)變化以及周期性偏移等原因導致成像時帶有噪聲。這樣的影像直接用于分割時，容易產生離散斑塊和孤立點現(xiàn)象，有時甚至會出現(xiàn)錯誤分割的情況。針對這種情況，本文提出了將改進的ICM算法用于遙感影像分割的方法： 首先利用BF對遙感影像進行預處理，以提高影像的質量； 然后用多閾值最大類間方差法(Otsu)代替K-means作為獲取初始標記的算法，用以克服K-means聚類算法類別數(shù)不易確定和算法復雜度不易控制等問題。

2.1 影像預處理

BF是一種保邊去噪效果良好的濾波器。首先對遙感影像進行雙邊濾波，這樣不但可以保留影像地物邊緣，還可以去除影像噪聲，增加影像對比度。BF由2個系數(shù)構成： 一個是由幾何空間距離決定的濾波器系數(shù)，另一個是由像元差值決定的濾波器系數(shù)。

BF中輸出像元的值，即對影像進行濾波操作后輸出的像元值，依賴于對鄰域像元值的加權，即

(12)

式中： (i,j)為影像中的任一像元； (i′,j′)為像元(i,j)鄰域系統(tǒng)內的像元；f(i,j)為輸入像元的值；g(i,j)為輸出像元的值；w(i,j,i′,j′)為高斯核函數(shù)，由鄰域系統(tǒng)內像元的核系數(shù)乘積構成，即

(13)

式中：hx和hy分別為由像元幾何空間距離決定的定義域核和由像元差值決定的值域核。這2個核均為未知參數(shù)，但因篇幅所限，本文不再贅述參數(shù)的估計方法，為了計算方便，直接使用文獻[14]中定義的參數(shù)。

2.2 初始標記獲取

ICM算法是計算最小后驗能量的常用方法。傳統(tǒng)的ICM算法使用K-means聚類算法獲取初始標記，該算法首先需要根據(jù)初始聚類中心來確定一個初始劃分，然后對初始劃分進行優(yōu)化。此外，K-means算法需要不斷地調整樣本分類和計算調整后的新的聚類中心，因此，當數(shù)據(jù)量較大時，該算法耗用的時間和復雜度是非常大的。

本文采用Otsu法作為獲取初始標記的算法，在劃分L個類別C1,C2, ...,CL的情況下，多閾值類間方差可表示為

(14)

此外，為了保證分割的準確性，本文在類間方差法的基礎上，引入類內方差，即

(15)

正確的分割分類結果應使類間樣本離散度大，同時保證類內樣本的聚集性好，這就要求在類間方差越大的基礎上保證類內方差最小。綜合考慮兩者的關聯(lián)度，本文采用可分性度量準則λ作為分割閾值。當λ達到最大值時，表示類間分離性和類內聚集性達到平衡，此時的閾值是最適合進行分割的最佳閾值。通過遍歷影像直方圖主峰的方式獲得多個分割閾值λ，并確定初始分割類別數(shù)為λ+1?？蓪⒍攘繙蕜tλ表示為

(16)

3 實驗結果與分析

為了驗證改進算法相對于傳統(tǒng)ICM算法的優(yōu)越性和通用性，本文分別采用航攝影像和多光譜衛(wèi)星遙感影像進行了實驗，并對分割效果進行了評價。分割效果評價首先從人眼視覺主觀感受方面分析改進ICM算法用于遙感影像分割時的表現(xiàn)，其次引入分割準確率和Kappa系數(shù)對傳統(tǒng)ICM算法、無濾波ICM和本文算法進行定量比較，最后從分割類數(shù)與Kappa系數(shù)的關系以及分割準確率與分割類數(shù)的關系對算法的魯棒性進行分析。之所以選擇分割準確率和Kappa系數(shù)作為定量評價指標，是因為研究過程中將分割看成是一種特殊的分類過程，因而可以用混淆矩陣(confusion matrix)對分割結果構建重疊矩陣(overlap matrix)進行評價。假設將影像分割為N類，將混淆矩陣的每一項Cij定義為同時屬于分割算法標記的第i類和參考分割結果標記的第j類的像元個數(shù)，則分割準確率Accuracy定義為

(17)

3.1 實驗一

選擇科羅拉多大峽谷的航攝影像進行分割實驗?？屏_拉多大峽谷航攝影像(圖2(a))由國際空間站2014年4月14日拍攝，該地區(qū)位于美國亞利桑那州西北部的凱巴布高原，總面積接近3 000 km2，影像大小為960像元×960像元。通過影像預處理，執(zhí)行程序遍歷影像直方圖主峰的方式獲得分割閾值個數(shù)為λ=7，因而確定初始聚類為8； 而目視判別影像區(qū)域可分為3個部分： 深色植被、淺色植被和其他地物，故確定分割類數(shù)為3。分割結果和分割精度分別如圖2和表1所示。

(a) 原始圖像(b) 傳統(tǒng)ICM (c) 無濾波ICM (d) 本文算法

圖2 科羅拉多大峽谷航攝影像不同算法分割結果對比Fig.2 Comparison of segmentation results of Grand Canyon aerial image by using different algorithms表1 不同算法的分割精度Tab.1 Segmentation accuracy of different algorithms

仔細對比圖2中的分割影像可以看出，傳統(tǒng)ICM算法的分割效果一般，并且出現(xiàn)了一些離散斑塊和孤立像元點； 無濾波ICM算法由于改進了初始標記的獲取方式，分割效果與傳統(tǒng)ICM算法相比有了一定的改善，但因原始圖像沒有經過濾波處理，分割效果與本文算法相比在邊緣和細節(jié)方面的表現(xiàn)差了很多。本文算法通過改進初始標記獲取方法和加入BF，使得分割效果明顯優(yōu)于前2種算法，并且在分割后的影像中邊緣更加平滑，細節(jié)信息更加突出。此外，從表1中的數(shù)據(jù)可以看出，本文算法相比于傳統(tǒng)ICM算法，在分割準確率上提高了18%，Kappa系數(shù)提高了14%，表明本文算法用于遙感影像分割時的結果更加準確。

3.2 實驗二

選擇印尼科摩多國家公園真彩色影像和阿拉胡埃拉湖多光譜影像進行分割實驗。印尼科摩多國家公園真彩色影像(圖3(a))是由Terra衛(wèi)星搭載的高級星載熱輻射和反射輻射計(advanced spaceborne thermal emission and reflection radiometer，ASTER)傳感器于2000年7月20日獲取的，影像大小為550像元×550像元。圖3(a)中藍色表示水，灰色表明裸露地面，綠色表示植被，白色表示云。初始聚類數(shù)為10，分割類數(shù)為4，實驗結果如圖3所示。

(a) 原始圖像 (b) 傳統(tǒng)ICM (c) 無濾波ICM(d) 本文算法

圖3科摩多國家公園影像不同算法分割結果對比

Fig.3ComparisonofsegmentationresultsofKomodoNationalParkimagebyusingdifferentalgorithms

阿拉胡埃拉湖多光譜影像(圖4(a))是由美國航天局Earth Observing-1(EO-1)衛(wèi)星搭載的高級陸地成像儀(advanced land imager，ALI)于2010年12月17日獲取的，影像大小為600像元×600像元。由于暴雨的影響，大量泥沙沖入該地區(qū)湖泊、河流中，導致水域顏色為黃色，圖4(a)中白色表示云，其余地物可看成是植被。初始聚類數(shù)為12，分割類數(shù)為3，實驗結果如圖4所示。

(a) 原始圖像(b) 傳統(tǒng)ICM (c) 無濾波ICM (d) 本文算法

圖4阿拉胡埃拉湖多光譜影像不同算法分割結果對比

Fig.4ComparisonofsegmentationresultsofAlajuelaLakemultispectralimagebyusingdifferentalgorithms

對比圖3和圖4的分割結果可以看出，用傳統(tǒng)的ICM算法進行多光譜影像分割時，出現(xiàn)的錯分現(xiàn)象比較嚴重，并且分割區(qū)域邊緣不平滑。反之，本文算法無論在分割效果方面的表現(xiàn)，還是在分割區(qū)域邊緣的表現(xiàn)，均優(yōu)于傳統(tǒng)ICM算法和無濾波ICM算法。

3.3 實驗三

選擇俄羅斯西伯利亞地區(qū)的Chara沙地Landsat 8多光譜衛(wèi)星遙感影像和新西蘭拉凱阿河航攝影像進行分割實驗。

俄羅斯西伯利亞地區(qū)Chara沙地的Landsat 8多光譜影像(圖5(a))于2013年9月獲取，由11個波段組成，本次研究使用Band3(R)、Band2(G)、Band1(B)波段組合合成模擬真彩色影像，大小為1 712像元×1 712像元。影像中可分辨出的地物主要包括沙地、植被、云、湖泊、河流及其他等6類，初始聚類為10，分割類數(shù)為6，實驗結果如圖5所示。

(a) 原始圖像(b) 傳統(tǒng)ICM (c) 無濾波ICM (d) 本文算法

圖5Chara沙地多光譜遙感影像不同算法分割結果對比

Fig.5ComparisonofsegmentationresultsofCharasandmultispectralremotesensingimagebyusingdifferentalgorithms

從圖5可以看出，本文算法的分割結果均優(yōu)于傳統(tǒng)ICM算法和無濾波ICM算法。

圖6示出實驗三影像通過3種分割算法獲得的Kappa系數(shù)與影像分割類別數(shù)目的關系。

圖6 Kappa系數(shù)與分割類別數(shù)的關系Fig.6 Relationship between Kappa coefficients and number of clusters

從圖6可以看出，當分割類數(shù)少于4時，3種分割算法均表現(xiàn)為低度一致性，說明3種算法對于簡單地物類別的遙感影像分割效果均一般； 當分割類數(shù)介于5到8之間時，本文算法相比于傳統(tǒng)ICM算法影像分割結果體現(xiàn)出了高度一致性； 而分割類數(shù)多于9時，由于地物類別復雜以及光譜特征相近等原因產生了一些誤分割的現(xiàn)象，導致Kappa系數(shù)降低，但本文算法仍然具有中度一致性的分割結果。

圖7(a)是選擇的新西蘭拉凱阿河航攝影像，大小為800像元×800像元。人眼從影像中可分辨出的地物包括耕地、人工林地、灌木林、草地、河邊灘涂、水域和未利用土地等類型。通過遍歷影像直方圖主峰，確定初始聚類數(shù)為16，分割類數(shù)為8。實驗結果如圖7所示。圖8示出圖7中拉凱阿河航攝影像取不同分割類數(shù)時獲得的分割準確率。

(a) 原始圖像(b) 傳統(tǒng)ICM (c) 無濾波ICM (d) 本文算法

圖7拉凱阿河航攝影像不同算法分割結果對比

Fig.7ComparisonofsegmentationresultsofRakaiaRiveraerialimagesbyusingdifferentalgorithms

圖8 拉凱阿河航攝影像分割準確率與分割類別數(shù)的關系Fig.8 Relationship between segmentation accuracy and number of clusters from aerial images of Rakaia River

從圖8可以看出，本文算法的分割類數(shù)與實際類數(shù)一致，并在分割類數(shù)與實際類數(shù)有較小偏差時獲得的分割準確率均高于傳統(tǒng)ICM算法和無濾波ICM算法。

4 結論

針對傳統(tǒng)的迭代條件模式(ICM)算法用于遙感影像分割時容易出現(xiàn)離散斑塊和孤立點的現(xiàn)象，本文改進了傳統(tǒng)的ICM算法。本文算法首先加入雙邊濾波器(BF)進行遙感影像預處理，在降噪的同時較好地保留了邊緣信息； 然后用多閾值最大類間方差法(Otsu)代替原有的K-means算法，通過遍歷影像直方圖主峰的方式獲取多個分割閾值，進而確定初始聚類數(shù)目，獲取初始標記。本文算法克服了傳統(tǒng)ICM算法中用K-means獲取初始標記時類別數(shù)不確定和算法復雜度不易控制的問題。實驗結果表明，改進的ICM算法用于遙感影像分割時分割效果優(yōu)于傳統(tǒng)ICM算法和無濾波ICM算法。

然而，對于地物類別復雜的遙感影像，利用本文算法進行分割時仍可能產生誤分割和未分割現(xiàn)象，這將是今后的研究重點。