实验四位场边缘识别程序设计实验汇总

1、实验四位场边沿鉴识程序设计实验汇总实验四位场边沿鉴识程序设计实验汇总20/20实验四位场边沿鉴识程序设计实验汇总重磁资料办理与解说实验四位场边沿鉴识程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本源理地质目标体边沿时指断裂结构线、不一样样地质体的界限限，其实是拥有必然密度或磁性差其余地质体的界限限，在地质体的边沿周边，重、磁异样变化率较大，故全部的边沿鉴识方法均利用这一特色进行设计。此刻有重、磁位场边沿鉴识方法分为数理统计、数值计算以及其余三大类。数值类边沿鉴识方法均利用极大值地点或零值地点确立地质体的边沿地点，其依

2、照的理论基础是二度体铅垂台阶模型的重力异样特色。在该模型边沿处重力异样总水平导数和分析信号振幅达到极大值、垂导游数达到零值。故能够利用这些特色地点来确立二度体铅垂台阶的边沿地点，确立倾斜二度体、不规则二度体及三度体边沿地点的理论均为二度体铅垂台阶模型理论的推行，但确立的边沿地点和真切的地点有必然的误差。该误差跟着地质体界限形状、埋深、水平尺寸及物性差别等的变化而变化。所以，边沿鉴识结果是一种定性或半定量解说结果，与定量解说结果有必然差别，鉴识结果可作为边沿地点定量反演的初值。(1)垂导游数：垂导游数方法利用零值地点确立地质体的边沿地点，重力异样能够直接使用，对磁力异样必然转变为磁源重力异样或化

3、极磁力异样才能够使用VDR(,xy,)zg(x,y,z)z(1.1)（2）分析信号振幅：分析信号振幅也是利用极大值地点来确立地质体的边沿地点合用于重、磁力异样ASMTHDR2VDR2(1.2)（3）总水平导数（THDR）g(x,y,z)2g(x,y,z)2THDR(x,y,z)y(1.3)x2、输入/输出数据格式设计依照上述原理，此刻对上述各样边沿鉴识方法进行程序设计。2.1输入数据格式设计本次实验给了正演的重力异样数据，为.GRD格式，均为实型变量。比方：DSAA201201-1000.0000001000.000000-1000.0000001000.0000005.549671E-012

4、3.5398465.549671E-015.634658E-015.721339E-015.808522E-015.897312E-015.987253E-016.078691E-016.171604E-012.2输出数据格式设计计算结果输出数据格式与输入格式对应，格式为.GRD格式，均为实型变量。比方：DSAA201201-1000.0001000.000-1000.0001000.000-0.14650840.3190881-3.6523044E-02-3.3485338E-02-3.3061244E-02-3.2748722E-02-3.2688729E-02-3.2654848E-02

5、-3.2723978E-02-3.2787599E-02-3.2927759E-02-3.3138681E-022.3参数文件数据格式设计将以上部重量保留在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超出80，全路径名。在该文件中保留的参数以下：输入数据文件名input_file，字符串变量，长度不超出80；输出vdr数据文件名output_file_vdr，字符串变量，长度不超出80；输出thdr数据文件名output_file_thdr，字符串变量，长度不超出80;输出asm数据文件名output_file_asm，字符串变量，长度不超出80factor_m：扩边比率因子，

6、实型变量（1）；3.整体设计此次程序采纳IPO结构设计，第一经过读取cmd文件，获得有关输入参数：输入数据文件名gravity.grd、输出vdr文件名field_vrd.grd、输出thdr文件名field_thdr.grd、输出asm文件名field_asm.grd、扩边比率因子factor_m；此后确立确立扩边网格的大小，扩边数据点号地点；再从观察面位场数据文件中读取数据。下一步，进行二维余弦扩边，将扩完边的数据进行迅速二维傅里叶变换，变换到频次域；接下来在频次域求出在x,y,z方向的导数并反变换；最后求出VDR、THDR、ASM数据。最后去除扩边部分后输出。整体设计见表1。输入参数：输

