实验二：三度直立六面体正演程序设计实验报告

1、 重磁资料处理与解释实验二直立六面体正演程序设计 专业名称： 地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-11-30目 录1 基本原理12 输入/输出数据格式设计12.1 场源参数数据格式12.2 计算点坐标数据格式设计22.3 计算结果输出数据格式设计22.4 参数文件数据格式设计23 总体设计24 测试结果34.1 测试参数34.2 测试结果45 结论及建议6附录：源程序代码7图1 直立六面体模型示意图1 基本原理在空间直角坐标系o-xyz中，形体（直立六面体）模型如图1所示。设该直立六面体x方向的坐标范围为，y方向坐标为，z方向（铅垂向下为正

2、）坐标为；又设该直立六面体剩余密度为，根据正演理论得知，其在空间任意一点处产生的重力异常为 （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)；第六列第七列

3、为x坐标的起点X1和终点X2(km)；第八列第九列为y坐标的起点Y 1和终点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 计算结果输出数据格式设计计算结果输出数

4、据格式与输入格式对应。非规则网时，采用一个计算点为一个记录的方式设计。第1列保存计算点x坐标x_coord，第2列保存计算点y坐标y_coord，第3列保存计算点计算结果field。以上各量均为实型变量。规则网时采用grd文件形式进行输出。2.4 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。在该文件中保存的参数如下：场源参数文件名filename_source，字符串变量，长度不超过80；计算点坐标文件名filename_obser，字符串变量，长度不超过80；计算结果输出文件名filename_field，字符串变量，

5、长度不超过80anglez：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）还是非

6、规则网（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

7、,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坐标的起点和终

8、点(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格式）。形式如下： DSAA27 27-26.0 26.0-26.0 26.0-5.3 1.7 0.15 0 -0.15 -0.3 -0.4 -0.5 -0.6 -0.7 -0.85 （3）有关参数保存在forward_3D.cmd文件中，如下:'场源参数文件名

9、' source.dat'计算面坐标文件名' BXYZL.GRD'计算异常结果输出文件名' output_field.grd'磁化方向与z轴夹角,与x轴夹角,与y轴夹角' 50,85,5type_net：测网类型（1：曲面非规则网；2：曲面规则网；3：平面非规则网；4：平面规则网）如果测网类型是平面网时Z坐标的值'-5.34.2 测试结果 图1. 测区地形图当测点数据为曲面规则网时，测试结果如下：1.重力异常：图3. 观测面重力异常（g.u）2.磁力异常：图4. 观测面磁力异常（nT）3.化极磁力异常：图5. 观测面化极磁力异常（

10、nT）4.结果分析：图3，重力异常在观测面上对异常体的位置定位不明显，等值线只有一个圈闭，显示一个高值异常；图4，磁力异常在观测面上能够对异常体位置大致判定，高值等值线显示3个圈闭，效果比由重力异常确定的异常体位置精确，但异常上方显示有低值异常；图5，化极磁力异常在观测面上能准确确定异常体位置，3个异常体大致在3个磁异常高值圈闭区域，同时低值异常消失，分辨效果比图4更好。故而，重磁异常需联合分析，紧靠重力或磁力异常可能不会很好确定场源形态。同时对比可看出，磁力异常比重力异常反演精度高，而化极磁力异常效果更好。5 结论及建议此次三度体直立六面体正演实验较好地得到了场源产生的异常场。并且我们可看出

11、，重力异常场对场源性质的判定不如磁异常场效果好。因而在实际反演工作中，可以联合重、磁资料综合解释，单一的重力异常可能不能正确得到反演结果。附录：源程序代码!*!功能：直立六面体正演法计算复杂三度体的重磁异常并输出!参数声明：!参数文件中的参数：! cmdfile：存放参数的文件(forward_3D.cmd)!filename_source：场源参数文件名!filename_coord：计算点坐标文件名!filename_field：计算点异常结果（重力、磁力或化极磁力异常）输出文件名!anglez：t0方向与Z轴夹角!anglex：t0方向与X轴夹角!angley：t0方向与Y轴夹角! an

12、glezz: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轴夹角的一维数组!

13、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：计算点的重力、磁力或化极磁异常的一维数组!*program forwar

14、d_3Dimplicit nonecharacter*80 cmdfilecharacter*80 filename_source,filename_coord,filename_fieldreal anglez,angley,anglex,anglezzinteger type_net,methodreal constant_z,iinteger N_sourcereal,allocatable : density_source(:),mag_source(:),Az_source(:),Ax_source(:),Ay_source(:),X_source(:,:),Y_source(:,:

15、),Z_source(:,:)integer N_point,N_line,N_coordreal xmin,xmax,ymin,ymaxreal,allocatable : x_coord(:),y_coord(:),z_coord(:)real,allocatable : field(:)cmdfile='forward_3D.cmd'call Read_cmd(cmdfile,filename_source,filename_coord,filename_field,anglez,anglex,angley,anglezz,type_net,method,constant

16、_z)call read_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)call read_source(filename_source,N_source,density_sour

17、ce,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source)call read_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)call read_coordinate(filename_coord,N_point,N_line,N_coord,X_coord,Y_coord,z_coord

18、,constant_z,Xmin,Xmax,Ymin,Ymax,type_net)if(method=1) thencall calculate_gra(N_source,X_source,Y_source,Z_source,density_source,N_coord,X_coord,Y_coord,Z_coord,field)else if(method=5) thencall calculate_mag(N_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source,N_coord,X_coord,

19、Y_coord,Z_coord,field,anglez,anglex,angley)else if(method=6) thencall calculate_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)endifcall Output_field(x_coord,y_coord,field,N_point,N_line,N_coord,xmin,xmax,y

20、min,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)end program!*! 子程序read_cmd!功能：读取参数文件!输入参数：!cmdfile：参数文件名!输出参数：!filename_source：场源参数文件名!filename_coord：计算点坐标文件名!filename_field：计算点异常结果（重力、磁力或化极磁

