流体是地球形成与演化中的主要因素,是物质运移和能量交换的重要途经(车用太等,2006)。研究表明,地震前后区域构造应力场的变化会改变地下流体的循环体系,并造成流体间的混合或流体和岩石间的相互作用,导致地下流体化学组分的异常变化(高小其等,2002)。前人针对异常变化的提取进行了探索,并提出了一些分析方法。例如,和宏伟等(1999)以云南地区5级以上地震为研究对象,采用一阶差分法对1987~1997年云南地区23个水氡观测点的水氡观测数据进行了处理,对应率达到70%左右;范雪芳等(2002)总结了地下流体中期和中短期前兆异常的4种判定方法;叶秀薇(2004)利用从属函数对粤闽地区地下流体与地震之间的关系进行了初步研究;杨兴悦等(2011)利用从属函数、变差率以及趋势速率方法提取甘肃东南部水氡资料的地震前兆异常,获得的异常信号很好地对应了600km以内5.0≤M≤8.0地震。此外,研究中常用到的方法还有高通滤波法、多项式拟合法、经验模态分解法、概率密度分布法等(邱泽华等,2010;钟伟等,2006;Huang et al,1998;Manshour et al,2009、2010)。
长期以来,国内外学者围绕地震孕育、发展、发生与地下流体的关系进行了探索,提出了不同理论、模式来进行解释。如,Scholz等(1973)提出了解释地震成因的膨胀扩容模式,米雅奇金(1983)提出了岩石破坏的雪崩不稳定裂隙形成模式,郭增建等(1973)提出了能量积累单元与调整单元构成的组合模式。上述模式均在一定程度上解释了某些大地震前出现的部分异常现象。
本文采用从属函数法和自适应阈值法,对江苏省盐城市纺织厂井(以下简称盐纺井)的Ca2+、Cl-浓度观测资料进行分析,尝试提取Ca2+和Cl-浓度在地震前的异常变化。
1 盐纺井简介 1.1 水文地质概况盐纺井位于江苏省盐城市大庆中路,测项为Ca2+、Cl-浓度。观测站周围为第四系覆盖层,覆盖层厚度在千米以上,无基岩出露。测点地处盐城断裂西端,向北毗邻NE向展布的洪泽-沟墩断裂,以西靠近泾口-沙沟断裂,向东为南黄海区内NW向沿海岸线的滨海大断裂。因此,台站受4条断裂形成的“口”字型构造复合控制(图 1)。地下水类型为自流抽水井,采样深度为600m以下,属热水井,水温最高可达48℃;水质类型为HCO3-Na型(表 1)。台站自1990年开始观测,原始值为日值。
根据《地震水文地球化学观测技术规范》(中国地震局,2014),Ca2+浓度采用EDTA容量法测定,Cl-浓度采用硝酸银容量法测定。
2 观测资料选取和目标地震的筛选 2.1 观测资料的选取本文主要分析该测点2007~2016年的月均值数据,Ca2+和Cl-浓度月均值原始曲线如图 2所示,两者含量具有一定的准同步性。
目标地震选取震中距200km范围内ML≥4.0的地震,若一年内在同一区域多次发生中强地震,则以震级较大的地震为主,最终选取的地震见表 2。
选取盐纺Ca2+和Cl-离子2007~2016年的月均值为数据序列,依次用从属函数和自适应阈值等方法进行异常的提取与分析。
3.1 利用从属函数法进行资料分析(1) 计算:从属函数法公式为
$ {\mu _i} = {\left({1 + \frac{\alpha }{{\left| {{k_i}} \right| + \left| {{\gamma _i}} \right|}}} \right)^{ - 1}} $ | (1) |
$ {V_i} = \left\{ {\begin{array}{*{20}{c}} {{\mu _i}\left({{k_i} > 0} \right)}\\ {0\left({{k_i} \le 0} \right)} \end{array}} \right. $ | (2) |
其中,μi为从属函数值;ki为观测时间序列M(t)的斜率;γi为M(t)与t的相关系数;α为经验函数,由小到大进行选取,直至异常表现突出;拟合窗长为n=5;由于盐纺井离子浓度上升为异常,因此只取ki>0的μi来描述,记作Vi。
(2) 异常指标:当盐纺离子浓度的从属函数值V≥0.5时,视为异常。
(3) 异常特征:在统计时间段内,盐纺井Ca2+离子浓度出现了6次异常,Cl-离子浓度出现了4次异常,异常特征及其对应地震见图 3和表 3。
盐纺井Ca2+离子浓度从属函数异常持续时间最长为11个月,最短为3个月;盐纺井Cl-离子浓度从属函数异常持续时间最长为13个月,最短为3个月。异常的持续时间与震级之间没有明显的数量关系。
(4) 对应地震情况:在统计时间段内,盐纺井Ca2+离子的6次异常有4次对应地震,对应率为60%,漏报率为20%,虚报率为40%。盐纺Cl-离子的4次异常均对应地震,对应率为100%,漏报率为20%,虚报率为0。
3.2 利用自适应阈值法进行异常提取(1) 计算:自适应阈值法公式为
$ {F_i} = \left| {{y_i} - \frac{{\sum\limits_{j = i - L}^{i - 1} {{y_i}} }}{L}} \right| - K{S_i}\;\;\;\left({i = L + 1, L + 2, \cdots, n} \right) $ | (3) |
其中,Fi为自适应阈值;yi为观测值;
(2) 异常指标:当盐纺离子浓度的自适应阈值超过其2倍均方差时,视为异常。
(3) 异常特征:在统计时间段内,盐纺井Ca2+离子、Cl-离子浓度均出现了5次异常,异常特征及其对应地震见图 4和表 4。
盐纺井Ca2+离子浓度自适应阈值异常持续时间最长为3个月,最短为2个月;盐纺井Cl-离子浓度从属函数异常持续时间最长为6个月,最短为1个月。异常的持续时间与震级之间没有明显的数量关系。
(4) 地震对应情况:在统计时间段内,盐纺井Ca2+、Cl-离子浓度均出现5次异常,有4次异常发生后对应地震,有1次异常发生在震后,对应率为80%,漏报率为20%,虚报率20%。
4 盐纺井Ca2+、Cl-映震机理讨论由于盐纺井仅仅开展了Ca2+和Cl-离子观测,没有其它测项可以协助开展机理分析,相邻地区也无其它流体观测手段,因此,本文以2012年高邮ML5.3地震为例,以盐城台(距盐纺井约10km)的波速比和地磁谐波振幅比为主要依据,使用膨胀扩容(DD)模式讨论盐纺井Ca2+和Cl-离子的映震机理。
4.1 盐城台波速比利用2009年1月~2013年12月盐城台120km范围内的地震目录资料,计算单台波速比,计算方法为
$ \frac{{{v_{\rm{p}}}}}{{{v_{\rm{s}}}}} = 1 + \frac{{n\sum\limits_{i = 1}^n {\Delta t_i^2 - {{\left({\sum\limits_{i = 1}^n {\Delta {t_i}} } \right)}^2}} }}{{n\sum\limits_{i = 1}^n {\Delta {t_i}{t_{{\rm{P}}i}} - \sum\limits_{i = 1}^n {{t_{{\rm{P}}i}}\sum\limits_{i = 1}^n {\Delta {t_i}} } } }} $ | (4) |
式中,tPi为P波走时,Δti为P波和S波的走时差,i=1,2,…,n(n为每次地震到时的数据个数)。
计算波速比的约束条件为:Δt≤15s,相关系数R≥0.95,波速比误差估计σ≤0.05,震级下限ML1.5,选择视窗18个月的窗口,逐月滑动时间窗,计算得到波速比(vP/vS)随时间的变化。
图 5为参与计算的地震分布情况,图 6为波速比计算结果。由图 6可见,在2012年7月20日高邮地震前,盐城台波速比出现了显著的下降变化,并于震后迅速恢复,具体过程可以分为3个阶段。第一阶段(2010年8月~2011年2月),随着应力的积累,新的微裂隙产生或原有裂隙扩张,岩体发生膨胀和扩容,波速比开始下降,岩体的扩容导致了本身的孔隙压降低,孔隙中的流体呈现不饱和状态,深部的卤水开始渗入到裂隙中;第二阶段(2011年3月~2012年6月),由于卤水的渗流速度较慢,因此波速比到达最低点后将在低值波动一段时间;第三阶段(2012年7月~2012年8月),随着深部卤水的渗入,岩体的孔隙压和弹性模量增大,波速比回升,岩体加速破裂进而发震。
地磁谐波振幅比YZHx(NS)和YZHy(EW)的定义为
$ {Y_{ZHx}}\left({{\rm{NS}}} \right) = \left| {\frac{{Z\left(\omega \right)}}{{{H_x}\left(\omega \right)}}} \right| $ | (5) |
$ {Y_{ZHy}}\left({{\rm{EW}}} \right) = \left| {\frac{{Z\left(\omega \right)}}{{{H_y}\left(\omega \right)}}} \right| $ | (6) |
其中,Z(ω)、Hx(ω)、Hy(ω)分别为地磁场的垂直分量、南北和东西向水平分量的频谱,ω为圆频率。
冯志生等(2013)指出,对于随时间周期变化的不均匀场源,在地球介质为均匀各向同性的平面导体的条件下YZHx(NS)和YZHy(EW)与地下介质的电阻率ρ呈正比,即YZHx(NS)和YZHy(EW)的变化可以近似反映地下介质电阻率ρ的变化。
利用谐波振幅比研究盐城台周边地下介质的电阻率变化,资料选取、处理和计算过程为:①选取盐城台地磁三分量F、H、D相对观测资料,并转换成Z、Hx、Hy;②利用富式拟合去噪声,加汉宁窗去高频干扰,窗长为6h;③计算每分钟的三分量富氏谱Z(ω)、Hx(ω)、Hy(ω);④计算10~60min共6个频带的地磁谐波振幅比,并使用年滑动平均法去年变化,计算结果如图 7所示。
盐城台谐波振幅比在多个周期内出现了下降—转折—恢复的变化特征。下降的起始时间约在2010年7月~2011年2月,转折时间均集中在2012年1~3月,恢复所用的时间不完全一致,跨度较短的于2012年10月恢复,跨度较长的在2014年才缓慢回升。
4.3 膨胀扩容(DD)模式根据Scholz等(1973)提出的膨胀扩容模式,孕震过程可以分为弹性应变、膨胀裂隙形成和地下水扩散3个阶段。如图 8所示,2012年高邮ML5.3地震孕震过程与这一模式较为符合,具体的过程如下:
(1) 弹性应变阶段(2010年8月之前),该阶段还未出现前兆现象。
(2) 膨胀裂隙形成阶段(2010年8月~2012年2月),因膨胀而新形成的裂隙中,由于流入了流体,导致震源区的介质处于不饱和的状态,引起盐城台的波速比和谐波振幅比(深部电阻率)出现明显的下降变化。
(3) 地下水扩散阶段(2012年2月~2012年7月),由于地下水的进一步扩散,膨胀岩石的物理性质恢复,使得波速比、电阻率等转折回升;当矿化度较高的深层卤水扩散至浅层,并与浅层地下水混合后,盐纺水化站的Ca2+、Cl-浓度上升;由于流体的渗入,破裂过程得以加速,进而发震。
5 结论与讨论 5.1 地震中短期指示意义通过对盐纺井Ca2+和Cl-离子浓度进行的异常提取,结果显示Ca2+和Cl-离子浓度异常出现时间多在震前半年左右,对台站200km范围内的地震具有较好的中短期指示意义。
5.2 映震灵敏性成因分析中国地震局地下流体学科技术管理组于2015年和2016年对盐纺井进行了水化学成分测试,获得了水质、氢氧同位素和氦同位素结果①。测试结果显示,盐纺井离子特征与该地区地下水特征基本一致,但是总溶解固体(TDS)特征低于区域地下水特征,表明该井属于较封闭的地下水环境;Na-K-Mg三角图显示该井水属于部分平衡水范畴(图 9),表明该井水的水-岩反应相对较为充分,水流系统相对稳定;地热温标显示该井水的热储温度介于71~150℃之间,热储深度大致在2.3~5km;氢氧同位素测值显著偏离全球大气降水线(图 10),其氧同位素富集可能表明该井水为大气降水补给,且水体经过了长距离的运移和蒸发作用,与该井水较低TDS、时间序列无明显年变特征一致,2015年测量的3 He/4 He结果也显示该井水壳源He同位素贡献率高达62%,大气He同位素贡献率为38%。
① 高小其, 2016, 2016年华东流体异常现场核实工作汇报,中国地震局地壳应力研究所内部资料.
地下水地球化学成因数据显示该井地下水补给循环较深,不易受浅层降水及环境干扰。波速比和地磁谐波振幅比的结果表明,盐纺井Ca2+、Cl-离子浓度在震前的快速上升可能与深部流体上涌有关。井水的水-岩反应程度适中,比较适合开展地震地下流体地球化学观测,较易获取地震前兆异常信息,这可能是该井Ca2+和Cl-具有较好映震灵敏性的原因。
然而,本文研究处于初步阶段,仍然存在不足之处,比如每一次地震前是否均会出现类似于2012年高邮ML5.3地震前的3个阶段,还有待进一步深入研究。
致谢: 中国地震局地壳应力研究所高小其教授在盐纺井水文地球化学资料收集过程中给予了支持和帮助,并对论文撰写给予指导,在此表示感谢。
车用太、鱼金子, 2006, 地震地下流体学, 北京: 气象出版社.
|
范雪芳、张淑亮、王吉易, 2002, 地下流体中期和中短期前兆异常的四种判定方法, 地震, 22(4): 136-139. |
冯志生、李鸿宇、张秀霞等, 2013, 地磁谐波振幅比异常与强地震, 华南地震, 33(3): 9-15. DOI:10.3969/j.issn.1001-8662.2013.03.002 |
高小其、迪里夏提、许秋龙等, 2002, 乌鲁木齐9、10号泉F-映震特征的研究, 地震学刊, 22(1): 5-11. DOI:10.3969/j.issn.1672-2132.2002.01.002 |
郭增建、秦保燕、徐文耀等, 1973, 震源孕育模式的初步讨论, 地球物理学报, 16(1): 43-48. |
和宏伟、和国文、李永莉等, 1999, 云南地区水氡前兆异常动态演化与地震关系研究, 地震研究, 22(4): 365-371. |
米雅奇金B И, 1983, 地震孕育过程, 冯德益、顾瑾平, 译, 北京: 地震出版社.
|
邱泽华、张宝红、池顺良等, 2010, 汶川地震前姑咱台观测的异常应变变化, 中国科学:地球科学, 40(8): 1031-1039. |
杨兴悦、王燕、缑亚江等, 2011, 甘肃东南部水氡中期异常特征研究, 地震研究, 34(1): 8-18. |
叶秀薇, 2004, 粤闽地区地下流体从属函数异常与地震关系的初步研究, 防灾减灾工程学报, 24(2): 195-201. |
中国地震局, 2014, 地震水文地球化学观测技术规范, 北京: 地震出版社.
|
钟伟、杨宝俊、张智, 2006, 多项式拟合技术在强噪声地震资料中的应用研究, 地球物理学进展, 21(1): 184-189. DOI:10.3969/j.issn.1004-2903.2006.01.027 |
Huang N E, Shen Z, Long S R, et al, 1998, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc Roy Soc A, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193 |
Manshour P, Ghasemi F, Matsumoto T, et al, 2010, Anomalous fluctuations of vertical velocity of Earth and their possible implications for earthquakes, Phys Rev E Stat Nonlin Soft Matter Phys, 82(3): 036105. DOI:10.1103/PhysRevE.82.036105 |
Manshour P, Saberi S, Sahimi M, et al, 2009, Turbulencelike behavior of seismic time series, Phys Rev Lett, 102(1): 014101. DOI:10.1103/PhysRevLett.102.014101 |
Scholz C H, Sykes L R, Aggarwal Y P, 1973, Earthquake prediction:a physical basis, Science, 181(4102): 803-810. DOI:10.1126/science.181.4102.803 |