整理实验二边缘识别实验报告正文.docx
《整理实验二边缘识别实验报告正文.docx》由会员分享，可在线阅读，更多相关《整理实验二边缘识别实验报告正文.docx（43页珍藏版）》请在冰豆网上搜索。
整理实验二边缘识别实验报告正文
（1）非煤矿矿山的建设项目（注：
对煤矿建设项目有单独特别规定）；
一、安全评价
4.将环境影响价值纳入项目的经济分析
2）应用环境质量标准时，应结合环境功能区和环境保护目标进行分级。
大纲要求
一、环境影响评价的发展与管理体系、相关法律法规体系和技术导则的应用
（一）规划环境影响评价的适用范围和责任主体
[例题-2005年真题]《中华人民共和国环境影响评价法》规定，建设项目可能造成轻度环境影响的，应当编制（ ）。
大纲要求
（一）安全评价的内涵
一、基本原理
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主要变量设计
●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*n3
2.2位场边缘识别程序数据格式设计
●cmd文件，输入参数文件，格式如下：
（cmd.txt）
input_field_filenameanomaly.grd
output_field_filenameVDR_THDR_anomaly.grd
num_area0
expan_2D__method1
Field_Deriv_method7
说明：
num_area:
取0：
表示一般扩边区域
取1：
表示计算区域相对于一般扩边再放大一倍
expan_2D__method:
取1：
表示余弦扩边，端点值取边界平均值
取2：
表示余弦扩边，端点值取全场区平均值
取3：
表示余弦扩边，端点值取常数0
取4：
表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取边界平均值
取5：
表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取全场区平均值
取6：
表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取常数0
Field_Deriv_method:
取1：
计算垂向导数VDR
取2:
计算总水平导数THDR
取3:
计算解析信号振幅ASM
取4:
计算倾斜角TA
取5:
计算垂向二阶导数VDR_VDR
取6:
计算倾斜角总水平导数TA_THDR
取7:
计算垂向导数总水平导数VDR_THDR
●原始位场文件，输入文件，grd格式（ASCII）
●延拓后位场文件，输出文件，grd格式（ASCII）
●
DSAA
mpointnline
XminXmax
YminYmax
ZminZmax
Z11Z12
......
Grd格式如下：
3、总体设计
四、测试结果
4.1测试用例
（1）观测面上的重力异常存放在“anomaly.grd”中，坐标单位为m。
（2）形体边界存放在“rectangle.bln”中，坐标单位
4.2测试结果
图4.2.3ASM（解析信号振幅）结果图
图4.2.5VDR_VDR（垂向二阶导）结果图
图4.2.4AT（倾斜角）边缘识别结果图
图4.2.6VDR_THDR结果图
（垂向导数总水平导数，余弦扩边方法）
图4.2.7TA_THDR（倾斜角总水平导数）结果图
图4.2.8VDR_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:
FrequencyPotentialfieldDerivative
!
Author:
gesang
!
Data:
2013-11
!
**********************************************************
PROGRAMFre_PoteField_Deriv
implicitnone
CHARACTER*（80）input_field_filename
CHARACTER*（80）output_field_filename
INTEGERexpan_2D__method,num_area,Field_Deriv_method
INTEGERmpoint,nline
INTEGERm0,m1,m2,m3
INTEGERn0,n1,n2,n3
REALXmin,Xmax
REALymin,Ymax
REALdx,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:
!
读取文件参数
CALLread_cmdinfo（input_field_filename,output_field_filename,&
expan_2D__method,Field_Deriv_method,num_area）
!
读取点数线数及测区范围
CALLget_mpoint_and_nline（input_field_filename,mpoint,nline,Xmin,Xmax,&
Ymin,Ymax）
!
计算扩边端点
CALLget_expan_num（mpoint,m0,m1,m2,m3,num_area）
CALLget_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））
!
从网格化文件读入初始场值
CALLRead_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）
CALLexpan_2D（m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,dx,dy,&
expan_2D__method）!
expan_2D__method:
扩边方法选择
!
计算导数因子
CALLDeriv_operator_sub（dx,dy,m3,n3,Deriv_operator_x,&
Deriv_operator_y,Deriv_operator_z）
!
求导运算
Field_Real=Ori_Field
Field_Image=0.0
CALLfield_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:
!
输出求导后位场文件
CALLoutput_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）
ENDPROGRAMFre_PoteField_Deriv
!
********************************开始****************
!
*********得到输入参数子程序开始***********
SUBROUTINEread_cmdinfo（input_field_filename,output_field_filename,&
expan_2D__method,Field_Deriv_method,num_area）
CHARACTER*（*）input_field_filename
CHARACTER*（*）output_field_filename
INTEGERexpan_2D__method,Field_Deriv_method,num_area
INTEGERunit
CHARACTER*80s
CALLget_unit（unit）
OPEN（unit,file='cmd.txt',status='old'）
READ（unit,*）s,input_field_filename
READ（unit,*）s,output_field_filename
READ（unit,*）s,num_area
READ（unit,*）s,expan_2D__method
READ（unit,*）s,Field_Deriv_method
CLOSE（unit）
ENDSUBROUTINEread_cmdinfo
!
******得到输入参数子程序结束******
!
****读入grd格式文件中点数和线数及范围开始*****
SUBROUTINEget_mpoint_and_nline（filename,mpoint,nline,&
Xmin,Xmax,Ymin,Ymax）
CHARACTER*（*）filename
INTEGERmpoint,nline
REALXmin,Xmax,Ymin,Ymax
INTEGERunit
CALLget_unit（unit）
open（unit,file=filename,status='old',err=999）
READ（unit,*）
READ（unit,*）mpoint,nline
READ（unit,*）Xmin,Xmax
READ（unit,*）Ymin,Ymax
close（unit）
RETURN
PRINT*,'
PAUSE
STOP
ENDSUBROUTINEget_mpoint_and_nline
!
****读入grd格式文件中点数和线数及范围结束********
!
**********得到扩边端点子程序开始************
SUBROUTINEget_expan_num（mpoint,m0,m1,m2,m3,flag）
IMPLICITNONE
INTEGERmtemp,factor_m,mu
INTEGERmpoint,m0,m1,m2,m3
INTEGERflag
mtemp=mpoint
factor_m=1
DOWHILE（（mod（mtemp,2）.eq.0）.and.（mtemp.ne.0））
mtemp=mtemp/2
Enddo
IF（mtemp.eq.1）THEN
m3=mpoint*2
ELSE
mu=int（log（float（mpoint））/0.693147+factor_m）
m3=2**mu
if（flag==1）then
m3=m3*2!
计算区域为原来倍
else
endif
ENDIF
m1=1+（m3-mpoint）/2
m2=m1+mpoint-1
m0=1
ENDSUBROUTINEget_expan_num
!
*******得到扩边端点子程序结束*************
!
*****读入grd格式文件中场值子程序开始********
SUBROUTINERead_field（filename,m0,m1,m2,m3,n0,n1,n2,n3,G）
PARAMETER（vial=1.701411e+38）
CHARACTER*（*）filename
INTEGERm0,m1,m2,m3,n0,n1,n2,n3
REALG（m0:
m3,n0:
n3）
INTEGERunit
CALLget_unit（unit）
G=vial
open（unit,file=filename）
DOi=1,5
READ（unit,*）
ENDDO
READ（unit,*）（（G（i,j）,i=m1,m2）,j=n1,n2）
CLOSE（unit）
ENDSUBROUTINERead_field
!
*****读入grd格式文件中场值子程序结束********
!
**************扩边子程序集开始********************************
!
*****扩边子程序调用主子程序开始*****
SUBROUTINEexpan_2D（m0,m1,m2,m3,n0,n1,n2,n3,&
Ori_Field,dx,dy,expan_2D__method）
INTEGERm0,m1,m2,m3
INTEGERn0,n1,n2,n3
REALOri_Field（m0:
m3,n0:
n3）
INTEGERexpan_2D__method
REALdx,dy
IF（expan_2D__method<=3）THEN
CALLexpan_cos2（m0,m1,m2,m3,n0,n1,n2,n3,&
Ori_Field,expan_2D__method）
ELSE
CALLMinCurv（Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,&
dx,dy,expan_2D__method）
ENDIF
ENDSUBROUTINEexpan_2D
!
*****扩边子程序调用主子程序结束*****
!
*********2D最小曲率扩边子程序开始***********
SUBROUTINEMinCurv（Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,dx,dy,flag）
INTEGERm0,m1,m2,m3,n0,n1,n2,n3
REALOri_Field（m0:
m3,n0:
n3）
REALdx,dy
INTEGERnum!
空白点点数
INTEGER,ALLOCATABLE:
:
P_POINT（:
）,P_LINE（:
）
INTEGERflag
REALeigval
eigval=1.701411e+38
CALLBlank_Point_number（Ori_Field,m3,n3,num,eigval）
ALLOCATE（P_POINT（1:
num））
ALLOCATE（P_LINE（1:
num））
CALLBlank_Point_position（Ori_Field,m3,n3,num,eigval,P_POINT,P_LINE）
CALLexpan_cos2（m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,flag-3）
!
余弦扩边作为空白部分初值及端点值
!
****2D网格数据无约束最小曲率迭代
CALLMinCurV_2D_net_sub（1E-4,10000,dx,dy,m3,n3,&
Ori_Field,num,P_POINT,P_LINE）
ENDSUBROUTINEMinCurv
!
*********2D最小曲率扩边子程序结束***********
!
****获得数据中的空白点数开始*****************
SUBROUTINEBlank_Point_number（FIELD,mpoint,line,bpn,eigval）
IMPLICITNONE
INTEGERmpoint,line,bpn,i,j
REALFIELD（1:
mpoint,1:
line）,eigval
bpn=0
DOi=1,mpoint,1
DOj=1,line,1
IF（FIELD（i,j）>=0.9*eigval）bpn=bpn+1
ENDDO
ENDDO
ENDSUBROUTINE
!
****获得数据中的空白点数结束*****************
!
****获得数据中的空白点位置开始********************
SUBROUTINEBlank_Point_position（FIELD,mpoint,line,bpn,eigval,P_POINT,P_LINE）
IMPLICITNONE
INTEGERmpoint,line,bpn,i,j,k
INTEGERP_POINT（1:
bpn）,P_LINE（1:
bpn）
REALFIELD（1:
mpoint,1:
line）,eigval
k=1
DOi=1,mpoint,1
DOj=1,line,1
IF（FIELD（i,j）>=0.9*eigval.AND.jP_POINT（k）=i
P_LINE（k）=j
k=k+1
ENDIF
ENDDO
ENDDO
ENDSUBROUTINE
!
****获得数据中的空白点位置结束********************
!
****2D无约束最小
升级会员 