基于Sentinel-2A與NPP-VIIRS夜間燈光數(shù)據(jù)的城市建成區(qū)提取

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

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

0 引言

近年來，中國的城市化水平大幅度提升。經(jīng)國家統(tǒng)計(jì)局統(tǒng)計(jì)(http://www.stats.gov.cn/tjsj/)，2017年中國的城市化率已達(dá)58.52%，相比2010年增長了11.93%。然而快速城市化在提升人類生活水平的同時(shí)，也給城市帶來了一系列環(huán)境及社會問題，如生態(tài)惡化、熱島效應(yīng)、環(huán)境污染、交通擁擠等[1]。城市建成區(qū)是一片具有建成環(huán)境的連續(xù)區(qū)域，包括城市用地(道路、建筑及工業(yè)設(shè)施)及其周邊的生態(tài)用地(公園綠地、行道樹、水體等)，作為衡量城市發(fā)展的一個(gè)重要指標(biāo)，其范圍及空間分布不僅可以反映城市發(fā)展的規(guī)模及趨勢，如城市擴(kuò)張水平、人口密度等[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個(gè)波段，空間分辨率最高可達(dá)10 m[9]，能夠準(zhǔn)確識別城市建成區(qū)的邊界，可為城市邊界提取提供良好的數(shù)據(jù)支持。但目前將Sentinel-2A數(shù)據(jù)應(yīng)用于城市建成區(qū)識別、提取的研究較少，因此，本文對VANUI指數(shù)加以改進(jìn)，基于Sentinel-2A與NTL數(shù)據(jù)提出基于建筑校正的城市夜間燈光指數(shù)(building adjusted NTL urban index，BANUI)，以包頭市南部為研究區(qū)，進(jìn)行城市建成區(qū)的提取。

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

1.1 研究區(qū)概況

包頭市位于內(nèi)蒙古自治區(qū)中西部，總面積約為27 768 km2，平均海拔為2 000 m，中部地勢較高，由中部向南北兩側(cè)地勢逐漸變低，橫貫其中部的陰山山脈將包頭市分成中部山區(qū)、北部高原草場與南部平原地區(qū)3個(gè)部分，其中南部地區(qū)包括包頭市中心城區(qū)與土默特右旗。近年來，包頭市經(jīng)濟(jì)發(fā)展迅速，南部平原地區(qū)城市規(guī)模不斷擴(kuò)張，在2000—2016年間，包頭市城區(qū)面積增長了71.10%。但當(dāng)?shù)貧夂蚋稍?、地形起伏較大，生態(tài)環(huán)境脆弱，而城市的不斷擴(kuò)張又加劇了對生態(tài)環(huán)境的破壞，由此引發(fā)了一系列環(huán)境問題。包頭市南部地區(qū)聚集了全市89.69%的人口，城市分布密集，具備進(jìn)行城市建成區(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個(gè)波段。

目前使用較多的NTL數(shù)據(jù)有DMSP-OLS和NPP-VIIRS 2種，其中，NPP-VIIRS數(shù)據(jù)的空間、時(shí)間分辨率都比較高，并已經(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ù)處理時(shí)將亮度值小于0的異常像元賦值為0。

研究區(qū)2016年土地利用調(diào)查數(shù)據(jù)來源于包頭市環(huán)境保護(hù)局，但其中缺少城市建成區(qū)的地面實(shí)況數(shù)據(jù)，因此基于此數(shù)據(jù)通過手動標(biāo)注來獲取城市建成區(qū)的范圍，作為本次研究的驗(yàn)證數(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é)果進(jìn)行對比分析和精度驗(yàn)證。技術(shù)路線如圖2所示。

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

2.2 VANUI

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

VANUI=(1-NDVI)NTL，

(1)

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

2.3 BANUI

2.3.1 NDBI的提取

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

(2)

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

處理后的Sentinel-2A數(shù)據(jù)包含12個(gè)波段，不同波段的組合會導(dǎo)致城市建成區(qū)范圍有所差異，因此分析不同地類(建筑、植被、水體、未利用地)上的光譜特征并確定最佳波段組合具有重要意義。本文利用感興趣區(qū)(region of interest，ROI)工具選取研究區(qū)各地類的代表性像元，并利用ArcGIS10.2軟件進(jìn)行分區(qū)統(tǒng)計(jì)，結(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波段中選出合適的波段進(jìn)行組合。由圖3可知，在第11和12波段上，建筑用地的像元亮度值較高，水體和植被的亮度值遠(yuǎn)小于建筑用地；在第8波段上，植被的亮度值達(dá)到極值，水體與建筑的亮度值相對較小；在第11、第8波段或第12、第8波段的亮度差值中，建筑為正數(shù)，植被為負(fù)數(shù)，水體為接近于0的負(fù)值。利用這3個(gè)波段的組合構(gòu)建了2個(gè)NDBI，分別記作NDBI11和NDBI12，公式為

(3)

(4)

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

為了直觀地對比NDBI11與NDBI12在不同地物上的表達(dá)差異，選取涵蓋了建筑、植被、水體與未利用地4種類型的矩形區(qū)域，并截取對應(yīng)的NDBI圖像，結(jié)果如圖4所示。圖4顯示，NDBI11和NDBI12都能夠較好地反映建筑用地的范圍，建筑的亮度值明顯高于圖中其他2種地類，水體的亮度值在圖中并不高，植被的亮度值也明顯低于建筑，且植被覆蓋程度越高的地方亮度值越低。此外，2幅圖植被亮度差異不明顯，NDBI12中建筑的亮度略高，因此NDBI12中建筑與植被更容易區(qū)分；但NDBI12中水體的DN值略高于NDBI11的水體，這在一定程度上會干擾建筑用地的提取。為定量分析2種波段組合所獲取的NDBI的優(yōu)劣，本文計(jì)算出各種地類的NDBI概率分布并繪制其正態(tài)分布曲線，結(jié)果如圖5所示。由圖5可知，2幅圖中建筑用地的NDBI值最高，水體次之，植被的NDBI值最低，其值遠(yuǎn)小于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)分布曲線以及色譜分離分析理論可以計(jì)算建筑與其余各種地類之間的分離度R，目前該指數(shù)被廣泛應(yīng)用于色譜柱中2個(gè)相鄰元素間分離程度的計(jì)算。一般當(dāng)R<1時(shí)，2波峰有部分重疊；R=1時(shí)，其分離程度可達(dá)98%；R=1.5時(shí)，分離程度接近99.7%，此時(shí)2個(gè)相鄰元素完全分離；R=2時(shí)，2波峰從基線起就完全分離。通常情況下，利用色譜峰的保留時(shí)間以及基線上的峰寬可計(jì)算分離度R，本文則借助各地類正態(tài)曲線峰值所對應(yīng)的NDBI值計(jì)算分離度，公式為

(5)

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

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

由表1計(jì)算結(jié)果可知，建筑與植被的分離度均大于2，表明2種地類的正態(tài)曲線完全分離，在NDBI中差異明顯，但NDBI12中兩者的分離度更高；建筑與水體在NDBI11中的分離度高于NDBI12中的分離度，因此NDBI11中建筑與水體差別明顯；建筑與未利用地的分離度都比較低、不易區(qū)分。由此可以推斷，除了建筑與未利用地，其他地類信息在NDBI11中都得到了不同程度的抑制，而NDBI12對水體信息的抑制作用并不明顯，NDBI11與NTL數(shù)據(jù)的結(jié)合可以有效消除噪聲的影響、并增強(qiáng)建筑信息的表達(dá)。因此本文選擇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。為檢驗(yàn)BANUI的提取效果，同樣利用Sentinel-2A與NPP-VIIRS數(shù)據(jù)計(jì)算得到VANUI，其空間分辨率為10 m。

選取E109.35°～110.85°范圍內(nèi)一橫截面，該橫截面可以最大程度地經(jīng)過包頭市南部城市建成區(qū)，信息豐富，基于此可檢驗(yàn)BANUI中建筑信息的表達(dá)。橫截面上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ū)較遠(yuǎn)，其面積只有153.47 km2(研究區(qū)總面積為4 965.78 km2)，對于城市建成區(qū)的提取影響甚微，因此未利用地的數(shù)據(jù)并未包含在其中。

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

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

