地下水位是流体学科中非常重要的监测手段之一，可反映某些地震前兆问题^[1-4]。20世纪70年代以来，国内外地下水位气压效应研究得到深入发展，为分析技术的提高提供了依据和帮助，但目前关于地下水超采对井-含水层影响的研究较少。

昌黎何家庄地震台(以下简称昌黎井)作为承压水位观测井，在近年的观测中出现了不同程度的环境干扰，经多次异常核实后发现主要与周围地下水开采有关。2016-06~2022-10受附近温泉馆抽取同层水位影响，水位大幅下降6.07 m。为克服断流影响，台站安装了稳流泵不间断抽水，以保证气氡、气氦的正常观测。本文采用小波分析、交叉小波及小波相干分析方法，在不消除地下水位其他自然周期影响(如固体潮、海潮、多年周期变化等)的情况下，针对超采影响期间前后对昌黎井水位、气压观测资料的相关性进行研究；同时结合一元线性回归方法分析气压系数和水位相位滞后(或超前)等情况；最后对重点现象进行讨论。

1 昌黎井概况

昌黎井成井于1975-12，于1976-01正式投入观测，属承压水井。井孔地面高程为15.0 m，属热水自流井，井深301.4 m，其中0~75.76 m深度为第四系粘土细砂，75.77~85.84 m深度为中粗砂夹砾石，85.84 m深度以下为前震旦系混合花岗岩和角闪斜长片麻岩，172~180 m深度为主要出水段。0~91.92 m深度井管直径为170 mm，91.92 m至井底井管直径为110 mm。完井后下直径为168 mm的套管至91.92 m，其下为裸眼出水，水泥止水段长15 m；井内下0.4寸塑料管至79.00 m，出水口由0.4寸塑料管引出^[5]。

2 分析方法

计算过程包括3个部分：小波连续变换、交叉小波分析、小波相干谱计算。其中，小波连续变换中的母波采用Morlet计算，计算过程参考文献[6]，该方法在地球物理学研究中应用广泛，主要用于提取水位和气压的周期能量。交叉小波分析是结合小波变换和交叉谱分析的计算过程，可以从不同时间尺度研究2个时间序列在时频域中的相关性，还可以计算2个时间序列的共有周期和相位变化、相位角的相对关系，计算过程参考文献[7]。小波相干谱分析可度量时间-频率空间中2个时间序列的局部相关的密切程度，即使是2个时间序列对应交叉小波能量谱中低能量的区域，在小波相干谱中也表现为显著相关，计算过程参考文献[6-7]。

3 水位与气压相关性

深井水位在气压的作用下发生变化，该现象被称为水位的气压效应，用气压系数或气压效率表示，二者相关性变化机理参考文献[8]。

昌黎井水位动态变化首先具有显著的固体潮日变，其次与气压呈负相关，这种负相关在日或月尺度的短期序列变化中明显，长周期年变化中并不明显，原因有二：一是人类活动或季节性的影响，如春季开始的抽水灌溉间接导致观测层水位降低，雨季降雨补给水位又得到回升，这种有规律的起伏变化在部分年份比较显著，并与气压年变呈一定正相关；二是该井水位具有与天体运行一致的10 a变化周期^[9]，当部分年份水位上升或下降变化显著于其他因素时，与气压也就不存在相关之说了。

3.1 相关性的小波相干估计

首先估计年尺度正相关变化，采用2012~2019年水位和气压月均值进行分析(2016~2017年温泉馆营业初期，地下水超采量大且不规律，导致水位年变畸变)，预处理过程保证了数据的连续和准确性，并剔除了错误数据，计算结果见图 1。其次估计日尺度的负相关，随机抽取无环境干扰的2012-01 03:00~05:00整点值，结果见图 2，各图均包含水位和气压周期、交叉小波(XWT)、小波相干(WTC)估计。图 1和2中粗实线表示通过95%显著性检验的临界值^[6-7]，细实线表示连续小波变换的数据边界效应影响较大的区域；右向箭头表示水位与气压同相位，反之表示反相位；垂直向下箭头表示水位相位超前气压90°，反之表示滞后90°。色条渐变中颜色数值越高说明周期性越强或相关系数越高。

图 1显示：1)2012~2013年和2018~2019年水位动态周期在256~512 d范围内最显著，均值接近1 a尺度，同时也存在其他周期性不强的短周期成分(图 1(a))；2)气压动态周期仅有1 a时间尺度，周期性非常强，数值在8.00左右(图 1(c))；3)交叉小波(XWT)、小波相干(WTC)反映2个序列1 a尺度的周期相关性很强，且二者同相位呈正相关变化，但2012年至2014上半年以水位相位超前气压相位为主，剩余年份以滞后气压相位为主，尤其2015年变化相对不规则(图 1(b)和1(d))，说明含水层水位受外界干扰较大，非自然环境或地球物理场短时间的变化所为；4)交叉小波和小波相干块状区域内，整体看相位不论超前还是滞后，箭头指向与横坐标偏差较小，说明水位相位和气压相位的同步性较好。

