CN109471190B - 一种不同高度重力数据和井中重力数据联合反演方法 - Google Patents

Abstract

一种不同高度重力数据和井中重力数据联合反演的方法，由于不同高度观测面的重力数据包含的波场信息不同，该方法通过延拓技术得到不同高度的重力数据，将延拓数据与观测数据联合反演改善反演结果的多解性问题。当研究区只有地面观测重力数据时，下延数据与地面测量数据联合反演能够提高反演结果分辨率；当研究区只有航空重力数据时，将航空重力数据下延到地面或地面以下合适平面，下延数据与航空重力数据联合反演也能够提高反演结果分辨率；反演过程中加入井中重力、密度数据做约束，能够进一步提高反演结果的垂直分辨率。本发明改进的深度加权函数具有更强的针对性和灵活性，应用在以上联合反演方法中更容易获得理想的反演效果。

Description

一种不同高度重力数据和井中重力数据联合反演方法

技术领域

本发明涉及一种位场数据反演方法，特别涉及一种不同高度重力数据和井中重力数据联合反演方法。

背景技术

重力数据反演面临严重的多解性问题，高斯定理指出其本质原因：有无数种不同的密度分布可以引起相同的重力异常；此外，异常体引起的重力异常是连续分布的，而我们采集的信息只是某一平面上的部分离散数据，这种信息的不足进一步加重了反演结果的多解性。利用已知的先验地质信息对密度进行合理约束是降低反演多解性的有效手段，常用的约束方法有最小质量约束、最小体积约束、最大光滑约束、最小惯量约束等；另外，增加观测数据信息也能够降低反演结果多解性，增加的观测数据可以是同种地球物理数据加密或不同空间采集，也可以是不同类地球物理数据，但实际工作中，由于经费或者测量条件限制，无法获得更多观测数据，通过延拓技术获取更多可靠信息加入到反演中是可行的。

垂直分辨率差是重力数据反演面临的另一严重问题，本质原因与多解性问题相同，重力数据的核函数随深度的增加逐渐衰减是密度反演结果分辨率差的直接原因，通过深度加权抵消核函数衰减速率，能够有效克服反演结果的趋肤效应，提高垂直分辨率，因此，深度加权函数的选取对反演结果的分布非常重要。深度加权函数中不同的参数选择能够带来不同程度的深度加权效果，根据先验信息交互选择适当的深度加权函数能够得到合适的密度分布。

不同高度观测面的重力数据与地面观测重力数据相比，包含的波场信息不同，为了充分挖掘观测数据中的信息，本发明通过向上延拓和向下延拓技术获得不同高度的观测重力数据，并将其与原始测量重力数据同时反演，以期通过增加观测数据信息降低反演结果的多解性。由于向上延拓可压制高波数成分，即突出低波数成分；向下延拓可放大高波数成分，即突出高波数成分，一般情况下，高波数成分包含更多的浅部信息，低波数成分包含更多的深部信息，故理论上，上延数据、下延数据与原始数据联合反演能够在一定程度上提高反演效果。另外，当存在不同高度测量数据时，可由实际观测数据代替延拓数据，直接联合不同高度观测面的实际数据进行反演，此时与联合不同高度延拓面数据反演的方法完全相同。同时，本发明引入更灵活的深度加权函数，通过交互式逐步调整修正，可有效提高垂直分辨率。我们也在反演过程中加入井中重力数据、井中密度数据做约束，进一步提高反演结果的垂直分辨率。最后，使用共轭梯度法求解目标函数最优解，并通过理论数据验证本发明技术的可行性。

发明内容

本发明的目的是提供一种不同高度重力数据和井中重力数据联合反演方法。

重力数据的密度反演问题是典型的线性问题：

g＝Aρ (1)

其中，M维向量ρ是为密度向量；N维向量g是观测重力异常数据；N×M维矩阵A是正演模型算子，反演过程中已知观测重力异常g和正演算子A，求解密度分布ρ。由于观测数据中不可避免的噪声干扰，使得反演问题病态，问题的解不稳定，吉洪诺夫正则化方法作为一种常规的处理手段，可有效的解决这一问题。利用该理论，求解(1)式所示的反演问题可转化为如下目标函数：

φ＝φ_g+αφ_ρ (2)

