CN112596106B - 一种球坐标系下重震联合反演密度界面分布的方法 - Google Patents

Abstract

本发明提供了一种球坐标系下重震联合反演密度界面分布的方法，该方法包括以下步骤：估算初始的密度界面深度，计算修正后的密度界面深度迭代之后重力异常矩阵值计算，计算已知地震约束点与反演计算点的差值Δh_i，均方根误差分析，多次迭代计算及迭代终点判断。本发明的有益效果在于：该方法针对于区域大尺度的密度界面反演研究，将地球近似为球体，用球单元体替代常规的棱柱体，克服了大尺度研究忽略地球曲率的问题，反演结果更接近真实情况，同时反演过程中引入了深度加权因子和地震约束，提高了反演精度，减少了密度界面反演的多解性和不确定性。

Description

一种球坐标系下重震联合反演密度界面分布的方法

技术领域

本发明属于地球物理勘探技术领域，具体涉及一种球坐标系下重震联合反演密度界面分布的方法。

背景技术

在地球物理勘探领域，密度界面反演在推断区域构造研究、获取结晶基底面和莫霍面起伏等方面具有重要意义。密度界面反演是指根据重力异常变化，反演密度界面深度。现阶段界面深度的反演方法有很多，大体上可分为空间域和频率域方法。空间域主要包含直接迭代法、压缩质面法和非线性反演方法，空间域反演界面一般具有不连续性且计算速度慢；频率域密度界面算法引进了傅里叶变换，提高了计算效率。

目前小范围区域Dec密度界面反演方法已经比较成熟，该方法采用的是在笛卡尔坐标系下将地下模型剖分为规则的棱柱体。但是，当涉及到区域性大尺度乃至全球尺度问题时(譬如区域莫霍面反演、区域构造地质解释等)，由于需要将地球曲率考虑进去(Anderson.1976；Grombein et al.,2013；Uieda at al.,2016，Aowei Hao et al.2019)，此时采用上述方法得到的结果与真实情况之间存在一定的偏差。而且重力界面反演存在多解性也是现有技术难以解决的一个技术问题。

发明内容

为了弥补现有技术的不足，本发明提供了一种球坐标系下重震联合反演密度界面分布的方法，该方法针对于区域大尺度的密度界面反演研究，将地球近似为球体，用球单元体替代常规的棱柱体，克服了大尺度研究忽略地球曲率的问题，反演结果更接近真实情况，同时反演过程中引入了深度加权因子和地震约束，提高了反演精度，减少了密度界面反演的多解性和不确定性。

一种球坐标系下重震联合反演密度界面分布的方法，包括以下步骤：

(1)估算初始的密度界面深度：根据地质或相关物性资料，给出初始的密度界面参考深度,即重力异常基准面的深度值，给出待反演的球坐标系下重力异常矩阵为g_obs(q)，同时给出已知地震约束点数据将地下密度层在经度和纬度上剖分为M和N个球单元体，其中/>为在经度方向第i个、纬度方向第j个球单元体在第k次迭代之后相对于重力异常基准面的深度值；球单元体观测矩阵共有M×N个点,将每个测点的重力异常近似为无限平板产生的，则可以估算初始密度界面深度/>

其中G为万有引力常量，Δρ为密度层上下界面的密度差值；

(2)计算修正后的密度界面深度根据估算的所述初始密度界面深度/>采用球坐标下起伏密度界面重力异常正演方法，获得球坐标系下每个测点的重力异常值，进而获得初始重力异常矩阵为/>此时重力的残差/>在所述初始密度界面深度的基础上，通过引入球坐标下界面反演深度加权函数w(Δr)，并基于地震数据约束方法引入已知的地震点约束信息/>分别对重力迭代修正项和地震数据修正项进行修正，得到第1次迭代后的密度界面深度/>

同理，第k次迭代的重力异常残差在对重力迭代修正项和地震数据修正项进行修正后，得到第k次迭代后的密度界面深度/>

