實(shí)驗(yàn)2磁性體磁場(chǎng)正演程序

1、應(yīng)用地磁學(xué)課程實(shí)驗(yàn)報(bào)告 應(yīng)用地磁學(xué)實(shí)驗(yàn)報(bào)告 姓 名： 張嘉琪 學(xué) 號(hào)： 1010112225 指導(dǎo)教師： 李淑玲 實(shí)驗(yàn)地點(diǎn)： 實(shí)驗(yàn)室319 實(shí)驗(yàn)日期： 2014-05-24第15頁 實(shí)驗(yàn)二：磁性體磁場(chǎng)正演 一、實(shí)驗(yàn)?zāi)康模?、通過球體、水平圓柱體磁場(chǎng)的正演計(jì)算，掌握簡(jiǎn)單規(guī)則磁性體正演磁場(chǎng)的計(jì)算方法；2、通過計(jì)算認(rèn)識(shí)球體與水平圓柱體磁場(chǎng)的一般分布規(guī)律，了解影響磁性體磁場(chǎng)的主要因素（如磁性體的形體、物性參數(shù)、走向或計(jì)算剖面的選擇等），培養(yǎng)學(xué)生實(shí)際動(dòng)手能力與分析問題的能力。二、實(shí)驗(yàn)內(nèi)容用matlab語言或c語言編程實(shí)現(xiàn)球體和水平圓柱體的磁場(chǎng)(包括za、ha、t)的正演計(jì)算。 三、實(shí)驗(yàn)要求假設(shè)地磁場(chǎng)方向

2、與磁性體磁化強(qiáng)度方向一致且均勻磁化的情況下,當(dāng)?shù)卮艌?chǎng)t=50000nt，磁傾角i=60，球體與水平圓柱體中心埋深r=30m，半徑r=10m，磁化率k=0.2（si），計(jì)算（觀測(cè)）剖面磁化強(qiáng)度水平投影夾角a=0時(shí)：1、正演計(jì)算球體的磁場(chǎng)（za、hax、hay、t），畫出對(duì)應(yīng)的平面等值線圖、曲面圖及主剖面異常圖；2、正演計(jì)算水平圓柱體的磁場(chǎng)（za、ha、t），畫出主剖面異常結(jié)果圖；3、通過改變球體與水平圓柱體的幾何參數(shù)、磁化強(qiáng)度方向（i）、計(jì)算剖面的方位角(a)，觀察主剖面磁場(chǎng)za的變化，分析磁化方向與計(jì)算剖面對(duì)磁性體磁場(chǎng)特征的影響。 四、實(shí)驗(yàn)原理球體與水平圓柱體磁場(chǎng)（za、ha、t）的計(jì)算公式是

3、以磁化強(qiáng)度傾角i、有效磁化傾角is和剖面與磁化強(qiáng)度水平投影夾角a來表達(dá)。1、球體磁場(chǎng)的正演公式： 2、水平圓柱體磁場(chǎng)的正演公式：3、有效磁化強(qiáng)度ms與有效磁化傾角is：五、計(jì)算程序代碼：1、球體matlab代碼：clc;clear;% 測(cè)點(diǎn)分布范圍dx=5; % x方向測(cè)點(diǎn)間距dy=5; % y方向測(cè)點(diǎn)間距nx=81; % x方向測(cè)點(diǎn)數(shù)ny=81; % y方向測(cè)點(diǎn)數(shù)xmin=-200; % x方向起點(diǎn)ymin=-200; % y方向起點(diǎn)x=xmin:dx:(xmin+(nx-1)*dx); % x方向范圍y=ymin:dy:(ymin+(ny-1)*dy); % y方向范圍x,y=meshgr

4、id(x,y); % 轉(zhuǎn)化為排列 % 球體參數(shù)i=pi/3; %磁化傾角ia=0; %剖面磁方位角r=10; % 球體半徑 mv=4/3*pi*r3u=4*pi*10(-7);%磁導(dǎo)率t=0.5*10(-4);%地磁場(chǎng)強(qiáng)度k=0.2;%磁化率m=k*t/u; %磁化強(qiáng)度 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.

5、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),co

6、ntourf(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);xlabel(x(m),ylabel(y(m),title(理論球體za異常);axis equal,axis(-50 5

7、0 -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),shading interp,xlabel(x(m),ylabel(y(m),zlabel(理論球體hax異常),colorbar;subplot(222),mesh(x,y,hay),shading interp,xlabel(x(m),yl

8、abel(y(m),zlabel(理論球體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);

9、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-,linewi

10、dth,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,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 on

11、plot(x,za2,r-,linewidth,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(z

12、a);legend(h,boxoff); 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;% 測(cè)點(diǎn)分布范圍dx=5; % x方向測(cè)點(diǎn)間

13、距dy=5; % y方向測(cè)點(diǎn)間距nx=81; % x方向測(cè)點(diǎn)數(shù)ny=81; % y方向測(cè)點(diǎn)數(shù)xmin=-200; % x方向起點(diǎn)ymin=-200; % y方向起點(diǎn)x=xmin:dx:(xmin+(nx-1)*dx); % x方向范圍y=ymin:dy:(ymin+(ny-1)*dy); % y方向范圍x,y=meshgrid(x,y); % 轉(zhuǎn)化為排列 % 水平圓柱體參數(shù)i=pi/3; %磁化傾角a=0;%剖面磁方位角is=(tan(tan(i)*sec(a)(-1);r=10; % 圓柱體橫截面半徑 ms=pi*r2; %圓柱體橫截面面積u=4*pi*10(-7); %磁導(dǎo)率t=0.5*1

14、0(-4);%地磁場(chǎng)強(qiáng)度k=0.2;%磁化率m=k*t/u; %磁化強(qiáng)度 a/mms=m*(cos(i)*cos(a)2+(sin(i)2);m=ms*s; %單位長(zhǎng)度的有效磁矩d=30; % 圓柱體中心點(diǎn)埋深 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*

15、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),ylabel(y(m),title(理論圓柱體ha異常);axis equa

16、l,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,za),shading interp,xlabel(x(m),ylabel(y(m),zlabel(理論圓柱體za異常),colorbar;subplot(222),mesh(x,y,ha),sh

17、ading interp,xlabel(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,

18、1.3),xlabel(x(m),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*

19、is-pi/2)-2*d*x*cos(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

20、(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); 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(磁力異常（半徑）),grid on;endh=legend(za);legend(h,boxoff);6、 實(shí)驗(yàn)結(jié)果：球體實(shí)驗(yàn)結(jié)果：平面等值線圖：曲面圖：主剖面視圖：球體參數(shù)改變后的主剖面視圖：