实验2磁性体磁场正演程序

1、应用地磁学实验报告 姓 名： 张嘉琪 学 号： 1010112225 指导教师： 李淑玲 实验地点： 实验室319 实验日期： 2014-05-24实验二：磁性体磁场正演 一、实验目的：1、通过球体、水平圆柱体磁场的正演计算，掌握简单规则磁性体正演磁场的计算方法；2、通过计算认识球体与水平圆柱体磁场的一般分布规律，了解影响磁性体磁场的主要因素（如磁性体的形体、物性参数、走向或计算剖面的选择等），培养学生实际动手能力与分析问题的能力。二、实验内容用matlab语言或c语言编程实现球体和水平圆柱体的磁场(包括za、ha、t)的正演计算。 三、实验要求假设地磁场方向与磁性体磁化强度方向一致且均匀磁化

2、的情况下,当地磁场t=50000nt，磁倾角i=60°，球体与水平圆柱体中心埋深r=30m，半径r=10m，磁化率k=0.2（si），计算（观测）剖面磁化强度水平投影夹角a=0°时：1、正演计算球体的磁场（za、hax、hay、t），画出对应的平面等值线图、曲面图及主剖面异常图；2、正演计算水平圆柱体的磁场（za、ha、t），画出主剖面异常结果图；3、通过改变球体与水平圆柱体的几何参数、磁化强度方向（i）、计算剖面的方位角(a)，观察主剖面磁场za的变化，分析磁化方向与计算剖面对磁性体磁场特征的影响。 四、实验原理球体与水平圆柱体磁场（za、ha、t）的计算公式是以磁化强度

