专利摘要

本发明公开了一种基于归一化异常权函数的高分辨率密度反演方法，包括：具体反演过程将地下三维空间分为有限个密度一定的长方体单元，由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加，当存在n个观测点和m个长方体模型单元时，重力异常正演表达式用矩阵相乘的方式来表示。本发明提出归一化异常权和深度权相结合的密度反演方法，能有效地提高反演结果的水平和垂直分辨率，且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行，因此反演过程更加稳定，能有效地压制噪声的干扰。

权利要求

1.一种基于归一化异常权函数的高分辨率密度反演方法，其特征在于，包括：

具体反演过程将地下三维空间分为有限个密度一定的长方体单元，由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加，当存在n个观测点和m个长方体模型单元时，重力异常正演表达式用矩阵相乘的方式来表示：

（1）

其中，代表n个观测点上的重力异常；为重力正演核函数矩阵；为m个长方体模型单元的密度；

核函数矩阵为：

（2）

为万有引力常数；为观测点坐标；为地下长方体单元的最小与最大角点的坐标；

通过以下形式来求解此欠定问题：

（3）

其中为观测异常与计算异常差值的二范数平方，代表数据的拟合泛函；为权函数，用于约束反演结果使结果更加符合实际；为正则化参数，起到平衡数据拟合函数与权函数的作用；

在权函数中引入深度加权函数：

（4）

（5）

为第i层长方体单元的埋深，为一个经验常数，用于改变深度权对浅部压制效果的大小，k为地下块体纵向上分割的个数；

所述的基于归一化异常权函数的高分辨率密度反演方法，还包括：

在公式（6）中将权函数中的原始异常g替换为异常的垂向导数，反演结果还会更加收敛：

（6）

（7）

其中，为相应i层深度对应高度的异常值，不同高度的值通过延拓的方式计算获得；为重力原始异常的最大值，作用是对整体数据进行归一化处理，避免不同的重力数据整体数值的大小对反演产生影响；a为常数，为了均衡异常权与其他权之间的数量级；的作用是将括号内的一维数组变为n×n的对角矩阵，只有对角线上元素与A中的元素一一对应，其余部分为0；k为地下块体纵向上分割的个数。

2.根据权利要求1所述的基于归一化异常权函数的高分辨率密度反演方法，其特征在于，还包括：

将（5）式、（7）式结合后代入到（3）式中的权函数里，得到具体的反演函数为：

（8）

（8）式的矩阵形式为：

（9）

其中为重力归一化异常与深度权相结合后的新型权函数，起到提高反演结果分辨率的作用；

三维密度反演最终为求解（9）式，获得的数值，即为地下空间块体的密度分布。

说明书

技术领域

本发明涉及密度反演方法，具体涉及一种基于归一化异常权函数的高分辨率密度反演方法。

背景技术

目前，三维密度反演是重力数据反演方法中反演目标最为全面的方法。因为密度反演是一个求解欠定方程的问题，也就无法避免地球物理反演中的多解性问题，如何通过重力场的物理性质以及地质条件等密度分布必须遵守的规律来限制反演结果，从而减小反演的多解性，获得更加符合实际的密度分布特征至关重要。

发明内容

本发明的主要目的在于提供一种基于归一化异常权函数的高分辨率密度反演方法。

本发明采用的技术方案是：一种基于归一化异常权函数的高分辨率密度反演方法，包括：

具体反演过程将地下三维空间分为有限个密度一定的长方体单元，由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加，当存在n个观测点和m个长方体模型单元时，重力异常正演表达式用矩阵相乘的方式来表示：

（1）

其中， 代表n个观测点上的重力异常； 为重力正演核函数矩阵； 为m个长方体模型单元的密度；

核函数矩阵为：

（2）

为万有引力常数； 为观测点坐标； 为地下长方体单元的最小与最大角点的坐标；

通过以下形式来求解此欠定问题：

（3）

其中 为观测异常与计算异常差值的二范数平方，代表数据的拟合泛函； 为权函数，用于约束反演结果使结果更加符合实际； 为正则化参数，起到平衡数据拟合函数与权函数的作用；

在权函数 中引入 深度加权函数：

（4）

（5）

为第i层长方体单元的埋深， 为一个经验常数，用于改变深度权对浅部压制效果的大小，k为地下块体纵向上分割的个数。

进一步地，基于归一化异常权函数的高分辨率密度反演方法还包括：

在公式（6）中将权函数中的原始异常g替换为异常的垂向导数，反演结果还会更加收敛：

（6）

（7）

其中， 为相应i层深度对应高度的异常值，不同高度的值通过延拓的方式计算获得； 为重力原始异常的最大值和最小值，作用是对整体数据进行归一化处理，避免不同的重力数据整体数值的大小对反演产生影响；a为常数，为了均衡异常权与其他权之间的数量级； 的作用是将括号内的一维数组 变为n×n的对角矩阵，只有对角线上元素与A中的元素一一对应，其余部分为0；k为地下块体纵向上分割的个数。

