第八讲 重力异常反演.docx
《第八讲 重力异常反演.docx》由会员分享，可在线阅读，更多相关《第八讲 重力异常反演.docx（53页珍藏版）》请在冰点文库上搜索。
第八讲重力异常反演
应用重力学第八讲
重力异常反演
④正问题（ForwardCalculation）
④已知地质体的形状、产状和剩余密度等，通过理论计算来求得异常的分布和规律。
md
正演
④反问题（Inversion）
④已知异常的分布特征和变化规律，求场源的赋存状态（如产状、形状和剩余密度等）。
?
④解正问题是解反问题的基础，解反问题是目的。
④仅从地质角度，解重力反演问题的目标
9矿体类问题:
寻找、研究或推断金属或非金属矿体；
9构造类问题:
研究地质构造，包括控矿构造，如含石油、天然气、煤的构造以及区域性的深部构造等。
④从地球物理角度，解重力反演问题的目标
9矿体类问题:
确定地质体的几何和物性参数；
9构造类问题:
确定物性分界面的深度及起伏；
9密度分布问题:
确定密度的分布。
一、计算地质模型体的几何及物性参数
（一）直接法
④直接利用由反演目标引起的局部异常，通过某种
积分运算和函数关系，求得与异常分布有关地质体的某些参量。
④三度体剩余质量的求法
⎡mn⎤
⎢m∑∑ijg.u.
mg.u.⎥
M=2.386
{⊗x⊗y}
{⊗g}+2ð{R2}{⊗g（R）}
⎣i=1j=1⎦
④三度体重心水平坐标的求法
2.386mn
{}⊗y[
x}{⊗g}]
x0≈⊗xm∑∑imijg.u.
{M}t
2.386
{⊗y}
i=1j=1
mn
[{y}{⊗g}]
y0≈⊗xm∑∑imijg.u.
{M}t
i=1j=1
④二度体横截面积的求法
n
⊗x∑⊗gi
⎣=i=1
4Gtg−1x
D
④二度体横截面重心水平坐标的求法
x0=
n
∑（x⋅⊗g）i
i=1
n
∑⊗gi
i=1
（二）特征点法
④根据异常曲线上的一些点或特征点（如极大值点、零值点、拐点）的异常值及相应的坐标求取场源体的几何或物性参数;
④仅适用于剩余密度为常数的几何形体。
异常曲线形态分类
④第一类是单峰异常，零值点在无穷远处如球体的Δg曲线、台阶的Vxz曲线等；
④第二类是具有极大值、极小值和一个零值点如球体的Vxz曲线、台阶的Vzz、Vzzz曲线；
④第三类是具有一个极大值、两个极小值和两个零值点如球体、水平圆柱体的Vzz和Vzzz曲线；
④第四类是台阶的Δg曲线，一边高一边低的形态
应用条件
④对异常作平滑处理，尽量准确确定原点的位置；
④对异常曲线作分离处理，获得单纯由研究对象引起的异常；
④对剩余（局部）异常进行分类，判明该异常的场源体接近于何种可能的几何形体，然后选用相应的反演公式。
球体Δg的反演（第一类曲线）
单峰异常，零值点在无穷远处
GMDGMD
⊗g1=
223/2
⊗g2=223/2
（x1+D）（x2+D）
⊗g（x2+D2）3/2
n=1=2
21
⊗g（x2+D2）3/2
1/2
⎛x2−x2n2/3⎞
D=⎜21⎟
n2/3−1
⎝⎠
D=1
2
x1/n−
n2/3
x1′/n
−1
当n=2时
D=1.305x1/2
=0.6524（x1/2
−x1′/2）
当n=3时
当n=4时
D=0.9622x1/3
D=0.8111x1/4
=0.4811（x1/3
=0.4056（x1/4
−x1′/3）
−x1′/4）
max
M=14.99D2⊗g
1/3
⎛⎞
1/3
⎡⎤
R3M
0.62
{M}t
=⎢⎥=⎜⎟
m⎣ð⎛′−⎛
⎦⎜⎛}⎟
4（0）
h=D−R
⎧⎛′⎫
⎝{g/cm3⎠
M实=⎨⎬
⎩⎛⎭g/cm3
{M}t
mm
{D}
=1.305{x1/2}
{M}=14.99{D}2{⊗g}
tmmax
1/3
g.u.
m
{R}=
0.62
⎛{M}⎞
⎜t⎟
⎜{⎛}⎟
⎝gc/m3⎠
D
R
⎛=1g/cm3
mm
{D}=1.305{x1/2}
=1.305⋅153=199.67m
{M}=14.99{D}2{⊗g}
200m
tmmax
g.u.
=14.99⋅199.67⋅199.67⋅6.98=4.1714⋅106t
1/3
⎛⎞
m
{R}=
0.62
{M}t
⎜⎟
⎜{⎛}⎟
⎝gc/m3⎠
=0.62⋅（4.1714⋅106/1）1/3=99.81m
100m
球体Vxz的反演（第二种类型的曲线）
具有极大值、极小值和一个零值点
11
xmax
=−D,xmin
=D
22
（V）=−（V）=
48GM
xzmax
xzmin
255D3
D=2xmin
=−2xmax
=xmin
−xmax
xz
M=0.01746D3（V）
max
球体Vzz的反演（第三种类型曲线）
具有一个极大值、两个极小值和两个零值点
（V）=
2GM
zzmaxD3
x0=2{D}m,x0=−2{D}m
zz
M=0.00749D3（V）
max
0zzmax
=0.00265x3（V）
球体Vzzz的反演（第三种类型曲线）
（V）=
6GM
zzz
maxD4
D=1.225x0
M=2.499D4（V）
⋅10−3=5.622x4（V）
zzz
max0
zzz
max
台阶曲线的反演（第四种类型曲线）
台阶的Δg曲线，一边高一边低的形态
⊗gmax
=⊗g（x）+⊗g（−x）=2ðG⎧
2ðG⎧（2−1
x1/n）
=G⎧ð+tg
nD
D=x1/n
tg（2−n）ð
2n
1
⎧=[⊗g（x）+⊗g（−x）]
2ðG
hx1/n
1
[⊗（x）+⊗g（−x）]
4ðG⎛
=−
（2−n）ð
tg
2n
Hx1/n
1
[⊗g（x）+⊗g（−x）]
4ðG⎛
=+
（2−n）ð
tg
2n
（三）选择法
④根据异常分布和变化特征，结合地质和其他地球物理和物性等资料，给出初始地质体模型;
④进行正演计算，将理论异常与实测异常对比;
m
∑kkkk12n
⎫=[⊗g
k=1
-f（x,y,z,b,b,b）]2
=min
④若两者偏差较大，对模型进行修改，重算其理论异常计算，再次进行对比……;
④如此反复进行，直至两种异常的偏差达到事前要求的误差范围为止，则这最后的理论模型就可作为所求的解答了。
修改模型参数
否
初始模型
计算理论重力异常
与实测异常对比
是否满足精度
是
反演结果
选择法特点
④异常可以是整条剖面或整个测区的数据，受个别点误差的影响较小，抗干扰的能力较强。
④所求的地质体可以是一个或几个复杂的不规则几何形体,密度分界面,或者密度的分布。
④需要重复而复杂的正演计算,可编制相应程序由计算机来自动完成。
④解释复杂重力资料时，能够考虑研究区已知的地质构造资料，在反演过程中利用这些资料，控制或约束计算结果，使得到的地质模型更接近实际的地质体。
（四）人机交互式反演方法
（姚长利--重磁异常正反演解释系统）
二、计算密度分界面的深度
④密度分界面与区域构造和储油构造有密切的关
系，因此计算密度分界面的起伏或深度的变化在区域构造研究和石油勘探中具有重要的意义。
（一）线性回归法
④如果界面起伏平缓，可以认为重力变化与界面的起伏近似呈线性关系。
h=a+b⋅⊗g
n
∑ii
⎫（a,b）=（h−h）2
i=1
=min
∂⎫∂⎫
令==0,=0,则
∂a∂b
2hn
a=∑=⊗gi∑hi−∑⊗gi∑⊗giib==∑⊗gihi−∑hi∑⊗gi
∑∑∑
n⊗g2−（⊗g）2
n⊗g2−（⊗g）2
iiii
（二）压缩质面法（刘云龙,1977）
条件：
界面起伏较小，埋藏深度较大。
基本原理
④将界面从最小深度h和最大深度H处向中间挤压，使之在界面平均深度D=（h+H）/2上压缩成一个面密度不均匀分布的物质面，
④将该物质面剖分成局部面密度均匀分布的水平物质带（二
维）或物质片，面密度
⎧j=⎛⊗hj
④计算物质带或物质片的正演，或迭代反演物质带或物质片的面密度，进而求出界面深度。
g（i）
2Gtg−12（
j−1）−i2（j−2）−1⎤
−1
⎧
⎡
n
⊗=∑⎢
−tg⎥j
j=1⎣
22⎦
（三）迭代法（Cordell,1968）
Δg
y
O
tσ
Q长方体
z
x
P测点
④由无限平板重力异常公式给出t初值
M
t1,q
=gobs（q）
2ðG⎛
④计算模型初值的重力异常
gcalc,n,p
=∑
q=1
g
Gf（P,Q,tn,q,⎛,D）
④计算t的下一个修改值
tn+1,q=
tn,q
（obs,q）
g
M
∑obs,q
calc,q
calc,n,q
④目标函数
⎫=（g−g）2=min
q=1
④判断计算结果满足要求否，不合要求则转到第2步继续计算；否则停止计算。
t1,q
=gobs（q）
2ðG⎛
gcalc,n,p
M
=∑
q=1
Gf（P,Q,tn,q,⎛,D）
n=n+1
tn+1,q=
tn,q
gobs,q
（）
g
calc,n,q
M
∑obs,q
calc,q
⎫=（g−g）2=min?
q=1
否是
反演结果
理论界面深度正演的重力异常
md=Gm
理论界面深度反演界面深度（迭代50次）
mm’=G-1d
（四）频率域反演法（Parker-Oldenburg,1973）
④假定在x-z直角坐标系中，重力异常用⊗g（x）表示，场源层
的上部边界为z=0，下部边界为z=h（x），这个边界显示界
k
面的起伏。
④根据帕克的二维傅里叶变换公式得到重力异常傅立叶变换
F[⊗g（x）]=
−2ð
G⎛e−k
∞
∞
z0∑
n=1
n−1
k
F[hn（x）]
n!
④从上式的无限和式中提出n=1的项，并重新排列，得到
F[h（
z0
x）]=
−=F[⊗=g（x）]e−∑
n−1
k
F[hn
（x）]
2ðG⎛
n=2n!
④假定已知地层与下部介质之间的密度差⎛及参考面深度z0已
知或给定，就可以应用上式进行下列迭代计算：
（1）给界面起伏h（x）的初值，例如h（x）=0，
（2）将h（x）的初值代入上式的右端项，计算右端项的傅立叶变换，
（3）右端项的傅立叶反变换即改进的界面起伏h（x）
（4）判断：
计算结果是否满足某个收敛标准，或是否达到给定的最大迭代次数，如果是，即停止计算；否则转到第2步，以本次迭代结果作为初值，继续迭代计算。
④级数的收敛与界面的起伏大小有关。
记H为界面相对于参考面的最大起伏值，当H/z0<1时，级数均匀收敛并与波数k无关；当H/z0=1时，即界面与观测面相交时，级数不收敛，而且，界面越接近观测面收敛速度越慢。
④上式右端第一项包含了一个指数项，这就会放大数据中的高频干扰，故在右端项乘以一个低通滤波器来加以限制。
反演
正演
三、计算地层密度的分布
④根据重力异常反演地下密度的不均匀分布。
Δg
X
X
σ
Z
④根据密度分布图，就可以判断是否存在高（低）密
度体及其大小、位置，或发现密度的分界面，再结合地质资料就可以作出地质上的结论。
④这个方法思路简单、清楚、直观，但是实现起来相当困难。
④这种方法目前尚处于理论研究阶段或限于推断深部地质构造问题，很少用于解决实际重力勘探问
题。
（一）二维密度成像
⊗g=−2G⎛[xln（x2+z2）+2−1x]
inn
⎩k+1⎛l+1=a⎛
ztg
z
⎩k⎛l
N
inn
gi=∑ain⎛n,i=1,2,...,m
n=1
⎪
⎧a11⎛1
⎪a21⎛1
⎨
+a12⎛1
+a22⎛1
++a1N⎛N=⊗g1
++a2N⎛N=⊗g2
⎪
⎪⎩am1⎛1+am2⎛1++amN⎛N=⊗gm
AP=G
（二）三维密度成像
Δg
y
Ox
σ
z
〉=（⎩2+ç2+⎛2）1/2
⊗g=−G⎛[⎩ln（ç+〉）+çln（⎩+〉）−⎛tg−1
⎩=ç]x2
y2z2
=
⊗g−G[
⎛
⎩ln（ç+〉）+çln（⎩+〉）−⎛tg−1
⎛〉
⎩=ç]
⎛〉
x1y1z1
x2y2z2
x1y1z1
=F
⊗g=F⎛
④下半空间剖分为M个长方体，观测面上有N个观测点
⊗gij=Fij⎛j
M
⊗gi=
∑
j=1
Fij⎛j,
i=1,2,...,N
⎡⊗g1
⎤⎡F1F1
12
F1M
⎤⎡⎛1⎤
⎢⎥⎢⎥⎢⎥
⎢
⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎥
⎢
⎥
⎢⊗g2⎥=
⎢F2F122
⎥⎢⎛2⎥
⎢⎣⊗gN
⎥⎦⎢⎣FN1
FNM⎥⎦⎢⎣⎛M⎥⎦
G=FP
（姚长利，2007）
（刘天佑，2007）
（刘天佑，2007）
重计算地质模型体的几何及物性参数
力直接法、特征点法、选择法、
异人机交互反演法
常
反计算密度分界面的深度
演线性回归法、压缩质面法、迭代法、
方频率域反演法
法计算地层密度的分布
二维、三维密度成像
四、反问题解的多解性
④处在不同深度的、具有不同起伏的界面，还包括在这些界面之间的许多界面，都能够引起在测量精度范围内的相同的重力异常。
（据Skeels，1947）
内特尔顿（1987）指出：
④与一定的异常宽度对应的场源的可能的最大深度，就是引起同样宽度异常的点源（球体）的深度。
④在这个最大深度和地面之间，存在一个可能源的锥形区。
④同一个异常可以由埋藏很浅的薄透镜体，或不太狭窄的较厚的物体，或球体引起，而且在球体上方较浅的深度上，却有无限多个可能的场源引起这个相同的异常。
④引起同一异常的不同的场源，有一个唯一的共同的特性，就是它们的剩余质量必须是相同的。
引起多解性的原因
④场的等效性；
④观测数据离散、有限;
④实测的异常包含一定误差的；
④数据整理带来的误差。
④反问题解的非唯一性是客观存在的。
限制多解性的方法
④除了应用重力资料外，还引用工区内的地质、钻井、物性和其他地球物理资料等，尽可能地增加已知条件和约束条件，则反问题的解答数目就会大大减少，甚至可以得到单一的解。
如球体反问题的求解，如果已知剩余密度值，则它的半径和顶部埋深也就唯一地被确定了。
④提高仪器的测量精度。
④改进各项校正的计算方法，使得能够更加精确地得到重力异常的分布。
（刘天佑，2007）
（刘天佑，2007）
五、重力归一化总梯度法
④原苏联地球物理学家别列兹金于1967年提出了重力归一化总梯度法，它是一种利用重力场中的特征点
（如极大点或极小点）或者解析函数的奇点，探测重力异常的场源体并估计其位置的方法。
④在物体的内部，引力位V满足泊松方程
∇2V=−4ðG⎛
④在物体的外部，引力位V满足拉普拉斯方程
∇2V
∂2V∂2V∂2V
=++=0
∂x2∂y2∂z2
④剩余质量的引力位及导数在场源体以外空间都是解
析函数，而在场源处则失去解析性。
④在解析函数中，失去解析性的点叫函数的奇点。
④确定场源的问题，就是通过对异常解析延拓来确定函数的奇点问题。
以二度水平圆柱为例
⊗g（x,z）=
2G⎣z
z2+x2
④在x=0处，把Δg（x,y）延拓到圆柱中心z=0处，在
圆柱中心处的Δg（0,0）将变为无穷大。
④圆柱中心是极点类型的奇点。
寻找重力场的奇点有一些困难：
④实测异常是由解析成分（有用信号）及随机干扰组成，要用实测异常求奇点首先必须把干扰滤掉；
④实测异常为离散形式，而且其观测剖面长度是有限的，在确定奇点时，需施加某些限制；
④不知道实测异常的解析式，无法求其奇点。
④二度水平圆柱重力异常
⊗gx（,z）=2G⎣
D
D2+x2
④向下延拓到任一点P（x,z）时
⊗g（x,z）=
2G⎣
D−z
④若沿Z轴（x=0）
⊗g（z）=
2G⎣
（D−
1
z）2
+x2
D−z
④在z=0的附近将上式展开成台劳级数为
g（z）
∞1n!
0）n
⊗=∑
2G⎣
n+1
z=0
（z−
n=0
n!
（D−z）
∞∞n
21n2G⎣
⎛z⎞
=∑=∑⎜⎟
G⎣n+1z
n=0DDn=0⎝D⎠
④当向下延拓趋于场源（z→D）时,级数发散，即D为Δg（z）的一个奇点。
④当级数不取无限项时，情况就不同了，这时
2G⎣
N⎛z⎞n
⊗g（z）
=∑⎜⎟
Dn=0⎝D⎠
④通过柱体中心（z≥D），因式是一有限项级数，只
有当z→∞时才发散，否则只要N为有限值，级数都是收敛的。
④所以向下延拓直接求与场源质量有关的奇点是很困难的。
④作如下改进：
2G⎣
n
N⎛z⎞
2G⎣
n
∞⎛z⎞
⊗g（z）
=∑⎜⎟
+∑
⎜⎟
Dn=0⎝D⎠
Dn=N+1⎝D⎠
④右边的第二项是级数的余项，也即舍去误差。
④令ζ=z/D
⊗g（⎛）=
⊗gN（⎛
）+rN（⎛）
④若N取适当值，则ζ<1时，有rN（ζ）<ΔgN（ζ），反之，在ζ>1时将有rN（ζ）>ΔgN（ζ），用它可以确定与场源有关的奇点。
④以rN（ζ）对ΔgN（ζ）进行归一化，则
g
升级会员 