位場邊緣識別程序設計實驗

1、精選優(yōu)質文檔-傾情為你奉上精選優(yōu)質文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質文檔-傾情為你奉上專心-專注-專業(yè)重磁資料處理與解釋實驗四位場邊緣識別程序設計實驗專業(yè)名稱：地球物理學學生姓名：學生學號：指導老師：王萬銀、紀新林、紀曉琳、邱世燦提交日期：2016-1-31、基本原理地質目標體邊緣時指斷裂構造線、不同地質體的邊界線，實際上是具有一定密度或磁性差異的地質體的邊界線，在地質體的邊緣附近，重、磁異常變化率較大，故所有的邊緣識別方法均利用這一特點進行設計。現(xiàn)在有重、磁位場邊緣識別方法分為數(shù)理統(tǒng)計、數(shù)值計算以及其他三大類。數(shù)值類邊緣識別方法均利用極大值位置或零值位置確定地質體的邊

2、緣位置，其依據(jù)的理論基礎是二度體鉛垂臺階模型的重力異常特征。在該模型邊緣處重力異?？偹綄?shù)和解析信號振幅達到極大值、垂向導數(shù)達到零值。故可以利用這些特征位置來確定二度體鉛垂臺階的邊緣位置，確定傾斜二度體、不規(guī)則二度體及三度體邊緣位置的理論均為二度體鉛垂臺階模型理論的推廣，但確定的邊緣位置和真實的位置有一定的偏差。該偏差隨著地質體邊界形狀、埋深、水平尺寸及物性差異等的變化而變化。因此，邊緣識別結果是一種定性或半定量解釋結果，與定量解釋結果有一定區(qū)別，識別結果可作為邊緣位置定量反演的初值。(1)垂向導數(shù)：垂向導數(shù)方法利用零值位置確定地質體的邊緣位置，重力異?？梢灾苯邮褂?，對磁力異常必須轉化為磁源

3、重力異?；蚧瘶O磁力異常才可以使用 ( 1.1) （2）解析信號振幅：解析信號振幅也是利用極大值位置來確定地質體的邊緣位置適用于重、磁力異常 (1.2) （3）總水平導數(shù)（THDR） (1.3)2、輸入/輸出數(shù)據(jù)格式設計依據(jù)上述原理，現(xiàn)在對上述各種邊緣識別方法進行程序設計。2.1 輸入數(shù)據(jù)格式設計本次實驗給了正演的重力異常數(shù)據(jù)，為.GRD格式，均為實型變量。例如：DSAA 201 201 -1000. 1000. -1000. 1000. 5.E-01 23. 5.E-01 5.E-01 5.E-01 5.E-01 5.E-01 5.E-01 6.E-01 6.E-012.2 輸出數(shù)據(jù)格式設計計

4、算結果輸出數(shù)據(jù)格式與輸入格式對應，格式為.GRD格式，均為實型變量。例如：DSAA 201 201 -1000.000 1000.000 -1000.000 1000.000 -0. 0. -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-02 -3.E-022.3 參數(shù)文件數(shù)據(jù)格式設計將以上部分量保存在一個文件中，該文件名變量為cmdfile，字符串變量，長度不超過80，全路徑名。在該文件中保存的參數(shù)如下：輸入數(shù)據(jù)文件名input_file，字符串變量，長度不超過80；輸出vdr數(shù)據(jù)文件名output_f

5、ile_vdr，字符串變量，長度不超過80；輸出thdr數(shù)據(jù)文件名output_file_thdr，字符串變量，長度不超過80;輸出asm數(shù)據(jù)文件名output_file_asm，字符串變量，長度不超過80factor_m：擴邊比例因子，實型變量（1）；總體設計此次程序采用IPO結構設計，首先通過讀取cmd文件，得到相關輸入?yún)?shù)：輸入數(shù)據(jù)文件名gravity.grd、輸出vdr文件名field_vrd.grd、輸出thdr文件名field_thdr.grd、輸出asm文件名field_asm.grd、擴邊比例因子factor_m；然后確定確定擴邊網(wǎng)格的大小，擴邊數(shù)據(jù)點號位置；再從觀測面位場數(shù)據(jù)

