刘序俨，方宏芳，张凯，张翊宁，马超文，黄声明，季颖锋

（1.福建省地震局，福建 福州 350003；2.福建省漳州市地震局，福建 漳州 363000；3.福建省地震局漳州地震台，福建 漳州 363000）

漳州台深井应变不变量同震响应分析

刘序俨1，方宏芳2，张凯3，张翊宁3，马超文2，黄声明1，季颖锋1

（1.福建省地震局，福建 福州 350003；2.福建省漳州市地震局，福建 漳州 363000；3.福建省地震局漳州地震台，福建 漳州 363000）

在解决钻孔应变观测加衬模型的待定常数的基础上，巧妙地利用不变量与坐标系选择无关的性质，推导出了5组观测组合计算不变量的简洁表达式，并以玉树地震为例，选取漳州台深井观测系统的探头上、下端面8分量应变观测单元的观测资料，采用加衬模型分别根据探头上下端面的5组观测组合公式计算了地震当天的最大主应变、最小主应变、最大主应变相应的主方向、最大剪应变以及面应变与体应变分钟值，分别给出了探头上端面234观测组合在2010年玉树7.1级地震发震当天7h30m～9h30m时间窗内的上述不变量分钟值同震响应图像，最后对同震响应性质进行了分析，并对最大、最小主应变以及最大主应变相应的主方向的可靠性进行了讨论。

多分量应变观测；加衬模型；应变不变量；同震响应；可靠性分析

引言

鉴于钻孔处的应变量不等于原有应变量，此时不能采用仅适用于各向同性均匀弹性介质中的摩尔圆定理，而要改用加衬模型来计算应变不变量。如何确定加衬模型的待定常数，如何根据加衬模型给出根据探头上下端面的5组观测组合计算不变量的表达式，如何对应变观测数据进行预处理，面应变与体应变又如何计算，应变不变量同震响应图像的实质是什么，又如何对应变不变量进行可靠性分析？以上问题正是本文要探讨的主要内容。

1 台点及仪器概况

漳州地震台位于龙海市九湖镇院后村石狮岩，地理坐标：E 117°37′54.9″，N 24°28′58″，海拔42.3 m。台站所在地为亚热带水果花卉种植区，年平均温度为22℃，年平均湿度为70%，年平均降雨量为1 400 mm，年平均雷电天气约20多天，海浪干扰轻度。台站位于著名的北东向长乐-诏安断裂带的南端，台站基岩为花岗闪长岩，覆盖层厚度为15～32 m。台站建有一条长100 m、走向近EW向的窿道，窿道主巷道总体走向为SWW266°，支巷为NNW356°，从巷道交叉口以西安装SS-Y型伸缩仪EW向，仪器长度24.9 m，支巷安装伸缩仪NS向仪器，长度17.8 m。

2008年12月，在台站一口253 m的钻井内安装了一套 “RZB-3D型地壳形变宽频带综合观测系统”，探头上下端面分别确均匀布设了 “米”字形的四分量水平应变测量单元，除此以外，还布设了系统参考元件、钻孔倾斜测量单元、应变地震波测量单元、精密地温测量单元及水位、气压与井温等辅助测量单元，共有17个观测量，其中水平应变测量单元相邻两条测线相隔45°。钻井的直径为Φ130 mm，井口以Φ146 mm套管下至孔口以下14.80 m处。整孔深度达到253 m，静止水位21.10 m，探头安装在钻孔247.30至251.83 m之间。

2 应变不变量表达式

在地球自由表面，任一质点的应变张量矩阵为一个2阶对称矩阵，矩阵中的元素仅有3个独立分量，因此，需要在通过该质点的3个不同方向进行应变观测[1]，才能确定该质点处的应变张量矩阵，由于在不同的正交应变曲线坐标系中，各坐标系中的应变张量矩阵都是相似矩阵，由不同坐标系下的应变张量矩阵皆能给出该质点的主应变及其相应的主方向[2]。

既然主应变及其主方向为该质点的几何不变量，与坐标系的选择无关，为方便起见，本文对上端面的2#测线方向作为Y轴，与2#测线成正交方向的4#测线方向作为X轴，探头端面的中心作为原点O，设εi为i#测线的线应变观测值 （i=2，3，4，1），设εmax与εmin（εmax＞εmin）为探头端面的主应变，设θ为εmax在该坐标系中的 “方位角”，即 “主方向”，且最大与最小主应变的主方向是相互正交的。表1为漳州台钻孔应变观测系统上下端面水平应变测量单元参数表。

