高亚洲主要包括青藏高原和周围的高山系(天山、帕米尔高原、兴都库什和喜马拉雅等,图1),有“亚洲水塔”之称。高亚洲冰川覆盖面积约为1.3×105km2,占全球山地冰川面积的43%[1]。在青藏高原地区,冻土面积达1.05×106km2[2],湖泊覆盖面积达36889km2[3],湿地面积为1.33×105km2(其中43.5%分布在长江、黄河和澜沧江-湄公河地区)[4]。高亚洲地区冰雪融化水和降雨会补给周边地区,根据水文平衡方程,高亚洲地区的径流既受降水的影响,也受冰川、积雪、土壤水、湖泊水、多年冻土和地下水等质量变化(或称物质平衡)的影响,同时还受蒸发作用的影响[5]。本文将分析气候驱动因素对高亚洲及其邻区2003～2017年质量变化的影响。

Yao等[6]将高亚洲的7个冰川区域划分为3个气候断面(图1、2):1)沿喜马拉雅中部-羌塘高原-青海高原东部-祁连山,印度季风的影响减弱;2)沿藏东南-羌塘高原-帕米尔高原,印度季风减弱,但西风带增强;3)沿藏东南-喜马拉雅,季风影响减弱。根据Yao等[6]的研究,高亚洲东南部和南部地区的冰川质量平衡受到印度季风和东亚季风环流的影响,在夏季带来丰沛的降雨,而青藏高原西部和帕米尔地区则受西风带的影响而干燥少雨。青藏高原内陆地区冰川物质平衡受环流影响较小,属于大陆性气候,东亚季风只出现在青藏高原东部(如祁连山东部)。巨大的地形地貌差异也是影响高亚洲地区质量平衡变化的重要因素。青藏高原及其周边地区冰川质量平衡有系统区域性差异:在喜马拉雅山脉地区,印度季风的减弱也减少了降水,使冰川质量平衡呈减少趋势;在帕米尔高原东部地区,西风带的增强会增加降水,从而产生一个弱的正质量平衡;在高亚洲中部,这2种环流的影响大大减弱。

图1高亚洲及邻近地区地形、大气环流(箭头)和冰川分布(白点)[6-7]

GRACE卫星可以观测每月的重力变化,并用来研究冰川和其他水文成分的质量变化[8,9,10]。Yi等[11]比较用GRACE得到的质量变化与季风指数认为,印度季风影响略强于西风带,印度季风也会影响帕米尔高原冬季和夏季的降水。Zhan等[7]利用2003～2015年GRACE月重力变化估计的质量变化进行复主成分分析(CPCA),揭示了前3个主成分对应的气候因子对质量变化的贡献:印度季风贡献最大,主要影响喜马拉雅东部、喀喇昆仑、帕米尔高原和印度西北部的物质平衡;厄尔尼诺贡献次之,主要影响青藏高原东部、喜马拉雅东部-羌塘高原-帕米尔和印度西北部;西风带和拉尼娜则主要影响青藏高原西部和西北部。Zhan等[7]认为,厄尔尼诺是第2主成分的驱动因素,与以往多数研究认为应是西风带不一致[6]。Zhan等[7]直接比较CPCA分析得到的时间模态时间系列与季风指数,这种比较值得商榷,因为时间模态时间系列是基于GRACE估计的质量平衡月系列的分析结果,其反映的是累计的质量变化,而印度季风指数和西风指数则可反映不同月份降水导致的质量变化,时间模态时间系列与季风指数没有直接的可比性。此外,Zhan等使用球谐表达的质量变化估计[12]认为,存在高阶截断和高斯平滑等导致的泄露误差影响[13]。因此,必须重新审视主要气候因素对高亚洲地区质量平衡的驱动。

本研究使用最新发布的GRACERL06质量模型数据集[14],直接得到研究区该卫星最长时间跨度(2003～2017年)的质量变化,有效减少了信号泄露的影响。对研究区时空质量变化数据进行复经验正交函数(CEOF,与CPCA相同)分析,得出3个主成分对应的气候驱动因素。此外,将厄尔尼诺和南方涛动(ENSO)同时考虑,新结果有助于解释Zhan等[7]的结果与之前研究不同的原因。