更进一步地，基于归一化异常权函数的高分辨率密度反演方法还包括：

将（5）式、（7）式结合后代入到（3）式中的权函数 里，得到具体的反演函数为：

（8）

（8）式的矩阵形式为：

（9）

其中 为重力归一化异常与深度权相结合后的新型权函数，起到提高反演结果分辨率的作用；

三维密度反演最终为求解（9）式，获得 的数值，即为地下空间块体的密度分布。

本发明的优点：

本发明针对现有重力密度反演在同一平面上权函数相同所导致的结果分辨率低的问题，提出归一化异常权和深度权相结合的密度反演方法，新型的权函数是依据异常与深度的对应关系来建立，因此所提出方法能有效地提高反演结果的水平和垂直分辨率，且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行，因此反演过程更加稳定，能有效地压制噪声的干扰。通过模型正演异常对所提出密度反演方法开展试验，反演结果表明该方法相对常规方法能获得更加准确和收敛的结果，且对于噪声具有良好的压制效果，并且利用重力垂向梯度异常进行加权时会取得更加高的分辨率。将该方法应用于山东地区的铁矿实测重力异常的反演，获得铁矿的空间分布特征，为下一步精细勘探提供了重要的支撑。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1是本发明的基于深度权的三维密度反演结果切片图；

图2是本发明的不同深度权函数的获得方式示意图；

图3是本发明的矿区实测重力异常图；

图4是本发明的总水平导数及断裂划分图；

图5是本发明的实测重力异常的平均对数功率谱图；

图6是本发明的位场分离后矿区重力异常图；

图7是本发明的三维密度反演结果高值三维成像图；

图8是本发明的估计铁矿范围解释图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明的具体反演过程是将地下三维空间分为有限个密度一定的长方体单元，由地下密度分布所产生的异常近似为所有地下块体单元分别对地表观测点产生的异常的叠加，当存在n个观测点和m个长方体模型单元时，重力异常正演表达式用矩阵相乘的方式来表示：

（1）

其中， 代表n个观测点上的重力异常； 为重力正演核函数矩阵； 为m个长方体模型单元的密度；

通过Holstein（1999）提出的长方体模型的正演计算知，核函数矩阵为：

（2）

为万有引力常数； 为观测点坐标； 为地下长方体单元的最小与最大角点的坐标；

三维密度反演就是求解式（1）中的 ，由于地下块体的数量远大于地表观测点个数，所以解方程问题为欠定问题；通过以下形式来求解此欠定问题：

（3）

其中 为观测异常与计算异常差值的二范数平方，代表数据的拟合泛函； 为权函数，用于约束反演结果使结果更加符合实际； 为正则化参数，起到平衡数据拟合函数与权函数的作用；

由于核函数矩阵G随着深度的增加快速衰减，所以反演结果在浅部较为集中，这经常与实际地质情况不符，在权函数 中引入 深度加权函数：

（4）

（5）

为第i层长方体单元的埋深， 为一个经验常数，用于改变深度权对浅部压制效果的大小，k为地下块体纵向上分割的个数。

图1为本发明基于深度权的三维密度反演结果切片图，（a）为重力异常剖面; （b）为深度权反演结果切片。

图1中为只用深度加权的三维密度反演结果的切片图，白色虚线为模型体的实际位置。从图中看出，由于核函数矩阵G随着深度的增加快速衰减，所以反演结果在深部具有较低的分辨率，引入下半空间的空间分布权函数，从而获得较高水平和垂直分辨率的结果。

图2为本发明不同深度权函数的获得方式示意图，（a）为原始异常向上延拓5m；（b）为原始异常向上延拓35m；（c）为原始异常向上延拓65m；（d）为原始异常向上延拓95m。

因为不同的重力异常的极值是不同的，所以为了去除异常极值大小带来的影响，将异常进行归一化处理，这样处理后的结果仅能体现异常高低的特征。又因为实测重力异常显示的是地下的综合结果，为了使权函数在深部时能主要体现深部的重力异常特征，利用向上延拓的方式来确定不同深度的权函数所利用的重力异常。如图2所示，

将向上延拓5m的异常数据的计算结果作为5m深处的权函数，按照这种方式形成一个和核函数相对应的对角矩阵。具体形式如下：采用不同高度的延拓异常进行归一化处理，从而反演结果将更加收敛；而不会像以往方法一样产生十分发散的结果，因此以往方法反演结果发散的区域异常权函数会使其无变化，接近零值，从而结果将更加收敛；

