基于GEE遥感大数据云平台的城市扩张研究——以南昌市为例

冉 玥，查夏凌云，赵丽红

(江西农业大学国土资源与环境学院，330045，南昌)

0 引言

城市是人类活动的聚集中心，社会、政治、经济、文化、自然等多因素共同作用聚合的产物。土地作为人类生存和发展的物质基础，土地利用分类信息在资源开发利用、城市扩张规划与研究、生态环境评价与保护等方面有着广泛的应用[1-2]。城市扩张在一定程度上反应了一段时期内地域经济的发展以及人口的增长和迁移，因此对城市扩张的研究对于社会发展及规划有着重要的启示意义。利用遥感图像对土地利用进行监测和分类是一种常用的研究方法。

遥感技术的发展使得遥感数据的数量和类型更多更复杂，在应用的过程中容易产生数据冗余，操作的过程也比较繁琐。基于云计算平台的遥感大数据研究方法有着高效率、低成本、高稳定性、易更新等特征。随着近年来国内外对地理云计算研究的深入，国内市场上也涌现了一批以PIE-Cloud遥感云平台、“四维地球”为代表的功能强大的产品。Google Earth Engine是由谷歌、卡内基梅隆大学和美国地质调查局联合开发的云计算平台，是目前世界上最先进的PB级别地理数据科学分析及可视化平台[3]。蒲东川[4]等使用GEE平台，综合利用光谱特征和独立成分分析后的Landsat8OLI影像特征对北京市城镇用地进行提取；戴声佩等[5]基于GEE平台提供的Landsat8 OLI影像，分别采用随机森林和支持向量机分类算法两种不同的算法对海南岛土地利用类型进行了分类研究，结果表明随机森林算法对各类地物信息的识别更加精确；方梦阳等[6]基于GEE平台提供的黄河流域近千景Landsat影像，采用随机森林分类算法对其进行监督分类并制作了黄河流域20年逐年土地覆盖数据。鉴于此，本研究以南昌市为研究区，基于GEE平台利用夜间灯光数据对城市建成区进行提取，再选取2000—2020年间5期历史影像数据监测城市土地利用动态变化，结合经济人口数据进行剖析，对城市化扩张的驱动力进行探析。

1 研究区概况和数据源

1.1 研究区概况

南昌位于我国的华东地区，江西中部偏北，位于东经E115°27′～116°35′、北纬N28°09′～29°11′之间。全境以平原为主，地形东南相对平坦，西北丘陵起伏，水网密布，气候湿润温和，属亚热带季风区。全市总面积7 194.98 km2，平原占35.8%，水面占29.8%，辖6区3县，设1个国家新区(赣江新区)。2020年全市生产总值(GDP)达到5 745.51亿元，常住人口达到538.29万人。

图1 研究区概况

1.2 数据来源与预处理

1.2.1 NPP/VIIRS夜间灯光数据 “索米”国家极轨卫星伙伴卫星(Suomi National Polar-obrating Partnership or Suomi NPP)于2011年发射，其搭载的VIIRS传感器上的DNB(Day Night Bands)波段能在500 m分辨率的尺度上对地表开展每天覆盖一次的高灵敏度的夜光观测，可以提供适用于地球系统的夜间可见光和近红外光的全球日常测量科学和应用。夜间灯光数据反映了夜晚地球表面的分区灯光强度，可以反映城市甚至小规模居民聚集地、车流等发出的微弱亮光，因此可以明显区分城区和非城区[7]。本文选取2021年1—12月12期数据，使用GEE平台的ImageCollection.median()函数对夜光数据集合创建中值图像，通过计算所有匹配波段堆栈中的每个像素的所有值的中值实现。

1.2.2 Landsat5/8 TOA数据 Landsat5 TOA和Landsat8 TOA数据集是分别由Landsat5 TM传感器和Landsat8 OLI/TIRS(operation land imager/ thermal infrared sensor)传感器经过大气校正和正射校正的地表反射率。本文选取GEE平台上的2000年、2005年、2010年、2015年、2020年Landsat5/8 TOA数据集进行研究，空间分辨率均为30 m。使用GEE提供的方法对影像进行时间和研究区的选择，进行去云处理后，使用从夜光数据中提取的建成区进行掩膜，得到研究区。

1.2.3 SRTM数据 航天飞机雷达地形任务(The Shuttle Radar Topography Mission，SRTM)数字高程模型产品由NASA提供，空间分辨率为30 m。

