频率域位场处理与转换(1225-)

1、重磁实验报告 实验名称：频率域位场处理与转换 学生姓名：刘国强专业名称：勘查技术与工程学生学号：201226020126指导老师：巨鹏 王万银目录一实验内容与要求2（1）实验内容：2（2）要求：2二基本原理2（1）空间域位场延拓：2（2）频率域位场延拓：3三.输入输出数据格式设计31.输入数据：32.重要变量名3四总体设计4五测试结果5六结论与建议6一实验内容与要求（1） 实验内容：对一个网格化数据进行向上延拓，得到延拓后的结果，并计算延拓后的均方误差。（2） 要求：进行软件设计、代码编写、软件测试、编写实验报告。二基本原理分别介绍空间域和频率域位场延拓的计算公式。（1）空间域位场延拓：延拓分

2、为向上延拓和向下延拓。向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。对于二度体（令z向下为正）： U(x,z)=-z-+U,0-x2+z2d （z<0） 其中U（，0）为剖面上各点的实测值。Z为延拓的高度，即转换后的平面和观测平面的距离，由于z坐标向下为正，因此z<0。空间域的延拓实际是积分的计算。 向下延拓是由实测磁场向磁源方向延拓。 （2）频率域位场延拓：设场源位于z=H平面以下（H>0），则磁场为在z=H平面以上对x、y、z连续的函数，若z=0观测平面的磁场T（x，y，0）为已知，则由外部的狄利克莱问题： T（x，y，z）=-z2-+-+T（，0）(x-)2

3、+(y-)2+z23/2dd 对其进行变量x，y进行傅里叶变换成 STu，v，z=-T（x，y，z）e-2i(ux+vy)dxdy 因此：STu，v，z= STu，v，0e2(u2+v2)1/2z 在频率域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。三.输入输出数据格式设计1.输入数据：（1）建立文件：wave.par（2）观测面高度之差：height=3.3m（3）低高度网格化数据：A20_mag.grd（4）特征值：eigval=1.701411*1038（6）输出文件名：expound_a20_2D.grd2.重要变量名（1）输入文件：input_filename(2) 输出文件：

4、out_filename(3) 点数：mpoint(4)线数：nline(5)异常值：filed四总体设计结束调用子程序，输出文件调用子程序，计算延拓因子调用子程序，对扩边数据进行快速傅里叶变换调用子程序，进行扩边调用程序，读入A20_mag.grd场值调用程序，计算扩边点位调用子程序，读入A20_mag.grd数据开始调用子程序，读入wave.par文件 五测试结果A53_mag.filenameA20_mag_filenameExpound_A20_mag_filenameme分析：将底高度向上延拓3.3m后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了。六结论与建议结

5、论：向上延拓可以压制浅部的重磁异常特征，突出深部异常。向下延拓是由实测磁源异常得到。建议：本次实验程序编写比较困难，以后要提高编程能力。最后得到了扩边后的数据，分析能力比较差，分析结果不够准确，多看资料，加强这方面能力。七源程序及代码PROGRAM wavenumber_continuationINTEGER mpoint,line,m0,m1,m2,m3,n0,n1,n2,n3,m,nREAL height,eigval,Xmin,Xmax,Ymin,Ymax,dx,dyCHARACTER*80 input_filename,output_filename,cmd_filename,Expo

6、und_a20_2D_filenameREAL,ALLOCATABLE:Field(:,:),Freal(:,:),Fimage(:,:)REAL,ALLOCATABLE:U(:),v(:),w(:,:),Fcount_Real(:,:),Fcount_image(:,:)cmd_filename='wave.par'CALL Read_CMD(cmd_filename,height,input_filename,Expound_a20_2D_filename,output_filename,eigval)CALL Read_mn(input_filename,mpoint,l

7、ine,Xmin,Xmax,Ymin,Ymax)CALL Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx)CALL Calculate_m1m2_dx(line,n0,n1,n2,n3,n,Ymin,Ymax,dy)ALLOCATE(U(m0:m3),V(n0:n3),W(m0:m3,n0:n3),Fcount_Real(m0:m3,n0:n3),Fcount_image(m0:m3,n0:n3)ALLOCATE(Field(m0:m3,n0:n3),Freal(m0:m3,n0:n3),Fimage(m0:m3,n0:n3) CALL

8、Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval)CALL Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line)CALL Putout_GRD_Expound_A20_2D(Field,Expound_a20_2D_filename,m0,m3,n0,n3,m,n,eigval,Xmin,Xmax,Ymin,Ymax) Freal=Field Fimage=0.CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,2) !2说明是正变换CALL WAVE2D

