0 引言

2019年6月17日22时55分四川省宜宾市长宁县(28.34°N, 104.90°E)发生M_S6.0地震, 震源深度16km。统计结果显示, 1900年以来, 震中50km范围内未发生6级以上地震；震中100km范围内发生过2次6级以上地震, 分别为1917年7月31日云南大关6.8级和1974年5月11日云南大关北7.1级地震, 其中云南大关北7.1级地震距本次地震震中约80km。震中100~200km范围内发生6.0级以上地震5次, 均属1935~1936年四川马边6级强震群序列。需要指出的是本次地震发生区自2018年12月以来5级及以上地震活动较为活跃, 其南侧分别于2018年12月16日和2019年1月3日发生兴文5.7级和珙县5.3级地震。

此次地震位于川东南地区, 距离最近的断裂带华蓥山断裂带, 实测为逆断层(何登发等, 2019)。从构造上看, 本次地震发生于四川盆地边缘的长宁背斜构造附近规模不大的次级断层上, 由于该地区活动断层较多, 地质构造复杂, 所以余震相对较为丰富(尹欣欣等, 2019)。

本次长宁6.0级地震前是否存在可探测的异常活动？不少学者围绕这一问题开展过研究。黄辅琼等(2019)观测到在本次长宁地震发生前的2天时间范围内, 川滇地区一些观测台站的水位水温数据存在可识别的显著异常现象；马亮等(2019)在对地震频谱进行系统跟踪的过程中也发现在本次长宁地震前存在明显的频谱偏移现象；而李艳娥等(2019)发现长宁地震前5年华蓥山断裂带存在视应力增强的过程；阮祥等(2008)认为可能是附近的盐矿井注水诱发地震, 张致伟等(2018)对地震发生时展现出的特征进行了分析。此外, 还有专家认为地震发生和区域的构造活动等因素有关(易桂喜等, 2019), 何登发等(2019)认为长宁背斜区的地震主要发生在背斜南翼, 属于沿基底断层下盘断坡部位活动形成的天然地震。而余怀忠等(2019)在2018年底, 利用基于加卸载响应比(LURR)(Yin et al, 1995)的多方法综合预测模型提出该地区存在明显异常(图 1), 相关结果提交至中国地震台网中心2019年度地震趋势研究报告, 说明该地区存在显著应力积累导致的地震危险性。

本文将分别采用地震释放的贝尼奥夫应变(BS)和地下水位(WL)作为响应量进行LURR的计算, 通过多重物理参量检验, 根据LURR演化过程对震前较长时段长宁地区的应力应变场进行分析, 在此基础上, 收集了震源区附近的地壳形变观测资料, 研究对应时段内可能存在的地壳变形过程, 并进一步将该算法应用到其他更多震例的检验中, 提高算法的可信度。在此研究之前, Zhu等(2019)开展地下流固耦合及渗流问题的研究工作, 探讨了地下流场特性及物理传递过程, 为建立震前异常的地下水位加卸载响应比识别方法提供了理论支撑。目前, 本论文所使用的水位LURR异常识别方法尚未见相关研究。

1 LURR基本原理

LURR是通过系统对外界响应的变化判断其应力状态, 从而实现地震预测的目的(尹祥础, 1987、2004)。在岩石实验中, 随着应力的不断增加, 岩石会经历弹性、损伤和失稳过程。当应力较低时, 加载和卸载的响应相当(图 2)；而当应力达到较高水平时, 岩石内部出现的损伤会致使加载和卸载的响应出现差异, 因此通过对比加卸载过程中的响应变化可以判断岩石的应力状态。

可以定义响应率(X)量化描述加卸载过程的响应变化。根据固体介质的本构关系(Wawersik et al, 1971), 假设载荷的增量为ΔP, 相应的响应增量为ΔR, 响应率可以定义为

$ X=\lim\limits_{\Delta p \rightarrow 0}\left(\frac{\Delta R}{\Delta P}\right) $

(1)

当应力较低时, 加载和卸载的响应率为线性；随着应力的增加, 由于介质内部损伤的出现, 导致加载和卸载时的响应率发生变化, 随着应力的增加, 这种差异会愈发明显, 直到最终破坏(Yu et al, 2006)。这样, 加卸载响应比被定义为

$ Y=\frac{X_{+}}{X_{-}} $

(2)