对于4条测线，任取3个方向作为一组，则有4种组合， 1#、2#、3#、4#本身也可以形成一个组合，则共有五种组合。每种组合皆可计算出相应组合的εmax、εmin与θ。

对于钻孔应变，因钻孔处的应变值不等于原有的应变值，此时不能采用摩尔圆定理以不变量 εmax、 εmin与 θ来描述 αi方向上的应变量 εi，因为摩尔圆定理仅适用于钻孔的各向同性的均匀弹性介质。对于钻孔应变观测则要采取所谓的加衬模型[3、4]

式中a与b为待定常数，它们与套管材料和岩石特性及其两者的耦合有关。εmax、εmin与θ与之所以称为应变不变量，是指这些量与坐标系的选择无关，按弹性力学，通常在剪应变为零的平面上，在相互正交的两个方向上线应变会取得最大值和最小值，这两个值称为主应变，与主应变相应的方向称为主方向。根据弹性力学，在各自同性均匀弹性介质 （未开孔）中，由εmax、εmin可以计算面应变s[5、6]，其公式为

表1 RZB-3D型地壳形变宽频带综合观测系统上下端面水平应变测Table 1 Azimuth of level strain observations on upper and lower surfaces of the detector in the borehole strain observation system

3 玉树地震深井应变观测资料的预处理

3.1 加衬模型的建立

如要采用加衬模型计算不变量εmax、εmin与θ，首先要确定a与b值，而欲确定a与b，则必须先要给出未开孔时的应变固体潮不变量理论值，但在有关固体潮文献[7～9]中仅给出了应变固体潮理论值公式。为此，本文假设（1）式中的a与b均等于1/2，此时（1）式回归到了摩尔圆定理[5、6、10]并以玉树地震前1个月，即2010年3月份的应变固体潮整点理论值替换观测值，则由（7）～（11）式分别可计算出5组观测组合的应变固体潮不变量的理论值。计算结果表明，5组不变量几乎完全相同，仅在小数点后2位存在计算误差，这是预料之中的，因为不管应变固体潮理论值，还是摩尔圆定理，它们都是根据弹性力学理论推导出来的。于是，不但证明了（7）～（11）式的正确性，也等价地证实了不变量的唯一性，因此，无论从哪种观测组合都可以获得不变量的理论值。

鉴于应变为一相对值，我们选择2010年3月1日00h作为起始时刻，把应变观测值与相应的不变量理论值分别减去此起始时刻的相应值，这样就分别取得了相对应变值εi与模型拟合值然后根据 3月份的整点观测值εi， 在的前提下，根据最小二乘法就可得到a与b的估计值。

3.2 玉树地震同震响应的不变量分钟值计算

2010年4月14日7 时49分，我国青海玉树地区发生7.1级强烈地震，地震造成严重人员伤亡和重大经济损失，距离震中2 300 km之外的漳州深井综合观测系统清晰的记录到了同震响应波形。

在本文3.1节，已利用玉树地震前1个月，即2010年3月份的应变观测资料确定了加衬模型的待定常数a与b，之所以这样做，是考虑到在该月份漳州深井应变系统未记录到较大的同震响应，可以避免同震响应对确定a与b的影响。我们可以把该月份计算得到的a与b代入（7）～（11）式，就可以分别获得上下端面的5组观测组合的不变量曲线图。这里要指出，在将a与b值代入（7）～（11）式之前，首先要对上下端面的4条测线的应变观测值作零漂消除，在消除零漂后，还要作基线改正，即选择2010年4月1日00h作为起始时刻，，把此时刻之后的应变观测值分别减去起始时刻的相应观测值，就可以获得该时刻相对于起始时刻的相对应变值，因为应变观测本身就是一种相对观测。

4 玉树地震同震响应分析

图1 上端面1#、2#、3#、4#测线2010年4月14日线应变分钟值变化曲线图Fig.1 The minute curve of strain variation in line on the upper surface 1#,2#,3#,4#at test line at April 14,2010

4.1 应变同震响应图像

为节省篇幅，本文仅给出上端面1#、2#、3#、4#测线2010年4月14日线应变分钟值变化曲线图。图2为1#测线在地震当日7：30～9：30时间窗的线应变分钟值变化曲线图。

