自動確定類別數的RJMCMC+SA圖像分割算法研究

張文坤，汪西原,2，韓佳雪

(1.寧夏大學物理與電子電氣工程學院，寧夏 銀川 750021；2.寧夏沙漠信息智能感知自治區重點實驗室，寧夏 銀川 750021)

全色遙感圖像分割和分割過程中類別數的確定是重點問題。圖像分割方法主要分為基于閾值、區域、邊緣、聚類的分割方法等[1]，從數學角度分析，其將數學函數離散化或者將數字圖像連續化，進行數學意義上的分割，在物理意義上無法解釋部分算法中一階或二階不可導等問題[2]；從圖像本質分析，基于閾值的分割方法忽略了圖像的空間特征，基于區域的分割方法需要人為指定種子點來獲取全局最優解，基于統計學的貝葉斯模型，其對數字圖像建模，將每個像素值看作一個隨機變量，求其最大后驗概率分布，實現穩定可靠的圖像分割。圖像分割中類別數的確定多采用聚類方法，基于近鄰間的關系找出不同特征族群，對數據進行劃分，但易受聚類中心的個數及初始聚類中心的影響。基于遍歷的徑向基函數(radial basis function，RBF)網絡能夠實現非線性映射，其隱含層單元將非線性可分的輸入空間變換到線性可分的特征空間，輸出節點數代表分類數，但隱含層基函數的中心選擇和寬度很難準確選擇且參數的數量影響模型的復雜度，很難自動確定類別數。

綜合分析RBF網絡與非線性回歸模型的特點，本文利用貝葉斯理論將非線性回歸模型中基函數的參數視為具有某種已知先驗分布的隨機變量，根據參數的先驗分布建立圖像分割的后驗概率分布模型。

為計算后驗概率，一般方法是貝葉斯置信傳播、圖切等[3]。貝葉斯置信傳播用于有向無環圖模型的概率推理，也可用于無向圖模型，如馬爾科夫隨機場(Markov random field，MRF)。貝葉斯網絡中的先驗知識可能屬于同一集群，而極大似然估計在圖像分割中沒有考慮數據分組時的概率和空間約束問題，因此將MRF模型和貝葉斯置信傳播相結合[4]，求解最大后驗概率估計中的交互勢函數，可實現分割結果區域一致性，并且通過學習局部間的聯系得到全局最優解，但會夸大待求節點的邊緣影響，為保證貝葉斯網絡中節點的完整性，將貝葉斯網絡與多尺度理論結合的多尺度貝葉斯網絡模型[5]更加有效，而在模型尺度選擇上無法準確固定樹的層數且算法復雜度較高。

為此，本文結合GREEN[6]提出的可逆跳馬爾科夫蒙特卡洛算法(reversible jump Markov Chain Monte Carlo，RJMCMC)模擬該后驗概率分布，構建靈活、全面的RJMCMC混合轉移核，完成后驗概率估計。根據雅克比行列式思想構建適當的接受概率，實現不同參數的“維度-匹配”，使混合轉移核在不同維度參數空間之間跳轉，并采用模擬退火理論(simulated annealing，SA)來約束優化轉移核跳出局部最優達到全局最優，進而計算出徑向基函數的個數和參數，完成圖像類別數的確定和分割。在圖像分割前利用高斯曲率濾波(Gauss curvature filtering，GC)原理[6]對圖像進行幾何平滑預處理，進一步避免了模型優化過程中陷入局部最優，取得很好的分割效果。

本文運用經典的信息準則(Akaike information criterion，AIC)、貝葉斯準則(Bayesian information criterion，BIC)、最小長度描述法(minimum description length，MDL)和H-Q信息準則(Hannan Quinn criter，HQC)選出4種基于遍歷的RBF分割模型和4種分割算法，分別與RJMCMC+SA算法進行對比實驗，結果表明本文算法具有很好的復雜度和精確度且能夠自動確定圖像類別數。

1 圖像預處理及模型參數估計思想

1.1 GC圖像平滑處理需求

