地震是自然界最具破坏力的自然地质灾害之一,是地球构造运动的一种表现形式。通常强震和大震是由于受到局部应力场以及周边活动断裂的影响而形成应力高度积累,从而在应力高度积累部位及其附近地区形成的明显构造差异运动,往往还伴有重力场显著变化(滕吉文等,2008; Chen et al,1979; Kuo et al,1999; Zhu et al,2010)。在地震地球物理观测与分析的学科领域中,一般会采用定期流动重力观测和连续重力重复观测的测量技术,来捕捉整个设定测网覆盖区内与场源变化相关的震前重力场变化信息(陈运泰等,1980; 祝意青等,2009a; 陈石等,2011)。地下介质变化过程在地表场兆反映中的重力场表现较为明显,且其物理意义明确,在构造运动中是一种较好的表现形式(刘芳等,2016; 高倩,陈石,2015)。祝意青等(2001,2003,2008,2009b,2013,2015)、陈石等(2011,2014)和Chen等(2016)将前人研究成果应用于我国川滇地区、新疆地区以及青藏高原东北部的中长期地震危险分析研究中,并捕捉到震前重力变化高梯度变化特征,为地震危险区的判定提供一定参考。

对于所获取的地表观测重力信号,如何开展其异常的定量化研究,一直以来都是研究人员的研究重点。焦新华和吴燕冈(2009)引入地球物理反演方法对重力场变化的场源特征进行研究,基于“以场溯源,场源结合”的分析思路对已知场源变化信息进行定量化解释。在反演过程中需要精确分离地球物理变化过程中蕴含的复杂信息,来获取有效观测和精确的先验信息以及参数,这一直以来也是重力场变化反演和解释的难点(高倩,陈石,2015; 刘芳等,2016)。

大地测量学的分析方法是地震重力场分析的常用方法,本文尝试采用勘探地球物理学中的欧拉反褶积方法,对2012—2018年辽宁地区观测的流动重力实测数据进行差分变化分析。结合研究区实际地质情况合理设计理论模拟的地表观测重力场模型,对观测到的流动重力场变化进行正演计算,选用欧拉反褶积方法反演场源参数,反复计算选取合适的构造指数和参数后,对2012—2018年实测数据进行反演和解释。

欧拉反褶积方法是由英国地球物理学家Reid等(1990)结合Thompson(1982)对欧拉齐次方程的研究总结后提出的方法,该方法可对重力位场数据进行快速反演,在先验信息较少的情况下,在自动、半自动扫描后确定场源位置,找出构造范围并推算其异常源信息(范美宁,2006),计算公式为:

(x-x₀)[SX(]T[]x[SX)]+(y-y₀)[SX(]T[]y[SX)]+(z-z₀)[SX(]T[]z[SX)]=N(B-T)(1)

式中:(x,y,z)为观测点位置;(x₀,y₀,z₀)为场源位置; T为位场异常; N为构造指数; B为背景场。

构造指数N反映的是位场异常伴随场源深度变化的衰减率,N参数的选取是场源位置确定的关键,一般与场源几何结构有关。结合前人研究成果分析重力异常发现,N一般选择0～2较为合适(范美宁,2006; 刘芳等,2016; 高倩,陈石,2015; 鲁宝亮等,2009)。

结合勘探地球物理研究中基于欧拉反褶积所建立的超定方程组,计算出最佳构造指数N,选取该构造指数进行欧拉反演计算成图。具体步骤如下:①对离散数据进行网格化处理; ②计算在x,y,z 3个方向的梯度的位场异常; ③根据异常体性质选择合适的N和滑动窗口的大小且窗口滑动大小,最小为3×3以上; ④解方程组求出异常场源参数等,并反复调整到最佳反演结果。

欧拉反褶积方法一直以来在勘探学中应用较广,它是在场源信息未知时就可对场源参数进行反演和解释。计算要点主要是在反演时,最佳构造指数N与在x,y,z 3个方向的导数组成线性方程组,计算求解即可。N在重力反演计算中与场源体几何形状相关,实际应用时要根据复杂的地质情况和场源形态多次试算出合适的N。再根据测点间的间距和异常形态特征等因素选取合理的滑动窗口,窗口大小能够覆盖场源异常,反复调整参数直至得到最佳反演结果(刘芳等,2016; 高倩,陈石,2015)。

