频率域位场处理和转换实验.doc

1、 重磁资料处理与解释实验二频率域位场处理和转换实验学 院：地测学院专业名称：勘查技术与工程学生姓名： 学生学号：指导老师：提交日期：2018 年 1 月 9 日二 0 一八年一月0目录1 基本原理 21.1 位场的方程 21.2 二维傅里叶变换及卷积性质 .2（1）傅里叶变换 .2（2）卷积性质 .21.3 频率域位场延拓原理 32 输入/输出数据格式设计 .32.1 输入数据格式设计 .32.2 输出数据格式设计 .32.3 参数文件数据格式设计 .33 总体设计 43.1 频率域位场处理与转换的一般步骤 .43.2 软件总体设计结果流程图 44 测试结果 54.1 测试参数 .5（1）向上

2、延拓 .5（2）向下延拓 .54.2 测试结果 .65 结论及建议 7附录：源程序代码 811 基本原理1.1 位场的方程由场论知识可知，位场方程分为两大类：有源的 Possion 方程 ，以及无源的 Laplace 方程 。02U02ULaplace 方程的第一边值问题 通常为 Dirichlet 问题，第二边值问题1|fS通常称为 Nueman 问题。若 P 点在 S 平面内称为内部问题，反之称2|fnUs为外部问题。由唯一性定理可知，Dirichlet 的内部和外部问题的解是唯一的，而 Nueman 内部问题的解不是唯一的，有一常数差，但其外部问题解是唯一的。外部问题的解的唯一性的原因：

