频率域磁异常处理与转换

磁法勘探上机实验报告姓名：王长阔 学号：1010173229 指导教师：李淑玲 日期：2020.4.20一、实验目的1、加深对磁性体磁异常在频率域处理转换原理与作用的认识；2、用Matlab语言编程实现频率域磁异常处理与转换，如向上延拓、导数计算、∆T磁异常化极等处理，培养程序开发与数据处理的动手能力。二、实验内容1、利用两个大小与埋深不同的球体产生的平面磁异常数据（如ΔT、Za），选择频率域向上延拓算子、导数计算算子、化磁极算子，通过频率域实现磁异常的处理与转换；2、磁异常数据准备：自行准备，利用正演实验计算球体的磁异常数据（ΔT、Za）；3、频率域处理转换算子：向上延拓算子、垂向二阶导数算子、化磁极算子。1）球体正演参数自行设定，也可参考以下参数设置：（1）正演2个球体的叠加异常（如ΔT、Za），用于向上延拓、垂向二阶导数计算、化磁极的数据源；（2）平面磁异常计算范围：x：-200m至200m，y：-200m至200m，地磁场总强度T=5*10-9（3）假设球体1的参数：半径r1=2m,球心坐标（50m,0m,5m）,磁化强度（4）假设球体2的参数：半径r1=15m,球心坐标（-50m,0m,30m）,磁化强度2）如上数据可以存储为grd文件，处理转换时直接调用即可。三、实验要求球体某分量磁异常上延计算、导数计算、化磁极可仟选其一，具体要求如下：1、上延计算：对ΔT或Za向上延拓不同高度，画出对应结果图；2、导数计算：利用ΔT或Za异常计算垂向二阶导数，画出对应结果图；3、化磁极：利用ΔT数据进行化极处理，画出对应结果图；4、观察转换前后的异常特征，分析这些处理转换的作用。四、实验原理在频率域实现磁异常处理转换计算简便，速度快。通过频率域处理转换的基本过程为：1）先将对磁异常进行傅里叶变换求磁异常频谱；2）将磁异常频谱与频域处理转换的算子相乘，得到转换后的频谱；3）对转换后的频谱进行反傅里叶变换得到磁异常。频率域向上延拓的计算公式：∆频率域磁异常垂向二阶导数的计算公式：∂∆频率域化磁极的计算公式：（不考虑剩磁，即磁化方向与地磁场方向一致）Zl0、m0、五．计算程序代码%测点分布范围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/4;%有效磁化倾角isa=0;%剖面磁方位角r1=2;%球体半径mr2=15;%球体半径mv1=(4*pi*r1^3)/3;v2=(4*pi*r2^3)/3;u=4*pi*10^(-7);%磁导率M1=0.1;%磁化强度A/mM2=0.2;%磁化强度A/mm1=M*v1;%磁矩m2=M2*v2;%磁矩D1=5;%球体1埋深mD2=30;%球体2埋深m%球体Za理论磁异常Za1=(u*m1*((2*D1.^2-(X-50).^2-Y.^2)*sin(i)-3*D1*(X-50).*cos(i)*cos(a)-3*D1*Y.*cos(i)*sin(a)))./(4*pi*((X-50).^2+Y.^2+D1.^2).^(5/2));Za2=(u*m2*((2*D2.^2-X.^2-Y.^2)*sin(i)-3*D2*X.*cos(i)*cos(a)-3*D2*Y.*cos(i)*sin(a)))./(4*pi*(X.^2+Y.^2+D2.^2).^(5/2));Za=Za1+Za2;figure(1),clf,pcolor(X,Y,Za);shadinginterp,xlabel('x(m)'),ylabel('y(m)'),title('理论球体Za异常');%磁异常向上延拓n=512;%傅立叶变换充零后的点数Hu=30;%向上延拓高度nn=n/2;%二维傅立叶变换Zaf=fft2(Za,n,n);%径向频率dr=1/((n-1)*dx);%基本频率、采样间隔r=0.0:dr:((nn-1)*dr);%频率范围%计算第一象限延拓滤波因子fori=1:nnrx=(i-1)*dr;forj=1:iry=(j-1)*dr;rr=(rx*rx+ry*ry).^0.5;w(j,i)=exp(2*pi*rr.*(-Hu));w(i,j)=w(j,i);endendfori=(nn+1):nforj=1:nnw(j,i)=w(j,n+1-i);%计算第三象限延拓滤波因子w(i,j)=w(n+1-i,j);%计算第二象限延拓滤波因子w(i,nn+j)=w(n+1-i,nn+1-j);%计算第四象限延拓滤波因子endend%乘以上延因子进行上延计算Zafu=Zaf.*w;%傅立叶反变换Za1=real(ifft2(Zafu));Zau=Za1(1:ny,1:nx);figure(2),clf,pcolor(X,Y,Zau);shadinginterp,xlabel('x(m)'),ylabel('y(m)'),title('向上延拓30米后Za异常');Hu=60;%向上延拓高度nn=n/2;%二维傅立叶变换Zaf=fft2(Za,n,n);%径向频率dr=1/((n-1)*dx);%基本频率、采样间隔r=0.0:dr:((nn-1)*dr);%频率范围%计算第一象限延拓滤波因子fori=1:nnrx=(i-1)*dr;forj=1:iry=(j-1)*dr;rr=(rx*rx+ry*ry).^0.5;w(j,i)=exp(2*pi*rr.*(-Hu));w(i,j)=w(j,i);endendfori=(nn+1):nforj=1:nnw(j,i)=w(j,n+1-i);%计算第三象限延拓滤波因子w(i,j)=w(n+1-i,j);%计算第二象限延拓滤波因子w(i,nn+j)=w(n+1-i,nn+1-j);%计算第四象限延拓滤波因子endend%乘以上延因子进行上延计算Zafu=Zaf.*w;%傅立叶反变换Za1=real(ifft2(Zafu));Zau=Za1(1:ny,1:nx);figure(3),pcolor(X,Y,Zau),shadinginterp,xlabel('x(m)'),ylabel('y(m)'),title('向上延拓60米后Za异常');Hu=90;%向上延拓高度nn=n/2;%二维傅立叶变换Zaf=fft2(Za,n,n);%径向频率dr=1/((n-1)*dx);%基本频率、采样间隔r=0.0:dr:((nn-1)*dr);%频率范围%计算第一象限延拓滤波因子fori=1:nnrx=(i-1)*dr;forj=1:iry=(j-1)*dr;rr=(rx*rx+ry*ry).^0.5;w(j,i)=exp(2*pi*rr.*(-Hu));w(i,j)=w(j,i);endendfori=(nn+1):nforj=1:nnw(j,i)=w(j,n+1-i);%计算第三象限延拓滤波因子w(i,j)=w(n+1-i,j);%计算第二象限延拓滤波因子w(i,nn+j)=w(n+1-i,nn+1-j);%计算第四象限延拓滤波因子endend%乘以上延因子进行上延计算Zafu=Zaf.*w;%傅立叶反变换Za1=real(ifft2(Zafu));Zau=Za1(1:ny,1:nx);figure(4),pcolor(X,Y,Zau);shadinginterp,xlab