實(shí)驗(yàn)四位場邊緣識別程序設(shè)計(jì)實(shí)驗(yàn)

1、重磁資料處理與解釋實(shí)驗(yàn)四位場邊緣識別程序設(shè)計(jì)實(shí)驗(yàn)專業(yè)名稱：地球物理學(xué)學(xué)生姓名：學(xué)生學(xué)號：指導(dǎo)老師：王萬銀、紀(jì)新林、紀(jì)曉琳、邱世燦提交日期：2016-1-31、基本原理地質(zhì)目標(biāo)體邊緣時(shí)指斷裂構(gòu)造線、不同地質(zhì)體的邊界線，實(shí)際上是具有一定密度或磁性差異的地質(zhì)體的邊界線，在地質(zhì)體的邊緣附近，重、磁異常變化率較大，故所有的邊緣識別方法均利用這一特點(diǎn)進(jìn)行設(shè)計(jì)?，F(xiàn)在有重、磁位場邊緣識別方法分為數(shù)理統(tǒng)計(jì)、數(shù)值計(jì)算以及其他三大類。數(shù)值類邊緣識別方法均利用極大值位置或零值位置確定地質(zhì)體的邊緣位置，其依據(jù)的理論基礎(chǔ)是二度體鉛垂臺階模型的重力異常特征。在該模型邊緣處重力異?？偹綄?dǎo)數(shù)和解析信號振幅達(dá)到極大值、垂向?qū)?shù)

2、達(dá)到零值。故可以利用這些特征位置來確定二度體鉛垂臺階的邊緣位置，確定傾斜二度體、不規(guī)則二度體及三度體邊緣位置的理論均為二度體鉛垂臺階模型理論的推廣，但確定的邊緣位置和真實(shí)的位置有一定的偏差。該偏差隨著地質(zhì)體邊界形狀、埋深、水平尺寸及物性差異等的變化而變化。因此，邊緣識別結(jié)果是一種定性或半定量解釋結(jié)果，與定量解釋結(jié)果有一定區(qū)別，識別結(jié)果可作為邊緣位置定量反演的初值。(1)垂向?qū)?shù)：垂向?qū)?shù)方法利用零值位置確定地質(zhì)體的邊緣位置，重力異?？梢灾苯邮褂?，對磁力異常必須轉(zhuǎn)化為磁源重力異?；蚧瘶O磁力異常才可以使用 ( 1.1) （2）解析信號振幅：解析信號振幅也是利用極大值位置來確定地質(zhì)體的邊緣位置適用于

3、重、磁力異常 (1.2) （3）總水平導(dǎo)數(shù)（THDR） (1.3)2、輸入/輸出數(shù)據(jù)格式設(shè)計(jì)依據(jù)上述原理，現(xiàn)在對上述各種邊緣識別方法進(jìn)行程序設(shè)計(jì)。2.1 輸入數(shù)據(jù)格式設(shè)計(jì)本次實(shí)驗(yàn)給了正演的重力異常數(shù)據(jù)，為.GRD格式，均為實(shí)型變量。例如：DSAA 201 201 -1000.000000 1000.000000 -1000.000000 1000.000000 5.549671E-01 23.539846 5.549671E-01 5.634658E-01 5.721339E-01 5.808522E-01 5.897312E-01 5.987253E-01 6.078691E-01 6.17

4、1604E-012.2 輸出數(shù)據(jù)格式設(shè)計(jì)計(jì)算結(jié)果輸出數(shù)據(jù)格式與輸入格式對應(yīng)，格式為.GRD格式，均為實(shí)型變量。例如：DSAA 201 201 -1000.000 1000.000 -1000.000 1000.000 -0.1465084 0.3190881 -3.6523044E-02 -3.3485338E-02 -3.3061244E-02 -3.2748722E-02 -3.2688729E-02 -3.2654848E-02 -3.2723978E-02 -3.2787599E-02 -3.2927759E-02 -3.3138681E-022.3 參數(shù)文件數(shù)據(jù)格式設(shè)計(jì)將以上部分量保

