用户名: 密码: 验证码:
探地雷达二维有限元正演模拟
详细信息    本馆镜像全文|  推荐本文 |  |   获取CNKI官网全文
摘要
探地雷达(Ground Penetrating Radar, GPR)是探测地下结构及其分布规律的一种重要浅层地球物理勘探方法,探地雷达正演模拟则是探地雷达研究的一个重要课题,其模拟结果对野外实际工作、室内数据处理与解释具有重要的指导意义。目前探地雷达正演模拟的方法主要有射线追踪法、时域有限差分法和有限单元法,有限单元法的优越性体现在它能模拟较复杂的模型,且有更好的稳定性。故本文选择有限单元法进行探地雷达正演模拟。
     论文研究得到了国家自然科学基金项目“基于小波有限元探地雷达正演及偏移处理研究”(40804027)、国家自然科学基金项目“复杂GPR模型无网格法正演及多参数混合智能优化反演(41074085)”、湖南省自然基金重点项目(09JJ3084)等科研项目的支持。其主要内容如下:
     1.从Maxwell方程组出发,推导了探地雷达有限元波动方程,给出了伽辽金法求解该微分问题的推导过程,并详细阐述了各主要步骤的细节处理方法,如:质量矩阵、衰减系数矩阵、刚度矩阵的具体求解过程,激励源的加载办法、差分代替微商将探地雷达有限元方程转化成对线性方程组的措施等。
     2.引入了大型稀疏矩阵的存储方法、不完全LU分解技术的大型线性方程组求解的Biconjugate Gradient Stabilized (BICGSTAB)算法,并结合差分代替微商,讨论了时间步长和网格步长须满足的数值稳定性条件,讨论了集中质量矩阵法对求解有限元方程的影响及其优缺点。
     3.阐述了基于透射机理的透射边界条件和基于衰减机理的Sarma边界条件的原理,详细推导这两种边界条件的探地雷达有限元理论公式;通过在Sarma边界条件衰减层内加入过渡带,压制了介质区和衰减层交界面处的人为反射。以中心电磁脉冲源模型为例,对比了两种边界条件对人工截断边界的处理效果。
     4.考虑到透射边界条件与Sarma边界条件不同的理论机制,提出了一种结合透射边界条件和Sarma边界条件的混合边界条件,主要思想是使探地雷达波经过Sarma边界条件的衰减吸收后再通过透射边界条件将剩余能量透射出去,集成了二者的优势。以二维均匀模型的全波场快照为例,证明了该混合边界条件对到达截断边界处的探地雷达波的处理优于单一边界条件。
     5.编制了带衰减项与不带衰减项的探地雷达二维有限元正演模拟程序,并自定义了衰减比的概念及计算公式,通过对衰减媒质和非衰减媒质的模拟计算,说明了GPR波在衰减媒质中传播时会产生频散和被媒质吸收。因此,对于常见的衰减媒质中进行正演模拟时,运用带衰减项的雷达波方程进行正演模拟,显然更符合雷达波在地下介质中传播的实际情况。最后,通过对一系列的典型地电模型进行探地雷达正演模拟,加深了对探地雷达波在地下介质中的传播规律的理解和掌握,有助于更好地指导探地雷达资料的地质解释。