图2 上端面1#测线2010年4月14日7：30～9：30线应变分钟值变化曲线Fig.2 The minute curve of strain variation in line on the upper surface 1#during 7:30 and 9:30 at April 14,2010.

从图1与图2可以看出，由于两者的时间窗长度不一样，时间窗短的同震响应细节比时间长的可以得到更充分的显示。因此，在后文所采用时间窗为玉树地震当天7：30-9：30这2个小时，因为玉树地震波到达漳州台大约在8 h左右。

4.2 不变量同震响应图像

采用公式（5）～（9）分别进行上、下端面五种组合的最大主应变、最小主应变、最大剪应变和最大主应变方位角的计算。鉴于篇幅，在这里仅以上端面2#3#4#组合为例给出地震当日7：30-9：30时间窗内上述不变量图，分别见图3～6。

4.3 面应变、体应变同震响应图像

随着深井应变观测技术的进展，我国许多台站都同时安装有均匀分布的4分量应变测量单元，可以形成两组相互正交的观测组合，由 （4）式可根据每组相互正交的观测的应变观测值之和求得面应变，在近地表处，体应变等于该点处的面应变的2/3。为节省篇幅，本文仅给出上端面其中一种正交组合的面应变与体应变同震响应曲线图。

图3 上端面2#3#4#组合2010年4月14日7：30～9：30最大主应变同震响应曲线图Fig.3 The co-seismic response curve of maximum principal strain on the upper surface 2#,3#,4#between 7:30 and 9:30 at April 14th,2010

图4 上端面2#3#4#组合2010年4月14日7：30～9：30最小主应变同震响应曲线Fig.4 The co-seismic response curve of minimum principal strain on the upper surface 2#,3#,4#between 7:30 and 9:30 at April 14,2010

图5 上端面2#3#4#组合2010年4月14日7：30～9：30最大剪应变同震响应曲线Fig.5 The co-seismic response curve of maximum shear strain on the upper surface 2#,3#,4#between 7:30 and 9:30 at April 14,2010

图6 上端面2#3#4#组合2010年4月14日7：30～9：30最大主应变方向角同震响应曲线Fig.6 The co-seismic response curve of maximum principal strain azimuth on the upper surface 2#,3#,4#between 7:30 and 9:30 at April 14,2010

如图7～8分别为地震当日上端面1#3#组合在7：30～9：30时间窗内面与体应变同震响应变化曲线图。

图7 上端面1#3#组合2010年4月14日7：30～9：30面应变同震响应曲线10-8Fig.7 The co-seismic response curve(10-8)of surface strain on the upper surface of 1#,3#,between 7:30 and 9:30 at April 14,2010

图8 上端面1#3#组合2010年4月14日7：30～9：30体应变同震响应曲线 10-8Fig.8 The co-seismic response curve(10-8)of bulk strain on the upper surface of 1#,3#,between 7:30 and 9:30 at April 14,2010

从同震响应曲线可以发现，在玉树地震波到达漳州台深井观测系统之前，上下端面的4条测线的线、面、体应变、最大与最小主应变及其主方向的分钟值观测曲线都是非常平稳的，这是因为上述诸量的应变固体潮主要是由低频日波和半日波组成的。而地震波到达后所记录到的变化是强震所引起的高频应变波，因此由玉树地震所引起的同震响应是一种叠加在低频应变固体潮曲线上的高频变化，也正是深井观测系统对地震的同震应变响应，所记录到的应变响应波形以一种振荡衰减的形式出现，而地震过后，曲线又恢复到这些量的应变固体潮应有的缓慢低频变化状态。

因限于篇幅，本文虽然仅给出了上端面其中一组观测组合的不变量响应曲线图，但笔者却对上下端面其余几种观测组合的不变量响应进行了计算。结果发现，除主方向响应出现极大差异外，其余几组组合的同一种不变量彼此之间虽然也有所差异，但远比主方向要小，且相似性也很高。除了对玉树地震同震响应进行计算外，笔者也曾用1个月的观测资料计算了上下端面共10种观测组合的不变量，且给出了他们的整点值曲线图，发现最大与最小主应变曲线大体相似，数值也很接近，对他们一个月观测资料的整点值进行两两相关分析，其相关系数大多在0.7～0.8以上，在自由度f高达700以上的情况下，可以认为两者几乎是完全相关的，对主应变值进行两两回归分析，发现回归系数大多位于0.7～0.8或1.2～1.4之间，大多数接近于1，基线差大约为10-8的数倍。但最大主应变的方向角曲线图很不相似，数值相差也很大，可以认为主方向是不可靠的，最大与最小主应变值则是较可靠的。

