基于GEE和GIS的黄河三角洲面积多尺度时间序列分析

苏志明,孙永福,2,3*,宋玉鹏,缪栋杰,4,宋丙辉,杜 星

(1.自然资源部 第一海洋研究所,山东 青岛 266061;2.国家深海基地管理中心,山东 青岛 266237;3.青岛海洋科学与技术试点国家实验室 海洋地质过程与环境功能实验室,山东 青岛 266235;4.南京大学 地理与海洋科学学院,江苏 南京 210023)

海岸带是海洋与陆地系统交叉作用、相互影响的地理单元,是研究滨海城市发展的重要平台[1]。河口三角洲地区是典型的河流动力和海洋动力交互控制地带,同时又是多尺度动力结构的耦合区[2]。三角洲岸线长度、形态和所围陆地面积的变化会直观地反映这种河流和海洋对河口三角洲的动力作用,海岸带的合理利用与开发对生态环境的可持续发展有重要的影响,一直以来是研究的热点区域[1-8]。改革开放以来,我国的大陆岸线在自然和人为等因素的影响下发生了巨大的变化。黄河三角洲是中国乃至世界各大河三角洲中海陆变迁最活跃的地区,黄河尾闾不断变迁、不断向海延伸,海岸线蚀退淤进交替演变,三角洲面积也在逐年发生变化。通过对黄河三角洲大陆岸线的变化检测与成因分析,可以为大陆岸线的变迁规律研究提供借鉴和思考,对于沿海生态环境保护、资源开发和可持续发展意义重大。

遥感具有宏观、高效、可重复观测的优点,结合GIS技术能够快速准确地分析大陆岸线长时间序列的变化信息,从而实现海岸线的动态监测[9]。王奎峰等[10]以黄河三角洲1976—2014年的MSS、TM、ETM+影像为数据源,采用一般高潮线法提取不同时期的海岸线,结合RS遥感与GIS技术,对黄河三角洲的岸线变化进行监测。李贺等[11]基于1976—2018年的LandSat长时间序列卫星影像,对黄河三角洲进行了冲淤演变的研究。王一鸣等[12]选取1973-2016年41景影像资料,借助ArcGIS软件得到研究区域面积,并结合利津站水沙资料,利用SPSS软件分析河口地区年造陆面积与入海年输沙量的关系。近年来新兴的遥感大数据处理云平台Google Earth Engine,凭借内置的大量免费数据集和超强的数据处理能力成为遥感分析有力的工具[13-16]。Kilian等[17]基于Python利用图像的监督分类和亚像元分类方法,结合GEE平台编写了CoastSat开源软件,使用户能够从30年来公开的卫星图像中获取全球任何沙质海岸线的位置。在2019年Deng等[18]利用GEE平台研究了长江流域开放水体的动态监测,并且提出了一种新的水体提取规则MIWDR进行长江流域的水体提取。2018年Wang等[19]利用现有的成熟地表水指数方法,提出一种快速估算年际地表水面积的方法,并将GEE平台将所提出的方法应用于1990—2017年长江中游流域地表水的变化分析中。

综上所述,利用软件处理遥感影像已经能够提取并分析黄河三角洲长时间序列的水边线,但容易产生数据质量和主观选取造成的数据差异性,而在大数据平台(GEE)支撑下的遥感水边线提取技术能够通过调用该区域所有的遥感影像进行批量处理,在不减少影像幅数的情况下剔除掉影响水边线提取的噪音像素,在所有的良好观测像素中采用统一的水边线提取规则进行水体提取,可以有效减少这种差异性并且提高处理的速度。因此,本文基于遥感大数据平台Google Earth Engine对1984—2018年35期遥感影像数据进行处理,分析近35 a黄河三角洲面积的变化及其原因。

1 研究区域与研究方法

1.1 研究区概况

