CN112114374B - 复杂地质体的多密度界面反演方法 - Google Patents

Abstract

本发明公开了一种复杂地质体的多密度界面反演方法，步骤为：构建与目标复杂地质体相符的密度反演模型；将地表处实测布格重力异常g代入基于位场分离目的的切割法计算公式中，以求出各密度层的底界面在地面产生的重力异常；应用位场大深度向下延拓方法将目标复杂地质体的各密度层的底界面在地面产生的重力异常延拓至基准面；将各密度层底界面在其基准面引起的重力异常依次进行界面反演，得到各界面相对基准面的深度，进而得到各密度层底界面的实际深度；该反演方法能够有效获得复杂地质体密度层实际分布状况，提高反演结果的准确率，相对于Parker‑Oldenburg反演方法，其反演精度可提升1％～6％。

Description

复杂地质体的多密度界面反演方法

技术领域

本发明涉及地球物理学领域技术领域，特别涉及一种复杂地质体的多密度界面反演方法。

背景技术

Parker-Oldenburg反演方法是Oldenburg在Parker公式的基础上提出的频率域密度界面迭代反演方法。由于该方法具有计算快速的优点，其应用得到了快速发展。但公式中存在一个向下延拓因子，使得公式在反演迭代过程中存在高频信号的震荡现象，且随着反演深度的增大，这种现象越明显，严重影响了反演结果的收敛性。为解决这一情况，通常在反演过程应用低通滤波器消除高频震荡现象，这样，虽然保证了反演结果能够收敛，但反演精度大大降低，且低通滤波器的参数设置也难以控制。

在设计密度反演模型时，为提高反演效果，许多学者采用变密度的密度反演模型，如：双曲线模型、抛物线模型、指数模型等。这些变密度的密度反演模型针对不同的地质环境很大程度上改善了反演效果，提高了精度，但主要是针对单一密度模型，当模型较复杂时，即：不能用一个函数表达密度与深度的关系，上述模型往往不能满足反演要求。

为实现复杂地质体的多密度界面反演，本文首先设计了多界面复杂密度体模型，模型中认为各密度层分界面起伏是引起地表重力异常的主要原因，并将各密度层顶界面平均深度处设为各密度层底界面的基准面；其次，应用切割法从布格重力异常中分离出各密度层底界面在地面引起的重力异常；然后，应用大深度向下延拓方法将各异常延拓至各基准面；最后应用Parker-Oldenburg反演方法实现各底界面的精确反演。完成复杂地质体的密度模型设计及提高Parker-Oldenburg反演方法的收敛性。

发明内容

本发明的目的是提供一种解决目前基于Parker-Oldenburg反演方法获得的复杂地质体的多密度层分界面的反演计算结果偏差较大的问题的复杂地质体的多密度界面反演方法。

为此，本发明技术方案如下：

一种复杂地质体的多密度界面反演方法，步骤如下：

S1、根据目标复杂地质体的地层密度分布资料，构建与目标复杂地质体相符的密度反演模型；

S2、将地表处实测的布格重力异常g代入基于位场分离目的的切割法计算公式中，计算出目标复杂地质体的各密度层的底界面及以上各密度层分界面在地面产生的重力异常，进而求出各密度层的底界面在地面产生的重力异常；

S3、应用位场大深度向下延拓方法将目标复杂地质体的各密度层的底界面在地面产生的重力异常延拓至基准面；其中，自地表向下，第一密度层的基准面为地表面，其他各密度层的基准面为其上层密度层底界面的平均深度处；

S4、将各密度层底界面在其基准面引起的重力异常依次通过Parker-Oldenburg反演方法实现界面反演，以得到各界面相对基准面的深度，进而求出各密度层底界面的实际深度。

进一步地，在步骤S1中，构建与复杂地质体相符的密度反演模型的具体方法为：

以地表为x-y平面，沿深度方向为z轴方向，从地表沿z轴方向依次构建第一密度层、第二密度层、第三密度层、…、第n密度层；其中，

对应每个密度层的密度依次为：ρ₁、ρ₂、ρ₃、…、ρ_n；