在非線性回歸模型中參數的選擇和優化始終存在局部極值問題，為使優化函數跳出局部最優解取得全局最優，KRIZHEVSKY等[7]提出了預訓練和微調理論。本文引入GC理論，在分割之前對圖像進行預處理。GC理論是一種具有邊緣保護的圖像平滑算法，充分利用了圖像的離散特征、微分幾何的連續理論，其假設圖像各塊曲面是分塊可展的，使用已知的幾何曲面來優化其對應的正則項，并利用圖像的離散特性來隱式地優化曲率，即減小曲率而無需計算曲率，避免了計算復雜的幾何流且濾波運算的復雜度大幅度降低，具有很高的執行效率。

圖像的GC處理中，利用3×3的滑動濾波窗口，平滑圖像中的像素點，如圖1所示。

圖1 高斯曲率濾波變分模型示意圖

GC根據鄰域像素構成的切平面中的最小距離修正調整像素值。對圖像中的極大值和極小值像素進行約束，降低像素的峰值、提高鞍點與谷值，保留其整體狀態，對圖像進行幾何優化，使分割優化過程中在學習率變化不大的情況下，盡可能避免陷入局部最優狀態。

1.2 后驗概率分布及其參數估計方法

一幅全色遙感圖像S，設為定義在圖像像素s上的標號場，其中，s為像素位置，Xs為定義于s上的隨機變量，為像素s的特征標號；設為定義在s上的特征場，其中，Hs為標號為xs的像素強度，分別為X，H的實現；定義為類別。設隸屬于類別k的像素標號集為如果Xk為空集，則對應類為空類；反之為實類。

多元參數映射模型采用Holmes和Mallick近似方案，構造非線性回歸模型M(model)，即

其中，D為參數矩陣化形式。尺度參數2σ服從共軛逆伽馬先驗分布系數α1:j服從方差為，期望為0的高斯分布；尺度參數δ2服從模糊共軛先驗參數超參數Λ服從無信息共軛先驗分布先驗模型階分布是一個有約束的截斷泊松分布。模型先驗分布為

對參數進行積分，得聯合后驗分布歸一化常數表達式為

2 圖像分割模型的建立過程

2.1 全色遙感圖像的標號場模型

假設Xs為獨立分布的，且滿足 Gaussian概率分布，則X的概率密度函數為

其中，aj為常數，表征s鄰域特征標號的作用強度；Ns為像素s的鄰域像素集合；I為指示函數，有

2.2 圖像模型的建立

對于全色遙感圖像，假設其像素強度Hs為獨立分布，且Hs滿足Gamma分布，則H的概率密度函數為

其中，()Γ·為Gamma函數；參數矢量式(8)刻畫出圖像中像素點的近似概率，結合式(2)，完備刻畫出圖像分割對象的服從概率，設整合參數向量則聯合后驗概率為

3 RJMCMC+SA算法原理及其實現

3.1 不同維度下參數空間的“維度-匹配”

聯合后驗概率(式(9))需對參數中非線性函數的高維積分進行評估，由于參數空間的維度不同，為解析地獲得不同維度下的參數空間計算，實現不同參數的“維度-匹配”，采用建議接受率在馬爾科夫蒙特卡洛算法(Markov chain Monte Carlo，MCMC)中構建一個非周期的馬爾科夫轉移核p(x,y)，定義一個特征分布π(x)，在細致平穩條件下采用Hastings算法[9]來解決轉移核和π(x)分布問題，則其建議接受率為

假設2個參數空間維度分別為m1和m2，在參數空間中存在一個奇異測度 π，則聯合平穩建議分布為為滿足細致平穩性和參數奇異自洽性，定義一個自由向量u和當前參數θ，在 2個子空間中不同子空間中式(10)改寫為

3.2 RJMCMC混合轉移核

RJMCMC混合轉移核由生成、刪除、分裂、組合、特征標號轉移組成。假設類別數k是變化不固定的，混合轉移都是自由移動，其分布分別為其中0≤k≤kmax。當k=0時，；k=1時，時，生成和刪除的接受率概率為

其中，p(k)為模型Mk的先驗分布，c是調整尺寸移動比例的常數。

同理，刪除概率移動為

接受率為

