白俊龍，王章瓊，張 明，李元松

武漢工程大學(xué)土木工程與建筑學(xué)院，湖北 武漢 430074

隨著無人機和傳感器技術(shù)的不斷成熟，無人機遙感系統(tǒng)逐漸成為空間數(shù)據(jù)獲取的一種重要手段［1-2］。與衛(wèi)星遙感和普通航空遙感技術(shù)相比，無人機遙感具有時效性強、機動靈活、獲取成本較低等優(yōu)勢［3］。此外，無人機還具備云下超低空飛行能力，彌補了傳統(tǒng)遙感技術(shù)因云層遮擋而無法獲取高質(zhì)量數(shù)字影像的不足，是傳統(tǒng)遙感手段的有力補充［4］，在災(zāi)害調(diào)查與解譯、監(jiān)測與評估、地質(zhì)勘察等方面發(fā)揮著極大的作用［5-6］。

現(xiàn)階段，無人機遙感數(shù)據(jù)的獲取能力已遠超其有效處理能力［7］，遙感圖像數(shù)據(jù)中所蘊含的大量有價值信息亟待充分挖掘與利用［8］。傳統(tǒng)的目視解譯往往需要專業(yè)的知識與豐富的經(jīng)驗，且因數(shù)據(jù)量巨大，導(dǎo)致解譯效率低下；而前沿的計算機解譯處理方式如機器學(xué)習(xí)，通常需要提供大量的訓(xùn)練樣本，否則難以保證后續(xù)的分類精度，導(dǎo)致信息提取難以完成［9-10］；此外，各相關(guān)領(lǐng)域能夠滿足模型訓(xùn)練條件的訓(xùn)練樣本數(shù)據(jù)，如滑坡樣本數(shù)據(jù)十分稀少等問題，也進一步增加了該方法的應(yīng)用難度。

K-means 算法是一種無監(jiān)督學(xué)習(xí)分類算法［11］，能對原始圖像進行聚類，再采用 Otsu 閾值法對聚類結(jié)果進行快速分割，結(jié)合形態(tài)學(xué)優(yōu)化等圖像處理，可實現(xiàn)無人機遙感圖像地物的精確提取與分類處理。鑒于此，本文根據(jù)實際工程需要，利用 MATLAB 強大的數(shù)值計算［12］、圖像處理功能［13］，借助其自帶的圖形用戶界面開發(fā)環(huán)境，設(shè)計了一個集圖像分割與分類等功能于一體的無人機遙感圖像分類系統(tǒng)，用于地物信息的快速識別與提取。

1 系統(tǒng)框架

圖形用戶界面可以內(nèi)嵌MATLAB 支持和用戶自定義的所有圖像處理函數(shù)，避免直接面對枯燥的代碼，通過具象化的圖形，如窗口、按鈕、菜單等簡便地操縱計算機［14］。MATLAB GUI 包含多個控件，如 Axes（坐標(biāo)軸）、Push Button（按鈕）、Pop-up Menu（彈出式菜單）等，用鼠標(biāo)拖動控件并雙擊即可對其屬性進行設(shè)置，與命令行界面相比，操作簡單，執(zhí)行效率高，因此，被廣泛應(yīng)用［15-16］。

基于MATLAB 的圖形用戶界面開發(fā)環(huán)境（GUIDE）設(shè)計無人機遙感圖像分類系統(tǒng)，該系統(tǒng)包含函數(shù)功能演示模塊和圖像系統(tǒng)處理模塊，其中函數(shù)功能演示模塊主要用于展示具體函數(shù)對圖像的處理效果，幫助用戶選擇合適的圖像處理函數(shù)，以期在處理不同圖像時均能達到較好的效果，在此對該模塊不進行過多的描述。

圖像系統(tǒng)處理模塊是整個系統(tǒng)的核心部分。該模塊用于對無人機遙感圖像進行系統(tǒng)化處理，最終實現(xiàn)圖像地物分類與信息提取。該部分主要包含圖像導(dǎo)入、圖像處理和圖像分類（圖1）。