对应每个密度层与下方密度层的分界面为其底界面，各密度层的底界面依次为：第一密度层底界面、第二密度层底界面、第三密度层底界面、…、第n-1密度层底界面；

对应每个密度层与上方密度层的分界面为其顶界面，各密度层的顶界面依次为：第一密度层顶界面、第二密度层顶界面、第三密度层顶界面、…、第n密度层底界面；其中，第一密度层顶界面为地表平面；

基于各密度层分界面起伏是造成地表重力异常的主要原因，设定地表x-y平面为第一密度层的基准面，第一密度层底界面的初始平均深度处为第二密度层的基准面，第二密度层底界面的初始平均深度处为第三密度层的基准面、…、第n-1密度层底界面的初始平均深度处为第n密度层的基准面；

对应每个密度层的地层平均厚度，即密度层的基准面与下方密度层的基准面之间的距离，依次为h₁、h₂、h₃、…、h_n-1；各密度层的地层平均厚度通过估算方法获得；

进一步地，在步骤S2中，基于位场分离目的的切割法计算公式如下：

式(1)中，α是加权系数；g(x,y)是待位场分离点的重力异常；点第一次切割的区域场；/>是待位场分离点周围四点重力异常值的算术平均值；其中，/>为g(x,y)点第一次切割的区域场；

α＝β+γ 式(3)，

式(2)和(3)中，r为切割半径，其取值为对应密度层至地表的各密度层的地层平均厚度之和；g(x+r,y),g(x-r,y),g(x,y+r)和g(x,y-r)是半径为r的周围4个点的重力异常；β和γ分别根据x轴方向和y轴方向的半二阶差分量确定；

同理，重复上述公式(1)～(3)，将所得一次区域场作为第二次切割的原始异常，得到第二次切割的区域场/>如此迭代，直至满足：

式(4)中，为获取的第n次区域场，ε₀为允许误差；

则区域场：

进而，得到近似深度为r以上的各地层分界面产生的重力异常

L_loc＝g-R_reg 式(6)，

式(6)中，g为地表处实测的布格重力异常；

基于此，通过该密度层顶界面的平均深度r_i-1和该密度层底界面的平均深度r_i，求出目标复杂地质体的各密度层的底界面在地面产生的重力异常：

u_i＝L_{loc_i}-L_{loc_i-1} 式(7)，

式(7)中，i为该密度层的编号，L_{loc_i}为第i密度层底界面及以上地层底界面在地面产生的重力异常。

进一步地，在步骤S3中，其具体转换步骤如下：

S301、假设第i密度层底界面在其基准面产生初始重力异常：

u₁(x,y，0)＝u(x,y) 式(8)，

式(8)中，u(x,y)为该密度层底界面在地面产生的重力异常；

S302、将初始异常值运用傅里叶变换法做向上延拓，延拓至地面，公式如下：

U₁(k_x,k_y，0)＝DFT(u₁(x,y，0)) 式(9)，

式(9)和(10)中，DFT表示傅里叶变换，IDFT表示傅里叶逆变换，k_x为x方向上的波数，k_y为y方向上的波数；

S303、对初始值u₁(x,y,0)进行校正，公式如下：

u₂(x,y，0)＝u₁(x,y，0)-(u₁(x,y,r_i)-u(x,y)) 式(11)，

式(11)中，u₂(x,y,0)为校正后的重力异常；

S304、重复S301～S303，直至上延计算值与密度层在地面上的异常值的差值小到可以忽略的程度，即：

|u_n(x,y,r_i)-u(x,y)＜ε 式(12)，

式(12)中，u_n(x,y,r_i)为第n次向上延拓后的重力异常，ε为允许误差；所得u_n(x,y,0)即为u(x,y)向下延拓r_i-1的重力异常，即第i密度层底界面在其基准面产生的重力异常。

与现有技术相比，该复杂地质体的多密度界面反演方法能够有效获得复杂地质体密度层实际分布状况，有效解决了Parker-Oldenburg反演方法存在的收敛性差的缺陷，提高反演结果的准确率，相对于Parker-Oldenburg反演方法，其反演精度绝大部分提升1％～6％。

附图说明

图1为本发明的复杂地质体的多密度界面反演方法中针对目标复杂地质体构建的密度反演模型的结构示意图；

