0 引言

松辽盆地作为东北地区的主要地震活动构造单元，近几年中强震频发。该地区历史上曾发生1119年前郭6 3/4级地震(吴戈等，1988)，是有记载以来松辽盆地内部发生的最大一次地震。近10年来，松原地区地震活动的频次和强度达到较高的水平，扶余—肇东断裂是东北地区地震活动性最为活跃的断裂构造，2013年在扶余—肇东断裂曾发生前郭M_S5.8震群。除此之外，自2017年7月23日起，该地区先后发生了4次M_S≥4.0地震及小震活跃现象。松原作为东北地区重要的石油开采和经济文化重点地区，研究本地区内中强地震与地震前地球物理观测异常之间的关系，对于减少人民生命财产损失、维护社会稳定具有重要意义。

我国地电阻率定点观测始于1967年，已经发展成为我国地震地球物理场观测的重要分支之一。通过记录测区范围内介质电阻率随时间的变化特征，有助于探索地震的孕育过程。在50多年的观测过程中，积累了上百次发生在观测台网内及附近5~8级地震前突出的中短期异常变化(钱家栋等，1985；国家地震局预测预防司，1998；汪志亮等，2002；杜学彬，2010；Lu et al，2016；解滔等，2018)。地震发生前，震中附近的视电阻率变化表现为：偏离之前多年背景观测值范围的年尺度持续性下降或上升，变化幅度约1%~7%，且通常伴有年变形态畸变，部分地震发生前3个月内出现过加速下降变化；不同方向观测的视电阻率，呈现出与主压应力方位有关的各向异性变化，且变化特征与实验结果相吻合(赵玉林等，1995；钱复业等，1996；杜学彬等，2007)；视电阻率的下降/上升变化形态，与震源机制解给出的台站所处位置的挤压增强/相对膨胀状态有关(赵玉林等，1996；解滔等，2020a)。

本文依据绥化台测区电性结构，采用影响系数理论分析地表浅层介质电阻率变化对观测的影响，讨论绥化台地电阻率反常年变形态的原因。采用断层虚位错模式，将地震同震位移按大小相等、方向相反的方式进行加载，获取地震前震中附近区域的相对介质变形空间分布特征，并结合地电阻率各向异性异常变化，分析绥化台地电阻率异常变化与松原地震之间的关系，为该区域的地电阻率资料分析提供一定的参考。

1 构造环境及台站简介

2017—2019年松原地区先后发生4次中强地震，分别为2017年7月23日M_S4.9、2017年8月15日M_S4.5、2018年5月28日M_S5.7和2019年5月18日M_S5.1地震。这4次地震发生于松辽盆地内部，位于NNE-NE向的扶余—肇东断裂和NW向第二松花江断裂的交汇部位(图 1)，哈尔滨、长春、白城等地震感明显，部分地区房屋被损毁，无人员伤亡(任静等，2020)。

绥化台地处绥化市西郊，海拔179m，位于吉黑褶皱系松辽凹陷东部隆起区，区内第四级盖层约10m，覆盖在白垩系基岩之上，白垩系四方台组为NNE向延伸的向斜，属东部隆起区的绥化凹陷部分(刘长生等，2018；高研等，2019)。该区域主要发震断裂为呼兰河断裂(图 1)，该断裂形成于燕山运动晚期，为正断层，与其距离不远的断层露头处疙瘩山山体为灰岩。

绥化台地电阻率观测于1980年投入使用，目前采用中国地震局研制生产的ZD-8BI数字地电仪器进行观测。测区内布设N15°E和N75°W两条共中心点测道(图 2)，采用四级对称观测装置，供电极距AB均为1000m，测量极距MN均为300m，两测道的中心点距观测室442m。外线路为地埋方式，埋设深度为2.5m，用绝缘程度较高的电线连接电极。观测场地内地势平坦，地形开阔，附近没有大型变电站、高压线路、通讯塔、金属水管、蔬菜大棚等可能对电阻率产生影响的因素，台站数据观测质量在全国评比中一直名列前茅，在2009年安达4.5级和2009年松原4.2级地震前分别出现显著下降和年变畸变异常(刘长生等，2012)。