7、入数据文件名gravity.grd、输出vdr文件名gravity_vrd.grd、输出thdr文件名gravity_thdr.grd、输出asm文件名gravity_asm.grd、扩边比率因子factor_m。确立扩边网格的大小m*n(m,n均为2的幂次方)从输入数据文件名中读取数据对原始数据进行二维余弦扩边对扩边后的数据进行迅速二维傅里叶正变换将傅里叶变换后的数据与导数因子相乘求出重力数据在x,y,z方向的导数对导数进行迅速二维傅里叶反变换分别求出VDR、THDR、ASM值。去除扩边部分后输出结果图4.1整体设计N-S图4.测试结果图4.2重力异样原始图图4.3垂导游数（VDR）图4.4

8、分析信号振幅（ASM）图4.5总水平导数（THDR）分析：由图4.3可看出，VDR方法的零值线可较正确鉴识模型体的界限；由图4.4可看出，ASM的极大值点界限可大概鉴识模型体界限，但精度不是很高。比较图4.2到4.5能够看出，THDR方法对模型界限的鉴识见效是最好的。结论及建议经测试，VDR与THDR对模型体的界限地点鉴识见效较好，而ASM对模型体界限鉴识见效较差。三种方法中，THDR见效最好。附录：边沿鉴识程序源代码!*!程序功能：实现频次域位场导数运算进行边沿鉴识cmd文件参数:cmdfile:寄存有关参数的文件名变量input_file:观察面位场数据文件output_file_vdr：

9、场值的水平导数数据文件output_file_thdr:场值垂导游数数据文件output_file_asm:场值的分析信号振幅数据文件!factor_m：扩边因子.grd文件参数：!N_point,N_line：点数(x方向)、线数(y方向)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值Ur:初始观察面场值扩边参数：!m1,m2：x方向实质数据起点和终点点号地点!1,m：x方向扩边后数据起点和终点点号地点!n1,n2：y方向实质数据起点和终点点号地点!1,n3：y方向扩边后数据起点和终点点号地点求导参数：field_re:初始观察面信号的实部field

10、_im：初始观察面信号的虚部Px_re:x方导游数信号的实部Px_im:x方导游数信号的虚部Py_re:y方导游数信号的实部Py_im:y方导游数信号的虚部Pz_re:z方导游数信号的实部Pz_im:z方导游数信号的虚部!W(m,n)：径向圆频次field_vdr:对场值作水平导数的结果field_thdr:对场值作垂导游数的结果field_asm:场值的分析信号振幅的结果!*programdeviationparameter(eigval=3.701411*1e5)character*(80)cmdfilecharacter*80input_file,output_file_vdr,outp

11、ut_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(:,:)integerN_point,N_lineintegerm,n,m1,m2,n1,n2realfact

12、or_mrealxmin,xmax,ymin,ymax,dx,dycmdfile=deviation.cmdcallread_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)callread_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)callcalculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(field_re(1:m,1:n),field_im(1:m,

13、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)callinput_grd(field_re,input_file,m1,m2,n1,n2,m,n)callexpand_cos_2D(m1,m2,m,n1,n2,n,field_re,field_

14、im)callFFT2(field_re,field_im,m,n,2)CALLcal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)callWAVE2D(m,n,dx,dy,U,V,W)callfactor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)callFFT2(px_re,px_im,m,n,1)callFFT2(py_re,py_im,m,n,1)callFFT2(pz_re,pz_im,m,n,1)calldeviration(m,n,px_re,py

15、_re,pz_re,field_vdr,field_thdr,field_asm)callOUTPUT_GRD(field_vdr,output_file_vdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)callOUTPUT_GRD(field_thdr,output_file_thdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)callOUTPUT_GRD(field_asm,output_file_asm,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)deall

