城市是人类在适应并且改造自然过程中建立起来的“自然-经济-社会”三者为一体的复杂生态系统。城市作为社会经济活动的中心区域，人口迁移和土地利用的快速扩张构成了全球性挑战。重庆市作为西南唯一的直辖市，有着特殊的地形地貌，与平原地区的演变相比，其城市化发展会有一定的特殊性。遥感数据资源和地理信息系统独特的空间分析能力，为城市建成区监测提供了充足的数据与技术支持^[1]。卫星遥感影像客观性强，国内外学者在运用遥感数据开展城市扩展研究方面有了许多积累。Masek等^[2]利用TM数据研究了美国华盛顿特区城市扩展的动态变化；Yang等^[3]利用多时相遥感影像分析了美国亚特兰大的土地利用/土地覆盖状况；蒲东川等^[4]利用独立成分分析和随机森林法对北京市城镇用地进行提取。但目前，中国城市研究区主要集中在东部经济发达的平原地区，相对而言，西部山地城市的研究案例较少。

基于长时间序列的制图必然要使用大量遥感影像，而传统的遥感影像处理软件（如ENVI、ERDAS等）需要将影像下载到本地后才能对其逐一进行处理，耗时费力。谷歌地球引擎（Google Earth Engine，GEE）平台是由Google公司、美国卡内基梅隆大学以及美国地质调查局联合开发的云端运算平台^[5]，具有对卫星影像数据和其他地球观测数据进行专业处理的能力^[6]。自2010年发布以来，GEE平台提供了海量的光学和雷达影像数据集，包括Landsat、MODIS、哨兵系列、陆地观测卫星等，此外还有一些现成的数据产品，常见的有植被指数（vegetation index，VI）、Hansen全球森林变化和数字高程模型（digital elevation model，DEM）数据。同时，GEE具备数据共享、编码存储、并行处理以及可视化的能力，支持JavaScript和Python编码语言的开发环境，包含了许多可供编程人员使用的应用程序接口（application program interface，API）库^[7]。目前，国内外学者已将GEE平台应用到多个领域。其中，作物制图相关研究数量最多，如农作物面积^[8]和农作物种类^{[9, 10]}等农业监测。随后也有相关的水体研究，如地表水^[11]、湖泊^[6]、冰川^[12]等监测。土地覆盖类型制图研究^{[13, 14]}、森林研究^{[15, 16]}、城市研究^[17]也逐渐成为近年来的热点。

因此，本文基于GEE云平台的Landsat TOA数据，以重庆市主城区为例，结合遥感分类和归一化建筑指数（normalized difference built-up index，NDBI）提取1999—2018年的城市建成区时间序列，通过探索性空间数据分析（exploratory spatial data analysis，ESDA）法进行相关分析，探究重庆市城市扩张的空间变化状况。

1 研究数据与方法 1.1 研究数据

重庆市地处长江上游，东临湘鄂，西靠四川，南接贵州，北连陕西。重庆市的山地面积超过全市总面积的3/4，坐拥长江、嘉陵江、乌江3条大江。地形由南北向长江河谷倾斜，构成“一岛、两江、三谷、四山”的城市格局。本文以重庆市主城区为研究区域，主城区总面积约为5 472.68 km²，截至2018年底，国内生产总值（gross domestic product，GDP）总值为8 208.39亿元，常住人口875.00万人，主城区城镇化率达到90.51%^[18]。

原始遥感影像为GEE平台提供的空间分辨率为30 m的Landsat系列1999—2018年的遥感数据。为保证研究区内云量最低，在GEE平台用Python语句Filter筛选出每年云量小于9% 的影像。2003年后Landsat 7 ETM+系列出现机载扫描行校正器故障，导致此后获取的影像出现数据条带缺失问题，故本文在2003年后综合选取Landsat 5和Landsat 8的遥感影像。其中，2010年因Landsat 8卫星未投入使用，故处理条带后仍使用Landsat 7卫星影像。

辅助数据包括重庆市主城区的矢量数据、SRTM（shuttle radar topography mission）空间分辨率为90 m的DEM数据、源于中国国家统计局的社会经济信息。训练样本通过Google Earth Pro对地物类型目视解译标定得到。

1.2 长时序城镇扩展演变研究方法

本文的技术路线如图 1所示。首先在GEE平台上利用质量评估（quality assessment，QA）波段的云量得分对影像进行云层和云阴影掩膜，并逐像元求取每年所有影像的中值；然后根据监督分类求取每年的土地覆被图；再通过计算求取NDBI，结合土地覆被中建成区和NDBI阈值获得初始的城市建成区提取结果。

1.2.1 城市建成区信息提取

NDBI突出了城市地区，在这些地区，短波红外的反射率通常比近红外区域的反射率高。NDBI负值代表水体，正值较高的代表建成区^[19]。NDBI计算公式如下：

