5 10 15 20 25 30 35 40 中国科技论文在线 探地雷达 FDTD 正演系统开发及图像处理 http://www.paper.edu.cn 徐荣华，崔伟雄* （中国矿业大学资源学院，江苏 徐州 221008） 摘要：本文以麦克斯韦方程组为基本原理，推导以 PML 为吸收边界的高阶有限差分（FDTD） 迭代公式，编制一套操作简便、精确的处理正演模拟系统，模拟得到二维地电剖面，准确反 映地下异常体位置，并能用作后序处理极大的提高了地质雷达研究的灵活性。同时，编制处 理程序对剖面增益（AGC）、二维滤波（FK）提取散射波特征，提高正演分辨率。 关键词：麦克斯韦方程组；正演模拟系统；增益；二维滤波；散射波 中图分类号：P631.4 The Finite-difference Time-domain(FDTD)Forward Simulation and image process for Ground Penetrating Radar(GPR) XU Ronghua, CU Weixiong (China Uinversity of Mining Technology, JiangSu XuZhou 221008) Abstract: In this paper, from the two vortices equation of Maxwell equation group,the author utilize the equation of PML as absorbed boundary to develop a convenient and accurate forward Simulation system. The section simulated by the soft can accurately response the position of the unnormal substace,and is also could be used for processing during the period of process. It highly improved the agility of the GPR research. In addition to above,the author make the process to improve the resolution of the image, such as automatic gain control(AGC)、two-dimensional filtering(F-K). The process program can separate the scattered wave from the section which has highly improved the resolution of exploration. Keywords:The two vortices equation of Maxwell equation; Forward Simulation system; AGC; Two-dimensional filtering; Scattered wave 0 引言 脉冲探地雷达(GPR)技术在浅层探测中的应用中，雷达天线接收到的信号时由直达波信 号，目标散射信号以及噪声信号几部分合成[1]。由于地下介质的不均匀性以及对波的不同程 度的吸收，使得接收天线获得的绕射脉冲振幅小，波形发生变异，由此得到的雷达剖面质量 较低，这样就影响了雷达剖面的准确解释，甚至导致错误的解释结果。因此，对接收信号进 行适当的处理是有必要的，从而改善雷达剖面质量，为后期的解释工作提供清晰可辨的图像。 同时由于直达波是直接耦合和经过地表传到接收仪器，相对于目标散射信号具有更大幅度， 如果直达波与目标散射信号在时间上能够分离，则可借时间划分来判定目标，在所关心目标 深度不是很深时，目标散射信号与直达波相互重叠，这样直达波对信号形成严重干扰，所以 抑制直达波对绕射波的提取是获得探地雷达有效信号从而提高探地雷达分辨率的有效手段。 本文利用自己编制雷达正演系统得到剖面，该软件采用 PML[2]方程为吸收边界，并采 用能量补偿和道内均衡方法对波场模拟得到的雷达剖面经行处理，使散射波能量明显增强， 并采用 FK 变换对剖面进行二维滤波，切除直达波和反射层能量提取散射波特征。 作者简介：徐荣华，（1986-），男，硕士在读，地震数据处理. E-mail: xuronghua86@163.com - 1 - 