16、ocate(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)endprogram!*!子程序：read_cmd功能：读取参数文件输入参数说明：!cmdfile：参数文件名输出参数说明：input_file:观察面位场数据文件output_file_vdr：对场值作水平导数办理后的数据文件output_file_thdr:对场值作垂导游数办理后的数据文件output_file_asm:对场值作总导数办理后的数据文件!factor_m：扩边因子!*Subroutin

17、eread_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)implicitnonecharacter*80strcharacter*(*)cmdfilecharacter*(*)input_file,output_file_vdr,output_file_thdr,output_file_asmrealfactor_mopen(10,file=cmdfile,status=old)read(10,*)str,input_fileread(10,*)str,output_file_

18、vdrread(10,*)str,output_file_thdrread(10,*)str,output_file_asmread(10,*)str,factor_mclose(10)endSubroutineread_cmd!*!子程序：read_grd功能：从原始观察.grd文件中读取有关参数输入参数说明：!filename_obser：输入文件名输出参数说明：!N_point,N_line：点数、线数!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!*subroutineread_grd(input_file,N_point,N_line,Xmi

19、n,Xmax,Ymin,Ymax)implicitnonecharacter*(*)input_fileintegerN_point,N_linerealXmin,Xmax,Ymin,Ymaxopen(10,file=input_file,status=old)Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)endsubroutineread_grd!*!子程序：calculate_mn功能：确立扩边数据点号地点输入参数说明：!factor_m:扩边比率因子（1.0）!a,b：点数

20、、线数输出参数说明：!m1,m2:x方向实质数据起点地点和终点地点点号!m:扩边后数据终点地点点号(起点地点点号为1)!n1,n2:y方向实质数据起点地点和终点地点点号!n:扩边后数据终点地点点号(起点地点点号为1)!*subroutinecalculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicitnoneintegera,b,m,n,m1,m2,n1,n2integermtemp,mu,nurealfactor_mmtemp=aDOWHILE(mod(mtemp,2).eq.0).and.(mtemp.ne.0)mtemp=mtemp/2EnddoIF

21、(mtemp.eq.1)THENm=a*2ELSEmu=int(log(float(a)/0.693147+factor_m)m=2*muENDIFm1=1+(m-a)/2m2=m1+a-1write(*,*)m,apausemtemp=bDOWHILE(mod(mtemp,2).eq.0).and.(mtemp.ne.0)mtemp=mtemp/2EnddoIF(mtemp.eq.1)THENn=b*2ELSEnu=int(log(float(b)/0.693147+factor_m)n=2*nuENDIFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,

22、m,npauseendsubroutinecalculate_mn!*!子程序：INPUT_GRD功能：读取grd文件中的数据输入参数说明：!filename_obser：输入文件名!m1,m2:x方向实质数据起点地点和终点地点点号!m:扩边后数据终点地点点号(起点地点点号为1)!n1,n2:y方向实质数据起点地点和终点地点点号!n:扩边后数据终点地点点号(起点地点点号为1)输出参数说明：!A：寄存输出数据的二维数组名!*SUBROUTINEINPUT_GRD(A,input_file,m1,m2,n1,n2,m,n)character*(*)input_fileintegerm1,m2,n1

23、,n2,m,nrealxmin,xmax,ymin,ymaxrealA(1:m,1:n)reali,j,kdoj=1,n,1doi=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=input_file,status=old)read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*)(A(i,j),i=m1,m2),j=n1,n2)Close(20)ENDSUBROUTINEINPUT_GRD!*!子程序：expand_cos_2D功能：二维扩边子程序

