综合地球物理作业

1、二度体及三度体重力正演姓名：XXX 学号：XXX前言模型正演是地球物理处理中的重要手段，同时也是一些反问题构建的基础，所以正演问题在地球物理问题中具有非常关键的地位。对一些地球物理问题的认识往往是从正演问题开始的。在重力勘探中，正演模拟的地位就非常重要。由于重力勘探采用的是相对测量的手段，所以正问题主要考虑的基本是地下密度异常所引起的相对重力异常。重力异常（重力加速度异常）是密度异常的一种体积效应。如果已知异常体的形状、位置等参数，理论上是可以解析计算出这一异常体引起异常的函数。不过实际实现过程中，这样做会耗费巨大的计算代价。通常的做法是将地下非规则的异常体近似成相近的规则异常体，求解规则异常

2、体的位置、几何参数以及密度；或者对地下空间进行网格剖分，求解每个网格内的密度参数。规则异常体正演的参数空间是异常体中心位置、几何规模（例如正方体模型的边长）、密度差。而基于网格剖分的模型正演的参数空间与所选择的网格大小有关，如果网格确定，那么参数空间就是每个网格里密度差。根据地质实际，地下构造往往在沿着走向上有一定的长度，在走向上构造的变化不大，例如向斜、背斜、断层等构造。对于这样的在走向上有一定延续度并且走向延续度远大于我们研究区域，这时我们称这样的模型为二度体，即我们只需要考虑垂直测线方向和垂直地面这两个空间尺度。但是对于块状地质异常体，例如花岗岩侵入体，这样没有明显构造走向的异常体，我们

3、必须考虑空间三个方向的异常体形态，对于这样的异常体，我们称为三度体。介于二者之间，我们仅考虑构造延伸尺度的异常体，我们称为二度半异常体。这次作业中，我主要考虑二度体和三度体的正演。对于二度体的情况，我考虑使用二度台阶异常的公式模拟一个褶皱界面所引起的重力异常。对于三度体，我主要考虑使用长方体重力正演公式模拟柱状花岗岩体所引起的重力异常。相关公式介绍二度台阶异常体公式：那么二度矩形异常公式：直立长方体模型重力异常公式：其中，i,j,k=1,2。计算代码正确性验证通过理论表达式计算坐标原点的理论异常值和程序数值计算所得结果，可以证明所获得的程序是正确的。二度体正演算例考虑灰岩和砂岩交界的起伏界面所

4、引起的密度异常。给定砂岩密度为，为下伏层，灰岩密度为，为上覆层。那么在界面处，相对水平层状介质有的密度异常。正演结果如下： 三度体正演算例：考虑花岗岩侵入大理岩所引起的密度异常，花岗岩密度：，大理岩密度：。结果如下：附录：二度体：正演函数forward：function g=forward(x,T) %¶þ¶È°å×´ÌåÒì³£×Ó³ÌÐò%UNTITLED Summary of this fu

5、nction goes here% Detailed explanation goes here x0=T(1);%Öеãºá×ø±ê a=T(2);%°ë¿í¶È H=T(3);%ÉÏµ× h=T(4);%ÏÂµ× Q=T(5);%ÃܶȲî f=6.67*10(-11);%Ò

6、53;Á¦³£Êý for i=1:length(x) %×ǫ́½×Òì³£ X(i)=x(i)-(x0-a); p1(i)=pi*(H-h); p21(i)=(X(i)2+H2)/(X(i)2+h2); p2(i)=X(i)*log(p21(i); p3(i)=2*H*atan(X(i)/H); p4(i)=-2*h*atan(X(i)/h); g1(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i);%À&#

7、237;ÂÛ¼ÆË㹫ʽ %ÓǪ̀½×Òì³£ X(i)=x(i)-(x0+a); p1(i)=pi*(H-h); p21(i)=(X(i)2+H2)/(X(i)2+h2); p2(i)=X(i)*log(p21(i); p3(i)=2*H*atan(X(i)/H); p4(i)=-2*h*atan(X(i)/h); g2(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i);%&

8、#192;íÂÛ¼ÆË㹫ʽ g(i)=g1(i)-g2(i);%¶þ¶È°å×´ÌåÒì³££¨×ǫ́½×ÓëÓǪ̀½×Òì³£²

9、 g(i)=106*g(i);%½«ÖØÁ¦Òì³£µÄµ¥Î»×ª»¯³É¸ß˹ endend调用程序s1：clear allclose allx=0:5:400; %¹Û²âµãºá×ø±

10、êz=0; %¹Û²âµã×Ý×ø±êT=200,20,400,10,0.3;%¸ø³öÕýÑÝÄ£ÐͲÎÊýTg=forward(x,T);%µ÷ÓÃÖØÁ¦Òì³£×

11、1;³ÌÐòfigure(1)plot(x,g,'.-r') %»³öÖØÁ¦Òì³£ylabel('Deltag/g.u.','FontSize',12)xlabel('²âµãºá×ø±ê/m','FontSize',12)mo=zeros(length(x),le

12、ngth(x);for j=(200-20)/5+1:(200+20)/5 for i=3:length(x) mo(i,j)=500; endendfigure(2)pcolor(x,x,mo)axis ij调用程序s2：clear allclose allx=-2000:20:2000; %¹Û²âµãºá×ø±êz=0; %¹Û²âµã×Ý×ø±êh0

13、=2000;A=1000;mx=-1000:1:1000;y=zeros(1,length(x);for i=1:length(mx) H(i)=h0+A*sin(i-1)/1000*pi); T(1)=mx(i); T(2)=0.5; T(3)=H(i); T(4)=h0; T(5)=-0.3; g=forward(x,T); y=y+forward(x,T); %if (mod(i,10)=0) %figure(i) %plot(x,g,'.-g') %ylabel('Deltag/g.u.','FontSize',12) %xlabel(&

14、#39;²âµãºá×ø±ê/m','FontSize',12) %endendfigure(10000)plot(x,y,'.-r') %»³öÖØÁ¦Òì³£ylabel('Deltag/g.u.','FontSize',12)xlabel('²âµãº&

15、#225;×ø±ê/m','FontSize',12)figure(100000)plot(mx,H,'.-r') %»³öÉî¶È½çÃæylabel('½çÃæÉî¶È/m','FontSize',12)xlabel('²âµãº

16、25;×ø±ê/m','FontSize',12)三度体：正演程序chang：function data_c=chang(dm,x1,x2,y1,y2,z1,z2)xx=zeros(2,1);yy=zeros(2,1);zz=zeros(2,1);xx(1)=x1;xx(2)=x2;yy(1)=y1;yy(2)=y2;zz(1)=z1;zz(2)=z2;sum=0;for ii=1:2 for jj=1:2 for kk=1:2 ro=(xx(ii)2+yy(jj)2+zz(kk)2)0.5; sum=sum+(-1)(ii+jj+kk)*(xx(ii)*log(yy(jj)+ro)+yy(jj)*log(xx(ii)+ro)+zz(kk)*(pi+atan2(zz(kk)*ro,xx(ii)*yy(jj); end endenddata_c=-6.67*dm*sum/1000;end调用程序changf：clearda=20;a=0:da:5000;middle_x=2500;middle_y=2500;dm=0.5;dl=20;h0=100;data=zeros(length(a)