1.2.4 样本数据 样本通过加载不同的遥感影像图层目视解译获取，共得到641个样本数据。2000—2010年的样本通过去云处理后的Landsat5 TOA影像，辅助以NDVI图像进行选取。2015年及以后的样本选取依据影像为Sentinel-2数据的Level-1C产品，该产品提供正射校正的大气顶层反射率，具有大范围、高分辨率和多光谱成像特征，其中QA60 波段可以对影像中的卷云和厚云进行消除，达到云掩膜的效果。分别使用GEE平台的ee.Filter.lt(name, value)和ee.Filter.lte(name, value)方法随机选择70%的样本作为训练样本，30%为验证验本。

1.2.5 经济人口数据 本文经济和人口数据来自于2021年《江西省统计年鉴》和《南昌市统计年鉴》。

2 研究方法

基于GEE云平台，通过夜光数据划分城市区域边界，考虑到耕地及其他地类占地较少，将城市土地利用地类分为建筑、植被、水体、裸地4类，依次进行特征构建、特征优化、样本选取后，利用随机森林的方法对城市建成区土地利用分类与精度评价。技术流程如图2。

图2 技术流程

2.1 城市建成区提取

利用夜间灯光数据提取城区的方法有很多种，如聚类阈值法、统计分析法，突变检测法，经验阈值法等。这些方法存在主观性较强、稳定性较差、工作量大和工作效率较低等问题。由于灯光的溢出性，灯光范围可能超出实际建成区边界，因此阈值的确定很大程度上影响了提取的结果[8-9]。

最大类间方差法是一种自适应阈值计算方法，又称为大津算法(OTSU算法)，由日本学者大津展之(Nobuyuki Otsu)于1979年提出[10]。这种算法的目的是在应用阈值对图像进行二值化分割后，使得目标与背景之间的类间方差最大。最大类间差算法的优越性在于它能够自动确定合适的阈值，且在GEE平台的实现简便。对于图像I(x,y),T为目标和背景的分割阈值，即使类间方差最大的阈值，由遍历的方法得到：

(1)

(2)

N0+N1=M×N

(3)

ω0=ω1=1

(4)

μ=ω0μ0+ω1μ1

(5)

g=ω0(μ0-μ)2+ω1(μ1-μ)2

(6)

g=ω0ω1(μ0-μ1)2

(7)

式中：ω0和ω1分别表示目标和背景的像素点数占整幅图像的比例，μ0和μ1分别为其平均灰度，图像的总平均灰度记为g，μ表示类间方差。

本文通过NPP/VIIRS数据的“avg_rad”波段，用最大类间方差法进行处理，最终将城市分为城市建成区和非建成区(图4)。

图3 夜间灯光数据(2021,南昌)

图4 建成区提取结果

2.2 特征构建

2.2.1 光谱指数特征 光谱指数特征在遥感影像分类中应用广泛，本文选取Landsat5/8 TOA影像中的6个原始波段(蓝波段、绿波段、红波段、近红外波段、短波红外1波段和短波红外2波段)和4个遥感专题指数作为光谱特征来提高分类精度。调用GEE平台的Image.normalizedDifference(bandnames)函数进行计算，表达式为：(B1-B2)/(B1+B2)。该函数通过原始光谱特征波段分别计算归一化植被指数(NDVI)，改进归一化水体指数(MNDVI)，归一化建筑指数(NDBI)；运用GEE的Image.expression(map,expression)函数计算增强型植被指数(EVI)。上述指数的计算公式如表 1，其中Rnir、Rr、Rb、Rg、Rswir分别为近红外、红、蓝、绿、短波红外波段的反射率。

表1 植被指数

2.2.2 纹理特征 纹理是图像的一个重要属性，是图像灰度在空间上以一定形式变化而产生的图案。灰度共生矩阵(Gray-Level Co-occurrence Matrix，GLCM)是一种被广泛应用的纹理统计分析方法[11]，稳定性较强，在影像分类中有着广泛的应用。18个纹理特征可以由GEE提供的glcmTexture(size, kernel, average)方法快速计算。由于归一化植被指数(NDVI)可以有效区分植被与非植被，利用Landsat影像计算归一化植被指数(NDVI)，基于NDVI构建纹理特征。为了避免冗余，选取研究中最常用的6个特征：角二阶矩(B8_asm)、对比度(B8_contrast)、相关性(B8_corr)、方差(B8_var)、熵(B8_ent)、逆差矩(B8_idm)。

2.2.3 地形特征 基于GEE平台ee.Terrain.products(input)方法计算海拔高度、坡度、坡向、山体阴影，以此构建地形特征。

2.3 特征优化方法

随机森林分类方法涉及到地物的光谱、形状、纹理和地形特征，可用于构建光谱特征指数的特征数量庞大。当使用全部特征输入分类器时，运算的效率和精度都会降低，在有限样本的情况下，会导致“维数灾难”[12]。因此需要进行特征优化，选择对分类精度更有效的特征。本文使用SEaTH算法对特征进行优化，该算法根据特征相关性，类间距离和类内距离对特征进行评价，然后综合利用多种评价结果获取最优的特征子集。在SEaTH算法中，利用Jeffries-Matudita距离(简称J-M距离，本文简称为J)来衡量类间可分性。分别设2个不同的类别为C1和C2，J-M距离可以表示为：

