实验三：频率域位场处理与转换实验报告

1、 重磁资料处理与解释实验三频率域位场处理和转换实验 专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-12-21 基本原理1.1位场的方程由场论知识可知，位场方程分为两大类：有源的Possion方程，以及无源的Laplace方程。Laplace方程的第一边值问题通常为Dirichlet问题，第二边值问题通常称为Nueman问题。若P点在S平面内称为内部问题，反之称为外部问题。由唯一性定理可知，Dirichlet的内部和外部问题的解是唯一的，而Nueman内部问题的解不是唯一的，有一常数差，但其外部问题解是唯一的。外部问题的解的唯一性的原因：。无

2、源区域位场可以表示为： （1-1） （1-2） 1.2二维傅里叶变换及卷积性质： （1）傅里叶变换： （1-3） （1-4）（2） 卷积性质： （1-5） （1-6） 1.3频率域位场延拓原理当已知实测平面的异常时，换算场源以外的异常称之为延拓，分为向上延拓和向下延拓。半空间狄利克莱问题解析解： （1-5）其中：称为延拓因子，为计算面Z坐标 ,Z轴向下为正方向， 为计算面频率域位场， 为延拓面的频率域位场。 2 输入/输出数据格式设计2.1 输入数据格式设计观测面位场数据保存在filename_obser文件中，为.grd格式。计算延拓误差时的精确场值文件保存在filename_obser2中

3、，为.grd格式。2.2 输出数据格式设计实际计算得到的数据保存在filename_output文件中，为.grd格式。2.3 参数文件数据格式设计将所要读取的参数保存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。在该文件中保存的参数如下： filename_obser:观测面位场数据文件名变量 filename_output：延拓后位场数据文件名变量 filename_obser2:延拓后准确位场数据文件名变量 factor_m：扩边因子 distance: 延拓的高度（m/z轴向下为正方向）3 总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，

4、得到相关输入参数：观测面位场数据文件名变量、延拓后位场数据文件名变量、延拓后准确位场数据文件名变量、扩边因子、延拓的高度（m/z轴向下为正方向）；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来计算延拓因子并且将扩完边的数据进行快速二维傅里叶变换后在频率域与延拓因子相乘；最后进行快速二维傅里叶反变换并且去除扩边部分后输出。总体设计见表1。图1 总体设计N-S图输入参数：filename_obser:观测面位场数据文件名变量、filename_output：延拓后位场数据文件名变量

5、、filename_obser2:延拓后准确位场数据文件名变量、factor_m：扩边因子、distance: 延拓的高度（m/z轴向下为正方向）确定扩边网格的大小m*n(m,n均为2的幂次方)从观测面位场数据文件中读取数据二维扩边子程序并为信号虚部赋值扩边后的数据进行二维快速傅里叶正变换计算延拓因子将扩完边并进行快速二维傅里叶正变换后的复数数组在频率域与延拓因子相乘对上步结果进行快速二维傅里叶反变换去除扩边部分后输出结果4 测试结果4.1 测试参数 （1）向上延拓 原始场值数据保存在a20_mag.grd”中，向上延拓3.3m，延拓后理论数据保存在“a53_mag.grd”中。网格大小为27

6、*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有关参数保存在filename.cmd文件中，如下:A20_mag.grd，存放初始观测面的数据文件A53_mag.grd，存放延拓后理论数据文件Output.grd，存放延拓观测面的数据文件1，扩边因子3.3 延拓的高度（m/z轴向下为正方向）（2）向下延拓 原始场值数据保存在a53_mag.grd”中，向下延拓3.3m，延拓后理论数据保存在“a20_mag.grd”中。网格大小为27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有关参数保存在fi

7、lename.cmd文件中，如下:A53_mag.grd，存放初始观测面的数据文件A20_mag.grd，存放延拓后理论数据文件Output.grd，存放延拓观测面的数据文件1，扩边因子-3.3，延拓的高度（m/z轴向下为正方向）4.2 测试结果 对比图2，图3，可以看出，向上延拓的结果与理论磁异常基本一致，均方误差仅为152.1404；对比图1和图3可以看出，向上延拓对地下3个浅部异常体的的磁异常进行了压制；而对比图1，图4，向下延拓的磁异常平面图不太光滑，与理论磁异常平面图相比有较大误差，为2488.719；对比图2和图4可以看出，向下延拓突出了地下三个浅部异常体的磁异常特征。5 结论及建