$ B=\frac{b_{\mathrm{SWIR}}-b_{\mathrm{NIR}}}{b_{\mathrm{SWIR}}+b_{\mathrm{NIR}}} $

(1)

式中，b_SWIR为短波红外波段；b_NIR为近红外波段。

1.2.2 ESDA

ESDA适用于检测空间单元间是否存在空间相关性^[20]，全局空间自相关指标能够利用单一值来描述观测范围的自相关程度，也可以用来分析空间数据在局部子区域所表现出来的空间异质性的局部空间自相关^{[21, 22]}。

1）全局空间自相关通常以全局Moran’s I为全局空间自相关指标^[23]，它体现空间邻接或邻近的区域单元属性值的相似水平，计算公式为：

$ \text { Moran's } I=\frac{\sum\limits_{i=1}^n \sum\limits_{j=1}^n C_{i j}\left(x_i-\bar{x}\right)\left(x_j-\bar{x}\right)}{S^2 \sum\limits_{i=1}^n \sum\limits_{j=1}^n C_{i j}} $

(2)

式中，x_i为第i地区的属性值；C_ij代表空间权重矩阵 C第i行第j列的元素；$ \bar{x}=\frac{1}{n} \sum\limits_{i=1}^n x_i$；方差S² =$\frac{1}{n} \sum\limits_{i=1}^n\left(x_i-\bar{x}\right)^2 $。

Moran’s I的取值范围通常为[-1, 1]，趋向-1则代表单元间具有“低-高”“高-低”聚集的属性；趋向1则代表单元间的关系越密切，性质越相似（高值聚集或低值聚集）；趋向0则表示单元间不相关，满足随机分布的要求^[24]。

对于全局Moran’s I，标准化统计量Z (I) 能够较好地检验空间自相关的显著性水平，计算公式为：

$ Z(I)=\frac{(I-E(I))}{\sqrt{V(I)}} $

(3)

式中，V (I) 是Moran’s I的理论方差；E (I)= - 1/(n - 1)，为其理论期望。

2）局部空间自相关是通过计算每个空间单元与邻近单元的相关程度，探究整个区域观测值的变化是否存在空间异质现象。本文运用Moran散点图以及空间关联局域指标（local indications of spatial association，LISA）集聚图对重庆市主城区建成区比例进行局部自相关分析，其中，局部Moran’s I计算公式如下：

$ I_i=\frac{\left(x_i-\bar{x}\right)}{S^2} \sum\limits_j\left(x_j-\bar{x}\right) $

(4)

Moran散点图提供了城市化建设过程中每个格网区域城市建设的局域空间自相关分析。各点横坐标值代表各区域标准化后城市化建设得分，纵坐标值代表标准化后建成区占格网面积百分比的空间位置。

2 结果与分析 2.1 重庆主城区建成区的时间动态变化

20年间重庆市主城区建成区结果见图 2。为凸显其在空间尺度上各阶段的变化，将20年划分为4个阶段进行量化研究。由图 2可知，20年中，重庆市的城市化过程迅速，城市建成区面积逐步扩大，扩张效应明显。以两江交汇所在地的渝中区为城市发展中心，逐渐向外扩散，形成了如今囊括9个城区的重庆市中心主城区的建设现状。

由图 3可知，建成区的面积占比呈增长趋势，从1999年到2018年增长了612.12 km²，2018年的建成区面约为1999年的5.53倍，体现了重庆主城区城市化的快速发展。2009—2014年间，建成区面积的增长速度最快，增长速率达到35.97 km²/a。而2014年后增长速率降至最低，值为28.88 km²/a，经查证，这与2014年颁布的《国家新型城镇化规划》相关，该规划明确要求严格控制城市边界无序扩张，积极盘活可用的存量建设用地。

2.2 重庆市主城区建成区的空间动态变化 2.2.1 城市建成区扩展方向分析

为了反映市重庆主城区城市用地空间扩展在方向上的差异，对于主城区紧凑的城市形态，采用16方位分配法的“等扇分析法”进行探究。16方位包括东（E）、南（S）、西（W）、北（N）、东北（NE）、东南（SE）、西北（NW）、西南（SW）、东北偏北（NNE）、东北偏东（NEE）、东南偏东（SEE）、东南偏南（SSE）、西北偏北（NNW）、西北偏西（NWW）、西南偏西（SWW）、西南偏南（SSW）。以1999年重庆市渝中区的解放碑为圆心，以东西方向为横轴，南北方向为纵轴，北偏东11.25°方向为起始点，将主城区建成区分为16个夹角、面积均相等的扇形区域，并以64.45 km为半径，使之包含所有建成区，如图 4所示。

重庆市主城区的平行岭谷地貌，从根本上限制了早期老城区拓展的空间，形成了南北方向的多个限制带。因此，重庆市受制于山地阻隔效应，建成区沿谷地呈先南后北的蔓延趋势。

在1997年成为直辖市前，重庆市由四川省管辖，经济发展受限^[25]，列为直辖市后GDP不断增长，城市建成区明显扩展。截至2016年底，重庆市主城区建成区的面积与四川成都市辖区相近。

重庆市主城区薄弱的交通基础设施建设也使得组团之间联系不通畅。2001年后，地铁线、高铁线及外环高速的铺设，催生了一些新的组团或产业带，导致区域尺度建成区沿交通线的扩展具有很强的空间集中性。

因此，重庆市主城区建成区向西北扩展至沙坪坝、江北区是受当地自然环境、经济发展状况、交通建设等情况决定的。

2.2.2 城市建成区探索性空间数据分析

根据ESDA方法，考虑研究尺度与集聚效应，本文打破行政区边界，将研究区划分为边长为1.5 km的格网。首先计算全局Moran’s I以及相对应的标准化统计量Z (I)，探究空间邻近的区域单元建成区占格网面积比例的相似度。

从表 1中可以看出，1999—2018年间的Moran’s I值均大于0.75，且在0.82左右浮动，说明城市用地的空间分布呈正相关，且相关程度较大。Z (I)值远大于置信度为0.05时的临界值1.96，表明20年间城市用地的空间分布在总体上存在显著的空间正相关性，呈现出明显的趋同集聚现象。

重庆市城市化过程中，每个区域的城市化建设与周围区域都具有一定的相关性。为了全面揭示区域空间差异的变化趋势规律，采用局部自相关分析方法进行探究，即局部Moran’s I值。利用GeoDa软件^[23]对重庆市每个阶段的建成区占比进行空间自相关分析，得到Moran散点图，见图 5。

由图 5可知：①多数格网区域位于第一、三象限，即城市化较高的区域相互邻接，城市化较低的区域也相互邻接，表明在城市化建设过程中建成区集聚性较强，城市化建设快的区域带动周围区域发展，而城市化建设慢的区域一般处于郊区，周围发展慢，呈现“高-高”“低-低”聚集分布的特点。②第一象限的格网数目远大于第三象限的格网数目，说明城市化快的区域的建成区比城市化慢的区域的建成区更加聚集。③位于第一象限的格网数目逐年增多，说明建成区带动了周围区域的城市化发展，建成区面积逐渐扩大；反之，在远离城市的郊区，城市化发展落后，周围区域不能被其带动发展。

Moran散点图不能反映空间相关类型的显著性水平，需要引入LISA聚类图展示空间要素聚类或分散的地理位置。由于空间异质性是普遍存在的，全局性指标通常难以有效揭示空间内部的局部差异，通过LISA指数可以在局部空间自相关的基础上进一步探究局部空间聚集性。因此，利用GeoDa软件绘制5期LISA聚类图，结果见图 6。由图 6可知，“高-高”格网区域逐年增多，说明这些区域及其周围的城市化发展都很迅速。主城区北部、东部、南部的边缘区域主要是山地和丘陵地貌，不利于城市化发展，只有少数居民在大山里居住，表明20年的经济社会发展对主城区周边山区地带的城市化格局影响不大。在城市群中，出现部分地区处于“低-高”的情况，说明在城市化建设过程中存在时序性和不平衡性，周围区域发展快，但由于城市原因，当地城市化发展慢，形成当地落后于周围城市发展的格局。

3 结束语

基于GEE平台提供的多源、多时相的遥感影像进行探索性空间数据分析，重庆城市化在时间和空间上有以下特点：①1999—2018年，主城区建成区面积增加612.12 km²，扩张速度先快后慢。主城区建成区的扩张主要通过占用植被用地实现，2018年重庆市主城区建成区面积达到747.37 km²，约为1999年的5.53倍。从扩张速度来看，重庆市主城区建成区面积增长经历了先快后慢的过程，其中2009—2014年间，建成区面积增长速度最快，达到35.97 km²/a。2014—2018年间的增长速度放缓，这与2014年发布的《国家新型城镇化规划》密不可分。②主城区建成区沿西北方向的沙坪坝、江北区扩张明显，长江两岸存在极化趋势。渝中区因地势及历史原因，发展程度较周边区县快，以其为中心，因东南部地区地势相对较高，重庆市主城区主要沿NWW方向发展，其中紧邻渝中区的沙坪坝区、江北区扩张趋势明显。以渝中区为主的核心区在城市化建设过程中出现集聚效应，且表现出扩散作用，带动周边地区的发展。而远离城市的远郊、山区城市化发展缓慢，“高-高”格网区域连绵分布。以长江为轴线，长江北岸城区较南岸扩张速度快，未来越过西部的中梁山和东部铜锣山山脉这两条限制带，两侧区域将会成为重庆市主城区的扩展区域。

本文方法的可扩展性较强，基于GEE的城市扩展也适用于其他山地地区。