高分辨率全色遙感圖像多級閾值分割

楊 蘊，李 玉，趙泉華

(遼寧工程技術大學 測繪與地理科學學院 遙感科學與應用研究所，遼寧 阜新 123000)

1 引 言

近年來，隨著遙感傳感器技術的進步，已發射了如Landsat,SPOT,IKONOS,QuickBird和WorldView等衛星，其對地面目標的采樣距離和重訪周期在不斷減小，獲得了大量像素大小為十幾米到亞米之間的高分辨率全色遙感圖像數據[1]。這些數據可以提供同一地區更多的地表詳實數據，便于更精確地實現地表觀測和地形測量。但高分辨率全色圖像包含信息的高復雜性對遙感圖像分析方法提出了更高的要求，為了有效利用這些信息，通常將圖像分割成一系列同質區域來過濾同一類別地物內部的異質性，因此圖像分割是高分辨全色遙感圖像分析的第一步和關鍵步驟[2]。現有全色遙感圖像分割方法可大致分為三類：基于區域、邊緣和閾值的方法[3]。其中，閾值法是最常用的方法之一，具有簡單易行、性能穩定等優點[4]。在該類方法中，精準分割的關鍵取決于最優閾值的選擇，在實際使用中，閾值的選取有多種方法，如直方圖法[5]、大津法[4]、最大熵法[6]、模糊聚類法[7]等。其中，最大熵法以其分割效果好、適用范圍廣等特點，具有廣泛的應用前景。

基于最大熵的兩級閾值圖像分割方法由Kapur等[6]提出，其用窮舉法選取使分割后圖像熵最大的閾值將圖像分為目標和背景兩部分。由于真實圖像通常包含多個目標，因此提出了如文獻[8-10]中的基于最大熵的多級閾值圖像分割方法，當圖像直方圖分布呈現明顯雙峰或多峰時，其能較好地實現背景和多目標的分離。然而圖像往往具有模糊性，造成直方圖分布復雜，使得很難獲得最優閾值。對此，一些學者將一型模糊集中的隸屬度和熵結合，構造了最大模糊熵來量化圖像的模糊性。如宋歡歡等[11]提出的基于模糊熵的自適應多閾值圖像分割方法，Chakraborty等[12]提出的最大模糊熵閾值法的快速算法。基于最大模糊熵的分割精度有所提高，但由于像素隸屬度是一個確定值，其提供的模糊特性不足，使得存在較多類屬和決策不確定性的問題[13]。區間二型模糊集是一型模糊集的推廣，它將隸屬度由固定值變換為一個區間，提供的自由度使得比一型模糊集具有更好的不確定性處理能力。近幾年，一些學者已將二型模糊集與熵相結合，用于在文獻[14-15]中獲得最優閾值，但都限于兩級閾值圖像分割。本文提出了一種基于區間二型模糊熵的多級圖像分割算法，在多級全色遙感圖像分割場景中定義區間二型模糊熵，利用一種搜索全局最優值的優化方法對區間二型模糊熵中的參數集進行優化，自動確定全局最優參數組合，從而達到確定最優閾值的同時獲得比窮舉搜索算法更少的時間開銷。

常用的優化算法主要包括遺傳算法(Genetic Algorithm, GA)[16]、粒子群算法(Particle Swarm Optimization, PSO)[17]和差分進化算法(Differential Evolution Algorithm, DE)[18]等，它們都是通過群體內個體之間的相互合作與競爭產生的群體智能來指導和優化搜索的方向，但當求解問題復雜時需要較大種群規模，計算量隨之變大，不利于求解。同時，隨著進化代數的增加，其種群的多樣性會變小，過早的收斂到局部最優點，使算法的性能在進化的過程中變差。為了解決上述問題，研究者將量子計算理論與GA相結合，提出了量子遺傳算法(Quantum Genetic Algorithm, QGA)[19]。它將潛在的解利用多態量子比特編碼為染色體，其在概率上表示所有可能狀態的疊加，因此在進化過程中可生成不同狀態的個體，能保持解決方案的多樣性，由此可克服傳統優化算法早熟收斂的問題。同時，量子旋轉門的旋轉運算，可充分利用當前最佳個體的信息進行下一次搜索以更新個體并避免停滯。與傳統優化算法相比，QGA具有種群規模小、種群多樣性好、適應性強等優點。