根据上述对欧拉反褶积方法的认识和应用需求,结合辽宁地区重力观测网络形态和规模,设计符合研究区观测规模的理论模型,该模型覆盖上地壳、中地壳和下地壳不同深度位置场源的3个直面正六面体模型。该模型如图1所示,在图中矩形框为3个场源模型在地表观测网上的水平投影,理论异常范围在0～120×10^-8 m/s²,测网点间距设为5 km,测网长度在X,Y两方向上为180 km。其中M1模型位于中地壳,参数设为:X(50,60),Y(100,110),Z(10,13),中心埋深11.5 km; M2模型位于下地壳,参数设为:X(80,110),Y(40,70),Z(20,25),中心埋深22.5 km; M3模型位于位于上地壳,参数设为:X(80,110),Y(150,155),Z(5,6),中心埋深5.5 km,密度变化选取为3.0×10^-3g/cm³。3个场源模型深度不同,具体参数信息如表1所示。

图1 不同深度场源组合模型及正演异常
Fig.1 The combined model of multi-depth bodies and its forward anomaly

表1 理论场源模型参数
Tab.1 Model parameters of theory field source

通过上述设定模型计算分析后发现,当选取构造指数N=1、滑动窗口选5～8个测点距离时,反演误差深度在15%以内,效果较好,并得到如图2所示的最佳反演结果。图中共计216组解,实际模型位置用矩形框表示,反演模型的结果用黑色圆点表示,计算结果聚集程度较好,将不同深度场源特征明显区分出来。M2和M3模型在水平位置上误差较小,结果较为理想,而M1误差较大。由于上述模型试算结果较好,因此采用该参数对辽宁流动重力测网实测重力变化资料进行反演计算。

图2 基于欧拉反褶积的反演结果三维空间分布
Fig.2 The 3D inversion results using Euler Deconvolution method

辽宁地区处于太平洋板块和欧亚大陆板块交汇的构造复杂部位,属华北地台北部、赤峰—开原深大断裂(内蒙地轴)以南的地台区。著名的郯庐地震带呈NW向、沿沈阳—海城—辽东湾—渤海东部—莱州湾—山东潍坊一线贯穿整个研究区。华北地台、蒙吉—兴安造山带、郯庐断裂带基本控制着研究区的构造背景。辽宁地区的华北地台部分自东向西依次分布有胶辽台隆、下辽河断陷、燕山台褶带。辽宁地区流动重力监测是我国东北部地区对郯庐地震带中北段前兆监测的常规地球物理方法之一。

2012年以来,辽宁地区重力综合测网每年观测2期,观测时间为每年5月和9月,承接101个重力测点、111个重力测段。根据学科发展的规划,对辽宁地区地震重力区域测网进行调整、优化和改造,2018年以来的2期测量包括69个重力测点和84个重力观测测段,测点分布见图3。

图3 辽宁地区测点及断裂分布图
Fig.3 Distribution map of surveys points and fault in Liaoning area

利用2012年9月至2018年9月辽宁地区共计13期流动重力观测数据进行分析,图4给出了1年尺度重力变化差分结果。由于重力测点分布不均匀,笔者对其进行了空间差值。

(1)2012年9月至2013年9月重力分布图显示(图4a),辽宁地区整体以重力正变化分布为主,在辽宁南部岫岩—海城老震区和辽宁西部朝阳—北票地区出现明显重力四象限现象,NW向和NE向有2个明显的重力高梯度带展布,辽宁东部地区尤其在沈阳—海城附近重力有自西向东由负转正的趋势变化。

(2)2013年9月至2014年9月的重力变化显示(图4b),与上一年度相比全区发生正负重力异常反向变化,辽宁西部地区朝阳一带重力异常由正变负尤为明显,该区朝阳—北漂断裂、阜新—锦州断裂梯度变化较大; 辽东地区重力负值面积扩大。

