重磁,最新边缘识别实验报告

位场边缘识别与程序设计位场边缘识别与程序设计 实验实验 学生姓名 专业名称 学生学号 指导老师 一 基本原理 2 1 1 位场边缘识别方法研究进展概述 2 1 2 数值计算类边缘识别方法的研究及计算 2 1 2 1 垂向导数 VDR 2 1 2 2 总水平导数 THDR 2 1 2 3 解析信号振幅 ASM 2 1 2 4 倾斜角 TA 3 1 2 5 二阶导数类边缘识别 3 二 输入 输出数据格式设计 3 2 1 主要变量设计 3 2 2 位场边缘识别程序数据格式设计 4 三 总体设计 5 四 测试结果 5 4 1 测试用例 5 4 2 测试结果 6 4 3 结果分析 8 五 结论与建议 9 5 1 结论 9 5 2 建议 9 附录 边缘识别程序源代码 10 一 基本原理 1 1 位场边缘识别方法研究进展概述位场边缘识别方法研究进展概述 地质目标体边缘是指断裂构造线 不同地质体的边界线 实际上是具有一 定密度或磁性差异的地质体的边界线 在地质体的边缘附近 重 磁异常变化 率较大 故所有的边缘识别方法均利用这一特点进行设计 现有重 磁位场边 缘识别方法分为数理统计 数值计算和其他三大类 1 2 数值计算类边缘识别方法的研究及计算数值计算类边缘识别方法的研究及计算 数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置 其依据的理论基础是二度体铅垂台阶模型的重力异常特征 在该模型边缘处重 力异常总水平导数和解析信号振幅达到极大值 垂向导数达到零值 故可以利 用这些特征位置来确定二度体铅垂台阶的边缘位置 确定倾斜二度体 不规则 二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广 但确定 的边缘位置与真实位置有一定偏差 该偏差随着地质体边界形状 埋深 水平 尺寸及物性差异等的变化而变化 因此 边缘识别结果是一种定性或半定量解 释结果 与定量解释结果有一定区别 识别结果可作为边缘位置定量反演的初 值 1 2 1 垂向导数 VDR 垂向导数方法利用零值位置确定地质体的边缘位置 重力异常可直接使用 对磁力异常必须转换成磁源重力异常 假重力异常 或化极磁力异常才可以使 用 垂向导数方法研究历史较早 方法较成熟 应用频率较高 通过我们的试 验和研究认为 随着地质体埋深的增大 垂向导数所确定的边缘位置偏离实际 位置的距离越大 1 2 2 总水平导数 THDR 总水平导数是利用其极大值位置来确定地质体的边缘位置 适用于重力异 常 对磁力异常必须换算成磁源重力异常或化极磁力异常才可以使用 1 2 3 解析信号振幅 ASM 解析信号振幅也是利用极大值位置来确定地质体的边缘位置 适用于重力 异常和磁力异常 1 2 4 倾斜角 TA 倾斜角实质上是垂向导数和总水平导数的比值 由于倾斜角为一阶导数的 比值 所以能很好地平衡高幅值异常和低幅值异常 起到边缘增强的效果 1 2 5 二阶导数类边缘识别 为了增强边缘分辨能力 和克服当总水平导数等于 时倾斜角存在 解析 奇点 会使得计算结果不稳定 2 输入 输出数据格式设计 依据上述原理 现对上述各种边缘识别方法实现进行程序设计 2 1 主要变量设计主要变量设计 cmd file 参数文件名 input field filename 输入位场文件名 output field filename 输出转换后位场文件名 expan 2D method 2D 扩边方法选择 Field Deriv method 边缘识别方法选择 num area 计算区域大小选择 mpoint nline 点数 线数 m0 m1 m2 m3 扩边后 X 方向上各端点 n0 n1 n2 n3 扩边后 Y 方向上各端点 Xmin Xmax X 坐标最小 最大值 ymin Ymax Y 坐标最小 最大值 Ori Field 原始场值 二维数组 m3 n3 Deriv Field 转换后后场值 二维数组 m3 n3 Deriv operator x x 方向转换因子 二维数组 m3 n3 Deriv operator y y 方向转换因子 二维数组 m3 n3 Deriv operator z z 方向转换因子 二维数组 m3 n3 2 2 位场边缘识别程序数据格式设计位场边缘识别程序数据格式设计 cmd 文件 输入参数文件 格式如下 cmd txt input field filename anomaly grd output field filename VDR THDR anomaly grd num area 0 expan 2D method 1 Field Deriv method 7 说明 num area 取 0 表示一般扩边区域 取 1 表示计算区域相对于一般扩边再放大一倍 expan 2D method 取 1 表示余弦扩边 端点值取边界平均值 取 2 表示余弦扩边 端点值取全场区平均值 取 3 表示余弦扩边 端点值取常数 0 取 4 表示最小曲率扩边 空白处初值取余弦扩边值 端点初 值取边界平均值 取 5 表示最小曲率扩边 空白处初值取余弦扩边值 端点初 值取全场区平均值 取 6 表示最小曲率扩边 空白处初值取余弦扩边值 端点初 值取常数 0 Field Deriv method 取 1 计算垂向导数 VDR 取 2 计算总水平导数 THDR 取 3 计算解析信号振幅 ASM 取 4 计算倾斜角 TA 取 5 计算垂向二阶导数 VDR VDR 取 6 计算倾斜角总水平导数 TA THDR 取 7 计算垂向导数总水平导数 VDR THDR 原始位场文件 输入文件 grd 格式 ASCII 延拓后位场文件 输出文件 grd 格式 ASCII Grd 格式如下 DSAA mpoint nline Xmin Xmax Ymin Ymax Zmin Zmax Z11 Z12 3 总体设计 四 测试结果 4 1 测试用例测试用例 1 观测面上的重力异常存放在 anomaly grd 中 坐标单位为 m 2 形体边界存放在 rectangle bln 中 坐标单位 图 3 1 重 磁异常边缘识别程序设计流程图 4 2 测试结果测试结果 图 4 2 2 THDR 总水平导数 结果图 1000 800 600 400 20002004006008001000 1000 800 600 400 200 0 200 400 600 800 1000 0 0 002 0 004 0 006 0 008 0 01 0 012 0 014 0 016 0 018 0 02 0 022 0 024 0 026 0 028 0 03 0 032 0 034 0 036 0 038 0 04 0 042 0 044 0 046 0 048 0 05 0 052 1000 800 600 400 20002004006008001000 1000 800 600 400 200 0 200 400 600 800 1000 0 006 0 01 0 014 0 018 0 022 0 026 0 03 0 034 0 038 0 042 0 046 0 05 0 054 0 058 0 062 0 066 图 4 2 3 ASM 解析信号振幅 结果图 图 4 2 6 VDR THDR 结果图 垂向导数总水平导数 余弦扩边方法 1000 800 600 400 20002004006008001000 1000 800 600 400 200 0 200 400 600 800 1000 0 4E 005 8E 005 0 00012 0 00016 0 0002 0 00024 0 00028 0 00032 0 00036 0 0004 0 00044 4 3 结果分析结果分析 1 从图 4 2 1 图 4 2 7 可以看出 VDR 垂向导数 VDR VDR 垂向二阶 导 AT 倾斜角 实验结果零值线与实际模型 白线框 基本重合 与实验原理 一致 说明这几种方法程序实现是正确的 2 同时从图 4 2 1 图 4 2 7 可以也可看出 THDR 总水平导数 VDR THDR 垂向导数总水平导数 AT THDR 倾斜角总水平导数 ASM 解析信号振幅 实验结果极大值线与实际模型 白线框 基本重合 与 实验原理一致 说明这几种方法程序实现也是正确的 3 该实验模型 从以上各方法边缘识别效果来看 ASM 解析信号振幅 分 辨能力较低 二阶导数类方法识别效果比一阶导数类识别效果好 就一阶导数 类来看 AT 倾斜角 由于其他三种方法 二阶导数类方法中 AT THDR 倾斜角总 水平导数 识别最精确 但是从 图 4 2 7 看出 其稳定性较差 故相比较而言 VDR THDR 垂向导数总水平导数 VDR VDR 垂向二阶导 效果是比较理想 的 4 对比图 4 2 6 和图 4 2 8 可以看出在频率域处理选用扩边方法不同对计 算结果有一定影响 其中最小曲率扩边计算效果为佳 五 结论与建议 5 1 结论结论 此方法在一定程度上是有效的 且从实验结果可以垂向导数总水平导数和 垂向二阶导识别的效果较好 5 2 建议建议 没有学好位场边缘识别方法没有可以能提出的建议 附录 边缘识别程序源代码 PROGRAM Fre PoteField Deriv PURPOSE Frequency Potential field Derivative Author gesang Data 2013 11 PROGRAM Fre PoteField Deriv implicit none CHARACTER 80 input field filename CHARACTER 80 output field filename INTEGER expan 2D method num area Field Deriv method INTEGER mpoint nline INTEGER m0 m1 m2 m3 INTEGER n0 n1 n2 n3 REAL Xmin Xmax REAL ymin Ymax REAL dx dy REAL ALLOCATABLE Ori Field REAL ALLOCATABLE Field Real REAL ALLOCATABLE Field Image REAL ALLOCATABLE Deriv Field REAL ALLOCATABLE Deriv operator x REAL ALLOCATABLE Deriv operator y REAL ALLOCATABLE Deriv operator z z方向导数因子 INPUT 读取文件参数 CALL read cmdinfo input field filename output field filename expan 2D method Field Deriv method num area 读取点数线数及测区范围 CALL get mpoint and nline input field filename mpoint nline Xmin Xmax Ymin Ymax 计算扩边端点 CALL get expan num mpoint m0 m1 m2 m3 num area CALL get expan num nline n0 n1 n2 n3 num area ALLOCATE Ori Field m0 m3 n0 n3 ALLOCATE Field Real m0 m3 n0 n3 ALLOCATE Field Image m0 m3 n0 n3 ALLOCATE Deriv Field m0 m3 n0 n3 ALLOCATE Deriv operator x m0 m3 n0 n3 ALLOCATE Deriv operator y m0 m3 n0 n3 ALLOCATE Deriv operator z m0 m3 n0 n3 从网格化文件读入初始场值 CALL Read field input field filename m0 m1 m2 m3 n0 n1 n2 n3 Ori Field PROCESS 扩边 dx Xmax Xmin mpoint 1 dy Ymax ymin nline 1 CALL expan 2D m0 m1 m2 m3 n0 n1 n2 n3 Ori Field dx dy expan 2D method expan 2D method 扩边方法选择 计算导数因子 CALL Deriv operator sub dx dy m3 n3 Deriv operator x Deriv operator y Deriv operator z 求导运算 Field Real Ori Field Field Image 0 0 CALL field Deriv Field Real Field Image Deriv operator x Deriv operator y Deriv operator z m0 m1 m2 m3 n0 n1 n2 n3 Field Deriv method Deriv Field Field Real OUTPUT 输出求导后位场文件 CALL output grd output field filename m0 m1 m2 m3 n0 n1 n2 n3 mpoint nline Deriv Field Xmin Xmax Ymin Ymax DEALLOCATE Ori Field DEALLOCATE Field Real DEALLOCATE Field Image DEALLOCATE Deriv Field DEALLOCATE Deriv operator x DEALLOCATE Deriv operator y DEALLOCATE Deriv operator z END PROGRAM Fre PoteField Deriv 开始 得到输入参数子程序开始 SUBROUTINE read cmdinfo input field filename output field filename expan 2D method Field Deriv method num area CHARACTER input field filename CHARACTER output field filename INTEGER expan 2D method Field Deriv method num area INTEGER unit CHARACTER 80 s CALL get unit unit OPEN unit file cmd txt status old READ unit s input field filename READ unit s output field filename READ unit s num area READ unit s expan 2D method READ unit s Field Deriv method CLOSE unit END SUBROUTINE read cmdinfo 得到输入参数子程序结束 读入grd格式文件中点数和线数及范围开始 SUBROUTINE get mpoint and nline filename mpoint nline Xmin Xmax Ymin Ymax CHARACTER filename INTEGER mpoint nline REAL Xmin Xmax Ymin Ymax INTEGER unit CALL get unit unit open unit file filename status old err 999 READ unit READ unit mpoint nline READ unit Xmin Xmax READ unit Ymin Ymax close unit RETURN PRINT PAUSE STOP ENDSUBROUTINE get mpoint and nline 读入grd格式文件中点数和线数及范围结束 得到扩边端点子程序开始 SUBROUTINE get expan num mpoint m0 m1 m2 m3 flag IMPLICIT NONE INTEGER mtemp factor m mu INTEGER mpoint m0 m1 m2 m3 INTEGER flag mtemp mpoint factor m 1 DO WHILE mod mtemp 2 eq 0 and mtemp ne 0 mtemp mtemp 2 End do IF mtemp eq 1 THEN m3 mpoint 2 ELSE mu int log float mpoint 0 693147 factor m m3 2 mu if flag 1 then m3 m3 2 计算区域为原来倍 else endif END IF m1 1 m3 mpoint 2 m2 m1 mpoint 1 m0 1 END SUBROUTINE get expan num 得到扩边端点子程序结束 读入grd格式文件中场值子程序开始 SUBROUTINE Read field filename m0 m1 m2 m3 n0 n1 n2 n3 G PARAMETER vial 1 701411e 38 CHARACTER filename INTEGER m0 m1 m2 m3 n0 n1 n2 n3 REAL G m0 m3 n0 n3 INTEGER unit CALL get unit unit G vial open unit file filename DO i 1 5 READ unit ENDDO READ unit G i j i m1 m2 j n1 n2 CLOSE unit ENDSUBROUTINE Read field 读入grd格式文件中场值子程序结束 扩边子程序集开始 扩边子程序调用主子程序开始 SUBROUTINE expan 2D m0 m1 m2 m3 n0 n1 n2 n3 Ori Field dx dy expan 2D method INTEGER m0 m1 m2 m3 INTEGER n0 n1 n2 n3 REAL Ori Field m0 m3 n0 n3 INTEGER expan 2D method REAL dx dy IF expan 2D method 0 9 eigval bpn bpn 1 ENDDO ENDDO ENDSUBROUTINE 获得数据中的空白点数结束 获得数据中的空白点位置开始 SUBROUTINE Blank Point position FIELD mpoint line bpn eigval P POINT P LINE IMPLICIT NONE INTEGER mpoint line bpn i j k INTEGER P POINT 1 bpn P LINE 1 bpn REAL FIELD 1 mpoint 1 line eigval k 1 DO i 1 mpoint 1 DO j 1 line 1 IF FIELD i j 0 9 eigval AND jeps abs and iteration iteration max eps error 0 CALL MinCurV 2D boundary U mpoint line alph DO k 1 number eigval 1 i M eigval k j N eigval k temp U i 2 j U i 2 j alph1 U i j 2 U i j 2 alph2 U i 1 j 1 U i 1 j 1 U i 1 j 1 U i 1 j 1 alph3 U i 1 j U i 1 j alph4 U i j 1 U i j 1 temp temp alph0 eps error MAX ABS temp U i j eps error U i j temp ENDDO iteration iteration 1 ENDDO print iteration eps error DO j 1 line 1 DO i 1 mpoint 1 FIELD i j U i j ENDDO ENDDO deALLOCATE U STAT ierr ENDSUBROUTINE 2D网格数据无约束最小曲率迭代结束 2D余弦扩边开始 SUBROUTINE expan cos2 m1 m2 m3 m4 n1 n2 n3 n4 G flag PARAMETER pi 3 141592654 INTEGER flag INTEGER m1 m2 m3 m4 INTEGER n1 n2 n3 n4 REAL G m1 m4 n1 n4 REAL G1 m1 m4 n1 n4 REAL G2 m1 m4 n1 n4 REAL sum INTEGER num sum 0 0 num 0 IF flag 1 then 扩边端点值取边界平均值 DO i m2 m3 sum sum G i n2 G i n3 num num 2 END DO DO i n2 1 n3 1 sum sum G m2 i G m3 i num num 2 END DO DO i m1 m4 G i n1 sum num G i n4 sum num END DO DO i n1 1 n4 1 G m1 i sum num G m4 i sum num END DO ELSE IF flag 2 then 扩边端点值取全部数据平均值 DO i m2 m3 DO j n2 n3 sum sum G i j num num 1 END DO END DO DO i m1 m4 G i n1 sum num G i n4 sum num END DO DO i n1 1 n4 1 G m1 i sum num G m4 i sum num END DO ELSE IF flag 3 then 扩边端点值取常数 DO i m1 m4 G i n1 0 G i n4 0 END DO DO i n1 1 n4 1 G m1 i 0 G m4 i 0 END DO END IF 扩边 G1 G G2 G DO j n2 n3 DO i m1 1 m2 1 G1 i j G1 m2 j G1 m1 j cos m2 i 1 0 m2 m1 pi 2 0 G1 m1 j END DO DO i m3 1 m4 1 G1 i j G1 m4 j G1 m3 j cos m4 i 1 0 m4 m3 pi 2 0 G1 m3 j END DO END DO DO i m1 m4 DO j n1 1 n2 1 G1 i j G1 i n2 G1 i n1 cos n2 j 1 0 n2 n1 pi 2 0 G1 i n1 END DO DO j n3 1 n4 1 G1 i j G1 i n4 G1 i n3 cos n4 j 1 0 n4 n3 pi 2 0 G1 i n3 END DO END DO DO i m2 m3 DO j n1 1 n2 1 G2 i j G2 i n2 G2 i n1 cos n2 j 1 0 n2 n1 pi 2 0 G2 i n1 END DO DO j n3 1 n4 1 G2 i j G2 i n4 G2 i n3 cos n4 j 1 0 n4 n3 pi 2 0 G2 i n3 END DO END DO DO j n1 n4 DO i m1 1 m2 1 G2 i j G2 m2 j G2 m1 j cos m2 i 1 0 m2 m1 pi 2 0 G2 m1 j END DO DO i m3 1 m4 1 G2 i j G2 m4 j G2 m3 j cos m4 i 1 0 m4 m3 pi 2 0 G2 m3 j END DO END DO G G1 G2 2 0 END SUBROUTINE expan cos2 2D余弦扩边开始 扩边子程序集结束 计算导数因子开始 SUBROUTINE Deriv operator sub dx dy m n Deriv operator x Deriv operator y Deriv operator z IMPLICIT NONE REAL dx dy REAL dz INTEGER m n i j REAL W 1 m 1 n U 1 m V 1 n REAL Deriv operator Z 1 m 1 n REAL Deriv operator x 1 m 1 n REAL Deriv operator y 1 m 1 n CALL WAVE2D sub U V W m n dx dy 计算垂向导数因子 实数 Deriv operator z W 计算X和y导数因子 虚数 do j 1 n do i 1 m Deriv operator x i j U i Deriv operator y i j V j end do end do END SUBROUTINE Deriv operator sub SUBROUTINE WAVE2D sub U V W m n dx dy REAL W 1 m 1 n U 1 m V 1 n INTEGER m n REAL dx dy REAL i j CALL Wave1D sub n dy V CALL Wave1D sub m dx U DO j 1 n 计算园频率 DO i 1 m W i j SQRT U i U i V j V j END DO END DO END SUBROUTINE WAVE2D sub Subroutine Wave1D sub N dy V REAL dy INTEGER N REAL V 0 N 1 REAL deltf deltf 2 0 3 141592654 N dy Do j 0 N 2 1 V j j deltf End do Do j N 2 1 N 1 1 V j j n deltf End do END subroutine Wave1D sub 计算导数因子子程序结束 求导子程序集开始 求导子程序总模块开始 SUBROUTINE field Deriv Field Real Field Image Deriv operator x Deriv operator y Deriv operator z m0 m1 m2 m3 n0 n1 n2 n3 Field Deriv method IMPLICIT NONE INTEGER m0 m1 m2 m3 n0 n1 n2 n3 REAL Field Real m0 m3 n0 n3 Field Image m0 m3 n0 n3 REAL Deriv operator x m0 m3 n0 n3 Deriv operator y m0 m3 n0 n3 Deriv operator z m0 m3 n 0 n3 INTEGER Field Deriv method IF Field Deriv method 1 THEN 计算垂向导数 CALL VDR Field Real Field Image m3 n3 Deriv operator z ELSE IF Field Deriv method 2 THEN 计算总水平导数 call THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y ELSE IF Field Deriv method 3 THEN 计算解析信号振幅 call ASM meth Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operat or z ELSE IF Field Deriv method 4 THEN 计算倾斜角 CALL TA Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operator z ELSE IF Field Deriv method 5 THEN 计算垂向二阶导数 call VDR VDR Field Real Field Image m3 n3 Deriv operator z ELSE IF Field Deriv method 6 THEN 计算倾斜角总水平导数 CALL TA THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operat or z ELSE 计算垂向导数总水平导数 call VDR THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv oper ator z END IF END SUBROUTINE field Deriv 延拓子程序总模块结束 VDR SUBROUTINE VDR Field Real Field Image m3 n3 Deriv operator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 CALL FFT2 Field Real Field Image m3 n3 2 CALL multiply sub Field Real Field Image m3 n3 Deriv operator z CALL FFT2 Field Real Field Image m3 n3 1 END SUBROUTINE VDR VDR 相乘运算开始 SUBROUTINE multiply sub Field Real Field Image m3 n3 Deriv operator IMPLICIT NONE INTEGER m3 n3 i j REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator 1 m3 1 n3 DO i 1 m3 DO j 1 n3 Field Real i j Field Real i j Deriv operator i j Field Image i j Field Image i j Deriv operator i j END DO END DO END SUBROUTINE multiply sub 相乘运算结束 THDR SUBROUTINE THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y INTEGER m3 n3 I J REAL Deriv operator X 1 m3 1 n3 Deriv operator y 1 m3 1 n3 REAL Field Real 1 m3 1 n3 REAL Field Image 1 m3 1 n3 REAL Field Real x 1 m3 1 n3 Field Image x 1 m3 1 n3 REAL Field Real y 1 m3 1 n3 Field Image y 1 m3 1 n3 CALL FFT2 Field Real Field Image m3 n3 2 do i 1 m3 do j 1 n3 Field Real x i j Deriv operator X i j Field Image i j Field Image x i j Deriv operator X i j Field Real i j Field Real y i j Deriv operator y i j Field Image i j Field Image y i j Deriv operator y i j Field Real i j end do end do CALL FFT2 Field Real x Field Image x m3 n3 1 CALL FFT2 Field Real y Field Image y m3 n3 1 do i 1 m3 do j 1 n3 Field Real i j sqrt Field Real x i j 2 Field Real y i j 2 end do end do END SUBROUTINE THDR THDR ASM meth SUBROUTINE ASM meth Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Field Real THDR 1 m3 1 n3 Field Image THDR 1 m3 1 n3 REAL Field Real VDR 1 m3 1 n3 Field Image VDR 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 REAL Deriv operator X 1 m3 1 n3 Deriv operator y 1 m3 1 n3 Field Real THDR Field Real Field Image THDR Field Image CALL THDR Field Real THDR Field Image THDR m3 n3 Deriv operator x Deriv operator y Field Real VDR Field Real Field Image VDR Field Image CALL VDR Field Real VDR Field Image VDR m3 n3 Deriv operator z DO i 1 m3 do j 1 n3 Field Real i j sqrt Field Real THDR i j 2 Field Real VDR i j 2 end do end do END SUBROUTINE ASM meth ASM meth TA SUBROUTINE TA Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 REAL Deriv operator X 1 m3 1 n3 Deriv operator y 1 m3 1 n3 REAL Field Real THDR 1 m3 1 n3 Field Image THDR 1 m3 1 n3 REAL Field Real VDR 1 m3 1 n3 Field Image VDR 1 m3 1 n3 Field Real THDR Field Real Field Image THDR Field Image CALL THDR Field Real THDR Field Image THDR m3 n3 Deriv operator x Deriv operator y Field Real VDR Field Real Field Image VDR Field Image CALL VDR Field Real VDR Field Image VDR m3 n3 Deriv operator z DO i 1 m3 do j 1 n3 Field Real i j atand Field Real VDR i j Field Real THDR i j end do end do END SUBROUTINE TA TA VDR VDR SUBROUTINE VDR VDR Field Real Field Image m3 n3 Deriv operator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 CALL FFT2 Field Real Field Image m3 n3 2 CALL multiply sub Field Real Field Image m3 n3 Deriv operator z CALL multiply sub Field Real Field Image m3 n3 Deriv operator z CALL FFT2 Field Real Field Image m3 n3 1 END SUBROUTINE VDR VDR VDR VDR TA THDR SUBROUTINE TA THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv operator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 REAL Deriv operator X 1 m3 1 n3 Deriv operator y 1 m3 1 n3 REAL Field Real THDR 1 m3 1 n3 Field Image THDR 1 m3 1 n3 REAL Field Real TA 1 m3 1 n3 Field Image TA 1 m3 1 n3 Field Real TA Field Real Field Image TA Field Image CALL TA Field Real TA Field Image TA m3 n3 Deriv operator x Deriv operator y Deriv oper ator z Field Real THDR Field Real TA Field Image THDR Field Image CALL THDR Field Real THDR Field Image THDR m3 n3 Deriv operator x Deriv operator y Field Real Field Real THDR END SUBROUTINE TA THDR TA THDR VDR THDR SUBROUTINE VDR THDR Field Real Field Image m3 n3 Deriv operator x Deriv operator y Deriv oper ator z INTEGER m3 n3 REAL Field Real 1 m3 1 n3 Field Image 1 m3 1 n3 REAL Deriv operator z 1 m3 1 n3 REAL Deriv operator X 1 m3 1 n3 Deriv operator y 1 m3 1 n3 REAL Field Real THDR 1 m3 1 n3 Field Image THDR 1 m3 1 n3 REAL Field Real VDR 1 m3 1 n3 Field Image VDR 1 m3 1 n3 Field Real VDR Field Real Field Image VDR Field Image CALL VDR Field Real VDR Field Image VDR m3 n3 Deriv operator z Field Real THDR Field Real VDR Field Image THDR Field Image CALL THDR Field Real THDR Field Image THDR m3 n3 Deriv operator x Deriv operator y Field Real Field Real THDR END SUBROUTINE VDR THDR VDR THDR 离散Fourier变换子程序集开始 c c c c c 2 D FAST FOURIER TRANSFORM FOR COMPLEX FUNCTION c c c c DREAL The real part of the function to be transformed c DIMAGE The imaginary part of the function to be transformed c M The number of points M must be the nu th power of 2 c N The number of points N must be the nu th power of 2 c NF 1 reverse transform c 2 normal transform c c Author gesang c TREAL Real part of the working array c TIMAG Imaginary part of the working array c c SUBROUTINE FFT2 DREAL DIMAG M N NF real dreal m n dimag m n real ALLOCATABLE treal timag nmmax max m n allocate treal nmmax timag nmmax STAT ierr DO i 1 m IF n ne 1 THEN do j 1 n treal j dreal i j timag j dimag i j end do call fft treal timag n nf do j 1 n dreal i j treal j dimag i j timag j end do END IF END DO DO j 1 n IF m ne 1 THEN do i 1 m treal i dreal i j timag i dimag i j end do call fft treal timag m nf