潜水位微动态固体潮影响及主要影响分波识别

史昊鑫, 许 模, 葛建宏, 郭 健

(成都理工大学 地质灾害防治与地质环境保护国家重点实验室, 成都 610059)

0 引 言

固体潮与地下水之间的联系是水文地质领域的热门研究课题, 20世纪以来, 不少学者对固体潮与地下水之间的关系进行了探索。Bredehoeft[1]研究了各种要素对井水位变化的影响, 并建立了地下水位与固体潮的函数关系, 成为研究地下水对气压效应及固体潮的理论基石; Melchior[2]、Rhoads等[3]给出了固体潮没有滞后的解释; Rojstaczer等[4]给出了地下水位对固体潮和气压效应的响应方程; 王勇[5]分析指出, 地层年代越老, 岩性越坚硬, 对固体潮响应能力越强; 王铁城等[6]验证了在研究水位与固体潮的关系时, 预先消除气压影响的必要性; Merritt[7]通过分析固体潮、海洋潮汐和气压效应, 估算了美国佛罗里达州含水层系统的水力特性; 呼晶磊等[8]通过对马17井的观测, 确定了井水位存在固体潮效应，且水位与固体潮效应之间存在正相关关系; 史浙明等[9]对不同深度承压井的潮汐进行对比分析, 探索了微动态与地震的关系; 刘阳等[10]通过对潮汐效应的分析, 解释了潮汐效应机理, 并得出了潮汐效应没有滞后性; 苏乔等[11]采用时间序列分析法和互相关分析法分析了潮汐对滨海地下水位的影响并证实了地下水位与潮位频率相同; 贺前钱等[12]利用井水位和降雨数据模拟计算土壤含水率变化, 估计了地下水变化导致的重力效应; Kazama等[13]通过建立的水文模型, 导出了地下水引起的重力扰动, 并可用来纠正重力数据中的水文干扰; Wang等[14]通过分析潜水含水层中地下水的潮汐效应实现了对俄克拉荷马州一地区的含水层渗漏监测。

以往对固体潮的研究主要集中在承压含水层, 针对固体潮对潜水含水层影响的研究较少。本研究以四川省德阳市中江县垮梁子乡UG01监测井的地下水位高精度连续监测数据为基础, 采用频谱分析法和互相关分析方法, 探讨固体潮对潜水位的影响, 并利用调和分析法进一步识别影响水位动态的主要固体潮分波。

1 研究区及井概况

研究区位于四川省德阳市中江县冯店镇垮梁子乡, 丘陵地貌, 地势由东北向西南逐渐降低(图1)。区内除第四系全新统堆积层覆盖外, 基岩主要为侏罗系蓬莱镇组砂泥岩互层(图2), 岩层产状趋于水平, 裂隙较发育, 近水平的层面和两组垂向裂隙构成了地下水运移的主要通道。

图2 研究区地层Fig.2 Stratum in the study area

图1 冯店镇地貌图Fig.1 Geomorphic map of Fengdian Town

研究区主要含水岩组为厚度为1～13 m的第四系全新统(Q4)粉质黏土夹少量块石、碎石的耕植土覆盖层， 以及下伏侏罗系蓬莱镇组上段(J3p2)的厚层—巨厚层的砂岩与泥岩、粉砂岩互层, 因区内地下水位多在20 m深度以下, 故地下水类型主要为裂隙水。

监测井UG01井位于山丘斜坡中部, 位置坐标为104°53′53″N、30°38′56″E, 井口高程530 m, 初见水位506.78 m。钻井完成后, 在34.58 m深度砂岩地层中埋入渗压计, 并用粗砂和膨润土进行分层回填。

根据研究区2007—2017年的降雨量数据, 降雨主要集中在5—9月(图3), 占全年总雨量的82.47%; 10月至次年4月的降雨量仅占全年的17.53%; 10 mm以上的降雨日数年均22.5 d。为减少降雨对UG01井潜水位的影响, 此次数据采集选在秋冬季节(2018-08-20—2019-01-28), 该时段平均降雨频率为0.1 cpd(circle per day, 日循环次数)。

图3 2007—2017年月平均降雨量Fig.3 Average monthly rainfall from 2007 to 2017

2 数据来源和研究方法

2.1 数据来源

本次地下水位数据采自RST振弦式渗压计(加拿大产), 振弦式压力传感器输出为频率信号, 该频率与施加在振弦换能器膜片上的压力成比例且不受导线阻抗和接触电阻的影响, 这类的振弦式传感器可以安装在钻孔中或打入软土地基中。数据采集时间间隔设置为1 h, 精度为0.1 cm, 选取2018-08-20T00:00—2019-01-28T00:02时间段的监测数据进行研究, 地下水位变化曲线如图4(图中灰色时间段为选取的特征数据段, 具体见3.2节)。

