干旱指数作为监测和评价旱情的重要工具，可以直观、准确、定量地表征干旱强度、时空分布等特征^[1]。传统干旱指数的计算依赖于地面单基站观测与遥感反演等手段，如以降水与地表径流为基础的代表性指标Z指数、标准化降雨指数SPI(standardized precipitation index)和标准化径流指数SRI(standardized runoff index)等^[2-3]。但基于单基站观测数据计算的干旱指数受数据质量和站点覆盖面限制，难以监测大尺度的干旱事件；基于传统光学或微波遥感数据反演的干旱指数仅包含地表水、树冠水与浅层地下水变化导致的气象干旱事件，难以表征陆地水总储量变化引起的水文干旱事件。

GRACE/GRACE-FO系列卫星可观测地球重力场变化，监测地表水资源储量变化情况，估计陆地水储量盈亏。因此，可利用GRACE卫星数据反演等效水高并进一步计算干旱指数，识别并定量表征水文干旱特征^[4-5]。然而，由于GRACE实际空间分辨率较低(300 km×300 km)，对小区域的陆地水储量变化不够敏感，不适合小区域水文干旱监测；同时，GRACE任务结束和GRACE-FO任务开始之间存在数据空白，会影响对水文干旱的评估效果。因此，需要寻找一种新的干旱监测手段作为补充。

GNSS垂向位移时间序列包含环境负载(非潮汐大气负载、非潮汐海洋负载、水文负载等)引起的地球弹性形变信息，并且对小区域的环境负载变化敏感^[6]。扣除非潮汐大气负载和非潮汐海洋负载等影响得到仅包含陆地水文负载信号的GNSS垂向位移时间序列，可以用来反演陆地水储量变化及监测因陆地水储量变化引起的干旱事件^[7]。国内学者在GNSS垂向位移监测干旱方面开展大量研究^[8-9]。但这些研究仅将GNSS垂向位移与GRACE-DSI、CI(comprehensive meteorological drought index)、PDSI(Palmer drought severity index)等干旱指数进行对比分析，未将其转化为干旱指数并提供干旱严重程度的定量表征。

随着中国大陆构造环境监测网络的建立，10多年来连续运行的GNSS基准站为水文大地测量提供了宝贵数据。本文从GNSS原始时间序列中提取仅包含水文负载信号的垂向位移时间序列，然后计算GNSS垂向位移的标准化异常数据GNSS-DSI，并将GNSS-DSI与GRACE-DSI进行对比，最后利用GNSS-DSI分析云南地区的干旱特征，为干旱监测与评估提供新的工具。

1 研究区域与数据来源 1.1 研究区域

云南省位于我国西南地区，地势北高南低，海拔落差悬殊，气候兼具低纬气候、季风气候、山原、高原气候等特点。全省包含长江、珠江、元江、澜沧江、怒江、大盈江6大水系，年降水量充沛。虽然云南水资源量很大，但在时空分布上极不均衡，并且受大气环流异常、水汽输送减少、青藏高原和社会环境等因素影响，云南省干旱频发^[10]。因此，开展干旱监测对主动抗旱减灾、旱情风险管理与生态环境保护具有重要作用。

1.2 数据来源 1.2.1 GNSS坐标时间序列

从国家地震科学数据中心获取云南地区26个GNSS台站2011~2019年的时间序列进行分析，站点分布如图 1所示，其结果由GAMIT/GLOBK软件解算得到，该产品已在数据解算过程中扣除固体潮、海潮、极潮等因素的影响^[11]。GNSS时间序列中包含由构造运动引起的趋势性变化、由环境负载引起的季节性变化和地壳运动与仪器更换引起的阶跃。其中，季节性变化主要体现在垂直方向，主要影响因素为大气、非潮汐海洋和陆地水储量负载^[12]。因此，本文仅采用GNSS时间序列中的垂直分量进行分析。

1.2.2 GRACE及降雨数据

