實(shí)驗(yàn)三：頻率域位場處理與轉(zhuǎn)換實(shí)驗(yàn)報告

1、 重磁資料處理與解釋實(shí)驗(yàn)三頻率域位場處理和轉(zhuǎn)換實(shí)驗(yàn) 專業(yè)名稱：地球物理學(xué)學(xué)生姓名：學(xué)生學(xué)號：指導(dǎo)老師：王萬銀、紀(jì)新林、紀(jì)曉琳、邱世燦提交日期：2016-12-21 基本原理1.1位場的方程由場論知識可知，位場方程分為兩大類：有源的Possion方程，以及無源的Laplace方程。Laplace方程的第一邊值問題通常為Dirichlet問題，第二邊值問題通常稱為Nueman問題。若P點(diǎn)在S平面內(nèi)稱為內(nèi)部問題，反之稱為外部問題。由唯一性定理可知，Dirichlet的內(nèi)部和外部問題的解是唯一的，而Nueman內(nèi)部問題的解不是唯一的，有一常數(shù)差，但其外部問題解是唯一的。外部問題的解的唯一性的原因：。無

2、源區(qū)域位場可以表示為： （1-1） （1-2） 1.2二維傅里葉變換及卷積性質(zhì)： （1）傅里葉變換： （1-3） （1-4）（2） 卷積性質(zhì)： （1-5） （1-6） 1.3頻率域位場延拓原理當(dāng)已知實(shí)測平面的異常時，換算場源以外的異常稱之為延拓，分為向上延拓和向下延拓。半空間狄利克萊問題解析解： （1-5）其中：稱為延拓因子，為計算面Z坐標(biāo) ,Z軸向下為正方向， 為計算面頻率域位場， 為延拓面的頻率域位場。 2 輸入/輸出數(shù)據(jù)格式設(shè)計2.1 輸入數(shù)據(jù)格式設(shè)計觀測面位場數(shù)據(jù)保存在filename_obser文件中，為.grd格式。計算延拓誤差時的精確場值文件保存在filename_obser2中

3、，為.grd格式。2.2 輸出數(shù)據(jù)格式設(shè)計實(shí)際計算得到的數(shù)據(jù)保存在filename_output文件中，為.grd格式。2.3 參數(shù)文件數(shù)據(jù)格式設(shè)計將所要讀取的參數(shù)保存在一個文件中，該文件名變量為cmdfile，字符串變量，長度不超過80，全路徑名。在該文件中保存的參數(shù)如下： filename_obser:觀測面位場數(shù)據(jù)文件名變量 filename_output：延拓后位場數(shù)據(jù)文件名變量 filename_obser2:延拓后準(zhǔn)確位場數(shù)據(jù)文件名變量 factor_m：擴(kuò)邊因子 distance: 延拓的高度（m/z軸向下為正方向）3 總體設(shè)計此次程序采用IPO結(jié)構(gòu)設(shè)計，首先通過讀取cmd文件，

4、得到相關(guān)輸入?yún)?shù)：觀測面位場數(shù)據(jù)文件名變量、延拓后位場數(shù)據(jù)文件名變量、延拓后準(zhǔn)確位場數(shù)據(jù)文件名變量、擴(kuò)邊因子、延拓的高度（m/z軸向下為正方向）；然后確定確定擴(kuò)邊網(wǎng)格的大小，擴(kuò)邊數(shù)據(jù)點(diǎn)號位置；再從觀測面位場數(shù)據(jù)文件中讀取數(shù)據(jù)。下一步，進(jìn)行二維余弦擴(kuò)邊，將擴(kuò)完邊的數(shù)據(jù)進(jìn)行快速二維傅里葉變換，轉(zhuǎn)換到頻率域；接下來計算延拓因子并且將擴(kuò)完邊的數(shù)據(jù)進(jìn)行快速二維傅里葉變換后在頻率域與延拓因子相乘；最后進(jìn)行快速二維傅里葉反變換并且去除擴(kuò)邊部分后輸出。總體設(shè)計見表1。圖1 總體設(shè)計N-S圖輸入?yún)?shù)：filename_obser:觀測面位場數(shù)據(jù)文件名變量、filename_output：延拓后位場數(shù)據(jù)文件名變量

