CN100349014C - 重力勘探数据处理变密度地形校正方法 - Google Patents

Abstract

一种重力勘探数据处理变密度地形校正方法，用重力勘探的常规方法采集测量目标地区地球表面的重力场强度，以表层最大密度进行常密度地形校正，然后再采用表层变密度补充地形校正的方法消除过度校正所造成的影响，所采用的模型是同时考虑沿地表和垂直于地表两个方向变化的变密度模型。本发明在地形起伏较大的地区，在经过常密度地形校正的布格重力异常后仍与地形呈负相关的虚假异常明显的情况下，经表层变密度补充地形校正，可以使重力异常的走向与区域地质走向一致，是一种稳妥、易行而有效的重力勘探数据处理变密度地形校正方法。

Description

重力勘探数据处理变密度地形校正方法

技术领域

本发明涉及石油地质勘探开发中的地球物理重力勘探数据处理技术，具体是一种重力勘探数据处理变密度地形校正方法。

技术背景

地球的引力场称为重力场，地球表面的重力场强度因地壳密度的变化而有微量差异，称为重力异常。测量地球表面的重力场强度，计算重力异常，并根据重力异常反算地壳密度分布，进而探测地下矿产和地质构造的地球物理方法称为重力勘探方法，简称重力勘探。

重力勘探包括重力(场强度)数据的观测、处理和推断解释三个环节。重力勘探数据处理是指将观测重力数据通过地形校正、中间层校正和正常场校正(含高度校正)获得重力异常，并对重力异常进行分解、转换的整个过程。其中的地形校正是一项比较基础而又复杂的校正。

所谓地形校正就是：过重力观测点作水平面，计算水平面以上实际存在的地形质量对观测点的重力效应，并从测点观测重力值中减去，然后以与地形质量等密度的物质充填水平面到地形面之间的空间，计算这些充填物对观测点的重力效应，并加到测点观测重力值中。地形校正和中间层校正合二为一的地球物理意义是消除大地水准面与地形面之间的物质对测点观测重力值的影响。

常规地形校正和中间层校正采用的是均匀密度模型。这一模型对地形体质量分布的模拟整体上是可以的，但对表层的模拟程度较差，因为表层密度受出露地层岩性和风化程度的影响，变化较大。远处表层密度变化对地形校正的影响不大，但测点数公里内地表密度变化对地形校正有明显的影响，尤其是山区。

地表密度不均匀的情况下对地形校正将产生异常影响，如果采用常密度进行校正就会出现两种情况：真实密度大于校正密度的地方，校正不足；真实密度小于校正密度的地方，校正过度。校正不足的地方就会产生与地形正相关的假异常(山形异常)，校正过度的地方就会产生与地形负相关的假异常(镜像异常)。假异常达到一定的量级，在布格重力异常图上就会表现为等值线的同向扭曲(图6)，重力垂直二次导数处理后就会形成圈闭异常(图8)，与地下局部构造在重力垂直二次导数图上形成的圈闭异常混在一起，真假难辨，导致错误地地质解释。因此，在山区进行高精度重力勘探，必须进行变密度地形校正。

发明内容

本发明的目的在于提供一种表层变密度补充地形校正与常密度地形校正相结合，稳妥、易行而有效的重力勘探数据处理变密度地形校正方法。

本发明在常密度地形校正的基础上进行表层补充地形校正的变密度地形校正，简称为“常规校正+表层补充校正”的变密度校正方法，是稳妥而有效的方案。所谓“表层补充校正”就是计算常密度地形校正后表层残余质量对测点的重力效应，并加以消除。

本发明的“常规校正+表层补充校正”的变密度校正方法，是在常密度地形校正的基础上进行表层三维变密度补充地形校正。

本发明通过以下技术方案实现：

采用重力勘探的常规方法采集测量目标地区地球表面的重力场强度，然后，

(1)实测或根据目标地区现有的地面密度数据，确定地表最大密度σ₀和下界残余密度δ，并以σ₀进行常密度地形校正和布格校正，得到布格重力异常图