21、力异常）输出文件名!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坐标的值!*subroutine Read_cmd(cmdfile,filename_source,filename_coord,filename_field,anglez,anglex,angley,a

22、nglezz,type_net,method,constant_z)implicit nonecharacter *(*) cmdfilecharacter *(*) filename_source,filename_coord,filename_fieldreal anglez,anglex,angley,anglezzinteger type_net,methodreal constant_zopen(10,file=cmdfile,status='old')read(10,*) filename_sourceread(10,*) filename_coordread(10

23、,*) filename_fieldread(10,*) anglez,anglex,angley,anglezzread(10,*) type_netread(10,*) methodread(10,*) constant_zclose(10)end subroutine Read_cmd!*! 子程序：read_N_source!功能：读取场源点的数据个数!输入参数：!filename:文件名!输出参数：!N_source:场源个数!*subroutine read_N_source(filename,N_source)implicit nonecharacter *(*) filenam

24、einteger N_sourceopen(10,file=filename,status='old')N_source=0do while(.not.eof(10)read(10,*,end=100,ERR=100)N_source=N_source+1100 end doclose(10)end subroutine read_N_source!*! 子程序：read_source!功能：读取场源点参数!输入参数：!filename_source:场源参数文件名变量!N_source：场源个数!输出参数：!density_source：场源剩余密度的一维数组!mag_sou

25、rce：场源磁化强度的一维数组!Az_source：磁化方向与Z轴夹角的一维数组!Ax_source：磁化方向与X轴夹角的一维数组!Ay_source：磁化方向与Y轴夹角的一维数组!X_source：X坐标的起点和终点的二维数组!Y_source：Y坐标的起点和终点的二维数组!Z_source：Z坐标的起点和终点的二维数组!*subroutine read_source(filename_source,N_source,density_source,mag_source,Az_source,Ax_source,Ay_source,X_source,Y_source,Z_source)impli

26、cit nonecharacter *(*) filename_sourceintegerN_sourcereal 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)integer kopen(10,file=filename_source,status='ol

27、d')do k=1,N_source,1read(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)enddoclose(10)end subroutine read_source!*! 子程序：read_N_coord!功能：读取计算点个数!输入参数：!filename：文件名!输出参数：!m：每条线上的点数（规则网）

28、!n：测线数（规则网）!N_coord：计算点个数!*subroutine read_N_coord(filename,m,n,N_coord,type_net)implicit nonecharacter *(*) filenameinteger m,n,N_coordinteger type_netif(type_net.eq.2).or.(type_net.eq.4) thenopen(10,file=filename,status='old')read(10,*)read(10,*) m,nclose(10)N_coord=m*nelse if(type_net.eq.

29、1).or.(type_net.eq.3) then open(10,file=filename,status='old') N_coord=0 do while(.not.eof(10)read(10,*,end=100,ERR=100)N_coord=N_coord+1100 end doclose(10)end ifwrite(*,*)m,n,N_coordpauseend subroutine read_N_coord!*! 子程序：read_coordinate!功能：读取计算点数据!输入参数：!filename：文件名!m：每条线上的点数（规则网）!n：测线数（规则

30、网）!N_coord：计算点数!constant_z：当测网类型是平面网时Z坐标的值!type_net：测网类型（1：曲面非规则网；2：曲面规则网；3：平面非规则网；4：平面规则网）!输出参数：!x_coord：计算点x坐标的一维数组!y_coord：计算点y坐标的一维数组!z_coord：计算点z坐标的一维数组!*subroutine read_coordinate(filename_obser,N_point,N_line,N_coord,X_coord,Y_coord,z_coord,constant_z,Xmin,Xmax,Ymin,Ymax,type_net)implicit non

31、echaracter*(*)filename_obserinteger type_netinteger N_coord,N_point,N_linereal Xmin,Xmax,Ymin,Ymax,constant_zreal X_coord(1:N_coord),Y_coord(1:N_coord),z_coord(1:N_coord)real i,j,kif (type_net=1) thenopen(12,file=filename_obser,status='old') do i=1,N_coord,1 read(12,*)X_coord(i),Y_coord(i),z

32、_coord(i) enddoclose(12) else if(type_net=3) then open(12,file=filename_obser,status='old') do i=1,N_coord,1 read(12,*)X_coord(i),Y_coord(i) enddo do i=1,N_coord,1 z_coord(i)=constant_z enddo close(12) elseif(type_net=2) then open(10,file=filename_obser,status='old') read(10,*) read(

33、10,*) read(10,*) xmin,xmax read(10,*) ymin,ymax read(10,*) write(*,*)'haha' write(*,*)xmin,xmax,ymin,ymax pause do J=1,N_line,1 read(10,*)(Z_coord(j-1)*N_point+i),i=1,N_point,1) enddo close(10) do j=1,N_line,1 do i=1,N_point,1 k=(j-1)*N_point+i X_coord(k)=Xmin+(j-1)*(Xmax-Xmin)/(N_point-1) Y

34、_coord(k)=Ymin+(i-1)*(Ymax-Ymin)/(N_line-1) enddo enddo elseif(type_net=4) then do j=1,N_line,1 do i=1,N_point,1 k=(j-1)*N_point+i X_coord(k)=Xmin+(j-1)*(Xmax-Xmin)/(N_point-1) Y_coord(k)=Ymin+(i-1)*(Ymax-Ymin)/(N_line-1) enddo enddo do i=1,N_coord,1 z_coord(i)=constant_z enddo !xmin=huge(xmin)!xmax

35、=-huge(xmax)!ymin=huge(ymin)!ymax=-huge(ymax)!do i=1,N_coord,1!xmin=min(xmin,x_coord(i)!xmax=max(xmax,x_coord(i)!ymin=min(ymin,y_coord(i)!ymax=max(ymax,y_coord(i)!end doend ifend subroutine read_coordinate!*! 子程序：calculate_gra!功能：计算场源在计算点处产生的重力异常!输入参数：!N_source：场源点个数!X_source：场源点X坐标的一维数组!Y_source：场源

36、点Y坐标的一维数组!Z_source：场源点Z坐标的一维数组!density_source：场源点剩余密度!N_coord：计算点数!X_coord：计算点X坐标的一维数组!Y_coord：计算点Y坐标的一维数组!Z_coord：计算点Z坐标的一维数组!输出参数：!field：计算点的重力异常的一维数组!*subroutine calculate_gra(N_source,X_source,Y_source,Z_source,density_source,N_coord,X_coord,Y_coord,Z_coord,field)implicit noneintegerN_sourcereal

37、X_source(1:N_source,1:2),Y_source(1:N_source,1:2),Z_source(1:N_source,1:2),density_source(1:N_source)integer N_coordrealX_coord(1:N_coord),Y_coord(1:N_coord),Z_coord(1:N_coord),field(1:N_coord)integer l,N,i,j,kreal delt_gdo l=1,N_coord,1 field(l)=0.0 do N=1,N_source,1 do i=1,2,1 do j=1,2,1do k=1,2,1

38、field(l)=field(l)+(-1)*(i+j+k)*delt_g(density_source(N),X_source(N,i),Y_source(N,j),Z_source(N,k),X_coord(l),Y_coord(l),Z_coord(l)end do end do end do end doend doend subroutine calculate_gra!*! 外部函数：delt_g!功能：计算单个直立六面体所产生的重力异常delt_g!输入参数说明：!density_source：场源点剩余密度!X12：场源点的X坐标!Y12：场源点的Y坐标!Z12：场源点的Z坐标

39、!X：计算点的X坐标!Y：计算点的Y坐标!Z：计算点的Z坐标!G：万有引力常数!*real function delt_g(density_source,X12,Y12,Z12,X,Y,Z)implicit nonereal density_source,X12,Y12,Z12real X,Y,Zreal G,rG=66.72r=sqrt(X12-X)*2+(Y12-Y)*2+(Z12-Z)*2)delt_g=G*density_source*(X12-X)*log(Y12-Y)+r)+(Y12-Y)*log(X12-X)+r)-(Z12-Z)*ATAN2(X12-X)*(Y12-Y),(Z1

40、2-Z)*r)*(-1)end function delt_g!*! 子程序：calculate_mag!功能：计算场源在计算点处产生的磁异常或化极磁力异常!输入参数：!N_source：场源个数!mag_source：场源磁化强度!Az_source：磁化方向与Z轴夹角的一维数组!Ax_source：磁化方向与X轴夹角的一维数组!Ay_source：磁化方向与Y轴夹角的一维数组!X_source：X坐标的起点和终点的二维数组!Y_source：Y坐标的起点和终点的二维数组!Z_source：Z坐标的起点和终点的二维数组!N_coord：计算点个数!x_coord：计算点的x坐标的一维数组!y

41、_coord：计算点的y坐标的一维数组!z_coord：计算点的z坐标的一维数组!anglez：t0方向与Z轴夹角!anglex：t0方向与X轴夹角!angley：t0方向与Y轴夹角!输出参数：!field：计算点的磁力异常的一维数组!*subroutine calculate_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)implicit noneinteger

42、 N_sourcereal 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)integer N_coordreal X_coord(1:N_coord),Y_coord(1:N_coord),Z_coord(1:N_coord),field(1:N_coord)real anglez,anglex,angleyinteg