5 认识与讨论

要采用加衬模型，首先要给出固体潮不变量理论值，本文在采用摩尔圆定理计算出了不变量理论值并由最小二乘原理确定了加衬模型待定常数的基础上，给出了玉树地震所引起的漳州台深井应变不变量同震响应图像，并对同震响应及其可靠性进行了分析。发现最大与最小主应变曲线大体相似，数值也很接近，但最大主应变的方向角曲线图很不相似，数值相差也很大，可以认为主方向是不可靠的，最大与最小主应变值则是较可靠的。出现这种情况可能归于观测误差与加衬模型的不完善性，虽然该模型比摩尔圆定理更符合实际情况，但可能造成上下端面4条测线的影响是各不相同的，同时应变测量误差也会给主方向计算造成极大的误差，且不是一个常数，甚至非常之大。

我国在上世纪90年代末到本世纪初，陆续研制了4分量应变观测系统[11～13]。在相近时期，日本深井观测点的深度达1 000～3 800 m，深井应变观测精度达10-8，而在美国则推出了板块边界观测 （PBO）计划[14、15]，该计划使用了大量的深井钻孔应变仪。对这些深井应变观测资料，据笔者所知，国外学者大多是把注意力放在应变观测资料对Parkfield地震的最大同震信号的分析上[16]，或放在利用Parkfield近场的高精度连续钻孔应变与孔隙压力观测资料对断层成核、断层响应与地震预报的问题上[17]。蒋靖祥等[18、19]，对新疆分量式钻孔应变阶变异常的统计特征及其在地震预测中的应用进行了大量的研究，使用RZB-1型压容式钻孔应变仪的长期观测资料和震例分析，结果表明钻孔应变仪在一定程度上记录大震和台站附近的中、强地震的异常变化。如：乌什钻孔应变仪在多次伽师地震前都记录到张、压性交替变化。2003年2月4日伽师6.8级地震前，出现趋势、短期、临震异常。本文采用应变不变量分析了玉树地震所引起的同震响应，但未提取到震前信息，有待积累更长的观测时间后再进行更深一步的分析研究。

钻孔应变分量变观测系统除了可计算出应变不变量外，同时提供了如此之多可对观测资料进行校核的方法，这是地面观测系统无法比拟的，这既是钻孔应变观测系统的长处，同时也给观测本身提出了更高的质量要求，因为有多余观测，就必须在逻辑上满足许多约束条件。由于有观测误差，我们不能奢求多余观测给出的多组不变量结果都完全一样，而只能要求他们之间的差异局限在一个允许的误差范围内就行，可以说是对钻孔应变观测系统内在质量的一项苛刻要求，仅利用相对应变量照样可以对震前、震时、震后应变状态进行分析。笔者认为，本文是对漳州台深井多分量观测资料是否可用来提取应变不变量所做的一次尝试，从中也反映了目前多分量应变观测中存在的一些问题，虽然深井观测可以大大降低外界干扰，但又产生了一些新的问题，诸如应变测量的定向精度，测量单元的标定，探头与岩壁的耦合以及模型等问题，都是值得认真考虑的，以上一些看法，可供仪器研制与使用者参考。

[1] （美）地壳运动测量专业委员会，构造形变测量 [M]，刘序俨，邹其嘉译，陈鑫连校，北京：地震出版社，1984.

[2] 刘序俨，季颖锋，黄声明，等.地形变应变张量矩阵的不变量分析 [J].大地测量与地球动力学，2011， 31 （4）： 66-70.

[3]潘立宙.与地应力测量有关的几个公式的推导和讨论[A].中国地质科学院，地质力学研究所，国家地震局地震地质大 队.地应力测量的原理和应用[C].北京：地质出版社，1981.

[4]骆鸣津，池顺良，马鸿钧.地层绝对应力测量与钻孔应变测量 [M].北京：地震出版社，2008.

[5]Jaeger JC，ERasticity.Fracture and Flow[M].London Methuen&co.1964.

