基于数字化水位的琼东北地区构造应力场时序特征分析＊

解晓静 张艺潇 张 帆 李 盛 郑在壮

1) 海南省地震局，海南海口 570203

2) 中国地震局地球物理研究所，北京 100081

引言

区域地壳构造应力场的分析是一项较困难复杂的事情，通常通过测量应变实现应力的测量，其耗资较高。研究表明，承压含水层井水位的变化反映了含水层孔隙压力的变化，而含水层孔隙压力的变化与含水层所受的应力状态有密切的关系[1]。测量地壳深部弹性应力应变的方法有观测记录断层活动区附近钻井内流体的压力变化，而测量钻井水位的变化是研究流体压力变化的方法之一[2]。流动性和难压缩性是地壳结构中活跃组分之一的地下水所存在的主要性状，而当地下水在含水层中形成一个封闭条件较好的承压系统时，其水位的升降变化就能较客观地反映出周边区域地壳应力应变的状态。构造应力变化将导致孔隙体积的变化，从而造成孔隙压力的改变；井孔水位的变化量与岩石中的孔隙压力有关，进而与构造应力有关。因此，利用已有承压井水位连续观测资料，研究区域现今构造应力场状态将是一条既有效又经济的途径[2]。

数字化观测以来，积累了较丰富的井水位观测资料，而利用承压井水位观测资料进行区域地壳应力应变状态的反演，从而进行地震形势的分析研究等，前人已做了较多定性的分析工作[3]：李永善[4]研究了构造应力引起地下水位变化的特征，并推导出了含水层水位变化与其垂直向应力变化之间的关系式；张昭栋等[5-6]利用含水层参数、固体潮效应和气压效应等3 种方法反演含水层应力应变，并进行了分析比较；黄辅琼等[2]利用华北地区40 多口深井水位动态变化资料，定性地分析了大华北地区的构造应力场状态；孙小龙等[7]利用华北地区63 口井的水位变化和井含水层水文地质参数资料，反演出华北地区构造应力场变化图像；而丁风和等[8]和其他研究人员[9-13]利用气压系数和维尼迪科夫潮汐调和分析结果（O1、M2波潮汐因子）等参数，反演得到井孔含水层部分介质参数，进而在水平层状含水层模式下，定量地计算出观测井的垂直向应力场的动态变化过程。

鉴于前人的研究，既有定性的分析也有定量的计算。本文拟借鉴丁风和等[8]定量计算的方法，选择资料可靠性较好的数字化井水位资料，利用气压系数和O1、M2波潮汐因子，反演井孔含水层部分介质参数，定量分析海南省区域构造应力场的时序变化过程，从而分析探寻与地震有关的部分地球物理信息。

1 资料的分析与处理

1.1 观测井基本概况

本文选取了位于琼东北地区的5 口数字水位观测井，其中，海口ZK26 井、向荣村井、火山流体井，这3 口井均沿马袅—铺前断裂带分布，且为活断裂，井距离海岸较近，井水位效应表现为海潮潮汐显著、气压潮汐次之、固体潮汐最弱的复合潮汐形态；文昌潭牛井和琼海加积井，主要沿文昌—琼海—三亚断裂带分布，这2 口井水位潮汐效应表现为以固体潮汐显著、气压潮汐次之、海潮潮汐最弱的复合潮汐形态（图1）。各井承压性均较好，观测资料相对较长且连续稳定。表1 为各井经纬度、井深、含水层岩性、地下水类型、水位埋深等基本情况。

表1 开展研究的琼东北地区5 口井的基本情况Table 1 Basic situation of 5 wells in northeast Hainan

图1 研究的5 口观测井空间位置分布图Fig.1 Spatial distribution of 5 observation wells studied

1.2 数据的整理与处理

本研究从海南省地球物理台网中心数据库下载井孔数字化水位和气压的整点值数据，选取研究资料为2008—2020 年，其中，文昌潭牛井因数据库故障原因，2008—2012 年气压资料缺失，故该井研究起始时间为2013 年。首先将水位、气压两测项数据进行逐月检查，其存在的突跳或错误数据需做删除处理，并按999999 进行缺数标记。水位数据整理成以米为单位的数据文件。