3、 。0;rrnU无源区域位场可以表示为：（1-1 ）dsnGWpW41)(（1-2 ）zyxhdzzy ,*,2,x 23221.2 二维傅里叶变换及卷积性质（1）傅里叶变换（1-3 ）dxyyxgyxgFvuGevui)(2),(),(),( （1-4 ）Gvuyx vyuxi)(21),(),(),( （2）卷积性质2（1-5 ）vuPGyxpF,*,*,g（1-6 ） yxpvuPg11.3 频率域位场延拓原理当已知实测平面的异常时，换算场源以外的异常称之为延拓，分为向上延拓和向下延拓。半空间狄利克莱问题解析解：（1-7 ）zvueyxWFhzyxF2,其中： 称为延拓因子， 为计算面

4、Z 坐标 ,Z 轴向下为正方向，ezvu)(22为计算面频率域位场， 为延拓面的频率域位场。 ,yxWFzyx,2 输入/输出数据格式设计2.1 输入数据格式设计观测面位场数据保存在 filename_obser 文件中，为.grd 格式。计算延拓误差时的精确场值文件保存在 filename_obser2 中，为.grd 格式。2.2 输出数据格式设计实际计算得到的数据保存在 filename_output1 和文件 filename_output2 中，为.grd 格式。2.3 参数文件数据格式设计将所要读取的参数保存在一个文件中，该文件名变量为 cmdfile，字符串变量，长度不超过 80

5、，全路径名。在该文件中保存的参数如下：filename_obser:低高度观测面位场数据文件filename_output1：向上延拓后位场数据文件filename_output2：向下延拓后位场数据文件filename_obser2:高高度观测面数据文件3factor_m：扩边因子distance1:向上延拓的高度（m/z 轴向下为正方向）distance2:向下延拓的高度（m/z 轴向下为正方向）3 总体设计3.1 频率域位场处理与转换的一般步骤 即 可 。去 掉 扩 边 部 分 ， 输 出步第 ；进 行 反 变 换 得 到步第 ；计 算步第 ；以 及 方 向 转 换 因 子 、 延 拓

6、因 子计 算 导 数 因 子步第 ；对 扩 边 后 的 数 据 进 行步第 进 行 扩 边 处 理 ；对 网 格 化 的 数 据步第 ),(:6 ),(,5,),(F:4),( )Y(,D:3 ,WFT2),(:1121zyxQzyxfzyxMtuvf zu,vyxyxn 3.2 软件总体设计结果流程图此次程序采用 IPO 结构设计，首先通过读取 cmd 文件，得到相关输入参数：观测面位场数据文件名变量、延拓后位场数据文件名变量、延拓后准确位场数据文件名变量、扩边因子、延拓的高度（m/z 轴向下为正方向） ；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。下一步

7、，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来计算延拓因子并且将扩完边的数据进行快速二维傅里叶变换后在频率域与延拓因子相乘；最后进行快速二维傅里叶反变换并且去除扩边部分后输出。总体设计见表 1。44 测试结果4.1 测试参数（1）向上延拓原始场值数据保存在a20_mag.grd”中，向上延拓 3.3m，延拓后理论数据保存在“a53_mag.grd”中, 延拓后的数据保存在 output1.grd 中。网格大小为27*27（m） ，Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.（2）向下延拓原始场值数据保存在a53_mag.gr

8、d”中，向下延拓 3.3m，延拓后理论数据保存在“a20_mag.grd”中延拓后的数据保存在 output2.grd 中。网格大小为27*27（m ） ，Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有关参数保存在 filename.cmd 文件中，如下:A20_mag.grd A53_mag.grdoutput1.grdoutput2.grd13.3-3.3（1）输入有关参数（2）计算有关参数I：（3）从文件输入有关数据(1)进行扩边处理(2)Fourier 正变换(3)计算延拓因子(4)进行乘积运算P:(5)Fourier 反变换O： 输出计算结果表

9、1 总体设计 N-S 图54.2 测试结果-2501-5012-0512-350467589-2501-501-05128648图 1 低高度观测面(a20_mag)位场等值线图 图 2 高高度观测面(a53_mag)位场等值线图-2501-5012-50-8640185-105125-35046758图 3 高高度观测面(a53_mag)向下延拓位场等值线图 图 4 低高度观测面(a20_mag)向上延拓位场等值线图65 结论及建议由测试结果可知，向上延拓可以对浅部异常进行压制，且误差较小；向下延拓可突出浅部异常，且延拓结果误差较大。通过本次频率域位场处理与转换设计实验，我收获颇丰，同时也感

10、触很深！首先，刚开始我对空间域转频率域以及位场延拓的概念与处理是有不理解的，但是经过老师耐心细致的讲解和同学们之间的帮助，最终克服了这些困难。其次，加深了我对于位场延拓的作用和具体延拓步骤的理解。感谢老师孜孜不倦的教导与同学不厌其烦的讲解。谢谢！图 5 运行结果图7附录：源程序代码!*! 程序功能：实现频率域位场延拓! cmd 文件参数:! cmdfile:存放有关参数的文件名变量! filename_obser:低高度观测面位场数据文件! filename_output1：向上延拓后位场数据文件! filename_output2：向下延拓后位场数据文件! filename_obser2:高

11、高度观测面数据文件! factor_m：扩边因子! distance1:向上延拓的高度（m/z 轴向下为正方向）! distance2:向下延拓的高度（m/z 轴向下为正方向）! .grd 文件参数：! N_point,N_line：点数 (x 方向)、线数(y 方向)! x_min,x_max：x 的最小值和最大值! y_min,y_max：y 的最小值和最大值! Ur:初始观测面场值! 扩边参数：! m1,m2：x 方向实际数据起点和终点点号位置! 1,m：x 方向扩边后数据起点和终点点号位置! n1,n2：y 方向实际数据起点和终点点号位置! 1,n3：y 方向扩边后数据起点和终点点号位

12、置! 延拓参数：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! Fr：延拓观测面的信号的实部! Fi：延拓观测面的信号的虚部! U：延拓后准确场值 ! W(m,n)：径向圆频率! Y(m,n)：延拓因子! 误差参数：! error:延拓后场值与真实场值的均方误差!*8*Program plyytimplicit nonereal eigvalparameter(eigval=3.701411*1e5)character*(80)cmdfile,filename_obser,filename_output1,filename_output2,filename_obser2real

13、,allocatable:Ur(:,:),Ui(:,:),y(:,:),Fr(:,:),Fi(:,:),U(:,:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,n1,n2integer factor_mreal xmin,xmax,ymin,ymax,dx,dyreal distance1,distance2,errorcmdfile=filename.cmdcall read_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,2:正变换 )! 输出参数说

14、明：! Freal：信号的实部! Fimage：信号的虚部 (对于实信号而言，赋值为 0)! 对应频率分布说明：! 数据 Freal(m,n)和 Fimage(m,n)对应的频率分布位置为:! m 方向： 0,1,.,m/2-1,m/2,-(m/2-1),-1! n 方向： 0,1,.,n/2-1,n/2,-(n/2-1),-1!*SUBROUTINE FFT2(Freal,Fimage,m,n,nf)implicit noneINTEGER m,n,nfREAL Freal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE:Treal(:),Timage(:)

15、integer nmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOP DO i=1,m,1IF (n.ne.1) THENdo j=1,n,1Treal(j)=Freal(i,j)18Timage(j)=Fimage(i,j)end docall FFT(Treal,Timage,n,nf)do j=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)end doEND IFEND DODO j=1,n,1IF(m.n

16、e.1) THENdo i=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)end docall FFT(Treal,Timage,m,nf)do i=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)end doEND IFEND DOdeallocate(Treal,Timage,STAT=ierr)END SUBROUTINE FFT2!*! 子程序： FFT! 功能：复数组 1-D 快速 Fourier 变换! 输入参数说明：! Xreal(n): 输入数据实部 ! Ximage(n): 输入数据虚部 1

