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

《实验2磁性体磁场正演程序.doc》由会员分享，可在线阅读，更多相关《实验2磁性体磁场正演程序.doc（15页珍藏版）》请在三一办公上搜索。

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）的计算公式是以磁化强度倾角I、有效磁化倾角is

3、和剖面与磁化强度水平投影夹角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,y); % 转化为排列

4、% 球体参数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-D.2)*cos(i)

5、*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),contourf(X,Y,Hax);x

6、label(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);xlabel(X(m),ylabel(Y(m),title(理论球体Za异常);axis equal,axis(-50 50 -50 50),colorba

7、r;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),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(

8、理论球体Hay异常),colorbar;subplot(223),mesh(X,Y,Za),shading interp,xlabel(X(m),ylabel(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-

9、D.2)*cos(i)*cos(a)-3*D*x.*sin(i)./(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-,linewidth,1.3);xlabel(X

10、(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异常); %绘制异常剖面图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-,lin

11、ewidth,1.3),xlabel(X(m),ylabel(磁力异常(磁倾角改变),grid on;endh=legend(Za);legend(h,boxoff); figure(5),clf,for a=0: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,boxof

12、f); figure(6),clf,for i=pi/3;a=0; R=10:5:20 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方向测点间距dy=5; % Y方向测点间距n

13、x=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*10(-4);%地磁场强度k=0.2

14、;%磁化率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*D*X.*sin(2*Is-pi/

15、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),ylabel(Y(m),title(理论圆柱体Ha异常);axis equal,axis(-200 200 -

16、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,Za),shading interp,xlabel(X(m),ylabel(Y(m),zlabel(理论圆柱体Za异常),colorbar;subplot(222),mesh(X,Y,Ha),shading interp,xlab

17、el(X(m),ylabel(Y(m),zlabel(理论圆柱体Ha异常),colorbar;subplot(223),surf(X,Y,T),shading interp,xlabel(X(m),ylabel(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),

18、ylabel(圆柱体Za异常);endsubplot(312)for x=-200:5:200 Ha1=(-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*co

19、s(2*Is-pi/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,Za1,r-*),xlabel(Y(m),ylabel(磁力异常(磁

20、倾角),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); hold on plot(X,Za1,r-*),xlabel(Y(m),ylabel(磁力异常（方位角）),grid on;endh=legend(Za);legend(

21、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(磁力异常（半径）),grid on;endh=legend(Za);legend(h,boxoff);六、 实验结果：球体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：（半径改变）磁倾角改变后的主剖面视图：剖面方位角改变后的剖

22、面图：圆柱体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：（半径改变）磁倾角改变后的主剖面视图：剖面方位角改变的异常图：七、结果分析：球体分析：平面图特征：球体的磁场T不仅与地磁场方向有关，还与观测剖面有关。由球体设置参数及观测剖面可知，球体相当于斜磁化，Za和T不相等，剖面异常曲线不对称，平面异常为正负伴生的近等轴状异常，其中T受斜磁化影响更大，在相同的磁化倾角其负值较大，异常极大值和极小值的连线与磁化强度矢量的水平投影方向一致。剖面特征：球体沿特定方向磁化，该方向在地面投影即为主剖面方位，主剖面视图中的Hay为一直线，不随x变化，斜磁化时，主剖面Za为两边有负值的

23、非对称曲线，当磁倾角由pi/2至0变化时，Za极大值减小，极大值开始增大；当磁方位角由pi/2至0变化时，当球体半径增大时，Za异常曲线形态不变，异常幅值明显变大。圆柱体分析： 平面图特征：Ha异常为近等轴状，中间出现极大值点，Za和T不相等，T受斜磁化影响比Za更大（正值更小，负值更大），异常极大值和极小值的连线与磁化强度矢量的水平投影方向一致。剖面特征：磁倾角为pi/2时，圆柱体相当于垂直磁化，异常关于原点对称，圆柱体半径增大时，异常曲线形态基本保持不变，幅值变大。八、实验小结： 通过此次实验熟悉了解了利用matlab对基本的规则球体与水平圆柱体的平面等值线与剖面图的绘制，并了解了磁异常的基本分布规律；，了解影响磁性体磁场的主要因素（如磁性体的形体、物性参数、走向或计算剖面的选择等），同时培养了实际动手能力与分析问题的能力。