蓟县台小辛庄山洞周边短时地下水抽取对应变变化的影响

张磊 邵永新 马建英 李恩建 金大利 王静

1）中国矿业大学（北京），北京市海淀区学院路丁11号 100083

2）天津市地震局，天津河西区友谊路19号 300201

3）山东省地震局，山东济南 250014

0 引言

伸缩仪是精密测量地壳岩体两点间水平距离相对变化的仪器，适用于观测地壳应变和固体潮水平分量的连续变化，为研究地震孕育过程中水平应变的变化规律提供数据，也为地球弹性研究提供重要数据（熊仲华，2006）。抽水开采是洞体应变观测中的一种重要干扰源（金克俭等，1986），也是伸缩应变观测较为常见的干扰因素（李国斌等，2004；李国斌，2006；汪翠枝等，2009）。已有研究表明（裴晓峰，1995），蓟县地震台定点形变观测所受干扰程度最大、发生最为频繁的是机井抽水的干扰。SSY型伸缩仪安装于蓟县台东约4.5km处的小辛庄山洞内，距洞口220～260m处，自2001年5月开始正式观测，基线长39.7m，NS测向方位角144.87°，EW测向方位角59.94°。小辛庄山洞周边存在多口用于日常饮水和农田灌溉的抽水井，水井距观测山洞最近约300m，井深一般在100～200m，多分布在山洞的南北两侧，伸缩应变受抽水干扰严重，表现为相对拉张变化，尤其以NS向更为敏感。

蓟县地震台位于首都圈中东部，“燕山沉降带”的“京东台凹”的北部，是全国定点地壳形变基本台。台站附近断裂构造较为复杂，存在一系列断裂和褶皱，台站附近有蓟县山前断裂、杨庄断裂和黄崖关断裂。小辛庄山洞山体为震旦系雾迷山组灰岩（含燧石），岩层层面走向335°，倾角20°～33°，测洞南距蓟县断裂1.5km。仪器室内的日温差小于0.03℃，年温差小于0.05℃，相对湿度达98%。桑梓井位于小辛庄山洞西南约37km处的蓟县桑梓，自2004年开始正式观测，井深300m左右，含水层属承压水，其北部为蓟县山前断裂（图1）。

图1 天津地倾斜台站与断裂分布

为了研究蓟县台小辛庄山洞伸缩应变受周边短时地下水抽取的影响，本文利用水位的变化获取与地下水抽取有关的信息，但是小辛庄山洞周围没有水位观测井，故文中利用天津地震前兆台网中蓟县桑梓井水位数据，通过分析小辛庄伸缩应变与桑梓井水位的年动态变化特征，判断二者对短时抽取地下水的响应是否存在一致性，进而建立水位变化量与应变变化量的关系，研究伸缩应变受周边短时抽水干扰变化幅度，以期对在季节性抽水时段伸缩应变出现短期加速拉张变化的异常判定提供帮助。

1 年动态特征分析

本文整理了蓟县台小辛庄山洞伸缩应变和蓟县台桑梓井水位2004年1月1日～2014年6月30日以来日均值数据（图2），由图可知，二者均表现出了较好的年动态变化特征：

（1）小辛庄伸缩应变自2004年以来NS、EW两方向年变形态呈现出一定的规律性变化。NS向每年年初变化基本平稳，4月前后受季节性抽水的影响呈现拉张变化，并在此基础上波动变化，6～8月份雨季到来受降雨影响呈现压缩变化，此后逐渐恢复呈相对拉张变化；EW向每年年初整体上呈现压缩变化，4月前后受地下水抽取影响呈现相对拉张变化，之后受降雨影响呈现压缩变化，此后逐渐恢复呈现相对拉张变化。

（2）桑梓水位自2004年开始观测以来具有较好的年周期变化特征。每年3月前后为水位的高值期，之后转降；6月前后降到年度最低值，下降幅度为5m左右，之后转升；到11月前后达到年度第2次高水位状态，上升幅度为3m左右，之后再次下降；到12月初达到年度第2次低水位期，下降幅度为1m左右，之后又开始上升，一直到转年的高水位期。

图2 蓟县台桑梓井水位和小辛庄伸缩应变NS向、EW向日均值曲线