1、数据和方法

1.1GRACE等效水柱高数据

如图1所示,研究区为24°～45°N、65°～105°E的高亚洲及其邻区,比Zhan等[7]研究范围更大,可更全面地了解高亚洲质量变化背后的气候驱动因素。使用NASA喷气推进实验室提供的GRACERL06质量模型数据[14],时间跨度为2003-01～2017-06,12个缺失月份的数据通过线性插值得到,质量变化在0.5°×0.5°的3360个网格上以等效水柱高(EWH)给出。值得注意的是,0.5°网格大小并不代表质量模型的实际空间分辨率,实际分辨率应该是3°[14,15]。使用Peltier等[16]的ICE-6G_D模型,该质量模型已进行冰川均衡调整的改正。

1.2气候指数数据

相关气候指数能较好地反映印度季风、西风带和ENSO的爆发时间、强度以及引起的质量变化。印度季风指数定义为2个区域(40°～80°E、5°～15°N和70°～90°E、20°～30°N)850hPa高度的纬向切变差,该指数的下载地址为http://apdrc.soest.hawaii.edu/projects/monsoon/daily-data.html。西风指数定义为研究区经度范围内2个不同纬度(24°N和45°N)500hPa高度差的平均值,该指数的下载地址为https://www.esrl.noaa.gov/psd/data/gridded/index.html。ENSO指数时间序列是利用经验正交函数方法分析热带太平洋盆地区域(30°S～30°N,100°E～70°W)5个典型变量的结果,这5个变量是海平面气压、海表面温度、纬向风、经向风和长波辐射,ENSO指数的下载地址为https://www.esrl.noaa.gov/psd/enso/mei/。

1.3复经验正交函数法(CEOF)