图4 UG01井潜水位变化曲线Fig.4 Curve of groundwater level in Well UG01

研究中的固体潮监测数据来自中国科学院精密测量科学与技术创新研究院位于丽江的OSG-066超导重力仪, 采用的数据时间从2018-08-20—2019-01-28, 时间间隔1 h, 得到的固体潮变化曲线如图5所示(图中灰色时间段为选取的特征数据, 具体见3.2节)。

图5 当地固体潮变化曲线Fig.5 Curve of local earth tide

对研究区的大气压强进行实时监测, 选取的监测设备为Baro-Diver气压计, 气压数据每隔1 h采集一次, 精度为0.01 cmH2O(0.009 806 hPa)。选用的2018-08-20T00:00—2019-01-28T15: 00的当地气压变化曲线如图6所示。

图6 当地气压变化曲线Fig.6 Curve of local barometric pressure

2.2 研究方法

1)高通滤波法。 高通滤波器是通过频率筛选并过滤信号的装置, 其原理是当频率低于f0的信号输入这一滤波器时,由于容抗很大而受阻,输出减小,频率越低输出越小。 当频率高于f0的信号输入这一滤波器时,由于容抗很小,对信号无衰减作用。频率f0由下式决定

f0=1/(2πR1C1),

(1)

式中:R1为电阻;C1为电容。

2)频谱分析法。由于离散的数据信号在时域图上特征不明显, 可通过快速傅里叶变换(FFT)将时域信号转换成频域信号, 凸显数据特征。因此, 本文在研究潜水位与固体潮相关周期性时借助快速傅里叶变换, 将离散的时域信号变换为频域信号分析水位和固体潮的周期性。

傅里叶正变换公式

(2)

傅里叶逆变换

(3)

式中:ω—角频率(Hz);t—时间(s)。式(2)为时间域的函数转换为频率域的函数的积分公式。

对潜水位和固体潮进行FFT分析可以简单理解为: 把水位或固体潮变化看作是多个简谐波叠加后的结果, 而FFT则是将水位或固体潮变化的一个整体叠加波拆分成多个简谐波, 进而从频域的角度去分析, 由此便可以看到时域上不明显的周期性特征。

3)互相关分析法。互相关函数(cross-correlation function, CCF)是一种表示两个变量之间相关程度和方向的函数, 给定两个时间序列Xt和Yt(t=0,±1,…,±n),如果不是平稳的时间序列,可以经过d阶差分将其转换为平稳的时间序列,将差分后的序列记为(x1,y1),(x2,y2),…,(xn,yn)。

样本的互协方差函数

(4)

(5)

4)调和分析法。固体潮和水位可以看作是由一系列分波叠加的结果,在实际分析中一般选取有限个较为主要的分波。调和分析法可以将各个分波的振幅和迟角从数据中分解出来。

潮位表达式

(6)

式中:S0为余水位;fj为交点因子;μj为交点订正角;σj是第j个分波频率;hj、gj是分波的振幅和迟角；v0j为第j个分潮的初相位;t为时间。

2.3 数据预处理

受固体潮影响的地下水位波动为典型的微动态变化[15], 监测的水位数据中除了拥有需要的周期性变化的信号, 也具有非周期性变化的或不需要研究的周期性信号(如降雨)。一次降雨引起的地下水位变动远大于固体潮的影响, 要找到固体潮对水位的影响, 就需要对数据进行预处理, 消除潜在因素的影响。监测时间内的降雨频率约为0.1 cpd, 通过高通滤波器, 将频率为0～0.1 cpd(周期大于10 d)的低频信号消除, 去除了水位中的非周期性变化后, 进一步对水位数据中的周期性变化(气压效应)进行消除[16-17], 得到处理后的水位(图7)。

图7 预处理后的潜水位变化曲线Fig.7 Curve of groundwater level after pretreatments

3 固体潮对潜水位影响的基本特征

3.1 周期性

地球固体潮是在太阳、月亮等天体的引潮力作用下, 固体地球产生的周期性的弹性变形现象。地球周期性的膨胀和压缩变形, 导致岩石空隙的扩大和缩小, 从而使井水位出现周期性的下降和上升现象, 其变化形态与理论固体潮形态相似, 称为地下水位的固体潮效应[18]。用精密仪器可以观测到地球固体表层也有和海洋潮汐相似的周期性升降现象。固体潮具有1/2日、日、半月、年、11年、18.6年的周期性变化[19]。

研究区固体潮的周期性分析采用频谱分析法对其固体潮和地下水位数据进行快速傅里叶变换, 并绘制于半对数坐标中, 得到固体潮频谱图和潜水位频谱图。

从固体潮频谱分布图(图8a)可见, 在1、2和3 cpd处有明显的峰值, 且1、2 cpd振幅高于其余的潮峰, 表明固体潮在12、24 h的周期比8 h的周期更显著。因此, 可以判定固体潮存在明显的1日、1/2日、1/3日周期性。