图3表明，蓟县台小辛庄伸缩应变观测对抽水响应敏感，尤其是NS向更为突出。受季节性抽水影响，蓟县台小辛庄伸缩应变在4月前后加速拉张变化显著（图4），水位在3月前后持续高值平稳变化后转降。对比2004年以来历年年初桑梓井水位与小辛庄伸缩应变数据（图4）发现，伸缩两分向加速拉张与水位加速下降在时间上具有很好的同步性。虽然蓟县台小辛庄山洞与桑梓井空间上还存在一定的距离，但二者对季节性抽水的响应表现出很好的同步性变化。杨耀栋等（2011）发现，天津平原区深层地下水受开采影响表现为径流-越流-开采型动态特征，因受灌溉强开采影响，低水位多出现于5～6月份，有时出现在9～10月份；丰水期8月停采后，水位逐渐回升，有时出现高水位，大多至翌年1～3月水位最高，高水位期较最大降水期滞后5～7个月，年水位变幅4～6m。桑梓水位自2004年观测以来的年变化形态与上述区域深层地下水动态变化特征一致。另外，陈鹏等（2004）在蓟县地震台开展抽水实验研究发现，小辛庄伸缩应变观测受抽水干扰的幅度既与抽水井到形变监测仪器距离有关，也与抽水井水位的最大降深有关。

2 回归分析

蓟县台小辛庄伸缩应变加速拉张变化与桑梓井水位加速下降变化同步性良好，本文利用最小二乘线性回归分析方法研究二者之间的关系，建立伸缩应变变化量与水位变化量之间的线性回归模型，并利用t检验方法对模型进行显著性检验，最后将2014年水位实测变化量代入回归模型计算应变量的理论变化区间，验证2014年蓟县台小辛庄伸缩应变加速拉张变化与地下水抽取的关系。

2.1 回归分析方法（李博纳等，2006）

对于给定数据点 （xi，yi），i＝1，2…n，满足如下关系

图3 2012年3月12～13日蓟县台小辛庄伸缩应变NS、EW向受抽水干扰数据变化

图4 蓟县台小辛庄伸缩应变NS向、EW向与多年桑梓井水位日均值曲线的对比

其最小二乘回归直线可以表示为

方差σ2的无偏估计

相关系数为

利用t检验方法进行回归模型显著性检验。假设：H0：β1＝0，H1：β1≠ 0，当H0成立时

对于给定显著性水平α，记t为统计量（7）的观测值，若

则拒绝H0：β1＝0，即认为自变量x对因变量y有显著性影响。

给定显著性水平α，在x＝x0处，标准差为

2.2 回归分析结果

2.2.1 回归模型的建立

由于小辛庄伸缩应变两方向加速拉张与桑梓井水位加速下降变化在时间上具有很好的同步性，本文统计了2004～2013年年初桑梓水位最高时水位和伸缩应变NS向、EW向观测值以及伸缩应变两方向加速拉张变化结束时的应变最大值和水位值，计算水位最高至伸缩两方向拉张结束时段内水位和应变观测的变化量（表1、表2）。

根据表1、表2可以得到水位变化量与伸缩应变变化量的对应关系（表3）。表3中xNi、yNi分别为研究水位变化量与伸缩NS向应变变化关系的样本点，xNi为蓟县台桑梓井水位变化量，yNi为蓟县台小辛庄伸缩NS向应变变化量；xEi、yEi为研究水位变化量与伸缩 EW向应变变化关系的样本点，xEi为蓟县台桑梓井水位变化量，yEi为蓟县台小辛庄伸缩EW向应变变化量。xNi，yNi与xEi、yEi的分布如图5所示，两组数据的相关系数分别为 0.8342、0.5884，可见蓟县台桑梓井水位变化量与小辛庄伸缩应变变化量相关性良好。

表1 蓟县台桑梓井水位最高值和小辛庄伸缩应变NS向拉张最大值时水位与应变的变化

表2 蓟县台桑梓井水位最高值和小辛庄伸缩应变EW向拉张最大值时水位与应变的变化

表3 桑梓井水位变化量与小辛庄伸缩应变变化量的对应关系

由表3数据，利用2.1中的最小二乘线性回归方法，根据式（2～4）计算得到最小二乘回归直线相关参数，如表4所示。

表4 最小二乘线性回归结果参数表

由表4可得蓟县台桑梓井水位变化量与小辛庄伸缩NS向、EW向应变变化量之间的最小二乘回归直线（图5）表示为