組合操作隨機選擇基函數中的一個μ1和其相鄰的μ2，根據Euclidean距離公式計算出2個參數的勒貝格測度保證其可逆性，分裂操作為

分裂概率移動為

同理，組合概率移動為

接受率為

特征標號轉移也就是RBF中心的更新，根據式(5)得特征場是一個基函數的全概率分布，即

3.3 RJMCMC+SA算法的實現

RJMCMC算法[10-11]能夠在參數空間中尋找到像素鄰域模型解空間的最優解，SA結合概率突跳特性在解空間中隨機尋找類別目標函數的全局最優解。

SA模擬非齊次馬爾可夫鏈[12]，狀態z在迭代i處的不變分布為為溫度下降進程且。在弱假設正則化下，根據Metropolis-Hastings(MH)算法，馬爾科夫鏈從狀態z轉移到狀態z′的過程中，建議分布為

假設可逆跳躍的齊次轉移核為,κ（z′ z），其狀態分布滿足，則RJMCMC+ SA的建議接受率為

RJMCMC+SA算法的實驗流程如下：

步驟1.遙感圖像初始化，利用GC算法對圖像進行平滑處理，設定圖像Hs,xs。

步驟2.初始化混合貝葉斯概率模型中的參數迭代次數i；設定均勻分布u：U[0,1]，溫度下降進程Ti。

步驟3.

則執行“生成”(式 13，15)；

則執行“刪除”(式 14，15)；

則執行“分裂”(式 16，18)；

則執行“組合”(式 17，18)；

else 更新徑向基函數的中心(圖像特征標號，式(20))；

end if。

步驟4.執行MH算法，SA運算(式(22))。

步驟5.迭代次數i←i+1，返回步驟3。

步驟6.計算徑向基函數的系數。

步驟7.通過步驟3～6計算出分割模型中圖像每個像素Hs符合類別X的徑向基函數；優化徑向基函數及其參數，取得X屬于Xk全局最優解空間，得出K個徑向基函數。

4 RJMCMC+SA算法評價

4.1 RJMCMC+SA算法的數據擬合性和復雜度

采用AIC，BIC，MDL和HQC選擇的4種基于高斯核函數訓練得到的 RBF圖像分割模型和GC預處理下的RJMCMC算法、RJMCMC+SA算法進行質量檢測對比，分析其復雜度和數據擬合能力(即似然函數)。模型評價策略依賴最大似然估計，對于一個模型Mk，其最優評價是似然估計與懲罰項的和為

其中，Θ為模型的估計參數；P為模型的懲罰項。AIC解決模型的擬合數據優良性和復雜性之間的平衡問題；BIC的懲罰項比AIC大，避免出現維度災難；MDL判別力求在模型精度和復雜度之間尋找平衡；HQC可避免小因素而忽略最優估計概率，其在貝葉斯方法中的漸進性比較好。定義ζ為模型的參數數量，4種模型選擇方法可表示為

實驗流程如下：

(1)初始化，即設定迭代次數5 000，kmax=50；

(2)迭代，即

隨機選擇中心點μ1:k，在測試集數據上計算丟失率(Loss)，更新模型，在訓練集數據上計算丟失率。

(3)根據不同的選擇方法來確定模型，即

丟失率(Loss)=均方誤差(mean-square error,MSE)

實驗選取伯克利大學實驗室圖像分割數據集BSDS500和 GF-2號全色遙感影像(譜段范圍0.45～0.90 μm、空間分辨率1 m)，不同分割方法下的實驗結果見表1和表2。分析表中數據，在普通圖像數據集上RJMCMC+SA算法相比AIC選出的最優RBF模型和RJMCMC算法，RBF的個數、數據損失量差別不大，且相對較低。在信息量大、細節豐富的GF-2號全色遙感圖像上，RJMCMC+SA算法相比HQC選出的最優RBF模型，RBF個數減少 250個，測試數據損失下降 3.09個點，且比RJMCMC算法在測試數據Loss和收斂速度方面得到進一步提升，說明在處理全色遙感圖像時RJMCMC+SA算法在數據擬合方面優于傳統 RBF網絡模型，同時表明SA在約束優化RJMCMC混合轉移核過程中，在基函數的中心點和個數方面取得優異效果。相比AIC、BIC、MDL方法通過增加懲罰項來實現數據擬合優良與復雜性的平衡，RJMCMC+SA算法在貝葉斯推斷中，其轉移核不僅在高低維參數空間跳變不斷調整后驗概率參數，而且進行了SA優化，在擬合基函數參數和基函數的選擇方面更加準確。對比分析表1和表2，在收斂速度上，RJMCMC+SA算法明顯快于其他5種算法，因為應用到模型集的核函數中心產生方式不同，雖都需要遍歷所有中心點，但前5種算法基函數中心的選擇是在μmin:μmax之間隨機產生，RJMCMC+SA算法是在μj:μj+1之間通過Euclidean距離計算得出，因此 RJMCMC+SA算法比其他算法收斂快，計算復雜度低。

