解滔, 薛艳, 卢军
中国地震台网中心, 北京 100045
地震是构造应力持续积累,并最终超过断层强度导致破裂失稳的结果,对这一过程中断层的力学状态以及可能伴随的地球物理和地球化学现象的研究,有助于认识地震晚期孕育过程这一复杂的科学难题和开展地震中短期预测.对断层面力学状态演化的研究,旨在揭示地震成核过程与滑动速率-状态的依从关系(马瑾等,2012;马瑾和郭彦双,2014;Scholz,1998;Lay and Kanamori,1981;McCaffrey et al.,2008;唐荣江和朱守彪,2020).当前地震晚期孕育过程的研究不再局限于发震断层面本身,已经逐步扩展至研究整个孕震区与高应力-应变水平积累有关的介质变化(Scholz et al.,1973;Mjachkin et al.,1975;Roeloffs,1988;Hartmann and Levy,2005).许多文献和研究报到了与地震有关的地球物理和地球化学观测数据变化,但存在较大的争议.争议的核心问题之一在于,这些数据的变化,与地震前震源区及附近区域的应力积累或介质变形之间,并未建立明晰的联系.如果地震前观测数据变化的原因是应力-应变的变化,那么在分析这些数据变化与地震孕育过程之间的关系时,需要有应力-应变变化的背景作为参照(吴忠良等,2009).
电阻率是地下岩土介质最基本的物理属性之一,应力作用下介质内裂隙率和微裂隙结构会发生改变,并引起电阻率发生变化.20世纪50—80年代,日本、苏联、美国先后开展了用于地震预测研究的视电阻率实验观测,并报道了多次中等地震前突出的视电阻率下降变化(Barsukov and Sorokin,1973;Nersesov et al.,1979;Mazzella and Morrison,1974;Morrison et al.,1977,1979;Park,1991;Madden et al.,1993),以及地震发生前后观测值的“准同震”阶跃变化,位于海边山洞内的观测站记录到了与潮汐同步的变化(Rikitake and Yamazaki,1967,1969,1976;Yamazaki,1974,1975).全球唯有中国在全国主要地震活动区域内,持续开展了规模化和规范化的连续观测,在百余次5~8级地震前,震中附近观测站记录到了持续数月至两年左右的异常变化及其恢复过程(钱复业等,1982;国家地震局预测预防司,1998;钱家栋等,1985;赵玉林等,2001;汪志亮等,2002;杜学彬,2010;Lu et al.,2016).
本文将利用我国自1967年以来的视电阻率连续观测数据,结合发生在观测站网内及附近的16次/组MS≥7.0地震,采用断层虚位错模式(赵玉林等,1996;解滔等,2020),将地震同震位移按大小相等但方向相反的方式进行加载,获取地震前产生这部分同震位移所需的应力-应变积累及其空间分布范围和特征,分析视电阻率下降/上升变化与挤压/膨胀变形区域之间的关系;并以此为纽带,从“介质变形-电阻率变化”的角度,尝试将应力作用下岩土介质电阻率变化的细观机制,与地震前视电阻率变化的宏观现象相联系.
1966年3月河北邢台7.2级地震之后,原中国科学院兰州地球物理研究所在邢台震区架设了第一个地电观测站,并于1967年4月投入观测,正式开始了视电阻率方法在地震监测预报中的应用研究.观测采用对称四级装置,每个观测站布设两个或三个方向的观测(图1),多数观测站供电极距为500~3000 m不等.截至目前,全国共有89个视电阻率观测站,主要分布在南北地震带、华北、东北和新疆等地震重点监视地区,其中23个观测站具有井下观测装置.
图1 我国视电阻率观测采用的对称四级装置及布极方式(a) 两个方向的观测装置; (b) 三个方向的观测装置.Fig.1 The Schlumberger resistivity arrays and its arrangements used in apparent resistivity monitoring stations in China(a) Arrays along two directions; (b) Arrays along three directions.
观测主要采用人工直流供电的方式,多数观测站供电电流1~3 A,每小时观测一次.目前采用的地电仪测量人工电位差分辨力为0.01 mV,在不考虑背景环境噪声的情况下,可分辨约万分之几的视电阻率变化.定期采用标准电阻或标准电源对仪器测量系统进行检定,确保观测数据的真实性;观测系统具有长期稳定性,在背景噪声较低的观测站,观测精度优于1‰.
杜学彬(2010)对中国192次MS≥4.0地震前674次视电阻率变化进行了分析,认为地震前与地震孕育过程有关的视电阻率变化的最大空间范围约400 km左右,其中300 km范围内居多.自1967年开始观测以来,中国共计16次/组MS≥7.0地震发生在观测站约400 km范围内(图2).通过文献调研和历史数据分析,这些地震前的视电阻率变化特征示于表1.
表1 观测站网内16次/组MS≥7.0地震前的视电阻率变化Table 1 Theapparent resistivity changes before the 16 times/groups earthquakes of MS≥7.0
续表1
图2 发生在视电阻率观测站网内16次/组MS≥7.0地震(图中含1990年景泰MS6.2地震)Fig.2 The 16 times/groups earthquakes of MS≥7.0 within the monitoring networks (The Jingtai MS6.2 earthquake in 1990 is included)
环渤海地区共发生3次/组MS≥7.0地震,分别为1969年渤海MS7.4地震,1975年辽宁海城MS7.3地震,1976年河北唐山—滦县MS7.8和MS7.1地震.1970年之前,视电阻率实验观测限于邢台—河间地震区,距离渤海地震约400 km,地震前7—10个月开始四个观测站出现下降变化.自1970年开始,华北地区开始建设规模化的观测站网,北京附近地区的站网密度相对较大.从1974年开始,华北东部地区较大范围内9个观测站先后出现下降变化,1975年海城地震发生之后下降变化仍然持续,直至1976年唐山地震之后才逐渐恢复.而辽宁地区的台安观测站则在海城地震前出现上升变化.唐山地震前约半个月,唐山、马家沟和昌黎站出现加速下降形态的临震变化,直至地震发生.
1970年云南通海MS7.8地震之后开始在南北地震带进行视电阻率观测,之后云南地区共发生5次/组MS≥7.0地震,分别为1974年大关MS7.1地震,1976年龙陵MS7.3和MS7.4地震,1988年澜沧—耿马MS7.6和MS7.2地震,1995年孟连MS7.3地震、1996年丽江MS7.0地震.1974年大关地震前15个月开始,会理站开始出现上升变化,4个月后西昌和米易站出现下降变化;雅安站1973年4月至10月存在缺数,恢复观测后数据呈上升变化,但与之前的数据存在阶跃下降,地震前的年尺度变化形态不易判断.1976年龙陵地震前楚雄站出现上升变化,而腾冲站1974年7月至1976年3月受干扰严重,无法分析龙陵地震前的变化情况.1988年云南澜沧—耿马地震前腾冲站出现持续20个月的下降变化;同期通海站也出现下降变化,但后续分析发现主要为测区内城镇化建设的影响(杜家伦,1990);楚雄和元谋站未观测到地震前的视电阻率异常变化.1995年孟连地震和1996年丽江地震前未出现视电阻率异常变化.
四川地区共发生5次/组MS≥7.0地震,分别为1973年炉霍MS7.6地震,1976年松潘平武两次MS7.2地震,2008年汶川MS8.0地震,2013年芦山MS7.0地震,2017年九寨沟MS7.0地震.1973年炉霍地震前约半年内,松潘、雅安和甘孜站相继出现下降变化,而康定站则出现上升变化. 1976年松潘—平武地震14个月前,武都和松潘站相继出现下降变化,且震前半月内出现临震加速下降变化(Lu et al.,2016).2008年汶川地震前22个月开始,成都、甘孜和江油站同步出现下降变化,而武都站出现上升变化.2013年芦山地震前成都站出现下降变化,甘孜站未出现异常变化,江油站在汶川地震后停测,于2015年12月恢复观测.2017年九寨沟地震前玛曲站出现下降变化,武都站于2014年改建为井下观测,垂直方向于2016年出现上升变化,天水站则表现为高频扰动变化.
青海地区共发生3次MS≥7.0地震,分别为1990年共和MS7.0地震、2010年玉树MS7.1地震和2021年玛多MS7.4地震.1990年共和地震前,武威站NS方向观测在1988年年底因不明原因导致年变形态消失,并持续多年;EW方向观测未受影响,自1989年下半年开始出现上升变化,而1990年10月20日距武威站约130 km发生景泰MS6.2地震.甘孜站距离2010年玉树地震约365 km,2009年4月—9月期间出现仪器故障,恢复观测后存在阶跃型变化,难以确认玉树地震前的真实变化形态.2021年玛多地震前,400 km范围内有甘孜、玛曲、白水河、金银滩和拦隆口5个观测站,其中玛曲站距离最近,约356 km,地震前均未出现明显的异常变化.
地震是构造应力在断层闭锁段长期作用,最终导致断层失稳错动的结果.地震前发震断层及附近区域以介质变形的形式积累应变能,部分应变能以断层错动的方式予以释放,并产生同震滑动.将地震同震位错按大小相等但方向相反的方式进行加载(图3),获取地震前能够产生这部分同震滑动所需的应力-应变积累的空间分布特征,这构成断层虚位错模式的基本思想(赵玉林等,1996;解滔等2020).
图3 断层同震位错与虚位错模式(a)逆冲断层;(b)正断层;(c)走滑断层.Fig.3 Coseismic dislocation of fault and virtual dislocation model(a) Thrust fault; (b) Normal fault; (c) Strike-slip fault.
虚位错模式的计算结果,仅表示这部分虚位移所能引起的应力-应变的变化量.构造区域内的真实应力-应变水平,应为地震发生之后的绝对应力-应变水平加上这部分变化量.但是,在地震发生之前及之后,构造区域内的绝对应力-应变水平是难以获取的.在整体为挤压环境的构造区域,逆冲、走滑、或逆冲兼走滑型地震计算结果中的挤压区域,是地震前挤压增强的区域;而计算结果中的拉张区域,并不能区分其是绝对的拉张或挤压区域,但可认为是相对膨胀的区域.由于区域内绝对的应力-应变水平难以获取,且应力水平与岩土介质微裂隙活动之间的定量关系也未明晰,目前仅能根据体应变变化量的分布,获取地震前挤压增强区和相对膨胀区域,并在此基础上开展分析工作.本文采用断层滑动位错模型分析软件Coulomb3.3(Lin and Stein,2004;Toda et al.,2005)进行断层虚位错模式的计算.
16次/组MS≥7.0地震的震源机制解和同震滑动模型示于表2.1995年孟连地震、1996年丽江地震、2010年玉树地震和2021年玛多地震前无视电阻率异常变化,本文未对其进行断层虚位错模式计算;2008年汶川地震采用Shen(2009)给出的精细同震滑动模型.武威站1989年至1990年的数据变化可能同时受1990年共和地震和景泰地震影响,计算时同时考虑了这两次地震.表2中除部分地震外,其余采用Chen等(2017)基于发生在中国的地震统计给出的面波震级和矩震级之间的经验关系进行估计,断层破裂参数通过以下经验关系计算(Wells et al.,1994).
表2 观测站网内16次/组MS≥7.0地震的震源机制解和同震滑动模型Table 2 Focal mechanism solutions and coseismic slip models of the 16 times/groups earthquakes of MS≥7.0 within apparent resistivity monitoring networks
走滑型地震:
(1)
式中LR为断层破裂长度(单位:km);WR为断层破裂宽度(单位:km);DR为断层平均滑动位移(单位:m);MW为矩震级.
逆冲型地震:
(2)
正断型地震:
(3)
同时含有逆冲和走滑、或正断和走滑分量的混合型地震:
(4)
环渤海地区3次/组地震断层虚位错模式体应变变化的计算结果示于图4.1969年渤海地震前,大柏舍、柳竹、大曹庄和牛家桥4个观测站呈下降变化,均位于震前挤压增强区域(图4a).1975年海城地震前,台安站呈上升变化,位于震前相对膨胀区域;昌黎站呈下降变化,位于挤压增强区域(图4b).1976年唐山地震前,昌黎、马家沟、唐山、青光、宝坻、西集、八里桥、小汤山和忠兴庒9个观测站为下降变化,位于地震前挤压增强区域;徐庄子站为上升变化,位于相对膨胀区域(图4c),与赵玉林等(1996)给出的结果是一致的.
图4 环渤海地区3次/组MS≥7.0地震断层虚位错模式计算结果(挤压为负) (a) 1969年渤海MS7.4地震; (b) 1975年海城MS7.3地震; (c) 1976年唐山MS7.8和滦县MS7.1地震.Fig.4 Results from fault virtual dislocation model for 3 times/groups earthquakes of MS≥7.0 occurred in Bohai Rim area (Compression is negative)(a) The Bohai MS7.4 earthquake in 1969; (b) The Haicheng MS7.3 earthquake in 1975; (c) The Tangshan MS7.8 and Luanxian MS7.1 earthquakes in 1976.
云南地区3次/组地震断层虚位错模式体应变变化的计算结果示于图5.1974年大关地震前西昌和米易观测站为下降变化,位于地震前挤压增强区域;会理站为上升变化,位于相对膨胀区域(图5a).1976年龙陵地震前,楚雄站位于相对膨胀区域(图5b),呈现上升变化. 1988年澜沧—耿马地震前,腾冲站位于挤压增强区域(图5c),表现为下降变化.
图5 云南地区3次/组MS≥7.0地震断层虚位错模式计算结果(挤压为负)(a) 1974年大关MS7.1地震; (b) 1976年龙陵MS7.4、MS7.3地震;(c) 1988年澜沧-耿马MS7.6、MS7.2地震.Fig.5 Results from fault virtual dislocation model for 3 times/groups earthquakes of MS≥7.0 occurred in Yunnan province (Compression is negative)(a) The Daguan MS7.1 earthquake in 1974; (b) The Longling MS7.4 and MS7.3 earthquakes in 1976; (c) The Lancang MS7.6 and Gengma MS7.2 earthquakes in 1988.
四川地区5次/组地震断层虚位错模式体应变变化的计算结果示于图6.1973年炉霍地震前,甘孜、松潘和雅安站呈现下降变化,位于挤压增强区域;康定站则出现上升变化,位于相对膨胀区域(图6a).1976年松潘—平武地震前,松潘和武都站位于挤压增强区域,呈现下降变化(图6b).2008年汶川地震前,成都、江油和甘孜站位于挤压增强区域,呈现下降变化;武都站位于相对膨胀区域,呈上升变化(图6c).2013年芦山地震前,成都站为下降变化,位于挤压增强区域(图6d).2017年九寨沟地震前,玛曲站位于挤压增强区域,呈下降变化,武都站位于相对膨胀区域,呈上升变化(图6e);天水站在九寨沟地震前的高频扰动变化,可能反映中短期阶段观测区电场信号的不稳定性(杜学彬等,2017).
图6 四川地区5次/组MS≥7.0地震断层虚位错模式计算结果(挤压为负)(a) 1973年炉霍MS7.6地震; (b) 1976年松潘—平武两次MS7.2地震; (c) 2008年汶川MS8.0地震; (d) 2013年芦山MS7.0地震; (e) 2017年九寨沟MS7.0地震.Fig.6 Results from fault virtual dislocation model for 5 times/groups earthquakes of MS≥7.0 occurred in Sichuan province (Compression is negative)(a) The Luhuo MS7.6 earthquake in 1973; (b) The Songpan-Pingwu double MS7.2 earthquakes in 1976; (c) The Wenchuan MS8.0 in 2008; (d) The Lushan MS7.0 earthquakes in 2013; (e) The Jiuzhaigou MS7.0 earthquakes in 2017.
图7为1990年共和地震和景泰地震断层虚位错模式体应变变化量的叠加结果.如果对两次时间相隔半年的地震单独进行分析,武威站位于共和地震前的挤压增强区域,但位于景泰地震前的相对膨胀区域.两次地震的叠加结果显示,武威站位于地震前的相对膨胀区域,且观测数据呈上升变化.由此可见,因武威站距离景泰地震更近,其观测数据变化可能和景泰地震的晚期孕育过程之间关系更为紧密.
图7 1990年共和MS7.0和景泰MS6.2地震断层虚位错模式计算结果(挤压为负)Fig.7 Results from fault virtual dislocation model for the Gonghe MS7.0 and Jingtai MS6.2 earthquakes in 1990 (Compression is negative)
以上分析结果显示,12次/组MS≥7.0地震前,视电阻率出现下降变化的观测站均位于地震前的挤压增强区域,而出现上升变化的观测站则位于相对膨胀区域.武都站在三次/组地震前分别出现了下降和上升变化,且变化形态与观测站所处区域在地震前的挤压增强和相对膨胀的变形特征相对应.
为了分析地震前视电阻率变化的原因,国内外学者从实验和理论模型方面开展了大量的研究,试图从“介质变形-电阻率变化”的角度,建立视电阻率变化与地震晚期孕育过程之间力学机制上的联系.实验室内应力加卸载实验结果显示,含水岩石在压应力加载过程中电阻率呈现下降变化,卸载过程中呈现上升恢复变化,多数岩石在临近破裂时加速下降,岩石破裂后出现快速回返(图8a),且横向和纵向两个方向测量时电阻率变化幅度存在差异(图8b);无水岩石在压应力作用下则呈现小幅度上升变化(Brace et al.,1965;Yamazaki,1966;Morrow and Brace,1981;张金铸和陆阳泉,1983;Jouniaux et al.,2006).野外原地实验结果显示,压应力加载时视电阻率出现下降变化,应力卸载过程中呈恢复变化;在地表不同方向观测时视电阻率呈现出各向异性变化,即垂直于应力加载方向观测的变化幅度最大,平行方向观测时变化幅度最小,斜交方向介于二者之间;变化幅度随应力加载强度的增加而增加,随力源距离的增加而减小(赵玉林等,1983;国家地震局预测预防司,1998).
图8 实验室内岩石样本应力加载时电阻率变化(a)轴向应力加载至岩石破裂时应变-电阻率变化(张金铸等,1983);(b)轴向应力加载时横向和纵向电阻率变化(修改自Brace el al.,1965).Fig.8 Resistivity changes of rock samples under stress loading in laboratory(a) Strain-resistivity changes from axial stress loading to rock fracture (Zhang et al., 1983); (b) Changes of transverse and longitudinal resistivity under axial stress loading (modified from Brace el al., 1965).
在实验中同步测量岩石样品的体积变化,发现在应力加载超过破裂应力强度的50%之后,岩石的体积开始增大而出现扩容,并伴随有声发射现象,反映出新裂隙的不断产生(Brace et al.,1966; Brace and Orange,1968;Brace,1975;Scholz,1967).基于该实验现象,学者提出了膨胀扩散模式(DD模式)和裂隙雪崩扩展模式(IPE模式)用以解释地震前的地球物理和地球化学异常变化(Scholz et al.,1973;Mjachkin et al.,1975).DD模式和IPE模式均强调在高应力承载条件下,介质内部存在微裂隙活动. 在这两种模式中,地震前介质电阻率均呈现持续性下降变化(图9).这些实验研究基本上明确了在实验样品尺度下,岩土介质在应力作用下诱发的微裂隙活动是电阻率变化的重要原因.
图9 地震异常变化的两种物理模式(a) 膨胀扩散模式(修改自Scholz et al.,1973); (b) 裂隙雪崩扩展模式(修改自Mjachkin et al.,1975).Fig.9 Two physical models for earthquake precusors(a) The dilatancy and diffusion model (modified from Scholz et al., 1973); (b) The fracture avalanche expansion model (modified from Mjachkin et al., 1975).
岩土力学实验和理论研究表明,对于低围压环境下初始含裂隙的介质,在压应力持续加载到一定程度之后,新生裂隙不断产生;无论初始裂隙如何分布,最终形成的新裂隙系统的优势展布方向将大致沿最大主压应力方向(李新平等,2002;张恒等,2015;张志强等,2020).在高围压条件下,引起完整岩石产生微破裂和扩容所需的应力水平将大于引起断层滑动所需的力. 但是,处于低围压水平的浅部松散介质,在相对较低的应力作用下,介质微裂隙的产生和扩容可通过松散物质颗粒的再移动实现(马瑾,1982).解滔和卢军(2020a)推导出了含裂隙介质的等效电阻率公式,以及电阻率变化与裂隙变化之间的微分关系.在微裂隙系统沿最大主压应力方向优势排列和扩展的假设下,裂隙的微小变化将引起较大幅度的电阻率变化,且沿裂隙优势展布方向的最小电性主轴变化幅度最大(Xie et al.,2020).
我国定点观测站的视电阻率观测,其深度探测范围与供电极距相当(赵和云和钱家栋,1982;杜学彬等,2008),主体探测范围位于低围压的地壳浅表地层内,且位于潜水位之下,观测值是千米尺度范围介质电阻率的综合反映.一般而言,在数月至两年左右的时间段内,非地热地区数十米至千米深度范围内,地下介质温度基本保持不变;水溶液有足够的时间进出裂隙,含水率也基本保持不变.因此,分析认为介质裂隙结构的变化,是地震前视电阻率变化的主要原因(钱家栋等,1985;国家地震局预测预防司,1998;杜学彬等,2007).大地震前近震中区域的视电阻率观测呈现出与主压应力方位有关的各向异性变化,表现为垂直于最大主压应力方位的测道变化幅度最大,平行方向最小,斜交方向介于二者之间(杜学彬,2010),比如汶川地震前成都站和江油站不同方向的视电阻率变化,与主震分段震源机制解给出的主压应力方位之间存在较好的对应关系(解滔和卢军,2020a).据不同学者的研究结果,1976年唐山MS7.8地震的主压应力方位介于71°~85°之间(中国震例,1976—1980),滦县MS7.1地震主压应力方位为75°(USGS).昌黎站位于滦县地震以东,距离约35 km,取75°作为昌黎站附近的主压应力方位;其余观测站位于唐山主震发震断裂以西,这里取平均值78°作为主压应力方位.各观测站测道与P轴之间的夹角示于表3,结合表1中的变化幅度可见,对于震中距140 km范围内的9个观测站,与P轴夹角越大的测道变化幅度越大.采用对称四级装置在地表进行观测时,视电阻率各向异性变化与水平电性主轴真电阻率变化之间存在π/2的角度差.因此,唐山地震的近震中区域范围内,各观测站视电阻率的各向异性变化特征与实验室内实验、野外原地实验、以及含裂隙介质电阻率模型给出的结果一致.此外,在构造应力作用下地层电阻率发生变化的上界面深度对视电阻率各向异性变化也产生影响,上界面距离地表越近,各向异性变化特征越显著(解滔和卢军,2020b).
表3 唐山地震周围观测站各测道与主压应力方位的夹角Table 3 Angles between the measuring channels of the stations around the Tangshan earthquake and the principal compressive stress
实验结果和含裂隙介质电阻率公式,从细观尺度呈现出电阻率变化与裂隙变化之间的关系,而震例分析则从宏观尺度呈现出地震前视电阻率变化的统计特征.二者之间需要有“介质变形-电阻率变化”这一中间过程的联系,即视电阻率出现变化的观测站,其所在的区域在地震前受到震源区高应力-应变水平积累的影响特征.以上对12次/组MS≥7.0地震的虚位错模式分析结果显示,视电阻率的下降/上升变化,与地震前的挤压增强区/相对膨胀区之间具有一致的对应性,这揭示出视电阻率变化和地震之间可能存在的联系过程:在地震晚期孕育阶段,震源区处于较高的应力-应变水平,并引起附近区域介质也产生不同程度和形态的变形,诱发地下浅层范围内介质的微裂隙活动,进而引起介质电阻率变化;这一过程中,震中周围视电阻率的下降/上升变化与观测站所在区域的介质变形特征有关.此外,对1976年唐山地震前震中周围视电阻率观测数据的分析发现,由震中向外围方向,变化幅度整体呈现减弱的空间分布,这也符合震源区应力-应变积累程度较高、向外围方向逐渐衰减的分布特征(赵玉林等,1996).
通过文献调研和历史数据分析,自观测以来我国共计16次/组MS≥7.0地震发生在视电阻率观测站网约400km范围,其中12次/组地震前出现视电阻率变化.采用断层虚位错模式计算了这12次/组地震前应变变化量的空间分布,视电阻率出现下降变化的观测站位于地震前的挤压增强区域,而出现上升变化的观测站则位于相对膨胀区域.已有实验结果和理论分析表明,岩土介质在应力作用下诱发的微裂隙活动会引起电阻率发生变化.应力作用下介质电阻率变化的微观机制,可通过地震前孕震区及附近区域的相对变形特征,与震中周围的视电阻率变化联系起来.因此,地震前的视电阻率变化,可能与地震晚期孕育过程之间存在“介质变形-电阻率变化”机制上的联系.
致谢四川省地震局和云南省地震局提供了部分观测站视电阻率的历史观测数据,两位审稿老师提出了中肯的修改建议,对文章的完善有很大的帮助,在此表示衷心的谢意.