為了提高傳統模糊熵的分割精度，同時實現更好的時間性能，本文采用自適應QGA(Adaptive QGA, AQGA)對在多級全色遙感圖像分割場景中定義的區間二型模糊熵參數進行自動選擇。利用量子比特將模糊參數編碼為量子染色體，設置若干量子染色體構成初始模糊參數種群，并以定義的區間二型模糊熵作為評價函數，對種群中的個體進行適應度評價，保留和記錄最優個體。利用量子旋轉門實現量子染色體向最大適應度方向進化，并根據父代和子代量子染色體的差異和進化代數自適應地調整量子旋轉角的大小，以最終進化的種群中適應度最大的個體為最優參數，據此，由最大模糊性原則得到多級閾值，實現高分辨全色遙感圖像的最優多級閾值分割。

2 算法描述

2.1 圖像的概率劃分

設z={z(m,n), (m,n)S}為定義在圖像域S、大小為M×N的待分割圖像，其中，(m,n)表示像素位置，z(m,n)為像素(m,n)的光譜測度值，z(m,n){0,1,,L-1}，L為圖像的光譜測度量化級，zmin和zmax分別為圖像z的像素光譜測度的最小值和最大值。對于光譜測度為i，i[zmin,zmax]的像素點，其在圖像中出現的概率pi為：

pi=ni/M×N，

(1)

其中ni為圖像像素光譜測度為i的個數，滿足如下關系：

(2)

假設T1,T2, ... ,TC為z的C個閾值，滿足zmin=T0

2.2 多級區間二型模糊熵

因嶺型分布隸屬函數由升和降兩分布組成，因此對圖像C級閾值處理，需用C個閾值T1,T2,,TC來估計Sk,k{1,2,,C+1}的隸屬度μk，具有2×C個參數，即模糊參數集P={a1,b1,,aC,bC}，其中，0≤a1≤b1≤≤aC≤bC≤L-1，如圖1所示。

圖1 多級嶺型隸屬函數Fig.1 Multi-level ridge membership function

第k個隸屬函數可表示為：

(3)

其中，k{1,2,,C+1}，a0=b0=0，aC+1=bC+1=L-1。以μk(i)為基，構造區間二型模糊集，可表示為：

Fk=

(4)

其中：k{1,2,,C+1}，2πr和r分別為Fk的上和下隸屬函數，其用μk(i)構造，可表示為：

(5)

其中，λ>1為指數變化系數，根據經驗取λ(1, 5]，本文取λ= 3。圖2為λ=3時，μ1(i)的降嶺型分布變換示意圖，原隸屬度曲線如圖2(a)所示，變換后的隸屬度曲線如圖2(b)所示。

圖2 區間二型模糊隸屬度Fig.2 Interval 2 fuzzy membership degree

(6)

由像素光譜值及其概率pi，對Sk進行的概率劃分，其概率分布為：

(7)

其中:k{1,2,,C+1}，T0=zmin，TC+1=zmax。則Sk的區間二型模糊熵可定義為：

Ek=E(Sk)=

(8)

其中:k{1,2,,C+1}。對圖像C級閾值處理將圖像域S分成C+1個同質區域，因此，對含有C個閾值的區間二型模糊總熵的判別函數可表示為：

(9)

以圖像的區間二型模糊總熵最大為準則，搜索參數的最佳組合，確定最優參數集P={a1,b1,,aC,bC}，使真正代表目標和背景的信息量最大：

