直立长方体重力异常及其一阶偏导数

2016年1月24日 354 次阅读 0 条评论 0 人点赞

这是我的一份实验报告 

 

 

 

本科生实验报告

 

 

填写说明

 

  1. 适用于本科生所有的实验报告(印制实验报告册除外);
  2. 专业填写为专业全称,有专业方向的用小括号标明;
  3. 格式要求:
    1. A4纸双面打印(封面双面打印)或在A4大小纸上用蓝黑色水笔书写。
    2. 打印排版:正文用宋体小四号,1.5倍行距,页边距采取默认形式(上下2.54cm,左右2.54cm,页眉1.5cm,页脚1.75cm)。字符间距为默认值(缩放100%,间距:标准);页码用小五号字底端居中。
      1. 具体要求:

      题目(二号黑体居中);

      摘要("摘要"二字用小二号黑体居中,隔行书写摘要的文字部分,小4号宋体

      关键词(隔行顶格书写"关键词"三字,提炼3-5个关键词,用分号隔开,小4号黑体)

  4. 正文部分采用三级标题;
  5. 1
    ××(小二号黑体居中,段前0.5)
  6. 1.1 ×××××小三号黑体×××××(段前、段后0.5行)
  7. 1.1.1小四号黑体(段前、段后0.5行)
    1. 参考文献(黑体小二号居中,段前0.5行),参考文献用五号宋体,参照《参考文献著录规则(GB/T 77142005)》。

 

直立长方体重力异常及其一阶偏导数

  1. 理论公式

模型图

 


(xs,ys,zs)为长方体上任意点,(x,y,z)为测点。r    长方体上某点到测点的距离。

 




 



所以剩余密度为σ直立长方体重力异常计算公式为

    

导数计算公式采用差分近似计算的方法,公式为

  1. 程序代码

1.主函数

clear all;

clc;

X=0:10:1000;

Y=0:10:1000;%自变量

X1=X-1;

Y1=Y-1;%求导公式差分项变量

uiCols=size(X,2);

uiRows=size(Y,2);%自变量大小

G=zeros(uiRows,uiCols);%定义重力异常数组

G1=zeros(uiRows,uiCols);

G2=zeros(uiRows,uiCols);%垂向导数差分项

Gx=zeros(uiRows,uiCols);%X方向一阶倒数数组

Gy=zeros(uiRows,uiCols);%Y方向一阶倒数数组

Gz=zeros(uiRows,uiCols);%垂向一阶倒数

GG=zeros(uiRows,uiCols);

dX1=400; dX2=600;

dY1=400; dY2=600;

dZ1=10; dZ2=dZ1+100;

dDensity=1.0;

for i=1:1:uiRows;

for j=1:1:uiCols;

G(i,j)=Kernel(X(j),Y(i),0,dX2,dY2,dZ2,dDensity)...

-Kernel(X(j),Y(i),0,dX2,dY2,dZ1,dDensity)...

-Kernel(X(j),Y(i),0,dX2,dY1,dZ2,dDensity)...

+Kernel(X(j),Y(i),0,dX2,dY1,dZ1,dDensity)...

-Kernel(X(j),Y(i),0,dX1,dY2,dZ2,dDensity)...

+Kernel(X(j),Y(i),0,dX1,dY2,dZ1,dDensity)...

+Kernel(X(j),Y(i),0,dX1,dY1,dZ2,dDensity)...

-Kernel(X(j),Y(i),0,dX1,dY1,dZ1,dDensity);

end;

end;%重力异常

for i=1:1:uiRows;

for j=1:1:uiCols;

Gx(i,j)=Kernel(X1(j),Y(i),0,dX2,dY2,dZ2,dDensity)...

-Kernel(X1(j),Y(i),0,dX2,dY2,dZ1,dDensity)...

-Kernel(X1(j),Y(i),0,dX2,dY1,dZ2,dDensity)...

+Kernel(X1(j),Y(i),0,dX2,dY1,dZ1,dDensity)...

-Kernel(X1(j),Y(i),0,dX1,dY2,dZ2,dDensity)...

+Kernel(X1(j),Y(i),0,dX1,dY2,dZ1,dDensity)...

+Kernel(X1(j),Y(i),0,dX1,dY1,dZ2,dDensity)...

-Kernel(X1(j),Y(i),0,dX1,dY1,dZ1,dDensity);

end;

end;%x方向一阶倒数

for i=1:1:uiRows;

for j=1:1:uiCols;

Gy(i,j)=Kernel(X(j),Y1(i),0,dX2,dY2,dZ2,dDensity)...

-Kernel(X(j),Y1(i),0,dX2,dY2,dZ1,dDensity)...

-Kernel(X(j),Y1(i),0,dX2,dY1,dZ2,dDensity)...

+Kernel(X(j),Y1(i),0,dX2,dY1,dZ1,dDensity)...

-Kernel(X(j),Y1(i),0,dX1,dY2,dZ2,dDensity)...

+Kernel(X(j),Y1(i),0,dX1,dY2,dZ1,dDensity)...

+Kernel(X(j),Y1(i),0,dX1,dY1,dZ2,dDensity)...

-Kernel(X(j),Y1(i),0,dX1,dY1,dZ1,dDensity);

end;

end;%y方向一阶倒数

for i=1:1:uiRows;

for j=1:1:uiCols;

G1(i,j)=Kernel(X(j),Y(i),1,dX2,dY2,dZ2,dDensity)...

-Kernel(X(j),Y(i),1,dX2,dY2,dZ1,dDensity)...

-Kernel(X(j),Y(i),1,dX2,dY1,dZ2,dDensity)...

+Kernel(X(j),Y(i),1,dX2,dY1,dZ1,dDensity)...

-Kernel(X(j),Y(i),1,dX1,dY2,dZ2,dDensity)...

+Kernel(X(j),Y(i),1,dX1,dY2,dZ1,dDensity)...

+Kernel(X(j),Y(i),1,dX1,dY1,dZ2,dDensity)...

-Kernel(X(j),Y(i),1,dX1,dY1,dZ1,dDensity);

end;

end;

for i=1:1:uiRows;

for j=1:1:uiCols;

G2(i,j)=Kernel(X(j),Y(i),-1,dX2,dY2,dZ2,dDensity)...

-Kernel(X(j),Y(i),-1,dX2,dY2,dZ1,dDensity)...

-Kernel(X(j),Y(i),-1,dX2,dY1,dZ2,dDensity)...

+Kernel(X(j),Y(i),-1,dX2,dY1,dZ1,dDensity)...

-Kernel(X(j),Y(i),-1,dX1,dY2,dZ2,dDensity)...

+Kernel(X(j),Y(i),-1,dX1,dY2,dZ1,dDensity)...

+Kernel(X(j),Y(i),-1,dX1,dY1,dZ2,dDensity)...

-Kernel(X(j),Y(i),-1,dX1,dY1,dZ1,dDensity);

end;

end;%垂向导数差分计算

Gz=(G1-G2)/2;%垂向导数

GG=sqrt((G-Gx).^2+(G-Gy).^2);

figure

subplot(2,2,1);%划分画板

contour(X,Y,G,50);%重力异常

subplot(2,2,2);

contour(X,Y,(G-Gx),50);%X方向倒数

subplot(2,2,3);

contour(X,Y,(G-Gy),50);%Y方向倒数

subplot(2,2,4);

contour(X,Y,Gz,50);%垂向导数

figure (2);

contour(X,Y,GG,50);

2.子函数

function dOutput=Kernel(dXo,dYo,dZo,... % 观察点坐标

dXs,dYs,dZs,... % 场源坐标

dDensity) % 密度

GC=6.67E-11;%万有引力常量

CGS2SI=1E3;

SI2MGAL=1E5;

dX=dXs-dXo;

dY=dYs-dYo;

dZ=dZs-dZo;

dR=sqrt(dX^2+dY^2+dZ^2);

dDensity=dDensity*CGS2SI;

dOutput=-dX*log(dR+dY)-dY*log(dR+dX)+dZ*atan(dX*dY/(dZ*dR));

dOutput=dOutput*GC*dDensity*SI2MGAL;

  1. 计算结果

 



 

 

 

 

 

 

 

  1. 结果分析

通过观察所得图像可以看出因为立方体为规则图形,所以重力异常图呈对称分布。X方向和Y方向一阶倒数变化最快的地方是地质体的边界.。在水平方向,地质体的边界一共有四个。X方向的一阶倒数的两个极值点分别是地质体在水平X方向的两个边界,正异常为左边界负异常为右边界。Y方向水平一阶导数的两个极大值分别是水平Y方向的两个边界,正异常为下边界,负异常为上边界。垂向导数也可以判断地质体的边界,也可以描述地质体的形态。

 

学生实验 心得

这次试验根据直立长方体重力异常计算公式和利用差分法近似计算直立长方体水平导数和垂向导数,运用matlab编写程序,计算重力异常和导数,并绘制图像。差分法在这个程序中可以简单地实现,只需在原来的计算重力异常的程序进行简单的修改即可。另外自变量间隔不宜取得太小,会使程序无法继续运行下去。Matlab功能强大,但是计算效率不算太高,无法应用于计算量非常大的任务。

 

 

 

学 生(签 名):

指导

教师

评语

 

 

 

 

 

 

 

成 绩 评 定:

指 导 教 师 (签 名):

年 月 日

这个人太懒什么东西都没留下

文章评论(0)