基于局部空間信息的可變類模糊閾值光學遙感圖像分割

楊 蘊 李 玉 趙泉華

隨著遙感圖像分辨率的提高,其在人們生活中的使用比例越來越大,如何快速、精確地從其中提取目標信息已成為遙感領域的一個重要研究課題[1].分辨率的提高在帶來信息豐富的同時也給目標分析帶來巨大困難,而圖像分割技術通過給每個像素分配標簽,使同一標簽的像素具有相同的視覺特性,可提高圖像分析的效率[2].因此執行圖像分割用于更好描述圖像是由圖像處理到圖像分析的關鍵步驟,對于遙感圖像目標的高效分析具有重要意義.

現有圖像分割方法大致可分為四類:基于區域[3]、邊緣[4]、學習分類[5]和閾值[6]的方法.其中,區域法的優點是對噪聲不敏感,但需要制定合適且高效的相似性準則,易造成圖像過分割;而單純基于邊緣檢測的方法往往不能提供較好的區域結構;學習分類法利用大量的訓練樣本能較好地解決圖像中的噪聲和均勻問題,但目前存在網絡模型難以確定,易出現欠學習和局部最優等問題[7-8];閾值法具有簡單易行、性能穩定等優點[9],目前流行的閾值方法有Otsu 法[10]、最大熵法[11]和聚類法[12]等,這些方法基于不同準則來選取閾值,都具有各自的優點,但也有自身的缺點.Otsu和最大熵法分別以分割后圖像的類間方差和熵最大為目標函數,窮舉選取使目標函數獲得最大值的閾值,但其僅考慮像素本身對分割結果的影響,對噪聲和異常值敏感[13].為此,人們引入鄰域像素進行二維閾值化,在一定程度上減弱了噪聲和異常值對分割結果的影響,但存在信息損失嚴重、耗時過長等缺點,特別是在多閾值圖像分割中更為突出[14].聚類法是一種不依賴于分類的先驗知識,遵循同質性準則對像素劃分的方法,能夠有效分離具有全局特征的數據[15],其中最具代表是Kmeans和模糊C 均值(Fuzzy C-means,FCM)算法.Kmeans 算法以歐氏距離作為相似度測度,以誤差平方和為聚類準則函數,迭代計算使得準則函數收斂為止[16].它具有簡單、計算快速等優點,但其強制圖像的各個像素屬于特定的類別,使得對具有復雜場景的遙感圖像的分割結果很差.FCM 方法利用模糊集理論中的隸屬函數對每個像素點進行劃分,考慮了模糊性和不確定性,可對紋理和背景簡單的圖像取得良好的效果[17],但FCM 僅考慮像素光譜信息,對噪聲較敏感.為此常將圖像局部空間信息引入目標函數,如文獻[18-19] 中,通過迭代計算像素鄰域和聚類中心之間的距離來提高FCM 算法對圖像分割的魯棒性.然而,局部空間信息的引入會導致算法復雜度升高,降低了該算法的實際應用價值.

為了克服傳統閾值法需人為設定類別數、對噪聲敏感、計算復雜度高等局限性,提出了一種基于局部空間信息的可變類模糊閾值遙感圖像分割方法.首先,以圖像中待定類像素點的光譜一階矩和二階矩為初始類中心和自適應閾值,利用二分法原理來初始化類別數目.對于任意待定類別,當相鄰兩次迭代得到類別的像素光譜均值差小于規定值時,該類查找結束,并以此均值作為該類的中心.并在確定各類的迭代過程中,計算與其二階鄰域系統間的相似度,以最大相似度準則進行類別區域合并.由確定的類別數和類中心為輸入,通過嶺形模糊隸屬函數計算圖像中各像素點對于不同類別的隸屬程度.考慮圖像中屬于同一類別的像素具有一定程度的連通性,即各像素與其鄰域之間關系緊密,不可能存在孤立的像素,而隸屬度和標號場是各像素屬性的映射,因此也滿足此局部空間關系.為了利用此局部空間信息同時降低計算復雜度,在隸屬度域中定義一個加權模糊濾波器對每一類別中像素點的隸屬度進行濾波,以濾波后的隸屬度集合為依據,按照最大隸屬原則確定圖像的標號場.最后,對標號場中局部異常標號進行替換,將修正后的標號場由對應的類中心賦色得到分割圖像.