17、9! N: 点数（N 必须为 2 的幂次方）! NF: 正、反变换标志量 (1:反变换;2:正变换 )! 输出参数说明：! Xreal(n): 变换后频谱实部 ! Ximage(n): 变换后频谱虚部 ! 对应频率分布说明：! 数据 Xreal(n)和 Ximage(n)对应的频率分布位置为:! 0,1,.,n/2-1,n/2,-(n/2-1),-1!*SUBROUTINE FFT(Xreal,Ximage,n,nf)implicit noneINTEGER n,nfREAL Xreal(1:n),Ximage(1:n)integer nu,n2,nu1,k,k1,k1n2,l,i,ibitr

18、real f,p,arg,c,s,treal,timagenu=int(log(float(n)/0.693147+0.001)n2=n/2nu1=nu-1f=float(-1)*nf)k=0DO l=1,nu,1DO while (k.lt.n)do i=1,n2,1p=ibitr(k/2*nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*s20Xreal(k1n2)=X

19、real(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k) thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageend ifEND DOIF(nf

20、.ne.1) THENdo i=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)end doEND IFEND SUBROUTINE FFT21FUNCTION IBITR(J,NU)implicit noneinteger ibitr,j,nuinteger j1,itt,i,j2j1=jitt=0do i=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doEND FUNCTION IBITR!*! 子程序： cal_dxdy! 功能：计算 x,y 方向的点距! 输入参

21、数说明：! x_min,x_max：x 的最小值和最大值! y_min,y_max：y 的最小值和最大值! N_point,N_line：点数 (x 方向)、线数(y 方向)! 输出参数说明：! dx,dy：x,y 方向的点距!*subroutine cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicit nonereal xmin,xmax,ymin,ymaxinteger N_POINT,N_LINEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINE22end subrou

22、tine cal_dxdy!*! 子程序： WAVE2D! 功能：计算 2D 径向圆频率 W! 输入参数说明：! dx: x 方向点距 ! dy: y 方向线距! m: 点数（M 必须为 2 的幂次方）! n: 线数（N 必须为 2 的幂次方）! 输出参数说明： ! W(m,n): 径向圆频率 !*SUBROUTINE WAVE2D(m,n,dx,dy,W)implicit noneINTEGER m,nREAL dx,dyREAL W(1:m,1:n),U(1:m),V(1:n)real pi,delx,delyinteger midm,midn,i,j,xx,yypi=3.1415926m

23、idm=m/2+1midn=n/2+1delx=pi*2.0/float(m)/dxdely=pi*2.0/float(n)/dydo j=1,n,1yy=jif(yy.gt.midn) yy=yy-nv(j)=dely*(yy-1)end dodo i=1,m,1xx=i23if(xx.gt.midm) xx=xx-mu(i)=delx*(xx-1)end dodo j=1,n,1do i=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j)end doend doEND SUBROUTINE WAVE2D!*! 子程序： calculate_Y! 功能：计算延拓因子 Y