圖1 圖像系統(tǒng)處理模塊結(jié)構(gòu)框架Fig.1 Structure framework of image processing module

1.1 圖像導(dǎo)入部分

使用imread 函數(shù)實現(xiàn)單張圖像的導(dǎo)入，通過imshow 函數(shù)實現(xiàn)數(shù)字圖像在系統(tǒng)中的實時顯示；若導(dǎo)入的為圖像文件夾，則先使用detectSURFF eatures 函數(shù)檢測和提取每張圖像上的Surf 特征點，根據(jù)特征點匹配實現(xiàn)圖像拼接，采用MSAC 算法對拼接之后的圖像進行幾何變換，生成最終的全景圖像用于后續(xù)處理。

若點擊“色彩轉(zhuǎn)換”按鈕可將RGB 圖像轉(zhuǎn)換到HSV、Lab 和YCbCr 色彩空間，通過改變特征顏色的深度和亮度實現(xiàn)圖像數(shù)據(jù)的增強，有利于后續(xù)圖像分割處理。

1.2 圖像處理部分

圖像分區(qū)處理使用了K-means 算法，根據(jù)顏色特征對圖像進行聚類分割，綜合程序運行時間和圖像聚類效果，在此將圖像聚類數(shù)設(shè)置為5。因此，分區(qū)處理可將圖像分為5 個不同區(qū)域，單個區(qū)域內(nèi)可能只包含一到兩種地物也可能包含多種地物。對于地物較少的區(qū)域，例如對絕大部分為“耕地”包含少部分“土路”的區(qū)域，可選擇“摳圖”操作，直接手工圈除“土路”部分，使該區(qū)僅保留“耕地”，然后對該區(qū)域圖像進行圖像增強、分割、去噪和形態(tài)學(xué)處理等一系列操作，優(yōu)化圖像成像效果，并最終生成和保存“耕地”區(qū)域的標(biāo)記，從而完成“耕地”區(qū)域的提取，下文會有具體的操作演示；若一個區(qū)域內(nèi)包含地物較多，則表明這幾類地物僅依靠顏色特征不易區(qū)分，此時可選擇手動獲取ROI（region of interest）感興趣區(qū)域，并生成相應(yīng)區(qū)域的標(biāo)記，即先用ginput 函數(shù)手動選取目標(biāo)區(qū)域，然后使用roipoly 函數(shù)實現(xiàn)感興趣區(qū)域的提取，最后用label2rgb 函數(shù)完成區(qū)域的標(biāo)記。

1.3 圖像分類部分

經(jīng)上述步驟完成各類地物的分割與提取后，通過對保存的標(biāo)記數(shù)據(jù)進行疊加，實現(xiàn)各類地物的合并，并在原圖上通過顏色區(qū)分展示具體的分類效果。此外，對每類地物進行標(biāo)注后，能夠更加清楚每類顏色具體代表哪類地物。各類標(biāo)記數(shù)據(jù)疊加部分的程序如下：

[filename, pathname, q] = uigetfile({'*.*'; '*.mat'},'select file','MultiSelect','on');

matFile=cell(length(filename),1);

for k=1:length(filename)

input=load(fullfile(pathname,filename{k}));

BJ=cell2mat(struct2cell(input));

matFile{k}=BJ*k*30;

L1=sum(cat(3,matFile{:}),3);

End

Lrgb=label2rgb(L1,jet,'w','shuffle');

end

2 模塊化設(shè)置

無人機遙感圖像分類系統(tǒng)采用自上而下的層次化設(shè)計思路，界面包含主界面和子界面。用戶通過主界面調(diào)用相應(yīng)的操作子界面，實現(xiàn)無人機遙感圖像分類處理。