黄河三角洲位于黄河入海口处,北临渤海,东靠莱州湾,面积约为5 400 km2,主要位于东营境内,地势西南高、东北低;年均降水量530～630 mm,雨量主要集中在夏季。历史上黄河三角洲分为古代、近代和现代黄河三角洲,古黄河三角洲形成于1855年以前。1855年6月黄河回流入渤海后,尾闾决口50多次,大的改道11次,分别于1855—1934年和1934年至今形成以宁海和渔洼为顶点的近代和现代三角洲[20]。本文中的研究区主要为现代黄河三角洲陆上部分,地理范围为(118°30′～119°24′E,37°30′～38°18′N)。除了黄河河道曾经多次变道,水沙通量也发生显著的变化。历史上黄河年均入海沙量约占全球入海沙量的6%,居世界第二位;含沙量达25 kg/m3,居世界第一位[21]。但是受自然因素和人类活动的共同影响,黄河入海径流量和输沙量自20世纪50年代以来表现出明显减少的趋势,尤其近年来急剧减少[3]。从2002年7月开始,黄河水利委员会实施了黄河调水调沙工程,利用上游的水库调节水沙,从而减轻下游河道淤积和河床升高。黄河三角洲是黄河从上游携带泥沙而沉积形成,所以入海口区域多为粉砂泥质岸线,远离入海口区域多为人工岸线。综合前人研究[22-24]和历史资料,将研究区分为4个岸段:刁口河岸段、人工岸段、河口岸段和莱州湾岸段(图1)。

图1 研究区域 Fig.1 Study area

1.2 数据来源

遥感影像的数据来源于GEE数据库,这些数据分别是LandSat 5的TM传感器、LandSat 7的ETM+传感器和LandSat 8的OLI传感器经过大气校正后的地表反射率数据集,覆盖周期为16 d,空间分辨率为30 m。

径流和泥沙数据采用1984—2018年黄河利津水文站实测的年径流量和年输沙量数据,其中1984—2002年数据来自《东营市水利志》[25],2003—2018年数据来自于《中国河流泥沙公报》[26]。黄河三角洲降水数据来源于国家气象科学数据中心的中国地面气候资料日志数据集(V 3.0)(http:∥data.cma.cn),利用Matlab软件剔除气象数据集中的异常值,按照年份提取出山东省内黄河流经地区的气象站降水数据,包括章丘站、惠民站、济南站和垦利站。

1.3 研究方法

1.3.1 岸线提取

目前提取海岸线的方法主要有人机互译法[27-29]、同月同潮位法、平均低潮位法、一般高潮线法[9]、平均高潮线法[5-6,30]和最大潮高潮线法[7]等,在实际应用中都取得了不错的效果。本文采用平均高潮线法进行岸线提取,取瞬时水边线与一般高潮线的中线作为岸线,即两线所围面积的一半为最终面积。岸线提取流程见图2。

图2 岸线提取流程Fig.2 The process of shoreline extraction

1)瞬时水边线的提取

本文通过修正归一化水体指数(MNDWI)提取瞬时水边线,MNDWI是徐涵秋[31]针对McFeeters[32]提出的归一化水体指数(NDWI)改进的。MNDWI的数学表达式为:

式中,LGreen和LMIR分别为绿波段和中红外波段的像元亮度值。

基于GEE利用MNDWI法提取水体的方法为:①将研究区在每一年内所有Landsat影像进行去云处理,只保留每景影像的有效观测像元;②经过裁剪和镶嵌合成该年份研究区的无云影像;③计算MNDWI水体指数;④利用阈值分割法提取水体,提取面积最大的要素转化为矢量;⑤将数据导入至Arc GIS中,参考Google Earth历史影像进行目视解译。

2)一般高潮线的提取

一般高潮线定义为高潮滩和中潮滩的界线,在提取一般高潮线时参考樊彦国等[5]的提取方法,区别在于处理过程的前半部分在GEE中完成,相比于ENVI等遥感图像处理平台,GEE的图像处理速度显著提升。提取过程为:①在GEE中将经过预处理的影像进行缨帽变换,得到亮度、绿度和湿度三个波段合成的影像;②非监督分类,根据分类的结果在Arc Map中选取训练样本,通过直方图和散点图评估训练样本,使不同分类的颜色尽量不重叠;③监督分类,将影像分为高潮滩、中低潮滩、水体和陆地四类;④分类后在Arc Map中目视解译提取出一般高潮线。

1.3.2 Mann-Kendall趋势及突变检验

在时间序列趋势分析中,最初由Mann和Kendall提出了Mann-Kendall(M-K)非参数检验方法,此方法不需要样本遵循一定的分布,也不受少数异常值的影响,被许多学者用于径流、降水和气温等非正态分布的数据。