2 年变化分析

图 3为绥化台两个测道2014—2019年的地电阻率日均值观测曲线，两测道具有同步的趋势变化和清晰的年变化形态。年变特征整体上表现为“夏高冬低”形态，每年2—3月为观测最低值，7—8月降水量增加时观测值处于高值。N15°E测道年变化平均幅度约3.6%(图 3(a))，N75°W测道年变化平均幅约2.5%(图 3(b))。

2.1 地电阻率影响系数理论

地电阻率观测值变化为测区内不同区域介质电阻率变化的加权和。因此，根据测区实际的电性结构，可以计算不同供电极距情况下各层介质对观测的影响系数分布，计算结果有助于综合评估观测系统对测区浅层干扰源的抑制能力以及对深部岩层介质电阻率变化的敏感程度。如果将地电阻率测区划分为任意的n块区域，假设每个小区域的介质电阻率为ρ_i (i=1，2，…，n)，当测区电性结构、观测装置的位置和布极方式确定时，地电阻率ρ_a可表示为各个小区域介质电阻率的函数(钱家栋等，1985)，即

$ {\mathrm{d}}\left(\ln \rho_{\mathrm{a}}\right)=\sum\limits_{i=1}^n \frac{\partial \ln \rho_{\mathrm{a}}}{\partial \ln \rho_i} {\mathrm{~d}}\left(\ln \rho_i\right) $

(1)

在实际应用时，各区域介质电阻率在一定时间内的相对变化幅度较小，即Δρ_i/ρ_i远小于1，因此可以将式(1)进行泰勒级数展开，大于一阶项的部分可忽略不计，最终地电阻率的相对变化可以近似地表示为各个小区域介质电阻率相对变化的加权和，即

$ \frac{\Delta \rho_{\mathrm{a}}}{\rho_{\mathrm{a}}}=\sum\limits_{i=1}^n B_i \frac{\Delta \rho_i}{\rho_i} $

(2)

其中，B_i为影响系数，且满足条件$\sum\limits_{i=1}^n B_i=1$，其表达式为

$ B_i=\frac{\partial \ln \rho_{\mathrm{a}}}{\partial \ln \rho_i}=\frac{\rho_i}{\rho_{\mathrm{a}}} \frac{\partial \rho_{\mathrm{a}}}{\partial \rho_i} $

(3)

2.2 影响系数分析

绥化台已有的电测深资料(张亚江等，2005)显示，测区内NS和EW两条电测深实测曲线宏观上呈现下降趋势，两个方向差异较小(图 4，黑线)，图 4中的红色曲线是以水平均匀层状模型对两条电测深实测曲线进行反演得到的电性结构，NS和EW向电测深曲线大致为Q型，反演结果与实测值基本吻合。根据该反演结果，将测点下介质化简为三层地电断面，如表 1所示。

绥化台实际测道方位与电测深曲线方位存在一定差异，本文中N15°E测道电性结构采用NS方向电测深曲线的结果，N75°W测道采用EW方向电测深曲线的结果。两个测道各层介质影响系数分布情况如图 5所示，在供电极距AB/2小于10m时，观测变化几乎全部受浅层介质的影响，不反映深部电阻率的变化；随着供电极距的增加，地电阻率观测数据开始反映深部介质电阻率的变化信息，在AB/2大于10m且小于100m的区间内，浅层的影响系数逐渐变小，第2和第3层影响系数逐渐增大，主要受第2层介质的电性结构影响；在AB/2大于150m之后，第2层介质影响系数开始随着极距的增加而变小，第3层影响系数逐渐占据主导地位并趋近于1。第1层介质电阻率的影响系数在AB/2大于100m后变为负值，受季节性降雨和温度变化的影响，浅层数米以内介质电阻率在春季之后随温度升高和降雨量的增加呈现下降变化，在秋季之后随温度降低和降雨量的减小呈现上升变化，尤其是浅表土层进入冻土阶段后，电阻率将呈现大幅上升。由于绥化台两个测道在AB/2=500m时，第一层介质的影响系数均为负值，因此年变化表现为“夏高冬低”的反常形态。