其中，α和β分别为重力和地震的收敛约束项，在实际测试中的一般范围为α∈(11.5),β∈(0.751.25)。

(3)第k次迭代后的重力异常矩阵值计算：采用球坐标下起伏密度界面重力异常正演方法，将第k次迭代后的所述密度界面深度/>带入球坐标下起伏密度界面重力异常正演公式，更新得到第k次迭代后的重力异常矩阵值/>

(4)计算已知地震约束点与反演计算点的差值Δh_i：将所述已知的地震约束点和反演计算的深度点/>带入下述公式进行差值Δh_i计算：

其中和/>分别为已知的地震约束点和反演计算的深度点，i∈n；其中反演计算的深度点/>可通过每次迭代的密度界面深度/>获取；

(5)均方根误差分析：每次迭代之后计算重力异常矩阵值与待反演的球坐标系下重力异常矩阵g_obs(q)之间的重力拟合均方根误差RMS_gra，并计算已知地震约束点与反演计算点的差值Δh_i的地震约束点拟合均方根误差RMS_seis，公式如下

(6)多次迭代计算及迭代终点判断：当所述重力拟合均方根误差RMS_gra和所述地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限时，迭代终止，输出密度界面反演结果，此时的可以认为是在经度方向第i个，纬度方向第j个球单元体相对于重力基准面的密度界面深度值；否则，重复步骤(2)-(5)进行迭代计算，直至重力拟合均方根误差RMS_gra和地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限，迭代终止，输出密度界面反演结果。

进一步的，步骤(2)中所述球坐标下起伏密度界面重力异常正演方法为：将地下密度层在经度和纬度上剖分为M和N个球单元体，测点P为观测面上任意测点，由地下岩层在测点P引起的球坐标下起伏密度界面重力异常g_p为：

其中，Forward(Δr_ij)为在经度方向第i个，纬度方向第j个球单元体引起的球坐标下起伏密度界面重力异常。

进一步的，步骤(2)中所述球坐标下界面反演深度加权函数w(Δr)的计算原理及方法为：

重力异常随着深度的加深迅速衰减，重力信号对深部不敏感，引入深度加权函数是为了平衡不同深度对于地表重力观测的权重，从而使反演结果更加接近于真实地下模型空间，纠正反演不同深度引起的畸变，使反演能够适应界面起伏剧烈的情况，球坐标系下界面偏离基准面的球坐标下界面反演深度加权函数：

其中η为深度加权的收缩因子，一般为1或2，H为球坐标下相对于观测面的深度，r₀为参考界面深度半径,r为实际界面起伏半径，Δr为实际界面起伏半径相对于参考界面深度半径的差值，沿着半径方向为正值，其中r＝r₀+Δr；在经度和纬度方向相同尺度的单元体，其只存在半径方向的差异，单元体的体积可近似为球单元体各陵边的乘积，即公式可简化为

当Δr为0是，实际界面与参考界面重合时，w(Δr)＝1；当Δr＞0时，实际界面位于参考界面上方，此时w(Δr)＜1，对浅部的界面隆起进行压制；当Δr＜0时，实际界面位于参考界面下方，此时w(Δr)＞1，对深部的凹陷进行突出。

进一步的，步骤(2)中所述地震数据约束方法为：

假定界面模型的差异和重力异常的差异存在线性的关系：

其中为实际界面和实际计算模型之间的修正值，Δg^M×N为理论重力异常和计算重力异常的残差，a和b为线性拟合的系数；

假定区域内共有n个地震测点作为深度约束点，通过最小二乘法线性回归可知：

其中和/>分别为第i个地震约束点的观测重力异常和计算重力异常，/>和分别为已知的地震深度点和反演计算的深度点，其中i∈n。

本发明的有益效果在于：

①本发明提供的一种球坐标系下重震联合反演密度界面分布的方法通过引入球坐标下界面反演深度加权函数的方法，平衡不同深度对于地表重力观测的权重，从而使反演结果更加接近于真实地下模型空间，纠正反演不同深度引起的畸变，使反演能够适应界面起伏剧烈的情况。

