基于GEE的东北三省城市建设用地扩张研究

王明常，刘 鹏，陈学业，王凤艳，宋玉莲，刘瀚元

1.吉林大学地球探测科学与技术学院，长春 130026 2.自然资源部城市国土资源监测与仿真重点实验室，广东 深圳 518000 3.中国地质调查局西安矿产资源调查中心，西安 711500

0 引言

我国城市化发展迅速。东北作为老工业基地，20世纪80年代就有较高的城市化水平。卫星遥感影像能够直观准确地反映出过去几十年我国城市建设用地的扩张情况。GEE(Google Earth Engine)遥感云计算平台自免费开放以来得到了广泛的应用[1-2]，包括：水体变化监测[3-4]、土地利用变化监测[5]、畜牧业动态变化监测[6]、气候变化对农作物的影响研究[7]、城市水文评价多尺度不透水遥感监测[8]、洪泛平原湿地长期变化监测[9]和土地覆被监测[10]等。机器学习算法和云计算平台的发展为城市用地制图的发展带来了新机遇[11]。传统的大范围城市建设用地面积测绘费时费力，更新速度较慢，随着遥感云平台的出现和遥感影像数据的不断增加，基于GEE和遥感影像的城市建设用地提取研究逐渐趋于成熟，可快速提供大范围长时间序列更为可靠的数据，为制定应对气候变化、国土资源利用以及城市规划与管理政策发挥重要作用。

Liu等[12]利用GEE对全球建设用地分布信息进行提取，根据夜间灯光和Landsat数据综合多个遥感指数，构建归一化城市综合指数，基于大量地面样本数据，用最佳阈值分割对城市范围进行提取，实现城市用地的自动化识别。王博帅等[13]利用珞珈一号和NPP/VIIRS(suomi national polar-orbiting partnership/visible infrared imaging radiometer suite)夜光遥感影像比较了人类居住指数、植被覆盖和建筑共同校正的城市夜光指数、支持向量机3种提取城市建成区的方法，结果表明使用珞珈一号夜光遥感影像通过植被覆盖和建筑共同校正的城市夜光指数方法提取的城市建成区效果最好。Liu等[14]利用GEE和Landsat影像结合城市用地综合指数方法，绘制30 m分辨率1985—2015年全球城市动态图。Gong等[15]利用2018年的10 m分辨率哨兵卫星数据、OpenStreetMap、珞珈一号夜间灯光数据、腾讯移动定位和高德兴趣点(POI)等社会大数据，实现全国范围城市土地利用制图，制作了地块尺度面向对象的高分辨率城市土地利用图。Zhang等[16]提出了分层语义认知作为城市功能识别的一般认知结构，利用高分辨率卫星影像和兴趣点数据对北京城市功能区进行识别。以上研究都取得较高的精度和较好的效果，但大多数采用一次分类，对长时间序列大范围城市建设用地扩张研究较少。

为研究城市建设用地(路面、混凝土、砖、石及其他人造不透水覆盖层类型)扩张情况，本文以东北三省为研究区，基于GEE云平台，采用二次分类方法提取建设用地。首先利用夜间灯光数据和阈值分割法提取城市区域；其次构建特征数据并进行优化，利用人工解译和随机森林算法结合最优特征数据，对1989、1995、2001、2007、2013和2019年的Landsat去云月度合成影像中的建设用地进行提取；最后进行精度评价与比较，研究建设用地扩张时空变化情况，分析推力、弹力、阻力和支持力等驱动城市建设用地扩张的机制。

1 研究区概况和数据源

1.1 研究区概况

东北三省为辽宁省、吉林省和黑龙江省(图1)，地处中国中高纬度地区(118°50′E—135°05′E，38°43′N—53°33′N)，总面积约为91.57万km2，地形以平原、山地为主，主要属湿润、半湿润温带大陆性季风气候，年降水量300～1 000 mm，平均海拔高度为323 m。

图1 研究区概况

1.2 数据源与预处理

1.2.1 夜间灯光数据

NPP/VIIRS夜间灯光数据是由Suomi-NPP卫星的可见光/红外辐射成像仪(VIIRS)获取的夜间灯光遥感影像(DNB波段)，为全球每月的夜间灯光数据，该数据集有背景噪声、火光和极光等干扰。NPP/VIIRS夜间灯光数据亮度值每年呈增加趋势，选取2019年1—12月12期数据，采用平均值法对12期数据求年平均值，得到研究区年度平均夜间灯光数据。

1.2.2 多光谱数据