5、、filename_obser2:延拓后準(zhǔn)確位場數(shù)據(jù)文件名變量、factor_m：擴(kuò)邊因子、distance: 延拓的高度（m/z軸向下為正方向）確定擴(kuò)邊網(wǎng)格的大小m*n(m,n均為2的冪次方)從觀測面位場數(shù)據(jù)文件中讀取數(shù)據(jù)二維擴(kuò)邊子程序并為信號虛部賦值擴(kuò)邊后的數(shù)據(jù)進(jìn)行二維快速傅里葉正變換計算延拓因子將擴(kuò)完邊并進(jìn)行快速二維傅里葉正變換后的復(fù)數(shù)數(shù)組在頻率域與延拓因子相乘對上步結(jié)果進(jìn)行快速二維傅里葉反變換去除擴(kuò)邊部分后輸出結(jié)果4 測試結(jié)果4.1 測試參數(shù) （1）向上延拓 原始場值數(shù)據(jù)保存在a20_mag.grd”中，向上延拓3.3m，延拓后理論數(shù)據(jù)保存在“a53_mag.grd”中。網(wǎng)格大小為27

6、*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有關(guān)參數(shù)保存在filename.cmd文件中，如下:A20_mag.grd，存放初始觀測面的數(shù)據(jù)文件A53_mag.grd，存放延拓后理論數(shù)據(jù)文件Output.grd，存放延拓觀測面的數(shù)據(jù)文件1，擴(kuò)邊因子3.3 延拓的高度（m/z軸向下為正方向）（2）向下延拓 原始場值數(shù)據(jù)保存在a53_mag.grd”中，向下延拓3.3m，延拓后理論數(shù)據(jù)保存在“a20_mag.grd”中。網(wǎng)格大小為27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有關(guān)參數(shù)保存在fi

7、lename.cmd文件中，如下:A53_mag.grd，存放初始觀測面的數(shù)據(jù)文件A20_mag.grd，存放延拓后理論數(shù)據(jù)文件Output.grd，存放延拓觀測面的數(shù)據(jù)文件1，擴(kuò)邊因子-3.3，延拓的高度（m/z軸向下為正方向）4.2 測試結(jié)果圖1 z=-2.0km平面理論磁異常圖圖2 z=-5.3km平面理論磁異常圖圖3 -5.3km平面延拓磁異常圖（均方誤差152.1404）圖4 -2.0km平面延拓磁異常圖（均方誤差2488.719） 對比圖2，圖3，可以看出，向上延拓的結(jié)果與理論磁異?；疽恢?，均方誤差僅為152.1404；對比圖1和圖3可以看出，向上延拓對地下3個淺部異常體的的磁異

8、常進(jìn)行了壓制；而對比圖1，圖4，向下延拓的磁異常平面圖不太光滑，與理論磁異常平面圖相比有較大誤差，為2488.719；對比圖2和圖4可以看出，向下延拓突出了地下三個淺部異常體的磁異常特征。5 結(jié)論及建議由測試結(jié)果可知，向上延拓可以對淺部異常進(jìn)行壓制，且誤差較?。幌蛳卵油乜赏怀鰷\部異常，且延拓結(jié)果誤差較大。經(jīng)此次試驗(yàn)，對位場延拓基本步驟有了進(jìn)一步了解，完成程序編寫且測試正確。附錄：程序源代碼!*! 程序功能：實(shí)現(xiàn)頻率域位場延拓! cmd文件參數(shù):! cmdfile:存放有關(guān)參數(shù)的文件名變量! filename_obser:觀測面位場數(shù)據(jù)文件! filename_output：延拓后位場數(shù)據(jù)文件!

9、 filename_obser2:延拓后準(zhǔn)確位場數(shù)據(jù)文件!factor_m：擴(kuò)邊因子! distance:延拓的高度（m/z軸向下為正方向）! .grd文件參數(shù)：!N_point,N_line：點(diǎn)數(shù)(x方向)、線數(shù)(y方向)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值! Ur:初始觀測面場值!擴(kuò)邊參數(shù)：!m1,m2：x方向?qū)嶋H數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號位置!1,m：x方向擴(kuò)邊后數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號位置!n1,n2：y方向?qū)嶋H數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號位置!1,n3：y方向擴(kuò)邊后數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號位置!延拓參數(shù)：! Ur:初始觀測面信號的實(shí)部!Ui：初始觀測面信號的虛部