1 算法描述

設I={Ip(i,j):(i,j)∈S,i=1,2,···,M,j=1,2,···,N}為定義在圖像空域S、大小為M×N的待分割遙感影像,其中,為像素光譜測度矢量,p為像素索引,(i,j)表示像素p的格點位置,l為光譜測度矢量的維度.對于全色遙感圖像,l=1;對于多光譜遙感圖像,l ≥3.

1.1 類別及其中心初始化

利用圖像的全局信息快速地對分割圖像的類別數及其中心粗略估計,其包含以下步驟:類中心初始化、劃分閾值選取、待定類識別、類別區域融合和查找結束.

1.1.1 類中心初始化

為了避免隨機選擇初始值而降低類別查找的動態性,本文以圖像光譜的一階矩作為初始類中心,即

其中,s∈{1,2,···,l}為波段索引,s為波段 的初始類中心.

1.1.2 劃分閾值選取

以圖像中待定類各像素點與類中心的二階中心矩的算術平方根為自適應閾值,較之各待定類像素點與類中心的距離從而確定類屬劃分,波段s中第k個類別的自適應閾值為

其中,s∈{1,2,···,l},k ∈{1,2,···,K}為類別索引,K為確定的類別數目,M′×N′為確定第k類中心時圖像中待定類像素點個數,則第k個類別的自適應閾值為

1.1.3 待定類識別

當查找第k類時,以歐氏距離為準則,計算波段s中各待定類像素點與類中心數值特征的絕對差異:

1.1.4 類別區域融合

考慮到類別查找時,由于遙感圖像的局部非均勻性,在區域上會出現局部最優現象,為此通過計算相鄰區域的相似度[20]來指導合并過程,確定類別數.其將圖像每一維度均勻量化為16 個等級,然后在特征空間中計算每個類別區域的歸一化直方圖,使用Bhattacharyya 系數[21]來度量相鄰類別區域間的相似性,可表示為:

其中,l為光譜測度矢量的維度,H istR和H istQ分別是類別區域R和Q的歸一化直方圖,上標u代表它們的第u個元素.設以R為中心的二階鄰域系統內所有類別區域為Qk,k ∈{1,2,···,8},計算相似度ρ(R,Qk),取其最大值對應的類別區域索引,可表示為:

其中,z∈{1,2,···,8}.若ρ(R,Qz) 大于閾值,將R與Qz合并,并更新類別及其中心,其中閾值越大,合并條件越嚴格,反之亦然,本文取0.85.

1.1.5 查找結束

當第k類確定后,在待定類的像素點中重復執行上面步驟,直到所有像素完成標記,類別查找結束.

1.2 像素隸屬度分配

定義一個模糊隸屬函數μk(·),由初始化的類別數和類中心使各像素與每個類相關聯,得到各像素所屬類別的隸屬度.嶺形分布隸屬函數(見圖1)采用正弦函數作為核心,具有主值區間寬、過渡平穩的特點[22],能很好反映各像素與類別隸屬程度的模糊關系,定義為:

圖1 多級嶺形隸屬函數Fig.1 Multilevel ridge membership function

其中,a,b,c為模糊參數,將其與類中心和圖像相關聯,當k=K=1 時,a=min(Is),b=CK;當k=K時,a=Ck-1,b=Ck,c=max(Is);當 1＜k ＜K時,a=Ck-1,b=Ck,c=Ck+1 .

其中,K為類別數.圖像中所有像素的隸屬度向量組成一個隸屬度集μ(I).

1.3 隸屬度域濾波

