格尔木长时间序列遥感生态指数变化监测
2022-08-15汪东川孙志超王海庆王鸿艺
汪东川,陈 星,孙志超,辛 燕,王海庆,柴 华,王鸿艺
1 天津城建大学地质与测绘学院,天津 300384 2 天津城建大学/天津市土木建筑结构防护与加固重点实验室,天津 300384 3 天津绿茵景观生态建设股份有限公司,天津 300384
资源型城市依托资源产业的开发而兴建,由于资源的不可再生性,过度开发在带来社会贡献与经济效益的同时对生态环境造成破坏[1],监测生态环境质量分布与变化对生态环境保护有重要作用。生态系统为人类生产生活提供了重要的物质基础[2],也影响着居民的身心健康、社会的文明进步、人才的吸引和经济发展的长期稳定性[3]。对生态环境进行准确、快速的评价与分析,能为城市环境治理与绿色发展提供可靠的科学依据[4]。
遥感技术的发展使长时序对地观测研究得以实现[5—7],利用多时相遥感数据获取土地利用类型或地表覆盖因子等信息,对地表动态变化特征进行量化表达,快速、准确、多尺度的分析变化规律。在长时间序列遥感的实际应用中,学者常选择3—10年相同间隔年份进行分析[8],时间选择缺少指向性,容易遗漏发生重要变化的年份。有研究通过观察多年数据规律进行人为分段[9]或基于梯度下降法自动分段拟合划分断点位置[10],利用主观判断或数学统计分析方法选取特征时间节点,对长时间序列进行分段。本研究引用Mann-Kendall突变检验法对格尔木RESI均值进行非线性检验,根据突变点对长时间序列数据分段。
目前利用遥感数据研究生态环境多通过地表覆盖类型[8]或植被覆盖指数[11]、不透水面指标[14]、地表温度[12]、地表湿度[13]等生态指标来表示或量化生态环境的优劣情况[15]。然而生态环境是一个综合动态系统,单一因子无法全面反映其变化,因此如何构建综合评价指标是监测生态环境质量研究的重要问题。近年来,有研究者计算多种指数均值[16]或以面积权重法构建综合模型[17],也有利用变异系数法[18]、层次分析法[19]等地学统计方法构建综合模型,但结果往往受人为主观因素影响。2013年徐涵秋[20]综合多方面生态因子提出新型遥感生态指数(Remote Sensing Ecological Index,RSEI),以主成分分析法综合绿度、湿度、热度和干度4个指标,根据各指标对主分量的贡献值确定权重,综合、客观的定量评价生态环境状况。近年来,利用RSEI对不同空间尺度的生态环境进行监测与评价的研究证明RSEI具有可扩展性、可视化和可比性[21—23]。由于RSEI指标选取全面,在长时间序列研究中数据量多,处理工作重复繁琐。谷歌地球引擎(Google Earth Engine,GEE)强大的数据运算功能可以批量、快速的处理影像数据集[24],因此使用GEE云平台预处理并计算各项指标,不仅提高了效率,还可以进行最小云量筛选以获取质量最优的影像进行研究[10,25]。
格尔木拥有丰富的矿产资源,以资源开采为主的人类活动造成生态环境严重破坏。近年来,学者们对格尔木地区的生态环境进行了各方面的研究[26—27],但对格尔木全区综合生态指数的研究较为少见。本文以格尔木为研究区,通过GEE计算1990—2019年遥感生态指数,进一步分析格尔木生态环境时空变化并探讨其影响因素。
1 研究区概况、数据来源
1.1 研究区概况
格尔木地处青藏高原,位于青海省西部,是青藏高原的第三大城市,由柴达木盆地中南部和唐古拉山镇两块互不相连的区域组成,唐古拉山区域内有一级河流长江源头沱沱河(图1)。随着全球气候变化与人类活动的影响,区内土地荒漠化日益严重,生态系统脆弱敏感。受地质构造作用,格尔木拥有众多矿产资源,对矿产资源的粗犷式开采加剧了资源的破坏和环境的恶化[28]。格尔木是新疆、甘肃、四川、西藏、青海五省(自治区)的中心枢纽地带,是青海省两个国家级交通枢纽站之一[29],该地区生态环境可持续发展关系到青藏地区的经济繁荣和社会稳定。
图1 研究区位置
1.2 数据来源及预处理
研究区矢量数据来自中国科学院资源环境科学与数据中心(https://www.resdc.cn/)。遥感数据利用GEE平台编程(JavaScript API)在线调用经过辐射校正、大气校正的Landsat 5/8地表反射数据集,对6—9月影像进行云量筛选,镶嵌合成格尔木近30年(1990—2019年)最小云量影像。通过归一化水体指数(Normalized Difference Water Index,NDWI)阈值法对水体掩膜处理,以减少水域对湿度指标的影响[30]。在GEE平台进行波段运算,得到各年份绿度、干度、热度与湿度四种指数,通过主成分分析得到遥感生态指数(RSEI),再利用ArcGIS与MATLAB等软件进行统计分析,本文技术路线(图2)。
图2 技术路线
2 研究方法
2.1 指数计算
根据徐涵秋提出的遥感生态指数(RSEI)[20],选择归一化差值植被指数(Normalized Difference Vegetation Index,NDVI)表示绿度指标,反映区域内植被生长情况及其分布密度[11]；建筑指数(Index-based Built-up Index,IBI)与土壤指数(Soil Index,SI)求均值得到干度指标(Normalized Difference Soil Index,NDSI),反映土地荒漠化与城镇化光谱响应[31]；单窗算法反演得到地表温度(Land Surface Temperature,LST)表示热度指标[32],反映不透水面热平衡[33]；缨帽变化计算湿度指标(Wet),通过土壤含水量反映土地退化情况[34—35]。各因子计算公式见表1：
表1 指标计算方法
由于各指标在数值单位和大小上有差异,因此要进行标准化处理[20],再通过主成分分析计算初始遥感生态指数(RSEI0)。标准化公式如下：
(1)
式中，N表示标准化后的指标值,I表示该指标在某像元的值,Imin为该指标的最小值,Imax为该指标的最大值。
基于主成分分析法计算初始遥感生态指数(RSEI0)：
RSEI0=1-PC1[f(NDVI,NDSI,LST,Wet)]
(2)
式中，PC1为第一主成分；f是对各指标进行标准化处理。同样为便于指标的度量和比较也要对RSEI0进行归一化处理,得到遥感生态指数(RSEI),值越大表示生态质量越好,反之生态质量越差。
2.2 生态指数动态变化分析
2.2.1Mann-Kendall突变检验
Mann-Kendall突变检验广泛应用于水文序列数据与气象站点数据,是常用的非参数统计方法,不受异常值干扰[36]。区域遥感生态指数均值常用来判断区域生态总体变化趋势,因此本研究对格尔木30年遥感生态指数均值序列进行非参数检验,以获取生态环境突变点对长时间序列数据进行分段研究。
在Mann-Kendall检验中,对长度为n的时间序列R(x1,x2,…,xn)构造秩序列Sk：
(3)
假定时间序列随机独立的情况下,定义标准化统计量UFk：
(4)
式中,E(Sk)与var(Sk)为Sk的期望和方差。
再构建逆时间序列,计算统计量定义为UBk,满足UBk=-UFk。当UFk>0时,序列为上升趋势；UFk<0,序列为下降趋势,UFk、UBk在置信水平上有交点,则该时刻为发生突变时间[37]。
2.2.2变化轨迹法
变化轨迹法利用代码完整的表征时间维度的变化过程[37],有利于深入探讨某一属性动态变化规律,在土地利用变化的研究中有十分广泛的应用[38]。本研究基于遥感生态指数等级,引用变化轨迹法对格尔木RSEI变化轨迹进行可视化。对变化轨迹编码计算公式进行简化：
(5)
式中,Tt为变化轨迹代码,研究节点共有n个,xi为第i个节点对应的RSEI等级代码(1为差、2为较差、3为中等、4为良、5为优)；从左至右,变化轨迹编码首位数字为初始年份RSEI等级,第二位数字为第二节点年份RSEI等级,之后以此类推。
2.2.3Sen+Mann-Kendall趋势分析
Sen+Mann-Kendall显著性趋势检验是长时间序列趋势判断的重要方法[39],由于Sen氏趋势分析本身不能判断序列趋势的显著性,可以通过Mann-Kendall检验Sen氏趋势分析结果实现显著性分析。Sen氏趋势分析通过计算时间序列的中值,能较好的消除异常值的干扰,计算公式如下[40]：
(6)
式中,QRSEI表示时间序列数据趋势的变化；mean为取中值函数；1
Mann-Kendall检验统计量S的计算过程如下：
(7)
(8)
当n≥10时,S近似服从标准正态分布,使用检验统计量Z进行趋势检验：
(9)
其中var表示方差。采用显著性水平α=0.05进行显著性检验,即|Z|≥1.96时发生显著变化,|Z|<1.96时发生不显著变化。
3 结果与分析
3.1 生态环境整体分析
根据格尔木RSEI主成分分析结果(表2)可知：(1)第一主成分(PC1)平均贡献率达到73%以上,说明PC1集中了4个指标的大部分特征信息,能够表示区域内生态环境情况。(2)绿度与湿度的载荷值为正,干度与热度的载荷值为负。符合绿度和湿度越高,地表植被覆盖越高,土壤水分越充足,生态环境质量越好；干度和热度越高,土壤沙化、岩石裸露等问题越严重的自然特征[41]。(3)对比PC1载荷值的变化,LST的绝对值远大于其他指标,且NDSI与LST载荷值的绝对值之和始终大于NDVI与Wet之和,说明格尔木土地荒漠化与不透水面热平衡对生态环境的破坏作用大于植被和土壤湿度的优化作用。综上所述,RSEI指数可以综合反应格尔木生境质量。
表2 主成分分析结果(第一主成分的载荷和贡献率)
图3为1990—2019年各指标及RSEI均值变化情况,RSEI在0.31与0.43之间波动变化,RSEI由0.38下降至0.33,整体降低13%。从单一指标来看绿度变化保持平稳,湿度明显下降30%,干度下降17%,而热度上升12%。综合可知各指标对格尔木RSEI的影响呈下降趋势,这与格尔木市RSEI变化趋势相一致。
图3 生态指标及遥感生态指数均值变化
3.2 格尔木市遥感生态指数突变分析
对1990—2019年的RSEI均值进行Mann-Kendall突变检验,结果如图4。UFk曲线从2001年之后持续低于0,表明格尔木生态质量从2001年开始有变差趋势；UFk曲线与UBk曲线在2006年存在有效突变点,且统计量小于0；UFk曲线在2015年之后超出显著性水平。这与年均RSEI变化趋势基本相同,因此根据突变分析的结果,本研究决定以1990、2001、2006、2015、2019年为时间节点,对格尔木长时间序列数据进行分段。
图4 RSEI突变检验(1990—2019)
对RSEI进行等级划分,并统计各生态等级占比(图5)。由图可知：等级为优的区域所占比例较为稳定,且主要集中在区内冰川积雪地带；等级为良好的部分占比较少,从1990年的2.60%降至2019年的1.50%,主要位于唐古拉山区域,从均匀分布到只集中在唐古拉山地区西北部,中部的分布明显减少；中等级的区域主要分布在柴达木盆地南部山区和唐古拉山区域,在2001年大幅度减少后占比基本稳定；格尔木RSEI等级以较差为主,分别达到48.96%、46.91%、49.07%、59.04%和45.72%；而生态环境等级差的区域主要分布在柴达木盆地北部,1990年仅占12.27%,之后基本稳定在25%上下浮动。
图5 RSEI等级分布(1990—2019)
通过变化轨迹法得到轨迹代码,表示RSEI等级变化过程。RSEI等级未发生变化的区域占格尔木总面积的61.94%,其中稳定不变的轨迹占总面积的35.93%。根据轨迹代码面积,保留发生变化的前15种轨迹代码,占格尔木总面积的30.26%(图6)。其中RSEI等级向好发展的区域仅占总面积的3.55%,主要为较差转为中等,其分布范围在柴达木盆地区域南部和唐古拉山地区东北部。RSEI等级降低的区域占总面积的27.25%,其中2001年发生变化后保持稳定的区域占14.41%,结合等级变化百分比(表3),说明格尔木RSEI等级在2001年下降最为明显；由较差转向差的区域占13.50%,主要分布在柴达木地区北部；由中等转为较差的区域占13.19%,主要分布在唐古拉山地区。根据相关研究,格尔木平均年气温在研究时段内有增长趋势[42],导致冰川融化速度加快,并且降水量上升。短期内河流径流量增加,周边植被覆盖增加,长期发展导致冻土地区植被退化,荒漠化加重,生境质量下降[43]。工业发展对环境的污染也不可忽视,但结合表3结论,RSEI在2015年之后下降幅度有所减缓,说明格尔木城市发展在改善区域生态环境工作中做出了积极的贡献。
表3 RSEI等级变化百分比/%
图6 RSEI等级变化轨迹
3.3 格尔木市遥感生态指数趋势分析
通过Sen+Mann-Kendall趋势分析良好的反映了1990—2019年格尔木RSEI变化趋势的空间分布特征。表4统计可知,格尔木生态环境质量以下降为主,占区域总面积的89.73%,上升区域仅占3.24%,稳定不变即没有发生显著变化的区域为7.03%。
表4 1990—2019年RSEI变化状况统计
根据图7变化趋势分布特点选择典型区域,对其变化原因展开讨论。格尔木RSEI显著上升区域仅有0.29%。近年来国家加快柴达木循环经济试验区建设,对格尔木老工业基地进行高新技术改造,促进产业绿色发展,这一措施减缓了工业发展对地表温度的影响,区域热度指标有所下降。另一方面,三北防护林工程、退耕还林、防沙治沙等项目的加速发展对土地荒漠化治理已有成效,干度指标下降。在城市内部推进林草间作治理、建设农田网格防护林等措施在显著提高区域植被覆盖度的同时也降低了土壤蒸散量[44],使绿度与湿度指标同时上升。因此格尔木市区及乌图美仁乡RSEI主要受人类活动的积极影响显著上升。RSEI不显著上升区域占2.95%。玛曲乡位于可可西里自然保护区,根据李睿等[45]的研究,这一地区气温升高、降水量增加使气候向暖湿化转变,有利于植被生长,区域绿度指标上升明显。而人类活动,如牧业生产和青藏公路、铁路修建等,对该地区土地景观也有一定影响,但整体影响较弱,干度指标微弱下降。在气候变化与人类活动的共同影响下该区域RSEI呈不显著上升。大灶火沟位于东昆仑造山成矿带,主要以黑色及有色金属冶炼产业发展为主,2005年国务院批复建立柴达木地区循环经济试验区[46],促进资源高效、循环利用实现区域可持续发展,有利于热度指标下降,对RSEI有积极影响,使得该区域RSEI不显著上升。RSEI不显著下降区域占43.05%。唐古拉山镇处于长江源头沱沱河流域,气候暖湿化使得植被覆盖逐步改善,减缓草地退化趋势,但绿度指标依然在下降。加之该地区土壤保持功能较差,降雨量增加也加剧了这一区域的土壤侵蚀[47],导致干度指标上升,因此气候变化是唐古拉山镇RSEI不显著下降的主要原因。RSEI显著下降区域占46.68%。近年来,察尔汗盐湖工业区发展迅速,并提出了一系列环保措施,但盐田开采面积增大影响植被生长[48],不可避免地导致绿度下降。柴达木盆地区域近30年降雨量变化不明显,但降水日数略有增加,年平均温度有升高趋势,导致蒸发量增加,土地沙漠化面积增加,整体干度上升[49]。在气候与人类活动共同的负面影响下造成了察尔汗盐湖工业区RSEI显著下降。因气候变暖,格拉丹东冰川地区导致冰川消融增加[50],易引发山洪等地质灾害,植被受到严重破坏,绿度下降。气候变暖导致次生灾害,使该地区RSEI显著下降。
图7 RSEI变化趋势 (1990—2019年)
4 结论
本研究基于长时间序列遥感影像,构建格尔木市1990—2019年遥感生态指数(RSEI)。通过检测突变点对时间序列数据进行分段研究,结合趋势分析对格尔木生态环境进行时空可视化分析,主要结论包括：
(1)1990—2019年间,各数据的第一主成份(PC1)贡献率平均值达到73%以上,绿度与湿度指标对生态环境有积极影响,干度与热度指标有消极影响,与自然环境特征一致,且热度指标对RSEI的影响最大,载荷值关系与RSEI均值变化趋势相一致,说明RSEI可以综合评价格尔木生态环境质量。30年间格尔木RSEI均值从0.38下降至0.33,表明格尔木生态环境整体有退化趋势。
(2)根据突变检验结果,以1990、2001、2006、2015、2019年为节点对格尔木长时间序列数据分段研究。将RSEI划分为5个等级,发现格尔木遥感生态指数等级分布以较差为主,占总面积的45%以上,主要分布在柴达木盆地北部地区。1990—2019年格尔木RSEI等级未发生变化的区域占61.94%,发生变化的轨迹中以等级下降为主,其中2001年下降最为明显,而2015年之后下降幅度有所减缓。
(3)受不同程度的气候与人类活动影响,格尔木遥感生态指数变化以下降为主。显著上升区域集中分布在人类活动密集的格尔木市区及乡镇所在地；不显著上升区域分布于唐古拉山中部玛曲乡及柴达木盆地南部成矿带；不显著下降区域大面积分布在唐古拉山镇以及柴达木盆地南部地区；显著下降区域主要分布在柴达木北部盐湖工业区与唐古拉山南部冰川地带。
