实验3-磁异常处理与转换(张嘉琪)

《应用地磁学》实验报告

姓名：张嘉琪

学号： 1010112225

指导教师：李淑玲

实验地点：实验室319

实验日期： 2014-05-31

实验三：磁异常处理与转换

一、实验目的：

1、加深对磁性体磁异常在空间域处理转换原理与作用的认识

2、用Matlab 语言编程实现水平圆柱体磁异常(包括Za 、Ha 、Δt)的向上延拓和分量转换，培养学生数据处理的实际动手能力。

二、实验内容

利用两个大小与埋深不同的水平圆柱体产生的磁异常（ΔT 、Za 、Ha ），进行上延计算与分量转换计算，异常数据利用实验一中水平圆柱体的正演程序计算而来。

三、实验要求

球体或水平圆柱体磁异常上延计算与分量计算可任选其一，但均要求对计算结果误差进行分析，具体要求如下：

1、上延计算：对ΔT 、Za 、Ha 各分量向上延拓5m 、10m ，画出对应结果图；

2、分量计算：进行Za →Ha ，Ha →Za 分量转换，画出对应结果图；

3、观察向上延拓、分量转换前后的异常特征，分析向上延拓与分量转换的作用；

4、要求对计算结果误差进行分析，如有可能提出改进措施。

5

152535455565758595105115125135145155 -30

-20

-10

0 10 20

30

(n T ) 0

10 20

30

(m )

图1、水平圆柱体模型及其产生的磁异常分量Za 和Ha

四、实验原理

在空间域内讨论磁异常的转换和处理的基础是磁异常的位函数，它具有调和函数的性质。因此可根据某观测面上的实测磁异常，换算成场源以外其他空间位置的磁异常，也可以将某种实测分量换算成其他的分量，增加解释信息。

1、向上延拓：换算平面位于实测平面之上。主要用途是削弱局部异常干扰，

反映深部异常。设坐标原点位于计算点下方实测剖面上，延拓高度为一个点距

h ，则原点的向上延拓公式：

1()21222()2(,0)14(0,)(,0)arctan 43n h n h n n Za nh h Za h d Za nh h n ξπξπ∞∞+-=-∞

=-∞-=

=++∑∑? 1()21222()2

