《实验四-位场边缘识别程序设计实验汇总.doc》由会员分享，可在线阅读，更多相关《实验四-位场边缘识别程序设计实验汇总.doc（20页珍藏版）》请在三一办公上搜索。
1、重磁资料处理与解释实验四位场边缘识别程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本原理地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数
2、达到零值。故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。(1)垂向导数：垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用 ( 1.1) （2）解析信号振幅：解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于
3、重、磁力异常 (1.2) （3）总水平导数（THDR） (1.3)2、输入/输出数据格式设计依据上述原理，现在对上述各种边缘识别方法进行程序设计。2.1 输入数据格式设计本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。例如：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 输出数据格式设计计算结果输出数据格式与输入格式对应，格式为.GRD格式，均为实型变量。例如：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 参数文件数据格式设计将以上部分量保
5、存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。在该文件中保存的参数如下：输入数据文件名input_file，字符串变量，长度不超过80；输出vdr数据文件名output_file_vdr，字符串变量，长度不超过80；输出thdr数据文件名output_file_thdr，字符串变量，长度不超过80;输出asm数据文件名output_file_asm，字符串变量，长度不超过80factor_m：扩边比例因子，实型变量（1）；3. 总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入数据文件名gravity.grd、输出vdr文件
6、名field_vrd.grd、输出thdr文件名field_thdr.grd、输出asm文件名field_asm.grd、扩边比例因子factor_m；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来在频率域求出在x,y,z方向的导数并反变换；最后求出VDR、THDR、ASM数据。最后去除扩边部分后输出。总体设计见表1。输入参数：输入数据文件名gravity.grd、输出vdr文件名gravity_vrd.grd、输出thdr文件名gravity_thdr.grd、输出asm文件
7、名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.2 重力异常原始图图4.5 总水平导数（THDR）图4.4 解析信号振幅（ASM）图4.3 ASM（解析信号振幅）分析：由图4.3可看出，VD
8、R方法的零值线可较准确识别模型体的边界；由图4.4可看出，ASM的极大值点边界可大致识别模型体边界，但精度不是很高。对比图4.2到4.5可以看出，THDR方法对模型边界的识别效果是最好的。5 结论及建议 经测试，VDR与THDR对模型体的边界位置识别效果较好，而ASM对模型体边界识别效果较差。三种方法中，THDR效果最好。附录：边缘识别程序源代码!*! 程序功能：实现频率域位场导数运算进行边缘识别! cmd文件参数:! cmdfile:存放有关参数的文件名变量! input_file:观测面位场数据文件! output_file_vdr：场值的水平导数数据文件! output_file_thd
9、r:场值垂向导数数据文件! 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_im：初始观测面信号的虚部
10、! 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:场值的解析信号振幅的结果!*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!功能：读取参数文件!输入参数说明：!cmdfile：参数文件名!输出参数说明：! input_file:观测面位场数据文件! output_file_vdr：对场值作水平导数处理后的数据文件! output_file_thdr:对场值作垂向导数处理后的数据文件! output_file
17、_asm:对场值作总导数处理后的数据文件!factor_m：扩边因子!*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文件中读取相关参数!输入参数说明：!filename_obser：输入文件名!输出参数说明：!N_point,N_line：点数、线数!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!功能：确定扩边数据点号位置!输入参数说明：!factor_m: 扩边比例因子（1.0）!a,b：点数、线数!输出参数说明：!m1,m2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!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文件中的数据!输入参数说明：!filename_obser：输入文件名!m1,m2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!n:扩边后数据终点位置点号(起点位置点号为1
23、)!输出参数说明：!A：存放输出数据的二维数组名!*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!功能：二维扩边子程序并为信号虚部赋值!输入参数说明：!m1,m2: x方向实际数据起点位置和终点位置点号!m:扩边后数据终点位置点号(起点位置点号为1)!n1,n2: y方向实际数据起点位置和终点位置点号!n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部!Ui：初始观测面信号的虚部!输出参数说明：! Ur:初始观测面
25、信号的实部!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!功能：一维扩边子程序!输入参数说明：!n0,n3：扩边后数据起点位置和终点位置!n1,n2：实际数据起点位置和终点位置!feild
27、(i),(i=n1,n1+1,.,n2)：实际数据!输出参数说明：!field(i),(i=n0,.,n1-1)：扩边后左边的数据!field(i),(i=n2+1,.,n3)：扩边后右边的数据!*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!功能：复数组2-D快速Fourier变换!输入参数说明：!m0,m3：x方向的起点和终点!n0,n3：y方向的起点和终点!field：输入信号(需要赋值给Freal，实部)!m,n：x,y方向扩边后数据终点点号位置（起始点号为1）!NF: 正、反变
29、换标志量(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!* 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!功能：复数组1-D快速Fourier变换!输入参数
32、说明：!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!*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!功能：计算x,y方向的点距!输入参数说明：!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!N_point,N_line：点数(x方向)、线数(y方向)!输出参数说明：!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!*
