一种地震面波走时和重力异常联合反演方法与系统

1.本发明涉及地震数据处理技术领域，特别是涉及一种地震面波走时和重力异常联合反演方法与系统。


背景技术：

2.面波层析成像从数据来源上，可以分为地震面波和背景噪声面波层析成像，这主要是依据面波是从一个真实的地震震源处激发，还是从一个虚拟的台站处激发来考虑的。对于后一种情况，其中的“面波”更准确来说是一种在特定条件下面波占据主要成分的叠加信号。在射线类方法中，又可以按照是否将反演群相速度图作为中间步骤分为两步法和一步法。同时，由于面波本质上是一个(粘)弹性介质的本征值问题，它在低频下也常常通过地球自由震荡的方式表现出来，因而可以用来研究地球的大尺度结构。尽管分类比较繁多，但这也从侧面说明了目前已经对面波开展了多手段的研究工作。
3.重力异常数据反映了地球真实的密度分布相对于参考椭球体的扰动。考虑到重力大小随距离平方衰减的性质，它通常被认为反映了地壳浅部和莫霍面附近的密度异常。因此，重力异常数据通常用来研究莫霍面起伏或者地壳内的密度变化。目前已经有很多文章利用重力异常来反演莫霍面结构，反演浅地表和地壳内密度结构等。重力反演研究的主要思路是首先利用特定的场分离方法从观测重力异常中提取需要研究的异常体对应的重力异常信号，之后再利用反演方法获取该异常体的密度结构。
4.由于地震面波具有频散特征，各个周期的面波传播速度不同。由于这种频散效应是由地下介质唯一决定的，因此利用频散信息可以约束地下速度结构信息。重力异常数据反映了地球真实的密度分布相对于参考椭球体的扰动，考虑到重力大小随距离平方衰减的性质，其反映了浅层地质结构的密度异常。由于每一种类型的数据仅仅对某些区域或者某些尺度的物理量敏感，不能完全地刻画地下介质的精细结构，而联合反演可以通过取长补短的方式弥补这一缺陷。通过利用岩石物理实验中的速度
‑
密度经验关系，联合地震面波频散和重力异常数据可以改善单一数据速度结构成像的分辨率，提高反演结果的可靠性。


技术实现要素：

