2017, Vol. 32
地质雷达技术是用于研究水文、环境与气候变化等重大科学问题研究的最有效地球物理方法之一，并广泛应用于检测工程、道路和建筑内部结构，探测冻土层、沙漠沙丘演化特征，具有探测速度快、分辨率高和探测结果直观成像的优势.工程道路检测 (蔡晖, 2014; 廖红建等, 2016) 需要准确地获取病害体的空间位置，如空洞体的大小，顶底界面的深度.但由于受到浅地表介质不均匀，空洞顶底界面会在雷达剖面上会形成复杂的反射波和绕射波，难以准确地确定空洞体的大小、深度，无法实现高分辨率工程探测.地质雷达探测冻土层的目的是研究冻土层随着气候和海拔的时空变化规律 (Arcone et al., 1998, 2002；Wainstein et al., 2014；刘澜波和钱荣毅, 2015; Merz et al., 2015)，需要准确地获取冻土层内部横向变化规律，季节性融化层和多年冻土层的深度，但由于受到冻土层横向速度变化和内部孤立块体绕射的影响，雷达探测难以实现冻土层内部结构的高分辨率成像.地质雷达广泛应用于探测沙漠演化规律 (Bristow et al., 1996, 2000; Doolittle et al., 2006; Livingstone et al., 2007)，但由于受到沙丘内部含水量横向变化及内部隔层的影响，仅仅在雷达剖面上难以清晰地反映沙漠内部的真实结构和潜水面的起伏形态.
由于受到近地表层横向速度变化，孤立块体的绕射波和倾斜界面的复杂反射波等影响，使得雷达剖面不能真实反映被探测目标体的空间位置，目标体的分辨率较低，严重影响探测效果.偏移处理是将地下反射点的信息归位到它真实的位置处，提高成像分辨率的最重要技术.应用较好的偏移处理方法处理雷达数据，能从偏移剖面中准确地获取冻土层内部时空变化规律，工程建筑物内部真实结构以及沙漠内部层理的演化特征，同时对研究水文、环境和气候变化具有重要意义.弹性波逆时偏移方法 (Guitton et al., 2007) 较Kirchhoff偏移方法 (Sun et al., 2006; Sun and McMechan, 2011; Li and Fomel, 2013)，能有效地利用全波场信息，较单程波偏移方法 (Moran et al., 2000; Sakamoto et al., 2015)，无传播方向和倾角限制，不需要对波场做上、下行波分离，并能较好的利用回转波、多次波成像，能适应速度横向变化的介质，是一种成像精度较高的偏移方法.电磁波与弹性波有着相似的运动学特征，当介质的导电率较低 ( < 10 mS/m)，频率范围在105~109Hz时，地质雷达数据就可用反射地震资料处理方法进行处理 (Fisher et al., 1992).常规的Kirchhoff偏移、有限差分偏移和频率-波数域偏移已经广泛应用于雷达数据成像处理 (Feng and Sato, 2004; Sun et al., 2006; Brown et al., 2009; Sun and McMechan, 2011)，有限差分逆时偏移在地质雷达上的应用，最早由Fisher于1992年提出，他利用声波算法对GPR数据做叠后偏移.随后经Sanada和Ashida (1999)、Leuschen和Plumb (2001)、Zhou等 (2005)、Bradford (2015)和Liu等 (2014)等人的发展相继完善，现今的有限差分逆时偏移理论是完全基于麦克斯韦方程组，并且能适应多种天线模式.
本文为解决常规Kirchhoff偏移方法对地下目标体成像精度低，绕射不完全收敛的难题，借鉴前人有限差分逆时偏移的理论基础，系统推导出麦克斯韦方程组TM波方程的高阶有限差分格式、稳定性条件、PML边界条件，并编程实现有限差分逆时偏移算法.利用两层模型和空洞模型的逆时偏移剖面，论证了偏移算法的准确性.通过对比两层模型和空洞模型的Kirchhoff偏移和逆时偏移剖面，论证得出逆时偏移方法能有效地提高垂直界面和空洞底界面成像分辨率，逆时偏移剖面能清晰地反映两层模型和空洞模型的空间形态和内部结构信息.
1 逆时偏移方法逆时偏移由正向延拓、逆时延拓和成像三个步骤组成 (如图 1):1) 正向延拓是波场沿着时间轴的正方向传播到最大时刻并将运算结果保存下来；2) 逆时延拓是从接收波场的最后一个采样点开始，向着负时间方向延拓至零时刻，并读取正向延拓中相应时刻的波场值；3) 应用互相关成像条件得到地下真实的构造信息 (Yan and Liu, 2013)，公式为
|
图 1 正向延拓示意图 (a) 和逆时延拓示意图 (b) Figure 1 The diagram of forward continuation (a) and reverse continuation (b) |
|
(1) |
其中image (x, z) 是点 (x, z) 的成像结果：S(x, z, t) 表示时间域的正向延拓波场；R(x, z, t) 表示时间域的逆时延拓波场.
2 二维TM波正演逆时偏移算法的关键在于正演问题，本文中正演算法是从麦克斯韦方程组推导出的二维TM波方程 (葛德彪和闫玉波, 2002).逆时延拓的算法和正演算法相似，不同之处在于差分离散时，正演算法是用t时刻的波场值推导计算t+1时刻的波场值，而逆时延拓算法是用t+1时刻的波场值推导计算t时刻的波场值.
2.1 基本理论公式对二维问题，可设麦克斯韦旋度方程组 (式2、3) 中 (葛德彪和闫玉波, 2002) 所有物理量均与z坐标无关，公式为
|
(2) |
|
(3) |
其中：E为电场强度，单位为伏特/米 (V/m)；D为电通量密度，单位为库伦/米2(C/m2)；H为磁场强度，单位为安培/米 (A/m)；B为磁通量密度，单位为韦伯/米2(Wb/m2)；J为电流密度，单位为安培/米2(A/m2)；Jm为磁流密度，单位为伏特/米2(V/m2).即
|
(4) |
|
(5) |
|
(6) |
式中的ε表示介质介电系数，单位为法拉/米 (F/m)，μ表示磁导系数，单位为亨利/米 (H/m)，σ表示电导率，单位为西门子/米 (S/m)，σm为导磁率，单位为欧姆/米 (Ω/m)，Hx、Hy分别为磁场的x和y分量，Ez为电场分量.
二维TM波FDTD离散时，电场值 (Ez) 根据周围4个磁场值 (Hx、Hy) 完成更新，磁场值 (Hx、Hy) 根据周围2个电场值 (Ez) 完成更新 (图 2).在FDTD中E、H各分量空间节点与时间步长的取值为网格节点的整数和半整数倍，如表 1所示.
|
图 2 二维TM波的Yee元胞 Figure 2 The Yee cell of 2-Dimension TM wave |
|
表 1 TM波的YEE元胞中E、H各分量节点位置 Table 1 YEE cell the node position of each component |
对上面3、4、5式中偏导用中心差商来代替，可得到：
|
(6) |
式中
|
(7) |
式中
|
(8) |
式中m=(i, j).
|
真空中的σ=0，σm=0，ε=ε0=8.85×10-12 F/m，μ=μ0= 4π×10-7 H/m，CA、CB、CP和CQ均为引入的辅助参数，Δt为时间步长，Δx和Δy分别为x方向和y方向的空间步长，下标i和j表示空间网格位置，Ck(N)表示高阶交错网格有限差分系数 (Liu et al., 1998).
2.3 数值稳定性条件FDTD方法在执行如上式FDTD算法时，随着时间步长的增加，时间步长Δt与空间步长Δx、Δy之间必须要满足稳定性条件，按照Cournant稳定性条件，二维情况下FDTD算法中空间和时间间隔之间应满足的条件为
|
(9) |
若Δx=Δy=δ, 上式为
|
(10) |
数值色散对空间离散间隔的要求为
|
(11) |
其各个物理参数见表 2.
|
|
表 2 物理参数及数值大小 Table 2 Physical parameter and numerical magnitude |
相比其他边界吸收条件，PML边界条件 (葛德彪, 2002) 存在以下几点优势：
(1) 计算PML边界区域的值时，仅仅需要更新边界区域的电场和磁场值，而不需要更新整个网格区域的值 (图 3).
|
图 3 全空间计算示意图 Figure 3 The diagram of the whole space calculation |
(2) PML边界处可灵活的设置边界衰减层，即在x方向只有水平方向衰减 (图 4)，y方向只有垂直方向衰减 (图 4)，角点区域x方向和y方向都有衰减.
|
图 4 x方向和y方向衰减示意图 Figure 4 The diagram of attenuation of x and y direction |
如图 5a和6a所示，分别设置两层模型和空洞结构模型用于检验逆时偏移算法的准确性，应用正演算法得出速度模型的零偏移距雷达剖面，并用正确的速度模型计算得出Kirchhoff偏移剖面和逆时偏移剖面，比较两种偏移方法的优缺点.横向和纵向网格间距为0.02 m，边界层为50个网格单元，雷克子波主频为900 MHz，时间步长选取0.02 ns，采样点数为2500.逆时偏移方法不可避免的在近地表处产生低频噪声，选用拉普拉斯滤波方法滤除低频噪声 (Zhang and Sun, 2009).
|
图 5 两层模型的偏移成像 (a) 速度模型; (b) 切除直达波的零偏移距雷达剖面; (c) Kirchhoff偏移剖面; (d) 零偏移距逆时偏移剖面. Figure 5 The migration imaging of two layer model (a) Velocity model; (b) Zero-offset radar section after mute the direct wave; (c) Kirchhoff migration section; (d) Zero-offset RTM section. |
|
图 6 空洞模型的偏移成像 (a) 速度模型; (b) 切除直达波的零偏移距雷达剖面; (c) Kirchhoff偏移剖面; (d) 零偏移距逆时偏移剖面. Figure 6 The migration imaging of cavitas model (a) Velocity model; (b) Zero-offset radar section after mute the direct wave; (c) Kirchhoff migration section; (d) Zero-offset RTM section. |
如图 5a所示，设置两层速度模型，雷达波速度分别为0.20 m/ns和0.15 m/ns，利用有限差分正演算法计算得出零偏移距雷达剖面 (图 5b)，分别应用Kirchhoff偏移和零偏移距逆时偏移算法得出相应的偏移剖面 (图 5c和5d).如图 5b所示，原速度模型中界面的尖端点处 (图 5a中A、B、C) 会引起明显的绕射波 (A、B、C箭头).对底界面C处的成像效果显示，逆时偏移剖面与原速度剖面完全对应 (图 5d中B)，而Kirchhoff剖面存在偏移假象 (图 5c中B).如图 5c和5d所示，Kirchhoff偏移和逆时偏移剖面都存在偏移画弧的假象，相比Kirchhoff偏移剖面 (图 5c中C和D)，逆时偏移剖面中偏移画弧的假象明显较弱 (图 5d中C和D).逆时偏移能有效地利用回转波，多次波成像，在垂直界面位置处 (图 5a中A)，逆时偏移方法的成像效果明显要好 (图 5c和5d中A)，能对垂直界面清晰成像.
3.2 空洞模型如图 6a所示，设置空洞模型，空洞外层介质为混凝土结构，空洞内部为空气，雷达波速度分别为0.15 m/ns和0.3 m/ns，正演计算得出零偏移距雷达剖面 (图 6b)，并用Kirchhoff偏移和逆时偏移算法计算得出对应的偏移剖面 (图 6c和6d).空洞的顶底界面是明显的波阻抗界面，会形成明显的反射和绕射波 (图 6b中A和B).Kirchhoff偏移和逆时偏移都能对空洞的顶界面清晰成像 (图 6c和6d中A)，但Kirchhoff偏移对空洞底界面的成像较差，不能清晰地反映空洞底界面的深度位置 (图 6c中红箭头)，而逆时偏移剖面能清晰地反映空洞底界面的深度 (图 6d中B).在空洞体两侧，Kirchhoff偏移和逆时偏移剖面都存在偏移假象，相比与Kirchhoff偏移剖面 (图 6c中B、C、E)，逆时偏移剖面中偏移假象的能量明显较弱 (图 6d中B、C、E).
4 总结 4.1相比于应用声波方程的地质雷达有限差分正演算法，麦克斯韦方程组TM波方程的地质雷达有限差分正演算法，考虑了电导率等参数，真实地还原了引起电磁波反射的主要因素是介电常数的差异，使得正演模拟过程更加接近电磁波在地下的传播规律，并且算法简单，易于实现.
4.2Kirchhoff偏移成像易受陡倾界面的影响，难以实现高分辨率偏移成像.相比于Kirchhoff偏移，逆时偏移存在以下两点优势：1) 能较好的利用回转波成像，可实现陡倾界面的清晰成像；2) 利用正向延拓和逆时延拓波场的互相关成像，能使反射波准确归位，绕射波完全收敛，可实现地下目标体的高分辨率成像.
4.3应用有限差分零偏移距正演模拟算法计算得出两层模型和空洞模型的零偏移距雷达剖面，并用有限差分零偏移距逆时偏移算法计算得出对应的的逆时偏移剖面.结果表明两个模型的逆时偏移剖面与对应的原速度模型完全相同，有效地验证了有限差分零偏移距逆时偏移算法的准确性.
4.4两层模型和空洞模型的逆时偏移结果表明，逆时偏移方法能有效地提高垂直界面和空洞底界面的成像分辨率，逆时偏移剖面能清晰地反映两层模型和空洞模型的空间形态和内部结构信息.为改善陡倾角地层、沙漠与冻土内部结构及裂缝探测分辨率与可靠性，区分地下空洞高顶、底界面提供了新的技术保证.
致谢 感谢国家基金委重点项目 (项目号2145803) 和教育部基科研费优秀教师基本项目委员会 (项目号53200859531) 的支持，感谢给予我实验帮助的宋翱、刘皓男同学，向你们致以最崇高的敬意和最真挚的谢意.| [] | Arcone S A, Lawson D E, Delaney A J, et al. 1998. Ground-penetratinng radar reflection profiling of groundwater and bedrock in an area of discontinuous permafrost[J]. Geophysics, 63(5): 1573–1584. DOI:10.1190/1.1444454 |
| [] | Arcone S A, Prentice M L, Delaney A J. 2002. Stratigraphic profiling with ground-penetrating radar in permafrost:A review of possible analogs for Mars[J]. Journal of Geophysical Research, 107(E11): 5108. |
| [] | Bradford J. 2015. Reverse-time prestack depth migration of GPR data from topography for amplitude reconstruction in complex environments[J]. Journal of Earth Science, 26(6): 791–798. DOI:10.1007/s12583-015-0596-x |
| [] | Bristow C, Pugh J, Goodall T. 1996. Internal structure of aeolian dunes in Abu Dhabi determined using ground-penetrating radar[J]. Sedimentology, 43(6): 995–1003. DOI:10.1111/sed.1996.43.issue-6 |
| [] | Bristow C S, Bailey S D, Lancaster N. 2000. The sedimentary structure of linear sand dunes[J]. Nature, 406(6791): 56–59. DOI:10.1038/35017536 |
| [] | Brown J, Nichols J, Steinbronn L, et al. 2009. Improved GPR interpretation through resolution of lateral velocity heterogeneity:Example from an archaeological site investigation[J]. Journal of Applied Geophysics, 68(1): 3–8. DOI:10.1016/j.jappgeo.2008.08.014 |
| [] | Cai H. 2014. Application of ground penetrating radar in detection of tunnel karst[J]. Subgrade Engineering(1): 175–178. |
| [] | Doolittle J A, Jenkinson B, Hopkins D, et al. 2006. Hydropedological investigations with ground-penetrating radar (GPR):Estimating water-table depths and local ground-water flow pattern in areas of coarse-textured soils[J]. Geoderma, 131(3-4): 317–329. DOI:10.1016/j.geoderma.2005.03.027 |
| [] | Feng X, Sato M. 2004. Pre-stack migration applied to GPR for landmine detection[J]. Inverse Problems, 20(6): S99–S115. DOI:10.1088/0266-5611/20/6/S07 |
| [] | Fisher E, McMechan G A, Annan A P, et al. 1992. Examples of reverse-time migration of single-channel, ground-penetrating radar profiles[J]. Geophysics, 57(4): 577–586. DOI:10.1190/1.1443271 |
| [] | Ge D B, Yan Y B. 2002. Finite-difference Time-domain Method for Electromagnetic Waves (in Chinese)[M]. Xi'an:Xi'an University of Electronic Science and Technology Press. |
| [] | Guitton A, Kaelin B, Biondi B. 2007. Least-squares attenuation of reverse-time-migration artifacts[J]. Geophysics, 72(1): S19–S23. DOI:10.1190/1.2399367 |
| [] | Leuschen C J, Plumb R G. 2001. A matched-filter-based reverse-time migration algorithm for ground-penetrating radar data[J]. IEEE Transactions on Geoscience and Remote Sensing, 39(5): 929–936. DOI:10.1109/36.921410 |
| [] | Li S W, Fomel S. 2013. Kirchhoff migration using eikonal-based computation of traveltime source derivatives[J]. Geophysics, 78(4): S211–S219. DOI:10.1190/geo2012-0375.1 |
| [] | Liao H J, Zhu Q N, Zan Y W, et al. 2016. Detection of ballastless track diseases in high-speed railway based on ground penetrating radar[J]. Journal of Southwest Jiaotong University (in Chinese), 51(1): 8–13. |
| [] | Liu L B, Qian R Y. 2015. Ground penetrating radar:A critical tool in near-surface geophysics[J]. Chinese Journal of Geophysics (in Chinese), 58(8): 2606–2617. DOI:10.6038/cjg20150802 |
| [] | Liu S X, Lei L L, Fu L, et al. 2014. Application of pre-stack reverse time migration based on FWI velocity estimation to ground penetrating radar data[J]. Journal of Applied Geophysics, 107: 1–7. DOI:10.1016/j.jappgeo.2014.05.008 |
| [] | Liu Y, Li C C, Mou Y G. 1998. Finite-difference numerical modeling of any even-order accuracy[J]. Oil Geophysical Prospecting, 33(1): 1–10. |
| [] | Livingstone I, Wiggs G F S, Weaver C M. 2007. Geomorphology of desert sand dunes:A review of recent progress[J]. Earth Science Reviews, 80(3-4): 239–257. DOI:10.1016/j.earscirev.2006.09.004 |
| [] | Merz K, Maurer H, Buchli T, et al. 2015. Evaluation of ground-based and helicopter ground-penetrating radar data acquired across an alpine rock glacier[J]. Permafrost and Periglacial Processes, 26(1): 13–27. DOI:10.1002/ppp.v26.1 |
| [] | Moran M L, Greenfield R J, Arcone S A, et al. 2000. Multidimensional GPR array processing using Kirchhoff migration[J]. Journal of Applied Geophysics, 43(2-4): 281–295. DOI:10.1016/S0926-9851(99)00065-8 |
| [] | Sakamoto T, Sato T, Aubry P J, et al. 2015. Ultra-wideband radar imaging using a hybrid of Kirchhoff migration and Stolt F-K migration with an inverse boundary scattering transform[J]. IEEE Transactions on Antennas and Propagation, 63(8): 3502–3512. DOI:10.1109/TAP.2015.2431725 |
| [] | Sanada Y, Ashida Y. 1999. An imaging algorithm for GPR data[C].//Symposium on the Application of Geophysics to Engineering and Environmental Problems 1999. Tulsa, OK:Environment and Engineering Geophysical Society, 565-573. |
| [] | Sun R, McMechan G A. 2011. Prestack 2D parsimonious Kirchhoff depth migration of elastic seismic data[J]. Geophysics, 76(4): S157–S164. DOI:10.1190/1.3581359 |
| [] | Sun R, McMechan G A, Lee C S, et al. 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data[J]. Geophysics, 71(5): S199–S207. DOI:10.1190/1.2227519 |
| [] | Wainstein P, Moorman B, Whitehead K. 2014. Glacial conditions that contribute to the regeneration of Fountain Glacier proglacial icing, Bylot Island, Canada[J]. Hydrological Processes, 28(5): 2749–2760. DOI:10.1002/hyp.v28.5 |
| [] | Yan H Y, Liu Y. 2013. Acoustic prestack reverse time migration using the adaptive high-order finite-difference method in the time-space domain[J]. Chinese Journal of Geophysics, 56(3): 181–195. |
| [] | Zhang Y, Sun J. 2009. Practical issues in reverse time migration:True amplitude gathers, noise removal and harmonic source encoding[J]. First Break, 26: 29–35. |
| [] | Zhou H, Sato M, Liu H J. 2005. Migration velocity analysis and prestack migration of common-transmitter GPR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 43(1): 86–91. DOI:10.1109/TGRS.2004.839920 |
| [] | 蔡晖. 2014. 探地雷达在隧道岩溶探测中的应用[J]. 路基工程(1): 175–178. |
| [] | 葛德彪, 闫玉波. 2002. 电磁波时域有限差分方法[M]. . |
| [] | 廖红建, 朱庆女, 昝月稳, 等. 2016. 基于探地雷达的高铁无砟轨道结构层病害检测[J]. 西南交通大学学报, 51(1): 8–13. |
| [] | 刘澜波, 钱荣毅. 2015. 探地雷达:浅表地球物理科学技术中的重要工具[J]. 地球物理学报, 58(5): 2606–2617. DOI:10.6038/cjg20150802 |