对一个序列通过计算得到检验统计量Z,并将其与α显著水平下的Z1-α/2相比较,如果|Z|＞Z1-α/2,就可以认为序列中存在相应的变化趋势。文中取95%显著度,Z1-α/2=1.96。在突变检测中,统计量UFk为时间序列X按照顺序计算的秩序列标准化参数,UFk为标准正态分布,对于给定的显著性水平α,若|UFk|＞U,表明序列存在明显的趋势变化。UBk为按照时间序列X逆序计算的秩序列标准化参数。若统计量UFk大于0,则表明序列呈上升趋势,若小于0则为下降趋势。若UFk和UBk存在交点,且交点位于临界点之间,那么交点对应的时间为发生突变的时刻[33-34]。

1.3.3 交叉小波变换(Cross Wave Transform,XWT)

与传统的傅里叶变换相比,小波变换能够有效地提取信号中的时频特性,有利于分析长时间序列的数据[35]。本文通过小波变换分析黄河三角洲面积和泥沙量径流量以及降水量的周期性变化特征,然后通过交叉小波和小波相干进一步分析数据之间的响应关系。

交叉小波可以对2个时间序列在不同时频域中的相互关系进行分析研究[36-37],设WX(s)、WY(s)分别为给定的两个时间序列X和Y的交叉小波变换,其交叉小波谱为,其中的复共轭,交叉小波功率谱密度为,其值越大表示两个时间序列的相关程度越高[38]。

1.3.4 小波相干(Wavelet Transform Coherence,WTC)

小波相干谱可以反映2个时间序列在时频空间上的局部相关程度[39],定义为:

式中:S为平滑器;s为伸缩尺度;和分别为X、Y的小波变换;为交叉小波谱。对上述小波分析的结果显著性都采用理论红噪声谱进行检验[39]。

2 结果与分析

2.1 黄河入海径流量、输沙量和降水量结果

为了分析黄河三角洲面积和黄河入海水沙量(包含输沙量、径流量和降水量)的关系,首先分析水沙量趋势:对1984—2018年水沙时间序列数据做五年滑动平均处理,并用M-K非参数检验法对水沙量进行趋势检验和突变检测。

1)采用滑动平均法,取k=5绘制出年输沙量、年径流量和年降水量的趋势图:做五年滑动平均后可看出年输沙量(图3)呈显著的下降趋势;做五年滑动平均后可看出年径流量(图4)在1984—2001年呈显著下降趋势、2001—2008年呈上升趋势、2009—2018年总体呈下降趋势;做五年滑动平均后可看出年降水量(图5)总体呈波动上升趋势。黄河三角洲入海口1984—2018年的年输沙量、年径流量和年降水量的M-K趋势检验统计值Z值和斜率见表1。年输沙量和年径流量的Z值与斜率均小于0,二者呈增大的趋势,年径流量的减小速率大于年输沙量,但年径流量|Z|＜1.96,趋势不显著;年降水量的Z值与斜率均＞0,呈增加的趋势,|Z|＜1.96,趋势不显著。

表1 趋势性检验Table 1 Trend test

2)用M-K检验法对年输沙量、年径流量和年降水量进行趋势分析和突变检验、年输沙量的统计量UF(图3b)和年径流量的UF(图4b)整体均为负值,表明二者一直是减少的趋势,年输沙量突变发生在1996年和1998前后、年径流量突变发生在1985和2012年前后,两者突变点的时间存在偏差,这与黄河三角洲的“水沙异源”特点以及一些其他非自然因素的影响有关[40];年降水量的统计量UF和UB(图5b)虽然有交点,但均在±1.96之间,变化趋势不显著,故无突变点。

图3 年输沙量滑动平均图和M-K检验Fig.3 Moving average chart and M-K test of annual sediment discharge

图4 年径流量滑动平均图和M-K检验Fig.4 Moving average chart and M-K test of annual runoff

图5 年降水量滑动平均和M-K检验Fig.5 Moving average chart and M-K test of annual precipitation discharge

2.2 岸线提取结果分析

首先提取出瞬时水边线(图6)和一般高潮线(图7),并基于平均高潮线法计算得出黄河三角洲在1984—2018年的陆地面积,然后对计算结果开展分析。