同一類別區域內部在分割過程由于隨機噪聲產生的孤立像素可認為是誤分割,在隸屬度矩陣中表現為異常的隸屬度值.為了利用此局部空間信息同時降低計算復雜度,將模糊加權函數[23] 應用于均值濾波器移動窗口內的隸屬度值,在隸屬度域中定義一個模糊加權均值濾波器來過濾隨機噪聲,以μk(i,j),k ∈{1,2,···,K}為二維模糊濾波器的輸入,其輸出為

其中,(m,n)∈ηp,μmax(i,j)、μmin(i,j)和μmav(i,j) 分別是ηp內輸入值的最大、最小和平均值.

1.4 標號域濾波和賦色

由濾波后的隸屬函數,按最大隸屬原則得到標號場,如式(10)所示

其中,k∈{1,2,···,K},Ls={Ls(i,j)=k:i=1,2,···,M,j=1,2,···,N}.為了消除標號場Ls中由于誤判的類別標號(假設表現為脈沖噪聲)并保持其邊緣,以移動窗口內的中值構建窗口加權函數WM(·) 來定義模糊加權中值濾波器,如式(11)所示

1.5 算法流程

算法流程圖如圖2 所示.

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

綜上所述,可對本文方法的流程總結如下.

S1 利用二分法原理和區域間最大相似度準則對圖像的類別數及其中心初始化;

S2 將類中心和各像素與嶺形分布隸屬函數相關聯,得到各像素所屬類別的隸屬度;

S3 在隸屬度域中定義一個模糊加權均值濾波器來過濾由于噪聲引起的異常隸屬度值;

S4 對濾波后的隸屬度集按照最大隸屬原則確定圖像的標號場,并由模糊加權中值濾波器來消除其中由于類別初始化時誤判的類別標號;

S5 對修正后的標號場由對應的類中心賦色得到分割圖像.

2 實驗結果與討論

選擇模擬圖像、全色遙感圖像和多光譜遙感圖像對本文方法進行驗證,其中研究區的真彩色圖像是由WorldView-2 的紅、綠、藍三波段(波段5=紅,波段3=綠,波段2=藍)合成.

2.1 實驗設置

實驗是在具有Intel(R) Core (TM) CPU,i7-6 700,3.4 GHz,16 GB RAM 的DELL 計算機上,應用Windows 8 環境下的MATLAB R2017a 實現的.為了驗證提出方法的優越性,將所提出方法與一些著名的閾值算法的性能進行比較,選擇了Kmeans、FCM和文獻[19]的方法作為對比算法,其中Kmeans和FCM 算法可在文獻[16]和文獻[17]中得到詳細描述,對其分割結果從目視分析、監督評價、非監督評價和運行時間來進行比較分析.對于模擬圖像,由于已知參考圖像,選擇監督評價中應用最廣泛的混淆矩陣來計算用戶精度、產品精度、總精度和Kappa 系數對分割結果進行定量評價.對于全色遙感圖像,利用文獻[24] 提出的面積加權方差(Areaweighted variance,WV)和Jeffries-Matusita (JM)距離對分割結果進行非監督評價.對多光譜遙感圖像,其光譜信息增加,為了較全面評價分割質量,除了利用WV和JM 外,增加文獻[25]提出的可見色差E對分割結果進行非監督評價.WV用于測量全局分割同質區域內的優度,并按每個同質區域的面積進行加權.其定義如下

其中,i∈{1,2,···,n}為同質區域索引,n是同質區域數,ai為區域i的面積,vi是區域i所有波段的平均方差,l是遙感圖像的波段數,vib是b波段中第i個區域的方差.較低的WV 表明分割結果具有較高的區域內同質性.JM 用于測量分割同質區域間的整體優度,并按每個同質區域的面積進行加權,其定義如下

其中,Ji(b) 是波段b中區域i的JM 距離,可表示為

其中,Li為區域i的邊界長度,Lik為區域i及其相鄰區域k∈{1,2,···,Ω}的共同邊界長度,Bik是根據區域鄰接圖計算區域i和k的Bhattacharyya 距離,可表示為