表1 伯克利大學實驗室數據集上的實驗結果

表2 GF-2號全色遙感圖像上的實驗結果

4.2 RJMCMC+SA算法自動確定圖像類別數

采用本文算法對伯克利大學實驗室數據庫的一幅圖像和GF-2號全色遙感圖像進行分割，在圖像類別數穩定后，取算法前100次迭代數據，分析其對應的圖像類別數變化情況，如圖2所示。伯克利大學實驗室數據庫一幅圖像中迭代次數在 1～30時，類別數在0～7之間跳變，從30次以后穩定為4類且與人工標記的圖像類別數一致。GF-2號全色遙感圖像中迭代次數在1～50時，類別數在0～8之間跳變，從50次以后穩定為5類且與目視結果相同。綜合圖可知，圖像類別數可很快收斂到實際類別數。

圖2 不同迭代次數下圖像類別數變化

圖像分割完成前，在不同的迭代次數下，類別數不確定、分割區域不明確。RJMCMC+SA算法在確定類別數的同時還對分割結果進行了優化。對比圖中紅色實線區域分割結果，RJMCMC+SA算法分割下，伯克利大學實驗數據庫一幅圖像(圖(a1)～(b1))迭代次數分別為 10和15，全色遙感圖像(圖(a4)～(c4))迭代次數分別為10、15和30時，圖像類別數不穩定、分割區域不明確，而分別在30次和50次迭代之后，類別數穩定且分割結果明顯優于之前。因為 RJMCMC+SA算法在確定類別數的過程中，在SA溫度下降進程的約束下，執行分裂和合并操作，對RBF中心點進行特征標號移動操作，合并相似幅值的基函數，不斷更新RBF的中心點，加強了類別分割精確度。

4.3 RJMCMC+SA算法與閾值分割算法比較

選用迭代自組織數據分析算法(iterative selforganizing data analysis techniques algorithm，ISODATA)、GC+ISODATA算法、最大類間距算法、Color Slices算法和RJMCMC+SA算法進行比較，證明本文算法的優越性。

在圖2所示的實驗中，得出在迭代50次后分割結果趨于穩定，因此，在實驗中5種算法均迭代50次，完成全色遙感圖像分割，并對比分析，如圖4所示。ISODATA算法、GC+ISODATA算法分割結果模糊，抗噪聲差；最大類間距算法分割明顯，但是無法確定圖像類別數，不同類別的地物劃分為了一類；Color Slices算法能夠分割出多種不同地物，類別數過多，相似地物劃分為了多個類，導致分割雜亂，出現過分割現象，分割結果不精確；RJMCMC+SA算法相對以上算法效果明顯優異，地物目標分割清晰，噪聲少，分割精確度高，適合全色遙感圖像分割。

圖3 不同迭代次數下的分割結果

圖4 5種算法的分割結果

4.4 RJMCMC+SA算法定量評價

選擇ISODATA算法和本文算法進行實驗并求其混淆矩陣，作定量評價。ISODATA算法在徑向基函數中心點的選擇上：當圖像中屬于某個目標類別的像素點數過少時刪除此類別，當屬于此類別的像素點數過多且分散程度較大時則分裂這個類為兩個子類別。其在基函數中心選擇上與本文算法具有類似思想。