8、议由测试结果可知，向上延拓可以对浅部异常进行压制，且误差较小；向下延拓可突出浅部异常，且延拓结果误差较大。经此次试验，对位场延拓基本步骤有了进一步了解，完成程序编写且测试正确。附录：程序源代码!*! 程序功能：实现频率域位场延拓! cmd文件参数:! cmdfile:存放有关参数的文件名变量! filename_obser:观测面位场数据文件! filename_output：延拓后位场数据文件! filename_obser2:延拓后准确位场数据文件!factor_m：扩边因子! distance:延拓的高度（m/z轴向下为正方向）! .grd文件参数：!N_point,N_line：点数(

9、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方向扩边后数据起点和终点点号位置!延拓参数：! Ur:初始观测面信号的实部!Ui：初始观测面信号的虚部!Fr：延拓观测面的信号的实部!Fi：延拓观测面的信号的虚部!U：延拓后准确场值 !W(m,n)：径向圆频率!Y(m,n)：延拓因子! 误差参数：! error:延拓后场值与真实场值的均方误差!*pr

10、ogram prolongationparameter(eigval=3.701411*1e5)character*(80)cmdfile,filename_obser,filename_output,filename_obser2real,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 distance,error!real

11、*8 eigvalcmdfile='filename.cmd'call read_cmd(cmdfile,filename_obser,filename_output,factor_m,filename_obser2,distance)call read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax) call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,

12、1:n),Fi(1:m,1:n),U(1:m,1:n),W(1:m,1:n)call input_grd(Ur,filename_obser,m1,m2,n1,n2,m,n)call input_grd(U,filename_obser2,m1,m2,n1,n2,m,n)call expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)call FFT2(Ur,Ui,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,W)call calculate_Y(m,n,W

13、,y,distance)call calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)call FFT2(Fr,Fi,m,n,1)call output_grd(filename_output,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax,Ymin,Ymax)call calculate_error(Fr,U,m1,m2,n1,n2,m,n,error)deallocate(Ur,Ui,y,Fr,Fi,u,w)end program!*! 子程序：read_cmd!功能：读取参数文件!输入参数说明：!cmdfile：参数文

14、件名!输出参数说明：!filename_obser：存放初始观测面的数据文件!filename_output：存放延拓观测面的数据文件!filename_obser2：存放参照观测面的数据文件! distance:延拓的高度（m/z轴向下为正方向）!factor_m：扩边因子!*subroutine read_cmd(cmdfile,filename_obser,filename_output,factor_m,filename_obser2,distance)implicit nonecharacter*(*)cmdfile,filename_obser,filename_output,fi

15、lename_obser2integer factor_mreal distance!real*8 eigvalcharacter*80 stropen(10,file=cmdfile,status='old') read(10,*)filename_obser read(10,*)filename_obser2 read(10,*)filename_output read(10,*)factor_m read(10,*)distanceclose(10)end subroutine read_cmd!*! 子程序： read_grd!功能：从原始观测.grd文件中读取相关参数

16、!输入参数说明：!filename_obser：输入文件名!输出参数说明：!N_point,N_line：点数、线数!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!*subroutine read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicit nonecharacter*(*)filename_obserinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=filename_obser,status='

17、old') Read(10,*) Read(10,*)N_line,N_point Read(10,*)Xmin,Xmax Read(10,*)Ymin,YmaxClose(10)end subroutine read_grd!*! 子程序：calculate_mn!功能：确定扩边数据点号位置!输入参数说明：!factor_m: 扩边比例因子（>1.0）!a,b：点数、线数!输出参数说明：!m1,m2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!n:扩边后数据终点位置点号(起点位置点号

18、为1)!*subroutine calculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicit noneinteger a,b,m,n,m1,m2,n1,n2integer mtemp,mu,nuinteger factor_mmtemp=a DO WHILE (mod(mtemp,2).eq.0).and.(mtemp.ne.0) mtemp=mtemp/2End doIF (mtemp.eq.1) THEN m=a*2ELSE mu=int(log(float(a)/0.693147+factor_m) m=2*muEND IFm1=1+(m-a)/2

19、m2=m1+a-1write(*,*)m,apausemtemp=b DO WHILE (mod(mtemp,2).eq.0).and.(mtemp.ne.0) mtemp=mtemp/2End doIF (mtemp.eq.1) THEN n=b*2ELSE nu=int(log(float(b)/0.693147+factor_m) n=2*nuEND IFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,m,npauseend subroutine calculate_mn!*! 子程序：INPUT_GRD!功能：读取grd文件中的数据!输入参数说明：

20、!filename_obser：输入文件名!m1,m2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!n:扩边后数据终点位置点号(起点位置点号为1)!输出参数说明：!A：存放输出数据的二维数组名!*SUBROUTINE INPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)character*(*)filename_obserinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,k

21、do j=1,n,1 do i=1,m,1 A(i,j)=3.701411*1e10enddo enddo Open(20,file=filename_obser,status='old') read(20,*) read(20,*) read(20,*)xmin,xmax read(20,*)ymin,ymax read(20,*) read(20,*) (A(i,j),i=m1,m2),j=n1,n2) Close(20)END SUBROUTINE INPUT_GRD!*! 子程序：expand_cos_2D!功能：二维扩边子程序并为信号虚部赋值!输入参数说明：!m1,m

