地震约束下的密度反演

地震约束下的密度反演
摘要：广义线性联合反演对其模型的先验了解依赖性很大，只有在给定初始猜测接近真实地球物理模型时才可以求出有意义的解估计。

地震波旅行时可以认为是目前最精确的地球物理探测量，本文在地震波旅行时的限制和纵波速度的约束下，利用常规最小二乘线性反演进行密度反演。

利用该方法得到的反演结果比较接近真实值，作为先验信息可初始化重震迭代联合反演的模型参数。

并且在速度参数接近准确的情况下，该反演结果比单由速度密度关系式得到的结果更接近真实值。

关键字：密度反演；旅行时限制；速度约束
Density Inversion under the Restriction of Seismic Data
Abstract: Joint inversion depends on much of the prior understanding of the model. Only when the given initial guess is close to the real physical model we can find out meaningful solution. We think that seismic waves travel time is the most accurate geophysical detecting quantity. Based on the restrictions of seismic waves travel time and the constraints of p-wave velocity, we can finish the density inversion using conventional least square linear inversion. The inversion results we get from the method as the priori information can initialize the parameters of the model used in iterative joint inversion of seismic and gravity data. When the velocity parameters are close to the accurate, the inversion results are close to the real value.
Key Words: density inversion; restrictions of travel time; constraints of velocity
0 引言
地球物理反演是在地球物理学中利用地球表面观测到的物理现象推测地球内部介质物理状态的空间变化及物性结构的一个分支[1]。

地球物理反演可分为单一地球物理现象的反演和多种地球物理现象的联合反演[2]。

把单一地球物理资料的反演变为反映各种不同地球物理场或地球物理资料的联合反演, 是当今地球物理资料反演的总趋势, 也是发展的必然[3]。

联合反演是地球物理综合研究的主要方法之一[4]。

对于缺乏测井信息的深部探测，地震信息与非震信息的联合更是显得尤为重要。

广义线性反演是目前联合反演中应用最广泛的反演方法，但是广义线性联合反演对其模型的先验信息依赖性很大，只有在给定初始猜测接近真实地球物理模型时才可以求出有意义的解估计[1]。

要想使模型参数良好的收敛于真实值，必须取得接近真实的模型参数。

地震波旅行时可以认为是目前最精确的地球物理探
测量，本文利用地震波旅行时的限制对地下介质进行时域区块结构划分，结合用资料处理过程中获得的速度信息得到空间域区块结构分布，并在速度-密度关系约束下进行密度反演。

本方法可作为重震联合反演的重要组成部分，也可作为联合反演问题研究的重要入口。

刘皓等人也曾在模型分块处理的基础上做了地震约束的重力反演，但他是在由速度-密度关系得出重力模型的基础上进行的参数调整[5]。

1 多边形重力正演
图1 多边形截面
Fig.1 Polygon section
人们通常使用Parker法进行异常正演计算，但是Parker法依赖于接近层状分布的地下介质模型，而不适用于块体或复杂结构异常计算。

因此，本方法利用多边形截面法计算横截面为任意形状二度体的重力异常，如图1所示。

首先计算点O与多边形各边组成的多个三角形，，···，所引起的重力异常，并以顺时针给各三角形赋以正、负号。

(1)
式中：。

整个二度体在O点引起的重力异常，即横截面，，···，等N个二度体在O点引起的重力异常总和，即
(2)
2 地震走时限制下的密度反演
地震波旅行时可以认为是目前最精确的地球物理探测量，可用其限制密度反演。

我们假设密度反演结果与旅行时完全一致，在旅行时的限制下进行密度反演。

这时的密度反演问题是一个完全线性问题。

首先，基于地震时域解释资料获得地下结构信息进行区块模型建立。

应用地下界面和断层等先验信息将地下介质划分成有限个介质区块，如图2，Zone-1~Zone-6是由地震波同相轴划分得到的。

再应用资料处理过程中获得的有效的均方根速度信息进行时深转换，将时域模型映射到深度域模型，同时也得到模型层速度。

垂直射线时深转换常用均方根速度Dix转换，这时我们获得了层速度截面。

利用层速度图进行时间层的垂直射线或成像射线转换，得到深度层位并组成了时深转换得到的初始模型。

Dix转换的基础公式为：
(4)
式中，为第（n-1）个地层界面和第n个地层界面限制的地层速度。

和为双程零偏移距时间, 和为相应的均方根速度。

图2 简单模型图
Fig.2 Simple Model
假设模型中基岩剩余密度为0，只考虑Zone-1~Zone-6的各个区块的反演问题。

则模型参数和观测数据分别为：
，(5)
(6)
其中，为各区块的密度，为各地面观测点的重力异常。

且有多边形截面正演公式可得：
(7)
即有：(8)
我们可用阻尼最小二乘法解这个线性反问题，即
(9)
其中阻尼系数的选择主要靠实验，它起到分辨率和方差之间的调节作用。

