CN112147709A - 一种基于部分光滑约束的重力梯度数据三维反演方法 - Google Patents

Abstract

本发明涉及一种基于部分光滑约束的重力梯度数据三维反演方法，首先计算长方体单元引力位在三维空间x，y，z三个方向上的二阶偏导数得到长方体元的正演公式，给定模型后根据正演公式正演计算得到观测数据；再依据一定的先验信息对粗糙度约束矩阵进行人为赋值，得到粗糙度矩阵，并采用积分灵敏度的加权方式克服“趋肤效应”；最后写出含有积分灵敏度矩阵，粗糙度矩阵，正则化参数，观测数据等参数的目标方程，并采用共轭梯度法解目标方程，得到反演结果。该技术适用于在有一定地质或地球物理先验条件下的三维重力和重力梯度数据反演。

Description

一种基于部分光滑约束的重力梯度数据三维反演方法

技术领域

本发明涉及地球物理勘探领域中三维重力梯度数据反演问题，具体说，涉及一种基于部分光滑约束的重力梯度数据三维反演方法。

背景技术

如何快速准确的确定地质体位置和物性是找矿地质勘探中最重要的问题之一，而重力或重力梯度的物性反演技术是解决该问题的关键技术环节。在实际工作中，因为 采集的有效数据有限，测量精度受各种因素干扰大，特别是场源的等效性使得三维物 性反演存在严重的多解性，故反演问题一般是欠定问题。解决这种问题目前最常用的 办法就是加入地质和其他地球物理方法提供的先验信息或者合理适当的约束条件来处 理。Li和Oldenburg(1998)提出的粗糙度矩阵给出了一种反演的光滑模型，并通过 增加深度加权矩阵压制了“趋肤效应”，Li(2000)进一步将光滑反演方法推广到重力 梯度全张量数据的反演中，Martinez和Li(2013)采用该方法对航空重力梯度张量数 据的不同组合进行反演，证明了全张量数据联合反演结果能够更好的反演地质体结构。 Silva等(2001)讨论了现有物性反演方法中的约束条件，认为应根据实际的地质问题 选择约束条件。光滑反演一般采取的是全局光滑，得到的反演结果边界模糊且范围扩 大，模型分辨率较低。

发明内容

本发明所要解决的问题是提高反演结果的分辨率，进而较真实反映地下地质体的位置和物性。

本发明提供一种基于部分光滑约束的重力梯度数据三维反演方法，包括如下步骤：

1)分别计算引力位在三维空间x，y，z三个方向上的二阶偏导数，得到重力梯 度正演公式，再利用上述公式，计算给定模型正演重力梯度值作为反演的观测数据。

2)依据先验信息对粗糙度约束矩阵进行人为赋值，得到粗糙度矩阵；

3)采用积分灵敏度的加权方式形成积分灵敏度矩阵。

4)基于步骤1)中观测数据，步骤2)中的粗糙度矩阵，步骤3)中的积分灵敏 度矩阵，选取正则化参数，得到目标方程，

5)采用共轭梯度法解目标方程，得到反演结果。

所述步骤1)具体为：以(ξ,η,ζ)为长方体内体积微元的坐标，长方体单元坐标值最小和最大的两个点的坐标分别为(ξ₁,η₁,ζ₁)，(ξ₂,η₂,ζ₂)，点(x,y,z)为观测点，舒晴 等人推导得到的无解析奇点的重力梯度全张量正演计算公式如下所示：

其中，G为万有引力常量，r＝[(x-ξ)²+(y-η)²+(z-ζ)²]^1/2，ρ为长方体的剩余 密度，s＝(-1)^i+j+k。

所述步骤1)还包括：建立一个三维地质理论模型，将地下空间剖分为J行K列 L层整齐排列的长方体单元，对于不同的正演模型，各长方体单元除密度值外完全相 同，以靠近地面那一层每个单元的顶部中心位置为观测点，根据正演公式计算观测点 重力梯度值。

所述步骤2)中，因为重力及其梯度数据对于异常体深度信息不敏感，故一般需 要其它地质或地球物理方法提供异常体的深度信息，采用深度信息作为先验信息进行 约束。同时，异常体倾向也可以作为先验信息添加到粗糙度矩阵中。

上述先验信息可以由其它地质或地球物理信息提供，如异常体埋深，异常体倾向等，先验信息越多则反演结果越准确。

所述步骤2)中的粗糙度矩阵为：

式中M为光滑约束的条件数，N为地下模型单元的个数。

所述步骤3)具体为：观测数据d与模型m变化的关系可表示为

