0 引言

SS-Y型伸缩仪是中国自行研制的适用于观测地壳应变和固体潮水平分量连续变化的仪器，可为研究地震孕育过程中水平应变的变化规律提供数据(胡卫建等，2002)。目前，在日常分析及震情会商中虽然提出大量异常现象，但极少得到地震“支持”(马栋等，2011)。异常无震原因比较复杂(邱泽华，2004)，加之对地震孕育与发生的各种因素及地震前兆异常信息知之甚少，尤其是在前兆异常信息核实过程中，各种观测环境干扰及观测技术不完善等问题的存在，对识别与提取地震前兆信息造成更多困惑。因此，如何利用有效的数学物理方法进行核实工作的论证(牛安福等，2011)，便成为识别地震前兆异常信息的首要问题。

易县地震台(以下称易县台)为形变综合台站(国家基本台)，位于太行山隆起带与华北平原拗陷区分界部位，台址地处太行山北段低山丘陵地区，紫荆关断裂带东侧，位于京津冀交汇地区，是首都“外八线”工程重要改造台站。台基岩性为中生代侵入闪长岩花岗基底岩与震旦系沉积岩和闪长岩小侵入体接合部，山洞岩石构造不均一。现配备SS-Y型伸缩仪、DSQ型水管倾斜仪、VP型垂直摆倾斜仪、TJ-Ⅱ型钻孔体应变仪器、YRY-4型分量钻孔应变仪、台站跨断层短水准6种形变观测手段，曾连续5年获得全国形变综合观测资料评比前3名，其中SS-Y型伸缩仪自2001年数字化仪器改造后，数据连续、可靠，符合规范(中国地震局，2004)要求。

2012年8月起，易县台伸缩仪NE向数据由压缩状态转为拉张状态，异常核实发现，2012年起该台观测山洞周边采石场开采频繁且开采量大，荷载变化是否为引起异常的主要原因？在此选取发生转折变化的2012年伸缩仪NE向观测数据进行分析，利用Ansys软件Solid45单元进行建模，核实荷载变化对观测数据的干扰。

1 异常核实

SS-Y型伸缩仪用于精密测量地壳表面两点间水平相对变化(杨婕等，2009)，易县台于2001年6月开始采用该仪器观测地壳应变固体潮水平分量，2013年11月29日— 30日更换仪器，将原电涡流传感器更换为差动变压传感器，测量位移长度变短。

自2006年2月以来，易县台伸缩仪NE向观测数据处于下降趋势，2006年7月4日文安M_L 5.1地震后出现转折，2007年2月至2008年9月底，持续下降幅度近4×10-6，2008年汶川地震后出现转折，2009年10月至2012年8月观测曲线仍呈下降趋势，但速率有所减小，2012年8月起出现转折，转为上升变化(图 1)。

对伸缩仪各项性能、场地环境干扰、供电设施运行及周围观测环境调查发现：①伸缩仪输出数据与数采面板数据一致；②固体潮曲线记录清晰，无人为干扰，气温、气压均季节性周期变化；③观测场地环境(采石场、抽水井等)发生变化，其中：1^#采石场位于观测山洞以东1.58 km，2012年开始进行大规模开采，主要采取爆破形式。2^#采石场位于山洞西南7 km处，2000年至今一直进行开采。3^#采石场位于山洞东南6.67 km处，2012年至今一直进行开采；抽水井在易县台西南0.5 km处，自2012年1月起每日6:30抽水，约3.5小时，抽水量约80 m³。因抽水干扰主要对体应变、分量应变数据影响较大，对洞体应变其他观测手段影响较小，基本排除抽水干扰。

2 荷载模型建立

易县台伸缩仪NE向观测数据2012年起由下降趋势转为上升变化，即由压性变化转为张性变化。核实发现，2012年易县台观测洞室(伸缩仪NE向)周围3个采石场每天有20 — 30辆翻斗车运送石头，每年约开采6个月，采石量约1.08×10⁴ t，伸缩仪NE向观测数据异常是否与采石场造成的荷载变化有关，在此建立荷载模型进行论证。

2.1 有限元模型Solid45单元建立

有限元法的基本思想是，把复杂形体拆分为若干形状简单的单元，利用单元节点变量，对单元内部变量进行插值，实现对总体结构的分析，将连续体进行离散化，即称网格划分。

采用层状地基模型，依据地质勘察资料及土工常规试验测定的土性指标，将地基土分为3层，每层厚度与实际厚度一致，在水平面上的长度和宽度取塔底换算半径的5倍，该塔横截面为正方形，换算半径在此取为1/2边长，并约束地基边界面的各自由度，竖向深度取地基土的勘察深度。

