改实验一：二度直立柱体正演程序设计实验报告(202,24李龙坤,20161118)

1、 重磁资料处理与解释实验一二度直立柱体正演程序设计 专业名称：勘查技术与工程学生姓名：李龙坤学生学号：201326020224指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-11-18目 录1 基本原理12 输入/输出数据格式设计12.1 场源参数数据格式设计12.2 计算点坐标数据格式设计22.3 计算结果输出数据格式设计22.4 参数文件数据格式设计23 总体设计24 测试结果34.1 测试参数34.2 测试结果4五 结论及建议4附录：源程序代码51 基本原理12 输入/输出数据格式设计12.1 场源参数数据格式12.2 计算点坐标数据格式设计22.3 计算结果输出数据格式设计

2、22.4 参数文件数据格式设计23 总体设计24 测试结果34.1 测试参数34.2 测试结果45 结论及建议5附录：源程序代码6图1 二度直立柱体示意图1 基本原理在二度直立柱体示意图，但不太确定。直角坐标系o-xz中，形体（二度直立柱体）模型 如图1所示。设该二度直立柱体x方向的坐标范围为， z方向（铅垂向下为正）坐标为；又设该二度直立柱体剩余密度为。根据正演理论得知，其在空间任意一点处产生的重力异常为 （1-1） （1-1）式中，为万有引力常数，在国际单位制中其值为； 。若有n个二度直立柱体，其引起的重力异常符合叠加原理，由此得到 （1-2） 2 输入/输出数据格式设计2.1 场源参数数

3、据格式设计场源参数按照一个二度直立柱体为一个记录进行设计，在数据文件中占一行。第一列为剩余密度density_source(g/cm3)；第二列第三列为x坐标的起点和终点起点X1和终点X2(m)；第四列第五列为z坐标的起点和终点Z1和终点Z2)(m，向下为正)。以上各量均为实型变量，各量的意义见图1所示。例如：0.2 -100 -50 50 2000.3 -50 50 50 2002.2 计算点坐标数据格式设计 计算点坐标数据格式设计为规则网测线，采用一个计算点为一个记录的方式设计。第1列保存计算点x坐标x_coord（m），第z列保存计算点z坐标z_coord（m）。以上各量均为实型变量。保

4、存在Input_file_source文件名变量中。例如：用一个模型数据在此示意一下会更好。以下数据格式设计相同。已经知道了 -400 -20-398 -202.3 计算结果输出数据格式设计计算结果输出数据格式与输入格式对应，设计为非规则网规则测线，采用一个计算点为一个记录的方式设计。第1列保存计算点x坐标x_coord（m），第2列保存计算点z坐标z_coord（m），第3列保存计算点计算结果field（mGal）。以上各量均为实型变量。保存在Output_file_source文件名变量中。例如：用一个模型数据在此示意一下会更好。以下数据格式设计相同。-400 -20 0.24-398 -

5、20 0.352.4 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。在该文件中保存的参数如下：场源参数个数： N_SOURCE，整型数据。计算点参数个数：N_COORD，整型数据。场源参数文件名input_file_source，字符串变量，长度不超过80；计算点坐标文件名input_file_coordinate，字符串变量，长度不超过80；计算结果输出文件名output_file_field，字符串变量，长度不超过803 总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入场源参数个数、计

6、算点参数个数、场源文件名、计算点坐标文件名、输出结果文件名；再通过场源参数文件，得到场源密度及其X、Y坐标；然后用过计算点参数文件，得到计算点X、Y坐标。然后利用重力异常计算公式得到每个计算点所对应的重力异常，然后由于流程图可能会出现乱码现象，因此这次更改为N-S图。将其输出到指定文件，利用GRAFER画图。输入参数文件名、输入场源文件名、计算点坐标文件名、输出结果文件名输入场源个数以及参数（ N_SOURCE,N_COORD:场源点与计算点个数、density(1:N_SOURCE):场源密度、x_source(1:N_SOURCE,1:2):场源x点坐标、 z_source(1:N_SOU

7、RCE,1:2):场源z点坐标、万有引力常数G）输入计算点坐标（ x_coord(1:N_COORD):计算点x坐标、z_coord(1:N_COORD):计算点y坐标）利用重力异常公式计算重力异常将重力异常输出到指定文件中图2 总体设计N-S图K=1,开始N_SOURCE,N_COORD,density,x_source(:,:),z_source(:,:),x_coord(:),z_coord(;)Field(k)=0.0N=1 i=1 j=1Field(k)=field(k)+(-1)*(i+j)*delt()j>2 是 否j=j+1i>2 是 否i=i+1N>N_SO