函數(shù)功能演示模塊屬于輔助操作模塊，主要用于展示和對比各類圖像處理函數(shù)的具體效果，界面左側(cè)用于顯示各個函數(shù)的具體處理效果，右側(cè)為控制面板，可以選擇具體要查看和對比的函數(shù)。

圖像系統(tǒng)處理模塊屬于系統(tǒng)核心部分，主要用于對遙感圖像進行分類。界面包含圖像處理效果顯示區(qū)和操作面板2 部分，分別位于界面的上側(cè)區(qū)域和下側(cè)區(qū)域，下文即對該模塊的實現(xiàn)進行詳細描述。

2.1 圖像系統(tǒng)處理模塊

圖像系統(tǒng)處理模塊通過集成圖像拼接、圖像分區(qū)、分割、去噪、形態(tài)學(xué)優(yōu)化等算法，對無人機遙感圖像進行系統(tǒng)化處理、可視化操作，實現(xiàn)圖像的分割與分類。由圖1 可知，無人機遙感圖像的處理步驟主要包含圖像導(dǎo)入、圖像處理和圖像分類。

2.1.1 圖像導(dǎo)入部分 首先，即對圖像導(dǎo)入的選擇與操作，該部分可以選擇導(dǎo)入已經(jīng)拼接好的單張遙感圖像，節(jié)省操作時間，也可以直接選擇遙感圖像文件夾，進行拼接處理；圖像導(dǎo)入成功后可以根據(jù)成像效果和處理經(jīng)驗，選擇對圖像進行色彩轉(zhuǎn)換、幾何變換或顏色分量計算等操作。

選擇色彩轉(zhuǎn)換操作可將圖像由RGB 空間變換到 HSV 空間，生成H、S、V這 3 個分量，色調(diào)H代表顏色屬性，決定光譜主波長；飽和度S表示光譜主波長在強度中的比例，即顏色的鮮艷度；強度V表示光譜的亮度大小。圖像經(jīng)HSV 變換后，能夠提高結(jié)果圖像的紋理特征，增強其空間細節(jié)表現(xiàn)能力［17］，有助于遙感圖像理解和圖像分析，如沖溝［18］、沙地［19］等地物的識別與提取。此外還可將圖像由RGB 空間轉(zhuǎn)換到Y(jié)CbCr 空間或Lab 空間，用于特殊地物的提取如建筑物陰影［20］、棉花植被覆蓋度［21］的提取等。

幾何變化可將導(dǎo)入圖像進行平移、旋轉(zhuǎn)和鏡像處理，有助于用戶多角度查看地物，便于目視解譯。顏色分量的計算主要用于不同植被類型的識別，相對于衛(wèi)星與航空遙感圖像可以使用基于近紅外波段與紅波段的植被指數(shù)（如歸一化植被指數(shù)，NDVI）來指示綠色植物，而由傳統(tǒng)數(shù)碼相機拍攝得到的無人機遙感圖像往往由于缺少近紅外波段的數(shù)據(jù)，在植被提取方面只能另辟蹊徑，即借助可見光的紅光波段與綠光波段構(gòu)建顏色指數(shù)來提取植被信息［22-23］。本系統(tǒng)可計算的顏色指數(shù)主要包括過綠指數(shù)（excess green index，ExG）、歸一化綠-紅差值指數(shù)（normalized green-red difference index，NGRDI）和 綠 - 藍 差 值 指 數(shù)（green-blue difference index，GBDI）。

由此可知，色彩轉(zhuǎn)換、幾何變換和顏色分量計算這3 種操作，主要針對特殊地物的識別與提取，在對于絕大部分地物類型進行識別處理時，用戶可選擇跳過這些步驟直接對原始輸入圖像進行分區(qū)處理，提高處理效率。

2.1.2 圖像處理部分

1）K-means 聚類

K-means 算法是一種基于劃分聚類的無監(jiān)督學(xué)習(xí)算法，常以歐式距離作為衡量數(shù)據(jù)對象間相似度的指標(biāo)，對象間距離越小，相似度越大，核心思想為［24］：

