重磁实验三.docx

资源描述

重磁实验三.docx

《重磁实验三.docx》由会员分享，可在线阅读，更多相关《重磁实验三.docx（22页珍藏版）》请在冰豆网上搜索。

重磁实验三.docx

重磁实验三

《重、磁资料处理与解释》上机实验报告

实验三：

位场边缘识别程序设计实验

姓名：

学号：

专业：

地球物理学

指导教师：

王万银、纪晓琳

完成时间：

2017.1.10

1基本原理3

1.1水平导数因子和垂向导数因子3

1.2垂向导数3

1.3总水平导数3

1.4解析信号振幅3

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

2.1输入/输出数据文件名3

2.2重要变量名4

3总体设计4

4测试结果5

5结论及建议7

附录：

源程序代码8

1基本原理

1.1水平导数因子和垂向导数因子

1.2垂向导数

1.3总水平导数

1.4解析信号振幅

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

2.1输入/输出数据文件名

输入数据和输出数据文件名均保存在“parameter.txt”中。

第一行为输入的重力异常文件；第二行～第四行依次为输出的垂向一阶导数、总水平导数以及解析信号振幅文件；第五行为输入的扩边比例因子。

gravity.grd

VDR.grd

THDR.grd

ASM.grd

1.5

2.2重要变量名

filename_Field:

重力异常文件

filename_VDR:

垂向一阶导数文件

filename_THDR:

总水平导数文件

filename_ASM:

解析信号振幅文件

Field（m,n）:

重力异常

VDR:

垂向一阶导数

THDR:

总水平导数

ASM:

解析信号振幅

Xreal,Ximage:

x方向导数实部、虚部

Yreal,Yimage:

y方向导数实部、虚部

Zreal,Zimage:

z方向导数实部、虚部

factor_x:

扩边比例因子（>1.0）

point:

点数

line:

线数

扩边以后总点数（满足2的幂次方）

扩边以后总线数（满足2的幂次方）

3总体设计

输入有关参数：

filename_Field,filename_VDR,filename_THDR,filename_ASM,factor_x

从文件输入有关参数：

point,line

计算有关参数：

m,m1,m2,n,n1,n2

从文件输入有关数据：

xmin,xmax,ymin,ymax,dx,dy,Field

进行扩边处理

傅里叶正变换

计算导数因子

进行乘积运算

傅里叶反变换

输出计算结果：

VDR,THDR,ASM

4测试结果

5结论及建议

（1）由重力异常垂向导数（图2）可以看出，VDR方法的零值线可大致识别异常体的边界；由重力异常解析信号振幅（图4）可以看出，ASM方法的极大值可大致识别异常体的边界。

以上两种方法的精度都不是很高，特征值识别的异常体边界与实际边界都存在偏差。

（2）由重力异常总水平导数（图3）可以看出，THDR方法的极大值可较准确地识别异常体的边界。

对比图2到图3，THDR方法对异常体边界的识别效果最好。

（3）本次程序设计过程中的一些认识：

由于有了实验二位场延拓的基础，本次位场边缘识别实验的程序编写就得心应手了许多，只需把实验二中延拓因子的乘积运算替换为导数因子的乘积运算即可。

要注意的一点是乘以x,y方向导数的时候，由于导数因子为虚数

，因此重力异常x,y方向导数的实部为负的重力异常的虚部乘以导数因子的虚部

，重力异常x,y方向导数的虚部为重力异常的实部乘以导数因子的虚部

。

对VDR、THDR及ASM与重力异常x,y,z方向导数的关系没有理解，导致程序编写有问题；x,y方向对应的圆频率U,V是一维数组，在主程序中定义成了二维数组导致出错；程序中相似的变量有很多，一个字母写错就会造成变量混乱，导致输出结果差之千里。

附录：

源程序代码

******************************************************************c

功能：

位场边缘识别

文件参数说明：

filename_Field:

重力异常文件

filename_VDR:

垂向一阶导数文件

filename_THDR:

总水平导数文件

filename_ASM:

解析信号振幅文件

输入参数说明：

factor_x:

扩边比例因子（>1.0）

point:

点数

line:

线数