24、! 输入参数说明：! distance:延拓的高度（m/z 轴向下为正方向）! m,n： x,y 方向扩边后数据终点点号位置（起始点号为 1）! W：径向圆频率! 输出参数说明：! Y：延拓因子!*subroutine calculate_Y(m,n,W,y,distance)implicit noneinteger m,nreal pi,distancereal y(1:m,1:n),W(1:m,1:n)real i,jpi=3.1415926do i=1,m,1do j=1,n,1Y(I,J)=1*exp(-1*w(i,j)*distance)enddoenddo24end subrout

25、ine calculate_Y!*! 子程序： calculate_FmulY! 功能：将延拓因子与原信号做乘积! 输入参数说明：! m,n： x,y 方向扩边后数据终点点号位置（起始点号为 1）! Y：延拓因子! Ur：初始信号的实部! Ui：初始信号的虚部! 输出参数说明：! Fr：延拓信号的实部! Fi：延拓信号的虚部!*subroutine calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)implicit noneinteger m,nreal Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n)rea

26、l i,jdo i=1,mdo j=1,nFr(i,j)=Ur(i,j)*y(i,j)Fi(i,j)=Ui(i,j)*y(i,j)enddoenddoend subroutine calculate_FmulY!*25*! 子程序：output_grd! 功能：按照 grd 格式输出延拓后的场值! 输入参数说明：! N_point,N_line：点数 (x 方向)、线数(y 方向)! filename_output：输出文件的文件名! m1,m2：x 方向实际数据起点和终点点号位置! 1,m：x 方向扩边后数据起点和终点点号位置! n1,n2：y 方向实际数据起点和终点点号位置! 1,n3：y

27、 方向扩边后数据起点和终点点号位置! x_min,x_max,y_min,y_max：x/y 的最小值、最大值! A：输出数组!*subroutine output_grd(filename_output,N_point,N_line,m,n,m1,m2,n1,n2,eigval,A,Xmin,Xmax,Ymin,Ymax)character*(*)filename_outputinteger N_point,N_lineinteger m,n,m1,m2,n1,n2real Xmin,Xmax,Ymin,Ymaxreal A(1:m,1:n)real amin,amax,i,jDO j=n1

28、,n2,1do i=m1,m2,1if(ABS(A(i,j).lt.eigval) thenamin=MIN(amin,A(I,j) amax=MAX(amax,A(I,j)endifEND DOenddoOPEN(20,file=filename_output,status=unknown,form=formatted)26write (20,(A)DSAAwrite (20,*)m2-m1+1,n2-n1+1write (20,*)Xmin,Xmaxwrite (20,*)Ymin,Ymaxwrite(20,*)amin,amaxDo j=n1,n2,1write(20,*) (A(i,j

29、),i=m1,m2)End doClose(20)end subroutine output_grd!*! 子程序：calculate_error! 功能：计算均方误差! 输入参数说明：! distance:延拓的方向! m1,m2：x 方向实际数据起点和终点点号位置! 1,m：x 方向扩边后数据起点和终点点号位置! n1,n2：y 方向实际数据起点和终点点号位置! 1,n3：y 方向扩边后数据起点和终点点号位置! Fr,U：待求均方误差的两数组! 输出参数说明：! error:延拓后场值与真实场值的均方误差!*subroutine calculate_error(distance,Fr,U,

30、m1,m2,n1,n2,m,n,error)implicit noneinteger m1,m2,n1,n2,m,nreal Fr(1:m,1:n),U(1:m,1:n)real distancereal errorreal i,j27error=0.0do j=n1,n2do i=m1,m2error=error+(U(i,j)-Fr(i,j)*2enddoenddoerror=sqrt(error)if(distance0.0)thenwrite(*,*)请输入向上延拓的均方误差elsewrite(*,*)请输入向下延拓的均方误差end ifwrite(*,*)errorend subroutine