Numerical and theoretical analyses of seismic wave response in non-persistent jointed rock mass
-
摘要:
为研究地震产生的应力波在断续节理岩体中的传播规律和应力分布趋势，首先，采用数值模拟方法分析应力波通过贯通节理的传播规律，并与已有理论研究结果进行对比，验证数值分析的准确性和适用性；然后，对应力波在断续节理岩体中的传播进行数值模拟，分析透射系数在水平方向的分布趋势以及不同节理连续性对波传播的影响，并结合波的衍射原理，给出定性的理论解释。结果表明：应力波通过断续节理时，节理的透射作用会使应力波振幅减小，引起波的衰减，岩桥的衍射作用则会使波阵面由平面变为曲面，波的传播方向发生改变，从而导致应力波振幅在水平方向的分布发生变化；应力波通过断续节理的透射系数与岩桥尺寸Lr和衍射角μ相关，当衍射角比较小时，透射系数主要受岩桥尺寸Lr的影响，当衍射角较大时，岩桥尺寸Lr和衍射角μ共同影响应力波在岩体中的传播。
Abstract:In order to reveal the regulation of seismic stress wave propagation and the distribution of wave amplitude in jointed rock masses, the wave propagation across persistent joints is simulated firstly, and the results are verified to be consistent with analytical results. Then, the horizontal distribution of transmission coefficient and the effect of joint persistency on wave propagation are researched by simulating seismic stress wave propagation across rock masses with a single non-persistent joint. Moreover, the theoretical analysis is proposed to explain the simulation results. The result shows that the joint segments induce wave attenuation, and rock bridges cause the change of the wave propagation direction. On the other hand, the transmission coefficient of stress wave through non-persistent joints is related to the size of rock bridge （Lr） and the diffraction angle （μ）. When the diffraction angle is small, the transmission coefficient is mainly affected by the size of rock bridge; otherwise, the size of rock bridge and the diffraction angle both affect the propagation of stress wave across non-persistent joints.
-
Keywords:
- seismic stress wave /
- non-persistent joints /
- numerical simulation /
- wave diffraction
-
引言
地球重力场是研究地球结构与性质的重要物理参数，其变化与地球介质密度和地球动力学过程（固体地球潮汐、内部热流、固体和液体之间质量的交换、表面负荷和地震构造运动等）关系密切（孙和平，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 ChinaThe 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 methods3.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 ChinaThe 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 methods3.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 method4. 华北地区实际重力测点试验
本文以现有重力测点进行模拟试验，基于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 method5. 讨论与结论
本文采用了最小二乘配置法研究华北地区的局部重力场建模问题，比较了最小二乘配置法与反距离加权法、最小曲率法恢复局部重力场的效果，讨论了噪声对最小二乘配置结果的影响以及协方差函数拟合问题，主要研究结论如下：
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 应力波通过贯通节理的数值计算模型
Figure 1. The numerical model of stress wave propagation across persistent joints
![]()
图 2 在x-t平面上特征线相交于j和n（引自Cai，Zhao，2000）
图中j和n取整数，α表示波速，$\varDelta$表示一定时段
Figure 2. Characteristics conjunction points at integer values of j and n in the x-t plane （after Cai，Zhao，2000）
j and n are integer，α is velocity，$\varDelta $ is time interval
![]()
图 3 波通过单条节理的透射系数T1与节理法向刚度Kn的关系
Figure 3. Transmission coefficient T1 versus normal stiffness of joints Kn for normally incident wave transmission across a single joint
![]()
图 4 波通过节理组的透射系数Tn与无量纲节理间距ξ的关系
Figure 4. Transmission coefficient Tn versus nondimensional joint space ξ for normally incident wave transmission across joint set
![]()
图 6 波在断续节理岩体中传播的数值模型
Figure 6. The numerical model of wave propagationacross non-persistent jointed rock mass
![]()
图 7 针对不同节理连续性k情形时透射系数TP的水平分布规律
Figure 7. The horizontal distribution of transmission coefficient TP for different persistency k
![]()
图 8 不同测点处P波的透射系数TP与节理连续性k的关系曲线
Figure 8. The transmission coefficient TP of P wave versus joint persistency k at different measuring points
-
Cai J G,Zhao J. 2000. Effects of multiple parallel fractures on apparent attenuation of stress waves in rock masses[J]. Int J Rock Mech Min Sci,37(4):661–682. doi: 10.1016/S1365-1609(00)00013-7
Francon M. 1986. Recent developments in monochromatic birefringent filters[G]//Trends in Quantum Electronics. Heidelberg: Springer: 93−98.
Jiao Y Y,Fan S C,Zhao J. 2005. Numerical investigation of joint effect on shock wave propagation in jointed rock masses[J]. J Test Evaluat,33(3):197–203.
Pyrak-Nolte L J,Myer L R,Cook N G W. 1990. Transmission of seismic waves across single natural fractures[J]. J Geophys Res,95(B6):8617–8638.
Schoenberg M. 1980. Elastic wave behavior across linear slip interfaces[J]. J Acoust Soc Am,68(5):1516–1521. doi: 10.1121/1.385077
Wasantha P L P,Ranjith P G,Xu T,Zhao J,Yan Y L. 2014. A new parameter to describe the persistency of non-persistent joints[J]. Eng Geol,181:71–77. doi: 10.1016/j.enggeo.2014.08.003
Zhao J,Zhao X B,Cai J G. 2006. A further study of P-wave attenuation across parallel fractures with linear deformational behaviour[J]. Int J Rock Mech Min Sci,43(5):776–788. doi: 10.1016/j.ijrmms.2005.12.007
Zhao X B,Zhao J,Cai J G,Hefny A M. 2008. UDEC modelling on wave propagation across fractured rock masses[J]. Computers Geotechn,35(1):97–104. doi: 10.1016/j.compgeo.2007.01.001
Zhu J B,Deng X F,Zhao X B,Zhao J. 2013. A numerical study on wave transmission across multiple intersecting joint sets in rock masses with UDEC[J]. Rock Mech Rock Eng,46(6):1429–1442. doi: 10.1007/s00603-012-0352-9
-
期刊类型引用(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)
下载:
计量
- 文章访问数: 1642
- HTML全文浏览量: 456
- PDF下载量: 31
- 被引次数: 5
