基金项目:中国地震局地球物理勘探中心青年基金项目(YFGEC2023008).
第一作者简介:李勇江(1988-),工程师,主要从事相对重力数据采集分析工作.E-mail:651451568@qq.com.
(中国地震局地球物理勘探中心,河南 郑州 450002)
(Geophysical Exploration Center,China Earthquake Administration,Zhengzhou 450002,Henan,China)
the Shanxi Fault Subsidence Zone; the Bayesian adjustment method; the classical adjustment method; time-varying gravity field; mobile gravity
DOI: 10.20015/j.cnki.ISSN1000-0666.2025.0066
重力场是地球的基本物理场,其静态和动态特征反映了地球表层及内部物质的分布变化和运动状态。地震的孕育和发生经常伴随着地球重力场的变化(祝意青等,2015; Montagner et al,2016; 陈石等,2015; 陈兆辉等,2018; 王青华等,2020; 赵云峰等,2023; 李树鹏等,2024; 杨雄等,2024),其变化的范围与量级和地震震级存在密切联系。因此,定期复测高精度陆地重力场是一种有效且重要的中长期前兆监测手段,可以应用于地震孕育、发展、发生过程的监测与研究中。其工作流程一般为:在测区内每年进行1~2期流动重力观测,之后在绝对重力点的控制下,通过平差计算得到每个测点的重力值,从而得到半年尺度、一年尺度及多年尺度的重力场时空变化,进而分析提取地震前兆异常。因此,合适的平差方法是处理流动重力观测数据,得到可靠结果的关键。
目前最常使用的是经典平差法(Classical adjustment method,CLS),这是一种以最小二乘拟合为原理的线性平差方法。然而,目前用于流动重力观测的主流仪器为LCR-G型、Burris型及CG-5、CG-6型相对重力仪,这些仪器均存在不同程度的非线性漂移,由此产生的误差是经典平差法无法有效消除的。针对这一问题,Chen等(2019)引入了贝叶斯平差方法(Bayesian adjustment method,BAY)替代传统的经典平差法。不同于经典平差法将仪器的漂移率作为常量,贝叶斯平差法将仪器的漂移率光滑假设为已知的先验信息,以测网中若干绝对重力点值作为约束条件,将影响仪器读数的主要参数作为未知参数待定求解,通过贝叶斯原理以及赤池-贝叶斯信息准则(ABIC)选取最优的参数值,在得到平差值的同时获取每台重力仪的漂移特性(Wang et al,2023),此方法已在华北地区及川滇地区多期观测数据处理中得到实际应用(徐伟民等,2021; 黄江培等,2022; Wang et al,2023; Han et al,2024)。
本文使用经典平差法和贝叶斯平差法对2015—2020年第一期山西地区流动重力观测数据进行处理,并对其点值精度、段差残差进行对比,检验贝叶斯平差方法在山西地震重力监测网的实际效果; 以2015-01期数据为基准,即以各测期平差结果减去2015-01期平差结果,得到山西地区1~5年期的重力场时空变化图像,结合2种方法得到的结果,分析总结研究区重力场时空变化特征。
山西地区位于鄂尔多斯块体东缘,主体构造为山西断陷带,整体走向NE-NNE向,总体呈S状,内部由一系列的构造盆地和盆地间的隆起构成,东西两侧分别是太行山与吕梁山,南部与汾渭盆地、秦岭构造带相连,其北部的晋冀蒙交界区是阴山构造带、张家口—渤海地震带的交会部位,构造背景复杂。研究表明,鄂尔多斯块体在青藏高原、四川盆地以及华北地块相对运动的作用下,整体产生了逆时针的转动(李延兴等,2005),这种运动造成了山西断陷带地震活动频繁发生,该区历史上曾发生过1303年侯马8级强震,1695年临汾73/4级强震以及1626年灵丘7级强震。然而,自1999年大同—阳高MS5.6地震后,该区并无5级以上地震发生。
山西地震重力监测网始建于1991年,至1997年测网覆盖山西地区主要构造(冯建林等,2020),2015年测网内测点总数为149个(图1a),2018年将测点优化至55个(图1b),测网网型基本保持不变。自2022年开始,测网由每年2期观测调整为每年1期观测。优化后测网内测点平均间距超过40 km,无法有效监测5级以下地震。因此,本文计算结果反映地更多是研究区地震平静期内的重力场时空演化特征,该结果可作为区域重力场背景变化特征。
本文使用2015—2020年山西地震重力监测网记录的重力观测数据,记录仪器为LCR-G型金属弹簧相对重力仪,时间跨度覆盖测网优化调整前后。在测网内的大同、太原、夏县和长治等4个绝对重力点值控制下,对各期数据分别进行经典平差和贝叶斯平差计算,得到各测期测点重力值,并对2种平差方法进行点值精度和段差残差分析。
其中,经典平差计算使用中国地震局推广的LGADJ程序完成,贝叶斯平差计算使用中国地震局地球物理研究所重力研究团队提供的开源Python库Geoist程序完成。
点值精度是判断平差结果是否准确可靠的重要指标。经典平差和贝叶斯平差方法得到的平均点值精度见表1,由表1可见,各测期2种平差方法的最大平均点值精度分别为10.5×10-8 m/s2和11.3×10-8 m/s2,均达到了较高的精度。除平均点值精度外,还得到各测期测网内各测点的点值精度,但由于篇幅所限,仅列出2016-01期(测网优化前)与2019-02期(测网优化后)使用2种平差方法得到的测网内各测点的点值精度分布。从图2可见,经典平差法(图2a)与贝叶斯平差法(图2b)各测点的点值精度分布整体相似。
图2 2016-01期、2019-01期使用经典平差法(a)与贝叶斯平差法(b)得到的各测点点值精度
Fig.2 Accuracy of the observed values at the monitoring points in 2016-01 termand 2019-01 term obtained through the classical adjustmentmethod(a)and the Bayesian adjustment method(b)
表1 各测期2种计算方法的平均点值精度
Tab.1 Average accuracy of the observed values at the monitoring points in 6 measurement periods obtained through two adjustment methods
为进一步对比2种平差方法在点值精度上的差异,本文统计了2016-01期、2019-01期各测点的点值精度在不同数值区间的测点数量(图3)。
图3显示,2种平差方法的点值精度呈正态分布。2016-01期经典平差法和贝叶斯平差法的点值精度分别为(7~9)×10-8 m/s2(图3a-1)和(8~10)×10-8 m/s2(图3b-1),两者相差约(1~2)×10-8 m/s2,表明2种方法的点值精度均较高,经典平差法略优于贝叶斯平差法。2019-01期经典平差法和贝叶斯平差法的点值精度分别为(6~8)×10-8 m/s2(图3a-2)和(6~7)×10-8 m/s2(图3b-2),贝叶斯平差法略优于经典平差法。从点值精度分布区间和分布形态上来看,2种平差方法的点值精度水平整体相当。
图3 2016-01期、2019-01期经典平差法(a)与贝叶斯平差法(b)的点值精度统计
[JZ]Fig.3 Statistics of the accuracy of the observed values at the monitoring points in 2016-01 term and 2019-01term obtained through the classical adjustment method (a)and the Bayesian adjustment method(b)
由外部因素造成的观测误差是随机的,因此经平差计算得出的段差残差序列应具备随机分布特征。为验证2种平差方法计算结果的非线性成分,分别对2种平差方法计算得到的段差残差序列进行二项式拟合,如图4中蓝色线和红色线所示。从图中可以看出,经典平差法计算结果在各测期中都有非线性成分,特别是在测网优化调整前,由于测点数量较多,非线性特征更为明显; 而贝叶斯平差法计算结果更符合随机分布特征,即更好地消除了由仪器非线性漂移引起的系统误差。
通过对比2种平差方法的点值精度和段差残差,可以得出:①2种平差方法计算结果的平均精度均较高,其主要原因是山西地震重力监测网观测所使用的LCR-G型金属弹簧相对重力仪,相较于CG-5型石英弹簧重力仪,具有更低的漂移率,由仪器非线性漂移带来的误差相对更小,因此2种平差方法计算结果精度相当; ②贝叶斯平差法计算的段差序列更符合随机分布特征,这说明该方法能更好地拟合仪器的非线性漂移,而经典平差法在计算时可能将这一漂移作为误差进行了处理。另外,黄江培等(2023)在云南地区的研究结果显示,由于最小二乘平差的边缘效应,经典平差法计算结果的点值精度在测网边缘和远离控制点的位置相对较差,而贝叶斯平差法的点值精度虽然在空间分布态势上与经典平差法一致,但在上述位置点值精度明显更优。因此在整体点值精度差别不大的情况下,贝叶斯平差法使整个测网平差计算结果的误差分布更均匀,从而提高测网整体特别是在测网边缘及远离控制点位置的测点精度。
基于2种平差方法计算得到的研究区各测期重力点值,再以2015-01期数据为基准进行作差运算,即以各测期平差结果减去2015-01期平差计算结果,得到使用2种平差方法得出的山西地区1~5年尺度重力场动态变化,如图5所示。从图中可见,2种平差方法得到的不同时间尺度重力场变化在形态上保持一致,只是在局部细节上有所不同,且在部分地区重力变化量存在一定差别,如1年期临汾至长治之间区域(图5a-1、b-1)和4年期忻州至朔州之间区域的重力场变化量差异(图5a-4,b-4)较大,最大为4年期累积变化差别达到20×10-8 m/s2。分析其原因为:2种平差方法计算结果之间的差别在中长期重力变化中被累积增大,加之2018年测网优化后测点点距加大,造成测网内可用数据量减少,在测点点值差别不大的情况下,仍可能使部分地区由插值得到的变化量大于其实际变化量。
从图5可见,1年尺度重力场变化特征为:整体以吕梁以北—太原—阳泉为分界线,呈北正南负分布,且变化量在北东、南西方向增大,这与山西断陷带整体构造方向一致,正负异常之间存在重力梯度带,但变化量较小,为(30~40)×10-8 m/s2(图5a-1、b-1); 2年尺度重力场变化特征为:研究区内重力场负变化区域较大,主要集中在吕梁山、运城盆地及临汾盆地区域,正变化区主要分布于太行山区、朔州至大同之间的区域(图5a-2、b-2)。0等值线沿长治向西经临汾后偏折向北,至太原—忻州—朔州一线,整体走向NNE向,沿断陷带内主要构造盆地西边界展布; 3~5年尺度重力场变化特征为:重力场以负变化为主,整体形态相似,但变化量不大,其中太原盆地附近区域重力变化量最小,2种方法所得的重力变化量级均为(10~20)×10-8 m/s2,临汾盆地附近重力变化量小于30×10-8 m/s2,区域内负变化最为显著的区域是运城盆地,5年尺度重力变化量可达-70×10-8 m/s2。研究区内正变化区域零散分布于长治北部的太行山区以及朔州至大同之间,但重力变化量较小,最大变化存在于朔州北部,5年尺度重力变化量为30×10-8 m/s2(图5a-3~a-5,b-3~b-5)。
从不同时间尺度重力场动态变化结果可以看出,山西地区内部重力场变化量整体较小,这说明研究区内重力场中长期变化相对稳定,也符合研究区处于5级以上地震活动平静期这一实际情况。3~5年尺度的重力场变化特征显示,研究区内不同构造单元呈现出不同的变化特征。以北部的忻定盆地和南部的临汾盆地为界,可以将山西断陷带分为北、中、南3个部分,其重力场整体变化特征为中部小、两端大且变化态势相反,即中部的太原盆地、临汾盆地重力场变化较小,特别是太原盆地,在各时间尺度的重力场变化中几乎都处于0值区附近; 北部的朔州至大同之间,以及南部的运城盆地附近区域重力场变化相对明显,且南北呈现出相反的变化态势。郭良迁等(2010)对山西断陷带内应变率和位移场的研究发现,山西断陷带内部各构造单元呈现出了不同的运动特征,即南北两端的位移速度较大,中段较小,且南北两端位移方向相反,断陷带整体围绕太原盆地北部进行逆时针旋转,这一结论与本文得到的重力场变化特征相吻合。
根据地下物质迁移假说,地壳内部物质会在构造应力的作用下发生迁移,从而导致壳内密度发生改变,进而引发该区域重力场发生相应变化(祝意青等,2009)。研究区内南北两端是山西断陷带与其周边构造单元的交会部位,即北部与阴山构造带、张家口—渤海地震带交会,南部与汾渭盆地、秦岭构造带相连,加之山西断陷带是鄂尔多斯块体的东部边界,在不同构造单元的互相作用下,这些区域的构造活动较强,相比而言,断陷带内部则相对稳定,这反映在区域重力场变化上就形成了中部变化小、南北变化大的特征。
本文利用经典平差法和贝叶斯平差法对2015—2020年山西地震重力监测网的流动重力观测数据进行了处理,对2种方法的计算结果进行了对比,并对研究区的重力场时空变化特征进行了分析总结,可以得到以下结论:
(1)2种平差方法的计算结果精度相当,但贝叶斯平差法对仪器非线性漂移的处理效果更好。2种方法所得的不同时间尺度的重力场变化在形态上基本相同,只是在部分地区出现了数值上的差别,这一差别可能是由测点分布稀疏导致的。
(2)山西地区重力场时空变化整体平静,与研究区整体处于5级以上地震平静期这一现状相符。在山西断陷带内部,重力场变化呈现出中间小、两端大且变化趋势相反的态势,反映了研究区两端与周边构造单元交会部位构造运动较强,内部相对稳定的特征。
贝叶斯平差法的提出主要是为了降低由于仪器的非线性漂移造成的系统误差对观测精度产生的影响,对于今后可能出现的具备不同漂移特性的重力仪的观测数据处理有实际应用价值。但受研究区用于流动重力观测的仪器型号所限,本文研究得到的结果仅对LCR-G型金属弹簧重力仪具备一定的参考价值,但此类重力仪本身的非线性漂移并不显著,因此得到的结果差异并不显著。对于如CG-5等有着较为明显非线性漂移的石英弹簧重力仪,这一方法的处理效果本研究并没有涉及。如果能结合不同型号、不同漂移率的重力仪在同一测区内的数据,使用2种平差方法分别进行处理分析,可能会得到更有价值的结果,今后将做进一步研究。
本文使用的贝叶斯平差计算程序Geoist由中国地震局地球物理研究所陈石研究员及重力学研究团队开发并开源,中国地震局第一监测中心陈兆辉高级工程师和云南省地震局黄江培高级工程师在数据处理和论文撰写过程中提供了指导和宝贵意见,在此一并表示衷心感谢。