实验2磁性体磁场正演程序.docx
- 文档编号:25830883
- 上传时间:2023-06-15
- 格式:DOCX
- 页数:20
- 大小:176.21KB
实验2磁性体磁场正演程序.docx
《实验2磁性体磁场正演程序.docx》由会员分享,可在线阅读,更多相关《实验2磁性体磁场正演程序.docx(20页珍藏版)》请在冰豆网上搜索。
实验2磁性体磁场正演程序
《应用地磁学》实验报告
姓名:
张嘉琪
学号:
1010112225
指导教师:
李淑玲
实验地点:
实验室319
实验日期:
2014-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;%磁化倾角i
a=0;%剖面磁方位角
R=10;%球体半径m
v=4/3*pi*R^3
u=4*pi*10^(-7);%磁导率
T=0.5*10^(-4);%地磁场强度
k=0.2;%磁化率
M=k*T/u;%磁化强度A/m
m=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/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));
holdon
plot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁倾角改变)'),gridon;
end
h=legend('Za');legend(h,'boxoff');
figure(5),clf,
fora=0:
pi/6:
pi
A=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));
holdon
plot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(磁方位改变)'),gridon;
end
h=legend('Za');legend(h,'boxoff');
figure(6),clf,
fori=pi/3;a=0;
R=10:
5:
20
v=4/3*pi*R^3
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));
holdon
plot(x,Za2,'r-','linewidth',1.3),xlabel('X(m)'),ylabel('磁力异常(球体半径)'),gridon;
end
h=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;%圆柱体横截面半径m
S=pi*R^2;%圆柱体横截面面积
u=4*pi*10^(-7);%磁导率
T=0.5*10^(-4);%地磁场强度
k=0.2;%磁化率
M=k*T/u;%磁化强度A/m
Ms=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:
200
Za1=(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异常');
end
subplot(312)
forx=-200:
5:
200
Ha1=(-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异常');
end
subplot(313)
forx=-200:
5:
200
T1=(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/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);
holdon
plot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常(磁倾角)'),gridon;
end
h=legend('Za');legend(h,'boxoff');
figure(5),clf,
fora=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);
holdon
plot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常(方位角)'),gridon;
end
h=legend('Za');legend(h,'boxoff');
figure(6),clf,
forR=10:
5:
20
i=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);
holdon
plot(X,Za1,'r-*'),xlabel('Y(m)'),ylabel('磁力异常(半径)'),gridon;
end
h=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对基本的规则球体与水平圆柱体的平面等值线与剖面图的绘制,并了解了磁异常的基本分布规律;,了解影响磁性体磁场的主要因素(如磁性体的形体、物性参数、走向或计算剖面的选择等),同时培养了实际动手能力与分析问题的能力。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验 磁性 磁场 程序