劉智麗，張啟斌，岳德鵬，郝玉光，蘇 凱

(1.北京林業(yè)大學(xué)林學(xué)院，北京 100083；2.中國林業(yè)科學(xué)研究院沙漠林業(yè)實驗中心，巴彥淖爾 015200)

0 引言

近年來，中國的城市化水平大幅度提升。經(jīng)國家統(tǒng)計局統(tǒng)計(http://www.stats.gov.cn/tjsj/)，2017年中國的城市化率已達58.52%，相比2010年增長了11.93%。然而快速城市化在提升人類生活水平的同時，也給城市帶來了一系列環(huán)境及社會問題，如生態(tài)惡化、熱島效應(yīng)、環(huán)境污染、交通擁擠等[1]。城市建成區(qū)是一片具有建成環(huán)境的連續(xù)區(qū)域，包括城市用地(道路、建筑及工業(yè)設(shè)施)及其周邊的生態(tài)用地(公園綠地、行道樹、水體等)，作為衡量城市發(fā)展的一個重要指標(biāo)，其范圍及空間分布不僅可以反映城市發(fā)展的規(guī)模及趨勢，如城市擴張水平、人口密度等[2]，還對城市發(fā)展規(guī)劃起著重要作用。因此，對城市建成區(qū)的提取得到越來越多的關(guān)注[3]。

隨著遙感技術(shù)的不斷發(fā)展，城市建成區(qū)提取的自動化程度與精度都在不斷提升[4-5]。夜間燈光數(shù)據(jù)(nighttime light data，NTL)能夠直接反映城市建成區(qū)的夜間亮度，分布連續(xù)[6]，但受城市邊緣亮度影響較大，因而提取準(zhǔn)確度不高。因此，許多研究將NTL與光學(xué)遙感影像結(jié)合起來提取城市建成區(qū)，例如Zhang等[7]將NTL與歸一化植被指數(shù)(normalized difference vegetation index，NDVI)相結(jié)合，得出基于植被校正的城市夜間燈光指數(shù)(vegetation adjusted NTL urban index，VANUI)，有效豐富了城市邊緣信息[8]，但VANUI會放大水體的反射率，混淆城市邊緣地帶的水體與城市建成區(qū)，且空間分辨率不夠高，提取結(jié)果仍無法滿足更高精度的研究需求。

歐洲航天局發(fā)射的Sentinel-2A光學(xué)遙感衛(wèi)星覆蓋了可見光、近紅外以及短波紅外范圍內(nèi)的13個波段，空間分辨率最高可達10 m[9]，能夠準(zhǔn)確識別城市建成區(qū)的邊界，可為城市邊界提取提供良好的數(shù)據(jù)支持。但目前將Sentinel-2A數(shù)據(jù)應(yīng)用于城市建成區(qū)識別、提取的研究較少，因此，本文對VANUI指數(shù)加以改進，基于Sentinel-2A與NTL數(shù)據(jù)提出基于建筑校正的城市夜間燈光指數(shù)(building adjusted NTL urban index，BANUI)，以包頭市南部為研究區(qū)，進行城市建成區(qū)的提取。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

包頭市位于內(nèi)蒙古自治區(qū)中西部，總面積約為27 768 km2，平均海拔為2 000 m，中部地勢較高，由中部向南北兩側(cè)地勢逐漸變低，橫貫其中部的陰山山脈將包頭市分成中部山區(qū)、北部高原草場與南部平原地區(qū)3個部分，其中南部地區(qū)包括包頭市中心城區(qū)與土默特右旗。近年來，包頭市經(jīng)濟發(fā)展迅速，南部平原地區(qū)城市規(guī)模不斷擴張，在2000—2016年間，包頭市城區(qū)面積增長了71.10%。但當(dāng)?shù)貧夂蚋稍?、地形起伏較大，生態(tài)環(huán)境脆弱，而城市的不斷擴張又加劇了對生態(tài)環(huán)境的破壞，由此引發(fā)了一系列環(huán)境問題。包頭市南部地區(qū)聚集了全市89.69%的人口，城市分布密集，具備進行城市建成區(qū)遙感監(jiān)測的基礎(chǔ)，而中部和北部地區(qū)大多為高原山地和草原，人口稀疏、城市數(shù)量屈指可數(shù)，缺乏城市遙感監(jiān)測的條件，因此將包頭市南部平原地區(qū)作為本次研究的研究區(qū)(圖1)。