其中φ_g是数据拟合函数，φ_ρ是模型目标函数，α代表正则化参数。数据拟合函数定义为：

φ_g＝||g-Aρ||²＝(g-Aρ)^T(g-Aρ) (3)

模型目标函数为最小模型约束：

φ_ρ＝||ρ||²＝ρ^Tρ (4)

通常，在目标函数中加入深度加权函数W，可有效克服重力反演结果中的趋肤问题。反演结果中异常体的位置很大程度上受深度加权函数影响，较弱的深度加权函数容易获得较浅的异常体埋深，而较强的深度加权则可以更容易获得深部密度变化特征。交互调节深度加权函数的强度，使反演的密度参数分布在更接近真实情况的深度范围内。

引入深度加权函数后，目标函数可以展开写成：

φ＝(g-Aρ)^T(g-Aρ)+α(Wρ)^T(Wρ) (5)

当仅有地面重力数据时，通过上延技术和下延技术能够得到不同高度面(高度或深度为z)的观测重力数据g^z：上角标z>0表示地面以下观测面；z<0表示地面以上观测面。将延拓数据g^z作为观测数据加入到目标函数的数据拟合函数中

φ_g＝λ₁||g-Aρ||²+λ₂||g^z-A^zρ||² (6)

A^z是与延拓数据g^z对应的正演算子。此外，还可以在数据拟合项中加入井中重力数据g^b、井中岩石密度数据ρ^d，得到目标函数如下：

φ_g＝λ₁||g-Aρ||²+λ₂||g^z-A^zρ||²+λ_b||g^b-A^bρ||+λ_d||ρ^d-A^dρ|| (7)

其中，g^b代表井中重力数据，A^b是与井中重力数据g^b对应的正演算子；ρ^d代表井中岩石剩余密度数据，A^d可以看成是与井中岩石密度ρ^d对应的正演算子，它是一个稀疏矩阵，矩阵每行仅有一个与取样岩石位置有关的元素为1，其他元素均为0。通过公式(7)实现了不同数据的联合，通过调整系数λ₁、λ₂、λ_b和λ_d的大小能够调整不同高度观测面重力数据、井中重力数据、井中岩石密度数据对目标函数的作用大小：系数越小，此类数据对目标函数影响越小；系数越大，此类数据对目标函数影响越大；当系数为0时，说明某一类数据对目标函数不起作用，即此类数据不参与反演过程。

反演过程中，较常用的深度加权函数如下

W＝(A^TA)^β/4 (8)

该深度加权函数中A为地面观测时模型正演算子，β为权重参数。常规反演中仅通过调节参数β来调节深度加权函数。在此基础上，通过加入其他观测面正演算子对其进行改进，改进后深度加权函数如下所示：

上式中，A_r表示改进后的重力数据正演算子，Az表示不同观测面的正演算子：上角标z>0表示地面以下观测面；z<0表示地面以上观测面。使用时可根据需求选择A、A^z或两者任意组合。通过选择不同观测面正演算子A^z和权重系数β，可以获得适用性更强、更灵活和更有针对性的深度加权函数。调整幂次β时，β越大，深度方向的加权强度越大，β越小，深度方向的加权强度越小，β＝0时，W为单位矩阵，没有深度加权效果；选择不同高度观测面的正演算子A_r时，观测面越高，加权强度越小，观测面越低，加权强度越大。当不同高度的观测数据同时反演时，我们可以根据需要选择深度加权函数中的正演算子A_r。尤其是当异常埋深较深时，我们可以同时交互调节深度加权函数中的参数A_r和β，以使密度分布在合适的深度范围内。

附图说明

图1是深度加权函数分布的中心竖直切片(a)A_r＝A，β＝1；(b)A_r＝A，β＝2；(c)A_r＝A^-200，β＝2；(d)A_r＝A²⁰⁰，β＝2；(e)β＝2；(f)β＝2。

图2(a)正演模型中心切片及(b)模型正演地面重力数据。

图3地面重力数据单独反演结果切片(a)未加深度加权；(b)加深度加权。

图4延拓重力数据单独反演结果，加深度加权(a)上延100m；(b)下延100m；(c)下延200m；(d)下延300m。

