引言

地下水具有分布广、易流动和不可压缩等特征，当井-含水层系统处于封闭性良好的承压体系中时，井水位系统即为一个天然体应变计(Hsieh et al, 1987).外力作用会引起含水层介质体应变的微小变化，而井水位可对此作出灵敏的响应，其响应灵敏度可用水位与体应变的振幅比来表示(刘序俨等，2009).在诸多外力作用中，地震波是比较常见的一种，其周期性的张、压作用可引起含水层介质的弹性变形，进而导致井水位出现周期性的振荡现象，通常将水位仪所记录到的这种振荡称为水震波或水文地震图.

对水震波的研究已近百年，水震波的记录几乎与地震波的记录具有同样悠久的历史(Blanchard，Byerly, 1935).水震波响应机理的研究一直是热点问题，早期研究发现水震波的产生是长周期瑞雷波作用的结果(Eaton，Takasaki, 1959; Rexin et al, 1962)，并得到了井水位对地震波作用的理论解析解(Cooper et al, 1965)，进一步研究表明水震波幅值与含水层有效水柱高度和含水层厚度等密切相关(Liu et al, 1989).近期的研究成果显示，地震波不仅能引起井水位水震波，还可改变井-含水层系统的水文参数，例如导水系数和储水系数等(Brodsky et al, 2003; Elkhoury et al, 2006; Manga et al, 2012)，该结论得到了岩石物理实验的验证(Elkhoury et al, 2011).

普遍认为，在地震波作用于井-含水层系统时，其周期性张、压作用既可引起地壳介质的弹性变形和地面振动，也可使岩体介质出现裂隙(Montgomery，Manga, 2003; Pride et al, 2004)，或使裂隙中的固结物重新排列(Wang et al, 2004; Mays, 2010)，而上述作用过程均可在井水位信息中得以反映.尽管如此，随着观测资料的不断积累和研究的继续深入，远场地震引起的井水位同震响应逐渐呈现出复杂性，其响应机理也并非唯一.例如，单一机理并不能解释一次地震所引起的多口观测井的水位响应现象(Yan et al, 2014; Ma，Huang, 2017)，也不能解释一口井对多次地震的响应(Weingarten，Ge, 2014).水震波的响应特征和幅度，除了受地震震级和震中距的影响外，还与井-含水层的介质参数、构造位置等密切相关(Barbour, 2015; Shi et al, 2015; Yan et al, 2016).

除了机理研究取得进展外，水位数据的采集技术也有所改进，从最早的模拟记录，发展到当前的数字化分钟值及秒值记录.众所周知，水震波的周期一般在十几秒至二十几秒之间，相比传统的模拟记录或数字化分钟值采样，秒采样率的数字化水位记录更能清晰完整地记录到远场地震所引起的水震波(舒优良等，2006；He et al, 2017)，为水震波机理及影响因素的分析提供更为丰富的资料(Shih et al, 2013；舒优良等，2014；Sun et al, 2015).

2016年11月25日新疆阿克陶发生M_S6.7地震，新10井数字化高频采样水位仪记录到了该地震所引起的水震波.本文拟对比分析其水震波与地震波地表垂向运动之间的相关性特征及其与井-含水层系统介质参数的关系，以期进一步理解地下水在地震孕育和发生过程中的动态变化.

1.   新10井构造环境与井孔条件

新10井位于乌鲁木齐市南部，其地理坐标为(43.70°N，87.62°E)，海拔高度为1056 m，位于柳树沟—红雁池逆冲断裂及其派生断裂的交汇部位(图 1a)，沿断裂挤压破碎带广泛出露泉水，多数泉水流量常年稳定，涌水量达72—108 m³/h.

图  1  新10井构造环境(a)及其井孔结构(b)图

F₁：柴窝堡盆地南缘断裂；F₂：红雁池断裂；F₃：雅玛克里断裂；F₄：西山断裂；F₅：二道沟断裂；F₆：阜康断裂; F₇：兴地断裂；F₈：柯坪断裂；F₉：西昆仑北缘断裂；F₁₀：布伦口断裂

Figure  1.  Geotectonic environment (a) and borehole structure (b) of Xin10 well

