一种起伏地层的频率域快速高精度正演方法与流程

本发明涉及地球物理重力勘探技术领域，具体涉及一种起在伏地层条件下重力位场的频率域正演方法；该方法通过将场源体的上、下界面分别在其对应的平均水平面上泰勒展开，并在傅里叶计算时引入高斯算法，显著提升了正演计算的数值精度。

背景技术：

重力勘探作为一项最基本的地球物理方法，已被广泛运用于地形校正、水文地质勘察、大地水准面测量、固体矿产资源勘探、岩石圈壳幔结构研究等技术领域，通过对重力勘探方法的研究可以开展相关领域的地质工作，例如对划分地球内部圈层结构进行、对区域构造的深入研究、以及对储油构造等地下矿藏的探明等。重力勘探相比于其他地球物理勘探方在勘查成本、勘查深度以及获取重力场信息效率方面具有显著优势，同时，重力勘探对于横向不均匀地质体的勘测，特别是在划分异常体边界方面也具有独特的优势。综上所述，用于研究地球内部构造的密度界面反演方法，一直以来都是属于重力勘探的重点研究内容；而反演的可行性及可信度又是由正演的速度和精度所决定的，因此在一个高效的反演程序中，快速高精度正演方法占有极其重要的地位。

在用重力勘探方法进行研究时首先需要对重力数据进行校正，为了消除地形起伏产生的影响，在得到地形高程数据和密度信息之后需要准确而快速的正演出地形产生的重力异常，完成地形改正。因为重力勘探深度大，加之重力卫星技术的快速发展，世界各国的科学家们热衷于使用卫星重力数据研究地球的内部构造和圈层划分，而通过卫星重力数据研究内部构造时需要对地壳密度进行校正以消除密度不均造成的影响，在获取地壳各层密度分布和地壳各层上、下界面的深度信息时需要一个准确而高效的正演方法作为支撑，以计算出地壳密度不均产生的重力异常。

现有的重力正演方法按照计算域可以分为空间域和频率域两大类，其中，常规的地形校正或者地壳密度校正通常采用空间域正演方法，虽然该方法正演精度高，但是正演计算式特别复杂，空间域中场源体的几何参数与正演公式对应繁琐，计算效率非常低。而频率域方法由于采用了快速傅里叶变换，正演计算式大为精简，异常谱和场源体的几何特性对应很简单，使得计算效率得到显著提升；但是存在的不足是相比于空间域方法其在计算精度上大打折扣，这是由于快速傅里叶算法是采用梯形法则计算连续傅里叶变换的数值方法，属于连续型积分计算，且连续傅里叶变换属于震荡积分问题，故传统的快速傅里叶数值算法不可避免会产生边界震荡效应、强制周期化等一系列频率域正演问题，影响频率域正演的数值精度。目前，针对频率域计算精度的常用解决方法是将积分区间离散化，使得连续型积分计算转换为矩形积分计算，但仍存在进一步优化和提高的空间。

技术实现要素：

本发明的目的在于提供一种兼顾高计算效率和高计算精度的用于起伏地层频率域计算的正演方法，以解决背景技术中提出的问题。

为实现上述目的，本发明提供了一种起伏地层频率域快速高精度正演方法，包括如下步骤，

步骤1)建立关于起伏地层的场源体模型，根据parker正演推导公式，可得到该场源体模型在高度为z0的水平观测面上产生的重力位的频率域表达式：

其中，u表示重力位，g为万有引力常量，ρ为场源体密度，h1和h2分别代表场源体上、下界面的高度，k表示频率域波数：

步骤2)将场源体的上、下界面分别在其对应的平均水平面l1和l2上按泰勒级数展开，使得上、下界面的绝对起伏值转换为相对(平均水平面的)起伏量，则公式(1)中代表起伏面的两项指数项可以写为：

步骤3)将公式(2)和(3)代入公式(1)中并对重力位u求z向偏导，得到起伏地层模型在水平观测面z0上产生的重力异常的频率域表达式：

步骤4)在对公式(4)做傅里叶变换计算时，将积分区间离散化，使得连续型积分计算转换为矩形积分计算，具体说明如下，