图 2显示，日尺度的水位和气压变化负相关非常明显，二者均有12 h最显著周期，同时水位24 h周期也较显著。在小波交叉和小波相干粗黑实线区域内，12 h周期的箭头指向与横坐标平行或有较小向下偏角，表明二者无相位差或水位相位略超前，24 h周期信号颜色强度比12 h高，但只有箭头垂直向上的相位变化通过95%显著性检验，说明水位相位滞后气压90°的变化周期也偶有出现，但并不是主要变化特征，可能与天气变化有关。

3.2 相关性的一元线性回归计算

该方法主要提取不同年份(有干扰和无干扰时期)、不同季节的气压系数、水位滞后(或超前)时间。目前较普遍的求取气压系数的方法有回归分析法、差分法及作图法等，其中最常用的是一元回归、二元回归与一阶差分法^[10]。昌黎井水位动态变化与气压相关性较明显，属于显著相关，可采用一元回归进行分析。

参与分析的数据采用整点值，筛选原则为：1)井孔附近地震活动相对平静；2)不论水位监测有无背景干扰，数据都要相对稳定；3)气压变化幅度较大，且包含上升和下降过程。

以2016-01数据为例，假设水位滞后时间分别为－4 h、－3 h、－2 h、－1 h、0 h、1 h、2 h、3 h、4 h、5 h，依次作井水位对气压的响应特征分析(图 3)，再分别求出对应时间水位与气压的相关系数，取相关系数最大时对应时间作为滞后的最佳时间，并由此确定一元回归方程和气压系数。

由于气压变化具有随机性，本文依季节分别抽取冬季(2016-01共29 d)、春季(2015-04共27 d)、夏季(2015-07~08共31 d)、秋季(2015-10共25 d)4个时段数据，计算结果见表 1中Ⅰ类对应数据。可以看出，2016-01水位与气压最大相关系数R²为0.483，超前时间为2 h；2015-04最大相关系数R²为0.487，超前时间为3 h；2015-07~08最大相关系数R²为0.196，超前时间为2 h；2015-10最大相关系数R²为0.478，超前时间为3 h。由此推测，昌黎井水位对气压的响应特征为：超前2~3 h，夏季相关程度最差，R²平均为0.180，其余3个季节R²均值在0.40~0.49之间。根据计算结果确定，上述4个时段内水位与气压的一元线性回归方程依次为：冬季y=－0.003x+ 4.125，气压系数a_p为－0.003 m/hPa；春季y=－0.003x+3.605，气压系数a_p为－0.003 m/hPa；夏季y=－0.002x+3.056，气压系数a_p为－0.002 m/hPa；秋季y=－0.003x+3.862，气压系数a_p为－0.003 m/hPa。虽然4个季节的R²不同，但气压系数基本一致，可看作背景干扰时期的气压系数。该结果是原始整点值线性拟合的特征，原始数据中带有趋势性变化成分，为对比去趋势后的结果，表 1计算了不同滞后时间下水位与气压整点值一阶差分的相关系数，见表中Ⅱ类对应数据。可以看出，除夏季R²在水位相位超前气压相位2 h时最高外，其他季节均在超前3 h时最高，与Ⅰ类对应数据的结果一致。

表 2统计了2017-04地下水超采严重时期的参数变化。可以看出，Ⅰ类R²降至0.2左右，a_p变为－0.011±0.002 m/hPa，相关性最强时水位相位超前4 h，结果与上述研究时段差异很大；Ⅱ类仍显示水位相位超前气压相位3 h时二者相关性最强，气压系数在不同的相位变化中浮动范围很大，与Ⅰ类的气压系数差距更大。

对2012年无背景干扰时期的气压效应进行分析，同样抽取各季节50~60 d左右时段数据。表 3显示，Ⅰ类R²大多在0.2~0.3左右，a_p为－0.004 ~－0.001 m/hPa，最高R²在水位相位超前1~4 h中均有出现；Ⅱ类显示4个季节中水位相位超前气压相位3 h时R²均最高，气压系数与Ⅰ类很接近。可以看出，利用昌黎井水位和气压的原始整点值及整点值一阶差分进行线性拟合，相位变化均有显著特征，即在不同季节和不同场地环境背景下，只要动态变化处于稳定状态，水位相位就一直超前气压相位，但气压系数a_p和相关系数R²会随监测环境的变化出现较大浮动。其中，采用整点值一阶差分得到的结果更确切，基本为水位相位超前3 h，而原始整点值的计算结果浮动偏大；场地环境无干扰时，不同季节Ⅰ、Ⅱ类数据得到的a_p基本一致，但当场地环境有显著超采干扰时，计算结果没有参考价值。

4 讨论

分析显示，无论场地环境是否存在抽水干扰，昌黎井地下水位日尺度短期分析的气压效应均具有超前3 h的变化特征；同时水位和气压12 h和24 h的周期性显著，但年周期变化中水位气压效应易受超采干扰，相位由超前变为滞后。

1) 气压系数偏低。