5.本发明的目的是提供一种地震面波走时和重力异常联合反演方法与系统，以解决单一数据集对地下结构约束能力较低的问题。
6.为实现上述目的，本发明提供了如下方案：
7.一种地震面波走时和重力异常联合反演方法，包括：
8.获取待反演的地震面波频散数据和布格重力异常数据；
9.对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
10.计算敏感核；
11.根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
12.基于岩石的速度
‑
密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
13.根据修正后的线性化反演模型建立联合反演目标函数；
14.对所述反演目标函数进行求解得到最优解；
15.利用所述最优解对待反演的地震面波进行反演成像。
16.优选的，所述对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据，包括：
17.采用模型：
18.f(m)＝f(α,β,ρ)＝d
19.对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；其中，m表示真实模型，d表示理论观测数据，f表示表示真实模型到观测数据的映射，α表示p波速度，β表示s波速度，ρ表示介质密度。
20.优选的，所述线性化反演模型，包括：
[0021][0022]
其中，ω表示面波频率，表示在面波频率为ω下的第i个观测面波走时，t
i
(ω)表示第i个理论面波走时，δ0t
i
(ω)表示第i个观测面波走时和理论面波走时之间的残差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δc
j
(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，k
ij
(ω)是相应的敏感核。
[0023]
优选的，所述基于岩石的速度
‑
密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型，包括：
[0024]
获取岩石的速度
‑
密度关系；
[0025]
将根据所述速度
‑
密度关系将s波速度结构转换成密度和p波速度公式；所述密度和p波速度公式为：
[0026]
α＝f
α
(β),ρ＝f
ρ
(β)
[0027]
其中，α表示p波速度，β表示s波速度，ρ表示密度，f
α
表示s波速度模型到p波速度的映射，f
ρ
表示s波速度到密度的映射；
[0028]
根据密度和p波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型。
[0029]
优选的，所述修正后的线性化反演模型为：
[0030][0031]
d/dβ是对s波速度的全偏导数，t
i
(ω)表示第i个理论面波走时，δt
i
(ω)表示修正后的第i个观测面波走时和理论面波走时之间的残差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的s波速度差值。
[0032]
优选的，所述根据修正后的线性化反演模型建立联合反演目标函数，包括：
[0033]
采用公式：
[0034][0035]
建立联合反演目标函数；其中，g
s
表示面波数据理论算子，g
g
表示重力数据理论算子，σ是二阶吉洪诺夫正则化系数，s是二阶导数矩阵，γ是阻尼系数，i是单位矩阵，w
s
是面波走时权重矩阵，w
g
是重力数据的权重矩阵，是重力数据的权重矩阵，其中p是一个0
‑
1之间的常数，n
s
是面波的数据个数，n
g
是重力的数据个数，表示面波数据中理论模型与真实模型之间的标准差，表示重力数据中理论模型与真实模型之间的标准差。
[0036]
本发明还提供了一种地震面波走时和重力异常联合反演系统，包括：
[0037]
数据集获取模块，用于获取待反演的地震面波频散数据和布格重力异常数据；
[0038]
正演模块，用于对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
[0039]
计算模块，用于计算敏感核；
[0040]
反演模块，用于根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
[0041]
线性化反演模型修正模块，用于基于岩石的速度
‑
密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
[0042]
联合反演目标函数构建模块，用于根据修正后的线性化反演模型建立联合反演目标函数；
[0043]
求解模块，用于对所述反演目标函数进行求解得到最优解；
[0044]
反演成像模块，用于利用所述最优解对待反演的地震面波进行反演成像。
[0045]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0046]
本发明提供了一种地震面波走时和重力异常联合反演方法与系统，本发明提供的一种地震面波走时和重力异常联合反演方法包括：首先对待反演的地震面波频散数据和布格重力异常数据进行正演得到理论观测数据并计算敏感核；然后根据理论观测数据和敏感核进行反演得到线性化反演模型；基于岩石的速度
‑
密度关系对线性化反演模型进行修正得到修正后的线性化反演模型；最后根据修正后的线性化反演模型建立联合反演目标函数；对反演目标函数进行求解得到最优解；利用最优解对待反演的地震面波进行反演成像。本发明通过利用地震面波频散数据和布格重力异常数据对地震面波进行联合反演成像，不仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性。
附图说明
[0047]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施
例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0048]
图1是依照本发明实施例的速度模型的3次均匀b样条插值结果；
[0049]
图2是依照本发明实施例的地震面波走时和重力异常联合反演方法流程图；
[0050]
图3是依照本发明实施例的检测板测试分析图，(a
‑
d)为真实模型，(e
‑
h)为面波单独反演结果，(i
‑
l)为联合反演结果；
[0051]
图4是依照本发明实施例的异常体测试分析图，(a
‑
b)为真实模型，(c
‑
d)为面波单独反演结果，(e
‑
f)为联合反演结果；
[0052]
图5是依照本发明实施例的4条剖面上的速度反演结果，第一列为真实模型，第二列为面波单独反演结果，第三列为联合反演结果。
具体实施方式
[0053]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0054]
本发明的目的是提供一种地震面波走时和重力异常联合反演方法与系统，以解决单一数据集对地下结构约束能力较低的问题。
[0055]
为实现上述目的，本发明提供了如下方案：
[0056]
一种地震面波走时和重力异常联合反演方法，包括：
[0057]
步骤1:获取待反演的地震面波频散数据和布格重力异常数据；
[0058]
步骤2:对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
[0059]
其中，步骤2具体包括：
[0060]
采用了高斯
‑
勒让德数值积分方法计算重力场的径向分量，然后基于径向分量采用模型：
[0061]
f(m)＝f(α,β,ρ)＝d
[0062]
对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；其中，m表示真实模型，d表示理论观测数据，f表示表示真实模型到观测数据的映射，α表示p波速度，β表示s波速度，ρ表示介质密度。
[0063]
步骤3:计算敏感核；
[0064]
步骤4:根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；在本发明中，所述线性化反演模型，包括：
[0065][0066]
其中，ω表示面波频率，表示在面波频率为ω下的第i个观测面波走时，t
i
(ω)表示第i个理论面波走时，δ0t
i
(ω)表示第i个观测面波走时和理论面波走时之间的残
差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δc
j
(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，k
ij
(ω)是相应的敏感核。
[0067]
步骤5:基于岩石的速度
‑
密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
[0068]
其中，步骤5具体包括：
[0069]
获取岩石的速度
‑
密度关系；
[0070]
将根据所述速度
‑
密度关系将s波速度结构转换成密度和p波速度公式；所述密度和p波速度公式为：
[0071]
α＝f
α
(β),ρ＝f
ρ
(β)
[0072]
其中，α表示p波速度，β表示s波速度，ρ表示密度，f
α
表示s波速度模型到p波速度的映射，f
ρ
表示s波速度到密度的映射；
[0073]
根据密度和p波速度公式对所述线性化反演模型进行修正得到修正后的线性化反演模型；所述修正后的线性化反演模型为：
[0074][0075]
d/dβ是对s波速度的全偏导数，t
i
(ω)表示第i个理论面波走时，δt
i
(ω)表示修正后的第i个观测面波走时和理论面波走时之间的残差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的s波速度差值。
[0076]
步骤6:根据修正后的线性化反演模型建立联合反演目标函数；
[0077]
其中，步骤6具体包括：
[0078]
采用公式：
[0079][0080]
建立联合反演目标函数；其中，σ是二阶吉洪诺夫正则化系数，s是二阶导数矩阵，γ是阻尼系数，i是单位矩阵，w
s
是面波走时权重矩阵，w
g
是重力数据的权重矩阵，是重力数据的权重矩阵，其中p是一个0
‑
1之间的常数，n
s
是面波的数据个数，n
g
是重力的数据个数，表示面波数据中理论模型与真实模型之间的标准差，表示重力数据中理论模型与真实模型之间的标准差。
[0081]
步骤7:对所述反演目标函数进行求解得到最优解；
[0082]
步骤8:利用所述最优解对待反演的地震面波进行反演成像。
[0083]
下面结合具体的实施例对本发明提供的一种地震面波走时和重力异常联合反演方法进行进一步的说明：
[0084]
每一种类型的数据仅仅对某些区域或者某些尺度的物理量敏感，它们不能完全地
刻画地下介质的精细结构。而本发明中的联合反演可以通过取长补短的方式弥补这一缺陷。因此，本发明需要用两套数据来构建一个联合反演问题，包括：模型参数化：通过一系列参数作为模型空间的“基矢量”，使得模型空间内每一点都能够用这些参数表示(即给定插值方式),同时给定这些参数的初始值。一般来说，表征一个模型通常包括模型类的参数化方法和节点类的参数化方法。他们的区别是在前者中这些参数表征的是模型基函数的系数，而在后者中表示的是网格节点上的模型物理量的具体数值。当然在某些情况下两种表述是完全一致的。正演计算：通过这些参数计算出该模型对应的理论观测数据，在本发明的方法中，在正演建模时考虑了面波射线的几何形态以及地球的曲率对面波和重力数据的影响。
[0085]
请一并参阅图1
‑
5，为使本发明的目的和技术方案更加的便于理解，下面将参照附图更为详细的说明本发明的实施例。
[0086]
本发明采用了高斯
‑
勒让德数值积分方法计算重力场的径向分量：
[0087][0088]
其中，(r,θ,φ)表示球坐标系的三个方向，对应径向分量、水平分量和纵向分量，i,j,k表示每个方向对应的积分阶数，g
r
是在接收点(r0,θ0,φ0)处的重力场的径向分量，ω
ijk
表示权重系数，n
i,j,k
表示使用的高斯
‑
勒让德积分的最高阶数。
[0089]
由于本发明在面波走时正演中使用了快速多极算法方法，该方法要求利用球面坐标系的均匀网格节点来进行有限差分型的计算，因此需要对这种均匀分布的节点网格选择一个比较好的插值方法。对于均匀分布的网格，一个比较好的参数化方法是利用3次均匀b样条来对速度模型插值。图1给出了经三次均匀b样条插值后的一个速度模型的分布特征，可以看到有间断的模型在参数化后具有了很好的平滑特征。
[0090]
同时，为了更准确地计算面波走时，在模型参数化过程中有一些小技巧。第一，由于实际问题使用的初始模型的尺度可能不足以对速度模型进行足够精确的采样，这通常会导致走时场计算不准确。因此，在做正演计算时，首先需要构建一套计算网格。这套网格是用模型网格在更密的网格节点上插值得到的。
[0091]
另一个问题是震源附近的速度结构对走时的影响非常大，而且震源位置并不一定就位于网格节点上。这样会造成初始几个临近节点中的走时计算不准确，从而影响最终的结果。因此先要利用模型网格将震源附近的网格局部加密，然后将这个局部网格上的计算出来的走时插值到计算网格上。最终，为了执行一次正演计算，本发明在模型网格节点之外，还需要另外构建两套用于计算的网格。
[0092]
对于面波一维模型的正演，本发明仍然采用节点式的模型，但是在深度方向让网格的间距可以自由变化以适应速度模型的跳变。在局部化假设中，为了利用汤姆森
‑
哈斯克尔矩阵方法计算频散曲线，需要将网格节点模型转换为层状模型：这里本发明假设垂直方向的速度线性变化，在两个网格节点之间插入几个均匀层，层间的弹性参数用其中点位置的值来代替。在这样的假设下，就可以利用三维节点模型，首先计算出频散图，之后用快速多极算法方法计算各个周期的面波走时。
[0093]
同时，为了计算的方便，本发明让重力正演计算的网格与面波正演的网格位于同一套网格节点上。最终就可以利用这套三维网格节点上的模型参数来合成理论观测数据,正演问题也可以进一步表示为：
[0094]
f(m)＝f(α,β,ρ)＝d#(2.2)
[0095]
其中，m表示真实模型，d表示理论观测数据,f表示真实模型到观测数据的映射，α,β,ρ表示真实模型的三个模型参数，分别表示p波速度，s波速度，介质密度。
[0096]
从之前的论述可以看出，面波的正演分为两个步骤。因此，面波走时的反演问题也可以相应地分成两步:首先要利用面波走时反演每个周期的相速度、群速度频散图，之后再用一系列的一维反演来获得地下三维结构。
[0097]
第一步线性化反演可以写成：
[0098][0099]
其中，ω表示面波频率，是该周期(频率)下的第i个观测面波走时，t
i
(ω)是该周期下的第i个理论面波走时，δ0t
i
(ω)表示第i个观测面波走时和理论面波走时之间的残差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δc
j
(ω)表示自由表面第j个网格节点处观测频散速度值和理论频散速度值之间的差值，k
ij
(ω)是敏感核函数(模型参数的广义导数)。
[0100]
由于面波频散对s波的速度结构更敏感，可以用岩石物理中的经验关系将s波速度结构转换成密度和p波速度公式：
[0101]
α＝f
α
(β),ρ＝f
ρ
(β)#(2.4)
[0102]
其中，α表示p波速度，β表示s波速度，ρ表示密度，f
α
表示s波速度模型到p波速度的映射，f
ρ
表示s波速度到密度的映射。
[0103]
传统的两步法利用(2.3)式来反演频散图，之后再用(2.4)式来反演s波速度。但第一步反演可能会将误差传递到第二步反演中。为了克服上述问题，可以将(2.3)
‑
(2.4)式结合起来，得到：
[0104][0105]
其中，d/dβ是对s波速度的全偏导数，t
i
(ω)是该周期下的第i个理论面波走时，δt
i
(ω)表示第i个观测面波走时和理论面波走时之间的残差，c
j
(ω)是在自由表面第j个网格节点处的频散速度值，δβ表示真实模型和理论模型之间的s波速度差值。
[0106]
需要说明的是，在当前的局部化假设下，每一个走时只依赖于两个台站的之间的面波射线路径之下的弹性参数。因此，面波走时敏感核也是稀疏矩阵。利用(2.4)式，可以把密度的扰动量转换为对s波速度的扰动量：
[0107]
[0108]
其中，β表示s波速度，ρ表示密度，d/dβ是对s波速度的全偏导数，δβ表示真实模型和理论模型之间的s波速度差值。
[0109]
之后利用(2.5)式，可以将每一次迭代过程中面波和重力数据组装到一个稀疏线性系统中：
[0110][0111]
其中，g
s
表示面波数据理论算子，g
g
表示重力数据理论算子，δβ
i
表示真实模型和理论模型之间的s波速度差值，δt表示面波走时残差，δg表示重力异常残差。
[0112]
上式的解可以看作一个优化问题：
[0113][0114]
其中，σ和s分别是二阶吉洪诺夫正则化系数和二阶导数矩阵，γ是阻尼系数，i是单位矩阵，w
s
和w
g
分别是面波走时和重力数据的(对角)权重矩阵：
[0115][0116][0117]
其中，p是一个0
‑
1之间的常数，需要根据实际问题的特点来调整，n
s，g
分别是面波和重力的数据个数，σ
s，g
是两套数据集中理论模型与真实模型之间的标准差。(2.8)式可以在2范数下通过最小二乘算法高效地求解。之后，可以通过：
[0118]
β
i+1
＝β
i
+δβ
i
#(2.10)
[0119]
来更新模型，直到误差降到合理的范围内。
[0120]
图3示出了中国四川盆地地区检测板测试分析图，利用初始模型加上
±
10％的速度扰动构成的，在水平方向上的尺度大约是1.5
°×
1.5
°
，在垂直方向上，考虑到两种数据敏感核的有效范围，异常体仅放置在60km以内，利用检测板模型正演得到所有台站对走时和所有观测点的布格重力异常后，在正演的面波数据中加入标准差1.4s的高斯噪声，重力异常数据中加入10mgal的高斯噪声作为观测数据，之后利用和实际数据同样的方法进行反演。
[0121]
从该结果可以看出可以看出，在60km之内的各个深度范围内，面波单独反演总体上恢复了速度结构，但是在射线路径之外的分辨率有所下降。而联合反演在射线路径内部得到了和单独反演类似的结果，同时在浅层射线路径之外的模型分辨率也有了显著提高。
[0122]
为了研究联合反演方法相比单独反演在减少反演多解性上的优势，本发明进行了一个异常体合成测试。本发明在一个一维模型上的不同深度的位置加入了4个异常值为
±
0.4km/s的矩形异常体(图5)。之后用了和检测板测试同样的方法来恢复模型，反演结果如图4
‑
5所示。
[0123]
图4和图5给出了反演结果在两个深度上的水平剖面和在4条垂直剖面上的反演结果。可以看出，面波单独反演恢复了模型的主要特征，但同时也引入了很多假异常。这些异常可能是由于对带噪声的走时数据过拟合引起的。而联合反演的结果明显更加“干净”，更接近真实模型。
[0124]
本发明实施例提供了一种地震面波走时和重力异常联合反演方法。其中，方法包括：获取待反演的地震面波频散数据和布格重力异常数据；基于程函方程的直接面波层析成像和球坐标系下的自适应高斯
‑
勒让德数值积分来进行两套数据集相应正演计算和反演敏感核；利用岩石物理实验中的速度
‑
密度经验关系，将面波直接反演方法以及球坐标系下的高斯
‑
勒让德重力积分方法相结合；对观测到的面波走时数据和布格重力异常直接建模进行反演成像。在本发明的方法中，在正演建模时考虑了面波射线的几何形态以及地球的曲率对面波和重力数据的影响。理论合成测试表明，和面波单独反演相比，尤其是在重力有较强约束能力的区域内，本发明能获得更合理的结果。最后，将该方法在中国四川盆地区域进行测试分析，实施例验证了新方法反演地下速度结构的可靠性。从反演误差上看，联合反演的速度模型相对于单独反演大大减小了重力数据的残差，说明引入重力数据可以提高联合反演可靠性，反演更为准确的地下速度结构信息，可为深部矿产资源勘查和地震灾害预警提供先验的速度模型。
[0125]
本发明还提供了一种地震面波走时和重力异常联合反演系统，包括：
[0126]
数据集获取模块，用于获取待反演的地震面波频散数据和布格重力异常数据；
[0127]
正演模块，用于对所述待反演的地震面波频散数据和所述布格重力异常数据进行正演得到理论观测数据；
[0128]
计算模块，用于计算敏感核；
[0129]
反演模块，用于根据所述理论观测数据和所述敏感核进行反演得到线性化反演模型；
[0130]
线性化反演模型修正模块，用于基于岩石的速度
‑
密度关系对所述线性化反演模型进行修正得到修正后的线性化反演模型；
[0131]
联合反演目标函数构建模块，用于根据修正后的线性化反演模型建立联合反演目标函数；
[0132]
求解模块，用于对所述反演目标函数进行求解得到最优解；
[0133]
反演成像模块，用于利用所述最优解对待反演的地震面波进行反演成像。
[0134]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0135]
本发明提供了一种地震面波走时和重力异常联合反演方法与系统，本发明提供的一种地震面波走时和重力异常联合反演方法包括：首先对待反演的地震面波频散数据和布格重力异常数据进行正演得到理论观测数据并计算敏感核；然后根据理论观测数据和敏感核进行反演得到线性化反演模型；基于岩石的速度
‑
密度关系对线性化反演模型进行修正得到修正后的线性化反演模型；最后根据修正后的线性化反演模型建立联合反演目标函数；对反演目标函数进行求解得到最优解；利用最优解对待反演的地震面波进行反演成像。本发明通过利用地震面波频散数据和布格重力异常数据对地震面波进行联合反演成像，不
仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性。
[0136]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上，本说明书内容不应理解为对本发明的限制。