CN108845353B

CN108845353B - 一种重震联合反演的方法及装置

Info

Publication number: CN108845353B
Application number: CN201810952725.9A
Authority: CN
Inventors: 赵崇进; 于鹏; 张罗磊
Original assignee: Tongji University
Current assignee: Tongji University
Priority date: 2018-08-21
Filing date: 2018-08-21
Publication date: 2020-01-03
Anticipated expiration: 2038-08-21
Also published as: CN108845353A

Abstract

本发明公开了一种重震联合反演的方法及装置，方法包括以下步骤：获取研究区的重力数据、地震观测数据、物性统计数据；建立研究区密度和速度的初始模型；建立一种新型结构耦合方式，构建重震联合反演的目标函数；并对目标函数进行优化并求解目标函数，进而获得研究区的密度和速度分布结果；图像绘制。装置包括数据采集器、处理器、显示器、输入设备；数据采集器内设置有数据获取模块；处理器内设置有数据计算模块、数据分析模块；显示器内设置有数据解释模块、数据成图模块；输入设备用于输入数据和参数。本发明克服了密度和速度的经验关系等人为因素的影响，获得数据可靠，适用性强，具有快速、高效、经济的优势。

Description

一种重震联合反演的方法及装置

技术领域

本发明涉及一种方法及装置，尤其涉及一种重震联合反演的方法及装置。

背景技术

地球物理反演的多解性以及单一地球物理方法解决问题的局限性，使得综合地球物理研究具有非常明显的重要性优势性，其中，联合反演是综合地球物理研究的重要定量解释工具。联合反演就是应用多种地球物理观测数据，通过地质体的岩石物性和几何参数之间的相互关系获得同一个地下地质以及地球物理模型。由于要推测的地球物理模型只有一个，它所产生的各种地球物理场应与地表观测到的多种地球物理数据相符合，从而才能达到改善反演的非唯一性，因此，联合反演长期以来一直是地球物理学科关注的焦点和难点问题。

然而，复杂的地表条件和复杂的地下结构使地球物理当前所面对的勘探目标愈加复杂和困难，已有的地球物理方法获得的反演结果都是依据某种岩石物性参数从某一个角度来认识地质地球物理问题，单凭一种地球物理方法解决问题表现出局限性，必须以其他方法作必要补充，并且已有的联合反演未能够建立合理的的耦合机制并未能减少反问题的多解性，具有着局限性、不稳定性和众多缺陷。此外，已有的联合反演对人为因素的依赖性强，获得的数据可靠性较差，对于复杂地质构造条件的适用性也较差，无法获得广阔的市场前景和良好的推广使用价值。

发明内容

为了解决上述技术所存在的不足之处，本发明提供了一种重震联合反演的方法及装置。

为了解决以上技术问题，本发明采用的技术方案是：一种重震联合反演的方法，包括以下步骤：

步骤一、获取研究区的重力数据、地震观测数据，收集和统计研究区物性统计数据；

步骤二、根据研究区的先验信息，建立研究区密度和速度的初始模型；

步骤三、基于研究区的区域模型矢量点积约束建立一种新型结构耦合方式，将密度和速度这两种物性连接，再将新型结构耦合方式引入到一般联合反演的目标函数中，然后采用交替耦合的方式重新构建重震联合反演的目标函数；

步骤四、采用交替耦合的方式和共轭梯度算法并对目标函数进行优化并求解目标函数，进而获得研究区的密度和速度分布结果；其中，共轭梯度算法是利用C++语言实现的；

步骤五、通过等值线涂色方式进行图像绘制，表现地下密度和速度的高低分布情况。

进一步地，步骤三中所述新型结构耦合方式是一种区域模型梯度矢量余弦相似度平方的新型结构耦合方式，如公式Ⅰ所示：

τ(m⁽¹⁾,m⁽²⁾)＝1-cos²θ₁₂ Ⅰ

其中，τ为重力和地震的耦合约束泛函，m⁽¹⁾表示密度模型矢量，m⁽²⁾表示速度模型矢量；cosθ₁₂为区域余弦相似度，cos²θ₁₂表示cosθ₁₂的平方，cosθ₁₂如公式Ⅱ所示：