本文研究的5 口水位井地处海南，属于我国沿海地区，其水位除了受固体潮、气压潮的影响，同时还受海潮的影响，表现出复合潮汐效应。井水位受海潮荷载影响的程度，与井孔距海岸的距离大小的关系十分明显，一般距离在5 km 以内时，海潮的影响显著；超过5 km 后，海潮的影响明显减弱。井水位可以分离成363 个潮汐分波，主要有Q1、O1、M1、S1K1、J1、OO1、2N2、N2、M2、L2、S2K2、M3。其中，O1、S1K1分波主要受海潮影响，S2K2分波主要受温度和气压影响，M2分波主要受固体潮影响。琼北地区的海口ZK26 井、向荣村井、火山流体井，这3 口井距离海岸均约5 km，再根据各井水位潮汐分波振幅特征的分析（图2），可看出该3 口井振幅较大的分波主要是O1、S1K1波，其次是S2K2波，因此，主要受海潮潮汐影响，进而选择其中的主要分波O1波进行维尼迪科夫潮汐调和分析得出3 口井的潮汐因子。而琼海加积井和文昌潭牛井，距离海岸约为20 km，由图3 可看出，该2 口井水位的潮汐分波主要是M2波，受固体潮汐影响显著，因此，选择M2分波进行维尼迪科夫潮汐调和分析计算2 口井的潮汐因子。图3 为各井潮汐因子月值曲线。

图2 琼东北地区 5 口井水位潮汐分波振幅特征图Fig.2 Characteristics of tidal wave amplitude of water level of 5 wells in northeast Hainan

图3 琼东北地区 5 口井水位潮汐因子月值曲线图Fig.3 Monthly value curve of water level tidal factor of 5 wells in northeast Hainan

各井气压系数的主要获取步骤为：首先对水位和气压数据进行别尔采夫滤波，其次对滤波后的数据进行一阶差分处理，最后根据差分结果，构建一元回归模型来拟合气压系数[8]。其大小因井而异，变化范围一般为1—10 mm/hPa，个别井特殊时段下气压系数值会出现负值或较大值，本文获取的各井气压系数值基本处于此合理变化范围之内。图4 为各井气压系数月值曲线。

图4 琼东北地区5 口井水位气压系数月值曲线图Fig.4 Monthly value curve of water level pressure coefficient of 5 wells in northeast Hainan

因多数井水位观测资料的趋势性下降变化，直接受区域地下水的开采等影响，且降水量的年周期变化也对井水位观测资料的趋势性变化具有干扰影响，为此需采用去趋势的分析方法，剔除地下水开采和降水补给等长周期干扰对水位的影响[8]。但长周期干扰（趋势性变化）剔除后，降雨和气压的年周期变化对大部分水位的影响依然存在。接着对剔除趋势性变化后的水位进行一般矩平去周期分析，进一步剔除降水和气压对其的影响，然后再进行井水位变化量的计算[8]。地下水开采及降水量的补给，是绝大多数地区多数井孔水位的主要影响因素，因此，本文采用多次剔除这些干扰因素的方法，再进行水位变化量的计算，结果才更可靠合理。图5 为各井水位变化量月值曲线。

图5 琼东北地区5 口井水位变化量月值曲线图Fig.5 Monthly value curve of water level change of 5 wells in northeast Hainan

经过计算得出的潮汐因子、气压系数和含水层参数等，均为月值数据，基本能满足地震的中短期预测对资料间隔的需求。

2 基本理论及参数求取过程

2.1 含水层参数的反演计算

依据丁风和等[8-9]研究结果，在理想的承压含水层模式下，井水位的气压系数和潮汐因子可分别表示为：

上两式联合可得到：

式中，n为含水层的孔隙度，α为固体骨架的体积压缩系数，β为含水层内水的体积压缩系数。ρg为水的重度，在研究过程中取ρg=9.8×103hPa/m。

