地磁作业(转换、延拓)

应用地磁学实验报告实验2——磁异常转换计算学号： 10105218 姓名：朱占升一、实验目的1、掌握水平圆柱体磁场异常分布；2、用Matlab实现水平圆柱体的磁异常场正演计算；3、利用正演结果进行磁异常分量之间的换算；4、通过程序换算认知测点间距即采样点数对换算效果的影响；5、加深对磁法勘探的理解认识；二、程序代码%磁法异常换算%剖面为北向A’=0度，则有I=is。

%所测数据均在同一水平面，柱体深30m，半径8m,测点数为100,%测点间距依次选取2 4 6 8mclcclearfprintf('\n柱体深30m，半径8m,测点数为100,点距依次选取2 4 6 8m\n');for b=2:2:8fprintf('\n点距取%dm\n',b);figure('color','w','NumberTitle','off ','name','za-→ha');x2=1:b:100*b; %点距为bmx1=x2-100*b/2;h=30; %柱深R=8; %柱体半径 ms=pi*(R^2);%柱体截面积k=0.2; %磁铁矿磁化率u=4*pi*10^(-7); %磁导率B=50000; %nT磁感应强度H=B/u ; %磁化场强度M=k*H; %磁化强度m=M*s; %磁矩a=0; %剖面为北向A’=0度.I=90/180*pi; %倾斜角is=atan(tan(I)*csc(pi/2-a));hold onza=u*m*((h.^2-x1.^2)*sin(is)-2*h*x1.*cos(is))./(2*pi*(x1.^2+h.^2).^2);hax=-1*u*m*((h.^2-x1.^2)*cos(is)+2*h*x1.*sin(is))./(2*pi*(x1.^2+h.^2).^2);plot(x2, za,'.-m');plot(x2, hax,'.-g');title('za转换为hax')xlabel(' X剖面走向/m');ylabel('磁异常nT');c=[0.4268 0.1749 0.1103 0.0813 0.06450.0536 0.0458 0.0400 0.0355 0.1759];n=length(c); %转换系数个数m=length(za);for i=(n+1):(m-n)haxz=0;for j=1:nhaxz=haxz+c(j)*(za(i+j)-za(i-j));endhaxzh(i)=haxz;endfor i=1:80haxzh1(i)=haxzh(i+10);endx=((n+1)*b):b:(100*b-n*b);plot(x,haxzh1,'.-')legend('za','hax','za→hax');end三、输出图形及解释通过图形可以知道当点距大于6m是转换出的H a与正演的H a有明显的图形失真（假异常），随着测点间距的加大，失真越厉害。

分析其原因：随着测点间距加大，单位距离内的采样点数逐渐减少，即采样密度逐渐降低。

获取的有效信息不足造成转换后的图形失真，而且推测当采样间距达到一定值之后，转换将会产生不稳定的情况。

测点间距即采样密度与地下的场源密切相关，它与场源体的几何参数、物性参数有着直接的关系。

本次实验中选用了水平圆柱体，其几何物性参数如程序中所设，在这种条件下他的最佳采样点距为2～5m。

当超过这个范围则出现失真、假异常。

四、实验小结本次试验通过Matlab程序实现对水平圆柱体正演后，对其引起的磁异常场进行了分量变换。

通过设定不同的测量点距，认识了采样密度不同的情况下，磁异常场转换失真的特点。

本次实验加深了对空间域磁异常的处理与转换的理解，对公式有了更深入的理解。

对于通过运用计算机实现磁场问题的模拟有了更深刻的理解，也加强了对于模型设置等程序语言的锻炼。

通过实验自己分析问题解决问题的能力不断加强，在实践中知识得到了不断的丰富。

应用地磁学实验报告实验3——磁异常频率域延拓学号： 10105218 姓名：朱占升一、实验目的1、掌握水平圆柱体磁场异常三维分布；2、用Matlab实现水平圆柱体的磁异常场延拓计算；3、利用程序分析延拓效果与高度的关系；4、通过柱体三维延拓了解延拓的失真问题；5、加深对磁法勘探的理解认识；二、程序代码1、水平圆柱体二维向上延拓（浅部异常在右侧）% 10105218 朱占升%剖面为北向A’=0度，则有I=is。

