一种基于集成学习的城市新增建设用地快速提取方法
杜培军, 梁昊, 王欣, 栗云峰     
南京大学地理与海洋科学学院,自然资源部国土卫星遥感应用重点实验室,江苏省地理信息技术重点实验室,江苏  南京  210023
摘要: 城市新增建设用地是热岛效应调控、大气污染和扬尘管控、生态服务变化的重要指示信息。为准确提取新增建设用地,提出一种运用多时相Sentinel遥感数据和集成学习算法的城市新增建设用地快速提取方法。基于多时相Sentinel-1和Sentinel-2遥感影像提取光谱特征、纹理特征和后向散射特征,进行面向对象分割、全域均值滤波和归一化后得到组合特征集,运用随机森林、旋转森林、支持向量机和极限学习机多分类器集成学习进行分类来提取新增建设用地。提取了南京市2017年4月至10月间的新增建设用地并统计了各行政区分布的面积,提取整体精度达0.95,Kappa系数达0.88。相比与基于像素方法,面向对象技术可有效降低“椒盐现象”,提高斑块完整性; 与分类后提取方法比较,直接变化提取方法可减少误差产生环节及误差累积,从而降低系统误差,提高提取精度。
关键词: 遥感影像    集成学习    面向对象图像分析    建设用地    变化检测    
A New Urban Built-up Land Extraction Method Based on Ensemble Learning
DU Pei-jun, LIANG Hao, WANG Xin, LI Yun-feng     
School of Geography and Ocean Science, Key Laboratory for Land Satellite Remote Sensing Applications of Ministry of Natural Resources, Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, Nanjing, Jiangsu 210023, China
Abstract: New built-up land is an important indicator to regulate urban heat island, air pollution and dust, and ecological service change. In order to accurately extract new urban built-up land, an object-oriented method of extracting new built-up land using multi-temporal Sentinel data and ensemble learning is proposed. Based on Sentinel-1 and Sentinel-2 remotely sensed images, spectral, texture and backscatter features are extracted, then object-oriented segmentation, global mean filtering and normalization are executed so that combined feature set is produced. Finally ensemble learning based on random forest, rotation forest, support vector machine and extreme learning machine is used to directly extract new built-up land. New built-up land of Nanjing city from April to October, 2017 is extracted, and the area in each district is calculated. The overall accuracy can up to 0.95 while the Kappa coefficient can up to 0.88. Comparing to pixel-based method, Object-oriented technology can effectively reduce the "pepper and salt" phenomenon and improve the integrity of plaque. Comparing to the method of change detection after classification, this method can reduce error generation and error accumulation, so that systematic error can be reduced and extraction accuracy can be improved.
Key words: Remote sensing image    Ensemble learning    Object-oriented image analysis    Built-up land    Change detection    
0 引言

建设用地指由人类活动形成的非农业生产建设地段,包括房屋、道路和构筑物等,具有变化频繁迅速、显著影响生态环境等特点[1]。新增建设用地为非建设用地转化为建设用地的区域,其准确提取可协助城市规划和管理机关掌握城市发展趋势与发展速度,实现地表覆盖动态监测,为预测和评估潜在的城市热岛效应、大气污染、建设工地扬尘、生态服务变化等提供动态信息。基于生态保护红线、永久基本农田、城镇开发边界3条控制线,综合利用遥感信息快速检测出新增和违章建筑,服务于城市环境监测与调控精准,具有重要的实用价值。其中,高效的遥感变化检测和提取方法是最为关键的技术支撑。