9、_sub(U,V,W,m,n,dx,dy)CALL Factor_conti(height,m0,m3,n0,n3,U,V,W,Fcount_Real,Fcount_image) !计算延拓因子CALL Multi_sub(Freal,Fimage,Fcount_Real,Fcount_image,m0,m3,n0,n3)CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,1)CALL Putout_GRD(Freal,output_filename,m0,m1,m2,m3,n0,n1,n2,n3,eigval,Xmin,Xmax,Ymin,Ymax)Dealloca

10、te(Field,Freal,Fimage,U,V,W,Fcount_Real,Fcount_image)END PROGRAM!搜寻未打开的文件的通道号!SUBROUTINE SEARCH_UNIT(NUNIT_START,NUNIT_IN)! logical unit_open! integer nunit_in,nunit_start! nunit_in=nunit_start! unit_open=.TRUE. ! do while (unit_open)! nunit_in=nunit_in+1! inquire(unit=nunit_in,opened=unit_open) ! e

11、nd do!END SUBROUTINE!读入数据文件!SUBROUTINE Read_CMD(cmd_filename,height,input_filename,Expound_a20_2D_filename,output_filename,eigval) real height,eigval character*80 input_filename,output_filename,cmd_filename,Expound_a20_2D_filename! call search_unit(10,nunit_in) OPEN(11,file=cmd_filename,status='

12、old') read(11,*) height read(11,*) input_filename read(11,*) Expound_a20_2D_filename read(11,*) output_filename read(11,*) eigval close(11)END SUBROUTINE!读入源文件中的值!SUBROUTINE Read_mn(input_filename,mpoint,line,Xmin,Xmax,Ymin,Ymax) integer mpoint,line real Xmin,Xmax,Ymin,Ymax character*80 input_fi

13、lename ! call search_unit(10,nunit_in) OPEN(12,file=input_filename) READ(12,*) READ(12,*) mpoint,line READ(12,*) Xmin,Xmax READ(12,*) Ymin,Ymax ! read(12,*) close(12)END SUBROUTINE!计算扩边点位!SUBROUTINE Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx) integer mpoint,mtemp,m0,m1,m2,m3,m real Xmin,Xma

14、x,dx,factor_m,mu m0=1 factor_m=1.5 dx=(Xmax-Xmin)/(mpoint-1) mtemp=mpoint do while(mod(mtemp,2).eq.0).and.(mtemp.ne.0) mtemp=mtemp/2 end do if(mtemp.eq.1) then m=mpoint*2 else mu=int(log(float(mpoint)/0.693147+factor_m) m=2*mu end if m1=1+(m-mpoint)/2 m2=m1+mpoint-1 m3=mEND SUBROUTINE!输入源文件的重磁异常值!SU

15、BROUTINE Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval) character*(*) input_filename integer m0,m1,m2,m3,n0,n1,n2,n3 real eigval real Field(m0:m3,n0:n3) Field=eigval ! call search_unit(10,nunit_in) open(13,file=input_filename,form='formatted') read(13,*) read(13,*) read(13,*)

