吳亞玲 黃義忠 彭秋志

摘要：為了提高從海量遙感影像中解譯建設用地的效率，設計一種基于Python從Landsat8 OLI遙感影像中自動識別與提取建設用地的軟件系統(tǒng)。利用時域中值濾波去除云噪聲，并通過z值標準化獲得可適用統(tǒng)一閾值劃分的穩(wěn)定像元；結(jié)合多種地物指數(shù)開展閾值劃分以初步識別可能的建設用地像元；結(jié)合像元鄰域組合關系進一步精化提取結(jié)果，并自動計算對應的Kappa系數(shù)。研究表明，與傳統(tǒng)基于大型專業(yè)遙感影像解譯平臺相比，該系統(tǒng)不僅能保持較高的解譯精度，而且具有體積小、速度快、自動化程度高、可遷移性和可擴展性強等特點。

關鍵詞：Python；遙感影像；自動識別；建設用地

DOIDOI：10.11907/rjdk.181322

中圖分類號：TP319

文獻標識碼：A 文章編號：1672-7800（2018）010-0161-04

英文摘要Abstract：In order to improve the efficiency of interpreting construction land from massive remote sensing images， the software system based on Python which can automatically identify and extract construction land from Landsat 8 OLI remote sensing images was designed. Firstly， median filtration was used to remove cloud interference， the stable pixels of uniform threshold value partitioning that was obtained through normalization of z value. Then，the threshold partitioning combined with various features indices are employed to identify possible pixels of construction land. Finally， the meta-neighborhood combination relationship is further refined to extract the result， and the corresponding Kappa coefficient is automatically calculated. The results show that the proposed system can not only maintain high accuracy of interpretation， but also has the advantages of small size， high speed and high degree of automation compared with the traditional method of recognition and extraction of construction land based on large professional remote sensing images interpretation system.

英文關鍵詞Key Words：Python；remote sensing image；automatic identification；construction land

0 引言

社會經(jīng)濟的持續(xù)發(fā)展以及城市化進程中伴隨的城市人口、社會發(fā)展要素重組等，對城市建設用地利用方式產(chǎn)生了巨大影響[1-2]。城市建設用地動態(tài)監(jiān)測空間基礎數(shù)據(jù)中遙感影像分辨率的提高、分類方法的改進，對實現(xiàn)城市建設用地的高效利用具有重要意義[3]。

隨著遙感影像大范圍普及，高效、快捷、自動化地獲取建設用地信息已經(jīng)成為一種趨勢[3-4]。近年來，為了獲得更高的分類精度，眾多從遙感影像中提取建設用地信息的方法相繼被提出。這些方法就分類單元而言可分為基于像元的分類和基于對象的分類[5]，就分類途徑而言可分為監(jiān)督分類和非監(jiān)督分類[6]，就分類結(jié)果表現(xiàn)形式而言可分為確定性分類和模糊分類[7-8]，就分類所用信息類別而言可分為僅基于光譜信息的分類和以光譜信息為主的多信息組合分類等[9-10]。然而，日益多樣的方法雖然帶來更多選擇可能，但同時也增加了學習成本、增大了選擇困難。在實際工程應用中，不僅需要考慮解譯精度，而且需要兼顧簡單易學、便捷高效、成本低廉等要求。因此，有必要針對特定的專題解譯需求，開發(fā)出一系列自動化軟件系統(tǒng)，以減輕對大型專業(yè)遙感影像處理平臺的依賴，大幅降低成本投入。

本文借助面向?qū)ο蟮慕忉屝陀嬎銠C程序設計語言——Python，綜合運用多波段指數(shù)法[11]和歸一化閾值控制法[12 -14]，提出一種節(jié)約時間、更新速度快的建設用地識別與提取技術，以期自動化和批量化從海量遙感影像中解譯建設用地。

1 系統(tǒng)原理與流程設計