F₁: Southern margin fault of Chaiwopu basin; F₂: Hongyanchi fault; F₃: Yama-Kerrey fault; F₄: Xishan fault; F₅: Erdaogou fault; F₆: Fukang fault; F₇: Xingdi fault; F₈: Kalpin fault; F₉: Northern margin fault of West Kunlun; F₁₀: Bulunkou fault

下载: 全尺寸图片 幻灯片

新10井于1980年9月成井，其井孔柱状图如图 1b所示，井深为28 m，0—9.2 m深处井孔直径为146 mm，9.2—13.04 m深处井孔直径为130 mm，13.04 m以下井孔直径为110 mm.井孔在16—24 m深度处岩芯破碎，为主要出水段，故在该处设置滤水管.含水层为二叠系红雁池组硅质砂岩和砾岩，厚度约为8 m，井水位埋深为1 m左右. 2016年10月对该井井水位观测仪器系统进行升级改造，改用中国地震局地壳应力研究所研制生产的SWY-Ⅱ型数字水位仪，分辨率为1 mm，采样率由原来的1分钟提升为1秒，信息记录能力得到明显提高.

2.   新10井水位同震响应分析

新10井记录到了2016年11月25日22点24分新疆阿克陶M_S6.7地震的数字化水震波波形(图 2a)，震中距为1234 km.为了对比分析该井水震波对地震波的响应，本文收集了距新10井东北侧约10 km处的水磨沟台SLJ-100型力平衡加速度计强震仪记录到的垂向速度波形数据(图 2b).对比二者的波形图可以看出，新疆阿克陶M_S6.7地震引起的新10井水震波的形态与水磨沟台地震波垂向分量的形态整体基本一致，局部形态略有差异，且可清晰地看到P波、S波和面波.

图  2  新10井水震波(左)与水磨沟台地震波(右)的形态(a, b)及时频特征(c, d)对比

Figure  2.  Waveforms (a, b) and time-frequency characteristics (c, d) of hydroseismogram in the Xin10 well (left) and those of the seismic waves in Shuimogou station (right)

下载: 全尺寸图片 幻灯片

为了更加清楚地比对新10井水震波与水磨沟台地震波信号的相似性，本文利用S变换法提取了二者的时频特征. S变换类似于连续小波变换，是局部谱的一种表示，通过对局部谱在时间上进行简单的平均操作得到傅里叶谱，其不仅保留了每个频率的绝对相位特征，而且与傅里叶变换一样，具有无损可逆性(Stockwell et al, 1996).图 2c, d分别为利用S变换得到的新10井水震波和水磨沟台地震波垂向分量的时频特征图.从图中可以看出，水震波和地震波均存在两个明显的周期信号，即6—10 s (0.1—0.15 Hz)和15—30 s(0.03—0.07 Hz).其中，15—30 s长周期信号仅出现在面波到达的初始时段，这也是面波的主要频段，而6—10 s短周期信号一直持续约15分钟.从图 2c, d水震波和地震波的变化幅度对比来看，6—10 s周期段内水震波和地震波的幅度值均较高，但15—30 s周期段内，地震波的幅值相对较低，而水震波的幅值相对较大，说明新10井水位对该周期段的体应变具有较高的放大效能.

2.1   时间域响应特征

从水震波与地震波的对比可知，二者的持续时间、信号周期及变化幅度均存在一定的相似性，尤其在面波到达后的20分钟内.图 3a为在面波时段内各时刻水震波(水位)和地震波垂向速度的幅度散点图.从统计关系来看，水震波的变化幅度与地震波垂向分量的变化幅度具有简单的线性相关性，即地震波振幅越大，水震波的振幅也越大，反之则越小，且当水震波幅度较大时(|h|>13 mm)，二者的振幅比(k=138)明显高于整体的(|h|＜30 mm)平均振幅比(k=79).水震波与地震波的振幅比反映了井水位对地震波引起的地壳介质体应变的放大作用，可见当地震波振幅较大时，井水位对其放大作用也较大.

图  3  时间域水位与地表运动速度垂向分量的幅度对比