其中, X代表响应率, “+”和“-”分别表示加载和卸载。当处于弹性范围时, X₊=X_-, 因此Y=1;当介质出现损伤时, X₊>X_-, 则Y>1。根据Y值的变化可以判断固体介质的应力状态。

2 长宁地震前的LURR 2.1 以贝尼奥夫应变(BS)作响应量

在地震研究实践时, 通常取一定时间窗和空间窗内地震所释放出的能量作为响应量, LURR被定义为

$ Y_{m}=\frac{\left(\sum_{i=1}^{N+} E_{i}^{m}\right)_{+}}{\left(\sum_{i=1}^{N-} E_{i}^{m}\right)_{-}} $

(3)

其中, E_i为第i个地震事件所释放出的能量, N+和N-分别代表加载和卸载的事件总数。m可以取0~1之间的任意数字, 特别当m=1时, E^m代表释放的地震能量；当m=1/2时, E^m代表BS。为了避免由于地震记录完备性导致的比值强烈波动, 计算时间窗中通常包含了多个加卸载过程。

首先, 利用传统的地震释放的BS作为响应量计算了震前的LURR时间序列, 计算的起止时间为2010年1月1日~2020年1月1日, 选择区域为半径150km内, 计算时间窗为1年, 滑动步长1个月。数据来源于中国地震台网中心(CENC)的地震目录, 将中小震(1.0≤M≤4.0)释放的BS作为加卸载响应, 预测震级M≥5.0(刘月等, 2012)。并通过G-R关系拟合发现2010年以来该区的地震完备震级在1.0左右, 如图 3所示, 表明长宁地区监测能力较好。震源机制取自全球震源机制解^①, 选取公认测定的1组数据作为详细的主断层面参数, 计算库仑破裂应力的内摩擦系数为0.4。详细的LURR演化曲线见图 4。结果显示, 2010~2019年间, 每当LURR高值出现, 均会伴随一定的显著地震事件, 其中在2018年前的几次LURR高值对应的实际震级分别5.2、5.2、5.6、5.1。而在本次长宁M_S6.0地震前, LURR异常尤为显著, 异常值更高, 持续时间更长。自2018年初开始, LURR值由1.2逐渐升高, 并在2019年初达到1.7左右的峰值, 此后异常一直持续至本次地震发生, 之后LURR逐渐发生回落。其中达到峰值时的LURR水平对应于在此期间相继发生的珙县M_S5.6、兴文M_S5.7、长宁M_S6.0以及其他多次M_S≥5.0地震, 而对于本次长宁地震震级水平与对应的LURR高值水平有所偏差的原因将在下文详细讨论。

① http://www.globalcmt.org/.

2.2 以地下水位(WL)为响应量

为了检验上述震前变化的可信度, 进一步利用WL作为响应量计算LURR, 并定义为WL-LURR。通过对比两组LURR曲线, 实现多重物理参量检验, 提升信度。

固体潮加载会使得地下岩石产生变形、破裂, 由于深井水位具有的流动性, 这一过程会被地下井水位记录。当构造应力较小时, 地下岩石发生的是弹性变形, 加载和卸载的水位变化没有差异。随着应力的不断积累, 当超过弹性极限, 岩石出现裂隙并导致容水体积发生变化, WL可能会因此发生变化。随着裂隙不断增加, 差异会愈加明显, 直至灾变发生。另一方面, 根据Kaiser效应(Li et al, 1993), 岩石在加卸载阶段产生的破裂程度不同, 通常加载阶段产生的裂隙要多于卸载阶段(Kurita et al, 1979), 这就会导致在加卸载过程中的水位变化出现差异。所以将平均水位作为响应量引入到LURR的计算, 并定义加载和卸载阶段的平均水位之比为WL-LURR

$ Y=\frac{\left(\sum_{i=1}^{N+} H_{i}\right) /(N+)}{\left(\sum_{I=1}^{N-} H_{i}\right) /(N-)} $

(4)

其中, H_i表示第i个水位值, N+和N-分别代表处于加载和卸载阶段的水位测值总数。