图2为本发明的复杂地质体的多密度界面反演方法中基于位场分离目的的切割法的算法示意图；

图3(a-1)为本发明的实施例1的复杂地质体模型的第一密度层底界面等深图；

图3(a-2)为本发明的实施例1的复杂地质体模型的第二密度层底界面等深图；

图3(a-3)为本发明的实施例1的复杂地质体模型的第三密度层底界面等深图；

图3(b)为本发明的实施例1的复杂地质体模型的各密度层底界面产生的重力总异常图；

图3(c-1)为本发明的实施例1的复杂地质体模型的第一密度层底界面在地面产生的重力异常示意图；

图3(c-2)为本发明的实施例1的复杂地质体模型的第二密度层底界面及以上各密度层分界面在地面产生的重力异常示意图；

图3(c-3)为本发明的实施例1的复杂地质体模型的第三密度层底界面及以上各密度层分界面在地面产生的重力异常示意图；

图3(d-1)为本发明的实施例1的复杂地质体模型的第一密度层底界面在地面产生的重力异常示意图；

图3(d-2)为本发明的实施例1的复杂地质体模型的第二密度层底界面在地面产生的重力异常示意图；

图3(d-3)为本发明的实施例1的复杂地质体模型的第三密度层底界面在地面产生的重力异常示意图；

图3(e-1)为本发明的实施例1的复杂地质体模型的第一密度层底界面在基准面产生的重力异常的示意图；

图3(e-2)为本发明的实施例1的复杂地质体模型的第二密度层底界面在基准面产生的重力异常的示意图；

图3(e-3)为本发明的实施例1的复杂地质体模型的第三密度层底界面在基准面产生的重力异常的示意图；

图3(f-1)为本发明的实施例1的复杂地质体模型的第一密度层底界面的反演等深示意图；

图3(f-2)为本发明的实施例1的复杂地质体模型的第二密度层底界面的反演等深示意图；

图3(f-3)为本发明的实施例1的复杂地质体模型的第三密度层底界面的反演等深示意图；

图4为本发明的实施例2的天津市团泊-大寺地区地质构造及其上各类钻孔位置的示意图

图5(a)为本发明的实施例2的研究区的第一密度层(第四系)的切割重力异常；

图5(b)为本发明的实施例2的研究区的第二密度层(新近系)的切割重力异常；

图5(c)为本发明的实施例2的研究区的第三密度层(古近系)的切割重力异常；

图5(d)为本发明的实施例2的研究区的第四密度层(中生界)的切割重力异常；

图5(e)为本发明的实施例2的研究区的第五密度层(古生界～蓟县系)的切割重力异常；

图5(f)为本发明的实施例2的研究区的第六密度层(长城系)的切割重力异常；

图6(a)为本发明的实施例2的研究区的第一密度层(第四系)在第一密度层基准面的重力异常；

图6(b)为本发明的实施例2的研究区的第二密度层(新近系)在第二密度层基准面的重力异常；

图6(c)为本发明的实施例2的研究区的第三密度层(古近系)在第三密度层基准面的重力异常；

图6(d)为本发明的实施例2的研究区的第四密度层(中生界)在第四密度层基准面的重力异常；

图6(e)为本发明的实施例2的研究区的第五密度层(古生界～蓟县系)在第五密度层基准面的重力异常；

图6(f)为本发明的实施例2的研究区的第六密度层(长城系)在第六密度层基准面的重力异常；

图7(a)为本发明的实施例2的研究区的第一密度层(第四系)在第一密度层底界面的密度等深图；

图7(b)为本发明的实施例2的研究区的第二密度层(新近系)在第二密度层底界面的密度等深图；

图7(c)为本发明的实施例2的研究区的第三密度层(古近系)在第三密度层底界面的密度等深图；

图7(d)为本发明的实施例2的研究区的第四密度层(中生界)在第四密度层底界面的密度等深图；

图7(e)为本发明的实施例2的研究区的第五密度层(古生界～蓟县系)在第五密度层底界面的密度等深图；

图7(f)为本发明的实施例2的研究区的第六密度层(长城系)在第六密度层底界面的密度等深图。