24、并为信号虚部赋值输入参数说明：!m1,m2:x方向实质数据起点地点和终点地点点号!m:扩边后数据终点地点点号(起点地点点号为1)!n1,n2:y方向实质数据起点地点和终点地点点号!n:扩边后数据终点地点点号(起点地点点号为1)Ur:初始观察面信号的实部!Ui：初始观察面信号的虚部输出参数说明：Ur:初始观察面信号的实部!Ui：初始观察面信号的虚部!*subroutineexpand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicitnoneintegerm,n,m1,m2,n1,n2realUr(1:m,1:n),Ui(1:m,1:n)real,allocatable:

25、u(:),r(:)integerj,i,kallocate(u(1:m)doj=n1,n2,1doi=1,m,1u(i)=Ur(i,j)enddocallexpand_cos_1d(1,m1,m2,m,u(1)doi=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n)doi=1,m,1doj=1,n,1r(j)=Ur(i,j)enddocallexpand_cos_1d(1,n1,n2,n,r(1)doj=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)doi=1,mdoj=1,nUi(i,j)=0en

26、ddoenddoendsubroutineexpand_cos_2D!*!子程序：expand_cos_1d功能：一维扩边子程序输入参数说明：!n0,n3：扩边后数据起点地点和终点地点!n1,n2：实质数据起点地点和终点地点!feild(i),(i=n1,n1+1,.,n2)：实质数据输出参数说明：!field(i),(i=n0,.,n1-1)：扩边后左侧的数据!field(i),(i=n2+1,.,n3)：扩边后右侧的数据!*Subroutineexpand_cos_1d(n0,n1,n2,n3,Field)RealField(n0:n3)pi=3.141592654Field(n0)=(F

27、ield(n1)+Field(n2)/2.0Field(n3)=Field(n0)i1=n0i2=n1DOi=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1)Enddoi1=n2i2=n3DOi=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1)EnddoEndsubroutineexpand_cos_1d!*!功能：FFT2功能：复数组2-D迅速Fourier变换输入参数说明：!m0,m3：x方向

28、的起点和终点!n0,n3：y方向的起点和终点!field：输入信号(需要赋值给Freal，实部)!m,n：x,y方向扩边后数据终点点号地点（初步点号为1）!NF:正、反变换标记量(1:反变换;2:正变换)输出参数说明：!Freal：信号的实部!Fimage：信号的虚部(关于实信号而言，赋值为0)对应频次散布说明：!数据Freal(m,n)和Fimage(m,n)对应的频次散布地点为:!m方向：0,1,m/2-1,m/2,-(m/2-1),-1!n方向：0,1,n/2-1,n/2,-(n/2-1),-1!*SUBROUTINEFFT2(Freal,Fimage,m,n,nf)implicitno

29、neINTEGERm,n,nfREALFreal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE:Treal(:),Timage(:)integernmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0)STOPDOi=1,m,1IF(n.ne.1)THENdoj=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)enddocallFFT(Treal,Timage,n,nf)doj=1,n,1Fr

30、eal(i,j)=Treal(j)Fimage(i,j)=Timage(j)enddoENDIFENDDODOj=1,n,1IF(m.ne.1)THENdoi=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)enddocallFFT(Treal,Timage,m,nf)doi=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)enddoENDIFENDDOdeallocate(Treal,Timage,STAT=ierr)ENDSUBROUTINEFFT2!*!子程序：FFT功能：复数组1-D迅速Fourier变

31、换输入参数说明：!Xreal(n):输入数据实部!Ximage(n):输入数据虚部!N:点数（N必然为2的幂次方）!NF:正、反变换标记量(1:反变换;2:正变换)输出参数说明：!Xreal(n):变换后频谱实部!Ximage(n):变换后频谱虚部对应频次散布说明：!数据Xreal(n)和Ximage(n)对应的频次散布地点为!0,1,n/2-1,n/2,-(n/2-1),-1:!*SUBROUTINEFFT(Xreal,Ximage,n,nf)implicitnoneINTEGERn,nfREALXreal(1:n),Ximage(1:n)integernu,n2,nu1,k,k1,k1n2