5、存在一個(gè)文件中，該文件名變量為cmdfile，字符串變量，長度不超過80，全路徑名。在該文件中保存的參數(shù)如下：輸入數(shù)據(jù)文件名input_file，字符串變量，長度不超過80；輸出vdr數(shù)據(jù)文件名output_file_vdr，字符串變量，長度不超過80；輸出thdr數(shù)據(jù)文件名output_file_thdr，字符串變量，長度不超過80;輸出asm數(shù)據(jù)文件名output_file_asm，字符串變量，長度不超過80factor_m：擴(kuò)邊比例因子，實(shí)型變量（1）；3. 總體設(shè)計(jì)此次程序采用IPO結(jié)構(gòu)設(shè)計(jì)，首先通過讀取cmd文件，得到相關(guān)輸入?yún)?shù)：輸入數(shù)據(jù)文件名gravity.grd、輸出vdr文件

6、名field_vrd.grd、輸出thdr文件名field_thdr.grd、輸出asm文件名field_asm.grd、擴(kuò)邊比例因子factor_m；然后確定確定擴(kuò)邊網(wǎng)格的大小，擴(kuò)邊數(shù)據(jù)點(diǎn)號位置；再從觀測面位場數(shù)據(jù)文件中讀取數(shù)據(jù)。下一步，進(jìn)行二維余弦擴(kuò)邊，將擴(kuò)完邊的數(shù)據(jù)進(jìn)行快速二維傅里葉變換，轉(zhuǎn)換到頻率域；接下來在頻率域求出在x,y,z方向的導(dǎo)數(shù)并反變換；最后求出VDR、THDR、ASM數(shù)據(jù)。最后去除擴(kuò)邊部分后輸出。總體設(shè)計(jì)見表1。輸入?yún)?shù)：輸入數(shù)據(jù)文件名gravity.grd、輸出vdr文件名gravity_vrd.grd、輸出thdr文件名gravity_thdr.grd、輸出asm文件

7、名gravity_asm.grd 、擴(kuò)邊比例因子factor_m。 確定擴(kuò)邊網(wǎng)格的大小m*n(m,n均為2的冪次方)從輸入數(shù)據(jù)文件名中讀取數(shù)據(jù)對原始數(shù)據(jù)進(jìn)行二維余弦擴(kuò)邊對擴(kuò)邊后的數(shù)據(jù)進(jìn)行快速二維傅里葉正變換將傅里葉變換后的數(shù)據(jù)與導(dǎo)數(shù)因子相乘求出重力數(shù)據(jù)在x,y,z方向的導(dǎo)數(shù)對導(dǎo)數(shù)進(jìn)行快速二維傅里葉反變換分別求出VDR、THDR、ASM值。去除擴(kuò)邊部分后輸出結(jié)果圖4.1 總體設(shè)計(jì)N-S圖4.測試結(jié)果圖4.2 重力異常原始圖圖4.3 垂向?qū)?shù)（VDR）圖4.2 重力異常原始圖圖4.5 總水平導(dǎo)數(shù)（THDR）圖4.4 解析信號振幅（ASM）圖4.3 ASM（解析信號振幅）分析：由圖4.3可看出，VD

8、R方法的零值線可較準(zhǔn)確識別模型體的邊界；由圖4.4可看出，ASM的極大值點(diǎn)邊界可大致識別模型體邊界，但精度不是很高。對比圖4.2到4.5可以看出，THDR方法對模型邊界的識別效果是最好的。5 結(jié)論及建議 經(jīng)測試，VDR與THDR對模型體的邊界位置識別效果較好，而ASM對模型體邊界識別效果較差。三種方法中，THDR效果最好。附錄：邊緣識別程序源代碼!*! 程序功能：實(shí)現(xiàn)頻率域位場導(dǎo)數(shù)運(yùn)算進(jìn)行邊緣識別! cmd文件參數(shù):! cmdfile:存放有關(guān)參數(shù)的文件名變量! input_file:觀測面位場數(shù)據(jù)文件! output_file_vdr：場值的水平導(dǎo)數(shù)數(shù)據(jù)文件! output_file_thd

9、r:場值垂向?qū)?shù)數(shù)據(jù)文件! output_file_asm:場值的解析信號振幅數(shù)據(jù)文件!factor_m：擴(kuò)邊因子! .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)號位置!求導(dǎo)參數(shù)：! field_re:初始觀測面信號的實(shí)部! field_im：初始觀測面信號的虛部