(2)将地面密度数据σ(ξ，η)网格化，获得地面密度网格化数据σ(m，n)，依公式Δσ(m，n)＝σ(m，n)-σ₀①，计算残余密度Δσ(m，n)，并绘制等值线图；

(3)比较布格重力异常图

、地形图h₂(m，n)、残余密度分布图Δσ(m，n)，确定是否存在山形异常；

(4)根据实际条件确定残余密度随深度的衰减率的初始值k₀(m，n)；

(5)依据以下公式计算地形残余密度体的下界k₁₀(m，n)，

$h_{10}(m,n)=h_2(m,n)+\frac{1}{k_0(m,n)}\·\ln \frac{\mathrm{\δ}}{\mathrm{\Δ\σ}(m,n)}$    ②

(6)计算上顶为h₂(m，n)、下底为h₁₀(m，n)、密度模式为： $\mathrm{\Δ\σ}(m,n,h)=\mathrm{\Δ\σ}(m,n)e^{-k(m,n)\·[h_2(m,n)-h]}$ 的方柱体的重力效应，并从布格重力异常

中减去，得到

以上计算分别采用以下公式：

(a)z＞h₂＞h₁，即柱体质量全部在观测点之上：

$\mathrm{\Δ}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\πG\Δ\σ}\frac{\sin 2\mathrm{\πua}}{\mathrm{\πu}}\·\frac{\sin 2\mathrm{\πvb}}{\mathrm{\πv}}\·e^{-i2\mathrm{\πu}{\mathrm{\ξ}}_0}\·e^{-i2\mathrm{\πv}{\mathrm{\η}}_0}$

$\·\frac{1}{k+2\mathrm{\π}\sqrt{u^2+v^2}}[e^{-2\mathrm{\π}\sqrt{u^2+v^2}(z-h_2)}-e^{-2\mathrm{\π}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]$    ③

(b)z＜h₁＜h₂，即柱体质量全部在观测点之上：

$\mathrm{\Δ}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=-2\mathrm{\πG\Δ\σ}\frac{\sin 2\mathrm{\πua}}{\mathrm{\πu}}\·\frac{\sin 2\mathrm{\πvb}}{\mathrm{\πv}}\·e^{-i2\mathrm{\πu}{\mathrm{\ξ}}_0}\·e^{-i2\mathrm{\πv}{\mathrm{\η}}_0}$

$\·\frac{1}{k-2\mathrm{\π}\sqrt{u^2+v^2}}[e^{-2\mathrm{\π}\sqrt{u^2+v^2}(h_2-z)}-e^{-2\mathrm{\π}\sqrt{u^2+v^2}(h_1-z)-k(h_2-h_1)}]$    ④

(c)h₁≤z≤h₂，即柱体顶面在观测点之上，底面在观测点之下：

$\mathrm{\Δ}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\πG\Δ\σ}\frac{\sin 2\mathrm{\πua}}{\mathrm{\πu}}\·\frac{\sin 2\mathrm{\πvb}}{\mathrm{\πv}}\·e^{-i2\mathrm{\πu}{\mathrm{\ξ}}_0}\·e^{-i2\mathrm{\πv}{\mathrm{\η}}_0}$

$\·\left\{\frac{1}{k+2\mathrm{\π}\sqrt{u^2+v^2}}\right[e^{-k(h_2-z)}-e^{-2\mathrm{\π}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]-\frac{1}{k-2\mathrm{\π}\sqrt{u^2+v^2}}[e^{-2\mathrm{\π}\sqrt{u^2+v^2}(h_2-z)}-e^{-k(h_2-z)}]\}$    ⑤

(7)对 、h₂(m，n)作相关分析，根据相关性的正负调整衰减率k₀(m，n)，得到k₁(m，n)，重复以上(5)至(7)，直至布格异常与实际地形的相关性最小；

公式③④⑤是负指数型变密度模式的方柱体模型的波数域正演公式，在计算方柱体重力效应时将它们离散化后，再进行反傅里叶变换数值计算；

公式③④⑤采用傅里叶变换数值计算的偏移抽样技术计算，为：

$\mathrm{\Δg}(x,y,z)$

$=d_1{d}_2\underset{m_1=-N1/2}{\overset{N1/2}{\mathrm{\Σ}}}\underset{m_2=-N2/2}{\overset{N2/2}{\mathrm{\Σ}}}\mathrm{\Δ}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\η}}^{*})d_1,(-\frac{N_1}{2}+m_2+{\mathrm{\η}}^{*})d_2,z]$

