0 引言

承德地区地质构造活动相对平静，2008年以来该地区M_L＞3.0地震仅有4次。2021年6月16日17时57分承德市双桥区发生M_L 2.9地震，6月17日8时7分发生M_L 3.6地震，这2次地震震中位于小营—张营断裂、双峰寺—红石砬断裂附近，震中附近有新建的双峰寺水库，地上载荷的变化会直接影响地壳形变量（孙伶俐等, 2010, 2011；吕琳等，2019），地震是否由水库蓄水引起的地下断层错动诱发的，这引起笔者的思考。地下断层错动产生的地形变化可通过形变观测仪器显示出来。形变观测中的地倾斜量能直接反映地壳介质的微动态变化。常用的观测仪器有宽频带倾斜仪、水管倾斜仪、宽频带石英水平摆倾斜仪等。承德地震台现有形变观测仪器为VP型宽频带倾斜仪，本文利用宽频带倾斜仪8时整点测值，研究承德地区地形变与双峰寺水库蓄水载荷间的相关性。

1 台址、水库概述

承德地震台位于承德市双桥区狮子沟镇殊像寺村，从地质构造看，地处燕山沉陷带、承德台凸中部，南与中生代承德盆地相邻，北靠燕窝铺—周台子背斜。台站所处位置为太古界深变质片麻岩所覆盖。北有NE走向中更新世晚期活动的双塔山—双峰寺断裂，南有NE走向东河套—红石砬第四纪活动断裂，东有控制武烈河流向近SN走向武烈河第四纪活动断裂。承德地震台形变观测山洞岩性为变质片麻岩、砾岩，洞内湿度为75%，洞内温度为12.5℃，高程420 m。自架设以来，承德地震台宽频带倾斜仪观测精度优良，主要干扰因素有标定、调零、进洞人为干扰等。

2019年5月双峰寺水库主体工程建成。水库拦河坝设计最大坝高51.1 m，最大坝长533 m。库岸总体呈弯曲带状，最窄处不足1 km，最宽处达3.75 km（郭建礼等，2010）。总库容1.373亿m³，校核洪水位为395.11 m，设计洪水位为392.50 m。双峰寺水库位于承德地震台东北方向，二者间直线距离约为7 km（图 1）。双峰寺水库地处内蒙古高原向华北平原的过渡带，为燕山北麓低山河谷地貌，地势总体北高南低，区域内岩相复杂多变，岩浆岩、变质岩、沉积岩均有分布，其中，岩浆岩、变质岩分布最广泛，沟谷、河岸多分布松散沉积物。双峰寺水库主要是以城市防洪为主，不同年份水库蓄水变化主要取决于降雨量及防汛要求。

利用双峰寺水库2020年、2021年蓄水量数据与承德地震台VP型宽频带倾斜仪测值，计算二者间的皮尔逊相关系数并建立载荷模型，从而计算由水库蓄水量变化所引起的地倾斜变化。

2 库水位变化与观测资料的一致性分析

双峰寺水库蓄水量为当日早8时统计的数据。将水库蓄水量变化看作一个整体，结合水库形态、水位、蓄水量将水库蓄水量变化近似为长方体。2020年双峰寺水库最小蓄水量为3 700万m³，最大为8 486万m³，可把2020年蓄水变化形状近似于1 000 m×4 300 m×8 m的长方体。2021年双峰寺水库最小蓄水量为1 010.57万m³，最大为4 597.05万m³，把2021年蓄水变化形状近似于2 000 m×3 500 m×7 m的长方体。2020年、2021年水库蓄水量变化如图 2所示，2020年水库蓄水最大变化量出现在6月24日至11月10日，变化量为3 586.48万m³，水位上升8.26 m；2021年出现在6月29日至8月25日，变化量为4 786万m³，水位上升7.11 m。

双峰寺水库2019年建成。为研究VP型宽频带倾斜仪的实时变化，绘制了2014—2021年承德地震台宽频带倾斜仪日均测值曲线（数据已年归零）及一阶差分速度曲线（图 3）。由图 3可见，EW、NS向测值基本上都存在有规律的零漂，且2014—2021年漂移方向大体一致，并没有出现明显的不同，并且零漂一阶差分速度曲线表明漂移速率整体稳定。结合GNSS数据分析河北区域地形变变化特征，研究表明河北区域整体运动较平稳（屈曼等，2021），故笔者认为零漂可能是由地壳运动规律形变所致。

