泰安地震台洞体应变与钻孔体应变参数解算结果对比分析

孟建国 蒋其峰 卢双苓 李惠玲 杨彬 周均太

1）泰安基准地震台，山东省泰安市罗汉崖路2号 271000

2）山东省地震局，济南 250014

3）邹城地震台，山东济宁 273500

0 引言

将应变连续观测数据应用于地震活动预测是自邢台地震后在李四光倡导和组织下开展起来的。国内对洞体应变资料的研究始于20世纪80年代，经过近30年的发展，洞体应变理论研究基础不断完善。目前，洞体应变伸缩仪在我国很多台站都有安装，分2测向和3测向2种，得到的资料不仅可用于各方向的线应变分析，还可经换算后用于主应变及剪切应变的分析。对地表附近水平面内的面应变变化进行监测的测项除了洞体应变，还有钻孔应变。两者具有很强的可比性。同时，区域内地下物质的应变积累与释放直接影响了地下物质密度的改变，可以通过重力观测进行验证（马栋等，2013）。本文研究的内容及方法参考了杜慧君等（1993）用Venedikov方法对应变观测资料的调和分析，计算结果表明，应变观测 NS向固体潮幅度及均方差明显优于EW向。本研究中还应用了蒋骏等（1993、1994）根据应变花理论（尹祥础，1991）对多方向潮汐线应变组合观测确定平面应变状态及提取潮汐面应变、体应变、平面剪切应变等信息的方法和公式。张雁滨等（1996、1997）对中国部分潮汐应变观测点的理论应变状态和观测应变状态的计算为本文的研究提供了很好的理论依据。

1 观测点概况

泰安基准地震台地处泰山南麓、莱芜弧形断裂带的北侧附近，台址基岩为花岗片麻岩，岩体完整，致密均匀，地脉动水平低，波导性能良好，测量信噪比高。台站形变观测仪器和测震观测仪器位于台站的专用山洞内。山洞洞室覆盖厚度29m以上，总进深130m，其中主洞进深76m，高差33m，避免了洞室内空气流动对仪器观测造成的影响，室温年变幅约0.6℃。泰安台的洞体应变观测仪器为SS-Y短基线铟钢伸缩仪，属首批“九五”数字化试验仪器。观测试验运行始于1998年2月，1999年通过验收并正式运行，于2008年12月进行了“十五”仪器升级改造。其中伸缩仪改造升级部分为标定器、传感器以及数据采集器，原铟钢基线未做改动。伸缩仪布设在原水管仪的基线墩上，与水管基线并行。伸缩仪基线EW向长10.00m，NS向长 30.03m；方位角 EW 向为 124°37′，NS向为 39°48′。2008年 9月，在山洞旁边距离伸缩仪直线距离10m处，安装了 TJ-2型体积式应变仪，探头底部的实际埋深为74.80m。探头与地层耦合采用水泥固结方式。

2 观测资料评价

本文采用伸缩仪的整点值数据进行计算，数据处理范围从2009年1月～2014年6月。

调和分析方法采用Venedikov调和分析公式（张雁滨等，2000），固体潮的表达为各种谐波的综和。即

式中，y（t）为t时的固体潮观测值；Hn为某潮汐波的观测振幅；ωn为某潮汐波的角频率；φn为某潮汐波的初位相；y0（t）为固体潮在t时刻的零点漂移值。Venedikov调和分析是把日波、半日波群分开。分离出不同的日波波群和半日波群，然后分别求出它们各自的振幅比和相位差以及每个潮汐波计算结果的均方差和总均方差。

公式（1）可改写成

式中，tj为自Tj起算的以小时表示的时间；Tj为tj为零时的时间；φn（T）j为某潮汐波在Tj时初相位。

根据滤波理论，利用偶数滤波C求出与Hncosφn（Ti）有关的数Mi，利用奇数滤波S求出与 Hnsinφn（Ti）有关的数 N，而后得到未知数为 xk，yk的线性方程组

