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

这是我的一份实验报告 

本科生实验报告

填写说明

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 7714－2005）》。
      

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

1. 理论公式
   

模型图

(x_s,y_s,z_s)为长方体上任意点，(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方向的两个边界，正异常为下边界，负异常为上边界。垂向导数也可以判断地质体的边界，也可以描述地质体的形态。