物探数据处理实验报告

本科生实验报告

实验课程物探数据处理

学院名称地球物理学院

专业名称勘查技术与工程

学生姓名00000

学生学号0000000000000

指导教师李勇

实验地点5417

实验成绩

二〇一七年九月二〇一七年十月

物探数据处理实验报告

实验一

1.1实验目的

本次实验的目的是将课本上给出的地质体情况，和计算参量的代数公式用c 语言通过程序表现出来。在得知已知参量，例如形状大小、空间位置的情况下将其带入编写好的程序，从而得出地下板状体的物性值数据，并将所得数据其通过成图数据成图后在用反演原理与提前已知的异常体情况对比，判断所得结果的正误，并分析最终的实验结果。

1.2实验原理

下图1-1为实验原理示意图。

图1-1 二维板模型

利用如下公式进行c语言程序的实现：

1.3实验参数

中点横坐标x0(m)：1000 中点深度坐标z0(m)：1000 板高度(m)：400 板长度(m)：200 磁化强度：2000 二维板角度：45、90 磁化角度：90 测点k：100

}

1.4 实验结果图示

图1-2 磁异常X分量

图1-3 磁异常Z分量

图1-4 重力异常

1.5 实验结果分析

根据实验数据结果图可知，通过观察x方向的磁分量可以看出：图中存在着明显的两个相互对称的正负异常，可知在地质体的上部积累了一定的负磁荷，而在地质体的下部积累了同样大小的正磁荷，实验中磁化角度取90度，二维板倾角取不同的角度。通过观察z方向上的磁异常可知：地质体呈现出中部正异常的对称图形，总体呈现先增后减的趋势。通过观察重力异常的图像可知：得出的结果基本与已知相符，所以通过重力异常、磁场x方向异常、磁场z方向异常数据曲线，就可对所求数据做相应的分析。

实验二

2.1 实验目的

通过给定的已知相关物性参数，建立截面为多边形的水平柱体模型，带入所编写好的c 语言程序中得到对应的正演结果。主要考查了c 语言程序进编程，和对所得结果的理解，分析，判断能力。

2.2 实验原理

图2-1为实验原理图

图2-1 多边形截面水平柱体截面形状

根据下列公式进行编程得到实验数据，作图分析。

；

]1

111[)1(2)1(2)