其中，τ＝1，2；I＝1，2，…，m；k＝1，2，…p；m ＞p

上式中未知数的个数为2p个，方程的个数为2m个，ak，βk为第k波群的起止波号。每连续48h的观测数据，经调和分析后，对日波和半日波群分别得出2个方程。最后用最小二乘法求解以上线性方程组，得到潮汐因子δk和相位差σk

通过MATLAB软件编程，将2013年全年洞体应变2分项整点值及体应变全年整点值带入运算方程，分别求出各参数项。表1、表2分别为洞体应变的2分项调和分析结果。

计算结果表明，泰安台伸缩仪NS向潮汐因子保持在0.7001、EW向保持在0.3008范围内，应变观测NS向固体潮幅度及均方差明显优于EW向。在张渤带及其邻区洞体应变仪测项潮汐因子及其误差（马栋等，2013）的对比分析中，36个台数据统计发现72%的台站潮汐因子取值为0～1和1～2的占27%、2～3的占11%；充分说明了泰安台仪器响应灵敏度、洞室条件等因素都符合形变台站仪器运行管理规范所需要求。

表1 洞体应变NS项应变调和分析表

表2 洞体应变EW项应变调和分析表

3 应变参数提取及对比

目前，泰安台安装的洞体应变为2个分量的伸缩仪，研究人员将通过构建计算模型提取应变组合观测信息（蒋骏等，1993；张雁滨等，1997），进而对应变参数进行提取解算。为分析和验证该方法的可信度，我们将通过实际观测资料与解算结果进行对比。首先根据实际观测的伸缩仪资料将2个分量潮汐线应变观测的方位角定为α1、α2，对地表潮汐实际的剪切应变eθλ可用理论值进行模拟响应系数（潮汐因子）δk由2个不同方向的潮汐线应变实测资料用调和分析拟合获得，已由式（5）得出。将泰安台伸缩仪观测NS、EW向2个不同方位（方位角分别为 α1＝39.8000°、α2＝124.6167°）组成方程式（6）。

根据式（6）中两 2个布设潮汐线应变观测值 e11、e22，（因作为已知，）即可求得另外 2个未知数 eθθ、eλλ（刘序俨，1994；蒋骏等，1993；张雁滨等，1996）。

通过求得的 3个已知量 eθθ、eλλ及 eθλ进而可以通过式（7）（8）计算得到其它的应变矢量

求得最大主应变e1、最小主应变e2。最大剪应变τ和最大主应变方位角α依式（9）、（10）求得

通过MATLAB软件编程，导入观测数据，其中伸缩仪2分项数值提取为2013年全年日均值。通过运算得出最大主应变e1、最小主应变e2、最大剪应变τ和最大主应变方位角α，并分别示于图1～4。

通过对比可清晰看到，主应力及剪应力自年初开始呈正弦波趋势转变，应力开始不断积累，应力方位角一致，到105天（3月15日）附近到达应力积累顶点，后又不断进行释放，一直到220天（8月10日）附近应力释放完全，最大主应力恢复到年初水准；而最小主应力、最大剪应力略高于年初水平，此时最大主应力方位角突然变为反向，接续了应力的积累到释放的过程。

值得注意的是，在第195天（7月14日）附近时，主应力及剪应力出现了4天左右的相对迟滞，随后应力加速释放。从最大主应变方位角时序图中可以明显看出，195～202天（7月14日～7月21日）附近的方位角出现了一次停滞，这应该是某一个不明作用力对原年变趋势的一个附加作用。从方位角及主应力变化形态上判断，该作用力为压应力。从第195天（7月14日）一直到2013年年底，整个洞体应变拟合后的最大主应变一直维持在呈相对固定斜率的加速积累状态，积累量及变化速率均高于上半年；最小主应变的变化趋势则表现为相对迟滞，变化量及变化速率明显低于上半年；而与之对应的主应变方位角则变化较为突出，这可能是新应变作用力与原正常应变之间相互作用的结果（牛安福等，2003）。