建設用地自動識別與提取系統(tǒng)流程如圖1所示?；驹硎牵米詣踊募闅v和迭代方法串聯(lián)起每一個解譯步驟，實現(xiàn)“進去的是原始數(shù)據(jù)，出來的是最終產(chǎn)品”之目的。該系統(tǒng)主要包括6個關鍵解譯步驟：①對下載的原始遙感影像進行數(shù)據(jù)預處理，以便實現(xiàn)數(shù)據(jù)輸入自動化，主要模塊包括解壓縮文件、重命名文件、提取關鍵波段文件和獲取影像最小重疊區(qū)；②疊合同一時期多期影像，利用時域中值濾波法實現(xiàn)去云降噪；③通過z值標準化獲得可適用統(tǒng)一閾值劃分的穩(wěn)定像元；④結(jié)合多種已被證明效果良好且穩(wěn)定的地物指數(shù)開展閾值劃分，以初步識別可能的建設用地像元；⑤結(jié)合像元鄰域組合關系進一步精化提取結(jié)果；⑥自動計算Kappa系數(shù)，并根據(jù)該系數(shù)對前兩步進行迭代式調(diào)整，保證解譯結(jié)果的質(zhì)量。

2 系統(tǒng)實現(xiàn)

系統(tǒng)實現(xiàn)基于Python語言進行編碼，數(shù)據(jù)預處理和質(zhì)量控制部分通過自行編碼實現(xiàn)，主要空間分析函數(shù)調(diào)用自ArcGIS中的空間分析模塊。

2.1 數(shù)據(jù)預處理

（1）對下載的原始數(shù)據(jù)進行批量解壓。對海量原始影像壓縮文件夾進行遍歷，識別并逐次解壓。文件解壓關鍵代碼如下：

targzfile.extract（file， unzipedfolder） #從壓縮文件targzfile中把解壓出的文件file放到文件夾unzipedfolder里

（2）重命名。在系統(tǒng)中，為了后續(xù)流程的完整性，對已解壓文件名不統(tǒng)一的波段重新命名。文件重命名關鍵代碼如下：

arcpy.Rename_management（name_fr， name_to） #從將文件name_fr重命名為name_to

（3）提取關鍵波段。海量影像的所有波段中，依據(jù)建設用地識別方法，對波段進行篩選，單獨提取出關鍵的波段數(shù)，并進行另存。

（4）獲取最小重疊區(qū)。去除空值和影像邊界，其留下的最小邊界作為后面研究范圍。以關鍵波段中某一波段為代表，利用Python腳本對其影像進行批量化掩摸，計算最小重疊區(qū)。關鍵代碼如下：

arcpy.gp.CellStatistics_sa（in， out， "MINIMUM"； "DATA"） #獲取所有輸入圖層的最小值

剔除外圍以0值表示的空值黑邊，只保留有值區(qū)，裁剪出最小重疊區(qū)。關鍵代碼如下：

arcpy.gp.Con_sa（in ， 1， out， ""， "Value > 0"） #剔除外圍黑邊

2.2 云噪聲去除

為了消除時相干擾，去除云噪聲，需進行時域中值濾波。其原理是云的反射率在各波段都位于異常極大值區(qū)，而云影大多位于異常極小值區(qū)，取中值具有很好的去云和去云影效果。關鍵代碼如下：

arcpy.gp.CellStatistics_sa（in， out， "MEDIAN"， "DATA"） #獲取所有輸入圖層的中值

2.3 穩(wěn)定像元獲取

不同影像中波段值存在時相差異，需對生成的波段中值進行z值標準化，并計算影像的z值標準差，再以0.5倍標準差為閾值，替換不穩(wěn)定像元，最后反z值標準化，即基于統(tǒng)一的平均值和標準差系數(shù)反算回去，獲取每個單波段的穩(wěn)定像元圖層。該步驟能保證在后續(xù)步驟中影像均采用相同閾值進行分割。關鍵代碼如下：