11)(1(22121

2ln )1(2)1(2)11)(1(21{12x i z

i tg x i z i tg x i x i z i z i z i x i z i x i x i x i z i x i z i x i x i x i z i z i z i x i z i x i z i z i n i G g ++---?-++-++-+-++++++?-++-++-+-+?=∑=?σ；)(21

Q M z P M x X +=?π ；)(21

P M z Q M x Z -=

?π

式中：

；

z i

x i z i x i x i x i z i z i z i z i x i z i tg x i z i tg x i x i z i z i x i x i z i z i n i Q 2212

12ln )1(2)1(2)1(221]1111[)1(2)1(2)1)(1({1++++?+-+-+-+?

-++---?+-+-++--+=∑=；

z i

x i z i x i x i x i z i z i x i x i z i z i x i z i tg x i z i tg x i x i z i z i z i z i n i P 2212

12ln )1(2)1(2)1

()1(21]1111[)1(2)1(2)1(2{1++++?+-+-++--+?+++---?+-+-+-+=∑=

2.3 实验参数

多边形半径r:200m ；磁化强度M:2000 拟合边数:10 磁化角度ct:90,120 测点k:100 中点深度z0(m):1000 中点横坐标x0(m):1000

2.3.1 讨论a

b

tg 1-的值

a

b

实际为由坐标原点对多边形第i 边的夹角，在具体计算中应当进行讨论： 当0.115->a 时，如果a b >0.0；当a>0.0,b>0.0,则θ=a b

tg 1-；当a<0.0，b<0.0，

则θ=—π+a

b

tg 1-；

如果a b <0.0，则θ=π+a

b tg 1-；

当0.115-≤a 时，若b<0.0，则θ=—2π；若b>=0,则θ=2π

2.4 实验结果

图2-2 重力异常

图2-3 磁异常X分量

图2-4 磁异常Z分量

2.5 实验结果分析

通过重力异常图可知，通过观察x方向的磁分量可以看出：图中仍然存在着明显的两个相互对称的正负异常，故我们仍然可以认为地质体的上部积累了一定的负磁荷，而在地质体的下部积累了同样大小的正磁荷。所以左端出现了较大的正异常，右端出现了较大的正异常。通过观察z方向上的磁异常可知：地质体在中线（x=50）呈现出中部正异常的对称图形，总体呈现先增后减的趋势。通过观察重力异常的图像可知：得出的结果基本与已知相符，出现了一个重力异常值曲线，最大值点对应的测线坐标与地质体实际位置一致。

实验三

3.1 实验目的

本次实验的目的是通过给定的地下异常体参量，对地下的正长方体进行正演计算，然后将计算参量的代数公式用c 语言通过程序表现出来。在得知已知参量，例如形状大小、空间位置的情况下将其带入编写好的程序，并通过积分来实现正演。

从而得出地下板状体的物性值数据，并将所得数据其通过成图数据成图后在用反演原理与提前已知的异常体情况对比，判断所得结果的正误，并分析最终的实验结果。将其运用的主要方法有通过划分单元来进行体积分，即8个体积元叠加即为立方体所对应的积分值，偶数为正，奇数为负，则可以计算出相应重磁异常。

3.2 实验原理

利用如下公式进行编程：

；2

02020202

02

01ln ln ({c z c z b y b y a x a

x }

zk)R

(z yk)xk)(y (x tg zk)(z R]xk)[(x yk)(y R]yk)[(y xk)x G Δg -

+

-+-

+

----?-++--++---=σ；2

020********yk)]}

-(y ln[R M z zk)]-(z ln[R M y x k)R

-(x zk)-yk)(z -(y tg 1-M x {-41c z c

z b y b y a x a x Δx -

+

-+-+

++++=π；2

020********x k)]}

-(x ln[R M z ]yk)R

-(y zk)-x k)(z -(x tg 1-M y -zk)]-(z ln[R M x {41c z c z b y b y a x a x Δy -

+

-+-+

++++=π

；2

020********}

zk)R

-(z yk)-x k)(y -(x tg 1-M z x k)]-(x ln[R M y -yk)]-(y ln[R M x {41c z c z b y b y a x a x Δz -

+

-+-+

-++=π式中：])(2)(2)(2[2

/1zk z yk y xk x R -+-+-=；

3.3 实验参数

列数:101 行数: 101 磁化强度:2000 磁倾角:45 磁偏角:45 立方体上坐标:（300,300,10） 立方体下坐标:（700,700,410） 密度:1 起终坐标:（0,1000）

3.4 实验结果图示

图3-1 重力异常平面图 图3-2 磁异常x 分量平面图

3-3 磁异常y 分量平面 图3-4 磁异常z 分量平面图

3.5实验结果分析

根据实验数据结果图可知，通过观察x方向的磁分量可以看出：图中存在着两个对称的正负异常，且异常呈现出相似性形态幅值基本一致。通过观察z方向上的磁异常可知：地质体呈现出中部正异常的对称图形，切呈现出对角线对称的形态上方为正下方为负的负异常。通过观察重力异常的图像可知：得出的结果基本与已知相符，得到一个明显的类似同心圆的重力异常，边界较为清晰。

实验程序1

#include

#include

#define PI 3.1415926

#define G 6.67e-6

void main()

{

FILE *fp1,*fp2,*fp3;

fp1=fopen("g.txt","w");

fp2=fopen("deltax.txt","w");

fp3=fopen("deltaz.txt","w");

double r1,r2,r3,r4; //定义板截面角点A,B,C,D到P点之间的距离

double f1,f2,f3,f4; //定义r1.r2.r3.r4与X轴正向夹角

int k,i; //定义P点的坐标xk,zk;

double xk,zk;

double cgm=2.67e+3;

double x0=1000.0; //中点横坐标

double z0=1000.0;

double l=400.0; //板高度

double b=200.0; //板长度

double I,M=2000.0;

double z1,z2,g,deltax,deltaz;

double J=90.0*PI/180.0; //计算二维板角度

I=90.0*PI/180.0; //计算磁化角度

double a[100]={0.0};double b1[100]={0.0};double ct[100]={0.0};

for(k=0;k<101;k++)

{

xk=k*20.0;

zk=0.0;

b1[0]=z0-zk-l*sin(J);

a[0]=xk-x0+b+l*cos(J);

b1[1]=z0-zk+l*sin(J);

a[1]=xk-x0+b-l*cos(J);

b1[2]=z0-zk-l*sin(J);

a[2]=xk-x0-b+l*cos(J);

b1[3]=z0-zk+l*sin(J);

a[3]=xk-x0-b-l*cos(J);

r1=sqrt(pow(b1[0],2.0)+pow(a[0],2.0));

r2=sqrt(pow(b1[1],2.0)+pow(a[1],2.0));

r3=sqrt(pow(b1[2],2.0)+pow(a[2],2.0));

r4=sqrt(pow(b1[3],2.0)+pow(a[3],2.0));

for(i=0;i<4;i++)

{

if(fabs(a[i])>1.0e-15)

{

if((b1[i]/a[i])>0.0)

{

if(a[i]>0.0&&b1[i]>0.0)

{

ct[i]=atan(b1[i]/a[i]);

}

else if(a[i]<0.0&&b1[i]<0.0)

{

ct[i]=atan(b1[i]/a[i])-PI;

}

}

else if((b1[i]/a[i])<0.0)

{

ct[i]=PI+atan(b1[i]/a[i]);

}

}

else if(fabs(a[i])<=1.0e-15)

{

if(b1[i]<0.0)

{

ct[i]=-PI/2.0;

}

else

{

ct[i]=PI/2.0;

}

}

f1=PI-ct[0];

f2=PI-ct[1];

f3=PI-ct[2];

f4=PI-ct[3];

z1=z0-l*sin(J);

z2=z1+l*sin(J);

g=2.0*G*cgm*((z2*(f2-f4)-z1*(f1-f3))+xk*(sin(J)*sin(J)*log((r2*r3)/(r 1*r4))+cos(J)*sin(J)*(f1-f2-f3+f4))+2.0*b*(sin(J)*sin(J)*log(r4/r3)+c os(J)*sin(J)*(f3-f4)));

deltax=((M*sin(J))/(2.0*PI))*(log((r2*r3)/(r1*r4))*cos(J-I)-(sin(J-I) )*(f1-f2-f3+f4));

deltaz=((M*sin(J))/(2.0*PI))*((sin(J-I))*log((r2*r3)/(r1*r4))+(cos(J-I))*(f1-f2-f3+f4));

fprintf(fp1,"%f\n",g); //将所得数据写入文件

fprintf(fp2,"%f\n",deltax);

fprintf(fp3,"%f\n",deltaz);

}

}

fclose(fp1);

fclose(fp2);

fclose(fp3);

实验程序2

include"stdio.h"

#include"math.h"

#define PI 3.1415926 //圆周率

#define G 6.67e-11 // 万有引力常数

#define r 200.0 //多边形半径

#define M 2000.0 //磁化强度

#define thegama 2.67e+3

void main()

{

double ct11(double a,double b); //声明子函数

FILE *fp1,*fp2,*fp3;

fp1=fopen("deltag.txt","w");

fp2=fopen("deltax1.txt","w");

fp3=fopen("deltaz1.txt","w");

double x0=1000.0;

double z0=1000.0;

int n=10,N=100;

double R=200.0;

double M=2000.0;

double cgm=2.67e+3;

double ct=PI/4.0;

double Mx=M*cos(ct);

double Mz=M*sin(ct);

double x[10000]={0.0},x1[1000]={0.0};

double z[10000]={0.0},z1[1000]={0.0};

double xk[10000]={0.0},zk[1000]={0.0};

double sum1=0.0;

double Q=0.0;

double P=0.0;

double

a1,a2,a3,a4,a5,a6,a7,aa,bb,ct1,ct2,deltag,deltax1,deltaz1,a,b;

int i,k;

for(i=0;i

{

x1[i]=x0+R*cos(i*2*PI/n);

z1[i]=z0-R*sin(i*2*PI/n);

}

for(k=0;k

{

xk[k]=20.0*k;///////////////////平面坐标

zk[k]=0.0;

}

for(k=0;k

{

sum1=0.0;

Q=0.0;

P=0.0;

for(i=0;i

{

x[i]=xk[k]-x1[i];

z[i]=zk[k]-z1[i];

}

for(i=0;i

{

a1=(z[i+1]-z[i])*(x[i]*z[i+1]-x[i+1]*z[i])/(pow(z[i+1]-z[i],2.0)+ pow(x[i+1]-x[i],2.0));

a2=(x[i+1]-x[i])*(x[i]*z[i+1]-x[i+1]*z[i])/(pow(z[i+1]-z[i],2.0)+pow( x[i+1]-x[i],2.0));

a3=(z[i+1]-z[i])*(x[i]-x[i+1])/(pow(z[i+1]-z[i],2.0)+pow(x[i]-x[i +1],2.0));

a4=pow(z[i+1]-z[i],2.0)/(pow(z[i+1]-z[i],2.0)+pow(x[i]-x[i+1],2.0));

aa=pow(x[i+1],2.0)+pow(z[i+1],2.0);

bb=pow(x[i],2.0)+pow(z[i],2.0);

ct1=ct11(z[i],x[i]); //调用子函数

ct2=ct11(z[i+1],x[i+1]);

a5=(a1/2.0)*log(aa/bb)+a2*(atan(ct1)-atan(ct2));

sum1+=a5;

a6=a3*(atan(ct1)-atan(ct2))-(a4/2.0)*log(aa/bb);

Q+=a6;

a7=a4*(atan(ct1)-atan(ct2))+(a3/2.0)*log(aa/bb);

P+=a7;

}

deltag=2.0*G*cgm*sum1;

deltax1=(Mx*P+Mz*Q)/(2.0*PI);

deltaz1=(Mx*Q-Mz*P)/(2.0*PI);

fprintf(fp1,"%lf\n",deltag);

fprintf(fp2,"%lf\n",deltax1);

fprintf(fp3,"%lf\n",deltaz1);

}

fclose(fp1);

fclose(fp2);

fclose(fp3);

}

///////////////////////////////////////////

double ct11(double a,double b) //判断子函数{

double ct;

if(fabs(a)>1.0e-15)

{

if((b/a)>0.0)

{

if(a>0.0&&b>0.0)

ct=atan(b/a);

else

ct=atan(b/a)-PI;

}

if((b/a)<0.0)

{

ct=PI+atan(b/a);

}

}

else

{

if(b<0.0)

{

ct=-PI/2.0;

}

else

{

ct=PI/2.0;

}

}

return(ct);

}

实验程序3

#include

#include

#include

#define l 101 // 列数

#define h 101 // 行数

#define M 2000.0 //磁化强度

#define II 60 // 磁倾角

#define aa 60 //磁偏角

#define PI 3.1415926 //圆周率

#define G 6.674E-11 // 万有引力常数

#define SI 1e6 // 将m/s^2换算为g.u.的比例因子

#define DENSITY 1E3 // 将g/cm^3换算为kg/m^3的比例因子

//计算atan值

double batan(double m,double n,double x,double r);

//Kernel函数,用途：为Cuboid函数准备

double Kernel(double dXk, double dYk, double dZk,

//dXk, dYk, dZk: 观测点坐标

double dXs, double dYs, double dZs);

//dYs, dZs, dZs: 立方体坐标，s取1,2 //Cuboid函数,用途: 计算立方体重力异常,参数说明

double Cuboid(double dXk, double dYk, double dZk,

//dXk, dYk, dZk: 观测点坐标

double dX1, double dY1, double dZ1,

// dX1, dY1, dZ1: 立方体坐标1

double dX2, double dY2, double dZ2,

//dX2, dY2, dZ2: 立方体坐标2

double dDens); //dDens: 立方体密度

// Cix函数

double Cix(double dXk, double dYk, double dZk,

double dXs, double dYs, double dZs,

double Mx, double My, double Mz);

// Cixxx函数: 用于计算立方体△x异常

double Cixxx(double dXk, double dYk, double dZk,

double dX1, double dY1, double dZ1,

double dX2, double dY2, double dZ2,

double Mx, double My, double Mz );

// Ciy函数

double Ciy(double dXk, double dYk, double dZk,

double dXs, double dYs, double dZs,

double Mx, double My, double Mz);

// Ciyyy函数: 用于计算立方体△Y磁异常

double Ciyyy(double dXk, double dYk, double dZk,