应用GPS垂向位移定量分析2011—2020年云南省极端干旱时空特征

杨兴海， 袁林果*， 姜中山， 冯显杰， 汤苗， 李长海

1 西南交通大学地球科学与环境工程学院， 成都 611756 2 云南大学国际河流与生态安全研究院， 昆明 650091

0 引言

干旱是一种渐发性的自然灾害，根据影响干旱的自然因素将其分为气象干旱、农业干旱、水文干旱与社会经济干旱，持续的气象干旱是引发农业干旱、水文干旱与社会经济干旱的主要因素.干旱严重制约我国社会经济的发展，近年来极端干旱事件的爆发给我国造成了巨大的经济损失.分析干旱的严重程度与时空分布特征，对政府制定抗旱与灾后生态环境修复决策具有重要意义.干旱指数是评估干旱特征的重要指标(Zhao et al.，2017a；Palmer，1965)，传统的气象干旱指数通常基于地面观测站和遥感卫星所观测的气候变量、水文变量，如标准化降水指数SPI(Standardized Precipitation Index，McKee et al.，1993)、标准化径流指数(Shukla and Wood，2008)、标准化地下水位指数(Bloomfield and Marchant，2013)、帕尔默干旱指数PDSI(Palmer Drought Severity Index，Palmer，1965)、自校正帕默尔干旱指数SCPDSI(Self-Calibrating Palmer Drought Severity Index，Wells et al.，2004；van der Schrier et al.，2013)、标准化降水蒸散指数SPEI(Standardized Precipitation Evapotranspiration Index，Vicente-Serrano et al.，2010)与植被条件指数(Kogan，1995)等.但传统地面观测站存在空间分布不均、区域覆盖率低等制约，难以反映大尺度干旱特征，遥感卫星能观测到地表浅层(土壤湿度、植被覆盖等)的水文成分，但无法监测更深层次的水文成分.随着大地测量技术的发展，GRACE系列重力卫星(GRACE / GRACE-FO)提供了全球重力场变化的观测数据(Zhou et al.，2018；Duan et al.，2022)，该数据能有效记录全球范围内陆地水储量的变化(Syed et al.，2008；Zhao et al.，2018；Rao and Sun.，2022).大量研究表明，GRACE数据反演的陆地水储量与降水量、土壤水及地下水的时空相关性较好(李琼等，2013；Han et al.，2019)，因此GRACE陆地水储量可作为研究干旱时空分布的补充数据(Ma et al.，2017).为更准确定量评估干旱严重程度，Zhao等(2017b)提出GRACE干旱指数GRACE-DSI(GRACE Drought Severity Index).GRACE-DSI能反映大尺度陆地水储量亏损情况，有更广泛的适用性，补充了传统气象干旱指数的不足(Liu et al.，2020)，同时与气象干旱指数相关性好(冉艳红等，2021；Satish Kumar et al.，2021).但GRACE的空间分辨率为250～300 km，因此GRACE-DSI更适用于评估全球范围的干旱事件.此外，GRACE时间分辨率为1个月，还存在较多数据缺失，这限制了GRACE-DSI监测短期干旱事件.

陆地水的质量变迁会引起地球的弹性响应.中国陆态网应用GPS大地测量技术近实时地高精度监测地表位移(Jin et al.，2019；Zhang et al.，2019)，为中国区域陆地水储量变化的相关研究提供了新的数据来源.大量研究表明，陆态网提供的GPS垂向位移时间序列包含线性趋势信号、阶跃信号、季节性信号及白噪声.其中，季节性信号包含大气、海洋与陆地水等环境负荷的变化、基岩的热膨胀效应(Yan et al.，2009)及共模误差等.Blewitt和Lavallée(2002)与姜卫平等(2018)对GNSS时间序列的线性项、季节项与阶跃项进行建模，并利用最小二乘方法求解时间序列模型的各项参数，使季节性信号得以分离.此外，Yuan等(2018a，b)与Ming等(2017)分别使用主成分分析与独立成分分析(ICA, Independent Component Analysis)时空滤波，来降低陆态网GPS时间序列中的共模误差等有色噪声信号，闫昊明等(2010)对基岩热膨胀效应引起的位移进行建模，使得季节性信号中的环境负荷位移得以分离.在此基础上扣除环境负荷位移中的大气负荷位移与海洋负荷位移后，陆地水负荷位移信号得以恢复(Fu and Freymueller，2012；Fu et al.，2012；Chen et al.，2013；Chew and Small，2014；胡顺强等，2021).由于陆地水负荷与地表形变的关系可通过格林函数(Farrell，1972)或球谐函数(沈迎春等，2017)建立，因此国内外学者基于负荷位移理论、使用GPS垂向位移反演陆地水储量(Argus et al.，2014；Jiang et al.，2021a；何思源等，2018；成帅等，2021).GPS陆地水储量为单天解，相比GRACE陆地水储量，GPS陆地水储量的时间分辨率更高，同时GPS垂向位移对小尺度的水文负荷变化很敏感，因此更适用于研究小尺度的水文干旱事件.有学者基于GPS陆地水储量计算了水储量异常(Argus et al.，2017；Jiang et al.，2021b)，以反映陆地水储量的亏损情况，进而评估干旱事件的时空分布特征.Jiang等(2021a)参考GRACE-DSI提出GPS干旱指数GPS-DSI(GPS Drought Severity Index)，并将其作为分析干旱严重程度的指标.

