CN101285896B - 一种地球物理勘探中的重磁数据处理方法 - Google Patents

Abstract

一种地球物理勘探中的重磁数据处理方法，其包括步骤：在波数域将重磁资料频谱分别与任意方向的n阶导数或积分频率响应函数、延拓回返响应函数及差分法下延算子频率响应函数相乘，得到下延后的频谱；再经过快速傅氏反变换，获得向下延拓后的三维重磁数据体。利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。该方法不仅有更大的向下延拓深度并且有较好的抗噪能力，从而获得稳定的下半空间三维重磁资料分布图像。

Description

一种地球物理勘探中的重磁数据处理方法

技术领域

本发明涉及地球物理勘探技术领域，特别是涉及一种地球物理勘探中的重磁数据处理方法，可以获得地下空间重磁异常的分布，因而可以有效的对地下地质体进行精细描述。

背景技术

重力勘探与磁力勘探均属于地球物理勘探，重力勘探是利用重力仪在野外观测地下物质密度差异引起的重力异常以查明地下的地质构造和岩性异常体。磁法勘探则是利用磁力仪在野外观测地下物体磁性差异引起的磁力异常以查明地下的地质构造和磁性异常体如火山岩等。

重磁勘探是一个系统工程，可以分为三个环节：野外重磁资料采集、室内重磁资料处理以及室内重磁资料分析解释。

野外采集的重磁数据是不规则的数据，如当在地面上进行重磁勘探时，如果地形十分复杂，或者在某些点位上无法实测，这时实际测点的分布可能是不规则的。在航空磁测中由于航空磁测定位精度的提高，测线往往按实际航迹来恢复，这时实际测点的分布也是不规则的。然而，在对重磁测资料作数据处理时，总是要求数据按规则网格分布的。因此就需要由不规则网格上的实际重磁数据换算出规则网格节点上的重磁数据，这个过程就是数据网格化。显然，数据网格化的问题实际上是插值问题，用不规则分布的插值节点上的值来计算规则网格节点上的值，以适应数据处理的需要。

由于重磁及其导数都是解析函数，因此可以把它们从已知区域延拓到场源以外区域而仍然保持解析性，这就是重磁异常的向上延拓。但是，重磁异常的向下延拓不同于向上延拓，它是一个典型的不适定问题，重磁异常的向下延拓不能过场源，在接近场源处发散。

重磁异常的向上延拓可以利用积分公式比较精准的计算，但是却不能用来计算向下延拓，因为向下延拓的计算是对异常按指数关系放大。异常的“频率”越高，放大的越强烈。因此在利用积分公式进行向下延拓计算时，放大最厉害的是异常中包含的数据误差及一些“高频”干扰，结果异常曲线出现强烈的“振荡”现象，计算出的下延值无法利用。

为了解决这个问题，Henderson(1960)提出了一个有效的、得到普遍应用的方法。这个方法基于以下思路：地面以下、场源体之上重磁场值在很大程度上可以根据地面以上连续的、不同高度水平面上的重磁场值预测。因此，如果观测面上的重磁异常垂直向上延拓到四个连续的、相距一个格点距的水平面上，根据这四个向上延拓值及观测面上的重磁值，应用Lagrangian插值公式可以外推出向下延拓值。郑志毅(1974)提出了一个“积分插值”法解析延拓。这个方法与Henderson方法非常相似。该方法首先通过积分，计算出上延不同高度的重磁异常值，然后利用这些上延值与原观测面上的实测值，形成一个插值多项式，并外推出(外插)某个深度处的下延值。

向上延拓具有“低通滤波”的特性，它的主要作用是使重磁异常变得更为平滑，相对突出了区域异常的特征。向下延拓则是向上延拓的逆过程，具有“高通滤波”的特性，其作用是相对突出了重磁局部异常，分解在水平方向叠加的异常，定性估计场源的深度，以及由于下延使延拓面更接近场源，异常等值线圈闭的形状与场源体水平截面形状更为接近，因而可用来了解复杂异常源的平面轮廓。

下延计算属于不适定问题，引力位在场源体外和场源体内分别满足拉普拉斯方程和泊松方程，场源深度又属未知的；因而在理论上未能解决其计算方法，只能用插值公式进行外推。外推的深度越大，下延值的误差越大。由于下延计算的高通滤波性质，局部干扰和误差会被放大，使下延计算可能失败，因而在处理实际资料时向下延拓深度不能太大；而且每下延一次要对这个深度的下延值进行平滑处理，然后再继续下延。

