频率域位场处理和转换实验.docx
- 文档编号:17041371
- 上传时间:2023-07-21
- 格式:DOCX
- 页数:31
- 大小:329.04KB
频率域位场处理和转换实验.docx
《频率域位场处理和转换实验.docx》由会员分享,可在线阅读,更多相关《频率域位场处理和转换实验.docx(31页珍藏版)》请在冰点文库上搜索。
频率域位场处理和转换实验
《重磁资料处理与解释》实验二
频率域位场处理和转换实验
学院:
地测学院
专业名称:
勘查技术与工程
学生某某:
学生学号:
指导教师:
提交日期:
2018年1月9日
二0一八年一月
1根本原理
1.1位场的方程
由场论知识可知,位场方程分为
两大类:
有源的Possion方程
,以与无源的Laplace方程
。
Laplace方程的第一边值问题
通常为Dirichlet问题,第二边值问题
通常称为Nueman问题。
假如P点在S平面内称为内部问题,反之称为外部问题。
由唯一性定理可知,Dirichlet的内部和外部问题的解是唯一的,而Nueman内部问题的解不是唯一的,有一常数差,但其外部问题解是唯一的。
外部问题的解的唯一性的原因:
无源区域位场可以表示为:
〔1-1〕
〔1-2〕
1.2二维傅里叶变换与卷积性质
〔1〕傅里叶变换
〔1-3〕
〔1-4〕
〔2〕卷积性质
〔1-5〕
〔1-6〕
1.3频率域位场延拓原理
当实测平面的异常时,换算场源以外的异常称之为延拓,分为向上延拓和向下延拓。
半空间狄利克莱问题解析解:
〔1-7〕
其中:
称为延拓因子,
为计算面Z坐标,Z轴向下为正方向,
为计算面频率域位场,
为延拓面的频率域位场。
2输入/输出数据格式设计
2.1输入数据格式设计
观测面位场数据保存在filename_obser文件中,为.grd格式。
计算延拓误差时的准确场值文件保存在filename_obser2中,为.grd格式。
2.2输出数据格式设计
实际计算得到的数据保存在filename_output1和文件filename_output2中,为.grd格式。
2.3参数文件数据格式设计
将所要读取的参数保存在一个文件中,该文件名变量为cmdfile,字符串变量,长度不超过80,全路径名。
在该文件中保存的参数如下:
filename_obser:
低高度观测面位场数据文件
filename_output1:
向上延拓后位场数据文件
filename_output2:
向下延拓后位场数据文件
filename_obser2:
高高度观测面数据文件
factor_m:
扩边因子
distance1:
向上延拓的高度〔m/z轴向下为正方向〕
distance2:
向下延拓的高度〔m/z轴向下为正方向〕
3总体设计
3.1频率域位场处理与转换的一般步骤
3.2软件总体设计结果流程图
此次程序采用IPO结构设计,首先通过读取cmd文件,得到相关输入参数:
观测面位场数据文件名变量、延拓后位场数据文件名变量、延拓后准确位场数据文件名变量、扩边因子、延拓的高度〔m/z轴向下为正方向〕;然后确定确定扩边网格的大小,扩边数据点号位置;再从观测面位场数据文件中读取数据。
下一步,进展二维余弦扩边,将扩完边的数据进展快速二维傅里叶变换,转换到频率域;接下来计算延拓因子并且将扩完边的数据进展快速二维傅里叶变换后在频率域与延拓因子相乘;最后进展快速二维傅里叶反变换并且去除扩边局部后输出。
总体设计见表1。
表1总体设计N-S图
I:
〔1〕输入有关参数
〔2〕计算有关参数
〔3〕从文件输入有关数据
P:
(1)进展扩边处理
(2)Fourier正变换
(3)计算延拓因子
(4)进展乘积运算
(5)Fourier反变换
O:
输出计算结果
4测试结果
4.1测试参数
〔1〕向上延拓
原始场值数据保存在’’a20_mag.grd〞中,向上延拓3.3m,延拓后理论数据保存在“a53_mag.grd〞中,延拓后的数据保存在output1.grd中。
网格大小为27*27〔m〕,Xmin=-26m,Xmax=26m,Ymin=-26m,Ymax=26m.
〔2〕向下延拓
原始场值数据保存在’’a53_mag.grd〞中,向下延拓3.3m,延拓后理论数据保存在“a20_mag.grd〞中延拓后的数据保存在output2.grd中。
网格大小为27*27〔m〕,Xmin=-26m,Xmax=26m,Ymin=-26m,Ymax=26m.
有关参数保存在filename.cmd文件中,如下:
A20_mag.grd
A53_mag.grd
output1.grd
output2.grd
1
3.3
-3.3
4.2测试结果
5结论与建议
由测试结果可知,向上延拓可以对浅部异常进展压制,且误差较小;向下延拓可突出浅部异常,且延拓结果误差较大。
通过本次频率域位场处理与转换设计实验,我收获颇丰,同时也感触很深!
首先,刚开始我对空间域转频率域以与位场延拓的概念与处理是有不理解的,但是经过教师耐心细致的讲解和同学们之间的帮助,最终克制了困难。
其次,加深了我对于位场延拓的作用和具体延拓步骤的理解。
感谢教师孜孜不倦的教诲与同学不厌其烦的讲解。
谢谢!
附录:
源程序代码
!
******************************************************************************
!
程序功能:
实现频率域位场延拓
!
cmd文件参数:
!
cmdfile:
存放有关参数的文件名变量
!
filename_obser:
低高度观测面位场数据文件
!
filename_output1:
向上延拓后位场数据文件
!
filename_output2:
向下延拓后位场数据文件
!
filename_obser2:
高高度观测面数据文件
!
factor_m:
扩边因子
!
distance1:
向上延拓的高度〔m/z轴向下为正方向〕
!
distance2:
向下延拓的高度〔m/z轴向下为正方向〕
!
.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方向扩边后数据起点和终点点号位置
!
延拓参数:
!
Ur:
初始观测面信号的实部
!
Ui:
初始观测面信号的虚部
!
Fr:
延拓观测面的信号的实部
!
Fi:
延拓观测面的信号的虚部
!
U:
延拓后准确场值
!
W(m,n):
径向圆频率
!
Y(m,n):
延拓因子
!
误差参数:
!
error:
延拓后场值与场值的均方误差
!
******************************************************************************
Programplyyt
implicitnone
realeigval
parameter(eigval=3.701411*1e5)
character*(80)cmdfile,filename_obser,filename_output1,filename_output2,filename_obser2
real,allocatable:
:
Ur(:
:
),Ui(:
:
),y(:
:
),Fr(:
:
),Fi(:
:
),U(:
:
),W(:
:
)
integerN_point,N_line
integerm,n,m1,m2,n1,n2
integerfactor_m
realxmin,xmax,ymin,ymax,dx,dy
realdistance1,distance2,error
cmdfile='filename.cmd'
callread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,&
filename_obser2,distance1,distance2)
callread_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)
callcalculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)
allocate(Ur(1:
m,1:
n),Ui(1:
m,1:
n),y(1:
m,1:
n),Fr(1:
m,1:
n),Fi(1:
m,1:
n),U(1:
m,1:
n),W(1:
m,1:
n))
callinput_grd(Ur,filename_obser,m1,m2,n1,n2,m,n)
callinput_grd(U,filename_obser2,m1,m2,n1,n2,m,n)
!
低高度向上延拓
callexpand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)
callFFT2(Ur,Ui,m,n,2)
CALLcal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)
callWAVE2D(m,n,dx,dy,W)
callcalculate_Y(m,n,W,y,distance1)
callcalculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)
callFFT2(Fr,Fi,m,n,1)
calloutput_grd(filename_output1,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax,Ymin,Ymax)
callcalculate_error(distance1,Fr,U,m1,m2,n1,n2,m,n,error)
!
高高度向下延拓
callexpand_cos_2D(m1,m2,m,n1,n2,n,U,Ui)
callFFT2(U,Ui,m,n,2)
CALLcal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)
callWAVE2D(m,n,dx,dy,W)
callcalculate_Y(m,n,W,y,distance2)
callcalculate_FmulY(m,n,U,Ui,y,Fr,Fi)
callFFT2(Fr,Fi,m,n,1)
calloutput_grd(filename_output2,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax,Ymin,Ymax)
callcalculate_error(distance2,Fr,Ur,m1,m2,n1,n2,m,n,error)
deallocate(Ur,Ui,y,Fr,Fi,u,w)
endprogram
!
****************************************************************************
!
子程序:
read_cmd
!
功能:
读取参数文件
!
输入参数说明:
!
cmdfile:
参数文件名
!
输出参数说明:
!
filename_obser:
存放低高度观测面的数据文件
!
filename_output1:
存放向上延拓观测面的数据文件
!
filename_output2:
存放向下延拓观测面的数据文件
!
filename_obser2:
存放高高度观测面的数据文件
!
distance1:
向上延拓的高度〔m/z轴向下为正方向〕
!
distance2:
向下延拓的高度〔m/z轴向下为正方向〕
!
factor_m:
扩边因子
!
****************************************************************************
subroutineread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,&
filename_obser2,distance1,distance2)
implicitnone
character*(*)cmdfile,filename_obser,filename_output1,filename_output2,filename_obser2
integerfactor_m
realdistance1,distance2
character*80str
open(10,file=cmdfile,status='old')
read(10,*)filename_obser
read(10,*)filename_obser2
read(10,*)filename_output1
read(10,*)filename_output2
read(10,*)factor_m
read(10,*)distance1
read(10,*)distance2
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(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)
implicitnone
character*(*)filename_obser
integerN_point,N_line
realXmin,Xmax,Ymin,Ymax
open(10,file=filename_obser,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
integerfactor_m
mtemp=a
mtemp=mtemp/2
Enddo
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
pause
mtemp=b
mtemp=mtemp/2
Enddo
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
pause
endsubroutinecalculate_mn
!
*************************************************************************
!
子程序:
INPUT_GRD
!
功能:
读取grd文件中的数据
!
输入参数说明:
!
filename_obser:
输入文件名
!
m1,m2:
x方向实际数据起点位置和终点位置点号
!
m:
扩边后数据终点位置点号(起点位置点号为1)
!
n1,n2:
y方向实际数据起点位置和终点位置点号
!
n:
扩边后数据终点位置点号(起点位置点号为1)
!
输出参数说明:
!
A:
存放输出数据的二维数组名
!
*************************************************************************
SUBROUTINEINPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)
character*(*)filename_obser
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=filename_obser,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
DOi=i1+1,i2,1
Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1))
Enddo
Endsubroutineexpand_cos_1d
!
************************************************************************
!
功能:
FFT2
!
功能:
复数组2-D快速Fourier变换
!
输入参数说明:
!
m0,m3:
x方向的起点和终点
!
n0,n3:
y方向的起点和终点
!
field:
输入信号(需要赋值给Freal,实部)
!
m,n:
x,y方向扩边后数据终点点号位置〔起始点号为1〕
!
NF:
正、反变换标志量(1:
反变换;2:
正变换)
!
输出参数说明:
!
Freal:
信号的实部
!
Fimage:
信号的虚部(对于实信号而言,赋值为0)
!
对应频率分布说明:
!
数据Freal(m,n)和Fimage(m,n)对应的频率分布位置为:
!
m方向:
0,1,.......,m/2-1,m/2,-(m/2-1),.....
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 频率 域位场 处理 转换 实验