徐文复杂地质体重力异常计算及分析研究.docx
- 文档编号:25033653
- 上传时间:2023-06-04
- 格式:DOCX
- 页数:50
- 大小:814.64KB
徐文复杂地质体重力异常计算及分析研究.docx
《徐文复杂地质体重力异常计算及分析研究.docx》由会员分享,可在线阅读,更多相关《徐文复杂地质体重力异常计算及分析研究.docx(50页珍藏版)》请在冰豆网上搜索。
徐文复杂地质体重力异常计算及分析研究
防灾科技学院
毕业论文
题 目
复杂地质体重力异常计算及分析研究
学生姓名
徐文
学 号
115043331
系别
地震科学系
专业
地球物理学
班 级
1150114
开题时间
2015年3月14日
答辩时间
2015年6月8日
指导教师
武晔
职称
教授
复杂地质体重力异常计算及分析研究
作者:
徐文
指导教师:
武晔
摘要:
本文详细推导了规则地质体重力异常的计算公式,并且对规则地质体和复杂地质体模型进行了正演数值计算,进而得到重力异常曲线,并采用延拓方法对复杂地质体重力异常数据处理及分析研究。
关键词:
MATLAB、重力异常分析解释
引言1
1重力异常正演计算的基本原理2
1.1规则地质体的正演计算3
1.1.1密度均匀的球体3
1.1.2均匀密度的水平圆柱体3
1.1.3均匀密度的倾斜台阶4
1.2复杂地质体的正演计算5
1.2.1无限长水平圆柱体的重力异常叠加5
1.2.2圆柱体与倾斜台阶叠加5
2重力异常的数值模拟6
2.1规则地质体的数值模拟6
2.1.1球体重力异常模拟分析6
2.1.2均匀密度的水平圆柱体异常9
2.2.3均匀密度的倾斜台阶12
2.2复杂地质体的数值模拟14
2.2.1无限长水平圆柱体的重力异常叠加14
2.2.2圆柱体与倾斜台阶异常的叠加24
3重力异常的延拓31
4结论33
感谢词34
附录34
重力勘探的历史34
我国在重力勘探方面的发展36
重力勘探的应用36
参考文献37
引言
重力异常的正演就是根据观测重力异常求取观测异常的场源体,首先我们必须了解不同形状、大小、产状、和场源密度等的场源体或地质体所引起的重力异常的特征、大小、分布等。
重力异常的正演就是求解这个问题。
在目前流行的反演问题方法中,模型体的正演计算就是反演过程的重要组成。
正演计算中,首先我们应该研究一些简单规则的物体引起的重力异常,例如球体、圆柱体、台阶及半平面等引起的重力异常。
研究这些简单地质体的正演的目的,一方面这些简单地质体接近于一些实际情形下的一些地质体;另一方面,复杂地质体的异常可以用多个简单的地质体异常叠加得到,所以为了研究重力异常和异常场源体的对应关系我们必须研究重力异常的正问题。
MATLAB是当今最优秀的应用科技软件之一,它以强大的科学计算与可视化功能、简单方便、开放式可扩展环境,特别是所附带的三十多种面向不同领域的工具箱支持,使它在许多科学领域当中成为了辅助设计、算法研究和应用开发的基本工具和首选平台。
MATLAB具有其他高级语言难以比拟的一些优点,编写简单,编程效率高,易学易懂。
在地球物理正反演计算、时频信号处理等方面,MATLAB都被广泛地使用,已经被人认可为能够有效提高工作效率、改善设计手段的工具软件。
许多学者对重力异常正演研究,看起来太繁琐,没有运用MATLAB使得文章写得比较空洞,我使用MATLAB对重力异常进行正演,就是为了把重力异常正演做得更加直观,让我们能直观的认识到各种地质体对地面的重力异常曲线规律。
1重力异常正演计算的基本原理
地球表面的任何物体都受到重力的影响,即收到离心力和地球万有引力的合力的作用。
地球的重力值随地球表面观测点变化而变化。
重力的变化与地下物质密度分布不均匀密切相关;而地下物质密度与地质构造与地下矿产分布息息相关。
因此,研究地下物质密度分布不均匀引起的重力变化,即重力异常,便可以了解和推断地下的结构、地壳的构造,以及某些矿产资源等。
但是要了解与地下物质变化引起的重力值变化绝不是一个简单的事。
首先,这种变化是非常小的相对于重力。
例如一个局部地质体或者矿床引起的重力值变化不到整个地球的1000万分之一。
所以,要观测到这个微小的变化,首先必须采用,精度高,灵明度高、稳定性好、适合野外条件工作、便于携带的仪器。
其次重力仪器测量的值不全是重力异常值,还有其他外界因素的影响,例如空气,温度及轻微的振动都会对仪器读数产生影响。
再者,根据仪器读数计算的重力值,不完全是由地下地质体引起的,它包含了地形起伏、测点的高程变化、地球并非球体以及自转引起的重力变化,只有消除这些影响,我们才能得到与由地下物质密度不均匀引起的重力变化相接近。
至此,对于重力勘探的工作并没有完成。
一个测点的重力异常是由地下地质体或者所有的密度分布不均匀引起的叠加异常,要得到地下某个地质体,例如一个可能的矿床,潜在的储油气构造等的位置、产状、大小等信息,必须从叠加异常中分离出单纯由这些勘探目标引起的异常。
把勘探目标的异常分离出来,求出或反演引起这个异常的地质体。
在重力勘探中最困难的问题也有重力异常的分离和反演。
重力解释存在多解性,因为重力勘探中对地下物质的了解是通过反演得出的地质体,我们不能清晰的看出地下到底有什么东西,假如在一些有利条件中,也许我们能得到比较可靠的重力解释。
重力勘探中不能脱离重力场的基本依据,即重力的物理场及场源之间的关系,重力勘探就是依据他们之间的关系来推断场源的情况。
1.1规则地质体的正演计算
1.1.1密度均匀的球体
在我们实际勘探中,一些近似等轴状的地质体,如岩株、矿巢、旷囊,隆状构造等,都可以把它们的异常当做球体异常来计算,假如地质体的水平尺寸小于它的埋藏深度,效果更好。
在以下的地质体模型中,为了简便计算,总把地面当做水平面,即是XOY坐标面,Z轴的方向铅垂向下,代表重力的方向。
对于均匀的球体,它的剩余质量全部集中于重心也就是原点,假设埋深为D,球体的半径为R,剩余密度为
,它的剩余质量就是
。
以下为球体重力异常
(三度球体)
(二度球体)
原点处的重力异常值最大,其最大值为
式子里面的
,D的单位为m,M的单位为t(吨)
1.1.2均匀密度的水平圆柱体
在对于一些横截面近似于圆形、沿着水平方向延伸较长的地质体,如两边比较陡峭的长轴背斜及向斜构造、扁豆状矿体等,在一定精度范围内,它们的异常可以近似于水平圆柱体异常来看待
假设圆柱体沿着y方向延伸,半径为R,中轴线埋藏深度为D,剩余密度为
,则剩余线密度
。
若取坐标原点处在中轴线中点在地面的投影处,且让Y轴于中轴线平行,则重力异常在任意位置的表达式为:
(二度水平圆柱体)
(三度水平圆柱体)
1.1.3均匀密度的倾斜台阶
倾斜台阶是一类非常,常见的构造,它近似表示地层的超覆、倾斜的接触带以及倾斜断裂等。
计算它的重力异常时,我们选定下图1-1作为参考图,原点在倾斜面的延长线与地面的交点上。
X轴为其走向。
假设台阶的顶面与底面埋藏深度分别为h、H,剩余密度为
,倾斜角为
,则计算异常的理论公式如下。
1.2复杂地质体的正演计算
为了了解和掌握重力场和场源,即重力异常和异常场源的对应关系,我们必须研究重力正演问题,所谓正演问题就是给定地下地质体的产状,密度等通过理论知识求取地表的异常大小、特征和变化规律等。
研究复杂地质体,就必须研究以上比较简单的地质体类型,因为某些复杂地质体就是,多个简单地质体的叠加得到。
所以以下的地质体是由多个简单地质体的叠加得到的地质体。
1.2.1无限长水平圆柱体的重力异常叠加
无限长水平的叠加,就是通过对多个水平圆柱体的研究,分析研究其重力异常规律,对比与一个水平圆柱体异常的区别。
叠加用的公式基本和单个水平圆柱体公式一样公式如下
(三度体)
(二度体)
1.2.2圆柱体与倾斜台阶叠加
在实际生活中,地下物质不是一层不变的,随着地理位置的变化还存在着地下物质的变化,更重要的是地下构造也会产生变化,这种变化让人难以捉摸,有时会带领我们走向误区。
所以对圆柱体和倾斜台阶的叠加,是要让我们掌握一种处理,物质界面的方法,更重要的是让我们了解地下物质变化的某些规律
叠加就是两个重力异常的相加,原理与水平圆柱体相同。
公式也一样如下
然后将得到的重力异常相加,我们就可以看到叠加异常了。
2重力异常的数值模拟
前面讲述了,各个地质体模拟的主要步骤和方法,以下是阐述怎样用这些方法来用MATLAB实现重力异常的正演。
2.1规则地质体的数值模拟
2.1.1球体重力异常模拟分析
以下为一个地质体模型,为均匀球体,取R=50m,D=100m,
其模型如下图2-1所示。
D=150m
图2-1
R=60m
以下为二度球体重力异常的计算程序:
clearall
G=6.67*1e-2;
D=150;
R=60;
sigma=1;
x=-200:
10:
200;
M=(4/3)*pi*(R^3)*sigma;
g=(G*M*D)./((x.^2+D^2).^1.5);
plot(x,g,'*-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
程序执行过后,便出现以下图形(图2-2):
图2-2
二度球体重力异常
三度球体重力异常计算程序如下:
clearall
G=6.67*1e-2;
R=60;
D=150;
sigma=1;
x=-200:
10:
200;
y=-200:
10:
200;
M=(4/3)*pi*(R^3)*sigma;
fori=1:
size(x,2)
forj=1:
size(x,2)
g(j,i)=(G*M*D)/((x(i)^2+y(j)^2+D^2).^1.5);
end
end
%三维曲面图
surf(x,y,g);
colorbar;
legend('R=60D=150sigma=1');
boxon;
xlabel('X(m)');
ylabel('Y(m)');
zlabel('\Deltag/g.u.');
%等值线图
contourf(x,y,g,12);
colorbar;
legend('R=60D=150sigma=1');
boxon;
xlabel('X(m)');
ylabel('Y(m)');
在执行上述程序时,可以先将等值线图注释,再执行程序,然后将三维曲面图注释,再执行程序,这样在对话框出来的图片就比较清晰,方便观察。
按照上面方法执行程序过后效果图如以下图2-3、图2-4所示。
图2-3
图2-4
在上述图中,我们可以清晰的看到,球体异常是沿着中心点向四周逐渐减小,而且它的变化率都相同。
2.1.2均匀密度的水平圆柱体异常
1、下面举实际例子,来分析无限长水平圆柱体的重力异常,模型参数为:
半径R=20m,高度D=60m,
=1t/
。
模型图如下图2-5所示
D
图2-5
(1)、按照二度圆柱体重力异常计算公式如下:
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
lamda=pi*(R^2)*sigma;
g=(2*G*lamda*D)./(x.^2+D^2);
plot(x,g,'*-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
执行程序过后,我们可以看到下二度体水平圆柱体异常图(图2-6):
二度水平圆柱体异常图(图2-6)
(2)、三度圆柱体重力异常分析按照一下程序运行:
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
y=-200:
10:
200;
lamda=pi*(R^2)*sigma;
fori=1:
size(x,2)
forj=1:
size(x,2)
g(j,i)=(2*G*lamda*D)/(x(i)^2++D^2);
end
end
%三维曲面图
surf(x,y,g);
colorbar;
boxon;
xlabel('X(m)');
ylabel('Y(m)');
zlabel('Deltag/g.u.');
%等值线图
contourf(x,y,g,12);
colorbar;
axissquare;
boxon;
xlabel('X(m)');
ylabel('Y(m)');
在执行此程序时,我先注释掉后面的等值线图,然后再注释掉曲面图分别运行,运行后,效果图如下图2-7、图2-8所示:
等值线图
三维曲面图
图2-8
图2-7
(3)、无限长水平圆柱体的基本特征
按照下面程序执行:
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
lamda=pi*(R^2)*sigma;
g=(2*G*lamda*D)./(x.^2+D^2);
Vxz=(4*G*lamda*D*x)./(x.^2+D^2).^2;
Vzz=(2*G*lamda*(D^2-x.^2))./(x.^2+D^2).^2;
Vzzz=4*G*lamda*D*(D^2-3*x.^2)./(x.^2+D^2).^3;
subplot(2,2,1);
plot(x,g,'-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
subplot(2,2,2);
plot(x,Vxz,'-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('Vxz/E.');
subplot(2,2,3);
plot(x,Vzz,'-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('Vzz/E.');
subplot(2,2,4);
plot(x,Vzzz,'-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('Vzzz/E.');
执行程序后出下以下图2-9,在图中我们可以看到,重力异常和高阶导数的变化规律。
图2-9
2.2.3均匀密度的倾斜台阶
1、分析下面台阶重力异常。
给定的参数模型数据是:
H=60m,h=20m,倾斜角
为45度,90度,135度,
=1t/
。
模型图如图2-10所示
h
H
图2-10
二度体倾斜台阶模型
重力异常分析图的计算程序如下面所示:
clearall
G=6.67*1e-2;
H=60;
h=20;
sigma=1;
theta=[pi/4pi/23*pi/4];
x=-200:
10:
200;
fori=1:
size(x,2)
forj=1:
size(theta,2)
a1(j,i)=pi*(H-h)+2*H*atan((x(i)+H*cot(theta(j)))/H)-2*h*atan((x(i)+h*cot(theta(j)))/h);
a2(j,i)=x(i)*(sin(theta(j))^2)*log(((H+x(i)*sin(theta(j))*cos(theta(j)))^2+x(i)^2*(sin(theta(j)))^4)/((h+x(i)*sin(theta(j))*cos(theta(j)))^2+x(i)^2.*(sin(theta(j)))^4));
a3(j,i)=-2*x(i)*sin(theta(j))*cos(theta(j))*atan((x(i)*(H-h)*(sin(theta(j)))^2)/(x(i)^2*(sin(theta(j)))^2+(H+h)*x(i)*sin(theta(j))*cos(theta(j))+H*h));
g(j,i)=G*sigma*(a1(j,i)+a2(j,i)+a3(j,i));
end
end
plot(x,g(1,:
),'r:
','LineWidth',2);
holdon
plot(x,g(2,:
),'b-','LineWidth',2);
plot(x,g(3,:
),'g-.','LineWidth',2);
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
h=legend('\alpha=45\circ','\alpha=90\circ','\alpha=135\circ',2);
执行程序过后,结果图如下图2-11所示
图2-11
在图中我们可以看出,倾斜台阶的重力异常是从小到大的增加,然后到达一个稳定值,在每一个角度的最大斜率都相差不大。
但是拐点不同,角度增加拐点的距离变大。
2.2复杂地质体的数值模拟
上面为简单的地质体模型的模拟,下面是稍微复杂一点的重力异常的叠加,在叠加之中我们可以看到许多特点
2.2.1无限长水平圆柱体的重力异常叠加
1、两个无限长水平圆柱体,下面举实际例子,来分析两个无限长水平圆柱体的重力异常,模型参数为:
半径R=20m,高度D=60m,
=1t/
。
参考模型图2-12
D
图2-12
R
(1)、三维曲面图分析,在运行以下程序:
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
y=-200:
10:
200;
lamda=pi*(R^2)*sigma;
fori=1:
size(x,2)
forj=1:
size(x,2)
g1(j,i)=(2*G*lamda*D)/((x(i)+R)^2++D^2);
g2(j,i)=(2*G*lamda*D)/((x(i)-R)^2++D^2);
g=g1+g2;
end
end
surf(x,y,g);
colorbar;
axissquare;
boxon;
xlabel('X(m)');
ylabel('Y(m)');
zlabel('\Deltag/g.u.');
在运行以上程序过后,可以看到如下重力异常三维曲面图2-13:
图2-13
(2)、两个无限长水平圆柱体等势图,在执行以下程序
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
y=-200:
10:
200;
lamda=pi*(R^2)*sigma;
fori=1:
size(x,2)
forj=1:
size(x,2)
g1(j,i)=(2*G*lamda*D)/((x(i)+R)^2++D^2);
g2(j,i)=(2*G*lamda*D)/((x(i)-R)^2++D^2);
g=g1+g2;
end
end
contourf(x,y,g,12);
colorbar;
axissquare;
boxon;
xlabel('X(m)');
ylabel('Y(m)');
运行程序后,效果图片(2-14)如下所示:
图2-14
(3)、两个无限长圆柱体二度体变化曲线
按照下面程序运行。
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
lamda=pi*(R^2)*sigma;
g1=(2*G*lamda*D)./((x+R).^2+D^2);
g2=(2*G*lamda*D)./((x-R).^2+D^2);
g=g1+g2;
plot(x,g,'*-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
运行结果如下图2-15所示:
图2-15
(4)、三个无限长水平圆柱体二度重力异常叠加
按照如下程序:
clearall
G=6.67*1e-2;
R=20;
D=60;
sigma=1;
x=-200:
10:
200;
lamda=pi*(R^2)*sigma;
g3=(2*G*lamda*D)./(x.^2+D^2);
g1=(2*G*lamda*D)./((x+2*R).^2+D^2);
g2=(2*G*lamda*D)./((x-2*R).^2+D^2);
g=g1+g2+g3;
plot(x,g,'*-','MarkerEdgeColor','r');
boxon;
xlabel('X(m)');
ylabel('\Deltag/g.u.');
按照上程序,运行结果图2-16如下:
图2-16
(5)、三个无限长水平圆柱体等势图
在执行以下程序
clearall
G=6.67*1e-2;
R=20;
D1=120;
D2=120
D3=120;
sigma1=1;
sigma2=1;
sigma3=1;
x=-200:
10:
200;
y=-200:
10:
200;
lamda1=pi*(R^2)*sigma1;
lamda2=pi*(R^2)*sigma2;
lamda3=pi*(R^2)*sigma3;
fori=1:
size(x,2)
forj=1:
size(x,2)
g1(j,i)=(2*G*lamda1.*D1)/(x(i)^2++D1.^2);
g2(j,i)=(2*G*lamda2.*D2)/((x(i)-3/2*R)^2++D2.^2);
g3(j,i)=(2*G*lamda3.*D3)/((x(i)+3/2*R)^2++D3.^2);
g=g1+g2+g3;
end
end
contourf(x,y,g,12);
colorbar;
axissquare;
boxon;
xlabel('X(m)');
ylabel('Y(m)');
在执行程序后会出现如下等势图-图2-17
图2-17
(6)、二度体圆柱体的对比,图2-18为一个圆柱体,图2-19为二个圆柱体,图2-20为三个圆柱体。
图2-19
图2-18
图2-20
在上面三幅图中,我们可以清晰的看到从一个到三个圆柱体的二度体重力异常变化,异常最大值在逐渐增大,但是变化率在逐渐减小,这就可以推论,假如有无限个圆柱体,那么它的重力异常变化曲线应该趋近于一条直线,就相当于一个水平地层的重力异常。
(7)、三种等势图的对比,图2-21为一个圆柱,图2-22为两个圆柱。
图2-23为三个圆柱体
图2-21
图2-23
图2-22
从三个等势图片也可以看出,最大值在增加,变化更加缓慢,也就是和上面的斜率减小相互呼应。
3、不同高度的两圆柱之间的叠加重力异常。
与上面的地质模型相似,改变高度
clearall
G=6.67*1e-2;
R=20;
D1=60;
D2=120
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 复杂 质体 重力 异常 计算 分析研究