云南省位于中国版图的西南部(97.31°—106.11°E，21.80°—29.15°N)，地处青藏高原东南侧，属于低纬度地区，海拔落差明显，地势自西北向东南呈阶梯状递减，各阶梯的温差较大.区域内分布着众多水系，以金沙江、澜沧江与怒江为主要代表.同时，该地区属于典型的亚热带高原季风气候，受到南亚季风与东亚季风的影响，具有干季、湿季分明的特点.全区地表水资源丰富、年平均降水量充沛，但降水的时间和空间分布不均，且常出现季风活动和其他气候条件的异常变化，导致旱灾时有发生.例如，在全球变暖、厄尔尼诺与异常反气旋的影响下，2019—2020年云南省发生了严重干旱，期间降水量不足历史平均值的一半，气温达到历史新高(Ding and Gao，2020).此次干旱是近十年来最严重的一次，其干旱覆盖范围广、干旱持续时间长，期间主要河流湖泊水位急剧下降、多数水源干涸，使得饮用水短缺、农产品大幅减产、经济损失严重，也是导致森林火灾频繁发生的重要原因(Du et al.，2021)，因此研究该地区干旱的时空分布特征具有重要意义.此前，已有学者基于GRACE时变重力场数据，分析西南地区干旱事件中陆地水储量的变化(李琼等，2013；Long et al.，2014；Wu et al.，2020)，同时西南地区密集的陆态网GPS站也为分析干旱事件提供了有力的数据基础.Jiang等(2017)基于GPS水汽反演技术获得可降水量，结合GPS垂向位移研究了云南省的干旱事件，结果表明云南省的GPS可降水量、GPS垂向位移与GRACE等效水高GRACE-EWH(GRACE Equivalent Water Height)的时间变化趋势一致，验证了利用GPS研究干旱的可行性.也有学者结合GPS垂向位移与多种气象干旱指数，分析西南地区干旱事件(姚朝龙等，2019)，研究表明该地区的GPS垂向位移对干旱事件有显著的响应、与气象干旱指数一致性较高.Ma等(2021)利用GPS可降水量计算标准化降水换算指数，结合SPEI及综合气象干旱指数，定量评估云南省的干旱事件.然而，GPS可降水量无法反映地面实际的水储量变化，因此在研究水文干旱时具有不确定性.何思源等(2018)与成帅等(2021)分别使用格林函数法与Slepian球谐系数法，利用GPS垂向位移反演云南地区的陆地水储量，但没有评估干旱事件对陆地水储量的影响.Jiang等(2021b)基于FastICA(Fast Independent Component Analysis)反演模式，利用GPS垂向位移分析了云南地区2011—2019年的陆地水储量亏损情况，但没有研究2020年极端干旱事件及定量分析水文干旱特征(严重程度、时空分布).

本文旨在利用GPS垂向位移反演云南省时变陆地水储量，结合GRACE产品和SCPDSI产品，研究云南省2011—2020年干旱的时空分布趋势，并定量分析2019—2020年极端干旱的严重程度及时空分布.首先，利用FastICA反演模式，反演每日的GPS-EWH(GPS Equivalent Water Height).其次，将云南省按照地形地貌特征划分为五个区域，根据GPS-EWH、GRACE-EWH和SCPDSI计算各区域的GPS-DSI、GRACE-DSI与归一化SCPDSI(SCPDSI-Z)，并比较其相关性.最后，根据GPS-EWH计算五个区域的GPS-EWH异常时间序列，评估干旱事件中陆地水的亏损情况，同时结合GPS-DSI、GRACE-DSI与SCPDSI-Z，分析云南省2019—2020年极端干旱事件的时空分布.

1 数据来源及预处理

1.1 GPS时间序列数据