Pearson[17]提出的经验正交函数(EOF)法能有效提取大气时间序列[18]和地球物理数据[19,20]中的主要成分。然而,由于缺乏相对应的相位信息,EOF法只能检测驻波而不能探测到行波。为此,Rasmusson等[21]引入Hilbert复经验正交函数(CEOF)法,CEOF又称为复主成分分析法(CPCA)[7],能提取随时间空间变化的非平稳模式。本研究进行CEOF分析,相关公式详见附录A(http://www.jgg09.com/CN/volumn/current.shtml)。分析过程主要有5个步骤:1)计算xj(t)——研究区域内以等效水柱高表示的任意网格单元j(j=1,2,3,…,3306)的月观测向量;2)计算Xj(t)——任意单元网格j的EWH复时间序列,其中实部和虚部分别由xj(t)及其Hilbert变换给出;3)对每一网格单元,对复观测向量Xj(t)去平均处理,并用标准偏差进行标准化处理;4)计算任意2个网格复观测向量Xj(t)和Xk(t)的相关系数,构成复相关系数矩阵rj,k;5)分解rj,k,得到所有特征值λn和复特征向量ejn(n=1,2,3…,174),并计算复特征向量在原始观测向量Xj(t)的投影,得到相应的主成分Pn(t)。空间模态由复特征向量的实部给出,相位由复特征向量的实部和虚部给出,时间模态由主成分的实部给出。

使用增量时间模态去除时间模态中累计的质量变化,可与印度季风指数、西风指数直接比较,任何月份的增量时间模态等于原始时间模态的该月值减去前月值。此外,还将对增量时间模态结果及其相应的季风指数时间系列进行小波周期振幅谱分析[11],对原始时间模态及其对应的ENSO指数时间系列进行相同的分析。

2、结果

对GRACE月质量数据系列进行常规的信号分析[10,22,23],对每个网格用最小二乘回归拟合线性趋势、周年和半年的变化。图2给出了研究区2003-01～2017-06月质量变化的趋势速率,其主要异常特征为:在喜马拉雅东部和藏东南(A)、喜马拉雅中部(B)、喜马拉雅西部(C)和印度西北部(D)出现大的负信号,表明质量变化呈下降趋势;其次,在羌塘高原(F)、西昆仑(G)和东昆仑(H)出现明显的正信号,表明质量变化呈上升趋势;此外,帕米尔(E)和天山(I)出现负信号,分别反映弱的和显著的质量变化衰减趋势。

图2RACERL06月质量模型给出的研究区

质量变化趋势(以每年等效水柱高表示)

A:藏东南和喜马拉雅东部;B:喜马拉雅中部;C:喜马拉雅西部;D:印度西北部;E:帕米尔;F:羌塘高原;G:西昆仑;H:东昆仑;I:天山

Zhan等[7]给出藏东南、喜马拉雅东部和整个喜马拉雅的质量变化速率分别为-4.6Gt/a、-4.1Gt/a和-14Gt/a;而图2显示,对应地区的结果分别为-7.6Gt/a、-5.5Gt/a和-16.3Gt/a,均偏大。在印度西北部,图2显示质量变化速率为-17.9Gt/a,与Rodell等[24]的-17.7Gt/a以及Yi等[11]的-20.2Gt/a相近,但是大于Zhan等[7]的-13.6Gt/a。这些差异可能归于信号的泄露[7],而本文直接使用GRACE质量模型可较好地减少信号泄露的影响。此外,这些差异也与研究时段不同有关。

从表1看出,前4个特征值对质量变化的贡献率分别为53%、27%、6%和3%,累计贡献率分别为53%、80%、86%和89%。图3～5给出了前3个CEOF主成分的相应结果,它们对质量变化的贡献累计超过85%。

表1前4个主成分特征值及对应的质量贡献率

图3显示了第1主成分对应的CEOF分析结果,分别用矢量箭头的方向和长度表示质量变化的方向顺序和大小(图3(a))。在印度北部和孟加拉国较大的矢量指向东北方向,然后穿过喜马拉雅、藏东南和兴都库什-帕米尔地区。这说明质量变化是沿着矢量方向依次发生的,与印度季风环流的前进方向基本一致。空间模态显示,印度北部、喜马拉雅、藏东南和天山等地出现正信号,而兴都库什-帕米尔地区出现负信号。空间模态信号的大小与相位矢量的长度相近,因为均包含复特征向量的实部。时间模态(图3(b))呈下降趋势,速率为-1.7cm/a。当滞后1个月时,增量时间模态(图3(c))与印度季风指数具有很强的相关性,相关系数达到0.92±0.16,对应于95%的置信区间。从小波分析结果(图3(d))发现,增量时间模态与印度季风都具有相同的周期(1a)。因此,与第1主成分相对应的质量变化很可能是由印度季风引起的,其对研究区的质量变化贡献占53%。

图3质量变化CEOF分析的第1主成分结果及其与2003-01～2015-12印度季风指数的对比

第2主成分较大的相位矢量基本指向东方向(图4(a)),主要分布在印度西北部、喜马拉雅、藏东南、兴都库什-帕米尔和天山地区,表明质量变化向东依次发生,这与西风带前进的特征吻合较好。空间模态(图4(a))表明,在印度西北部、喜马拉雅、藏东南、兴都库什-帕米尔、喀喇昆仑和天山地区有正信号,而在西昆仑和东昆仑出现负信号,这些信号通常对应于较大长度的相位矢量。时间模态(图4(b))的质量变化有3段趋势速率:-3.5cm/a(2003-01～2007-12)、-1.3cm/a(2008-01～2012-12)和1.9cm/a(2013-01～2017-06)。当滞后1个月时,第2主成分的月增量时间模态与西风指数显示出很强的相关性(图4(c)),相关系数达到0.70±0.15,对应于95%的置信区间。小波分析结果表明(图4(d)),第2主成分月增量时间模态与西风指数具有相同的1a周期。因此,与第2主成分相对应的研究区的质量变化应该是由西风带引起的,西风带对研究区总质量变化的影响为27%。

第3主成分较大相位矢量和显著的空间模态信号主要出现在天山、帕米尔、喜马拉雅和印度西北部地区(图5(a))。相位矢量在印度西北部和喜马拉雅指向东北,在天山和帕米尔地区则转为指向西北,说明第3主成分的质量变化按这些方向依次发生。时间模态在2003～2017年显示了较大的年际变化(图5(b)),但不存在长期趋势。当滞后1个月时,第3主成分的时间模态和ENSO指数显示中等相关(图5(c)),相关系数为0.42±0.15,对应于95%的置信区间。分析第3主成分的时间模态和ENSO指数的小波周期(图5(d))发现,两者均存在3a和6a周期。因此,ENSO可被认为是第3主成分对应的研究区内质量变化的驱动因素。第3主成分的空间模态和相位矢量(图5(a))与第1主成分的结果(图3(a))十分相似,尤其在帕米尔、喜马拉雅和印度西北部地区,意味着ENSO可能通过印度季风(即印度季风被ENSO扰动)影响这些地区的质量变化,但仅占整个研究区总质量变化的6%。

图4质量变化CEOF分析的第2主成分结果及其与2003-01～2017-06西风指数的对比

图5质量变化CEOF分析的第3主成分结果及其与2003-01～2017-06ENSO指数的对比

3、讨论

3.1气候因素对整个研究区的影响

本研究使用的数据时间跨度为2003-01～2017-06,比同类研究更长,因此估计得到的冰川及其他水文因素贡献的总陆地水储量结果会不同[7,11,25]。类似地,冰川实地测量[6]、卫星测高[26,27]和光学遥感[28]的结果也可能有差异。例如,Yao等[6]研究了青藏高原冰川在过去30a里的质量变化,发现质量变化速率(每年等效水柱)在藏东南-西藏西部为-400mm/a,西藏北部为-550mm/a,青藏高原东北部为-450mm/a,但本文基于GRACE的估计结果在羌塘高原(F)、西昆仑(G)和东昆仑(H)出现明显的增加速率信号(图2)。

研究区内复杂的质量时空变化基本来自前3个主成分的贡献,其累计贡献率达到86%。前3个主成分相对应的空间模态、相位信息和时间模态与气候因素影响的特征基本一致,印度季风、西风带和ENSO分别占整个研究区内总质量变化的53%、27%和6%。尽管本研究与先前的研究[6,7]时段不完全相同,但都认为印度季风是引起研究区内质量变化最重要的气候因素。西风带为影响研究区内质量变化的第2重要气候因素,比ENSO重要,这支持了先前大多数的研究[6,11],但与Zhan等[7]的结果不一致。与Zhan等[7]的研究相比,本文的研究区域分别向西和向北延伸了4°和5°,因此完整地包括了阿富汗高地和天山地区。将研究区和Zhan等[7]保持一致,再进行CEOF分析,基本可以重现他们的研究结果,但限于篇幅,未具体展示。因此,在本研究中,西风带对研究区域质量变化的重要性高于ENSO。ENSO通过作用于印度季风对帕米尔、喜马拉雅和印度西北部地区的质量变化产生影响,这也证实了前人在青藏高原南部的发现[29,30]。

Zhan等[7]利用原始时间模态与相应季风指数进行比较发现,第1主成分原始时间模态与印度季风指数的滞后相关系数为0.83±0.23。本研究则使用第1主成分月增量时间模态与印度季风指数相比较,由于理论上的可比性,相关性得到显著增强,相关系数高达0.92±0.16(图3(c));对于第2主成分,本文的月增量时间模态与西风指数的相关系数高达0.70±0.15(图4(c))。

3.2气候因素对特定区域的影响

风力强但随时间减弱的印度季风主要作用于印度北部、孟加拉国、喜马拉雅、藏东南、兴都库什-帕米尔和天山地区(图3(a)),意味着越来越少的降雨导致这些地区的质量变化呈下降趋势。这也证实了Yao等[6]的研究结果,即沿藏东南-羌塘高原-帕米尔和沿藏东南-喜马拉雅地区印度季风存在减弱效应。

风力较大但随时间逐渐加强的西风带主要出现在印度西北部、喜马拉雅、藏东南、兴都库什-帕米尔、喀喇昆仑、天山、西昆仑和东昆仑等地区,带来越来越多的降水,使质量变化呈增加的趋势(图4(a))。这也与Yao等[6]认为沿藏东南-羌塘高原-帕米尔地区西风带效应增强的结论一致。

较强的ENSO主要出现在印度西北部、喜马拉雅、天山和帕米尔地区(图5(a)),ENSO的效应是通过作用于印度季风来实现的[29,31]。

一些地区同时受到多个气候驱动因素的影响,本文将印度季风、西风带和ENSO分别缩写为I、W、E,总结这些地区影响质量变化的驱动因素和顺序:印度西北部(I、W、E)、喜马拉雅(I、E、W)、藏东南(I、E、W)、帕米尔(W、I、E)、天山(W、I、E)、孟加拉国(I)、兴都库什(W、I)、西昆仑(W)、东昆仑(W)。除了大气环流(如印度季风和西风带)和海-气相互作用(如ENSO),巨大的地形地貌差异和上升的气温对现存的冰川质量变化[6]和总陆地水储量变化也有重要影响。以印度季风为主要影响因素的区域质量变化一般呈下降趋势,以西风带为主要影响区域的质量变化一般呈上升趋势。在帕米尔和天山地区,虽然受到西风带的影响,但也受印度季风减弱和与之相关的ENSO的影响,夏季气温升高也会影响冰川的融化[10],因此质量变化呈下降趋势。值得注意的是,印度西北部质量下降大部分是由于地下水过度开采造成的[24]。

4、结语

利用2003-01～2017-06的GRACERL06月质量模型,对高亚洲及邻区质量变化进行CEOF分析,研究质量变化背后的气候驱动因素,主要结论如下:

1)前3个主要成分对研究区内质量变化的贡献率分别为53%、27%和6%,分别归因于印度季风、西风带和ENSO,3个主成分与相应气候因素的滞后相关系数依次为0.92±0.16、0.70±0.15和0.42±0.15,对应于95%的置信区间。使用月增量时间模态有效改善了印度季风和西风指数分别与第1主成分、第2主成分的相关性。研究还发现,ENSO通过作用于印度季风对特定区域(帕米尔、喜马拉雅和印度西北部)的质量变化产生影响。