6、文件中讀取數(shù)據(jù)。下一步，進行二維余弦擴邊，將擴完邊的數(shù)據(jù)進行快速二維傅里葉變換，轉換到頻率域；接下來在頻率域求出在x,y,z方向的導數(shù)并反變換；最后求出VDR、THDR、ASM數(shù)據(jù)。最后去除擴邊部分后輸出?？傮w設計見表1。輸入?yún)?shù)：輸入數(shù)據(jù)文件名gravity.grd、輸出vdr文件名gravity_vrd.grd、輸出thdr文件名gravity_thdr.grd、輸出asm文件名gravity_asm.grd 、擴邊比例因子factor_m。 確定擴邊網(wǎng)格的大小m*n(m,n均為2的冪次方)從輸入數(shù)據(jù)文件名中讀取數(shù)據(jù)對原始數(shù)據(jù)進行二維余弦擴邊對擴邊后的數(shù)據(jù)進行快速二維傅里葉正變換將傅里葉變

7、換后的數(shù)據(jù)與導數(shù)因子相乘求出重力數(shù)據(jù)在x,y,z方向的導數(shù)對導數(shù)進行快速二維傅里葉反變換分別求出VDR、THDR、ASM值。去除擴邊部分后輸出結果圖4.1 總體設計N-S圖4.測試結果圖4.2 重力異常原始圖圖4.3 垂向導數(shù)（VDR）圖4.2 重力異常原始圖圖4.5 總水平導數(shù)（THDR）圖4.4 解析信號振幅（ASM）圖4.3 ASM（解析信號振幅）分析：由圖4.3可看出，VDR方法的零值線可較準確識別模型體的邊界；由圖4.4可看出，ASM的極大值點邊界可大致識別模型體邊界，但精度不是很高。對比圖4.2到4.5可以看出，THDR方法對模型邊界的識別效果是最好的。5 結論及建議 經(jīng)測試，VD

8、R與THDR對模型體的邊界位置識別效果較好，而ASM對模型體邊界識別效果較差。三種方法中，THDR效果最好。附錄：邊緣識別程序源代碼!*! 程序功能：實現(xiàn)頻率域位場導數(shù)運算進行邊緣識別! cmd文件參數(shù):! cmdfile:存放有關參數(shù)的文件名變量! input_file:觀測面位場數(shù)據(jù)文件! output_file_vdr：場值的水平導數(shù)數(shù)據(jù)文件! output_file_thdr:場值垂向導數(shù)數(shù)據(jù)文件! output_file_asm:場值的解析信號振幅數(shù)據(jù)文件!factor_m：擴邊因子! .grd文件參數(shù)：!N_point,N_line：點數(shù)(x方向)、線數(shù)(y方向)!x_min,x_

9、max：x的最小值和最大值!y_min,y_max：y的最小值和最大值! Ur:初始觀測面場值!擴邊參數(shù)：!m1,m2：x方向實際數(shù)據(jù)起點和終點點號位置!1,m：x方向擴邊后數(shù)據(jù)起點和終點點號位置!n1,n2：y方向實際數(shù)據(jù)起點和終點點號位置!1,n3：y方向擴邊后數(shù)據(jù)起點和終點點號位置!求導參數(shù)：! field_re:初始觀測面信號的實部! field_im：初始觀測面信號的虛部! Px_re:x方向導數(shù)信號的實部! Px_im:x方向導數(shù)信號的虛部! Py_re:y方向導數(shù)信號的實部! Py_im:y方向導數(shù)信號的虛部! Pz_re:z方向導數(shù)信號的實部! Pz_im:z方向導數(shù)信號的虛部

