基于谷歌地球引擎的开封城区2010—2019年水体分布变化研究

周 珂, 柳 乐, 程承旗, 苗 茹, 杨 阳

(1.河南大学计算机与信息工程学院, 开封 475004； 2.北京大学工学院, 北京 100871; 3.河南大学, 河南省大数据分析与处理重点实验室, 开封 475004; 4.河南大学， 空间信息处理工程实验室, 开封 475004)

水是生命之源,任何城市的发展都离不开水域的变化。纵观历史,开封市的发展与水系工程建设息息相关。开封市作为八大古都之一、中国31个重点防洪城市之一、中原经济区以及郑州大都市圈核心城市之一,城市发展受到文物保护以及生态保护等多种要求。目前随着遥感技术的发展,遥感技术以其不受地理空间限制、覆盖面广、获得数据速度快等优点[1]已成为监测水体、检测水质以及研究水体变化的重要手段[2-4],通过遥感影像去分析开封市城区水体分布的年际变化趋势,对开封未来的城市发展以及生态环境的保护具有重要意义。

在遥感技术中常用的提取城市间水体的方法有单波段阈值法、多波段谱间关系法、水体指数法等方法[5]。随着遥感技术的不断发展,机器学习分类方法以其自动化程度高,分类精度较准确以及效率高等优点在遥感科学领域取得广泛应用。目前在大数据背景下海量遥感数据集的快速处理变得尤为重要,所以利用高性能云计算平台使用机器学习分类方法对城市水体信息的提取研究很有必要。

谷歌地球引擎(google earth engine,GEE)是谷歌旗下的1款产品,是1个行星尺度的地球科学数据与分析平台,具有强大的计算能力,将数PB(数据存储单位)级的卫星图像和地理空间数据集与行星尺度分析功能结合在一起,使科学家、研究人员和开发人员可以使用它检测地球表面的变化、绘制趋势图并量化差异[6]。在数据方面,GEE云平台在云端存储大量多源影像数据供用户使用,包括Landsat系列、MODIS系列以及Sentinel等常用的影像数据集,且用户可以上传自己的数据在GEE进行计算分析等操作。在语言方面,GEE提供了JavaScript API和Python API, 只需了解JavaScript和python语言的基本语法和数据类型即可上手操作。目前,有许多科学工作者运用GEE云平台进行了大量的地理科学研究工作。Luo等[7]研究了GEE在考古和文化遗产方面的应用；Liu等[8]使用GEE结合Landsat系列影像数据研究了越南玉显县1993—2017年的红树林年际变化检测,得出结论:红树林和红树林-虾塘面积呈现出负相关关系；裴杰等[9]使用GEE利用近2 000景30 m分辨率Landsat-NDVI长时间序列数据研究植被覆盖度变化长时间序列遥感监测。以上研究使用GEE云平台都取得了较好的效果。

结合以上分析,利用GEE高性能云计算平台,结合2010—2019年的Landsat系列影像数据，利用机器学习的方法对开封市城区的水体分布变化进行研究。

1 研究区概况和数据来源

1.1 研究区域概况

选取河南省开封市城区作为研究区(图1)。开封市位于113°5′15″～115°15′42″ E、34°11′45″～35°01′20″N,地处河南省东部。开封市与水有着千丝万缕的联系,境内河流湖泊众多,水资源丰富,分属黄河和淮河2大水系,虽属北城,却拥有“一城宋韵半城水”的赞誉。

图1 研究区Fig.1 Research area

1.2 数据来源及数据处理

1.2.1 Landsat系列影像数据

开封地处华北平原,属温带季风气候,四季分明,4—8月份的影像数据含云量较少,能够较清晰地显示各地物的地貌特征。为了更好地研究开封市城区的水体分布变化,选取的影像为4—8月份云量极少(均小于等于10%)的影像数据。根据研究对象的空间特性,研究所使用的影像数据来自Landsat系列的Landsat 5、Landsat 7以及Landsat 8影像数据。在2010—2011年选取Landsat 5的TM影像作为数据源；在2012年选取Landsat 7 ETM+影像作为数据源；在2013—2019年选取Landsat 8 OLI影像作为数据源。所有操作都在GEE上通过编写代码实现,所获得数据结果存储在Google云盘中。Landsat 5 TM影像和Landsat 8 OLI影像各波段信息如表1、表2所示。