其中，M⁽¹⁾＝m⁽¹⁾-m⁽¹⁾ _apr，表示密度模型参数的变化，M⁽²⁾＝m⁽²⁾-m⁽²⁾ _apr，表示速度模型参数的变化；m⁽¹⁾ _apr表示先验密度模型，m⁽²⁾ _apr表示先验速度模型，η⁽¹⁾表示密度模型的非零项，η⁽²⁾表示速度模型的非零项。

进一步地，步骤三中所述重震联合反演的目标函数如公式Ⅲ所示：

其中，P^α为目标函数，m⁽¹⁾表示密度模型矢量，m⁽²⁾表示速度模型矢量，表示重力的数据拟合泛函，

表示地震的数据拟合泛函，s₁表示重力的模型稳定泛函，s₂地震的模型约束泛函，γ₁表示重力的数据误差权重，γ₂表示地震数据误差的权重，α₁为重力的正则化因子，α₂为地震的正则化因子，τ为重力和地震的耦合约束泛函，β表示模型矢量耦合泛函的权重因子。

进一步地，步骤四中优化后的目标函数如公式Ⅳ所示：

将目标函数划分为多个交替进行的子反演系统，相对独立地进行目标函数优化，

其中，n-1表示第n-1次迭代，n表示当前模型的优化更新。

一种重震联合反演的装置，包括数据采集器、处理器、显示器、输入设备；数据采集器内设置有数据获取模块；处理器内设置有数据计算模块、数据分析模块；显示器内设置有数据解释模块、数据成图模块；输入设备用于输入数据和参数；

数据获取模块用于获取研究区的重力数据、地震观测数据、研究区物性统计数据；

数据计算模块用于重震联合反演的计算，以得到研究区的密度和速度分布结果；数据分析模块用于对重震联合反演计算结果进行聚类分析，以得到相应的分类结果；

解释模块用于根据研究区的物性统计数据将聚类结果中的不同类划分到不同的地层，进行聚类均值统计，以得到重震联合反演的地质综合解释结果；数据成图模块用于显示研究区的密度和速度的分布情况。

进一步地，数据获取模块采用网格差值算法对重力数据和地震观测数据进行分析处理成为重震联合反演所需的重力数据、地震观测数据；

数据成图模块采用了等值线的方式对地下密度和速度的分布进行显示。

本发明提供了一种重震联合反演的方法及装置，克服了密度和速度的经验关系等人为因素的影响，利用有限的地质和地震约束条件，得到可靠的密度和速度模型，实现对多数据的联合反演，特别是对于复杂地质构造条件具有很强的适用性。该发明快速、高效、经济，具有广阔的市场前景和良好的推广使用价值。

附图说明

图1为本发明装置的组成框图。

图2为本发明研究区的计算区域示意图。

图3本发明实施例中岩丘模型密度分布图。

图4本发明实施例中岩丘模型速度分布图。

图5为本发明实施例中岩丘模型获得的重力异常观测数据。

图6为本发明实施例中岩丘模型获得的初至走时观测数据。

图7为本发明实施例中岩丘模型获得的密度初始模型。

图8为本发明实施例中岩丘模型获得的速度初始模型。

图9为本发明实施例中岩丘模型对应的密度分布单一反演结果。

图10为本发明实施例中岩丘模型对应的速度分布单一反演结果。

图11为本发明实施例中岩丘模型基于重震联合反演的密度分布结果。

图12为本发明实施例中岩丘模型基于重震联合反演的速度分布结果。

图13为本发明实施例中岩丘模型的理论模型点积图。

图14为本发明实施例中岩丘模型的单一反演模型点积图。

图15为本发明实施例中岩丘模型的重震联合反演模型点积图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。

一种重震联合反演的方法，包括以下步骤：

步骤二、根据研究区的先验信息，建立研究区密度和速度的初始模型；其中先验信息一般包括地震资料信息、钻井信息、地质信息、岩石物性信息等；

步骤五、通过等值线涂色方式进行图像绘制，表现地下密度和速度的高低分布情况；其中，等值线涂色方式是基于软件surfer进行绘制的，采用的网格化方式是最近邻距(Nearest Neighber)或克里金(Kriging)的网格化方式。