10、!Fr：延拓觀測面的信號的實(shí)部!Fi：延拓觀測面的信號的虛部!U：延拓后準(zhǔn)確場值 !W(m,n)：徑向圓頻率!Y(m,n)：延拓因子! 誤差參數(shù)：! error:延拓后場值與真實(shí)場值的均方誤差!*program 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

11、,N_lineinteger m,n,m1,m2,n1,n2integer factor_mreal xmin,xmax,ymin,ymax,dx,dyreal distance,error!real*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 calc

12、ulate_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,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)C

13、ALL 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,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

14、,n2,m,n,error)deallocate(Ur,Ui,y,Fr,Fi,u,w)end program!*! 子程序：read_cmd!功能：讀取參數(shù)文件!輸入?yún)?shù)說明：!cmdfile：參數(shù)文件名!輸出參數(shù)說明：!filename_obser：存放初始觀測面的數(shù)據(jù)文件!filename_output：存放延拓觀測面的數(shù)據(jù)文件!filename_obser2：存放參照觀測面的數(shù)據(jù)文件! distance:延拓的高度（m/z軸向下為正方向）!factor_m：擴(kuò)邊因子!*subroutine read_cmd(cmdfile,filename_obser,filename_output,f

15、actor_m,filename_obser2,distance)implicit nonecharacter*(*)cmdfile,filename_obser,filename_output,filename_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,

16、*)factor_m read(10,*)distanceclose(10)end subroutine read_cmd!*! 子程序： read_grd!功能：從原始觀測.grd文件中讀取相關(guān)參數(shù)!輸入?yún)?shù)說明：!filename_obser：輸入文件名!輸出參數(shù)說明：!N_point,N_line：點(diǎn)數(shù)、線數(shù)!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*(*)

17、filename_obserinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=filename_obser,status='old') Read(10,*) Read(10,*)N_line,N_point Read(10,*)Xmin,Xmax Read(10,*)Ymin,YmaxClose(10)end subroutine read_grd!*! 子程序：calculate_mn!功能：確定擴(kuò)邊數(shù)據(jù)點(diǎn)號位置!輸入?yún)?shù)說明：!factor_m: 擴(kuò)邊比例因子（>1.0）!a,b：點(diǎn)數(shù)、線數(shù)!輸出參數(shù)

18、說明：!m1,m2: x方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為1)!n1,n2: y方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為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/2

19、End doIF (mtemp.eq.1) THEN m=a*2ELSE mu=int(log(float(a)/0.693147+factor_m) m=2*muEND IFm1=1+(m-a)/2m2=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+

20、b-1write(*,*)m1,m2,n1,n2,m,npauseend subroutine calculate_mn!*! 子程序：INPUT_GRD!功能：讀取grd文件中的數(shù)據(jù)!輸入?yún)?shù)說明：!filename_obser：輸入文件名!m1,m2: x方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為1)!n1,n2: y方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為1)!輸出參數(shù)說明：!A：存放輸出數(shù)據(jù)的二維數(shù)組名!*SUBROUTINE INPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)

21、character*(*)filename_obserinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,k 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

22、,m2),j=n1,n2) Close(20)END SUBROUTINE INPUT_GRD!*! 子程序：expand_cos_2D!功能：二維擴(kuò)邊子程序并為信號虛部賦值!輸入?yún)?shù)說明：!m1,m2: x方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為1)!n1,n2: y方向?qū)嶋H數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)位置點(diǎn)號(起點(diǎn)位置點(diǎn)號為1)! Ur:初始觀測面信號的實(shí)部!Ui：初始觀測面信號的虛部!輸出參數(shù)說明：! Ur:初始觀測面信號的實(shí)部!Ui：初始觀測面信號的虛部!*subroutine expand_cos_2D(m1,m2,m,n1,n

23、2,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(:)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

24、) 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 j=1,n Ui(i,j)=0 enddoenddoend subroutine expand_cos_2D!*! 子程序：expand_cos_1d!功能：一維擴(kuò)邊子程序!輸入?yún)?shù)說明：!n0,n3：擴(kuò)邊后數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置!n1,n2：實(shí)際數(shù)據(jù)起點(diǎn)位置和終點(diǎn)位置!feild(i),(i=n1,n1+1,.,n2)：實(shí)際數(shù)據(jù)!輸出參數(shù)說明：!field(i),(i=n0,.,n1-1)：擴(kuò)