分析表3和表4所得的混淆矩陣[13]，RJMCMC+SA算法比ISODATA算法的用戶精度、產品精度高，且 Kappa系數和總精度分別高 0.49和36.7%，說明RJMCMC+SA算法在分類精度方面優于 ISODATA算法。采用 GC進行圖像平滑后ISODATA算法分類的Kappa系數和總精度分別提高0.18和13.9%，見表3和表5，因ISODATA算法在徑向基函數中心點(聚類中心點)的選擇上忽略了像素的空間特征，在類邊界處像素強度變化劇烈且被動擴散性影響下，容易陷入局部極值，GC平滑了像素點且抑制被動擴散，提高了算法的精度，說明在GC平滑濾波后提高了聚類中心點的選擇準確性。分析表4和表5數據，在經過GC處理后，RJMCMC+SA算法分割的 Kappa系數和總精度分別比ISODATA算法提高0.31和22.8%，表2實驗數據已然表明 RJMCMC+SA算法在基函數中心點選擇上相對4種模型的優越性，再次對圖像全局像素點峰值和谷值進行GC平滑預處理后，產品精度和用戶精度均大于85%且Kappa系數大于0.9，說明GC+RJMCMC+SA算法的分割精度更高。

表3 ISODATA算法分割的混淆矩陣

表4 GC算法預處理后RJMCMC+SA算法分割的混淆矩陣

表5 GC算法預處理后ISODATA算法分割的混淆矩陣

5 RJMCMC+SA算法主觀評價

AIC、BIC、MDL和HQC選擇的4種模型及GC+RJMCMC+SA算法分割結果對比。

觀察分析圖5(a1)～(f1)，針對細節較少、圖像復雜度低的伯克利大學數據圖像，4種模型和RJMCMC+SA算法均能實現顯著目標區域的分割，由于4種選擇方法的懲罰項不同，圖(b1)中AIC選擇的模型的分割效果相對于其他3種分割模型最優，藍色框中區域的噪聲最少，紅色框中的細節分割明顯，RJMCMC+SA算法亦能達到優異分割效果。在GF-2全色遙感圖像實驗中，4種模型的分割結果差異明顯，出現混分、錯分現象，如圖5(b2)～(e2)中飛機場停站樓及飛機(紅色實線邊框)區域，本文算法相比4種模型的分割結果具有低噪聲，清晰度高的優勢，更適合處理全色遙感圖像。

圖5 4種分割模型和RJMCMC+SA算法的圖像分割結果

圖(a2)～(d2)與圖(a4)～(d4)為 ISODATA 算法和RJMCMC+SA算法在10～50次迭代次數下對同一全色遙感圖像的分割結果，在同樣的迭代次數下，ISODATA算法相比RJMCMC+SA算法出現錯分和過分割現象，當迭代到 50次時本文算法分割效果明顯趨于穩定，分割區域清晰、準確，能夠清晰地分割出 3個油罐、廠房和道路，而 ISODATA算法的分割類別數混淆，精度差。在細節信息更為復雜的全色遙感圖像下，如圖(a3)～(d3)所示，ISODATA分割結果模糊無法分辨，且收斂緩慢、類別混亂不準確，細節分割錯誤；迭代次數同為50時，如圖3(d3)與圖5(f2)所示，RJMCMC+SA算法可清晰分辨出目標類別，說明本文算法分割的有效性和準確性高。

6 結束語

本文算法構建一種基于非線性回歸模型的RJMCMC+SA圖像分割算法。首先對圖像進行GC平滑處理，在幾何上避免了局部極值問題，然后利用 RJMCMC+SA算法實現貝葉斯形式化后的后驗概率分布，進而確定徑向基函數的參數和個數，完成全色遙感圖像中地物目標類別數的自動確定和分割。解決了非線性回歸模型中不同參數維度空間的跳變、模型冗余、模型可行性和復雜性不平衡、分割不穩定性、計算收斂速度慢等問題。相比傳統的分割算法，本文算法綜合像素強度、位置特征、參數相關性、空間特征，在處理信息量大、復雜度高的全色遙感圖像方面取得更大的優勢，更適合全色遙感圖像的分割。