②本发明提供的一种球坐标系下重震联合反演密度界面分布的方法通过基于地震数据约束方法引入已知的地质等相关信息得约束可以有效地降低多解性的存在。相对重力而言，地震结果可以很好的反演地下界面信息，提高反演准确性，常规的已知地震深度点约束采用“硬约束”，假定已知点附近的范围内都为已知深度，当地震深度数据不准确时容易长生局部的形变，采用“软约束”，不追求某一个或者某几个约束点的正确性，而追求整体反演模型的最优化。

③本发明提供了一种球坐标系下重震联合反演密度界面分布的方法，该方法针对于区域大尺度的密度界面反演研究，将地球近似为球体，用球单元体替代常规的棱柱体，克服了大尺度研究忽略地球曲率的问题，反演结果更接近真实情况，同时反演过程中引入了深度加权因子和地震约束，提高了反演精度，减少了密度界面反演的多解性和不确定性。

附图说明

图1是本发明实施例1中球坐标下起伏密度界面重力异常正演方法中起伏密度层剖分为球单元体的二维剖面示意图；

图2是本发明实施例1中的一种球坐标系下重震联合反演密度界面分布的方法的步骤流程图；

图3是本发明实施例2中一种球坐标系下重震联合反演密度界面分布方法的测试结果示意图；其中(a)为理论模拟的莫霍界面分布图；(b)为球坐标下起伏密度界面重力异常正演方法得到的理论重力异常分布图；(c)为迭代次数与重力拟合均方根误差RMS_gra的收敛曲线图；(d)为迭代次数和地震约束点的RMSseis收敛曲线图；(e)迭代7次后地震约束点在界面结果和中异常结果的差异比对图；(f)为为已知重力计算异常值和理论重力异常之间的残差δg_k示意图；(g)为理论密度界面分布结果和反演密度界面结果的残差图；(h)为迭代7次后密度界面反演结果分布结果。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

1.球坐标下起伏密度界面重力异常正演原理和方法

将地下密度层在经度和纬度上剖分为M和N个球单元体，测点P为观测面上任意测点(起伏密度层剖分为球单元体的二维剖面示意图参见图1)，由地下岩层在测点P引起的球坐标下起伏密度界面重力异常g_p为：

其中，Forward(Δr_ij)为在经度方向第i个，纬度方向第j个球单元体引起的球坐标下起伏密度界面重力异常。

2.球坐标下界面反演深度加权函数w(Δr)的计算原理及方法

重力异常随着深度的加深迅速衰减，重力信号对深部不敏感，常规的反演会长生“趋肤效应”，引入深度加权函数就是为了平衡不同深度对于地表重力观测的权重，从而使反演结果更加接近于真实地下模型空间，纠正反演不同深度引起的畸变，使反演能够适应界面起伏剧烈的情况，球坐标系下界面偏离基准面的球坐标下界面反演深度加权函数：

其中η为深度加权的收缩因子，一般为1或2，H为球坐标下相对于观测面的深度，r₀为参考界面深度半径,r为实际界面起伏半径，Δr为实际界面起伏半径相对于参考界面深度半径的差值，沿着半径方向为正值，其中r＝r₀+Δr；在经度和纬度方向相同尺度的单元体，其只存在半径方向的差异，单元体的体积可近似为球单元体各陵边的乘积，即公式可简化为

当Δr为0是，实际界面与参考界面重合时，w(Δr)＝1；当Δr＞0时，实际界面位于参考界面上方，此时w(Δr)＜1，对浅部的界面隆起进行压制；当Δr＜0时，实际界面位于参考界面下方，此时w(Δr)＞1，对深部的凹陷进行突出，通过球坐标系下深度加权函数可以压制反演时的浅部隆起和深部凹陷，提高界面反演的准确性。

3.地震数据约束的原理和方法

