研究表明，GNSS时间序列中的周期性信号可用来探测地表负荷变化及变形^[1-3]，其中地表水文负荷质量的变化是引起地表负荷变形的主要原因^[4]。水文质量的时空分布变化使岩石圈变形^[5]，影响活动断裂带的应力积累状态^[6]，这种应力变化可达数kPa^[7]，在一定程度上会影响断裂带的活动性。四川地区降水丰富，年降水量在900 mm以上，夏季降雨量为冬季降雨量10倍以上，季节性水文负荷变化显著^[8]；该地区位于南北断裂带中部，毗邻菱形块体边界，构造运动强烈且具有明显的分区性，地震活动频繁^[9]，可为季节性水文负荷对地壳变形及断层活动影响的研究提供天然实验场地。

目前，国内外学者对水文负荷对地壳变形及活动断层的影响进行了大量研究^[10-16]。由于四川地区复杂的断层运动模式与季节性水文负荷具有很强的时空相关性，本文通过调查该区域季节性负荷位移与季节性负荷应力应变特征，揭示季节性水文负荷对地壳形变的影响及负荷应力扰动对活动断裂带应力积累的调节模式，以研究不同构造区域季节性负荷应力的扰动模式，从而更深入地认识地壳深部的动力学过程，更细致地了解活动构造的变形机制。

1 数据与方法 1.1 季节性信号提取

GNSS时间序列中的季节性信号主要由环境负载所引起，包括大气负荷、海洋潮汐负荷、非潮汐海洋负荷、非潮汐大气负荷和水文负荷等^[17]，本文将着重分析水文负荷的影响，因此需尽可能扣除其他环境负载的影响。四川及周边地区GNSS连续运行参考站时间序列数据可从中国地震局GNSS数据产品服务平台(ftp://ftp.cgps.ac.cn/products/position/)获取，大气负荷和海洋潮汐负荷在GAMIT解算过程中分别使用atl.grid大气潮汐模型和FES2004海洋潮汐模型进行剔除；非潮汐海洋负荷和非潮汐大气负荷使用ESMGFZ产品中的非潮汐海洋负荷(NTOL)模型和非潮汐大气负荷(NTAL)模型进行剔除。为保证结果的准确性，本文剔除观测时间小于2.5 a的数据，得到四川及周边地区共52个GNSS站的观测数据，GNSS坐标时间序列函数模型可用式(1)表示:

$ \begin{array}{l} y\left( {{t_i}} \right) = {y_0} + v{t_i} + A{\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{t_i}} \right) + B{\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{t_i}} \right) + \\ C{\rm{sin}}(4{\rm{ \mathsf{ π} }}{t_i}) + D{\rm{cos}}(4{\rm{ \mathsf{ π} }}{t_i}) + \\ \sum\limits_{j = 1}^k {{g_j}} H({t_i} - {T_{gj}}) + {\varepsilon _i} \end{array} $

(1)

式中，t_i为以a为单位的历元；y₀为测站的正常位置；v为长期速率项；A、B、C和D分别为谐波函数振幅；H(t_i-T_gj)为由于地震和仪器更换导致的跳变项；ε为观测噪声。

本文基于式(1)构建GNSS坐标时间序列方程，利用最小二乘法对待估参数进行求解。其中，Asin(2πt_i)+Bcos(2πt_i)和Csin(4πt_i)+Dcos(4πt_i)为采用谐波函数表示的周期项，分别为周年和半周年变形信号，主要由呈周期变化的信号源(如季节性负荷变化)所引起。由于半周年变形约为周年变形的1/3^[16]，相较于周年变形，半周年变形太小，因此本文利用式(2)对周年变形信号进行分析，可求得任意时刻地壳的季节性位移：

$ \begin{array}{l} y({t_i}) = A{\rm{sin}}(2{\rm{ \mathsf{ π} }}{t_i}) + B{\rm{cos}}(2{\rm{ \mathsf{ π} }}{t_i}) = \\ {\rm{Amp}} \times {\rm{sin}}(2{\rm{ \mathsf{ π} }}{t_i} + {\rm{Pha}}) \end{array} $

(2)

式中，$ {\rm{Amp}} = \sqrt {{\rm{ }}{A^2} + {B^2}} ;{\rm{Pha}} = {\rm{arctan}}\left( {B/A} \right)$。