联合反演中，前人对不同物性参数的耦合方式做了大量的研究，具有代表性的是具有物性之间经验关系的耦合方式，但该方式对于先验经验关系的依赖性太强，因此交叉梯度和梯度点积的耦合方式得到了比较广泛地应用，但由于其需要计算模型梯度且计算区域较小会导致反演的不稳定性等问题，针对目前耦合方式研究存在的问题，本发明基于研究区的区域模型矢量点积约束提出了一种区域模型梯度矢量余弦相似度平方的新型结构耦合方式，如公式Ⅰ所示：

τ(m⁽¹⁾,m⁽²⁾)＝1-cos²θ₁₂ Ⅰ

其中，τ为重力和地震的耦合约束泛函，m⁽¹⁾表示密度模型矢量，m⁽²⁾表示速度模型矢量；cosθ₁₂为区域余弦相似度，cos²θ₁₂表示cosθ₁₂的平方，cosθ₁₂又如公式Ⅱ所示：

其中，M⁽¹⁾＝m⁽¹⁾-m⁽¹⁾ _apr，表示密度模型参数的变化，M⁽²⁾＝m⁽²⁾-m⁽²⁾ _apr，表示速度模型参数的变化；m⁽¹⁾ _apr表示先验密度模型，m⁽²⁾ _apr表示先验速度模型；为避免分母为零出现奇异性的情况，对分母进行了|η⁽ⁱ⁾|＝max{|M⁽ⁱ⁾|,η_i ^min},η_i ^min>0的处理，当i＝1时，η⁽¹⁾表示密度模型的非零项，i＝2时，η⁽²⁾表示速度模型的非零项，η_i ^min以梯度模均值乘以估计噪音值作为小值的标准。

区域模型梯度矢量余弦相似度表示归一化后指定模型区域空间内两个矢量的点积，若两个矢量方向愈接近(同向或反向)，其余弦值趋于1或-1，可来表征两个模型空间方向接近或相似的程度。

针对计算区域的指定，只要两种物性在指定的模型空间区域内满足相关或无关条件即可，即可以将整个计算模型空间划分为几个不同变化的区域，而这几个区域只需要保证物性变化的相关或无关条件。在二维情况下，假设计算模型区域如图2所示，横向范围网格数为a～b，纵向范围网格数为c～d。

将区域余弦相似度新型结构耦合方式引入到一般联合反演的目标函数中，得到重震联合反演的目标函数，如公式Ⅲ所示：

其中，P^α为目标函数，m⁽¹⁾表示密度模型矢量，m⁽²⁾表示速度模型矢量，

表示重力的数据拟合泛函，

表示地震的数据拟合泛函，s₁表示重力的模型稳定泛函，s₂地震的模型约束泛函，γ₁表示重力的数据误差权重，γ₂表示地震数据误差的权重，α₁为重力的正则化因子，α₂为地震的正则化因子，τ为重力和地震的耦合约束泛函，β表示模型矢量耦合泛函的权重因子，按照高精度约束低精度的原则，权重因子β可以作为控制不同物性模型之间的重要性大小的指标，来体现不同约束对联合反演的贡献。

考虑到除了耦合项以外，与各类方法有关的目标函数项仍是相互独立的，因此在实现具体最优化过程中利用交替耦合的方式对目标函数进行优化分步，即将目标函数按各方法划分为多个交替进行的子反演系统，各方法相对独立地进行目标函数优化，如公式Ⅳ所示：

其中，n-1表示第n-1次迭代，n表示当前模型的优化更新。

该优化方案是使各方法的数据误差拟合项和模型稳定泛函项相互独立，而模型耦合项则是多方法之间的纽带，交替式同步进行计算，即各方法的初始模型或第n-1次模型参数需要依次输入到各自的子反演系统进行基于结构耦合的反演优化以更新第n次的新模型，然后按收敛条件判断是否满足迭代结束条件，若满足则输出结果，若不满足则返回更新。该优化方式的特点是避免了公式Ⅲ优化过程的复杂性和不确定性，每个反演子系统对不同方法可以选取不同的权重因子，这样更适合各方法按自身的反演能力反映其对整体反演的贡献影响，在具体实施联合反演过程中也容易操作和实现。