重力界面反演存在多解性，加入已知的地质等相关信息得约束可以有效地降低多解性的存在。相对重力而言，地震结果可以很好的反演地下界面信息，提高反演准确性。常规的已知地震深度点约束采用“硬约束”，假定已知点附近的范围内都为已知深度，当地震深度数据不准确时容易长生局部的形变，本次我们采用“”软约束”，不追求某一个或者某几个约束点的正确性，追求整体反演模型的最优化。

假定界面模型的差异和重力异常的差异存在线性的关系：

其中为实际界面和实际计算模型之间的修正值，Δg^M×N为理论重力异常和计算重力异常的残差，a和b为线性拟合的系数；

假定区域内共有n个地震测点作为深度约束点，通过最小二乘法线性回归可知：

其中和/>分别为第i个地震约束点的观测重力异常和计算重力异常，/>和分别为已知的地震深度点和反演计算的深度点，其中i∈n。

4.一种球坐标系下重震联合反演密度界面分布的方法，包括以下步骤(步骤流程图参见图2)：

(1)估算初始的密度界面深度：根据地质或相关物性资料，给出初始的密度界面参考深度,即重力异常基准面的深度值，给出待反演的球坐标系下重力异常矩阵为g_obs(q)，同时给出已知地震约束点数据将地下密度层在经度和纬度上剖分为M和N个球单元体，其中/>为在经度方向第i个、纬度方向第j个球单元体在第k次迭代之后相对于重力异常基准面的深度值；球单元体观测矩阵共有M×N个点,将每个测点的重力异常近似为无限平板产生的，则可以估算初始密度界面深度/>

其中G为万有引力常量，Δρ为密度层上下界面的密度差值；

(2)计算修正后的密度界面深度根据估算的初始密度界面深度/>采用球坐标下起伏密度界面重力异常正演方法，获得球坐标系下每个测点的重力异常值，进而获得初始重力异常矩阵为/>此时重力的残差/>在初始密度界面深度/>的基础上，通过引入球坐标下界面反演深度加权函数w(Δr)，并基于地震数据约束方法引入已知的地震点约束信息/>分别对重力迭代修正项和地震数据修正项进行修正，得到第1次迭代后的密度界面深度/>

同理，第k次迭代的重力异常残差在对重力迭代修正项和地震数据修正项进行修正后，得到第k次迭代后的密度界面深度/>

其中，α和β分别为重力和地震的收敛约束项，在实际测试中的一般范围为α∈(11.5),β∈(0.751.25)。

(3)第k次迭代后的重力异常矩阵值计算：采用球坐标下起伏密度界面重力异常正演方法，将第k次迭代后的密度界面深度/>带入球坐标下起伏密度界面重力异常正演公式，更新得到第k次迭代后的重力异常矩阵值/>

(4)计算已知地震约束点与反演计算点的差值Δh_i：将已知的地震约束点和反演计算的深度点/>带入下述公式进行差值Δh_i计算：

其中和/>分别为已知的地震约束点和反演计算的深度点，i∈n；其中反演计算的深度点/>可通过每次迭代的密度界面深度/>获取；

(5)均方根误差分析：每次迭代之后计算重力异常矩阵值与待反演的球坐标系下重力异常矩阵g_obs(q)之间的重力拟合均方根误差RMS_gra，并计算已知地震约束点与反演计算点的差值Δh_i的地震约束点拟合均方根误差RMS_seis，公式如下

(6)多次迭代计算及迭代终点判断：当重力拟合均方根误差RMS_gra和地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限时，迭代终止，输出密度界面反演结果，此时的可以认为是在经度方向第i个，纬度方向第j个球单元体相对于重力基准面的密度界面深度值；否则，重复步骤(2)-(5)进行迭代计算，直至重力拟合均方根误差RMS_gra和地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限，迭代终止，输出密度界面反演结果。

5.采用上述一种球坐标系下重震联合反演密度界面分布方法的进行密度界面分布测试

