实验四 位场边缘识别程序设计实验.docx
《实验四 位场边缘识别程序设计实验.docx》由会员分享，可在线阅读，更多相关《实验四 位场边缘识别程序设计实验.docx（30页珍藏版）》请在冰点文库上搜索。
实验四位场边缘识别程序设计实验
《重磁资料处理与解释》实验四
位场边缘识别程序设计实验
专业名称：
地球物理学
学生姓名：
学生学号：
指导老师：
王万银、纪新林、纪晓琳、邱世灿
提交日期：
2016-1-3
1、基本原理
地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。
现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。
数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。
在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值。
故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。
该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。
因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。
（1）垂向导数：
垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用
（1.1）
（2）解析信号振幅：
解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常
（1.2）
（3）总水平导数（THDR）
（1.3）
2、输入/输出数据格式设计
依据上述原理，现在对上述各种边缘识别方法进行程序设计。
2.1输入数据格式设计
本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。
例如：
DSAA
201201
-1000.0000001000.000000
-1000.0000001000.000000
5.549671E-0123.539846
5.549671E-015.634658E-015.721339E-015.808522E-01
5.897312E-015.987253E-016.078691E-016.171604E-01
…
2.2输出数据格式设计
计算结果输出数据格式与输入格式对应，格式为.GRD格式，均为实型变量。
例如：
DSAA
201201
-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，字符串变量，长度不超过80
factor_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效果最好。
附录：
边缘识别程序源代码
!
****************************************************************************************************
!
程序功能：
实现频率域位场导数运算进行边缘识别
!
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:
场值的解析信号振幅的结果
!
****************************************************************************************************
programdeviation
parameter（eigval=3.701411*1e5）
character*（80）cmdfile
character*80input_file,output_file_vdr,output_file_thdr,output_file_asm
real,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_line
integerm,n,m1,m2,n1,n2
realfactor_m
realxmin,xmax,ymin,ymax,dx,dy
cmdfile='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_vdr：
对场值作水平导数处理后的数据文件
!
output_file_thdr:
对场值作垂向导数处理后的数据文件
!
output_file_asm:
对场值作总导数处理后的数据文件
!
factor_m：
扩边因子
!
****************************************************************************
Subroutineread_cmd（cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm）
implicitnone
character*80str
character*（*）cmdfile
character*（*）input_file,output_file_vdr,output_file_thdr,output_file_asm
realfactor_m
open（10,file=cmdfile,status='old'）
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_m
close（10）
endSubroutineread_cmd
!
***************************************************************************
!
子程序：
read_grd
!
功能：
从原始观测.grd文件中读取相关参数
!
输入参数说明：
!
filename_obser：
输入文件名
!
输出参数说明：
!
N_point,N_line：
点数、线数
!
x_min,x_max：
x的最小值和最大值
!
y_min,y_max：
y的最小值和最大值
!
***************************************************************************
subroutineread_grd（input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax）
implicitnone
character*（*）input_file
integerN_point,N_line
realXmin,Xmax,Ymin,Ymax
open（10,file=input_file,status='old'）
Read（10,*）
Read（10,*）N_line,N_point
Read（10,*）Xmin,Xmax
Read（10,*）Ymin,Ymax
Close（10）
endsubroutineread_grd
!
**************************************************************************
!
子程序：
calculate_mn
!
功能：
确定扩边数据点号位置
!
输入参数说明：
!
factor_m:
扩边比例因子（>1.0）
!
a,b：
点数、线数
!
输出参数说明：
!
m1,m2:
x方向实际数据起点位置和终点位置点号
!
m:
扩边后数据终点位置点号（起点位置点号为1）
!
n1,n2:
y方向实际数据起点位置和终点位置点号
!
n:
扩边后数据终点位置点号（起点位置点号为1）
!
**************************************************************************
subroutinecalculate_mn（a,b,m,n,m1,m2,n1,n2,factor_m）
implicitnone
integera,b,m,n,m1,m2,n1,n2
integermtemp,mu,nu
realfactor_m
mtemp=a
DOWHILE（（mod（mtemp,2）.eq.0）.and.（mtemp.ne.0））
mtemp=mtemp/2
Enddo
IF（mtemp.eq.1）THEN
m=a*2
ELSE
mu=int（log（float（a））/0.693147+factor_m）
m=2**mu
ENDIF
m1=1+（m-a）/2
m2=m1+a-1
write（*,*）m,a
pause
mtemp=b
DOWHILE（（mod（mtemp,2）.eq.0）.and.（mtemp.ne.0））
mtemp=mtemp/2
Enddo
IF（mtemp.eq.1）THEN
n=b*2
ELSE
nu=int（log（float（b））/0.693147+factor_m）
n=2**nu
ENDIF
n1=1+（n-b）/2
n2=n1+b-1
write（*,*）m1,m2,n1,n2,m,n
pause
endsubroutinecalculate_mn
!
*************************************************************************
!
子程序：
INPUT_GRD
!
功能：
读取grd文件中的数据
!
输入参数说明：
!
filename_obser：
输入文件名
!
m1,m2:
x方向实际数据起点位置和终点位置点号
!
m:
扩边后数据终点位置点号（起点位置点号为1）
!
n1,n2:
y方向实际数据起点位置和终点位置点号
!
n:
扩边后数据终点位置点号（起点位置点号为1）
!
输出参数说明：
!
A：
存放输出数据的二维数组名
!
*************************************************************************
SUBROUTINEINPUT_GRD（A,input_file,m1,m2,n1,n2,m,n）
character*（*）input_file
integerm1,m2,n1,n2,m,n
realxmin,xmax,ymin,ymax
realA（1:
m,1:
n）
reali,j,k
doj=1,n,1
doi=1,m,1
A（i,j）=3.701411*1e10
enddo
enddo
Open（20,file=input_file,status='old'）
read（20,*）
read（20,*）
read（20,*）xmin,xmax
read（20,*）ymin,ymax
read（20,*）
read（20,*）（（A（i,j）,i=m1,m2）,j=n1,n2）
Close（20）
ENDSUBROUTINEINPUT_GRD
!
*************************************************************************
!
子程序：
expand_cos_2D
!
功能：
二维扩边子程序并为信号虚部赋值
!
输入参数说明：
!
m1,m2:
x方向实际数据起点位置和终点位置点号
!
m:
扩边后数据终点位置点号（起点位置点号为1）
!
n1,n2:
y方向实际数据起点位置和终点位置点号
!
n:
扩边后数据终点位置点号（起点位置点号为1）
!
Ur:
初始观测面信号的实部
!
Ui：
初始观测面信号的虚部
!
输出参数说明：
!
Ur:
初始观测面信号的实部
!
Ui：
初始观测面信号的虚部
!
*************************************************************************
subroutineexpand_cos_2D（m1,m2,m,n1,n2,n,Ur,Ui）
implicitnone
integerm,n,m1,m2,n1,n2
realUr（1:
m,1:
n）,Ui（1:
m,1:
n）
real,allocatable:
:
u（:
）,r（:
）
integerj,i,k
allocate（u（1:
m））
doj=n1,n2,1
doi=1,m,1
u（i）=Ur（i,j）
enddo
callexpand_cos_1d（1,m1,m2,m,u
（1））
doi=1,m,1
Ur（i,j）=u（i）
enddo
enddo
deallocate（u）
allocate（r（1:
n））
doi=1,m,1
doj=1,n,1
r（j）=Ur（i,j）
enddo
callexpand_cos_1d（1,n1,n2,n,r
（1））
doj=1,n,1
Ur（i,j）=r（j）
enddo
enddo
deallocate（r）
doi=1,m
doj=1,n
Ui（i,j）=0
enddo
enddo
endsubroutineexpand_cos_2D
!
**************************************************************************
!
子程序：
expand_cos_1d
!
功能：
一维扩边子程序
!
输入参数说明：
!
n0,n3：
扩边后数据起点位置和终点位置
!
n1,n2：
实际数据起点位置和终点位置
!
feild（i）,（i=n1,n1+1,...,n2）：
实际数据
!
输出参数说明：
!
field（i）,（i=n0,...,n1-1）：
扩边后左边的数据
!
field（i）,（i=n2+1,...,n3）：
扩边后右边的数据
!
**************************************************************************
Subroutineexpand_cos_1d（n0,n1,n2,n3,Field）
RealField（n0:
n3）
pi=3.141592654
Field（n0）=（Field（n1）+Field（n2））/2.0
Field（n3）=Field（n0）
i1=n0
i2=n1
DOi=i1,i2-1,1
Field（i）=Field（i1）+cos（pi/2.0*（i2-i）/（i2-i1））*（Field（i2）-Field（i1））
Enddo
i1=n2
i2=n3
升级会员 