MEAN = arcpy.GetRasterProperties_management（in， "MEAN"， "Band_1"） #計算多期影像平均值

STD = arcpy.GetRasterProperties_management（in， "STD"， "Band_1"） #計算多期影像標準差

arcpy.gp.Minus_sa（in， MEAN， numerator） #獲得每期影像與多期影像平均值之差，以作為標準化的分子numerator

arcpy.gp.Divide_sa（numerator， STD， out） #獲得每期影像與numerator之比，得到該期影像的標準化值

arcpy.gp.Con_sa（in ，1， out， ""， "Value < 0.5"） #獲得均值兩側(cè)0.5倍標準差范圍內(nèi)的像元，即穩(wěn)定像元

2.4 建設用地初步識別

建設用地提取過程實質(zhì)是遙感影像信息的處理過程[15] ，遙感技術手段具有客觀、宏觀、快速、準確、動態(tài)等特點[16]。不同類型的地物反射光譜強度在不同波段上具有不一致性，分析和找出建設用地與其它背景地物間的差異，充分利用遙感影像光譜信息、時相信息、地形信息和鄰域信息，逐步排除確定性無關地物。運用波段差異提取出其特征信息的NDVI（歸一化植被指數(shù)）[17]、MNDWI（歸一化水體指數(shù)）[18]，其值介于（-1，1）之間，經(jīng)過其處理后比值運算用于消除地形差異的影響。

通過地物指數(shù)計算閾值分割方法，在已選用的地物指數(shù)參數(shù)中，構(gòu)建各參數(shù)的閾值范圍，如：NDVI、MNDWI。同時在提取城市建設用地等地物中，僅依據(jù)指數(shù)去除植被和水體是完全不夠的，在系統(tǒng)中根據(jù)每兩個波段間的歸一化差值取閾值去除高亮地物（旱裸土、部分灘涂、部分高亮人工地面、部分高亮屋頂）、水體及陰坡植被、裸地及植被、海島湖泊及灘涂等其它不屬于建設用地的地物類型。其中地物參數(shù)選取和閾值分割均在Python中批量化處理，用到的波段組合原理為：

NDVI=NIR-RedNIR+Red（1）

MNDWI=Green-SWIR1Green+SWIR1（2）

D73=SWIR2×Red（3）

ND75=SWIR2-SWIR1SWIR2-SWIR1（4）

ND42=NIR-GreenNIR+Green（5）

ND75ND42=ND75-ND42（6）

式（1）-式（6）均表示在Python中計算用到的指數(shù)。針對實驗區(qū)現(xiàn)狀，反復試驗后，優(yōu)先構(gòu)建出一套完整的建設用地提取流程，用來提取建設用地信息。其中，部分關鍵代碼如下：

（1）根據(jù)波段計算出NDVI、MNDWI指數(shù)。

（2）采用歸一化閾值法確定植被的閾值范圍為0.25，粗略去除植被。

arcpy.gp.Con_sa（MNDWI ，1， Pout， 0， 'Value > 0.25'） #獲取像元大于0.25的NDVI，用于識別植被

（3）采用步驟（2）中的歸一化閾值法，確定水體的閾值范圍為0.1，去除水體。

arcpy.gp.Con_sa（MNDWI ，1， Pout， 0， 'Value > 0.10'） #獲取像元大于0.1的MNDWI，用于識別水體

（4）在前3步的基礎上，通過光譜和紋理分析，去除試驗區(qū)內(nèi)比較明顯的植被、水體以后，通過光譜差異和借用實驗區(qū)范圍內(nèi)的Google Earth影像對比，進一步識別出需要識別去除的高亮地物（含裸土、部分灘涂、部分高亮人工地面、部分高亮屋頂）、水體及偏陰坡植被、裸地及植被。

arcpy.gp.Con_sa（D73 ， 1， Pout， 0， 'Value > 0.06'） #獲取像元大于0.06的D73，用于識別高亮地物