(1)数值计算模型。计算单元采用空间8结点六面体单元(图 2)，亦可退化为无中间结点的空间四面体单元和三棱柱单元，各结点沿其坐标x、y、z共3个平移自由度，分布式荷载可作用于单元的各个侧面，采用该单元可分析大变形、大应变、塑性和屈服等问题。

若假定单元材料为各向同性，则其弹性刚度矩阵为

$ \left[ D \right]=\frac{E}{(1+\mu)(1-2\mu)}\left[ \begin{matrix} 1-\mu & v & v & {} & {} & {} \\ v & 1-\mu & v & {} & {} & {} \\ v & v & 1-\mu & {} & {} & {} \\ 0 & 0 & 0 & \frac{1-2\mu }{2} & {} & {} \\ 0 & 0 & 0 & 0 & \frac{1-2\mu }{2} & {} \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\mu }{2} \\ \end{matrix} \right] $

(1)

式中，E、µ为材料的弹性模量和泊松比。弹性刚度矩阵可由塑料理论推导得出。

(2) Solid45单元构造三维实体结构。采用Solid45单元用于三维实体结构构造，通过8个节点来定义，各节点有3个沿着x、y、z方向平移的自由度，单元具有塑性、蠕变、膨胀、应力强化、大变形和大应变能力，用于沙漏控制的缩减积分选项。Solid45单元几何描述见图 3。

2.2 模型建立

易县台伸缩仪NE向数据2012年出现转折上升趋势，结合该台观测山洞周围环境，清晰可见：1^#采石场2012年开始开采，开采量大，距伸缩仪所在洞室最近，对其NE向数据影响应最大，2^#、3^#采石场距山洞较远，开采量较少，对伸缩仪荷载影响量级较小，因此利用Ansys中的Solid45单元对1^#采石场荷载进行建模。根据实际情况，建模时需注意：①开采的石头基本为花岗岩，密度、成分并非完全一致，但开采用料来源于一个山体，相差不多，由此设定开采的石料为均匀、各向同性的弹性物质；②采石位置主要是山体，主要采用铲车开采，开采面并非规则形状，建模时将开采面作规则锥形。

根据实际情况，对1^#采石场进行建模，建立一个底面积为50 m²、高为30 m的锥形，假设开采的岩石为均匀、各向同性的花岗岩，选取其弹性模量为1.5×109 g/cm³，泊松比为0.15，密度为2.6 g/cm³，进行模拟，计算应力大小。

3 模拟计算

通过建立的采石场三维模型进行模拟计算，得到1^#采石场采石量载荷变化对于整个建模区域等效应变的三维等值线，具体结果见图 4，其中：模型输出结果为：STEP = 1，SUB = 1，TIM E= 1，RSYS = 1，SXY产出的应力为均匀的。

由图 4可知，1^#采石场锥形模型应力最大点主要体现在顶部，所产生的应力最大值为392 751；最小应力主要在锥形内部，所产生的应力为-0.436×10^-7。通过应力与应变公式计算应变值，即

$ \sigma =E\times \varepsilon $

(2)

式中，σ为应力，E为弹性模量，ε为应变。

将花岗岩弹性模量E = 1.5×10⁹ g/cm³代入式(2)，得到采石场2012年开采造成的最大主应变应为ε = 2.62×10^-7，伸缩仪NE向2012年数据曲线的年变幅为4.08×10^-7，二者之间对比可知，NE向实际应变幅度是模拟计算结果的约2倍。由此可知，易县台NE向2012年起出现的转折上升变化，应该与台站周边采石场开采造成的荷载变化有关。

4 结论

通过荷载变化对伸缩仪观测产生的实际影响量进行计算，可知：①从时间上判定荷载对地应变产生的影响及应变量值；②因地壳结构复杂，荷载距应变观测场地较近、变化量较大时应变量值较为准确。

综上所述，可知：易县台2012年伸缩仪NE向应变年变幅达4.08×10^-7，是1^#采石场荷载模型模拟计算最大应变值2.62×10^-7的2倍。由于只选取开采量较大的1^#采石场建模计算，若增加2^#、3^#采石场，那么荷载对伸缩仪的影响将更大。因此，易县台伸缩仪NE向分量2012年张性转折变化与采石场载荷卸载有较大相关性。

本次建模未确定具体采石量，且未考虑伸缩仪所在位置与采石场具体方位，只考虑主应力大小，未考虑剪应力，所得结果不够精确，需要进一步搜集准确资料，进行荷载建模。