高级搜索

震源机制一致性的显著性检验方法

郭祥云, 陈学忠, 李艳娥, 隗永刚, 陈丽娟

郭祥云, 陈学忠, 李艳娥, 隗永刚, 陈丽娟. 2019: 震源机制一致性的显著性检验方法. 地震学报, 41(6): 709-722. DOI: 10.11939/jass.20190045
引用本文: 郭祥云, 陈学忠, 李艳娥, 隗永刚, 陈丽娟. 2019: 震源机制一致性的显著性检验方法. 地震学报, 41(6): 709-722. DOI: 10.11939/jass.20190045
shu
Guo Xiangyun, Chen Xuezhong, Li Yan’e, Wei Yonggang, Chen Lijuan. 2019: Significance test of focal mechanism consistency：Taking the foreshock sequence of the MS5.4 Xiuyan earthquake on November 29，1999 as an example. Acta Seismologica Sinica, 41(6): 709-722. DOI: 10.11939/jass.20190045
Citation: Guo Xiangyun, Chen Xuezhong, Li Yan’e, Wei Yonggang, Chen Lijuan. 2019: Significance test of focal mechanism consistency：Taking the foreshock sequence of the MS5.4 Xiuyan earthquake on November 29，1999 as an example. Acta Seismologica Sinica, 41(6): 709-722. DOI: 10.11939/jass.20190045
shu

震源机制一致性的显著性检验方法

  • 1.

    中国北京 100081 中国地震局地球物理研究所

  • 2.

    中国重庆 401147 重庆市地震局

基金项目: 中国地震局地球物理研究所基本业务专项（DQJB19B10）和国家重点研发计划（2018YFC1503405）共同资助
详细信息
    通讯作者:

    陈学忠: e-mail：cxz8675@163.com

  • 中图分类号: P315.3＋3

Significance test of focal mechanism consistency：Taking the foreshock sequence of the MS5.4 Xiuyan earthquake on November 29，1999 as an example

  • 1.

    Institute of Geophysics，China Earthquake Administration，Beijing 100081，China

  • 2.

    Chongqing Earthquake Agency，Chongqing 401147，China