10、 !W(m,n)：徑向圓頻率! field_vdr:對場值作水平導數(shù)的結果! field_thdr:對場值作垂向導數(shù)的結果! field_asm:場值的解析信號振幅的結果!*program deviationparameter(eigval=3.*1e5)character*(80)cmdfilecharacter*80 input_file,output_file_vdr,output_file_thdr,output_file_asmreal,allocatable: field_re(:,:),field_im(:,:)real,allocatable: Px_re(:,:),Px_im

11、(:,:),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,n1,n2real factor_mreal xmin,xmax,ymin,ymax,dx,dycmdfile=deviation.cmdcall read_cmd(cmdfile,factor_m,input_file,out

12、put_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,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

13、: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,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)c

14、all 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)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

15、,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,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,f

16、ield_asm,u,v,w)end program!*! 子程序：read_cmd!功能：讀取參數(shù)文件!輸入?yún)?shù)說明：!cmdfile：參數(shù)文件名!輸出參數(shù)說明：! input_file:觀測面位場數(shù)據(jù)文件! output_file_vdr：對場值作水平導數(shù)處理后的數(shù)據(jù)文件! output_file_thdr:對場值作垂向導數(shù)處理后的數(shù)據(jù)文件! output_file_asm:對場值作總導數(shù)處理后的數(shù)據(jù)文件!factor_m：擴邊因子!*Subroutine read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_th

17、dr,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) read(10,*)str,input_file read(10,*)str,output_file_vdr read(10,*)str,output_file_thdr read(10,*)str,output_file_as

18、m read(10,*)str,factor_mclose(10)end Subroutine read_cmd!*! 子程序： read_grd!功能：從原始觀測.grd文件中讀取相關參數(shù)!輸入?yún)?shù)說明：!filename_obser：輸入文件名!輸出參數(shù)說明：!N_point,N_line：點數(shù)、線數(shù)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!*subroutine read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicit nonecharacter*(*)input_fil

19、einteger 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_grd!*! 子程序：calculate_mn!功能：確定擴邊數(shù)據(jù)點號位置!輸入?yún)?shù)說明：!factor_m: 擴邊比例因子（1.0）!a,b：點數(shù)、線數(shù)!輸出參數(shù)說明：!m1,m2: x方向實際數(shù)據(jù)起點位置和終點位置點號!m

20、:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2: y方向實際數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為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,nureal factor_mmtemp=a DO WHILE (mod(mtemp,2).eq.0).and.(mtemp.ne.0) mtemp=mtemp/2End doIF (mtemp.eq.1) THEN m=a*2EL

21、SE mu=int(log(float(a)/0.+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.+factor_m) n=2*nuEND IFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,m,npauseend subroutin

22、e calculate_mn!*! 子程序：INPUT_GRD!功能：讀取grd文件中的數(shù)據(jù)!輸入?yún)?shù)說明：!filename_obser：輸入文件名!m1,m2: x方向實際數(shù)據(jù)起點位置和終點位置點號!m:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2: y方向實際數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!輸出參數(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

23、 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.*1e10enddo enddo Open(20,file=input_file,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!功能：二維

24、擴邊子程序并為信號虛部賦值!輸入?yún)?shù)說明：!m1,m2: x方向實際數(shù)據(jù)起點位置和終點位置點號!m:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2: y方向實際數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)! Ur:初始觀測面信號的實部!Ui：初始觀測面信號的虛部!輸出參數(shù)說明：! 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)

25、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) enddo call expand_cos_1d(1,n1,n2,n,r(1) do j=1,n,1 Ur(i,j)=r(j) enddoend

26、dodeallocate(r)do i=1,m do j=1,n Ui(i,j)=0 enddoenddoend subroutine expand_cos_2D!*! 子程序：expand_cos_1d!功能：一維擴邊子程序!輸入?yún)?shù)說明：!n0,n3：擴邊后數(shù)據(jù)起點位置和終點位置!n1,n2：實際數(shù)據(jù)起點位置和終點位置!feild(i),(i=n1,n1+1,.,n2)：實際數(shù)據(jù)!輸出參數(shù)說明：!field(i),(i=n0,.,n1-1)：擴邊后左邊的數(shù)據(jù)!field(i),(i=n2+1,.,n3)：擴邊后右邊的數(shù)據(jù)!*Subroutine expand_cos_1d(n0,n1,n2

27、,n3,Field) Real Field(n0:n3) pi=3. 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 s

