长安大学重磁报告三

1、重磁实验报告（三）复杂形体正演 姓名： 学号： 专业：勘查技术与工程 指导教师：鲁宝亮、王万银 完成日期：2013年12月18日1目录一、基本原理：1 1、重力异常计算公式12、磁力异常计算公式23、化极磁力异常计算公式3二、输入/输出数据格式设计：31、输入数据文件名的格式设计32、输出数据文件名的格式设计43、重要变量的名称4三、总体设计：6四、测试结果：7五、结论及建议：9附录：源程序代码9 2一、基本原理：（1） O X A(x,y,z) Y dv() Z 地质体重力异常的计算 重力异常计算公式： 式中：V为地质体的剩余质量对A点的单位质量产生的引力位；v为地质体的体积。我们还可以推导

2、出计算重力异常垂向梯度或重力垂向梯度异常的基本公式为：计算重力异常水平梯度或重力水平梯度异常的基本公式为：一阶导数类边缘识别计算公式：垂向导数：总水平导数：解析信号振幅：（2） O X (x,y,z) a b Y （） c () Z 直立长方体示意图均匀磁化直立长方体磁场各分量表达式： （1） （2） （3） 把（1）、（2）、（3）式代入磁异常计算公式中得： 其中，为异常体的坐标，、 、为地磁场的方向余弦；为总磁化强度，、为其在X、Y、Z坐标方向的分量。设M为的模，I为磁化强度倾角，D为磁化强度磁偏角,且满足关系式 （3）化极磁力异常磁极因子为： 式中：、为原磁化方向的方向余弦；为新磁化方向

3、的方向余弦（磁化方向即为磁化强度方向）。该公式适合于转换到任意磁化方向。二、输入/输出数据格式设计： 1、输入数据文件名格式设计： （1）输入数据及文件名对应如下： input_source_filename ：source.dat input_plane_terrain_filename ：AXYZL.GRD input_rugged_terrain_filename ：BXYZL.GRD body_num ：3 input_source_filename对应source.dat文件，存放直立六面体相关参数；input_plane_terrain_filename对应AXYZL.GRD文件，

4、存放平面地形数据；input_rugged_terrain_filename对应BXYZL.GRD文件，存放曲面地形数据；body_num赋值为3，表示地质体数量。 （2）场源参数保存在“source.dat”中。第一列为剩余密度(g/cm3)；第二列为磁化强度(10-6CGSM)；第三列为磁化方向倾角(DEG)；第四列为磁化方向与x轴的夹角(DEG)；第五列为磁化方向与y轴的夹角(DEG)；第六列第七列为x坐标的起点和终点(km)；第八列第九列为y坐标的起点和终点(km)；第十列第十一列为z坐标的起点和终点(km，向下为正)。 0.7,20000,50,85,5,-9,-3,4,8,2,7

5、0.8,34000,50,85,5,0,7,0,5,3,7 0.9,17000,50,85,5,-4,8,-8,-3,2,7 （3）起伏地形数据保存在“bxyzl.grd”文件中（GRD格式）；平面坐标范围与曲面坐标范围相同，但其z坐标值为-5.3km，形式如下： DSAA 27 27 -26 26 -26 26 -5.3 1.7 0.15 0 -0.15 -0.3 -0.4 -0.5 -0.6 -0.7 -0.85 可知，网格化数据按GRD格式保存，有27条线、每条线27个点；点（线）坐标从-26m26m，点（线）距为2m；Zmin ，Zmax分别表示Z坐标的最小值与最大值。 2、输出数据文

6、件名格式设计： （1）输出数据及文件名对应如下：plane_gravity_anomaly_filename ：plane_gravity_anomaly.grdplane_magnetic_anomaly_filename ：plane_magnetic_anomaly.grdrugged_surface_gravity_anomaly_filename ：rugged_surface_gravity_anom aly.grdrugged_surface_magnetic_anomaly_filename ：rugged_surface_magnetic_an omaly.grdplane_