图6 1984—2018年瞬时水边线变化Fig.6 Changes of instantaneous waterline during 1984—2018

图7 1984—2018年一般高潮线变化Fig.7 Changes of mean high tide line during 1984—2018

2.2.1 拟合结果分析

绘制了黄河三角洲总体(图8a)和各岸段(图8b)近35 a的面积时间序列变化散点图,然后通过非参数方法局部加权回归(LOWESS)[41]拟合。从黄河三角洲总面积(图8a)可知:陆地面积表现出了明显的先增加、后缓慢减少的趋势。1984—1997年为黄河三角洲总面积的增长期,在1993年达到峰值3 443.31 km2,比1984年增加了144.09 km2,河口岸段面积的增长为总面积的增长做出了主要贡献;1997—2018年为黄河三角洲总面积的缓慢减少期,2018年比1993年减少了66.01 km2,减少速率随时间逐渐降低趋于平缓,主要是由河口岸段面积增长速度小于北部刁口河岸段面积的侵蚀速度所致。

图8b明显地反映出黄河三角洲各分区面积的变化趋势,其中刁口河岸段面积呈显著的减少趋势,河口岸段呈显著的增加趋势,1984—1996年大于1996—2018年面积的增长速率;莱州湾岸段和人工岸段没有显著的变化趋势,莱州湾岸段近35 a面积极差为23.11 km2,人工岸段为25.51 km2。

图8 面积的时间序列变化Fig.8 Time series changes of area

2.2.2 河口岸段面积与水沙量交叉小波和小波相干结果分析

利用交叉小波和小波相干方法分析黄河三角洲河口岸段面积和水沙量的时间序列在不同时间尺度的变化和相互关系。交叉小波变换利用2个时间序列的连续小波变换分析共同高能量区和相位关系,小波相干分析2个时间序列在时频空间中局部的相关关系,可以进一步解释在交叉小波能量谱中低能量区的相关性。

图9给出了河口岸段面积与降水量、泥沙量以及径流量的交叉小波和小波相干功率谱图。由河口岸段面积和降水量的交叉小波功率谱(图9a)可知,两者在1989—1991年存在一年的显著的共振周期,通过小波相干图(图9b)可知在2003—2010年存在周期为5～7 a的显著的负相位关系;由河口岸段面积和泥沙量与径流量的交叉小波和小波相干功率谱(图9c～图9f)可知,面积和泥沙量与径流量不存在显著的相关关系。

图9 各因子和河口岸段面积的交叉小波和小波相干功率谱Fig.9 CWT and WTC power spectrum between factors and estuary coastal area

3 讨 论

3.1 黄河三角洲总面积和各岸段面积变化特征分析

黄河三角洲主要是由黄河挟带泥沙冲积而形成的粉砂淤泥质海岸,其面积稳定性较差,受入海水沙、河流改道以及海洋动力条件等因素的影响。从1978—1998年黄河三角洲的面积总体呈现增长的趋势,随后直到2016年起黄河三角洲面积开始呈现减少的趋势[12]。本研究也得到相似的结果:1984—1997年黄河三角洲总面积呈增长趋势,且增长速度逐渐变缓;1998—2018年黄河三角洲面积呈缓慢减少趋势。1976年黄河入海口由原来的北部刁口河河口改为南部清水沟河口,这导致了自1976年后刁口河岸段面积急剧减少、现行河口岸段面积增加,是影响黄河三角洲总面积的主要原因。

本文根据岸段的特点,将整体岸段按照发展特点分为刁口河岸段、河口岸段、莱州湾段和人工岸段四个区域。与牛明香[22]划分的刁口段、东营港及附近岸段、入海口段及莱州湾岸段四个区域相比,虽然岸线提取方法有差异,两者的数据不完全一致,但整体趋势大致相同。

莱州湾岸段在研究时段内没有显著的变化趋势,1984—2007年面积大致稳定,2007年后有轻微的减少趋势,这种变化一方面是由于黄河入海泥沙减少,渤海冷流南下[42]输运到莱州湾岸段的泥沙量相应减少;另一方面由于防潮堤的建设减缓了岸线的扩张。河口岸段面积受黄河入海泥沙影响最显著,在研究时间段内呈显著的增长趋势,按照增长速率可分为1984—1996年和1996—2018年两个时间段。由于胜利油田孤东围堤之间有稳固的防潮堤,所以该人工岸段相对稳定。刁口河岸段从1976年黄河改道后一直呈蚀退趋势,虽然人工固化海岸逐渐增多,但仍然不及海水动力作用的侵蚀速度,故总体面积呈减少趋势。