(3)2014年9月至2015年9月重力场分布显示(图4c),辽宁地区整体出现大范围负变化,最

图4 2012—2018年辽宁地区1年尺度重力变化
Fig.4 One-year scale gravity chaige in Liaoning from 2012 to 2018

大负值区出现在朝阳—北票断裂附近,可达-40×10^-8 m/s²重力异常区; 辽东地区由负转正,发生正负变化较大说明断裂断层活跃有所增强。

(4)2015年9月至2016年9月,辽宁北部沈阳及其周边重力场由负变正,辽西地区出现明显重力四象限现象,变化最大处在朝阳—北票断裂附近,达30×10^-8m/s²,辽南海城—岫岩老震区亦出现重力四象限现象(图4d),且其周边地震活动增加,表明老震区附近断裂带有所活动。

(5)2016年9月至2017年9月,与上一期相比重力负值区增大(图4e),总体重力变化较为平稳,说明断裂活动稳定。

(6)2017年9月至2018年9月,在辽宁地区测网整体改造后,测点间距变大等值线稀疏,重力整体变化平稳,辽东出现重力高值变化区,测区重力变化等值线走向大体与上期一致(图4f)。

本文选用欧拉反褶积方法对辽宁地区7期年尺度重力变化做反演,计算得到如图5所示的重力变化的场源形态特征。图中圆圈表示反演的场源位置,圆圈大小与场源深度一一对应,彩色实心圆圈直径从小到大依次分别代表场源深度0～5 km,5～10 km,10～15 km,空心彩色圆圈代表场源深度15～25 km。选取红、粉、黄、粉红、蓝、紫这6种颜色表示年尺度变化反演解,在反演计算

图5 辽宁地区基于2012—2018年重力变化反演的场源特征空间分布
Fig.5 Spatial distribution of field source characteristics based on 2012—2018 gravity inversion in Liaoning area

时,以20 km为网格间距,以8倍间距为窗长,对2012—2018年共计7期年尺度流动重力数据计算,得到反演后位于地壳内部的场源深度为10～50 km。2012—2018年研究区内共发生5级以上地震1次,4级以上地震7次,3级以上地震48次。从反演场源位置分布整体特征上看,引起重力场变化的场源体位置落在6次差分结果上,其中重合度较好的在朝阳—北票断裂和锦州—阜新断裂带区域、场源深度在20～40 km,和营口—海城断裂(郯庐断裂带中北段)和庄河断裂带的条状区域、其场源深度集中在20～50 km。

通过欧拉反褶积方法的计算,在理论模型分析基础上,选用最优解参数对辽宁地区流动重力数据重力场变化进行反演计算,得到与场源位置相关的几何参数,获得如下结论:

(1)欧拉反褶积方法在位场资料解释中应用较广,经过试算和分析发现该方法可反演重力场异常的场源三维空间位置,对多期重复观测重力资料反演,发现持续聚集的位置区域与区域构造活动性增强区域一致。

(2)根据构建的理论模型和实验结果,欧拉反褶积可以选择合适的构造指数和滑动窗口大小,不需要给出准确场源的物性参数即可以得到较合理的反演结果。借助欧拉反褶积反演,得到的场源解位置集中的区域,往往更加接近真实解区域位置。通过计算辽宁地区7期年尺度重力场差分变化发现,重合度较好的地区集中在朝阳—北票断裂附近NE向、海城—岫岩断裂周围NE向区域,场源深度集中在20～40 km,这个位置与近年来辽南及沿海地区和辽蒙交界地震活动频繁活动地区相对应。

(3)辽宁临海海域内并无实测点,在数据差值计算和反演时会产生虚假解,因此对于海域内反演结果深度及位置一致性较差。由于场源变化具有随机性,进而部分地区集中度较低,虽然短期内存在一定的场源运动(物质运移),但连续性较差。

综上认为,采用欧拉反褶积方法对辽宁地区重力场变化进行分析,所识别的场源特征对于区域潜在危险性分析是有意义的。