頻率域位場處理與轉(zhuǎn)換(1225-)

1、重磁實驗報告 實驗名稱：頻率域位場處理與轉(zhuǎn)換 學生姓名：劉國強專業(yè)名稱：勘查技術(shù)與工程學生學號：201226020126指導老師：巨鵬 王萬銀目錄一實驗內(nèi)容與要求2（1）實驗內(nèi)容：2（2）要求：2二基本原理2（1）空間域位場延拓：2（2）頻率域位場延拓：3三.輸入輸出數(shù)據(jù)格式設(shè)計31.輸入數(shù)據(jù)：32.重要變量名3四總體設(shè)計4五測試結(jié)果5六結(jié)論與建議6一實驗內(nèi)容與要求（1） 實驗內(nèi)容：對一個網(wǎng)格化數(shù)據(jù)進行向上延拓，得到延拓后的結(jié)果，并計算延拓后的均方誤差。（2） 要求：進行軟件設(shè)計、代碼編寫、軟件測試、編寫實驗報告。二基本原理分別介紹空間域和頻率域位場延拓的計算公式。（1）空間域位場延拓：延拓分

2、為向上延拓和向下延拓。向上延拓是只將實測平面上的數(shù)據(jù)換到平面之上的平面上的計算。對于二度體（令z向下為正）： U(x,z)=-z-+U,0-x2+z2d （z<0） 其中U（，0）為剖面上各點的實測值。Z為延拓的高度，即轉(zhuǎn)換后的平面和觀測平面的距離，由于z坐標向下為正，因此z<0?？臻g域的延拓實際是積分的計算。 向下延拓是由實測磁場向磁源方向延拓。 （2）頻率域位場延拓：設(shè)場源位于z=H平面以下（H>0），則磁場為在z=H平面以上對x、y、z連續(xù)的函數(shù)，若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 在頻率域中，延拓變成了對觀測數(shù)據(jù)的傅里葉變換乘以延拓因子。三.輸入輸出數(shù)據(jù)格式設(shè)計1.輸入數(shù)據(jù)：（1）建立文件：wave.par（2）觀測面高度之差：height=3.3m（3）低高度網(wǎng)格化數(shù)據(jù)：A20_mag.grd（4）特征值：eigval=1.701411*1038（6）輸出文件名：expound_a20_2D.grd2.重要變量名（1）輸入文件：input_filename(2) 輸出文件：

4、out_filename(3) 點數(shù)：mpoint(4)線數(shù)：nline(5)異常值：filed四總體設(shè)計結(jié)束調(diào)用子程序，輸出文件調(diào)用子程序，計算延拓因子調(diào)用子程序，對擴邊數(shù)據(jù)進行快速傅里葉變換調(diào)用子程序，進行擴邊調(diào)用程序，讀入A20_mag.grd場值調(diào)用程序，計算擴邊點位調(diào)用子程序，讀入A20_mag.grd數(shù)據(jù)開始調(diào)用子程序，讀入wave.par文件 五測試結(jié)果A53_mag.filenameA20_mag_filenameExpound_A20_mag_filenameme分析：將底高度向上延拓3.3m后結(jié)果與高高度原場值圖形近似，但是場值變小了，且異常部位相對不明顯了。六結(jié)論與建議結(jié)

5、論：向上延拓可以壓制淺部的重磁異常特征，突出深部異常。向下延拓是由實測磁源異常得到。建議：本次實驗程序編寫比較困難，以后要提高編程能力。最后得到了擴邊后的數(shù)據(jù)，分析能力比較差，分析結(jié)果不夠準確，多看資料，加強這方面能力。七源程序及代碼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!讀入數(shù)據(jù)文件!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