潮汐因子Bg通过维尼迪科夫（Venedikov）潮汐调和分析获取；气压系数Bp主要利用水位数据和同井观测的气压数据经过别尔采夫滤波和一阶差分，再通过一元回归模型进行拟合后获取；最后经过滑动推算得到含水层孔隙度n与水的体积压缩系数β，然后再利用气压系数Bp和潮汐因子Bg的公式进行推导计算得出固体骨架的体积压缩系数α。

2.2 井水位变化量与含水层垂直向应力变化量的关系

依据前人的研究结果[8-9]，在理想的水平层状含水层（一维）模式下，井孔含水层垂直向的应力变化量Δσz与含水层部分介质参数、井水位变化量之间存在以下定量关系，并且能在一定程度上反映出该区域的构造活动特征，即

式中，Δσz为含水层垂直向应力变化量，E为含水层固体骨架的杨氏模量（E=1/α），ΔH指剔除地下水开采、降雨影响后的含水层应力变化引起的压力水头变化量，即井水位变化量。

从物理意义及异常性质判定来讲，当井孔含水层系统所受应力增强，即Δσz＞0 时，井水位上升，水位埋深值H变小，其变化量ΔH＜0；当井孔含水层系统所受应力减弱，即Δσz＜0 时，井水位下降，水位埋深值H变大，其变化量ΔH＞0[9]。

3 应力变化及地震形势分析

3.1 基于井水位变化量和含水层垂直向应力变化量的定量分析

本文利用琼东北地区5 口观测井的气压系数和O1、M2波潮汐因子等参数，反演得到井孔含水层孔隙度、水的体积压缩系数、含水层固体骨架的体积压缩系数等参数，在理想承压含水层水平层状模式下，定量地计算了2008 年以来5 口观测井的垂直向应力场的动态变化过程（图6）。图6 中各井井水位变化量和含水层垂直向应力变化量的时序月值曲线显示，井孔含水层系统所受应力增强时即Δσz＞0，ΔH＜0。5 口井水位变化量ΔH＜0 时，其相应的各井含水层垂直向应力变化量Δσz＞0 对应地很好。

图6 海口ZK26 井、琼海加积井、向荣村井、火山流体井、文昌潭牛井水位变化量与含水层垂直向应力变化量月值曲线Fig.6 Monthly curve of water level variation and vertical stress variation of aquifer Haikou ZK26 well，Qionghai Jiaji well，Xiangrong village well，Volcanic fluid well and Wenchang Tanniu well

其中，海口ZK26 井、琼海加积井2 口深井的含水层垂直向应力变化量曲线形态相似，且2017 年以来两井—含水层系统应力明显较2017 年前减弱；向荣村井、火山流体井深度相差不大的2 口较深井的含水层垂直向应力变化量曲线形态也相似，井—含水层系统以压应力为主的变化相较海口ZK26 井、琼海加积井明显较弱，且2016 年以来两井—含水层系统应力变化同样呈减弱状态；文昌潭牛井，2018 年以前，井—含水层系统以压应力为主的应力变化量相对集中，而2018 年以后较零散减弱。综合分析，基于数字化水位的井水位变化量与含水层垂直向应力变化量的研究结果表明，琼东北地区近几年来应力变化活动水平相对减弱，需继续跟踪其转折变化等情况的发生。

3.2 基于重力场累积动态及GPS 基线的区域构造活动分析

基于井水位变化量与含水层垂直向应力变化量的应力场时序变化特征研究结果，表明琼东北地区近几年来应力变化活动水平相对较弱，下面将结合本区重力场累积动态变化以及基于GPS 基线观测数据的区域构造活动的多手段综合分析，进一步增强对基于井水位变化量与含水层垂直向应力变化量研究结果的可参考性。