3、倾角i、有效磁化倾角is和剖面与磁化强度水平投影夹角a来表达。1、球体磁场的正演公式： 2、水平圆柱体磁场的正演公式：3、有效磁化强度ms与有效磁化倾角is：五、计算程序代码：1、球体matlab代码：clc;clear;% 测点分布范围dx=5; % x方向测点间距dy=5; % y方向测点间距nx=81; % x方向测点数ny=81; % y方向测点数xmin=-200; % x方向起点ymin=-200; % y方向起点x=xmin:dx:(xmin+(nx-1)*dx); % x方向范围y=ymin:dy:(ymin+(ny-1)*dy); % y方向范围x,y=meshgrid(x,

4、y); % 转化为排列 % 球体参数i=pi/3; %磁化倾角ia=0; %剖面磁方位角r=10; % 球体半径 mv=4/3*pi*r3u=4*pi*10(-7);%磁导率t=0.5*10(-4);%地磁场强度k=0.2;%磁化率m=k*t/u; %磁化强度 a/mm=m*v; %磁矩d=30; % 球体埋深 m% 球体za理论磁异常za=(u*m*(2*d.2-x.2-y.2)*sin(i)-3*d*x.*cos(i)*cos(a)-3*d*y.*cos(i)*sin(a)./(4*pi*(x.2+y.2+d.2).(5/2);% 球体hax理论磁异常hax=(u*m*(2*x.2-y.2

5、-d.2)*cos(i)*cos(a)-3*d*x.*sin(i)+3*x.*y.*cos(i)*sin(a)./(4*pi*(x.2+y.2+d.2).(5/2);%球体hay理论磁异常hay=(u*m*(2*y.2-x.2-d.2)*cos(i)*sin(a)-3*d*y.*sin(i)+3*x.*y.*cos(i)*cos(a)./(4*pi*(x.2+y.2+d.2).(5/2);%球体t理论异常t=hax*cos(i)*cos(a)+hay*cos(i)*sin(a)+za*sin(i);%绘平面异常等值线图(二维)figure(1),clf,subplot(221),contour

6、f(x,y,hax);xlabel('x(m)'),ylabel('y(m)'),title('理论球体hax异常');axis equal,axis(-50 50 -50 50),colorbar;subplot(222),contourf(x,y,hay);xlabel('x(m)'),ylabel('y(m)'),title('理论球体hay异常');axis equal,axis(-50 50 -50 50),colorbar;subplot(223),contourf(x,y,za);xl

7、abel('x(m)'),ylabel('y(m)'),title('理论球体za异常');axis equal,axis(-50 50 -50 50),colorbar;subplot(224),contourf(x,y,t);xlabel('x(m)'),ylabel('y(m)'),title('理论球体t异常');axis equal,axis(-50 50 -50 50),colorbar;%绘制曲面图（三维）figure(2),clf,subplot(221),mesh(x,y,hax)

8、,shading interp,xlabel('x(m)'),ylabel('y(m)'),zlabel('理论球体hax异常'),colorbar;subplot(222),mesh(x,y,hay),shading interp,xlabel('x(m)'),ylabel('y(m)'),zlabel('理论球体hay异常'),colorbar;subplot(223),mesh(x,y,za),shading interp,xlabel('x(m)'),ylabel('

9、y(m)'),zlabel('理论球体za异常'),colorbar;subplot(224),surf(x,y,t),shading interp,xlabel('x(m)'),ylabel('y(m)'),zlabel('理论球体t异常'),colorbar;%绘制主剖面异常等值线za1=(u*m*(2*d.2-x.2)*sin(i)-3*d*x.*cos(i)*cos(a)./(4*pi*(x.2+d.2).(5/2);hax1=(u*m*(2*x.2-d.2)*cos(i)*cos(a)-3*d*x.*sin(i)

10、./(4*pi*(x.2+d.2).(5/2);hay1=(u*m*(-x.2-d.2)*cos(i)*sin(a)./(4*pi*(x.2+d.2).(5/2);t1=hax1*cos(i)*cos(a)+hay1*cos(i)*sin(a)+za1*sin(i);figure(3),clf,subplot(221)plot(x,za1,'g-','linewidth',1.3);xlabel('x(m)'),ylabel('理论球体za异常');subplot(222)plot(x,hax1,'k-','

11、;linewidth',1.3);xlabel('x(m)'),ylabel('理论球体hax异常');subplot(223)plot(x,hay1,'r-','linewidth',1.3);xlabel('x(m)'),ylabel('理论球体hay异常');subplot(224)plot(x,t1,'b-','linewidth',1.3);xlabel('x(m)'),ylabel('理论球体t异常'); %绘制异常

12、剖面图figure(4),clf,for i=0:pi/6:pi/2 za2=(u*m*(2*d.2-x.2)*sin(i)-3*d*x.*cos(i)*cos(a)./(4*pi*(x.2+d.2).(5/2); hold onplot(x,za2,'r-','linewidth',1.3),xlabel('x(m)'),ylabel('磁力异常(磁倾角改变)'),grid on;endh=legend('za');legend(h,'boxoff'); figure(5),clf,for a=0

13、:pi/6:pia=pi/3; za2=(u*m*(2*d.2-x.2)*sin(i)-3*d*x.*cos(i)*cos(a)./(4*pi*(x.2+d.2).(5/2); hold onplot(x,za2,'r-','linewidth',1.3),xlabel('x(m)'),ylabel('磁力异常(磁方位改变)'),grid on;endh=legend('za');legend(h,'boxoff'); figure(6),clf,for i=pi/3;a=0; r=10:5:20

14、v=4/3*pi*r3 m=m*v; za2=(u*m*(2*d.2-x.2)*sin(i)-3*d*x.*cos(i)*cos(a)./(4*pi*(x.2+d.2).(5/2); hold onplot(x,za2,'r-','linewidth',1.3),xlabel('x(m)'),ylabel('磁力异常(球体半径)'),grid on;endh=legend('za');legend(h,'boxoff');圆柱体程序代码：clc;clear;% 测点分布范围dx=5; % x方向测点

15、间距dy=5; % y方向测点间距nx=81; % x方向测点数ny=81; % y方向测点数xmin=-200; % x方向起点ymin=-200; % y方向起点x=xmin:dx:(xmin+(nx-1)*dx); % x方向范围y=ymin:dy:(ymin+(ny-1)*dy); % y方向范围x,y=meshgrid(x,y); % 转化为排列 % 水平圆柱体参数i=pi/3; %磁化倾角a=0;%剖面磁方位角is=(tan(tan(i)*sec(a)(-1);r=10; % 圆柱体横截面半径 ms=pi*r2; %圆柱体横截面面积u=4*pi*10(-7); %磁导率t=0.5*

16、10(-4);%地磁场强度k=0.2;%磁化率m=k*t/u; %磁化强度 a/mms=m*(cos(i)*cos(a)2+(sin(i)2);m=ms*s; %单位长度的有效磁矩d=30; % 圆柱体中心点埋深 m % 圆柱体za理论磁异常za=(u*m*(d.2-x.2)*sin(is)-2*d*x.*cos(is)./(2*pi*(x.2+d.2)2);% 圆柱体ha理论磁异常ha=(-u*m*(d.2-x.2)*cos(is)+2*d*x.*sin(is)./(2*pi*(x.2+d.2)2);%圆柱体t理论异常t=(u*m*sin(i)*(d.2-x.2)*cos(2*i-pi)-2

17、*d*x.*sin(2*is-pi/2)./(sin(is)*(d.2-x.2)*sin(2*is-pi/2)-2*d*x.*cos(2*is-pi/2);%绘平面异常等值线图(二维)figure(1),clf,subplot(221),contourf(x,y,za);xlabel('x(m)'),ylabel('y(m)'),title('理论圆柱体za异常');axis equal,axis(-200 200 -200 200),colorbar;subplot(222),contourf(x,y,ha);xlabel('x(m)&

18、#39;),ylabel('y(m)'),title('理论圆柱体ha异常');axis equal,axis(-200 200 -200 200),colorbar;subplot(223),contourf(x,y,t);xlabel('x(m)'),ylabel('y(m)'),title('理论圆柱体t异常');axis equal,axis(-200 200 -200 200),colorbar;%绘制曲面图（三维）figure(2),clf,%clf清除图形subplot(221),mesh(x,y,z

19、a),shading interp,xlabel('x(m)'),ylabel('y(m)'),zlabel('理论圆柱体za异常'),colorbar;subplot(222),mesh(x,y,ha),shading interp,xlabel('x(m)'),ylabel('y(m)'),zlabel('理论圆柱体ha异常'),colorbar;subplot(223),surf(x,y,t),shading interp,xlabel('x(m)'),ylabel('

20、y(m)'),zlabel('理论圆柱体t异常'),colorbar;%主剖面视图figure(3),clf,subplot(311) for x=-200:5:200 za1=(u*m*(d.2-x.2)*sin(is)-2*d*x.*cos(is)./(2*pi*(x.2+d.2)2); hold on; plot(x,za1,'b-*','linewidth',1.3),xlabel('x(m)'),ylabel('圆柱体za异常');endsubplot(312)for x=-200:5:200 h

21、a1=(-u*m*(d2-x2)*cos(is)+2*d*x*sin(is)/(2*pi*(x2+d2)2); hold on; plot(x,ha1,'b-*','linewidth',1.3),xlabel('x(m)'),ylabel('圆柱体ha异常');endsubplot(313)for x=-200:5:200 t1=(u*m*sin(i)*(d2-x2)*cos(2*i-pi)-2*d*x*sin(2*is-pi/2)/(sin(is)*(d2-x2)*sin(2*is-pi/2)-2*d*x*cos(2*is-p

22、i/2); hold on; plot(x,t1,'b-*','linewidth',1.3),xlabel('x(m)'),ylabel('圆柱体t异常');end%绘制异常剖面图figure(4),clf,for i=0:pi/6:pi/2 is=(tan(tan(i)*sec(a)(-1); ms=m*(cos(i)*cos(a)2+(sin(i)2); m=ms*s; za1=(u*m*(d.2-x.2)*sin(is)-2*d*x.*cos(is)./(2*pi*(x.2+d.2)2); hold on plot(x,z

23、a1,'r-*'),xlabel('y(m)'),ylabel('磁力异常(磁倾角)'),grid on;endh=legend('za');legend(h,'boxoff'); figure(5),clf,for a=0:pi/6:pi/2 i=pi/3; is=(tan(tan(i)*sec(a)(-1); ms=m*(cos(i)*cos(a)2+(sin(i)2); m=ms*s; za1=(u*m*(d.2-x.2)*sin(is)-2*d*x.*cos(is)./(2*pi*(x.2+d.2)2);

24、hold on plot(x,za1,'r-*'),xlabel('y(m)'),ylabel('磁力异常（方位角）'),grid on;endh=legend('za');legend(h,'boxoff'); figure(6),clf,for r=10:5:20 i=pi/3; a=0; s=pi*r2; m=ms*s; za1=(u*m*(d.2-x.2)*sin(is)-2*d*x.*cos(is)./(2*pi*(x.2+d.2)2); hold on plot(x,za1,'r-*'),xlabel('y(m)'),ylabel('磁力异常（半径）'),