28、ubroutine expand_cos_1d!*! 功能：FFT2!功能：復數(shù)組2-D快速Fourier變換!輸入?yún)?shù)說明：!m0,m3：x方向的起點和終點!n0,n3：y方向的起點和終點!field：輸入信號(需要賦值給Freal，實部)!m,n：x,y方向擴邊后數(shù)據(jù)終點點號位置（起始點號為1）!NF: 正、反變換標志量(1:反變換;2:正變換)!輸出參數(shù)說明：!Freal：信號的實部!Fimage：信號的虛部(對于實信號而言，賦值為0)!對應頻率分布說明：!數(shù)據(jù)Freal(m,n)和Fimage(m,n)對應的頻率分布位置為:!m 方向： 0,1,.,m/2-1,m/2,-(m/2-1)

29、,.,-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,n) allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOP DO i=1,m,1

30、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 Treal(i)=Freal(i,j) Timage(i)=Fimage(i,j) end do call FFT(Treal,Timage,m,nf) do i=1

31、,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!功能：復數(shù)組1-D快速Fourier變換!輸入?yún)?shù)說明：!Xreal(n): 輸入數(shù)據(jù)實部 !Ximage(n): 輸入數(shù)據(jù)虛部 !N: 點數(shù)（N 必須為 2 的冪次方）!NF: 正、反變換標志量(1:反變換;2:正變換)!輸出參數(shù)說明：!Xreal(n): 變換后頻譜實部 !Ximage(n): 變換后頻譜虛部 !對應頻率分布說

32、明：!數(shù)據(jù)Xreal(n)和Ximage(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.+0.001)n2=n/2nu1=nu-1f=float(-1)*nf)k=0 DO l=1,nu,1 DO

33、while (k.lt.n) do i=1,n2,1 p=ibitr(k/2*nu1,nu) arg=6.*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 en

34、d 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 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

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

36、y方向)!輸出參數(shù)說明：!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!輸入?yún)?shù)說明：!dx: x 方向點距 !dy: y 方向線距!m: 點數(shù)（M 必須為 2 的

37、冪次方）!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(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*(y

38、y-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!* ! 子程序：factor! 功能：計算x,y,z方向導數(shù)的實部和虛部! 輸入?yún)?shù)說明：!m: 點數(shù)（M 必須為 2 的冪次方）!n: 線數(shù)（N 必須為 2 的冪次方）! field_re:初始觀測面信號的實部! field_im：初始觀測面信號的虛部 !W(m,n)：徑向圓頻

39、率! 輸出參數(shù)說明：! Px_re:x方向導數(shù)信號的實部! Px_im:x方向導數(shù)信號的虛部! Py_re:y方向導數(shù)信號的實部! Py_im:y方向導數(shù)信號的虛部! Pz_re:z方向導數(shù)信號的實部! Pz_im:z方向導數(shù)信號的虛部 !* subroutine factor(m,n,field_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:

40、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. do i=1,m,1 do j=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

41、.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! 功能：計算異常的水平導數(shù)、垂向導數(shù)及總導數(shù)! 輸入?yún)?shù)說明：!m: 點數(shù)（M 必須為 2 的冪次方）!n: 線數(shù)（N 必須為 2 的冪次方）! Px_re:x方向導數(shù)信號的實部! Px_im:x方向導數(shù)信號的虛部! Py_re:y方向導數(shù)信號的實部! Py_im:y方向導數(shù)信號的虛部! Pz_re:z方向導數(shù)信號的實部! Pz_im:z方向導數(shù)信號的虛部 ! 輸出參數(shù)說明：! field_vdr:對場值