本文采用国家地震科学数据中心(https:∥data.earthquake.cn/)提供的中国陆态网GPS时间序列，其结果由GAMIT/GLOBK软件解算得到，解算时扣除了固体潮、海潮与极潮的影响(Davis，2004).为保证观测数据能够覆盖整个研究区域，采用2011年1月1日至2020年12月31日云南省内及其周边共44个观测站(如图1所示)的位移时间序列.GPS时间序列包含构造运动引起的长期线性趋势项、地球物理因素引起的季节项及仪器更换或地震引起的阶跃项(姜卫平等，2018)，其中季节性信号包含由大气、海洋与陆地水等物质的质量负荷变化引起的位移.由于陆地水负荷引起的垂向位移比水平位移更显著(Argus et al.，2014)，因此本文使用GPS垂向位移时间序列进行研究，后文中的位移时间序列均指垂向位移的时间序列.

为消除GPS位移时间序列中由其他地球物理现象引起的位移信号、改正地震与仪器更换引起的阶跃信号，利用德国地学研究中心GFZ(GeoForschungsZentrum，https:∥www.gfz-potsdam.de/en/home/)提供的非潮汐大气负荷与非潮汐海洋负荷产品(Dill and Dobslaw，2013)扣除GPS位移时间序列中的大气负荷与海洋负荷位移；采用欧洲中期气象预报中心ECMWF(European Centre for Medium-Range Weather Forecasts)提供的第五代全球大气再分析日平均地表温度产品(https:∥cds.climate.copernicus.eu/cdsapp#!/home)，使用Yan等(2009)与闫昊明等(2010)提出的GPS基岩热膨胀模型，改正GPS时间序列中基岩热膨胀效应引起的位移.最后，使用最小二乘法与GPS时间序列模型(姜卫平等，2018)对各项参数进行建模，并扣除趋势项及阶跃项.采用的拟合模型如下：

(1)

其中，y(ti)为ti时刻的位移观测值；y0为测站的标准位置；a0为线性项速度；Ck与Sk为季节项的振幅，k=1时对应周年项，k=2时对应半周年项；gj为阶跃项，H为海维西特阶梯函数，阶跃前为0，阶跃后为1，Tgj阶跃的时间；e(ti)为随机噪声.

数据处理后得到垂向水文负荷位移时间序列，为保证反演精度还需筛选站点数据.移除数据缺失时间超过所研究时间域30%的GPS站数据，同时比较邻近站点的时间序列信号特征，去除含有异常信号的站点数据.筛选过程中发现MMMZ站缺失数据超过5年、YNGM站含有异常信号，因此移除这两个站的数据，得到42个GPS站的位移时间序列.

为保证反演的时变陆地水储量的完整性，期望选择的GPS站在每个历元都有数据，但实际部分历元的数据会存在缺失.因此，本文使用Liu等(2018)开发的GNSS缺失数据插值软件(GMIS)，对数据缺失较少的GPS站进行插值.该软件基于克里金-卡尔曼插值方法(Mardia et al.，1998)，将克里金空间域插值方法与卡尔曼滤波时间域最优估计结合.如图2所示，云南地区所有的GPS站点都有明显的周年振幅，表明该地区的水文负荷季节性变化显著.GPS站的周年振幅平均值为6.9 mm，最大可达10.1 mm，最小为3.3 mm，周年振幅较大的站主要集中在云南西南部.GPS站最大沉降量对应的时间主要集中在9月下旬与10月初，该时间比最大降水量对应的月份(7月与8月)滞后约2个月.

图2 GPS的垂向水文负荷位移周年振幅.箭头方向为最大沉降量对应的时间，颜色图为振幅内插结果Fig.2 Annual amplitude of GPS vertical hydrological loading displacements. The direction of arrows represents the time corresponding to the maximum settlement, and the color plot represents the amplitude interpolation result

1.2 GRACE数据

GRACE地球时变重力场数据可以反映地表质量负荷的时空变化，在陆地水储量监测中得到广泛的应用.美国德克萨斯大学空间研究中心CSR(Center for Space Research，http:∥www2.csr.utexas.edu/grace/)提供了时变重力场数据，其中2017年6月前的数据由GRACE重力卫星观测得到，2018年6月之后的数据由GRACE-FO重力卫星观测所得.为验证GPS陆地水储量结果的可靠性，本文采用2011年1月至2020年12月的GRACE时变重力场Mascon解产品(RL06 v02版本)与之对比.Mascon产品解算时将原始的C20系数替换为卫星激光测距的C20系数(Loomis et al.，2019)，同时还将GRACE-FO的C30系数替换为卫星激光测距的C30系数，基于Richard Peltier等(2018)提出的ICE6G-D模型改正冰后回弹.相较于球谐系数产品而言，Mascon产品不需要进行滤波、平滑或者缩放处理，已被广泛应用于陆地水储量变化相关研究(Save et al.，2016).Mascon解被估计在1°×1°的六边形格网中，以0.25°×0.25°的格网形式提供，但其实际的分辨率由卫星轨道高度决定，约为250～300 km(Save et al.，2016).

1.3 降水及干旱指数

降水是陆地水补给的一个重要途径，直接关系着气象干旱的发生与否.本文采用国家气象科学数据中心(http:∥data.cma.cn/)提供的月平均降水量格网数据，其时间范围与GPS位移时间序列保持一致.该格网数据由气象观测站数据空间内插所得，云南省内共均匀分布109个气象站点，所得格网空间分辨率为0.5°×0.5°.

干旱指数可直观地反映复杂的干旱情况，有利于研究干旱严重程度及时空分布情况.SCPDSI由Wells等(2004)提出，是PDSI的改进方案.在SCPDSI计算过程中将参考作物替换成实际的植被覆盖，并考虑了潜在蒸散量与积雪的动态变化，使干旱的空间特征更突出，因此SCPDSI被广泛应用于气象干旱的相关研究.为研究云南省水文干旱与气象干旱的相关性，本文采用东英吉利大学气候研究组CRU(Climatic Research Unit，https:∥lr1.uea.ac.uk/cru/data)提供的SCPDSI，时间分辨率为1个月，空间分辨率为0.5°×0.5°.

2 研究方法

2.1 时变陆地水储量反演方法

固体地球是弹性体，陆地水的质量变化会引起地球表面的弹性形变，Farrell(1972)使用格林函数确定地表质量变化与形变之间的关系.Argus等(2014)基于圆盘负荷理论，提出利用GPS垂向位移反演陆地水储量，其最小二乘反演模型如下：

‖Ax-b‖2+β2‖L(x)‖2→min,

(2)

其中，A为格林函数系数矩阵，基于Wang等(2012)提出的负荷勒夫数(PREM地球模型)计算得到；x为待估计的EWH向量；b为GPS垂向位移观测向量；β为平滑因子；L为拉普拉斯平滑矩阵.式(2)理论上应考虑反演区域内与区域外所有质量负荷的影响，何思源等(2018)研究表明，基于格林函数反演方法、利用GPS垂向位移反演云南地区陆地水储量变化时，将边界向外扩充3°，场外负荷源的影响可以忽略不计.

Jiang等(2021b)采用FastICA反演模式，反演了云南省的时变陆地水储量.该方法不同于单个历元反演模式，无论观测数据的历元多少，其进行反演的次数都等于独立分量IC(Independent Component)个数，因此能有效地提高反演的时间效率.并且，GPS垂向位移时间序列含有许多噪声信号，FastICA方法可以有效过滤噪声，具有收敛速度快、稳定性强的优点(Hsu et al.，2020；Jiang et al.，2021b；Hyvärinen and Oja，2000).其反演步骤如下：

(1)利用GPS垂向位移时间序列构建观测矩阵，行号为历元的序号，列号为GPS站的序号.

(2)使用FastICA算法将观测矩阵分解为空间函数与时间函数IC，并根据每个IC对滤波后数据的方差贡献进行排序，以确定IC的个数，分解公式为

Xr×p=Vr×nWn×p,

(3)

其中，X为GPS垂向位移观测矩阵；V为时间函数，垂向位移与反演的EWH共用一套时间函数；W为垂向位移空间函数；r为历元个数；n为IC个数；p为站点个数.

(3)根据式(2)将每个垂向位移空间函数IC分别反演为EWH空间函数IC，反演过程中将边界向外扩充3°，并采用交叉验证方法(Jiang et al.，2021b；Matthews and Segall，1993)确定最优的平滑因子为0.012.

(4)各时间函数IC与对应的EWH空间函数IC乘积之和为等效水高总量：

Yr×m=Vr×nUn×m,

(4)

其中，Y为EWH向量；V为时间函数；U为EWH空间函数；m为等效水高格网个数.

本文沿用上述反演模式，使用2个IC来分析，共减少90%的数据方差，IC1与IC2对滤波后数据的方差贡献率(S1,S2)分别为63.4%与36.6%.如图3b—c所示，85.7%的GPS站IC1垂向位移空间函数绝对值大于10 mm，IC2为61.9%.其中，IC2的空间函数在滇西南较大，最大可达29 mm；而IC1则相对分布均匀，稳定在10～20 mm之间.此外，IC1的时间函数在每年第228天附近达到最小值，对应了最大降水月份(7月)；IC2的时间函数则在每年第318天附近达到最小值，滞后于最大降水月份约两个月.两个IC的时间函数都有短期下降或抬升的趋势，可能与水文极端事件(极端干旱、极端降水)有关.

图3 IC1与IC2的时间函数 (a) 与垂向位移空间函数 (b—c)Fig.3 Temporal functions (a) and vertical displacements spatial functions (b—c) of IC1 and IC2

2.2 干旱指数

用干旱指数反映水储量的亏损情况，能直观地评估干旱的发展过程及严重程度(Van Loon，2015).常用的PDSI、SPI等气象干旱指数通常仅考虑了少数水文或气象变量，不能完整地反映陆地水储量的亏损情况，同时各种干旱指数考虑的水文与气象变量不完全相同，因此各种干旱指数存在差异(Shi et al.，2021；Wang et al.，2015).Zhao等(2017a，b)计算的干旱指数GRACE-DSI，可以反映陆地水储量的亏损情况，但受限于时空分辨率低与数据缺失.Jiang等(2021a)基于GPS-EWH计算的干旱指数GPS-DSI，同样能反映陆地水储量的亏损情况，且不存在数据缺失，可以作为评估干旱的一种有效补充手段.GPS-DSI计算方式如下：

(5)

为更准确地评估干旱，本文采用与式(5)相同的标准化方法，将GRACE-EWH和SCPDSI计算为GRACE-DSI与SCPDSI-Z，结合GPS-DSI、GRACE-DSI和SCPDSI-Z分析干旱特征.参考Zhao等(2017b)的GRACE-DSI分类标准，将文本中使用的干旱指数划分为表1中的11个类别.

表1 干旱程度分类Table 1 Classification of drought

3 结果分析

3.1 陆地水储量变化时空分布

本文估计了研究区2011—2020年的GPS-EWH日变化量，图4展示了GPS-EWH与GRACE-EWH的周年振幅及多年平均降水量，三种数据周年变化的空间分布趋势一致，即滇西南的周年变化最大，逐渐向滇东北方向减小.其中滇西南的GPS-EWH周年振幅最大值约为300 mm，GRACE-EWH约为220 mm；滇东北的GPS-EWH周年振幅最小值约为20 mm，GRACE-EWH约为40 mm.GPS-EWH与GRACE-EWH存在差异，其主要原因可能是GRACE的空间分辨率低，对局部区域的陆地水储量变化不够敏感(成帅等，2021；何思源等，2018；Argus et al.，2014)，及GPS时间序列中还含有未建模的非水文负荷位移信号(Fu et al.，2015).图4c中多年平均降水量在滇西南达到最大值，约1750 mm；在滇东北减小到约900 mm；在滇西北达到最小值，约600 mm.多年平均降水量的总体空间分布趋势与GPS-EWH和GRAEC-EWH周年振幅一致，表明该地区陆地水储量变化受降水驱动.此外，滇西北降水量最小，但GPS-EWH的周年振幅较大，说明该地区陆地水储量可能受其他因素(径流、蒸散)影响较大.

本文为研究云南地区的陆地水储量时间变化特征，比较了研究区域内GPS-EWH、GRACE-EWH与降水量的月度时间序列.如图5a所示，GPS-EWH与GRACE-EWH的相关系数为0.82，两者的时间变化趋势具有高度一致性.但GPS-EWH与GRACE-EWH的时间序列周年振幅分别约为127 mm与64 mm，GRACE-EWH比GPS-EWH低估了50.4%的陆地水储量.

GPS-EWH与GRACE-EWH都是陆地水的状态量，而降水属于通量，两者是不同的物理量，陆地水储量与降水的关系可由下式简单表示：

(6)

其中，t为时间，TWSA为陆地水储量，P为降水，R

为地表径流量，E为蒸发量.由式(6)可知陆地水储量由降水、径流与蒸发量累积得到，因此图5中GPS-EWH与GRACE-EWH的相位滞后于降水2个月，在其他相关研究中有相同结论(Jiang et al.，2017；何思源等，2018；成帅等，2021).

本文将水储量等效水高进行微分，以此统一陆地水与降水的物理量，如图5b所示，GPS-EWH微分、GRACE-EWH微分与降水相关性较高、时间变化趋势一致.但三者也存在差异：(1)GRACE-EWH微分比GPS-EWH微分的变化幅度小，其原因可能是GRACE空间分辨率低；(2)GPS-EWH微分与降水的波谷不一致，结合公式(6)分析，可能是干季时陆地水储量主要受蒸发与地表径流影响.

图5 陆地水储量的时间变化.EWH时间序列(a)，EWH微分时间序列(b)，EWH异常时间序列(c).R、R1、R2与R3分别为GPS-EWH与GRACE-EWH、GPS-EWH微分与降水量、GPS-EWH微分与GRACE-EWH微分及GRACE-EWH微分与降水量的相关系数Fig.5 Temporal variations of terrestrial water storage. EWH (a), EWH differentials (b), and EWH anomalies (c). R, R1, R2, and R3 are the correlation coefficients between GPS-EWH and GRACE-EWH, GPS-EWH differential and precipitation, GPS-EWH differential and GRACE-EWH differential, as well as GRACE-EWH differential and precipitation, respectively

3.2 等效水高与干旱指数及其相关性分析

本文按照地形地貌特征将研究区域分为图1中的五个区域(滇东北、滇东南、滇中、滇西北与滇西南)，并探究各区域的GPS-DSI、GRACE-DSI与SCPDSI-Z之间的相关性.如图6a所示，GPS-DSI与GRACE-DSI在滇西南与滇西北地区相关性最高，两者在滇东的相关性最低.对比图6a与图6b可知，在滇中与滇西地区GPS-DSI与SCPDSI-Z的相关性更高.如图6b所示，GPS-DSI与SCPDSI-Z在滇中和滇西南地区相关性最高，相关系数约为0.85，表明在该地区GPS-DSI能作为评估干旱的补充数据.但在滇东地区GPS-DSI与SCPDSI-Z呈现出负相关，对比图2中GPS站点分布，可能是滇东地区GPS站空间分布稀疏导致GPS-DSI不确定性较大.

图7a展示了五个区域GPS-EWH与GRACE-EWH的时间序列及其相关系数.GPS-EWH与GRACE-EWH表现出高度的一致性，两者在滇西北、滇西南与滇中地区的相关系数最高，分别为0.87、0.86与0.79.在滇东南与滇东北GPS-EWH与GRACE-EWH的相关系数分别为0.56与0.42，呈现出弱相关性，同时GPS-EWH与GRACE-EWH的时间变化趋势略有区别.结合图2分析，可能是GPS站空间分布稀疏导致了滇东南与滇东北地区的GPS-EWH具有较大不确定性.

图7b展示了云南省五个区域GPS-DSI、GRACE-DSI与SCPDSI-Z的时间序列及它们间的相关系数R1、R2与R3.滇中、滇西南与滇西北的R2分别为0.72、0.71与0.56，大于R3和R1，表明GPS-DSI对小区域的干旱事件更敏感.此外，五个区域的GRACE-DSI与SCPDSI-Z相关性较好，R3均约为0.5，表明水文干旱与气候干旱的一致性较好.但滇东地区的R2明显小于R3，表明该地区的GPS-DSI可能具有较大不确定性，这与图6b的结果相同.除滇东地区的GPS-DSI外，三种干旱指数在2019年出现大幅度下降，一直持续到2020年，表明期间有严重的干旱事件发生.

图6 GPS-DSI与GRACE-DSI (a) 及GPS-DSI与SCPDSI-Z (b) 的相关系数Fig.6 Correlation coefficients between GPS-DSI and GRACE-DSI (a), as well as GPS-DSI and SCPDSI-Z (b)

图7 云南省五个区域的EWH (a) 与干旱指数 (b) 时间序列.R、R1、R2和R3分别为GPS-EWH与GRACE-EWH、GPS-DSI与GRACE-DSI、GPS-DSI与SCPDSI-Z及GRACE-DSI与SCPDSI-Z的相关系数Fig.7 Time series of EWH (a) and drought severity indices (b) in the five regions of Yunnan Province. R, R1, R2, and R3 are the correlation coefficients between GPS-EWH and GRACE-EWH, GPS-DSI and GRACE-DSI, GPS-DSI and SCPDSI-Z, as well as GRACE-DSI and SCPDSI-Z, respectively

综上可知，三种干旱指数存在一定的差异，其主要原因可能是：(1)滇东地区GPS站点空间分布稀疏，导致估计的GPS-EWH与GPS-DSI具有较大不确定性；(2)GPS时间序列中含有未建模的非水文负荷位移信号；(3)GRACE的空间分辨率较低，因此在小区域内GPS-DSI比GRACE-DSI更敏感；(4)GRACE卫星在任务最后几个月数据质量不佳，导致该时期GRACE-DSI与GPS-DSI、SCPDSI-Z的差异较大；(5)SCPDSI-Z属于气象干旱指数，与水文干旱指数考虑的因素不同，使得它们表现出时间差异、干旱程度差异.因此，结合多种干旱指数，使其优势互补，能更好地评估干旱事件.

3.3 干旱的时空分布特征

图8展示了干旱指数分别在DSI≤-0.8(中度干旱及以上)、DSI≤-1.3(重度干旱及以上)与DSI≤-1.6(极端干旱及以上)干旱区间的频率.如图8a—c所示，云南省大部分地区发生中度及以上干旱事件的频率超过0.2.GPS-DSI、GRACE-DSI与SCPDSI-Z在中度及以上干旱区间且频率超过0.2的格网分别有87%、43%与74%.如图8d—f所示，GPS-DSI与SCPDSI-Z属重度及以上干旱区间且频率大于0.1的格网分别有53%与38%，该类型干旱主要发生在滇中与滇南.图8g—i表明GPS-DSI与SCPDSI-Z在极端干旱区间的频率较低，且极端干旱主要集中在滇西南地区，GRACE-DSI在该种干旱区间的频率最低.

图8 DSI≤-0.8 (a—c)、DSI≤-1.3 (d—f)与DSI≤-1.6 (g—i)的频率空间分布Fig.8 Spatial distributions of frequency of DSI≤-0.8 (a—c), DSI≤-1.3 (d—f) and DSI≤-1.6 (g—i)

图9 云南省五个区域GPS-EWH (a)、GPS-DSI异常与降水异常时间序列(b)Fig.9 Time series of GPS-EWH (a), GPS-DSI anomalies, and precipitation anomalies (b) in the five regions of Yunnan Province

图9b展示了五个区域的降水异常与GPS-DSI时间序列，降水异常用于评价降水量的亏损程度.滇西南、滇西北、滇中的GPS-DSI与降水异常整体变化趋势相似，降水异常显著上升时GPS-DSI会出现上升趋势，降水异常显著下降后GPS-DSI也会下降.但滇东南2018年的GPS-DSI存在异常的波谷，滇东北GPS-DSI与降水异常的变化趋势一致性较差，因此分析这两个地区的水文干旱时需参考GRACE-DSI.

根据图9与图7b中的GPS-DSI与GRACE-DSI时间序列，云南省在2011—2020年间发生了三次明显的水文干旱事件.第一次水文干旱开始于2011年8月份左右，该次水文干旱持续时间短、干旱程度较轻，在2012年4月左右恢复为无旱.第二次开始于2012年10月左右，此次水文干旱程度与第一次相似，并于2013年12月左右恢复为无旱.第三次水文干旱开始于2019年4月左右，持续到2020年10月左右出现缓解，而后又加剧.Ma等(2017)与Jiang等(2017)分别利用GRACE陆地水储量和GPS可降水量分析云南省干旱，均发现第一次与第二次水文干旱.Shi等(2021)通过分析滇中、滇西南与滇西北的GRACE陆地水储量异常，发现了这三次水文干旱事件.在这三次干旱事件中，第三次水文干旱事件在10年内持续时间最长、干旱程度最严重.此次干旱事件中降水异常大多数月份为亏损，虽然有少许月份的降水异常为正值，但降水持续时间短，并没有缓解水文干旱，同时在滇西南降水量亏损的月份最多，其亏损量在2019年8月达到最大值.持续的降水亏损导致GPS-DSI与GRACE-DSI大幅度下降，2020年5月大部分地区都达到异常干旱.

综上所述，云南省时常有中度及以上的干旱发生，偶尔有重度干旱与极端干旱发生.2011—2020年共有3次显著的干旱事件发生，其中极端干旱开始于2019年4月，持续到2020年底.

3.4 2019—2020年极端干旱的空间分布

为进一步评估2019—2020年云南省极端干旱事件的空间分布，图10展示了此次干旱最严重的时期(2019年6—11月、2020年3—8月)的GPS-DSI、GRACE-DSI与SCPDSI-Z的空间分布.如图10所示，三种干旱指数整体的干旱程度相似，2019年云南大部分地区属于中度干旱，少部分地区达到重度干旱，2020年云南大部分地区达到重度干旱与异常干旱，直到2020年8月才出现消退的迹象.

图10 2019年(左3列)与2020年(右3列)GPS-DSI、GRACE-DSI与SCPDSI-Z空间分布Fig.10 Spatial distributions of GPS-DSI, GRACE-DSI and SCPDSI-Z in 2019 (left 3 columns) and 2020 (right 3 columns)

具体分析图10中GPS-DSI空间分布可知，2019年6—7月滇西北与滇中地区的干旱程度最严重，可达到重度干旱与极端干旱.2019年雨季后，8—10月极端干旱消退，11月重度干旱也消退，但大部分地区仍属于中度干旱.2019年冬季与次年春季降水量偏小，陆地水没有得到足够的补充，导致2020年春季干旱再次加重.2020年3月滇中、滇西北和滇西南绝大部分地区发展为极端干旱，4—7月干旱范围继续扩大，使云南大部分地区发展为极端干旱与异常干旱.2020年8月，得到降水补充后，干旱程度有所减轻，滇西北大部分地区降低为中度干旱，仅滇西南存在部分地区属于异常干旱.

由于滇东地区的GPS-DSI不确定性较大，因此研究此区域的干旱还需参考GRACE-DSI与SCPDSI-Z.如图10中的GRACE-DSI所示，2019年6—8月干旱的覆盖面积最广，大部分地区为中度干旱，9—11月滇西北与滇东北逐渐减轻为无旱，而滇西南出现干旱加重趋势，此结论与冉艳红等(2021)的研究结果一致.2020年3月，滇南的GRACE-DSI仍为干旱，4—7月干旱范围扩大、干旱程度加重，大部分地区发展为极端干旱与异常干旱，8月干旱由南向北逐渐减弱.此外，根据SCPDSI-Z可知，2019年6—11月大部分地区都为中度干旱，6—10月小部分地区出现重度干旱或极端干旱.除滇东的GPS-DSI外，2020年3—7月SCPDSI-Z的空间分布与GPS-DSI、GRACE-DSI相似，滇西北地区为无旱，而其他地区干旱覆盖范围广、干旱程度高，直到8月干旱有所缓解.

综上可知，2019—2020年云南省经历了严重的干旱，其覆盖范围广、持续时间长.2019年6—11月云南全区以中度干旱为主，2020年3—8月云南大部分地区以极端干旱为主，其中滇西南2020年可达异常干旱.此外，GPS-DSI、GRACE-DSI与SCPDSI-Z存在差异：(1)滇东地区的GPS-DSI比GRACE-DSI和SCPDSI-Z更湿润；(2)2019年8—11月滇西南地区的GRACE-DSI有极端干旱和异常干旱，但GPS-DSI与SCPDSI-Z都没有出现此类干旱；(3)2020年SCPDSI-Z的干旱程度相比GPS-DSI和GRACE-DSI低.本文3.2节已详细讨论上述差异的原因.

4 结论

本文反演了云南省2011—2020年的GPS-EWH，并结合GPS-EWH、GRACE Mascon产品和SCPDSI产品计算了GPS-DSI、GRACE-DSI和SCPDSI-Z，从而定量分析云南省干旱的时空分布特征.为分析干旱事件中水储量亏损情况，本文基于GPS-EWH和降水产品计算了云南省五个区域的GPS-EWH异常和降水异常，并分析其时间变化趋势.此外，展示了云南省2019年6—11月及2020年3—8月GPS-DSI、GRACE-DSI以及SCPDSI-Z的空间分布，评估了此次干旱事件的时空变化.本文的主要研究结论如下：

(1)GPS-EWH与GRACE-EWH周年振幅的空间分布一致，均在滇西南最大，并逐渐向滇东北减少.GPS-EWH与GRACE-EWH时间序列的周期性变化趋势一致，均在每年8—10月最大、4—6月最小，但GRACE-EWH时间序列的周年振幅比GPS-EWH小50.4%.此外，GPS-EWH微分、GRACE-EWH微分与降水量相关性较高，表明陆地水储量受降水的驱动.

(2)滇西北、滇西南与滇中GPS-DSI与SCPDSI-Z的相关系数分别为0.56、0.71与0.72，高于GPS-DSI与GRACE-DSI、GRACE-DSI与SCPDSI-Z的相关系数，表明在局部地区GPS-DSI比GRACE-DSI对干旱更敏感.但滇东南与滇东北地区GPS-DSI与SCPDSI-Z的相关性较低，对比GPS的站点空间分布，其原因可能与滇东南与滇东北地区GPS站点分布稀疏有关.

(3)分析GPS-DSI、GRACE-DSI与SCPDSI-Z的频率结果可知，云南省时常发生中度干旱，偶尔发生重度与极端干旱.同时三种干旱指数、GPS-EWH异常与降水异常结果表明，2011—2020年共有3次显著的干旱事件发生，其中极端干旱事件发生在2019—2020年.极端干旱期间滇中、滇西南与滇西北的GPS-EWH亏损值达10年内最大值，同时GPS-DSI、GRACE-DSI与SCPDSI-Z均大幅下降.

(4)结合GPS-DSI、GRACE-DSI和SCPDSI-Z分析2019—2020年极端干旱，2019年6—11月滇中及滇西大部分地区为中度干旱，少部分地区为重度干旱，2020年3—8月干旱加剧，云南大部分地区达到重度干旱与异常干旱.此次干旱在这10年内持续时间最长、干旱程度最严重.

实验结果表明，GPS-DSI具有较高的时空分辨率，可以作为传统干旱监测手段的补充.但GPS-DSI受限于GPS站的空间分布密度，因此结合GPS-DSI、GRACE-DSI与SCPDSI-Z可以更有效地评估干旱的严重程度与时空分布.此外，随着我国北斗导航定位系统的全面建成及其他GNSS导航卫星系统的不断发展，将来有望联合GNSS多星座观测数据，进行陆地水储量及水文干旱研究.

致谢感谢国家地震科学数据中心提供的GPS坐标时间序列，感谢美国德克萨斯大学空间研究中心提供的时变重力场Mascon解产品，感谢德国地学研究中心提供的负荷位移产品，感谢中国气象局提供的降水数据，感谢东英吉利大学气候研究组提供的SCPDSI数据.文中部分图使用GMT6.1开源软件绘制 (Wessel et al., 2019).