K-means聚類引導的無人機遙感圖像閾值分類方法

白俊龍，王章瓊，閆海濤

(1.武漢工程大學土木工程與建筑學院，武漢 430073；2.中交第二公路勘察設計研究院有限公司，武漢 430052)

0 引言

無人機遙感技術以其時效性強、機動靈活、獲取成本較低等優(yōu)勢得以快速發(fā)展[1]，應用范圍也日趨廣泛，是傳統(tǒng)遙感手段的有力補充[2-3]。近年來，為充分利用無人機遙感影像豐富的信息數(shù)據(jù)，避免“數(shù)據(jù)豐富，信息貧乏”[4-5]現(xiàn)象的出現(xiàn)，越來越多的研究者開始了高分辨率遙感圖像的處理研究[6]。基于機器學習的圖像分類方法是現(xiàn)階段的研究熱點與重點[7]，但此類方法通常需要提供大量的訓練樣本來保證后續(xù)的分類精度[8-9]，且目前各相關領域能夠滿足模型訓練條件的訓練樣本數(shù)據(jù)十分稀少，導致該方法還難以用于實際工程應用。而無監(jiān)督聚類方法不需要訓練樣本，僅憑地物地磁輻射強弱在遙感圖像上所反映的光譜信息，即可將數(shù)據(jù)按其自然分布特性進行聚類，自動判別地物類別，并能克服人工解譯所帶來的主觀性因素影響[10]。

K-means算法是聚類算法中最為流行的算法，具有思想簡單、處理快速及效果好等優(yōu)點。自1967年MacQueen[11]首次使用后，K-means算法在聚類算法中得到了廣泛的應用。但在使用K-means算法時，其聚簇個數(shù)k需要事先確定，若k值選取不當會使最終的聚類結果陷入局部最優(yōu)。而在實際操作中，往往因為數(shù)據(jù)量過大和缺乏經(jīng)驗導致k值很難確定；且該算法對初始值十分敏感[12]，若k值選取得過小，則會導致同一簇內(nèi)數(shù)據(jù)對象差異很大，若k值選擇過大，則會導致不同簇間差異很小[13]。故在使用K-means算法進行聚類時，首先應該確定一個最優(yōu)聚類數(shù)目k。

鑒于K-means算法的優(yōu)點與缺點，結合無人機高分辨率遙感圖像的特點，本文提出一種K-means聚類引導的閾值分類方法。首先采用K-means聚類對圖像進行聚類處理，實現(xiàn)圖像的初始分割，其中利用Average Silhouette指標值來定量確定最優(yōu)聚類數(shù)目k；然后對初始分割結果進行閾值分割和圖像優(yōu)化處理完成對象的提?。辉賹λ刑崛〗Y果進行合并實現(xiàn)全圖的識別與分類。并基于本文提出的分類方法處理步驟，采用MATLAB/GUI平臺，設計了適用于無人機遙感圖像分類處理系統(tǒng)。

1 K-means聚類引導的閾值分類方法

1.1 K-means聚類

K-means算法是一種基于劃分聚類的無監(jiān)督學習算法[14]，具有簡單、快速、易于解釋、性能高，伸縮性好等多項優(yōu)點[15]。K-means算法通常以歐式距離作為衡量數(shù)據(jù)對象間相似度的指標，相似度與數(shù)據(jù)對象間的距離成反比，即對象間距離越小，則相似度越大[16]，其核心思想為[13]：①從給定的數(shù)據(jù)集中隨機選取k個初始聚類中心；②計算其余數(shù)據(jù)對象與聚類中心的歐氏距離，根據(jù)歐式距離，將數(shù)據(jù)對象分配到與其相距最近的聚類中心所在的簇中；③計算每個簇中數(shù)據(jù)對象的平均值作為新的聚類中心；④進行下一次迭代，直到聚類中心不再變化或達到最大的迭代次數(shù)停止。

1.2 聚類數(shù)目k的確定