圖7為包頭市城市建成區(qū)邊緣地帶，包含了水體、植被以及連續(xù)分布建筑等信息，因此將其作為檢驗(yàn)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)中建筑與水體距離更遠(yuǎn)，二者更容易區(qū)分。

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

2.4 分水嶺分割算法

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

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

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

2.5 對比分析與精度驗(yàn)證

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

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

LSI可衡量整個(gè)景觀的形狀與面積相同的正方形或圓形之間的偏差程度，以此來表示該景觀整體形態(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)部斑塊無限分散時(shí)，AI的值趨向于0；隨著斑塊的逐漸集聚，AI的值也隨之增加。

3 成果分析與檢驗(yàn)

3.1 提取結(jié)果

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

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

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

Fig.9Extractedurbanbuilt-upareasfromdifferentdata

3.2 對比分析與精度驗(yàn)證

通過上文所述精度驗(yàn)證方法得到不同數(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與驗(yàn)證數(shù)據(jù)的誤差。由表4可知，BANUI提取結(jié)果的景觀指數(shù)與驗(yàn)證數(shù)據(jù)十分接近，提取的城市建成區(qū)總面積比驗(yàn)證數(shù)據(jù)多了4.91%，LSI和AI的差距也較小，提取效果最佳；而其余2個(gè)數(shù)據(jù)源所提取結(jié)果與驗(yàn)證數(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進(jìn)行改進(jìn)，提出了用于提取城市建成區(qū)的BANUI指數(shù)。

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

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

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