由潜水位频谱分布图(图8b)可见, 潜水位在5个不同的周期内存在明显的波峰, 其中1、2和3 cpd的周期性波动特征与固体引力潮重叠, 从变化周期的角度来看, 潜水位存在与固体潮相似的1日、1/2日和1/3日周期性关联。

图8 固体潮(a)及潜水位(b)频谱分布图Fig.8 Spectrogram of earth tide(a) and groundwater level(b)

3.2 负相关性

为充分研究反映潜水位与固体潮的相关性, 从时间序列中抽取2018-11-06—13和2019-01-02—09两周的固体潮和地下水位数据进行曲线对比。水位与固体潮变化曲线如图9、图10。从图形的形态特征来看, 固体潮的峰(谷)值与潜水位的谷(峰)值对应关系明确, 固体潮与潜水位之间表现出明显的负相关性。

图9 2018-11-06—13固体潮和潜水位周变化曲线Fig.9 Curves of earth tide and water level from 6 to13 in November of 2018

图10 2019-01-02—09固体潮和潜水位周变化曲线Fig.10 Curves of earth tide and water level from 2 to 9 in January of 2019

利用式(4)、(5)对固体潮和地下水位作互相关函数计算, 得到分析结果如图11所示。地下水位与固体潮在Lag=0时的相关性系数最高, 为-0.864, 反映固体潮与潜水位呈负相关且不存在明显的时间滞后。即当固体潮增大的时候, 地下水位随即下降, 水位埋深增大; 当固体潮减小的时候, 地下水位随即上升, 水位埋深减小。呼晶磊等[8]在研究了马17井承压含水层地下水位后得出水位与固体潮呈负相关, 并解释其原因是当潮汐应力增大时, 含水岩层发生体膨胀, 含水层的孔隙水压降低, 井水位下降; 当潮汐应力变小时, 含水岩层发生体压缩, 含水层的孔隙水压升高, 井水位上升。本研究也在潜水含水层中得出了相同结果, 但是由于本研究的井结构特殊, 采用的传感器也非常规的压电式压力传感器, 故在非承压含水层中这种负相关关系的物理原因仍需进一步探索。

图11 固体潮与潜水位互相关系数分布图Fig.11 Cross-relation coefficients between earth tide and groundwater level

4 主要影响分波的识别

从前文可知固体潮对潜水位存在明显影响, 为进一步分析影响地下水位的固体潮分波, 采用调和分析法得到地下水位和固体潮的主要分波的振幅, 进而结合各波的重叠频率域来识别影响水位波动的固体潮分波[20]。调和分析法采用MTALAB实现, 从结果中提取出日波、半日波和1/3日波的振幅, 结合各波对应的频率, 得到各分波的频率振幅数据(图12)。

图12 固体潮(a)及潜水位(b)主要分波频谱图Fig.12 Main frequency components of earth tide(a) and groundwater level(b)

固体潮和潜水位的频谱分析显示, 固体潮频谱分析图中的主要成分为O1(频率 0.928 8 cpd, 振幅 31.74 μGal)、K1(频率 1.003 2 cpd, 振幅 43.70 μGal)日波, N2(频率1.896 0 cpd, 振幅12.65 μGal)、M2(频率1.932 0 cpd, 振幅64.86 μGal)、S2(频率1.999 2 cpd, 振幅30.65 μGal)半日波, M3(频率2.899 2 cpd, 振幅1.02 μGal)1/3日波; 而潜水位频谱分析图中的主要成分为K1(频率1.003 2 cpd, 振幅9.1 mm)日波, S2(频率1.999 2 cpd, 振幅9 mm)半日波, SK3(频率3.002 4 cpd, 振幅1.6 mm)1/3日波, 因此, 从两者间重叠的峰值部分可知, 固体潮对水位的主要影响分波为K1波和S2波。

5 结 论

垮梁子乡的UG01监测井为一封闭井, 其潜水位实时监测数据不受气压效应影响, 固体潮作用是潜水位微动态变化的主要原因, 主要结论如下:

(1)通过频谱分析法, 得到固体潮存在明显的8、12和24 h的周期变化, 研究区潜水位也表现出了相同的周期性波动。

(2)时间序列和互相关分析均表明, 固体潮和潜水位标高呈负相关关系, 即固体潮增大时, 水位下降, 固体潮减小时, 水位上升; 且固体潮与潜水位之间不存在明显的滞后性, 在非承压含水层中这种负相关关系的物理原因仍需进一步探索。

(3)利用调和分析法, 得出影响地下水位的主要固体潮分波为K1日波、S2半日波。

致谢: 中国科学院精密测量科学与技术创新研究院的陈晓东研究员提供了固体潮数据, 在此表示衷心的感谢!