水井含水层的水文地质参数具有明确的物理基础与机理。一般情况下,获取这些参数需要水文地质工作者花费较大的人力、物力和财力。而基于气压系数和潮汐因子获取既简便易行,又可掌握水位动态变化过程信息。这对于有效地利用水位记录资料分析地震前兆异常和机理具有十分重要的意义。利用井潮、气压等资料分析井-含水层系统的固体潮效应、气压效应、以及它们与含水层参数(孔隙度和渗透系数等)等,一直以来是国内外地震工作者研究的热点(Bredehoeft,1967; Kamp,Gale,1983; Narasinmhan et al,1984; 田竹君,谷园珠,1985; 张昭栋等,1989,1995; 李春洪等,1990; Erskine,1991; 丁风和等,2007)。如在不排水状态下,一般都是给出气压系数或是潮汐因子和孔隙度、固体骨架的体积压缩系数、水的体积压缩系数间的定量关系,相应的贮水率和渗透系数等含水层参数也可获得。在川滇地区,方慧娜(2013)也仅以四川的南溪、邛崃井等为研究对象,利用井水位的气压效应和固体体积模量的实验值等,反演了孔隙度、储水率、导水系数等含水层参数,但也是取固体体积模量的实验值来完成的。而本文可滑动拟合在不排水状态下孔隙度和渗透系数的定量结果,研究时序特征分析,对动态捕捉和跟踪异常具有非常重要的意义。

本文剔除川滇地区12口井的水位资料的观测干扰因素后,利用气压系数和潮汐因子,滑动拟合得到不排水状态下,各井含水层的孔隙度、水和固体骨架的体积压缩系数,在此基础上,进一步获得各井含水层的渗透系数等,定量分析区域井孔孔隙度和渗透系数的时序变化过程。

针对水位与气压间的非线性及考虑滞后的卷积回归法,近年来一直是气压校正的主要方法之一。其核心思想是:引入井水位对气压的阶跃响应函数的概念,通过利用相对应的水位和气压数据拟合出阶跃响应函数的最佳值,再由该最佳阶跃响应函数校正水位。如果考虑固体潮因素(一般由传感器测量或理论固体潮代替),该方法可同时校正观测井水位数据中的气压和固体潮组分(Rasmussen,Crawford,1997; Toll,Rasmussen,2007; 王丽亚等,2012)。

(1)气压系数的获取:假设将水位的影响因素分解为气压、潮汐、趋势成分和随机成分4个分项。首先,用卷积回归法剔除潮汐成分,再用一般多项式分段拟合去趋势和随机成分后,得到剩余水位(只有气压效应项的水位),且趋势明显的气压数据也要事先进行线性去趋势。最后利用剩余水位和校正后的气压采用1阶差分来获取气压系数。

(2)潮汐因子等的获取:首先用卷积回归法剔除井水位的气压成分(气压校正),再用一般多项式分段拟合去趋势和随机成分,得到的剩余水位只有潮汐响应项,最后利用维尼迪柯夫调和分析程序就可获得各波群的潮汐因子等参数。

值得注意的是,由于大部分井点原始水位的趋势性上升变化,直接受区域降雨量的影响,区域地下水开采的年际变化对井水位的趋势变化有扰动作用。因此,利用卷积回归法进行水位数据中的气压和固体潮组分校正后,还要进一步采用去趋势分析方法,剔除地下水开采和降水补给等长周期干扰对水位的影响。这样可以为下一步进行含水层参数的定量分析等提供可靠的数据。

依据前人的研究结果(Bredehoeft,1967; 张昭栋等,1989,1995; 李春洪等,1990; 丁风和等,2015a,b),在不排水状态下,井水位的气压系数B_P和潮汐因子B_g可分别表示为:

B_P=(nβ)/(α+nβ)(1)

B_g=-(1-n)/(ρ_g[(1-n)α+nβ])(2)

(1)(2)式联合可得到:

B_g=-(1-n)/(ρ_gnβ[((1-n)(1-β_P))/(B_P)+1])(3)

式中:α为固体骨架的体积压缩系数; β为水的体积压缩系数; n为含水层的孔隙度; ρ_g为水的重度; 且ρ_g=0.098 hPa/mm。所以含水层的孔隙度n、水的体压缩系数β依据(3)式可滑动得到。最后利用(1)式或(2)式,可求出固体骨架的体积压缩系数α。

相应的贮水率Ss也可按照下式求出:

Ss=ρ_g[(1-n)α+nβ](4)

动态响应条件下,含水层的渗透系数K、潮汐波频率ω、导压系数a和相位滞后φ等满足以下3个公式(李春洪等,1990):

tgφ=(r²₀)/(2T)ωk₀(5)

k₀=ln(1.12)/((ω/a)^1/2·r₀)(6)

α=K/(Ss)(7)

