重磁实验报告（三）复杂形体正演姓名：学号：专业：勘查技术与工程指导教师：鲁宝亮、王万银完成日期：2013年12月18日目录一、基本原理： (1)1、重力异常计算公式 (1)2、磁力异常计算公式 (2)3、化极磁力异常计算公式 (3)二、输入/输出数据格式设计： (3)1、输入数据文件名的格式设计 (3)2、输出数据文件名的格式设计 (4)3、重要变量的名称 (4)三、总体设计： (6)四、测试结果： (7)五、结论及建议： (9)附录：源程序代码 (9)一、基本原理：（1Z地质体重力异常的计算 重力异常计算公式：⎰⎰⎰-+-+--==∂∂=∆vz z y x d d d z G V z V g 2/3222])()()[()(ζηξζηξζσ式中：V 为地质体的剩余质量对A 点的单位质量产生的引力位；v 为地质体的体积。
我们还可以推导出计算重力异常垂向梯度或重力垂向梯度异常的基本公式为：ζηξζηξηξζσd d d z y x y x z G V z g vzz ⎰⎰⎰-+-+------==∂∆∂2/5222222])()()[()()()(2 计算重力异常水平梯度或重力水平梯度异常的基本公式为：ζηξζηξξζσd d d z y x x z G V x g v xz ⎰⎰⎰-+-+---==∂∆∂2/5222])()()[())((3 ζηξζηξηζσd d d z y x y z G V y g vyz ⎰⎰⎰-+-+---==∂∆∂2/5222])()()[())((3 一阶导数类边缘识别计算公式： 垂向导数：zz y x g z y x VDR ∂∆∂=),,(),,(总水平导数：22)),,(()),,((),,(yz y x g x z y x g z y x THDR ∂∆∂+∂∆∂= 解析信号振幅：22VDR THDR ASM +=（2） O X •(x,y,z)a b Y （000,,z y x ） c(ζηξ,,)Z直立长方体示意图均匀磁化直立长方体磁场各分量表达式：222222220000000)]}(ln[)](ln[)()()())((arctan {4),,(cz cz b y b y a x a x z y x ax y r M z r M z z r x y x M z y x H +-+-+--++-++-+-+----=ηζζζξηξπμ （1）222222220000000)]}(ln[)()()())((arctan )](ln[{4),,(cz c z b y b y a x a x z y x ay x r M z z r y y x M z r M z y x H +-+-+--++-+-+-----+=ξζζηηξζπμ （2）2222220000000})())((arctan )](ln[)](ln[{4),,(cz c z b y b y ax ax Z y x a z r y x M x r M y r M z y x Z +-+-+------++-+=ζηξξηπμ （3）把（1）、（2）、（3）式代入磁异常T ∆计算公式中得：22222262252243210000000})())((arctan )()()())((arctan )()()())((arctan )](ln[)](ln[)](ln[{4),,(cz c z b y b y ax ax z r y x k z z r y y x k z z r x y x k z r k y r k x r k M z y x T +-+-+----•+-+-+---•+-+-+---•+-+•+-+•+-+•=∆ζηξζζηηξζζξηξζηξπμ 其中21222])()()[(z y x r -+-+-=ζηξ，ζηξ,,为异常体的坐标，、 、为地磁场的方向余弦；为总磁化强度，、、为其在X 、Y 、Z 坐标方向的分量。
设M 为的模，I 为磁化强度倾角，D 为磁化强度磁偏角,且满足关系式D I M M x cos cos = D I M M y sin cos = I M M z sin =（3）化极磁力异常磁极因子为：])([2])([22/1221112/12222212）（）（v u v u i v u v u i q q l l ++++++=γβαπγβαπ 式中：1α、1β、1γ为原磁化方向的方向余弦；222γβα、、为新磁化方向的方向余弦（磁化方向即为磁化强度方向）。
该公式适合于转换到任意磁化方向。
二、输入/输出数据格式设计：1、输入数据文件名格式设计：（1）输入数据及文件名对应如下： input_source_filename ：source.dat input_plane_terrain_filename ：AXYZL.GRD input_rugged_terrain_filename ：BXYZL.GRD body_num ：3input_source_filename 对应source.dat 文件，存放直立六面体相关参数；input_plane_terrain_filename 对应AXYZL.GRD 文件，存放平面地形数据；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 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个点；点（线）坐标从-26m～26m，点（线）距为2m；Zmin ，Zmax分别表示Z坐标的最小值与最大值。
2、输出数据文件名格式设计：（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_pole_magnetic_anomaly_filename :plane_pole_magnetic_anomaly.grd rugged_surface_pole_magnetic_anomaly_filename : rugged_surface_pole_m agnetic_anomaly.grd （2）对给定三维规则地质体，计算其在给定平面上的重力异常plane_gravit y_anomaly，曲面上的重力异常rugged_surface_gravity_anomaly；平面上的磁力异常plane_magnetic_anomaly，曲面上的磁力异常rugged_surface_magneti c_anomaly；平面化极磁力异常plane_pole_magnetic_anomaly，曲面化极磁力异常。
输出格式如下：DSAA27 27-26 26-26 26field_min field_max 可知，网格化数据按GRD格式保存，有27条线、每条线27个点；点（线）坐标从-26m～26m，点（线）距为2m，field min ，field max分别表示网格化数据的最小值与最大值。
3、重要变量的名称：INTEGER m0,m1,m2,m3 m0,m3为扩边后的起点和终点点位m1,m2为原始数据的起点和终点点位INTEGER n0,n1,n2,n3 n0,n3为扩边后的起点和终点线位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_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_anomaly 曲面磁力异常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 新磁化方向方向余弦的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 线距三、总体设计：A 段流程图：四、测试结果：平面重力异常曲面重力异常平面磁力异常曲面磁力异常平面化极磁力异常曲面化极磁力异常五、结论及建议：结论：（1）在计算重力和磁力异常时，我们使用了三个直立六面体，在用数组存放六面体的X、Y、Z坐标起点与终点时，我们应该使用三个二维数组，这样在编程时才更便于计算场点到源点的距离R。