步驟1：從給定的數(shù)據(jù)集中隨機選取k個初始聚類中心；

步驟2：計算其余數(shù)據(jù)對象與聚類中心的歐氏距離，根據(jù)歐式距離，將數(shù)據(jù)對象分配到與其相距最近的聚類中心所在的簇中；

步驟3：計算每個簇中數(shù)據(jù)對象的平均值作為新的聚類中心；

步驟4：進行下一次迭代，直到聚類中心不再變化或達到最大的迭代次數(shù)停止。

在圖像分區(qū)部分，K-means 算法可對圖像進行初分割。在使用K-Means 算法時，其聚簇個數(shù)k需要提前確定。因此，綜合程序運行時間和圖像聚類效果，本文將圖像聚類數(shù)k值設(shè)置為5，將圖像粗略分為5 個不同區(qū)域，其中單個區(qū)域內(nèi)可能包含不同數(shù)量地物類別。根據(jù)圖像處理效果即可以任選其中1 個區(qū)域也可選擇5 個區(qū)域之間的任意組合。如圖2 所示，區(qū)域4 和區(qū)域5 組合區(qū)域大部分為“植被”，包含小部分“水體”，表明該部分“水體”與“植被”僅依靠顏色特征已無法分割。由于“水體”區(qū)域相對較小，為提高處理效率可使用ginput函數(shù)，在“水體”周圍選擇多個點，用imfill 函數(shù)將此區(qū)域填充成塊，最終實現(xiàn)將“水體”部分去除，得到完整的“植被”區(qū)域。其中使用K-means 聚類實現(xiàn)圖像初分割的核心代碼如下：

圖2 圖像處理流程與效果Fig.2 Image processing flow and effect

global input;

Lab = applycform(input, makecform('srgb2lab'));

I = reshape(double(Lab(:,:,2:3)),size(ab,1)*size(ab,2),2);

nColors=5;

[cluster_idx,cluster_center] =kmeans(I,nColors,'distance','sqEuclidean','Replicates',2);

pixel_label=reshape(cluster_idx,nrows,ncols);

2）Otsu 閾值分割

閾值分割是一種基于區(qū)域的圖像分割技術(shù)，因其實現(xiàn)簡單、運算快捷、性能較穩(wěn)定而被廣泛應(yīng)用。

Otsu 閾值法是實際應(yīng)用中最為廣泛的閾值分割法之一［25］。

基本原理：假設(shè)一副圖像由前景色和背景色組成，首先根據(jù)計算方法（最大類間方差法）得到一個閾值T，若點（x，y）的灰度值f（x，y）＞T，則將其歸為對象點；反之，歸為背景點。由此可將圖像的前景色和背景色盡可能分開，實現(xiàn)圖像的二值分割。

經(jīng)過初分割得到的“植被”區(qū)域圖像，此時只包含目標(biāo)對象——“植被”和背景，使用MATLAB自帶的Otsu 函數(shù)能有效地將目標(biāo)對象與背景進行分割。

3）結(jié)果優(yōu)化與標(biāo)記

進行Otsu 閾值分割之后得到的二值圖像往往包含背景噪聲［圖3（a）］。中值濾波是一種非線性操作，使用medfilt2 函數(shù)對二值圖像進行中值濾波，可以濾除圖中的椒鹽噪聲，同時保留邊緣；使用imfill 函數(shù)對圖中的密集獨立的的孔洞區(qū)域進行填充使其連成一整塊區(qū)域，實現(xiàn)二值圖像的優(yōu)化處理［圖3（b）］。再使用bwlabel 函數(shù)對分割后的二值圖像進行處理，生成“植被”區(qū)域的標(biāo)簽數(shù)據(jù)，實現(xiàn)該區(qū)的信息提?。蹐D3（c）］。