[6]米恩斯.丁中一译，王仁校.应力与应变 [M].北京：科学出版社 ，1982.

[7]P.梅尔基奥尔.杠品仁，吴庆鹏，陈益惠，等译.行星地球的固体潮[M].北京：科学出版社，1984.

[8]北京大学地球物理系，武汉测绘学院大地测量系，中国科学技术大学地球物理教研室编.重力与固体潮[M].北京：地震出版社，1982.

[9]方俊.固体潮[M].北京：科学出版社，1982.

[10]W.杨，R.布迪纳斯.岳珠峰，高竹山，王峰会译.应力应变公式手册[M].北京：科学出版社，2005.

[11]池顺良.深井宽频钻孔应变地震仪与高频地震学-地震预测观测技术的发展方向，实现地震预报的希望[J].地球物理学进展，2007，22（4）：1164-1170.

[12]欧阳祖熙.RZB-1型电容器式地应变计 [A].第一届全国地应力会议论文选编 [C].北京：地震出版社，1977.

[13] 欧阳祖熙.地壳形变深井综合观测技术的新进展[J].国际地震动态，2009， （11）:1-13.

[14]李海亮，李宏.钻孔应变观测现状与展望[J].地质学报，2010，84（6）：895-900.

[15]陈群策，李宏，廖椿庭，等.地应力测量与监测技术实验研究[J].地球学报，2011，32（supp.I）：113-124.

[16]J.Longbein， J.R.Marray， H.A.Snyder.Coseismic and initial postseismic deformation from the 2004 Parkfield,Californica,Earthquake,observed by global positioning system,electronic distance meter，greepmeter,and bore hole strainmeters[J].Bulletin of the Seismological Society of America， 2006， 96(4b)： 5304-5320

[17]J.S.Johston， R.D.Borcherdt， A.T.Linde， etal.Continuous borehole strain and pore pressure in the near field of the 28 September 2004 M6.0 Parkfield， California， Earthquake， implications of nucleation， fault response,earthquake prediction and tremor[J].Bulletion of the seismological Society of America， 2006， 96(4B)：556-572.

[18]蒋靖祥，光华，王在华，等.新疆分量式钻孔应变阶变异常的统计特征及其在地震预测中的应用研究[J].岩石力学与工程学报，2004，23：4072-4078．

[19] 蒋靖祥，赖爱京.2003年12月1日中哈交界MS6.1级地震钻孔应变异常分析 [J].四川地震，2005， 114(1)： 24-28.

Analysis on Deep-well Strain Invariant's Co-seismic Response to Yushu Earthquake Stimulation at Zhangzhou Seismic Station

LiIU Xuyan1， FANG Hongfang2， ZHANG Kai3，ZHANG Yining3，MA Chaowen2，HUANG Shengming1， JI Yingfeng1

(1.Earthquake Administration of Fujian Province,Fuzhou 350003， China； 2.Earthquake Administration of Zhangzhou City， Zhangzhou 363000， China； 3.Zhangzhou Seismic Sstation of Earthquake Administration of Fujian Province，Zhangzhou,363000,China)

Based on how to solve for the constants determined in a borehole strain observations lining model,this paper uses one of the invariant′s properties to derive five sets of the invariant′s expression using observation data combination.The property is that the invariant has nothing to do with the coordination system selection.Then the paper applies it to the case in Yushu earthquake by selecting better observed datum at the upper and lower surfaces of a detector in an 8-component strain observing system.Then,this paper calculates the maximum and the minimum principal strains,principal strain directions,the maximum shear strain,the maximum strain on a surface and volumetric strain at minute interval.Briefly,this paper selects the datum at 2#3#4#from the upper surfaces of a detector within time window 7:30～9:30 when earthquake occurred resulting in co-seismic response for strain invariants,surface-strain and bulk-strain.Finally,this paper furthermore discusses properties of these co-seismic responses,as well as reliability of the result of maximum,minimum strain,and principal directions of maximum strain given in this paper.

Multi-component strain observation;Lining model;Strain invariant;Co-seismic response;Reliability analysis

P315.727

A

1001-8662（2011）04-0051-11

2011-09-15

2010年度福建省地震局基金课题资助

刘序俨，男，1939年生，研究员.主要从事从事固体潮与地形变研究.E-mail:xuyanliu@126.com.