xmin,xmax:

点坐标的最大值，最小值

ymin,ymax:

线坐标的最大值，最小值

dx:

x方向点距

dy:

y方向线距

扩边参数说明：

m1,m2:

实际数据起点位置和终点位置（m2-m1+1=point）

扩边以后总点数（满足2的幂次方）

n1,n2:

实际数据起点位置和终点位置（n2-n1+1=line）

扩边以后总线数（满足2的幂次方）

输入数据说明：

Field（m,n）:

重力异常

Freal（m,n）:

重力异常实部

Fimage（m,n）:

重力异常虚部

导数参数说明：

U（m）:

x方向对应的圆频率

V（n）:

y方向对应的圆频率

W（m,n）:

径向圆频率

Xreal,Ximage:

x方向导数实部、虚部

Yreal,Yimage:

y方向导数实部、虚部

Zreal,Zimage:

z方向导数实部、虚部

输出数据说明：

VDR:

垂向一阶导数

THDR:

总水平导数

ASM:

解析信号振幅

------------------------------------------------------------------c

ProgramPotenEdge

implicitnone

character*（80）filename_Field,filename_VDR,filename_THDR,filename_ASM

realfactor_x

integerpoint,line

realxmin,xmax,ymin,ymax,dx,dy

integerm,m1,m2,n,n1,n2

real,allocatable:

Field（:

）,Freal（:

）,Fimage（:

）

real,allocatable:

VDR（:

）,THDR（:

）,AASM（:

）

real,allocatable:

Xreal（:

）,Yreal（:

）,Zreal（:

）

real,allocatable:

Ximage（:

）,Yimage（:

）,Zimage（:

）

real,allocatable:

U（:

）,V（:

）,W（:

）

读取文件名及参数

callGet_parameter（filename_Field,filename_VDR,filename_THDR,filename_ASM,factor_x）

读取点线数

callGet_point_line（filename_Field,point,line）

计算扩边数据

callGET_M1M2m（point,factor_x,m,m1,m2）

callGET_M1M2m（line,factor_x,n,n1,n2）

开辟数组

allocate（Field（1:

m,1:

n）,Freal（1:

m,1:

n）,Fimage（1:

m,1:

n））

allocate（VDR（1:

m,1:

n）,THDR（1:

m,1:

n）,AASM（1:

m,1:

n））

allocate（Xreal（1:

m,1:

n）,Yreal（1:

m,1:

n）,Zreal（1:

m,1:

n））

allocate（Ximage（1:

m,1:

n）,Yimage（1:

m,1:

n）,Zimage（1:

m,1:

n））

allocate（U（1:

m）,V（1:

n）,W（1:

m,1:

n））

读取重力异常

callInput_Field（filename_Field,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,dx,dy,Field）

二维扩边

callExpend_cos_2D（m,m1,m2,n,n1,n2,Field）

2DFFT正变换

Freal=Field

Fimage=0.0

callFFT2（Freal,Fimage,m,n,2）

计算二维圆频率

callWAVE2（m,n,dx,dy,U,V,W）

乘以导数因子

callProduct_Factor（m,n,U,V,W,Freal,Fimage,Xreal,Ximage,Yreal,Yimage,Zreal,Zimage）

2DFFT反变换

callFFT2（Xreal,Ximage,m,n,1）

callFFT2（Yreal,Yimage,m,n,1）

callFFT2（Zreal,Zimage,m,n,1）

输出垂向一阶导数、总水平导数和解析信号振幅

VDR=Zreal

callOutput_grd（filename_VDR,point,line,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,VDR）

THDR=sqrt（Xreal**2+Yreal**2）

callOutput_grd（filename_THDR,point,line,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,THDR）

AASM=sqrt（VDR**2+THDR**2）

callOutput_grd（filename_ASM,point,line,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,AASM）

EndProgramPotenEdge

读取文件名及参数子程序

SubroutineGet_parameter（filename_Field,filename_VDR,filename_THDR,filename_ASM,factor_x）

implicitnone

character*（*）filename_Field,filename_VDR,filename_THDR,filename_ASM

realfactor_x

open（10,file='parameter.txt',status='old'）