表1 Landsat 5 TM 各波段信息Table 1 Information of each band of Landsat 5 TM

表2 Landsat 8 OLI 各波段信息Table 2 Information of each band of Landsat 8 OLI

1.2.2 样本数据

样本数据均从Google Earth上选取。Google Earth由谷歌公司开发的一款虚拟地球平台,可以为用户提供免费高度清晰的遥感影像[10],用户可以使用Google Earth查看全球范围内近10年来的高分辨率影像,很适合用来选取长时间序列地物分类的样本。根据各地物在卫星影像上的可见光反射率不同以及纹理信息以目译的方法进行选取样本,同时结合外业工作,在2010—2019年共选取1 227份样本数据,分布均匀,涵盖了研究区的所有地物类型,用于训练数据以及验证分类结果。选取的所有样本保存为kml格式，通过Arcgis软件转换成矢量文件上传到GEE云平台。样本选取依据如表3所示。

表3 样本选取依据Table 3 Sample selection basis

1.2.3 数据预处理

(1)影像数据处理。从GEE中下载的Landsat 影像数据均为一级校正的TOA遥感影像数据,在GEE中选择影像的过程如图2所示。选取的2010—2019年影像信息如表4所示,选取的影像数据均为4—8月份的最佳影像数据。其中2012年遥感影像使用的是Landsat 7 ETM+影像数据,据USGS官方公布,自2003年6月以来,由于扫描线矫正器(SCL)出故障导致Landsat 7影像数据出现空白区域,下载的影像缺失区域呈条带状，如图3(a)所示,因此借助ENVI的Landsat_gapfill插件使用三角插值方法对2012年的影像数据进行填充修复,修复后的影像能够对缺失的条带空白地区进行填充,清晰展示各地物的地貌特征,如图3(b)所示。

表4 2010—2019年影像信息Table 4 2010—2019 Image information

图2 GEE中选取影像流程Fig.2 Select image flow chart in GEE

图3 2012年影像修复前后对比Fig.3 Comparison before and after image restoration in 2012

(2)样本数据处理。对上传到GEE的样本数据添加土地覆盖标签landcover属性:水体为0,不透水面为1,植被为2,裸地为3。

2 研究方法

2.1 工作流程

工作流程如图4所示。第1步,在GEE云平台中获取Landsat系列最佳时相数据,并进行数据预处理。第2步,在Google Earth中制作训练样本以及验证样本。第3步,进行光谱特征构建。第4步,使用分类回归树(classification and regression tree,CART)、随机森林(random forest)、支持向量机(Support vector machine,SVM)3种分类器对样本数据集进行训练并对训练后的结果进行精度验证。最后,根据分类结果对开封市城区10年间水体分布变化情况进行研究。

图4 工作流程Fig.4 Work flow chart

2.2 光谱特征构建

在GEE云平台计算经过数据预处理后得到的研究区2010—2019年每景影像的归一化水指数(normalized difference water index,NDWI)[11]、改进的归一化水指数(modified normalized difference water index,MNDWI)[12]、归一化植被指数(normalized vegetation index,NDVI)、归一化建筑指数(normalized difference building index, NDBI)。将计算出的特征指数作为独立的光谱波段添加进原始影像中作为分类特征。4种特征指数的计算方法如表5所示,RRed为传感器红光波段反射率；RGreen为绿光波段反射率；RNIR为近红外波段反射率；RSWIR为短波红外波段反射率。

表5 各指数计算公式Table 5 Calculation formula of each index

2.3 水体提取方法

利用机器学习分类的方法对开封市城区的水体进行提取,其中CART分类、SVM分类、RF分类在遥感领域具有广泛的应用[13-15]。因此使用CART、SVM、RF 3种分类方法在GEE中以样本数据的landcover属性为标签在图像上叠加点以进行训练, 使用与训练相同的光谱特征对影像进行地物分类进而进行水体信息的定量提取,然后对3种分类器的提取水体的精度进行评价,选择最优的分类结果统计开封市水体面积及其他地物的面积，并对水体的提取结果结合空间分布进行分析。

2.3.1 CART分类