其中,μi,μk和vi,vk分別是相鄰區域i和k的均值和方差.較低的JM 表明分割結果具有較高的區域間異質性.

在可見色差的情況下,E結合區域內和區域間視覺誤差來評價分割結果的質量,定義為

其中,Eintra表示區域內視覺誤差,其旨在評估欠分割的程度,可表示為

其中,M×N為圖像大小,I和I′分別表示原始和分割圖像,‖·‖Lab表示CIELab 顏色空間中的色差,th表示可見色差的閾值,u(·) 表示階躍函數,

其中,th為可調閾值,本文th=6 .Einter為區域間視覺誤差,用來評估過分割程度,定義為

其中,K表示分割的同質區域數量,wij表示區域i和區域j之間的連接長度,如果區域i和j未連接,則等于零,c表示歸一化因子,本文選擇c=1/6.由于Eintra與Einter之間存在著相互關系,因此當E達到較小的值時,被視為實現較好的分割.

2.2 模擬圖像分割

圖3 (a1)和圖3 (a2)是尺寸為 1 28×128 像素的模板,用于生成標記為I-V 的五個均勻區域的模擬圖像(見圖3 (b1)和圖3 (b2)),其各波段中均勻區域的像素強度服從表1 所示的高斯分布.

表1 各同質區域的高斯分布參數Table 1 Gaussian distribution parameters of homogeneous regions

圖3 模擬圖像Fig.3 Simulated images

分別利用對比算法和本文方法對圖3 (b1)和3 (b2)進行分割實驗,分割結果如圖4 所示,其中第一行和第二行分別為圖3 (b1)和圖3 (b2)的分割結果.由圖4 可以看出,Kmeans 的分割結果中含有大量分割噪聲,這是因為其僅考慮像素的測量,對光譜的變化很敏感,因此很難區分不同的均勻區域.對于FCM,它利用隸屬度對像素模糊劃分,隸屬度提供的自由度使得比傳統硬閾值劃分具有更好的不確定性處理能力,從而減弱噪聲和異常值對分割結果的影響,它可以成功地區分區域I、III和IV,但區域II和V 的分割結果比較模糊.文獻[19] 通過定義中心像素與鄰域像素隸屬度間的制約關系改善了FCM 得到的分割結果,但對部分區域不能有效分割,如圖4 (c2)的區域II.本文方法在正確識別類別數的情況下,考慮了隸屬度空間中各像素隸屬度之間的空間信息,從而加強了鄰域系統中像素之間的連接,此外在標號域中建立的模糊中值濾波剔除了局部異常標號,取得了更好的分割結果.

圖4 模擬圖像分割結果Fig.4 Simulated image segmentation results

為了進行定量精度評價,計算不同同質區域的誤差矩陣,從而計算出產品精度、用戶精度、總體精度和Kappa 系數[26]來描述分割結果的準確性(見表2).從表2 可以看出,本文方法各同質區域的用戶精度和產品精度均在96 % 以上,總體精度和Kappa 系數平均為98.8 %和99.1 %.在對比算法中,Kmeans 的總體精度和Kappa 系數平均分別為54.3 %和51.4 %;FCM 的總體精度和Kappa 系數平均分別為78.1 %和77.3 %,其優于Kmeans;文獻[19] 的總體精度和Kappa 系數分別為89.1 %和87.5 %,優于Kmeans和FCM,但精度仍低于本文方法.定量結果表明,本文所提出的方法可取得良好的分割效果,具有廣闊的應用前景.

表2 模擬圖像分割的定量評價結果Table 2 Quantitative evaluation results of simulated image segmentation

2.3 全色遙感圖像分割