read（10,*）filename_Field

read（10,*）filename_VDR

read（10,*）filename_THDR

read（10,*）filename_ASM

read（10,*）factor_x

close（10）

EndSubroutineGet_parameter

读取点线数子程序

SubroutineGet_point_line（filename_Field,point,line）

implicitnone

character*（*）filename_Field

integerpoint,line

open（20,file=filename_Field,status='old'）

read（20,*）

read（20,*）point,line

close（20）

EndSubroutineGet_point_line

******************************************************************c

功能：

扩边数据位置确定子程序

输入参数说明：

mpoint:

实际数据个数

factor_x:

扩边比例因子（>1.0）

输出参数说明：

m1,m2:

实际数据起点位置和终点位置（m2-m1+1=mpoint）

扩边以后总点数（满足2的幂次方）

------------------------------------------------------------------c

SUBROUTINEGET_M1M2m（mpoint,factor_x,m,m1,m2）

implicitnone

INTEGERmpoint,m,m1,m2

REALfactor_x

integermtemp,mu

mtemp=mpoint

dowhile（（mod（mtemp,2）.eq.0）.and.（mtemp.ne.0））

mtemp=mtemp/2

enddo

if（mtemp.eq.1）then

m=mpoint*2

else

mu=int（log（float（mpoint））/0.693147+factor_x）

m=2**mu

endif

m1=1+（m-mpoint）/2

m2=m1+mpoint-1

ENDSUBROUTINEGET_M1M2M

读取低高度观测面数据子程序

SubroutineInput_Field（filename_Field,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,dx,dy,Field）

implicitnone

character*（*）filename_Field

integerpoint,line

integerm,m1,m2,n,n1,n2

realxmin,xmax,ymin,ymax,dx,dy

realX_Field（1:

m）,Y_Field（1:

n）,Field（1:

m,1:

n）

integeri,j

open（30,file=filename_Field,status='old'）

read（30,*）

read（30,*）xmin,xmax

read（30,*）ymin,ymax

read（30,*）

doj=n1,n2,1

read（30,*）（Field（i,j）,i=m1,m2,1）

enddo

close（30）

dx=（xmax-xmin）/（m2-m1）

dy=（ymax-ymin）/（n2-n1）

EndSubroutineInput_Field

二维扩边子程序

SubroutineExpend_cos_2D（m,m1,m2,n,n1,n2,Field）

implicitnone

integerm,m1,m2,n,n1,n2

realField（1:

m,1:

n）

realField_m（1:

m）,Field_n（1:

n）

integeri,j

doj=n1,n2,1

doi=m1,m2,1

Field_m（i）=Field（i,j）

enddo

callExpend_cos_1D（1,m1,m2,m,Field_m）

doi=1,m,1

Field（i,j）=Field_m（i）

enddo

doi=1,m,1

doj=n1,n2,1

Field_n（j）=Field（i,j）

enddo

callExpend_cos_1D（1,n1,n2,n,Field_n）

doj=1,n,1

Field（i,j）=Field_n（j）

enddo

EndSubroutineExpend_cos_2D

一维扩边子程序

SubroutineExpend_cos_1D（n0,n1,n2,n3,Field）

implicitnone

real,parameter:

pi=3.141592654

integern0,n1,n2,n3

RealField（n0:

n3）

integeri,i1,i2

Field（n0）=（Field（n1）+Field（n2））/2.0

Field（n3）=Field（n0）

i1=n0

i2=n1

DOi=i1+1,i2-1,1

Field（i）=Field（i1）+cos（pi/2.0*（i2-i）/（i2-i1））*（Field（i2）-Field（i1））

Enddo

i1=n2

i2=n3

DOi=i1+1,i2-1,1

Field（i）=Field（i1）+cos（pi/2.0*（i2-i）/（i2-i1））*（Field（i2）-Field（i1））

Enddo

EndSubroutineExpend_cos_1D

******************************************************************c

功能：

复数组2-D快速Fourier变换

输入参数说明：

Freal（m,n）:

输入数据实部

Fimage（m,n）:

输入数据虚部

点数（M必须为2的幂次方）

线数（N必须为2的幂次方）

NF:

正、反变换标志量（1:

反变换;2:

正变换）

输出参数说明：

Freal（m,n）:

变换后频谱实部

Fimage（m,n）:

变换后频谱虚部

对应频率分布说明：

数据Freal（m,n）和Fimage（m,n）对应的频率分布位置为:

m方向：

0,1,.......,m/2-1,m/2,-（m/2-1）,......,-1

n方向：

0,1,.......,n/2-1,n/2,-（n/2-1）,......,-1

------------------------------------------------------------------c

SUBROUTINEFFT2（Freal,Fimage,m,n,nf）

implicitnone

INTEGERm,n,nf

REALFreal（1:

m,1:

n）,Fimage（1:

m,1:

n）

real,ALLOCATABLE:

Treal（:

）,Timage（:

）

integernmmax,ierr,i,j

nmmax=max（m,n）

allocate（Treal（1:

nmmax）,Timage（1:

nmmax）,STAT=ierr）

if（ierr.ne.0）STOP'数组分配错误，STOP!

DOi=1,m,1

IF（n.ne.1）THEN

doj=1,n,1

Treal（j）=Freal（i,j）

Timage（j）=Fimage（i,j）

enddo

callFFT（Treal,Timage,n,nf）

doj=1,n,1

Freal（i,j）=Treal（j）

Fimage（i,j）=Timage（j）

enddo

ENDIF

ENDDO

DOj=1,n,1

IF（m.ne.1）THEN

doi=1,m,1

Treal（i）=Freal（i,j）

Timage（i）=Fimage（i,j）

enddo

callFFT（Treal,Timage,m,nf）

doi=1,m,1

Freal（i,j）=Treal（i）

Fimage（i,j）=Timage（i）

enddo

ENDIF

ENDDO

deallocate（Treal,Timage,STAT=ierr）

ENDSUBROUTINEFFT2

******************************************************************c

功能：

复数组1-D快速Fourier变换

输入参数说明：

Freal（n）:

输入数据实部

Fimage（n）:

输入数据虚部

点数（M必须为2的幂次方）

NF:

正、反变换标志量（1:

反变换;2:

正变换）

输出参数说明：

Freal（n）:

变换后频谱实部

Fimage（n）:

变换后频谱虚部

对应频率分布说明：

数据Freal（n）和Fimage（n）对应的频率分布位置为:

0,1,.......,n/2-1,n/2,-（n/2-1）,......,-1

------------------------------------------------------------------c

SUBROUTINEFFT（Xreal,Ximage,n,nf）

implicitnone

INTEGERn,nf

REALXreal（1:

n）,Ximage（1:

n）

integernu,n2,nu1,k,k1,k1n2,l,i,ibitr

realf,p,arg,c,s,treal,timage

nu=int（log（float（n））/0.693147+0.001）

n2=n/2

nu1=nu-1

f=float（（-1）**nf）

k=0

DOl=1,nu,1

DOwhile（k.lt.n）

doi=1,n2,1

p=ibitr（k/2**nu1,nu）

arg=6.2831853*p*f/float（n）

c=cos（arg）

s=sin（arg）

k1=k+1

k1n2=k1+n2

treal=Xreal（k1n2）*c+Ximage（k1n2）*s

timage=Ximage（k1n2）*c-Xreal（k1n2）*s

Xreal（k1n2）=Xreal（k1）-treal

Ximage（k1n2）=Ximage（k1）-timage

Xreal（k1）=Xreal（k1）+treal

Ximage（k1）=Ximage（k1）+timage

k=k+1

enddo

k=k+n2

ENDDO

k=0

nu1=nu1-1

n2=n2/2

ENDDO

DOk=1,n,1

i=ibitr（k-1,nu）+1

if（i.gt.k）then

treal=Xreal（k）

timage=Ximage（k）

Xreal（k）=Xreal（i）

Ximage（k）=Ximage（i）

Xreal（i）=treal

Ximage（i）=timage

endif

ENDDO

IF（nf.ne.1）THEN

doi=1,n,1

Xreal（i）=Xreal（i）/float（n）

Ximage（i）=Ximage（i）/float（n）

enddo

展开阅读全文