变化检测是通过观察不同时期遥感影像鉴别差异的过程[2-3]。变化检测有多种分类体系[4],根据是否需要先分类,可分为分类后比较法和直接比较检测法两类[5]。分类后比较法通过对分类后不同时相对应单元类属的比较来检测变化。如:对美国陆地卫星4~5号专题制图仪(Thematic Mapper, TM)影像进行面向对象分割和决策树分类,然后进行变化分析[6]; 运用先分类后变化检测的方法分析滑坡变化[7]; 结合像元级分类结果和面向对象分割得到对象级的变化检测结果等[8]。直接比较检测法对不同时相遥感影像中的对应单元直接进行分析来检测变化区域,如:将多分辨率变化影像进行决策级融合得到不同的变化强度[9]; 利用纹理、几何、形状等特征构建信息融合策略的变化检测方法[10]; 利用不同特征的差异影像和集成学习直接提取变化目标[11]; 运用Dempster-Shafer(D-S)理论融合多种变化检测结果得到最终的变化检测结果等[12]。研究表明,没有哪一种变化检测方法最优、可适用于所有的应用,应具体问题具体分析[4, 13]

城市建设用地等非自然生态物具有鲜明的回波反射/散射/偏振等信号特征,光学遥感难以精准提取,且易受云雨遮蔽等因素影响。可借助合成孔径雷达(Synthetic Aperture Radar,SAR)获取地物不同极化方式的回波信息,表征目标不同的物理散射特性,丰富建设用地观测和提取的手段[14]

现基于“哨兵1号”(Sentinel-1)SAR和“哨兵2号”(Sentinel-2)光学遥感数据,运用面向对象技术,在中分辨率尺度上建立了一套运用集成学习算法进行变化检测来提取新增建设用地的方法。

1 研究区和数据 1.1 研究区概况

研究区位置见图 1(a)(b)。江苏省南京市是长三角城市群中仅次于上海的唯一特大城市,面积6 622.45 km2,2017年常住人口达833万人,是中国东部地区重要的中心城市、全国重要的科研教育基地和综合交通枢纽,城市化水平高。其地理坐标北纬31°14′—32°37′,东经118°22′—119°14′,属北亚热带湿润气候区,四季分明,雨水充沛。

图 1 研究区位置

1.2 数据

选用中分辨率的Sentinel-1A和Sentinel-2A数据,既可以保证地物内部的异质性较小及不同地物的可分性较高,又不会产生分辨率过高带来的同种地物光谱可分性(类内可分性)增加和不同地物的光谱可分性(类间可分性)降低等问题[15-16]。综合运用SAR和光学数据,波段范围较广,光谱特征丰富; 结合辅助数据,提高实验可靠性。

选用2017年4月21日和10月18日的Sentinel-1A Level-1级别干涉宽幅模式(TOPS Mode)的地距多视产品S1A IW GRD(Ground Range Detected)数据,方位向和距离向的分辨率均为10 m,已进行多视处理,只有幅度信息没有相位信息。以及2017年4月2日和10月9日的Sentinel-2A Level-1C级(L1C)多光谱数据(MSI),L1C级数据是经过几何精校正的正射影像,没有进行辐射定标和大气校正。选用10 m分辨率的2,3,4,8波段,即蓝、绿、红和近红外波段,数据无云层覆盖,影像质量较好。2种数据均在ESA网站上(https://scihub.copernicus.eu/dhus/#/home)下载。

参考第一次全国地理国情普查分类体系,结合实地调研采集样本点。根据研究区和数据集特点,制定本研究的分类体系:将房屋建筑(区)、道路和构筑物(包括水泥等硬化地表)定为建设用地,其余定为非建设用地。于2017年10月在整个研究区范围内展开实地调研,采集样本362个,其中新增建设用地样本130个,非新增建设用地样本232个。

此外,基于房天下网站(http://nanjing.fang.com/)等房地产家居网络平台,获取部分楼盘的开工时间、建成时间等信息作为验证参考。

2 新增建设用地提取方法

对两期Sentinel-1A和Sentinel-2A数据进行预处理并提取光谱特征、纹理特征和后向散射特征,进行多尺度分割、全域均值滤波和归一化得到组合特征集,通过随机森林、旋转森林、支持向量机和极限学习机多分类器集成学习提取新增建设用地。该方法将土地利用分类与变化检测有机结合,基于多时相特征信息,通过二值分类直接提取特定变化区域,简化了技术流程,提高了效率。同时流程的简化可以减少误差的产生及累积,从而降低系统误差,提高提取精度。整体技术流程见图 2

图 2 技术流程

2.1 哨兵数据特征提取

基于欧洲空间局(ESA)提供的SNAP6.0开源平台(https://step.esa.int/main/download/)对Sentinel-1A数据进行预处理,包括辐射校正、热噪声去除、Refined Lee滤波、SAR Mosaic和裁剪,得到后向散射特征VH和VV。

Sentinel-2A L1C是经过几何精校正的大气表观反射率产品,需要进行辐射定标和大气校正。处理基于ENVI5.5试用版进行,运用FLAASH模型进行大气校正,得到地表反射率。

精确配准对变化检测非常重要,可有效减少因配准误差而产生的伪变化区域[17]。基于ArcMap10.3平台对2个时期的影像进行配准,运用Georeferencing工具,以4月影像为基准影像对10月影像进行校正,在全图均匀选取17个控制点,运用二阶多项式进行校正,使总体误差控制在半个像元以内。预处理后的2017年4月和10月的Sentinel-2A数据B2、B3和B4波段真彩色合成图见图 3(a)(b)

图 3 Sentinel-2A真彩色合成图

利用ENVI5.3平台获取Sentinel-2A数据的光谱特征和纹理特征。利用Sentinel-2A的B2、B3、B4和B8波段计算遥感指数,包括归一化植被指数(Normalized Difference Vegetation Index,NDVI)、归一化水体指数(Normalized Difference Water Index,NDWI)和增强植被指数(Enhanced Vegetation Index,EVI); 根据波长范围对应关系,选择相应波段应用Landsat TM的变换系数进行缨帽变换(Tasselled Cap Transformation)得到土壤亮度指数(Soil Brightness)、绿度指数(Greenness)和黄度指数(Yellow Stuff); 进行主成分分析(Principal Component Analysis,PCA)得到第一主成分PCA1,包含最大的数据方差百分比和丰富的纹理信息,计算其灰度共生矩阵(Gray Level Co-occurrence Matrix,GLCM)得到8个纹理算子:均值(Mean)、方差(Variance)、协同性(Homogeneity)、反差性(Contrast)、非相似性(Dissimilarity)、熵(Entropy)、角二阶矩(Second Moment)和相关性(Correlation)。在此基础上,结合Sentinel-1A的后向散射特征VH和VV,所有特征见图 4

图 4 选用特征

基于ArcMap10.3平台,运用样本数据对所有特征进行采样,分析建设用地和非建设用地各个特征的取值规律,发现建设用地取值较高的特征有可见光、部分纹理特征和后向散射特征等,取值较低的特征有近红外、植被指数、水体指数和部分纹理特征等。建设用地和其他地类在所有特征中均存在差异,但部分特征具有相关性,且难以建立精细有效的判别规则,需要机器学习等方法进行深入分析。

2.2 耦合面向对象和集成学习的新增建设用地提取 2.2.1 面向对象分割

基于eCognition8.9平台,运用多尺度分割算法(Multiresolution Segmentation)进行面向对象分割[18]。经过反复试验,确定最优分割尺度为30,形状指数为0.1,紧密度为0.5,效果见图 5。可见同一地物大多被划分在同一对象之内,对象内部异质性较低,不同对象间差异性较大。

图 5 图像分割结果

由于2期影像获取时的太阳高度角、大气条件等不同,同类地物的光谱特性也会存在一定的差异,经过辐射校正后同一地物的光谱亮度也不尽相同,即使分割参数相同分割结果也不会完全一致,从而导致伪变化[18]。对此采用一种分割策略,即对第一期影像进行分割,将分割结果直接用于第二期影像,在此基础上进行变化检测和新增建设用地提取。

2.2.2 均值滤波

基于ArcMap10.3平台,对分割对象的所有特征逐个进行全域均值滤波,用滤波结果代替原特征进行分析。对于同一地物,滤波后的分割对象相比原始像元可以显著降低地物内部的异质性,使地物成为无差异整体。

2.2.3 集成学习

集成学习是一种多分类器集成算法[19]。分类器在特征空间中不同区域的性能存在差异,单一分类器难以适用于所有情况,易造成错分漏分现象。可结合其他分类器实现互补,通过多种分类器的综合分析提高效果[20]。为体现不同分类器的适用性和互补性,采用精度加权投票法进行集成[11],即先独立设计各个分类器,训练并检验其分类精度,以此作为各分类器的权重进行决策级集成。这样既可以体现各分类器的适用性,也可以降低直接投票法票数相同的概率。采用分类效果较优的随机森林(Random Forest,RF)[21-23]、旋转森林(Rotation Forest,RoF)[24-25]、支持向量机(Support Vector Machine,SVM)[26-27]和极限学习机(Extreme Learning Machine,ELM)[28-29]作为基分类器,算法差异大,可实现优势互补。

3 结果与分析

集成学习算法基于Matlab R2017a平台实现,随机选取50%的样本用于训练,剩余用作测试。为了分析比较方法的性能,共设置4种方法进行实验:(1)为提出的面向对象直接变化提取的方法; (2)基于像素尺度,不做均值滤波,直接变化提取; (3)对两期数据分别进行面向对象分类,通过变化分析提取新增建设用地,即分类后提取; (4)为基于像素的分类后提取。4种方法用于交叉比较面向对象与基于像素以及直接变化提取与分类后提取的差异。提取效果用整体精度(Overall Accuracy,OA)和Kappa系数来衡量,取值0~1,数值越大精度越高。最终结果见表 1(图 5)。

表 1 实验结果统计
方法 1 2 3 4
OA 0.950 3 0.922 7 0.828 7 0.850 8
Kappa 0.887 4 0.833 6 0.630 1 0.659 8

由OA和Kappa系数可知,本研究提出的方法1提取精度最高,其次是方法2和方法4,精度最低的是方法3。说明多时相直接提取的方法优于分类后比较的方法,可以有效降低系统误差,提高提取精度。直接提取方法中,面向对象的方法优于基于像元的方法,说明面向对象方法可有效降低同一地物内部的异质性和像元误分类的概率,从而提高精度。分类后比较的方法中,面向对象反而不如基于像元的方法,因为任一时相错分时,后续变化检测中整个地物都会错分,导致误差累积,而基于像元的方法不影响地物内部部分正确分类的像元。

基于方法1结果统计南京市2017年4—10月的新增建设用地面积共约117.038 km2,约占南京市总面积的1.767 3%。各区新增建设用地面积统计见表 2

表 2 南京市各区新增建设用地面积统计
行政区划 行政区面积/km2 新增建设用地面积/km2 新增建设用地占全市新增建设用地面积比例/% 新增建设用地占行政区面积比例/ %
江宁区 1 572.9 32.191 8 27.505 3 2.046 7
浦口区 912.3 24.518 6 20.949 2 2.687 6
溧水区 1 067.3 24.315 2 20.775 4 2.278 2
六合区 1 485.5 15.006 5 12.821 9 1.010 2
栖霞区 381.88 14.384 4 12.290 4 3.766 7
高淳区 802 7.591 4 6.486 3 0.946 6
雨花台区 134.6 5.545 2 4.737 9 4.119 8
建邺区 82.7 2.534 8 2.165 8 3.065 1
秦淮区 49.2 1.359 9 1.161 9 2.764
玄武区 80.97 0.898 1 0.767 4 1.109 2
鼓楼区 53.1 0.831 6 0.710 5 1.566 1

表 2可见,新增建设用地主要分布在江宁区、浦口区、溧水区、六合区和栖霞区,江宁区最多,达32.191 8 km2,鼓楼区、玄武区、秦淮区、建邺区、雨花台区和高淳区相对较少,鼓楼区最少,只有0.831 6 km2。考虑到行政区面积,从新增建设用地占行政区面积比例分析,雨花台区、栖霞区、建邺区、秦淮区和浦口区比例较大,高淳区、六合区、玄武区、鼓楼区、江宁区和溧水区相对较小。

4 结论

为准确提取新增建设用地,提出了一种基于多时相Sentinel-1A和Sentinel-2A数据的新增建设用地提取方法,在面向对象尺度上运用集成学习算法直接分类来提取特定变化类型。所提方法简化了变化检测流程,既不需要各个时相单独分类,也不需要生成差异数据等,而是基于两时相原始特征规律直接进行特定变化类型提取,减少了误差的产生环节和误差累积,从而降低系统误差,提高精度。

基于南京市的试验表明所提方法提取新增建设用地的整体精度可达0.95,Kappa系数达0.88,斑块提取完整,轮廓清晰。

后续研究一方面将扩展地物提取类型,增加时间节点,尝试进行多种变化类型的时间序列分析; 另一方面将综合分析城市建设用地变化与热岛效应、大气污染、工地扬尘、生态服务变化等的相互关联与作用,以期为城市生态环境监测与治理提供动态的信息支持。

参考文献
[1]
谈明洪, 李秀彬, 吕昌河. 20世纪90年代中国大中城市建设用地扩张及其对耕地的占用[J]. 中国科学(D辑), 2004, 34(12): 1157-1165.
[2]
SINGH A. Review article digital change detection techniques using remotely-sensed data[J]. International Journal of Remote Sensing, 2010, 10(6): 989-1003.
[3]
BRUZZONE L, BOVOLO F. A Novel Framework for the design of change-detection systems for very-high-resolution remote sensing images[J]. Proceedings of the IEEE, 2013, 101(3): 609-630. DOI:10.1109/JPROC.2012.2197169
[4]
佟国峰, 李勇, 丁伟利, 等. 遥感影像变化检测算法综述[J]. 中国图象图形学报, 2015, 20(12): 1561-1571. DOI:10.11834/jig.20151201
[5]
AUTHOR P C C, JONCKHEERE I, NACKAERTS K, et al. Review article digital change detection methods in ecosystem monitoring: a review[J]. International Journal of Remote Sensing, 2004, 25(9): 1565-1596. DOI:10.1080/0143116031000101675
[6]
张雨霁, 李海涛, 顾海燕. 基于决策树的面向对象变化信息自动提取研究[J]. 遥感信息, 2011(2): 91-94. DOI:10.3969/j.issn.1000-3177.2011.02.016
[7]
HUANG F, CHEN L, YIN K, et al. Object-oriented change detection and damage assessment using high-resolution remote sensing images, tangjiao landslide, three gorges reservoir, China[J]. Environmental Earth Sciences, 2018, 77(5): 183. DOI:10.1007/s12665-018-7334-5
[8]
HAN M, ZHANG C, ZHOU Y. Object-wise joint-classification change detection for remote sensing images based on entropy query-by fuzzy ARTMAP[J]. GIScience & Remote Sensing, 2018, 55(2): 265-284.
[9]
柳思聪, 杜培军, 陈绍杰. 决策级融合的多分辨率遥感影像变化检测[J]. 遥感学报, 2011, 15(4): 846-862.
[10]
杜培军, 柳思聪. 融合多特征的遥感影像变化检测[J]. 遥感学报, 2012, 16(4): 663-677.
[11]
WANG X, LIU S, DU P, et al. Object-based change detection in urban areas from high spatial resolution images based on multiple features and ensemble learning[J]. Remote Sensing, 2018, 10(2): 276. DOI:10.3390/rs10020276
[12]
LUO H, LIU C, WU C, et al. Urban change detection based on dempster-shafer theory for multitemporal very high-resolution imagery[J]. Remote Sensing, 2018, 10(7): 980. DOI:10.3390/rs10070980
[13]
周启鸣. 多时相遥感影像变化检测综述[J]. 地理信息世界, 2011, 9(2): 28-33. DOI:10.3969/j.issn.1672-1586.2011.02.007
[14]
TOUZI R, BOERNER W, LEE J, et al. A review of polarimetry in the context of synthetic aperture radar:concepts and information extraction[J]. Canadian Journal of Remote Sensing, 2004, 30(3): 380-407. DOI:10.5589/m04-013
[15]
BRUZZONE L, CARLIN L. A multilevel context-based system for classification of very high spatial resolution images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(9): 2587-2600. DOI:10.1109/TGRS.2006.875360
[16]
CARLEER A P, DEBEIR O, WOLFF E. Assessment of very high spatial resolution satellite image segmentations[J]. Photogrammetric Engineering and Remote Sensing, 2005, 71(11): 1285-1294. DOI:10.14358/PERS.71.11.1285
[17]
李德仁. 利用遥感影像进行变化检测[J]. 武汉大学学报(信息科学版), 2003(S1): 7-12.
[18]
董璐, 惠文华, 胡琰. 一种面向对象遥感变化检测的影像分割策略[J]. 遥感信息, 2016, 31(2): 80-85. DOI:10.3969/j.issn.1000-3177.2016.02.015
[19]
CHI M, KUN Q, BENEDIKTSSON J A, et al. Ensemble classification algorithm for hyperspectral Remote sensing data[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(4): 762-766. DOI:10.1109/LGRS.2009.2024624
[20]
刘培, 杜培军, 谭琨. 一种基于集成学习和特征融合的遥感影像分类新方法[J]. 红外与毫米波学报, 2014, 33(3): 311-317.
[21]
BELGIU M, DRAGUT L. Random forest in remote sensing: A review of applications and future directions[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 114: 24-31. DOI:10.1016/j.isprsjprs.2016.01.011
[22]
DU P, SAMAT A, WASKE B, et al. Random forest and rotation forest for fully polarized SAR image classification using polarimetric and spatial features[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 105: 38-53. DOI:10.1016/j.isprsjprs.2015.03.002
[23]
LI J, HUANG X, GAMBA P, et al. Multiple feature learning for hyperspectral image classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(3): 1592-1606. DOI:10.1109/TGRS.2014.2345739
[24]
XIA J, DU P, HE X, et al. Hyperspectral remote sensing imagec lassification based on rotation forest[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(1): 239-243. DOI:10.1109/LGRS.2013.2254108
[25]
XIA J, CHANUSSOT J, DU P, et al. Spectral-spatial classification for hyperspectral data using rotation forests with local feature extraction and markov random fields[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(5): 2532-2546. DOI:10.1109/TGRS.2014.2361618
[26]
KUO B, HO H, LI C, et al. A Kernel-based feature selection method for SVM with RBF kernel for hyperspectral image classification[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(1): 317-326. DOI:10.1109/JSTARS.2013.2262926
[27]
XIA J, CHANUSSOT J, DU P, et al. Rotation-based support vector machine ensemble in classification of hyperspectral data with limited training samples[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(3): 1519-1531. DOI:10.1109/TGRS.2015.2481938
[28]
HUANG G, ZHU Q, SIEW C. Extreme learning machine: Theory and applications[J]. Neurocomputing, 2006, 70(1-3): 489-501. DOI:10.1016/j.neucom.2005.12.126
[29]
HUANG G B, ZHOU H, DING X, et al. Extreme learning machine for regression and multiclass classification[J]. IEEE Trans Syst Man Cybern B Cybern, 2012, 42(2): 513-529. DOI:10.1109/TSMCB.2011.2168604