(,0)14(0,)(,0)arctan 43n h n h n n Ha nh h Ha h d Ha nh h n ξπξπ∞∞+-=-∞=-∞-==++∑∑? 344arctan 1)0,()0,(),0(2)21()21(22+?=+?=

-??∑∑+-∞-∞

=∞-∞=n nh T d h h nh T h T h n h n n n πξξπ 2、磁异常分量间的换算：原点处的磁异常的分量换算公式

Za →Ha ：)]0,()0,([)0,0(1i a i a N

i i a Z Z a H ξξ--=∑=，

Ha →Za ：)]0,()0,([)0,0(1i a i a N

i i a H H a Z ξξ---=∑=， 式中，)ln 211(1121ξξ+=πa ，i i i a ξξπ1ln 22+=，)ln 1(211

-+=N N N a ξξπ。 五、计算程序代码

水平圆柱体向上延拓：

% 两个水平圆柱体剖面异常向上延拓

% 重点理解向上延拓的作用和计算误差来源。

clc;

clear;

% 测点分布范围

dx=5; % X 方向测点间距 m

nx=101; % X 方向测点数

xmin=-250; % X 方向起点 m

x=xmin:dx:(xmin+(nx-1)*dx); % X 方向范围 m

% 两个水平圆柱体参数

i=pi/2; %有效磁化倾角is

R1=5; % 水平圆柱体1半径 m

R2=50; % 水平圆柱体2半径 m

v1=pi*R1^2;

v2=pi*R2^2;

u=4*pi*10^(-7); %磁导率

Ms=0.2; %磁化强度 A/m

m1=Ms*v1; %磁矩

m2=Ms*v2; %磁矩

D1=25; % 水平圆柱体1埋深 m

D2=100; % 水平圆柱体2埋深 m

% 园柱体Za理论磁异常

Za1=(u*m1*((D1.^2-(x-50).^2)*sin(i)-2*D1*(x-50).*cos(i)))./(2*pi*((x-50).^2+D1. ^2).^2);

Za2=(u*m2*((D2.^2-(x+50).^2)*sin(i)-2*D2*(x+50).*cos(i)))./(2*pi*((x+50).^2+D2. ^2).^2);

Za=Za1+Za2;

% 圆柱体Ha理论磁异常

Ha1=(-u*m1*((D1.^2-(x-50).^2)*cos(i)+2*D1*(x-50).*sin(i)))./(2*pi*((x-50).^2+D1 .^2).^2);

Ha2=(-u*m2*((D2.^2-(x+50).^2)*cos(i)+2*D2*(x+50).*sin(i)))./(2*pi*((x+50).^2+D2 .^2).^2);

Ha=Ha1+Ha2;

%圆柱体ΔT理论异常

T1=(u*m1*((D1.^2-(x-50).^2)*cos(2*i-pi)-2*D1*(x-50).*sin(2*i-pi/2)))./(((D1.^2-(x-50).^2)*sin(2*i-pi/2)-2*D1*(x-50).*cos(2*i-pi/2)));

T2=(u*m2*((D2.^2-(x+50).^2)*cos(2*i-pi)-2*D2*(x+50).*sin(2*i-pi/2)))./(((D2.^2-(x+50).^2)*sin(2*i-pi/2)-2*D2*(x+50).*cos(2*i-pi/2)));

T=T1+T2;

% Za 向上延拓

Zau1=ones(1,nx)*(0/0); % 上延后异常值，先赋空值

h=3; % 延拓高度,必须为点距dx的倍数

n=10; % 级数，即参与积分的区间段（积分点数-1）/2

for i=(h*n+1):(nx-h*n)

tmp=0;

for j=(i-h*n):h:(i+h*n)

k=(j-i)/h % 参与积分的区间段（积分点数-1）/2

tmp=tmp+Za(j)*atan(4/(4*k*k+3))/pi;

end

Zau1(i)=tmp; % 上延后异常值

end

figure(1),clf,plot(x,Za,'b',x,Zau1,'r:'),xlabel('X (m)'),ylabel('Za磁异常(nT.)'),

legend('原始异常','上延3m异常'),title('向上延拓Za磁异常');

% Ha 向上延拓

Hau1=ones(1,nx)*(0/0);

h=3;

n=10;

for i=(h*n+1):(nx-h*n)

tmp=0;

for j=(i-h*n):h:(i+h*n)

k=(j-i)/h

tmp=tmp+Ha(j)*atan(4/(4*k*k+3))/pi;

end

Hau1(i)=tmp;

end

figure(2),clf,plot(x,Ha,'b',x,Hau1,'r:'),xlabel('X (m)'),ylabel('Ha磁异常(nT.)'),

legend('原始异常','上延3m异常'),title('向上延拓Ha磁异常');

% T 向上延拓

Tu1=ones(1,nx)*(0/0);

h=3;

n=10;

for i=(h*n+1):(nx-h*n)

tmp=0;

for j=(i-h*n):h:(i+h*n)

k=(j-i)/h

tmp=tmp+T(j)*atan(4/(4*k*k+3))/pi;

end

Tu1(i)=tmp;

end

figure(3),clf,plot(x,T,'b',x,Tu1,'r:'),xlabel('X (m)'),ylabel('ΔT磁异常(nT.)'),

legend('原始异常','上延3m异常'),title('向上延拓ΔT磁异常');

2、水平圆柱体Za—Ha互算

clc，clear，

fprintf('\n柱体1深25m，半径5m,柱体2深100m，半径50m，测点数为101,点距选取5m\n'); b=2;

fprintf('\n点距取%dm\n',b);

figure('color','w','NumberTitle','off','name','za-→ha');

x2=1:5:505; %点距为bm

x1=x2-101*5/2;

%柱体参数

D1=25; %柱深

D2=100;

R1=5; %柱体半径 m

R2=50;

s1=pi*(R1^2);%柱体截面积

s2=pi*(R2^2);

u=4*pi*10^(-7); %磁导率

Ms=0.2; %磁化强度

m1=Ms*s1; %磁矩

m2=Ms*s2;

is=pi/2;

hold on

% 园柱体Za理论磁异常

Za1=(u*m1*((D1.^2-(x1-50).^2)*sin(i)-2*D1*(x1-50).*cos(i)))./(2*pi*((x1-50).^2+ D1.^2).^2);

Za2=(u*m2*((D2.^2-(x1+50).^2)*sin(i)-2*D2*(x1+50).*cos(i)))./(2*pi*((x1+50).^2+ D2.^2).^2);

Za=Za1+Za2;

% 圆柱体Ha理论磁异常

Ha1=(-u*m1*((D1.^2-(x1-50).^2)*cos(i)+2*D1*(x1-50).*sin(i)))./(2*pi*((x1-50).^2 +D1.^2).^2);

Ha2=(-u*m2*((D2.^2-(x1+50).^2)*cos(i)+2*D2*(x1+50).*sin(i)))./(2*pi*((x1+50).^2 +D2.^2).^2);

Ha=Ha1+Ha2;

%Za向Hx转换

c=[0.4268 0.1749 0.1103 0.0813 0.0645 0.0536 0.0458 0.0400 0.0355 0.1759];

n=length(c); %转换系数个数

m=length(Za);

for i=(n+1):(m-n)

Haz=0;

for j=1:n

Haz=Haz+c(j)*(Za(i+j)-Za(i-j));

end

Hazh(i)=Haz;

end

for i=1:80

Hazh1(i)=Hazh(i+10);

end

x=((n+1)*5):5:(100*5-n*5);

figure(1),clf,

plot(x2, Za,'.-m',x2, Ha,'.-g',x,Hazh1,'.-'),title('Za转换为Ha'), xlabel(' X剖面走向/m'), ylabel('磁异常nT');

legend('Za','Ha','Za转换后的Ha');

六、实验结果

1、Za向上延拓3m：

2、Ha向上延拓3m：

3、ΔT向上延拓3m：

4、Za转换Ha

七、结果分析

1、二维向上延拓

在50m处25米深有一水平柱体，半径为5m，-50m处100米深有一柱体，半径为50m。由出的图可知，向上延拓可以突出深部异常趋势，压制浅部矿体的异常。但是随着延拓高度的增加，异常值逐渐减小。向上延拓的高度越高则异常曲线越接近深部场源的异常，异常峰值较好的反映了深部异常体的位置，但是异常值明显减小。

2、磁异常分量转换

由图可知实际的Ha异常曲线与转换后的磁异常曲线并不重合，但大体趋势一致，主要是跟点距及测点的选择有关，点距越小，测点越多，则转换出来的Ha曲线越接近真实值。

3、小结：通过这次上机实验，我加深对磁性体磁异常在空间域处理转换原理与作用的认识，同时可以使用Matlab语言编程实现水平圆柱体磁异常(包括Za、Ha、Δt)的向上延拓和分量转换。