气压对承压水位的影响不同于潜水位，是由于气压对含水层和井孔水面施加了附加应力，引起含水层中孔隙水压的增(减)量小于井孔中水所承受的大气压强的增(减)量，从而引起井孔中的水与含水层的水互相渗流，造成水位变化。通常灰岩井气压系数最大，砂岩次之，角砾岩、大理岩、玄武岩、砂砾岩居中，第四系松散砂层最低^[11]。昌黎井地表以下至75.76 m深度覆盖第四系粘土细沙，中间层为中粗砂夹砾石-混合花岗岩和角闪斜长片麻岩，主要出水段172~180 m处为破碎岩石。从整体岩性结构来看，不利于气压系数的增强，故本文分析的R²、a_p均较低。

2) 日尺度分析水位相位超前。

观测井对大气载荷的响应可分为高频、中频和低频，由于气压对含水层系统作用的是面积力，与固体潮相比其能量集中在较低频段^[12]，低频响应与含水层渗透率无关，而依赖于封闭层和非饱和带的扩散系数。在这个频段内(即承压性不理想或半承压水状态)，信号的衰减和放大及相位的滞后或超前都有可能出现^[13]。

首先，利用2种方法判断昌黎井地下水的承压性，一为傅里叶谱分析^[14]，提取井水位潮汐波中占全部引潮力95%以上的5个潮汐分波(M₂、K₁、O₁、N₂、S₂)，如果5个分波均有显示，则为承压性地下水；二为水-岩平衡状态法(Piper三角图)^[15]，该方法可生成三线图和三角图2个图件，其中三线图可用来判断地下水水质类型、成因及来源等信息，三角图可分析水-岩平衡、地下水循环深度等，还可判断含水层体系的开放和闭合条件、时间及运移过程。2种方法结合使用可提高结果的准确性，当水-岩平衡状态分析方法判定结果为未成熟水或部分平衡水，而傅里叶谱分析判定为承压水时，则可综合判定该井为混合水井^[16]。

图 4显示，5个潮汐分波在傅里叶谱中都有显示，且M₂波谱值最显著，仅N₂波谱值偏小，说明昌黎井地下水为承压水。图 5为该井2017~2019年及2022~2023年水质化验结果Piper三角图，可以看出：①昌黎井水样为Cl-Na·Ca型，与成井时类型基本一致^[17]，表明水的来源基本稳定；②水样为部分平衡水，说明有一定程度的水-岩反应，但反应不充分；③地下水中不仅有大气降水补给，还有较深层地下水的混入^[16]。综合判定，该井地下水为混合水，并非完全承压水或承压性不理想。

其次，对相位超前进行理论判定。当考虑到非饱和带气压效应不可忽略或隔水层承压性不理想时，相位出现超前情况。依据大气压对半承压井水位作用的响应关系，以井孔内水位上升为正，气压加载振幅为A，气压响应的相位θ(ω)计算公式^[13]为：

$ \theta(\omega)=\arg \left(\frac{x_0 \rho g}{A}\right) $

式中，ω为大气加载频率，x₀为井孔中水位波动幅度，ρ为水密度，g为重力加速度。当给定一定的参数值时，相位超前现象是一定存在的^[10]。

3) 月尺度分析超采对气压效应的影响。

通过小波相干分析和日常监测环境调查，本文认为有两起超采干扰事件存在：①以2014年上下半年为分界线，相位发生转折变化，由相位超前转为滞后。原始曲线显示，2014全年水位趋势下降，下半年变化趋势与气压相反，打破正常年变，全年降幅0.74 m。此事件应判定为超采干扰，因为相位发生转折且一直未恢复，意味着井-含水层参数发生变化，其中导水系数是相位滞后的主要因素，在无震例出现的前提下，只有超采干扰才能带来如此影响。②2016-06调查发现，附近温泉馆抽取同层水位，此事件实际在2015-12已经开始，但用水量不大，2016-05~2017-02为水位快速下降期，整个干扰期间水位降幅达6.07 m，这个时段的相位同样滞后，与2014~2015年有很好的衔接，说明干扰源的本质是相同的。

由Piper图分析可知，昌黎井地下水位年动态变化主要受降水影响，其次还有其他径向水流的渗入，包括导水系数极小的弱透水层和导水系数较大的含水层。附近温泉馆每天不间断抽水，水位下降速度远大于含水层补给速度，打破了原来稳定补给的平衡状态，水均衡为负，意味着上部含水层的疏干和压密。这一疏干层位属于第四系黄土或粘土层，渗透性远低于其他层位，当有降雨或其他径流渗入，弱透水层的补给和压力传导非常缓慢，导致水位年变相位滞后于气压相位。

5 结语

随着经济社会的发展，沿海地区地热开采逐年增多，造成地下水位趋势性下降及上部水位下移，对水位多年周期性变化和水化学组分的分析带来很大影响，甚至造成前兆异常的错误研判。正常情况下，井-含水层各参数的变化具有动态特征，因井而异，因震而异。虽然超采带来的水量变化难以计算，但可根据具体参数的变化特征总结出一些规律，为正确识别地球物理场异常提供帮助。本文仅对超采前后水位对气压的响应进行了研究，其他影响参数有待进一步分析。