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

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

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

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

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

1 系統(tǒng)框架

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

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

圖像系統(tǒng)處理模塊是整個(gè)系統(tǒng)的核心部分。該模塊用于對(duì)無(wú)人機(jī)遙感圖像進(jìn)行系統(tǒng)化處理，最終實(shí)現(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ù)實(shí)現(xiàn)單張圖像的導(dǎo)入，通過(guò)imshow 函數(shù)實(shí)現(xiàn)數(shù)字圖像在系統(tǒng)中的實(shí)時(shí)顯示；若導(dǎo)入的為圖像文件夾，則先使用detectSURFF eatures 函數(shù)檢測(cè)和提取每張圖像上的Surf 特征點(diǎn)，根據(jù)特征點(diǎn)匹配實(shí)現(xiàn)圖像拼接，采用MSAC 算法對(duì)拼接之后的圖像進(jìn)行幾何變換，生成最終的全景圖像用于后續(xù)處理。

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

1.2 圖像處理部分

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

1.3 圖像分類部分

經(jīng)上述步驟完成各類地物的分割與提取后，通過(guò)對(duì)保存的標(biāo)記數(shù)據(jù)進(jìn)行疊加，實(shí)現(xiàn)各類地物的合并，并在原圖上通過(guò)顏色區(qū)分展示具體的分類效果。此外，對(duì)每類地物進(jìn)行標(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è)置

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

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

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

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

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

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

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

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

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

2.1.2 圖像處理部分

1）K-means 聚類

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

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

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

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

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

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

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

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

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

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

進(jìn)行Otsu 閾值分割之后得到的二值圖像往往包含背景噪聲［圖3（a）］。中值濾波是一種非線性操作，使用medfilt2 函數(shù)對(duì)二值圖像進(jìn)行中值濾波，可以濾除圖中的椒鹽噪聲，同時(shí)保留邊緣；使用imfill 函數(shù)對(duì)圖中的密集獨(dú)立的的孔洞區(qū)域進(jìn)行填充使其連成一整塊區(qū)域，實(shí)現(xiàn)二值圖像的優(yōu)化處理［圖3（b）］。再使用bwlabel 函數(shù)對(duì)分割后的二值圖像進(jìn)行處理，生成“植被”區(qū)域的標(biāo)簽數(shù)據(jù)，實(shí)現(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ù)尋找二值圖像的輪廓邊界線及對(duì)應(yīng)的坐標(biāo)，并用shapewrite 函數(shù)將對(duì)應(yīng)的坐標(biāo)信息和“Line”元素寫(xiě)入形文件，對(duì)結(jié)果進(jìn)行矢量化輸出，便于導(dǎo)入其它軟件進(jìn)行處理，實(shí)現(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 步驟，對(duì)初分割得到的其它區(qū)域進(jìn)行處理，完成圖中其它地物的提取，得到對(duì)應(yīng)地物的標(biāo)簽數(shù)據(jù)。選擇“多標(biāo)記疊加”按鈕會(huì)出現(xiàn)一個(gè)新的操作子界面（圖5），可對(duì)此前處理得到的各類地物標(biāo)記數(shù)據(jù)進(jìn)行合并與類別標(biāo)注。

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

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

圖5 系統(tǒng)子界面：（a）地物標(biāo)記疊加，（b）面積計(jì)算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 開(kāi)發(fā)的無(wú)人機(jī)遙感圖像分類系統(tǒng)，包含函數(shù)功能演示模塊和圖像系統(tǒng)處理模塊兩部分。函數(shù)功能演示模塊為輔助模塊；圖像系統(tǒng)處理模塊為核心模塊，可對(duì)圖像進(jìn)行拼接、分區(qū)、分割、去噪、形態(tài)學(xué)優(yōu)化等操作，最終實(shí)現(xiàn)無(wú)人機(jī)遙感圖像分類。

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

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