7、pole_magnetic_anomaly_filename :plane_pole_magnetic_anomaly.grdrugged_surface_pole_magnetic_anomaly_filename : rugged_surface_pole_m agnetic_anomaly.grd （2）对给定三维规则地质体，计算其在给定平面上的重力异常plane_gravity_anomaly，曲面上的重力异常rugged_surface_gravity_anomaly；平面上的磁力异常plane_magnetic_anomaly，曲面上的磁力异常rugged_surface_magn

8、etic_anomaly；平面化极磁力异常plane_pole_magnetic_anomaly，曲面化极磁力异常。输出格式如下： DSAA 27 27 -26 26 -26 26 field_min field_max 可知，网格化数据按GRD格式保存，有27条线、每条线27个点；点（线）坐标从-26m26m，点（线）距为2m，field min ，field max分别表示网格化数据的最小值与最大值。 3、重要变量的名称： INTEGER m0,m1,m2,m3 m0,m3为扩边后的起点和终点点位 m1,m2为原始数据的起点和终点点位 INTEGER n0,n1,n2,n3 n0,n3为扩

9、边后的起点和终点线位 n1,n2为原始数据的起点和终点线位 INTEGER mpoint,nline mpoint表示点数 nline表示线数 REAL remaining_density 剩余密度(g/cm3) REAL magnetic_intensity 磁化强度(10-6CGSM) REAL magnetic_inclination 磁化方向倾角(DEG) REAL magnetic_x_angle 磁化方向与x轴的夹角(DEG) REAL magnetic_y_angle 磁化方向与y轴的夹角(DEG) REAL X_source 二维数组，x坐标的起点和终点(km) REAL Y_

