白俊龍，王章瓊，閆海濤

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

0 引言

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

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

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

1 K-means聚類引導(dǎo)的閾值分類方法

1.1 K-means聚類

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

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

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

對于樣本點(diǎn)i，其Silhouette值[17]Si的計(jì)算公式如下：

(1)

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

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

1.3 分類步驟

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

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

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

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

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

具體分類流程參見圖2。

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

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

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

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

2 K-means聚類引導(dǎo)的圖像閾值分類

2.1 原始圖像與k值確定

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

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

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

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

2.2 K-means聚類初分割

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

3.1 精度評價(jià)指標(biāo)

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

1)生產(chǎn)者精度(producer’s accuracy,PA)。表示在分類中某一類別被正確分類的像元數(shù)與此類別真實(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)。指被正確分類的像元個(gè)數(shù)與總像元個(gè)數(shù)的比值，公式為：

(4)

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

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

(5)

3.2 精度評價(jià)

利用449個(gè)驗(yàn)證樣本點(diǎn)數(shù)據(jù)，分別建立4種土地利用類型數(shù)據(jù)集在整體區(qū)域的混淆矩陣(表1)，計(jì)算出的OA為91.09%，Kappa系數(shù)為0.88?？梢姡褂帽疚姆椒▽o人機(jī)遙感圖像分類得到的OA和Kappa系數(shù)都很高，說明分類結(jié)果與地表真實(shí)信息的總體一致性程度高，該方法適用于無人機(jī)遙感圖像的分類處理。

表1 混淆矩陣Tab.1 Confusion matrix

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

4 結(jié)論

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

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

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