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

实验2--磁性体磁场正演程序?应用地磁学?实验报告姓名：张嘉琪学号：1010112225指导教师：李淑玲实验地点：实验室319实验日期：2021-05-24实验二：磁性体磁场正演一、实验目的：1、通过球体、水平圆柱体磁场的正演计算，掌握简单规那么磁性体正演磁场的计算；2、通过计算认识球体与水平圆柱体磁场的一般分布规律，了解影响磁性体磁场的主要因素〔如磁性体的形体、物性参数、走向或计算剖面的选择等〕，培养学生实际动手能力与问题的能力。二、实验内容用Matlab语言或C语言编程实现球体和水平圆柱体的磁场(包括Za、Ha、Δt)的正演计算。三、实验要求假设地磁场方向与磁性体磁化强度方向一致且均匀磁化的情况下,当地磁场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和剖面与磁化强度水平投影夹角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);%转化为排列%球体参数i=pi/3;%磁化倾角ia=0;%剖面磁方位角R=10;%球体半径mv=4/3*pi*R^3u=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)*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);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Hax异常');axisequal,axis([-5050-5050]),colorbar;subplot(222),contourf(X,Y,Hay);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Hay异常');axisequal,axis([-5050-5050]),colorbar;subplot(223),contourf(X,Y,Za);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体Za异常');axisequal,axis([-5050-5050]),colorbar;subplot(224),contourf(X,Y,T);xlabel('X(m)'),ylabel('Y(m)'),title('理论球体ΔT异常');axisequal,axis([-5050-5050]),colorbar;%绘制曲面图〔三维〕figure(2),clf,subplot(221),mesh(X,Y,Hax),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Hax异常'),colorbar;subplot(222),mesh(X,Y,Hay),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Hay异常'),colorbar;subplot(223),mesh(X,Y,Za),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Za异常'),colorbar;subplot(224),surf(X,Y,T),shadinginterp,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)))./(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(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,fori=0:pi/6:pi/2Za2=(u*m*((2*D.^2-x.^2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.^2+D.^2).^(5/2));holdonplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁倾角改变)'),gridon;endh=legend('Za');legend(h,'boxoff');figure(5),clf,fora=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));holdonplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁方位改变)'),gridon;endh=legend('Za');legend(h,'boxoff');figure(6),clf,fori=pi/3;a=0;R=10:5:20v=4/3*pi*R^3m=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));holdonplot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(球体半径)'),gridon;endh=legend('Za');legend(h,'boxoff');圆柱体程序代码：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);%转化为排列%水平圆柱体参数i=pi/3;%磁化倾角a=0;%剖面磁方位角Is=(tan(tan(i)*sec(a)))^(-1);R=10;%圆柱体横截面半径mS=pi*R^2;%圆柱体横截面面积u=4*pi*10^(-7);%磁导率T=0.5*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*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异常');axisequal,axis([-200200-200200]),colorbar;subplot(222),contourf(X,Y,Ha);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体Ha异常');axisequal,axis([-200200-200200]),colorbar;subplot(223),contourf(X,Y,T);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体ΔT异常');axisequal,axis([-200200-200200]),colorbar;%绘制曲面图〔三维〕figure(2),clf,%clf去除图形subplot(221),mesh(X,Y,Za),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体Za异常'),colorbar;subplot(222),mesh(X,Y,Ha),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体Ha异常'),colorbar;subplot(223),surf(X,Y,T),shadinginterp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论圆柱体ΔT异常'),colorbar;%主剖面视图figure(3),clf,subplot(311)forx=-200:5:200Za1=(u*m*((D.^2-x.^2)*sin(Is)-2*D*x.*cos(Is)))./(2*pi*(x.^2+D.^2)^2);holdon;plot(x,Za1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Za异常');endsubplot(312)forx=-200:5:200Ha1=(-u*m*((D^2-x^2)*cos(Is)+2*D*x*sin(Is)))/(2*pi*(x^2+D^2)^2);holdon;plot(x,Ha1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Ha异常');endsubplot(313)forx=-200:5:200T1=(u*m*sin(i)*((D^2-x^2)*cos(2*i-pi)-2*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)));holdon;plot(x,T1,'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体T异常');end%绘制异常剖面图figure(4),clf,fori=0:pi/6:pi/2Is=(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);holdonplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常(磁倾角)'),gridon;endh=legend('Za');legend(h,'boxoff');figure(5),clf,fora=0:pi/6:pi/2i=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);holdonplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常〔方位角〕'),gridon;endh=legend('Za');legend(h,'boxoff');figure(6),clf,forR=10:5:20i=pi/3;a=0;S=pi*R^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);holdonplot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常〔半径〕'),gridon;endh=legend('Za');legend(h,'boxoff');实验结果：球体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：〔半径改变〕磁倾角改变后的主剖面视图：剖面方位角改变后的剖面图：圆柱体实验结果：平面等值线图：曲面图：主剖面视图：球体参数改变后的主剖面视图：〔半径改变〕磁倾角改变后的主剖面视图：剖面方位角改变的异常图：七、结果分析：球体分析：平面图特征：球体的磁场ΔT不仅与地磁场方向有关，还与观测剖面有关。由球体设置参数及观测剖面可知，球体相当于斜磁化，Za和ΔT不相等，剖面异常曲线不对称，平面异常为正负伴生的近等轴状异常，其中ΔT受斜磁化影响更大，在相同的磁化倾角其负值较大，异常极大值和极小值的连线与磁化强度矢量的水平投影方向一致。剖面特征：球体沿特定方向磁化，该方向在地面投影即为主剖面方位，主剖面视图中的Hay为一直线，不随x变化，斜磁化时，主剖面Za为两边有负值的非对称曲线，当磁倾角由pi/2至0变化时，Za极大值减小，极大值开始增大；当磁方位角由pi/2至0变化时，当球体半径增大时，Za异常曲线形态不变，异常幅值明显变大。圆柱体分析：平面图特征：Ha异常为近等轴状，中间出现极大值点，Za和ΔT不相等，ΔT受斜磁化影响比Za更大〔正值更小，负值更大〕，异常极大值和极小值的连线与磁化强度矢量的水平投影方向一致。剖面特征：磁倾角为pi/2时，圆柱体相当于垂直磁化，异常关于原点对称，圆柱体半径增大时，异常曲线形态根本保持不变，幅值变大。八、实验小结：通过此次实验熟悉了解了利用matlab对根本的规那么球体与水平圆柱体的平面等值线与剖面图的绘制，并了解了磁异常的根本分布规律；，了解影响磁性体磁场的主要因素〔如磁性体的形体、物性参数、走向或计算剖面的选择等〕，同时培养了实际动手能力与分析问题的能力。