中国科技论文在线 1 时域有限差分正演模拟系统 1.1 理论推导 FDTD 公式 本系统以麦克斯韦方程组为理论推导 FDTD 迭代公式即 http://www.paper.edu.cn ∇× = − E B ∂ t ∂ D ∂ t ∂ ∇× H j = + (1.1.1) B 0 D ρ ∇ • = ∇ • = 式中，ρ为电荷密度（C/m3）；J 为电流密度（A/m2）；D 为电位移(C/m2)；H 为磁场 强度（A/m）；B 为磁感应强度(T)；E 为电场强度（V/m）。式 1 为微分形式的法拉第电磁 感应定律；式 2 称为安培电流环路定律。采用有限差分法，首先要把变量空间进行离散化， 也就是建立合适的网格剖发体系。本文采用的差分近似格式[3]如下表 1. 微分式 du dx i 表 1 不同微分格式的差分近似 Tab.1 different difference format 有限差分近似 u u −+ i 1 1 + h 2 i d u 2 dx 2 i u i 1 + − + u i 1 − u 2 i h 2 示意图 45 50 55 在二维 TM 波情况下，根据表 1 差分格式采用 Yee 氏[1]网格对麦克思维方程组离散可得 n 用于地质雷达正演的 FDTD 格式： j E i j CD CB i j ( , ) ( , ) × + × j H i i j , ) ( ( , − − 1 2 CA i ( , ) j H , ) − i j ( , ) i + ( + E z H [ = 1 2 z 1 + 2 + 1 + + n n n n 1 2 1 2 y y x n + H x 1 2 n + H y i ( , 1 2 i ( j + + 1 2 ) 1 2 = H x n − 1 2 j ), = H y n i ( , 1 2 − j + i ( + ) 1 2 1 2 + CD E i × [ n z 1 2 j ( , ) − ) − H n + 1 2 x E i ( , z n (1.1.2) i ( , j + )] 1 2 j + 1)] (1.1.3) j ), + CD × [ n iE ( z + j ),1 − iE ,( z j )] (1.1.4) 上式中： CA i j ( , ) = 1 − 1 + i j t ( , ) σ Δ i j 2 ( , ) ε t i j ( , ) σ Δ i j 2 ( , ) ε ; CD = 1 t Δ s εμ Δ 0 0 CB i j ( , ) ; = i j ( , ) ε ε 0 σ + i j ( , ) 2 , tΔ 、 sΔ 分 Δ t 60 别为波长模拟时间步长和空间步长。 1.2 数值频散压制和稳定性条件 数值频散是由于用有限差分离散代替电磁波传播方程连续微商造成的部分，它往往导致 脉冲波形的畸变和虚假折射[4]，由理论推导可得二维空间中 TM 波差分方程的数值频散关系 - 2 - 

中国科技论文在线 http://www.paper.edu.cn 如式（2.2.1），当时间和空间步长无限减小时[4][5][6]即 tΔ 时数值频散关系满足解析 方法得到的二维 TM 波频散关系如公式（2.2.2）。在二维 Yee 氏网格中，为了使波场模拟 稳定，必须满足式（2.2.3）。因此只要选择合适的网格参数，可以实现理想的频散关系和差 分格式稳定。 0→Δs ) sin ( 2 2 ω Δ 2 t t 2 ) = 1 x Δ 2 sin ( 2 K x Δ x 2 2 ) + 1 y Δ 2 sin ( 2 K y Δ y 2 ) (1.2.1) 1 ( υ Δ 2 ω 2 υ = 2 K x + 2 K y (1.2.2) t Δ = s Δ 2 υ max (1.2.3) 式中：ω-角频率； xK , yK -波矢量在 x 和 y 方向的分量；υ-波速， maxυ 为模型最大速 度。 1.3 正演模拟实例 通 过一 个具 体实 例， 阐述 正演 模拟 的效 果。 模型参 数为 ：磁 导率 为真 空介 质 下 75 0μ =4π×10-7F/m，其它参数如下表 2。 表 2 三层模型参数 Tab.2 The parameter of three layers of model 模型层位 6 15 30 相对介电常数 电导率(s/m) 0.0035 0.008 0.05 0.1224 0.077 0.054 速度（m/ns） 顶界面深度（米） 到达界面时间（ns） 0 1.5 4 6 15 30 0 29 80 图 a 三层地电模型 图 b 雷达合成记录 Fig.a Three layers of terrestrial electricity model Fig.b The synthetical record of radar Fig.2 Three layers of terrestrial electricity model and its synthetical record 图 2 三层地电模型及合成记录 图 1 为正演模拟中的各个时刻波场快照图，图 2.a 是三层正演模型，图 2.b 是正演模拟 后的合成记录。该模型是中间震源激发，两边接收信号。通过激励源位置下界面处顶界面深 65 70 80 85 90 图 1 各个时刻波场快照图 Fig.1 the shot snap of each time - 3 - 