进行上延计算时，由浅部场源体引起的“高频”异常随高度增加而衰减的同时，所求的由深部场源体引起的宽缓的“低频”异常也随高度衰减；进行下延计算时，由浅部场源体引起的“高频”异常随深度增加而增大的同时，由深部场源体引起的“低频”异常也得到放大。为了解决这个问题，Polowski(1995)提出了一个优选延拓方法，理论模型计算和实际应用已显示出它的效果。

前苏联学者B·M·别列兹金首先提出归一化总梯度法，并已较好地应用于油气勘探领域，国内一些学者肖一鸣、王家林等首先将该方法引入国内，并对此进行了改进研究，最近，肖鹏飞等又提出了稳定算法。徐世浙、陈生昌从不同的角度提出了向下延拓的方法。杨辉等提出重力异常二维、三维视深度滤波的方法，实现了下延过场源并定量确定矿体深度。

中国专利公开号CN1877366A，公开了一种重磁延拓回返垂直导数目标优化处理技术，原理是在观测面上进行资料的目标处理，获得观测面上的处理后的二维数据。步骤方法是通过延拓回返压制噪音，针对不同的勘探深度的目标，通过不同次数的垂直导数来提高分辨率。资料处理的目的是突出不同深度的勘探目标，实现资料的目标处理，与本方法适用的范围、目的及效果都不尽相同，本方法最终获得观测面及以下半空间的3D数据体。中国专利公开号CN1667433A，公开了一种位场高分辨率视深度滤波方法，该专利原理是将位场数据下延到地面以下研究地下不同深度的位场数据的分布，采用了水平一次导数以及垂直一次导数、垂直二次导数方法处理，以提高分辨率，目的是确定局部构造及断裂的分布及质心深度，该专利①每次处理只考虑了一个方向的导数；②导数的阶数只为1或为2；③没有考虑通过积分以增强区域异常；而本方法①综合考虑各个方向的导数增强局部异常；②导数的阶数n为有理数，由于导数的阶数可以连续变化，使其能够适应更为复杂的地质情况；③同时综合考虑各个方向的积分增强区域异常；其适用的范围、目的及效果都不尽相同。

虽然这些方法在重磁资料的应用中取得了好的效果，但是我们现在面对的地质问题更为复杂。因此，需要根据资料的观测精度以及针对不同的地质问题提出合适分辨率的方法以满足不同的地质需求。

发明内容

为了解决现有技术中存在更为复杂的地质问题以及为了满足不同的地质需求针对不同的地质问题提出合适分辨率的方法，本发明提供一种地球物理勘探中的重磁数据处理方法。

本发明的地球物理勘探中的重磁数据处理方法包括以下步骤：

(1)首先利用重力仪或磁力仪在野外采集重磁资料，或者在室内对原有的重磁异常图进行数字化而得到重磁资料，然后对野外采集或室内数字化的重磁资料进行各种常规校正及改正；

(2)利用Surfer软件中的Kriging插值方法将重磁资料网格化，形成规则网重磁资料f_a(x，y)，其中x、y分别为重磁资料的平面网格点坐标；

(3)进行镶边处理：利用低精度的重磁资料进行补充镶边；将规则网重磁资料f_a(x，y)转变为其线数和点数扩大为2的整数幂重磁资料f_b(x，y)；

(4)利用快速傅氏正变换(FFT)将镶边后的重磁资料f_b(x，y)变换到波数域，获得f_b(x，y)的频谱F_b(u，v)，F_b(u，v)＝FFT[f_b(x，y)]，其中u、v分别为x和y方向的波数；

(5)在波数域，将频谱F_b(u，v)与重磁场任意方向的n阶导数或积分频率响应函数 $q^n(u,v)={[2\mathrm{\πi}(\mathrm{\αu}+\mathrm{\βv})+\mathrm{\γ}\sqrt{u^2+v^2}]}^n$ 相乘，增强局部或区域重磁异常，

即：

F_c(u，v)＝F_b(u，v)×qⁿ(u，v)

其中n为有理数，i为虚数；α、β、γ分别为所求导数或积分方向的三个方向余弦；n为导数或积分的阶数，即分辨率参数；

