1、课程设计-长方体的重磁异常正演计算学 院：地球物理学院专 业：地球物理学指导老师：姓 名：学 号：一、 序言野外重磁测量的最后成果是重磁异常的等值线平面图和剖面平面图,这些图件反映了待查的地下目标物与围岩之间的密度及磁性差异所产生的重磁异常特征。为了解各种地质现象及目标物与重磁异常的对应规律和本质联系以及重磁异常特征与各种地质(目标物)体形状、产状等的定性和定量关系，以便根据测得的重磁异常推断出地下的地质情况及目标物的分布,因此,要根据实际地质模型简化出的规则密度体及磁性体的重磁场进行数学物理的解析，通过分析对比计算出的各种规则形体的重磁异常剖面和平面上的特征,从中找出其规律，为所获资料的解释
2、推断建立感性认识,寻求数理依据。这种根据已知质体及磁性体的形态、质量及磁性、空间等分布来计算其重磁场分布的过程，称为重磁场正演问题二、 原理设有 x、y、z 坐标系，其 x、y 平面与地面重合，z 轴锤子向下，原点为O，用三组平行于坐标面的截平面去切某形体，就得到形状规则的长方体，今研究其长方体在地面引起的重磁异常。它是一个典型的三度体，由计算三度体重磁异常的一般公式出发。其中三度体的 g 异常的一般公式出发，g=G3222()()zdxyz现计算 xy 平面内测点 P（ x、y、o）的 g 值，则只有 z=0，如果将计算点P 作为坐标原点（即移动的坐标原点） ，则有 x=y=z=0，代入上式
3、有32221()xyzdgG经过积分，计算长方体 g 异常值得解析表达式为12ln()l()|xyzVRtgz R 同理由函数的对称性可得：21ln()l()|xyzGtgx 21l()l()|xyzVRty R式中： 22Rkkyxz一个体积为 v、密度均匀的物体之引力位为1vVGdr式中 G 为引力常数， 为密度差。同一均匀磁化物体的磁位为114pvUMgradA综合上两式，可得：14pgradVGA上式即为磁位与引力位之间的泊松公式。该式表明，同一个既均匀磁化有密度均匀物体的磁位，可由其引力位来计算。对于正长方体的情况，它的磁位可表示为21()4zyxUdzrA磁异常的各分量就是磁位在各
4、分量的导数，即：211lnln|4 xyzxyzXMtgRMx 211ll|xyzyxzUYty211lnln|4 xyzzxyZtgRz 式中： 22Rkkyxz注意：1.式中凡长度均用 m 作单位。 为剩余密度； g 以 mgal 为单位；（ g/c） G 为万有引力常数，其值 G=0.0067， 为磁化强度在 x 方向分量， 为磁化MxMy强度在 y 方向分量， 为磁化强度在 z 方向分量。角度取值范围 ，而不z 0是取主值 。22.公式推导是把计算点作为原点，即计算点 x=y=z=0，而实际上，模型和计算点是在统一坐标系下，计算点 xp、 yp、 zp 均是有值的，所以，实际计算公式中
5、的代限值应为：x1=Xp-X1 x2=Xp-X2y1=Yp-Y1 y2=Yp-Y2z1=Yp-Z1 z2=Yp-Z2其中 X1,X2,Y1,Y2,Z1,Z2 为长方体各顶点坐标。 Xp， Yp， Zp 为计算点坐标。三、 程序本程序采用的模型参数为1000.0,1000.0,1000.0,100.0,0.0,100.0,500.0,500.0,500.0其中正长方体的中心点坐标为（1000,1000,1000） ，三个磁化方向的分量为Mx=100， My=0， Mz=100，本程序采用的是正方体，其长宽高均为 500 米。这些数据通过文件输入的方式进行；剩余密度取 2.67g/cm,测线距与测
6、点距均设置为 20 m；本程序设置了 101 条测线，每条测线上有 101 个测点。其输出均为文件输出，均为 txt 格式文件。程序代码：#include#include#define pi 3.1416#define G 6.67E-3float atg(float z,float x) /定义反正切子函数float c; if(fabs(x)(1E-15) if(x0if(x0)c=(pi/2);else if(z=0)c=(-pi/2);return c; void main()FILE *fp,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;fp=fopen(“正长方体重力
7、异常.txt“,“w“);fp1=fopen(“正长方体磁异常在 X 方向分量.txt“,“w“);fp2=fopen(“正长方体磁异常在 Y 方向分量.txt“,“w“);fp3=fopen(“正长方体磁异常在 Z 方向分量.txt“,“w“);fp4=fopen(“正长方体重力位在 X 方向导数.txt“,“w“);fp5=fopen(“正长方体重力位在 Y 方向导数.txt“,“w“);fp6=fopen(“模型体参数.txt“,“r“);float x0,y0,z0,a,b,c,q,R,dg,dGx,dGy,DX,DY,DZ,x1,y1,z1,Mx,My,Mz,x2,y2,z2,g10
8、1,X101,Y101,Z101,Gx101,Gy101;float p,dx,dy,xk,yk; fscanf(fp6,“%f,%f,%f,%f,%f,%f,%f,%f,%f“,printf(“请输入模型的剩余密度，测点距以，测线距n“); /输入模型数据scanf(“%f,%f,%f“,x0=x0-a/2; /计算长方体各个顶点的坐标x1=x0+a/2;y0=y0-b/2;y1=y0+b/2;z0=z0-c/2;z1=z0+c/2;int i,j,k,m,n;for(i=0;i101;i+)for(j=0;j101;j+)gj=0.0;Xj=0.0;Yj=0.0;Zj=0.0;Gxj=0.
9、0;Gyj=0.0;for(n=0;n2;n+) /嵌套三重循环语句for(m=0;m2;m+)for(k=0;k2;k+)q=pow(-1.0,k+m+n); /判断每一变量代限结果的正负xk=j*dx;yk=i*dy;x1=xk-xk;y1=ym-yk;z1=zn;R=sqrt(x1*x1+y1*y1+z1*z1);dg=(-G*p*q*(x1*log(y1+R)+y1*log(x1+R)+(z1*atg(x1*y1),(z1*R);gj=gj+dg;dGx=(-G*p*q*(z1*log(y1+R)+y1*log(z1+R)+(x1*atg(z1*y1),(x1*R);Gxj=Gxj+d
10、Gx;dGy=(-G*p*q*(z1*log(x1+R)+x1*log(z1+R)+(y1*atg(z1*x1),(y1*R);Gyj=Gyj+dGy;DX=(q/(4*pi)*(-Mx*atg(y1*z1),(x1*R)+My*log(R+z1)+Mz*log(R+y1);Xj=Xj+DX;DY=(q/(4*pi)*(-My*atg(x1*z1),(y1*R)+Mx*log(R+z1)+Mz*log(R+x1);Yj=Yj+DY;DZ=(q/(4*pi)*(-Mz*atg(x1*y1),(z1*R)+Mx*log(R+y1)-My*log(R+x1);Zj=Zj+DZ;fprintf(fp,
11、“%f %f %fn“,xk,yk,gj);fprintf(fp1,“%f %f %fn“,xk,yk,Xj);fprintf(fp2,“%f %f %fn“,xk,yk,Yj);fprintf(fp3,“%f %f %fn“,xk,yk,Zj);fprintf(fp4,“%f %f %fn“,xk,yk,Gxj);fprintf(fp5,“%f %f %fn“,xk,yk,Gyj);fclose(fp);fclose(fp1);fclose(fp2);fclose(fp3);fclose(fp4);fclose(fp5);fclose(fp6);四、 结果图图 1、正方体的重力异常图图 2、重力位 X 方向导数异常图图 3、重力位 Y 方向导数异常图图 4、磁异常在 X 方向分量图 5、磁异常在 Y 方向的分量图 6、磁异常在 Z 方向的分量