43、er i,j,k,l,nreal Mx(1:N_source),My(1:N_source),Mz(1:N_source)real Hax(1:N_coord),Hay(1:N_coord),Za(1:N_coord)real delt_Hax,delt_Hay,delt_Za,delt_x,delt_y,delt_z,rdo i=1,N_source,1Mx(i)=mag_source(i)*sind(Az_source(i)*cosd(Ax_source(i)My(i)=mag_source(i)*sind(Az_source(i)*cosd(Ay_source(i)Mz(i)=mag_s

44、ource(i)*cosd(Az_source(i)end dodo l=1,N_coord,1Hax(l)=0.0Hay(l)=0.0Za(l)=0.0do n=1,N_source,1do i=1,2,1do j=1,2,1do k=1,2,1delt_x=X_source(n,i)-X_coord(l)delt_y=Y_source(n,j)-Y_coord(l)delt_z=Z_source(n,j)-Z_coord(l)r=sqrt(delt_x*2+delt_y*2+delt_z*2)delt_Hax=(-Mx(n)*atan2(delt_y*delt_z,r*delt_x)+My(n)*log(r+delt_z)+Mz(n)*log(r+delt_y)*0.1delt_Hay=(Mx(n)*log(r+delt_z)-My(n)*atan2(delt_x*delt_z,r*delt_y)+Mz(n)*log(r+delt_z)*0.1delt_Za=(Mx(n)*log(r+delt_y)+My(n)*log(r+delt_x)-Mz(n)*at