10、! Px_re:x方向?qū)?shù)信號的實(shí)部! Px_im:x方向?qū)?shù)信號的虛部! Py_re:y方向?qū)?shù)信號的實(shí)部! Py_im:y方向?qū)?shù)信號的虛部! Pz_re:z方向?qū)?shù)信號的實(shí)部! Pz_im:z方向?qū)?shù)信號的虛部 !W(m,n)：徑向圓頻率! field_vdr:對場值作水平導(dǎo)數(shù)的結(jié)果! field_thdr:對場值作垂向?qū)?shù)的結(jié)果! field_asm:場值的解析信號振幅的結(jié)果!*program deviationparameter(eigval=3.701411*1e5)character*(80)cmdfilecharacter*80 input_file,output_file_v

11、dr,output_file_thdr,output_file_asmreal,allocatable: field_re(:,:),field_im(:,:)real,allocatable: Px_re(:,:),Px_im(:,:),Py_re(:,:),Py_im(:,:),Pz_re(:,:),Pz_im(:,:)real,allocatable: field_vdr(:,:),field_thdr(:,:),field_asm(:,:)real,allocatable: U(:),V(:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,

12、n1,n2real factor_mreal xmin,xmax,ymin,ymax,dx,dycmdfile=deviation.cmdcall read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)call read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax) call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(field_re(1:m

13、,1:n),field_im(1:m,1:n)allocate(Px_re(1:m,1:n),Px_im(1:m,1:n),Py_re(1:m,1:n),Py_im(1:m,1:n),Pz_re(1:m,1:n),Pz_im(1:m,1:n)allocate(field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n)allocate(U(1:m),V(1:n),W(1:m,1:n)call input_grd(field_re,input_file,m1,m2,n1,n2,m,n)call expand_cos_2D(m1,m2,m,n1

14、,n2,n,field_re,field_im)call FFT2(field_re,field_im,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,U,V,W)call factor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)call FFT2(px_re,px_im,m,n,1)call FFT2(py_re,py_im,m,n,1)call FFT2(pz_re,pz_im,m,n,1

15、)call deviration(m,n,px_re,py_re,pz_re,field_vdr,field_thdr,field_asm)call OUTPUT_GRD(field_vdr,output_file_vdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_thdr,output_file_thdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_asm,output_file_asm,m1,m2,n1,n2,m,n,

16、eigval,xmin,xmax,ymin,ymax)deallocate(field_re,field_im,px_re,px_im,py_re,py_im,pz_re,pz_im,field_vdr,field_thdr,field_asm,u,v,w)end program!*! 子程序：read_cmd!功能：讀取參數(shù)文件!輸入?yún)?shù)說明：!cmdfile：參數(shù)文件名!輸出參數(shù)說明：! input_file:觀測面位場數(shù)據(jù)文件! output_file_vdr：對場值作水平導(dǎo)數(shù)處理后的數(shù)據(jù)文件! output_file_thdr:對場值作垂向?qū)?shù)處理后的數(shù)據(jù)文件! output_file

17、_asm:對場值作總導(dǎo)數(shù)處理后的數(shù)據(jù)文件!factor_m：擴(kuò)邊因子!*Subroutine read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)implicit nonecharacter*80 strcharacter*(*)cmdfilecharacter*(*)input_file,output_file_vdr,output_file_thdr,output_file_asmreal factor_mopen(10,file=cmdfile,status=old)

18、read(10,*)str,input_file read(10,*)str,output_file_vdr read(10,*)str,output_file_thdr read(10,*)str,output_file_asm read(10,*)str,factor_mclose(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_ma

19、x：y的最小值和最大值!*subroutine read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicit nonecharacter*(*)input_fileinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=input_file,status=old) Read(10,*) Read(10,*)N_line,N_point Read(10,*)Xmin,Xmax Read(10,*)Ymin,YmaxClose(10)end subroutine read

20、_grd!*! 子程序：calculate_mn!功能：確定擴(kuò)邊數(shù)據(jù)點(diǎn)號位置!輸入?yún)?shù)說明：!factor_m: 擴(kuò)邊比例因子（1.0）!a,b：點(diǎn)數(shù)、線數(shù)!輸出參數(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)!*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,nu

21、real 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)/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

22、 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文件中的數(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

23、)!輸出參數(shù)說明：!A：存放輸出數(shù)據(jù)的二維數(shù)組名!*SUBROUTINE INPUT_GRD(A,input_file,m1,m2,n1,n2,m,n)character*(*)input_fileinteger 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=input_file,status=old) read(20,*) read(20,*) read(20,*)xmin,x

