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

《重磁资料处理与解释》实验四位场边缘识别程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本原理地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值。故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。(1)垂向导数：垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用(,,)(,,)gxyzVDRxyzz(1.1)（2）解析信号振幅：解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常22ASMTHDRVDR(1.2)（3）总水平导数（THDR）22(,,)(,,)(,,)gxyzgxyzTHDRxyzxy(1.3)2、输入/输出数据格式设计依据上述原理，现在对上述各种边缘识别方法进行程序设计。2.1输入数据格式设计本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。例如：DSAA201201-1000.0000001000.000000-1000.0000001000.0000005.549671E-0123.5398465.549671E-015.634658E-015.721339E-015.808522E-015.897312E-015.987253E-016.078691E-016.171604E-01…2.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-3.2723978E-02-3.2787599E-02-3.2927759E-02-3.3138681E-02…2.3参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为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文件名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文件名gravity_asm.grd、扩边比例因子factor_m。确定扩边网格的大小m*n(m,n均为2的幂次方)从输入数据文件名中读取数据对原始数据进行二维余弦扩边对扩边后的数据进行快速二维傅里叶正变换将傅里叶变换后的数据与导数因子相乘求出重力数据在x,y,z方向的导数对导数进行快速二维傅里叶反变换分别求出VDR、THDR、ASM值。去除扩边部分后输出结果图4.1总体设计N-S图4.测试结果分析：由图4.3可看出，VDR方法的零值线可较准确识别模型体的边界；由图4.4可看出，ASM的极大值点边界可大致识别模型体边界，但精度不是很高。对比图4.2到4.5可以看出，THDR方法对模型边界的识别效果是最好的。5结论及建议经测试，VDR与THDR对模型体的边界位置识别效果较好，而ASM对模型体边界识别效果较差。三种方法中，THDR效果最好。图4.2重力异常原始图图4.3垂向导数（VDR）图4.4解析信号振幅（ASM）图4.3ASM（解析信号振幅）图4.5总水平导数（THDR）图4.2重力异常原始图附录：边缘识别程序源代码!****************************************************************************************************!程序功能：实现频率域位场导数运算进行边缘识别!cmd文件参数:!cmdfile:存放有关参数的文件名变量!input_file:观测面位场数据文件!output_file_vdr：场值的水平导数数据文件!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_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,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(:,:)integerN_point,N_lineintegerm,n,m1,m2,n1,n2realfactor_mrealxmin,xmax,ymin,ymax,dx,dycmdfile='deviation.cmd'callread_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,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_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_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)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)endprogram!****************************************************************************!子程序：read_cmd!功能：读取参数文件!输入参数说明：!cmdfile：参数文件名!输出参数说明：!input_file:观测面位场数据文件!output_file_vd