相应的式（10）回归直线方差的无偏估计分别为

图5 蓟县台小辛庄伸缩应变变化量与桑梓井水位变化量线性拟合分析结果

2.2.2 回归模型显著性检验

本文利用t检验方法对回归模型（10）进行显著性检验，依据式（7～9）分别用表示回归模型式（10）中直线xE的t统计量的观测值。则有：

（1）对于蓟县台桑梓井水位变化量与小辛庄伸缩NS向应变变化量的回归直线：44.6858＋328.7437xN，有n＝10，假设：H0：β1＝0；H1：β1≠0，对于给定显著性水平α＝0.05，t0.025（n-2）＝t0.025（8）＝2.3060，t统计量的观测值，则拒绝H0，认为自变量xN对因变量yN有显著性影响即水位变化对伸缩NS向应变变化影响显著，置信水平达到了95%。

（2）对于蓟县台桑梓井水位变化量与小辛庄伸缩EW向应变变化量的回归直线：有n＝10，给定显著性2.0582＜t0.025（8），不满足显著性水平为α＝0.05的假设检验；于是取α＝0.1，有t0.05（n-2）＝1.8595，），则在给定显著性水平α＝0.1条件下，自变量xE对因变量yN有显著性影响即水位变化对伸缩EW向应变变化影响显著，置信水平达到了90%。

上述t检验结果表明，水位变化对伸缩应变变化影响显著，并且水位变化量与应变变化量的回归直线方程NS向置信水平高于EW向，这与实际观测中伸缩应变NS向对抽水响应更为敏感这一实测现象也是一致的。

2.2.3 应变量变化区间分析

蓟县台小辛庄伸缩应变观测受抽水干扰影响严重，一般在每年4月前后开始加速拉张变化（图6），但2014年伸缩应变两方向自3月中旬开始加速拉张变化，时间较往年同期有所提前，且变化较为突出（图7）。震例研究表明（汪翠枝等，2009），蓟县伸缩应变震前一般表现为短临前兆异常，NS向在文安地震前、EW向在卢龙地震前均出现短期加速拉张变化。针对这种异常变化，首先对比小辛庄伸缩应变和桑梓井水位日均值曲线，发现水位加速转降与伸缩两方向加速拉张变化开始的时间基本一致。另外，调查发现，此期间仪器工作正常，周边也并未发现新的干扰源，故初步判定此异常变化可能与季节性抽水有关。

图6 蓟县台小辛庄伸缩应变NS、EW向日均值多年对比曲线

本文依据表1、表2，统计了2014年水位最高时水位和伸缩NS向、EW向观测值、伸缩应变两方向加速拉张变化结束时的应变最大值和水位值，分别计算出2014年该时段内蓟县台桑梓井水位和小辛庄伸缩 NS向、EW向的实测应变变化量（表5），分别用（xN0，yN0）和（xE0，yE0）表示，相应的标准差分别用δN和δE表示。如表5所示，依据式（10）和式（9）有：

（1）当蓟县台桑梓井水位变化量xN0＝2.8841时，代入式得到2014年伸缩NS向应变变化量的理论估计值标准差δN＝408.8522，2014年蓟县台小辛庄伸缩应变NS向的实测变化量yN0＝674.4375，落入区间即（583.9633，1401.6677）内，置信水平为 95%；

图7 2014年蓟县台小辛庄伸缩应变NS向、EW向与桑梓井水位日均值对比曲线

表5 2014年桑梓井水位最高值和小辛庄伸缩应变拉张最大值时水位与应变变化

（2）当蓟县台桑梓井水位变化量xE0＝2.1929时，代入式得到2014年伸缩EW向应变变化量的理论估计值、标准差δE＝472.3435，2014年蓟县台小辛庄伸缩EW向的实测应变变化量yE0＝481.05，落入区间即（-93.7461，850.9409）内，置信水平为 90%。

对于给定任意xN0和xE0，yN0和yE0的变化区间理论上下限曲线如图5中虚线所示。2014年蓟县台桑梓井水位的实测变化量和小辛庄伸缩NS向、EW向的实测应变变化量（2.8841，674.4375）、（2.1929，481.05）在理论控制线内。因此，伸缩应变在2014年3月中旬以来出现的加速拉张变化与季节性抽水有关，作为地震前兆异常的信度不高。