(a)水位与垂向速度散点图; (b)不同周期τ的水位幅值变化; (c)不同周期τ的垂向速度变化; (d)不同周期τ的水位与垂向速度振幅比m

Figure  3.  Amplitude comparison between water level and vertical component of ground motion velocity in time domain

(a) Scatter diagram of water level versus vertical velocity; (b) Amplitude of water level in different periods τ; (c) Vertical velocity changes in different periods τ; (d) Amplitude ratio m of water level to vertical velocity in different periods τ

下载: 全尺寸图片 幻灯片

为了进一步分析二者的相关性，本文将面波阶段的水震波和地震波的垂向速度在每个周期内的波峰、波谷值及对应的峰值时间逐一进行统计.若以波峰、波谷幅值作为响应振幅，两次峰值的时间差作为响应周期，则可得到不同周期信号的振幅变化分布.图 3b，c分别为水震波和地震波的时序数据中不同周期信号的振幅统计结果.可以看出，时序数据信号中水震波和地震波在不同周期下的振幅变化较为离散，但二者的振幅平均值的变化存在一定的相似特征，即整体上均表现为周期越长振幅越大，仅8—11 s周期信号的振幅随周期的增大略有减小.

图 3d为不同周期下水震波与地震波的振幅比，由于统计得到的周期并非一一对应，因此仅显示了二者共有周期的振幅比.可以看出，由时序曲线统计得到的水震波与地震波的幅值比虽然较为离散，但比值多为100左右，此结果与图 3a中利用所有数据统计得到的结果基本一致.

2.2   频率域响应特征

利用时序数据统计得到的振幅比，虽然在一定程度上能反映水震波与地震波之间的相关性特征，但难免有遗漏之处.这是由于为水震波或地震波是多种周期成分的波形信号相互叠加的综合表现，单纯从峰值及其对应时间统计得到的振幅并不能完整地表示其叠加信号的分项特征.图 2中水震波和地震波的S变换结果，清晰地显示了其不同周期成分的叠加现象.为此，本文利用快速傅立叶变换对新10井水位和水磨沟台垂向速度分量数据进行了时频转换，其结果如图 4所示.从图中可以看出，水位和垂向速度在5—50 s周期(频率为0.05—0.20 Hz)中，频谱幅值最大的频段均介于0.10—0.15 Hz (7—10 s周期)范围内，但在频率低于0.08 Hz(周期大于12.5 s)的频段，水位整体高于垂向速度.

图  4  水位(a)与地表运动速度垂向分量(b)的频谱对比

Figure  4.  Spectrum comparison between water level (a) and vertical component of ground motion velocity (b)

下载: 全尺寸图片 幻灯片

对图 4中不同频段的水位与垂向速度进行比较，得到与时间域统计结果相同的认识.虽然二者的对应关系较为离散(图 5a)，但整体呈正相关，即水位越大, 垂向速度也越大，且相比时间域的统计结果，频率域得到的振幅比更为全面和完整，可以呈现出不同周期段振幅比的变化特征(图 5b)：当频率大于0.08 Hz (周期小于12.5 s)时，水位和垂向速度的振幅比随着频率的减小而增大；当频率小于0.08 Hz时，振幅比则基本转平或有所下降(在100—150倍之间波动).说明在0.08 Hz频段(周期12.5 s)左右，新10井水位对地震波引起的地壳变形的放大作用最强.

图  5  频率域水位与地表运动速度垂向分量的散点图(a)及其振幅比m (b)

Figure  5.  Scatter diagram (a) and amplitude ratio m (b) of water level to vertical component of ground motion velocity in frequency domain

下载: 全尺寸图片 幻灯片

结合时间域(图 3)和频率域(图 4，5)的分析结果可知，新疆阿克陶M_S6.7地震引起的新10井水震波与地表垂向运动速度分量之间具有一定的相似性，二者数据信号的周期基本一致，响应幅度呈正相关，振幅比介于100—150之间，且在高频阶段(大于0.08 Hz)振幅比随着频率的减小而增大.

2.3   同震响应机理分析