arcpy.gp.Con_sa（D73 ，1， Pout，0， 'Value < 0.01'） #獲取像元小于0.01的D73，用于識別水體及偏陰坡植被

arcpy.gp.Con_sa（ND75ND42 ，1， Pout， 0， 'Value < -0.40'） #獲取像元大于-0.04的ND75ND42，識別裸地及植被

2.5 建設用地精化提取

初步識別可能的建設用地后，系統(tǒng)用于精化提取建設用地的方法主要是采用空間鄰域關系計算。分別對水體、植被的邊界像元通過鄰域均值計算后，采用條件函數(shù)逐一判別消除其它地物。其關鍵代碼如下：

arcpy.gp.Con_sa（Vm ， 1， Vm01 ， 0， "VALUE > 0.25"） #獲取鄰域內(nèi)均值Vm超過0.25的植被

arcpy.gp.Con_sa（Wm ， 1， Wm01 ， 0， "VALUE > 0.25"） ##獲取鄰域內(nèi)均值Wm超過0.25的水體

arcpy.gp.CellStatistics_sa（Vm01&Wm01;， Vm01sumWm01， "SUM"， "DATA"） #獲取鄰域內(nèi)水體和植被Vm01& Wm01均超過0.25的像元的和

arcpy.gp.Con_sa（Vm01sumWm01 ，1 ， Pout， 0， "VALUE = 2"） #獲取上一步中VALUE = 2的值

2.6 質(zhì)量控制

對遙感分類精度進行評價，需要參考數(shù)據(jù)作為評價基準。系統(tǒng)選擇Google Earth影像作為參考數(shù)據(jù)，在針對完整流程提取建設用地的基礎上，于所選實驗區(qū)內(nèi)隨機生成1 000個參考點，對建設用地提取結(jié)果進行精度評價。采用的評價指標為Kappa系數(shù)，其計算在Python中進行。

3 實例驗證

3.1 實驗區(qū)選擇與數(shù)據(jù)下載

實驗選取昆明中心城區(qū)（五華區(qū)、盤龍區(qū)、西山區(qū)、官渡區(qū)、呈貢新區(qū)和空港經(jīng)濟區(qū)）作為實驗區(qū)，如圖2所示。采用2015年的Landsat 8 OLI遙感影像數(shù)據(jù)，數(shù)據(jù)下載自https：//glovis.usgs.gov/。為保障影像質(zhì)量，避免云量干擾，下載的影像數(shù)據(jù)云量均控制在10%以內(nèi)。對地面同一地點，下載同時區(qū)影像數(shù)量確保大于或等于10，本文實驗區(qū)所在的影像行列號為129043。以上條件限制主要為了在后續(xù)自動提取建設用地過程中，能利用盡量多的時相信息消除云干擾影響，減少季相差異，從而提高分類精度。其Landsat 8 OLI遙感影像數(shù)據(jù)均經(jīng)過Level1標準地形校正，即經(jīng)過系列輻射校正、地面控制點幾何校正和DEM地形校正，并以單波段輸出格式為GeoTIFF的形式封裝在壓縮包文件內(nèi)，取樣方式均為三次卷積，地圖投影均采用UTM-WGS84投影。本文中運用到的波段為：Band3（Green）、Band4（Red）、Band5（NIR）、Band6（SWIR1）、Band7（SWIR2）5個波段，分辨率均為30m。

3.2 建設用地提取結(jié)果

圖2顯示了研究區(qū)Landsat 8 OLI影像的RGB波段合成效果，圖3是利用本文所述方法自動提取的建設用地效果。直觀對比，本文所述方法能夠很好地將建設用地與其它地物區(qū)別開來，特別是能很好地剔除耕地和灘涂；就精度定量評估而言，建設用地提取結(jié)果的Kappa系數(shù)為0.812 5，與大多數(shù)現(xiàn)有方法提取建設用地結(jié)果的精度水平相當，可以滿足大部分相關研究和應用的精度要求。

