戴 勇 高立新 杨彦明 魏建民 格 根
(内蒙古自治区地震局,呼和浩特 010010)
1966年邢台7.2级地震发生后,物探电阻率法作为重要的前兆监测手段被引入中国地震监测和预测领域。在1976年唐山7.8级、1976年松潘-平武7.2级、1988年澜沧-耿马7.6级和2008年汶川8.0级等一系列强震发生前均监测到电阻率的异常变化(钱复业等,1980; 赵和云等,1982; 钱家栋等,1985; Zhaoetal.,1996; 汤吉等,1998; Luetal.,2004,2016; 张学民等,2009; 杜学彬,2010; 中国地震局监测预报司,2010; Zhaoetal.,2011)。中国地震台网的地电观测主要采用对称四极装置(fixed-electrode quasi-Schlumberger arrays),供电极距一般约为1km,按均匀介质中常规物探电法估算其探测深度在0.705km以内。由于设置于地震台站用于监测前兆信息的电阻率反映了布极区下方探测范围内整体的电性特征,故该物理量又被称为视电阻率或地电阻率(赵和云等,1982; 姚文斌,1989; 杜学彬等,2008)。由于供电极距和探测深度较小,受浅层介质的状态影响,通常地电阻率的观测曲线存在年变、日变和阶跃等,这些干扰使得地震前的异常变化往往不明显,甚至被湮没(汪志亮等,2002; 王兰炜等,2011; 解滔等,2013,2015,2016; 石富强等,2014; 张国苓等,2015)。因此,总结地电阻率的典型变化并分析其产生原因,对于甄别异常特征、提取前兆异常、评价预测效能和震情判定等具有十分重要的意义。
图1 宝昌台地电布极区的剖面电阻率分布(魏建民等,2019)Fig. 1 Cross-sectional resistivity distribution at Baochang station survey area(after WEI Jian-min et al.,2019).
宝昌台位于内蒙古锡林郭勒盟太仆寺旗宝昌镇,地理坐标为(41.9°N,115.3°E),地质构造位于内蒙地轴东段,属四级构造单元,该区的主要断裂为华北断块区的北部边界断裂——赤峰-开原断裂(内蒙古自治区地震局,2006)。自1980年起宝昌台开始地电阻率观测,其观测曲线存在长期变化、年变和日变等典型变化,并在1989年大同-阳高6.1级及1998年张北6.2级等地震前出现异常变化,在中国地电阻率观测数据中具有代表性(高立新等,1999; 汪志亮等,2002; 戴勇等,2013; 徐锡泉等,2014)。本文通过反演方法获得了该台地电布极区地下的水平层状模型和三维有限元模型,在此基础之上获得了影响系数,同时结合数字信号处理及数值模拟等方法对上述变化的成因进行分析,所得结论可为地电阻率异常分析提供参考依据,并对形成异常核实的工作思路具有借鉴作用。
图1 为基于高密度电法反演得到的深度为2.5~101m的电阻率剖面结果。由图1 可见,地电布极区地下的电阻率基本呈现水平分布,大致可分为3层结构(魏建民等,2019)。
沿NS和EW测向对布极区开展垂向直流电测深,结果显示,电测深曲线类型属于KH型(图2),基于高密度电法和电测深结果,结合由位于布极区附近ZK42号钻孔得到的水文综合地质图(魏建民等,2019),采用尝试法(王家映,1998)进行反演,获得了布极区呈水平层状的地下电性结构模型(表1)。
图2 宝昌台的电测深曲线(魏建民等,2019)Fig. 2 Electrical sounding curve at Baochang station(after WEI Jian-min et al.,2019).
表1 水平层状电性结构模型Table1 Horizontal layered electrical structure model
宝昌台地电阻率观测的电流场可视为稳恒电流场(解滔等,2013),根据地电实测情况确定有限元数值模拟的参数: 设置载荷为2A的直流电流,模型边界的电位为0V,单元类型为二维热-电耦合面单元和三维热-电耦合体单元。
三维几何模型: 固定NS、EW测向的电性参数和极距,当水平尺寸为4i000m×4i000m时,地电阻率的有限元模拟值ρs在最底层厚度h4>3i500m后不随h4的增加而变化; 当h4=4i000m时,ρs在水平宽度d>4i000m后基本不随d的增加而变化。由此,确定三维几何模型的尺寸为4i000m×4i000m×4i071.5m。上述模型是基于2018年4月的观测结果通过反演建立的,在具体研究中将根据需要对层厚、层电阻率等参数进行调整,建立适合研究所需的水平层状模型和三维有限元模型。
本文采用电位分布解析表达式和地电阻率滤波器算法计算宝昌台对称四极装置相应的一维影响系数(O′Neilletal.,1984; 姚文斌,1989; 解滔等,2015),结果如图3 所示。B1、B2、B3、B4分别为模型第1、2、3、4层的影响系数,随着供电极距AB增加,B1、B2整体减小,B3、B4整体增大,并在AB达到一定距离后均趋于平缓。由此可知,增大供电极距可有效抑制第1层和第2层的干扰,同时可放大来自深部的电阻率变化的信息,这可为确定地电阻率观测装置参数提供借鉴。
图3 宝昌台各层介质的影响系数Fig. 3 The medium influence coefficients at Baochang station.
宝昌台现有地电观测装置的供电极距AB=560m、测量极距MN=80m,则NS测向各层的影响系数为B1=0.007、B2=0.016、B3=0.955、B4=0.022,EW测向各层的影响系数为B1=0.003、B2=0.012、B3=0.959、B4=0.026。由影响系数结果(图3)可知: 1)第1层和第2层影响系数均为正,浅层电阻率变化与地电阻率变化具有同向变化特征; 2)第3层影响系数达0.9,比其他3层的影响系数均大1个数量级,说明7~71m深度范围内的电阻率变化能够有效地反映在地电阻率的变化中。
表2列出了模型的第1层电阻率、地电阻率和影响系数B1。由此可见,第1层电阻率下降,则其影响系数将增大,地电阻率下降。其中,NS测向变化较快,EW测向变化较缓慢。
表2 宝昌台第1层电阻率、地电阻率、影响系数B1Table2 First layer resistivity,geo-resistivity and influencing coefficient B1 at Baochang station
宝昌台NS、EW测向的地电阻率自1993年至今一直存在长期的下降变化,且其变化速率存在显著的各向异性(图4)。针对地电阻率的长期变化,赵和云(1994)总结了26a以来全国近百个台站的地电阻率观测资料,认为斜型趋势变化一般不对应地震。当前在地电阻率存在的长期单调变化现象的异常信度判定以及是否将其纳入会商判定依据等方面,仍存在着不同的意见。可见,对地电长期变化的成因进行分析具有十分重要的意义。
图4 宝昌台地电阻率的年均值曲线Fig. 4 Annual mean curves of geo-resistivity at Baochang station.
本文对宝昌台1993—2016年NS、EW测向的地电阻率数据进行线性拟合,得到的变化速率分别为 -0.76Ω·m/a 和 -0.19Ω·m/a,24a中累积变化的幅度分别为 -18.24Ω·m 和 -4.56Ω·m。若上述变化由4层电阻率各自贡献,则由影响系数公式计算获得了各层的变化幅度(表3)。
表3 各层电阻率的理论变化幅度Table 3 The theoretical variation of resistivity in each layer
1993—2016年气温、降雨量的累积变化幅度分别为0.20℃和72mm,上述二者变化导致模型的第1层和第2层电阻率的实际变化幅度达不到引起地电阻率长期变化应有的幅度值。另外,由于第4层为石英斑岩层,该层的密度为2.54~2.66g/cm3(鄢泰宁,2014),孔隙度小、含水率低、温度变化小,电阻率的实际变化也达不到引起地电阻率长期变化应有的幅度值。由此判定,第3层的电阻率变化是引起2个测向的地电阻率长期变化的主要贡献源。
使第3层电阻率出现变化的可能因素为应力场作用和层内水位变化。第3层作为含水层,其水位呈现逐年下降趋势,导致层内的电阻率和地电阻率呈逐年上升趋势,这与实测的地电阻率长期下降的变化相反。因此,第3层电阻率的长期下降主要是受台站所在区域应力持续作用的结果。宝昌台所在区域长期受到近EW向主压应力的持续作用(徐菊生等,1999),使介质内部的裂隙走向逐渐沿主压应力方向优势排列,并出现导电通道连通、导电流体(如水)进入或重新分布,导致以垂直主压应力方向的地电阻率NS测向持续下降为主的地电阻率各向异性(杜学彬等,2001,2007)。
图5 宝昌台地电阻率的月均值曲线Fig. 5 Monthly mean curves of geo-resistivity at Baochang station.
图6 宝昌台测区周边的水位数据Fig. 6 Groundwater level data around Baochang station.
宝昌台NS、EW测向的地电阻率均存在冬春高、夏秋低的正向年变形态(图5),本节将对其成因进行逐层分析: 1)图6 为锡林郭勒盟太仆寺旗生态与农业气象监测站提供的2013年1月—2017年8月该区的水位资料,水位埋深主要位于16~17m,呈现冬春埋深小、夏秋埋深大的年周期波动特征。李丹等(2016)基于太仆寺旗2005—2015年的地下水位监测数据分析了该地区地下水位的动态变化特征,结果显示,太仆寺旗地下水位的年动态变化主要受降水及蒸发影响: 春季气候干旱少雨、蒸发作用强烈,地下水位下降迅速; 5月以后降雨逐渐增加,蒸发量开始减少,但地下水用于蒸发的消耗量大于降雨的补给量,地下水位进一步下降; 8—10月,由于汛期降雨的补给及气温降低,蒸发量进一步减小,水位迅速回升; 11月—次年2月为冻结期,蒸发强度很弱,地下水位相对稳定。由反映岩石电阻率与孔隙率、含水性之间关系的阿奇尔公式(张国民等,2001)可知,模型第3层的电阻率由于受到冬春埋深小、夏秋埋深大的地下水位波动影响,其值呈现冬春低、夏秋高的反向年变形态,而由第3层导致的地电阻率年变也应呈现冬春低、夏秋高的反向年变形态,但这与实际观测结果不符。2)第2层的主体位于冻土层之下、含水层之上,温度、水位等对该层的电阻率影响较小,故第2层的电阻率年变幅度小。3)将宝昌台监测的1993年1月——2008年12月的地电阻率月均值数据和同期温度、降雨量数据进行相关性分析,结果显示该台的地电阻率与温度、降雨量之间均存在较为显著的年相关性,且以负相关为主(戴勇等,2013)。宝昌台模型的第1层为含砾粉细砂层,该层是冻土层,受温度、降雨影响大,而一般地区的含砂黏土、砂土的电阻率变化范围为100~1i000Ω·m(刘刚等,2011),由影响系数公式计算得出NS、EW测向地电阻率的变化幅度分别为11.15Ω·m 和2.05Ω·m,而该值与NS测向8.01Ω·m、EW测向3.96Ω·m 的实测年变幅度相近。综合分析认为,由温度、降雨量季节性变化引起的第1层电阻率年变是地电阻率年变的主要贡献源。
采用小波方法进行去噪后的小时值曲线显示,宝昌台的地电阻率存在显著的日变现象,该日变形态与地电场、地磁形态相差甚远,与该台固体潮汐理论值之间不具有相关性,与布极区温度呈显著的负相关(戴勇等,2013)。温度一般通过影响室内观测仪器、室外观测系统(包括外线路和电极)、位于布极区及附近的金属导体和布极区地下土层等,进而引起地电阻率的观测值变化。1)ZD8B系列地电观测仪器的内置取样电阻在0~40℃范围内年的变化范围 ≤0.06% (熊仲华,2006),宝昌台观测室的温度年变化在18~22℃范围内,故认为观测仪器不受温度影响。2)宝昌台地电观测系统的外线路主要采用架空方式,所用电缆为铠装电缆,并用瓷瓶隔离固定电缆的钢绞线与电线杆,且定期对瓷瓶进行清理。在近年的历次检查中,外线路的绝缘性均符合规范要求(中国地震局,2001),温度无法通过外线路对地电阻率的观测值造成影响。3)宝昌台地电观测的NS和EW测向供电电极和测量电极均为铅电极,埋深为3.5m,电极位于冻土层之下的第2层,该层的温度变化较小。4)钢、铝等金属的电阻温度系数数量级为10-3~10-5℃-1(西安交通大学物理教研组,1974),而布极区的温度日变化幅度约为20℃,故铁丝网等金属电阻率的日变幅度接近0Ω·m。5)布极区地下电性结构模型呈4层水平模型,仅第1层(冻土层)受温度影响较大,其电阻率具有与温度呈负相关的特征,若宝昌台的地电阻率日变仅由第1层电阻率的日变所致,则由影响系数公式计算得到的地电阻率NS、EW测向的变化幅度分别为1.24Ω·m 和0.23Ω·m,而2个测向的地电阻率实测日变幅度均在0.5Ω·m 以内,可见理论结果与实测结果基本一致。综合分析认为,由温度日变引起的第1层电阻率日变是地电阻率日变的主要贡献源。
由于观测环境复杂,宝昌台的地电阻率观测值易受降雨、钢绞线埋设、抽水等干扰因素影响形成阶跃等。本节统计了2012—2018年间阶跃频次与降雨量的对应情况(图7, 8),结果显示,累计出现116次阶跃,且冬春频次低、夏秋频次高,其中49次阶跃与测区降雨在时间上吻合。
图7 宝昌台地电阻率的小时值曲线Fig. 7 Hourly value curves of geo-resistivity at Baochang station.
图8 2012—2018年阶跃的月频次分布Fig. 8 Monthly frequency distribution of step variations from 2012—2018.a 总阶跃频次; b 与降雨对应的阶跃频次
图9 宝昌台的地电阻率观测值与降雨量的时序曲线Fig. 9 Time series curves of geo-resistivity and precipitation at Baochang station.
本节以2017年7月1日降雨引起的地电阻率变化案例进行异常成因分析。2017年7月1日12时宝昌台所在区域开始降雨,至当日18时结束。12时,宝昌台NS测向的地电阻率观测值即出现突降,下降幅度为0.32Ω·m; 至当日18时,持续的快速下降终止,累计变化幅度为0.79Ω·m(图9)。
降雨对宝昌台地电阻率观测值的影响是通过测区铁丝网、悬空钢绞线等金属导体与地接触性加强,雨水大范围渗透和小区域汇集3个途径施加的。地电布极区积水坑位于三维影响系数为负的NS测向供电极与测量极之间,雨水汇集并不是引起本次地电阻率下降的主要因素。金属导体与地的接触性加强,相当于低阻体嵌入第1层,造成地电阻率无滞后的突降。测区第1层为含砾粉细砂层,《工程流体力学(第4版)》中给出的细砂的渗透系数为1×10-3~6×10-3cm/s,降雨引起的渗透的水力坡度为1(孔珑,2014),由达西定律可估算出6h内渗透的最大深度为1.2m,这说明自12时降雨开始至18时,雨水正好渗透至第1层的底面。第1层(表层)在降雨后6h内含水量持续增加,导致该层的电阻率持续降低,引起的地电阻率下降幅度为0.20Ω·m。NS测向第1层的影响系数B1为0.007,电阻率为72Ω·m,由影响系数公式计算可知,第1层电阻率的下降幅度达16.15Ω·m。
宝昌台EW测向的地电阻率观测值由2017年8月20日10时的 146.27Ω·m 下降至11时的 143.69Ω·m,下降幅度为1.76%。数值下降期间,电信公司在EW向西供电极附近埋设通信光缆和钢绞线,光缆埋深为1.2m,钢绞线埋深为0.9m,钢绞线的总体走向为NS向,钢绞线距西供电极13m。2017年8月21日西供电极附近约1i200m的钢绞线被抽出,之后EW向观测数据恢复至埋设前的水平。该次突降和突升与埋设和抽出钢绞线在时间上吻合。
属于低阻体的钢绞线位于第1层,由一维影响系数B1可知,其导致的地电阻率下降,初步分析结果与实际观测吻合。由于上述影响系数反映的某层电阻率变化对地电阻率的影响对于干扰定性分析仍显得较为粗略,故在此采用三维影响系数的分布特征进行干扰源影响的定性分析。8月20日电信公司埋设的钢绞线正好位于影响系数为正的区域(解滔等,2015),由此可判定属于低阻体的钢绞线引起的地电阻率变化是下降的,这与实际观测吻合。依据测区水平层状模型构建了三维有限元模型,并通过数值模拟方法计算了钢绞线埋设前后宝昌台地电阻率EW向地电阻率的变化幅度(表4)。结果显示,理论下降幅度为1.37%,略小于实测下降幅度。理论值与实测值之间出现较小差异主要是由于三维有限元模型仅是对实际情况进行简化后的模型,且计算过程出现的误差等因素也会导致两者存在不同。数值模拟结果进一步证明了8月20日地电阻率EW测向的突降变化是由钢绞线引起的。
表4 埋设钢绞线前后地电阻率的实测值及数值模拟值Table4 Geo-resistivity measurements and numerical simulations before and after steel strand installation
2018年10月16日9—17时,在距EW测向西测量极以西10m处开展了抽水实验,其间进行了水位、地电阻率的同步观测。
10月12—16日,布极区温度的逐日间差异性较小且无降雨,其间地电阻率NS测向值无明显变化,抽水对NS测向地电阻率观测没有影响。12—15日,地电阻率EW测向在每天9—17时的平均变化幅度为 0.11Ω·m,16日9—17时变化幅度达 0.21Ω·m,表明抽水对EW测向影响较大。在EW测向测量极与供电极之间且靠近测量极处(三维影响系数为负的区域)抽水时,位于机井下方的水位下降,且快速形成围绕抽水机井的漏斗状水面,局部电阻率快速上升,造成地电阻率快速下降,没有滞后效应。
(1)高密度电法反演结果显示,宝昌台地电布极区地下的电阻率基本呈现水平分布,电测深曲线类型属于KH型。通过尝试法反演确定,该台地电布极区的地下共分为4层,其中第3层为含水层,层深为6.5~71.5m,较上覆和下伏地层整体表现为高导层,当供电极距AB=560m、测量极距MN=80m时,第3层NS、EW向的影响系数均超过0.9,比其他3个层的影响系数均大1个数量级,这说明第3层的电阻率变化能够有效地反映于地电阻率变化中,该层的电阻率能有效携带应力变化、震前异常变化等信息。
(2)宝昌台NS、EW测向的地电阻率自1993年至今一直存在长期下降变化,且变化速率存在显著的各向异性,这主要是由第3层受台站所在区域的应力持续作用引起的。
(3)宝昌台NS、EW测向的地电阻率均存在冬春高、夏秋低的正向年变形态,由温度、降雨量的季节性变化引起的第1层电阻率年变是地电阻率年变的主要贡献源; 位于第3层的水位年周期波动虽然对宝昌台地电阻率年变存在影响,但并非主要因素,这是对高立新等(2019)研究结果的深化和重新认识。宝昌台地电阻率两测向均存在凌晨及上午高、下午及晚间低的正向日变形态,这主要是由温度对表层电阻率的影响所引起的,类似的日变形态也存在于西昌小庙、甘孜、乌加河、青光、宝坻新台等台站地电阻率数据的小时值曲线中。
(4)宝昌台地电阻率阶跃存在冬春频次低、夏秋频次高的特点,且多与测区内降雨、短期抽水、埋设钢绞线等吻合,通过数值模拟、三维影响系数计算等方法对典型阶跃进行定性和定量角度分析的结果也印证了这一统计结果。
致谢中国地震局兰州地震研究所杜学彬研究员在成文过程中给予了悉心指导; 中国地震台网中心解滔副研究员为本研究提供了部分计算程序; 内蒙古自治区宝昌地震台贾昕晔高级工程师、高昌志高级工程师为本研究提供了基础资料; 审稿专家为本文提出了宝贵的修改建议。在此一并表示感谢!