地震波引起的水震波，其响应幅度与井孔条件、含水层参数及地震波周期等密切相关.地震波除了会引起地表介质的上下波动之外，还会引起含水层内孔隙压力的波动变化.井水位的振荡由地表垂向运动和孔隙压力变化共同造成，且井水位对二者均有一定的放大作用.因此，水震波响应幅度既受到井-含水层系统内孔隙压力波动的影响，也受到地表垂向运动的影响.在承压性的井-含水层系统中，井水位对孔隙压力和地表垂向运动的响应可用放大因子(或振幅响应比)来表示(Cooper et al, 1965)，即

式中，A为水位对含水层孔隙压力波动的放大因子，A′为水位对地表垂向位移的放大因子，r_w为井孔内径，H_e为观测井有效水柱高度，T为导水系数，S为储水系数，τ为地震波周期，ω为地震波角频率，g为重力加速度，Kei和Ker分别为零阶开尔文函数的虚部和实部.

若含水层孔隙压力波动引起的水位变化为h、地表垂向运动引起的水位变化为s，则水位整体波动为x=h+s，放大因子A=x/h.令R为含水层孔隙压力波动引起的水位变化与地表垂向运动引起的水位变化之比，m为水位整体波动幅值与地表垂向位移幅值之比，则有

式中：E_w为水的弹性模量，E_w=22×10⁸ N/m²；ρ为密度，ρ=10³ kg/m³；v为地震波波速；n为含水层孔隙度.

依据上述公式，利用实测的水震波和地震波在各周期段内的最大幅值，可计算得到相应的R和M.本文中地表运动(地震波垂向分量)以速度表示，其质点位移可由速度转化得到，即s_max=v_maxτ/2π.

令放大因子比R′=A/A′，则水位变化振幅比R与地震波周期τ成反比关系，而放大因子比R′与周期τ成二次线性正比关系，如图 6所示.图中计算R实测值所需的波速v(2.82 km/s)由震中距除以面波到时与发震时刻之差获得，周期τ可依据时序曲线各峰值与谷值的时间差来计算或通过频谱图读取，孔隙度n取0.15；计算R′实测值只需知道有效水柱高度H_e，即可依据H_e=H+3d/8计算得到，H为隔水顶板以上的井内水柱高度，d为含水层厚度.新10井含水层埋深为16—24 m，则含水层厚度d为8 m，隔水顶板埋深为16 m，同时水位埋深为1 m，则井内水柱高15 m，因此H_e=H+3d/8=18 m.

图  6  水震波变化幅值比R (a)和放大因子比R′ (b)

Figure  6.  Amplitude ratio R of water level changes (a) and amplification factor ratio R′ (b) of hydroseismogram

下载: 全尺寸图片 幻灯片

计算结果显示：新10井面波时段水震波的水位变化幅度比R介于30—300之间，可见孔隙压力波动引起的水位变化幅度远大于地面垂向运动所引起的水位变化幅度，因此该时段内的水震波主要是由地震波作用下的含水层孔隙压力波动所致；新10井放大因子比R′介于0.1—50之间，10 s周期以下的R′＜1，即水位对含水层孔隙压力波动的放大作用小于对地表垂向运动的放大作用，10 s周期以上的R′>1，即水位对含水层孔隙压力波动的放大作用大于对地表垂向运动的放大作用.

由上述分析可知，依据时序图和频谱图可逐一统计得到各周期τ的m值，依据其R值可计算得到各自的A和A′，结果如图 7所示.计算中井径r_w=54 mm，含水层厚度d=8 m，有效水柱高度H_e=18 m，储水系数S=8×10^-6，渗透系数K值拟合结果为4×10^-2 cm/s.

图  7  不同渗透系数K下新10井水震波放大因子的理论值和实测值

(a)水位对含水层孔隙压力波动的放大因子A; (b)水位对地表垂向位移的放大因子A′

Figure  7.  Theoretical values and measured values of amplification factor for the Xin10 well under different permeability coefficients K

(a) Amplification factor A of water level fluctuation to aquifer pore pressure; (b) Amplification factor A′ of water level to vertical displacement of ground surface

下载: 全尺寸图片 幻灯片