HTML全文

    摘要
  • 摘要: 本文提出用四维空间的欧氏距离DFM来表示不同地震震源机制之间的一致性，并以1975年辽宁海城ML7.3地震序列和1999年辽宁岫岩MS5.4地震序列为例分析了主震与前震和余震的震源机制一致性与DFM值之间的关系，其结果显示，当欧氏距离DFM＜50时，两次地震的震源机制接近。为了对若干次地震组成的一组地震的震源机制一致性进行判定，引入了显著性检验方法。根据陈颙提出的震源机制一致性参数K，以符号检验法和统计检验量Z值检验法对岫岩MS5.4地震前小震的震源机制一致性进行了分析，其结果表明，在临近岫岩MS5.4地震前所发生地震的震源机制的一致性显著，置信度可达98%。
    Abstract: This paper selected the Euclidean distance DFM to measure the focal mechanism consistency between different earthquakes, and analyzed the relationship between the consistency of focal mechanism and DFM value taking the 1975 ML7.3 Haicheng earthquake sequence and 1999 MS5.4 Xiuyan earthquake sequence as examples. The result shows that when DFM is below 50 the two focal mechanisms become extremely consistent. In order to determine consistency of a number of focal mechanisms the significance test method is drawn into. According to the consistency parameter K proposed by Chen Y, we investigated the consistency of focal mechanism for small earthquakes prior to the MS5.4 Xiuyan earthquake using sign test and Z test. The result suggests that the focal mechanisms of small earthquakes occurred near the MS5.4 Xiuyan main shock are significantly consistent with each other with confidence of 98%. Therefore the method proposed in this paper can be applied to practical earthquake prediction effectively.
    • HTML全文

  • 引言

    地球重力场是研究地球结构与性质的重要物理参数，其变化与地球介质密度和地球动力学过程（固体地球潮汐、内部热流、固体和液体之间质量的交换、表面负荷和地震构造运动等）关系密切（孙和平，2004）。地球重力场的深入研究，为解决地球动力学、地震学、大地测量学等多学科的相关问题提供了必要的物理基础信息（Bayoud，Sideris，2003；Hosse et al，2014；王武星等，2014）。

    目前获取地球重力场信息手段趋于多元化，包括地表重力测量、航空重力测量、海洋重力测量、卫星重力测量。各种重力观测手段累积了丰富的重力数据，充分合理利用这些重力数据，构建出一个精确可靠的全球或局部重力场模型是十分重要且颇有意义的工作（Sadiq et al，2010；吴怿昊等，2016；陈石等，2017；Sobh et al，2019）

    全球地球重力场模型通常用球谐函数表达，其解算需要重力观测数据在全球范围内均匀分布。而在局部区域内，因观测环境的限制，地面重力测点在空间上往往呈不均匀分布，经典的球谐函数方法并不适用于局部重力场模型的构建，因此地面局部重力场建模一般通过对离散重力点进行插值来实现。

    空间插值法在重构局部重力场方面已广泛应用（向文，1997；徐遵义等，2010；徐伟民等，2016），常见的空间插值法有最小曲率法、Shepard插值法、反距离加权法、多面函数法等。空间插值方法的选取不仅与观测的物理量有关，还需要能够反映出测区数据的总体特征。由于地球内部的密度异常呈随机性分布等因素，简单的数学插值方法并不能满足高精度地球重力场建模的需求。同时，对含有噪声的重力异常信号，采用空间插值法推估未知测点位置的函数值时，会不可避免地引入噪声影响。以上问题均可制约局部重力场模型构建的精度。研究一种适合局部重力场模型的重构方法是重力数据资料处理需要面对的重要问题。

    最小二乘配置法是二十世纪六十年代提出而后逐渐发展起来的一种函数逼近理论，广泛应用于地表全球卫星定位系统、卫星重力等资料的数据处理。最小二乘配置法应用的关键在于协方差函数的选择和相关参数的确定。局部重力场的协方差函数表达了不同距离重力场量间的相关性，它反映局部重力场特性，这是一些数学类插值方法不具备的优势（文汉江，2000；彭泽辉等，2010）。另外，最小二乘配置方法从原理上可以对噪声信号滤波，按照噪声影响最小原则来确定待估信号（姚道荣等，2008；Wu et al，2017）。

    最小二乘配置是解决噪声影响的地壳形变场拟合和GPS高程拟合问题的一种有效方法（沙月进，2000；Wu et al，2011）。噪声是影响重力场建模效果的一个重要因素。针对局部重力场建模问题中的噪声因素影响，本文利用不同插值方法进行建模，并考虑不同噪声水平下重力场建模效果，以寻求一种更可靠有效的局部重力场建模方法。

    引言

    地球重力场是研究地球结构与性质的重要物理参数，其变化与地球介质密度和地球动力学过程（固体地球潮汐、内部热流、固体和液体之间质量的交换、表面负荷和地震构造运动等）关系密切（孙和平，2004）。地球重力场的深入研究，为解决地球动力学、地震学、大地测量学等多学科的相关问题提供了必要的物理基础信息（Bayoud，Sideris，2003；Hosse et al，2014；王武星等，2014）。

    目前获取地球重力场信息手段趋于多元化，包括地表重力测量、航空重力测量、海洋重力测量、卫星重力测量。各种重力观测手段累积了丰富的重力数据，充分合理利用这些重力数据，构建出一个精确可靠的全球或局部重力场模型是十分重要且颇有意义的工作（Sadiq et al，2010；吴怿昊等，2016；陈石等，2017；Sobh et al，2019）

    全球地球重力场模型通常用球谐函数表达，其解算需要重力观测数据在全球范围内均匀分布。而在局部区域内，因观测环境的限制，地面重力测点在空间上往往呈不均匀分布，经典的球谐函数方法并不适用于局部重力场模型的构建，因此地面局部重力场建模一般通过对离散重力点进行插值来实现。

    空间插值法在重构局部重力场方面已广泛应用（向文，1997；徐遵义等，2010；徐伟民等，2016），常见的空间插值法有最小曲率法、Shepard插值法、反距离加权法、多面函数法等。空间插值方法的选取不仅与观测的物理量有关，还需要能够反映出测区数据的总体特征。由于地球内部的密度异常呈随机性分布等因素，简单的数学插值方法并不能满足高精度地球重力场建模的需求。同时，对含有噪声的重力异常信号，采用空间插值法推估未知测点位置的函数值时，会不可避免地引入噪声影响。以上问题均可制约局部重力场模型构建的精度。研究一种适合局部重力场模型的重构方法是重力数据资料处理需要面对的重要问题。

    最小二乘配置法是二十世纪六十年代提出而后逐渐发展起来的一种函数逼近理论，广泛应用于地表全球卫星定位系统、卫星重力等资料的数据处理。最小二乘配置法应用的关键在于协方差函数的选择和相关参数的确定。局部重力场的协方差函数表达了不同距离重力场量间的相关性，它反映局部重力场特性，这是一些数学类插值方法不具备的优势（文汉江，2000；彭泽辉等，2010）。另外，最小二乘配置方法从原理上可以对噪声信号滤波，按照噪声影响最小原则来确定待估信号（姚道荣等，2008；Wu et al，2017）。

    最小二乘配置是解决噪声影响的地壳形变场拟合和GPS高程拟合问题的一种有效方法（沙月进，2000；Wu et al，2011）。噪声是影响重力场建模效果的一个重要因素。针对局部重力场建模问题中的噪声因素影响，本文利用不同插值方法进行建模，并考虑不同噪声水平下重力场建模效果，以寻求一种更可靠有效的局部重力场建模方法。

    1.   最小二乘配置法原理

    1.   最小二乘配置法原理

    1.1   方法原理

    最小二乘配置法（least squares collocation）的数学模型可以表示为

    $$ {{L}} {\text{＝}} {{AX}} {\text{＋}} { t} {\text{＋}} { n}{\text{，}} $$ (1)

    式中：L表示观测值向量；X代表非随机的参数向量，A为系数矩阵，AX表示观测向量的系统部分；t表示随机信号部分；n表示噪声部分。式（1）为一般的观测模型，观测值L包含了固有参数确定的系统部分、偏离这种系统的信号部分以及观测噪声。

    估值遵循高斯给出的最小二乘和最小方差，并利用拉格朗日乘数法则，可推导出参数向量X和信号部分t为 （Ruffhead，1987）

    $$ {{X}} {\text{＝}} {\left[ {{{{A}}^{\rm{T}}}{\rm{cov}}{{\left({{{L}}{\text{，}} {{L}}} \right)}^{ - 1}}{{A}}} \right]^{ - 1}}{{{A}}^{\rm{T}}}{\rm{cov}}{\left({{{L}}{\text{，}} {{L}}} \right)^{ - 1}}{{L}}{\text{，}} $$ (2)
    $$ { t} {\text{＝}} {\rm{cov}}\left({{ t}{\text{，}} { L}} \right){\rm{cov}}{\left({{{L}}{\text{，}} {{L}}} \right)^{ - 1}}\left({{{L}} {\text{－}} {{AX}}} \right){\text{，}} $$ (3)

    式中cov（L，L）＝cov（t，t）＋cov（n，n），cov（t，t）为观测值的信号间的协方差矩阵，cov（n，n）为噪声协方差矩阵，cov（t，L）为信号与观测值间的协方差。首先由式（2）求得模型的参数向量X，类似于最小二乘平差求解参数。再由式（3）即可求解信号t。

    在特殊情况下，当观测模型中不存在非随机的参数向量时，即A＝0，此时即为最小二乘预估。本文研究的局部重力场重构属于模型推估问题，可将局部重力场信号看成空间平稳随机场的随机量，仅考虑含有观测噪声的情况下，由式（3）可推导

    $$ { t} {\text{＝}} {\rm{cov}}\left({{ t}{\text{，}} {{L}}} \right){\left[ {{\rm{cov}}\left({{ t}{\text{，}} { t}} \right) {\text{＋}} {\rm{cov}}\left({{ n}{\text{，}} { n}} \right)} \right]^{ - 1}}{{L}}. $$ (4)

    假定计算区域的重力异常的协方差函数只与距离有关，若已知该区域协方差函数，则可以计算出各量间的协方差值cov（t，L）和cov（t，t），再利用式（4）进行推估。由于重力信号部分t可在局部区域任意取点并给出最优推估结果，因此可求得研究区域的连续部分的重力场模型。

    1.1   方法原理

    最小二乘配置法（least squares collocation）的数学模型可以表示为

    $ {{L}} {\text{＝}} {{AX}} {\text{＋}} { t} {\text{＋}} { n}{\text{，}} $

    		(1)

    式中：L表示观测值向量；X代表非随机的参数向量，A为系数矩阵，AX表示观测向量的系统部分；t表示随机信号部分；n表示噪声部分。式（1）为一般的观测模型，观测值L包含了固有参数确定的系统部分、偏离这种系统的信号部分以及观测噪声。

    估值遵循高斯给出的最小二乘和最小方差，并利用拉格朗日乘数法则，可推导出参数向量X和信号部分t为 （Ruffhead，1987）

    $ {{X}} {\text{＝}} {\left[ {{{{A}}^{\rm{T}}}{\rm{cov}}{{\left({{{L}}{\text{，}} {{L}}} \right)}^{ - 1}}{{A}}} \right]^{ - 1}}{{{A}}^{\rm{T}}}{\rm{cov}}{\left({{{L}}{\text{，}} {{L}}} \right)^{ - 1}}{{L}}{\text{，}} $

    		(2)

    $ { t} {\text{＝}} {\rm{cov}}\left({{ t}{\text{，}} { L}} \right){\rm{cov}}{\left({{{L}}{\text{，}} {{L}}} \right)^{ - 1}}\left({{{L}} {\text{－}} {{AX}}} \right){\text{，}} $

    		(3)

    式中cov（L，L）＝cov（t，t）＋cov（n，n），cov（t，t）为观测值的信号间的协方差矩阵，cov（n，n）为噪声协方差矩阵，cov（t，L）为信号与观测值间的协方差。首先由式（2）求得模型的参数向量X，类似于最小二乘平差求解参数。再由式（3）即可求解信号t。

    在特殊情况下，当观测模型中不存在非随机的参数向量时，即A＝0，此时即为最小二乘预估。本文研究的局部重力场重构属于模型推估问题，可将局部重力场信号看成空间平稳随机场的随机量，仅考虑含有观测噪声的情况下，由式（3）可推导

    $ { t} {\text{＝}} {\rm{cov}}\left({{ t}{\text{，}} {{L}}} \right){\left[ {{\rm{cov}}\left({{ t}{\text{，}} { t}} \right) {\text{＋}} {\rm{cov}}\left({{ n}{\text{，}} { n}} \right)} \right]^{ - 1}}{{L}}. $

    		(4)

    假定计算区域的重力异常的协方差函数只与距离有关，若已知该区域协方差函数，则可以计算出各量间的协方差值cov（t，L）和cov（t，t），再利用式（4）进行推估。由于重力信号部分t可在局部区域任意取点并给出最优推估结果，因此可求得研究区域的连续部分的重力场模型。

    1.2   协方差函数

    最小二乘配置法中常用的协方差函数有高斯协方差函数

    $$ C\left(d \right) {\text{＝}} C\left(0 \right){\rm{exp}}\left({ - {k^2}{d^2}} \right){\text{，}} $$ (5)

    西尔沃宁协方差函数

    $$ C(d) {\text{＝}} \dfrac{{C\left(0 \right)}}{{1 {\text{＋}} {k^2}{d^2}}}{\text{，}} $$ (6)

    式中C（0）为信号的方差，k为拟合参数，d为两点之间的距离。在应用最小二乘配置法建立区域应变场与速度场问题中，通常利用这两函数作为协方差函数（武艳强等，2009；江在森，刘经南，2010；管守奎等，2015）。

    利用此两协方差函数表示局部重力场协方差函数时，具有一定的局限性。高斯函数和西尔沃宁函数求出的协方差值恒为正值，而实际计算重力场量间的协方差也可能为负值。随着距离增大，重力场量间的协方差值趋于零，即重力场量间更加独立，重力异常因受局部质量分布不均和地区因素影响，协方差值在零值附近以小幅度正负来回摆动。高斯函数和西尔沃宁函数总是从正值趋于零，不能很好地表示实际重力场协方差函数的一些特征。由于不同局部区域因地形起伏变化不同，所反映出的局部协方差函数细节特征也不尽相同；高斯函数和西尔沃宁函数待拟合的参数只有一个，不能很好地拟合实际局部重力场协方差函数的细节特征。另外上述两种协方差函数不能简单地调和扩展到外部空间，不能当作空间协方差函数使用，限制了空间重力数据的使用。

    最小二乘配置法构建局部重力场时，常用Moritz协方差函数和Tscherning/Rapp 协方差函数（以下简称T-R协方差函数）（张皞等，2006；Yildiz，2012；Abd-Elmotaal，Kühtreiber，2013）。Moritz协方差函数用方差、相长度、曲率参数表达局部重力场的函数特征；T-R协方差函数用球谐函数模型表示重力场参数协方差函数，并引进参数拟合实际重力位场观测数据的特性，通过β，R，A，n等四个参数能充分拟合出局部区域重力场协方差函数，还能够表达出协方差函数与径向距离的关系。

    采用最小二乘配置法重构重力场模型的关键在于能够拟合出充分反映研究区域内重力场函数特征的协方差函数参数。刘敏等（2013）研究结果表明，T-R协方差函数适合于重构重力场模型问题，该协方差函数详细推导可参考Moritz （1980），基本推导过程如下：

    重力扰动位协方差函数表示为：

    $$ K\left({P{\text{，}}Q} \right) {\text{＝}} {{M}}\left[({T\left(P \right){\text{，}}T\left(Q \right)} \right]{\text{，}} $$ (7)

    式中T（P）和T（Q）为球面上两点的扰动位，M为球面上协方差值的平均值。由于算子的各向同性，K（P，Q）只与P，Q两点间的距离ψ有关，则

    $$ K\left({P{\text{，}}Q} \right) {\text{＝}} K\left(\psi \right) {\text{＝}} \frac{1}{{4{\text{π} ^2}}}\mathop \int \nolimits_{\lambda {\text{＝}} 0}^{2\text{π} } \mathop \int \nolimits_{\theta {\text{＝}} 0}^{\text{π}} \frac{1}{{2\text{π} }}\mathop \int \nolimits_{\alpha {\text{＝}} 0}^{2\text{π} } T\left({\theta{\text{，}}\lambda } \right)T\left({\theta '{\text{，}}\lambda '} \right)\sin\theta {\rm d}\theta {\rm d}\lambda {\rm d}\alpha {\text{，}}$$ (8)

    式中θ，λ为球面坐标，α为方位角。

    由于扰动位可用球谐函数级数表示，因此K（ψ）也可表示为球谐函数级数式

    $$ K\left({P{\text{，}}Q} \right) {\text{＝}} \mathop \sum \limits_{n {\text{＝}} 2}^\infty {\sigma _n}{\left({\frac{{{R^2}}}{{r{r{\text{′}\!\!\!\!}}}}} \right)^{n {\text{＋}} 2}}{{\rm P}_n}\left({{\rm{cos}}\psi } \right){\text{，}} $$ (9)

    式中，r，r′分别为P，Q的地心向径，σn为阶方差，R为比亚哈马（Bjerhammar）球半径，Pn（cos ψ）为n阶勒让德多项式。由球谐展开理论并顾及球谐函数的正交关系，有

    $$ {\sigma _n} {\text{＝}} \mathop \sum \limits_{m {\text{＝}} 0}^\infty \left({c_{nm}^2 {\text{＋}} s_{nm}^2} \right). $$ (10)

    以上是全球协方差模型。当进行局部重力场建模时，协方差函数需能反映局部重力场特性。式（9）含高阶次的阶方差求和项，而现有的重力位模型展开式均截断至某一阶次。因此，常用阶方差模型表示高于重力位模型的阶方差。本文采用T-R阶方差模型

    $$ {\sigma _n}\left({T{\text{，}}T} \right) {\text{＝}} \frac{A}{{\left({n {\text{－}} 1} \right)\left({n {\text{－}} 2} \right)\left({n {\text{－}} B} \right)}}{\text{，}} $$ (11)

    式中，A为重力异常场的方差。

    当全球协方差函数运用于局部时，采用 “移去−恢复” 方法，同时顾及参考位模型系数误差，则局部扰动位协方差函数模型为

    $$ K\left({P{\text{，}}Q} \right) {\text{＝}} \mathop \sum \limits_2^M {\varepsilon _n}\left({T{\text{，}}T} \right){\left({\frac{{{R^2}}}{{rr'}}} \right)^{n {\text{＋}} 2}}\!\!\!{{\rm P}_n}\left({{\rm{cos}}\psi } \right) {\text{＋}}\!\!\! \mathop \sum \limits_{n {\text{＝}} M {\text{＋}} 1}^\infty {\sigma _n}(T{\text{，}}T){\left({\frac{{{R^2}}}{{rr'}}} \right)^{n {\text{＋}} 2}}\!\!\!{{\rm P}_n}\left({{\rm{cos}}\psi } \right){\text{，}} $$ (12)

    式中，εn（T，T）为参考模型误差阶方差的相关量。通过平滑因子β与位模型误差阶方差联系

    $$ {\varepsilon _n}\left({T{\text{，}}T} \right) {\text{＝}} \beta \mathop \sum \limits_{m {\text{＝}} 0}^n \left[ {{{\left({\delta {c_{nm}}} \right)}^2} {\text{＋}} {{\left({\delta {s_{nm}}} \right)}^2}} \right]{\text{，}} $$ (13)

    式中，δcnm和δsnm分别为重力位模型相应阶次位系数的中误差。

    局部协方差函数的拟合过程中，可由经验协方差按最小二乘法规则拟合得到R，β，A和n参数，从而得知重力场量的协方差，根据式（4）求得未知信号值。

    1.2   协方差函数

    最小二乘配置法中常用的协方差函数有高斯协方差函数

    $ C\left(d \right) {\text{＝}} C\left(0 \right){\rm{exp}}\left({ - {k^2}{d^2}} \right){\text{，}} $

    		(5)

    西尔沃宁协方差函数

    $ C(d) {\text{＝}} \dfrac{{C\left(0 \right)}}{{1 {\text{＋}} {k^2}{d^2}}}{\text{，}} $

    		(6)

    式中C（0）为信号的方差，k为拟合参数，d为两点之间的距离。在应用最小二乘配置法建立区域应变场与速度场问题中，通常利用这两函数作为协方差函数（武艳强等，2009；江在森，刘经南，2010；管守奎等，2015）。

    利用此两协方差函数表示局部重力场协方差函数时，具有一定的局限性。高斯函数和西尔沃宁函数求出的协方差值恒为正值，而实际计算重力场量间的协方差也可能为负值。随着距离增大，重力场量间的协方差值趋于零，即重力场量间更加独立，重力异常因受局部质量分布不均和地区因素影响，协方差值在零值附近以小幅度正负来回摆动。高斯函数和西尔沃宁函数总是从正值趋于零，不能很好地表示实际重力场协方差函数的一些特征。由于不同局部区域因地形起伏变化不同，所反映出的局部协方差函数细节特征也不尽相同；高斯函数和西尔沃宁函数待拟合的参数只有一个，不能很好地拟合实际局部重力场协方差函数的细节特征。另外上述两种协方差函数不能简单地调和扩展到外部空间，不能当作空间协方差函数使用，限制了空间重力数据的使用。

    最小二乘配置法构建局部重力场时，常用Moritz协方差函数和Tscherning/Rapp 协方差函数（以下简称T-R协方差函数）（张皞等，2006；Yildiz，2012；Abd-Elmotaal，Kühtreiber，2013）。Moritz协方差函数用方差、相长度、曲率参数表达局部重力场的函数特征；T-R协方差函数用球谐函数模型表示重力场参数协方差函数，并引进参数拟合实际重力位场观测数据的特性，通过β，R，A，n等四个参数能充分拟合出局部区域重力场协方差函数，还能够表达出协方差函数与径向距离的关系。

    采用最小二乘配置法重构重力场模型的关键在于能够拟合出充分反映研究区域内重力场函数特征的协方差函数参数。刘敏等（2013）研究结果表明，T-R协方差函数适合于重构重力场模型问题，该协方差函数详细推导可参考Moritz （1980），基本推导过程如下：

    重力扰动位协方差函数表示为：

    $ K\left({P{\text{，}}Q} \right) {\text{＝}} {{M}}\left[({T\left(P \right){\text{，}}T\left(Q \right)} \right]{\text{，}} $

    		(7)

    式中T（P）和T（Q）为球面上两点的扰动位，M为球面上协方差值的平均值。由于算子的各向同性，K（P，Q）只与P，Q两点间的距离ψ有关，则

    $ K\left({P{\text{，}}Q} \right) {\text{＝}} K\left(\psi \right) {\text{＝}} \frac{1}{{4{\text{π} ^2}}}\mathop \int \nolimits_{\lambda {\text{＝}} 0}^{2\text{π} } \mathop \int \nolimits_{\theta {\text{＝}} 0}^{\text{π}} \frac{1}{{2\text{π} }}\mathop \int \nolimits_{\alpha {\text{＝}} 0}^{2\text{π} } T\left({\theta{\text{，}}\lambda } \right)T\left({\theta '{\text{，}}\lambda '} \right)\sin\theta {\rm d}\theta {\rm d}\lambda {\rm d}\alpha {\text{，}}$

    		(8)

    式中θ，λ为球面坐标，α为方位角。

    由于扰动位可用球谐函数级数表示，因此K（ψ）也可表示为球谐函数级数式

    $ K\left({P{\text{，}}Q} \right) {\text{＝}} \mathop \sum \limits_{n {\text{＝}} 2}^\infty {\sigma _n}{\left({\frac{{{R^2}}}{{r{r{\text{′}\!\!\!\!}}}}} \right)^{n {\text{＋}} 2}}{{\rm P}_n}\left({{\rm{cos}}\psi } \right){\text{，}} $

    		(9)

    式中，r，r′分别为P，Q的地心向径，σn为阶方差，R为比亚哈马（Bjerhammar）球半径，Pn（cos ψ）为n阶勒让德多项式。由球谐展开理论并顾及球谐函数的正交关系，有

    $ {\sigma _n} {\text{＝}} \mathop \sum \limits_{m {\text{＝}} 0}^\infty \left({c_{nm}^2 {\text{＋}} s_{nm}^2} \right). $

    		(10)

    以上是全球协方差模型。当进行局部重力场建模时，协方差函数需能反映局部重力场特性。式（9）含高阶次的阶方差求和项，而现有的重力位模型展开式均截断至某一阶次。因此，常用阶方差模型表示高于重力位模型的阶方差。本文采用T-R阶方差模型

    $ {\sigma _n}\left({T{\text{，}}T} \right) {\text{＝}} \frac{A}{{\left({n {\text{－}} 1} \right)\left({n {\text{－}} 2} \right)\left({n {\text{－}} B} \right)}}{\text{，}} $

    		(11)

    式中，A为重力异常场的方差。

    当全球协方差函数运用于局部时，采用 “移去−恢复” 方法，同时顾及参考位模型系数误差，则局部扰动位协方差函数模型为

    $ K\left({P{\text{，}}Q} \right) {\text{＝}} \mathop \sum \limits_2^M {\varepsilon _n}\left({T{\text{，}}T} \right){\left({\frac{{{R^2}}}{{rr'}}} \right)^{n {\text{＋}} 2}}\!\!\!{{\rm P}_n}\left({{\rm{cos}}\psi } \right) {\text{＋}}\!\!\! \mathop \sum \limits_{n {\text{＝}} M {\text{＋}} 1}^\infty {\sigma _n}(T{\text{，}}T){\left({\frac{{{R^2}}}{{rr'}}} \right)^{n {\text{＋}} 2}}\!\!\!{{\rm P}_n}\left({{\rm{cos}}\psi } \right){\text{，}} $

    		(12)

    式中，εn（T，T）为参考模型误差阶方差的相关量。通过平滑因子β与位模型误差阶方差联系

    $ {\varepsilon _n}\left({T{\text{，}}T} \right) {\text{＝}} \beta \mathop \sum \limits_{m {\text{＝}} 0}^n \left[ {{{\left({\delta {c_{nm}}} \right)}^2} {\text{＋}} {{\left({\delta {s_{nm}}} \right)}^2}} \right]{\text{，}} $

    		(13)

    式中，δcnm和δsnm分别为重力位模型相应阶次位系数的中误差。

    局部协方差函数的拟合过程中，可由经验协方差按最小二乘法规则拟合得到R，β，A和n参数，从而得知重力场量的协方差，根据式（4）求得未知信号值。

    2.   其它插值法原理

    反距离加权法是常见的插值方法之一。该方法假设距未采样点越近的观测点对未采样值大小影响越大；距离越远，影响越小，影响与距离成反比，即

    $$ H {\text{＝}} \frac{{\displaystyle\mathop \sum \limits_{i {\text{＝}} 1}^n \frac{1}{{{{{{D_i^p}} }}}}{H_i}}}{{\displaystyle\mathop \sum \limits_{i {\text{＝}} 1}^n \frac{1}{{{{{D_i^p}}}}}}}{\text{，}} $$ (14)

    式中，H为推估值，Hi为第i个样本值，Di为该样本距离，p是距离的幂。反距加权法主要依赖于反距离的幂值。H幂越高，内插结果越平滑。反距离加权法的权重公式与实际物理过程无关，无法确定其幂值大小是否合适。

    最小曲率法广泛用于地球科学，其插值基准是生成一个具有最小曲率（即弯曲度最小）且到各样点的Z值距离最小的曲面。在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面。使用最小曲率法时依据最大残差和最大循环次数两参数来控制最小曲率的收敛标准。最大残差通常介于10%—1%之间；最大循环次数与栅格大小有关，通常设置为生成栅格数量的1—2倍。

    2.   其它插值法原理

    反距离加权法是常见的插值方法之一。该方法假设距未采样点越近的观测点对未采样值大小影响越大；距离越远，影响越小，影响与距离成反比，即

    $ H {\text{＝}} \frac{{\displaystyle\mathop \sum \limits_{i {\text{＝}} 1}^n \frac{1}{{{{{{D_i^p}} }}}}{H_i}}}{{\displaystyle\mathop \sum \limits_{i {\text{＝}} 1}^n \frac{1}{{{{{D_i^p}}}}}}}{\text{，}} $

    		(14)

    式中，H为推估值，Hi为第i个样本值，Di为该样本距离，p是距离的幂。反距加权法主要依赖于反距离的幂值。H幂越高，内插结果越平滑。反距离加权法的权重公式与实际物理过程无关，无法确定其幂值大小是否合适。

    最小曲率法广泛用于地球科学，其插值基准是生成一个具有最小曲率（即弯曲度最小）且到各样点的Z值距离最小的曲面。在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面。使用最小曲率法时依据最大残差和最大循环次数两参数来控制最小曲率的收敛标准。最大残差通常介于10%—1%之间；最大循环次数与栅格大小有关，通常设置为生成栅格数量的1—2倍。

    3.   模型试验

    3.   模型试验

    3.1   噪声对最小乘配置方法的影响

    EGM2008模型是目前国际公认的高分辨率地球重力场模型（阶次分别为2 190，2 159）。EGM2008模型在我国大陆的精度与在全球范围内的精度相当，在我国大陆的总体精度为10.5×10−5 m/s2 （章传银等，2009），可以基本反映中国大陆地区的实际重力异常场。在本文的模型试验中，自由空气重力异常数据均基于EGM2008重力场模型系数进行计算，模型截断阶次为180阶。将其作为仿真测试数据。本文选择的研究区域范围为华北地区（29°N—43°N，109°E—124°E），区域网格大小为0.1° （图1）。

    图  1  华北地区自由重力异常场
    图中地震事件为中国地震台网中心记录的2012年5月至2019年1月期间的MS≥4.0地震
    Figure  1.  Free gravity anomaly field in North China
    The figure shows the MS≥4.0 events recorded by the China Earthquake Networks Center from May 2012 to January 2019
    下载: 全尺寸图片 幻灯片

    华北地区地质构造复杂，地震灾害频发，该区内的地震风险一直受到高度关注。地球重力场的动态变化是地震监测预报研究的重要信息，地震的孕育发生与活动断层紧密相关，断层运动既会影响深部物质的迁移运动，又会引起浅部物质的位移运动，而这二者均使局部重力场发生变化。重复重力测量能获取局部重力场的动态变化，利用重力测量结果可反演断层的活动参数（申重阳等，2002）。2009年中国地震局对华北地区主要活动构造带的地震动态重力区域网进行调整、优化、改造，优化后的华北地区相关省（自治区、直辖市）地震局自成体系的重力测网进行了有效连接，这使得该区域大尺度的重力场变化研究成为可能。深入研究华北地区重力场，有助于提高华北地区的地震监测预报能力。图1给出了华北地区180阶自由空气重力异常背景场。

    在研究区范围内加上一个满足正态分布的随机数据，模拟观测背景噪声，其均值为0，标准差为3×10−5 m/s2。实际重力测量中，由于地形和测量条件限制，重力测网测点在空间常呈离散不均匀分布，为模拟实际测点的分布情况，从加噪后的重力异常场网格节点中随机选取50%为模拟观测点，将这些点的重力异常值作为仿真值。利用最小二乘配置法、最小曲率法以及反距离加权法对离散点网格化，恢复区域局部重力场。各方法插值结果如图2所示。

    图  2  不同插值法网格化结果对比
    （a） 仿真重力异常场；（b） 反距离加权法插值结果；（c） 最小曲率法插值结果；（d） 最小二乘配置法插值结果
    Figure  2.  Comparison of meshing results of different interpolation methods
    （a） Simulated gravity anomaly field；（b） Inverse distance weighting interpolation results；（c） Minimum curvature interpolation results；（d） Least squares collocation interpolation results
    下载: 全尺寸图片 幻灯片

    对比图2c与图2d可以看出，最小二乘配置法和最小曲率法计算的结果与理论信号相近，能够较好地恢复实际的重力异常变化趋势。反距离加权法在恢复的模型中拟合了过多的噪声干扰，插值结果较差。为进一步对比不同方法结果之间的差异特征，分别将图2b，c，d所示结果与图2a相减，得到了如图3所示的拟合残差空间分布。

    图  3  不同插值法残差影像图
    （a） 最小二乘配置法；（b） 最小曲率法；（c） 反距离加权法
    Figure  3.  Images of residuals for different interpolation methods
    （a） Least squares collocation method；（b） Minimum curvature method；（c） Inverse distance weighting method
    下载: 全尺寸图片 幻灯片

    由图3可以看出，最小二乘配置法表现出的高频信号最少，剩余残差信号能力最弱，最小曲率法次之，反距离加权法最多。剩余残差信号主要反映了各方法本身局限性和噪声引起的误差大小。统计各插值方法的残差标准差：反距离加权法残差标准差为3.3×10−5 m/s2，最小曲率法标准差为2.7×10−5 m/s2，最小二乘配置法标准差为1.5×10−5 m/s2。因此判断各插值方法建模效果：最小二乘配置法最优，最小曲率法较好，反距离加权法次之。

    对比不同噪声水平下各插值方法的适用性。在图2a的重力异常网格节点中，随机选择50%节点。添加零均值，标准差分别为0，1×10−5，2×10−5，3×10−5，4×10−5和5×10−5 m/s2的满足正态分布噪声信号。图4为不同水平噪声下各插值法恢复重力场的残差标准差。可以看出：在不同大小的噪声下，最小二乘配置法的残差标准差均小于最小曲率法和反距离加权法，网格化后结果更精确，体现出该方法抑制噪声具有良好的效果；噪声水平越大，相较于最小曲率法，最小二乘配置法对噪声的抑制效果也更明显。

    图  4  不同插值法不同水平噪声对应的残差标准差
    Figure  4.  Residual standard deviation corresponding to different levels of noise for three interpolation methods
    下载: 全尺寸图片 幻灯片

    3.1   噪声对最小乘配置方法的影响

    EGM2008模型是目前国际公认的高分辨率地球重力场模型（阶次分别为2 190，2 159）。EGM2008模型在我国大陆的精度与在全球范围内的精度相当，在我国大陆的总体精度为10.5×10−5 m/s2 （章传银等，2009），可以基本反映中国大陆地区的实际重力异常场。在本文的模型试验中，自由空气重力异常数据均基于EGM2008重力场模型系数进行计算，模型截断阶次为180阶。将其作为仿真测试数据。本文选择的研究区域范围为华北地区（29°N—43°N，109°E—124°E），区域网格大小为0.1° （图1）。

    图  1  华北地区自由重力异常场
    图中地震事件为中国地震台网中心记录的2012年5月至2019年1月期间的MS≥4.0地震
    Figure  1.  Free gravity anomaly field in North China
    The figure shows the MS≥4.0 events recorded by the China Earthquake Networks Center from May 2012 to January 2019
    下载: 全尺寸图片 幻灯片

    华北地区地质构造复杂，地震灾害频发，该区内的地震风险一直受到高度关注。地球重力场的动态变化是地震监测预报研究的重要信息，地震的孕育发生与活动断层紧密相关，断层运动既会影响深部物质的迁移运动，又会引起浅部物质的位移运动，而这二者均使局部重力场发生变化。重复重力测量能获取局部重力场的动态变化，利用重力测量结果可反演断层的活动参数（申重阳等，2002）。2009年中国地震局对华北地区主要活动构造带的地震动态重力区域网进行调整、优化、改造，优化后的华北地区相关省（自治区、直辖市）地震局自成体系的重力测网进行了有效连接，这使得该区域大尺度的重力场变化研究成为可能。深入研究华北地区重力场，有助于提高华北地区的地震监测预报能力。图1给出了华北地区180阶自由空气重力异常背景场。

    在研究区范围内加上一个满足正态分布的随机数据，模拟观测背景噪声，其均值为0，标准差为3×10−5 m/s2。实际重力测量中，由于地形和测量条件限制，重力测网测点在空间常呈离散不均匀分布，为模拟实际测点的分布情况，从加噪后的重力异常场网格节点中随机选取50%为模拟观测点，将这些点的重力异常值作为仿真值。利用最小二乘配置法、最小曲率法以及反距离加权法对离散点网格化，恢复区域局部重力场。各方法插值结果如图2所示。

    图  2  不同插值法网格化结果对比
    （a） 仿真重力异常场；（b） 反距离加权法插值结果；（c） 最小曲率法插值结果；（d） 最小二乘配置法插值结果
    Figure  2.  Comparison of meshing results of different interpolation methods
    （a） Simulated gravity anomaly field；（b） Inverse distance weighting interpolation results；（c） Minimum curvature interpolation results；（d） Least squares collocation interpolation results
    下载: 全尺寸图片 幻灯片

    对比图2c与图2d可以看出，最小二乘配置法和最小曲率法计算的结果与理论信号相近，能够较好地恢复实际的重力异常变化趋势。反距离加权法在恢复的模型中拟合了过多的噪声干扰，插值结果较差。为进一步对比不同方法结果之间的差异特征，分别将图2b，c，d所示结果与图2a相减，得到了如图3所示的拟合残差空间分布。

    图  3  不同插值法残差影像图
    （a） 最小二乘配置法；（b） 最小曲率法；（c） 反距离加权法
    Figure  3.  Images of residuals for different interpolation methods
    （a） Least squares collocation method；（b） Minimum curvature method；（c） Inverse distance weighting method
    下载: 全尺寸图片 幻灯片

    由图3可以看出，最小二乘配置法表现出的高频信号最少，剩余残差信号能力最弱，最小曲率法次之，反距离加权法最多。剩余残差信号主要反映了各方法本身局限性和噪声引起的误差大小。统计各插值方法的残差标准差：反距离加权法残差标准差为3.3×10−5 m/s2，最小曲率法标准差为2.7×10−5 m/s2，最小二乘配置法标准差为1.5×10−5 m/s2。因此判断各插值方法建模效果：最小二乘配置法最优，最小曲率法较好，反距离加权法次之。

    对比不同噪声水平下各插值方法的适用性。在图2a的重力异常网格节点中，随机选择50%节点。添加零均值，标准差分别为0，1×10−5，2×10−5，3×10−5，4×10−5和5×10−5 m/s2的满足正态分布噪声信号。图4为不同水平噪声下各插值法恢复重力场的残差标准差。可以看出：在不同大小的噪声下，最小二乘配置法的残差标准差均小于最小曲率法和反距离加权法，网格化后结果更精确，体现出该方法抑制噪声具有良好的效果；噪声水平越大，相较于最小曲率法，最小二乘配置法对噪声的抑制效果也更明显。

    图  4  不同插值法不同水平噪声对应的残差标准差
    Figure  4.  Residual standard deviation corresponding to different levels of noise for three interpolation methods
    下载: 全尺寸图片 幻灯片

    3.2   协方差函数的拟合

    预选一个能够反映局部趋势的函数作为协方差函数，根据观测值采用拟合方法求得所选协方差函数中的待定参数，从而确定这个协方差函数，此过程即为协方差函数拟合。

    选定协方差函数后，协方差模型的参数需用实测数据拟合出来。利用已知数据计算出协方差即为经验协方差。实际观测值通常在离散点上给出，可由经验协方差值由数值积分计算。设每个观测值li代表一个区域Ai，观测值lj代表一个区域Aj。可得到经验协方差：

    $$ C\left({{d}} \right) {\text{＝}} \frac{{\displaystyle\mathop \sum \limits_{i{\text{,}}j} {A_j}{A_i}{l_i}{l_j}}}{{\displaystyle\mathop \sum \limits_{i{\text{,}}j} {A_i}{A_j}}} {\text{，}} $$ (15)

    其判定条件为$d {\text{－}} {{\Delta d}}/{2} {\text{＜}} {d_{ij}} {\text{＜}} d {\text{＋}} {{\Delta d}}/{2}$，其中Δd为选取经验协方差函数离散点间隔。

    本文最小二乘配置法计算在Gravsoft软件（Nielsen et al，2012）上实现。首先利用EMPCOV子程序，实现对局部重力场的经验协方差函数计算；然后利用COVFIT子程序，实现对经验协方差函数的拟合（图5），得到式（12）中T-R协方差函数的最优参数β，R，A，n，从而确定这个局部地区重力场的协方差函数；最后利用式（4）对研究区域任意一点的重力场进行滤波和推估计算，在GEOCOL子程序可实现对式（4）的计算。

    图  5  采用T-R协方差函数所得的最小二乘配置参数拟合结果
    （a） 协方差函数拟合结果图；（b） 参数A对协方差函数的影响；（c） 参数r′对协方差函数的影响；（d） 参数n对协方差函数的影响
    Figure  5.  Least-squares collocation parameter fitted by T-R covariance function
    （a） Result by fitting covariance function；（b） The effect of parameter A on the covariance function；（c） The effect of parameter r′ on the covariance function；（d） The effect of parameter n on the covariance function
    下载: 全尺寸图片 幻灯片

    图5a是协方差函数拟合的结果图，可见其拟合效果较好。拟合出各参数分别为：r′＝−8.17，A＝246，n＝10，β＝2。各参数变化对协方差函数曲线拟合的影响如图5b，c和d所示，对比图5a中拟合出的协方差函数，图5b中分别取A为230，246，260，可以看出，重力异常场的方差A影响着协方差函数的起始值，A值越大，起始值也越大，且随着距离增大，协方差值都趋于相同；图5c中分别取r′为−6，−8.17，−12，可以看出，r′影响协方差函数下降幅度的快慢，基本不影响协方差函数的起始值和最终趋向值；图5d将重力异常阶方差的截断阶数n分别取为7，10，15，可以看出，n越小，协方差下降幅度越小，下降速度也越慢，因为n表示的是重力异常阶方差的截断阶次，n越小，意味着模型拟合过程中有更多的低频信号，重力异常间的相关性更强，即协方差值更大。由于低频信号的影响，随着距离增大，协方差下降变化幅度也较慢。平滑因子β是与重力位系数的误差阶方差相关的量，在本拟合过程中，β基本不影响协方差函数图像的变化。

    应用最小二乘配置进行局部重力场建模，须逼近出符合实际的局部协方差函数，但这需具备密集且均匀分布的实测重力异常，这显然与实际研究问题相矛盾。由经验协方差拟合出的协方差函数，只是对实际局部重力场协方差函数的逼近，很难通过有限重力异常点精确地推估出局部重力场的真实协方差函数，这也是最小二乘配置法误差产生的一个重要原因。另外，当拟合出的协方差函数图形收敛性一致或接近时，最小二乘配置法计算的结果基本相同。

    4.   华北地区实际重力测点试验

    本文以现有重力测点进行模拟试验，基于EGM2008地球重力场模型提取各点180阶自由空气重力异常值作为模型输入测试数据，利用不同插值方法重构华北区域局部重力场。各方法恢复结果的对比如图6所示。在测网内部，各方法均基本反映出华北局部区域重力异常变化，相较其它两种方法，最小二乘配置恢复的局部重力场更舒缓平缓，符合实际重力场变化。

    图  6  基于不同插值法的重力场模型对比
    （a） 自由空气重力异常场；（b） 最小曲率法插值结果；（c） 反距离加权法插值结果；（d） 最小二乘配置法插值结果
    Figure  6.  Comparison of gravity field models bassed on different interpolation method
    （a） Free air gravity anomaly field；（b） Interpolation result of minimum curvature method；（c） Interpolation result of inverse distance weighting method；（d） Interpolation result of least square collocation method
    下载: 全尺寸图片 幻灯片

    为对比不同方法结果的差异特征，分别对图6b，c，d结果与图6a做差，得到各插值方法残差信号影像图，可见最小二乘配置法表现出最佳的插值结果，在测网内部几乎没有残差信号。统计各插值方法的残差标准差：最小二乘配置法为0.8×10−5 m/s2，最小曲率法为1.4×10−5 m/s2，反距离加权法为4.1×10−5 m/s2。根据现有测网分布，理论上利用最小二乘配置法可以恢复出较佳的华北区域180阶重力场模型。值得注意的是，图7a和7b中局部地区（33°N—34°N，111°E—114°E）的重力残差信号较为明显。对比图6a相同区域，易看出此区域重力测点空白，缺失重力数据。空白区域重力信号是由外部重力数据推估得到，由于内部重力场信号有较大变化，外部重力数据插值的结果无法真实反映内部实际重力值，因此存在较大误差。

    图  7  不同插值方法的残差影像图
    （a） 最小二乘配置法；（b） 最小曲率法；（c） 反距离加权法
    Figure  7.  Images of residuals of different interpolation methods
    （a） Least square collocation method；（b） Minimum curvature method；（c） Inverse distance weighting method
    下载: 全尺寸图片 幻灯片

    进一步探究信号受噪声污染时各插值方法对华北区域建模的影响。测网各测点信号施加零均值，标准差分别为0，1×10−5 ，2×10−5 ，3×10−5 ，4×10−5 ，5×10−5 m/s2的满足正态分布噪声信号。同理该模型试验步骤，各插值方法重构后的局部重力场与局部重力异常背景场相减，得到各插值方法的局部残差结果，统计各插值方法的残差标准差，如图8所示。当噪声标准差为4×10−5 m/s2时（噪声与信号比例为4%），最小二乘配置法的残差标准差为2.4×10−5 m/s2，最小曲率法为3.6×10−5 m/s2，反距离加权法为4.8×10−5 m/s2。华北区域局部重力场即使在受到噪声污染的情况下，最小二乘配置法依然具有很好的局部重力场构建能力。

    图  8  各插值法不同水平噪声对应的残差标准差
    Figure  8.  Residual standard deviation corresponding to different levels of noise in each interpolation method
    下载: 全尺寸图片 幻灯片

    4.   华北地区实际重力测点试验

    本文以现有重力测点进行模拟试验，基于EGM2008地球重力场模型提取各点180阶自由空气重力异常值作为模型输入测试数据，利用不同插值方法重构华北区域局部重力场。各方法恢复结果的对比如图6所示。在测网内部，各方法均基本反映出华北局部区域重力异常变化，相较其它两种方法，最小二乘配置恢复的局部重力场更舒缓平缓，符合实际重力场变化。

    图  6  基于不同插值法的重力场模型对比
    （a） 自由空气重力异常场；（b） 最小曲率法插值结果；（c） 反距离加权法插值结果；（d） 最小二乘配置法插值结果
    Figure  6.  Comparison of gravity field models bassed on different interpolation method
    （a） Free air gravity anomaly field；（b） Interpolation result of minimum curvature method；（c） Interpolation result of inverse distance weighting method；（d） Interpolation result of least square collocation method
    下载: 全尺寸图片 幻灯片

    为对比不同方法结果的差异特征，分别对图6b，c，d结果与图6a做差，得到各插值方法残差信号影像图，可见最小二乘配置法表现出最佳的插值结果，在测网内部几乎没有残差信号。统计各插值方法的残差标准差：最小二乘配置法为0.8×10−5 m/s2，最小曲率法为1.4×10−5 m/s2，反距离加权法为4.1×10−5 m/s2。根据现有测网分布，理论上利用最小二乘配置法可以恢复出较佳的华北区域180阶重力场模型。值得注意的是，图7a和7b中局部地区（33°N—34°N，111°E—114°E）的重力残差信号较为明显。对比图6a相同区域，易看出此区域重力测点空白，缺失重力数据。空白区域重力信号是由外部重力数据推估得到，由于内部重力场信号有较大变化，外部重力数据插值的结果无法真实反映内部实际重力值，因此存在较大误差。

    图  7  不同插值方法的残差影像图
    （a） 最小二乘配置法；（b） 最小曲率法；（c） 反距离加权法
    Figure  7.  Images of residuals of different interpolation methods
    （a） Least square collocation method；（b） Minimum curvature method；（c） Inverse distance weighting method
    下载: 全尺寸图片 幻灯片

    进一步探究信号受噪声污染时各插值方法对华北区域建模的影响。测网各测点信号施加零均值，标准差分别为0，1×10−5 ，2×10−5 ，3×10−5 ，4×10−5 ，5×10−5 m/s2的满足正态分布噪声信号。同理该模型试验步骤，各插值方法重构后的局部重力场与局部重力异常背景场相减，得到各插值方法的局部残差结果，统计各插值方法的残差标准差，如图8所示。当噪声标准差为4×10−5 m/s2时（噪声与信号比例为4%），最小二乘配置法的残差标准差为2.4×10−5 m/s2，最小曲率法为3.6×10−5 m/s2，反距离加权法为4.8×10−5 m/s2。华北区域局部重力场即使在受到噪声污染的情况下，最小二乘配置法依然具有很好的局部重力场构建能力。

    图  8  各插值法不同水平噪声对应的残差标准差
    Figure  8.  Residual standard deviation corresponding to different levels of noise in each interpolation method
    下载: 全尺寸图片 幻灯片

    5.   讨论与结论

    本文采用了最小二乘配置法研究华北地区的局部重力场建模问题，比较了最小二乘配置法与反距离加权法、最小曲率法恢复局部重力场的效果，讨论了噪声对最小二乘配置结果的影响以及协方差函数拟合问题，主要研究结论如下：

    1） 反距离加权法由于存在 “牛眼” 效应，易在离散数据周围引入较大的误差，拟合出的重力场曲线不符合实际特征，因此反距离加权法不适合局部重力场模型构建。当噪声与信号比例小于2%时，最小二乘配置法的建模效果略优于最小曲率法；当噪声于信号比例大于2%时，最小二乘配置法的建模精度明显高于最小曲率法。若用曲线斜率k表示噪声对各插值方法误差影响大小，即每10−5 m/s2噪声引入的误差大小，最小二乘配置法k＝0.4，最小曲率法为k＝0.9。因此基于不同水平的背景噪声的重力场异常场，相较于反距离加权法和最小曲率法，最小二乘配置法能够更可靠地恢复出局部重力场模型。

    2） 在局部重力场建模中，T-R协方差函数是可靠的选择。T-R协方差函数中的多个参数使拟合出的协方差函数能展现出更丰富的曲线特征。相较于高斯协方差函数和西尔沃宁协方差函数通过单一参数调节曲线特征，T-R协方差函数能更好地拟合实际经验协方差。良好的拟合结果保证了最小二乘配置最终计算结果的准确性。

    3） 基于现有华北区域重力测网测点，最小二乘配置法理论上能构建出较佳的局部重力场。相较于反距离加权法与最小曲率法，即使重力信号受到噪声污染时，最小二乘配置法依然能有效地抑制噪声，恢复出平缓流畅的重力场曲线。用曲线斜率k表示每10−5 m/s2噪声引入的误差大小，最小二乘配置法k＝0.4，最小曲率法k＝0.6。

    本文最小二乘配置法处理的数据仅是自由空气重力异常数据。最小二乘配置法的另一优势是能同时对多种数据进行联合处理。随着重力观测手段的日益丰富，观测技术的不断提升，局部地区积累了越来越多的各种重力数据，如卫星重力数据、航空重力数据、地面重力数据。最小二乘配置法对多源数据的综合处理能有效地提高局部地球重力场模型的精度，这也是以后的研究方向之一。

    5.   讨论与结论

    本文采用了最小二乘配置法研究华北地区的局部重力场建模问题，比较了最小二乘配置法与反距离加权法、最小曲率法恢复局部重力场的效果，讨论了噪声对最小二乘配置结果的影响以及协方差函数拟合问题，主要研究结论如下：

    1） 反距离加权法由于存在 “牛眼” 效应，易在离散数据周围引入较大的误差，拟合出的重力场曲线不符合实际特征，因此反距离加权法不适合局部重力场模型构建。当噪声与信号比例小于2%时，最小二乘配置法的建模效果略优于最小曲率法；当噪声于信号比例大于2%时，最小二乘配置法的建模精度明显高于最小曲率法。若用曲线斜率k表示噪声对各插值方法误差影响大小，即每10−5 m/s2噪声引入的误差大小，最小二乘配置法k＝0.4，最小曲率法为k＝0.9。因此基于不同水平的背景噪声的重力场异常场，相较于反距离加权法和最小曲率法，最小二乘配置法能够更可靠地恢复出局部重力场模型。

    2） 在局部重力场建模中，T-R协方差函数是可靠的选择。T-R协方差函数中的多个参数使拟合出的协方差函数能展现出更丰富的曲线特征。相较于高斯协方差函数和西尔沃宁协方差函数通过单一参数调节曲线特征，T-R协方差函数能更好地拟合实际经验协方差。良好的拟合结果保证了最小二乘配置最终计算结果的准确性。

    3） 基于现有华北区域重力测网测点，最小二乘配置法理论上能构建出较佳的局部重力场。相较于反距离加权法与最小曲率法，即使重力信号受到噪声污染时，最小二乘配置法依然能有效地抑制噪声，恢复出平缓流畅的重力场曲线。用曲线斜率k表示每10−5 m/s2噪声引入的误差大小，最小二乘配置法k＝0.4，最小曲率法k＝0.6。

    本文最小二乘配置法处理的数据仅是自由空气重力异常数据。最小二乘配置法的另一优势是能同时对多种数据进行联合处理。随着重力观测手段的日益丰富，观测技术的不断提升，局部地区积累了越来越多的各种重力数据，如卫星重力数据、航空重力数据、地面重力数据。最小二乘配置法对多源数据的综合处理能有效地提高局部地球重力场模型的精度，这也是以后的研究方向之一。

  • 图  1   海城ML7.3地震序列的震源机制图像和欧式距离DFM值

    Figure  1.   Focal mechanisms and Euclidean distance DFM values of the ML7.3 Haicheng earthquake sequence

    下载: 全尺寸图片 幻灯片

    图  1   海城ML7.3地震序列的震源机制图像和欧式距离DFM值

    Figure  1.   Focal mechanisms and Euclidean distance DFM values of the ML7.3 Haicheng earthquake sequence

    下载: 全尺寸图片 幻灯片

    图  2   岫岩MS5.4地震前小震震源机制解空间分布

    Figure  2.   Spatial distribution of focal mechanisms of small earthquakes occurred prior to the MS5.4 Xiuyan earthquake

    下载: 全尺寸图片 幻灯片

    图  3   1999年岫岩MS5.4地震前各次小震的欧式距离DFM值变化

    Figure  3.   Variation of Euclidean distance DFM values of small earthquakes prior to the MS5.4 Xiuyan earthquake in 1999

    下载: 全尺寸图片 幻灯片

    图  4   1999年岫岩MS5.4地震前欧式距离DFM平均值随时间的变化

    Figure  4.   Temporal variation of average Euclidean distance DFM value prior to the MS5.4 Xiuyan earthquake in 1999

    下载: 全尺寸图片 幻灯片

    图  5   1999年岫岩MS5.4地震前K值（a）和Z值（b）随时间的变化

    水平线分别表示通过α＝1%，2%和5%的显著性水平检验的临界值

    Figure  5.   Temporal variation of K value （a） and Z value （b） prior to the MS5.4 Xiuyan earthquake in 1999

    Horizontal lines show the critical values at the significance level α＝1%，2% and 5%

    下载: 全尺寸图片 幻灯片

    图  6   以表3中第41号（a）和第36号（b）地震的震源机制分别作为主震的震源机制所得的岫岩MS5.4地震前K值随时间的变化

    Figure  6.   Temporal variation of K value prior to the MS5.4 Xiuyan earthquake taking the focal mechanism of No.41 （a） and No.36 （b） earthquake in Table 3 as that of the main shock

    下载: 全尺寸图片 幻灯片

    表  1   海城ML7.3地震ML≥4.0前震和余震的震源机制（引自顾浩鼎等，1976）及其与主震之间的欧氏距离DFM值

    Table  1   Focal mechanisms for ML≥4.0 foreshocks and aftershocks of the ML7.3 Haicheng earthquake （after Gu et al，1976） and their Euclidean distance DFM value

    序号发震时刻地理坐标ML深度
    /km    		节面Ⅰ节面ⅡP 轴T 轴B 轴　DFM
    　 年-月-日 时：分：秒东经/°北纬/°走向/°倾角/°滑动角/°走向/°倾角/°滑动角/°方位角/°俯角/°方位角/°俯角/°方位角/°俯角/°
    11975-02-0407：50：47.0122.7540.674.7171108525.1 1765174.5 6115 1572230064 18.85
    21975-02-0410：35：35.0122.7840.674.3151098622.1 1768175.7 6113 1551829768 15.66
    31975-02-0419：36：06.0122.8040.657.31229081−15.2 2375−170.7 6617.5157 410072.50
    41975-02-0501：01：45.0122.9340.704.410298884.020886178.02522 343 414386 17.67
    51975-02-0502：56：29.0122.8240.674.510112886.0 2284178.0 663 157 631083 14.64
    61975-02-0512：33：00.0122.7740.684.1101266111.522180150.526828 17113 5959 29.52
    71975-02-0523：52：54.0122.6340.704.610125645.621785153.926422 16814 4764 23.77
    81975-02-0605：43：42.0122.9040.625.22329287−2.0 2288−177.0 673 338 023186 15.11
    91975-02-0612：24：57.0122.5040.805.4171628034.625956167.929531 351614954 132.71
    101975-02-0613：56：16.0122.8340.754.01011680−128.9 1440−15.716941 562530438 147.66
    111975-02-0802：30：23.0122.4740.824.01216178−31.826259−166.029338 371814848 131.25
    121975-02-1220：42：46.0122.7840.704.0 714352−118.2 346−58.917168 73 334122 143.64
    131975-02-1521：08：02.0122.7840.705.4121118426.2 1864173.3 6214 1582230363 18.79
    141975-02-1622：01：26.0122.8040.685.311123626.8 3084151.826015 16324 2062 25.26
    151975-02-1818：51：49.0122.6540.774.2171144414.5 1480133.125223 14239 542 38.94
    161975-02-2215：45：14.0122.7340.704.0121138536.2 1954173.8 6021 1622930054 26.42
    171975-02-2405：07：20.0122.8840.784.4 713257−52.925848−132.927860 17 511030 149.77
    181975-02-2504：52：10.0122.6240.734.41413474−4.222686−164.027614 180 8 6274 159.93
    191975-02-2605：09：53.0122.8240.674.3 81219020.021070180.025414 3471412170 16.62
    201975-03-2111：32：59.0122.9540.774.01126072−29.6 062−159.5 3834 132 723156 41.11
    211975-03-2923：16：36.0122.6040.774.1 61118524.1 1966174.5 6313 1582030365 16.92
    221975-04-1003：55：37.0122.4840.724.610118859.0 2781184.92533 1631035780 18.20
    231975-04-2100：17：06.0122.4540.774.0 81098430.2 1660173.1 5917 1572529959 22.14
    241975-07-0407：06：29.0122.6740.724.11013756−19.423974−144.428336 18512 8052 157.73
    下载: 导出CSV

    表  2   岫岩MS5.4地震MS≥4.0前震和余震的震源机制解（引自张萍，蒋秀琴，2001）及其与主震之间的欧氏距离DFM值

    Table  2   Focal mechanisms for MS≥4.0 foreshocks and aftershocks of the MS5.4 Xiuyan earthquake （after Zhang，Jiang，2001） and their Euclidean distance DFM value

    序号发震日期地理坐标MS深度
    /km    		节面Ⅰ节面ⅡP 轴T 轴B 轴DFM
    　 年-月-日 时：分：秒东经/°北纬/°走向/°倾角/°滑动角/°走向/°倾角/°滑动角/°方位角/°俯角/°方位角/°俯角/°方位角/°俯角/°
    11999-11-0907：01：40.6123.0340.534.19320707.522883159.8261263581120870159.33
    21999-11-0907：07：21.2123.0240.534.2813780−7.122883−169.9 6218331 3 8477158.79
    31999-11-2520：47：48.5123.0040.554.081374429.5 2770130.0 5625150 8 113839.46
    41999-11-2520：55：4.2123.0040.554.491394230.7 2770127.72711216746 123642.50
    51999-11-2623：34：01.0123.0240.534.491107010.6 2480159.72461515212 466717.00
    61999-11-2912：10：39.2123.0340.535.4929684−20.1 2870−173.6 624119136280700
    71999-11-2912：45：50.4123.0340.535.1930070173.62068420.1 1240262221936811.36
    81999-11-2916：16：47.6123.0340.535.081118020.3 1770169.424820341 83246716.52
    91999-11-3007：52：55.8123.0340.534.091116813.0 1578157.526232162153486524.78
    101999-11-3013：58：17.7123.0340.535.2929180−16.3 2674−169.627015180 2260709.33
    111999-11-3014：06：55.1122.9840.554.982927210.519780161.7 4122133101667018.25
    121999-11-3014：09：36.6123.0240.554.391125811.8 1680147.4 7728350 33575720.22
    131999-11-3020：19：44.2123.0340.534.3929583−20.2 2870−172.529325 3725279704.69
    141999-12-0101：47：1.9123.0340.534.292908010.219880169.8 6722338 61547817.03
    151999-12-0104：33：0.5123.0340.534.491108030.5 1460168.4 428274 03076020.71
    161999-12-0112：45：30.8123.0340.554.392997010.620580159.729014 18 21176816.64
    171999-12-1305：49：30.5123.0840.534.111 1376023.3 3670147.9 7942346 2 85238.18
    181999-12-2719：27：15.4123.0240.534.081127010.6 1980159.7 9130 0 93566614.73
    192000-01-1207：43：55.4123.0540.535.5814680−40.725650−166.9 70253331313950140.51
    202000-01-1213：00：31.9123.0340.554.391128826.0 2064177.826913359 12976210.20
    下载: 导出CSV

    表  3   1999年岫岩MS5.4地震前小震震源机制

    Table  3   Focal mechanisms for small earthquakes before the MS5.4 Xiuyan earthquake in 1999

    序号发震日期地理坐标 MS深度
    /km    		节面Ⅰ节面ⅡP 轴T 轴B 轴
    　 年-月-日 时：分：秒东经/°北纬/°走向/°倾向倾角/°走向/°倾向倾角/°方位角/°俯角/°方位角/°俯角/°方位角/°俯角/°
    11999-01-1503：05：18.6122.6240.65 3.411128WS80 42NW62261263581111059
    21999-01-0605：43：10.9122.8240.65 2.910105NE76 17ES80 6218331 323372
    31999-01-2007：37：43.4122.5340.70 3.212101NE80 13ES70 5625150 825965
    41999-04-0916：31：48.0121.0042.00 3.015138WS52 32ES702711216746 1243
    51999-04-2321：13：38.4122.8240.70 2.9 7107WS70 20NW862461515212 1368
    61999-04-2715：14：14.4124.5741.12 2.813125NE85 40ES30 62411913630430
    71999-04-3021：16：13.3122.8539.68 2.7 5142NE78 60NW42 12402622215239
    81999-05-1112：29：25.9122.7840.65 3.414114WS80 27NW7024820341 8 9068
    91999-05-1501：45：46.6122.2239.38 3.110118WS55 35NW802623216215 5052
    101999-05-2220：40：45.1122.4740.73 3.4 9133WS82 46NW8027015180 2 9378
    111999-05-2415：58：35.1122.2041.70 2.528 86NW80 0NE68 41221331024667
    121999-05-2903：13：30.9121.0042.00 3.321119NE70 34ES72 7728350 325364
    131999-06-0301：46：50.2122.2541.67 3.4 8 75NW64164NE9029325 372516464
    141999-06-0616：51：47.4122.8240.65 2.814110NE70 26ES80 6722338 623166
    151999-08-2112：56：25.9121.7041.25 2.721142NE70 45NW70 428274 018361
    161999-08-3115：11：21.0122.2339.35 2.914153WS80 64NW8029014 18 210875
    171999-09-0605：05：30.3122.6740.70 3.1 6114NE58 41ES64 7942346 225247
    181999-09-2013：22：57.0122.6740.68 2.9 5136NE62 52ES76 9130 0 925658
    191999-10-1007：21：15.5122.9540.65 2.8 9109NE59 24ES80 70253331321956
    201999-10-1216：04：20.3122.1840.48 3.414132WS80 46NW8026913359 1 8774
    211999-11-0414：46：26.9122.6039.27 3.5 7118WS57 14ES7015239252 534850
    221999-11-0903：34：1.3121.5738.52 3.9 9114NE62 22ES902703713937 2225
    231999-11-0908：21：36.2123.0040.53 2.9 8118NE70 32ES80 7720345 724067
    241999-11-0917：44：44.6123.0040.53 3.6 8121NE70 34ES80 8020347 823867
    251999-11-1618：57：48.5123.0040.53 3.0 9110NE90 20ES60 61211612129160
    261999-11-1703：09：31.4123.0040.53 2.7 9112WS70 24NW8224920156 9 4370
    271999-11-1723：59：50.8123.0040.55 2.5 9102WS80 14NW8423911148 4 4677
    281999-11-1808：15：41.2123.0040.53 2.8 9111NE80 29ES50 62371672028048
    291999-11-2100：50：1.9123.0040.55 3.0 8122NE72 40ES70 8128172 226564
    301999-11-2500：59：19.2123.0040.55 3.0 8112NE80 32ES42 60391722528340
    311999-11-2521：17：50.8123.0040.53 3.5 8138NE70 55ES70 9630187 127861
    321999-11-2522：08：12.4123.0040.53 3.2 8147WS76 61NW8028518194 2 9273
    331999-11-2523：19：17.5123.0040.55 3.2 8110WS80 24NW7024620338 6 8467
    341999-11-2603：18：23.3123.0040.55 2.7 9123NE72 39ES70 8027171 126361
    351999-11-2623：36：20.7123.0040.53 3.1 9112NE60 38ES60 7740344 325346
    361999-11-2701：33：26.6123.0040.53 2.6 8148WS80 53ES6819221 99 835065
    371999-11-2702：29：43.7123.0040.55 2.8 8122WS74 38NW6826023350 2 9263
    381999-11-2715：48：2.6123.0040.55 3.2 8118WS80 36NW50250373542010747
    391999-11-2808：15：3.6122.9840.55 3.3 8108NE86 20ES62 60211561727961
    401999-11-2905：56：59.9123.0040.55 3.2 9121NE80 31ES47 61371672028146
    411999-11-2909：25：51.2123.0040.53 2.7 8126NE38 44ES86 99373462922660
    下载: 导出CSV
    • 参考文献(53)
  • 陈颙. 1978. 用震源机制一致性作为描述地震活动性的新参数[J]. 地球物理学报,21(2):142–159.

    Chen Y. 1978. Consistency of focal mechanism as a new parameter in describing seismic activity[J]. Acta Geophysica Sinica,21(2):142–159 (in Chinese).
    陈颙,刘杰,杨文. 2015. 前震序列的图像特征研究[J]. 中国地震,31(2):177–187. doi: 10.3969/j.issn.1001-4683.2015.02.001

    Chen Y,Liu J,Yang W. 2015. Pattern characteristics of foreshock sequences[J]. Earthquake Research in China,31(2):177–187 (in Chinese).
    程万正,阮祥,张永久. 2006. 川滇次级地块震源机制解类型与一致性参数[J]. 地震学报,28(6):561–573. doi: 10.3321/j.issn:0253-3782.2006.06.001

    Cheng W Z,Ruan X,Zhang Y J. 2006. Types of focal mechanism solutions and parameter consistency of the sub-blocks in Sichuan and Yunnan Provinces[J]. Acta Seismologica Sinica,28(6):561–573 (in Chinese).
    崔子健,李志雄,陈章立. 2015. 云南景谷MS6.6、云南沧源MS5.5地震谱振幅相关系数特征分析[J]. 地震研究,38(4):535–540. doi: 10.3969/j.issn.1000-0666.2015.04.003

    Cui Z J,Li Z X,Chen Z L. 2015. An analysis of correlation coefficient characteristic of spectral amplitude of Jinggu MS6.6 earthquake and Cangyuan MS5.5 earthquake[J]. Journal of Seismological Research,38(4):535–540 (in Chinese).
    刁桂苓,于利民,李钦祖. 1992. 震源机制解的系统聚类分析:以海城地震序列为例[J]. 中国地震,8(3):86–92.

    Diao G L,Yu L M,Li Q Z. 1992. Hierarchical clustering analysis of the focal mechanism solution:Taking the Haicheng earthquake sequences for example[J]. Earthquake Research in China,8(3):86–92 (in Chinese).
    刁桂苓,于利民,李钦祖. 1994. 强震前后震源区应力场变化一例[J]. 地震学报,16(1):64–69.

    Diao G L,Yu L M,Li Q Z. 1994. One case of variation in the focal stress prior to a strong earthquake[J]. Acta Seismologica Sinica,16(1):64–69 (in Chinese).
    刁桂苓,赵英萍,啜永清,王勤彩,高景春,曹肃朝,王焱,朱振兴. 2004. 大同晚期强余震前震源机制解的一致性特征[J]. 内陆地震,18(3):202–206. doi: 10.3969/j.issn.1001-8956.2004.03.002

    Diao G L,Zhao Y P,Chuo Y Q,Wang Q C,Gao J C,Cao S C,Wang Y,Zhu Z X. 2004. Coherence characteristics of focal mechanism solutions of later-period strong aftershocks[J]. Inland Earthquake,18(3):202–206 (in Chinese).
    顾浩鼎,陈运泰,高祥林,赵毅. 1976. 1975年2月4日辽宁省海城地震的震源机制[J]. 地球物理学报,19(4):270–285.

    Gu H D,Chen Y T,Gao X L,Zhao Y. 1976. Focal mechanism of Haicheng,Liaoning Province,earthquake of February 4,1975[J]. Acta Geophysica Sinica,19(4):270–285 (in Chinese).
    郭增建,秦保燕,徐文耀,汤泉. 1973. 震源孕育模式的初步讨论[J]. 地球物理学报,16(1):43–48.

    Guo Z J,Qin B Y,Xu W Y,Tang Q. 1973. Preliminary study on a model for the development of the focus of an earthquake[J]. Acta Geophysica Sinica,16(1):43–48 (in Chinese).
    郭增建,秦保燕,张远孚,黎在良. 1977. 从水平力和垂直力的相互作用讨论我国境内地震的孕育和发生[J]. 地球物理学报,20(3):242–250.

    Guo Z J,Qin B Y,Zhang Y F,Li Z L. 1977. A discussion on the mutual action of horizontal and vertical stresses in the development of earthquake sources in China[J]. Acta Geophysica Sinica,20(3):242–250 (in Chinese).
    韩晓明,荣代潞. 2015. 美国南加州地区1981—2011年MW≥6.0地震前发震应力场与构造应力场趋于一致现象研究[J]. 地震学报,37(6):948–958. doi: 10.11939/jass.2015.06.006

    Han X M,Rong D L. 2015. Consistency of seismogenic stress field of preshocks to the tectonic stress field before eight earthquakes (MW≥6.0) in southern California of United States from 1981 to 2011[J]. Acta Seismologica Sinica,37(6):948–958 (in Chinese).
    金严,赵毅,陈颙,鄢家全,卓钰如. 1976. 辽宁省海城地震前震震源错动方式的一个特点[J]. 地球物理学报,19(3):156–164.

    Jin Y,Zhao Y,Chen Y,Yan J Q,Zhuo Y R. 1976. A characteristic feature of the dislocation model of foreshocks of the Haicheng earthquake,Liaoning Province[J]. Acta Geophysica Sinica,19(3):156–164 (in Chinese).
    李金,周龙泉,龙海英,聂晓红,郭寅. 2015. 天山地震带(中国境内)震源机制一致性参数的时空特征[J]. 地震地质,37(3):792–803. doi: 10.3969/j.issn.0253-4967.2015.03.010

    Li J,Zhou L Q,Long H Y,Nie X H,Guo Y. 2015. Spatial-temporal characteristics of the focal mechanism consistency parameter in Tianshan (within Chinese territory) seismic zone[J]. Seismology and Geology,37(3):792–803 (in Chinese).
    李丽,宋美琴,刘素珍,扈桂让. 2015. 山西地区震源机制一致性参数时空特征分析[J]. 地震,35(2):43–50. doi: 10.3969/j.issn.1000-3274.2015.02.005

    Li L,Song M Q,Liu S Z,Hu G R. 2015. Spatial-temporal characteristics of the consistency parameter of focal mechanisms in Shanxi area[J]. Earthquake,35(2):43–50 (in Chinese).
    刘方斌,曲均浩,李国一,田兆阳. 2018a. 山东长岛震群震源机制解一致性参数时空演化特征[J]. 地震工程学报,40(5):1034–1041.

    Liu F B,Qu J H,Li G Y,Tian Z Y. 2018a. Spatial-temporal characteristics of the focal mechanism solutions consistency parameter of Changdao earthquake swarm[J]. China Earthquake Engineering Journal,40(5):1034–1041 (in Chinese).
    刘方斌,曲均浩,李亚军,范晓易,苗庆杰. 2018b. 山东乳山地震序列震源机制解一致性参数特征[J]. 地震地质,40(5):1086–1099.

    Liu F B,Qu J H,Li Y J,Fan X Y,Miao Q J. 2018b. Research on characteristics of the focal mechanism solutions consistency of Rushan earthquake sequence,Shandong Province[J]. Seismology and Geology,40(5):1086–1099 (in Chinese).
    荣代潞. 2014. 研究中强地震前中小地震震源机制变化的一种方法[J]. 地震工程学报,36(2):286–291. doi: 10.3969/j.issn.1000-0844.2014.02.0286

    Rong D L. 2014. Analyzing changes of focal mechanism of small-moderate earthquakes before a moderate-strong earthquake[J]. China Earthquake Engineering Journal,36(2):286–291 (in Chinese).
    孙丽娜,李皓,齐玉妍,温超,刁桂苓. 2017. 2004年12月26日印度尼西亚MW9.0大震前后震源机制一致性变化特征研究[J]. 中国地震,33(3):424–431. doi: 10.3969/j.issn.1001-4683.2017.03.008

    Sun L N,Li H,Qin Y Y,Wen C,Diao G L. 2017. A study on forecast method of focal mechanism consistency before the MW9.0 Indonesia earthquake on December 26,2004[J]. Earthquake Research in China,33(3):424–431 (in Chinese).
    万永革. 2008. 美国Landers地震和Hector Mine地震前震源机制与主震机制一致现象的研究[J]. 中国地震,24(3):216–225. doi: 10.3969/j.issn.1001-4683.2008.03.003

    Wan Y G. 2008. Study on consistency of focal mechanism of mainshock and that of preshocks in Landers and Hector Mine earthquake in United States[J]. Earthquake Research in China,24(3):216–225 (in Chinese).
    泽仁志玛,刁桂苓,李志雄,王晓山,冯向东. 2009. 千岛岛弧2006年MW8.3地震前震源机制解的一致性变化[J]. 地震学报,31(4):467–470. doi: 10.3321/j.issn:0253-3782.2009.04.014

    Zeren Z M,Diao G L,Li Z X,Wang X S,Feng X D. 2009. Variation of focal mechanism consistency before 2006 Kuril Islands arc MW8.3 earthquake[J]. Acta Seismologica Sinica,31(4):467–470 (in Chinese).
    泽仁志玛,刁桂苓,李志雄,王晓山,冯向东. 2010. 大震前显示的地震震源机制趋于一致的变化[J]. 地震,30(1):108–114. doi: 10.3969/j.issn.1000-3274.2010.01.012

    Zeren Z M,Diao G L,Li Z X,Wang X S,Feng X D. 2010. Consistent distribution of stress before strong earthquake from focal mechanism[J]. Earthquake,30(1):108–114 (in Chinese).
    张敏强. 2010. 教育与心理统计学[M]. 北京: 人民教育出版社: 253−256.

    Zhang M Q. 2002. Educational and Psychological Statistics[M]. Beijing: People’s Education Press: 253−256 （in Chinese）.
    张萍,蒋秀琴. 2001. 岫岩—海城MS5.4地震序列的震源机制解及应力场特征[J]. 地震地磁观测与研究,22(2):76–82. doi: 10.3969/j.issn.1003-3246.2001.02.015

    Zhang P,Jiang X Q. 2001. The focal mechanism solutions and the crust stress field characteristics in Xiuyan-Haicheng (MS5.4)earthquake sequence[J]. Seismological and Geomagnetic Observation and Research,22(2):76–82 (in Chinese).
    张萍,于龙伟,李涯,迮安民,刘天阁,吴野,杨红艳. 2003. 岫岩—海城5.4级地震前小震震源机制解与记录特征分析[J]. 地震地磁观测与研究,24(1):29–38. doi: 10.3969/j.issn.1003-3246.2003.01.005

    Zhang P,Yu L W,Li Y,Ze A M,Liu T G,Wu Y,Yang H Y. 2003. An analysis on the focal mechanism solutions of small earthquakes and the record characteristics before the Xiuyan-Haicheng earthquake of 5.4[J]. Seismological and Geomagnetic Observation and Research,24(1):29–38 (in Chinese).
    张致伟,周龙泉,龙锋,阮祥. 2015. 汶川8.0和芦山7.0级地震序列应力场时空特征[J]. 地震地质,37(3):804–817. doi: 10.3969/j.issn.0253-4967.2015.03.011

    Zhang Z W,Zhou L Q,Long F,Ruan X. 2015. Spatial and temporal characteristic of stress field for Wenchuan MS8.0 and Lushan MS7.0 earthquake sequence[J]. Seismology and Geology,37(3):804–817 (in Chinese).

    Lund B,Bӧðvarsson R. 2002. Correlation of microearthquake body-wave spectral amplitudes[J]. Bull Seismol Soc Am,92(6):2419–2433. doi: 10.1785/0119990156

    Michael A J. 1991. Spatial variations in stress within the 1987 Whittier Narrows,California,aftershock sequence:New techniques and results[J]. J Geophys Res,96(B4):6303–6319. doi: 10.1029/91JB00195

    Wiemer S,Gerstenberger M,Hauksson E. 2002. Properties of the aftershock sequence of the 1999 MW7.1 Hector Mine earthquake:Implications for aftershock hazard[J]. Bull Seismol Soc Am,92(4):1227–1240. doi: 10.1785/0120000914
    • 相关文章
    • 施引文献(5)

  • 期刊类型引用(5)

    1. 田唯熙，张永仙，张盛峰，张小涛. 区域选取对图像信息法可预测性的影响. 地震学报. 2024(02): 208-225 . 本站查看
    2. 宋程，张永仙，夏彩韵，毕金孟，张小涛，吴永加，徐小远. 基于PI方法的华北2019年以来3次M_S≥5.0地震回溯性预测研究. 地震. 2024(02): 120-134 . 百度学术
    3. 宋程，张永仙，周少辉，毕金孟，徐小远. 2021年玛多M_S7.4地震的PI热点特征回溯性预测研究. 地震研究. 2023(02): 226-236 . 百度学术
    4. 田唯熙，张永仙. 基于图像信息方法的南北地震带地震预测研究. 地震. 2023(03): 159-177 . 百度学术
    5. 余娜，张晓清，袁伏全，杨晓霞. 基于PI法的门源M_S6.4地震前震中附近地震热点图像异常变化研究. 浙江大学学报(理学版). 2020(06): 724-729+742 . 百度学术

    其他类型引用(0)

    • 资源附件(0)
图(7)  /  表(3)
计量
  • 文章访问数:  3578
  • HTML全文浏览量:  692
  • PDF下载量:  91
  • 被引次数: 5
出版历程
  • 收稿日期:  2019-04-11
  • 修回日期:  2019-05-12
  • 网络出版日期:  2019-12-16
  • 发布日期:  2019-10-31

目录

摘要
HTML全文

  • 引言

  • 引言

  • 1.   最小二乘配置法原理

  • 1.   最小二乘配置法原理

  • 1.1   方法原理

  • 1.1   方法原理

  • 1.2   协方差函数

  • 1.2   协方差函数

  • 2.   其它插值法原理

  • 2.   其它插值法原理

  • 3.   模型试验

  • 3.   模型试验

  • 3.1   噪声对最小乘配置方法的影响

  • 3.1   噪声对最小乘配置方法的影响

  • 3.2   协方差函数的拟合

  • 4.   华北地区实际重力测点试验

  • 4.   华北地区实际重力测点试验

  • 5.   讨论与结论

  • 5.   讨论与结论

参考文献(53)
相关文章
施引文献(5)
资源附件(0)

/

返回文章
返回

本系统由 北京仁和汇智信息技术有限公司 开发  技术支持：info@rhhz.net 百度统计