因为在不排水状态下,相位滞后的正切值为零,也及开尔文函数的实部为零。所以含水层导压系数a可写为:

e^k₀=(1.12)/((ω/α)^1/2·r₀)=1(8)

a=(ω·r²₀)/(1.254 4)(9)

渗透系数K可写为:

K=(ω·r²₀)/(1.254 4)·Ss(10)

其中:φ为M₂波相位滞后; k₀为开尔文函数实部; r₀为井径; ω为潮汐波中某一波群的固定频率(选择振幅大、干扰少的M₂波频率,ω=1.932 4); a为含水层导压系数; K为渗透系数。最终可得到各井不排水状态下的含水层孔隙度和渗透系数。

图2～3给出了各井的孔隙度和渗透系数的时序月值曲线。各井孔隙度和渗透系数在大部分历史强震前有明显异常特征。如保山井、德宏法帕井、南溪井、腾冲井、泸州13井、姚安井、永胜井、开远井和丽江井等在汶川地震和芦山地震前

后,孔隙度和渗透系数同步性较好(孔隙度和渗透系数增大,则以构造张应力为主; 反之亦然),反应了区域构造应力场的变化(张应力或压应力)。2015年以来,大部分井孔的孔隙度和渗透系数变化相对平稳,后续需进一步跟踪分析。

针对上述大部分井孔出现的孔隙度和渗透系数同步变化现象,从应力变化(拉张和压缩)方面,构建了水动力学模型,分析其变化的成因机制。

自然状态下,含水层介质一般是由固相的骨架颗粒、相互连通的孔隙和其中的流体物质耦合在一起而形成的固、液、气三相物体,本文所研究的流体物质主要指水。含水层介质一般有4个体积:总体积(外观体积)、岩石骨架体积、孔隙体积、水的体积,且外观总体积等于岩石骨架体积、孔隙体积与水的体积之和。在含水层受力状态变化引起的孔隙率(度)和孔隙压力动态过程中(图4):当含水层介质受压力变形时,孔隙压力随之增大,则孔隙率(度)减小; 反之,当含水层介质受张力变形时,孔隙压力亦减小,则孔隙率

图2 12口井孔隙度月值曲线
Fig.2 Monthly value curves of porosity for the 12 wells

图3 12口井渗透系数月值曲线
Fig.3 Monthly value curves of permeability for the 12 wells

图4 含水层受力状态变化引起的孔隙率(度)和孔隙压力动态过程示意
Fig.4 The dynamic state and porosity of groundwater caused by the change of stress state in aquifer

(度)增加。另外,孔隙率(度)和渗透系数间存在明显的幂函数关系,孔隙率(度)在逐渐增大或减小的同时,相应地渗透系数也同步增大或减小(丁风和等,2015c)。

川滇地区作为我国最主要的活动构造带和地震活动区域之一,因其地理位置的重要性、特殊性和震情形势的严峻性而被中国地震局确定为地震强化监视地区。本文利用该区域12口井的气压系数和各井水位的M₂波潮汐因子等,滑动拟合得到不排水状态下,各井的孔隙度和渗透系数等含水层参数,并定量分析这些参数的时序变化特征,主要得出以下结论:

(1)各井孔隙度和渗透系数在大部分历史强震前有明显异常特征,孔隙度和渗透系数同步性较好,反应了区域构造应力场的变化。

(2)在含水层受力状态变化引起的孔隙率(度)和孔隙压力动态变化过程中,当含水层介质受压力变形时,孔隙压力随之增大,则孔隙率(度)和渗透系数同步减小; 反之亦然。

(3)与传统的通过现场抽水试验和室内实验等不同,本文研究结果表明,利用数字化水位等资料,结合气压系数和维尼迪科夫潮汐调和分析结果,而获取含水层介质的孔隙度和渗透系数是简便易行的。

(4)与前人所不同的是,本文得到的各井含水层的孔隙度和渗透系数是动态变化的,而前人对这些参数都是定性的分析或者是取其一般的可能值或平均值。

(5)卷积回归法(针对水位与气压间的非线性及考虑滞后)是进行气压校正和剔除潮汐成分的有效方法。

(6)进行井孔含水层的孔隙度和渗透系数的动态时序定量分析,在有效利用数字化水位记录资料来甄别地震前兆异常和机理分析具有重要意义。上述研究都是在假设含水层介质是线性、均质和各向同性的弹性体,井中的水为理想流体等情况下完成的。在实际应用过程中还需通过实践进行深入研究。

在本文研究过程中,得到了中国地震局地壳应力研究所刘耀伟研究员、中国地震台网中心黄辅琼研究员以及中国地震台网中心晏锐部长等的悉心指导和帮助。另外,中国地震台网中心张晶研究员还提供了潮汐调和分析程序。在此一并表示衷心的感谢。