实验四 位场边缘识别程序设计实验讲解Word下载.docx

实验四 位场边缘识别程序设计实验讲解Word下载.docx

《实验四 位场边缘识别程序设计实验讲解Word下载.docx》由会员分享，可在线阅读，更多相关《实验四 位场边缘识别程序设计实验讲解Word下载.docx（28页珍藏版）》请在冰豆网上搜索。

垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用

（1.1）

（2）解析信号振幅：

解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常

（1.2）

（3）总水平导数（THDR）

（1.3）

2、输入/输出数据格式设计

依据上述原理，现在对上述各种边缘识别方法进行程序设计。

2.1输入数据格式设计

本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。

例如：

DSAA

201201

-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格式，均为实型变量。

-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.2重力异常原始图

图4.3垂向导数（VDR）

图4.5总水平导数（THDR）

图4.4解析信号振幅（ASM）

图4.3ASM（解析信号振幅）

分析：

由图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（:

）

Px_re（:

）,Px_im（:

）,Py_re（:

）,Py_im（:

）,Pz_re（:

）,Pz_im（:

field_vdr（:

）,field_thdr（:

）,field_asm（:

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:

n））

allocate（Px_re（1:

n）,Px_im（1:

n）,Py_re（1:

n）,Py_im（1:

n）,Pz_re（1:

n）,Pz_im（1:

allocate（field_vdr（1:

n）,field_thdr（1:

n）,field_asm（1:

allocate（U（1:

m）,V（1:

n）,W（1:

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）

call