磁性体磁场正演
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
磁法勘探
上
机
实
验
报
告
*名:***
学号: ********** 指导教师:***
日期:2020.4.15
实验二:磁性体磁场正演
一、实验目的:
1、通过球体、水平圆柱体磁场的正演计算,掌握简单规则磁性体正演磁场的计算方法;
2、通过计算认识球体与水平圆柱体磁场的一般分布规律,了解影响磁性体磁场的主要因素(如磁性体的形体、物性参数、走向或计算剖面的选择等),培养学生实际动手能力与分析问题的能力。
二、实验内容
用Matlab 语言或C 语言编程实现球体和水平圆柱体的磁场(包括Za 、Ha 、Δt)的正演计算。
三、实验要求
假设地磁场方向与磁性体磁化强度方向一致且均匀磁化的情况下,当地磁场T=50000nT ,磁倾角I=60°,球体与水平圆柱体中心埋深R=30m ,半径r=10m ,磁化率k=0.1(SI ),计算(观测)剖面磁化强度水平投影夹角A ′=0°时:
1、正演计算球体的磁场(Za 、Hax 、Hay 、ΔT ),画出对应的平面等值线图、曲面图及主剖面异常图;
2、正演计算水平圆柱体的磁场(Za 、Ha 、ΔT ),画出主剖面异常结果图;
3、通过改变球体或水平圆柱体的几何参数、磁化强度方向(I )、计算剖面的方位角(A ′),观察主剖面磁场Za 的变化,分析磁化方向与计算剖面对磁性体磁场特征的影响。
四、实验原理
球体与水平圆柱体磁场(Za 、Ha 、ΔT )的计算公式是以磁化强度倾角I 、有效磁化倾角is 和剖面与磁化强度水平投影夹角A ′来表达。
1、球体磁场的正演公式: [[[⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪
⎪⎬⎫'-'---++='+-'--++='+-'--⋅++=
]sin cos 3cos cos 3sin )2( )(4]cos cos 3sin 3sin cos )2( )(4]sin cos 3sin 3cos cos )2( )(42222/52220
2222/522202222/52220A I Ry A I Rx I y x R R y x m Z A I xy I Ry A I R x y R y x m H A I xy I Rx A I R y x R y x m H a ay ax πμπμπμ
()]sin 2sin 32sin cos 3cos 2sin 3sin cos )2(cos cos )2(sin )2[(42222222
222222222/52
220A I yR A I xy A I xR A I R x y A I R y x I y x R R y x m
T '-'+'
-'--+'--+--++=∆πμ2、水平圆柱体磁场的正演公式:
⎪⎪⎭
⎪⎪⎬⎫+-+-=--+=]sin 2cos )[()(12]cos 2sin )[()(12222220222220s s s a s s s a i Rx i x R R x m H i Rx i x R R x m Z πμπμ()()()()[]
902cos 2902sin sin sin 2222220----+=
∆s s s s
i Rx i x R i I R x m T πμ 3、有效磁化强度Ms 与有效磁化倾角is :
⎪⎭
⎪⎬⎫'==+'=+=--)sec ()sin cos (cos )(112222/122A tgI tg M M tg i I A I M M M M x z s z x s 五、实验报告
(内容包括实验目的、实验内容、实验原理、计算程序代码、实验结果、结果分析或小结)
计算程序代码:
球体:dx=5;
dy=5;
nx=81;
ny=81;
xmin=-200;
ymin=-200;
x=xmin:dx:(xmin+(nx-1)*dx);
y=ymin:dy:(ymin+(ny-1)*dy);
[X,Y]=meshgrid(x,y);
T=5*10^(-5);
a=0;
i=pi/3;
R1=10;
v1=4/3*pi*R1^3;
u=4*pi*10^(-7);
k=0.1;
M=k*T/u;
m=M*v1;
D=30;
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=(u*m*((2*X.^2-Y.^2-D.^2)*cos(i)*cos(a)-3*D*X.*sin(i)+3*D*X.*cos(i )*sin(a)))./(4*pi*(X.^2+Y.^2+D.^2).^(5/2));
Hay=(u*m*((2*X.^2-Y.^2-D.^2)*cos(i)*sin(a)-3*D*X.*sin(i)+3*D*X.*cos(i )*cos(a)))./(4*pi*(X.^2+Y.^2+D.^2).^(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异常');
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]),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),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.^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);