本研究依据的原理为：当应力加载到较高的水平, 系统进入扩容膨胀阶段, 此时, 任何微小的应力变化都可能会造成裂隙的产生；另一方面, 结合Kaiser效应可知, 微破裂更趋向由增加的库仑破坏应力造成, 即在循环的潮汐应力作用下, 加载阶段比卸载阶段引起更多的水位变化。整个演化过程可以简化为图 5, 当构造应力较低时, 岩石处于弹性阶段, 无裂隙产生, 此时未发生扩容现象, 加卸载过程体积无明显差异, 因此无水位差异出现, WL-LURR值接近于1, 表示为图 5 (a)；随着构造应力积累到较高水平, 岩石开始有裂隙产生, 此时发生扩容现象, 导致加卸载过程的水位响应出现差异, WL-LURR>1, 表示为图 5 (b)。

本文选择震中距最小的3口井的水位资料进行WL-LURR分析, 分别是距震中76km、85km和105km的南溪井、泸州井和荣昌华江井, 其深度从100m到300m不等, 在距本次地震震中100km范围内, 除上述南溪井和泸州井外无其他水位观测井分布, 100~150km范围内除荣昌华江井无其他观测井分布, 这3口井的分布情况见图 6, 资料取自于中国地震数据中心^②。

② http://data.earthquake.cn/.

需要说明的是在自然条件下井水位观测会受到气压、仪器和人为活动等影响, 这些均可能会造成一定的误差, 所以通过以下步骤对原始水位数据进行预处理以保证能够清晰观测到固体潮的水位响应, 具体如下。

(1) 去除数据极值：通过依次计算连续的2个水位测点的水位差值, 如果差值大于平均差值的2倍, 则将其判定为异常极值并去除。其原因在于在岩石发生破裂时, 其体积为弹性体积的0.2~2倍(Brace et al, 1966), 对应的水位变化幅度也不会超过这一范围, 如果超出则认为是干扰并予以去除。

(2) 补全数据：将缺失的水位测点数值用线性插值技术进行补全。

(3) 滤波：将补全后的数据用Butterworth滤波器进行带通滤波, 保留频段为12~24h潮汐相关过程的数据。

图 7为预处理前后的水位资料, 并在表 1中给出了3处水位观测井的详细信息。

图 8为使用WL计算的WL-LURR时间序列, 计算时间窗为60d, 滑动步长为30d, 采用的震源参数与使用BS计算LURR时相同。

对比3口井WL-LURR的时间序列, 可以看出在长宁地震发生前的较长一段时间(即2018年前)比值时间序列演化平稳, 均在1.02范围以下波动；而从2018年年初开始比值序列的演化打破之前的平稳状态并开始出现较大幅度的攀升, 超过1.02, 并在震前2个月达到峰值。虽然南溪水位观测井由于仪器改造的缘故, 水位数据只到2019年4月10日, 从图中我们可以看出南溪井水位LURR正处于明显的上升趋势。

必须指出的是, 3口井的WL-LURR异常存在准同步性, 但年变信息明显仍存在一定影响, 特别是荣昌华江的水位LURR异常的年变信息最明显。就本研究而言, 虽然WL-LURR存在年变信息且不可免除, 但是对比其幅度对LURR的趋势判定影响不大。也就是说, 从历年来的WL-LURR变化幅度对比可观测到长宁地震发生前的明显异常。

3 讨论

LURR是一个反映介质本构关系动态变化的灵敏参数, 计算的理论依据通过潮汐应力触发地震来完成, 当构造应力处于较低水平时(例如在弹性阶段内), 微小的应力变化(例如潮汐应力)不足以触发产生地震, 加载和卸载释放的地震BS没有显著差异, 这时的LURR值维持在1左右；当构造应力积累到相对较高水平时, 较小的应力变化就可以触发产生地震。从图 4的LURR演化曲线可以看出, 在大部分时间里LURR处于相对较低状态, 而在地震发生前的1~2年开始攀升并在震前数月达到峰值, 这与Yin等(2000)的观点相一致。事实上, 在以往多次显著地震前由BS得到的LURR时间序列均能观测到一定的异常变化过程(图 4)。LURR反映了区域的应力积累过程, 其异常的幅度、持续的时间尺度应该与地震震级存在一定的关联性。然而, 仔细对比不同地震前的LURR值, 发现本次长宁M_S6.0地震前的异常幅度虽比以往地震要高、持续的时间也更长, 但对比震级却高了很多, 即异常与对应的地震似乎不匹配。这暗示着本次长宁M_S6.0地震可能并非完全由构造应力积累所致, 其他外界因素应该对该地震能够达到如此规模有相当贡献, 具体原因还有待深入研究。