3.2 水沙量对河口岸段面积的影响

因为水沙量对现行河口岸段的影响最大,故单独选出河口岸段分析其面积的变化与利津站年输沙量、年径流量和黄河三角洲年降水量的关系。通过河口岸段面积的时间序列图(图8b)可知,1984—1996年河口岸段面积呈显著增长的趋势,同时期的输沙量和径流量通过5 a滑动平均(图3a、图4a)亦可得到为1984—2018年的最大值时间段,年降水量的5 a滑动平均也为增长的趋势,说明1984—1996年水沙量是河口面积增长的重要驱动力;1997—2003年河口岸段面积呈缓慢的减少趋势,这也与水沙量有密切的关系。由图3b,图4b和图5b可知年输沙量突变发生在1997年,年输沙量、年径流量和年降水量1997—2002年均为下降的趋势。根据利津水文站的观测数据可知,1996年黄河改道清8汊后,原清水沟流路失去水沙的供应,1997年黄河断流202 d,岸线不断被侵蚀,一直到2004年黄河处于严重的断流状态,导致清8汊流路的造陆面积变小,河口岸段面积整体呈减少趋势。2004—2011年由于调水调沙工程河口岸段面积呈增加趋势;2011—2018年北部沙嘴朝向(图6c和图7c)逐渐由东向变为北向,河口岸段面积呈减少趋势。此时间段输沙量、径流量和降水量与面积的变化关系不显著,可能是因为黄河2002年调水调沙之后,入海水沙模式与之前不同,具体原因尚需进一步探讨。

通过交叉小波和小波相干分析了河口岸段面积与水沙量在不同时间尺度的变化和相互关系,由河口岸段面积与泥沙径流量交叉小波谱和小波相干谱(图9c～图9f)可知,面积与泥沙径流不存在显著的相关关系,这进一步佐证了河口岸段的面积不仅受入海水沙的影响,还受潮汐、海流、泥沙粒径及调水调沙工程等因素的影响,未通过显著性检验的低能量区均以负相位为主,也可以说明河口面积的变化滞后于水沙的变化;由小波交叉谱图(图9a)可知,河口岸段面积与降水量只在1989—1991年存在1 a的显著共振周期,其他未通过显著性检验的低能量区均以负相位为主,由小波相干谱(图9b)可知在2003—2010年存在周期为5～7 a的显著的负相位关系,说明河口岸段面积变化滞后于黄河三角洲年降水量,这与黄河入海泥沙对黄河三角洲造陆存在明显的“后效”影响[22]一致。

4 结 论

基于Landsat卫星数据,利用遥感大数据平台Google Earth Engine和GIS技术分析了黄河三角洲地区岸线变迁趋势,并研究了河口岸段面积与黄河三角洲年输沙量、年径流量、年降水量的关系。得到主要结论:

1)1984—2018年黄河三角洲总面积先增加后缓慢减少,变化主要发生在北部刁口河岸段和现行河口岸段,因人工固化和防潮堤的修建,人工岸段和莱州湾岸段面积保持稳定,对黄河三角洲总面积的影响较小。

2)输沙量和径流量总体趋势基本一致,丰枯交替,整体呈下降趋势。河口岸段面积与水沙量存在密切的关系,并且输沙量、径流量和降水量对面积的变化存在明显的“后效”影响,比面积变化提前一年。此外,海洋中的潮汐和海流等海水动力作用以及人工固岸和调水调沙等非自然因素也对河口岸段的变化起着至关重要的影响。

本文使用的岸线提取方法相对方便快捷,具有一定的普适性,平均高潮线法提取岸线受潮汐影响较小,能够满足宏观分析的精度。LandSat系列卫星的TM、ETM+和OLI数据量大,可以筛选出少云清晰的影像,而GEE平台的高性能计算能力使分析处理数据更高效,特别是在对全球各地岸线的大尺度和长时间序列分析。所以,此方法未来可在其他河口地区进行试验。