Landsat5地表反射率数据集是由Landsat5 TM传感器在大气中校正且经过正射校正的地表反射率，包含4个可见光、近红外波段和2个短波红外波段；Landsat8地表反射率数据集是由Landsat8 OLI/TIRS(operational land imager/thermal infrared sensor)传感器经大气校正和正射校正的地表反射率，包含5个可见光、近红外波段和2个短波红外波段以及2个热红外波段。本研究所使用的3 142景Landsat表观反射率(TOA)产品均来自GEE云平台，分别选取了1989、1995、2001、2007年477、462、560、534景Landsat5 TM影像和2013、2019年510、599景Landsat8 OLI影像，空间分辨率均为30 m。对所有TOA产品进行去云处理[17]，选取云量值小于30%的影像，其余影像被掩膜删除。此方法的优点是不会破坏波段原始值。为了提高分类精度并避免积雪对训练样本的影响，选取4—11月份影像，对去云影像进行镶嵌、裁切处理并生成东北三省月度合成影像，用所有月度合成数据生成影像数据集。

1.2.3 地形数据

SRTM(shuttle radar topography mission)数字高程模型产品是由美国“奋进”号航天飞机上搭载的SRTM系统采集后经过数据处理制成的数字高程模型。本研究使用了30 m分辨率的SRTM1产品。

1.2.4 经济和人口数据

本文中经济和人口数据来自《中国统计年鉴》[18]、《辽宁统计年鉴》[19]、《吉林统计年鉴》[20]和《黑龙江统计年鉴》[21]。

2 研究方法

通过借鉴国内外城市土地利用信息提取的研究经验，基于GEE云平台，利用阈值分割法提取城市区域，构建特征数据并进行特征优化，合理选取训练样本，将最优特征作为输入对象，再利用随机森林算法提取建设用地，并进行精度评价与比较。技术流程如图2所示。

图2 技术流程

2.1 城市区域提取

基于NPP/VIIRS夜间灯光数据，利用阈值分割法提取城市区域。2019年度平均夜间灯光数据(图3)灰度值(Gdn)的范围为(0，168)，根据观察可知灯光存在溢出现象。经实验分析与比较，选取Gdn=0.8作为分割阈值划分城市区域边界，对非城市区域进行掩膜处理得到城市区域。先确定城市区域再提取建设用地，避免了林地、裸地和盐碱地等其他无关地物类别对提取建设用地的影响。

a.东北三省；b.沈阳；c.长春；d.哈尔滨。

2.2 建设用地提取

2.2.1 样本数据选取

Collect Earth是联合国粮农组织开发的开源目视解译工具，可以通过时间序列信息调用Landsat影像数据，利用关键物候期的彩色合成图像目视解译出更准确的土地利用信息[22]。将研究区划分为建设用地和非建设用地，基于Collect Earth利用目视解译的方法得到5 280个样本数据，2/3样本用于训练，1/3样本用于精度评价。

2.2.2 特征数据构建

光谱特征是地物分类的主要特征，在Landsat月度合成影像中选取蓝波段、绿波段、红波段、近红外波段、短波红外1波段和短波红外2波段作为光谱特征。

通过构建遥感专题指数可以有效提高分类精度[23]。经过处理后的Landsat 数据为月度合成影像，通过原始光谱波段构建指数特征：归一化植被指数(normalized difference vegetation index, NDVI)、增强型植被指数(enhanced vegetation index, EVI)[24]、改进归一化水体指数(modified normalized difference water index, MNDWI)[25]和归一化建筑指数(normalized difference built-up index, NDBI)[26]，计算公式分别为：

(1)

(2)

(3)

(4)

式中：Indvi、Ievi、Imndwi、Indbi分别为NDVI、EVI、MNDWI、NDBI的指数值；Rnir、Rr、Rb、Rg、Rmir分别为近红外、红、蓝、绿、中红外波段的反射率。

纹理作为图像的重要属性应用于图像分类中。由于短波红外波段用于分辨道路、水和植被等更有效，因此使用Landsat短波红外波段构建灰度共生矩阵(gray-level co-occurrence matrix, GLCM)，利用GEE云平台中的glcmTexture函数可计算出18个纹理特征。本研究选取最常用的6个纹理特征：角二阶矩、对比度、相关性、方差、逆差距和熵[27]。

本研究使用SRTMGL1_003数据构建地形特征：海拔高度、坡度、坡向和山体阴影。

2.2.3 特征数据优化

太多的特征数据会导致运算效率和分类精度的降低。SEaTH算法是一种有效的特征优化算法，可筛选出最优分类特征[28]。SEaTH算法根据地物特征的正态分布，利用Jeffries-Matusita距离(简称J-M距离，用J表示)衡量两类地物C1和C2间的可分性。J的计算公式如下：