Ground penetrating radar (GPR) is an important shallow layer geophysical exploration method which detects the underground structure and the distribution law. The GPR numerical simulation is one of the most significant subjects in the GPR research. The GPR numerical simulation result can help us to work outside and data processing and geologic interpretation of radar data. Currently, the methods of GPR forward simulation are ray tracing method, finite difference time domain method (FDTD) and finite element method (FEM). The finite element method superiority manifests in the simulation of the complex model. So the FEM is used to GPR numerical simulation in this paper.
     The paper was funded by a joint support of the National Natural Science Fund "The research of GPR' s forward modeling and migration based on finite element of wavelet" (40804027),"The complex GPR element-free method forward simulation and multi-parameter hybrid intelligent optimization inversion"(41074085) and the Hunan province natural fund important project(09JJ3084). The main content of this paper is as follows:
     1. The FEM wave equation method was derived from the Maxwell system of equations. The process of the derivation about the differential coefficient equation which was solved by the Galerkin method was elaborated, such as how to calculate the mass matrix, the attenuation coefficient matrix, the stiffness matrix, the method of disposing the source, in order to transform the GPR finite element equation to system of linear equations uesing difference to instead of the differential coefficient, etc.
     2. The method to storage the large-scale sparse matrix was introduced. The Biconjugate Gradient Stabilized (BICGSTAB) algorithm based on the LU incompletely decomposition technology was introduced to solve the large linear systems. The rules of the time step and gridding step which are influence the numerical value stability were given out, discussed the influence of the lumped mass matrix to solve the FEM equation and its advantages and disadvantages.
     3. The principle of the transmitting boundary condition which is based on transmission mechanism and the Sarma boundary condition which is based on attenuation mechanism were described in this paper. The formula of the two kinds of boundary condition which are used in GPR FEM forward simulation were derivated in detail. Through joined the transitional layer in the Sarma boundary condition to suppress the articicial reflection which is engendered from the interface between the medium area and the damping region. Taking the center electromagnetic pulse source model as example, compared the results of the two kinds of boundary condition.
     4. Considered the different theoretical mechanism of the transmitting boundary condition and the Sarma boundary condition, a mixed boundary condition which combined with the two kinds of boundary condition was given out in this paper. The main idea is the power of the GPR waves were attenuated by the Sarma boundary condition and then transmitted the residual energy by the transmitting boundary condition. And take the two-dimensional equal model entire wave field snapshot as the example, proved the effect of the mixed boundary condition better than just one kind of boundary condition.
     5. Compiled the 2D GPR forward simulation program which including the attenuation term or not. Defined the concept and the formula of the attenuation ratio. Through to simulate the attenuation medium and non attenuation medium models, explained the GPR waves transmit in the attenuation medium were absorbed by the medium and the frequency dispersion. Therefore, when the medium is attenuation medium, forward simulation must use the equation which including the attenuation term, obviously conforms to the actual situation of the GPR waves which disseminates in the underground medium. Finally, the results of the typical geoelectric structures can help us to understand the propagation law of GPR waves, so as to guide geologic interpretation of GPR data.
