专利名称：一种地球物理勘探中的重磁数据处理方法
技术领域：
本发明涉及地球物理勘探技术领域，特别是涉及一种地球物理勘探中的重 磁数据处理方法，可以获得地下空间重磁异常的分布，因而可以有效的对地下 地质体进行精细描述。
背景技术：
重力勘探与磁力勘探均属于地球物理勘探，重力勘探是利用重力仪在野外 观测地下物质密度差异引起的重力异常以査明地下的地质构造和岩性异常体。 磁法勘探则是利用磁力仪在野外观测地下物体磁性差异引起的磁力异常以査 明地下的地质构造和磁性异常体如火山岩等。
重磁勘探是一个系统工程，可以分为三个环节野外重磁资料采集、室内 重磁资料处理以及室内重磁资料分析解释。
野外采集的重磁数据是不规则的数据，如当在地面上进行重磁勘探时，如 果地形十分复杂，或者在某些点位上无法实测，这时实际测点的分布可能是不 规则的。在航空磁测中由于航空磁测定位精度的提高，测线往往按实际航迹来 恢复，这时实际测点的分布也是不规则的。然而，在对重磁测资料作数据处理 时，总是要求数据按规则网格分布的。因此就需要由不规则网格上的实际重磁 数据换算出规则网格节点上的重磁数据，这个过程就是数据网格化。显然，数 据网格化的问题实际上是插值问题，用不规则分布的插值节点上的值来计算规 则网格节点上的值，以适应数据处理的需要。
由于重磁及其导数都是解析函数，因此可以把它们从已知区域延拓到场源 以外区域而仍然保持解析性，这就是重磁异常的向上延拓。但是，重磁异常的 向下延拓不同于向上延拓，它是一个典型的不适定问题，重磁异常的向下延拓不能过场源，在接近场源处发散。
重磁异常的向上延拓可以利用积分公式比较精准的计算，但是却不能用来 计算向下延拓，因为向下延拓的计算是对异常按指数关系放大。异常的"频率" 越高，放大的越强烈。因此在利用积分公式进行向下延拓计算时，放大最厉害 的是异常中包含的数据误差及一些"高频"干扰，结果异常曲线出现强烈的"振 荡"现象，计算出的下延值无法利用。
为了解决这个问题，Henderson (I960)提出了一个有效的、得到普遍应用 的方法。这个方法基于以下思路地面以下、场源体之上重磁场值在很大程度 上可以根据地面以上连续的、不同高度水平面上的重磁场值预测。因此，如 果观测面上的重磁异常垂直向上延拓到四个连续的、相距一个格点距的水平面 上，根据这四个向上延拓值及观测面上的重磁值，应用Lagmngian插值公式可 以外推出向下延拓值。郑志毅(1974)提出了一个"积分插值"法解析延拓。 这个方法与Henderson方法非常相似。该方法首先通过积分，计算出上延不 同高度的重磁异常值，然后利用这些上延值与原观测面上的实测值，形成一个 插值多项式，并外推出(外插)某个深度处的下延值。
向上延拓具有"低通滤波"的特性，它的主要作用是使重磁异常变得更 为平滑，相对突出了区域异常的特征。向下延拓则是向上延拓的逆过程，具有 "高通滤波"的特性，其作用是相对突出了重磁局部异常，分解在水平方向叠 加的异常，定性估计场源的深度，以及由于下延使延拓面更接近场源，异常等 值线圈闭的形状与场源体水平截面形状更为接近，因而可用来了解复杂异常源 的平面轮廓。
下延计算属于不适定问题，引力位在场源体外和场源体内分别满足拉普拉 斯方程和泊松方程，场源深度又属未知的；因而在理论上未能解决其计算方法，只能用插值公式进行外推。外推的深度越大，下延值的误差越大。由于下延计 算的高通滤波性质，局部干扰和误差会被放大，使下延计算可能失败，因而在
处理实际资料时下延深度不能太大；而且每下延一次要对这个深度的下延值进 行平滑处理，然后再继续下延。
进行上延计算时，由浅部场源体引起的"高频"异常随高度增加而衰减的 同时，所求的由深部场源体引起的宽缓的"低频"异常也随高度衰减；进行下 延计算时，由浅部场源体引起的"高频"异常随深度增加而增大的同时，由 深部场源体引起的"低频"异常也得到放大。为了解决这个问题，Polowski (1995) 提出了一个优选延拓方法，理论模型计算和实际应用已显示出它的效果。
前苏联学者B'M 别列兹金首先提出归一化总梯度法，并己较好地应 用于油气勘探领域，国内一些学者肖一鸣、王家林等首先将该方法引入国内， 并对此进行了改进研究，最近，肖鹏飞等又提出了稳定算法。徐世浙、陈生昌 从不同的角度提出了向下延拓的方法。杨辉等提出重力异常二维、三维视深度 滤波的方法，实现了下延过场源并定量确定矿体深度。
中国专利公开号CN1877366A，公开了一种重磁延拓回返垂直导数目标优 化处理技术，原理是在观测面上进行资料的目标处理，获得观测面上的处理后 的二维数据。步骤方法是通过延拓回返压制噪音，针对不同的勘探深度的目标， 通过不同次数的垂直导数来提高分辨率。资料处理的目的是突出不同深度的勘 探目标，实现资料的目标处理，与本方法适用的范围、目的及效果都不尽相同， 本方法最终获得观测面及以下半空间的3D数据体。中国专利公开号 CN1667433A，公开了一种位场高分辨率视深度滤波方法，该专利原理是将位场 数据下延到地面以下研究地下不同深度的位场数据的分布，采用了水平一次导 数以及垂直一次导数、垂直二次导数方法处理，以提高分辨率，目的是确定局部构造及断裂的分布及质心深度，该专利①每次处理只考虑了一个方向的导 数；②导数的阶数只为1或为2;③没有考虑通过积分以增强区域异常；而本 方法①综合考虑各个方向的导数增强局部异常；②导数的阶数n为有理数，由 于导数的阶数可以连续变化，使其能够适应更为复杂的地质情况；③同时综合 考虑各个方向的积分增强区域异常；其适用的范围、目的及效果都不尽相同。 虽然这些方法在重磁资料的应用中取得了好的效果，但是我们现在面对的 地质问题更为复杂。因此，需要根据资料的观测精度以及针对不同的地质问题 提出合适分辨率的方法以满足不同的地质需求。
发明内容
为了解决现有技术中存在更为复杂的地质问题以及为了满足不同的地质 需求针对不同的地质问题提出合适分辨率的方法，本发明提供一种地球物理勘 探中的重磁数据处理方法。
本发明的地球物理勘探中的重磁数据处理方法包括以下步骤
(1) 首先利用重力仪或磁力仪在野外采集重磁资料，或在室内对重磁异常图 进行数字化，然后对野外采集或室内数字化的重磁资料进行各种常规校正及改 正；
(2) 利用Surfer中的Kriging插值方法将重磁资料网格化，形成规则网重 磁资料fa(x， y)，其中，x， y为重磁资料的平面网格点坐标；
(3) 进行镶边处理。为了使实际重磁资料的线数和点数满足快速傅氏变换的 要求，将实际重磁资料fa(x， y)的线数和点数扩大为2的整数幂的重磁资料 fb (x， y)。为了达到最好的镶边效果，利用低精度的重磁资料进行补充镶边， 这是一种最为合理、简单、实用的扩边方法；
(4) 利用快速傅氏正变换(FFT)，将fb (x， y)变换到波数域，获得实测重磁资料的频谱Fb (u， v)， Fb (u， v) =FFT [fb (x， y)]，其中u， v分别为x 和y方向的波数；
(5) 在波数域，将镶边后的重磁资料频谱Fb (u， v)与重磁场任意方向的n
阶导数或积分频率响应函数qn (u， v) =[2力(咖+ /^) + ^/"2+一]"相乘，增强局 部或区域重磁异常，艮P:
F。 (u， v) 二Fb (u， v) X qn (u， v)
其中i为虚数；a、 P、 Y分别为所求导数或积分方向的三个方向余弦； n为导数或积分的阶数，即分辨率参数；
(6) 在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F。 (u， v)与延拓回返频域响应函数
—,v) = [exp(-2;zWM2 +v2 )]p(6_exp(_2;zWM2 + v2 - 2cos(2;z^)-2cos(2;^v)]p相乘，抑
制高频成分，使其满足向下延拓的条件，艮口 Fd (u， v)(u， v) X —,v) 其中h为下延深度，p为延拓回返次数；
(7) 在波数域，在抑制高频成分的基础上，与差分法向下延拓频域响应函数 d (u， v) =(6-exp(—2;^Vm2 +v2 —2cos(2油m) —2cos(2Wv)))相乘，《寻到下延后的频 谱Fe (u， v)
Fe (u， v) =Fd (u， v) X d (u， v);
(8) 再经过快速傅氏反变换(IFFT)，获得向下延拓后的重磁资料f;(x， y)， fe (x， y) 二IFFT [Fe (u， v)];
(9) 进行反镶边处理，即从fs (x， y)中提取原数据大小的向下延拓深度h 后的重磁数据ff (x， y);
(10) 重复步骤(6)，从而获得下半空间三维重磁数据体fg (x， y， z)，其中，x， y， Z为重磁异常的空间网格点坐标；
(ll)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对 比、分析并提交下一步勘探的有利区带及目标。
其中n、 p、 h均为有理数，并且—3《n《3、 0《p《30、 0《h《500Km (千米)。
本发明较原有技术有较大的改进，使其完全符合油气勘探生产的需要并产 生了好用及实用的效果。是一项新颖、进步、实用的新的重磁异常处理方法。 使用该方法提高了重磁勘探确定局部异常及区域异常的可靠性，进而提高了油 气构造及深层构造预测的准确性，节省了勘探成本，降低了勘探风险。
下面结合附图和具体实施例对本发明作进一步详细说明。
附图简要说明
图1是依据本发明的JY坳陷某地区向下延拓3Km重力异常立体图2是JY坳陷某地区三维地震Tgl-Tg构造平面图3是JY坳陷某地区过CG18构造东西向地震剖面图4是依据本发明的JY坳陷某地区向下延拓3Km航磁异常平面图5是JY坳陷某地区过L151井地震剖面图。
具体实施例方式
采用野外采集的重磁数据，使用本发明的方法，可以具体提高重磁异常的 分辨率，以满足不同勘探目的的需要。
实施例1:参阅附图1、附图2、附图3。利用LCR型重力仪在JY拗陷某 地区采集1 / 5万高精度重力资料。按照发明内容所述的11个步骤对重力资料 进行处理，艮口
(l)首先利用重力仪在野外采集重力资料，然后对野外采集的重力资料进行 常规校正及改正；(2) 利用Surfer软件中的Kriging (克里金)插值方法将观测的重力资料 网格化，形成规则网重力资料fa (x， y)，其中，x， y为重力资料的平面网格 点坐标；首先选择菜单命令"Grid (网格)I Data (数据)"，在"Open (打开)" 对话框里选定数据文件后，单击"Open (打开)"，出现"离散数据插值"对话 框；从"Gridding Method (网格化方法)"组中选择"Kriging (克里金插值 法)"单击"Ok"即可；网格化数据大小为152行x210列；
(3) 进行镶边处理。为了使实际重力资料的线数和点数满足快速傅氏变换的 要求，将实际重力资料fa (x， y)的线数和点数扩大为2的整数幂的重力资料 fb (x， y)。为了达到最好的镶边效果，禾U用1/20万低精度的重力资料进行 补充镶边，这是一种最为合理、简单、实用的扩边方法。镶边后网格化数据大 小为256行x 256列；
(4) 利用快速傅氏正变换(FFT)，将fb (x， y)变换到波数域，获得实测重 力资料的频谱Fb (u， v)， Fb (u， v) =FFT [fb (x， y)]，其中u， v分别为x 和y方向的波数；
(5) 在波数域，将镶边后的重力资料频谱Fb (u， v)与重力场任意方向的n 阶导数或积分频率响应函数qn (u， v) =[2范'(^ +側+ ^/"2+一]"相乘，增强局 部或区域重磁异常，艮P:
F。 (u， v) =Fb (u， v) X q11 (u， v)
其中i为虚数；a、 P、 Y分别为所求导数或积分方向的三个方向余弦； n为导数或积分的阶数，即分辨率参数，n为有理数；
(6) 在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F。 (u， v)与延拓回返频域响应函数
pO,v) = [exp(-2;rWM2 +v2)]p(6 _ exp(-2;zWw2 +v2 — 2cos(2Wm) — 2cos(2;^v)]p相乘，抑制高频成分，艮P:
Fd (u， v)(u， v) X —,v) 其中h为下延深度，p为延拓回返次数；
(7) 在波数域，在抑制高频成分的基础上，与差分法向下延拓频域响应函数 d (u， v) =(6 — exp(—2^Vm2 + v2 — 2cos(2油m) — 2cos(2Wv)))相乘，得到下延后的频 谱Fe (u， v)
Fe (u， v) =Fd (u， v) X d (u， v);
(8) 再经过快速傅氏反变换(IFFT)，获得向下延拓后的重力资料fs(x， y)， fe (x， y) =IFFT [Fe (u， v)];
(9) 进行反镶边处理，即从f^ (x， y)中提取原数据大小的向下延拓深度h 后的重力数据ff (x， y)。 ff (x， y)数据大小为152行X210列；
(抑重复步骤(6)，从而获得下半空间三维重磁数据体fg (x， y， z)，其中， x， y， z为重磁异常的空间网格点坐标；
(ll)利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对 比、分析并提交下一步勘探的有利区带及目标。
图1为JY坳陷某地区利用本方法获得的重力资料向下延拓3Km立体图， 图2为JY坳陷某地区利用三维地震资料经过处理解释后得到的Tgl-Tg构造 图，图3为JY坳陷某地区过CG18构造东西向地震剖面图。从图1及图2中可 以看到，地震构造走向与重力异常走向、高点位置一致，验证了重力资料的结 果。利用处理后的重力资料发现了CG18异常(见图l)，地震资料构造图(图 2)及地震剖面图(图3)证实该异常为古潜山，取得了明显的地质效果。
实施例2:参阅附图4、附图5。利用机载磁力仪在JY坳陷某地区采集1 / 5万高精度航磁资料。按照发明内容所述的11个步骤对航磁资料进行处理，即
(1) 首先利用磁力仪在野外采集航磁资料，然后对野外采集的航磁资料进行
常规校正及改正；
(2) 利用Surfer软件中的Kriging (克里金)插值方法将观测的航磁资料 网格化，形成规则网航磁资料fa (x， y)，其中，x， y为航磁资料的平面网格 点坐标；首先选择菜单命令"Grid (网格)I Data (数据)"，在"Open (打开)" 对话框里选定数据文件后，单击"Open (打开)"，出现"离散数据插值"对话 框；从"Gridding Method (网格化方法)"组中选择"Kriging (克里金插值 法)"单击"Ok"即可；网格化数据大小为692行x 854列；
(3) 进行镶边处理。为了使实际航磁资料的线数和点数满足快速傅氏变换的 要求，将实际航磁资料fa (x， y)的线数和点数扩大为2的整数幂的航磁资料 fb (x， y)。为了达到最好的镶边效果，利用1/20万低精度的航磁资料进行 补充镶边，这是一种最为合理、简单、实用的扩边方法。镶边后网格化数据大 小为1024行x 1024列；
(4) 利用快速傅氏正变换(FFT)，将fb (x， y)变换到波数域，获得实测航 磁资料的频谱Fb (u， v)， Fb (u， v) =FFT [fb (x， y)]，其中u， v分别为x 和y方向的波数；
(5) 在波数域，将镶边后的航磁资料频谱Fb (u， v)与航磁场任意方向的n 阶导数或积分频率响应函数q11 (u， v) = [2瓜'(6m + + n/"2 + v2]"相乘，增强局 部或区域重磁异常，艮P:
Fc (u， v) =Fb (u， v) X q11 (u， v)
其中i为虚数；a、 0、 Y分别为所求导数或积分方向的三个方向余弦； n为导数或积分的阶数，即分辨率参数，n为有理数；(6) 在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常F。
(U， V)与延拓回返频域响应函数 <formula>formula see original document page 14</formula>相乘，抑 制高频成分，艮口
<formula>formula see original document page 14</formula> 其中h为下延深度，p为延拓回返次数；
(7) 在波数域，在抑制高频成分的基础上，与差分法向下延拓频域响应函数 d (u， v) =(6-exp(-2油V"2 + v2 — 2cos(2油m)-2cos(2^v)》相乘，f导到下延后的步贞 谱Fe (U， v)<formula>formula see original document page 14</formula>
(8) 再经过快速傅氏反变换(IFFT)，获得向下延拓后的航磁资料<formula>formula see original document page 14</formula>
(9) 进行反镶边处理，即从L (x， y)中提取原数据大小的向下延拓深度h 后的航磁数据ff (x， y)。 ff (x, y)数据大小为692行x 854列；
(10) 重复步骤(6)，从而获得下半空间三维重磁数据体fg (x， y， z)，其中，
x， y， z为重磁异常的空间网格点坐标；
(11) 利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对 比、分析并提交下一步勘探的有利区带及目标。
图4为JY坳陷某地区利用本方法获得的航磁资料向下延拓3Km平面图， 图5为JY坳陷某地区过L151井地震剖面图，从图5中可以看到由于火成岩厚 度较薄，在地震剖面上火成岩的地震反射与泥岩反射相同，只有一个强反射同 相轴，难于识别区分，利用钻井标定处理后的航磁异常较好的预测了该区火山 岩的分布，后得到数口井的验证。
上述离散数据网格化、镶边处理、反镶边处理的过程，本专业技术领域分 析人员能熟练完成。如上述，已经清楚详细地描述了本发明提出的地球物理勘探中的重磁数据 处理方法。尽管本发明的优选实施例详细描述并解释了本发明，但是本领域普 通的技术人员可以理解，在不背离所附权利要求定义的本发明的精神和范围的 情况下，可以在形式和细节中做出多种修改。
权利要求
1. 一种地球物理勘探中的重磁数据处理方法，该方法包括以下步骤(1)首先利用重力仪或磁力仪在野外采集重磁资料，或者在室内对原有的重磁异常图进行数字化而得到重磁资料，然后对野外采集或室内数字化的重磁资料进行各种常规校正及改正；(2)利用Surfer中的Kriging插值方法将重磁资料网格化，形成规则网重磁资料fa(x，y)，其中x、y分别为重磁资料的平面网格点坐标；(3)进行镶边处理将实际重磁资料fa(x，y)的线数和点数扩大为2的整数幂的重磁资料fb(x，y)；利用低精度的重磁资料进行补充镶边；(4)利用快速傅氏正变换(FFT)将fb(x，y)变换到波数域，获得实测重磁资料的频谱Fb(u，v)，Fb(u，v)＝FFT[fb(x，y)]，其中u、v分别为x和y方向的波数；(5)在波数域，将镶边后的重磁资料频谱Fb(u，v)与重磁场任意方向的n阶导数或积分频率响应函数<math-cwu><![CDATA[<math> <mrow><msup> <mi>q</mi> <mi>n</mi></msup><mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo></mrow><mo>=</mo><msup> <mrow><mo>[</mo><mn>2</mn><mi>πi</mi><mrow> <mo>(</mo> <mi>αu</mi> <mo>+</mo> <mi>βv</mi> <mo>)</mo></mrow><mo>+</mo><mi>γ</mi><msqrt> <msup><mi>u</mi><mn>2</mn> </msup> <mo>+</mo> <msup><mi>v</mi><mn>2</mn> </msup></msqrt><mo>]</mo> </mrow> <mi>n</mi></msup> </mrow></math>]]></math-cwu><!--img id="icf0001" file="S2008101105005C00011.gif" wi="75" he="5" top= "164" left = "99" img-content="drawing" img-format="tif" orientation="portrait" inline="no"/-->相乘，增强局部或区域重磁异常，即Fc(u，v)＝Fb(u，v)×qn(u，v)其中n为有理数，i为虚数；α、β、γ分别为所求导数或积分方向的三个方向余弦；n为导数或积分的阶数，即分辨率参数；(6)在波数域，根据向下延拓深度h，将增强后的局部或区域重磁异常Fc(u，v)与延拓回返频域响应函数id="icf0002" file="S2008101105005C00012.gif" wi="143" he="5" top= "242" left = "29" img-content="drawing" img-format="tif" orientation="portrait" inline="no"/>相乘，抑制高频成分，使其满足向下延拓的条件，即id="icf0003" file="S2008101105005C00021.gif" wi="75" he="5" top= "32" left = "37" img-content="drawing" img-format="tif" orientation="portrait" inline="no"/>其中h为下延深度，p为延拓回返次数；h，p为正有理数；(7)在波数域，在抑制高频成分的基础上，与差分法向下延拓频域响应函数
2、 根据权利要求l所述的重磁资料向下延拓的方法，其特征是所述 的n阶导数或积分，n为有理数，一3《n《3。
3、 根据权利要求1所述的重磁资料向下延拓的方法，其特征是所述 的延拓回返次数P， P为正有理数，0《p《30。
4、根据权利要求1所述的重磁资料向下延拓的方法，其特征是所述向下延拓深度h， h为正有理数，0《h《500Km。
全文摘要
一种地球物理勘探中的重磁数据处理方法，其包括步骤在波数域将重磁资料频谱分别与任意方向的n阶导数或积分频率响应函数、延拓回返响应函数及差分法下延算子频率响应函数相乘，得到下延后的频谱；再经过快速傅氏反变换，获得向下延拓后的三维重磁数据体。利用向下延拓之后的三维重磁数据体结合地震、地质、钻井资料进行对比、分析并提交下一步勘探的有利区带及目标。该方法不仅有更大的向下延拓深度并且有较好的抗噪能力，从而获得稳定的下半空间三维重磁资料分布图像。
文档编号G01V7/00GK101285896SQ200810110500
公开日2008年10月15日 申请日期2008年6月13日 优先权日2008年6月13日
发明者辉 杨 申请人:辉 杨