由此可以看出，聚類個數(shù)k值選取是否得當會直接影響最終的聚類效果，故使用K-means算法進行聚類首先需要確定一個合適的k值。傳統(tǒng)K-means聚類過程中k值多依靠經(jīng)驗確定，聚類結果容易陷入局部最優(yōu)。本文利用Average Silhouette指標值來定量確定最優(yōu)聚類數(shù)目k，從而提高初分割精度與處理效率、減少后續(xù)手動處理步驟。輪廓圖(Silhouette)是一種用來刻畫聚類效果的度量，每個點的Silhouette值都表示相比于其他簇，這個點與本簇內(nèi)的其他點的相似程度。

對于樣本點i，其Silhouette值[17]Si的計算公式如下：

(1)

Average Silhouette是計算所有對象的Silhouette值的平均值，反映了分配聚類對象的整體效果，即聚類質(zhì)量[18]，根據(jù)Average Silhouette值可以評估分類結果的準確度或合適度。在聚類數(shù)目與Average Silhouette指標值的關系曲線圖中，Average Silhouette指標值往往會隨聚類數(shù)目的增加出現(xiàn)先增大后減小的變化趨勢，圖中Average Silhouette指標值最大的點所對應的聚類數(shù)目，則表示該數(shù)據(jù)集的最優(yōu)聚類數(shù)目k(圖1)。

圖1 聚類數(shù)目與Average Silhouette值關系示意圖Fig.1 Schematic diagram of the relationship between the number of clusters and the Average Silhouette value

1.3 分類步驟

本文提出K-means聚類引導的無人機遙感圖像閾值分類方法，主要包括以下步驟：

1)根據(jù)已有的無人機遙感圖像數(shù)據(jù)計算Average Silhouette指標值，確定當前圖像數(shù)據(jù)集的最優(yōu)聚類數(shù)目k值；對待處理圖像進行K-means聚類初分割；根據(jù)K-means聚類結果，目視判定對象類別，將同類別對象進行組合并手工剔除非目標區(qū)域，得到新對象。

2)對新對象進行閾值分割，將得到的二值圖像進行形態(tài)學等優(yōu)化處理。

3)重復步驟1)和2)直至完成所有區(qū)域的處理。

4)最后對處理完的所有區(qū)域進行合并，得到完整分類結果。

具體分類流程參見圖2。

圖2 K-means聚類引導的閾值分類流程Fig.2 The process of threshold classification under K-means clustering guidance

基于MATLAB/GUI平臺，對本文提出的分類方法處理步驟進行集成，開發(fā)了無人機遙感圖像分類處理系統(tǒng)，后續(xù)使用該系統(tǒng)完成無人機遙感圖像分類的全部處理，系統(tǒng)操作界面如圖3所示。

另外，我國有15個崩塌、滑坡和泥石流多發(fā)區(qū)，分別是：橫斷山區(qū)、黃土高原地區(qū)、川北陜南地區(qū)、川西北龍門山地區(qū)、金沙江中下游地區(qū)、川滇交界地區(qū)、漢江安康—白河地區(qū)、川東大巴山地區(qū)、三峽地區(qū)、黔西六盤水地區(qū)、湘西地區(qū)、贛西北地區(qū)、贛東北上饒地區(qū)、北京懷柔—密云地區(qū)、遼東岫巖—鳳城地區(qū)。

圖3 無人機遙感圖像處理系統(tǒng)界面Fig.3 UAV remote sensing image processing system interface

2 K-means聚類引導的圖像閾值分類

2.1 原始圖像與k值確定

采用高精度無人機對桂林至柳州高速公路進行拍攝，將獲得的高分辨率可見光遙感圖像作為原始數(shù)據(jù)源，航片的地面分辨率為0.04 m。工作區(qū)地勢西北高東南低，地表起伏較大，區(qū)域內(nèi)可分為構造侵蝕剝蝕低山丘陵、構造侵蝕剝蝕丘陵、侵蝕堆積河谷階地、溶蝕堆積孤峰平原4個地貌區(qū)。根據(jù)《土地利用現(xiàn)狀分類》(GB/T 21010—2017)[19]，該區(qū)地物類型主要包括林地、草地、耕地、交通運輸用地、水域及其他土地等多種土地利用類型。選取該工作區(qū)所拍攝的部分圖像進行拼接，生成待分類的無人機遙感圖像，如圖4所示。