3 线性相关性分析

对承德地震台宽频带倾斜仪整点测值和双峰寺水库蓄水数据进行线性相关分析（王飞等，2019），线性相关性分析主要是计算皮尔逊相关系数，其通常用r表示，取值范围是[-1，1]。计算公式为

$ r=\frac{\sum(x-\bar{x})(y-\bar{y})}{\sqrt{\sum(x-\bar{x})^2 \sum(y-\bar{y})^2}}=\frac{x \text {与} y \text {的协方差}}{x \text {的标准差与} y \text {的标准差的乘积}} $

(1)

其中，|r|＜0.3时代表低度线性相关；0.3≤|r|＜0.5时代表中低度线性相关；0.5≤|r|＜0.8时代表中度线性相关；0.8＜|r|≤1时代表高度线性相关。

由原始数据计算出的皮尔逊相关系数（表 1）可见，水库蓄水量与宽频带倾斜仪8时整点测值间呈中高度线性相关（图 4）。为了验证其相关性，又将2020年、2021年承德地震台宽频带倾斜仪测值去趋势去固体潮后计算皮尔逊相关系数（表 2），结果表明水库蓄水量与宽频带倾斜仪8时整点测值间线性相关很低（图 5）。2020年6—11月水库蓄水量增加，宽频带倾斜仪8时整点测值没有明显变化。2021年6—8月水库蓄水量激增，2021年7月承德地震台宽频带倾斜仪测值漂移量变大，在7月7日、7月28日、8月7日、8月27日对NS向测值调零处理，在7月31日、8月22日对EW向测值调零处理，调零后日均值逐渐恢复至原有水平。笔者认为，2021年7月、8月宽频带倾斜仪两分向测值漂移量大有2方面原因：其一，由于6月承德2次地震震中距仪器较近，虽然震级不大，但宽频带倾斜仪需要一定的稳定期；其二，水库蓄水量激增对仪器产生影响。因此，水库蓄水量的突变与承德地震台宽频带倾斜仪测值间存在线性相关，但相关性很低。

4 载荷模型建立计算 4.1 质点载荷模型

为分析双峰寺水库蓄水载荷变化所引起的近地表应力影响，将仪器布设位置和水库区域作为一个整体，设其为均匀的半无限空间弹性体（李祖宁等，2007；闫伟等，2015）。把双峰寺水库蓄水量的变化设为集中载荷变化，其中，p为垂向正应力；M（x，y，z）为任意点。受力示意图如图 6所示（张英杰等, 2020, 2021）。

该问题最早由Boussinesq求解（胡卫建等，2002）。M（x，y，z）点的竖向正应力σ_z和竖向位移w的表达式为

$ \sigma_z=\frac{3 p}{2 \pi R^2} \cos ^3 \theta=\frac{3 p}{2 \pi} \cdot \frac{z^3}{R^5} $

(2)

$ w=\frac{p(1+\mu)}{2 \pi E}\left[\frac{z^2}{R^3}+(1-2 \mu) \frac{1}{R}\right] $

(3)

其中，E为弹性模量；μ为泊松比；p为作用力；R为M点到集中作用力点的距离

$ R=\sqrt{x^2+y^2+z^2} $

(4)

其中，x、y方向的倾斜量为T_x、T_y

$ T_x=\frac{\partial w}{\partial x}=\frac{\partial w}{\partial R} \cdot \frac{\partial R}{\partial x}=\left[-3 \frac{p(1+\mu)}{2 \pi E} z^2 R^{-4}-\frac{p(1+\mu)}{2 \pi E}(1-2 \mu) R^{-2}\right] R^{-1} x $

(5)

$ T_y=\frac{\partial w}{\partial y}=\frac{\partial w}{\partial R} \cdot \frac{\partial R}{\partial y}=\left[-3 \frac{p(1+\mu)}{2 \pi E} z^2 R^{-4}-\frac{p(1+\mu)}{2 \pi E}(1-2 \mu) R^{-2}\right] R^{-1} y $

(6)