P*=Argmax(E).

(10)

由于模糊集的隸屬值越接近0或1，模糊性越小，若越接近0.5，模糊性越大，則按分割區域的最大模糊性原則，可用模糊參數由式(11)的方式獲得閾值：

(11)

2.3 基于AQGA的區間二型模糊熵算法

QGA在函數尋優方面具有計算量小，收斂速度快等特點，因此本文將AQGA用于式(9)中最優參數組合的確定，以實現圖像多級閾值分割。在尋優過程中，包含以下步驟：量子染色體編碼，量子染色體測量和解碼，適應度評價，量子染色體進化。

2.3.1 量子染色體編碼

QGA中的染色體用量子比特(Q-bit)進行編碼，一個量子比特可表示|0態和|1態及其之間的任意中間態，可表示為：

|ψ〉=α|0〉+β|1〉，

(12)

其中，α和β分別是態|0和|1的概率幅，且滿足以下歸一化條件：

(13)

其中，|α|2和|β|2分別表示態|0和態|1的概率幅。

本文的目的是將AQGA作為優化工具，尋求最優參數集P={a1,b1,,aC,bC}，以實現圖像多級閾值分割。因此對P進行編碼。假設第t代種群為,D}，其中D是種群規模，代表第t代種群中P的第K個潛在解個體，其在量子染色體的存儲和表達如下：

(14)

(15)

對應的染色體表示為：

(16)

則該染色體可分別以概率1/8,3/8,1/8和3/8表示4種狀態000,0001,100和101。由此可看出，通過制備不用的概率幅可以使得一條染色體同時描述多種狀態，使得潛在解具有多樣性，有利于最優解的尋找。

2.3.2 量子染色體測量和解碼

(17)

2.3.3 適應度評價

(18)

在最大適應度值所對應的量子染色體中隨機選擇一條作為進化目標，記為R的第h列的元素，h可表示為：

h=round(rand(1)×R)，

(19)

其中，rand(1)表示生成[0, 1]之間的一個隨機數，若rand(1) = 0，取h= 1。

2.3.4 量子染色體進化

本文采用量子旋轉門G，其使量子比特在單位圓內發生轉動，從而導致相位發生變化，由此增加種群的潛在多樣性。記第t代的量子門為G(t)，G(t) = {GK(t),K= 1, 2,,D}，其中，D為種群規模，,表示第t代第K個染色體的量子門，K{1, 2,,D}，{1, 2,, 2ZC}為作用于量子比特的旋轉門，可表示為：

(20)

其中，φK,q為種群中第K個染色體中第q個量子比特的旋轉角，其更新過程如下：

(21)

圖3 量子旋轉門調整量子比特Fig.3 Quantum bits are adjusted by quantum rotation

因此，第t代種群的更新過程可描述為：

Q(t+1)=G(t)×Q(t)，

(22)

其中，G(t)為第t代的量子門，Q(t)和Q(t+ 1)分別為第t代及其進化后第t+ 1代的種群。

(23)

2.4 算法流程

圖4為本文方法的流程圖，為了更好地理解其中的過程，對流程總結如下。

圖4 算法流程圖Fig.4 Algorithm flow chart

S1 統計待分割圖像中各個像素點的光譜測度值，計算光譜測度為i，i[zmin,zmax]的像素點在出現的概率pi，如式(1)所示；

S2 以一型模糊集中的多級嶺型模糊隸屬度函數(見式(3))為基，對原隸屬度曲線(見圖2)進行變換，構建區間二型模糊集，如式(4)所示；

S3 由構造的區間二型模糊集和閾值個數，在多級分割場景中定義模糊熵準則，如式(8)所示；

S4 對區間二型模糊總熵(見式(9))的參數集P={a1,b1,,aC,bC}進行編碼，初始化種群Q(t)，其存儲和表達如式(14)所示；