10、source 二维数组，y坐标的起点和终点(km) REAL Z_source 二维数组，z坐标的起点和终点(km） REAL Z_coordinate 二维数组，GRD数据的Z坐标 REAL plane_terrain_data 平面地形数据 REAL rugged_terrain_data 曲面地形数据 REAL plane_gravity_anomaly 平面重力异常 REAL plane_magnetic_anomaly 平面磁力异常 REAL rugged_surface_gravity_anomaly 曲面重力异常 REAL rugged_surface_magnetic_anom

11、aly 曲面磁力异常 REAL plane_pole_magnetic_anomaly 平面化极磁力异常 REAL rugged_surface_pole_magnetic_anomaly 曲面化极磁力异常 REAL x_original_unit_vector 原磁化方向方向余弦的X坐标 REAL y_original_unit_vector 原磁化方向方向余弦的Y坐标 REAL z_original_unit_vector 原磁化方向方向余弦的Z坐标 REAL x_new_unit_vector 新磁化方向方向余弦的X坐标 REAL x_new_unit_vector 新磁化方向方向余弦的

12、Y坐标 REAL x_new_unit_vector 新磁化方向方向余弦的Z坐标 REAL expand_rugged_surface_magnetic_anomaly(:,:) 余弦扩边后的曲面磁力异常场 REAL field_real(:,:) FFT变换后磁力异常场的实部 REAL field_image(:,:) FFT变换后磁力异常场的虚部 REAL pole_factor_real(:,:) 磁极因子实部 REAL pole_factor_image(:,:) 磁极因子虚部 dx 点距 dy 线距3、 总体设计：读取场源数据 读取文件名与文件参数定义变量参数，定义输入输出文件名 读

13、取地形数据 开始 计算平面重力异常A段流程图： 开始 计算曲面重力异常定义变量参数，定义输入输出文件名 计算平面磁力异常 计算曲面磁力异常 二维FFT正变换 结束 输出平面重力异常 计算磁极因子输出平面化极磁力异常输出曲面重力异常 乘积运算 A输出平面化极磁力异常 二维FFT反变换输出平面磁力异常 二维余弦扩边 结束 A 输出曲面磁力异常 读入场值 确定扩边点位 确定扩边点位4、 测试结果： 平面重力异常 曲面重力异常 平面磁力异常 曲面磁力异常 平面化极磁力异常 曲面化极磁力异常五、结论及建议：结论： （1）在计算重力和磁力异常时，我们使用了三个直立六面体，在用数组存放六面体的X、Y、Z坐标

14、起点与终点时，我们应该使用三个二维数组，这样在编程时才更便于计算场点到源点的距离R。 （2）平面重力异常与曲面重力异常相比较，等位线显得更加圆滑；曲面重力异常更能突出局部异常。平面磁力异常与曲面磁力异常相比较，异常的范围更大；曲面磁力异常更能突出局部异常。 （3）对磁力异常化极后，我们可以清楚地看到高值异常周围有四块低值区域，而且分布比较对称；曲面化极磁力异常与平面化极磁力异常相比，有三个分开的高值区域，更能清楚的突出三个异常体的存在。建议： 我们每次做完实验，交完报告后，都不会得到任何来自老师的反馈，那我们就很难知道自己在写实验报告，编程序过程中的不足，建议老师能指出我们一些实验方面的不足和

15、改正的办法。附录：源程序代码!复杂形体正演程序设计!-PROGRAM COMPLEX_BODY_FORWARDINTEGER mpoint,nlineINTEGER body_numREAL xmax,xmin,ymax,ymin,zmax,zmin,A,B,C,eigvalCHARACTER *120 input_source_filename,input_rugged_terrain_filenameCHARACTER *120 plane_gravity_anomaly_filename,plane_magnetic_anomaly_filename,rugged_surface_gra

16、vity_anomaly_filenameCHARACTER *120 rugged_surface_magnetic_anomaly_filenameCHARACTER *120 plane_pole_magnetic_anomaly_filename,rugged_surface_pole_magnetic_anomaly_filenameREAL,ALLOCATABLE : remaining_density(:)REAL,ALLOCATABLE : magnetic_intensity(:)REAL,ALLOCATABLE : magnetic_inclination(:)REAL,A

17、LLOCATABLE : magnetic_x_angle(:)REAL,ALLOCATABLE : magnetic_y_angle(:)REAL,ALLOCATABLE : X_source(:,:)REAL,ALLOCATABLE : Y_source(:,:)REAL,ALLOCATABLE : Z_source(:,:)REAL,ALLOCATABLE : Z_coordinate(:,:)REAL,ALLOCATABLE : plane_terrain_data(:,:)REAL,ALLOCATABLE : rugged_terrain_data(:,:)REAL,ALLOCATA

18、BLE : plane_gravity_anomaly(:,:)REAL,ALLOCATABLE : plane_magnetic_anomaly(:,:)REAL,ALLOCATABLE : rugged_surface_gravity_anomaly(:,:)REAL,ALLOCATABLE : rugged_surface_magnetic_anomaly(:,:)REAL,ALLOCATABLE : plane_pole_magnetic_anomaly(:,:)REAL,ALLOCATABLE : rugged_surface_pole_magnetic_anomaly(:,:)RE

19、AL,ALLOCATABLE : field_real(:,:)REAL,ALLOCATABLE : field_image(:,:)REAL,ALLOCATABLE : pole_factor_real(:,:)REAL,ALLOCATABLE : pole_factor_image(:,:)REAL,ALLOCATABLE : random_terrain(:,:)call read_parameter(input_source_filename,input_rugged_terrain_filename,plane_gravity_anomaly_filename,& plane

20、_magnetic_anomaly_filename,rugged_surface_gravity_anomaly_filename,rugged_surface_magnetic_anomaly_filename,& plane_pole_magnetic_anomaly_filename,rugged_surface_pole_magnetic_anomaly_filename,body_num) !读取文件名与文件参数 call read_source(input_source_filename,body_num,remaining_density,magnetic_intens

21、ity,magnetic_inclination,magnetic_x_angle,magnetic_y_angle,& X_source,Y_source,Z_source) !读取场源数据call read_parameter_grd(input_rugged_terrain_filename,mpoint,nline,Xmin,Xmax,Ymin,Ymax,Zmin,Zmax)call read_terrain_grd(input_rugged_terrain_filename,mpoint,nline,Z_coordinate) !读取曲面地形数据dx=(xmax-xmin)/

22、(mpoint-1) !确定间隔dy=(ymax-ymin)/(nline-1) ALLOCATE (remaining_density(1:body_num)ALLOCATE (magnetic_intensity(1:body_num)ALLOCATE (magnetic_inclination(1:body_num)ALLOCATE (magnetic_x_angle(1:body_num)ALLOCATE (magnetic_y_angle(1:body_num)ALLOCATE (X_source(1:2,1:body_num)ALLOCATE (Y_source(1:2,1:bod

23、y_num)ALLOCATE (Z_source(1:2,1:body_num)ALLOCATE (Z_coordinate(1:mpoint,1:nline)ALLOCATE (plane_terrain_data(1:mpoint,1:nline)ALLOCATE (rugged_terrain_data(1:mpoint,1:nline)ALLOCATE (plane_gravity_anomaly(1:mpoint,1:nline)ALLOCATE (plane_magnetic_anomaly(1:mpoint,1:nline)ALLOCATE (rugged_surface_gra

24、vity_anomaly(1:mpoint,1:nline)ALLOCATE (rugged_surface_magnetic_anomaly(1:mpoint,1:nline)plane_terrain_data=Zminrugged_terrain_data=Z_coordinatecall gravity_anomaly_calculate(remaining_density,X_source,Y_source,Z_source,Xmin,Xmax,Ymin,Ymax,plane_terrain_data,mpoint,nline,body_num,& plane_gravity

25、_anomaly) !计算平面重力异常call gravity_anomaly_calculate(remaining_density,X_source,Y_source,Z_source,Xmin,Xmax,Ymin,Ymax,rugged_terrain_data,mpoint,nline,body_num,& rugged_surface_gravity_anomaly) !计算曲面重力异常call magnetic_anomaly_calculate(magnetic_intensity,magnetic_inclination,magnetic_x_angle,magneti

26、c_y_angle,X_source,Y_source,Z_source,& Xmin,Xmax,Ymin,Ymax,plane_terrain_data,mpoint,nline,body_num,plane_magnetic_anomaly,A,B,C) !计算平面磁力异常 call magnetic_anomaly_calculate(magnetic_intensity,magnetic_inclination,magnetic_x_angle,magnetic_y_angle,X_source,Y_source,Z_source,& Xmin,Xmax,Ymin,Ym

27、ax,rugged_terrain_data,mpoint,nline,body_num,rugged_surface_magnetic_anomaly,A,B,C) !计算曲面磁力异常eigval=1.701411e+38call OUTPUT_grd(plane_gravity_anomaly,plane_gravity_anomaly_filename,1,mpoint,1,nline,mpoint,nline,eigval,xmin,xmax,ymin,ymax) !输出平面重力异常call output_grd(rugged_surface_gravity_anomaly,rugge

28、d_surface_gravity_anomaly_filename,1,mpoint,1,nline,mpoint,nline,eigval,xmin,xmax,ymin,ymax) !输出曲面重力异常 call output_grd(plane_magnetic_anomaly,plane_magnetic_anomaly_filename,1,mpoint,1,nline,mpoint,nline,eigval,xmin,xmax,ymin,ymax) !输出平面磁力异常call output_grd(rugged_surface_magnetic_anomaly,rugged_surf

29、ace_magnetic_anomaly_filename,1,mpoint,1,nline,mpoint,nline,eigval,xmin,xmax,ymin,ymax) !输出曲面磁力异常call EX_EDGE(mpoint,m0,m1,m2,m3,1) !确定扩边点位call EX_EDGE(nline,n0,n1,n2,n3,1)ALLOCATE (field_real(m0:m3,n0:n3)ALLOCATE (field_image(m0:m3,n0:n3)ALLOCATE (pole_factor_real(m0:m3,n0:n3)ALLOCATE (pole_factor_

30、image(m0:m3,n0:n3)ALLOCATE (random_terrain(m0:m3,n0:n3)ALLOCATE (plane_pole_magnetic_anomaly(m0:m3,n0:n3)ALLOCATE (rugged_surface_pole_magnetic_anomaly(m0:m3,n0:n3)call Read_field(plane_magnetic_anomaly_filename,m0,m1,m2,m3,n0,n1,n2,n3,plane_magnetic_anomaly) !读入场值call Read_field(rugged_surface_magn

31、etic_anomaly_filename,m0,m1,m2,m3,n0,n1,n2,n3,rugged_surface_magnetic_anomaly) call EXP_2D_cos_sub(m1,m2,m3,m4,n1,n2,n3,n4,plane_magnetic_anomaly,1) !二维余弦扩边call EXP_2D_cos_sub(m1,m2,m3,m4,n1,n2,n3,n4,rugged_surface_magnetic_anomaly,1)!平面化极磁力异常field_real=plane_magnetic_anomaly !确定实部虚部field_image=0.0c

32、all random_pole_calculate(Field_Real,Field_Image,m0,m3,n0,n3,dx,dy,A,B,C,0.0,0.0,1.0,pole_factor_real,pole_factor_image,random_terrain)plane_pole_magnetic_anomaly=random_terraincall output_grd_sub(plane_pole_magnetic_anomaly_filename,mpoint,nline,Xmin,Xmax,Ymin,Ymax,m0,m1,m2,m3,n0,n1,n2,n3,plane_pol

33、e_magnetic_anomaly) !曲面化极磁力异常field_real=rugged_surface_magnetic_anomaly !确定实部虚部field_image=0.0call random_pole_calculate(Field_Real,Field_Image,m0,m3,n0,n3,dx,dy,A,B,C,0.0,0.0,1.0,pole_factor_real,pole_factor_image,random_terrain)rugged_surface_pole_magnetic_anomaly=random_terraincall output_grd_sub

34、(rugged_surface_pole_magnetic_anomaly_filename,mpoint,nline,Xmin,Xmax,Ymin,Ymax,m0,m1,m2,m3,n0,n1,n2,n3,& rugged_surface_pole_magnetic_anomaly)DEALLOCATE (remaining_density)DEALLOCATE (magnetic_intensity)DEALLOCATE (magnetic_inclination)DEALLOCATE (magnetic_x_angle)DEALLOCATE (magnetic_y_angle)D

35、EALLOCATE (X_source)DEALLOCATE (Y_source)DEALLOCATE (Z_source)DEALLOCATE (Z_coordinate)DEALLOCATE (plane_terrain_data)DEALLOCATE (rugged_terrain_data)DEALLOCATE (plane_gravity_anomaly)DEALLOCATE (plane_magnetic_anomaly)DEALLOCATE (rugged_surface_gravity_anomaly)DEALLOCATE (rugged_surface_magnetic_an

36、omaly) DEALLOCATE (plane_pole_magnetic_anomaly) DEALLOCATE (rugged_surface_pole_magnetic_anomaly) DEALLOCATE (field_real)DEALLOCATE (field_image)DEALLOCATE (pole_factor_real)DEALLOCATE (pole_factor_image)DEALLOCATE (random_terrain)END PROGRAM COMPLEX_BODY_FORWARD!-读取文件名与文件参数-SUBROUTINE read_paramete

37、r(input_source_filename,input_rugged_terrain_filename,plane_gravity_anomaly_filename,& plane_magnetic_anomaly_filename,rugged_surface_gravity_anomaly_filename,rugged_surface_magnetic_anomaly_filename,& plane_pole_magnetic_anomaly_filename,rugged_surface_pole_magnetic_anomaly_filename,body_nu

38、m) CHARACTER *(*) input_source_filename,input_rugged_terrain_filename CHARACTER *(*) plane_gravity_anomaly_filename,plane_magnetic_anomaly_filename,rugged_surface_gravity_anomaly_filename,& rugged_surface_magnetic_anomaly_filename CHARACTER *(*) plane_pole_magnetic_anomaly_filename,rugged_surfac

39、e_pole_magnetic_anomaly_filename INTEGER body_num OPEN(11,file='parameter.txt',status='old') READ(11,*) input_source_filename READ(11,*) input_rugged_terrain_filenameREAD(11,*) plane_gravity_anomaly_filename READ(11,*) plane_magnetic_anomaly_filename READ(11,*) rugged_surface_gravity

40、_anomaly_filename READ(11,*) rugged_surface_magnetic_anomaly_filename READ(11,*) plane_pole_magnetic_anomaly_filename READ(11,*) rugged_surface_pole_magnetic_anomaly_filename READ(11,*) body_numCLOSE(11)END SUBROUTINE read_parameter!-读取场源数据子程序开始-SUBROUTINE read_source(filename,body_num,remaining_den

41、sity,magnetic_intensity,magnetic_inclination,magnetic_x_angle,magnetic_y_angle,& X_source,Y_source,Z_source)INTEGER body_num CHARACTER *(*) filenameREAL remaining_density(1:body_num),magnetic_intensity(1:body_num),magnetic_inclination(1:body_num),magnetic_x_angle(1:body_num),& magnetic_y_ang

42、le(1:body_num),X_source(1:2,1:body_num),Y_source(1:2,1:body_num),Z_source(1:2,1:body_num)OPEN (20,file=filename,status='old') DO i=1,body_num READ (20,*) remaining_density(i),magnetic_intensity(i),magnetic_inclination(i),magnetic_x_angle(i),magnetic_y_angle(i),& X_source(1,i),X_source(2,

43、i),Y_source(1,i),Y_source(2,i),Z_source(1,i),Z_source(2,i) END DOCLOSE(20)END SUBROUTINE read_source!-读取场源数据子程序结束-!-读取地形数据子程序开始-SUBROUTINE read_parameter_grd(filename,mpoint,nline,Xmin,Xmax,Ymin,Ymax,Zmin,Zmax)CHARACTER *(*) filenameINTEGER mpoint,nlineREAL Xmin,Xmax,Ymin,Ymax,Zmin,ZmaxOPEN (30,file

44、=filename,status='old') READ(30,*) READ(30,*) mpoint,nline READ(30,*) Xmin,Xmax READ(30,*) Ymin,Ymax READ(30,*) Zmin,ZmaxCLOSE(30)END SUBROUTINE read_parameter_grdSUBROUTINE read_terrain_grd(filename,mpoint,nline,Z_coordinate)CHARACTER *(*) filenameINTEGER mpoint,nlineREAL Z_coordinate(1:mpo

45、int,1:nline)OPEN (40,file=filename,status='old') DO j=1,nline DO i=1,mpoint READ(40,*) Z_coordinate(i,j)END DO END DOCLOSE(40)END SUBROUTINE read_terrain_grd!-读取地形数据子程序结束-!-计算重力异常子程序开始-SUBROUTINE gravity_anomaly_calculate(remaining_density,X_source,Y_source,Z_source,Xmin,Xmax,Ymin,Ymax,Z_coo

46、rdinate,mpoint,& nline,body_num,gravity_anomaly)IMPLICIT NONEREAL(8) GPARAMETER(G=6.672e-11)INTEGER mpoint,nline,body_num,I,J,K,L,M,NREAL remaining_density(1:body_num),X_source(1:2,1:body_num),Y_source(1:2,1:body_num),Z_source(1:2,1:body_num),& Z_coordinate(1:mpoint,1:nline),X_coordinate(1:mpoint),Y_coordinate(1:nline),gravity_anomaly(1:mpoint,1:nline)REAL Xmin,Xmax,Ymin,Ymax,sum,r,dx,dygravity_anomaly=0.0dx=(xmax-xmin)/(m