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

白俊龍，王章瓊，閆海濤

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

0 引言

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

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

鑒于K-means算法的優點與缺點，結合無人機高分辨率遙感圖像的特點，本文提出一種K-means聚類引導的閾值分類方法。首先采用K-means聚類對圖像進行聚類處理，實現圖像的初始分割，其中利用Average Silhouette指標值來定量確定最優聚類數目k；然后對初始分割結果進行閾值分割和圖像優化處理完成對象的提??；再對所有提取結果進行合并實現全圖的識別與分類。并基于本文提出的分類方法處理步驟，采用MATLAB/GUI平臺，設計了適用于無人機遙感圖像分類處理系統。

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

1.1 K-means聚類

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

1.2 聚類數目k的確定

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

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

(1)

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

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

1.3 分類步驟

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

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

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

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

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

具體分類流程參見圖2。

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

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

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

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

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

2.1 原始圖像與k值確定

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

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

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

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

2.2 K-means聚類初分割

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

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

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

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

2.3 閾值分割與優化

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

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

圖9 優化二值圖像Fig.9 Optimized binary image

2.4 標記與合并

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

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

圖10 圖像處理系統子界面Fig.10 Sub interface of image processing system

3 結果與分析

3.1 精度評價指標

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

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

PA=xii/x+i，

(2)

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

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

UA=xii/xi+，

(3)

式中，xi+為第i行的總像元數。

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

(4)

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

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

(5)

3.2 精度評價

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

表1 混淆矩陣Tab.1 Confusion matrix

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

4 結論

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

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

3)對無人機遙感圖像進行分類處理，取得的OA為91.09%，Kappa系數為0.88，表明了該方法用于無人機遙感圖像分類處理，能夠實現地物的精確分類與信息提取。