另外，对于地下结构分布复杂性不均匀时，还需考虑各探测段数据拟合情况的权重问题。

对于区块分布较为复杂，分块较多的探测段需要拟合权重有必要增加。

这时需要引入加权矩阵
(11)
3 速度约束下密度反演
我们知道，地震波信息对纵向介质差异反应比较敏感，而重力异常信息对横向介质差异反应比较敏感。

当模型区块分布偏重于纵向分布时，由重力数据所得到反演结果与真实值会有比较明显的偏差。

这时地震波信息对密度反演的约束显得尤为重要。

纵横波速度与介质密度存在着正相关关系，尤其是纵波速度更为明显。

人们也经常利用速度密度经验关系式来估算介质密度。

因此，利用模型速度参数约束
密度反演是个比较好的方法。

速度-密度关系式通常使用Birch线性关系式和Gardner(1974)指数关系式，本次我们选择Gardner指数关系式：
(12)
式中为密度，单位；为纵波速度，单位。

但是密度与速度经验公式有地区差异性，因此在应用之前应确定适合的经验关系式，即确定适合的公式参数a，b
(13)
结合式(8)，可得数据与模型的关系式
(14)
其中，
(15)
(16)
式中，为M阶单位矩阵。

和分别为重力数据与速度数据的加权矩阵，且均为对角阵。

加权矩阵的确定是个比较困难的问题，它取决于多个因素[10]。

首先要求对模型参数进行无量纲化或规格化，以及不同数据之间的相对加权。

另外，还要考虑模型响应对模型参数改变的灵敏度等先验信息。

4密度反演的实现
4.1 简单模型的密度反演
以图2所示简单模型为例，其正演重力异常如图3所示，该正演结果带有一定量的随机噪音。

测线设计在模型两端分别增加一段测点，以便更好的反映地下模型的分布情况。

图3 简单模型重力异常
Fig.3 Gravity anomaly of simple model
表1为模型反演结果数据，由反演结果可知，单一重力数据反演情况下Zone-2和Zone-5较小且较深的区块密度偏差较大。

而速度约束有力的限制了较小区块的的反演结果偏差。

但有些区块，速度约束并没有起到太大的作用，这与
速度加权矩阵的设计有着很大的关系。

表1 简单模型反演结果
Table 1 Inversion of simple model
区块编号模型参数
() 单一反演
() 速度约束
()
Zone-1 2.32 2.3356 2.3164
Zone-2 2.35 2.3918 2.3593
Zone-3 2.33 2.338 2.321
Zone-4 2.45 2.4541 2.4466
Zone-5 2.49 2.3919 2.4569
Zone-6 2.46 2.4977 2.4916
图4为密度反演结果与模型参数的对比分布图，其中虚线直线表示的是关系式，即模型参数与反演结果相等的直线。

统计点越靠近该直线就说明反演结果越接近真实值。

由图可知，速度约束使大部分统计点更接近其真实值。

图4 反演结果分布图
Fig.4 Distribution of Inversion
4.2 随机干扰对反演的影响
下面讨论随机干扰对反演的影响。

一般情况下随机噪音对较小区块的影响较大。

表2和图5为不同干扰情况下的反演结果对比，噪音对Zone-4和Zone-6较大区块影响很小。

因此，可以认为随机噪音是造成该方法多解性的重要原因之一。

在干扰信号较明显时，初始反演阶段的区块划分应倾向于大块划分，且反演前有必要对随机噪音进行过滤，平滑等处理，已获得最优的反演结果。

表2 随机干扰下反演结果
Table 2 Inversion under random noise
区块编号模型参数
() 低干扰
() 高干扰
()
Zone-1 2.32 2.3164 2.4092
Zone-2 2.35 2.3593 2.41
Zone-3 2.33 2.321 2.4442
Zone-4 2.45 2.4466 2.4422
Zone-5 2.49 2.4569 2.346
Zone-6 2.46 2.4916 2.4339
图5 反演结果分布图
Fig.5 Distribution of Inversion
4.3 速度误差对反演的影响
本方法是基于地震时域走时和纵波速度实现的，因此，纵波速度误差对反演结果的影响是一个十分重要的讨论点。

本文只对其做了简单的数字统计。

表3为当计算速度比真实速度均偏大2%时的反演结果，图6为其统计分布图。

由此可知，计算速度偏大时反演结果明显偏离模型参数，且大部分反演结果偏小。

速度误差使各区块空间分布偏离实际。

有关速度误差影响还有待进一步讨论。

表3 速度误差下反演结果
Table 3 Inversion under velocity error
区块编号模型参数
() 单一反演
() 速度约束
()
Zone-1 2.32 2.3082 2.2793
Zone-2 2.35 2.3878 2.3373
Zone-3 2.33 2.3455 2.3173
Zone-4 2.45 2.4105 2.3984
Zone-5 2.49 2.296 2.3981
Zone-6 2.46 2.4217 2.4133
图6 反演结果分布图
Fig.6 Distribution of Inversion
4.4 多区块模型密度反演
下面我们做一个较为复杂的区块模型密度反演。

