《重磁资料处理与解释》实验一二度直立柱体正演程序设计专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-11-29目录1 基本原理 (1)2 输入/输出数据格式设计 (1)2.1 场源参数数据格式设计 (1)2.2 计算点坐标数据格式设计 (1)2.3 计算结果输出数据格式设计 (2)2.4 参数文件数据格式设计 (2)3总体设计 (2)4测试结果 (3)4.1 测试参数 (3)4.2 测试结果 (4)5结论及建议 (4)附录：源程序代码 (5)1 基本原理在空间直角坐标系o-xyz 中，形体（二度体）模型如图1所示。

设该直立六面体x 方向的坐标范围为21~ξξ，z 方向（铅垂向下为正）坐标为21~ζζ；又设该直立六面体剩余密度为σ，根据正演理论得知，其在空间任意一点),,(z y x 处产生的重力异常为()()()()2222ln 2arctan 11z x g G V G x x z z z ξζξσσξξζζξζζ⎧⎫⎛⎫-⎡⎤∆==--+-+-⎨⎬ ⎪⎣⎦-⎝⎭⎩⎭ （1-1） 式中，G 为万有引力常数，在国际单位制中其值为()2311-m 10676s kg ⋅⨯/.。

2 输入/输出数据格式设计2.1 场源参数数据格式设计场源参数按照一个二度体为一个记录进行设计，在数据文件中占一行。

第一列为剩余密度density(g/cm 3)；第二列～第三列为x 坐标的起点1ξ和终点2ξ(m)；第四列～第五列为z 坐标的起点ζ1和终点2ζ(m ，向下为正)。

以上各量均为实型变量，各量的意义见图1所示。

2.2 计算点坐标数据格式设计计算点坐标数据格式设计为非规则网，采用一个计算点为一个记录的方式设计。

第1列保存计算点x 坐标x_coordinate(m)，第2列保存计算点z 坐标z_coordinate(m)。

以上各量均为实型变量。

图1 直立二度体模型示意图2.3 计算结果输出数据格式设计计算结果输出数据格式与输入格式对应，设计为非规则网，采用一个计算点为一个记录的方式设计。

第1列保存计算点x坐标x_coordinate，第2列保存计算点z坐标z_coordinate，第3列保存计算点计算结果field（mgal）。

以上各量均为实型变量。

2.4 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为cmd_file，字符串变量，长度不超过80，全路径名。

在该文件中保存的参数如下：场源参数文件名Input_file_source，字符串变量，长度不超过80；计算点坐标文件名Input_file_coordinate，字符串变量，长度不超过80；计算结果输出文件名output_file_field，字符串变量，长度不超过803总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入场源文件名、计算点坐标文件名、输出结果文件名；从场源文件中读取输入的场源个数及场源参数。

下一步，输入计算点坐标。

然后计算重力异常，最后，输出计算结果。

总体设计见表1。

图2 总体设计N-S图4测试结果4.1 测试参数（1）场源参数保存在“2D_source.dat”中。

第一列为剩余密度(g/cm3)；第二列～第三列为x坐标的起点和终点(m)；第四列～第五列为z坐标的起点和终点(m，向下为正)。

模型如图3所示，模型参数如下：0.2,-100,-50,50,2000.3,-50,50,50,2000.2,50,100,50,200（2）计算点为平面非规则网（图2），数据保存在“2D_coordinate.dat”文件中,形式如下：-400 -20-398 -20-396 -20... ...396 -20398 -20400 -20Zmin=Zmax=-200m;Xmin=-400m,Xmax=400m,N_coordinate=401（3）有关参数保存在forward_2D.cmd文件中，如下:'场源个数',3'场源参数文件名','2D_source.dat''计算点个数',401'计算面坐标文件名','2D_coordinate.dat''计算结果输出文件名','2D_output.dat'（3）'绘图所用特征值（绘图时，大于eigval的数据才能画出）' 1.70141E+38 4.2 测试结果根据上述参数通过处理计算得到直立二度体观测剖面上重力异常，如图4：由图4得，在x=0处重力异常值为最大，约为0.65mgal，在x=0左右两侧对称分布且逐渐减小。

由于三个异常体关于x=0对称分布，且中间异常体剩余密度高于左右两侧，因而重力异常呈现中间为最高值向两侧逐渐减小。

所得图像与实际异常场的分布一致。

由图4，可看出重力异常曲线变化斜率变化最大处大约在x=-100m以及x=100m处，可大致推出异常体边界所在位置，与实际情况相符。

5结论及建议经测试，此次程序设计是正确的。

由测试结果所得图可知，由多边形截面法对二度直立柱体正演的效果是较好的。

附录：源程序代码!***************************************************************************************** ************************!程序功能：直立柱体法计算二度体重力异常与输出!参数说明：! cmd_file:存放文件名变量、场源个数和计算点个数的文件名变量。

! input_file_source：存放场源参数的文件名变量。

! input_file_coordinate：存放计算点参数的文件名变量。

! output_file_field：存放输出结果的文件名变量。

! N_source:场源个数。

! N_coordinate：计算点个数。

program forward_2dcharacter*80 cmd_filecharacter*80 input_file_source,input_file_coordinate,output_file_fieldinteger N_source,N_coordinatereal,allocatable::density(:),X_source(:,:),Z_source(:,:)real,allocatable::X_coordinate(:),Z_coordinate(:),field(:)cmd_file='forward_2D.cmd'callread_cmd(cmd_file,N_source,input_file_source,N_coordinate,input_file_coordinate,output_fil e_field)allocate(density(1:N_source),X_source(1:N_source,1:2),Z_source(1:N_source,1:2))allocate(X_coordinate(1:N_coordinate),Z_coordinate(1:N_coordinate),field(1:N_coordinate)) call input_source_2d_vertical(input_file_source,N_source,density,X_source,Z_source)callinput_coordinate_2d_vertical(input_file_coordinate,N_coordinate,X_coordinate,Z_coordinate) callforward_2d_vertical(N_source,density,X_source,Z_source,N_coordinate,X_coordinate,Z_coordin ate,field)calloutput_file_2d_vertical(output_file_field,N_coordinate,X_coordinate,Z_coordinate,field) deallocate(density,X_source,Z_source,X_coordinate,Z_coordinate,field)end program!***************************************************************************************** ************************! 子程序read_cmd!功能：读取文件名变量、场源个数和计算点个数。

!形参说明：! cmd_file:存放文件名变量、场源个数和计算点个数的文件名变量。

! N_source:场源个数。

! input_file_source：存放场源参数的文件名变量。

! N_coordinate：计算点个数。

! input_file_coordinate：存放计算点参数的文件名变量。

! output_file_field：存放输出结果的文件名变量。

! X_source：存放场源X方向边界坐标的二维数组。

! Z_source：存放场源Z方向边界坐标的二维数组。

! X_coordinate：存放计算点X坐标的一维数组。

! Z_coordinate：存放计算点Z坐标的一维数组。

! field：存放输出结果的三维数组（X,Z,deltg）。

!***************************************************************************************** ************************subroutineread_cmd(cmd_file,N_source,input_file_source,N_coordinate,input_file_coordinate,output_fil e_field)character*(*)cmd_fileinteger N_source,N_coordinatecharacter*(*)input_file_source,input_file_coordinate,output_file_fieldreal stropen(10,file=cmd_file,status='old')read(10,*)str,N_sourceread(10,*)str,input_file_sourceread(10,*)str,N_coordinateread(10,*)str,input_file_coordinateread(10,*)str,output_file_fieldclose(10)end subroutine read_cmd!***************************************************************************************** ************************! 子程序input_source_2d_vertical!功能：读取各场源（二度体）形状参数（X、Z坐标）及剩余密度。