圖1 研究區(qū)遙感影像Fig.1 Image of study area

1.2 數(shù)據(jù)源

本文所使用的數(shù)據(jù)有Sentinel-2A遙感影像數(shù)據(jù)、NPP-VIIRS NTL數(shù)據(jù)以及研究區(qū)2016年土地利用數(shù)據(jù)。遙感影像選用處理后的2016年6月份研究區(qū)Sentinel-2A Level-2A數(shù)據(jù)，來源于歐洲航天局網(wǎng)站(https://scihub.copernicus.eu/)，經(jīng)Sen2Cor模型處理，包含12個波段。

目前使用較多的NTL數(shù)據(jù)有DMSP-OLS和NPP-VIIRS 2種，其中，NPP-VIIRS數(shù)據(jù)的空間、時間分辨率都比較高，并已經(jīng)過輻射定標(biāo)，更適宜應(yīng)用于分析與研究[12]。因此，本文選用2016年6月份包頭NPP-VIIRS數(shù)據(jù)，下載自美國國家海洋大氣局網(wǎng)站(https://www.ngdc.noaa.gov/eog/viirs/download_dnb_composites.html)。數(shù)據(jù)預(yù)處理時將亮度值小于0的異常像元賦值為0。

研究區(qū)2016年土地利用調(diào)查數(shù)據(jù)來源于包頭市環(huán)境保護局，但其中缺少城市建成區(qū)的地面實況數(shù)據(jù)，因此基于此數(shù)據(jù)通過手動標(biāo)注來獲取城市建成區(qū)的范圍，作為本次研究的驗證數(shù)據(jù)。

本文中使用的所有地理空間數(shù)據(jù)投影均為WGS1984 UTM投影。

2 研究方法

2.1 技術(shù)路線

本文提出BANUI指數(shù)用于研究區(qū)城市建成區(qū)的提取。首先，基于Sentinel-2A數(shù)據(jù)得到研究區(qū)的歸一化建筑指數(shù)(normalized difference build-up index，NDBI)，將其與NPP-VIIRS數(shù)據(jù)結(jié)合得到BANUI數(shù)據(jù)；然后，基于BANUI利用分水嶺分割算法提取包頭市城市建成區(qū)；最后，對提取結(jié)果進行對比分析和精度驗證。技術(shù)路線如圖2所示。

圖2 技術(shù)路線Fig.2 Technology roadmap

2.2 VANUI

為了將城市與農(nóng)村更好地區(qū)分開來，增大兩者亮度值的差異，Zhang等[14]根據(jù)DMSP/OLS數(shù)據(jù)與NDVI提出VANUI指數(shù)，計算公式為

VANUI=(1-NDVI)NTL，

(1)

式中NDVI為歸一化植被指數(shù)。當(dāng)VANUI接近0時，區(qū)域植被覆蓋較好；VANUI的值較大時，植被相對稀疏。但是計算NDVI時使用的MODIS數(shù)據(jù)空間分辨率僅為500 m，NTL數(shù)據(jù)中的細節(jié)信息無法加以利用；且在VANUI中，非植被覆蓋地區(qū)的亮度值會被過度放大，水體反射與其他異常亮度點也會對提取精度有所影響。

2.3 BANUI

2.3.1 NDBI的提取

由TM影像中提取的NDBI可以準(zhǔn)確描述城市建筑用地信息，計算公式為

(2)

式中TM4和TM5分別為TM影像的近紅外和短波紅外波段的反射率。一般情況下，NDBI=0時，為耕地或林地；NDBI<0時，為水體；NDBI>0時，為建筑用地。

處理后的Sentinel-2A數(shù)據(jù)包含12個波段，不同波段的組合會導(dǎo)致城市建成區(qū)范圍有所差異，因此分析不同地類(建筑、植被、水體、未利用地)上的光譜特征并確定最佳波段組合具有重要意義。本文利用感興趣區(qū)(region of interest，ROI)工具選取研究區(qū)各地類的代表性像元，并利用ArcGIS10.2軟件進行分區(qū)統(tǒng)計，結(jié)果如圖3所示。

圖3 Sentinel-2A數(shù)據(jù)各波段的平均亮度值Fig.3 Mean brightness of each band in Sentinel-2A

參照式(2)中所使用波段的波長，本文從Sentinel-2A數(shù)據(jù)的第8、第11與第12波段中選出合適的波段進行組合。由圖3可知，在第11和12波段上，建筑用地的像元亮度值較高，水體和植被的亮度值遠小于建筑用地；在第8波段上，植被的亮度值達到極值，水體與建筑的亮度值相對較?。辉诘?1、第8波段或第12、第8波段的亮度差值中，建筑為正數(shù)，植被為負數(shù)，水體為接近于0的負值。利用這3個波段的組合構(gòu)建了2個NDBI，分別記作NDBI11和NDBI12，公式為

(3)

(4)

式中：S8,S11,S12分別為Sentinel-2A數(shù)據(jù)中第8、第11與第12波段的反射率。

為了直觀地對比NDBI11與NDBI12在不同地物上的表達差異，選取涵蓋了建筑、植被、水體與未利用地4種類型的矩形區(qū)域，并截取對應(yīng)的NDBI圖像，結(jié)果如圖4所示。圖4顯示，NDBI11和NDBI12都能夠較好地反映建筑用地的范圍，建筑的亮度值明顯高于圖中其他2種地類，水體的亮度值在圖中并不高，植被的亮度值也明顯低于建筑，且植被覆蓋程度越高的地方亮度值越低。此外，2幅圖植被亮度差異不明顯，NDBI12中建筑的亮度略高，因此NDBI12中建筑與植被更容易區(qū)分；但NDBI12中水體的DN值略高于NDBI11的水體，這在一定程度上會干擾建筑用地的提取。為定量分析2種波段組合所獲取的NDBI的優(yōu)劣，本文計算出各種地類的NDBI概率分布并繪制其正態(tài)分布曲線，結(jié)果如圖5所示。由圖5可知，2幅圖中建筑用地的NDBI值最高，水體次之，植被的NDBI值最低，其值遠小于0。2幅圖像中建筑用地與植被的差別都比較明顯，但在圖5(b)中，建筑與水體有較多重疊區(qū)域，此外建筑用地與未利用地在2幅圖中都難以區(qū)分。

(a)NDBI11(b)NDBI12

圖4 NDBI11和NDBI12的代表性區(qū)域

Fig.4RepresentativeareasofNDBI11andNDBI12

(a)NDBI11(b)NDBI12

圖5 不同地類在NDBI中的概率分布及其正態(tài)分布曲線

Fig.5ProbabilitydistributionandnormaldistributioncurvesofdifferentlandtypesinNDBI

基于上圖中正態(tài)分布曲線以及色譜分離分析理論可以計算建筑與其余各種地類之間的分離度R，目前該指數(shù)被廣泛應(yīng)用于色譜柱中2個相鄰元素間分離程度的計算。一般當(dāng)R<1時，2波峰有部分重疊；R=1時，其分離程度可達98%；R=1.5時，分離程度接近99.7%，此時2個相鄰元素完全分離；R=2時，2波峰從基線起就完全分離。通常情況下，利用色譜峰的保留時間以及基線上的峰寬可計算分離度R，本文則借助各地類正態(tài)曲線峰值所對應(yīng)的NDBI值計算分離度，公式為

(5)

式中：BR1和BR2分別為2條正態(tài)曲線的峰值所在NDBI的值；W1和W2分別為2條正態(tài)曲線在基線上的峰寬。計算結(jié)果如表1所示。

表1 NDBI11與NDBI12中建筑與各地類之間的分離度RTab.1 R between buildings and the other land types in NDBI11 and NDBI12

由表1計算結(jié)果可知，建筑與植被的分離度均大于2，表明2種地類的正態(tài)曲線完全分離，在NDBI中差異明顯，但NDBI12中兩者的分離度更高；建筑與水體在NDBI11中的分離度高于NDBI12中的分離度，因此NDBI11中建筑與水體差別明顯；建筑與未利用地的分離度都比較低、不易區(qū)分。由此可以推斷，除了建筑與未利用地，其他地類信息在NDBI11中都得到了不同程度的抑制，而NDBI12對水體信息的抑制作用并不明顯，NDBI11與NTL數(shù)據(jù)的結(jié)合可以有效消除噪聲的影響、并增強建筑信息的表達。因此本文選擇Sentinel-2A數(shù)據(jù)中第8波段和第11波段來提取NDBI。

2.3.2 BANUI的提取

在VANUI中，水體的值相對偏離正常范圍，且空間分辨率較低。為解決這一問題，本文利用NDBI與NPP-VIIRS數(shù)據(jù)構(gòu)建了BANUI指數(shù)，公式為

BANUI=(1+NDBI)VNTL，

(6)

式中VNTL為NPP-VIIRS數(shù)據(jù)的亮度值。經(jīng)過重采樣后的NPP-VIIRS數(shù)據(jù)的空間分辨率為20 m，最終BANUI的分辨率同樣為20 m。為檢驗BANUI的提取效果，同樣利用Sentinel-2A與NPP-VIIRS數(shù)據(jù)計算得到VANUI，其空間分辨率為10 m。

選取E109.35°～110.85°范圍內(nèi)一橫截面，該橫截面可以最大程度地經(jīng)過包頭市南部城市建成區(qū)，信息豐富，基于此可檢驗BANUI中建筑信息的表達。橫截面上VANUI與BANUI的值如圖6所示。由圖6可知，兩者變化趨勢比較相似；橫截面經(jīng)過包頭市一公園，公園內(nèi)部植被覆蓋度較高，且因NDVI對植被更敏感，導(dǎo)致VANUI的值要比BANUI更低些。

圖6 橫截面上VANUI與BANUI的值Fig.6 Values of VANUI and BANUI on the cross section

表2數(shù)據(jù)為研究區(qū)內(nèi)VANUI和BANUI的像元亮度值范圍、像元亮度均值、像元亮度值標(biāo)準(zhǔn)差，以及水體、植被、建筑的像元亮度均值。由于未利用地距離城區(qū)較遠，其面積只有153.47 km2(研究區(qū)總面積為4 965.78 km2)，對于城市建成區(qū)的提取影響甚微，因此未利用地的數(shù)據(jù)并未包含在其中。

表2 VANUI和BANUI的統(tǒng)計數(shù)據(jù)Tab.2 Statistical data of VANUI and BANUI

從表2可以看出，BANUI像元亮度值覆蓋范圍較廣、標(biāo)準(zhǔn)差更高，其變化幅度更大，所涵蓋的信息更多；水體在BANUI中的像元亮度均值較低，建筑在BANUI中的像元亮度均值較高，表明BANUI更好地凸顯了建筑信息。

圖7為包頭市城市建成區(qū)邊緣地帶，包含了水體、植被以及連續(xù)分布建筑等信息，因此將其作為檢驗BANUI區(qū)分建筑用地效果的典型區(qū)域。

(a)典型區(qū)域位置 (b)典型區(qū)域影像

圖7 典型區(qū)域

Fig.7Typicalarea

該區(qū)域水體、建筑、植被的概率分布以及正態(tài)分布情況如圖8所示。

(a)VANUI (b)BANUI

圖8 典型區(qū)域內(nèi)不同地類在VANUI與BANUI中的概率分布及正態(tài)分布曲線

Fig.8ProbabilitydistributionandnormaldistributioncurvesofdifferentlandtypesofVANUIandBANUIinthetypicalarea

圖8顯示，BANUI中建筑與植被的像元亮度值要高于VANUI中建筑與植被的亮度值，水體的亮度值在兩幅圖像中相差不大，導(dǎo)致圖8(b)中建筑與水體距離更遠，二者更容易區(qū)分。

計算VANUI和BANUI中建筑與其他地類的分離度，結(jié)果顯示，建筑與植被的分離度都超過了1.5(1.76和1.52)；而BANUI中建筑與水體的分離度(0.63)明顯高于VANUI中二者的分離度(0.29)，表明BANUI相比VANUI能夠更好地將建筑與水體區(qū)分開來。由此可見，BANUI可以豐富NPP-VIIRS的信息、更好地抑制水體反射造成的影響，并增強建筑信息的表達，因此將其視作提取城市建成區(qū)的有效數(shù)據(jù)。

2.4 分水嶺分割算法

為檢驗BANUI在城市建成區(qū)提取應(yīng)用中的優(yōu)勢，本文利用分水嶺分割算法從NPP-VIIRS與VANUI數(shù)據(jù)中提取城市建成區(qū)。

分水嶺分割算法是一種基于區(qū)域進行影像分割的有效方法，可用于城市建成區(qū)范圍的提取[15]。BANUI數(shù)據(jù)的像元被視作地形表面，每個像元的亮度值為地形表面上相應(yīng)像元的高程值，亮度值高的像元為高地，亮度值低的像元為低地，城市建成區(qū)的邊界為地形中集水盆地的邊界即分水嶺[16]。

主要步驟如下：①首先，對影像進行平滑處理，去除BANUI影像中椒鹽噪聲的影響，自適應(yīng)加權(quán)中值濾波方法可在去除噪聲的同時保留更多信息，而高斯濾波、中值濾波等常用方法則無法滿足這一需求；②利用形態(tài)重構(gòu)提高影像中同類物體的同質(zhì)性，并降低不同物體之間的異質(zhì)性，由此得到形態(tài)重構(gòu)影像；③計算影像的形態(tài)梯度，通常情況下像元亮度值變化程度較高的地區(qū)，影像梯度較大；④借助形態(tài)重構(gòu)影像求得局部極大值并標(biāo)記為前景對象，通過歐氏距離倒數(shù)法標(biāo)記背景影像；⑤利用分水嶺分割算法劃分標(biāo)有前景與背景對象的待處理影像，提取出城市建成區(qū)。

2.5 對比分析與精度驗證

為了檢驗本文提取城市建成區(qū)所采取方法的有效性，選取包頭市2016年的土地利用調(diào)查數(shù)據(jù)作為對照，通過人機交互方式在其基礎(chǔ)上繪制出城市建成區(qū)的范圍作為驗證數(shù)據(jù)，此范圍具有高精度、連續(xù)性高等特點。

通過實驗結(jié)果與驗證數(shù)據(jù)之間逐像元的對比建立混淆矩陣，進而計算總體精度、Kappa系數(shù)、生產(chǎn)者精度和用戶精度等指標(biāo)進行精度檢驗[17]。此外，需借助景觀指數(shù)來描述結(jié)果空間結(jié)構(gòu)上的特征，包括總面積(total area，TA)、景觀形狀指數(shù)(landscape shape index，LSI)和聚集指數(shù)(aggregation index，AI)等。

LSI可衡量整個景觀的形狀與面積相同的正方形或圓形之間的偏差程度，以此來表示該景觀整體形態(tài)的復(fù)雜程度[18]，取值范圍為[1，+∞)。LSI值越大，該景觀形狀與正方形的偏差越大。公式為

(7)

式中：0.25為正方形校正常數(shù)；E為景觀中所有斑塊的邊緣總長度，100 m；A為景觀總面積，hm2。

AI表示景觀內(nèi)部斑塊之間相臨的格網(wǎng)單元數(shù)目與在這些數(shù)中的最大值的比值，它反映了景觀的聚集程度與連續(xù)性[19]。公式為

(8)

式中g(shù)i為景觀內(nèi)斑塊i相鄰的格網(wǎng)單元數(shù)。AI的取值范圍為(0,100%]，當(dāng)景觀內(nèi)部斑塊無限分散時，AI的值趨向于0；隨著斑塊的逐漸集聚，AI的值也隨之增加。

3 成果分析與檢驗

3.1 提取結(jié)果

利用分水嶺分割算法由不同數(shù)據(jù)中提取的城市建成區(qū)(圖9)。由圖9可知，利從BANUI中提取的城市建成區(qū)與驗證數(shù)據(jù)相似度較高，而從VANUI和NPP-VIIRS中提取的城市建成區(qū)相比于驗證數(shù)據(jù)范圍較大。由BANUI中提取的城市建成區(qū)包括城市用地，以及水體、公園等其他地類，空間分布連續(xù)、成分完整，與上文對城市建成區(qū)的定義相吻合。

(a)NTL (b)VANUI (c)BANUI

圖9 由不同數(shù)據(jù)提取的城市建成區(qū)

Fig.9Extractedurbanbuilt-upareasfromdifferentdata

3.2 對比分析與精度驗證

通過上文所述精度驗證方法得到不同數(shù)據(jù)的提取精度，如圖10所示。從圖10可以看出，由BANUI提取出的城市建成區(qū)總體精度最高，為93.61%，Kappa系數(shù)為0.793 4，城市建成區(qū)的生產(chǎn)者精度與用戶精度分別為85.34%和81.34%，總體上提取的精度較高；VANUI位列其次，NPP-VIIRS精度最低。

圖10 不同數(shù)據(jù)源提取的城市建成區(qū)的精度對比Fig.10 Accuracy comparison of urban built-up areas extracted by different data sources

表3為不同數(shù)據(jù)源提取的城市建成區(qū)景觀格局指數(shù)對比，其中，TA.Diff，LSI.Diff，AI.Diff分別表示TA，LSI，AI與驗證數(shù)據(jù)的誤差。由表4可知，BANUI提取結(jié)果的景觀指數(shù)與驗證數(shù)據(jù)十分接近，提取的城市建成區(qū)總面積比驗證數(shù)據(jù)多了4.91%，LSI和AI的差距也較小，提取效果最佳；而其余2個數(shù)據(jù)源所提取結(jié)果與驗證數(shù)據(jù)相差較大，面積分別多出了24.12%和16.55%，精度較低。

表3 不同數(shù)據(jù)源提取的城市建成區(qū)景觀格局指數(shù)對比Tab.3 Comparison of landscape pattern indexes of urban built-up areas extracted by different data sources

4 結(jié)論

1)基于Sentinel-2A與NPP-VIIRS數(shù)據(jù)，對傳統(tǒng)的VANUI進行改進，提出了用于提取城市建成區(qū)的BANUI指數(shù)。

2)通過對Sentinel-2A數(shù)據(jù)各個波段的特征進行分析，選擇Sentinel-2A影像的第8和第11波段計算研究區(qū)NDBI，從而構(gòu)建BANUI指數(shù)

3)實驗結(jié)果表明，由BANUI所提取的研究區(qū)城市建成區(qū)的面積為971.25 km2，相比驗證數(shù)據(jù)僅相差4.91%，總體精度高達93.61%，Kappa系數(shù)為0.793 4，明顯優(yōu)于由VANUI和NPP-VIIRS數(shù)據(jù)的提取結(jié)果。

4)BANUI在城市建成區(qū)的提取上具有結(jié)果完整、覆蓋范圍廣、提取精度高等優(yōu)點，應(yīng)用前景廣闊，可為進一步研究NTL數(shù)據(jù)對城市建成區(qū)的提取提供思路，也可用于城市發(fā)展的規(guī)劃與監(jiān)測。