(6)在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F_c(u，v)与延拓回返频率响应函数

相乘，抑制高频成分，使其满足向下延拓的条件，即：

其中h为向下延拓深度，p为延拓回返次数；h、p为正有理数；

(7)在波数域，在抑制高频成分的基础上，与差分法向下延拓频率响应函数 $d(u,v)=(6-\exp (-2\mathrm{\πh}\sqrt{u^2+v^2}-2\cos \left(2\mathrm{\πhu}\right)-2\cos \left(2\mathrm{\πhv}\right)))$ 相乘，得到向下延拓后的频谱F_e(u，v)

F_e(u，v)＝F_d(u，v)×d(u，v)；

(8)再经过快速傅氏反变换(IFFT)，获得向下延拓后的重磁资料f_e(x，y)，f_e(x，y)＝IFFT[F_e(u，v)]；

(9)进行反镶边处理，即从f_e(x，y)中提取原数据大小的向下延拓深度h后的重磁数据f_f(x，y)；

(10)重复步骤(6)，从而获得下半空间三维重磁数据体f_g(x，y，z)，其中，x，y，z为重磁异常的空间网格点坐标；

(11)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。

其中，-3≤n≤3。

其中，延拓回返次数p，0＜p≤30。

另外地，向下延拓深度h，0＜h≤500Km。

本发明较原有技术有较大的改进，使其完全符合油气勘探生产的需要并产生了好用及实用的效果。是一项新颖、进步、实用的新的重磁异常处理方法。使用该方法提高了重磁勘探确定局部异常及区域异常的可靠性，进而提高了油气构造及深层构造预测的准确性，节省了勘探成本，降低了勘探风险。

下面结合附图和具体实施例对本发明作进一步详细说明。

附图简要说明

图1是依据本发明的JY坳陷某地区向下延拓3Km重力异常立体图；

图2是JY坳陷某地区三维地震Tg₁-Tg构造平面图；

图3是JY坳陷某地区过CG18构造东西向地震剖面图；

图4是依据本发明的JY坳陷某地区向下延拓3Km航磁异常平面图；

图5是JY坳陷某地区过L151井地震剖面图。

具体实施方式

采用野外采集的重磁数据，使用本发明的方法，可以具体提高重磁异常的分辨率，以满足不同勘探目的的需要。

实施例1：参阅附图1、附图2、附图3。利用LCR型重力仪在JY坳陷某地区采集1/5万高精度重力资料。按照发明内容所述的11个步骤对重力资料进行处理，即：

(1)首先利用重力仪在野外采集重力资料，然后对野外采集的重力资料进行常规校正及改正；

(2)利用Surfer软件中的Kriging(克里金)插值方法将观测的重力资料网格化，形成规则网重力资料f_a(x，y)，其中，x，y为重力资料的平面网格点坐标；首先选择菜单命令“Grid(网格)|Data(数据)”，在“Open(打开)”对话框里选定数据文件后，单击“Open(打开)”，出现“离散数据插值”对话框；从“Gridding Method(网格化方法)”组中选择“Kriging(克里金插值法)”单击“Ok”即可；网格化数据大小为152行×210列；

(3)进行镶边处理。为了使实际重力资料的线数和点数满足快速傅氏变换的要求，将实际重力资料f_a(x，y)的线数和点数扩大为2的整数幂的重力资料f_b(x，y)。为了达到最好的镶边效果，利用1/20万低精度的重力资料进行补充镶边，这是一种最为合理、简单、实用的扩边方法。镶边后网格化数据大小为256行×256列；

(4)利用快速傅氏正变换(FFT)，将f_b(x，y)变换到波数域，获得实测重力资料的频谱F_b(u，v)，F_b(u，v)＝FFT[f_b(x，y)]，其中u，v分别为x和y方向的波数；

(5)在波数域，将镶边后的重力资料频谱F_b(u，v)与重力场任意方向的n阶导数或积分频率响应函数 $q^n(u,v)={[2\mathrm{\πi}(\mathrm{\αu}+\mathrm{\βv})+\mathrm{\γ}\sqrt{u^2+v^2}]}^n$ 相乘，增强局部或区域重磁异常，即：

F_c(u，v)＝F_b(u，v)×qⁿ(u，v)

其中i为虚数；α、β、γ分别为所求导数或积分方向的三个方向余弦；n为导数或积分的阶数，即分辨率参数，n为有理数；