中国科技论文在线 http://www.paper.edu.cn 度和时间关系公式 v=2h /t 可以验证由剖面得到的速度符合模型参数。 2 图像处理 本文对雷达模拟剖面做了以下处理（本文重点介绍 FK 滤除直达波）： （1）去直流漂移，去除时间常数漂移。 （2）AGC，增强绕射波能量。 （3）抽取平均，抑制水平一致能量。 （4）巴特沃斯带通滤波，滤除低频和高频噪声。 2.1 二维傅氏变换 以 ( , ) y t x 表示多道雷达记录， ( Y f k 表示二维 FK 谱，则它的二维 FK 正反变换分别 , ) 95 100 为公式（2.1.1）(2.1.2)。 Y f k , ) ( y t x e ( , ) − i 2 ( π ft kx + ) dtdx (2.1.1) ∞ ∞ = ∫ ∫ −∞ −∞ ∞ ∞ y t x ( , ) = ∫ ∫ −∞ −∞ Y f k e , ) i ( 2 ( π ft kx + ) dfdk (3.1.2) 2.2 视速度的概念 105 由于波数 k 满足公式 k = 1 λ = 1 vt = ，且雷达探测中通常是沿测线观测的，因侧，用 f v 视 波长 （ λ∗ ） 、视 波数 （ k ∗ ） 和视 速度 （ v∗ ） 等来 描述所 观测 的雷 达波 的特 征 即 ∗ k = 1 ∗ λ = 1 v t ∗ = f v ∗ ，因此雷达资料构成了频率、波数、视速度的关系[7]。 2.3 视速度滤波（倾角滤波） 雷达在探测应用中通常为寻找异常体如管道、溶洞等，提取这些异常体所产生的散射波 110 能为后期解释提高精度和分辨率[9]。根据视速度和视波数的关系公式，可得 v ∗ = ，则可 f k ∗ 以根据视速度滤除不同的同向轴，下面图 2 为滤除直达波，从雷达资料提取散射特征。 图 a 圆形异常体模型 图 b 雷达合成原始记录 Fig.a The circle of unnormal structure Fig.b The synthetical record of radar - 4 - 

中国科技论文在线 http://www.paper.edu.cn Fk 滤波扇 形区域 115 120 125 130 135 图 c FK 扇形滤波器 图 d 切除直达波剖面 Fig.c F-K fan-shaped filtering Figure Fig.d the cut direct path wave section 图 3 圆形二度体波场模型、原始剖面、二维滤波器和切除直达波剖面 Fig.3 The round two-dimensional model、raw section、two-dimensional Filtering and direct path wave cut section 从图 3d 可以看出直达波已经被切除，但波形中心位置也受到损伤，因为波形中间位置 的视速度很大，刚好在扇形滤波器的范围内。因此 FK 切除直达波会损伤绕射弧的中间部分。 3 结论 本论文研制了一套实用、操作简便、绘制多种规则和不规则异常体、可视化的雷达正演 软件，为后序探地雷达研究提供了有利的条件，同时编制雷达数据处理软件极大的提高了处 理的灵活性和精确性，为地质雷达的各种应用性研究提供了理论依据。 [参考文献] (References) [1] 曾昭发,刘四新.探地雷达方法原理及应用[M]，科学出版社，2006.6 [2] 李宁.完美匹配层理论及其在地震波模拟中的应用[D]，中国地震局工程力学研究所，2006.5 [3] 王春燕.高阶交错网格有限差分地震波场计算[D]，成都理工大学，2007.4 [4] 何兵寿.地质雷达的模型计算及应用[D]，中国矿业大学，1999.5 [5] 冯德山,戴前伟.地质雷达二维时域有限差分正演[D]，勘察科学技术，2004 第 6 期 [6] 吴宝杰.探地雷达二维正演模拟及其工程实例[D]，浙江大学，2007.6 [7] 李振春,张军华.地震数据处理方法[M],中国石油大学出版社，2006.12 [8] 李昂,蒋延生.自适应对消在去除探地雷达信号直达波的应用[J]，电波科学学报 2004.4 - 5 -