J=2(1-e-B)

(8)

(9)

式中：B表示距离计算因子，mi和σ^2i(i=1，2)分别代表类C1和C2某特征分布的均值和方差。J的取值范围为[0-2]，2表示C1和C2在某个特征上可以完全分开；反之，0表示这2个类在这个特征上可以完全混淆。J越接近2，可分性越好；J越接近0，可分性越差[13-14]。

2.4 随机森林分类

随机森林(Random Forest)是集成学习中一种典型的Bagging算法，在此基础上综合了Cart决策树和特征随机选取的思想。从原始样本中采用Bootstrap重抽样技术随机抽选样本，再采用节点的随机分裂技术构建大量决策树。由于每一颗决策树都从数据集中随机采样进行训练，每颗决策树输出的结果同时具有相似性和差异性，使得最终得到的结果高效且稳定。除此之外，在抽样过程中未被抽取的样本称为袋外(Out of bag，OOB)数据，利用部分OOB数据可以对分类结果进行精度评价和特征重要性评价。随机森林算法需要定义的参数有2个：决策树的数量(Ntree)和构建分类树时从所有特征中随机选取的Mtry个特征。根据研究，随机抽取的特征变量个数通常为特征变量总数的平方根，当Ntree≥100时，各分类情况的OOB误差趋于稳定，Ntree一般不大于1 000。考虑到分类精度和运行速度，设置Ntree=100[15-17]。

2.5 精度评价

混淆矩阵(Confusion Matrix)又称为误差矩阵，是一种n行n列的矩阵(n为类别数)，在研究中常用来评价精度[18]。运用GEE平台的ImageCollection.errorMatrix(actual, predicted, order)函数可以利用验证样本计算混合矩阵，得到总体精度(Overall Accuracy，OA)和Kappa系数，作为精度评价指标。

P0=Tr/N

(10)

式中：总体分类精度用Po表示；Tr为正确分类的像元个数；N为研究区像元个数。

Pe=(a1×b1+a2×b2+…+ac×bc)/N×N

(11)

K=(P0-Pe)/1-Pe

(12)

式中：K为Kappa系数；Po代表总体分类精度。

2.6 城市土地利用变化和扩长度量

本文选用土地利用变化面积△S、土地利用变化速率v、土地利用扩张强度K、土地利用度LUD4个指标来衡量城市土地利用变化和扩张，计算公式分别为：

△S=Sb-Sa

(13)

(14)

(15)

(16)

其中：Sa和Sb分别为研究初期和末期的建筑用地面积，△T为a―b时刻的研究时段间隔，Dab为从a―b时刻新开发的建筑用地面积，Li是区内第i类土地利用强度分级指数(表 2)；Ai是第i类土地利用类型在区内所占的百分比[19]。依照扩张强度将城市扩长类型分为5类[20](表 3)。

表2 土地利用强度分级表

表3 城市扩张强度分级表

利用土地利用转移矩阵评估区域土地利用变化的结构特征与各用地类型变化的方向。该方法来源于系统分析中对系统状态与状态转移的定量描述，在研究中应用广泛[21]。

3 结果与分析

3.1 特征优化结果

将4种类别(C1建筑，C2裸地，C3植被，C4水体)两两组合，对特征构建得到的6个原始光谱特征、4个指数特征、6个纹理特征和4个地形特征通SEaTH算法进行优化，分别计算J-M距离。从20个特征(图5)中选取最接近2的12个最优特征(表 5)。本文主要对城市扩张进行研究，故只选取类别组合为建筑-裸地(C1-C2)、建筑-植被(C1-C3)、建筑-水体(C1-C4)计算J-M距离。结果如表 4。

图5 不同特征数据的均值

表4 不同类别组合的J-M距离

表5 特征优化结果

3.2 分类结果与精度评价

利用OA和Kappa系数对5年的分类结果(图 5)的精度进行评价(图 6)。结果表明5期不同年份的分类结果的精度较高，平均OA和Kappa系数分别为96.81%和94.56%。

图6 2000—2020年土地利用覆被分类结果

图7 分类结果精度评价

统计可知，2000年、2005年、2010年、2015年、2020年建筑用地面积分别占南昌市建成区总面积的31.41%、43.17%、51.43%、69.77%、79.42%。由图6可知，从2000—2020年，城市扩张趋势表现为以南昌市老城区(包括东湖区和西湖区)为中心，向西部地区(新建区)、南部地区(南昌县)、东北部地区(青山湖区)扩张，呈现“南延北控，西进东拓”的扩张特征。由于从灯光数据中提取到的进贤县和安义县建成区所占面积较小，且经济发展水平较其他区域较落后，本文不将上述二县纳入主要分析范围。