8、URCE 是 k>N_COORDN=N+1 否 是 否N_COORD=N_COORD+1输出field（k）结束 流程图1 2 计算二度直立柱体重力异常正演流程图 4 测试结果4.1 测试参数 （31）有关参数保存在cmdfile.txt文件中，如下:3401source.datcoord.txtfield.txt （12）场源参数保存在“source.dat”中。第一列为剩余密度(g/cm3)；第二列第三列为x坐标的起点和终点(m)；第四列第五列为z坐标的起点和终点(m，向下为正)。 0.2 -100 -50 50 2000.3 -50 50 50 2000.2 50 100 50 2

9、00 （23）计算点为规则网（图2）图2 是什么？根据王老师给的模板写的，但是没注意到这个小细节，已经删掉。，数据保存在“coord.txt”文件中（txt格式）。形式如下： -400 -20-398 -20-397 -20 398 20400 20 （3）有关参数保存在cmdfile.txt文件中，如下: 3 401source.datcoord.txtfield.txt4.2 测试结果 利用上面的测试参数，通过重力异常程序计算并利用GRAFER画图，如下：图3 重力异常剖面图通过上图3上图没有图名已添加可知，地下密度异常体的在关于0轴呈对称分布，最大异常在0.6mgal,最小异常在0.08

10、mgal。五 结论及建议 经测试，此次程序设计是正确的。由测试结果所得图可知，处理所得的重力异常的效果均较好。但是程序还是比较冗长，有待改进。进。附录：源程序代码! 功能：计算二度体在空间任意一点的重力异常值! 输入参数说明：cmdfile:输入参数文件名! input_file_source:场源参数文件名! input_file_coordinate：计算点参数文件名! output_file_field：输出参数文件名! N_SOURCE,N_COORD:场源点与计算点个数! density(1:N_SOURCE):场源密度! x_source(1:N_SOURCE,1:2):场源x点坐

11、标! z_source(1:N_SOURCE,1:2):场源z点坐标! x_coord(1:N_COORD):计算点x坐标! z_coord(1:N_COORD):计算点y坐标! field(1:N_COORD):异常值 !PROGRAM forward_2DIMPLICIT NONE character*80 cmdfile character*80 input_file_source character*80 input_file_coordinate character*80 output_file_field real,allocatable:x_source(:,:),z_sourc