4 結(jié)語

本文在回顧和總結(jié)遙感影像建設用地提取研究的基礎上，介紹了一種基于Python從Landsat8 OLI遙感影像中自動識別與提取建設用地的技術方法。與其它建設用地提取方法相比，其可以在大范圍遙感影像覆蓋區(qū)域?qū)崟r、高效、準確地自動提取建設用地信息，為快速監(jiān)測和評價建設用地變化提供新的途徑。在后續(xù)研究中，可將本文研究方法延伸到其它地物的獲取、監(jiān)測等領域。

參考文獻：

[1] 袁麗麗.城市化進程中城市土地可持續(xù)利用研究[D].武漢：華中農(nóng)業(yè)大學， 2005.

[2] 馬兵.揚州市城市化進程中城鄉(xiāng)建設用地動態(tài)變化研究[D].南京：南京農(nóng)業(yè)大學， 2009.

[3] 王宏勝，李永樹，吳璽，等.結(jié)合空間分析的面向?qū)ο鬅o人機影像土地利用分類[J]. 測繪工程， 2018（2）：57-61.

[4] 何少林，徐京華，張帥毅.面向?qū)ο蟮亩喑叨葻o人機影像土地利用信息提取[J].國土資源遙感，2013，25（2）：107-112.

[5] 張洋華，劉慧平，楊曉彤.基于像元與對象分類的景觀指數(shù)差異性分析[J].遙感技術與應用，2016，31（1）：119-125.

[6] 趙春霞，錢樂祥.遙感影像監(jiān)督分類與非監(jiān)督分類的比較[J].河南大學學報：自然科學版，2004，34（3）：90-93.

[7] 張卉，張超，趙冬玲.特征權重優(yōu)化高分辨率遙感影像模糊分類研究[J].遙感信息，2010（1）：94-98.

[8] 劉艷芳，蘭澤英，劉洋，等.基于混合熵模型的遙感分類不確定性的多尺度評價方法研究[J].測繪學報，2009，38（1）：82-87.

[9] 陳晨，張友靜.基于多尺度紋理和光譜信息的SVM分類研究[J].測繪科學，2009，34（1）：29-31.

[10] 厲小潤，朱潔爾，王晶，等.組合核支持向量機高光譜圖像分類[J].浙江大學學報：工學版，2013，47（8）：1403-1410.

[11] 樊彥國，李翔宇，張磊，等.基于多波段分析的無閾值自動光譜角制圖分類法[J].地理與地理信息科學，2010， 26（2）：38-41.

[12] VIOVY N， ARINO O， BELWARD A S. The best index slope ex-traction（BISE）：a method for reducing noise in NDVI time series[J].International Journal of Remote Sensing，1992，13：1585-1590.

[13] 徐影，丁一匯，趙宗慈.人類活動引起的我國西北地區(qū)21世紀溫度和降水變化情景分析[J].冰川凍土，2003，25（3）：327-330.

[14] 付忠良.一些新的圖像閾值選取方法[J].計算機應用，2000，20（10）：15-17.

[15] 李愛民.基于遙感影像的城市建成區(qū)擴張與用地規(guī)模研究[D].鄭州：解放軍信息工程大學，2009.

[16] 邢雅娟，劉東升，王鵬新.遙感信息與作物生長模型的耦合應用研究進展[J].地球科學進展，2009，24（4）：444-451.

[17] DEFRIES R， TOWNSHEND J. NDVI-derived land cover classification at global scales[J]. International Journal of Remote Sensing， 1994， 15（17）：3567-3586.

[18] 徐涵秋.福清市城鎮(zhèn)空間擴展規(guī)律及其驅(qū)動機制分析[J].遙感技術與應用， 2002， 17（2）：86-92.

（責任編輯：何 麗）