对于将LURR应用于WL中, 张昭栋等(1997)也曾做过类似研究, 然而其仅考虑到体应变加卸载在孕震过程中产生的变化, 对此进行了探讨(即利用井水位潮汐幅度的变化与固体潮体应变的比值作为响应量进行LURR的计算), 未考虑到扩容部分加卸载在孕震过程中所产生的差异。而实际地下岩石在周期应力作用下不仅发生弹性应变, 还会产生由偏应力导致的扩容变形(裂隙的变化与地震的孕育发生直接相关), 而张昭栋等(1997)也未针对扩容变化伴随的水位变化进行系统分析, 因而难以达到规律性的认识。本研究的创新之处在于岩石扩容膨胀(Brace et al, 1966；Scholz et al, 1973)引起的水位变化通过LURR加以解释, 且量化了预测目标。

从图 8中利用WL得出的WL-LURR演化曲线也可以观测到上述变化特征。相比于传统的BS计算LURR, WL具有更好的连续性, 所得WL-LURR波动幅度更小, 演化规律更为稳定, 因此能够更加准确地给出异常的临界水平。从进一步收集到的该地区地壳形变观测资料中发现, 在长宁地震震中距57km的水富水平摆倾斜仪观测中, EW和SN向在LURR出现异常的同期都发生明显转折变化(图 9)。从图 9可以看出倾斜仪南北和东西分量从2018年开始发生了明显转折变化, 这与上述LURR出现异常较为同步, 说明此时该地区地壳应力可能发生了较大程度的积累。结合BS-LURR以及SF倾斜异常的同期出现均在2018年后, 对比WL的WL-LURR时间序列, 其值在上述异常出现前一直维持在1.02以下运行, 而在异常出现后明显高于1.02, 所以可以将1.02作为水位LURR出现异常的阈值, 并认为当LURR>1.02时存在异常, 反之系统则处于稳定状态。

虽然不同井的WL-LURR在震前异常峰值仅为1.07左右, 仅超过稳定值0.07, 变化幅度只有7%, 但临界状态的阈值为1.02, 即变化幅度为0.02, 从变化幅度的角度分析, 其波动程度已经超出数倍(0.07/0.02=3.5), 在计算WL-LURR的过程中主要比较其相对变化情况与变化的同期性。换言之, 重点关注的不应该是WL-LURR变化的绝对值, 而应该是相对变化幅度, 这样, 可以将BS-LURR与WL-LURR以及SF倾斜变化结合起来, 通过多重物理参量检验提高异常的预测效率。

综上所述, 推测长宁地震产生的原因可能包含两个方面, 其一是构造应力积累, 当载荷超过介质的承载水平时, 触发微破裂或小地震的产生, 从而导致LURR的高值异常变化；另一方面, 由于工业活动也可能导致震源区介质本身承载能力降低, 使得同等应力下岩石更容易发生破坏。然而无论何种原因, 使用加卸载响应比方法均可能捕捉到有效的异常特征, 进一步通过多重物理参量检验能够提高预测效率。了解LURR方法的这一特点, 可为我们今后探测大地震发生的时间和位置提供更多的帮助。

4 结论

本文研究结果显示, 长宁地震前存在明显的加卸载响应比异常, 结合大地形变资料, 表明震源区介质存在明显的应力积累过程, 显示长宁地震的发生与构造活动有关。同时, 提出了新的地震前水位异常的加卸载响应比提取方法, 建立了库仑破裂应力、扩容、水位变化与加卸载响应比之间的动态联系, 通过不同种类参数的异常同期发生加深了对本次地震震前应力积累程度的认识。由图 3发现, 2018年前4次明显的LURR异常对应的震级分别为5.2、5.2、5.6、5.1, 根据异常幅度和震级的正比关系推算, 2018年的LURR异常对应的最大震级水平在5.7左右, 而实际震级为6.0, 考虑到震前LURR异常幅度与震级可能存在的不匹配性, 推测其一可能是估算模型本身的误差, 这需要进一步改进；其二是其他因素对该地震的发生也有重要贡献, 需要更多的依据加以支持。然而, 无论如何, 从加卸载响应比和大地形变观测时间序列中的异常可以获知显著的应力积累是从2018年年初开始, 其峰值的出现与长宁M_S6.0地震的发生对应, 而异常的持续发展与期间不断发生的M_S≥5.0地震相一致。而当前LURR的逐渐减小, 表明区域介质正在逐渐趋于稳定。