图5地面重力数据与下延重力数据联合反演结果(a)与下延100m重力数据；(b)与下延200m重力数据；(c)与下延300m重力数据；(d)与下延100m、200m、300m重力数据共同联合。

图6(a)模拟航空测量重力数据，飞行高度为300m；(b)模拟航空重力数据单独反演的密度分布切片。

图7(a)模拟航空重力数据与下延300m数据联合反演结果切片(b)模拟航空重力数据与下延400m数据联合反演结果切片。

图8(a)双长方体模型正演地面重力数据，含3％高斯噪声；(b)由(a)数据下延300m的重力数据。

图9是模拟井中测量数据图。(a)井A观测重力异常数据；(b)井A岩石剩余密度数据；(c)井B观测重力异常数据；(d)井B岩石剩余密度数据。含3％高斯噪声。

图10是双立方体模型反演结果。(a)地面重力数据结果，无深度加权；(b)地面重力数据结果，加深度加权；(c)地面重力数据和下延300m重力数据联合反演结果，加深度加权；(d)井中重力数据、井中岩石剩余密度数据约束下，地面重力数据和下延300m重力数据联合反演结果，加深度加权。

具体实施方式

首先，为了验证不同参数组合深度加权函数的效果，我们将地下三维空间剖分为38×42×15个大小为50m×50m×100m的长方体单元，42×38个观测点位于每个密度单元的中心，基于这个观测系统，我们调节深度加权函数的参数，对应的深度加权函数的分布如图1所示，归一化展示是为了方便比较不同参数对应的加权函数的强度。当A_r＝A，β＝1和β＝2时，对应的深度加权函数分布中心垂直切片如图1(a)和图1(b)，两者都能够有效克服反演结果的趋肤效应，β＝2对应的加权函数随深度衰减更快，加权强度更大，更容易得到较深的反演结果；当A_r＝A^-200，β＝2时，其深度加权函数分布如图1(c)，与图1(b)相比，随深度衰减速度慢，说明其加权强度小，容易获得较浅反演结果；当A_r＝A²⁰⁰，β＝2时，其深度加权函数分布如图1(d)，其随深度衰减速度快，得到的反演结果远离z＝200m平面；当β＝2时，其深度加权函数分布如图1(e)，地下0-200m部分几乎不衰减，200以下迅速衰减，此函数加入到反演中，反演结果容易在200m以下分布；当β＝2时，其深度加权函数分布如图1(f)，地下0-400m部分几乎不衰减，400m以下迅速衰减，此加权函数稳定，加入到反演中，反演结果将在400m以下分布，故容易得到相对较深的反演结果。综上，反演结果容易在所选正演算子所在观测面以下分布，即加权函数中的参数A_r所在观测面不宜低于异常体顶面。当已知密度异常体埋深的先验信息时，对加权函数的选择有一定意义，当埋深较浅时，β可取1，A_r可取地面重力数据与地面以上重力数据正演算子的组合；当埋深较深时，β可取2，A_r可取地面重力数据与地面以下重力数据正演算子的组合；选择几个不同深度观测面的重力数据正演算子与地面重力数据正演算子组合能够得到针对性更强、灵活性更好的深度加权函数。

为了验证不同类型数据约束对地面数据反演的改善效果，建立一个长方体的正演模型。将地下三维空间剖分为38×42×15个大小为50m×50m×100m的长方体单元，42×38个观测点位于每个密度单元的中心，背景密度为0，长方体异常大小为600m×1000m×400m，密度差为1g/cm³，顶面埋深为300m，中心位置为(950m，1050m，500m)，模型x方向中心竖直切片如图2(a)所示，其地面观测重力异常数据如图2(b)所示。

常规光滑反演方法处理地面数据，引入密度上下限约束0-1g/cm³，不加深度加权反演结果如图3(a)，密度分布在地表附近，趋肤效应严重。处理地面重力数据反演地下密度分布时，交互调节深度加权函数的参数A_r和β，最后选择其中，A表示地面重力数据的正演算子，A¹⁰⁰表示地面以下100m处观测重力数据的正演算子，反演结果如图3(b)所示，与图3(a)中结果对比，加入合适的深度加权函数后，异常体位置更准确，但密度异常分布仍然比较发散，且密度值远远小于理论密度值。这说明，引入深度加权函数能够在一定程度上改善反演效果，但地面重力数据单独反演效果并不理想。