以上公式中各个符号的意义如下：

x，y，z为计算点的三维空间坐标；

ξ₀，η₀为方柱体顶、底面中心的平面坐标，h₂，h₁分别为方柱体顶、底面的高程，a，b，分别为顶、底面的半长和半宽，Δσ为方柱体顶面密度，k为方柱体内密度随深度的指数衰减率；

为方柱体在x，y，z产生的重力效应，

为其波谱u，v分别为x，y方向上的波数；

$\mathrm{\Δ}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\η}}^{*})d_1,(-\frac{N_2}{2}+m_2+{\mathrm{\η}}^{*})d_2,z]$ 为波谱函数的

的抽样序列，其中；d₁，d₂分别为u，v方向上的抽样间隔，m₁，m₂分别为u，v方向上的抽样点序号，N₁，N₂分别为u，v方向上的抽样点数，η^*为u，v方向上的偏移抽样参数；G为万有引力常数。

本发明下界残余密度δ为-0.01g/cm³。

本发明的采用表层变密度补充地形校正与常密度地形校正相结合，在地形起伏较大的地区，在经过常密度地形校正的布格重力异常后仍与地形呈负相关的虚假异常明显的情况下，经表层变密度补充地形校正，可以使重力异常的走向与区域地质走向一致。本发明是一种稳妥、易行而有效的重力勘探数据处理变密度地形校正方法。

从图5、图6、图7这三张实际测量布格重力异常图上可以看出：在地形起伏较大的地区，它们的重力异常形态的差异是比较大的。未做地形校正的布格重力异常图(图5)，重力等值线零乱得很，与地形相关性较强；经过常密度地形校正的布格重力异常图，等值线圆滑了许多，但与地形呈负相关的虚假异常仍非常明显：对应于山脊是向高值方向同向扭曲的一组等值线，对应于山谷是向低值方向同向扭曲的一组等值线。在做过三维变密度补充地形校正的布格重力异常图上，这类与地形呈负相关的虚假异常不见了，而且等值线更加圆滑，重力异常的走向与区域地质走向是一致的。这就是三维变密度补充地形校正的效果，是在布格重力异常图展现出的明显效果。

三维变密度补充地形校正的效果在重力垂直二阶导数异常图上表现得更醒目。从图8、图9两张重力垂直二次导数异常图上可以看出：常密度地形校正后的重力垂直二阶导数异常图(图8)上幅度相当大的三排与地形相关的北东向虚假异常，在经过三维变密度补充地形校正后的重力垂直二阶导数异常图(图9)上不见了踪影。经过三维变密度补充地形校正后的重力垂直二阶导数异常图上显现出的重力垂直二次导数异常具有较明显的北西走向，与整体的构造走向相一致。这些局部重力异常是地下地质情况的反映，与其它物探资料的推断解释一致，也与已知的地质资料一致。

附图说明

图1是地表残余密度分布图；

图2是测区高程分布图；

图3是变密度地形校正底界高程图；

图4(a)是东西向变密度地形体密度剖面；

图4(b)是南北向变密度地形体密度剖面；

图5是未做地形校正的布格重力异常图；

图6是经过常密度地形校正的布格重力异常图；

图7是经过三维变密度补充地形校正的布格重力异常图；

图8是经过常密度地形校正后的重力垂直二阶导数异常图；

图9是经过三维变密度补充地形校正后的重力垂直二阶导数异常图。

具体实施方式

以下结合附图说明本发明具体实施方式。