32、,l,i,ibitrrealf,p,arg,c,s,treal,timagenu=int(log(float(n)/0.693147+0.001)n2=n/2nu1=nu-1f=float(-1)*nf)k=0DOl=1,nu,1DOwhile(k.lt.n)doi=1,n2,1p=ibitr(k/2*nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2

33、)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1enddok=k+n2ENDDOk=0nu1=nu1-1n2=n2/2ENDDODOk=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k)thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageendifENDDOIF(nf.ne.

34、1)THENdoi=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)enddoENDIFENDSUBROUTINEFFTFUNCTIONIBITR(J,NU)implicitnoneintegeribitr,j,nuintegerj1,itt,i,j2j1=jitt=0doi=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2enddoENDFUNCTIONIBITR!*!子程序：cal_dxdy功能：计算x,y方向的点距输入参数说明：!x_min,x_max：x的最小值和最大值!y

35、_min,y_max：y的最小值和最大值!N_point,N_line：点数(x方向)、线数(y方向)输出参数说明：!dx,dy：x,y方向的点距!*subroutinecal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicitnonerealxmin,xmax,ymin,ymaxintegerN_POINT,N_LINErealdx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEendsubroutinecal_dxdy!*!子程序：WAVE2D功能：计算2D径向圆频次W输入参数说明：!dx

36、:x方向点距!dy:y方向线距!m:点数（M必然为2的幂次方）!n:线数（N必然为2的幂次方）输出参数说明：!W(m,n):径向圆频次!*SUBROUTINEWAVE2D(m,n,dx,dy,U,V,W)implicitnoneINTEGERm,nREALdx,dyREALW(1:m,1:n),U(1:m),V(1:n)realpi,delx,delyintegermidm,midn,i,j,xx,yymidm=m/2+1midn=n/2+1delx=float(m)/dxdely=float(n)/dydoj=1,n,1yy=jif(yy.gt.midn)yy=yy-nv(j)=dely*(

37、yy-1)enddodoi=1,m,1xx=iif(xx.gt.midm)xx=xx-mu(i)=delx*(xx-1)enddodoj=1,n,1doi=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j)enddoenddoENDSUBROUTINEWAVE2D!*!子程序：factor功能：计算x,y,z方导游数的实部和虚部!输入参数说明：!m:点数（M必然为2的幂次方）n:线数（N必然为2的幂次方）field_re:初始观察面信号的实部field_im：初始观察面信号的虚部!W(m,n)：径向圆频次输出参数说明：Px_re:x方导游数信号的实部Px_im:x方导游数

38、信号的虚部Py_re:y方导游数信号的实部Py_im:y方导游数信号的虚部Pz_re:z方导游数信号的实部Pz_im:z方导游数信号的虚部!*subroutinefactor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)implicitnoneintegerm,nrealfield_re(1:m,1:n),field_im(1:m,1:n)realpx_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:

39、m,1:n)realU(1:m),V(1:n),W(1:m,1:n)integeri,jrealpipi=3.1415926doi=1,m,1doj=1,n,1px_re(i,j)=field_im(i,j)*u(i)*(-1)*pi*2.0px_im(i,j)=field_re(i,j)*u(i)*pi*2.0py_re(i,j)=field_im(i,j)*v(j)*(-1)*pi*2.0py_im(i,j)=field_re(i,j)*v(j)*pi*2.0pz_re(i,j)=field_re(i,j)*w(i,j)*pi*2.0pz_im(i,j)=field_im(i,j)*w(i,j)*pi*2.0enddoenddoendsubroutinefactor!*!子程序：deviration功能：计算异样的水平导数、垂导游数及总导数!输入参数说明：!m:点数（M必然为2的幂次方）n:线数（N必然为2的幂次方）Px_re:x方导游数信号的实部Px_im:x方导游数信号的虚部Py_re:y方导游数信号的实部Py_im:y方导游数信号的虚部Pz_re: