基于N阶多项式密度函数的重力异常反演方法与流程

本发明涉及一种基于n阶多项式密度函数的重力异常反演方法，属于地球物理勘探技术领域。

背景技术：

目前，重力反演主要分为两大类，一类是界面反演，该类方法一般假设地层的密度分布已知，通过观测重力异常与理论计算重力异常之间的拟合，确定某个地层界面的深度，此类方法应用较为广泛的案例是求取盆地的基底深度；另一类是物性反演，将地下进行网格化剖分，划分为许多矩形单元，每个矩形单元的密度值为一个常数。反演时，通过不断修改各个单元的密度值使得观测的重力异常与计算重力异常的残差在允许的范围内，不断迭代获得地下的密度分布。

近年来，许多学者关注到变密度体的重力异常计算问题。考虑密度随深度变化的情况，把地下密度分布表示为与深度相关的指数函数，抛物线函数，双曲线函数等不同的函数关系，并推导了相关的正演公式，利用变密度函数来研究盆地基底深度，取得了一些研究成果。但还存在以下问题：(1)密度随深度的变化关系多种多样，线性的或指数、抛物线、双曲线等函数关系不足以描述地下真实的密度分布特征；(2)由于密度函数形式的多样化使反演中选取某种密度-深度函数困难，且反演重力核函数不能统一；(3)目前多数反演都将反演界面和反演密度分布分开考虑，即反演界面时假设特定的密度函数，而事实上，地下真实密度分布与假定的密度函数是不相同的，这无疑增加了反演误差；(4)现存的物性反演方法，无论是规则网格,如矩形网格,还是不规则网格,如三角网格，都是在横向和纵向上均进行剖分，若网格单元数目过少，则反演不够精细，若网格单元数目越多，多解性问题就越严重。

技术实现要素：

针对现有技术存在的上述缺陷，本发明提出了一种新的基于n阶多项式密度函数的重力异常反演方法。在建立n阶多项式密度函数的目标体重力异常解析表达式的基础上，通过反演多项式密度函数的系数确定地下物质密度分布的方法技术。

本发明是采用以下的技术方案实现的：基于n阶多项式密度函数的重力异常反演方法，包括如下步骤：

(1)对实测重力异常进行预处理，得到研究目标体产生的重力异常；

(2)单元划分：将研究区划分为竖向并排的nr个矩形单元，各个矩形单元纵向长度相同，且大于或等于研究目标区的深度；

(3)设定密度多项式函数的最大阶数，构建重力核函数矩阵：将研究区划分为nr个竖向并排的矩形单元，每个直立矩形单元的密度表示为随深度z变化的函数：

其中，nz是z的最大阶数，aj表示多项式中各项的系数。那么，一个直立矩形单元在测点(xo,zo)处产生的重力异常可由下式计算：

则所有直立矩形单元在一个测点的重力异常为：

上式可简写为

其中，

若重力异常的测点数为nm个，重力异常与多项式密度函数的系数便可表示为如下矩阵形式：

其中，为重力核函数；

为要求解的多项式系数；

为实测重力异常值；

(4)构建反演目标函数，建立相应的重力反演线性方程组，求解方程，使计算的理论重力异常与实测重力异常残差在允许的范围内，得到相应各矩形单元密度函数的多项式系数；

目标函数的构建形式具体如下：

其中，

第一项为数据残差项；第二项为先验模型约束；第三项为横向平滑约束；第四项为纵向平滑约束；λ1,λ2,λ3分别为各项的权重系数；

p为对角分块矩阵，q1,q2,…,分别为各分块矩阵，且可统一为如下形式：h＝1,2,…nr；d为每个矩形单元中施加约束的深度点数，z1,z2,…,zd分别为各约束点的深度值，σ0为先验密度差模型，在密度模型趋势未知时，取σ0＝0；

h为密度差横向一阶差分矩阵，

v为密度差纵向二阶差分矩阵，w1,w2,…,分别为各分块矩阵，且可统一为如下形式：

(5)利用反演的多项式系数计算整个研究区的地下密度分布。通过以上具体步骤的处理，实现了对研究区重力异常的反演。

本发明的原理在于：利用重力异常反演地下的密度分布是重力勘探的主要目标之一。把密度随深度变化的关系用任意阶的多项式函数表示。多项式函数的最大阶数越高，就可以描述越复杂的地下密度分布情况。给定多项式系数便可以正演计算出相应产生的重力异常，使计算的理论重力异常与实测重力异常在一定约束下达到最小，即可获得地下的密度分布。本发明的核心在于：(1)地下物质密度的多项式函数表示；(2)以多项式密度函数的系数为未知量的重力反演的核函数矩阵的构建；(3)约束条件和先验信息的施加；(4)建立重力反演线性方程组并求解密度函数的多项式系数。