本发明的内容是重力勘探数据处理中的一项变密度地形校正技术，简称“常规校正+表层补充校正”的变密度校正方法，因为其基本思想是在常密度地形校正的基础上进行表层三维变密度补充地形校正，即以表层最大密度进行常密度地形校正，然后再采用表层变密度补充地形校正的方法消除过度校正所造成的影响。

表层变密度补充地形校正所采用的模型是同时考虑沿地表和垂直于地表两个方向变化的变密度模型。其具体形式如下：

$\mathrm{\Δ\σ}(\mathrm{\ξ},\mathrm{\η},h)=\mathrm{\Δ\σ}(\mathrm{\ξ},\mathrm{\η})e^{-k(\mathrm{\ξ},\mathrm{\η})\·[h_2(\mathrm{\ξ},\mathrm{\η})-h]}$

式中：Δσ(ξ，η)＝σ(ξ，η)-σ₀为地形体表面残余密度，简称地表残余密度；

σ(ξ，η)为地表密度；

k(ξ，η)为地形体残余密度的纵向衰减率；

h₂(ξ，η)为地形点高程；

h为地形体内p点的高程。

表层变密度模型的底界是由地表密度Δσ(ξ，η)和地形h₂(ξ，η)两个因素共同决定的。确定底界的具体公式是：

$h_1(\mathrm{\ξ},\mathrm{\η})=h_2(\mathrm{\ξ},\mathrm{\η})+\frac{1}{k(\mathrm{\ξ},\mathrm{\η})}\·\ln \frac{\mathrm{\δ}}{\mathrm{\Δ\σ}(\mathrm{\ξ},\mathrm{\η})}$

表层变密度模型的重力效应采用波数域正演技术计算。

“常规校正+表层补充校正”的变密度校正的具体做法如下：

(1)搜集地面密度资料(吸纳整理前人关于该区的地面密度数据，必要时进行实测)，确定地表最大密度σ₀和下界残余密度δ(一般-0.01g/cm³即可)，并以σ₀行常密度地形校正和布格校正，得到布格重力异常

(如图6)。

(2)将地面密度数据σ(ξ，η)网格化，获得地面密度网格化数据σ(m，n)，依据公式(1)

        Δσ(m，n)＝σ(m，n)-σ₀    ①

计算残余密度Δσ(m，n)，并绘制等值线图(如图1)。

(3)观察分析布格重力异常图

与地形h₂(m，n)(图2)和残余密度分布Δσ(m，n)(图1)的关系，确定是否存在山形异常(负相关)。

(4)确定(根据经验选择)残余密度随深度的衰减率的初始值k₀(m，n)。

(5)依据(2)式计算地形残余密度体的下界h₁₀(m，n)(如图3)：

$h_{10}(m,n)=h_2(m,n)+\frac{1}{k_0(m,n)}\·\ln \frac{\mathrm{\δ}}{\mathrm{\Δ\σ}(m,n)}$    ②

(6)根据以下(a)、(b)、(c)三种不同情况，分别采用公式(3)或(4)或(5)计算上顶为h₂(m，n)、下底为h₁₀(m，n)、密度模式(图4)为：

$\mathrm{\Δ\σ}(m,n,h)=\mathrm{\Δ\σ}(m,n)e^{-k(m,n)\·[h_2(m,n)-h]}$

的一组方柱体的重力效应，并从布格重力异常

中减去，得到

。(a)、

(b)、(c)三种不同情况对应的公式(3)、(4)、(5)的形式如下：

(a)z＞h₂＞h₁，即柱体质量全部在观测点之上：

$\mathrm{\Δ}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\πG\Δ\σ}\frac{\sin 2\mathrm{\πua}}{\mathrm{\πu}}\·\frac{\sin 2\mathrm{\πvb}}{\mathrm{\πv}}\·e^{-i2\mathrm{\πu}{\mathrm{\ξ}}_0}\·e^{-i2\mathrm{\πv}{\mathrm{\η}}_0}$