图8本发明复杂地质体的多密度界面反演方法流程图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

实施例1

以下首先通过模拟模型对该复杂地质体的多密度界面反演方法进行模拟，具体操作步骤如下：

S1、如图1所示，构建一个复杂地质体的理论模型：模型长度和宽度均为24km，包含四个密度层、三个地层分界面；

构建与目标复杂地质体相符的密度反演模型：以地表为x-y平面，沿深度方向为z轴方向，从地表沿z轴方向依次构建第一密度层、第二密度层、第三密度层和第四密度层；对应每个密度层的密度自上而下依次为：1.95g/cm³、2.15g/cm³、2.40g/cm³和2.70g/cm³；每个密度层与其下方密度层的分界面为该密度层的底界面，对应依次为：第一密度层底界面、第二密度层底界面和第三密度层底界面；每个密度层与其上方密度层的分界面为该密度层的顶界面，对应依次为：第一密度层顶界面、第二密度层顶界面、第三密度层顶界面和第四密度层顶界面；设定地表x-y平面为第一密度层的基准面，第一密度界面的平均深度处为第二密度层的基准面，第二密度界面的平均深度处为第三密度层的基准面；每个密度层的地层平均厚度为该密度层的基准面与下方密度层的基准面之间的距离，其自上而下依次为1.0km、1.0km和1.0km；

如图3(a-1)所示为第一密度层底界面等深图；如图3(a-2)所示为第二密度层底界面等深图；如图3(a-3)所示为第三密度层底界面等深图；其中，在图3(a-1)～图3(a-3)中，x轴方向和y轴方向上的网格间距均为0.2km；

S2、应用多密度分界面重力异常计算公式计算三个分界面在地面产生的重力总异常g，如图3(b)所示；接着，将该布格重力异常g代入本申请的基于位场分离目的的切割法计算公式中，计算出目标复杂地质体的各密度层的底界面及以上各密度层分界面在地面产生的重力异常；如图3(c-1)所示为第一密度层底界面在地面产生的重力异常示意图；如图3(c-2)所示为第二密度层底界面及以上各密度层分界面在地面产生的重力异常示意图；如图3(c-3)所示为第三密度层底界面及以上各密度层分界面在地面产生的重力异常示意图；

进而，求出各密度层的底界面在地面产生的重力异常；如图3(d-1)所示为第一密度层底界面在地面产生的重力异常示意图；如图3(d-2)所示为第二密度层底界面在地面产生的重力异常示意图；如图3(d-3)所示为第三密度层底界面在地面产生的重力异常示意图；

S3、应用本申请的位场大深度向下延拓方法将目标复杂地质体的各密度层的底界面在地面产生的重力异常延拓至基准面，即将各界面在地面产生的重力异常转换为在其密度层基准面引起的重力异常；其中，自地表向下，第一密度层的基准面为地表面，其他各密度层的基准面为其上层密度层底界面的平均深度处；

如图3(e-1)所示为第一密度层底界面在基准面产生的重力异常的示意图；如图3(e-2)所示为第二密度层底界面在基准面产生的重力异常的示意图；如图3(e-3)所示为第三密度层底界面在基准面产生的重力异常的示意图；

S4、将各密度层底界面在其基准面引起的重力异常以及各密度层的密度应用Parker-Oldenburg反演方法对进行反演，得到各密度层的底界面相对基准面的深度，再加上其对应基准面的深度得到各密度层底界面的实际深度；

如图3(f-1)所示为第一密度层底界面的反演等深示意图；如图3(f-2)所示为第二密度层底界面的反演等深示意图；如图3(f-3)所示为第三密度层底界面的反演等深示意图；

从图3(f-1)～图3(f-3)中可以看出，除第三密度层的底界面存在小幅度振荡现象，第一密度层和第二密度的底界面基本无震荡现象，即可以确定该反演方法能有效削弱Parker-Oldenburg反演方法的震荡现象。

将各密度层的底界面反演等深图中各网格点深度与模型真实界面中各网格点深度相减后求取绝对值，然后分别求取得到如下表1所示的各反演界面误差平均值、最大值和最小值进行比较。

表1：