分类与回归树(CART)由Breiman等[16]提出。CART是应用广泛的决策树学习方法,由特征选择、树的生成和剪枝组成,既可以用来分类也可以用来回归。CART分类时使用基尼系数(Gini)来选择最优的数据分割阈值划分决策树。假设样本数据共有n个类别，Gini的计算公式为

(1)

式(1)中：ck表示在数据集D中第k类的个数。

基于GEE云平台使用CART分类对样本数据根据选择的光谱特征进行训练自动选择合适的阈值构建决策树对地物进行分类,提取水体信息。

2.3.2 支持向量机(SVM)

支持向量机(SVM)由Cortes等[17]提出。支持向量机算法优化的目标是寻找最大化的边界[18-19]。边界定义为分离超平面与其最近的训练样本之间的距离。在GEE中对SVM分类选用径向基函数核(RBF kernel)[20]以及对应的其他参数使用样本数据对影像进行分类,实现对开封市城区水体分布信息的定量提取,径向基函数核定义为

(2)

2.3.3 随机森林(RF)

RF算法由Breiman[21]提出，是1种组合了多棵决策树的技术。因为该方法具有极好的准确率，能够评估各个特征在分类问题上的重要性以及能提供快速、可靠的分类结果,所以在影像分类上具有广泛的应用。在进行RF分类时参数的设置非常重要,几个重要的参数有树节点数n、最大叶节点数m、分割节点的特征数量v。在GEE中为了更好地定量提取水体信息,在进行RF分类时经过多次实验,将n设为50,将m设为无限制,将v设为0(输入特征总数的平方根)。

3 结果与分析

3.1 数据处理效率分析

与传统的处理遥感影像工具ENVI、ArcGIS相比,GEE集成多种遥感影像数据,对遥感影像进行分析只需编写少量代码就能实现,且GEE基本不占用本地资源。以2018年的影像处理为例,在GEE上加载2018年覆盖研究区域的影像进行裁剪到可视化需要5 s左右的时间,从影像获取到影像分类处理完毕并计算出结果的整个流程所需时间在15 s左右；而在本地环境(CPU:i5-4200H,内存:8 GB,带宽:100 MB)下,所用的2018年4月8日单景影像数据LC81230362018098LGN00,大小为914 MB,理论情况下,在百兆带宽的网络环境下下载时间为71.41 s左右,使用ENVI 5.3对影像进行裁剪、辐射校正、大气校正花费的时间在96 s左右(排除掉手动操作的时间),再进行后续分类操作需要的时间远比使用GEE处理的时间要长。使用GEE能够极大地提高运行效率以及计算能力,效率改善达到10个量级以上。

3.2 精度评价与分析

根据研究区域范围以及样本数量结合算法本身特征将样本数据按照7∶3的比例进行分配,从样本数据中随机选取70%的样本数据作为训练数据集进行分类的训练,剩下30%的样本数据作为验证数据集用于验证分类结果的精确度。由于用户精度、总体精度以及Kappa系数能够很好地反映地物的分类精度,所以使用基于混淆矩阵计算所得的水体提取精度、总体精度以及Kappa系数作为评价指标。

3种分类方法对开封市城区水体分类精度如图5所示。使用CART分类器提取水体的平均精度为96.4%,最小精度为93.3%；使用RF分类器提取水体的平均精度为97.1%,最小精度为94.6%；使用SVM分类提取水体的平均精度为95.9%,最小精度为92.1%。3种分类方法对开封市地物分类的总体精度(OA)和Kappa系数如图6、图7所示,使用CART分类方法的平均总体精度为96.6%,平均Kappa系数为0.950；使用RF分类方法的平均总体精度为96.9%,平均Kappa系数为0.954；使用SVM分类方法的平均总体精度为95.1%,平均Kappa系数为0.930。经过对比使用RF分类能够更好地提取开封市城区的水体信息。

图5 水体分类精度Fig.5 Classification accuracy of water body

图6 3种分类器的总体分类精度Fig.6 Overall classification accuracy of three classifiers

图7 3种分类器kappa系数Fig.7 Kappa coefficients of three classifiers

3.3 水体提取结果与分析

使用水体提取精度最好的RF分类方法分析开封市城区2010—2019年的水体分布变化情况。