δd_i＝F_ikδm_k

其中，F_ik为观测数据关于参数的灵敏度(导数)。数据灵敏度对模型m_k的积分 可表示为

故积分灵敏度矩阵S为对角矩阵；

S＝diag(F^TF)^1/2

其中，F为观测数据关于参数的灵敏度，F^T为F的转置，diag表示取对角矩阵。

所述步骤4)中的目标方程为：φ(m)＝(d-Gm)^T(d-Gm)+αm^TW_m ^TW_mm；式中，α为 正则化参数，W_m＝RS，S为积分灵敏度对角矩阵，R为粗糙度矩阵，m为要求解的模 型单元密度值。

所述步骤5)中，目标函数可以写成以下增广矩阵的形式：

采用共轭梯度法解上述目标函数，得到矩阵[ρ]的解。

本发明的有益效果在于：相对于现有技术中针对当进行全局光滑约束时，反演结果模型边界模糊且范围扩大，模型分辨率较低，本发明首次提出了根据已知先验信息 控制光滑范围的部分光滑约束矩阵，将部分光滑约束矩阵加入到重力梯度异常数据三 维反演中。与现有技术相比，能够显著提高反演结果的分辨率，提高了模型反演结果 的准确性。

附图说明

图1三维重力梯度反演方法流程图；

图2模型一模型示意图；

图3模型一切片图

图4模型一全局光滑反演结果图；

图5模型一部分光滑反演结果图；

图6模型二切片图；

图7模型二部分光滑反演结果图。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

本发明涉及一种基于部分光滑约束的重力梯度数据三维反演方法，包括以下步骤：

a、建立笛卡尔直角坐标系，具体为：以(ξ,η,ζ)为长方体内体积微元的坐标，长 方体单元坐标值最小和最大的两个点的坐标分别为(ξ₁,η₁,ζ₁)，(ξ₂,η₂,ζ₂)，以点(x,y,z) 为观测点，则一个剩余密度为ρ长方体的重力梯度正演公式如下：

其中，G为万有引力常量，r＝[(x-ξ)²+(y-η)²+(z-ζ)²]¹²,s＝(-1)^i+j+k。

b、输入观测数据d，这个观测数据由理论模型给出。将地下空间剖分为一系列21×21×10 的长方体单位元，其中每个单位元的大小为10m×10m×20m，观测点为最上层每个单位 元的顶部中心位置。模型一是一个地下存在两个异常体的三维地质模型。地质体位置如图2所示。其中左侧地质体剩余密度为0.8g/cm³，右侧地质体剩余密度为1g/cm³， 做三维切片图如图3所示，正演其重力梯度数据U_zz，得到观测数据。模型二为台阶模 型，剩余密度为0.8g/cm³，模型切片图如图6所示，正演其U_ZZ作为观测数据。

c、粗糙度矩阵按照以下示例的形式给出。

式中M为光滑约束的条件数，N为地下模型单元的个数。如上述矩阵的第一行 前两个元素分别为-1和1，该行其余元素全部为零，则表示第一个模型单元和第二 个模型单元之间进行了光滑约束，每多一条约束，则矩阵增加一行，最终形成粗糙 度矩阵。粗糙度矩阵的具体形式由其他地质地球物理信息给出，如在本示例模型一 中已知该异常体的埋深，模型二中除了埋深还把异常体倾向作为先验信息。

d、由于各个剖分单元与观测网格点之间几何关系，导致了距离观测点越远的 剖分单元的重力梯度张量导数越小，如果不对模型进行加权约束就会使反演结果出 现“趋肤效应”，即异常都集中在地表很难反演出真实构造。因此，本文引入积分灵 敏度矩阵S，对模型进行加权，矩阵S由于具有x,y,z三方向加权效果，使得反演 具有三方向的分辨能力，因此，反演算法的旁侧反演能力得到大大提高。

观测数据与模型变化的关系可表示为

δd_i＝F_ikδm_k

其中，F_ik为观测数据关于参数的灵敏度(导数)。数据灵敏度对模型m_k的积分 可表示为

故积分灵敏度矩阵S为对角矩阵；

S＝diag(F^TF)^1/2

其中，F为观测数据关于参数的灵敏度，F^T为F的转置，diag表示取对角矩阵。

e、写出含有约束矩阵W_m，正则化参数，观测数据等参数的目标方程,即 φ(m)＝(d-Gm)^T(d-Gm)+αm^TW_m ^TW_mm。式中，α为正则化参数，W_m＝RS，S为积分灵 敏度对角矩阵，R为粗糙度矩阵。本文正则化参数选取的方法是一种易于编程的方法， 即选用α的初值为一个较大的常数值，之后每次迭代过程中α总为之前的一半，从而 实现需要光滑约束的模型空间尽可能光滑，且数据拟合误差在可接受范围内。