假定研究区范围为0°-15°E，0°-10°N，上下密度界面的差值为280kg/m³，构建理论莫霍面界面深度见图3(a)，其变化范围在30.2km至50.6km之间，模型界面的起伏变化整体较大，在模型西部莫霍面较深，在东部莫霍面较浅。模型剖分和重力观测网的经度和纬度的网格间距为0.1°，计算网格共包含151×101个节点，重力异常的参考基准面选取莫霍面埋深的平均深度40km，重力观测面是高度为0的曲面，利用上述一种球坐标系下重震联合反演密度界面分布方法中的球坐标下起伏密度界面重力异常正演方法得到的理论重力异常见图3(b)，在莫霍面较浅的位置，底部高密度岩体抬升或侵入引起高的异常，相反地，在莫霍面较深的地方，莫霍面上部低密度体下降引起低的异常，其重力异常的范围为-93.43mGal至97.02mGal。

我们将理论的重力异常作为观测的重力异常进行空间域迭代反演，其迭代次数与重力拟合均方根误差RMS_gra的收敛曲线见图3(c)，迭代次数与地震约束点拟合均方根误差RMS_seis的收敛曲线见图3(d)。曲线图显示结果为：开始时快速收敛，在迭代6次之后趋于稳定。

本次实验取迭代7次后的结果(重力拟合均方根误差RMS_gra为0.032mGal,地震约束点拟合均方根误差RMS_seis为87.5m)；图3(e)展示的是迭代7次后地震约束点在界面结果和中异常结果的差异比对图，横轴为Δh_i，是界面反演结果和已知理论值的残差，其深度误差范围在250m之内，纵轴为Δg_i,它是已知地震约束点的重力理论值和计算值的残差结果，其范围在0.4mGal之内；图3(f)为已知重力计算异常值和理论重力异常之间的残差δg_k，其残差拟合结果整体较好，可控制在0.68mGal之内；图3(g)为理论密度界面分布结果和反演密度界面结果的残差图，即其残差可控制在650m之内；图3(h)为该方法迭代7次后的密度界面反演结果/>通过测试验证了上述一种球坐标系下重震联合反演密度界面分布方法的有效性，其密度界面反演结果精度高。在反演过程中以重力界面反演为主，同时加入已知地震点的“软约束”，降低反演的多解性。

应当理解，以上所描述的具体实施例仅用于解释本发明，并不用于限定本发明。由本发明的精神所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。

Claims (4)

1.一种球坐标系下重震联合反演密度界面分布的方法，其特征在于，所述方法包括以下步骤：

(1)估算初始的密度界面深度：根据地质或相关物性资料，给出初始的密度界面参考深度,即重力异常基准面的深度值，给出待反演的球坐标系下重力异常矩阵为g_obs(q)，同时给出已知地震约束点数据将地下密度层在经度和纬度上剖分为M和N个球单元体，其中为在经度方向第i个、纬度方向第j个球单元体在第k次迭代之后相对于重力异常基准面的深度值；球单元体观测矩阵共有M×N个点,将每个测点的重力异常近似为无限平板产生的，则可以估算初始密度界面深度/>

其中G为万有引力常量，Δρ为密度层上下界面的密度差值；

(2)计算修正后的密度界面深度根据估算的所述初始密度界面深度/>采用球坐标下起伏密度界面重力异常正演方法，获得球坐标系下每个测点的重力异常值，进而获得初始重力异常矩阵为/>此时重力的残差/>在所述初始密度界面深度/>的基础上，通过引入球坐标下界面反演深度加权函数w(Δr)，并基于地震数据约束方法引入已知的地震点约束信息/>分别对重力迭代修正项和地震数据修正项进行修正，得到第1次迭代后的密度界面深度/>

同理，第k次迭代的重力异常残差在对重力迭代修正项和地震数据修正项进行修正后，得到第k次迭代后的密度界面深度/>

其中，α和β分别为重力和地震的收敛约束项，在实际测试中的范围为α∈(1 1.5),β∈(0.75 1.25)；