本发明的有益效果是：采用本发明所述的基于n阶多项式密度函数的重力异常反演方法，(1)密度采用多项式函数表示，比常密度假设更符合实际情况，重力异常计算精度更高，反演结果更符合地质规律；(2)将地下划分为竖直并排的矩形单元，只需反演各个矩形单元的多项式系数，比常规的网格剖分单元数少，求解未知量少，降低了反演的多解性，也提高了反演的计算效率；(3)构建的目标函数灵活，可以方便的施加约束和由地质调查、测井和其他地球物理手段获得的先验信息。

附图说明

图1是本发明的流程框图。

图2(a)是采用2阶多项式表示密度随深度变化的盆地模型图。

图2(b)是盆地产生的负的重力异常值及反演误差图。

图2(c)是常规方法反演的密度分布图。

图2(d)是本发明方法(nz＝9)反演得到的地下密度分布图。

图3(a)矩形和平行四边形局部密度异常体模型图。

图3(b)是局部密度异常体模型产生的重力异常值及反演误差图。

图3(c)是常规方法反演的密度分布图。

图3(d)是本发明方法反演的密度分布图。

具体实施方式

为了使本发明目的、技术方案更加清楚明白，下面结合附图，对本发明作进一步详细说明。

该发明的详细技术操作简图如图1所示。该发明主要技术关键点为以下四个：(1)地下物质密度的多项式函数表示；(2)以多项式密度函数的系数为未知量的重力反演的核函数矩阵的构建；(3)约束条件和先验信息的施加；(4)建立重力反演线性方程组并求解密度函数的多项式系数。

实施例一：

如图2(a)所示，为一密度随深度变化的盆地模型，密度函数用2阶多项式表示δσ＝-400+274.8z-47.3z²，其中，z单位为km，δσ单位为kg/m³，浅层密度差大，在地表处最大达-0.4g/cm³。盆地的水平延伸长度为5.4km，盆地的最大基底深度为2.3km。图2(b)实线所示为此盆地产生的负的重力异常。根据本发明前述流程，进行重力反演：将研究区划分为60个竖向并排的矩形单元，每个矩形单元的密度分布均用随深度变化的9阶多项式表示，即nz＝9，构建重力核函数矩阵。取σmin＝-0.6g/cm³,σmax＝0.6g/cm³，σ0＝0，d＝300，权重λ1＝10^-4,λ2＝10^-3,λ3＝10^-3，计算出的各矩形单元的多项式系数汇总于表一。

表一:实施例一密度函数多项式系数计算结果表

图2(c)为假设密度线性变化，采用常规方法反演的密度分布，图2(d)为本发明方法，采用9阶多项式反演得到的地下密度分布。图中黑色线框代表密度异常体的真实位置，对比可见，本发明反演结果与真实密度吻合更好，呈现浅部密度差大、深部密度差小的趋势，且盆地基底与反演密度分布的整体轮廓相吻合，计算反演的密度分布与真实密度分布的平均绝对误差为0.042g/cm³。图2(b)中点线显示了图2(d)反演的密度分布计算的理论重力异常，可见观测重力异常与计算重力异常吻合较好，图2(b)中点划线显示了观测重力异常与理论计算重力异常的差值，两者的残差很小，在0.4mgal以下。

实施例二：

如图3(a)所示，为局部密度异常体模型，其中，矩形体为一正密度体，密度值为0.9g/cm³，平行四边形为一负密度体，密度值为-0.9g/cm³，两异常体的中心深度均在1.5km深处，负密度异常体随深度增加向正密度异常体一侧倾斜，图3(b)实线所示为此局部异常体产生的重力异常值。根据本发明前述流程，进行重力反演：将研究区划分为60个竖向并排的矩形单元，每个矩形单元的密度分布均用随深度变化的9阶多项式表示，即nz＝9，构建重力核函数矩阵。取σmin＝-1g/cm³,σmax＝1g/cm³，σ0＝0，d＝300，权重λ1＝10^-2,λ2＝10^-4,λ3＝10^-4，深度加权β＝2，z0＝0.5。计算出的各矩形单元的多项式系数汇总于表二。

表二：实施例二密度函数多项式系数计算结果表

图3(c)和图3(d)是采用常规平滑约束反演和采用本发明方法反演结果的对比，图中黑色线框代表密度异常体的真实位置，可以发现，常规平滑约束反演产生的密度分布更为平滑，虽然在异常体的中心显示密度差极值，但是，反演值与实际值相差很大，对局部异常体的边界刻画不清晰；而本发明方法通过反演高阶多项式的系数，可以有效的刻画局部异常体的真实密度分布和边界位置。图3(b)中点线显示了图3(d)反演的密度分布计算的理论重力异常，图3(b)中点划线显示了观测重力异常与理论计算重力异常的差值，两者的残差很小，在0.04mgal以下。

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