从实测值与理论计算值的对比结果来看(图 7)，新10井水位对面波引起的孔隙压力波动的放大作用较小(A值多小于1)，这与该井观测含水层的导水能力(与渗透系数K、含水层厚度d相关)或承压性较小有直接关系.另外，由图 7可见，在面波时段水震波放大因子(A或A′)的变化较为离散，有可能是因为在地震波作用过程中，由于周期性的张、压作用，岩体介质的裂隙形态及其中的填充物均会重新排列，进而引起含水层系统透水性质的变化(Manga et al, 2016)，即井-含水层系统透水性的波动变化.

3.   讨论与结论

2016年11月25日新疆阿克陶M_S6.7地震引起了新10井水位的同震响应，本文分别从时间域和频率域对比分析了该井水位仪所记录的水震波与地震波信号的相关性特征，并基于水位同震响应机理估算了新10井观测含水层的渗透系数.分析结果显示：

1) 新疆阿克陶M_S6.7地震引起的新10井水震波和地震波均存在两个显著的周期，即6—10 s和15—30 s，其中，长周期15—30 s的信号仅出现在面波到达的初始时段，而短周期6—10 s的信号一直持续约15分钟.

2) 新10井水震波响应幅度与地震波幅度整体呈正相关，且在高频阶段(频率大于0.08 Hz，周期小于12.5 s)二者的振幅比随着频率的减小而增大，表明该井水位对周期大于12 s的信号放大效能较高.

3) 利用水震波与地震波幅度比估算的新10井观测含水层渗透系数的量级为10^-2 cm/s，且在地震波作用过程中含水层的水文参数也存在波动变化.

虽然关于井水位同震响应的研究已很多，但多基于井水位分钟值记录数据来开展.由于地震波周期大多远小于1分钟，因此分钟值记录的水震波振幅必定不够完整，以此为基础统计得到的部分分析结果可能不够准确，因此，在水位同震响应研究中，采样率的高低对分析结果的影响不容忽视.前人的研究将井水位对地震波的响应形态分为振荡和阶变(Yan et al, 2014; Zhang et al, 2015)，如果提高数据采样率，使部分井水位记录更加完整，那么以前被视为阶变的形态就有可能呈现先振荡后升高(或下降)的形态，即水位在地震波作用阶段振荡, 振荡结束后水位值未恢复至振荡前水平.因此，采样率的提高也有助于更科学地对水位同震响应类型进行合理的分类，进而依据这些分类更有针对性地进行相关研究.

另外，在地下水动态研究中常常会用到渗透系数，该参数通常基于抽水试验和模型反演来获取，前者一般在成井时完成，这是由于考虑到地震观测数据的稳定性和连续性，在正式观测后不易开展抽水试验.如果要获取这类观测井的渗透系数，只能通过模型反演来实现.常用的方法是水位固体潮调和分析法(Hsieh et al, 1987)，即利用水位对潮汐的响应幅值和相位滞后来估算其水文参数.但是，由于井-含水层系统条件的差异，并非所有观测井(如新10井)的水位均能反映出明显的固体潮汐变化，有些井水位则根本无潮汐形态.因此，对于此类观测井，利用水位对面波(瑞雷波)的同震响应来反演其井-含水层的水文参数是一种不错的选择.由于远场地震引起的水震波变化只能为地震波作用所致，其结果的可信度较高.井水位对地震波的响应振幅也与井孔条件和地震波周期等密切相关(Cooper et al, 1965)，而影响井孔条件的一个重要因素即为井孔的共振频率ω，可表示为ω=，换算为周期τ=2π.新10井的共振周期约为8.5 s，可见其对该周期地震波的放大作用最为明显.

井水位对远场地震的响应较为复杂，其机理也不尽相同，仍需研究人员不断积累和探索.与本文类似，近年来也有水震波与地震波的对比研究(Shih et al, 2013; Shalev et al, 2015)，或逐步利用采样率较高的秒水位数据进行同震响应特征的分析(He et al, 2017).如何利用观测得到的宝贵数据来进一步探讨同震响应机理，并从中挖掘出更多的具有一定科学价值的信息，是当前研究的热点.本文在探讨新10井水位响应机理的过程中，反演得到了新10井观测含水层的渗透系数，这也是一种新的尝试，并期望在今后的工作中积累更多的资料对其进行进一步验证.