实验3-磁异常处理与转换

《应用地磁学》实验报告名：学号：指导教师：实验地点：实验室319实验日期：2014-05-31

实验三：磁异常处理与转换一、实验目的：1、加深对磁性体磁异常在空间域处理转换原理与作用的认识2、用Matlab语言编程实现水平圆柱体磁异常（包括Za、Ha.At）的向上延拓和分量转换，培养学生数据处理的实际动手能力。二、实验内容利用两个大小与埋深不同的水平圆柱体产生的磁异常（AT、Za、Ha），进行上延计算与分量转换计算，异常数据利用实验一中水平圆柱体的正演程序计算而来。三、实验要求球体或水平圆柱体磁异常上延计算与分量计算可任选其一，但均要求对计算结果误差进行分析，具体要求如下：1、上延计算：对AT、Za、Ha各分量向上延拓5m、10m，画出对应结果图；2、分量计算：进行Za今Ha，Ha»Za分量转换，画出对应结果图；3、观察向上延拓、分量转换前后的异常特征，分析向上延拓与分量转换的作用；4、要求对计算结果误差进行分析，如有可能提出改进措施。302010TT0,-20-30四、实验原理TT0,-20-30在空间域内讨论磁异常的转换和处理的基础是磁异常的位函数，它具有调和函数的性质。因此可根据某观测面上的实测磁异常，换算成场源以外其他空间位置的磁异常，也可以将某种实测分量换算成其他的分量，增加解释信息。1、向上延拓：换算平面位于实测平面之上。主要用途是削弱局部异常干扰，

反映深部异常。设坐标原点位于计算点下方实测剖面上，延拓高度为一个点距h,则原点的向上延拓公式：TOC\o"1-5"\h\z"小7、VZa(nh,0)[(n,1)hh化\14Za(0,-h)=乙J(2)hdg=乙Za(nh,0)—arctan兀(n--1)hg2+h2兀4n2+3n=-s2n=-syHa(nh,0)J(n+1)hh亦\14Ha(0,-h)=乙J(2)dg=乙Ha(nh,0)—arctan兀(n-1)hg2+h2兀4n2+3n=-s2n=-sAT(0,-h)=y^h0g工dg=:Eat(nh,0)Larctan—兀(n-1)hg2+h2兀4n2+3n=-sk2n=-s2、磁异常分量间的换算：原点处的磁异常的分量换算公式Za今Ha：H(0,0)=Ea.[Z(g.,0)-Z(-g.,0)]，i=1aiHa今Za：Z(0,0)=-yLa^[H(g,0)-H(-g,0)],i=1式中，a1n1-11glnTai式中，a1n1-11glnT2g1=Am2兀1g=(1+lnd)。2兀gN-1%两个水平圆柱体剖面异常向上延拓%重点理解向上延拓的作用和计算误差来源。clc;clear;%测点分布范围dx=5;%X方向测点间距mnx=101;%X方向测点数xmin=-250;%X方向起点mx=xmin:dx:(xmin+(nx-1)*dx);%X方向范围m%两个水平圆柱体参数i=pi/2;%有效磁化倾角isR1=5;%水平圆柱体1半径mR2=50;%水平圆柱体2半径mv1=pi*R「2;v2=pi*R2"2;u=4*pi*10”(-7);%磁导率Ms=0.2;%磁化强度A/mm1=Ms*v1;%磁矩m2=Ms*v2;%磁矩D1=25;%水平圆柱体1埋深mD2=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+D1ff;Ha2=(-u*m2*((D2."2-(x+50).”2)*cos(i)+2*D2*(x+50).*sin(i)))./(2*pi*((x+50)."2+D2ff;Ha=Ha1+Ha2;%圆柱体AT理论异常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)/2fori=(h*n+1):(nx-h*n)tmp=0;forj=(i-h*n):h:(i+h*n)k=(j-i)/h%参与积分的区间段(积分点数-1)/2tmp=tmp+Za(j)*atan(4/(4*k*k+3))/pi;endZau1(i)=tmp;%上延后异常值endfigure(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;fori=(h*n+1):(nx-h*n)tmp=0;forj=(i-h*n):h:(i+h*n)k=(j-i)/htmp=tmp+Ha(j)*atan(4/(4*k*k+3))/pi;endHau1(i)=tmp;endfigure(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;fori=(h*n+1):(nx-h*n)tmp=0;forj=(i-h*n):h:(i+h*n)k=(j-i)/htmp=tmp+T(j)*atan(4/(4*k*k+3))/pi;endTu1(i)=tmp;endfigure(3),clf,plot(x,T,'b',x,Tu1,'r:'),xlabel('X(m)'),ylabel('△T磁异常(nT.)'),legend('原始异常','上延3m异常'),title('向上延拓AT磁异常');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;%点距为bmx1=x2-101*5/2;%柱体参数D1=25;%柱深D2=100;R1=5;%柱体半径mR2=50;s1=pi*(R「2);%柱体截面积s2=pi*(R2"2);u=4*pi*10”(-7);%磁导率Ms=0.2;%磁化强度m1=Ms*s1;%磁矩m2=Ms*s2;is=pi/2;holdon%园柱体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.42680.17490.11030.08130.06450.05360.04580.04000.03550.1759];n=length(c);%转换系数个数m=length(Za);fori=(n+1):(m-n)Haz=0;forj=1:nHaz=Haz+c(j)*(Za(i+j)-Za(i-j));endHazh(i)=Haz;endfori=1:80Hazh1(i)=Hazh(i+10);endx=((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：

-250-200-160-100-50060X(m)32.6向上延拓及磁异常100150200260261H—0.60-0.5——原始异常一..上延3m异常3、^-250-200-160-100-50060X(m)32.6向上延拓及磁异常100150200260261H———原始异常一..上延3m异常3、^T向上延拓3m：_x1013向上延拓AT磁异常□.二IIIIIIIIII-原始异常3-上延3m异常2.6-2-1.5-TOC\o"1-5"\h\z1-.\o"CurrentDocument"0.6-:■::|—I....•:■-<-■■■'■-□5111111111-250-200-160-100-50060100150200260X(m)4、Za转换Ha七、结果分析1、二维向上延拓在50m处25米深有一水平柱体，半径为5m，-50m处100米深有一柱体，半径为50m。由出的图可知，向上延拓可以突出深部异常趋势，压制浅部矿体的异常。但是随着延拓高度的增加，异常值