$\·\frac{1}{k+2\mathrm{\π}\sqrt{u^2+v^2}}[e^{-2\mathrm{\π}\sqrt{u^2+v^2}(z-h_2)}-e^{-2\mathrm{\π}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]$    ③

(b)z＜h₁<h₂，即柱体质量全部在观测点之上：

$\mathrm{\&Delta;}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=-2\mathrm{\&pi;G\&Delta;\&sigma;}\frac{\sin 2\mathrm{\&pi;ua}}{\mathrm{\&pi;u}}\&CenterDot;\frac{\sin 2\mathrm{\&pi;vb}}{\mathrm{\&pi;v}}\&CenterDot;e^{-i2\mathrm{\&pi;u}{\mathrm{\&xi;}}_0}\&CenterDot;e^{-i2\mathrm{\&pi;v}{\mathrm{\&eta;}}_0}$

$\&CenterDot;\frac{1}{k-2\mathrm{\&pi;}\sqrt{u^2+v^2}}[e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_2-z)}-e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_1-z)-k(h_2-h_1)}]$    ④

(c)h₁≤z≤h₂，即柱体顶面在观测点之上，底面在观测点之下：

$\mathrm{\&Delta;}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\&pi;G\&Delta;\&sigma;}\frac{\sin 2\mathrm{\&pi;ua}}{\mathrm{\&pi;u}}\&CenterDot;\frac{\sin 2\mathrm{\&pi;vb}}{\mathrm{\&pi;v}}\&CenterDot;e^{-i2\mathrm{\&pi;u}{\mathrm{\&xi;}}_0}\&CenterDot;e^{-i2\mathrm{\&pi;v}{\mathrm{\&eta;}}_0}$

$\&CenterDot;\left\{\frac{1}{k+2\mathrm{\&pi;}\sqrt{u^2+v^2}}\right[e^{-k(h_2-z)}-e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]-\frac{1}{k-2\mathrm{\&pi;}\sqrt{u^2+v^2}}[e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_2-z)}-e^{-k(h_2-z)}\left]\right\}$    ⑤

式中：

(7)由于变密度地形校正方案是一个迭代过程，要反复修改残余密度模型中的参数，直至达到最佳结果，要对

、h₂(m，n)作相关分析，根据相关性的正负调整衰减率k₀(m，n)，得到k₁(m，n)，重复(5)---(7)，直至布格异常与地形的相关性最小。

需要特别说明的有以下两点：

1)公式③④⑤是一个负指数型变密度模式的方柱体模型的波数域正演公式，利用它们计算方柱体重力效应时，实际的做法是将它们离散化(即抽样)后，再进行反傅里叶变换数值计算。

2)因为公式③④⑤含有分式，所以必须采用傅里叶变换数值计算的偏移抽样技术，才能实现。具体计算公式是：

$\mathrm{\&Delta;g}(x,y,z)$

$=d_1{d}_2\underset{m_1=-N1/2}{\overset{N1/2}{\mathrm{\&Sigma;}}}\underset{m_2=-N2/2}{\overset{N2/2}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\&eta;}}^{*})d_1,(-\frac{N_1}{2}+m_2+{\mathrm{\&eta;}}^{*})d_2,z]$

公式③④⑤⑥中各个符号的意义如下：

x，y，z为计算点的三维空间坐标；

ξ₀，η₀为方柱体顶、底面中心的平面坐标，h₂，h₁分别为方柱体顶、底面的高程，α，b，分别为顶、底面的半长和半宽，Δσ为方柱体顶面密度，k为方柱体内密度随深度(离开顶面的距离)的指数衰减率；

为方柱体在x，y，z产生的重力效应，

为其波谱u，v分别为x，y方向上的波数；

$\mathrm{\&Delta;}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\&eta;}}^{*})d_1,(-\frac{N_2}{2}+m_2+{\mathrm{\&eta;}}^{*})d_2,z]$ 为波谱函数的

的抽样序列，其中；d₁，d₂分别为u，v方向上的抽样间隔，m₁，m₂分别为u，v方向上的抽样点序号，N₁，N₂分别为u，v方向上的抽样点数，η^*为u，v方向上的偏移抽样参数；

G为万有引力常数。

Claims (2)