圖5 (a1)～ (d1)為待分割的5 幅全色遙感圖像,其中,圖5 (a1) 是分辨率為0.5 m,尺度為322×288 像素的Pleiades-1 機場圖像,人為判讀類別數為4;圖5 (b1) 是分辨率為1 m,尺度為270×350 像素的IKONOS 郊區圖像,人為判讀類別數為5;圖5 (c1)和(d1)是分辨率為0.7 m 的EROS-B 圖像,其中,圖5 (c1)為 5 12×512 像素的碼頭圖像,人為判讀類別數為6,圖5(d1) 為256×256 像素的房屋圖像,人為判讀其類別數為7.采用本文算法和對比算法對全色遙感圖像進行分割實驗,其中圖5 (a2)～ (d2) 為采用Kmeans 方法進行分割實驗得到的結果;圖5 (a3)～ (d3)為采用FCM 方法進行分割實驗得到的結果;圖5 (a4)～ (d4)為采用文獻[19]方法進行分割實驗得到的結果;圖5 (a5)～ (d5)為采用本文方法進行分割實驗得到的結果.

由圖5 可以看出,通過各方法的分割,圖像表現為不同的同質區域,簡化了圖像表示,有利于圖像的后續分析.對Kmeans 法,各同質區域疊加在一起,使其含有大量的分割噪聲,分割結果較差;FCM 法分割精度有所提高,但仍有一些明顯的噪聲;文獻[19]的方法由于引入局部信息,具有一定的抗噪性.由于在隸屬度域和標號域進行了模糊加權濾波,本文方法在視覺上提供了比對比算法更好的分割結果.

表3 為利用WV和JM 對圖5 (a1)～ (d1)進行定量評價的結果.從表中可以看出,本文方法均具有較小的WV和JM 值,說明本文方法相比對比算法,分割結果的同質區域內部具有較低的方差,且各同質區域面積較為完整,相鄰同質區域間具有較大公共邊界,滿足圖像分割結果對空間連續且光譜均勻的要求.

表3 全色遙感圖像分割質量評價指標Table 3 Quality evaluation of panchromatic remote sensing image segmentation

圖5 全色遙感圖像和分割結果Fig.5 Panchromatic remote sensing images and segmentation results

由于場景的不同,屬于每一類的像素不盡相同,存在稀疏類、密集類和均勻類[27].為了便于分析,本文假設各場景為均勻類,其像素在各類別中均勻分布.基于此,在表4 中給出了對比算法和本文方法的計算復雜度,其中M×N是圖像的像素數,K是類別的數量,t是迭代數,ω是濾波窗口的大小,M ×N/K為每一類包含的像素數.

根據表4,由于每個類包含M×N/K個元素,Kmeans 在每次迭代中最多進行K次比較就可以確定適當的聚類,并進行M×N/K比較以確定適當的元素,具有低計算復雜度,文獻[19]的方法與迭代次數成線性關系且每次迭代時都需要計算像素鄰域值,具有較高的計算復雜度,本文方法考慮空間鄰域信息時利用在隸屬度和標號域濾波,降低了計算復雜度.此外,為了估計算法的實用性,我們比較了這些算法的運行時間,如表5 所示(以秒為單位).

表4 計算復雜度對比Table 4 Computational complexity comparison

由表5 可以看出,由于Kmeans和FCM 僅考慮像素光譜測度對分割結果的影響,忽略了局部空間信息,能夠在較短的時間內對圖像進行分割.文獻[19]在聚類的迭代過程中利用空間鄰域信息,需要在每次迭代中計算像素鄰域與聚類中心之間的距離,使其具有較高的計算復雜度,耗時較長且隨圖像尺度增大呈現非線性增長,相比之下,本文方法利用圖像的全局信息初始化類中心,耗時短,而在考慮局部空間信息時,本文在隸屬度域和標號域使用濾波方法,降低了計算復雜度,使本文方法整體計算時間較短.

表5 全色圖像分割時間對比(s)Table 5 Panchromatic images segmentation time comparison (s)

2.4 多光譜遙感圖像分割