从表1可以看出，各界面反演误差总体较小，证实本申请的基于密度反演模型实现的复杂地质体多密度界面反演方法结果与实际情况相符程度高，反演结果可信。

实施例2

以下以天津市团泊-大寺地区为例对本发明的多密度界面反演方法的具体实施方法进行进一步解释，其步骤如下：

S1、构建与目标复杂地质体相符的密度反演模型：

根据研究区内密度物性资料，天津市团泊-大寺地区重力场位于天津市中部重力高值区，区内的布格重力异常主要包括一个高的圈闭异常和一个低的圈闭异常，呈NE-SW向椭圆形展布，构造位于Ⅲ级构造单元沧县隆起上，区内划分4个Ⅳ级构造单元，分别为大成凸起、双窑凸起、白塘口凹陷和小韩庄凸起夹持在沧东断裂与天津断裂之间，区内断裂比较发育，如图4所示，其主要断裂为：天津断裂(F_Ⅰ-1)、大寺断裂(F_Ⅰ-3)和白塘口东断裂(F_Ⅰ-4)；区内断裂以NE走向为主，其中大寺断裂控制白塘口凹陷的形成演化，向下切入寒武系基底，向上断至明化镇组-第四系地层，走向NNE；其它断裂以新近系活动为主，后期(第四系)活动相对较弱。研究区内地表均被第四系覆盖，根据区内钻孔及周边资料，区内主要分布有太古宇、中元古界长城系、蓟县、中上元古界青白口系、下古生界寒武系、奥陶系、上古生界石炭系、二叠系、中生界侏罗系、白垩系、新生界的古近系、新近系和第四系等地层，除白塘口凹陷地层较全外，其余构造单元均有不同程度的地层缺失，且缺失地层均包含古近系和中生界。

根据研究区内密度物性资料，将研究区地层划分为7个密度层，分别为：第四系地层，密度值为2.00g/cm³；新近系地层，密度值为2.19g/cm³；古近系地层，密度值为2.38g/cm³；中生界地层，密度值为2.54g/cm³；古生界—中元古界蓟县系地层，密度值为2.68g/cm³；中元古界长城系地层，密度值为2.57g/cm³；太古宇地层，密度值为2.72g/cm³；

建立与研究区地层地质体情况相符的密度反演模型，具体如下：

以地表为x-y平面，沿深度方向为z轴方向，从地表沿z轴方向依次构建：第一密度层、第二密度层、第三密度层、第四密度层、第五密度层、第六密度层、第七密度层，对应密度依次为：2.00g/cm³、2.19g/cm³、2.38g/cm³、2.54g/cm³、2.68g/cm³、2.57g/cm³和2.72g/cm³；

各密度层与其上方密度层的分界面自上而下依次为：第一密度层顶界面、第二密度层顶界面、第三密度层顶界面、第四密度层顶界面、第五密度层顶界面、第六密度层顶界面；其中，第一密度层顶界面为地表平面；

各密度层与其下方密度层的分界面自上而下依次为：第一密度层底界面、第二密度层底界面、第三密度层底界面、第四密度层底界面、第五密度层底界面、第六密度层底界面；

设定x-y平面(即地面)为第一密度层的基准面，第一密度界面的初始平均深度处为第二密度层的基准面，第二密度界面的初始平均深度处为第三密度层的基准面、第三密度界面的初始平均深度处为第四密度层的基准面、第四密度界面的初始平均深度处为第五密度层的基准面、第五密度界面的初始平均深度处为第六密度层的基准面；

基于此，每个密度层的地层平均厚度，即各密度层的顶界面的初始平均深度处与其底界面的初始平均深度处之间的间距，自上而下依次为h₁、h₂、h₃、h₄、h₅、h₆；具体地，各密度层的地层平均厚度通过地质及钻孔资料初步估算得到，自上而下依次为：0.5km、1.5km、2.0km、2.5km、4.5km、6.0km。

S2、利用该研究区的布格重力异常g，并将其代入基于位场分离目的的切割法计算公式中，以计算出目标复杂地质体的各密度层的底界面及以上各密度层分界面在地面产生的重力异常，进而求出各密度层的底界面在地面产生的重力异常；

具体地，基于位场分离目的的切割法计算公式如下：

