直立六面体正演程序设计实验

长安大学长安大学 重力与磁法勘探实验报告重力与磁法勘探实验报告 直立六面体正演程序设计直立六面体正演程序设计 一 实验内容 一 实验内容 对给定的密度体和磁性体计算平面规则网上的重力异常 磁力异常 二 基本原理二 基本原理 1 1 模型示意及有关变量的意义 并分别给出重力异常和磁力异常的计算 模型示意及有关变量的意义 并分别给出重力异常和磁力异常的计算 公式公式 计算磁力异常计算磁力异常 T T 的计算公式 的计算公式 222 111 0 22 22 1 ln 2ln 3 ln u 4arctan 4 5arctan6arctan krxkrykrz xy T X Y ZMk xrzz xyxy kk yrzzrz 其中其中 L0 M0 N0以及以及 分别为地磁场及总磁化强度的方向余弦 分别为地磁场及总磁化强度的方向余弦 I0 A 0 I A 分别分别 为地磁场与总磁化强度方向的倾角和相对为地磁场与总磁化强度方向的倾角和相对 x 轴方向的偏角 有以下关系 轴方向的偏角 有以下关系 00 1kMN 00 2kLN 00 3kLM 0 4kL 0 5kM 0 6kN 000 coscosLIA 000 cossinMIA 00 sinNI coscosIA cossinIA sin I 计算重力异常计算重力异常 g g 的计算公式 的计算公式 222 111 222 111 3 g ln ln arctan z z V x y zGd d d r Gxyr yxr xy z z r 算法原理 算法原理 222 222 k 111 111 1 ij ijk ijk x y z If x y zf x y z x y z 注明 在注明 在 CGSM 中的磁导率中的磁导率 u0 1 磁化强度的单位为奥斯特 磁化强度的单位为奥斯特 Oe 换算关系换算关系 1nT 10 5Oe 10 5Gs 而实验场源数据给出的是而实验场源数据给出的是 10 6Oe 为单位 经程序设计使得输出的磁异常为为单位 经程序设计使得输出的磁异常为 单位单位 nT 计算重力异常的使用的计算重力异常的使用的 G 为为 6 672 10 8cm3 g s2 的单位为的单位为 g cm3 而而 xyz 坐标的单位为坐标的单位为 km 则计算出来的重力异常的单位为 则计算出来的重力异常的单位为 km s2 它可它可 以换算成以换算成 g u 1km s2 109g u 以上的工作在程序中以实现 磁异常的输出单位为以上的工作在程序中以实现 磁异常的输出单位为 nt 重力异常输出的单 重力异常输出的单 位为位为 g u 2 2 计算点输入坐标格式设计 计算点输入坐标格式设计 我们用的是平面规则网我们用的是平面规则网 三三 总体设计 总体设计 定义变量参数 及 文件名 调用子程序读入所 需的数据的子程序 四 测试结果四 测试结果 测试参数如下测试参数如下 地磁场倾角为地磁场倾角为 50 50 与 与 x x 轴的夹角为轴的夹角为 85 85 与 与 y y 轴夹角为轴夹角为 5 5 平面坐标范围与曲面坐标范围相同 但其平面坐标范围与曲面坐标范围相同 但其 z z 坐标值为坐标值为 5 3km 5 3km 1 1 命名文件 命名文件 地磁参数地磁参数 txt txt 调用子程序依次计算每个场源的重力异常 磁力异常 并用 dat 格式输出 调用子程序叠加计算总的场 源的重力异常 磁力异常并 用 dat 格式输出 得到结果 输出结 果 调用格式输出子程序 将叠加重力异常和叠加磁 异常以 grd 格式输出 2 2 场源参数文件 场源参数文件 source txtsource txt 3 3 平面参数文件 平面平面参数文件 平面 xyz txtxyz txt 4 4 输出的叠加重力异常文件 输出的叠加重力异常文件 g g 叠加叠加 grd grd 5 5 输出的叠加磁力异常文件 输出的叠加磁力异常文件 t t 叠加叠加 grd grd 6 6 输出的输出的 g1 dat g2 dat g3 dat gd dat g1 dat g2 dat g3 dat gd dat t1 datt1 dat t2 dat t2 dat t3 datt3 dat td dattd dat 数据文件就不一一展示了 数据文件就不一一展示了 7 7 叠加后的重力异常 磁力异常由 叠加后的重力异常 磁力异常由 surfersurfer 得到的图形得到的图形 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 0 05 0 1 0 15 0 2 0 25 0 3 0 35 0 4 0 45 0 5 0 55 0 6 0 65 0 7 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 1 1 1 2 1 3 1 4 1 5 1 6 1 7 g1 场源 场源 1 0 7 0 5 3 7 异常图 异常图 g2 场源 场源 2 4 8 8 3 2 7 异常图 异常图 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 0 04 0 08 0 12 0 16 0 2 0 24 0 28 0 32 0 36 0 4 0 44 0 48 0 52 0 56 g3 场源 场源 3 9 3 4 8 2 7 异常图 异常图 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 16 12 8 4 0 4 8 12 16 20 24 28 32 36 40 44 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 20 15 10 5 0 5 10 15 20 25 30 35 40 45 50 t1 场源 场源 1 0 7 0 5 3 7 异常图 异常图 t2 场源 场源 2 4 8 8 3 2 7 异常图 异常图 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 10 8 6 4 2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 t3 场源 场源 3 9 3 4 8 2 7 异常图 异常图 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 2 1 2 2 g1 g2 g3g1 g2 g3 叠加重力异常叠加重力异常 25 20 15 10 50510152025 25 20 15 10 5 0 5 10 15 20 25 30 25 20 15 10 5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 t1 t2 t3t1 t2 t3 叠加磁异常叠加磁异常 五 结论及建议五 结论及建议 本次实验了解了在场源背景为直立六面体的正演重磁异常的方法 绘制本次实验了解了在场源背景为直立六面体的正演重磁异常的方法 绘制 出了地下异常的图像 从出了地下异常的图像 从 t1 t2 t3 g1 g2 g3t1 t2 t3 g1 g2 g3 异常图中可看出 异常大的分异常图中可看出 异常大的分 布在场源的正上方 但三场源的重磁异常经过叠加后 高异常值的范围的位布在场源的正上方 但三场源的重磁异常经过叠加后 高异常值的范围的位 置改变了 磁异常有正负可以在以后的不规则体 无法用公式直接给出的重置改变了 磁异常有正负可以在以后的不规则体 无法用公式直接给出的重 磁异常 中磁异常的正演中 把不规则提划分成许多个直立六面体 然后运磁异常 中磁异常的正演中 把不规则提划分成许多个直立六面体 然后运 用此程序叠加计算该不规则体的产生的重磁异常 为不规则体的重磁异常正用此程序叠加计算该不规则体的产生的重磁异常 为不规则体的重磁异常正 演提供新的思路 对于异常的分布趋势以及为什么是这样分布的问题作了讨演提供新的思路 对于异常的分布趋势以及为什么是这样分布的问题作了讨 论 熟悉了物理场的单位换 加深了对地下异常分布的宏观认识 对于在课论 熟悉了物理场的单位换 加深了对地下异常分布的宏观认识 对于在课 本上的一些抽象的理解在结果图像中有了直观的认识 本上的一些抽象的理解在结果图像中有了直观的认识 编制程序中遇到的问题 编制程序中遇到的问题 1 1 物理场的单位换算 如不换算好的话 的出来的重磁异常都会为物理场的单位换算 如不换算好的话 的出来的重磁异常都会为 零 报告前面给出了换算思路 零 报告前面给出了换算思路 2 2 场源的读入问题 若有成千上百个场源 那又该怎么读入 场源的读入问题 若有成千上百个场源 那又该怎么读入 还有事先不知道文件里有多少个场源的话 又该怎么编制程序 还有事先不知道文件里有多少个场源的话 又该怎么编制程序 3 3 可以将平面拓展成曲面 即把在可以将平面拓展成曲面 即把在 zminzmin 和和 zmaxzmax 给定为不一样的值 给定为不一样的值 且且 zmin zmax zmin zmax 但这成为三维空间上某个点的重磁异常正演 用但这成为三维空间上某个点的重磁异常正演 用 surfersurfer 怎么怎么 画具有四维的重磁异常图像 画具有四维的重磁异常图像 4 4 将三重积分积出来的式子可以用三个循环累加的出结果 将三重积分积出来的式子可以用三个循环累加的出结果 附录 源程序代码 本程序按照自由格式编写 对于粘贴中难免附录 源程序代码 本程序按照自由格式编写 对于粘贴中难免 有一些问题 有一些问题 场源为直立六面体时的重力异常 dg 磁异常 dt 的正演 三个场源的数据存放在 source txt 中 地磁场参数 txt 存放地磁场的倾角 与 x 和 y 轴的夹角 平面 xyz txt 存放平面的参数 三个场源数据读到 s 3 11 中 t1 t2 t3 分别为三个场源产生的磁异常值存放的数组 g1 g2 g3 分别为三个场源产生的重力异常值存放的数组 td 和 gd 为三个场源叠加后产生的磁异常和重力异常存放的数组 is ix iy 和 is1 ix1 iy1 分别为地磁场和磁化强度的倾角 与 x 轴的夹角 与 y 轴的夹角 xmin xmax ymin ymax zmin zmax 为曲面的坐标 其中 zmin zmax 则在 z 给定时 xy 表示为平面 这里 给定 5 3km 其中的 u0 和 G 为常数注意单位 m n 分别为横测线数 纵测线数 都为 27 程序 program zhongci4 PARAMETER PI 3 1415926 m 27 n 27 character 80 filename1 filename2 filename3 filename4 filename5 filename6 filename7 filename8 filename9 file name10 real s 3 11 td 1 m 1 n t1 1 m 1 n t2 1 m 1 n t3 1 m 1 n gd 1 m 1 n g1 1 m 1 n g2 1 m 1 n g3 1 m 1 n real xmin xmax ymin ymax zmin zmax eigval is ix iy u0 G 控制网格里数的上限值不能超过 eigval eigval 10 9 filename1 t1 dat filename2 t2 dat filename3 t3 dat filename4 t 叠加 dat filename5 g1 dat filename6 g2 dat filename7 g3 dat filename8 g 叠加 dat filename9 t 叠加 grd filename10 g 叠加 grd call cdr is ix iy xmin xmax ymin ymax zmin zmax s call jisaunciyichang is ix m n t1 XMAX YMAX YMIN XMIN zmin zmax u0 filename1 s 1 call jisaunciyichang is ix m n t2 XMAX YMAX YMIN XMIN zmin zmax u0 filename2 s 2 call jisaunciyichang is ix m n t3 XMAX YMAX YMIN XMIN zmin zmax u0 filename3 s 3 call changyuandiejia t1 t2 t3 td m n filename4 call jszlyc m n g1 XMAX YMAX YMIN XMIN zmin zmax filename5 s 1 call jszlyc m n g2 XMAX YMAX YMIN XMIN zmin zmax filename6 s 2 call jszlyc m n g3 XMAX YMAX YMIN XMIN zmin zmax filename7 s 3 call changyuandiejia g1 g2 g3 gd m n filename8 调用场源叠加外部子例行程序 call OUTPUT GRD td filename9 1 27 1 27 27 27 eigval xmin xmax ymin ymax call OUTPUT GRD gd filename10 1 27 1 27 27 27 eigval xmin xmax ymin ymax end 所有的参数读入 subroutine cdr is ix iy xmin xmax ymin ymax zmin zmax s real is ix iy xmin xmax ymin ymax zmin zmax s 1 3 1 11 读入地磁场的参数 倾角与 x 轴的夹角 与 y 轴的夹角 open 12 file 地磁场参数 txt read 12 is ix iy close 12 读入曲面参数 当 zmin zmax 时则读入的为一平面 open 11 file 平面 xyz txt read 11 xmin xmax ymin ymax ZMIN ZMAX close 11 输入三个场源的参数 用数组 s 的一行表示一个场源的各个参数 open 12 file source txt do i 1 3 read 12 s i j j 1 11 enddo close 12 end subroutine GRD 格式数据输出子程序 将数组中的数据读入到网格中 SUBROUTINE OUTPUT GRD A filename m1 m2 n1 n2 m n eigval xmin xmax ymin ymax integer m1 n1 m2 n2 m n real A 1 m 1 n character 80 filename real amin amax eigval xmin xmax ymin ymax amin huge amin amax huge amax DO j n1 n2 1 do i m1 m2 1 if ABS A i j lt eigval then amin MIN amin A I j amax MAX amax A I j end if end do END DO OPEN 20 file filename form formatted write 20 A DSAA write 20 m2 m1 1 n2 n1 1 write 20 xmin xmax write 20 ymin ymax write 20 amin amax Do j n1 n2 1 write 20 A i j i m1 m2 End do Close 20 END subroutine OUTPUT GRD 计算磁异常 subroutine jisaunciyichang is ix m n t XMAX YMAX YMIN XMIN zmin zmax u0 filename s m3 m1 为 磁化强度 m3 从 1 变到三 parameter pi 3 141592 real M0 L0 N0 A B R K1 K2 K3 K4 K5 K6 character 80 filename integer m n m3 xz yz REAL i1 i2 i3 real is ix ix1 is1 real t 1 m 1 n x 1 2 y 1 2 z 1 2 s 1 3 1 11 real m1 u0 XMAX YMAX YMIN XMIN zmin zmax u0 1 查数据 M1 s m3 2 m1 为磁化强度 is1 ix1 iy1 分别为磁化强度的倾角 与 x 轴的夹角 与 y 轴的夹角 is1 s m3 3 ix1 s m3 4 x 1 s m3 6 x 2 s m3 7 y 1 s m3 8 y 2 s m3 9 z 1 s m3 10 z 2 s m3 11 z 应该可以控制 绘制曲面的时候 z 可以改变 i3 zmin open 30 file filename 得修改 K1 COSd IS SINd IX SINd IS1 SINd IS COSd IS1 SINd IX1 MO R N0 B K2 COSd IS COSd IX SINd IS1 SINd IS COSd IS1 COSd IX1 L0 R N0 A K3 COSd IS COSd IX COSd IS1 SINd IX1 COSd IS SINd IX COSd IS1 COSd IX1 LO B MO A K4 COSd IS COSd IX COSd IS1 COSd IX1 L0 A K5 COSd IS SINd IX COSd IS1 SINd IX1 M0 B K6 SINd IS Sind IS1 N0 R do i1 XMIN XMAX 2 用到坐标和数组序号之间的换算关系 xz 26 i1 2 1 do i2 ymin ymax 2 用到坐标和数组序号之间的换算关系 yz 26 i2 2 1 c 0 do i 1 2 do j 1 2 do k 1 2 r0 sqrt i1 x i 2 i2 y j 2 i3 z k 2 a1 x i i1 y j i2 b1 x i i1 2 r0 z k i3 z k i3 2 d1 y j i2 2 r0 z k i3 z k i3 2 e1 r0 z k i3 cc 1 i k j te cc U0 4 PI M1 k1 alog r0 x i i1 k2 alog r0 y j i2 k3 alog r0 z k i3 k4 atan a1 b1 k5 atan a1 d1 k6 atan a1 e1 c c te enddo enddo enddo t xz yz c 0 1 换算单位 使输出的单位为 nt write 30 i1 i2 t xz yz enddo enddo close 30 end subroutine 叠加 用到坐标和数组序号之间的换算关系 subroutine changyuandiejia t1 t2 t3 td m n filename character filename integer m n xz yz real td 1 m 1