第41卷第1期 2021年1月
地 震
E A R T H Q U A K E
Vol.41, No. 1
J a n.，2021
d o i：10. 12196/j. issn. 1000-3274. 2021. 01. 002
时变重力场球面模型反演算法和模拟实验
张贝1’2，陈石U2，李红蕾U2，杨锦玲韩建成U2，卢红艳K
(1.中国地震局地球物理研究所.北京100081; 2.北京白家疃地球科学国家野外观测
研究站，北京100095; 3.福建省地震局，福建福州3.50003)
摘要：时变重力场是研究地球内部介质物性变化的重要手段。本文提出丫 一种适用于地面流动
重力测量获得的时变重力信号的场源反演方法，该方法采用球坐标系下的六面体单元来模拟
场源介质，适合大尺度地璲流动重力测t t数据的等效源模墦构建。通过引人重力时变信号的•
阶光滑先验条件，压制f时变重力信号中的短周期高频分黾，H了用于提取~地震孕育相关的长
周期信号。通过理论和模型实验证明了本文算法的可靠性和稳定性，并使用南北地震带南段
2014—2017年的流动重力实测数据进行f反演解释.获得了地壳内部等效场源的视密度时变
信号，变化_t t级在正常地壳密度的±〇_ 7%。之间，其空间形态受川滇菱形块体边界控制。研究
成果可用于时变重力场模型解释和深部场源特征提取.可为地震重力前兆信号分析和相关研
究提供完备的方法保障。
关键词：时变重力场；位场反演；流动重力；等效源；密度变化
中图分类号：P315.7 文献标识码：A文章编号：1000-3274(2021)0卜0013-12
引言
时变重力场特别是高精度的微伽级时变重力信号，是研究地球内部物质迁移和变形过程的重要科学依据。21世纪以来，以GRACE E星为代表的空间对地观测手段不断发展，使得人类可以获取全球尺度的重力变化图像1.这些时变重力信号在全球陆地水储量变化’'、同震物质迁移1和两极冰川融化速率估计等方面取得了一系列重要的成果。除 此之外，通过陆基固定点的定期重复观测，也可以获得时变重力信号，相比卫星重力手段，陆基流动重力观测距离地壳内部场源更近、可观测的重力信号量级更大，但是由于测量受 到地形和道路限制，测点分布不均匀且容易受到局部近场源的环境变化干扰：'8:。
我们知道，时变重力场反映的是场源的属性变化。长期以来，科学家们一直期望通过陆基重力观测获取地壳深部孕震区介质变化的信息一1。但由于陆基流动重力测点分布不 *
*收稿日期：2020-07-16；修改回日期：2020-09-03
基金项目：科技部重点研发专项项H(2018YFC1504506; 2018YFC0603502)，国家自然科学基金面上项目(41774090)，国家自然科学基金地震联合基金项目（U1939205)，中国地震局地球物理研究所基本科
研业务费专项资助（DQJB20X09〉，中国地震科学试验场专项（2019CSES0105)联合资助 作者简介：张贝（1985-)，男，河南商丘人，助理研究员，主要从事计算地球动力学与数值模拟。
通讯作者：陈石，研究员。E-mail: ******************
14地 湓41卷
均、测量容易受到高频噪声污染等问题，如何从空间非均匀分布的有限重力测点数据中提 取地壳深部场源信号，一直是国内外研究的热点问题之一 除此之外，陆基流动重力测点还可能由于环境变化等原因产生位置变化，从而引起有效的重力时变信号采样点发生 改变。
本文以“以场求源、场源结合”的思想为指导，通过构建一个等效源模型11，以可控分 辨率的离散网格单元模拟场源介质，采用球坐标系下的六面体单元:18构建场源格林函数；通过引人时空光滑等先验条件[19:，给出了时变重力场正则化反演的核心算法。该方法对不 N期次的陆基时变重力数据测点变化+敏感，可以减小由于测点变化引起的时变重力信号 损失问题，同时，由于采用球坐标系下的六面体单元，不需要进行测点的坐标投影变换，特别适合于测点间距在十几至几十千米尺度的地震流动重力数据场源建模问题。
本文分以下几个部分开展研究，首先在第一部分给出了本文提出的等效源反演方法原 理和正则化约束定义及反演求解方法；在第二部分，通过设计的理论模型，测试了由于测 点变化、噪声干扰等情况下的算法稳定性和有效性；第三部分针对实际数据进行了测试，反演得到的重力时空场源变化特征与川滇菱形构造具有一致性；最后，详细讨论了该方法 的优缺点及可能在实践应用中遇到的问题，并对潜在的应用领域和方向进行了进一步分析 和探讨。
1方法原理
经典的重力位场反演方法1：同样适合时变重力场反演问题.两者之间的区別在于时变重力反演需要对多个离散时间节点上的观测重力信号进行反演和参数优化。由于重力位 场反演通常在没有合适先验约束的条件下，垂向分辨能力不佳。因此.本文在不引人过多假设条件的基础上，采用“等效源”的概念来描述时变重力场的场源特征，在诸如对时变卫 星重力信号的分析解释中，也有类似的用反演“等效水厚度”来描述时变重力场源特征的方 法1。在本节中，我们分三部分从基本方程、正则化约束和贝叶斯参数优化三个方面来展 开论述。
1.1基本方程
在时变重力场的场源观测方程可以表示为如式（1)的线性方程组形式。
/(x,y,z,t)=G m(1)其中，G为核函数矩阵，m为场源模型，/为在离散观测系统上的四维函数。如果分别用 M.N和K来表示场源模型数.观测点数和观测期/次数.则式（1)中的场源模型m可以写 为以下向量形式，m=[wM1，，…，];；观测场/的向M形式为/=[/.\_丨，/\'」，…，/、<]T。核函数G为全乂维矩阵，其中非零元素个数为I：吣X N t。
k-\ k - I it 1
通常，当观测系统和场源模型设计好后，可以确定核函数矩阵G，从场源模型m汁算 观测系统上的理论重力场/为正演问题。而从观测场/估计场源模型/»为反演问题。反演 问题是典型的不适定问题.一般需正则化处理后才可以求解。反演问题的求解可以转换为 式（2)所示的H标函数最小二乘问题：
<p=W|| Gm — /(x,y,z,l)|| L+<PAm) + <P, (m)(2)其中，等式右端第一项为观测部分.第二和三项为正则化约束部分，与模型先验假设有关。
1期
张贝等：时变重力场球面模型反演算法和模拟实验
1.2 正则化约束
由于重力位场反演问题在数学上属于典型的不适定问题，通常需要正则化约束才能定解。对于时变重力场的反演，正则化约束由空间项和时间项两部分构成。
(1)空间约束
在式(2)中的右端第二项少可以表示为如下形式：
其中，W,、和W；是与模型在空间各方向光滑程度相关的权系数（超参数）。具体物理含义是期望待求的场源模型在X，X Y和Y各个水平方向上都具有一阶光滑特征，具体光 滑程度与、W：和这三个待优化的未知超参数有关。
(2)时间约束
在式（2)中的右端第三项可以表示为如下形式：
<Pr(m) =J w,*(4)其中，是与模型在时间上光滑程度相关的权系数（超参数）。物理含义是期望待求的场源模型在时间上具有--阶光滑特征，光滑程度由超参数来控制。
1.3超参数优化
在反演计算中，式（2)〜式（4)中的W，至W，共5个超参数需要确定后才可以求解。为了合理确定这些超参数•本文引入贝叶斯方法中的A B I C准则来实现超参数的优化计算。具体实现是通过A B I C最小化实现。
ABIC =—21og(L) +2N(5)式（5)中，L是模型的似然函数.N是超参数个数。A B I C最小化问题是非线性问题，可以 采用单纯形方法或者牛顿法等具有多参数非线性优化能力的方法来求解。
2模型实验
本文研制的时变重力的场源反演方法，主要是为了解决阽地时变重力测M的场源模型 确定问题。为了验证方法的有效性，本文以南北地震带南段的实际重力测点分布为依据，通过设计两组不同特征的检测板模型，分别从无噪声和有噪声两个角度，i、J•论实际地表非 均匀分布重力测点的场源模型反演问题。
图1是南北地震带南段的实际流动重力测点分布，测点数为412个，如图中白圆点所示。从测点位置可以发现，测点间距不均，且几何形态与地形和交通状况密切相关。一 般认为，测点的非均匀分布程度直接影响场源的监测能力，对测点分布空区或显著稀疏分布位置的场源，其监测能力不会太高。在图1的重力测网屮.即使在最密集的西昌-攀枝花昭通地区，也存在显著的环形监控盲区。
本文基于球面六面体（Tesseroid)模型.统一采用以0.5°X0. 5°的模型单元尺度来实现场源模型的离散化，正负相间的场源体密度取±l X10 :i g/cm5，场源体埋深10 km，每个 场源模型的等效厚度为1km。表1给出了两组不同分辨率的检测板参数，以及地表重力测 网可观测的理论模型重力异常。下面我们分别从以下两个方面开展综合模型实验。
1641卷
地 震97° E 99° 101° 103° 105° 107°
A 程/m
-5000
-4000
-3000
|2〇〇〇
-1000
图1南北地震带南段流动重力观测网络
图中白圆点为陆地重力观测点位置.红实线为活动断裂构造，白实线为国界
Fig. 1 T h e  gravity obvservation netw ork in the southern section of the South-N orth seismic zone
th e  w h ite  d o ts sh o w  th e  location o f th e  g ra v ity  s ta tio n , th e  red solid lines show  th e
activ e fa u lts , and th e w h ite solid line is th e  b o rd e r
表1
检测板模型测网可观测重力异常 Table 1 The observable gravity anomalies for the checkerboard model
检测分辨率最大异常 /10_K m #s -最小异常
/10 8 m • s -
中数异常/10 8 m • s -平均异常/10-8 m • s -2模型11°X 1°32. 10-32. 06
-3. 77一 1. 61模型2()• 5〇X 0. 5。23.00
-23. 270. 200. 492. 1检测板实验
本文的检测板实验是通过设计一组正负相间的同尺度场源体，通过正演方法得到实际 测点位置的理论重力异常，然后利用反演方法获得场源模型参数，通过与已知场源参数对 比，来评价由于实际测网分布不均匀而导致的空间非均匀场源分辨能力的差异性问题。
图2是r x r 分辨率的检测板模型1在地表可观测的理论重力异常与地表实际测点的 空间分布情况.整个可观测的重力异常量级在±30X l (r s  m  • 变化，由于地表实际重力 测点分布不均匀，能通过现有陆地重力
测网观测到的异常如图3所示。图3中除了中甸以 西、百以南等没有重力测点分布的测网以外区域，在其测网内部也有多个类似于泸州至 文山之间的大部分空区存在，可见地表重力测点多以“环状”展布，而每个测环中心区域的 监控盲区也将十分不利于场源结构的反演。对比图2和图3的结果可以发现，如果仅通过 图3的地表观测重力异常场信息，很难与图2所示的理论重力异常特征联系起来，因此.
1期 张贝等：时变重力场球面模型反演算法和模拟实验 17由于地表测点的非均匀分布和观测局限导致的稀疏采样问题，如何影响反演？而本文设计 的正则化手段是否有利于恢复场源参数，有效压制噪声影响？下面我们进一步结合模型测 试结果进行讨论。
98° E 100° 丨 02° 104° 106° 1081
28。
N 26°24。22。图2检测板模型地表理论重力异常
Fig. 2 T h e  surface theoretical gravity anomaly of checkerboard model
98°E  100。 102。 104。 106。 108。
28。- N 26。-24。-22。-阁3地表重力测网可观测的重力异常
Fig. 3 T h e  surface observable gravity anomaly at the surface gravity ne tw o rk
^ mm
泸州
* P 义
攀枝花••贵
• • ••临沧• •• I f f ,
• _玉溪
百 河池
文山
f ：力异常/10*
-30 20 -10 0 10 20 30
■■ r s s w
鼸鏟罈鼉爾糴_8驄
«鼸激鑤麝隳癲___ _翁篇邐國钃爨_a r
»篇霤_垂鷂窸騸麴_
▲■钃■
_鶸繡_钃驂霾_正则化坐标
S 力异常/108m.s 2
, , , 'm m 10 2030
图4是基于图3所示的地表实际测网采样到的重力异常反演后的场源参数结果。图

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。