磁性体磁场正演

磁法勘探上机实验名：学号：指导教师：日期：2020.4.15实验二：磁性体磁场正演一、实验目的：1、通过球体、水平圆柱体磁场的正演计算，掌握简单规则磁性体正演磁场的计算方法；2、通过计算认识球体与水平圆柱体磁场的一般分布规律，了解影响磁性体磁场的主要因素（如磁性体的形体、物性参数、走向或计算剖面的选择等），培养学生实际动手能力与分析问题的能力。二、实验内容用Matlab语言或C语言编程实现球体和水平圆柱体的磁场（包括Za、Ha、At）的正演计算。三、实验要求假设地磁场方向与磁性体磁化强度方向一致且均匀磁化的情况下,当地磁场T=50000nT，磁倾角1=60°,球体与水平圆柱体中心埋深R=30m,半径r=10m,磁化率k=0.1（SI），计算（观测）剖面磁化强度水平投影夹角A=0°时：1、 正演计算球体的磁场（Za、Hax、Hay、AT）,画出对应的平面等值线图、曲面图及主剖面异常图；2、 正演计算水平圆柱体的磁场（Za、Ha、AT）,画出主剖面异常结果图；3、 通过改变球体或水平圆柱体的几何参数、磁化强度方向（I）、计算剖面的方位角（Az）,观察主剖面磁场Za的变化，分析磁化方向与计算剖面对磁性体磁场特征的影响。四、实验原理球体与水平圆柱体磁场（Za、Ha、AT）的计算公式是以磁化强度倾角I、有效磁化倾角is和剖面与磁化强度水平投影夹角A'来表达。1、球体磁场的正演公式：Li m \ fH= •L（2x2-y2-R2）cosIcosAax4兀（x2+y2+R2）5/2-3RxsinI+3xycosIsinA]imH=电 L（2y2-x2-R2）cosIsinAay4兀（x2+y2+R2）5/2 >-3RysinI+3xycosIcosA’]Z=电 [（2R2-x2-y2）sinIa4兀（x2+y2+R2）5/2-3RxcosIcosA‘-3RycosIsinA’]AT=Z m —[(2R2-x2-y2)sin21+(2x2-y2-R2)cos21cos2A4兀 2[2 2^5/2x2+y2+R2+(2y2一x2一R2)cos21sin2A-3xRsin21cosA+3xycos21sin2A一3yRsin21sinA]2、水平圆柱体磁场的正演公式：TOC\o"1-5"\h\zhm 1Z= _s [(R2一x2)sini一2Rxcosi]a 2兀(x2+R2)2 s shm1H=一 ―s [(R2一x2)cosi+2Rxsini]a 2兀(x2+R2)2 s s)Z2-x2)in(i―90)2Rxcos(i―90。2sini s ss3、有效磁化强度Ms与有效磁化倾角is：M=(M2+M2)1/2=M(cos21cos2A+sin21)s Mx z\=tg-1M=tg-1(tglsecA')五、实验报告（内容包括实验目的、实验内容、实验原理、计算程序代码、实验结果、结果分析或小结）计算程序代码：球体：dx=5dy=5nx=81ny=81xmin=200ymin=200x=xmin:dx:(xmin+(nx1)*dx)y=ymin:dy:(ymin+(ny1)*dy)[X,Y]=meshgrid(x,y)T=5*10A(5)a=0i=pi/3R1=10v1=4/3*pi*R"3u=4*pi*10P)k=0.1M=k*T/um=M*v1D=30Za=(u*m*((2*D.A2-X.A2-Y.A2)*sin(i)-3*D*X.*cos(i)*cos(a)-3*D*Y.*cos(i)*sin(a)))./(4*pi*(X「2+Y.A2+D.A2)「(5/2));Hax=(u*m*((2*X.A2-Y.A2-D.A2)*cos(i)*cos(a)-3*D*X.*sin(i)+3*D*X.*cos(i)*sin(a)))./(4*pi*(X.A2+Y.A2+D.A2).A(5/2));Hay=(u*m*((2*X.A2-Y.A2-D.A2)*cos(i)*sin(a)-3*D*X.*sin(i)+3*D*X.*cos(i)*cos(a)))./(4*pi*(X.A2+Y.A2+D.A2).A(5/2));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('理论球体At异常');axisequal,axis([-5050-5050]),colorbar;figure(2),clf,subplot(221),mesh(X,Y,Hax),xlabel('X(m1)'),ylabel('Y(m1)'),zlabel('球体Hax异常'),colorbar;subplot(222),mesh(X,Y,Hay),xlabel('X(m1)'),ylabel('Y(m1)'),zlabel('球体Hay异常'),colorbar;subplot(223),mesh(X,Y,Za),xlabel('X(m1)'),ylabel('Y(m1)'),zlabel('球体Za异常'),colorbar;subplot(224),mesh(X,Y,T),xlabel('X(m1)'),ylabel('Y(m1)'),zlabel('球体△t异常'),colorbar;Za1=(u*m*((2*D.A2-x.A2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.A2+D.八2)"(5/2));Hax1=(u*m*((2*x.A2-D.A2)*cos(i)*cos(a)-3*D*x.*sin(i)))./(4*pi*(x.A2+D.八2).八(5/2));Hay1=(u*m*((-x.A2-D.A2)*cos(i)*sin(a)))./(4*pi*(x.A2+D.A2).A(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('理论球体At异常');figure(4),clf;fori=0:pi/6:pi/2;Za2=(u*m*((2*D.A2-x.A2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.A2+D.八2).八(5/2));holdon;plot(x,Za2'b-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常:磁倾角改变'),gridon;endh=legend('Za');legend(h,'boxoff');figure(5),clf;fora=0:pi/6:pi;i=pi/3;Za2=(u*m*((2*D.A2-x.A2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.A2+D.八2).八(5/2));holdon;plot(x,Za2'g-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常：磁方位角改变'),gridon;endh=legend('Za');legend(h,'boxoff');figure(6),clf;forR1=10:5:20;v1=4/3*pi*R"3;m=M*v1;i=pi/3;a=0;Za2=(u*m*((2*D.A2-x.A2)*sin(i)-3*D*x.*cos(i)*cos(a)))./(4*pi*(x.A2+D.八2).八(5/2));holdon;plot(x,Za2'y-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常:球体半径改变'),gridon;endh=legend('Za');legend(h,'boxoff');El.2曲* -Hni zjlt：・m-Kiwt^r-r□ k 営疫棗-Or口自■口Knd吐sam ZiVi XBm*[>. MnkJdWk'、_、:SX-1□□•□«E谡诉aD«Aq:1l圆柱：dx=5dy=5nx=81ny=81xmin=200ymin=200x=xmin:dx:(xmin+(nx1)*dx)y=ymin:dy:(ymin+(ny1)*dy)[X,Y]=meshgrid(x,y)i=pi/3;a=0;Is=(tan(tan(i)*sec(a)))八(-1);R=10;S=pi*RA2;u=4*pi*10入(_7);T=0.5*10A(-4);k=0.2;M=k*T/u;Ms=((cos(i)*cos(a))A2+(sin(i))A2);m=Ms*S;D=30;Za=(u*m*((D.A2_X.A2)*sin(Is)_2*D*X.*cos(Is)))./(2*pi*(X.A2+D.A2)A2);Ha=(_u*m*((D.A2_X.A2)*cos(Is)+2*D*X.*sin(Is)))./(2*pi*(X.A2+D.A2)A2);T=(u*m*sin(i)*((D.A2_X.A2)*cos(2*i_pi)_2*D*X.*sin(2*Is_pi/2)))./(sin(Is)*((D.A2_X.A2)*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('理论圆柱体Hax异常');axisequal,axis([_200200_200200]),colorbar;subplot(223);contourf(X,Y,T);xlabel('X(m)'),ylabel('Y(m)'),title('理论圆柱体At异常');axisequal,axis([_200200_200200]),colorbar;figure(2),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('理论圆柱体At异常'),colorbar;figure(3),clf;subplot(311);forx=_200:5:200;Za1=(u*m*((D.A2_x.A2)*sin(Is)_2*D*x.*cos(Is)))./(2*pi*(x.A2+D.A2)A2);holdon;plot(x,Za1'b-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Za异常');endsubplot(312);forx=-200:5:200;Ha1=(-u*m*((DA2-xA2)*cos(Is)+2*D*x*sin(Is)))/(2*pi*(xA2+DA2)A2);holdon;plot(x,Ha1'y-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体Ha异常');endsubplot(313);forx=-200:5:200;Tl=(u*m*sin(i)*((DA2-xA2)*cos(2*i-pi)-2*D*x*sin(2*Is-pi/2)))./(sin(Is)*((DA2-xA2)*sin(2*Is-pi/2)-2*D*x*cos(2*Is-pi/2)));holdon;plot(x,T1,'g-*','linewidth',1.3),xlabel('X(m)'),ylabel('圆柱体T异常');endfigure(4),clf;fori=0:pi/6:pi/2;Is=(tan(tan(i)*sec(a)))A(-1);Ms=M*((cos(i)*cos(a)八2+sin(i)))八2;m=Ms*S;Za1=(u*m*((D.A2-X.A2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.A2+D.A2)A2);holdon;plot(X,Za1'r-*'),xlabel('Y(m)'),ylabel('磁力异常：磁倾角'),gridon;endh=legend('Za1');legend(h,'boxoff');figure(5),clf;fora=0：pi/6：pi/2;i=pi/3;Is=(tan(tan(i)*sec(a)))A(-1);m=Ms*S;Za1=(u*m*((D.A2-X.A2)*sin(Is)-2*D*X.*cos(Is)))./(2*pi*(X.A2+D.A2)A2);holdon;plot(X,Za1'y-*'),xlabel('Y(m)'),ylabel('磁力异常:方位角'),gridon;endh=legend('Za');legend(h,'boxoff');figure(6),clf;forR=10:5:20;—UB-K4Hst¥1tfs^sP9aFIT_h<w.-i~亠(-uoxoq\q)pu①6①H亠(-f