3 异常变化分析

从图 3所示观测曲线可以看出，自2017年5月开始N15°E和N75°W两个测道观测数据同步出现年变形态畸变，之后开始出现持续性的下降变化。通常为了分析异常产生的原因，首先要排除观测装置或环境变化的影响，再分析与地震的对应关系。黑龙江地震局和中国地震台网中心针对此次异常开展了多次异常核实，发现异常变化期间台站周围的观测环境无干扰源，观测系统运行正常，测量仪器检定结果合格。

排除装置和环境干扰后，对异常做进一步分析。采用傅里叶滑动方法(赵跃辰等，1984；杜学彬等，2017)对地电阻率原始曲线去年变，结果如图 6所示，截至2018年4月N15°E测道最大下降幅度约1.6%，N75°W测道最大下降幅度约1.4%，之后开始出现转折回升，2019年开始恢复正常年变化状态，且观测数据由之前的趋势下降转变为趋势上升变化(表 2)。

在此次绥化台地电阻率异常时段内，2017年7月23日松原M_S4.9地震发生在年变形态开始畸变后约2个月，之后1个月内又发生了2017年8月15日M_S4.5地震，2018年5月28日松原M_S5.7地震发生在异常达到最低值之后的转折回返阶段，而2019年5月18日松原M_S5.1地震则发生在异常恢复至另一背景变化阶段(图 6)。

4 讨论

为分析地震前地电阻率变化的原因，国内外学者从实验角度明确了应力作用下介质内微裂隙的变化是介质电阻率变化的重要原因(Mjachkin et al，1975)。岩石物理实验结果显示，在压应力作用下含水岩土介质的电阻率呈现下降变化，应力卸载时电阻率呈现恢复上升变化；在临近破裂时多数岩石的电阻率会加速下降，岩石破裂后电阻率则快速回返(Brace et al，1968；Jouniaux et al，2006)；垂直于最大主压应力方向观测的地电阻率变化幅度最大，平行方向变化幅度最小，斜交方向介于二者之间(陈大元等，1983；陈峰等，2013)；野外原地实验结果中地电阻率的变化特征与实验室内结果基本一致(赵玉林等，1983；国家地震局预测预防司，1998)。应力加载超过破裂应力强度的50%之后，岩石的体积开始增大而出现扩容，并伴随有声发射现象，反映出新裂隙的不断产生(Scholz，1967；Brace，1975)。

2017—2019年4次松原中强地震发生在松辽盆地内NW向第二松花江断裂和NNE-NE向扶余—肇东断裂的交汇处(图 1)，全球矩张量(GCMT)和吉林省地震局给出的震源机制解结果显示，4次松原中强地震的震源机制解具有明显的一致性(表 3)，节面Ⅰ为NW走向，节面Ⅱ为NE走向，均为高倾角走滑断层。李君等(2019)认为2018年松原M_S5.7地震的震源机制节面解与第二松花江断裂性质基本一致，由此推断第二松花江断裂为此次地震的发震断裂；而李永生等(2020)利用震源-矩心法的反演结果显示，与NE向的扶余—肇东断裂走向、倾角一致的界面是发震断层面。