2)西风带是引起研究区在研究时段内质量变化的第2重要气候因素,这支持了之前的许多研究[6]。当研究区不包括阿富汗高地和天山等西风带起重要作用的地区时,则第2重要气候因素是ENSO,第3是西风带。因此讨论气候因子对质量平衡的影响时,需要关注研究区域的差异性。

3)由于减弱的印度季风、与之相关的ENSO作用和增强的西风带的影响,印度西北部、喜马拉雅和藏东南的质量变化呈下降趋势。由于西风带的增强,兴都库什、西昆仑和东昆仑的质量变化呈上升趋势。兴都库什同时也受到减弱的印度季风影响。虽然受到增强的西风带的影响,但帕米尔和天山还受到减弱的印度季风及其施加的ENSO作用,同时还受气温上升的影响,因此质量变化出现下降趋势。

参考文献：

[3]谢高地,鲁春霞,冷允法,等.青藏高原生态资产的价值评估[J].自然资源学报,2003,18(2):189-196.

[4]牛亚菲.青藏高原生态环境问题研究[J].地理科学进展,1999,18(2):3-5.

[5]相龙伟,汪汉胜,贾路路.GRACE监测青藏高原及邻区陆地水储量变化结果的可变性[J].大地测量与地球动力学,2017,37(3):311-318.

[9]朱传东,陆洋,史红岭,等.高亚洲冰川质量变化趋势的卫星重力探测[J].地球物理学报,2015,58(3):793-801.

[22]罗志才,李琼,钟波.利用GRACE时变重力场反演黑河流域水储量变化[J].测绘学报,2012,41(5):676-681.

[23]周旭华,许厚泽,吴斌,等.用GRACE卫星跟踪数据反演地球重力场[J].地球物理学报,2006,49(3):718-723.

汪汉胜,张会同,相龙伟,江利明,沈强,STEFFENHolger.高亚洲及其邻区2003~2017年质量平衡的气候影响分析[J].大地测量与地球动力学,2020,40(12):1211-1218.

基金:国家重点研发计划(2017YFA0603103);国家自然科学基金(41431070,41974009);中国科学院前沿科学重点研究项目(QYZDB-SSW-DQC027,QYZDJ-SSW-DQC042).