22、2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部!Ui：初始观测面信号的虚部!输出参数说明：! Ur:初始观测面信号的实部!Ui：初始观测面信号的虚部!*subroutine expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicit noneinteger m,n,m1,m2,n1,n2real Ur(1:m,1:n),Ui(1:m,1:n)real,allocatable:u(:),r(:)

23、integer j,i,kallocate(u(1:m)do j=n1,n2,1 do i=1,m,1 u(i)=Ur(i,j) enddo call expand_cos_1d(1,m1,m2,m,u(1) do i=1,m,1 Ur(i,j)=u(i) enddoenddodeallocate(u)allocate(r(1:n)do i=1,m,1 do j=1,n,1 r(j)=Ur(i,j) enddo call expand_cos_1d(1,n1,n2,n,r(1) do j=1,n,1 Ur(i,j)=r(j) enddoenddodeallocate(r)do i=1,m do

24、 j=1,n Ui(i,j)=0 enddoenddoend subroutine expand_cos_2D!*! 子程序：expand_cos_1d!功能：一维扩边子程序!输入参数说明：!n0,n3：扩边后数据起点位置和终点位置!n1,n2：实际数据起点位置和终点位置!feild(i),(i=n1,n1+1,.,n2)：实际数据!输出参数说明：!field(i),(i=n0,.,n1-1)：扩边后左边的数据!field(i),(i=n2+1,.,n3)：扩边后右边的数据!*Subroutine expand_cos_1d(n0,n1,n2,n3,Field) Real Field(n0:n

25、3) pi=3.141592654 Field (n0)=(Field (n1)+Field (n2)/2.0 Field (n3)=Field (n0) i1=n0 i2=n1 DO i=i1,i2-1,1 Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1) End do i1=n2 i2=n3 DO i=i1+1,i2,1 Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1) End do End subroutine expand_

26、cos_1d!*! 功能：FFT2!功能：复数组2-D快速Fourier变换!输入参数说明：!m0,m3：x方向的起点和终点!n0,n3：y方向的起点和终点!field：输入信号(需要赋值给Freal，实部)!m,n：x,y方向扩边后数据终点点号位置（起始点号为1）!NF: 正、反变换标志量(1:反变换;2:正变换)!输出参数说明：!Freal：信号的实部!Fimage：信号的虚部(对于实信号而言，赋值为0)!对应频率分布说明：!数据Freal(m,n)和Fimage(m,n)对应的频率分布位置为:!m 方向： 0,1,.,m/2-1,m/2,-(m/2-1),.,-1!n 方向： 0,1,.

27、,n/2-1,n/2,-(n/2-1),.,-1!* SUBROUTINE FFT2(Freal,Fimage,m,n,nf) implicit none INTEGER m,n,nf REAL Freal(1:m,1:n),Fimage(1:m,1:n) real,ALLOCATABLE:Treal(:),Timage(:) integer nmmax,ierr,i,j nmmax=max(m,n) allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOP DO i=1,m,1 IF (n.ne.1) THEN

28、do j=1,n,1 Treal(j)=Freal(i,j) Timage(j)=Fimage(i,j) end do call FFT(Treal,Timage,n,nf) do j=1,n,1 Freal(i,j)=Treal(j) Fimage(i,j)=Timage(j) end do END IF END DO DO j=1,n,1 IF(m.ne.1) THEN do i=1,m,1 Treal(i)=Freal(i,j) Timage(i)=Fimage(i,j) end do call FFT(Treal,Timage,m,nf) do i=1,m,1 Freal(i,j)=T

29、real(i) Fimage(i,j)=Timage(i) end do END IF END DO deallocate(Treal,Timage,STAT=ierr) END SUBROUTINE FFT2!*! 子程序：FFT!功能：复数组1-D快速Fourier变换!输入参数说明：!Xreal(n): 输入数据实部 !Ximage(n): 输入数据虚部 !N: 点数（N 必须为 2 的幂次方）!NF: 正、反变换标志量(1:反变换;2:正变换)!输出参数说明：!Xreal(n): 变换后频谱实部 !Ximage(n): 变换后频谱虚部 !对应频率分布说明：!数据Xreal(n)和Xim

30、age(n)对应的频率分布位置为:!0,1,.,n/2-1,n/2,-(n/2-1),.,-1!*SUBROUTINE FFT(Xreal,Ximage,n,nf) implicit none INTEGER n,nf REAL Xreal(1:n),Ximage(1:n) integer nu,n2,nu1,k,k1,k1n2,l,i,ibitr real f,p,arg,c,s,treal,timage nu=int(log(float(n)/0.693147+0.001)n2=n/2nu1=nu-1f=float(-1)*nf)k=0 DO l=1,nu,1 DO while (k.lt

31、.n) do i=1,n2,1 p=ibitr(k/2*nu1,nu) arg=6.2831853*p*f/float(n) c=cos(arg) s=sin(arg) k1=k+1 k1n2=k1+n2 treal=Xreal(k1n2)*c+Ximage(k1n2)*s timage=Ximage(k1n2)*c-Xreal(k1n2)*s Xreal(k1n2)=Xreal(k1)-treal Ximage(k1n2)=Ximage(k1)-timage Xreal(k1)=Xreal(k1)+treal Ximage(k1)=Ximage(k1)+timage k=k+1 end do

32、 k=k+n2 END DO k=0 nu1=nu1-1 n2=n2/2 END DO DO k=1,n,1 i=ibitr(k-1,nu)+1 if(i.gt.k) then treal=Xreal(k) timage=Ximage(k) Xreal(k)=Xreal(i) Ximage(k)=Ximage(i) Xreal(i)=treal Ximage(i)=timage end if END DO IF(nf.ne.1) THEN do i=1,n,1 Xreal(i)=Xreal(i)/float(n) Ximage(i)=Ximage(i)/float(n) end do END

33、IF END SUBROUTINE FFT FUNCTION IBITR(J,NU) implicit none integer ibitr,j,nu integer j1,itt,i,j2j1=jitt=0 do i=1,nu,1 j2=j1/2 itt=itt*2+(j1-2*j2) ibitr=itt j1=j2 end do END FUNCTION IBITR!*! 子程序： cal_dxdy!功能：计算x,y方向的点距!输入参数说明：!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!N_point,N_line：点数(x方向)、线数(y方向)

34、!输出参数说明：!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_LINEend subroutine cal_dxdy!*! 子程序：WAVE2D!功能：计算2D径向圆频率W!输入参数说明：!dx: x 方向点距 !dy: y 方向线距!m: 点数（M 必须为 2 的幂次方）

35、!n: 线数（N 必须为 2 的幂次方）!输出参数说明： !W(m,n): 径向圆频率 !*SUBROUTINE WAVE2D(m,n,dx,dy,W) implicit none INTEGER m,n REAL dx,dy REAL W(1:m,1:n),U(1:m),V(1:n) real pi,delx,dely integer midm,midn,i,j,xx,yypi=3.1415926midm=m/2+1midn=n/2+1delx=pi*2.0/float(m)/dxdely=pi*2.0/float(n)/dydo j=1,n,1 yy=j if(yy.gt.midn) yy

36、=yy-n v(j)=dely*(yy-1)end dodo i=1,m,1 xx=i if(xx.gt.midm) xx=xx-m u(i)=delx*(xx-1)end dodo j=1,n,1 do i=1,m,1 w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j) end doend do END SUBROUTINE WAVE2D!*! 子程序：calculate_Y!功能：计算延拓因子Y!输入参数说明：! distance:延拓的高度（m/z轴向下为正方向）!m,n：x,y方向扩边后数据终点点号位置（起始点号为1）!W：径向圆频率!输出参数说明：!Y：延拓因子!*sub

37、routine 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,1 do j=1,n,1 Y(I,J)=1*exp(-1*w(i,j)*distance) enddoenddoend subroutine calculate_Y!*! 子程序：calculate_FmulY!功能：将延拓因子与原信号做乘积!输入参数说明：!m,n：x,y方向扩边后数据终点点号位置（起始点号为1）!Y：延拓因子!Ur

38、：初始信号的实部!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)real i,jdo i=1,m do j=1,n Fr(i,j)=Ur(i,j)*y(i,j) Fi(i,j)=Ui(i,j)*y(i,j) enddoenddoend subroutine calculate_FmulY!*! 子程序：output_grd!功能：按照grd格式输出延拓后的场值!输入参数说明：!N_point