本发明还提供了一种重震联合反演的装置，如图1所示，包括数据采集器、处理器、显示器、输入设备；数据采集器内设置有数据获取模块；处理器内设置有数据计算模块、数据分析模块；显示器内设置有数据解释模块、数据成图模块；输入设备用于输入数据和参数；

其中，数据获取模块采用网格差值算法对重力数据和地震观测数据进行分析处理成为重震联合反演所需的重力数据、地震观测数据；

下面结合具体的实施例，对本发明所提供的方法做进一步说明。

一、获取研究区的重力数据和地震观测数据

根据现有的研究区模型(一种岩丘模型)，该岩丘模型的密度分布如图3所示、速度分布如图4所示(对于各个密度分布或速度分布图，其右侧的标尺表示不同颜色表示的是不同密度或速度的值)，获得该研究区的重力异常和初至走时数据(走时数据为地震观测数据的一种)，分别如图5、图6所示。该模型的目标区域为30km×6km(横向为30km，纵向为6km)，网格个数为150×30，横纵向网格间距相同均为200m。

二、收集和统计研究区物性统计，根据研究区的先验信息，建立研究区密度和速度的初始模型，故根据已有的先验信息我们得到了密度和速度的初始模型，分别如图7、图8所示。可以看出在研究区的密度和速度整体呈由低到高的物性分布。

三、首先对该模型分别用重、震两种方法进行单一反演，重力采用了共轭梯度法反演、地震采用了基于快速扫描的伴随状态法初至层析成像，其结果分别如图9、10所示，并统计了三方法的数据拟合误差RMS和模型还原误差RMS，见表1。单一重力反演结果的整体密度结构和真实模型较为类似，但是对于x方向15km下方的岩丘，并没有任何的体现；单一磁法反演结果的磁性体位置和磁化强度均有较好地反映，但其形态却为一个椭圆形，这是由于磁法等位场方法是一个体积的效应对异常体的形态刻画能力不强，和真实的岩丘形态差距较大。对于速度反演结果而言，不仅各层的速度又较好的刻画，且岩丘的位置和形态也有较好的反映。

表1单一重震反演的数据拟合误差的RMS和模型差值的RMS

四、基于区域模型矢量点积约束建立新型结构耦合方式，采用交替耦合方式构建重震联合反演的目标函数，采用交替耦合方式和共轭梯度法优化并求解目标函数，进而获得研究区的密度和速度分布。

首先对目标函数中的重力反演目标函数求偏导，由公式II和公式IV可以得到公式Ⅴ，如下所示：

其中，δ表示对函数求偏导，M表示整个模型空间，d⁽¹⁾表示重力的观测数据，A表示重力的正演算子，A^*表示其伴随算子，W_d为数据灵敏度矩阵，W_m为综合灵敏度矩阵，且W_m＝diag(A^TA)^1/2，W_d ²表示W_d的平方，W_m ²表示W_m的平方；K⁽¹²⁾的具体表示如公式Ⅵ所示：

K⁽¹²⁾＝a((m⁽²⁾-m^(2)apr)-b(m⁽¹⁾-m^(1)apr)) VI

其中，a、b分别表示K⁽¹²⁾的系数，a、b的计算如下所示：

其中，j表示第j个模型网格，e＞0为一小值，N_m为模型网格的总数。

使目标函数，即公式IV达到最小，可以得到式子VII：

进而得到模型修正的梯度方向I⁽¹⁾，如公式VIII所示：

其中，

为指数稳定泛函。

运用共轭梯度的迭代流程，即可实现重震联合反演中重力方法目标函数的优化，具体流程如下：

其中，n表示第n次迭代，r表示数据残差，m表示模型矢量，β^α表示共轭梯度系数，

表示共轭梯度方向，

表示初始共轭梯度方向，I表示梯度方向，I₀表示初始梯度方向，k^α表示共轭梯度法中的步长。

对目标函数中的地震反演目标函数求偏导，又由公式Ⅱ和公式Ⅳ可以得到公式IX，如下所示：

其中，λ表示速度模型的拉格朗日算子。

进而得到速度模型修正的梯度方向I⁽²⁾，如公式x所示：

那么速度的修正公式为，如公式XI所示：