引文
[1]曾昭发,刘四新,冯晅.探地雷达原理与应用.北京:电子工业出版社,2010
    [2]杜军.盾构隧道壁后注浆探测图像识别及沉降控制研究:[博士学位论文].上海:同济大学,2006
    [3]郭山红,孙锦涛,谢仁宏,等.穿墙生命探测技术研究.南京理工大学学报,2005,29(2):186~192
    [4]高立兵,王赞,夏明军.GPR技术在考古勘探中的应用研究.地球物理学进展,2000,15(1):61~69
    [5]时玉彬,杨子杰.海洋环境监测高频地波雷达波形参数设计.电讯技术,2001,(6):1~4
    [6]郭高轩,吴吉春.应用GPR获取水文地质参数研究初探.水文地质工程地质,2005,(1)89~93
    [7]周峻峰.钻孔雷达探测金属矿的数值模拟:[硕士学位论文].吉林:吉林大学,2009
    [8]崔祥斌,孙波,张向培,等.极地冰盖冰雷达探测技术的发展综述.极地研究,2009,21(4):322~335
    [9]粟毅,黄春琳,雷文太.探地雷达理论与应用.北京:科学出版社,2006
    [10]曾昭发,高尔根.三维介质中探地雷达(GPR)波传播逐段迭代射线追踪方法研究和应用.吉林大学学报(地球科学版),2005,35(S1):119~136
    [11]喻振华,冯德山,戴前伟,等.复杂地电模型的探地雷达时域有限差分正演.物探化探计算技术,2005,27(4):11~15+5
    [12]王惠濂.探地雷达目的体物理模拟研究结果.地球科学-中国地质大学学报,1993,18(3):266~284
    [13]李大心.探地雷达方法与应用.北京:地质出版社,1994
    [14]Yee K S. Numerical Solution of Initial Boundary Value Problems Involving Maxwell' s Equation in Isotripic Media. IEEETrans.Antennas Propagate.1966, 14,302~307
    [15]Taflove A. Re-inventing electromagnetics:supercomputing solution of Maxwell' s equations via direct time integration on space grids. Computing Systems in Engineering,1992,3 (1):153~168
    [16]Bergmann T, Robertsson J O.A, Holliger K. Numerical properties of staggered finite-difference solutions of Maxwell' s equations for ground-penetrating radar modeling. Geophysical Research Letters,1996,23 (1),45~48
    [17]Carcione, Jose M. Radiation patterns for 2-D GPR forward modeling. Geophysics,1998,63 (2):424~430
    [18]Kowalsky M B, Dietrich P, Teutsch G, etal. Forward modeling of ground-penetrating radar data using digitized outcrop images and multiple scenarios of water saturation. Water Resources Research,2001,37 (6): 1615~1625
    [19]GoodMan D. Ground-penetrating radar simulation in engineering and archaeology. Geophysics.1994,59 (2):224~232
    [20]王长清,祝西里.电磁场计算中的时域有限差分法.北京:北京大学出版社,1994
    [21]高本庆.时域有限差分法.北京:国防工业出版社,1995
    [22]余文华,彭仲秋,任朗.探地雷达的时域有限差分模型.应用科学学报,1995,4(13):491~494
    [23]方广有,张忠治,汪文秉.脉冲探地雷达的模拟计算.微波学报,1998,14(4):288~295
    [24]何兵寿,张会星.地质雷达正演中的频散压制和吸收边界条件改进方法.地质与勘探,2000,36(3):59~63
    [25]刘四新,佐藤源之.时间域有限差分法(FDTD)对井中雷达的数值模拟.吉林大学学报(地球科学版),2003,33(4):545~550
    [26]戴前伟,冯德山,王启龙,等.时域有限差分法在地质雷达二维正演模拟中的应用.地球物理学进展,2004,19(4):898~902
    [27]冯德山.基于小波多分辨探地雷达正演及偏移处理研究:[博士学位论文].长沙:中南大学,2006
    [28]冯德山,戴前伟,翁晶波.时域多分辨率法在探地雷达三维正演模拟中的应用.中南大学学报(自然科学版),2007,38(5):901~906
    [29]刘四新,曾昭发,徐波.三维频散介质中地质雷达信号的FDTD数值模拟.吉林大学学报(地球科学版),2006,36(1):123~127
    [30]刘四新,曾昭发.频散介质中地质雷达波传播的数值模拟.地球物理学报,2007,50(1):320~326
    [31]仝传雪,刘四新,王春辉.时域有限差分法(FDTD)模拟探地雷达极化测量.吉林大学学报(地球科学版),2006,(S1):201~205
    [32]徐世浙.地球物理中的有限单元法.北京:科学出版社,1994
    [33]杜世通.变速不均匀介质中波动方程的有限元数值解.华东石油学院学报,1982,6(2):1-20
    [34]牟永光.弹性波问题有限元数值解科研成果论文集.北京:华东石油学院,1987:9~231
    [35]柯本喜.一种有限元新网格自动形成方法.石油地球物理勘探,1988,23(3): 294~300
    [36]邵秀民,蓝志凌.非均匀各向同性弹性介质中地震波传播的数值模拟.地球物理学报,1995,38(S1):39~55
    [37]张明,王润秋.有限元解三维弹性波方程的并行算法.计算数学,1995,5(2):127~135
    [38]王妙月,郭亚曦,底青云.二维线性流变体波的有限元模拟.地球物理学报,1995,36(4):499~506
    [39]刘晶波,吕彦东.结构-地基动力相互作用问题分析的一种直接方法.土木工程学报,1998,31(3):55-64
    [40]张美根,王妙月,李小凡,等.各向异性弹性波场的有限元数值模拟.地球物理学进展,2002,17(3):384~389
    [41]马德堂,朱光明.有限元法与伪谱法混合求解弹性波动方程.地球科学与环境学报,2004,26(1):61~64
    [42]薛东川,王尚旭,焦淑静.起伏地表复杂介质波动方程有限元数值模拟方法.地球物理学进展,2007,22(2):522~529
    [43]王月英.地震波有限元集中质量矩阵及并行算法模拟研究:[博士学位论文].北京:中国石油大学,2007
    [44]沈飚.探地雷达波波动方程研究及其正演模拟.物探化探计算技术,1994,16(1):29-33
    [45]底青云,王妙月.雷达波有限元仿真模拟.地球物理学报,1999,42(6):818~825
    [46]底青云,许琨,王妙月.衰减雷达波有限元偏移.地球物理学报,2000,43(2):257~263
    [47]Qingyun Di, Miaoyue Wang. Migration of ground-penetrating radar data with a finite-element method that considers attenuation and dispersion. Geophysics, 2004,69 (2):472~477
    [48]Fanning P J, Boothby T E. Three-dimensional modeling and full-scale testing of stone arch bridges. Comput&Structures,2001,79 (29):2645~2662
    [49]谢辉,钟燕辉,蔡迎春.电磁场有限元法在GPR正演模拟中的应用.河南科学,2003,21(3):295~298
    [50]Arias P, Armesto J, Capua D D, etal. Digital photogrammetry, GPR and computational analysis of structural damages in a mediaeval bridge. Engineering Failure Analysis,2007,14 (8):1444~1457
    [51]Lubowiecka I, Armesto J, Arias P, etal. Historic bridge modelling using laser scanning, ground penetrating radar and finite element methods in the context of structural dynamics. Engineering Structures,2009,31 (11):2667~2676
    [52]杨峰,彭苏萍.地质雷达探测原理与方法研究.北京:科学出版社,2010
    [53]RAO S S. The finite element method in engineering. New York:Oxford New York Toronto Sydney Paris Frankfurt,1982:50~300
    [54]杜世通.变速不均匀介质中波动方程的有限元法数值解.华东石油学院学报,1982,6(2):1~20
    [55]俞康胤,任甲祥.中心差分法时间步长选取原则——有限单元法合成地震记录制作中的稳定性.华东石油学院学报(自然科学版),1988,12(2):11-19
    [56]俞康胤.合成地震记录制作中的有限元法数值分析.华东石油学院学报,1982,4:11~28
    [57]朱建林,牟永光.有限元粘滞弹性波VSP地震模型.石油地球物理勘探,1994,29(6):724~732
    [58]廖振鹏.工程波动理论(第2版).北京:科学出版社,2002
    [59]Kriegsmann G., Morawetz C S. Numerical solutions of exterior problems with the reduced wave equations. Journal of Computational Physics,1978,28(2): 181~97
    [60]Bayliss A. Turkel E. Boundary conditions for the helmholtz equation in duct-like geometries. Int Assoc for Mathematics & Computers in Simulation, 1982,455~458
    [61]Engquist B. Majda A. Absorbing boundary conditions for numerical simulation of waves. Proceedings of the National Academy of Sciences of the United States of America,1977,74 (5):1765~1766
    [62]Mur, Gerrit. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations. IEEE Transactions on Electromagnetic Compatibility,1981,23 (4):377~382
    [63]Mei K K., Fang Jiayuan. Superabsorption--A method to improve absorbing boundary conditions. IEEE Transactions on Antennas and Propagation,1992, 40 (9):1001~1010
    [64]Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics,1994,114(2):185~200
    [65]熊章强,张大洲,秦臻,等.瑞雷波数值模拟中的边界条件及模拟实例分析.中南大学学报(自然科学版),2008,39(4):824~830
    [66]冯德山,陈承申,戴前伟.基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟.地球物理学报,2010,53(10):2484~2496
    [67]Claerbout J F. Imaging the Earth' s Interior. Blackwell Scientific Publications, 1985
    [68]牟永光.弹性波问题有限元数值解科研成果论文集.北京:华东石油学院,1987,9~231
    [69]Caughey T. Classical normal modes in damped linear systems. Journal of Applied Mechanics.1960,27:269~271
    [70]成谷,张宝金.反射地震走时层析成像中的大型稀疏矩阵压缩存储和求解.2008,23(3):674~680
    [71]陈德鹏.任意线源井地电法三维有限元正演研究:[硕士论文].长沙,中南大学,2009
    [72]David M. Young. Second-degree iterative methods for the solution of large linear systems. Journal of Approximation Theory,1972,5 (2):137~148
    [73]柳建新,蒋鹏飞,童孝忠,等.不完全LU分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用.中南大学学报(自然科学版),2009,40(2):484~491

© 2004-2018 中国地质图书馆版权所有 京ICP备05064691号 京公网安备11010802017129号

地址:北京市海淀区学院路29号 邮编:100083

电话:办公室:(+86 10)66554848;文献借阅、咨询服务、科技查新:66554700