%场源为不同深度的两水平圆柱体，延拓高度选取了50、80、100m三种高度clc;clear;% 测点分布范围b=5; % X方向测点间距X=-200:b:300;nx=length(X);a=0; %剖面为北向A’=0度.%浅部球体1参数I=90/180*pi; %倾斜角i=atan(tan(I)*csc(pi/2-a));%有效磁化倾角isR1=20; % 球体半径 ms1=pi*R1^2;u=4*pi*10^(-7); %磁导率k=0.2;B=50000;H=B/u;M=k*H ; %磁化强度 A/mm1=M*s1; %磁矩D1=30; % 球体埋深 m% 深部球体2参数I=90/180*pi; %倾斜角i=atan(tan(I)*csc(pi/2-a));%有效磁化倾角isR2=50; % 球体半径 mu=4*pi*10^(-7); %磁导率k=0.3;B=50000;H=B/u;M2=k*H ; %磁化强度 A/mm2=M2*s2; %磁矩D2=100; % 球体埋深 m%总Za磁异常x1=X-100;x2=X;Za1=u*m1*((D1.^2-x1.^2)*sin(i)-2*D1*x1.*cos(i))./(2*pi*(x1.^2+D1.^2).^2);Za2=u*m2*((D2.^2-x2.^2)*sin(i)-2*D2*x2.*cos(i))./(2*pi*(x2.^2+D2.^2).^2);Za=Za1+Za2;%总Hax磁异常Hax1=-1*u*m1*((D1.^2-x1.^2)*cos(i)+2*D1*x1.*sin(i))./(2*pi*(x1.^2+D1.^2).^2);Hax2=-1*u*m2*((D2.^2-x2.^2)*cos(i)+2*D2*x1.*sin(i))./(2*pi*(x1.^2+D2.^2).^2);Hax=Hax1+Hax2;figure('color','g','NumberTitle','off','name','单个柱体Za');plot(X,Za1,'.-',X,Za2,'.-');xlabel('x(m)'),ylabel('y(m)'),title('单个柱体Za'); xlabel('x剖面方向(m)'),ylabel('y(nT)');figure('color','g','NumberTitle','off','name','两柱体Za');plot(X,Za,'.-r'),title('总Za'); xlabel('x剖面方向(m)'),ylabel('y(nT)');%延拓部分figure('color','g','NumberTitle','off','name','延拓后Za');n=512; % 傅立叶变换充零后的点数Hu=50; % 向上延拓高度50mZayt=up(Za,b,nx,n,Hu); %延拓function函数upplot(X,Zayt,'.-');xlabel('x剖面方向(m)'),ylabel('y(nT)'),title('上延50、80、100m 后Za');Hu=80; % 向上延拓高度80mZayt=up(Za,b,nx,n,Hu);%延拓function函数uphold onplot(X,Zayt,'.-m');Hu=100; % 向上延拓高度100mZayt=up(Za,b,nx,n,Hu);%延拓function函数upplot(X,Zayt,'.-r');legend('上延50m','上延80m','上延100m');stem(0,-100)stem(100,-20)ylabel('Y坐标/m');x=[0 100];y=[-100 20];s=['柱体'];text(x,y,s,'color','g')2、水平柱体三维上延（浅部异常在左侧）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); % 转化为排列a=0;% 浅部柱体1参数I=90/180*pi; %倾斜角i=atan(tan(I)*csc(pi/2-a));%有效磁化倾角is R1=20; % 柱体半径 ms1=pi*R1^2;u=4*pi*10^(-7); %磁导率k=0.2;B=50000;H=B/u;M=k*H ; %磁化强度 A/mm1=M*s1; %磁矩D1=30; % 柱体埋深 mx1=X+100;x2=X;% 深部柱体2参数I=90/180*pi; %倾斜角i=atan(tan(I)*csc(pi/2-a));%有效磁化倾角is R2=50; % 柱体半径 ms2=pi*R2^2;u=4*pi*10^(-7); %磁导率k=0.3;B=50000;H=B/u;M2=k*H ; %磁化强度 A/mm2=M2*s2; %磁矩D2=100; % 柱体埋深 m% 两柱体Za总磁异常Za1=u*m1*((D1.^2-x1.^2)*sin(i)-2*D1*x1.*cos(i))./(2*pi*(x1.^2+D1.^2).^2);Za2=u*m2*((D2.^2-x2.^2)*sin(i)-2*D2*x2.*cos(i))./(2*pi*(x2.^2+D2.^2).^2);Za=Za1+Za2;% 两柱体Hax总磁异常Hax1=-1*u*m1*((D1.^2-x1.^2)*cos(i)+2*D1*x1.*sin(i))./(2*pi*(x1.^2+D1.^2).^2);Hax2=-1*u*m2*((D2.^2-x2.^2)*cos(i)+2*D2*x1.*sin(i))./(2*pi*(x1.^2+D2.^2).^2);Hax=Hax1+Hax2;subplot(221)surf(X,Y,Za),shading interp,xlabel('x(m)'),ylabel('y(m)'),title('两柱体总Za异常');colorbar;axis autosubplot(222)surf(X,Y,Hax),shading interp,xlabel('x(m)'),ylabel('y(m)'),title('两柱体总Hax异常');colorbar;axis auton=512; % 傅立叶变换充零后的点数Hu=80; % 向上延拓高度mZau=upcon(Za,dx,nx,ny,n,Hu); % Za向上延拓axis autosubplot(223)surf(X,Y,Zau),shading interp,xlabel('x(m)'),ylabel('y(m)'),title('上延80m两柱体Za 总磁异常');colorbar;axis autosubplot(224)surf(X,Y,Za2),shading interp,xlabel('x(m)'),ylabel('y(m)'),title('深部柱体Za磁异常');colorbar;axis auto三、输出图形及解释1、二维向上延拓（浅部异常在右侧）在0处100米深有一水平柱体，100处20米深有一柱体。