对于式(2)中时间t_i的选取，本文将1 a分为24个区间，每个区间约为15 d，这样既具有一定的时间跨度，又可保证数据在时间上的连贯性，可直观地看出四川地区受季节性水文负荷影响而产生的地壳运动变化。利用插值大地位移计算水平应变，由季节性负荷位移场得到季节性负荷应力应变场^[17]。

1.2 降水数据提取

为证明负荷位移及应力应变场模型的合理性，并方便对模型结果进行解释，利用四川地区降水数据进行验证。美国国家航空航天局GPM任务可提供全球雨雪观测数据，四川地区降水数据可从其网站(https://storm.pps.eosdis.nasa.gov/storm/)下载。本文选用四川地区2001-01~2019-08时间分辨率为月、空间分辨率为0.1°×0.1°的降水数据，原始数据格式为每格网当月平均每小时降水量，通过计算提取四川地区每月平均降水量，降水随时间的变化特征如图 1所示。四川地区既受东亚季风和印度洋季风影响，也受青藏高原环流系统影响^[18]，因此年内降雨干湿季节明显，年内降雨主要集中在5~9月，且均在7月出现降雨量峰值。

2 水文负荷季节性变化影响分析 2.1 季节性位移与降水相关性

本文基于式(2)计算四川地区1 a的季节性水文负荷位移，并对各站出现最大位移的月份进行统计，结果如图 2所示，图中箭头长短代表位移大小，箭头方向代表测站位移方向，不同颜色代表不同月份。

对于水平位移，上半年63%的测站在3~5月达到最大位移，下半年最大位移集中出现在9~11月；对于垂直位移，上半年73%的测站在2~4月达到最大位移，下半年则集中在8~10月。为更清楚地看出季节性位移与降水的关系，本文按照四川地区季节划分(春季为3~5月，夏季为6~8月，秋季为9~11月，冬季为12月至次年2月)，利用2001~2019年的降水数据得到月平均降水分布，结果如图 3所示。由图可知，四川地区春、夏、秋、冬四季降雨量的空间分布特征总体上相似，均呈东南向西北逐渐降低的趋势。其中，春季降雨主要集中在四川盆地，占全年降雨量的22%；夏季降雨量明显增多，主要集中在盆地东部和川西南山地地区，月平均降雨量达152 mm，占全年降雨量的51%；秋季降雨量高值区由川西南山地地区转移到盆地东部地区，占全年降雨量的23%；冬季降雨区域变化相较于其他季节尤为显著，高值区范围明显缩小，仅占全年降雨量的4%。四川地区年降水量在900 mm以上，夏季降雨量为冬季降雨量10倍以上，季节性水文负荷变化显著。

与降水数据对比分析发现，季节性位移与降水之间存在明显的相关性，最大位移集中出现在地面水文负荷逐渐达到最小值的春季和地面水文负荷逐渐达到最大值的秋季。图 4为四川地区平均水平位移、平均垂直位移与月平均降雨量的变化关系，从图中可以看出，平均水平位移在3~5月达到最大值，随着降雨量增多，位移逐渐减小，在9~11月达到最小值；平均垂直位移在3月达到最大值，4~7月随着降雨量增多，位移逐渐减小并在降雨量相对较大的9月达到最小值，10~12月随着降雨量逐渐减小，位移逐渐增大，地壳回弹上升，这与姜卫平等^[19]的研究结果一致，说明GPS观测值与降雨量在时间尺度上具有一致性。由此推测，降水可能是影响季节性位移周期性变化的主要原因。

从图 2可以看出，巴颜喀拉块体和川滇菱形块体的垂直位移在2~4月均达到最大，而四川盆地在6月才达到最大，这与四川地区降水分布和地形起伏有关，为此本文沿102°~106°E、31.5°N范围作地形剖面发现，龙门山两侧最大高差可达4 000 m，结合降水和地形分析可知，由于巴颜喀拉块体和川滇菱形块体在2~4月降雨量较小，地下水流失，蒸发作用增强，导致地面水文负荷减小，地壳向上回弹，从而达到垂向最大负荷位移。虽然四川盆地在2~4月降雨较少，但其地势较低，与西部山区最大高差可达4 000 m，且2~4月仍有冰雪融水和地下水汇入，因此四川盆地地壳上升的最大负荷位移出现在5~6月，进一步证明水文负荷是影响四川地区垂向变形的主要原因。

从图 2还可以看出，鲜水河断裂和安宁河-则木河断裂带西侧的川滇菱形块体附近的最大负荷位移量明显大于四川盆地与巴颜喀拉块体地区，水平最大位移达2.3 mm，垂直最大位移达9.4 mm，且两者达到最大负荷形变量的时间不同，这可能与四川地区地壳运动有关。龙门山断裂带右旋速率≤0.5 mm/a，逆冲断层收缩速率约为1~2 mm/a^[20]；四川盆地位于相对稳定的华南块体之上，变形量较小，但鲜水河断裂的滑动速率约为9~11 mm/a，安宁河-则木河断裂带的滑动速率约为7.9 mm/a^[21-22]。相较于四川盆地和龙门山断裂带，川滇菱形块体地区地壳运动活跃，因此本文推测，板块运动影响负荷位移大小，水文负荷变化调节最大位移出现的时间。图 2中水平季节性负荷位移较分散，但仍可以看出，四川盆地水平季节性负荷位移在2~3月达到最大，而巴颜喀拉块体和川滇菱形块体则在4~5月达到最大，且上半年四川地区地表季节性负荷位移整体由西南向东北运动。安宁河-则木河断裂带附近的地壳运动方向与断层走向近似垂直，从而在上半年加剧断裂带的逆冲活动；而龙门山断裂带附近的地壳运动方向与断层走向平行，在上半年加剧断裂带的走滑活动。

通过对比水平和垂向负荷位移发现，季节性水文负荷是导致四川地区地壳垂向负荷形变的主要原因，在地壳升降的过程中同时产生水平形变，并且季节性水文负荷对断层运动速率更快的川滇菱形块体地区的影响更大，使该地区地壳季节性负荷变形更明显。

2.2 水文负荷对活动断层应力应变场的扰动模式

地壳应力应变能直观反映地壳的相对运动趋势和变形程度，利用长期GNSS数据基于地壳变形理论模型结合区域地质构造背景计算区域应力应变场分析其时空分布特征，是研究地震活动与地壳变形关系的重要内容。本文由24份季节性负荷位移场计算得到24份水平应变场，分析该地区全年水平负荷应变的变化情况。面应变可反映地壳中应变能的状态及大小^[23]，本文计算四川地区各月负荷面应变平均值，并与降水数据进行对比，结果如图 5所示。由图可知，在降水较少的12月和1月负荷面应变达到最大值0.5×10^-9，在降水较多的6~7月负荷面应变达到最小值-0.5×10^-9，全年负荷面应变呈周期性变化，并与降雨量存在很强的相关性。另外，本文还对比4月和10月的负荷面应变、最大最小负荷主应变及经过内插处理的季节性负荷位移场，结果见图 6。从图中可以看出，4月四川地区受季节性水文负荷影响，地表负荷位移整体向东北方向运动，但川滇菱形块体地区季节性负荷位移大于巴颜喀拉块体地区和四川盆地，运动趋势受到巴颜喀拉块体和四川盆地阻挡，这与该地区负荷面应变为负值表示地壳收缩、负荷主应变表现出主压应变的情况一致。相反，10月地表负荷位移整体向西南方向运动，川滇菱形块体地区季节性负荷位移大于巴颜喀拉块体地区和四川盆地，面应变为正值表示地壳膨胀，主应变表现出主拉应变特征。由此可知，四川地区季节性负荷与地壳形变及断层应力应变变化存在很强的相关性。

为计算库仑应力，本文采用King等^[24]的计算方法，库仑破坏准则为：

$ {\sigma _f} = \tau - {\mu _f}{\sigma _n} $

(3)

式中，σ_n和τ为破坏面上的法向应力和剪切应力，μ_f为有效摩擦系数。

断层面上法向应力σ_n和剪切应力τ可分别由式(4)求得^[25]：

$ \begin{array}{l} {\sigma _n} = {\sigma _{xx}}{\rm{si}}{{\rm{n}}^2}\theta - {\tau _{xy}}{\rm{sin}}\theta {\rm{cos}}\theta + {\sigma _{yy}}{\rm{co}}{{\rm{s}}^{\rm{2}}}\theta \\ \tau = \frac{1}{2}({\sigma _{yy}} - {\sigma _{xx}}){\rm{sin}}2\theta + {\tau _{xy}}{\rm{cos}}2\theta \end{array} $

(4)

式中，θ为断层面与x轴(右旋剪切应力方向为正)的夹角；σ_xx、σ_yy和τ_xy为水平应力分量，可通过水平应变分量和逆胡克定律计算求得。

本文杨氏弹性模量取70 GPa，泊松比取0.25^[26]，通过逆胡克定律由水平应变分量计算得到水平应力分量σ_xx、σ_yy和τ_xy；再根据式(4)计算断层面上法向应力σ_n和剪切应力τ；最后取有效摩擦系数μ_f为0.4^[24]，并由式(3)计算得到断层面上的库仑应力。考虑到四川地区主要断裂带运动方式存在差异，因此本文进行分区计算，分区情况及各区域库仑应力如图 7所示，其中区域1为鲜水河断裂带，区域2为龙门山断裂带，区域3为安宁河-则木河断裂带。由图可知，3个区域的库仑应力值不尽相同，但都与季节性负荷具有很强的相关性。图 8为以龙门山断裂带北川-映秀段为例分析断层受荷载影响的结果，从图中可以看出，当夏季降水增多地面荷载增大时，沿断层面向下会产生与原应力方向相反的作用力，对断层原有的运动趋势产生一定抑制作用；相反，当冬季降水减少地面荷载减小时，沿断层面向上会产生与原应力方向相同的作用力，对断层原有的运动趋势产生一定促进作用。由图 7(b)~7(d)可知，3个区域的库仑应力与面积变化的相位非常接近，表明水文负荷导致地壳变形的同时会调节断层的库仑应力。在6~7月降水增多、地面水负荷增大时，库仑应力均为负值，代表应力得到释放，会减小断层长期积累的库仑应力；12月至次年1月降水减少，地面水负荷减小，库仑应力为正值，代表应力积聚，会增加断层周围的库仑应力。这也证明四川地区季节性水文负荷对断层长期积累的库仑应力具有一定的调节作用。

通过对比四川地区的水文负荷及断层应力应变场的变化情况可以看出，研究区水文负荷与断层应力应变之间存在很强的相关性，在降水增多、地面荷载增大的月份，断层周围的库仑应力减小，表明部分库仑应力得到释放，会降低断层大范围破裂的可能性；龙门山断裂带区域受水文负荷影响产生的库仑应力在±4 Pa之间变化，仅为鲜水河断裂带和安宁河-则木河断裂带区域库仑应力变化的1/10(±40 Pa)，季节性负荷对龙门山断裂带区域负荷应力积累的调节作用非常微弱。考虑到汶川地震的震后变形仍是不可忽略的因素，认为水文负荷对龙门山断裂带的影响相对较小。

3 结语

本文基于四川及周边地区52个GNSS站的观测数据，提取四川地区季节性负荷位移信号，构建季节性负荷应力应变场模型，研究四川地区季节性水文负荷对地壳形变及活动断层的影响，得出以下结论：

1) 四川地区季节性负荷位移受季节性水文负荷调节。GNSS站最大水平位移为2.3 mm，最大垂直位移为9.4 mm，全年呈周期性变化；巴颜喀拉块体和川滇菱形块体区域垂向负荷位移与水文负荷存在很强的相关性，四川盆地垂向负荷位移也呈现周期性变化，但与巴颜喀拉块体和川滇菱形块体相比，四川盆地垂向最大负荷位移出现的月份存在2~3个月滞后；此外，本文推测板块运动影响该区域最大负荷变形值，水文负荷会影响达到最大位移的时间。

2) 四川地区负荷应力应变场受季节性水文负荷调节。水平面应变、主应变与季节性水文负荷存在很强的相关性，在水文负荷较大的6~8月面应变为负值，地壳收缩，主应变为主压应变，这与季节性负荷位移场一致。

3) 龙门山断裂带、鲜水河断裂带和安宁河-则木河断裂带等区域的库仑应力均与季节性水文负荷具有很强的相关性，在水文负荷增大的月份，库仑应力为负值，表示断层部分长期积累的库仑应力得到释放，对断层活动具有一定的调节作用，但龙门山断裂带区域库仑应力的变化十分微弱。

致谢: 感谢中国陆态网提供四川省CORS站GPS观测数据，感谢美国航空航天局网站提供降水数据，本文部分图片由GMT软件生成，在此一并表示感谢。