S5 對Q(t)中各個個體進行解碼為二進制串，得到疊加態的確定解M(t)，其基本原理如式(17)所示；

S6 將S3定義的總熵準則作為適應度函數(見式(9))，根據適應度大小評價各個個體的優劣，保留適應度最高的個體作為進化目標；

S7 判斷是否滿足結束條件，如果滿足則按照當前最大適應度對應的最優參數集合確定多級閾值(見式(11))，否則繼續計算；

S8 利用第t-1代量子旋轉角調整策略(見式(23))更新G(t)，并將其作用于Q(t)中的所有個體(見式(21)和式(22))，得到進化后的新種群Q(t+1)；

S9 對Q(t+1)執行步驟S5，S6和S7，使得各個體向目標方向進化；

S10 將算法執行到最大代數或父子代的最大適應度的差值小于0.01，否則返回步驟S7。

3 實驗結果與討論

3.1 實驗設置

本文通過MATLAB r2015b編寫程序，其運行環境為Windows7 32位專業版操作系統的，具有4G內存，處理器為Intel(R)Core(TM)CPU 32 G的計算機。為了驗證本文提出的區間二型模糊熵方法的可行性，利用文獻[8]中基于最大熵的算法和文獻[12]中基于一型模糊熵的算法作為對比算法，對其分割結果從目視分析、監督和非監督評價來進行比較分析。同時為了說明本文AQGA的優越性，以本文定義的區間二型模糊熵為目標函數，選擇窮舉搜索算法(Exhaustive Algorithm, EA)[21]，GA，PSO和DE作為搜索對比算法。由于各方法具有相同的目標函數，其分割結果相差不大，因此在論文中不對其分割結果進行展示，從搜索速度、精度和魯棒性來進行比較分析，其中用EA方法驗證對區間二型模糊熵的參數集進行優化的必要性；GA，PSO和DE表明在優化搜索方法中選擇本文AQGA的優勢。對于各圖像的定量評價，由于模擬圖像已知參考圖像，選擇監督評價中應用最廣泛的混淆矩陣來計算用戶精度、產品精度、總精度和Kappa系數對分割結果進行定量評價。對于遙感圖像，利用文獻[22]提出的面積加權方差(Area-Weighted Variance, WV)和Jeffries-Matusita(JM)距離對分割結果進行非監督評價。WV用于測量全局分割同質區域內的優度，并按每個同質區域的面積進行加權。其定義如下：

(24)

其中：k{1,2,…,C+1}為同質區域索引，C+1是分割的同質區域數，ak為k區域的面積，vk是k區域的方差。較低的WV表明分割結果具有較高的區域內同質性。JM用于測量分割同質區域間的整體優度，并按每個區域的面積進行加權，其定義如下：

(25)

其中Jk是波段區域k的JM距離，可表示為：

(26)

其中：Lk為區域k的邊界長度，Lkb為k與相鄰區域b{1, 2, …,}的共同邊界長度，Bkb是根據區域鄰接圖計算區域k和b的巴氏距離，可表示為：

(27)

其中，μk，μb和vk，vb分別是相鄰區域k和b的均值和方差。較低的JM表明具有較高的區域間異質性。

此外，為了驗證本文利用AQGA優化模糊參數集的優勢，利用GA,PSO和DE作為對比優化算法，并從尋優時間和最優參數與本文方法進行比較。參考文獻中對優化算法參數的描述，為GA,PSO,DE和AQGA設置最佳參數，如表1所示。

表1 GA，PSO，DE和QGA的參數設置

3.2 模擬圖像分割

通過尺度為256×256像素的標準灰度圖像(見圖5(a))對所提方法進行測試，假定其編號為I-V的區域內各像素相互獨立，且服從均值和標準差如表2所示的Gaussian分布，生成如圖5(b)所示的模擬圖像，其光譜測度值的分布如圖5(c)所示。

圖5 模板圖像和模擬圖像Fig.5 Template and simulated images