从美国得克萨斯大学空间研究中心CSR(http://www2.csr.utexas.edu/grace/)下载GRACE/GRACE-FO RL06 V02 Mascon数据集。该数据产品已考虑椭球效应，并应用所有适当的修正，包含C₂₀、C₃₀替换、地心一阶项改正、冰川均衡改正(GIA)、加入大气海洋模型(GAD)、扣除冰后回弹、非潮汐海洋负载和大气负载影响等^[13]。本文直接应用该数据集，以研究区每个站点为圆心、300 km为半径的格网值计算站点处等效水高(EWH)时间序列，并将其转换为月尺度GRACE-DSI^[14]。除2017-06~2018-06外，缺失的GRACE数据采用样条插值方法进行补充，该方法通过对邻近点数值进行三次插值得到。降雨数据从全球降水气候学中心(http://gpcc.dwd.de)获取，时间分辨率为1个月，空间分辨率为0.25°×0.25°。

2 研究方法 2.1 GNSS水文负载信号提取

为提取仅包含陆地水文负载变化引起的地表弹性形变信息，使用德国地学研究中心GFZ(https://www.gfz-potsdam.de/en/home/)提供的环境荷载产品消除非潮汐海洋负载(non-tidal ocean loading, NTOL)和非潮汐大气负载(non-tidal atmospheric loading, NTAL)效应的影响。从图 2可以看出，GNSS垂向位移时间序列受非潮汐大气负载效应影响较大，周年振幅在2~3 mm之间；受非潮汐海洋负载效应影响较小，周年振幅最大为0.5 mm。最后对垂向位移时间序列进行建模，基于最小二乘法估计所有参数并去除时间序列中的长期线性趋势及跳变，以减小构造运动(如地震、火山喷发等)等其他因素的影响，得到以水文负载信号为主导的垂向位移时间序列^[15]。模型公式为：

$ \begin{aligned} y\left(t_i\right)= & a+b t_i+c \sin \left(2 \pi t_i\right)+d \cos \left(2 \pi t_i\right)+ \\ & e \sin \left(4 \pi t_i\right)+f \cos \left(4 \pi t_i\right)+ \\ & \sum\limits_{j=1}^{n_k} g_i H\left(t_i-T_{k j}\right)+v_i \end{aligned} $

(1)

式中，a为初始位置常数项(测站位置)；b为线性速率；c、d为年周期项系数；e、f分别为半年周期项系数；t_i为时间；$\sum\limits_{j=1}^{n_k} g_i H\left(t_i-T_{k j}\right)+v_i$为跳变改正项；g_i为跳变振幅(阶跃)；T_kj为跳变发生的历元；n_k为跳变个数；H为海维西特阶梯函数；v_i为观测噪声(残差序列)。

为保证计算干旱指数时数据的完整性，使用GMIS软件补齐缺失的GNSS数据。该软件可以很好地恢复GNSS时间序列中缺失的数据，并且减弱时间序列中白噪声的影响^[16]。以YNCX、YNHZ、YNMH站为例，使用GMIS软件得到滤波与插值补齐后的时间序列(图 3)。图中黑色为扣除长期趋势、阶跃、非潮汐海洋负载和非潮汐大气负载等影响的GNSS垂向位移时间序列，红色为滤波与插值补齐后的垂向位移时间序列。

2.2 水文干旱指数计算

GRACE对大尺度的陆地水储量变化十分敏感，该优势促进了其在干旱监测领域的发展。Zhao等^[17]通过GRACE数据反演等效水高(EWH)，计算EWH的标准化异常数据GRACE-DSI，结果表明，GRACE-DSI可以独立描述干旱事件。但GRACE存在长期数据空白，且受空间分辨率限制，会影响评估精度。受上述干旱评估方法启发，本文仅基于GNSS垂向位移数据计算干旱强度指数GNSS-DSI，计算方法与GRACE-DSI相同，且不存在数据缺失。

水文负载与其引起的地球弹性形变呈负相关关系，即水文负载增加，GNSS台站垂向位置向下运动；水文负载减小，GNSS台站垂向位置向上运动。为方便与GRACE-EWH、GRACE-DSI进行对比分析，本文对所有垂向位移数据进行反符号处理，利用处理后的垂向位移数据(vertical crustal displacement, VCD)进一步计算GNSS-DSI。GNSS-DSI定义为GNSS-VCD的标准化异常，计算方法是通过月内平均将每日VCD时间序列重新采样得到每月VCD时间序列，减去每月(1~12月)VCD平均值(气候平均值)得到VCD异常，再除以VCD的标准差。计算公式如下：

$ {\rm{GNSS - DS}}{{\rm{I}}_{i, j}} = \frac{{{\rm{VC}}{{\rm{D}}_{i, j}} - {{\overline {{\rm{VCD}}} }_j}}}{{{\sigma _j}}} $

(2)

式中，i为年(2011~2019年)；j为月(1~12月)；VCD_{i, j}为第i年第j个月的VCD；VCD_j为所有年份第j个月的垂向位移平均值；σ_j为所有年份第j个月的垂向位移的标准差；VCD异常为VCD_{i, j}－VCD_j，用于评价当月的VCD亏损。

为评估GNSS-DSI的可靠性，结合GRACE数据进行对比，按照上述方法计算GRACE-DSI，并将干旱指数等级划分成11个类别(表 1)。

3 结果与分析 3.1 VCD、EWH与DSI及其相关性分析

为评价GNSS量化水文干旱程度的能力，本文对GNSS和GRACE观测的均匀分布在云南不同地理位置(中部、东北部、东部、东南部、南部、西南部、西部、西北部)的8个台站的VCD、EWH以及进一步计算得到的DSI进行相关性分析。图 4为研究区26个GNSS站点处GNSS-VCD、GRACE-EWH和2种DSI的相关性系数以及位于不同相关性系数区间的台站数量统计(无相关[0~0.2]、弱相关(0.2~0.4]、中等相关(0.4~0.6]、强相关(0.6~0.8]、极强相关(0.8~1])。从图 4(a)可以看出，整个云南地区的GNSS-VCD与GRACE-EWH表现出高度的一致性，两者的相关性系数在0.58~0.94之间，位于中-极强相关区间，除东部地区存在2个台站呈现中等相关外，其他地区的24个台站均呈现强-极强相关。从图 4(b)、4(c)看出，云南地区整体的GNSS-DSI与GRACE-DSI具有较好的相关性，除呈无相关的3个台站和呈弱相关的2个台站外，其余21个台站均呈中-强相关。结合图 4和图 1中水文负载位移周年振幅可以看出，云南东部、东南部、东北部地区陆地水储量变化较弱，位于这3个区域的GNSS台站季节性变化还包含其他地球物理因素^[18]，导致这3个区域的GNSS-DSI与GRACE-DSI相关性较弱，并存在不稳定性。

图 5为8个区域的GNSS-VCD、GRACE-EWH和2种干旱指数的时间序列以及相关性系数(蓝色和黄色分别表示GRACE/GRACE-FO卫星在2016-11~2017-06期间的数据质量下降以及2017-07~2018-05期间的数据缺失)。从图中可以看出，云南中部、西南部地区GNSS-VCD与GRACE-EWH、GNSS-DSI与GRACE-DSI具有很强的相关性，GNSS-VCD与GRACE-EWH的相关系数R值分别为0.89、0.92(图 5(a)、5(f))，GNSS-DSI与GRACE-DSI的R值分别为0.78、0.74(图 5(i)、5(n))；陆地水储量年振幅较小的云南东北部、东部、东南部地区，GNSS- VCD与GRACE-EWH的R值分别为0.76、0.74、0.68(图 5(b)、5(c)、5(d))，GNSS-DSI与GRACE-DSI的R值分别为0.58、0.58、0.45(图 5(j)、5(k)、5(l))；其他3个地区的GNSS-VCD与GRACE-EWH均表现出良好的一致性，R值分别为0.83、0.93、0.91(图 5(e)、5(g)、5(h))，GNSS-DSI与GRACE-DSI也具有一定的相关性，R值分别为0.64、0.65、0.64(图 5(m)、5(o)、5(p))。从DSI变化趋势来看，8个区域的2种DSI分别在2012、2013和2019年出现较低值，最低达到-2以下，显示出3次较为严重的水文干旱事件。

由于GNSS台站在接收信号过程中会产生部分随机噪声和虚假的季节性信号，而GRACE存在长期的数据质量恶化及数据空白，并且GRACE空间分辨率较低(300 km)，对小区域水文负载变化不及GNSS-VCD与GNSS-DSI敏感，从而影响其对水文干旱事件的评估效果。因此，上述因素导致2种干旱指数存在一定差异，结合2种数据更有助于分析水文干旱事件。

3.2 GNSS-DSI表征水文干旱事件

为证明GNSS-DSI表征水文干旱事件的能力，本文调查云南地区的水文干旱特征。图 6为8个区域的GNSS-VCD、GNSS-VCD气候平均值以及GNSS-VCD异常和GNSS-DSI与降水(PRECI)异常时间序列的比较。表 2为云南地区2011~2019年通过GNSS-DSI探测到的水文干旱事件及其特征。结合图 6和表 2可以看出，GNSS-DSI探测到了云南地区发生的持续时间和严重程度不同的水文干旱事件，其中，2011-08~2012-03、2012-11~2013-04、2019-04~2019-12这3个时期的干旱事件最为显著，尤其是从2019-04开始持续9个月的特大干旱(D4)，各区域的VCD异常与PRECI异常均出现大幅下降，VCD平均亏损值和亏损峰值分别为-4.5 mm、-16.1 mm，DSI值最小为-2.43，期间东北部、东部、东南部区域从2019-08左右开始，西部、西北部从2019-10左右开始，干旱严重程度逐渐减弱。2011-08~2012-03干旱事件严重程度相对较弱，除2012-02左右DSI出现波谷，VCD亏损峰值最低为-14.0 mm，DSI值最小为-1.75，等级为极端干旱(D3)，其他时期均为中度-严重干旱(D1-D2)。2011-08~2011-11干旱主要集中在东北部、东部、东南部区域，2011-12左右开始扩散到整个云南地区，于2012-03逐渐恢复。曹言等^[19]基于TRMM降雨数据也探测到本次干旱事件。2012-11~2013-04干旱事件与2011年比较类似，但持续时间较短，DSI于2013-03左右出现波谷，VCD亏损峰值最低为-15.1 mm，DSI值最小为-1.96，等级为极端干旱(D3)，同时期南部、东南部的干旱强度相对其他区域较弱。金燕等^[20]也识别到本次干旱事件。

除3次比较严重的干旱外，GNSS-DSI还探测到4次持续时间短、强度较弱的干旱。第1次为2015-06~09左右，主要分布在云南东部、西部、西北部(图 6(c)、6(g)、6(h)、6(k)、6(o)、6(p))，VCD异常均为负值，DSI在2015-08左右出现波谷，3个区域的DSI最低值分别为-1.29、-1.4、-1.29，等级为中度-严重干旱(D1-D2)。第2次为2016-07~10左右，主要分布在云南东部和东南部(图 6(c)、6(d)、6(k)、6(l))，本次干旱与2015年相似，DSI在2016-09左右出现波谷，DSI最低值分别为-1.42、-1.11，等级为中度-严重干旱(D1-D2)，通过查阅中国气象局国家气候中心干旱事件监测信息，证实了本次干旱事件。第3次为2018-03~07左右，干旱范围较大，除云南西南部以外其他7个区域均有干旱，DSI在2018-04出现波谷，其中东南部强度最大，VCD平均亏损值和亏损峰值分别为-6.1 mm、-12.7 mm，DSI最低值为-1.81，强度为极端干旱(D3)，其他地区显示为中度-严重干旱(D1-D2)，于维等^[21]也发现了本次干旱事件。第4次为2018-10~12左右，主要分布在云南东北部，2018-11 DSI达到低值-1.41，等级为严重干旱(D2)，VCD平均亏损值与亏损峰值分别为-1.1 mm、7.4 mm，同时期东部、东南、西北也出现中度干旱(D1)。Jiang等^[22]利用GNSS垂向位移监测云南地区陆地水储量变化，同样发现了本文除2016年以外的6次干旱事件，并指出期间陆地水储量均呈亏损状态，说明该时期个别月份的PRECI异常虽为正值，但少量降水并未起到缓解干旱的作用，导致GNSS-DSI仍呈下降趋势。

上述研究结果表明，基于GNSS-VCD的标准化异常GNSS-DSI可成功识别和量化因陆地水储量亏损引起的水文干旱事件，说明GNSS-DSI可作为水文干旱程度的评价指标应用于实际水文干旱监测。

4 结语

本文将GNSS-VCD与GRACE反演的GRACE-EWH进行对比分析，发现GNSS-VCD可以很好地体现陆地水的季节性变化。基于此，本文利用GNSS-VCD进一步计算GNSS-DSI，然后对比分析GNSS-DSI与GRACE-DSI，并利用GNSS-DSI研究云南地区水文干旱事件，得到以下结论：

1) 云南省整体的GNSS-DSI与GRACE-DSI具有较好的相关性，共21个台站呈中强相关，中部、南部、西南部、西部、西北部5个区域的GNSS台站GNSS-DSI与GRACE-DSI相关性最好。东部、东南部、东北部个别GNSS台站季节性变化包含其他地球物理因素，导致这些区域的GNSS-DSI与GRACE-DSI相关性较弱。

2) GNSS-DSI探测到云南地区7次水文干旱事件并可提供定量表征，期间陆地水储量均呈亏损状态，最为严重的是2019年干旱事件，强度大且持续时间长。7次干旱期间降雨异常多为负值或呈下降趋势，个别月份降水异常虽为正值，但由于降水量不足、气候变化和人类活动等因素影响，陆地水储量仍处于亏损状态，导致GNSS-DSI为负值。

3) GNSS-DSI可以作为GRACE-DSI的补充，将GNSS-DSI与GRACE-DSI相结合，可以提高水文干旱监测的时空分辨率，更全面地了解干旱特征，改善对水文干旱的预测、评估与管理。