24、max 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!功能：二維擴(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:初始觀測面

25、信號的實(shí)部!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(:)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) enddoenddodea

26、llocate(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 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

27、(i),(i=n1,n1+1,.,n2)：實(shí)際數(shù)據(jù)!輸出參數(shù)說明：!field(i),(i=n0,.,n1-1)：擴(kuò)邊后左邊的數(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)/

28、(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_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: 正、反變

29、換標(biāo)志量(1:反變換;2:正變換)!輸出參數(shù)說明：!Freal：信號的實(shí)部!Fimage：信號的虛部(對于實(shí)信號而言，賦值為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,ALLOCATA

30、BLE: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 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) e

31、nd 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)=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ù)

32、說明：!Xreal(n): 輸入數(shù)據(jù)實(shí)部 !Ximage(n): 輸入數(shù)據(jù)虛部 !N: 點(diǎn)數(shù)（N 必須為 2 的冪次方）!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) inte

33、ger 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.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 tima

34、ge=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 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) Xim

35、age(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 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*j

36、2) ibitr=itt j1=j2 end do END FUNCTION IBITR!*! 子程序： cal_dxdy!功能：計(jì)算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_L

37、INEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEend subroutine cal_dxdy!*! 子程序：WAVE2D!功能：計(jì)算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,U,V,W) implicit none INTEGER m,n REAL dx,dy REAL W(1:m,1:n),U(1:m),V

38、(1:n) real pi,delx,dely integer midm,midn,i,j,xx,yymidm=m/2+1midn=n/2+1delx=float(m)/dxdely=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 SUBR

39、OUTINE WAVE2D!* ! 子程序：factor! 功能：計(jì)算x,y,z方向?qū)?shù)的實(shí)部和虛部! 輸入?yún)?shù)說明：!m: 點(diǎn)數(shù)（M 必須為 2 的冪次方）!n: 線數(shù)（N 必須為 2 的冪次方）! field_re:初始觀測面信號的實(shí)部! field_im：初始觀測面信號的虛部 !W(m,n)：徑向圓頻率! 輸出參數(shù)說明：! Px_re:x方向?qū)?shù)信號的實(shí)部! Px_im:x方向?qū)?shù)信號的虛部! Py_re:y方向?qū)?shù)信號的實(shí)部! Py_im:y方向?qū)?shù)信號的虛部! Pz_re:z方向?qū)?shù)信號的實(shí)部! Pz_im:z方向?qū)?shù)信號的虛部 !* subroutine factor(m,n,fi

40、eld_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im) implicit none integer m,n real field_re(1:m,1:n),field_im(1:m,1:n) real px_re(1:m,1:n),px_im(1:m,1:n),py_re(1:m,1:n),py_im(1:m,1:n),pz_re(1:m,1:n),pz_im(1:m,1:n) real U(1:m),V(1:n),W(1:m,1:n) integer i,j real pi pi=3.1415926 do i=1,m,1 do j=

41、1,n,1 px_re(i,j)=field_im(i,j)*u(i)*(-1)*pi*2.0 px_im(i,j)=field_re(i,j)*u(i)*pi*2.0 py_re(i,j)=field_im(i,j)*v(j)*(-1)*pi*2.0 py_im(i,j)=field_re(i,j)*v(j)*pi*2.0 pz_re(i,j)=field_re(i,j)*w(i,j)*pi*2.0 pz_im(i,j)=field_im(i,j)*w(i,j)*pi*2.0 end do end doend subroutine factor!* ! 子程序：deviration! 功能：計(jì)算異常的水平導(dǎo)數(shù)、垂向?qū)?shù)及總導(dǎo)數(shù)! 輸入?yún)?shù)說明：!m: 點(diǎn)數(shù)（M 必須為 2 的冪次方）!n: 線數(shù)（N 必須為 2 的冪次方）! Px_re:x方向?qū)?shù)信號的實(shí)部! Px_im:x方向?qū)?shù)信號的虛部! Py_re:y方向?qū)?shù)信號的實(shí)部! Py_im:y方向?qū)?shù)信號的虛部! Pz_re:z方向?qū)?shù)信號的實(shí)部! Pz_im:z方向?qū)?shù)信號的虛部 ! 輸出參數(shù)說明：! fi