表2 模擬圖像的Gaussian分布參數

以圖5(b)為輸入，利用AQGA尋優構造的多級區間二型模糊熵參數，得到最優參數集P={8, 76, 79, 114, 118, 152, 167, 202}。由此根據式(9)得到分割區域的總熵E(P)=6.56，根據式(11)定義的方式得到閾值T={42,96.5,135,184.5}。與表2的標準數據和圖5(c)相比，得到的T在各區域均值間，可滿足最佳分割要求。以T對圖5(b)進行分割，結果如圖6(a)所示。為了從視覺上評估所提方法的精度，將分割結果的輪廓線(見圖6(b))與原始圖像疊加，如圖6(c)所示。

由各區域參數(見表3)和圖6的結果可以看出，盡管區域間的均值(如Ⅱ和Ⅲ)和標準差(如Ⅰ和Ⅴ)較相近，本文方法仍能夠對其分割出來，具有較好的視覺效果。

圖6 模擬圖像實驗結果Fig.6 Experimental results of simulated images

為了對分割結果精度進行統計分析，以圖5(a)為標準分割數據，計算與分割結果的混淆矩陣[23]，如表3所示。

表3 分割結果的混淆矩陣

據此分別計算出產品精度、用戶精度、總體精度和Kappa值，如表4所示，其中各項指標越高，表明其分割精度越高。

表4 產品、用戶、總精度和Kappa值

從表4可以看出，各區域的產品精度和用戶精度均達到較大值，總精度和Kappa值也達到0.99以上，一般分類器Kappa值達到0.8以上就是優質分類器[24]，因此表明本文方法的有效性和準確性。

3.3 遙感圖像分割

選取4幅具有不同地物的高分辨全色遙感圖像，如圖7(a)～(d)所示。其中，圖7(a)是分辨率為2.5 m，尺度為370 pixel×370 pixel的CARTOSAT-1農田圖像，人為判讀含有道路、農田和裸地3類；圖7(b)是分辨率為0.5 m，尺度為500 pixel×500 pixel的Worldview-1碼頭圖像，人為判讀含有水體、房屋、道路、裸地和植被5類；圖7(c)是分辨率為0.7 m，尺度為435 pixel×435 pixel的房屋圖像，人為判讀含有房屋、道路和操場等5類，圖7(d)是分辨率為0.5 m，尺度為600 pixel×600 pixel的海岸圖像，人為判讀含有房屋、棧道道路和沙灘等7類。各圖像光譜測度值的分布如圖8所示，從原始圖像和其直方圖的對比可以看出，直方圖的峰谷數與圖像中同質區域數并不一致，使得僅基于直方圖的分割方法失效。

圖7 全色遙感圖像Fig.7 Panchromatic remote sensing images

圖8 全色遙感圖像的直方圖Fig.8 Histogram of panchromatic remote sensing image

以圖像統計直方圖中各像素光譜測度值出現的個數為輸入，區間二型模糊熵最大為準則，通過本文方法尋優最佳閾值組合，結果如表5所示。以表中閾值組合和對比算法對圖7的全色遙感圖像進行分割實驗，其中，圖9(a1)～(d1)、圖9(a2)～(d2)和圖9(a3)～(d3)分別為采用文獻[8]、文獻[12]和本文方法所得的分割結果。由對比可以看出，傳統最大熵法在目標和背景光譜測度相差較小時，易導致誤分割現象，使得如圖9(c1)的操場和圖9(d1)的沙灘內部含有較多的分割噪聲。模糊熵法將隸屬度和熵相結合來量化圖像的模糊性，其對分割結果的精度都有一定程度的提升，但由于使用的是一型模糊集，像素隸屬度是一個確定值，因此提供的模糊特性不足，使得對如圖9(c2)和9(d2)的目標區域有較多欠分割現象，精度仍需要提高。本文方法利用二型模糊集構建模糊熵，每個像素點隸屬度是一個函數，其可提供較大的自由度來更好地處理不確定性，這在一定程度上減少了同質區域內部的分割噪聲。

