0 引言

降雨是一种离散的不连续干扰因素，其对地形变的干扰主要由雨水从地表向地下渗透引起地块膨胀和收缩所致，为非线性干扰，干扰图像具有群体性质。降雨对地块作用必须有一个能产生水分转移的势梯度和一个进行水分转移的连续途径（Kundu et al，2016），降雨干扰与许多外界条件和初始条件有关，且与地形地貌、地下径流、渗透率、孔隙率以及初始含水量密切相关。法国Darcy通过实验提出，地下水渗流为一维的均匀运动，渗流速度与水力梯度成线性关系，即渗流速度与水力梯度的大小和方向沿流程不变，薛禹群（1997）对地下水动力学原理进行深入研究，将Darcy定律推广到三维空间。

我国对形变降雨干扰研究已有先例，如：刘序俨（1985）采用褶积回归法进行降雨干扰研究；刘序俨等（1991）建立数学物理模型，采用动态系统法对降雨进行干扰定量排除；张德涵（1994）基于Darcy定律，建立垂直形变观测受地下水运动支配的地面沉降模型，并采用有限元法给出数值解，由于台站自身样本数量较少，无法建立有效模型；张凌空（1995）提出体应变降雨干扰响应模型，将降雨和体应变的响应过程分为2个时间函数，即降雨造成的体应变积累时间函数和雨后体应变衰减时间函数。若要准确得到相应的时间函数，则测点的含水层渗透系数、系统时间常数等参数需为已知。这类参数专业性较强，如含水层渗透系数是一个动态参数，选取适合测点的含水层渗透系数值是建立降雨干扰数值模型的关键。本文采用降雨造成的体应变积累时间函数，通过徐州地震台体应变降雨干扰观测结果，运用数值模拟进行拟合，估算该测点含水层渗透系数范围，通过含水层渗透系数值的密度分布，选取较为符合测点的含水层渗透系数值，即可得到该台降雨干扰数值模型，明确降雨对体应变观测数据干扰的能力，厘清体应变干扰物理机制，为日后降雨期间地球物理异常核实以及消除降雨对体应变干扰提供基础，为区域地震监测预报提供参考。

1 台站简介及仪器参数

2006年12月12日徐州台安装TJ-2型体积式应变仪，探头深度60 m，2007年7月1日运行，2017年10月遭受雷击停测。2019年5月，井孔钻探加深，安装新的TJ-2型体应变仪，探头深度64.8 m，同年12月30日入网运行。

徐州台体应变观测井孔岩芯完整，岩性为寒武系中统厚层灰岩，地表为约厚2 m的粘土夹碎石层（图 1为观测井的钻孔柱状图）。井区地下水补给与径流条件良好。当地每年6—10月进入雨季，降雨量充沛。

2 体应变降雨响应特征分析

选取徐州台2014—2021年体应变观测数据，分析各年雨季所受降雨影响，发现该台体应变观测曲线呈正零漂变化。为免赘述，文中仅给出2020年7月20日—23日和2021年9月4日—5日降雨干扰曲线，见图 2。

基于部分降雨事件中徐州台体应变增量∆a（单位量：10^-9）及对应增量时间∆t（单位：h）、降雨量W（单位：mm），计算得到降雨系数β和降雨效率η，计算公式如下

$ \beta=\Delta a / W $

(1)

$ \eta=\beta / \Delta t $

(2)

统计徐州台降雨量及相应干扰参数数值，结果见表 1，可知：徐州台体应变对降雨的响应较为灵敏，平均趋动降雨量为80.5 mm，但最小趋动降雨量为41.8 mm（降雨序号8）。这是因为，2016年8月6日，当地短时集中降雨，导致雨水快速渗入地下形成地下径流，从而使得体应变对降雨响应变得更为灵敏，说明趋动降雨量随着地下介质的含水饱和度或降雨的缓急不同而不断发生变化。

基于表 1所示降雨参数，绘制2014—2021年徐州台体应变降雨效率柱状图，见图 3，可见：2020—2021年体应变降雨效率更高，其平均降雨效率为2014—2017年的5.5倍，查询表 1可知，该时段平均降雨系数为2014—2017年的4.0倍，表明随着时间的递进，体应变对降雨的响应时间在缩短，对降雨的响应幅度在增大，分析认为：①2020—2021年多为暴雨事件，即降雨量较大，降雨时间较短，雨水充分渗入地下，改变了介质密度以及重力分布，引起局部区域加速形变；②2019年安装的体应变仪，钻孔深度更深，钻孔直径更小，地下水结构可能有所改变，降雨后地下介质含水量更易饱和（薛禹群，1986）。

3 体应变降雨模型计算结果及分析

据张凌空（1995）提出的体应变降雨响应时间函数，曲线积累过程为

$ Y_1=\beta W\left(1-\mathrm{e}^{-\Delta t /\left(B_1 \cdot \tau\right)}\right) $

(3)

式中，Y₁为降雨产生的应变量；β为降雨系数；W为降雨量；∆t为降雨造成的体应变积累时间；B₁为常数；τ为系统时间常数。据Kundu等（2016）的研究，有

$ \tau=S / T $

(4)

式中，T为含水层导水系数；S为井孔面积，即S = 16 277 mm²。据薛禹群（1997）的研究，有