如图8所示,近10年研究区的景观格局发生了较大变化,开封市城区的水体主要分布在北部以及中部,开封市西北高东南低的地势造就了河流自西北向东南流的格局。通过表6的面积变化和图9的折线走势得出，开封市城区水体覆盖率先减少后升高,总体上是增加的趋势；植被与裸地总体上趋于减少；不透水面自2016年起较之前年份增长较多,总体上不透水面区域在不断增多。

图8 2010—2019研究区区各类土地分类结果Fig.8 Classification results of various types of land in the study area from 2010 to 2019

图9 2010—2019年研究区各类土地覆盖率Fig.9 Various land cover rates in the study area from 2010 to 2019

表6 2010—2019研究区各类土地面积Table 6 Area of various types of land in the study area from 2010 to 2019

2010—2011年水体和裸地面积在减少,植被和城镇的面积在上升,北部的裸地部分被植被覆盖,中部的裸地得到利用变成了不透水面,西北部的黄河流域由于水位下降部分转变成裸地,水体面积从28.17 km2下降到23.75 km2；2011—2012年黄河水位恢复,中部和东部水体区域增多,总体水体覆盖率基本恢复到2010年的水平,为26.65 km2;2012—2013年植被面积在增加,其他3类地物的面积均在减少,结合李肃秋等[22]的研究数据2013年平均降水量达到历年的最低点335.3 mm，分析此年受降水量影响导致水体面积降至最低,为17.75 km2；2014年降水量增加水体面积有所恢复；2014—2016年,开封市政府为了更好地将水系连贯起来提升城市形象以及改善居民生活修建完成开封西湖,中部的部分裸地变成了水体区域,水体面积增加。2016—2018年,受黄河水体面积影响以及城市在往西扩建发展使得西部地区的水体区域在不断减少。2018—2019年，城市在不断发展不透水面的覆盖率在提高，由44.93%上升到47.13%,城市的快速发展同时，环境以及人文生活都在发生改变。为了改善环境以及能更好地体现开封的“宋文化”,开封市通过“一渠六河”“十湖连通”等工程对城区的水系进行修整,水体面积增加从18.60 km2上升到31.16 km2。

总体上,随着城市的发展,城市水体面积也在不断发生变化,影响水体面积有2个主要原因,其一是由于每年的降水量不同导致研究区北部的黄河流域的水量每年在不断变化,其二开封市古城保护和城市双修建设,水体面积也在逐步增加。整体上,除掉黄河的影响,开封市的水体变化呈现好的发展态势,这是郑汴一体化期间开封市政府湖泊保护政策取得的显著成果。

4 结论

研究了开封市城区的年际水体动态变化,探索GEE云平台的年际水体变化监测能力,基于GEE云平台选取2010—2019年4—8月份的极少云(云量小于10%)Landsat系列影像数据,使用CART、SVM以及随机森林3种分类器提取开封市城区的水体分布信息,结论如下。

(1)GEE在云端存储海量数据,在进行科学研究时不需要在本地下载以及处理,只需将云端处理后的结果下载到本地即可,且其使用谷歌公司强大的计算能力,能够对开封市城区的年际水体信息的提取提供有力保障。

(2)使用CART、RF、SVM 3种分类器都能够对开封市城区年际水体信息进行有效提取,经过对比使用RF分类提取的精度最高,其提取的水体平均精度为97.1%,提取水体的最小精度为94.6%,平均总体分类精度为96.9%,平均Kappa系数为0.954。

(3)开封市城区水体的面积及空间处于不断变化中,水体面积总体趋势在增长,水体区域主要分布在北部和中西部。开封市近年来大力推进“郑汴一体化”的步伐,向西与郑州相向发展,水体区域在不断规划完善。总体上导致水体面积不断变化，主要因素既有北部黄河的水含量在不断变化,也有开封城市化脚步的加快使城区的景观格局发生改变。

由于使用的影像数据为中高分辨率影像数据,不可避免地出现了混合像元的现象致地物分类受到影响,且影响细微水体区域的提取。下一步研究方向将使用高分辨率卫星的影像数据进行水体的提取，将结合深度学习方法研究更加高效的水体提取算法,深入探索城市水体分布变化与城市生态环境变化之间的关系。