式(1)中，α是加权系数；g(x,y)是待位场分离点的重力异常；点第一次切割的区域场；/>是待位场分离点周围四点位场值的算术平均值，二者分别用如下公式表示：

α＝β+γ 式(3)，

式(2)和(3)中，r为切割半径，其取值为对应密度层底界面的初始平均深度至地表的地层平均厚度之和；g(x+r,y),g(x-r,y),g(x,y+r)和g(x,y-r)是半径为r的周围4个点的重力异常，如图2所示；β和γ分别根据x轴方向和y轴方向的半二阶差分量确定；

同理，重复上述公式(1)～(3)，将所得一次区域场作为第二次切割的原始异常，得到第二次切割的区域场/>如此迭代，直至满足：

式(4)中，为获取的第n次区域场，ε₀为允许误差；则区域场：

进而，得到近似深度为r以上的各地层分界面产生的重力异常：

L_loc＝g-R_reg 式(6)，

式(6)中，g为地表处实测的该研究区的布格重力异常；

运用上述理论，应用密度层上界面的初始平均深度h_i-1和下界面的初始平均深度h_i(i为密度层的编号)，求出目标复杂地质体的各密度层的底界面在地面产生的重力异常：

u_i＝L_{loc_i}-L_{loc_i-1} 式(7)，

式(7)中，L_{loc_i}为第i地层底界面及以上地层底界面产生的重力异常；

该步骤S2的求取过程自上而下进行，通过首先将切割半径r赋值为第一密度层的地层平均厚度，即r＝h₁，并代入至公式(1)～(6)中，得到第一密度层底界面在地面产生的重力异常L_{loc_1}，继而将切割半径r赋值为第一密度层的平均厚度和第二密度层的平均厚度之和，即r＝h₁+h₂，并代入至公式(1)～(6)中，得到第二密度层底界面及其以上各密度层分界面在地面产生的重力异常L_{loc_2}，以此类推，依次计算得到第三密度层底界面及其以上各地层分界面在地面产生的重力异常L_{loc_3}、第四密度层底界面及其以上各地层分界面在地面产生的重力异常L_{loc_4}、第五密度层底界面及其以上各地层分界面在地面产生的重力异常L_{loc_5}，和第六密度层底界面及其以上各地层分界面在地面产生的重力异常L_{loc_6}；进而，通过将第二密度层底界面及其以上各地层分界面在地面产生的重力异常L_{loc_2}减去第一密度层底界面在地面产生的重力异常L_{loc_1}，求得出第二密度层底界面在地面产生的重力异常u₂；依次类推，得到各密度层的底界面在地面产生的重力异常u₃、u₄、u₅和u₆；

S3、将目标复杂地质体的各密度层的底界面在地面产生的重力异常转换为在其密度层基准面引起的重力异常；具体地，

S301、假设第i密度层底界面在基准面产生初始重力异常：

u₁(x,y，0)＝u(x,y) 式(8)，

式(8)中，u(x,y)为密度层底界面在地面产生的重力异常；

S302、将初始异常值运用傅里叶变换法做向上延拓，延拓至地面，公式如下：

U₁(k_x,k_y，0)＝DFT(u₁(x,y，0)) 式(9)，

式(9)和(10)中，DFT表示傅里叶变换，IDFT表示傅里叶逆变换，k_x为x方向上的波数，k_y为y方向上的波数；

S303、对初始值u₁(x,y,0)进行校正，公式如下：

u₂(x,y，0)＝u₁(x,y，0)-(u₁(x,y,r_i)-u(x,y)) 式(11)，

式(11)中，u₂(x,y,0)为校正后的重力异常；

S304、重复S301～S303，直至上延计算值与密度层在地面上的异常值的差值小到可以忽略的程度，即：

|u_n(x,y,r_i)-u(x,y)|＜ε 式(12)，

式(12)中，u_n(x,y,r_i)为第n次向上延拓后的重力异常，ε为允许误差；所得u_n(x,y,0)即为u(x,y)向下延拓r_i-1的重力异常，即第i密度层底界面在其基准面产生的重力异常；

