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

一 、基本原理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附录：边缘识别程序源代码9 一 、基本原理 1.1 位场边缘识别方法研究进展概述地质目标体边缘是指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。现有重、磁位场边缘识别方法分为数理统计、数值计算和其他三大类。1.2 数值计算类边缘识别方法的研究及计算数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置与真实位置有一定偏差该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。1.2.1 垂向导数（VDR）垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可直接使用，对磁力异常必须转换成磁源重力异常（假重力异常）或化极磁力异常才可以使用。垂向导数方法研究历史较早，方法较成熟，应用频率较高通过我们的试验和研究认为：随着地质体埋深的增大，垂向导数所确定的边缘位置偏离实际位置的距离越大。1.2.2 总水平导数（THDR）总水平导数是利用其极大值位置来确定地质体的边缘位置，适用于重力异常，对磁力异常必须换算成磁源重力异常或化极磁力异常才可以使用。1.2.3 解析信号振幅（ASM）解析信号振幅也是利用极大值位置来确定地质体的边缘位置，适用于重力异常和磁力异常。1.2.4 倾斜角（TA）倾斜角实质上是垂向导数和总水平导数的比值。由于倾斜角为一阶导数的比值，所以能很好地平衡高幅值异常和低幅值异常，起到边缘增强的效果。1.2.5 二阶导数类边缘识别为了增强边缘分辨能力，和克服当总水平导数等于时倾斜角存在“解析奇点”，会使得计算结果不稳定。2、 输入/输出数据格式设计依据上述原理，现对上述各种边缘识别方法实现进行程序设计。2.1 主要变量设计l cmd_file:参数文件名l input_field_filename： 输入位场文件名l output_field_filename：输出转换后位场文件名l expan_2D_method： 2D扩边方法选择l Field_Deriv_method：边缘识别方法选择l num_area：计算区域大小选择l mpoint,nline：点数，线数l m0,m1,m2,m3：扩边后X方向上各端点l n0,n1,n2,n3：扩边后Y方向上各端点l Xmin,Xmax：X坐标最小，最大值l ymin,Ymax：Y坐标最小，最大值l Ori_Field：原始场值，二维数组，m3*n3l Deriv_Field：转换后后场值，二维数组，m3*n3l Deriv_operator_x：x方向转换因子，二维数组，m3*n3l Deriv_operator_y：y方向转换因子，二维数组，m3*n3l Deriv_operator_z：z方向转换因子，二维数组，m3*n32.2位场边缘识别程序数据格式设计l 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_THDRl 原始位场文件，输入文件，grd格式(ASCII)l 延拓后位场文件，输出文件，grd格式(ASCII)DSAAmpoint nlineXmin XmaxYmin YmaxZmin Zmax Z11 Z12 . . l Grd格式如下：3、 总体设计图 3.1 重、磁异常边缘识别程序设计流程图四、测试结果4.1 测试用例（1）观测面上的重力异常存放在“anomaly.grd”中，坐标单位为m。（2）形体边界存放在“rectangle.bln”中， 坐标单位4.2 测试结果图 4.2.2 THDR（总水平导数）结果图图 4.2.3 ASM(解析信号振幅)结果图图 4.2.4 AT(倾斜角)边缘识别结果图图 4.2.6 VDR_THDR结果图(垂向导数总水平导数，余弦扩边方法)图 4.2.5 VDR_VDR(垂向二阶导)结果图图 4.2.7 TA_THDR(倾斜角总水平导数)结果图图 4.2.8 VDR_THDR结果图(垂向导数总水平导数，最小曲率扩边)4.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，可以看出在频率域处理选用扩边方法不同对计算结果有一定影响，其中最小曲率扩边计算效果为佳。五、结论与建议5.1 结论此方法在一定程度上是有效的，且从实验结果可以垂向导数总水平导数和垂向二阶导识别的效果较好。5.2 建议没有学好位场边缘识别方法没有可以能提出的建议附录：边缘识别程序源代码!*! PROGRAM: Fre_PoteField_Deriv! PURPOSE: Frequency Potential field Derivative! Author : gesang! Data: 2013-11!* PROGRAM Fre_PoteField_Deriv implicit none CHARACTER *(80) input_field_filename CHARACTER *(80) output_field_filename INTEGER expan_2D_method,num_area,Field_Deriv_method INTEGER mpoint,nline INTEGER m0,m1,m2,m3 INTEGER n0,n1,n2,n3 REAL Xmin,Xmax REAL ymin,Ymax REAL dx,dy REAL,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) CALL 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_Field Field_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_filename CHARACTER *(*) output_field_filename INTEGER expan_2D_method,Field_Deriv_method,num_area INTEGER unit CHARACTER*80 s CALL get_unit(unit) OPEN(unit,file='cmd.txt',status='old') READ(unit,*) s ,input_field_filename READ(unit,*) s ,output_field_filenameREAD(unit,*) s ,num_area READ(unit,*) s ,expan_2D_method READ(unit,*) s ,Field_Deriv_method CLOSE(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)RETURN PRINT*,' PAUSE STOPENDSUBROUTINE get_mpoint_and_nline!*读入grd格式文件中点数和线数及范围结束*!*得到扩边端点子程序开始*SUBROUTINE get_expan_num(mpoint,m0,m1,m2,m3,flag) IMPLICIT NONE INTEGER 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) THEN m3=mpoint*2ELSE mu=int(log(float(mpoint)/0.693147+factor_m) m3=2*muif(flag=1) then m3=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,5 READ(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,m3 INTEGER n0,n1,n2,n3 REAL Ori_Field(m0:m3,n0:n3) INTEGER expan_2D_method REAL dx,dy IF(expan_2D_method<=3) THEN CALL expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,& Ori_Field,expan_2D_method) ELSE CALL MinCurv(Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,& dx,dy,expan_2D_method) END IFEND SUBROUTINE expan_2D!*扩边子程序调用主子程序结束*!*2D最小曲率扩边子程序开始*SUBROUTINE MinCurv(Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,dx,dy,flag) INTEGER m0,m1,m2,m3,n0,n1,n2,n3 REAL Ori_Field(m0:m3,n0:n3) REAL dx,dy INTEGER num !空白点点数 INTEGER ,ALLOCATABLE:P_POINT(:),P_LINE(:) INTEGER flag REAL eigval eigval=1.701411e+38 CALL Blank_Point_number(Ori_Field,m3,n3,num,eigval) ALLOCATE(P_POINT(1:num)ALLOCATE(P_LINE(1:num) CALL Blank_Point_position(Ori_Field,m3,n3,num,eigval,P_POINT,P_LINE) CALL expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,flag-3) !余弦扩边作为空白部分初值及端点值 !*2D网格数据无约束最小曲率迭代 CALL MinCurV_2D_net_sub(1E-4,10000,dx,dy,m3,n3,& Ori_Field,num,P_POINT,P_LINE) END SUBROUTINE MinCurv!*2D最小曲率扩边子程序结束*!*获得数据中的空白点数开始*SUBROUTINE Blank_Point_number(FIELD,mpoint,line,bpn,eigval)IMPLICIT NONEINTEGER mpoint,line,bpn,i,jREAL FIELD(1:mpoint,1:line),eigvalbpn=0DO i=1,mpoint,1 DO j=1,line,1 IF(FIELD(i,j)>=0.9*eigval) bpn=bpn+1 ENDDOENDDOENDSUBROUTINE!*获得数据中的空白点数结束*!*获得数据中的空白点位置开始*SUBROUTINE Blank_Point_position(FIELD,mpoint,line,bpn,eigval,P_POINT,P_LINE)IMPLICIT NONEINTEGER mpoint,line,bpn,i,j,kINTEGER P_POINT(1:bpn),P_LINE(1:bpn)REAL FIELD(1:mpoint,1:line),eigvalk=1DO i=1,mpoint,1 DO j=1,line,1 IF(FIELD(i,j)>=0.9*eigval.AND.j<bpn) THEN P_POINT(k)=i P_LINE(k)=j k=k+1 ENDIF ENDDOENDDOENDSUBROUTINE!*获得数据中的空白点位置结束*!*2D无约束最小曲率迭代边界处理开始*SUBROUTINE MinCurV_2D_boundary(U,mpoint,line,alph)REAL U(1-2:mpoint+2,1-2:line+2)DO j=1,line,1 U(1-1,j)=2*U(1,j)-U(1+1,j) U(mpoint+1,j)=2*U(mpoint,j)-U(mpoint-1,j)ENDDODO i=1,mpoint,1 U(i,1-1)=2*U(i,1)-U(i,1+1) U(i,line+1)=2*U(i,line)-U(i,mpoint-1)ENDDOU(1-1,1-1)=U(1+1,1-1)+U(1-1,1+1)-U(1+1,1+1)U(mpoint+1,1-1)=U(mpoint+1,1+1)+U(mpoint-1,1-1)-U(mpoint-1,1+1)U(1-1,line+1)=U(1+1,line+1)+U(1-1,line-1)-U(1+1,line-1)U(mpoint+1,line+1)=U(mpoint+1,line-1)+U(mpoint-1,line+1)-U(mpoint-1,line-1)i=1alph1=alph*2alph2=-2*(1+alph1)DO j=1,line pij=alph1*(U(i+1,j+1)-U(i-1,j+1)+U(i+1,j-1)-U(i-1,j-1)+alph2*(U(i+1,j)-U(i-1,j) U(i-2,j)=U(i+2,j)+pijENDDOi=mpointDO j=1,line pij=alph1*(U(i+1,j+1)-U(i-1,j+1)+U(i+1,j-1)-U(i-1,j-1)+alph2*(U(i+1,j)-U(i-1,j) U(i+2,j)=U(i-2,j)-pijENDDOj=1beta=1./alphbeta1=beta*betabeta2=-2*(1+beta1)DO i=1,mpoint,1 qij=beta1*(U(i+1,j+1)-U(i+1,j-1)+U(i-1,j+1)-U(i-1,j-1)+beta2*(U(i,j+1)-U(i,j-1) U(i,j-2)=U(i,j+2)+qijENDDOj=lineDO i=1,mpoint,1 qij=beta1*(U(i+1,j+1)-U(i+1,j-1)+U(i-1,j+1)-U(i-1,j-1)+beta2*(U(i,j+1)-U(i,j-1) U(i,j+2)=U(i,j-2)-qijENDDOENDSUBROUTINE MinCurV_2D_boundary!*2D无约束最小曲率迭代边界处理结束*!*2D网格数据无约束最小曲率迭代开始*SUBROUTINE MinCurV_2D_net_sub(eps_abs,iteration_max,dx,dy,& mpoint,line,FIELD,number_eigval,M_eigval,N_eigval)REAL eps_abs,dx,dyINTEGER iteration_max,mpoint,line,number_eigvalREAL FIELD(1:mpoint,1:line)INTEGER M_eigval(1:number_eigval),N_eigval(1:number_eigval)REAL,ALLOCATABLE:U(:,:)ALLOCATE(U(1-2:mpoint+2,1-2:line+2)DO j=1,line,1 DO i=1,mpoint,1 U(i,j)=FIELD(i,j) ENDDOENDDOeps_error=2*ABS(eps_abs)iteration=0alph=dx/dyalph0=-1./(2*(3+4*alph*2+3*alph*2) alph1=alph*4alph2=2*alph*2alph3=-4*(1+alph*2)alph4=-4*(1+alph*2)*alph*2DO WHILE(eps_error>eps_abs).and.(iteration<iteration_max) eps_error=0. CALL MinCurV_2D_boundary(U,mpoint,line,alph) DO k=1,number_eigval,1 i=M_eigval(k) j=N_eigval(k) temp=(U(i+2,j)+U(i-2,j)+alph1*(U(i,j+2)+U(i,j-2)+ & alph2*(U(i+1,j+1)+U(i-1,j+1)+U(i+1,j-1)+U(i-1,j-1)+ & alph3*(U(i+1,j)+U(i-1,j)+alph4*(U(i,j+1)+U(i,j-1) temp=temp*alph0 eps_error=MAX(ABS(temp-U(i,j),eps_error) U(i,j)=temp ENDDO iteration=iteration+1ENDDOprint*,iteration,eps_errorDO j=1,line,1 DO i=1,mpoint,1 FIELD(i,j)=U(i,j) ENDDOENDDOdeALLOCATE(U,STAT=ierr)ENDSUBROUTINE!*2D网格数据无约束最小曲率迭代结束*!*2D余弦扩边开始*SUBROUTINE expan_cos2(m1,m2,m3,m4,n1,n2,n3,n4,G,flag) PARAMETER (pi=3.141592654) INTEGER flag INTEGER m1,m2,m3,m4 INTEGER n1,n2,n3,n4 REAL G(m1:m4,n1:n4) REAL G1(m1:m4,n1:n4) REAL G2(m1:m4,n1:n4) REAL sum INTEGER num sum=0.0 num=0 IF(flag=1) then !扩边端点值取边界平均值 DO i=m2,m3 sum=sum+G(i,n2)+G(i,n3) num=num+2 END DO DO i=n2+1,n3-1 sum=sum+G(m2,i)+G(m3,i) num=num+2 END DO DO i=m1,m4 G(i,n1)=sum/num G(i,n4)=sum/num END DO DO i=n1+1,n4-1 G(m1,i)=sum/num G(m4,i)=sum/num END DO ELSE IF(flag=2) then !扩边端点值取全部数据平均值 DO i=m2,m3 DO j=n2,n3 sum=sum+G(i,j) num=num+1 END DO END DO DO i=m1,m4 G(i,n1)=sum/num G(i,n4)=sum/num END DO DO i=n1+1,n4-1 G(m1,i)=sum/num G(m4,i)=sum/num END DO ELSE IF(flag=3)