2.3 地下水抽取对小辛庄伸缩应变观测影响的机制分析

地下水抽取过程实质上为以井孔为中心的含水层的疏干过程，由于含水层失水产生地表沉降，形成以井孔为中心的局部收缩区域，受其影响井孔中心收缩区以外受到拉张作用，这就解释了小辛庄山洞附近短时抽水对伸缩应变的拉张影响。Terzaghi于1925年提出有效应力原理（王大纯等，2002）可以解释抽水引起地下水位动态变化进而引发地面沉降的问题。有效应力即为实际作用于砂层骨架上的应力，等于总应力减去孔隙水压力。抽水造成地下水水位下降，孔隙水压力降低，有效应力增加，砂层颗粒发生位移，颗粒接触面积增加，排列更趋紧密，孔隙度降低，砂层受到压密，从而产生地面沉降，含水层的孔隙度、给水度、渗透系数等参数均将变小。砂砾类岩土基本呈弹性形变，待砂层孔隙水压恢复后，砂层能大体恢复原状。

已有水文资料表明蓟县台小辛庄山洞附近地下水类型属于潜水，利用地下水潜水完整井稳定流Dupuit公式可以得到抽水井影响半径，该区域旱季和雨季井水位降深对形变干扰的影响范围在800～1000m（陈鹏等，2004；李恩建等，2006）。杨耀栋等（2011）的研究也表明，天津地区抽水开采是影响深层地下水位动态的主要因素，地下水开采与全淡水区的浅层地下水位动态具有良好的相关性。文中回归分析结果也表明，桑梓水位的短时下降变化与伸缩应变的短时加速拉张变化相关性良好，二者之间的回归模型置信水平达到了90%，同时也说明了利用桑梓水位作为区域地下水动态背景信息研究伸缩应变受抽水影响的短时变化是具有实际意义的。

3 结论与讨论

蓟县小辛庄伸缩应变观测受抽水干扰影响显著，每年年初受季节性抽水影响NS、EW两方向会表现出不同程度的短期加速拉张变化。桑梓水位多年来年动态特征良好且与区域深层地下水动态特征一致，在季节性抽水时段，水位加速下降与伸缩应变加速拉张变化时间上具有很好的同步性。本文利用天津地震前兆台网的观测数据开展了地下水抽取与伸缩应变短时变化关系的研究，得到如下结果：

（1）文中统计并计算了2004～2013年水位最高值至伸缩应变拉张结束时水位变化量和应变的变化量，利用最小二乘线性回归方法得到了水位变化量与应变变化量的回归直线（式（10）），二者呈正相关，即伸缩应变受抽水影响的程度随着水位变化量的增大而增大，这与陈鹏等（2004）抽水实验的研究结果也是一致的。

（2）水位变化量与伸缩NS向和EW向应变变化量的线性回归模型分别满足显著性水平为α＝0.05和α＝0.1的t检验，即回归模型的置信水平分别为95%和90%。回归模型置信水平较高，水位变化量与NS向应变变化量的回归模型置信水平更高，这与实际观测过程中伸缩应变NS向对抽水响应更为敏感的客观事实一致。

（3）伸缩应变两方向在2014年3月中旬开始加速拉张变化，考量年变形态发现，时间较往年有所提前，且变化较为突出。而调查中并未发现新的干扰源，且仪器工作正常。利用文中的回归分析方法计算得到2014年伸缩两方向应变变化量的理论变化区间，并且2014年应变的实测变化量在变化区间内，证明了该变化为季节性抽水所致，作为地震前兆异常的信度较低。

桑梓水位反映了区域地下水的动态背景，在季节性抽水时段桑梓水位变化量与小辛庄伸缩应变变化量的回归模型置信水平达到了90%以上，回归模型显示伸缩应变变化受水位变化影响显著。由于小辛庄山洞附近缺乏相应的水位观测，因此可以利用桑梓水位作为辅助资料研究伸缩应变受季节性抽水短时干扰的变化幅度。随着观测资料的积累，样本量逐渐增加，回归模型可以不断进行修正并进行实例验证，这有利于辅助判定在季节性抽水时段伸缩应变出现短期加速拉张变化的异常信度。另外，在受抽水干扰影响显著的形变观测台站，可以开展水位观测项目，以辅助研究抽水对地形变观测的短时影响。