25、邊后左邊的數(shù)據(jù)!field(i),(i=n2+1,.,n3)：擴(kuò)邊后右邊的數(shù)據(jù)!*Subroutine expand_cos_1d(n0,n1,n2,n3,Field) Real Field(n0:n3) 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,

26、i2,1 Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1) End do End subroutine expand_cos_1d!*! 功能：FFT2!功能：復(fù)數(shù)組2-D快速Fourier變換!輸入?yún)?shù)說明：!m0,m3：x方向的起點(diǎn)和終點(diǎn)!n0,n3：y方向的起點(diǎn)和終點(diǎn)!field：輸入信號(需要賦值給Freal，實(shí)部)!m,n：x,y方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號位置（起始點(diǎn)號為1）!NF: 正、反變換標(biāo)志量(1:反變換;2:正變換)!輸出參數(shù)說明：!Freal：信號的實(shí)部!Fimage：信號的虛部(對于實(shí)信號而言

27、，賦值為0)!對應(yīng)頻率分布說明：!數(shù)據(jù)Freal(m,n)和Fimage(m,n)對應(yīng)的頻率分布位置為:!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 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,

28、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 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

29、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)=Treal(i) Fimage(i,j)=Timage(i) end do END IF END DO deallocate(Treal,Timage,STAT=ierr) END SUBROUTINE FFT2!*! 子程序：FFT!功能：復(fù)數(shù)組1-D快速Fourier變換!輸入?yún)?shù)說明：!Xreal(n): 輸入數(shù)據(jù)實(shí)部 !Ximage(n): 輸入數(shù)據(jù)虛部 !N: 點(diǎn)數(shù)（N 必須為 2 的冪次

30、方）!NF: 正、反變換標(biāo)志量(1:反變換;2:正變換)!輸出參數(shù)說明：!Xreal(n): 變換后頻譜實(shí)部 !Ximage(n): 變換后頻譜虛部 !對應(yīng)頻率分布說明：!數(shù)據(jù)Xreal(n)和Ximage(n)對應(yīng)的頻率分布位置為:!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,t

31、image 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.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)-trea

32、l Ximage(k1n2)=Ximage(k1)-timage Xreal(k1)=Xreal(k1)+treal Ximage(k1)=Ximage(k1)+timage k=k+1 end do 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 EN

33、D 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 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_d

34、xdy!功能：計算x,y方向的點(diǎn)距!輸入?yún)?shù)說明：!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!N_point,N_line：點(diǎn)數(shù)(x方向)、線數(shù)(y方向)!輸出參數(shù)說明：!dx,dy：x,y方向的點(diǎn)距!*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_LINEen

35、d subroutine cal_dxdy!*! 子程序：WAVE2D!功能：計算2D徑向圓頻率W!輸入?yún)?shù)說明：!dx: x 方向點(diǎn)距 !dy: y 方向線距!m: 點(diǎn)數(shù)（M 必須為 2 的冪次方）!n: 線數(shù)（N 必須為 2 的冪次方）!輸出參數(shù)說明： !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.141592

36、6midm=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=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!功能：

37、計算延拓因子Y!輸入?yún)?shù)說明：! distance:延拓的高度（m/z軸向下為正方向）!m,n：x,y方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號位置（起始點(diǎn)號為1）!W：徑向圓頻率!輸出參數(shù)說明：!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,1 do j=1,n,1 Y(I,J)=1*exp(-1*w(i,j)*distance) enddoenddoend subrouti

38、ne calculate_Y!*! 子程序：calculate_FmulY!功能：將延拓因子與原信號做乘積!輸入?yún)?shù)說明：!m,n：x,y方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號位置（起始點(diǎn)號為1）!Y：延拓因子!Ur：初始信號的實(shí)部!Ui：初始信號的虛部!輸出參數(shù)說明：!Fr：延拓信號的實(shí)部!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格式輸出延拓后的場值!輸入?yún)?shù)說