二维傅里叶反变换可以写为：

将公式(5)在频率域离散且x、y方向离散点的数量分别为m个和n个，离散后的表达式可写为：

还可对于上述方法做进一步优选，在步骤4)的傅里叶变换计算中引入高斯节点，将积分区间细分，使得原来的矩形积分转换为高斯型积分，具体说明如下，

函数g(x)在[a，b]区间的一维高斯积分可以写为：

其中，k为高斯节点数，ak和tk分别表示高斯系数和区间[-1，1]上的节点，通过公式(7)可将积分区间由[a，b]转换到[-1，1]以实现变换积分的高斯积分；

将公式(7)的变换运用到公式(6)中，得到二维高斯型积分公式：

其中，(θix，αix)和(θiy，αiy)分别表示高斯系数和区间[-1，1]上的节点，因此和是在区间[0，1]上的高斯节点，对应的高斯系数分别为和

将公式(8)化简之后得：

由于公式(9)的计算精度与高斯节点的数量成正比，因此，可通过增加高斯节点的数量使公式(9)中的“≈”向“＝”逼近。

高斯节点数不同，其对应的高斯节点和高斯系数也不同，其中，高斯节点和高斯系数均可通过查表获得，表1中给出了若干个关于高斯节点及其对应系数的例子。

表1高斯型求积公式节点及系数表

从表1中可以看出，随着高斯节点数量的增加，其对应求积时的余项急剧减小，因此，节点数量越多则精度越高。即当高斯节点足够多时，公式(9)中的“≈”可视为“＝”。

本发明提供的技术方案至少具有如下有益效果：

1、该方法提供了一个更为稳定的正演计算式，通过将场源体的上、下界面分别在其对应的平均水平面上泰勒展开，使得计算时所用的场源体模型上、下界面的绝对深度数据转化为相对平均界面的起伏量数据，提高正演算法的精度。

2、该方法在傅里叶变换时通过引入高斯节点将矩形积分转换成了精度更高的高斯型积分，在保证计算效率的同时兼顾了数值精度，从而大幅度提高了频率域正演方法的可靠性，确保快速精准地完成地形改正和密度校正等目标；同时还可以根据正演精度的需要调整高斯节点的个数，高斯节点的个数越多则正演结果的精度越高。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图，其中：

图1是本发明中场源体模型和观测面的结构示意图；

图2是场源体上、下界面分别在其对应的平均水平面上展开的示意图；

图3是采用空间域方法得到的正演结果示意图(在本发明中视为理论值)；

图4是采用对比实施例1中传统的频率域方法得到的正演结果，且图4(a)和图4(b)分别为正演结果示意图以及与理论值对比后的误差示意图；

图5是采用实施例中本发明改进后的频率域方法得到的正演结果，且图5(a)和图5(b)分别为正演结果示意图以及与理论值对比后的误差示意图；

其中：1观测面，2上界面，3下界面。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

对比实施例1(采用传统的频率域正演方法)：

步骤1)建立关于起伏地层的场源体模型，根据parker正演推导公式，可得到该场源体模型在高度为z0的水平观测面上产生的重力位的频率域表达式：

其中，u表示重力位，g为万有引力常量，ρ为场源体密度，h1和h2分别代表场源体上、下界面的高度，k表示频率域波数：

步骤2)在对公式(1)做傅里叶变换计算时，将积分区间离散化，使得连续型积分计算转换为矩形积分计算，具体说明如下，

二维傅里叶反变换可以写为：

将公式(5)在频率域离散且x、y方向离散点的数量分别为m个和n个，离散后的表达式可写为：

实施例1(采用本发明的频率域正演方法)：

步骤1)建立关于起伏地层的场源体模型，根据parker正演推导公式，可得到该场源体模型在高度为z0的水平观测面上产生的重力位的频率域表达式：

其中，u表示重力位，g为万有引力常量，ρ为场源体密度，h1和h2分别代表场源体上、下界面的高度，k表示频率域波数：

步骤2)如图2所示，将场源体的上、下界面分别在其对应的平均水平面l1和l2上按泰勒级数展开，使得上、下界面的绝对起伏值转换为相对起伏量，则公式(1)中代表起伏面的两项指数项可以写为：

