实验二边缘识别实验报告正文

一 、基本原理 21.1 位场边缘识别方法研究进展概述 .21.2 数值计算类边缘识别方法的研究及计算 .21.2.1 垂向导数（VDR） .21.2.2 总水平导数（THDR） .21.2.3 解析信号振幅（ASM） .21.2.4 倾斜角（TA） .31.2.5 二阶导数类边缘识别 3二、 输入/输出数据格式设计 .32.1 主要变量设计 .32.2 位场边缘识别程序数据格式设计 4三、 总体设计 5四、测试结果 54.1 测试用例 .54.2 测试结果 .64.3 结果分析 .8五、结论与建议 95.1 结论 .95.2 建议 .9附录：边缘识别程序源代码 92一 、基本原理 1.1 位场边缘识别方法研究进展概述地质目标体边缘是指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。现有重、磁位场边缘识别方法分为数理统计、数值计算和其他三大类。1.2 数值计算类边缘识别方法的研究及计算数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置与真实位置有一定偏差该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。1.2.1 垂向导数（VDR）垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可直接使用，对磁力异常必须转换成磁源重力异常（假重力异常）或化极磁力异常才可以使用。垂向导数方法研究历史较早，方法较成熟，应用频率较高通过我们的试验和研究认为：随着地质体埋深的增大，垂向导数所确定的边缘位置偏离实际位置的距离越大。1.2.2 总水平导数（THDR）总水平导数是利用其极大值位置来确定地质体的边缘位置，适用于重力异常，对磁力异常必须换算成磁源重力异常或化极磁力异常才可以使用。31.2.3 解析信号振幅（ASM）解析信号振幅也是利用极大值位置来确定地质体的边缘位置，适用于重力异常和磁力异常。1.2.4 倾斜角（TA）倾斜角实质上是垂向导数和总水平导数的比值。由于倾斜角为一阶导数的比值，所以能很好地平衡高幅值异常和低幅值异常，起到边缘增强的效果。1.2.5 二阶导数类边缘识别为了增强边缘分辨能力，和克服当总水平导数等于时倾斜角存在“解析奇点” ，会使得计算结果不稳定。2、输入/输出数据格式设计依据上述原理，现对上述各种边缘识别方法实现进行程序设计。2.1 主要变量设计 cmd_file:参数文件名 input_field_filename： 输入位场文件名 output_field_filename：输出转换后位场文件名 expan_2D_method： 2D 扩边方法选择 Field_Deriv_method：边缘识别方法选择 num_area：计算区域大小选择 mpoint,nline：点数，线数 m0,m1,m2,m3：扩边后 X 方向上各端点 n0,n1,n2,n3：扩边后 Y 方向上各端点 Xmin,Xmax：X 坐标最小，最大值 ymin,Ymax：Y 坐标最小，最大值 Ori_Field：原始场值，二维数组，m3*n3 Deriv_Field：转换后后场值，二维数组， m3*n3 Deriv_operator_x：x 方向转换因子，二维数组，m3*n3 Deriv_operator_y：y 方向转换因子，二维数组，m3*n3 Deriv_operator_z：z 方向转换因子，二维数组，m3*n342.2 位场边缘识别程序数据格式设计 cmd 文件，输入参数文件，格式如下：（cmd.txt）input_field_filename anomaly.grdoutput_field_filename VDR_THDR_anomaly.grdnum_area 0expan_2D_method 1Field_Deriv_method 7说明：num_area: 取 0：表示一般扩边区域取 1：表示计算区域相对于一般扩边再放大一倍expan_2D_method:取 1：表示余弦扩边，端点值取边界平均值取 2：表示余弦扩边，端点值取全场区平均值取 3：表示余弦扩边，端点值取常数 0取 4：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取边界平均值取 5：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取全场区平均值取 6：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取常数 0Field_Deriv_method:取 1：计算垂向导数 VDR取 2: 计算总水平导数 THDR取 3: 计算解析信号振幅 ASM取 4: 计算倾斜角 TA取 5:计算垂向二阶导数 VDR_VDR取 6:计算倾斜角总水平导数 TA_THDR取 7:计算垂向导数总水平导数 VDR_THDR 原始位场文件，输入文件，grd 格式(ASCII) 延拓后位场文件，输出文件，grd 格式(ASCII) Grd 格式如下： DSAAmpoint nlineXmin XmaxYmin YmaxZmin Zmax Z11 Z12 . . 53、总体设计四、测试结果4.1 测试用例（1）观测面上的重力异常存放在“anomaly.grd”中，坐标单位为 m。（2）形体边界存放在“rectangle.bln”中， 坐标单位图 3.1 重、磁异常边缘识别程序设计流程图64.2 测试结果图 4.2.2 THDR（总水平导数）结果图-10864-20681-204681.4.30.648527-10864-2081-204681.643.5026图 4.2.3 ASM(解析信号振幅)结果图图 4.2.5 VDR_VDR(垂向二阶导)结果图-10864-2081-20-.9765.3-图 4.2.4 AT(倾斜角)边缘识别结果图-0-204681-81-9753809图 4.2.6 VDR_THDR 结果图(垂向导数总水平导数，余弦扩边方法)0-420681204E-5.3.84.3 结果分析1) 从图 4.2.1-图 4.2.7 可以看出，VDR（垂向导数） 、VDR_VDR(垂向二阶导)、AT( 倾斜角) 实验结果零值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现是正确的。2）同时从图 4.2.1-图 4.2.7 可以也可看出，THDR（总水平导数） 、VDR_THDR（垂向导数总水平导数） 、AT_THDR(倾斜角总水平导数) 、ASM(解析信号振幅)实验结果极大值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现也是正确的。3）该实验模型，从以上各方法边缘识别效果来看，ASM(解析信号振幅)分辨能力较低；二阶导数类方法识别效果比一阶导数类识别效果好；就一阶导数类来看 AT(倾斜角 )由于其他三种方法；二阶导数类方法中 AT_THDR(倾斜角总水平导数) 识别最精确，但是从（图 4.2.7）看出，其稳定性较差，故相比较而言VDR_THDR（垂向导数总水平导数） 、VDR_VDR( 垂向二阶导) 效果是比较理想的。4）对比图 4.2.6 和图 4.2.8，可以看出在频率域处理选用扩边方法不同对计算结果有一定影响，其中最小曲率扩边计算效果为佳。-10864-2081-204681.53679.2图 4.2.7 TA_THDR(倾斜角总水平导数)结果图-04-20681-04681E-5.4.306图 4.2.8 VDR_THDR 结果图(垂向导数总水平导数，最小曲率扩边)9五、结论与建议5.1 结论此方法在一定程度上是有效的，且从实验结果可以垂向导数总水平导数和垂向二阶导识别的效果较好。5.2 建议没有学好位场边缘识别方法没有可以能提出的建议附录：边缘识别程序源代码!*! PROGRAM: Fre_PoteField_Deriv! PURPOSE: Frequency Potential field Derivative! Author : gesang! Data: 2013-11!*PROGRAM Fre_PoteField_Derivimplicit noneCHARACTER *(80) input_field_filenameCHARACTER *(80) output_field_filenameINTEGER expan_2D_method,num_area,Field_Deriv_methodINTEGER mpoint,nlineINTEGER m0,m1,m2,m3INTEGER n0,n1,n2,n3REAL Xmin,XmaxREAL ymin,YmaxREAL dx,dyREAL,ALLOCATABLE: Ori_Field(:,:) REAL,ALLOCATABLE: Field_Real(:,:)REAL,ALLOCATABLE: Field_Image(:,:) REAL,ALLOCATABLE: Deriv_Field(:,:) REAL,ALLOCATABLE: Deriv_operator_x(:,:) REAL,ALLOCATABLE: Deriv_operator_y(:,:) REAL,ALLOCATABLE: Deriv_operator_z(:,:) !z方向导数因子!INPUT: !读取文件参数 CALL read_cmdinfo(input_field_filename,output_field_filename,&expan_2D_method,Field_Deriv_method,num_area)!读取点数线数及测区范围CALL get_mpoint_and_nline(input_field_filename,mpoint,nline,Xmin,Xmax,&Ymin,Ymax)!计算扩边端点CALL get_expan_num(mpoint,m0,m1,m2,m3,num_area) 10CALL get_expan_num(nline, n0,n1,n2,n3,num_area) ALLOCATE(Ori_Field(m0:m3,n0:n3)ALLOCATE(Field_Real(m0:m3,n0:n3)ALLOCATE(Field_Image(m0:m3,n0:n3)ALLOCATE(Deriv_Field(m0:m3,n0:n3)ALLOCATE(Deriv_operator_x(m0:m3,n0:n3)ALLOCATE(Deriv_operator_y(m0:m3,n0:n3)ALLOCATE(Deriv_operator_z(m0:m3,n0:n3)!从网格化文件读入初始场值CALL Read_field(input_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field) !PROCESS:!扩边dx=(Xmax-Xmin)/(mpoint-1)dy=(Ymax-ymin)/(nline -1)CALL expan_2D(m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,dx,dy,&expan_2D_method) !expan_2D_method:扩边方法选择!计算导数因子 CALL Deriv_operator_sub(dx,dy,m3,n3,Deriv_operator_x,&Deriv_operator_y,Deriv_operator_z)!求导运算Field_Real=Ori_FieldField_Image=0.0 CALL field_Deriv(Field_Real,Field_Image,Deriv_operator_x,Deriv_operator_y,&Deriv_operator_z,m0,m1,m2,m3,n0,n1,n2,n3,Field_Deriv_method) Deriv_Field=Field_Real! OUTPUT:!输出求导后位场文件CALL output_grd(output_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,&mpoint,nline,Deriv_Field,Xmin,Xmax,Ymin,Ymax)DEALLOCATE (Ori_Field)DEALLOCATE (Field_Real)DEALLOCATE (Field_Image)DEALLOCATE (Deriv_Field)DEALLOCATE (Deriv_operator_x)DEALLOCATE (Deriv_operator_y)DEALLOCATE (Deriv_operator_z)END PROGRAM Fre_PoteField_Deriv!*开始*! *得到输入参数子程序开始*SUBROUTINE read_cmdinfo(input_field_filename,output_field_filename,&expan_2D_method,Field_Deriv_method,num_area)CHARACTER *(*) input_field_filenameCHARACTER *(*) output_field_filenameINTEGER expan_2D_method,Field_Deriv_method,num_areaINTEGER unit11CHARACTER*80 s CALL get_unit(unit)OPEN(unit,file=cmd.txt,status=old)READ(unit,*) s ,input_field_filenameREAD(unit,*) s ,output_field_filenameREAD(unit,*) s ,num_areaREAD(unit,*) s ,expan_2D_methodREAD(unit,*) s ,Field_Deriv_methodCLOSE(unit)END SUBROUTINE read_cmdinfo! *得到输入参数子程序结束*!*读入 grd格式文件中点数和线数及范围开始*SUBROUTINE get_mpoint_and_nline(filename,mpoint,nline,&Xmin,Xmax,Ymin,Ymax) CHARACTER*(*) filenameINTEGER mpoint,nlineREAL Xmin,Xmax,Ymin,YmaxINTEGER unitCALL get_unit(unit)open(unit,file=filename,status=old,err=999)READ(unit,*)READ(unit,*) mpoint,nlineREAD(unit,*) Xmin,XmaxREAD(unit,*) Ymin,Ymax close(unit)RETURNPRINT*,PAUSESTOPENDSUBROUTINE get_mpoint_and_nline!*读入 grd格式文件中点数和线数及范围结束*!*得到扩边端点子程序开始*SUBROUTINE get_expan_num(mpoint,m0,m1,m2,m3,flag)IMPLICIT NONEINTEGER mtemp,factor_m,muINTEGER mpoint,m0,m1,m2,m3INTEGER flagmtemp=mpointfactor_m=1DO WHILE (mod(mtemp,2).eq.0).and.(mtemp.ne.0)mtemp=mtemp/2End doIF (mtemp.eq.1) THENm3=mpoint*2ELSEmu=int(log(float(mpoint)/0.693147+factor_m)m3=2*mu12if(flag=1) thenm3=m3*2 !计算区域为原来倍else endifEND IFm1=1+(m3-mpoint)/2m2=m1+mpoint-1m0=1END SUBROUTINE get_expan_num!*得到扩边端点子程序结束*!*读入 grd格式文件中场值子程序开始 *SUBROUTINE Read_field(filename,m0,m1,m2,m3,n0,n1,n2,n3,G)PARAMETER (vial=1.701411e+38)CHARACTER*(*) filenameINTEGER m0,m1,m2,m3,n0,n1,n2,n3REAL G(m0:m3,n0:n3)INTEGER unitCALL get_unit(unit)G=vialopen(unit,file=filename)DO i=1,5READ(unit,*)ENDDOREAD(unit,*) (G(i,j),i=m1,m2),j=n1,n2)CLOSE(unit)ENDSUBROUTINE Read_field!*读入 grd格式文件中场值子程序结束 *!*扩边子程序集开始*!*扩边子程序调用主子程序开始*SUBROUTINE expan_2D(m0,m1,m2,m3,n0,n1,n2,n3,&Ori_Field,dx,dy,expan_2D_method)INTEGER m0,m1,m2,m3INTEG