图1 泰安台洞体应变拟合最大主应变时序变化

图2 泰安台洞体应变拟合最小主应变时序变化

为了验证这一结论，我们比对了同样反映地表应变量变化的体应变数据。首先提取体应变2013年全年整时值，通过多项式曲线拟合对数据进行处理。相关步骤为：将观测序列划分为p个数据段，每段用m次多项式拟合。约束条件为：相邻数据段中拟合曲线在共同处的r阶导数相等。设所要求的曲线方程为

通过MATLAB软件实现公式的解算。将数据一次带入方程，运算得出体应变拟合数据并绘成图5。

图3 泰安台洞体应变拟合最大剪应变时序变化

图4 泰安台洞体应变拟合最大主应变方位角（弧度）时序变化

由图5可知，体应变对地表地应力观测值在170天（6月18日）左右因年变而开始逐渐反向恢复，但在第195天（7月14日）附近时突然加速恢复，之后地应力加速积累。该应力在体应变记录上表现为压应力。截止到205天（7月25日）附近，通过计算发现，应力变化量速率达到年变量速率的2倍。从250天（9月8日）起至2013年年底，拟合值变化速率放缓，明显低于上半年同期变化量。上述结果清晰证明了体应变应力的加速积累与洞体应变拟合解算结果是一致的。

综上所述，在2013年的195天（7月14日）起，泰安台所处区域出现一力源尚不明了的应力变化，导致地表观测应变量高于正常值。该应力为压应力，造成地表岩石挤压，其日作用力应变量为正常值的2倍。由于该应变作用力一直与原背景应力之间相互作用，放大了最大主应变的变化量，对最小主应变则起到抑制作用，同时也造成了最大主应变方位角的不断改变。

图5 体应变多项式曲线拟合时序变化

4 与重力观测结果的对比

通过对洞体应变2测向的拟合解算，并分析提取的相关参数后发现了自2013年第195天（7月14日）到年底出现了不明力源的应力加速积累，导致了地表应变呈加速受压状态。因地下物质受压会导致重力加速度的转变，结合泰安台PET重力仪监测数据，我们同样采取拟合法提取了重力参数，经分析发现，该时间段内重力加速度也发生相应的突然转变（见图6中的第200天（7月19日）附近出现的变化）。随后重力加速度出现了异常（如图6所示）。重力数据的整个变化趋势直接对应了洞体应变、体应变拟合提取的解算结果，这也验证了这一不明力源的作用。

图6 重力观测K-L法最佳直线拟合时序变化

5 结论

对洞体应变观测数据进行Venedikov调和分析计算，给出了洞体应变固体潮观测资料的维氏调和分析结果。计算结果表明应变观测NS向固体潮幅度及均方差明显优于EW向。泰安台洞体应变观测响应灵敏度、洞室条件等因素都符合前兆台站仪器运行管理规范所需要求。

通过对泰安台2个测向潮汐线应变组合观测进行了平面应变状态建模，提取了最大主应变、最小主应变、最大剪切应变及最大主应变方位角等信息，在对比研究中发现了该区域在2013年6月后出现了应变场地应力的异常，通过对比同样反映地表应变场变化的体应变拟合值，发现该应变场异常为客观存在，并且表现突出。综合数据判断：2013年自195天（7月14日）起，泰安台所处区域，出现一不明原因应变力，造成了地表一定的应变量高于正常值。主要表象为压应力，造成地表岩石挤压。日作用力应变量为正常应变量的2倍。该应变作用力一直与原正常形应变之间相互作用，放大了最大主应变的变化量，对最小主应变则起到抑制作用，同时也造成了最大主应变方位角的不断改变。通过分析应力场异常与重力场异常的对应关系，也验证了这一不明作用场的变化。