16、 Xmin,Xmax read(13,*) Ymin,Ymax read(13,*) read(13,*) (Field(i,j),i=m1,m2),j=n1,n2) Close(13)END SUBROUTINE!*2-D扩边*SUBROUTINE Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line) integer m0,m1,m2,m3,n0,n1,n2,n3 real Field(m0:m3,n0:n3) real:pi=3.141592654 do i=n1,n2 Field(m0,i)=(Field(m1,i)+Field(m2,i)/2 F

17、ield(m3,i)=Field(m0,i) end do do j=n1,n2 do i=m0+1,m1-1 Field(i,j)=Field(m0,j)+cos(pi/2.0*(m1-i)/(m1-m0)*(Field(m1,j)-Field(m0,j) end do do i=m2+1,m3-1 Field(i,j)=Field(m2,j)+cos(pi/2.0*(m3-i)/(m3-m2)*(Field(m3,j)-Field(m2,j) end do end do do i=m0,m3 Field(i,n0)=(Field(i,n1)+Field(i,n2)/2 Field(i,n3

18、)=Field(i,n0) end do do i=m0,m3 do j=n0+1,n1-1 Field(i,j)=Field(i,n0)+cos(pi/2.0*(n1-j)/(n1-n0)*(Field(i,n1)-Field(i,n0) end do do j=n2+1,n3-1 Field(i,j)=Field(i,n2)+cos(pi/2.0*(n3-j)/(n3-n2)*(Field(i,n3)-Field(i,n2) end do end doEND SUBROUTINE!*输出扩边Expound_a20_2D.grd*SUBROUTINE Putout_GRD_Expound_A

19、20_2D(Field,Expound_a20_2D_filename,m0,m3,n0,n3,m,n,eigval,Xmin,Xmax,Ymin,Ymax) real Field(m0:m3,n0:n3) character*(*) Expound_a20_2D_filename real Gmin,Gmax Gmin=Huge(Gmin) Gmax=-Huge(Gmax) do j=n0,n3 do i=m0,m3 Gmin=min(Gmin,Field(i,j) Gmax=max(Gmax,Field(i,j) end do ! call search_unit(10,nunit_in)

20、 end do open(14,file=Expound_a20_2D_filename,status='unknown') write(14,'(A)') 'DSAA' write(14,*) m3-m0+1,n3-n0+1 write(14,*) -(m3-m0),m3-m0 write(14,*) -(n3-n0),n3-n0 write(14,*) Gmin,Gmax do j=n0,n3 write(14,*) (Field(i,j),i=m0,m3) end do close(14)END SUBROUTINE!*傅里叶正变换*SUB

21、ROUTINE FFT2(DREAL,DIMAG,M,N,NF) real dreal(m,n),dimag(m,n) real,ALLOCATABLE:treal(:),timag(:) nmmax=max(m,n) allocate(treal(nmmax),timag(nmmax),STAT=ierr) DO i=1,m IF (n.ne.1) THEN do j=1,n treal(j)=dreal(i,j) timag(j)=dimag(i,j) end do call fft(treal,timag,n,nf) do j=1,n dreal(i,j)=treal(j) dimag(

22、i,j)=timag(j) end do END IF END DO DO j=1,n IF(m.ne.1) THEN do i=1,m treal(i)=dreal(i,j) timag(i)=dimag(i,j) end do call fft(treal,timag,m,nf) do i=1,m dreal(i,j)=treal(i) dimag(i,j)=timag(i) end do END IF END DO deallocate(treal,timag,STAT=ierr)END SUBROUTINESUBROUTINE FFT(XREAL,XIMAG,N,NF) real xr

23、eal(n),ximag(n) nu=int(log(float(n)/0.693147+0.001) n2=n/2 nu1=nu-1 DO l=1,nu DO while (k.lt.n) do i=1,n2 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+ximag(k1n2)*s timag=ximag(k1n2)*c-xreal(k1n2)*s xreal(k1n2)=xreal(k1)-treal ximag(k1n2)

24、=ximag(k1)-timag xreal(k1)=xreal(k1)+treal ximag(k1)=ximag(k1)+timag f=float(-1)*nf) k=0 k=k+1 end do k=k+n2 END DO k=0 nu1=nu1-1 n2=n2/2 END DO DO k=1,n i=ibitr(k-1,nu)+1 if(i.gt.k) then treal=xreal(k) timag=ximag(k) xreal(k)=xreal(i) ximag(k)=ximag(i) xreal(i)=treal ximag(i)=timag end if END DO IF

25、(nf.ne.1) THEN do i=1,n xreal(i)=xreal(i)/float(n) ximag(i)=ximag(i)/float(n) end do END IFEND SUBROUTINEFUNCTION IBITR(J,NU) j1=j itt=0 do i=1,nu j2=j1/2 itt=itt*2+(j1-2*j2) ibitr=itt j1=j2 end doend function!*计算延拓因子*v SUBROUTINE WAVE2D_sub(U,V,W,m,n,dx,dy)subroutine WAVE2D_sub(U,V,W,m,n,dx,dy) REA

26、L W(1:m,1:n),U(1:m),V(1:n) CALL Wave1D_sub(n,dy,V) CALL Wave1D_sub(m,dx,U) DO j=1,n DO i=1,m W(i,j)=SQRT(U(i)*U(i)+V(j)*V(j) END DO END DOEND subroutine WAVE2D_subSubroutine Wave1D_sub(N,dy,V) REAL dy INTEGER N REAL V(0:N-1) REAL deltf deltf=(2.0*3.141592654)/(N*dy) Do j=0,N/2,1 V(j)=j*deltf End do

27、Do j=N/2+1,N-1,1 V(j)=(j-n)*deltf End doEND subroutine Wave1D_sub!计算延拓因子!SUBROUTINE Factor_conti(height,m0,m3,n0,n3,U,V,W,Fcount_Real,Fcount_image) integer m0,m3,n0,n3 real height real U(m0:m3),V(n0:n3),W(m0:m3,n0:n3),Fcount_Real(m0:m3,n0:n3),Fcount_image(m0:m3,n0:n3) DO j=n0,n3 DO i=m0,m3 Fcount_Real(i,j)=EXP(-W(i,j)*height) END DO END DOEND SUBROUTINESUBROUTINE Multi_sub