圖4 拼接之后的無人機遙感圖像Fig.4 UAV remote sensing image after stitching

從采用無人機遙感系統(tǒng)獲取的高分辨率遙感圖像數(shù)據(jù)來看，針對同一批次的無人機遙感圖像數(shù)據(jù)集，無人機的飛行高度相對固定，單張圖像覆蓋的面積范圍有限，各張遙感圖像中所包含的地物類別數(shù)量相對較少且基本相等，故處理同一批次的遙感圖像數(shù)據(jù)，只需計算1個k值。通過對本次獲取的無人機遙感圖像數(shù)據(jù)集進行大量測試驗證，根據(jù)式(1)計算對應的Silhouette指標值，并根據(jù)聚類數(shù)目與Average Silhouette指標值的關系曲線圖(圖5)，可以確定一個最優(yōu)聚類數(shù)目k為4，故在使用K-means對圖像進行聚類時，不會陷入局部最優(yōu)，提高后續(xù)處理效率。

圖5 聚類數(shù)目與Average Silhouette值的關系Fig.5 Relationship between the number of clusters and the Average Silhouette value

2.2 K-means聚類初分割

首先，使用K-means聚類對原始圖像進行初分割。根據(jù)2.1節(jié)計算確定的最優(yōu)聚類數(shù)目，將初始聚類值設置為4，其初始聚類分割結果如圖6所示。通過目視判斷可知，初始分割可將部分地物較為完整地分離出來，如圖6(b)只包含“林地、草地”(本文將林地和草地按同一類型劃分)；絕大部分分割結果包含2類及以上地物，圖6(a)包含“其他土地”、“水域”和小部分“交通運輸用地”，圖6(c)包含“交通運輸用地”、“水域”和小部分“其他土地”、圖6(d)主要包含“林地、草地”和“水域”。

(a)分割結果1 (b)分割結果2 (c)分割結果3 (d)分割結果4

根據(jù)目視判定結果將圖6(b)和(d)進行組合，可得到新的處理對象。在系統(tǒng)中選用“摳圖”操作將新對象中的“水域”剔除，得到完整的“林地、草地”區(qū)域(圖7)，完成K-means聚類引導的初分割。

圖7 圖像組合與摳除Fig.7 Image combination and cutout

2.3 閾值分割與優(yōu)化

處理得到完整的“林地、草地”區(qū)域后，需對圖像進行閾值分割和進一步的優(yōu)化處理，生成二值圖像。使用可視化圖像閾值選擇GUI工具對“林地、草地”進行閾值分割，thresh_tool函數(shù)會先自動計算出一個自適應閾值，通常能對圖像進行很好地分割；若不滿意該結果，也可左右移動灰度直方圖中的閾值灰色豎線，調(diào)整閾值直至最佳分割效果。如圖8所示，對“林地、草地”進行分割的閾值為32。無人機遙感圖像分辨率高，在使用K-means聚類對兩種相互交融的地物(林地之間的其他土地)進行初分割時，難以得到精準的邊界，進行閾值分割之后得到的二值圖像往往包含背景噪聲。中值濾波是一種非線性操作，使用medfilt2函數(shù)對二值圖進行中值濾波，可以濾除圖中的椒鹽噪聲，同時保留邊緣。再使用imfill函數(shù)對圖中的密集獨立的的孔洞區(qū)域進行填充使其連成一整塊區(qū)域，實現(xiàn)二值圖像的優(yōu)化處理(圖9)。

圖8 圖像閾值選擇GUI工具Fig.8 GUI tool for image threshold selection

圖9 優(yōu)化二值圖像Fig.9 Optimized binary image

2.4 標記與合并

通過bwlabel函數(shù)對二值圖像中各個分離部分進行標注，返回一個和二值圖像大小相同的標記矩陣L，包含了二值圖像中每個連通區(qū)域的類別標簽。前面經(jīng)過K-means初分割和閾值分割之后得到的對象中只包含一種地物，故標記矩陣L中只有一種類別標簽，其對應的標簽值為1。使用label2rgb函數(shù)，可將標記矩陣L轉(zhuǎn)換為偽色彩圖像，至此即完成“林地、草地”的識別與標記提取。重復上述步驟，分別完成剩下K-means聚類結果的處理，最終從原始圖像中識別出“林地、草地”，“水域”，“交通運輸用地”和“其他土地”4類地物。