1、一种重力勘探数据处理变密度地形校正方法，用重力勘探的常规方法采集测量目标地区地球表面的重力场强度，其特征在于即以表层最大密度进行常密度地形校正，然后再采用表层变密度补充地形校正的方法消除过度校正所造成的影响，所采用的模型是同时考虑沿地表和垂直于地表两个方向变化的变密度模型，具体采用以下技术方案实现：

表层变密度补充地形校正采用沿地表和垂直于地表两个方向变化的变密度模型，具体如下：

$\mathrm{\&Delta;\&sigma;}(\mathrm{\&xi;},\mathrm{\&eta;},h)=\mathrm{\&Delta;\&sigma;}(\mathrm{\&xi;},\mathrm{\&eta;})e^{-k(\mathrm{\&xi;},\mathrm{\&eta;})\&CenterDot;[h_2(\mathrm{\&xi;},\mathrm{\&eta;})-h]}$

式中：Δσ(ξ，η)＝σ(ξ，η)-σ₀为地形体表面残余密度，简称地表残余密度；

σ(ξ，η)  为地表密度；

k(ξ，η)   为地形体残余密度的纵向衰减率；

h₂(ξ，η)  为地形点高程；

h           为地形体内p点的高程，

表层变密度模型的底界是由地表密度Δσ(ξ，η)和地形h₂(ξ，η)两个因素共同决定的，具体公式是：

$h_1(\mathrm{\&xi;},\mathrm{\&eta;})=h_2(\mathrm{\&xi;},\mathrm{\&eta;})+\frac{1}{k(\mathrm{\&xi;},\mathrm{\&eta;})}\&CenterDot;\ln \frac{\mathrm{\&delta;}}{\mathrm{\&Delta;\&sigma;}(\mathrm{\&xi;},\mathrm{\&eta;})}$

表层变密度模型的重力效应采用波数域正演计算；

(1)实测或根据目标地区现有的地面密度数据，确定地表最大密度σ₀。和下界残余密度δ，并以σ₀进行常密度地形校正和布格校正，得到布格重力异常图

(2)将地面密度数据σ(ξ，η)网格化，获得地面密度网格化数据σ(m，n)，依公式Δσ(m，n)＝σ(m，n)-σ₀①，计算残余密度Δσ(m，n)，并绘制等值线图；

(3)比较布格重力异常图

、地形图h₂(m，n)、残余密度分布图Δσ(m，n)，确定是否存在山形异常；

(4)根据实际条件确定残余密度随深度的衰减率的初始值k₀(m，n)；

(5)依据以下公式计算地形残余密度体的下界h₁₀(m，n)，

$h_{10}(m,n)=h_2(m,n)+\frac{1}{k_0(m,n)}\&CenterDot;\ln \frac{\mathrm{\&delta;}}{\mathrm{\&Delta;\&sigma;}(m,n)}$ ②

(6)计算上顶为h₂(m，n)、下底为h₁₀(m，n)、密度模式为： $\mathrm{\&Delta;\&sigma;}(m,n,h)=\mathrm{\&Delta;\&sigma;}(m,n)e^{-k(m,n)\&CenterDot;[h_2(m,n)-h]}$ 的方柱体的重力效应，并从布格重力异常 中减去，得到

以上计算分别采用以下公式：

(a)z＞h₂＞h₁，即柱体质量全部在观测点之上：

$\mathrm{\&Delta;}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\&pi;G\&Delta;\&sigma;}\frac{\sin 2\mathrm{\&pi;ua}}{\mathrm{\&pi;u}}\&CenterDot;\frac{\sin 2\mathrm{\&pi;vb}}{\mathrm{\&pi;v}}\&CenterDot;e^{-i2\mathrm{\&pi;u}{\mathrm{\&xi;}}_0}\&CenterDot;e^{-i2\mathrm{\&pi;v}{\mathrm{\&eta;}}_0}$

$\&CenterDot;\frac{1}{k+2\mathrm{\&pi;}\sqrt{u^2+v^2}}[e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(z-h_2)}-e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]$ ③

(b)z＜h₁＜h₂，即柱体质量全部在观测点之上：