為了體現本文方法的適用性,對多光譜遙感圖像進行了分割.圖6 (a1)～ (d1)顯示了0.5 米分辨率的WorldView-2 圖像,其含有不同土地利用和覆蓋類型.其中,圖6 (a1)～ (d1)的尺度分別為 4 76×432、397×356、4 45×342和5 29×435 像素,目視判讀分別含有5、6、6和7 個同質區域.利用本文方法和對比算法對其進行分割實驗,其中圖6 (a2)～ (d2)為利用Kmeans 方法所得的結果;圖6 (a3)～ (d3)為利用FCM 方法所得的結果;圖6 (a4)～ (d4)為利用文獻[19]方法 所得的結果;圖6 (a5)～ (d5)為利用本文方法所得的結果.

由圖6 可以看出,由于多光譜遙感圖像光譜信息豐富,具有更多的地物色彩,增加了分割的難度,Kmeans 算法不能對同質區域內部的異質性進行有效分離.FCM 由于考慮了模糊隸屬關系而能區分部分內部區域,但沒有克服其對噪聲的敏感性,分割仍然不能令人滿意.文獻[19]利用空間鄰域信息,計算每個迭代中像素鄰域與聚類中心之間的距離,提高了圖像分割效果,但固定的空間距離對不同局部信息的圖像不具有魯棒性,使得如圖6 (b4)的耕地內部出現誤分割現象.本文方法所采用的隸屬度空間濾波和相應的標號域濾波之間的相互作用可以增強同一鄰域系統中像素之間的聯系,可以在視覺上更好地捕獲同質區域,因此分割結果比文獻[19]方法更令人滿意.

表6 為利用MV、JM和E對圖6 中的4 幅圖像進行評價的結果.從表6 可以看出,本文方法的指標均是最優的,說明本文方法可有效分割出同質區域,且在各區域內部具有較高的純度,滿足人眼對各區域分離的色差.

表6 多光譜遙感圖像分割質量評價Table 6 Quality evaluation of multispectral remote sensing image segmentation

表7 為本文方法與對比算法對圖6 中的4 幅圖像的分割時間,并將其與表5 的全色遙感圖像運行時間進行對比.

表7 多光譜圖像分割計算時間對比(s)Table 7 Multispectral images segmentation time comparison (s)

圖6 多光譜遙感圖像和分割結果Fig.6 Multispectral remote sensing images and segmentation results

與全色遙感圖像分割的運行時間對比可以看出,由于多光譜遙感圖像數據維數和各均勻區域中的內部異質性的增加,使得計算復雜度增大,各方法對同尺度圖像的分割時間增加,較之不考慮局部空間信息的Kmeans和FCM,文獻[19]的耗時增加更明顯,而本文方法是在隸屬度域和標號域考慮局部空間信息,因此耗時短且魯棒性強,具有實用性.

3 結論

1)提出一種利用圖像全局信息進行類別及其中心初始化的方法,以圖像中待定類像素點的光譜一階矩和二階矩為初始類中心和自適應閾值,利用二分法原理和區域相似度準則對類別數及其中心粗略估計.

2)以估計的類別數和類中心為輸入,通過定義的嶺形模糊隸屬函數快速得到圖像中各像素點對于不同類別的隸屬程度,將圖像域轉化為隸屬度域.

3)考慮到隸屬度和標號場是各像素屬性的映射,首先在隸屬度域中定義一個模糊加權均值濾波器來過濾由于噪聲引起的異常隸屬度值,對濾波后的隸屬度集按照最大隸屬原則確定圖像的標號場.然后由移動窗口內的中值構建模糊加權中值濾波器,以消除標號場中由于誤判的類別標號.

4)在隸屬度域和標號域濾波時,由于原圖像分辨率、場景復雜度等的不同,造成濾波窗口鄰域大小難以自動獲取,其產生不同的平滑度,對最終分割結果中同質區域的邊界有較大影響.

綜上,本文利用圖像全局信息進行類別及其中心初始化可避免由于隨機選擇初始值而降低類別查找的動態性,后經嶺形模糊隸屬函數可快速對像素類屬劃分,同時由隸屬度域和標號域濾波,充分利用了圖像的局部空間信息以對類別數進行精確劃定,為可見光遙感圖像的閾值分割問題提供了一個比較實用的解決方案.