圖3 圖像優(yōu)化與標(biāo)記：（a）Otsu 閾值分割，（b）二值圖像優(yōu)化，（c）標(biāo)記疊加Fig.3 Image optimization and marking：（a）Otsu threshold segmentation，（b）binary image optimization，（c）marker overlay

4）圖形矢量化

使用bwboundaries 函數(shù)尋找二值圖像的輪廓邊界線及對應(yīng)的坐標(biāo)，并用shapewrite 函數(shù)將對應(yīng)的坐標(biāo)信息和“Line”元素寫入形文件，對結(jié)果進行矢量化輸出，便于導(dǎo)入其它軟件進行處理，實現(xiàn)軟件之間的數(shù)據(jù)交互。

圖4 是將提取得到的“植被”標(biāo)簽矢量文件導(dǎo)入ArcGIS 處理得到的效果圖，在此可以任意更改圖層的顏色信息，添加“字段”為圖層賦予更多的屬性等。

圖4 ArcGIS 軟件處理效果Fig.4 Processing effect of ArcGIS software

2.1.3 圖像分類部分 重復(fù)2.1.2 步驟，對初分割得到的其它區(qū)域進行處理，完成圖中其它地物的提取，得到對應(yīng)地物的標(biāo)簽數(shù)據(jù)。選擇“多標(biāo)記疊加”按鈕會出現(xiàn)一個新的操作子界面（圖5），可對此前處理得到的各類地物標(biāo)記數(shù)據(jù)進行合并與類別標(biāo)注。

式（1）～式（4）中，Cl、Cw分別表示傳感器的長和寬；Tl、Tw分別表示圖像的長和寬；H表示無人機的航拍高度；b表示傳感器尺寸與圖像尺寸的比值；d表示圖像的比例尺；f表示相機焦距；m表示比例換算系數(shù)；li為單個像素在地面上的實際投影長度；n表示所求區(qū)域的像素點個數(shù)；S表示所求區(qū)域的實際換算面積。

如圖5（a）所示，圖中主要包含“水體”、“植被”、“道路”、“土體”和“其他”5 部分，通過對 5 類地物進行疊加，得到一個彩色標(biāo)記矩陣，其中相同顏色區(qū)域代表同一種地物。此外，選擇其中一個地物類型，如“水體”，輸入無人機的航拍高度H=300 m、相機焦距f=35 mm、傳感器尺寸Cl×Cw=35.9 mm×24 mm 和圖像尺寸Tl×Tw=6 000×4 000，根據(jù)式（4）［26］換算出該區(qū)域的實際面積［圖5（b）］，由此實現(xiàn)各區(qū)域面積的快速調(diào)查。表1 為根據(jù)以上輸入?yún)?shù)，計算出圖中水體、植被、道路和土體的實際面積。

圖5 系統(tǒng)子界面：（a）地物標(biāo)記疊加，（b）面積計算Fig.5 System sub-interface：（a）feature marker overlay，（b）area calculation

表1 各區(qū)域?qū)嶋H面積Tab.1 Actual area of each region m2

3 結(jié) 論

1）基于MATLAB/GUI 開發(fā)的無人機遙感圖像分類系統(tǒng)，包含函數(shù)功能演示模塊和圖像系統(tǒng)處理模塊兩部分。函數(shù)功能演示模塊為輔助模塊；圖像系統(tǒng)處理模塊為核心模塊，可對圖像進行拼接、分區(qū)、分割、去噪、形態(tài)學(xué)優(yōu)化等操作，最終實現(xiàn)無人機遙感圖像分類。

2）K-means 聚類與 Otsu 閾值分割相結(jié)合的方法，實現(xiàn)了圖像的快速分類與信息提??；分類處理后的結(jié)果可以矢量化輸出，大大提高了本系統(tǒng)與其他軟件的交互性。

3）根據(jù)無人機的航拍高度、鏡頭焦距和傳感器尺寸等信息，可計算每一類地物的實際面積，實現(xiàn)區(qū)域面積大小的快速統(tǒng)計。