其中，k_n为步长，然后又由公式x可以得到，伴随状态下的层析成像的速度扰动的计算方程，如公式xii所示：

其中，▽为梯度运算，T为初至走时，n为▽T的法方向向量，T_obs为观测初至走时，而对于式中的权重因子β，利用自适应正则化因子的选取方式，依据数据拟合误差、模型约束误差以及联合反演耦合项这三项的值的大小，来确定相对的权重因子初始值以使这三项的大小相对平衡，然后根据反演的数据拟合误差变化情况来对权重因子做一定的调节。

基于本发明的重震联合反演方法对该模型进行试验，其反演结果如图11、12所示，并统计了其数据拟合误差RMS和模型还原误差RMS，见表2。

表2联合重震反演的数据拟合误差的RMS和模型差值的RMS

密度反演结果不仅显示明显密度分层，且对岩丘的形态也有较好的刻画；磁反演结果不仅对磁性体的位置和磁化强度有较好的反映，且盐丘的形态也有很好的刻画。速度反演结果在本次反演试验中模型还原误差RMS有所改善。

根据公式Ⅰ计算了盐丘模型的理论梯度点积，如图13、14、15所示，模型矢量耦合方式加强了对模型的刻画，从图中可以看出分界面处明显加强。

统计单一和联合反演结果的物性，并作出其物性相关图。结果显示，速度和密度的关系具有更强的线性变化，且更符合模型的变化趋势；而单一反演反应的却是随速度增大，磁化强度变小，不符合模型设计。

上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。

Claims

1.一种重震联合反演的方法，其特征在于：所述方法包括以下步骤：

步骤三、基于研究区的区域模型矢量点积约束建立一种新型结构耦合方式，将密度和速度这两种物性连接，再将新型结构耦合方式引入到一般联合反演的目标函数中，然后采用交替耦合的方式重新构建重震联合反演的目标函数；其中，新型结构耦合方式是一种区域模型梯度矢量余弦相似度平方的新型结构耦合方式，如公式Ⅰ所示：

τ(m⁽¹⁾,m⁽²⁾)＝1-cos²θ₁₂ Ⅰ

其中，M⁽¹⁾＝m⁽¹⁾-m⁽¹⁾ _apr，表示密度模型参数的变化，M⁽²⁾＝m⁽²⁾-m⁽²⁾ _apr，表示速度模型参数的变化；m⁽¹⁾ _apr表示先验密度模型，m⁽²⁾ _apr表示先验速度模型，η⁽¹⁾表示密度模型的非零项，η⁽²⁾表示速度模型的非零项；

步骤四、采用交替耦合的方式和共轭梯度算法并对目标函数进行优化并求解目标函数，进而获得研究区的密度和速度分布结果；其中，共轭梯度算法是利用C++语言实现的；其中，优化后的目标函数如公式Ⅳ所示：

其中，n-1表示第n-1次迭代，n表示当前模型的优化更新。

2.根据权利要求1所述的重震联合反演的方法，其特征在于：步骤三中所述重震联合反演的目标函数如公式Ⅲ所示：

3.根据权利要求1中所述的重震联合反演的方法，其特征在于：所述方法所使用的装置包括数据采集器、处理器、显示器、输入设备；所述数据采集器内设置有数据获取模块；所述处理器内设置有数据计算模块、数据分析模块；所述显示器内设置有数据解释模块、数据成图模块；所述输入设备用于输入数据和参数；

所述数据获取模块用于获取研究区的重力数据、地震观测数据、研究区物性统计数据；

所述数据计算模块用于重震联合反演的计算，以得到研究区的密度和速度分布结果；所述数据分析模块用于对重震联合反演计算结果进行聚类分析，以得到相应的分类结果；

所述解释模块用于根据研究区的物性统计数据将聚类结果中的不同类划分到不同的地层，进行聚类均值统计，以得到重震联合反演的地质综合解释结果；所述数据成图模块用于显示研究区的密度和速度的分布情况。

4.根据权利要求3所述的重震联合反演的方法，其特征在于：所述数据获取模块采用网格差值算法对重力数据和地震观测数据进行分析处理成为重震联合反演所需的重力数据、地震观测数据；

所述数据成图模块采用了等值线的方式对地下密度和速度的分布进行显示。