承德地震台VP型宽频带倾斜仪布设山洞的岩体为变质片麻岩、砾岩，取杨氏模量E = 5×10¹⁰ Pa；泊松比μ = 0.15；g = 9.8 m/s²；水的比重为1 t/m³。双峰寺水库与形变山洞间高程差为55 m，二者间直线距离约为7 km。计算得到2020年双峰寺水库蓄水量使宽频带倾斜仪测值产生的倾斜量变化为7.726 5×10^-5″、3.747 3×10^-10 rad；2021年的为5.789 9×10^-5″、2.808 1×10^-10 rad。2年的叠加倾斜量为1.351 6×10^-4″、6.555 4×10^-10 rad。

4.2 二维载荷模型

基于Boussinesq解，将质点模型优化成二维载荷模型。二维载荷模型是将水库载荷区域看成1个长方形，2020年水库蓄水看作1 000 m×4 300 m长方形，2021年看作2 000 m×3 500 m的长方形，对该水库载荷的总作用力p进行散点化，把水库区域等分为n个网格，作用力平均分布在网格上，各散点的作用力记为p_i，则每个网格作用力为

$ p_i=\frac{p}{n} $

(7)

然后，分别求取各散点对某点M的垂向位移W_i和倾斜量(T_xi，T_yi)，再求所有散点p_i对点M作用的矢量和，即可得到该点的垂向位移w和倾斜量(T_x，T_y)。网格个数为n，其他变量间的关系为

$ T_x=\sum T_{x i} ; T_y=\sum T_{y i} $

(8)

设x、y方向为水库近似长方形的长、宽方向，而非宽频带倾斜仪的NS、EW方向。将水库载荷按照相同尺度划分网格区域，利用式（3）、（5）、（6）计算出每个网格区域产生的倾斜变化量T_xi、T_yi，再利用式（8）求取M点的倾斜量。

计算得到2020年双峰寺水库蓄水量使宽频带倾斜仪测值产生的倾斜量变化为7.689 2×10^-5″、3.729 3×10^-10 rad；2021年的为5.701 8×10^-5″、2.765 4×10^-10 rad。2年的叠加倾斜量为1.339 1×10^-4″、6.494 6×10^-10 rad。与质点模型相比较，结果相差不大。

4.3 三维载荷模型

三维载荷模型是将水库载荷区域看成1个长方体，2020年水库蓄水近似于1 000 m×4 300 m×8 m长方体，2021年近似于2 000 m×3 500 m×7 m长方体。同样，在三维载荷模型中，设x、y方向为水库近似长方形的长、宽方向，而非宽频带倾斜仪的NS、EW方向。三维载荷模型是将网格划分为n个等体积的小长方体，利用式（3）、（5）、（6）分别计算各散点载荷对M点的影响，最后利用式（8）矢量化合成整个不规则载荷对该点的影响。

计算得到2020年双峰寺水库蓄水量使宽频带倾斜仪测值产生的倾斜量变化为7.609 0×10^-5″、3.690 4×10^-10 rad；2021年的为5.701 9×10^-5″、2.765 4×10^-10 rad。2年叠加为1.331 1×10^-4″、6.455 8×10^-10 rad。

通过对比3个模型的结果（表 3）可知，三者间相差不大，二维、三维载荷模型相对偏小，其收敛性也说明结果的真实可靠。形变学科异常核实中，邱泽华(2004a, 2004b)研究认为，大于3×10^-8 rad的倾斜量会对仪器产生影响。因此，认为双峰寺水库蓄水量变化对承德地震台VP型宽频带倾斜仪测值影响不大。

5 结论与讨论

通过计算双峰寺水库蓄水量与承德地震台VP型宽频带倾斜仪8时整点测值去趋势处理后数据间的皮尔逊相关系数可见，二者间相关性不大。

通过建立质点、二维、三维载荷模型计算载荷变化所引起的地形变量发现，3个结果间相差不大，且二维、三维模型结果趋向收敛，亦佐证了载荷模型计算的可靠性，倾斜量计算结果小于3×10^-8 rad，且远小于日变量和8时整点测值的变化。由此认为，双峰寺水库蓄水量对承德地震台VP型宽频带倾斜仪测值变化影响较小，不足以影响仪器的日常曲线变化趋势，进而认为2021年6月承德市2次地震为水库蓄水诱发地震的可能性较小。

因岩石岩性和参数均为估算得出，模型中假设均匀半空间弹性体，库区形状假设为长方体，但其实际形状变化较大等，所以，这些条件均会影响到计算结果，使结果可能存在一定的误差。后续可结合双峰寺水库附近测震类仪器产出数据进行更深入的分析。