12、e(:,:),density(:) real,allocatable:x_coord(:),z_coord(:),field(:) integer N_SOURCE,N_COORD cmdfile="cmdfile.txt" call read_cmd(cmdfile,N_SOURCE,input_file_source,N_COORD,input_file_coordinate,output_file_field) ALLOCATE(x_source(1:N_SOURCE,1:2),z_source(1:N_SOURCE,1:2),density(1:N_SOURCE)

13、ALLOCATE(x_coord(1:N_COORD),z_coord(1:N_COORD),field(1:N_COORD) CALL input_source_2D_vertical(input_file_source,x_source,z_source,N_SOURCE,density) call input_coord_2D_vertical(input_file_coordinate,N_COORD,x_coord,z_coord) call forward_2D_vertical(x_source,z_source,N_SOURCE,density,N_COORD,x_coord,

14、z_coord,field) call output_field_2D_vertical(output_file_field,N_COORD,x_coord,z_coord,field) DEALLOCATE(x_source,z_source,density) DEALLOCATE(x_coord,z_coord,field) end program forward_2D! 功能：利用参数文件读取各种参数! 输入参数说明：cmdfile:参数文件名! 输出参数说明：N_SOURCE,N_COORD:场源点与计算点个数! input_file_source:场源参数文件名! input_fil

15、e_coordinate：计算点参数文件名! output_file_field：输出参数文件名!SUBROUTINE read_cmd(cmdfile,N_SOURCE,input_file_source,N_COORD,input_file_coordinate,output_file_field) IMPLICIT NONE CHARACTER*(*) cmdfile integer N_SOURCE,N_COORD character*(*) input_file_source,input_file_coordinate,output_file_field open(10,file=c

16、mdfile,status='old') read(10,*) N_SOURCE read(10,*) N_COORD read(10,*) input_file_source read(10,*) input_file_coordinate read(10,*) output_file_field close(10)end SUBROUTINE read_cmd! 功能:场源点数据读入! 输入参数说明：input_file_source:场源参数文件名! N_SOURCE:场源点个数! 输出参数说明: density(1:N_SOURCE):场源密度! x_source(1:

17、N_SOURCE,1:2):场源x点坐标! z_source(1:N_SOURCE,1:2):场源z点坐标!subroutine input_source_2D_vertical(input_file_source,x_source,z_source,N_SOURCE,density) implicit none character*(*) input_file_source integer N_SOURCE,k real x_source(1:N_SOURCE,1:2),z_source(1:N_SOURCE,1:2),density(1:N_SOURCE) open(20,file=inp

18、ut_file_source,status='old') do k=1,N_SOURCE,1 read(20,*) density(k),x_source(k,1),x_source(k,2),z_source(k,1),z_source(k,2) WRITE(*,*) density(k),x_source(k,1),x_source(k,2),z_source(k,1),z_source(k,2) enddo close(10)end subroutine input_source_2D_vertical! 功能:计算点数据输入! 输入参数说明: N_COORD:计算点个数

19、! input_file_coordinate：计算点参数文件名! 输出参数说明: x_coord(1:N_COORD):计算点x坐标! z_coord(1:N_COORD):计算点y坐标!subroutine input_coord_2D_vertical(input_file_coordinate,N_COORD,x_coord,z_coord) implicit none character*(*) input_file_coordinate integer N_COORD,k real x_coord(1:N_COORD),z_coord(1:N_COORD) open(20,file

20、=input_file_coordinate,status='old') do k=1,N_COORD,1 read(20,*) x_coord(k),z_coord(k) end do close(20)end subroutine input_coord_2D_vertical! 功能:计算重力异常! 输入参数说明: N_SOURCE:场源点个数! density(1:N_SOURCE):场源密度! x_source(1:N_SOURCE,1:2):场源x点坐标! z_source(1:N_SOURCE,1:2):场源z点坐标! N_COORD:计算点个数! x_coord

21、(1:N_COORD):计算点x坐标! z_coord(1:N_COORD):计算点y坐标! 输出参数说明: field(1:N_COORD):异常值!subroutine forward_2D_vertical(x_source,z_source,N_SOURCE,density,N_COORD,x_coord,z_coord,field) implicit none integer N_SOURCE real x_source(1:N_SOURCE,1:2),z_source(1:N_SOURCE,1:2),density(1:N_SOURCE) integer N_COORD real

22、x_coord(1:N_COORD),z_coord(1:N_COORD) integer k,n,j,i real field(1:N_COORD) real Deltg do k=1,N_COORD,1 field(k)=0.0 do n=1,N_SOURCE,1 do i=1,2,1 do j=1,2,1 field(k)=field(K)+(-1)*(i+j)*Deltg(density(n),x_source(n,i),z_source(n,j),x_coord(k),z_coord(k) end doend doend do end doend subroutine forward_2D_vertical! 功能：外部子函数异常计算公式! 输入参数说明：DE:密度! XI:场源点x坐标! ZJ:场源点z坐标! Z:计算点x坐标! G:万有引力常数! 输出参数说明: Deltg:计算值!real function Deltg(DE,XI,ZJ,X,Z) REAL DE,XI,ZJ,X,Z REAL:G=6.672/1000这个值和你