圖9 分割結果Fig.9 Segmentation results

表5 最佳閾值組合和最大熵值

表6利用WV和JM對不同方法進行非監督評價的結果。從表中可以看出，本文方法在WV和JM方面均具有較小的值，說明本文方法相較于對比算法，分割結果的同質區域內部具有較小的差異程度，且各同質區域面積較為完整，相鄰同質區域間緊密相連，具有較大公共邊界，滿足圖像分割結果對空間連續且光譜均勻的要求。

表6 遙感圖像分割評價指標

為了能更加直觀觀察算法對評價指標的降低幅度，以表6中指標較大的文獻[8]的方法作為評價標準，計算文獻[12]和本文方法的定量評價指標降低幅度，如表7所示。其中，M0和J0為文獻[8]結果的MV和JM；M和J為文獻[12]和本文方法的MV和JM。由表中數據可看出，本文算法的兩個質量指標都得到較大的降低，MV平均降低提升了39.7%，JM平均降低了14.7%，表明了本文方法有較好性能，滿足算法魯棒性的要求。

表7 定量評價指標的降低幅度

以本文的區間二型模糊熵為目標函數，分別利用EA，GA，PSO，DE和本文AQGA對圖7(a)～(d)分割的模糊參數集進行優化，各算法的收斂曲線如圖10所示，其中各收斂圖的橫縱坐標分別為迭代次數和適應度值(區間二型模糊熵)。從圖10可以看出，本文AQGA相比其他三種算法，在收斂速度和全局收斂能力方面都具有一定的優越性。為了更客觀地評價算法的性能，對各算法運行20次的時間、精度和魯棒性進行比較，如表8所示，其中對應于以秒為單位的平均時間hm(s)，以比特/像素為單位的平均最佳目標函數值fm及其標準差std。

圖10 算法收斂曲線Fig.10 Algorithms convergence curve

表8 EA，GA，PSO，DE和AQGA的平均性能

由表8的對比可以看出，EA由于需要依次遍歷所有的元素，并判斷是否為可行解，因此用時較長，GA，PSO，DE和本文AQGA并不機械的檢查每一種情況，減少了搜索量，耗時相應減少。相比之下，AQGA的平均耗時最少，PSO與DE幾乎相同，GA耗時相對較多，對fm和std，AQGA對于所有圖像取得最大的fm和最小的std，這表明AQGA不僅能更好地優化問題，同時具有更強的魯棒性。

4 結 論

模糊熵法具有良好的性能，是一種主要的圖像閾值分割方法，但存在模糊特性不足、計算量大、自動性差等問題，使其在實際應用中受到限制。為了提高分割精度，本文利用嶺型模糊隸屬度函數構造的區間二型模糊集上定義圖像區間二型模糊熵，以增加其模糊特性，同時為了實現更好的時間性能，將自適應量子遺傳算法引入到區間二型模糊熵多級閾值圖像分割中，將其作為區間二型模糊熵中模糊參數集的優化工具，從而快速精確地確定圖像的多級閾值，實現圖像最優多級圖像閾值分割。通過對含有不同地物的高分辨率全色遙感圖像的實驗表明：與基于最大熵和模糊熵的方法相比，本文方法能在減少計算時間的同時獲得更好的分割結果，面積加權方差可降低39.7%，Jeffries-Matusita距離可降低14.7%，運行時間秒6.403 s。基本滿足高分辨全色遙感圖像分割結果對空間連續且光譜均勻的要求且具有較高的實時性。

對高分辨率全色圖像分割而言，由于地物種類繁多、背景環境復雜和實地考察困難等原因，造成人為確定其類別數比較困難。因此，在未來工作中，將對此問題進行研究，實現高分辨率全色遙感圖像的可變類分割。