(6)在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常Fc(u，v)与延拓回返频率响应函数

相乘，抑制高频成分，即：

其中h为向下延拓深度，p为延拓回返次数；

(7)在波数域，在抑制高频成分的基础上，与差分法向下延拓频率响应函数 $d(u,v)=(6-\exp (-2\mathrm{\πh}\sqrt{u^2+v^2}-2\cos \left(2\mathrm{\πhu}\right)-2\cos \left(2\mathrm{\πhv}\right)))$ 相乘，得到下延后的频谱F_e(u，v)

F_e(u，v)＝F_d(u，v)×d(u，v)；

(8)再经过快速傅氏反变换(IFFT)，获得向下延拓后的重力资料f_e(x，y)，f_e(x，y)＝IFFT[F_e(u，v)]；

(9)进行反镶边处理，即从f_e(x，y)中提取原数据大小的向下延拓深度h后的重力数据f_f(x，y)。f_f(x，y)数据大小为152行×210列；

(10)重复步骤(6)，从而获得下半空间三维重磁数据体f_g(x，y，z)，其中，x，y，z为重磁异常的空间网格点坐标；

(11)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。

图1为JY坳陷某地区利用本方法获得的重力资料向下延拓3Km立体图，图2为JY坳陷某地区利用三维地震资料经过处理解释后得到的Tg1-Tg构造图，图3为JY坳陷某地区过CG18构造东西向地震剖面图。从图1及图2中可以看到，地震构造走向与重力异常走向、高点位置一致，验证了重力资料的结果。利用处理后的重力资料发现了CG18异常(见图1)，地震资料构造图(图2)及地震剖面图(图3)证实该异常为古潜山，取得了明显的地质效果。

实施例2：参阅附图4、附图5。利用机载磁力仪在JY坳陷某地区采集1/5万高精度航磁资料。按照发明内容所述的11个步骤对航磁资料进行处理，即：

(1)首先利用磁力仪在野外采集航磁资料，然后对野外采集的航磁资料进行常规校正及改正；

(2)利用Surfer软件中的Kriging(克里金)插值方法将观测的航磁资料网格化，形成规则网航磁资料f_a(x，y)，其中，x，y为航磁资料的平面网格点坐标；首先选择菜单命令“Grid(网格)|Data(数据)”，在“Open(打开)”对话框里选定数据文件后，单击“Open(打开)”，出现“离散数据插值”对话框；从“Gridding Method(网格化方法)”组中选择“Kriging(克里金插值法)”单击“Ok”即可；网格化数据大小为692行×854列；

(3)进行镶边处理。为了使实际航磁资料的线数和点数满足快速傅氏变换的要求，将实际航磁资料f_a(x，y)的线数和点数扩大为2的整数幂的航磁资料f_b(x，y)。为了达到最好的镶边效果，利用1/20万低精度的航磁资料进行补充镶边，这是一种最为合理、简单、实用的扩边方法。镶边后网格化数据大小为1024行×1024列；

(4)利用快速傅氏正变换(FFT)，将f_b(x，y)变换到波数域，获得实测航磁资料的频谱F_b(u，v)，F_b(u，v)＝FFT[f_b(x，y)]，其中u，v分别为x和y方向的波数；

(5)在波数域，将镶边后的航磁资料频谱F_b(u，v)与航磁场任意方向的n阶导数或积分频率响应函数 $q^n(u,v)={[2\mathrm{\πi}(\mathrm{\αu}+\mathrm{\βv})+\mathrm{\γ}\sqrt{u^2+v^2}]}^n$ 相乘，增强局部或区域重磁异常，即：

F_c(u，v)＝F_b(u，v)×qⁿ(u，v)

其中i为虚数；α、β、γ分别为所求导数或积分方向的三个方向余弦；n为导数或积分的阶数，即分辨率参数，n为有理数；

(6)在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F_c(u，v)与延拓回返频率响应函数

相乘，抑制高频成分，即：

其中h为向下延拓深度，p为延拓回返次数；

(7)在波数域，在抑制高频成分的基础上，与差分法向下延拓频率响应函数 $d(u,v)=(6-\exp (-2\mathrm{\πh}\sqrt{u^2+v^2}-2\cos \left(2\mathrm{\πhu}\right)-2\cos \left(2\mathrm{\πhv}\right)))$ 相乘，得到下延后的频谱F_e(u，v)

