磁性體磁場(chǎng)正演

磁法勘探上機(jī)實(shí)驗(yàn)名：學(xué)號(hào)：指導(dǎo)教師：日期：2020.4.15實(shí)驗(yàn)二：磁性體磁場(chǎng)正演一、實(shí)驗(yàn)?zāi)康模?、通過(guò)球體、水平圓柱體磁場(chǎng)的正演計(jì)算，掌握簡(jiǎn)單規(guī)則磁性體正演磁場(chǎng)的計(jì)算方法；2、通過(guò)計(jì)算認(rèn)識(shí)球體與水平圓柱體磁場(chǎng)的一般分布規(guī)律，了解影響磁性體磁場(chǎng)的主要因素（如磁性體的形體、物性參數(shù)、走向或計(jì)算剖面的選擇等），培養(yǎng)學(xué)生實(shí)際動(dòng)手能力與分析問(wèn)題的能力。二、實(shí)驗(yàn)內(nèi)容用Matlab語(yǔ)言或C語(yǔ)言編程實(shí)現(xiàn)球體和水平圓柱體的磁場(chǎng)（包括Za、Ha、At）的正演計(jì)算。三、實(shí)驗(yàn)要求假設(shè)地磁場(chǎng)方向與磁性體磁化強(qiáng)度方向一致且均勻磁化的情況下,當(dāng)?shù)卮艌?chǎng)T=50000nT，磁傾角1=60°,球體與水平圓柱體中心埋深R=30m,半徑r=10m,磁化率k=0.1（SI），計(jì)算（觀測(cè)）剖面磁化強(qiáng)度水平投影夾角A=0°時(shí)：1、 正演計(jì)算球體的磁場(chǎng)（Za、Hax、Hay、AT）,畫(huà)出對(duì)應(yīng)的平面等值線圖、曲面圖及主剖面異常圖；2、 正演計(jì)算水平圓柱體的磁場(chǎng)（Za、Ha、AT）,畫(huà)出主剖面異常結(jié)果圖；3、 通過(guò)改變球體或水平圓柱體的幾何參數(shù)、磁化強(qiáng)度方向（I）、計(jì)算剖面的方位角（Az）,觀察主剖面磁場(chǎng)Za的變化，分析磁化方向與計(jì)算剖面對(duì)磁性體磁場(chǎng)特征的影響。四、實(shí)驗(yàn)原理球體與水平圓柱體磁場(chǎng)（Za、Ha、AT）的計(jì)算公式是以磁化強(qiáng)度傾角I、有效磁化傾角is和剖面與磁化強(qiáng)度水平投影夾角A'來(lái)表達(dá)。1、球體磁場(chǎng)的正演公式：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、水平圓柱體磁場(chǎng)的正演公式：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、有效磁化強(qiáng)度Ms與有效磁化傾角is：M=(M2+M2)1/2=M(cos21cos2A+sin21)s Mx z\=tg-1M=tg-1(tglsecA')五、實(shí)驗(yàn)報(bào)告（內(nèi)容包括實(shí)驗(yàn)?zāi)康?、?shí)驗(yàn)內(nèi)容、實(shí)驗(yàn)原理、計(jì)算程序代碼、實(shí)驗(yàn)結(jié)果、結(jié)果分析或小結(jié)）計(jì)算程序代碼：球體：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