累积动态变化图像，表示相对某一期或某一基准的重力变化，其可突出观测区域不同时期重力场动态演化的相对状态或累积信息。以2016 年9 月这一期为基准，计算2016-09—2017-10（累积1 年）、2016-09—2018-10（累积2 年）、2016-09—2019-08（累积3 年）和2016-09—2020-09（累积4 年）累积重力变化：海南岛陆地区4 年期间的重力场变化相对平稳，总体呈东北部负变化，西南部正变化特征；海南岛陆重力场变化较大的区域分布于琼西南。海口ZK26 井、向荣村井、火山流体井位于海南岛陆北部，琼海加积井、文昌潭牛井位于海南岛陆东部，2016 年以来海南岛陆北部和东部地区重力场累积变化比较平稳，无显著异常（图7）。

图7 海南岛陆重力场累积变化动态分析图像Fig.7 Dynamic analysis image of cumulative change of land gravity field in Hainan Island

GPS 观测数据可以反映地球内部应力场变化在地面的形变响应分布，及构造应力缓慢作用过程中地壳的某些形变和运动特征。因此，我们通过利用GPS 资料来分析研究区域地壳形变特征以及对区域地下流体变化的影响。

目前海南在琼北地区活动断裂带两侧布设有GPS 流动观测站（图8a），以期通过断裂两侧的GPS流动观测点对断裂活动进行形变观测，从而监测区域内火山及地震活动状态。海口ZK26 井、火山流体井、向荣村井、琼海加积井、文昌潭牛井均位于该监测区域周边，根据琼北地区流动GPS 的解算结果（图8b）可知，2008—2020 年以来参考站马鞍岭与各流动观测站之间的6 条基线逐年变化基本在15 mm以内，总体变化较小且稳定，即琼北及周边地区没有明显的拉张和挤压趋势，表明块体内部较稳定。

图8 （a） 琼北地区流动GPS 观测点分布；（b） 马鞍岭参考站与各流动站点基线逐年变化示意图Fig.8 （a） Distribution of mobile GPS observation points in northern Hainan；（b） Schematic diagram of annual change of baseline of Ma’anling reference station and mobile stations

4 结论与讨论

琼东北地区作为海南岛陆主要活动构造带的集中区域，是海南地震监测研究的重点区域。本文借鉴前人的研究结果，利用5 口井的气压系数和各井水位的主要分波O1、M2潮汐因子等，滑动拟合得到不排水状态下，各井部分含水层参数。接着在理想水平层状含水层模式下，利用这些参数与井水位变化量、含水层垂直向应力变化量的关系式，定量分析了琼东北地区构造应力场的时序变化特征。基于数字化水位的井水位变化量与含水层垂直向应力变化量研究结果，再结合区域重力场累积动态变化与基于GPS 基线观测数据构造活动的综合分析，表明琼东北地区近5 年来应力变化活动总体呈减弱状态，需跟踪其转折变化情况的发生。

同时发现琼东北地区井孔性质相当的井，主要表现为井深差不多的井孔，其井水位变化量与含水层垂直向应力变化量曲线形态相似。随着本区域井孔观测数据的积累，将继续对比其他井孔的应力变化定量分析结果，进行更为全面的论证。含水层介质构造应力变化的定量分析，对于有效利用数字化水位资料来甄别地震地球物理异常和机理分析具有重要的意义[8]，因此，本文的研究结果将积极运用到海南省日常的会商及异常核实等分析工作中，为本区域地震形势的分析研判提供有力依据。

与传统的通过现场抽水试验和室内实验等不同，利用数字化水位等资料，结合气压系数和维尼迪科夫潮汐调和分析结果，而获取含水层介质的孔隙度、固体骨架的体积压缩系数和水的体积压缩系数是简便易行的[8]。相较定性的分析，经参与了含水层垂直向应力变化量的分析而计算得到的各井含水层的孔隙度、固体骨架的体积压缩系数和水的体积压缩系数是动态变化的，因此，得到的琼东北地区构造应力场的变化结果是可靠、合理的[8]。另，本文均是在假定理想的水平层状含水层，介质是各向同性、均质的弹性体模型情况下完成的，实际应用过程中还需深入研究。

致谢

在本文的研究过程中，非常感谢内蒙古自治区地震局丁风和高级工程师提供的帮助！