$ T=K M $

(5)

式中，K为含水层渗透系数；M为含水层厚度，即M = 2 000 mm。

运用以上体应变降雨数学物理模型，由实际观测结果进行数值模拟计算，求得含水层渗透系数K，具体数值见表 1。对徐州台含水层渗透系数K进行拟合，结果见图 4，可知：含水层渗透系数变化较大，降雨量越大，降雨产生的应变量越大，含水层渗透系数越高。原因在于，地形变受降雨干扰主要是由雨水从地表向地下渗透引起地块的膨胀和收缩所致，降雨量及降雨效率越大，产生的地块膨胀及载荷效应越大，降雨干扰是一个复杂的非线性渗流问题，雨水渗入量通常受降雨强度、降雨持续时间、岩土初始含水量以及入渗面几何特征影响。

基于表 1的统计结果，得到徐州台含水层渗透系数K平均值为17.7 mm/h。然而，徐州地区常见短时强降雨，短时间内的极大降雨量，使得降雨效率出现偏离正常降雨事件的极大值，因此在计算过程中K值易出现极大值，若K值平均值将出现不适用普遍正常降雨的偏大现象，不能真实反映徐州台含水层渗透系数。为此，基于表 1中的降雨量统计结果，计算并分析该区K值核密度分布，核密度分布原理是，根据有限数据样本，运用核密度函数，对整体数据密度进行估计，即已知有限数据样本和一个核函数，输出整体数据的概率密度，其核密度图原理为

$ \hat{p}_n(x)=\frac{1}{n h} \sum\limits_{i=1}^n K\left(\frac{s_i-x}{h}\right) $

(6)

式中，x是属于整体数据中的任一变量；n是样本数据中数据点的个数；h是平滑带宽，控制平滑量，h值越大，数据越平滑；K是核密度函数，s_i是样本数据中的一个数据点称。由此对K值进行核密度研究，分析在有限降雨数据样本下，K值的概率密度分布，从而求解其核密度平均值，即在K值概率密度分布较为集中区域（图 5中上下2条红色线为取值区域）的平均值，以此作为体应变降雨积累模型计算中的K值，结果见图 5，可见：K值在6.4—24.1 mm/h区间，概率密度分布较大，其核密度平均值为15.5 mm/h。分析认为，含水层渗透系数K = 15.5 mm/h时，更接近真实含水层渗透系数，适用于降雨量为75—120 cm的降雨事件。因此，据K = 15.5 mm/h设计徐州台体应变降雨模型，模拟结果与体应变实际观测值更为接近。

将S = 16 277 mm²、K = 15.5 mm/h、M = 41 200 mm代入式（4）、式（5），求得系统时间常数τ，则τ = 0.026。B₁为无量纲常数，根据实际观测取值21.1，据式（3），计算徐州台2020—2021年降雨系数β，取其平均值，则该台降雨系数为β = 7.25。由此得到徐州台体应变降雨积累模型如下

$ Y_1=7.25 W\left(1-\mathrm{e}^{-\Delta t / 9.79}\right) $

(7)

基于式（7）所示体应变降雨积累模型，选取2022年7月14日和7月29日降雨事件，计算徐州台体应变，与实测值进行对比，结果见图 6，图中横轴表示降雨事件的时间序列，为了较好地对比固体潮降雨干扰曲线，前12个小时形变值取固体潮实测背景值，后续部分为模型计算值与实测值对比。

由图 6可知，在2次显著降雨事件中，利用该模型计算所得形变值与真实形变的相关系数均在0.8以上（表 2），模型响应效果较好。对于岩土初始含水量较低的首次降雨事件（如2022年7月22日降雨事件），二者的相关系数有所降低。原因在于，2022年徐州地区较为干旱，全年降雨量为225.5 mm，约为2019—2021年年均降雨量的25%，地下介质含水量较低，首次降雨事件中，2.2小时的集中降雨使得大量雨水从地表向地下渗透，地下介质含水量充分饱和后出现非线性渗流，进而导致降雨干扰变成更为复杂的地下水结构问题。同时，研究发现，2次降雨事件中拟合形变值均小于真实形变，可能是因为，观测场地周边山体特征、降雨后地下水的渗流、荷载变化等系列附加效应并未计算在内，因此降雨干扰产生的附加形变变化并未在模型中体现出来。

4 结束语

徐州台体应变受降雨干扰显著，钻孔加深安装新体应变仪后，体应变对降雨的响应特征表现为：降雨响应时间缩短、响应幅度增大，复杂的降雨方式导致体应变降雨干扰的不确定性，同时受到季节、地下介质含水饱和度的影响，体应变响应变化表现出差异性特征。

基于体应变降雨积累模型，由数值模拟可知，在降雨量75—120 cm的常见降雨事件中，含水层渗透系数K值在6.4—24.1 mm/h范围内分布密度较大，若K值取为核密度平均值，即K = 15.5 mm/h时，拟合结果更接近真实形变变化，二者的相关系数均在0.8以上，模型响应效果较好。由于该模型未考虑降雨导致的载荷变化影响，故模拟结果均小于实测值，若要进一步提高模型精度，则需根据测点所处山体的地形地貌特征、地质水文参数等建立合理的数学物理模型。