为了进一步开发地面重力数据中的信息，通过FFT上延技术得到不同高度观测面的重力数据，通过迭代下延技术得到不同深度观测面的重力数据。对比延拓数据与理论数据，FFT上延数据与理论数据的残差标准差较小，吻合较好；对于下延数据，当下延数据观测面在异常体上方，即观测面未穿过异常体时，下延数据与理论数据的残差标准差也较小，吻合较好；当下延数据的观测面穿过异常体后，下延数据与理论数据的差异随深度增加迅速变大。因此，穿过异常体及异常体以下观测面的延拓数据与实际模型信息差异较大，不再适合作为观测数据加入到反演程序。基于以上分析，我们尝试将上延数据以及异常体以上观测面的下延数据融合到反演中，以期改善反演效果。首先分别验证单独上延数据和单独下延数据反演的效果，深度加权函数不变，上延100m重力数据反演结果如图4(a)所示，异常体中心位置识别准确，与地面数据反演结果相比，密度分布发散更严重，异常体密度值最大不超过0.1g/cm³，仅为理论密度值得1/10，说明单独的上延数据反演结果分辨率低。下延100m重力数据反演结果如图4(b)，与地面数据反演结果相比，分辨率明显提高，异常体中心位置识别准确，边界较为清晰，异常体的最大密度值达到0.5g/cm^-3，大于理论模型密度值的一半。观察到了下延100m数据的良好反演效果后，为了进一步验证下延数据反演的规律和优势，我们又分别单独反演了下延200m数据和下延300m数据，结果如图4(c)、4(d)，与图5(b)共同说明，在下延数据合理的前提下，下延深度越深，反演结果的分辨率越高，反演的密度值越接近理论模型密度。

由于下延重力数据单独反演分辨率较高，我们将下延重力数据与原始地面观测重力数据联合反演地下密度分布，希望在拟合原始观测重力数据的前提下，提高反演结果分辨率。将地面重力数据分别与下延100m、200m、300m重力数据联合反演，目标函数中的数据拟合项分别为φ_g＝λ₁||g-Aρ||²+λ₂||g¹⁰⁰-A¹⁰⁰ρ||²、φ_g＝λ₁||g-Aρ||²+λ₂||g²⁰⁰-A²⁰⁰ρ||²、φ_g＝λ₁||g-Aρ||²+λ₂||g³⁰⁰-A³⁰⁰ρ||²，其中λ₁＝λ₂＝1，深度加权函数仍然不变，反演结果如图5(a)、5(b)、5(c)所示，参与联合反演的下延数据的下延深度越深，异常体边界越清晰，反演结果的分辨率越高，反演的密度值越接近理论异常值，这与下延数据单独反演的规律类似。另外，我们也计算了同时联合地面重力数据、下延100m、200m、300m重力数据等四个平面数据的反演，此时，目标函数中的数据拟合项为φ_g＝λ₁||g-Aρ||²+λ₂||g¹⁰⁰-A¹⁰⁰ρ||²+λ₃||g²⁰⁰-A²⁰⁰ρ||²+λ₄||g³⁰⁰-A³⁰⁰ρ||²，其中，λ₁＝λ₂＝λ₃＝λ₄＝1，反演结果如图5(d)，异常体分辨率也较高，但是，与图5(c)相比，效果提高并不明显，这说明，参与反演的延拓面越多，计算量越大，但反演效果不一定越好，选择一个合适的平面的下延数据与地面观测数据联合反演，可以在保证计算效率、拟合原始观测数据的前提下，提高单独地面重力数据反演的效果。

基于以上分析结果，我们进一步验证航空重力数据(测量平面在地面以上)与下延数据联合反演对比航空重力数据单独反演，是否也能改善反演效果。同样，基于以上密度异常体模型以及测量网格，我们假设仅存在航空测量重力数据，飞行测量高度为地面以上300m，采集到的重力数据如图6(a)，引入深度加权函数单独反演结果如图6(b)，可以看到，航空重力数据单独反演效果很差，难以收敛。将原始航空重力异常下延到地面，即下延高度为300m，得到的延拓重力数据与航空测量数据联合反演地下密度信息，引入深度加权函数得到的结果如图7(a)，能够看到，反演结果收敛，分辨率大大提高；保持深度加权函数不变，航空重力数据与下延400m数据联合反演，反演结果如图7(b)所示，密度分布进一步收敛，分辨率进一步提高。所以，航空重力数据与下延数据联合反演对比航空重力数据单独反演，效果明显提高。

