实验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、球体磁场的正演公式： [[[⎪⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪

⎪⎬⎫'-'---++='+-'--++='+-'--⋅++=

]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 五、计算程序代码：

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异常');

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),shading

interp,xlabel('X(m)'),ylabel('Y(m)'),zlabel('理论球体Hax异常'),colorbar; subplot(222),mesh(X,Y,Hay),shading

interp,xlabel('X(m)'),ylabel('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) );

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 ));