$\mathrm{\&Delta;}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=-2\mathrm{\&pi;G\&Delta;\&sigma;}\frac{\sin 2\mathrm{\&pi;ua}}{\mathrm{\&pi;u}}\&CenterDot;\frac{\sin 2\mathrm{\&pi;vb}}{\mathrm{\&pi;v}}\&CenterDot;e^{-i2\mathrm{\&pi;u}{\mathrm{\&xi;}}_0}\&CenterDot;e^{-i2\mathrm{\&pi;v}{\mathrm{\&eta;}}_0}$

$\&CenterDot;\frac{1}{k-2\mathrm{\&pi;}\sqrt{u^2+v^2}}[e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_2-z)}-e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_1-z)-k(h_2-h_1)}]$ ④

(c)h₁≤z≤2，即柱体顶面在观测点之上，底面在观测点之下：

$\mathrm{\&Delta;}\tilde{g}(u,v,z)=-\frac{d}{\mathrm{dz}}\tilde{V}(u,v,z)$

$=2\mathrm{\&pi;G\&Delta;\&sigma;}\frac{\sin 2\mathrm{\&pi;ua}}{\mathrm{\&pi;u}}\&CenterDot;\frac{\sin 2\mathrm{\&pi;vb}}{\mathrm{\&pi;v}}\&CenterDot;e^{-i2\mathrm{\&pi;u}{\mathrm{\&xi;}}_0}\&CenterDot;e^{-i2\mathrm{\&pi;v}{\mathrm{\&eta;}}_0}$

$\&CenterDot;\left\{\frac{1}{k+2\mathrm{\&pi;}\sqrt{u^2+v^2}}\right[e^{-k(h_2-z)}-e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(z-h_1)-k(h_2-h_1)}]-\frac{1}{k-2\mathrm{\&pi;}\sqrt{u^2+v^2}}[e^{-2\mathrm{\&pi;}\sqrt{u^2+v^2}(h_2-z)}-e^{-k(h_2-z)}\left]\right\}$

(7)对

、h₂(m，n)作相关分析，根据相关性的正负调整衰减率k₀(m，n)，得到k₁(m，n)，重复以上(5)至(7)，直至布格异常与实际地形的相关性最小；

公式③④⑤是负指数型变密度模式的方柱体模型的波数域正演公式，在计算方柱体重力效应时将它们离散化后，再进行反傅里叶变换数值计算；

公式③④⑤采用傅里叶变换数值计算的偏移抽样技术计算，为：

$=d_1{d}_2\underset{m_1=-N_1/2}{\overset{N_1/2}{\mathrm{\&Sigma;}}}\underset{m_2=-N_2/2}{\overset{N_2/2}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\&eta;}}^{*})d_1,(-\frac{N_1}{2}+m_2+{\mathrm{\&eta;}}^{*})d_2,z]$

以上公式中各个符号的意义如下：

x，y，z为计算点的三维空间坐标；

ξ₀，η₀为方柱体顶、底面中心的平面坐标，h₂，h₁分别为方柱体顶、底面的高程，a，b，分别为顶、底面的半长和半宽，Δσ为方柱体顶面密度，k为方柱体内密度随深度的指数衰减率；

为方柱体在x，y，z产生的重力效应，

为其波谱u，v分别为x，y方向上的波数；

$\mathrm{\&Delta;}\tilde{g}[(-\frac{N_1}{2}+m_1+{\mathrm{\&eta;}}^{*})d_1,(-\frac{N_2}{2}+m_2+{\mathrm{\&eta;}}^{*})d_2,z]$ 为波谱函数的

的抽样序列，其中；d₁，d₂分别为u，v方向上的抽样间隔，m₁，m₂分别为u，v方向上的抽样点序号，N₁，N₂分别为u，v方向上的抽样点数，η^*为u，v方向上的偏移抽样参数；G为万有引力常数。

2、根据权利要求1所述的一种重力勘探数据处理变密度地形校正方法，其特征在于下界残余密度δ为-0.01g/cm³。

Images