(3)第k次迭代后的重力异常矩阵值计算：采用球坐标下起伏密度界面重力异常正演方法，将第k次迭代后的所述密度界面深度/>带入球坐标下起伏密度界面重力异常正演公式，更新得到第k次迭代后的重力异常矩阵值/>

(4)计算已知地震约束点与反演计算点的差值Δh_i：将所述已知的地震约束点和反演计算的深度点/>带入下述公式进行差值Δh_i计算：

其中和/>分别为已知的地震约束点和反演计算的深度点，i∈n；其中反演计算的深度点/>可通过每次迭代的密度界面深度/>获取；

(5)均方根误差分析：每次迭代之后计算重力异常矩阵值与待反演的球坐标系下重力异常矩阵g_obs(q)之间的重力拟合均方根误差RMS_gra，并计算已知地震约束点与反演计算点的差值Δh_i的地震约束点拟合均方根误差RMS_seis，公式如下

(6)多次迭代计算及迭代终点判断：当所述重力拟合均方根误差RMS_gra和所述地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限时，迭代终止，输出密度界面反演结果，此时的可以认为是在经度方向第i个，纬度方向第j个球单元体相对于重力基准面的密度界面深度值；否则，重复步骤(2)-(5)进行迭代计算，直至重力拟合均方根误差RMS_gra和地震约束点拟合均方根误差RMS_seis达到可忽略的极小值或者迭代次数达到上限，迭代终止，输出密度界面反演结果。

2.根据权利要求1所述的一种球坐标系下重震联合反演密度界面分布的方法，其特征在于，步骤(2)中所述球坐标下起伏密度界面重力异常正演方法为：将地下密度层在经度和纬度上剖分为M和N个球单元体，测点P为观测面上任意测点，由地下岩层在测点P引起的球坐标下起伏密度界面重力异常g_p为：

其中，Forward(Δr_ij)为在经度方向第i个，纬度方向第j个球单元体引起的球坐标下起伏密度界面重力异常。

3.根据权利要求2所述的一种球坐标系下重震联合反演密度界面分布的方法，其特征在于，步骤(2)中所述球坐标下界面反演深度加权函数w(Δr)的计算原理及方法为：

重力异常随着深度的加深迅速衰减，重力信号对深部不敏感，引入深度加权函数是为了平衡不同深度对于地表重力观测的权重，从而使反演结果更加接近于真实地下模型空间，纠正反演不同深度引起的畸变，使反演能够适应界面起伏剧烈的情况，球坐标系下界面偏离基准面的球坐标下界面反演深度加权函数：

其中η为深度加权的收缩因子，值为1或2，H为球坐标下相对于观测面的深度，r₀为参考界面深度半径,r为实际界面起伏半径，Δr为实际界面起伏半径相对于参考界面深度半径的差值，沿着半径方向为正值，其中r＝r₀+Δr；在经度和纬度方向相同尺度的单元体，其只存在半径方向的差异，单元体的体积可近似为球单元体各陵边的乘积，即其中λ、/>分别为经度和纬度方向单元体的变化量，公式可简化为

当Δr为0是，实际界面与参考界面重合时，w(Δr)＝1；当Δr＞0时，实际界面位于参考界面上方，此时w(Δr)＜1，对浅部的界面隆起进行压制；当Δr＜0时，实际界面位于参考界面下方，此时w(Δr)＞1，对深部的凹陷进行突出。

4.根据权利要求3所述的一种球坐标系下重震联合反演密度界面分布的方法，其特征在于，步骤(2)中所述地震数据约束方法为：

假定界面模型的差异和重力异常的差异存在线性的关系：

其中为实际界面和实际计算模型之间的修正值，Δg^M×N为理论重力异常和计算重力异常的残差，a和b为线性拟合的系数；

假定区域内共有n个地震测点作为深度约束点，通过最小二乘法线性回归可知：

其中和/>分别为第i个地震约束点的观测重力异常和计算重力异常，/>和/>分别为已知的地震深度点和反演计算的深度点，其中i∈n。