F_e(u，v)＝F_d(u，v)×d(u，v)；

(8)再经过快速傅氏反变换(IFFT)，获得向下延拓后的航磁资料f_e(x，y)，f_e(x，y)＝IFFT[F_e(u，v)]；

(9)进行反镶边处理，即从f_e(x，y)中提取原数据大小的向下延拓深度h后的航磁数据f_f(x，y)。f_f(x，y)数据大小为692行×854列；

(10)重复步骤(6)，从而获得下半空间三维重磁数据体f_g(x，y，z)，其中，x，y，z为重磁异常的空间网格点坐标；

(11)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。

图4为JY坳陷某地区利用本方法获得的航磁资料向下延拓3Km平面图，图5为JY坳陷某地区过L151井地震剖面图，从图5中可以看到由于火成岩厚度较薄，在地震剖面上火成岩的地震反射与泥岩反射相同，只有一个强反射同相轴，难于识别区分，利用钻井标定处理后的航磁异常较好的预测了该区火山岩的分布，后得到数口井的验证。

上述离散数据网格化、镶边处理、反镶边处理的过程，本专业技术领域分析人员能熟练完成。

如上述，已经清楚详细地描述了本发明提出的地球物理勘探中的重磁数据处理方法。尽管本发明的优选实施例详细描述并解释了本发明，但是本领域普通的技术人员可以理解，在不背离所附权利要求定义的本发明的精神和范围的情况下，可以在形式和细节中做出多种修改。

Claims (4)

1.一种地球物理勘探中的重磁数据处理方法，该方法包括以下步骤：

(1)首先利用重力仪或磁力仪在野外采集重磁资料，或者在室内对原有的重磁异常图进行数字化而得到重磁资料，然后对野外采集或室内数字化的重磁资料进行各种常规校正及改正；

(2)利用Surfer软件中的Kriging插值方法将重磁资料网格化，形成规则网重磁资料f_a(x，y)，其中x、y分别为重磁资料的平面网格点坐标；

(3)进行镶边处理：利用低精度的重磁资料进行补充镶边；将规则网重磁资料f_a(x，y)转变为其线数和点数扩大为2的整数幂的重磁资料f_b(x，y)；

(4)利用快速傅氏正变换(FFT)将镶边后的重磁资料f_b(x，y)变换到波数域，获得f_b(x，y)的频谱F_b(u，v)，F_b(u，v)＝FFT[f_b(x，y)]，其中u、v分别为x和y方向的波数；

(5)在波数域，将频谱F_b(u，v)与重磁场任意方向的n阶导数或积分频率响应函数相乘，增强局部或区域重磁异常，即：

F_c(u，v)＝F_b(u，v)×qⁿ(u，v)

其中n为有理数，i为虚数；α、β、γ分别为所求导数或积分方向的三个方向余弦；n为导数或积分的阶数，即分辨率参数；

(6)在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F_c(u，v)与延拓回返频率响应函数

相乘，抑制高频成分，使其满足向下延拓的条件，即：

其中h为向下延拓深度，p为延拓回返次数；h、p为正有理数；

(7)在波数域，在抑制高频成分的基础上，与差分法向下延拓频率响应函数

相乘，得到向下延拓后的频谱F_e(u，v)

F_e(u，v)＝F_d(u，v)×d(u，v)；

(8)再经过快速傅氏反变换(IFFT)，获得向下延拓后的重磁资料f_e(x，y)，f_e(x，y)＝IFFT[F_e(u，v)]；

(9)进行反镶边处理，即从f_e(x，y)中提取原数据大小的向下延拓深度h后的重磁数据f_f(x，y)；

(10)重复步骤(6)，从而获得下半空间三维重磁数据体f_g(x，y，z)，其中，x，y，z为重磁异常的空间网格点坐标；

(11)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。

2.根据权利要求1所述的地球物理勘探中的重磁数据处理方法，其特征是：所述的n阶导数或积分，-3≤n≤3。

3.根据权利要求1所述的地球物理勘探中的重磁数据处理方法，其特征是：所述的延拓回返次数p，0＜p≤30。

4.根据权利要求1所述的地球物理勘探中的重磁数据处理方法，其特征是：所述向下延拓深度h，0＜h≤500Km。

Images