在公式（6）中将权函数中的原始异常g替换为异常的垂向导数，反演结果还会更加收敛：

（6）

（7）

其中， 为相应i层深度对应高度的异常值，不同高度的值通过延拓的方式计算获得； 为重力原始异常的最大值和最小值，作用是对整体数据进行归一化处理，避免不同的重力数据整体数值的大小对反演产生影响；a为常数，为了均衡异常权与其他权之间的数量级； 的作用是将括号内的一维数组 变为n×n的对角矩阵，只有对角线上元素与A中的元素一一对应，其余部分为0；k为地下块体纵向上分割的个数；将（5）式、（7）式结合后代入到（3）式中的权函数 里，得到具体的反演函数为：

（8）

（8）式的矩阵形式为：

（9）

其中 为重力归一化异常与深度权相结合后的新型权函数（空间分布权函数），起到提高反演结果分辨率的作用；

三维密度反演最终为求解（9）式，获得 的数值，即为地下空间块体的密度分布。采用预条件共轭梯度法来求解这一病态方程，共轭梯度法具有内存占用少和稳定性高的优势，比较适合如今数据量越来越大的重力反演。

实际数据处理：

为了验证的方法在实际数据中的效果，将提出的方法应用于山东铁矿矿区的实测重力数据中。该地区位于华北板块，地处鲁西与鲁中隆起的西北部。地表存在大面积的第四系覆盖，并且隐伏碳酸盐岩分布广泛。附近已经发现的矿床资料显示，该地区主要为矽卡岩型铁矿床。矿体主要在闪长岩体与碳酸盐岩的接触带上，所以侵入岩的接触带为重要的找矿位置，为高密度体。重力数据主要是通过寻找重力异常高值的侵入岩来判断铁矿位置。将地下划分为40×30×20个长方形块体，每个长方体的大小为750×1500×500m。对数据重新进行网格化，获得观测数据个数为40×30个，观测点坐标对应长方体中心位置。实测异常如图3所示。我们对实测重力数据进行三维密度反演，并筛选出反演结果中的高密度值，来判断铁矿位置与范围。

因为断裂可以为岩浆上涌提供通道，矿体所在的岩浆侵入的位置经常会存在断裂构造，首先对原始重力异常进行总水平导数的计算来划分出当地的断裂构造，如图4所示，图中黑色虚线为断裂位置，即总水平导数的高值区域。

因为重力异常为地下所有密度的综合相应，所以对于实测异常中存在的背景场为干扰反演的异常。通过匹配滤波法对实测异常进行位场分离，图5为异常的平均对数功率谱，将波数较小的区域滤掉，获得浅部的剩余重力异常，如图5所示。

从图6的重力异常中我们可以看出，重力高值区大致分为六块，其中Ⅰ、Ⅲ、Ⅵ处的重力异常较高，说明这三处的高密度体范围较大，图6中西南角的高值不完全在勘探区域中，所以无法反演出准确的结果，所以本次处理并不考虑此高值异常。利用提出的方法对此异常进行反演，获得的反演结果如图7所示。

图7为反演结果中高密度块体的三维成像，表面覆盖的为原始重力异常，并且按照图4所示用罗马数字标出了六块高值异常的位置，很好的对应了反演结果中的六块地下高密度体，并且矿体倾向与异常显示完全吻合，从图7中可以看出，总共有六个区域是高密度块体汇聚的地区，这些区域的铁矿分布不是相互连通的，通过对高密度块体的坐标可以获得各个区域铁矿的埋深与范围，最终得到图8。

图8展示的是矿体的水平范围，即该地区的三维密度反演结果解释图。通过图7的三维结果可以获得各处矿体的顶面埋深，其中图8中Ⅰ和Ⅴ的顶面埋深为800m；Ⅲ和Ⅵ的顶面埋深为1000m；Ⅱ的顶面埋深为1200m；Ⅳ的顶面埋深为1400m。

本发明针对现有重力密度反演在同一平面上权函数相同所导致的结果分辨率低的问题，提出归一化异常权和深度权相结合的密度反演方法，新型的权函数是依据异常与深度的对应关系来建立，因此所提出方法能有效地提高反演结果的水平和垂直分辨率，且由于在获取不同深度异常特征过程中主要依据向上延拓函数进行，因此反演过程更加稳定，能有效地压制噪声的干扰。通过模型正演异常对所提出密度反演方法开展试验，反演结果表明该方法相对常规方法能获得更加准确和收敛的结果，且对于噪声具有良好的压制效果，并且利用重力垂向梯度异常进行加权时会取得更加高的分辨率。将该方法应用于山东地区的铁矿实测重力异常的反演，获得铁矿的空间分布特征，为下一步精细勘探提供了重要的支撑。

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

一种基于归一化异常权函数的高分辨率密度反演方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。