f、则目标函数可以写成以下增广矩阵的形式：

采用共轭梯度法解上述目标函数，得到矩阵[ρ]的解。

按照上述过程模型一反演结果如图4、图5所示，全局光滑反演结果如图4所示， 部分光滑时，采用左半侧空间上五层仅z方向光滑，右侧空间上四层仅z方向光滑， 使该空间范围内单元体之间x,y方向无约束，其余空间位置仍然采用x,y,z三方向的全 局光滑，反演结果如图6所示。其中红色线框表示地质体实际位置。

对比图4和图5，两种光滑方式均可以确定地质体的水平位置，但部分光滑反演 得到的地质体边界位置更加清楚，这是因为在全局光滑中地质体与围岩之间存在光滑 约束，使地质体边界模糊。

对比图4和图5，部分光滑相比于全局光滑有明显的优势，在图4中右侧异常体 纵向位置不明显，而图5中右侧地质体纵向位置可以准确的的反演得到，同时对于左 侧的异常体，图5边界更加清晰且反演得到的剩余密度值更接近于所给定模型的剩余 密度值。

对于模型二，采用部分光滑约束为上五层没有x,y方向的光滑，若已知地下模型可能为倾斜构造及构造方向，则添加右倾对角单元之间的光滑关系，下五层仍然采用 x,y,z三方向的全局光滑，反演结果如图7所示。可以看出，在浅部可以很好的得到地 质体的水平位置及纵向位置，在较深部区域水平位置略有偏差，则部分光滑反演基本 可以得到台阶模型的轮廓位置。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，包括如下步骤：

1)分别计算引力位在三维空间x，y，z三个方向上的二阶偏导数，得到重力梯度正演公式，再利用上述公式，计算给定模型正演重力梯度值作为反演的观测数据。

2)依据先验信息对粗糙度约束矩阵进行人为赋值，得到粗糙度矩阵；

3)采用积分灵敏度的加权方式形成积分灵敏度矩阵；

4)基于步骤1)中观测数据，步骤2)中的粗糙度矩阵，步骤3)中的积分灵敏度矩阵，选取正则化参数，得到目标方程，

5)采用共轭梯度法解目标方程，得到反演结果。

2.根据权利要求1所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤1)具体为：以(ξ,η,ζ)为长方体内体积微元的坐标，长方体单元坐标值最小和最大的两个点的坐标分别为(ξ₁,η₁,ζ₁)，(ξ₂,η₂,ζ₂)，以点(x,y,z)为观测点，则一个剩余密度为ρ长方体的重力梯度正演公式如下：

其中，G为万有引力常量，r＝[(x-ξ)²+(y-η)²+(z-ζ)²]^1/2,s＝(-1)^i+j+k。

3.根据权利要求2所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤1)还包括：建立一个三维地质理论模型，将地下空间剖分为J行K列L层整齐排列的长方体单元，对于不同的正演模型，各长方体单元除密度值外完全相同，以靠近地面那一层每个单元的顶部中心位置为观测点，根据正演公式计算观测点重力梯度值。

4.根据权利要求3所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤2)中，所述先验信息为通过地质或地球物理方法提供异常体的深度信息或者异常体倾向。

5.根据权利要求3所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤2)中的粗糙度矩阵为：

式中M为光滑约束的条件数，N为地下模型单元的个数。

6.根据权利要求5所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤3)具体为：观测数据d与模型m变化的关系可表示为

δd_i＝F_ikδm_k

其中，F_ik为观测数据关于参数的灵敏度；数据灵敏度对模型m_k的积分可表示为

故积分灵敏度矩阵S为对角矩阵；

S＝diag(F^TF)^1/2

其中，F为观测数据关于参数的灵敏度，F^T为F的转置，diag表示取对角矩阵。

7.根据权利要求6所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤4)中的目标方程为：φ(m)＝(d-Gm)^T(d-Gm)+αm^TW_m ^TW_mm；式中，α为正则化参数，W_m＝RS，S为积分灵敏度对角矩阵，R为粗糙度矩阵，m为要求解的模型单元密度值。

8.根据权利要求6所述的基于部分光滑约束的重力梯度数据三维反演方法，其特征在于，所述步骤5)中，目标函数可以写成以下增广矩阵的形式：

采用共轭梯度法解上述目标函数，得到矩阵[ρ]的解。

Images