步骤3)将公式(2)和(3)代入公式(1)中并对重力位u求z向偏导，得到起伏地层模型在水平观测面z0上产生的重力异常的频率域表达式：

步骤4)在对公式(4)做傅里叶变换计算时，将积分区间离散化，使得连续型积分计算转换为矩形积分计算，具体说明如下，

二维傅里叶反变换可以写为：

将公式(5)在频率域离散且x、y方向离散点的数量分别为m个和n个，离散后的表达式可写为：

继续引入高斯节点，将积分区间细分，使得原来的矩形积分转换为高斯型积分，具体说明如下，

函数g(x)在[a，b]区间的一维高斯积分可以写为：

其中，k为高斯节点数，ak和tk分别表示高斯系数和区间[-1，1]上的节点，通过公式(7)可将积分区间由[a，b]转换到[-1，1]以实现变换积分的高斯积分；

将公式(7)的变换运用到公式(6)中，得到二维高斯型积分公式：

其中，(θix，αix)和(θiy，αiy)分别表示高斯系数和区间[-1，1]上的节点，因此和是在区间[0，1]上的高斯节点，对应的高斯系数分别为和

将公式(8)化简之后得：

由于公式(9)的计算精度与高斯节点的数量成正比，因此，可通过增加高斯节点的数量提高公式(9)的精度，当高斯节点足够多时，公式(9)中的“≈”可视为“＝”。

本发明方法用于正演模拟时，先在程序中设置好相关的正演参数，包括场源体上下界面数据、场源体密度数据和观测水平面高度等，然后根据界面起伏值确定界面深度的平均值(即对应的平均水平面)，根据界面起伏程度选取泰勒级数展开项数，若界面起伏较大则适当增加泰勒级数展开项数，一般来说泰勒级数的展开项数为6～10项。由于高斯节点和剖分参数越大则计算精度越高，但相应地计算时间也会随之增加，因此需根据正演模拟的精度要求选择合适的高斯节点的个数以及剖分参数nx和ny的大小。

为评价本发明方法的计算精度和速度，我们建立了如图1所示的模型：模型体的上下界面均为起伏面，上界面的平均深度l1＝-1km、下界面的平均深度l2＝-3km，x和y方向的长度均为256km，即0≤ξ≤256km、0≤η≤256km，模型密度为ρ＝800kg/m³，观测面为z0＝1.5km的水平面，x轴方向和y轴方向上的观测点数nx＝ny＝256，网格间距dx＝dy＝1km。

首先，采用1966年nagy提出的空间域正演算法对模型体在观测面上产生的重力异常进行正演计算，该算法的思路是将模型体剖分成多个棱柱体，分别计算单个棱柱体在观测点上的重力异常，然后累加求和，得到的正演结果如图3所示。由于棱柱体在空间某点处产生的重力异常具有解析解，因此空间域算法具有很高的精度，我们将该正演结果作为理论值。

然后，参照对比实施例1，我们采用传统的parker公式(1)对模型体进行正演计算并在做傅里叶变换时引入公式(6)，将所得结果与理论值作差，正演结果与误差分别如图4(a)和图4(b)所示；参照实施例1，我们采用优化后的公式(4)对模型体进行正演计算，同时在做傅里叶变换时引入公式(9)，采用4个高斯节点，泰勒展开项数为10项，将所得结果与理论值作差，正演结果与误差分别如图5(a)和图5(b)所示。

对比图4(b)和图5(b)，可知采用本发明改进后的算法可以大幅提高正演计算精度，这有力地证明了本发明方法在正演精度上的可靠性。另外，在耗时方面，空间域方法的正演时间为2165秒，而本发明改进后算法的正演时间仅需16秒，大大提高了计算效率。

以上所述仅为本发明的优选实施例，并非因此限制本发明的专利保护范围，对于本领域的技术人员来说，本发明可以有各种更改和变化。在本发明的精神和原则之内，凡是利用本发明说明书及附图内容所作的任何改进或等同替换，直接或间接运用在其它相关的技术领域，均应包括在本发明的专利保护范围内。