S4、将各密度层底界面在基准面引起的重力异常依次通过Parker-Oldenburg反演方法实现界面反演，求出底界面相对基准面的深度值(Δh)：

进而，将该深度值加上基准面深度值得出各地层底界面的深度值，即各界面的实际深度：

h＝r_i-1+Δh 式(13)。

为弥补区内钻孔的不足和分布的不均匀性，应用区内已有的大地电磁和地震成果资料制作8个虚拟钻孔(包括四个地震虚拟钻孔Vs1～Vs4和四个Mt虚拟钻孔Vmt1～Vmt4)，再加上如图4所示的该地区的7个真实钻孔(D1～D7)，共同验证各地层分界面反演深度的准确性。如下表2～6所示为研究区内钻孔资料对本发明的密度界面反演方法的结果精度评价。

表2：第四系底界面钻孔与反演成果对比表

表3：新近系底界面钻孔与反演成果对比表

表4：古近系底界面MT虚拟钻孔与反演成果对比表

表5：中生界底界面地震虚拟钻孔与反演成果对比表

表6：蓟县系底界面钻孔与反演成果对比表

从图7(c)、图7(d)、图7(e)和图7(f)中可以看出，白塘口凹陷内数据存在小幅度的震荡现象，且震荡程度由图7(c)到图7(f)逐渐增大；从表2～表6中可以看出，反演界面深度与钻孔揭露深度最大偏差为398m(表6，D3)，最小偏差9m，绝大部分偏差在100m附近，相对偏差最大值为18.07％(表6，D3)，最小值为1.91％，绝对大部分在5.0％附近，总体偏差较小，从各钻孔偏差可以看出，与地层偏差较大钻孔主要集中在D2、D3、D4、D6，分析认为：D2、D3、D4钻孔分布在F_Ⅰ-3断裂附近，该处布格重力异常的梯度带较陡，地层底界面坡度较大，由此造成了钻孔揭露深度与计算深度偏差较大，D6钻孔位于研究区的边界，由于边界效应影响，造成该处反演误差较大，因此，通过数据分析证实反演结果总体是可信的，证明反演方法是可靠的。

与目前常用的Parker-Oldenburg反演方法预测复杂地质体的密度层分布情况相比，除了对于第一密度层底界面(参见表2)来说，采用本申请的反演方法获得的反演深度与采用Parker-Oldenburg反演方法获得的反演深度相同；其余各密度层底界面的反演结果则通过采用本申请的反演方法均得到明显提升，其反演精度大部分绝对值提升50m～150m，提升1％～6％；其原因在于：反演第一密度层底界面时，两种方法反演的基准面均为地面，故反演结果相同；而反演其它密度层底界面时，本申请进行反演的基准面为各密度层顶界面初始平均深度处，而Parker-Oldenburg方法反演各密度层的基准面均为地面，故造成上述差异；综上所述，本申请的反演方法有效解决了Parker-Oldenburg反演方法存在的收敛性差的缺陷，提高反演结果的准确率，有效获得复杂地质体密度层实际分布状况。

Claims (4)

1.一种复杂地质体的多密度界面反演方法，其特征在于，步骤如下：

S1、根据目标复杂地质体的地层密度分布资料，构建与目标复杂地质体相符的密度反演模型；

S2、将地表处实测的布格重力异常g代入基于位场分离目的的切割法计算公式中，计算出目标复杂地质体的各密度层的底界面及以上各密度层分界面在地面产生的重力异常，进而求出各密度层的底界面在地面产生的重力异常；

S3、应用位场大深度向下延拓方法将目标复杂地质体的各密度层的底界面在地面产生的重力异常延拓至基准面；其中，自地表向下，第一密度层的基准面为地表面，其他各密度层的基准面为其上层密度层底界面的平均深度处；

S4、将各密度层底界面在其基准面引起的重力异常依次通过Parker-Oldenburg反演方法实现界面反演，以得到各界面相对基准面的深度，进而求出各密度层底界面的实际深度。

2.根据权利要求1所述的复杂地质体的多密度界面反演方法，其特征在于，在步骤S1中，构建与复杂地质体相符的密度反演模型的具体方法为：

以地表为x-y平面，沿深度方向为z轴方向，从地表沿z轴方向依次构建第一密度层、第二密度层、第三密度层、…、第n密度层；其中，