为分析绥化台地电阻率2017年开始的下降变化与震前区域应力应变分布的关系，本文以2018年松原M_S5.7地震为例，通过断层虚位错模式计算了震前产生同震位错所需的体应变的变化分布(图 7)。需要注意的是，构造区域内的绝对应力状态通常是难以获取的，采用虚位错模式计算时通常将区域内预应力水平设置为零，计算结果仅表示能引起这部分断层滑动量的应力应变的变化量(赵玉林等，1996；解滔等，2020b)。在整体为挤压构造环境的区域，计算结果中的挤压区域是地震前的挤压增强区域，而对于计算结果中的拉张区域，并不能区分其是真正的拉张区域还是挤压卸载区域，但可认为是相对膨胀区域。由图 7可见，地震前绥化台位于挤压增强区域，地震前的下降异常变化与测区的挤压变形特征吻合。此外，依据表 3中的震源机制解，4次地震的最大主压应力方位分别为81°、90°、264°和74°。N15°E测道与P轴方位的夹角大于N75°W测道，且N15°E测道地震前的变化幅度也大于N75°W测道，两个测道的各向异性变化特征与实验和震例给出结果一致(严玲琴等，2013；贾立峰等，2016；解滔等，2020b)。因此，绥化台地电阻率2017—2019年的异常变化可能与同期松原地区4次中强地震之间存在一定的联系。

图 8展示了松原地区周边部分地电阻率台站2014—2021年期间的观测资料。其中，四平台NW测道自2016年5月开始出现趋势下降变化(图 8(a))，下降幅度约0.5%，远远小于绥化台地电阻率的下降幅度；NE测道的年变形态和年变幅度与往年一致(图 8(b))；白城地电阻率EW和NE两个测道在2014—2021年年变趋势基本保持一致(图 8(c)、8(d))；榆树台观测环境较为复杂，观测数据多次出现因环境干扰造成的快速上升或下降变化，NE和NW两个测道年变形态不规律，2017—2019年4次中强地震前后未出现显著的异常变化(图 8(e)、8(f))；林甸台于2019年9月停测，NS测道自2015年起年变形态规律，变化幅度基本不超过0.7%，4次中强地震前后年变幅度、形态基本没有变化(图 8(g))，EW测道受到环境干扰，观测数据年变形态不规律，震前异常变化特征不清晰(图 8(h))。

综合分析松原地区周边部分地电阻率观测资料在4次中强地震前后的变化特征发现，绥化地电阻率两个测向出现约1.5%的趋势下降变化，四平地电阻率NW出现约0.5%的下降变化，变化幅度远小于绥化地电阻率，其他地电阻率观测数据的年变形态和年变幅度基本保持一致。导致绥化和四平2个台站地电阻率变化幅度存在差异的主要因素推测与区域主压应变率的大小有关：震源机制解(表 3)显示2017—2019年松原地区4次中强地震发震断层的优势走向比较一致，主压应力以近EW向的挤压为主，四平台位于绥化台西南方向，相距约400km；我国大陆及周边GPS观测数据显示，东北块体的主压应力方位角在73.9°N~N75.9°E之间，近EW方向，随着太平洋板块向欧亚板块边界的俯冲，主压应变率自东向西逐渐减小(李延兴等，2006)。

5 结论

本文通过地电阻率边界积分方程和断层虚位错方法，对2017—2019年松原地区4次中强地震前绥化地电阻率异常进行分析，得到以下初步结论：

(1) 绥化台地电阻率观测场地平坦，没有明显的干扰源，观测数据稳定，存在显著的年变形态，能够较好地记录地下介质地电阻率变化情况；

(2) 在供电极距AB/2大于330m后，绥化台测区内表层介质电阻率变化的影响系数变为负值，其是造成观测值“夏高冬低”反向年变的主要原因；

(3) 通过数据处理手段剔除年变化后，震前地电阻率快速下降异常变化更为突出，断层虚位错模式计算结果显示，绥化台位于震前的挤压区域，观测值的下降变化与区域构造应力有关，综合松原地区其他部分地电阻率观测数据，推断2017—2019年期间绥化地电阻率的快速下降变化与4次松原中强震可能存在联系。

致谢: 本文采用Coulomb3.3程序包计算断层虚位错模式的应变分布，审稿专家提出了宝贵的修改意见，在此一并表示感谢。