J=2(1-e-B)；

(5)

(6)

2.2.4 随机森林算法分类

随机森林(random forest, RF)算法是Bagging(bootstrap aggregating)的一个扩展变体，在以决策树为基学习器构建Bagging集成的基础上，进一步在决策树的训练过程中引入随机属性选择。RF算法的基学习器是决策树。RF学习器可以分类也可以回归，在最终的分类时可以输出概率。随着个体学习器数目的增加，RF算法通常会收敛到更低的泛化误差[30]。

RF分类器即通过构造多个训练集得到决策树，通过T轮训练得到分类器序列{h1(x),h2(x),…,hT(x)}，构成一个多决策树的分类器组合，通过学习器投票的置信度来确定输出结果[31]。采用Bagging抽样策略，使用2/3的样本用于分类器训练和1/3的样本用于分类精度评估，则有

(7)

式中：H(x)为样本x的分类结果；ht为单个决策树分类器；y为输出变量；Ⅱ(*)为指示函数。RF算法泛化误差的上界满足

(8)

决策树数目对分类精度有重要的影响[32]。通过综合考虑分类精度和运行速度，将决策树数目设置为n=400，利用RF算法在城市区域范围中提取出建设用地。

2.2.5 精度评价方法

利用验证样本数据构建混淆矩阵，评价建设用地与非建设用地的分类精度，具体评价指标包括总体精度(overall accuracy, OA)和Kappa系数[33]。

2.3 城市空间扩张度量

通过城市建设用地扩张面积(ΔS)、扩张速率(v)和扩张强度(K)3个时空特征测度指标来分析城市建设用地扩张时空变化，计算公式分别为：

ΔS=Sb-Sa；

(9)

(10)

(11)

式中：Sa和Sb分别为研究期间前期和后期的城市建设用地面积；ΔT为时间间隔。城市扩张强度指数是对扩张面积的年平均化，也是对扩张速率的标准化，消除了城市规模初始状态对扩张特征的影响。按照扩张强度可将城市扩张类型分为5类[34]：高速扩张型(K>20%)、快速扩张型(14%

3 结果与分析

3.1 特征优化结果

通过原始光谱波段和特征构建得到6个光谱特征、4个指数特征、6个纹理特征和4个地形特征共计20个特征数据(图4)，将3种类别(1.建设用地；2.水体；3.植被)进行两两组合，计算每种类别组合和每种特征的J-M距离(表1)。由于只提取建设用地，所以选择表1类别组合1-2、1-3中J>1的特征数据，最终得到12个最优特征数据(表2)。

表2 特征优化结果

海拔高度的单位为m；坡度、坡向的单位为(°);山体阴影为灰度值。

表1 类别组合和特征的J-M距离

3.2 建设用地提取精度评价与比较

利用OA和Kappa系数对不同年份建设用地提取结果进行评价(图5)，结果表明各年份提取结果都有较高的精度。东北三省平均OA和Kappa系数分别为95.88%和0.92，沈阳平均OA和Kappa系数分别为96.86%和0.94，长春平均OA和Kappa系数分别为96.08%和0.92，哈尔滨平均OA和Kappa系数分别为95.96%和0.92。Gong等[35]提供了1978—2017年中国建设用地制图产品，利用本文验证样本对其中2013年长春市建设用地进行精度评价，平均OA和Kappa系数分别为91.56%和0.88。本文的提取结果与Gong等提供的建设用地制图产品对比如图6所示，其中b为本研究建设用地提取结果轮廓叠加到Landsat影像上的显示效果。可以看出，与Gong等提取结果相比，本文所用方法提取结果精度较高而且效果较好。

图5 研究区精度评价

a.2013年长春市部分区域Landsat影像；b.本研究建设用地提取轮廓；c.本研究提取结果；d.Gong等[35]提取结果。

3.3 省会城市建设用地扩张时空变化

沈阳、长春和哈尔滨城市建设用地扩张时空变化如表3所示。可以看出：沈阳在2001—2007年为低速扩张型，长春和哈尔滨在2007—2013年为低速扩张型，其他时间3个城市都属于缓慢扩张型；1989—2019年沈阳建设用地扩张面积较大，而且扩张平均速率较快，长春次之，哈尔滨最慢。主城区30年的扩张情况如图7所示，沈阳主要向西、南方向扩张，长春向四周扩张，哈尔滨主要向西北方向扩张。

a.沈阳；b.长春；c.哈尔滨。

表3 扩张面积、扩张速率和扩张强度统计表

3.4 城市建设用地、经济和人口变化

1989、1995、2001、2007、2013和2019年东北三省所有城市以及长春、沈阳和哈尔滨建设用地提取结果如图8所示，城市建设用地面积、经济和人口变化如图9所示。东北三省所有城市30年建设用地面积总共扩张49.07%，经济增长25.93倍，人口增加11.12%；沈阳建设用地面积扩张2.25倍，经济增长28.74倍，人口增加47.42%；长春建设用地面积扩张2.09倍，经济增长63.39倍，人口增加20.30%；哈尔滨建设用地面积扩张1.80倍，经济增长31.85倍，人口增加9.02%。2013—2019年，东北三省GDP总量减少7.7%，总人口减少1.66%；哈尔滨人口减少4.4%且经济增长缓慢，沈阳、长春经济和人口虽呈增长趋势，但人口增长较慢，相应城市建设用地面积也增加缓慢。

a.东北三省；b.沈阳；c.长春；d.哈尔滨。

a.东北三省；b.沈阳；c.长春；d.哈尔滨。

3.5 城市建设用地扩张驱动机制

城市建设用地扩张驱动机制主要包括推力、支持力、弹力和阻力[36]。推力是推动城市扩张进程的主要因素，主要包括人口、经济和交通等因素。由上述东北三省城市建设用地、经济和人口变化情况可知，人口增加的速度和经济发展的速度与城市建设用地扩张速度变化规律基本一致。分析其原因主要为：人口增加会带来居住、就业、医疗、教育和娱乐等需求的增加，从而会加大对城市基础设施及相关工程的建设，经济的快速发展也会促进基础设施的扩建。而人口的减少和经济的下降则会导致城市建设用地扩张速率和强度的降低，甚至会导致城市建设用地收缩。人口因素和经济因素是促进城市扩张的主要因素。

支持力是指为城市扩张提供便捷条件的因素，包括自然环境和行政因素。丰富的矿产资源、水资源和肥沃的黑土地是东北三省城市发展的自然基础，三江平原、松嫩平原和辽河平原为城市发展建设提供了优越的地理条件。东北振兴战略是东北三省城市扩张重要的行政因素。

弹力作为推动城市扩张的动力作用于城市不同区域，包括开发区、高新区和大学城等的建设，影响城市化进程和城市扩张模式。

阻力是指阻碍城市扩张的因素，包括生态环境保护和文化遗址保护等。为了保护生态环境，东北三省一些城市放弃了很多盈利企业发展的机会；城市建设也要考虑到对文化遗址的保护。阻力因素的作用导致了城市扩张速率和强度的降低。

4 结论及展望

本研究基于GEE云平台，利用NPP/VIIRS年度平均夜间灯光数据和阈值分割法提取了城市区域；利用Landsat去云月度合成影像制作不同年份影像集，构建了指数特征、纹理特征和地形特征；利用SEaTH算法对特征数据进行了优化；利用RF算法分别提取了东北三省1989—2019年每6年时间间隔的城市建设用地，并进行精度评价和比较；分析了省会城市建设用地扩张时空变化情况，东北三省城市建设用地、经济和人口的变化以及城市建设用地扩张驱动机制。得到以下结论：

1)基于Landsat月度合成影像可以减少云和雪等的影响，获得较好的建设用地提取效果。利用夜间灯光和Landsat数据通过二次分类先提取城市区域再提取其中建设用地的方式，避免了林地、裸地和盐碱地等其他无关地物类别对分类的干扰，提高了分类效率和精度。

2)利用SEaTH算法进行特征优化，根据J-M距离的取值将特征数量从20个精简到12个，降低了运行载荷的同时提高了分类精度。本研究平均总体精度和Kappa系数分别为96.19%、0.92，有较高的精度和可靠性。

3)东北三省城市建设用地的扩张是复杂的动态变化过程，省会城市中沈阳扩张速率较快，长春次之，哈尔滨最慢。推力、支持力、弹力和阻力是驱动城市建设用地扩张的主要机制，人口因素和经济因素是推动城市建设用地扩张的主要因素，自然基础和行政因素提供了便捷条件，开发区、高新区和大学城等的建设影响城市扩张模式，生态环境和文化遗址保护会导致城市建设用地扩张速率和强度的降低。

利用遥感云平台可以进行长时间序列大范围城市建设用地提取，提取成果可为城市时空变化研究提供借鉴和参考。本研究仅讨论了城市建设用地的提取，对建设用地中更详细土地利用信息(如商业区、住宅区、工业区等)提取还需寻找更加有效的技术方案。各种遥感云平台和遥感数据的出现，为我们进行城市土地利用信息的提取提供了多种选择，以便获取高精度、高分辨率和更加丰富的城市土地利用信息数据产品。