为进一步验证下延数据与地面数据联合反演的适用性以及井数据的加入对反演效果的影响，本节设置不同埋深的双立方体模型，地下剖分网格和观测数据网格都与模型一相同，两个长方体密度异常与模型一中长方体形态、大小、密度差都相同，仅位置不同，左侧长方体顶面埋深300m，中心位置坐标为(550m，1050m，500m)，右侧长方体顶面埋深为500m，中心位置坐标为(1350m，1050m，700m)。在正演计算的理论数据中加入3％的高斯噪声，如图8(a)，将带噪声的地面数据下延300m得到的重力数据如图8(b)所示。建立两口模拟竖直井记为井A、井B，水平位置分别为(550m，1050m)、(1350m，1050m)，即位于两个密度异常体中心，井中模拟理论重力数据和岩石剩余密度数据采样间距为50m，同样填加3％高斯噪声，如图9所示。不引入深度加权的地面重力数据反演结果如图10(a)，趋肤效应明显，且不能明显区分两个密度异常体，引入深度加权函数即深度加权函数中A_r为地面重力数据正演算子和地下100m、200m重力数据正演算子的组合，引入深度加权后，地面重力数据单独反演结果如图10(b)，埋深相对较浅的异常体中心位置识别准确，埋深相对较深的异常体在反演结果的图像中呈现不明显，且无法区分两个异常体边界，反演结果整体发散严重，分辨率较低；深度加权函数不变，引入下延300m重力数据(图8(b))与地面数据联合反演，数据拟合项为φ_g＝λ₁||g-Aρ||²+λ₂||g³⁰⁰-A³⁰⁰ρ||²，λ₁＝λ₂＝1，结果如图10(c)，与图10(b)相比，埋深较浅的异常体不仅中心位置定位准确，而且异常体边界清晰，密度值的最大值达到理论异常值，埋深较浅的异常体反演效果并不理想，另外，由于下延300m数据中噪声很大，导致反演结果中也出现一些假异常；再引入两口井的井中重力、岩石剩余密度数据(图9)共同约束地面重力数据和下延300重力数据联合反演，数据拟合项为φ_g＝λ₁||g-Aρ||²+λ₂||g³⁰⁰-A³⁰⁰ρ||²+λ_b||g^b-A^bρ||+λ_d||ρ^d-A^dρ||，其中，λ₁＝λ₂＝λ_h＝λ_d＝1，结果如图10(d)，两个不同埋深的异常体中心位置都得到了准确的定位，边界清晰，分辨率高，两个异常体的最大异常值都达到了理论异常值，效果令人满意。

Claims (1)

1.一种不同高度重力数据和井中重力数据联合反演方法，该方法是：不同高度观测面的重力数据包含的波场信息不同，通过延拓技术得到不同高度的延拓重力数据，充分挖掘观测重力数据中的信息，并将延拓数据与观测数据联合反演改善反演多解性问题，并提高反演结果分辨率；在联合反演中引入改进的深度加权函数，改进的深度加权函数能够使反演结果更容易分布在合适的深度范围内；其特征在于：改进的深度加权函数表达式为：

上式中，A是正演模型算子，A_r表示改进后的重力数据正演算子，A^z表示不同观测面的正演算子：上角标z>0表示地面以下观测面；z<0表示地面以上观测面；使用时能根据需求选择A、A^z或两者任意组合；通过选择不同观测面正演算子A^z和权重系数β；调整幂次β时，β越大，深度方向的加权强度越大，β越小，深度方向的加权强度越小，β＝0时，W为单位矩阵，没有深度加权效果；选择不同高度观测面的正演算子A_r时，观测面越高，加权强度越小，观测面越低，加权强度越大；当不同高度的观测数据同时反演时，根据需要选择深度加权函数中的正演算子A_r；当异常埋深较深时，同时交互调节深度加权函数中的参数A_r和β，以使密度分布在合适的深度范围内。