3.3 城市土地利用与扩张分析

南昌市20年建筑用地变化面积、变化速率、扩张速率如表 6。在2000—2005年、2010—2015年属于低速扩张型，在2005—2010年、2015—2020年属于缓慢扩张型。

表6 建筑用地变化面积、变化速率与扩张速率

2000—2010年和2010—2020年的土地利用转移矩阵分别如表7和表 8，在提取的建成区范围内，由植被和裸地转为建筑用地的面积最多，在2000―2010年转换的面积分别占建成区总面积的17%和10%，在2010―2020年转换的面积则分别占建成区总面积的16%和14%。

表7 土地利用转移矩阵(2000—2010年)

表8 土地利用转移矩阵(2010—2020年)

2000年、2005年、2010年、2015年、2020年南昌市建成区的土地利用度依次为：2.41、2.65、2.85、3.31、3.55，总体呈逐年上升的趋势。大量的裸地、植被和水体转换为建筑用地。从2000—2005年，在青云谱区，位于南昌西郊的象湖附近区域、以瑶湖、艾溪湖为中心的青山湖区的东部地区，以及以南昌县为代表的北部地区为城市扩张的主要区域，表现为建筑用地增加面积大；从2005―2015年，南昌建城区逐渐开始向西部扩张，主要集中在红谷滩区和新建区；从2015年开始，建筑用地的开发主要集中在赣江新区和九龙湖区域。

随着社会的发展，居民生活水平的提高刺激了对第二产业和第三产业的需求，尤其是文化教育、娱乐饮食、社会保险等服务类行业的需求。学校、综合体、医疗服务、房地产的进驻和体系优化往往吸引了人口的流动，带来了经济的繁荣，人口与经济的相互作用共同推进了城市化。南昌市20年人口总体呈增长趋势，增加24.44%；经济发展迅猛，增长12.06倍，其中，从2000—2010年这10年的发展尤为显著，地区生产总值增长4.50倍，从2000―2020年经济稳步发展，地区生产总值增长2.67倍；相应的建筑面积随着经济人口的上涨而上升，增加2.52倍。

南昌市近10年各区县常住人口数变化如图 8，老城区(东湖区、西湖区)常住人口数呈下降趋势，其中，东湖区2020年的常住人口比2011年下降了11.41%；除这3个区域外其他区域均呈上涨趋势。其中，2020年红谷滩区的常住人口数比2011年增长了1.29倍，高新开发区和经济开发区的增长率则分别达到了82.08%和59.18%。整体看来，老城区的人口正在逐渐向郊区迁移，大学城、工业、高新产业和资本的引入在很大程度上引导了人口的流动。

图8 2000—2010年土地利用变化图

图9 2010—2020年土地利用变化图

图10 研究区建设用地、经济和人口直方图

图11 近10年南昌市各区县常住人口变化

4 结论与展望

本研究基于GEE平台，利用NPP/VIIRS年度夜间灯光数据对南昌市城市建成区提取，对Landsat5/8 TOA年度数据(2000年、2005年、2010年、2015年、2010年)构建光谱指数特征、纹理特征、地形特征等指数，基于SEaTH算法对特征进行优化，最后采用随机森林算法进行土地利用覆被分类，并进行精度评价，主要结论如下。

1)利用最大类间方差法和夜光数据可以提取较为精准的建成区范围，采用随机森林算法，基于优化后的特征进行分类可以得到较为准确的分类结果，平均总体精度和Kappa系数分别为96.81%和94.56%，与传统遥感分析操作方式相比效率更高。

2)南昌市建成区在2000—2010年、2010—2020年2个时间段分别经历了从低速到缓慢的扩张，扩张保持了与规划一致的“西进、东拓、北控、南延”趋势，土地利用度的逐年缓慢上升反映了土地利用规划对土地资源合理利用的控制效果，经济发展和人口的稳步上升促使建设用地的增加。

3)基于GEE平台的数据处理方式与传统模式相比在流程上省去了大量繁琐的重复操作，在保证精度的前提下提高了数据处理的效率，处理过程透明可视，代码可复用性强，在纠错和优化上有极大的优势。

本文仍存在不足，灯光数据提取精度有待提高，部分位于郊区的城区未被精确提取；考虑到从灯光数据中提取出的建成区所含耕地面积较小，没有将耕地加入分类，对土地利用度的计算有一定的影响；其次，随机森林分类方法仍存在“椒盐现象”，在后续的研究中可以采用其他方法解决该问题。