對4類地物標簽進行合并，得到彩色標記矩陣，將原始遙感圖像與彩色標記矩陣疊加即為最終的圖像分類結果。圖10為設計的無人機遙感圖像處理系統(tǒng)子界面。由圖可知，該界面除了對地物標簽進行合并與展示最終的分類結果外，還可以計算某類地物的像素面積大小及其占全圖總面積的比值。

圖10 圖像處理系統(tǒng)子界面Fig.10 Sub interface of image processing system

3 結果與分析

3.1 精度評價指標

混淆矩陣分析是進行圖像分類精度評價的常用方法[20]。本文使用混淆矩陣計算出的生產(chǎn)者精度、用戶精度、總體精度和Kappa系數(shù)作為評價指標，對無人機遙感圖像分類結果進行精度評價和分析。

1)生產(chǎn)者精度(producer’s accuracy,PA)。表示在分類中某一類別被正確分類的像元數(shù)與此類別真實參考像元之間的比率，公式為：

PA=xii/x+i，

(2)

式中：xii為i類地物正確分類的像元數(shù)目；x+i為第i列的總像元數(shù)量。

2)用戶精度(user’s accuracy,UA)。指正確分到某類的像元總數(shù)與所有被分為該類的像元總數(shù)之間的比率，公式為：

UA=xii/xi+，

(3)

式中，xi+為第i行的總像元數(shù)。

3)總體分類精度(overall accuracy,OA)。指被正確分類的像元個數(shù)與總像元個數(shù)的比值，公式為：

(4)

式中：n為總類別數(shù)；N為用于精度評估的像素總數(shù)。

4)Kappa系數(shù)。用于檢驗遙感圖像分類對于真實地物判斷的正確性程度，可以反映所分類別與遙感圖像之間的一致性。Kappa系數(shù)通常介于0～1之間，Kappa值越大表示分類精度越高。公式為：

(5)

3.2 精度評價

利用449個驗證樣本點數(shù)據(jù)，分別建立4種土地利用類型數(shù)據(jù)集在整體區(qū)域的混淆矩陣(表1)，計算出的OA為91.09%，Kappa系數(shù)為0.88。可見，使用本文方法對無人機遙感圖像分類得到的OA和Kappa系數(shù)都很高，說明分類結果與地表真實信息的總體一致性程度高，該方法適用于無人機遙感圖像的分類處理。

表1 混淆矩陣Tab.1 Confusion matrix

根據(jù)表1可知，用戶精度：水域>交通運輸用地>林地、草地>其他土地；生產(chǎn)者精度：交通運輸用地>水域>林地、草地>其他土地?？梢园l(fā)現(xiàn)其他土地利用類型的UA和PA均為最低，林地、草地次之，分析其原因，該類型分類精度低是因為緊挨林地、草地的其他土地對象中混入部分林地、草地所致；水域和交通運輸用地的兩種精度都很高，是因為水域多呈面狀或帶狀分布，交通運輸用地呈帶狀分布，區(qū)域內(nèi)地物分布均勻、單一，且這兩種形狀與其他地物接觸通常有較為明顯的邊界，易于識別區(qū)分。

4 結論

1)本文提出了K-means聚類引導的閾值分類方法用于無人機遙感圖像的分類處理，并基于Matlab/GUI平臺，對該方法處理步驟進行集成，開發(fā)了無人機遙感圖像分類處理系統(tǒng)。

2)根據(jù)Average Silhouette指標值與聚類數(shù)目之間的關系，確定出采用K-means聚類對無人機遙感圖像數(shù)據(jù)進行處理的最優(yōu)聚類數(shù)目為4。

3)對無人機遙感圖像進行分類處理，取得的OA為91.09%，Kappa系數(shù)為0.88，表明了該方法用于無人機遙感圖像分類處理，能夠?qū)崿F(xiàn)地物的精確分類與信息提取。