对应每个密度层的密度依次为：ρ₁、ρ₂、ρ₃、…、ρ_n；

对应每个密度层与下方密度层的分界面为其底界面，各密度层的底界面依次为：第一密度层底界面、第二密度层底界面、第三密度层底界面、…、第n-1密度层底界面；

对应每个密度层与上方密度层的分界面为其顶界面，各密度层的顶界面依次为：第一密度层顶界面、第二密度层顶界面、第三密度层顶界面、…、第n密度层底界面；其中，第一密度层顶界面为地表平面；

基于各密度层分界面起伏是造成地表重力异常的主要原因，设定地表x-y平面为第一密度层的基准面，第一密度层底界面的初始平均深度处为第二密度层的基准面，第二密度层底界面的初始平均深度处为第三密度层的基准面、…、第n-1密度层底界面的初始平均深度处为第n密度层的基准面；

对应每个密度层的地层平均厚度，即密度层的基准面与下方密度层的基准面之间的距离，依次为h₁、h₂、h₃、…、h_n-1；各密度层的地层平均厚度通过已知地质和/或钻孔资料获得。

3.根据权利要求1所述的复杂地质体的多密度界面反演方法，其特征在于，在步骤S2中，基于位场分离目的的切割法计算公式如下：

式(1)中，α是加权系数；g(x,y)是待位场分离点的重力异常；点第一次切割的区域场；/>是待位场分离点周围四点重力异常值的算术平均值；其中，/>为g(x,y)点第一次切割的区域场；

α＝β+γ 式(3)，

式(2)和(3)中，r为切割半径，其取值为对应密度层的基准面至地表的各密度层的地层平均厚度之和；g(x+r,y),g(x-r,y),g(x,y+r)和g(x,y-r)是半径为r的周围4个点的重力异常；β和γ分别根据x轴方向和y轴方向的半二阶差分量确定；

同理，重复上述公式(1)～(3)，将所得一次区域场作为第二次切割的原始异常，得到第二次切割的区域场/>如此迭代，直至满足：

式(4)中，为获取的第n次区域场，ε₀为允许误差；

则区域场：

进而，得到近似深度为r以上的各地层分界面产生的重力异常

L_loc＝g-R_reg 式(6)，

式(6)中，g为地表处实测的布格重力异常；

基于此，通过该密度层顶界面的平均深度r_i-1和该密度层底界面的平均深度r_i，求出目标复杂地质体的各密度层的底界面在地面产生的重力异常：

u_i＝L_{loc_i}-L_{loc_i-1} 式(7)，

式(7)中，i为该密度层的编号，L_{loc_i}为第i密度层底界面及以上地层底界面在地面产生的重力异常。

4.根据权利要求1所述的复杂地质体的多密度界面反演方法，其特征在于，在步骤S3中，其具体转换步骤如下：

S301、假设第i密度层底界面在其基准面产生初始重力异常：

u₁(x,y，0)＝u(x,y) 式(8)，

式(8)中，u(x,y)为该密度层底界面在地面产生的重力异常；

S302、将初始异常值运用傅里叶变换法做向上延拓，延拓至地面，公式如下：

U₁(k_x,k_y，0)＝DFT(u₁(x,y，0)) 式(9)，

式(9)和(10)中，DFT表示傅里叶变换，IDFT表示傅里叶逆变换，k_x为x方向上的波数，k_y为y方向上的波数；

S303、对初始值u₁(x,y,0)进行校正，公式如下：

u₂(x,y，0)＝u₁(x,y，0)-(u₁(x,y,r_i)-u(x,y)) 式(11)，

式(11)中，u₂(x,y,0)为校正后的重力异常；

S304、重复S301～S303，直至上延计算值与密度层在地面上的异常值的差值小到可以忽略的程度，即：

|u_n(x,y,r_i)-u(x,y)|＜ε 式(12)，

式(12)中，u_n(x,y,r_i)为第n次向上延拓后的重力异常，ε为允许误差；所得u_n(x,y,0)即为u(x,y)向下延拓r_i-1的重力异常，即第i密度层底界面在其基准面产生的重力异常。