模型图如图7所示。

其重力异常曲线如图8所示。

此曲线也带有一定量的随机干扰。

图7 多区块模型
Fig.7 Model with multi-zone
图8 多区块模型重力异常
Fig.8 Gravity anomaly of model with multi-zone
表4为该模型的密度反演结果。

由此可知，在这种模型情况下速度约束显得更加重要。

由于此模型多个区块均接近纵向分布，所以这些区块对重力异常反应不大。

图9为反演结果分布图，此次反演中，速度约束在权重比较小的情况下就可以很好的将反演结果约束在真实值附近。

表4 多区块模型反演
Table 4 Inversion of model with multi-zone
区块编号模型参数
() 单一反演
() 速度约束
()
Zone-1 2.23 2.2222 2.1948
Zone-2 2.35 2.5563 2.3934
Zone-3 2.52 2.1587 2.474
Zone-4 2.67 2.1475 2.6291
Zone-5 2.8 3.5271 2.8633
Zone-6 2.32 2.1345 2.2741
Zone-7 2.5 3.1833 2.5742
Zone-8 2.65 1.8588 2.6215
Zone-9 2.65 1.6387 2.5858
Zone-10 2.87 2.699 2.896
Zone-11 2.23 2.309 2.291
Zone-12 2.35 2.4812 2.2893
Zone-13 2.52 2.068 2.4333
Zone-14 2.67 1.9019 2.568
Zone-15 2.8 3.4612 2.7497
图9 反演结果分布图
Fig.9 Distribution of Inversion
5 结论
(1) 本方法的核心在于对于地震波走时探测认为是完全可靠的，而只对重力异常进行拟合。

本方法意在形成足够准确的模型参数进行下一步迭代联合反演，并且这种限制方法可以用于迭代联合反演。

(2) 不同类型的区块对重力异常有着不同程度的反应，因此速度约束在该方法中起着较大的作用。

加权矩阵的设计和阻尼系数的选择需要实验去选择。

参考文献
[1] 敬荣中，鲍光淑，陈绍裘. 地球物理联合反演研究综述[J].地球物理学进展，2003，18(3):535~540.
JING Rong-zhong, BAO Guang-shu, CHEN Shao-qiu. A review of the researches for geophysicalcombinative inversion[J]. Progress in Geophysics, 2003, 18(3):535~540.
[2] 杨文采. 地球物理反演和地震层析成象[M]. 北京: 地质出版社，1989.
YANG Wen-cai. Geophysical Inversion & Seismic Tomography[M]. Beijing: Geology Publishing House, 1989.
[3] 王家映. 地球物理反演理论[M]. 武汉: 中国地质大学出版社，1998.
WANG Jia-ying. Inversion Theory in Geophysics[M]. Wuhan: China University of Geosciences Press, 1998.
[4] 刘光鼎. 论地球科学[J]. 地学前缘. 1998，5(1-2): 1~8.
LIU Guang-ding. Theory of earth science[J]. 1998, 5(1-2).
[5] 刘皓，方盛明，嘉世旭. 地震-重力联合反演及效果——以天津-北京-赤诚地震探测剖面为例[J]. 地震学报，2011，33(4)：443~450.
LIU Hao, FANG Sheng-ming, JIA Shi-xu. Joint inversion of seismic and gravity data and its effect: Taking the Tianjin—Beijing—Chicheng seismic profile as an example[J]. Acta Seismologica Sinica, 2011, 33(4): 443-450.
[6] 孟令顺，杜晓娟. 勘探重力学与地磁学. 北京：地质出版社，2008.
MENG Ling-shun, DU Xiao-juan. Explorational Gravitation and magnetism. Beijing: Geology Publishing House, 2008.
[7] 渥·伊尔马滋. 地震资料分析——地震资料处理、反演和解释[M]. 北京：石油工业出版社，2006.
Yilmaz, O.. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data[M]. Beijing: Petroleum Industry Press, 2006.
[8] Gardner. G. H. F, Formation velocity and density——the diagnostic basics for stratigraphic traps. Geophysics. 1974, 39: 770-780.
[9] 朱广生，桂志先，熊新斌，等. 密度与纵横波速度关系[J]. 地球物理学报，1995,38(sup.1): 260~264.
ZHU Guang-sheng, GUI Zhi-xian, XIONG Xin-bin. Relationship between Density and P-wave Velocity[J]. Acta Geophysica Sinica, 1995, 38(sup.1):260~264.
[10] 李雄. 地震-重力联合反演不同的实现方式[J]. 天然气地球科学，1992，(5): 37~42.
LI Xiong. Joint Inversion of Seismic and Gravity Data: Different Implementation Ways[J]. Natural Gas Geoscience, 1992, (5): 37~42.。