实验二三度直立六面体正演程序设计实验报告.docx
- 文档编号:12273256
- 上传时间:2023-06-05
- 格式:DOCX
- 页数:29
- 大小:163.68KB
实验二三度直立六面体正演程序设计实验报告.docx
《实验二三度直立六面体正演程序设计实验报告.docx》由会员分享,可在线阅读,更多相关《实验二三度直立六面体正演程序设计实验报告.docx(29页珍藏版)》请在冰点文库上搜索。
实验二三度直立六面体正演程序设计实验报告
《重磁资料处理与解释》实验二
直立六面体正演程序设计
专业名称:
地球物理学
学生姓名:
学生学号:
指导老师:
王万银、纪新林、纪晓琳、邱世灿
提交日期:
2016-11-30
5结论及建议6
附录:
源程序代码7
图1直立六面体模型示意图
1基本原理
在空间直角坐标系o-xyz中,形体(直立六面体)模型如图1所示。
设该直立六面体x方向的坐标范围为
,y方向坐标为
,z方向(铅垂向下为正)坐标为
;又设该直立六面体剩余密度为
,
根据正演理论得知,其在空间任意一点
处产生的重力异常为
(1-1)
式中,
为万有引力常数,在国际单位制中其值为
;
为计算点
到场源点
的距离,可表示为
(1-2)
2输入/输出数据格式设计
2.1场源参数数据格式设计
场源参数按照一个直立六面体为一个记录进行设计,在数据文件中占一行。
第一列为剩余密度density_source(g/cm3);第二列为磁化强度mag_source(nT);第三列为磁化方向倾角Az_source(DEG);第四列为磁化方向与x轴的夹角Ax_source(DEG);第五列为磁化方向与y轴的夹角Ay_source(DEG);第六列~第七列为x坐标的起点X1和终点X2(km);第八列~第九列为y坐标的起点Y1和终点Y2(km);第十列~第十一列为z坐标的起点Z1和终点Z2)(km,向下为正)。
以上各量均为实型变量,各量的意义见图1所示。
2.2计算点坐标数据格式设计
计算点坐标数据格式设计为1:
曲面非规则网;2:
曲面规则网;3:
平面非规则网;4:
平面规则网;共四种形式。
其中非规则网采用一个计算点为一个记录的方式设计。
第1列保存计算点x坐标x_coord,第2列保存计算点y坐标y_coord,第3列保存计算点z坐标z_coord。
以上各量均为实型变量。
规则网采用grd文件形式进行设计。
2.3计算结果输出数据格式设计
计算结果输出数据格式与输入格式对应。
非规则网时,采用一个计算点为一个记录的方式设计。
第1列保存计算点x坐标x_coord,第2列保存计算点y坐标y_coord,第3列保存计算点计算结果field。
以上各量均为实型变量。
规则网时采用grd文件形式进行输出。
2.4参数文件数据格式设计
将以上部分量保存在一个文件中,该文件名变量为cmdfile,字符串变量,长度不超过80,全路径名。
在该文件中保存的参数如下:
场源参数文件名filename_source,字符串变量,长度不超过80;
计算点坐标文件名filename_obser,字符串变量,长度不超过80;
计算结果输出文件名filename_field,字符串变量,长度不超过80
anglez:
t0方向与Z轴夹角
angley:
t0方向与Y轴夹角
anglex:
t0方向与Z轴夹角
flag:
测网类型(1:
曲面规则网;2:
曲面非规则网;3:
平面规则网;4:
平面非规则网)
z_par:
如果测网类型是平面网时Z坐标的值
3总体设计
此次程序采用IPO结构设计,首先通过读取cmd文件,得到相关输入参数:
输入场源文件名、计算点坐标文件名、输出结果文件名,磁倾角、与X轴夹角和y轴夹角,给定计算异常类型,给定测网类型,绘图所以特征值;然后从场源文件中读取输入的场源个数及场源参数。
下一步,通过判断测网类型type_net的数值,选择测网类型是规则网(type_net=1或3)还是非规则网(type_net=2或4),选择好测网类型后,输入与测网类型相符的计算点坐标。
然后通过判断method的数值,选择计算重力异常(method=1),磁力异常(method=5)还是化极磁力异常(method=6)。
最后,根据不同的测网类型输出不同的计算结果。
总体设计见表1。
输入场源文件名、计算点坐标文件名、输出结果文件名,磁倾角、与X轴夹角和y轴夹角,给定计算异常类型,给定测网类型,绘图所以特征值。
输入场源个数及参数(
)
判断测网类型type_net的数值
type_net=1或3
曲面或平面非规则网
type_net=2或4
曲面或平面规则网
输入计算点坐标(x,y,z)
输入计算点坐标(x,y,z)
判断method的数值
判断method的数值
method=1计算重力异常
method=5计算磁力异常
method=6计算化极磁力异常
method=1计算重力异常
method=5计算磁力异常
method=6计算化极磁力异常
输出计算结果
输出计算结果
表1总体设计N-S图
4测试结果
4.1测试参数
(1)场源参数保存在“source.dat”中。
第一列为剩余密度(g/cm3);第二列为磁化强度(nT);第三列为磁化方向与z轴的夹角(DEG);第四列为磁化方向与x轴的夹角(DEG);第五列为磁化方向与y轴的夹角(DEG);第六列~第七列为x坐标的起点和终点(km);第八列~第九列为y坐标的起点和终点(km);第十列~第十一列为z坐标的起点和终点(km,向下为正)。
0.7,20000,50,85,5,-9,-3,4,8,2,7
0.8,34000,50,85,5,0,7,0,5,3,7
0.9,17000,50,85,5,-4,8,-8,-3,2,7
(2)计算点为曲面规则网,数据保存在“bxyzl.grd”文件中(GRD格式)。
形式如下:
DSAA
2727
-26.026.0
-26.026.0
-5.31.7
0.150-0.15-0.3-0.4-0.5-0.6-0.7-0.85
(3)有关参数保存在forward_3D.cmd文件中,如下:
'场源参数文件名'source.dat
'计算面坐标文件名'BXYZL.GRD
'计算异常结果输出文件名'output_field.grd
'磁化方向与z轴夹角,与x轴夹角,与y轴夹角'50,85,5
type_net:
测网类型(1:
曲面非规则网;2:
曲面规则网;3:
平面非规则网;4:
平面规则网)
如果测网类型是平面网时Z坐标的值'-5.3
4.2测试结果
图1.测区地形图
当测点数据为曲面规则网时,测试结果如下:
1.重力异常:
图3.观测面重力异常(g.u)
2.磁力异常:
图4.观测面磁力异常(nT)
3.化极磁力异常:
图5.观测面化极磁力异常(nT)
4.结果分析:
图3,重力异常在观测面上对异常体的位置定位不明显,等值线只有一个圈闭,显示一个高值异常;
图4,磁力异常在观测面上能够对异常体位置大致判定,高值等值线显示3个圈闭,效果比由重力异常确定的异常体位置精确,但异常上方显示有低值异常;
图5,化极磁力异常在观测面上能准确确定异常体位置,3个异常体大致在3个磁异常高值圈闭区域,同时低值异常消失,分辨效果比图4更好。
故而,重磁异常需联合分析,紧靠重力或磁力异常可能不会很好确定场源形态。
同时对比可看出,磁力异常比重力异常反演精度高,而化极磁力异常效果更好。
5结论及建议
此次三度体直立六面体正演实验较好地得到了场源产生的异常场。
并且我们可看出,重力异常场对场源性质的判定不如磁异常场效果好。
因而在实际反演工作中,可以联合重、磁资料综合解释,单一的重力异常可能不能正确得到反演结果。
附录:
源程序代码
!
*****************************************************************************************************************************
!
功能:
直立六面体正演法计算复杂三度体的重磁异常并输出
!
参数声明:
!
参数文件中的参数:
!
cmdfile:
存放参数的文件(forward_3D.cmd)
!
filename_source:
场源参数文件名
!
filename_coord:
计算点坐标文件名
!
filename_field:
计算点异常结果(重力、磁力或化极磁力异常)输出文件名
!
anglez:
t0方向与Z轴夹角
!
anglex:
t0方向与X轴夹角
!
angley:
t0方向与Y轴夹角
!
anglezz:
t0方向与Z轴夹角,为0
!
type_net:
测网类型(1:
曲面非规则网;2:
曲面规则网;3:
平面非规则网;4:
平面规则网)
!
method:
计算方法(1.计算重力异常;5.计算磁力异常;6.计算化极磁力异常)
!
constant_z:
测网类型为是平面网(type_net=3,4)时Z坐标的值
!
场源点的参数:
!
N_source:
场源点个数
!
density_source:
场源剩余密度的一维数组
!
mag_source:
场源磁化强度的一维数组
!
Az_source:
磁化方向与Z轴夹角的一维数组
!
Ax_source:
磁化方向与X轴夹角的一维数组
!
Ay_source:
磁化方向与Y轴夹角的一维数组
!
X_source:
X坐标的起点和终点的二维数组
!
Y_source:
Y坐标的起点和终点的二维数组
!
Z_source:
Z坐标的起点和终点的二维数组
!
计算点的参数:
!
N_point:
每条测线上的点数(规则网)
!
N_line:
测线数(规则网)
!
N_coord:
计算点数
!
xmin:
计算点x坐标的最小值
!
xmax:
计算点x坐标的最大值
!
ymin:
计算点y坐标的最小值
!
ymax:
计算点y坐标的最大值
!
x_coord:
计算点的x坐标的一维数组
!
y_coord:
计算点的y坐标的一维数组
!
z_coord:
计算点的z坐标的一维数组
!
field:
计算点的重力、磁力或化极磁异常的一维数组
!
****************************************************************************************************************************
programforward_3D
implicitnone
character*80cmdfile
character*80filename_source,filename_coord,filename_field
realanglez,angley,anglex,anglezz
integertype_net,method
realconstant_z,i
integerN_source
real,allocatable:
:
density_source(:
),mag_source(:
),Az_source(:
),Ax_source(:
),Ay_source(:
),X_source(:
:
),Y_source(:
:
),Z_source(:
:
)
integerN_point,N_line,N_coord
realxmin,xmax,ymin,ymax
real,allocatable:
:
x_coord(:
),y_coord(:
),z_coord(:
)
real,allocatable:
:
field(:
)
cmdfile='forward_3D.cmd'
callRead_cmd(cmdfile,filename_source,filename_coord,filename_field,anglez,anglex,angley,anglezz,type_net,method,constant_z)
callread_N_source(filename_source,N_source)
allocate(density_source(1:
N_source),mag_source(1:
N_source),Az_source(1:
N_source),Ax_source(1:
N_source),Ay_source(1:
N_source),X_source(1:
N_source,1:
2),Y_source(1:
N_source,1:
2),Z_source(1:
N_source,1:
2))
callread_source(filename_source,N_source,density_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source)
callread_N_coord(filename_coord,N_point,N_line,N_coord,type_net)
allocate(x_coord(1:
N_coord),y_coord(1:
N_coord),z_coord(1:
N_coord),field(1:
N_coord))
callread_coordinate(filename_coord,N_point,N_line,N_coord,X_coord,Y_coord,z_coord,constant_z,Xmin,Xmax,Ymin,Ymax,type_net)
if(method==1)then
callcalculate_gra(N_source,X_source,Y_source,Z_source,density_source,N_coord,X_coord,Y_coord,Z_coord,field)
elseif(method==5)then
callcalculate_mag(N_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source,N_coord,X_coord,Y_coord,Z_coord,field,anglez,anglex,angley)
elseif(method==6)then
callcalculate_mag(N_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source,N_coord,X_coord,Y_coord,Z_coord,field,anglezz,anglex,angley)
endif
callOutput_field(x_coord,y_coord,field,N_point,N_line,N_coord,xmin,xmax,ymin,ymax,type_net,filename_field)
deallocate(density_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source,x_coord,y_coord,z_coord,field)
endprogram
!
*****************************************************************************************************************************
!
子程序read_cmd
!
功能:
读取参数文件
!
输入参数:
!
cmdfile:
参数文件名
!
输出参数:
!
filename_source:
场源参数文件名
!
filename_coord:
计算点坐标文件名
!
filename_field:
计算点异常结果(重力、磁力或化极磁力异常)输出文件名
!
anglez:
t0方向与Z轴夹角
!
anglex:
t0方向与X轴夹角
!
angley:
t0方向与Y轴夹角
!
anglezz:
t0方向与Z轴夹角,为0
!
type_net:
测网类型(1:
曲面非规则网;2:
曲面规则网;3:
平面非规则网;4:
平面规则网)
!
method:
计算方法(1.计算重力异常;2.计算磁力异常;3.计算化极磁力异常)
!
constant_z:
如果测网类型是平面网时Z坐标的值
!
*****************************************************************************************************************************
subroutineRead_cmd(cmdfile,filename_source,filename_coord,filename_field,anglez,anglex,angley,anglezz,type_net,method,constant_z)
implicitnone
character*(*)cmdfile
character*(*)filename_source,filename_coord,filename_field
realanglez,anglex,angley,anglezz
integertype_net,method
realconstant_z
open(10,file=cmdfile,status='old')
read(10,*)filename_source
read(10,*)filename_coord
read(10,*)filename_field
read(10,*)anglez,anglex,angley,anglezz
read(10,*)type_net
read(10,*)method
read(10,*)constant_z
close(10)
endsubroutineRead_cmd
!
*****************************************************************************************************************************
!
子程序:
read_N_source
!
功能:
读取场源点的数据个数
!
输入参数:
!
filename:
文件名
!
输出参数:
!
N_source:
场源个数
!
*****************************************************************************************************************************
subroutineread_N_source(filename,N_source)
implicitnone
character*(*)filename
integerN_source
open(10,file=filename,status='old')
N_source=0
dowhile(.not.eof(10))
read(10,*,end=100,ERR=100)
N_source=N_source+1
100enddo
close(10)
endsubroutineread_N_source
!
**********************************************************************************************************************************
!
子程序:
read_source
!
功能:
读取场源点参数
!
输入参数:
!
filename_source:
场源参数文件名变量
!
N_source:
场源个数
!
输出参数:
!
density_source:
场源剩余密度的一维数组
!
mag_source:
场源磁化强度的一维数组
!
Az_source:
磁化方向与Z轴夹角的一维数组
!
Ax_source:
磁化方向与X轴夹角的一维数组
!
Ay_source:
磁化方向与Y轴夹角的一维数组
!
X_source:
X坐标的起点和终点的二维数组
!
Y_source:
Y坐标的起点和终点的二维数组
!
Z_source:
Z坐标的起点和终点的二维数组
!
**********************************************************************************************************************************
subroutineread_source(filename_source,N_source,density_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source)
implicitnone
character*(*)filename_source
integerN_source
realdensity_source(1:
N_source),mag_source(1:
N_source),Az_source(1:
N_source),Ax_source(1:
N_source),Ay_source(1:
N_source),X_source(1:
N_source,1:
2),Y_source(1:
N_source,1:
2),Z_source(1:
N_source,1:
2)
integerk
open(10,file=filename_source,status='old')
dok=1,N_source,1
read(10,*)density_source(k),mag_source(k),Az_source(k),Ax_source(k),Ay_source(k),X_source(k,1),X_source(k,2),Y_source(k,1),Y_source(k,2),Z_source(k,1),Z_source(k,2)
enddo
close(10)
endsubroutineread_source
!
*********************************************************************************************************************************
!
子程序:
read_N_coord
!
功能:
读取计算点个数
!
输入参数:
!
filename:
文件名
!
输出参数:
!
m:
每条线上的点数(规则网)
!
n:
测线数(规则网)
!
N_coord:
计算点个数
!
*********************************************************************************************************************************
subroutineread_N_coord(filename,m,n,N_coord,type_
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验 三度 直立 六面体 程序设计 报告