一种重力多界面反演方法
【技术领域】
[0001] 本发明涉及地球物理勘探技术领域，特别是涉及一种重力多界面反演方法，其为 重力勘探中的密度界面反演方法，可以得到地下多层界面的起伏分布，且具有较高的运算 效率。
【背景技术】
[0002] 地质勘探主要应用的地球物理勘探方法包括重力勘探、磁法勘探、电法勘探和地 震勘探等。重力勘探通过利用重力仪观测地下物质密度差异引起的重力异常以查明地下的 地质构造和岩性异常体。相较于地震等勘探方法而言，重力勘探具有横向分辨率好、成本低 廉等优点。重力勘探分为三个大的环节:野外重磁异常采集、重磁异常处理和重磁异常解 释。
[0003] 从地质角度，重力反演问题的目标在于寻找、研究或推断金属或非金属矿体和研 究地质构造。从地球物理角度，重力反演问题的目标在于确定地质体的几何和物性参数，确 定物性分界面的深度和起伏，以及密度的分布等等。
[0004] 总体来看，当前重力反演方法研究主要可分为两大类。第一类是密度分布研究，着 重于物性（即密度)在地下半空间中的分布情况的研究;第二类是界面反演研究，关注点在 于明显的密度分界面，一般用于处理密度已知(或密度变化规律已知）的情形下密度分界面 的起伏变化情况。
[0005] 重力界面反演方法的研究起步较早，自20世纪60年代起，Bott(1960)、Cordell et al. (1968)等即对重力异常在空间域的正反演问题进行了初步研究。这个时期的计算方法 都是基于空间域的，通过长方体模型进行拟合，逐次逼近消去剩余值，无法利用先验信息。 随着计算机技术的发展，快速傅里叶变换技术(FFT)投入使用，Parker(1973)借助FFT，给出 了频率域计算地下单层物质重力异常的正演公式，较以往的离散模型方法速度提高了至少 一个数量级(冯锐等，1986)，极大地提高了计算效率。随后由01denburg(1974)加以发展，建 立了反演问题迭代求解方法，后人称为Parker-Oldenburg反演法，目前已经成为了地球物 理领域最常用的重力反演方法之一，具有理论严谨、适应性强和快速有效的特点。除此以 外，也有其他科研人员提出了一些反演方法，如Spector et al. (1970)的功率谱分析方法、 Chavez et al. (1985)的Parker频谱展开线性反演法和Pederson(1977)的广义矩阵反演法 等。这些方法里有一些是对Parker-Oldenburg方法的改进或引申，另外一些方法则常常用 于某些特定情况，比如功率谱方法适用于对物质界面的大体预估。不过总体而言，以上方法 都存在同样的局限性:假定地下物质为单一界面模型，如此通过重力异常的反演只能得到 单一界面的深度变化情况。针对这一问题，王万银等（1993)给出了双界面反演方法，实质采 用了单层地壳双界面模型，方法上与Parker-Oldenburg方法非常类似，需要额外考虑的是 "调整深度"的合理选取以加快收敛速率。该方法理论较为简易，也可以引入控制点进行约 束(黄松，2010)，但模型仅针对单层物质的顶底界面，应用范围有限，对诸如俯冲带地区"双 莫霍"的反演问题无能为力，且只能反演两个界面，而且密度、顶界面、底界面三者的反演进 行中，对其中一项的反演总是相当于假定其他两个参数为已知确定值，使用上造成了一定 的局限性。为避开这种将未知界面"假定为已知确定值"的情况，一些基于双界面模型发展 的方法(冯旭亮等，2014)干脆以地形作为双界面模型的顶界面，实质又退化为了单界面反 演问题，且因模型本身受限，关注点转变为沉积盆地基底研究，方法无法用于深层界面反 演。其他的一些非线性全局最优化方法(朱自强等，1995; Snopek，2005;于鹏等，2007;李倩 等，2010;明圆圆等，2012)采用非线性全局最优化方法，包括模拟退火算法、遗传算法和人 工鱼群算法等，这类算法为避免线性/非线性方程组大规模矩阵求解问题而采用了零阶算 法，其优点是对初始模型的依赖程度低，但参数取值范围和选取参数存在随机性，采用启发 式算法导致每次的反演结果变动很大，且求解效率不高。
[0006] Martins et al. (2011)和冯旭亮等(2014)针对沉积盆地基底反演问题，采用过全 变差泛函与控制点信息相结合的方式进行研究。但需要指出，Martins et al. (2011)和冯 旭亮等（2014)的方法针对的仅是单一界面反演问题，而且反演对象为浅层物质，无法对深 层物质进行反演。
[0007] 综上所述，现有的地球物理勘探方法难以处理更为复杂的多界面反演情况，并且 运算效率不高。
【发明内容】
[0008] 本发明主要解决的技术问题是提供一种重力多界面反演方法，可以得到地下多层 界面的起伏分布，且具有较高的运算效率和反演精度。
[0009] 为解决上述技术问题，本发明采用的一个技术方案是：提供一种重力多界面反演 方法，该方法包括以下步骤：
[0010] (1)获取观测点的信息、重力观测值、初始模型范围、地下介质信息(地下介质的相 关数据）、待反演层位信息、与模型泛函各部分的正则化权重；
[0011] (2)对地下介质进行网格剖分；
[0012] (3)计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到 重力正演值；
[0013] (4)根据重力观测值与重力正演值，计算失配泛函？:^)，其中，X为待求的各柱状 体底界面深度；
[0014] (5)计算全变差函数F2(x);
[0015] (6)计算模型界面反演偏差函数F3(x);
[0016] (7)根据步骤(5)~(6)的计算结果构成模型泛函，并根据模型泛函各部分的正则 化权重λ2~λ 3与失配泛函Fl(X)-起构成目标函数：
[0017] F(x) =Fi(x)+A2F2(x)+hF3(x);
[0018] (8)采用迭代方法使目标函数最小，期间根据迭代的结果与给定的不等约束，不断 修改模型以使反演结果趋向真实值。
[0019] 其中，步骤(3)进一步包括利用长方体重力正演公式计算初始模型网格剖分后的 各柱状体对观测点的重力效应并进行累加，得到重力正演值。
[0020] 其中，步骤(1)进一步包括：
[0021] 获取用于约束反演模型的信息，其包括控制点信息、模型各层上下界信息、失配泛 函协方差矩阵以及模型界面反演偏差矩阵，其中，控制点信息包括控制点位置信息；
[0022] 步骤(4)和步骤(5)之间还包括:步骤(45)存在控制点信息时，计算控制点信息函 数F 4(x);
[0023] 步骤(7)包括:根据步骤(45)、（5)以及(6)计算结果构成模型泛函，并根据模型泛 函各部分的正则化权重λ2~λ 4与失配泛函FKx)-起构成目标函数：
[0024] F(X) =Fl(X)+A2F2(X)+X3F3(X)+X4F4(X)。
[0025] 其中，步骤(8)具体包括:不断重复步骤(3)~(8)，直至达到迭代退出条件为止，最 终求得各层界面深度反演结果，其中，迭代退出条件包括迭代次数达到预设的次数以及迭 代的结果误差小于预设的误差。
[0026] 其中，观测点的信息包括观测点数目和各观测点位置(x，y，z)，其中，用于表示X值 的X轴和用于表示y值的Y轴位于同一平面，用于表示z值的Z轴向下为正，X轴、Y轴以及Z轴三 者构成右手直角坐标系；
[0027]重力观测值包括重力异常观测值和/或重力梯度观测值，当进行模型实验研究时， 由标准模型计算出重力异常正演值和/或重力梯度正演值，若使用重力梯度观测值，则重力 梯度观测值至少包括梯度张量的其中一个分量；
[0028] 地下介质的相关数据包括地下介质层数、各层介质的密度和底界面初始深度；
[0029] 待反演层位信息包括要反演的层位编号，层位编号为由浅部到深部排列或自定义 的编号。
[0030] 其中，重力梯度张量的形式为：
[0032]其中G为万有引力常数，P为密度，v为体积微元，尸和f分别为观测点位置和场源 位置，算子▽作用于观测点，在χ-γ-ζ坐标系下的重力梯度为
[0036]其中，（114243) = (1，7,2)，.没为偏微分符号，&为在坐标1上的偏微分。
[0037] 其中，步骤(4)中，失配泛函形式为
[0039] 其中，fi(x)的取值如下：
[0040] 当无重力梯度观测值时：;1；'1(1) = 8瓜￥。1)；^)-8瓜￥。￡11(叉），
[0041 ]当有重力梯度观测值时：
7
[0042 ] 其中，gra Vobs (X)为重力观测值，gra Veal (X)为重力正演值，gra vgraobs (X)为重力梯 度观测值，gravgra^U)为重力梯度正演值，Cddg为失配泛函协方差矩阵，Μ为观测点数目。 [0043]其中，步骤(45)中，控制点信息函数为
[0045] 其中
[0046] f4(x)=k_Wx
[0047] 其中，W是一个BXN的稀疏矩阵，B为控制点数目，N为待求柱状体个数，稀疏矩阵里 面元素非0即1，保证在相应控制点i上的柱状体底界面深度^趋近于k里的对应元素值，而n 4 为归一化