跨孔雷达走时层析成像反演方法的研究
|
摘要
电磁波层析成像技术（Electromagnetic Tomography）是常用的跨孔测量数据的解释方法，其属性理论基于Radon变换与Radon变换的逆变换，即根据在物体外部的测量数据，依据一定的物理和数学关系反演物体内部物理量的分布，并由计算机以图像形式显示出来。根据观测数据及反演目的的不同，层析成像方法可分为走时层析成像和衰减层析成像两种。其中走时层析成像利用直达波的初至时数据来反演两个钻孔间的电磁波速度分布；而衰减层析成像是利用直达波的振幅或重心频率信息来反演两个钻孔间的衰减场。
层析成像方法能识别那些引起明显物理性质差异（如介电常数、电导率或磁导率）的地下结构（如分层、断裂和埋藏的公共设施等）。对于低电损耗介质，可以认为介电常数只与电磁波速度有关，而与衰减系数无关，介电常数可以通过单独进行走时层析成像而求得；而对于高损耗介质，介电常数则由电磁波速度和衰减系数共同决定。另一个进行衰减层析成像的原因是通过衰减系数来推导所探测介质的电导率结构。而进行跨孔雷达走时和衰减的联合层析成像则能够将介电常数以及电导率对电磁波传播的影响分开来研究，这能极大地增强对所探测介质的理解。
全波形反演算法由于计算量巨大及天线特性的限制，很难得到更广泛的应用；而衰减层析成像算法容易受到辐射模式、几何扩散等因素的影响，反演结果往往不理想；走时层析成像利用直达波波形的初至时进行反演，而初至时的位置与辐射模式、天线特性等无关，并且走时层析成像的计算效率相对全波形反演算法要高得多。因此，在本文中，我们主要进行了跨孔雷达走时层析成像理论的研究。
当速度场相对变化较小时，我们可以近似认为射线在地下介质中传播的路径是直射线。此时，走时方程是线性的，基于直射线追踪的走时层析成像的反演可通过最小二乘非迭代反演算法实现，实现过程中，需要通过初至时提取来得到走时的观测值，并用直射线追踪算法构建系数矩阵（雅可比矩阵）；当速度场相对变化较大时，射线在地下介质中传播的路径是弯曲射线。此时，走时方程是非线性的，基于弯曲射线追踪的走时层析成像的反演需要使用最小二乘迭代线性反演算法计算，实现过程中，首先需要通过初至时提取来得到走时的观测值；在每次迭代中，需要使用MSFM算法计算走时的计算值，并用基于MSFM和最速下降法的弯曲射线追踪算法构建系数矩阵（雅可比矩阵）；无射线追踪跨孔雷达走时层析成像可使用与基于弯曲射线追踪的走时层析成像相同的最小二乘迭代线性反演算法实现，在实现过程中，首先需要通过初至时提取来得到走时的观测值；在每次迭代中，需要使用MSFM算法计算走时的计算值，并用有限差分微扰的方式构建雅可比矩阵。
振幅比法衰减层析成像原理简单，但是需要考虑几何扩散及天线辐射模式等的影响；重心频率下移法衰减层析成像由频率域数据出发，无需考虑几何扩散、仪器响应、源/接收耦合、辐射模式、发射/透射系数以及由传播引起的相位累加等的影响。通过公式变形，振幅比法和重心频率下移法衰减层析成像具有与走时层析成像相似的反演方程，可使用最小二乘非迭代反演算法和迭代线性反演算法实现。
基于程函方程有限差分解的走时计算方法有FMM、HAFMM和MSFM等。其中FMM和HAFMM算法分别使用一阶和二阶精度的差分近似求解程函方程，而MSFM算法使用两个模板计算邻点走时，同时考虑了水平垂直方向及对角线方向上的信息，在理论上能提高走时的计算精度和计算效率。
在对射线追踪技术和FMM算法及其改进算法研究的基础上，本文提出了一种新的射线追踪方法。该算法将射线追踪分为用MSFM算法正向计算走时和用最速下降法反向追踪射线路径两个过程。与基于FMM和最速下降法、基于HAFMM和最速下降法的射线追踪方法相比，对于匀速模型和常速度梯度模型，基于MSFM和最速下降法的射线追踪方案计算的走时等值线与理论值差异最小，射线路径与理论路径最接近，射线路径传播时间误差最小，从而证明了本射线追踪方案的有效性和高精度性；而对于复杂随机介质模型，基于MSFM和最速下降法的射线追踪方法在不同网格距下计算的两组走时等值线以及两组射线路径的差异最小，证明了本射线追踪方案的稳定性和对于复杂速度模型的适用性。
数字图像分割法最初应用在折射地震波数据初至时提取中，在本文中，我们首次实现了使用数字图像分割法对跨孔雷达直达波数据的初至时提取。数字图像分割法的分割是在能量比图像上进行的，使用凸集投影技术（POCS）进行图像分割，并在分割时加入一些限制条件。分割后的能量比振幅图中初至区域的下缘的索引即为初至的位置。通过对简单模型、起伏界面模型和随机模型算例的分析可知，数字图像分割法在提取初至时的精度上要高于互相关法和信噪比最大法。
三种走时层析成像算法对于合成数据和实测数据的反演结果表明：直射线追踪走时层析成像的反演结果最差、误差最高；弯曲射线追踪走时层析成像的反演结果最接近真实模型，反演误差最小；无射线追踪走时层析成像的反演结果和弯曲射线追踪走时层析成像反演结果比较相近但稍差一些，反演误差也稍大一点。但是，无射线追踪走时层析成像提供了一种跨孔雷达走时层析成像的新方法和新尝试，并且其有效性和精度都能满足反演的需要。
在跨孔雷达走时层析成像中，加权因子、模型加权矩阵、走时计算方法、初至时提取方法、最小二乘反演算法、反演网格距、天线移动步长、射线覆盖角度等等都会对反演结果造成影响。
加权因子的最优值可以通过L曲线法确定，并且使用该最优值进行反演所得到的结果能够很好的平衡数据空间和模型空间之间的失配。综合考虑目标体重建效果和反演误差，拉普拉斯算子是最适合跨孔雷达走时层析成像反演的模型加权矩阵算子。
数字图像分割法提取初至时的精度要高于信噪比最大法和互相关法，并且当把数字图像分割法提取的初至时应用到跨孔雷达走时层析成像中时，所得到的反演结果也好于信噪比最大法和互相关法的情况。MSFM算法不但在走时计算精度上要高于FMM和HAFMM算法，使用其计算的走时进行走时层析成像的结果也是最好的。在跨孔雷达走时层析成像中，三种最小二乘算法（LSQR算法、GMRES算法和BICGSTAB算法）都可以很好的反演速度场，并且使用LSQR和BICGSTAB算法的反演结果更好一些。
减小反演网格距可以提高反演精度，但是反演精度的提高是以牺牲计算时间为代价的。但是，在基于二维弯曲射线追踪的走时层析成像算法中，这个计算时间的提高还是可以接受的。因此，在跨孔雷达走时层析成像中，我们可以通过减小反演网格距的尺寸来提高反演的精度。通过减小天线移动步长能获得更多的射线数（观测数据），从而能提高反演的精度；但是同时也会引入更多的大角度射线，从而降低了模型重建的效果。当我们减小天线移动步长并且只使用小角度射线进行反演时，我们即提高了反演精度，又能获得更好的反演结果。
Electromagnetic tomography technique is a commonly used interpretation method forcrosshole radar data. The theory of electromagnetic tomography is based on the Radontransform and inverse Radon transform, which inverse the distribution of the physicalparameters in the abnormal body according certain physical and mathematical relationship,and then the distribution of the physical parameters will be displayed by the computer in theform of image. According to the different observation data and inversion purposes,tomography technique can be divided into two types which are traveltime tomography andattenuation tomography, respectively. Traveltime tomography inverses the distribution ofvelocity between two boreholes using the first arrival data of direct wave. Attenuationtomography inverses the distribution of attenuation between two boreholes using theamplitude or centorid frequency data of direct wave.
Tomography method can identify the subsurface structure (i.e., bedding, fracture andburied utilities) that induce significant physical property contrasts (dielectric permittivity,electrical conductivity, or magnetic permeability). For medium with low dielectric loss, thedielectric permittivity can be considered only concerned with electromagnetic wave velocity,and has nothing to do with the attenuation coefficient; dielectric permittivity can be obtainedthrough the traveltime tomography alone. And for the high dielectric loss medium, thedielectric permittivity is decided by both the electromagnetic wave velocity and theattenuation coefficient. Another reason for attenuation tomography is to obtain theconductivity structure of the detecting medium through the attenuation coefficient. The jointcrosshole radar traveltime and attenuation tomography can study the effects of permittivityand conductivity on the propagation of electromagnetic wave separately, which can greatlyimprove the understanding to the detection medium.
The huge amount of computation and antenna characteristics limit the applications ofthe full wave form tomography techniques. The inversion process of attenuation tomographymay be effect by radiation patterns, geometric spreading and other factors which lead to a badresult. Traveltime tomography use the first arrival data of direct wave to inverse the velocityfield, which are not related with the radiation patterns, antenna characteristics, and so on. Andthe computational efficiency of traveltime tomography is much higher than full waveforminversion methods. Therefore, in this paper, we mainly studied the theory of crosshole radartraveltime tomography.
Straight rays may be adequate if a medium is characterized by smooth and negligiblysmall velocity variations. At this time, the equation of traveltime is linear, and the traveltimetomography based on straight raytracing may be realized through a least squares linearizednon-iteratively inversion scheme. During the inversion, the observation data of traveltimescan be obtained by the first arrival extraction, and the Jacobi matrix can be constructed byusing the straight raytracing technique. In case of strong variations, using curved rays isnecessary to obtain accurate results. At this time, the equation of traveltime is non-linear, andthe traveltime tomography based on curved raytracing may be realized through a leastsquares linearized iteratively inversion scheme. During the inversion, we should extract thefirst arrival time firstly to get the observation data of traveltime. In each iteration, thecalculated data of traveltime can be calculated using the MSFM algorithm, and the Jacobimatrix can be constructed by using the curved raytracing technique based on the MSFMalgorithm and the steepest descent technique. The traveltime tomography without raytracingmay be realized through a least squares linearized iteratively inversion scheme which are thesame as the one used in curved raytracing traveltime tomography. During the inversion, theobservation data of traveltimes can be obtained by the first arrival extraction. In eachiteration, the calculated data of traveltime can be calculated using the MSFM algorithm andthe Jacobi matrix can be constructed by using the finite difference approximate.
In the attenuation tomography based on amplitude ratio, amplitudes are easilycontaminated by factors such as scattering, geometric spreading, source and receivercoupling, radiation patterns, and transmission and reflection effects. It can be difficult toobtain reliable attenuation estimates from time-domain data. In the attenuation tomographyusing the centroid frequency downshift method, the frequency shift or pulse broadening of anEM pulse is not affected by far-field geometrical spreading or reflection losses, appears to bemore reliable than the time-domain amplitude-decay methods, and can be easily implementedby a tomography algorithm. The attenuation tomography can be realized through a leastsquares linearized iteratively inversion scheme which are the same as the one used intraveltime tomography.
The fast marching method (FMM) uses a first-order approximation to solve the eikonalequation, which makes the FMM algorithm having a low accuracy. The higher accuracy fastmarching method (HAFMM) which solves the eikonal equation using a second orderapproximation which can improve the accuracy of FMM. However, both FMM and HAFMMignore the information provided by diagonal points and may suffer from a large numericalerror along the diagonal direction. The multistencils fast marching method (MSFM) computes the solution at each grid point by solving the eikonal equation along two stencilsthat cover its entire neighbor points and then picks the solution that satisfies the upwindcondition, which can improve the accuracy of FMM and HAFMM greatly.
We presented a ray-based iteratively traveltime tomography algorithm for crossholeradar direct-arrival data using the multistencils fast marching method (MSFM) and thesteepest descent technique. The proposed scheme used MSFM to compute the traveltimesolution at each grid point by solving the traveltime eikonal equation along several stencilsand picked the solution that satisfies the upwind condition. In contrast to classical fastmarching method (FMM) and the higher accuracy fast marching method (HAFMM), MSFMis highly accurate for forward modeling of traveltime. Curved raypaths, which were neededfor the construction of the Jacobi matrix during inversion, were generated by following thesteepest descent direction through a computed traveltime field from each receiver back to thesource using the steepest descent technique. In order to verify the accuracy and efficiency ofthe new ray tracing method, we test the proposed scheme on three synthetic velocity models,and we compared our results with the one obtained by a FMM based and a HAFMM basedsteepest descend ray tracing methods. This comparison indicated that the suggested raytracing technique is efficiency and can achieve much better results both on accuracy andefficiency compared to the FMM based and the HAFMM based steepest descend ray tracingmethods.
Digital image segmentation is an accurate extraction method for first arrival times whichwas first used in refracted seismic data. Digital image segmentation method extracts firstarrival times by segment the color image of the energy ratio based on the projection ontoconvex sets (POCS). We first applied the digital image segmentation method into crossholeradar first arrival extraction. We employed three synthetic data sets to test the digital imagesegmentation method and the traditional signal to noise ratio method and correlation method.The contrast indicated that the digital image segmentation method is more accurate withsmaller residual which can provide more effective help for the traveltime tomography.
The comparison of synthetic data set and field data set inversion results using threetypes of traveltime tomography algorithms indicate that the traveltime tomography algorithmwithout raytracing is able to generate a solution as good as the one resulting from a curvedraytracing scheme.
In traveltime tomography, inversion results are easily affected by factors such asweighting operator, weighting matrix, traveltime calculation methods, first arrival extractionmethods, least square inversion algorithms, inversion grid size, antenna stepping, ray coverage angle, and so on.
The L-curve approach can be used to obtain an optimal value for the weighting factorwhich can also lead to a good misfit between the data space and model space during thetraveltime tomography. When considering both the reconstruction result and residual, theLaplace operator is the most suitable weighting operator for crosshole radar traveltimetomography.
The accuracy for the extracted first arrival using the digital image segmentation methodis higher than the one using the max signal-to-noise ratio method and the cross-correlationmethod, and the inversion result based on the extracted first arrival using the digital imagesegmentation method is the best during the three first arrival extraction methods. Theaccuracy for both traveltime calculating and traveltime tomography using the MSFMalgorithm is higher than the FMM and HAFMM algorithms. In crosshole traveltimetomography, the three kinds of least square algorithm (LSQR, GMRES and BICGSTAB) canget very good reconstructed velocity field.
Reducing inversion grid spacing can improve the retrieval accuracy, but the inversionaccuracy is achieved at the expense of computing time. However, in the two-dimensionaltraveltime tomography algorithm based on curved raytracing, the increase of the computingtime is acceptable. Therefore, we can improve the inversion accuracy by reducing the size ofthe inversion grid spacing. By reducing the moving step of antenna, we can get more numberof rays (observed data), which can improve the accuracy of inversion; but at the same time, itwill also introduce more high angle rays which can affect the inversion effect. When we bothreduce the antenna moving step and only using the small angle rays for inversion, we can notonly improve the inversion precision, but also obtain better inversion results.
层析成像方法能识别那些引起明显物理性质差异（如介电常数、电导率或磁导率）的地下结构（如分层、断裂和埋藏的公共设施等）。对于低电损耗介质，可以认为介电常数只与电磁波速度有关，而与衰减系数无关，介电常数可以通过单独进行走时层析成像而求得；而对于高损耗介质，介电常数则由电磁波速度和衰减系数共同决定。另一个进行衰减层析成像的原因是通过衰减系数来推导所探测介质的电导率结构。而进行跨孔雷达走时和衰减的联合层析成像则能够将介电常数以及电导率对电磁波传播的影响分开来研究，这能极大地增强对所探测介质的理解。
全波形反演算法由于计算量巨大及天线特性的限制，很难得到更广泛的应用；而衰减层析成像算法容易受到辐射模式、几何扩散等因素的影响，反演结果往往不理想；走时层析成像利用直达波波形的初至时进行反演，而初至时的位置与辐射模式、天线特性等无关，并且走时层析成像的计算效率相对全波形反演算法要高得多。因此，在本文中，我们主要进行了跨孔雷达走时层析成像理论的研究。
当速度场相对变化较小时，我们可以近似认为射线在地下介质中传播的路径是直射线。此时，走时方程是线性的，基于直射线追踪的走时层析成像的反演可通过最小二乘非迭代反演算法实现，实现过程中，需要通过初至时提取来得到走时的观测值，并用直射线追踪算法构建系数矩阵（雅可比矩阵）；当速度场相对变化较大时，射线在地下介质中传播的路径是弯曲射线。此时，走时方程是非线性的，基于弯曲射线追踪的走时层析成像的反演需要使用最小二乘迭代线性反演算法计算，实现过程中，首先需要通过初至时提取来得到走时的观测值；在每次迭代中，需要使用MSFM算法计算走时的计算值，并用基于MSFM和最速下降法的弯曲射线追踪算法构建系数矩阵（雅可比矩阵）；无射线追踪跨孔雷达走时层析成像可使用与基于弯曲射线追踪的走时层析成像相同的最小二乘迭代线性反演算法实现，在实现过程中，首先需要通过初至时提取来得到走时的观测值；在每次迭代中，需要使用MSFM算法计算走时的计算值，并用有限差分微扰的方式构建雅可比矩阵。
振幅比法衰减层析成像原理简单，但是需要考虑几何扩散及天线辐射模式等的影响；重心频率下移法衰减层析成像由频率域数据出发，无需考虑几何扩散、仪器响应、源/接收耦合、辐射模式、发射/透射系数以及由传播引起的相位累加等的影响。通过公式变形，振幅比法和重心频率下移法衰减层析成像具有与走时层析成像相似的反演方程，可使用最小二乘非迭代反演算法和迭代线性反演算法实现。
基于程函方程有限差分解的走时计算方法有FMM、HAFMM和MSFM等。其中FMM和HAFMM算法分别使用一阶和二阶精度的差分近似求解程函方程，而MSFM算法使用两个模板计算邻点走时，同时考虑了水平垂直方向及对角线方向上的信息，在理论上能提高走时的计算精度和计算效率。
在对射线追踪技术和FMM算法及其改进算法研究的基础上，本文提出了一种新的射线追踪方法。该算法将射线追踪分为用MSFM算法正向计算走时和用最速下降法反向追踪射线路径两个过程。与基于FMM和最速下降法、基于HAFMM和最速下降法的射线追踪方法相比，对于匀速模型和常速度梯度模型，基于MSFM和最速下降法的射线追踪方案计算的走时等值线与理论值差异最小，射线路径与理论路径最接近，射线路径传播时间误差最小，从而证明了本射线追踪方案的有效性和高精度性；而对于复杂随机介质模型，基于MSFM和最速下降法的射线追踪方法在不同网格距下计算的两组走时等值线以及两组射线路径的差异最小，证明了本射线追踪方案的稳定性和对于复杂速度模型的适用性。
数字图像分割法最初应用在折射地震波数据初至时提取中，在本文中，我们首次实现了使用数字图像分割法对跨孔雷达直达波数据的初至时提取。数字图像分割法的分割是在能量比图像上进行的，使用凸集投影技术（POCS）进行图像分割，并在分割时加入一些限制条件。分割后的能量比振幅图中初至区域的下缘的索引即为初至的位置。通过对简单模型、起伏界面模型和随机模型算例的分析可知，数字图像分割法在提取初至时的精度上要高于互相关法和信噪比最大法。
三种走时层析成像算法对于合成数据和实测数据的反演结果表明：直射线追踪走时层析成像的反演结果最差、误差最高；弯曲射线追踪走时层析成像的反演结果最接近真实模型，反演误差最小；无射线追踪走时层析成像的反演结果和弯曲射线追踪走时层析成像反演结果比较相近但稍差一些，反演误差也稍大一点。但是，无射线追踪走时层析成像提供了一种跨孔雷达走时层析成像的新方法和新尝试，并且其有效性和精度都能满足反演的需要。
在跨孔雷达走时层析成像中，加权因子、模型加权矩阵、走时计算方法、初至时提取方法、最小二乘反演算法、反演网格距、天线移动步长、射线覆盖角度等等都会对反演结果造成影响。
加权因子的最优值可以通过L曲线法确定，并且使用该最优值进行反演所得到的结果能够很好的平衡数据空间和模型空间之间的失配。综合考虑目标体重建效果和反演误差，拉普拉斯算子是最适合跨孔雷达走时层析成像反演的模型加权矩阵算子。
数字图像分割法提取初至时的精度要高于信噪比最大法和互相关法，并且当把数字图像分割法提取的初至时应用到跨孔雷达走时层析成像中时，所得到的反演结果也好于信噪比最大法和互相关法的情况。MSFM算法不但在走时计算精度上要高于FMM和HAFMM算法，使用其计算的走时进行走时层析成像的结果也是最好的。在跨孔雷达走时层析成像中，三种最小二乘算法（LSQR算法、GMRES算法和BICGSTAB算法）都可以很好的反演速度场，并且使用LSQR和BICGSTAB算法的反演结果更好一些。
减小反演网格距可以提高反演精度，但是反演精度的提高是以牺牲计算时间为代价的。但是，在基于二维弯曲射线追踪的走时层析成像算法中，这个计算时间的提高还是可以接受的。因此，在跨孔雷达走时层析成像中，我们可以通过减小反演网格距的尺寸来提高反演的精度。通过减小天线移动步长能获得更多的射线数（观测数据），从而能提高反演的精度；但是同时也会引入更多的大角度射线，从而降低了模型重建的效果。当我们减小天线移动步长并且只使用小角度射线进行反演时，我们即提高了反演精度，又能获得更好的反演结果。
Electromagnetic tomography technique is a commonly used interpretation method forcrosshole radar data. The theory of electromagnetic tomography is based on the Radontransform and inverse Radon transform, which inverse the distribution of the physicalparameters in the abnormal body according certain physical and mathematical relationship,and then the distribution of the physical parameters will be displayed by the computer in theform of image. According to the different observation data and inversion purposes,tomography technique can be divided into two types which are traveltime tomography andattenuation tomography, respectively. Traveltime tomography inverses the distribution ofvelocity between two boreholes using the first arrival data of direct wave. Attenuationtomography inverses the distribution of attenuation between two boreholes using theamplitude or centorid frequency data of direct wave.
Tomography method can identify the subsurface structure (i.e., bedding, fracture andburied utilities) that induce significant physical property contrasts (dielectric permittivity,electrical conductivity, or magnetic permeability). For medium with low dielectric loss, thedielectric permittivity can be considered only concerned with electromagnetic wave velocity,and has nothing to do with the attenuation coefficient; dielectric permittivity can be obtainedthrough the traveltime tomography alone. And for the high dielectric loss medium, thedielectric permittivity is decided by both the electromagnetic wave velocity and theattenuation coefficient. Another reason for attenuation tomography is to obtain theconductivity structure of the detecting medium through the attenuation coefficient. The jointcrosshole radar traveltime and attenuation tomography can study the effects of permittivityand conductivity on the propagation of electromagnetic wave separately, which can greatlyimprove the understanding to the detection medium.
The huge amount of computation and antenna characteristics limit the applications ofthe full wave form tomography techniques. The inversion process of attenuation tomographymay be effect by radiation patterns, geometric spreading and other factors which lead to a badresult. Traveltime tomography use the first arrival data of direct wave to inverse the velocityfield, which are not related with the radiation patterns, antenna characteristics, and so on. Andthe computational efficiency of traveltime tomography is much higher than full waveforminversion methods. Therefore, in this paper, we mainly studied the theory of crosshole radartraveltime tomography.
Straight rays may be adequate if a medium is characterized by smooth and negligiblysmall velocity variations. At this time, the equation of traveltime is linear, and the traveltimetomography based on straight raytracing may be realized through a least squares linearizednon-iteratively inversion scheme. During the inversion, the observation data of traveltimescan be obtained by the first arrival extraction, and the Jacobi matrix can be constructed byusing the straight raytracing technique. In case of strong variations, using curved rays isnecessary to obtain accurate results. At this time, the equation of traveltime is non-linear, andthe traveltime tomography based on curved raytracing may be realized through a leastsquares linearized iteratively inversion scheme. During the inversion, we should extract thefirst arrival time firstly to get the observation data of traveltime. In each iteration, thecalculated data of traveltime can be calculated using the MSFM algorithm, and the Jacobimatrix can be constructed by using the curved raytracing technique based on the MSFMalgorithm and the steepest descent technique. The traveltime tomography without raytracingmay be realized through a least squares linearized iteratively inversion scheme which are thesame as the one used in curved raytracing traveltime tomography. During the inversion, theobservation data of traveltimes can be obtained by the first arrival extraction. In eachiteration, the calculated data of traveltime can be calculated using the MSFM algorithm andthe Jacobi matrix can be constructed by using the finite difference approximate.
In the attenuation tomography based on amplitude ratio, amplitudes are easilycontaminated by factors such as scattering, geometric spreading, source and receivercoupling, radiation patterns, and transmission and reflection effects. It can be difficult toobtain reliable attenuation estimates from time-domain data. In the attenuation tomographyusing the centroid frequency downshift method, the frequency shift or pulse broadening of anEM pulse is not affected by far-field geometrical spreading or reflection losses, appears to bemore reliable than the time-domain amplitude-decay methods, and can be easily implementedby a tomography algorithm. The attenuation tomography can be realized through a leastsquares linearized iteratively inversion scheme which are the same as the one used intraveltime tomography.
The fast marching method (FMM) uses a first-order approximation to solve the eikonalequation, which makes the FMM algorithm having a low accuracy. The higher accuracy fastmarching method (HAFMM) which solves the eikonal equation using a second orderapproximation which can improve the accuracy of FMM. However, both FMM and HAFMMignore the information provided by diagonal points and may suffer from a large numericalerror along the diagonal direction. The multistencils fast marching method (MSFM) computes the solution at each grid point by solving the eikonal equation along two stencilsthat cover its entire neighbor points and then picks the solution that satisfies the upwindcondition, which can improve the accuracy of FMM and HAFMM greatly.
We presented a ray-based iteratively traveltime tomography algorithm for crossholeradar direct-arrival data using the multistencils fast marching method (MSFM) and thesteepest descent technique. The proposed scheme used MSFM to compute the traveltimesolution at each grid point by solving the traveltime eikonal equation along several stencilsand picked the solution that satisfies the upwind condition. In contrast to classical fastmarching method (FMM) and the higher accuracy fast marching method (HAFMM), MSFMis highly accurate for forward modeling of traveltime. Curved raypaths, which were neededfor the construction of the Jacobi matrix during inversion, were generated by following thesteepest descent direction through a computed traveltime field from each receiver back to thesource using the steepest descent technique. In order to verify the accuracy and efficiency ofthe new ray tracing method, we test the proposed scheme on three synthetic velocity models,and we compared our results with the one obtained by a FMM based and a HAFMM basedsteepest descend ray tracing methods. This comparison indicated that the suggested raytracing technique is efficiency and can achieve much better results both on accuracy andefficiency compared to the FMM based and the HAFMM based steepest descend ray tracingmethods.
Digital image segmentation is an accurate extraction method for first arrival times whichwas first used in refracted seismic data. Digital image segmentation method extracts firstarrival times by segment the color image of the energy ratio based on the projection ontoconvex sets (POCS). We first applied the digital image segmentation method into crossholeradar first arrival extraction. We employed three synthetic data sets to test the digital imagesegmentation method and the traditional signal to noise ratio method and correlation method.The contrast indicated that the digital image segmentation method is more accurate withsmaller residual which can provide more effective help for the traveltime tomography.
The comparison of synthetic data set and field data set inversion results using threetypes of traveltime tomography algorithms indicate that the traveltime tomography algorithmwithout raytracing is able to generate a solution as good as the one resulting from a curvedraytracing scheme.
In traveltime tomography, inversion results are easily affected by factors such asweighting operator, weighting matrix, traveltime calculation methods, first arrival extractionmethods, least square inversion algorithms, inversion grid size, antenna stepping, ray coverage angle, and so on.
The L-curve approach can be used to obtain an optimal value for the weighting factorwhich can also lead to a good misfit between the data space and model space during thetraveltime tomography. When considering both the reconstruction result and residual, theLaplace operator is the most suitable weighting operator for crosshole radar traveltimetomography.
The accuracy for the extracted first arrival using the digital image segmentation methodis higher than the one using the max signal-to-noise ratio method and the cross-correlationmethod, and the inversion result based on the extracted first arrival using the digital imagesegmentation method is the best during the three first arrival extraction methods. Theaccuracy for both traveltime calculating and traveltime tomography using the MSFMalgorithm is higher than the FMM and HAFMM algorithms. In crosshole traveltimetomography, the three kinds of least square algorithm (LSQR, GMRES and BICGSTAB) canget very good reconstructed velocity field.
Reducing inversion grid spacing can improve the retrieval accuracy, but the inversionaccuracy is achieved at the expense of computing time. However, in the two-dimensionaltraveltime tomography algorithm based on curved raytracing, the increase of the computingtime is acceptable. Therefore, we can improve the inversion accuracy by reducing the size ofthe inversion grid spacing. By reducing the moving step of antenna, we can get more numberof rays (observed data), which can improve the accuracy of inversion; but at the same time, itwill also introduce more high angle rays which can affect the inversion effect. When we bothreduce the antenna moving step and only using the small angle rays for inversion, we can notonly improve the inversion precision, but also obtain better inversion results.
引文
ALDRIDGE D F, OLDENBURG D W.1993. Two-dimensional tomographic inversion withfinite-difference traveltimes[J]. Journal of Seismic Exploration,2(3):257-274.
AL-GHAMDI A.2007. Automatic first arrival picking using energy ratios [D]. Saudi Arabia:King Fahd University of Petroleum&Minerals.
ALUMBAUGH D L, NEWMAN G A.2000. Image appraisal for2-D and3-Delectromagnetic inversion[J]. Geophysics,65(5):1455-1467.
AMMON C J, VIDALE J E.1993. Tomography without rays[J]. Bulletin of theSeismological Society of America,83(2):509-528.
ASAKAWA E, KAWANAKA T.1993. Seismic ray tracing using linear traveltimeinterpolation[J]. Geophysical prospecting,41(1):99-111.
BALANIS C A.1989. Advanced engineering electromagnetics[J].
BECHT A, TRONICKE J, APPEL E, ET AL.2004. Inversion strategy in crosshole radartomography using information of data subsets[J]. Geophysics,69(1):222-230.
BELINA F, IRVING J, ERNST J, ET AL.2012. Analysis of an iterative deconvolutionapproach for estimating the source wavelet during waveform inversion of crossholegeoradar data[J]. Journal of Applied Geophysics,78:20-30.
BERES M, HAENI F P.1991. Application of Ground-Penetrating-Radar Methods inHydrogeologie Studies[J]. Ground Water,29(3):375-386.
BREGMAN N D, BAILEY R C, Chapman C H.1989. Crosshole seismic tomography[J].Geophysics,54(2):200-215.
BROYDEN C G.1965. A class of methods for solving nonlinear simultaneous equations[J].Math. Comp,19(92):577-593.
CAI W, QIN F, SCHUSTER G T.1996. Electromagnetic velocity inversion using2-DMaxwell's equations[J]. Geophysics,61(4):1007-1021.
CAO S, GREENHALGH S.1993. Calculation of the seismic first-break time field and its raypath distribution using a minimum traveltime tree algorithm[J]. Geophysical JournalInternational,114(3):593-600.
CARDARELLI E.2002. Ray tracing in elliptical anisotropic media using the lineartraveltime interpolation (LTI) method applied to traveltime seismic tomography[J].Geophysical prospecting,50(1):55-72.
CARDIMONA S J, CLEMENT W P, KADINSKY-CADE K.1998. Seismic reflection andground-penetrating radar imaging of a shallow aquifer[J]. Geophysics,63(4):1310-1317.
CHENG N, HOUSE L.1996. Minimum traveltime calculation in3-D graph theory[J].Geophysics,61(6):1895-1898.
CLEMENT W P, BARRASH W.2006. Crosshole radar tomography in a fluvial aquifer nearBoise, Idaho[J]. Journal of Environmental&Engineering Geophysics,11(3):171-184.
CLEMENT W P, KNOLL M D.2000. Tomographic inversion of crosshole radar dataconfidence in results SAGEEP2000[C]. Environmental and Engineering GeophysicalSociety, Proceedings:553–562.
CONYERS L B.2006. Innovative ground-penetrating radar methods for archaeologicalmapping[J]. Archaeological Prospection,13(2):139-141.
COPPENS F.1985. First arrival picking on common-offset trace collections for automaticestimation of static corrections[J]. Geophysical Prospecting,33(8):1212-1231.
DANIELSSON P E, LIN Q.2003. A modified fast marching method[M]//Image Analysis.Springer Berlin Heidelberg,1154-1161.
DIJKSTRA E W.1959. A note on two problems in connexion with graphs[J]. Numerischemathematik,1(1):269-271.
DINES K A, LYTLE R J.1979. Computerized geophysical tomography[J]. Proceedings ofthe IEEE,67(7):1065-1073.
ENDRES A L, KNIGHT R.1991. The effects of pore‐scale fluid distribution on thephysical properties of partially saturated tight sandstones[J]. Journal of applied physics,69(2):1091-1098.
ENDRES A L, REDMAN J D.1996. Modelling the electrical properties of porous rocks andsoils containing immiscible contaminants[J]. Journal of Environmental and EngineeringGeophysics,1(B):105-112.
FENG S, SEN P N.1985. Geometrical model of conductive and dielectric properties ofpartially saturated rocks[J]. Journal of Applied Physics,58(8):3236-3243.
FISCHER R, LEES J M.1993. Shortest path ray tracing with sparse graphs[J]. Geophysics,58(7):987-996.
GAUTHIER O, VIRIEUX J, TARANTOLA A.1986. Two-dimensional nonlinear inversionof seismic waveforms: Numerical results[J]. Geophysics,51(7):1387-1403.
GIROUX B, BOUCHEDDA A, CHOUTEAU M.2009. Assisted traveltime picking ofcrosshole GPR data[J]. Geophysics,74(4): J35-J48.
GIROUX B, GLOAGUEN E, CHOUTEAU M.2007. bh_tomo—a Matlab borehole georadar2D tomography package[J]. Computers&Geosciences,2007,33(1):126-137.
GLOAGUEN E, CHOUTEAU M, MARCOTTE D, ET AL.2001. Estimation of hydraulicconductivity of an unconfined aquifer using cokriging of GPR and hydrostratigraphicdata[J]. Journal of Applied Geophysics,47(2):135-152.
GLOAGUEN E, MARCOTTE D, GIROUX B, ET AL.2007. Stochastic borehole radarvelocity and attenuation tomographies using cokriging and cosimulation[J]. Journal ofApplied Geophysics,62(2):141-157.
G KTüRKLER G, BALKAYA.2010. Traveltime tomography of crosshole radar datawithout ray tracing[J]. Journal of Applied Geophysics,72(4):213-224.
GREENHALGH S A, BING Z, GREEN A.2006. Solutions, algorithms and inter-relationsfor local minimization search geophysical inversion[J]. Journal of Geophysics andEngineering,3(2):101.
GREMAUD P A, KUSTER C M.2006. Computational study of fast methods for the eikonalequation[J]. SIAM journal on scientific computing,27(6):1803-1816.
HANAFY S, AL HAGREY S A.2005. Ground-penetrating radar tomography forsoil-moisture heterogeneity[J]. Geophysics,71(1): K9-K18.
HANSEN P C, O'LEARY D P.1993. The use of the L-curve in the regularization of discreteill-posed problems[J]. SIAM Journal on Scientific Computing,14(6):1487-1503.
HASSOUNA M S, FARAG A A.2007. Multistencils fast marching methods: A highlyaccurate solution to the eikonal equation on cartesian domains[J]. Pattern Analysis andMachine Intelligence, IEEE Transactions on,29(9):1563-1574.
HINKEL K M, DOOLITTLE J A, BOCKHEIM J G, ET AL.2001. Detection of subsurfacepermafrost features with ground‐penetrating radar, Barrow, Alaska[J]. Permafrost andPeriglacial Processes,12(2):179-190.
HOLE J A, ZELT B C.1995.3-D finite-difference reflection traveltimes[J]. GeophysicalJournal International,121(2):427-434.
HOLLENDER F.1999. Interprétation de la distorsion des signaux géoradar propagés etréflèchis. Développement d'une tomographie par bandes de fréquence[D].
HOLLIGER K, BERGMANN T.2002. Numerical modeling of borehole georadar data[J].Geophysics,67(4):1249-1257.
HOLLIGER K, MUSIL M, MAURER H R.2001. Ray-based amplitude tomography forcrosshole georadar data: A numerical assessment[J]. Journal of Applied Geophysics,47(3):285-298.
IRVING J D, KNOLL M D, KNIGHT R J.2007. Improving crosshole radar velocitytomograms: A new approach to incorporating high-angle traveltime data[J]. Geophysics,72(4): J31-J41.
JACKSON M J, TWEETON D R.1994. MIGRATOM-Geophysical tomography usingwavefront migration and fuzzy constraints[M]. US Department of Interior, Bureau ofMines.
JOHNSON T C, ROUTH P S, KNOLL M D.2005. Fresnel volume georadarattenuation-difference tomography[J]. Geophysical Journal International,162(1):9-24.
JULIAN B R, GUBBINS D.1977. Three-dimensional seismic ray tracing[J]. J. geophys,43(1):95-114.
KIM S, COOK R.1999.3-D traveltime computation using second-order ENO scheme[J].Geophysics,64(6):1867-1876.
KIM S.2001. An O(N) Level Set Method for Eikonal Equations[J]. SIAM journal onscientific computing,22(6):2178-2193.
KIRSCH A.1996. An Introduction to the Mathematical Theory of Inverse Problems [M].Basel: Springer.
KNIGHT R.2001. Ground penetrating radar for environmental applications[J]. AnnualReview of Earth and Planetary Sciences,29(1):229-255.
KURODA S, TAKEUCHI M, KIM H J.2005. Full waveform inversion algorithm forinterpreting cross-borehole GPR data[C]//2005SEG Annual Meeting. Society ofExploration Geophysicists.
LINES L R, TREITEL S.1984. Tutorial: A review of least-squares inversion and itsapplication to geophysical problems[J]. Geophysical prospecting,32(2):159-186.
LIU L, LANE J W, QUAN Y.1998. Radar attenuation tomography using the centroidfrequency downshift method[J]. Journal of Applied Geophysics,40(1):105-116.
LIU S, WU J, DONG H, ET AL.2012. The experimental results and analysis of a boreholeradar prototype[J]. Journal of Geophysics and Engineering,9(2):201.
MEBINGER J.2004. Effective automatic picking of traveltime data with highprecision[C]//Ground Penetrating Radar,2004. GPR2004. Proceedings of the TenthInternational Conference on. IEEE,1:91-94.
MELES G, GREENHALGH S, DER KRUK J, ET AL.2012. Taming the non-linearityproblem in GPR full-waveform inversion for high contrast media[J]. Journal of AppliedGeophysics,78:31-43.
MOGHADDAM M, CHEW W C, ORISTAGLIO M.1991. Comparison of the born iterativemethod and tarantola's method for an electromagnetic time‐domain inverse problem[J].International Journal of Imaging Systems and Technology,3(4):318-333.
MOLYNEUX J B, SCHMITT D R.1999. First-break timing: Arrival onset times by directcorrelation[J]. Geophysics,64(5):1492-1501.
MOSER T J.1991. Shortest path calculation of seismic rays[J]. Geophysics,56(1):59-67.
MOSER T J.1992. The shortest path method for seismic ray tracing in complicated media[J].Geologica Ultraiectina,83:1-180.
MOUSA W A, AL-SHUHAIL A A, AL-LEHYANI A.2011. A new technique for first-arrivalpicking of refracted seismic data based on digital image segmentation[J]. Geophysics,76(5): V79-V89.
MUSIL M, MAURER H R, GREEN A G.2003. Discrete tomography and joint inversion forloosely connected or unconnected physical properties: application to crosshole seismicand georadar data sets[J]. Geophysical Journal International,153(2):389-402.
NAKANISHI I, YAMAGUCHI K.1986. A numerical experiment on nonlinear imagereconstruction from first-arrival times for two-dimensional island arc structure[J].Journal of Physics of the Earth,34(2):195-201.
OLSSON O, FALK L, FORSLUND O, ET AL.1992. Borehole Radar Applied to theCharacterization of Hydraulically Conductive Fracture Zones in Crystalline ROCK1[J].Geophysical Prospecting,40(2):109-142.
PAIGE C C, SAUNDERS M A.1982. LSQR: An algorithm for sparse linear equations andsparse least squares[J]. ACM Transactions on Mathematical Software (TOMS),8(1):43-71.
PODVIN P, LECOMTE I.1991. Finite difference computation of traveltimes in verycontrasted velocity models: a massively parallel approach and its associated tools[J].Geophysical Journal International,105(1):271-284.
PRATT R G.1999. Seismic waveform inversion in the frequency domain, Part1: Theory andverification in a physical scale model[J]. Geophysics,64(3):888-901.
PROTHERO W A, TAYLOR W J, EICKEMEYER J A.1988. A fast, two-point,three-dimensional raytracing algorithm using a simple step search method[J]. Bulletin ofthe Seismological Society of America,78(3):1190-1198.
QIN F, LUO Y, OLSEN K B, ET AL.1992. Finite-difference solution of the eikonal equationalong expanding wavefronts[J]. Geophysics,57(3):478-487.
QUAN Y, HARRIS J M.1997. Seismic attenuation tomography using the frequency shiftmethod[J]. Geophysics,62(3):895-905.
SAAD Y, SCHULTZ M H.1986. GMRES: A generalized minimal residual algorithm forsolving nonsymmetric linear systems[J]. SIAM Journal on scientific and statisticalcomputing,7(3):856-869.
SABBIONE J I, VELIS D.2010. Automatic first-breaks picking: New strategies andalgorithms[J]. Geophysics,75(4): V67-V76.
SCHNEIDER JR W A.1995. Robust and efficient upwind finite-difference traveltimecalculations in three dimensions[J]. Geophysics,60(4):1108-1117.
SEN P N, SCALA C, COHEN M H.1981. A self-similar model for sedimentary rocks withapplication to the dielectric constant of fused glass beads[J]. Geophysics,46(5):781-795.
SETHIAN J A,1999. Leve Sets Methods and Fast Marching Methods, second ed[M].Cambridge Univ. Press.
SETHIAN J A, POPOVICI A M.1999.3-D traveltime computation using the fast marchingmethod[J]. Geophysics,64(2):516-523.
SETHIAN J A.1996. A fast marching level set method for monotonically advancing fronts[J].Proceedings of the National Academy of Sciences,93(4):1591-1595.
SLAWINSKI R A, SLAWINSKI M A.1999. On raytracing in constant velocity-gradientmedia: calculus approach[J]. Canadian journal of exploration geophysics,35:24-27.
SLEIJPEN G L G, FOKKEMA D R.1993. BiCGstab (l) for linear equations involvingunsymmetric matrices with complex spectrum[J]. Electronic Transactions on NumericalAnalysis,1(11):2000.
TARANTOLA A.1984. Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,49(8):1259-1266.
TIKHONOV A N, ARSENIN V I A, JOHN F.1977. Solutions of ill-posed problems[M].Washington, DC: Winston.
TRONICKE J, TWEETON D R, DIETRICH P, ET AL.2001. Improved crosshole radartomography by using direct and reflected arrival times[J]. Journal of AppliedGeophysics,47(2):97-105.
TRONICKE J.2007. The influence of high frequency uncorrelated noise on first-breakarrival times and crosshole traveltime tomography[J]. Journal of Environmental&Engineering Geophysics,12(2):173-184.
TSAI Y H R, CHENG L T, OSHER S, ET AL.2003. Fast sweeping algorithms for a class ofHamilton--Jacobi equations[J]. SIAM journal on numerical analysis,41(2):673-694.
UM J, THURBER C.1987. A fast algorithm for two-point seismic ray tracing[J]. Bulletin ofthe Seismological Society of America,77(3):972-986.
VAN TRIER J, SYMES W W.1991. Upwind finite-difference calculation of traveltimes[J].Geophysics,56(6):812-821.
VIDALE J E.1988. Finite-difference calculation of travel times[J]. Bulletin of theSeismological Society of America,78(6):2062-2076.
VIDALE J E.1990. Finite-difference calculation of traveltimes in three dimensions[J].Geophysics,55(5):521-526.
YATZIV L, BARTESAGHI A, SAPIRO G.2006. O (N) implementation of the fast marchingalgorithm[J]. Journal of computational physics,212(2):393-399.
ZHANG S, ZHOU Z, CHEN L, et al.2007. Seismic ray-tracing method of applyingparabolic Interpolation to travel-time[J]. Progress in Geophysics,1:004.
ZHANG Y T, ZHAO H K, QIAN J.2006. High order fast sweeping methods for staticHamilton–Jacobi equations[J]. Journal of Scientific Computing,29(1):25-56.
ZHAO H.2005. A fast sweeping method for eikonal equations[J]. Mathematics ofcomputation,74(250):603-627.
ZHOU B, FULLAGAR P K.2000. Borehole radar conductivity and permittivitytomography[C]//8th International Conference on Ground Penetrating Radar.International Society for Optics and Photonics,542-547.
ZHOU B, GREENHALGH S A.1998. Crosshole acoustic velocity imaging withfull-waveform spectral data:2.5-D numerical simulations[J]. Exploration Geophysics,29(3/4):680-684.
ZHOU C, LIU L, LANE JR J W.2001. Nonlinear inversion of borehole-radar tomographydata to reconstruct velocity and attenuation distribution in earth materials[J]. Journal ofApplied Geophysics,47(3):271-284.
ZHOU Z, ZHANG S, CHEN L.2004. Seismic ray-tracing calculation based on parabolictravel-time interpolation[J]. Journal of Central South University of Technology,11(2):199-205.
曾昭发,刘四新,王者江,等.2006.探地雷达方法原理及应用[M].北京：科学出版社.
成谷,张宝金,马在田.2007.浅谈反射地震走时层析中的正则化[J].地球物理学进展,22(6):1715-1721.
李灿苹,王南萍,李志宏.2010. VonKarman型自相关函数模拟随机介质[J].物探与化探,34(1):98-102.
李玉喜.2010.钻孔雷达层析成像软件系统的研究与开发[D].吉林：吉林大学.
佩特罗夫斯基.1981.地下无线电波法[M].北京:地质出版社:109-231.
宋文荣,吴仪芳,刘建达,等.1994.电磁波层析成像技术（CT）在地基勘探中的应用[J].地震学刊,2:62-64.
吴俊军.2012.跨孔雷达全波形层析成像方法的研究[D].吉林：吉林大学.
熊金良,李国发,王濮,等.2005.地震波走时计算的逆风差分算法[J].物探化探计算技术,26(4):295-298.
杨文采.1989.地球物理反演和地震层析成像[M].北京：地质出版社.
张东,谢宝莲,杨艳,等.2009.一种改进的线性走时插值射线追踪算法[J].地球物理学报,52(1):200-205.
张旻舳,师学明.2009.电磁波层析成像技术进展[J].工程地球物理学报,6(4):418-425.
张文生,何樵登.1996.有限差分法精确计算地震走时[J].物探化探计算技术,18(4):331-335.
AL-GHAMDI A.2007. Automatic first arrival picking using energy ratios [D]. Saudi Arabia:King Fahd University of Petroleum&Minerals.
ALUMBAUGH D L, NEWMAN G A.2000. Image appraisal for2-D and3-Delectromagnetic inversion[J]. Geophysics,65(5):1455-1467.
AMMON C J, VIDALE J E.1993. Tomography without rays[J]. Bulletin of theSeismological Society of America,83(2):509-528.
ASAKAWA E, KAWANAKA T.1993. Seismic ray tracing using linear traveltimeinterpolation[J]. Geophysical prospecting,41(1):99-111.
BALANIS C A.1989. Advanced engineering electromagnetics[J].
BECHT A, TRONICKE J, APPEL E, ET AL.2004. Inversion strategy in crosshole radartomography using information of data subsets[J]. Geophysics,69(1):222-230.
BELINA F, IRVING J, ERNST J, ET AL.2012. Analysis of an iterative deconvolutionapproach for estimating the source wavelet during waveform inversion of crossholegeoradar data[J]. Journal of Applied Geophysics,78:20-30.
BERES M, HAENI F P.1991. Application of Ground-Penetrating-Radar Methods inHydrogeologie Studies[J]. Ground Water,29(3):375-386.
BREGMAN N D, BAILEY R C, Chapman C H.1989. Crosshole seismic tomography[J].Geophysics,54(2):200-215.
BROYDEN C G.1965. A class of methods for solving nonlinear simultaneous equations[J].Math. Comp,19(92):577-593.
CAI W, QIN F, SCHUSTER G T.1996. Electromagnetic velocity inversion using2-DMaxwell's equations[J]. Geophysics,61(4):1007-1021.
CAO S, GREENHALGH S.1993. Calculation of the seismic first-break time field and its raypath distribution using a minimum traveltime tree algorithm[J]. Geophysical JournalInternational,114(3):593-600.
CARDARELLI E.2002. Ray tracing in elliptical anisotropic media using the lineartraveltime interpolation (LTI) method applied to traveltime seismic tomography[J].Geophysical prospecting,50(1):55-72.
CARDIMONA S J, CLEMENT W P, KADINSKY-CADE K.1998. Seismic reflection andground-penetrating radar imaging of a shallow aquifer[J]. Geophysics,63(4):1310-1317.
CHENG N, HOUSE L.1996. Minimum traveltime calculation in3-D graph theory[J].Geophysics,61(6):1895-1898.
CLEMENT W P, BARRASH W.2006. Crosshole radar tomography in a fluvial aquifer nearBoise, Idaho[J]. Journal of Environmental&Engineering Geophysics,11(3):171-184.
CLEMENT W P, KNOLL M D.2000. Tomographic inversion of crosshole radar dataconfidence in results SAGEEP2000[C]. Environmental and Engineering GeophysicalSociety, Proceedings:553–562.
CONYERS L B.2006. Innovative ground-penetrating radar methods for archaeologicalmapping[J]. Archaeological Prospection,13(2):139-141.
COPPENS F.1985. First arrival picking on common-offset trace collections for automaticestimation of static corrections[J]. Geophysical Prospecting,33(8):1212-1231.
DANIELSSON P E, LIN Q.2003. A modified fast marching method[M]//Image Analysis.Springer Berlin Heidelberg,1154-1161.
DIJKSTRA E W.1959. A note on two problems in connexion with graphs[J]. Numerischemathematik,1(1):269-271.
DINES K A, LYTLE R J.1979. Computerized geophysical tomography[J]. Proceedings ofthe IEEE,67(7):1065-1073.
ENDRES A L, KNIGHT R.1991. The effects of pore‐scale fluid distribution on thephysical properties of partially saturated tight sandstones[J]. Journal of applied physics,69(2):1091-1098.
ENDRES A L, REDMAN J D.1996. Modelling the electrical properties of porous rocks andsoils containing immiscible contaminants[J]. Journal of Environmental and EngineeringGeophysics,1(B):105-112.
FENG S, SEN P N.1985. Geometrical model of conductive and dielectric properties ofpartially saturated rocks[J]. Journal of Applied Physics,58(8):3236-3243.
FISCHER R, LEES J M.1993. Shortest path ray tracing with sparse graphs[J]. Geophysics,58(7):987-996.
GAUTHIER O, VIRIEUX J, TARANTOLA A.1986. Two-dimensional nonlinear inversionof seismic waveforms: Numerical results[J]. Geophysics,51(7):1387-1403.
GIROUX B, BOUCHEDDA A, CHOUTEAU M.2009. Assisted traveltime picking ofcrosshole GPR data[J]. Geophysics,74(4): J35-J48.
GIROUX B, GLOAGUEN E, CHOUTEAU M.2007. bh_tomo—a Matlab borehole georadar2D tomography package[J]. Computers&Geosciences,2007,33(1):126-137.
GLOAGUEN E, CHOUTEAU M, MARCOTTE D, ET AL.2001. Estimation of hydraulicconductivity of an unconfined aquifer using cokriging of GPR and hydrostratigraphicdata[J]. Journal of Applied Geophysics,47(2):135-152.
GLOAGUEN E, MARCOTTE D, GIROUX B, ET AL.2007. Stochastic borehole radarvelocity and attenuation tomographies using cokriging and cosimulation[J]. Journal ofApplied Geophysics,62(2):141-157.
G KTüRKLER G, BALKAYA.2010. Traveltime tomography of crosshole radar datawithout ray tracing[J]. Journal of Applied Geophysics,72(4):213-224.
GREENHALGH S A, BING Z, GREEN A.2006. Solutions, algorithms and inter-relationsfor local minimization search geophysical inversion[J]. Journal of Geophysics andEngineering,3(2):101.
GREMAUD P A, KUSTER C M.2006. Computational study of fast methods for the eikonalequation[J]. SIAM journal on scientific computing,27(6):1803-1816.
HANAFY S, AL HAGREY S A.2005. Ground-penetrating radar tomography forsoil-moisture heterogeneity[J]. Geophysics,71(1): K9-K18.
HANSEN P C, O'LEARY D P.1993. The use of the L-curve in the regularization of discreteill-posed problems[J]. SIAM Journal on Scientific Computing,14(6):1487-1503.
HASSOUNA M S, FARAG A A.2007. Multistencils fast marching methods: A highlyaccurate solution to the eikonal equation on cartesian domains[J]. Pattern Analysis andMachine Intelligence, IEEE Transactions on,29(9):1563-1574.
HINKEL K M, DOOLITTLE J A, BOCKHEIM J G, ET AL.2001. Detection of subsurfacepermafrost features with ground‐penetrating radar, Barrow, Alaska[J]. Permafrost andPeriglacial Processes,12(2):179-190.
HOLE J A, ZELT B C.1995.3-D finite-difference reflection traveltimes[J]. GeophysicalJournal International,121(2):427-434.
HOLLENDER F.1999. Interprétation de la distorsion des signaux géoradar propagés etréflèchis. Développement d'une tomographie par bandes de fréquence[D].
HOLLIGER K, BERGMANN T.2002. Numerical modeling of borehole georadar data[J].Geophysics,67(4):1249-1257.
HOLLIGER K, MUSIL M, MAURER H R.2001. Ray-based amplitude tomography forcrosshole georadar data: A numerical assessment[J]. Journal of Applied Geophysics,47(3):285-298.
IRVING J D, KNOLL M D, KNIGHT R J.2007. Improving crosshole radar velocitytomograms: A new approach to incorporating high-angle traveltime data[J]. Geophysics,72(4): J31-J41.
JACKSON M J, TWEETON D R.1994. MIGRATOM-Geophysical tomography usingwavefront migration and fuzzy constraints[M]. US Department of Interior, Bureau ofMines.
JOHNSON T C, ROUTH P S, KNOLL M D.2005. Fresnel volume georadarattenuation-difference tomography[J]. Geophysical Journal International,162(1):9-24.
JULIAN B R, GUBBINS D.1977. Three-dimensional seismic ray tracing[J]. J. geophys,43(1):95-114.
KIM S, COOK R.1999.3-D traveltime computation using second-order ENO scheme[J].Geophysics,64(6):1867-1876.
KIM S.2001. An O(N) Level Set Method for Eikonal Equations[J]. SIAM journal onscientific computing,22(6):2178-2193.
KIRSCH A.1996. An Introduction to the Mathematical Theory of Inverse Problems [M].Basel: Springer.
KNIGHT R.2001. Ground penetrating radar for environmental applications[J]. AnnualReview of Earth and Planetary Sciences,29(1):229-255.
KURODA S, TAKEUCHI M, KIM H J.2005. Full waveform inversion algorithm forinterpreting cross-borehole GPR data[C]//2005SEG Annual Meeting. Society ofExploration Geophysicists.
LINES L R, TREITEL S.1984. Tutorial: A review of least-squares inversion and itsapplication to geophysical problems[J]. Geophysical prospecting,32(2):159-186.
LIU L, LANE J W, QUAN Y.1998. Radar attenuation tomography using the centroidfrequency downshift method[J]. Journal of Applied Geophysics,40(1):105-116.
LIU S, WU J, DONG H, ET AL.2012. The experimental results and analysis of a boreholeradar prototype[J]. Journal of Geophysics and Engineering,9(2):201.
MEBINGER J.2004. Effective automatic picking of traveltime data with highprecision[C]//Ground Penetrating Radar,2004. GPR2004. Proceedings of the TenthInternational Conference on. IEEE,1:91-94.
MELES G, GREENHALGH S, DER KRUK J, ET AL.2012. Taming the non-linearityproblem in GPR full-waveform inversion for high contrast media[J]. Journal of AppliedGeophysics,78:31-43.
MOGHADDAM M, CHEW W C, ORISTAGLIO M.1991. Comparison of the born iterativemethod and tarantola's method for an electromagnetic time‐domain inverse problem[J].International Journal of Imaging Systems and Technology,3(4):318-333.
MOLYNEUX J B, SCHMITT D R.1999. First-break timing: Arrival onset times by directcorrelation[J]. Geophysics,64(5):1492-1501.
MOSER T J.1991. Shortest path calculation of seismic rays[J]. Geophysics,56(1):59-67.
MOSER T J.1992. The shortest path method for seismic ray tracing in complicated media[J].Geologica Ultraiectina,83:1-180.
MOUSA W A, AL-SHUHAIL A A, AL-LEHYANI A.2011. A new technique for first-arrivalpicking of refracted seismic data based on digital image segmentation[J]. Geophysics,76(5): V79-V89.
MUSIL M, MAURER H R, GREEN A G.2003. Discrete tomography and joint inversion forloosely connected or unconnected physical properties: application to crosshole seismicand georadar data sets[J]. Geophysical Journal International,153(2):389-402.
NAKANISHI I, YAMAGUCHI K.1986. A numerical experiment on nonlinear imagereconstruction from first-arrival times for two-dimensional island arc structure[J].Journal of Physics of the Earth,34(2):195-201.
OLSSON O, FALK L, FORSLUND O, ET AL.1992. Borehole Radar Applied to theCharacterization of Hydraulically Conductive Fracture Zones in Crystalline ROCK1[J].Geophysical Prospecting,40(2):109-142.
PAIGE C C, SAUNDERS M A.1982. LSQR: An algorithm for sparse linear equations andsparse least squares[J]. ACM Transactions on Mathematical Software (TOMS),8(1):43-71.
PODVIN P, LECOMTE I.1991. Finite difference computation of traveltimes in verycontrasted velocity models: a massively parallel approach and its associated tools[J].Geophysical Journal International,105(1):271-284.
PRATT R G.1999. Seismic waveform inversion in the frequency domain, Part1: Theory andverification in a physical scale model[J]. Geophysics,64(3):888-901.
PROTHERO W A, TAYLOR W J, EICKEMEYER J A.1988. A fast, two-point,three-dimensional raytracing algorithm using a simple step search method[J]. Bulletin ofthe Seismological Society of America,78(3):1190-1198.
QIN F, LUO Y, OLSEN K B, ET AL.1992. Finite-difference solution of the eikonal equationalong expanding wavefronts[J]. Geophysics,57(3):478-487.
QUAN Y, HARRIS J M.1997. Seismic attenuation tomography using the frequency shiftmethod[J]. Geophysics,62(3):895-905.
SAAD Y, SCHULTZ M H.1986. GMRES: A generalized minimal residual algorithm forsolving nonsymmetric linear systems[J]. SIAM Journal on scientific and statisticalcomputing,7(3):856-869.
SABBIONE J I, VELIS D.2010. Automatic first-breaks picking: New strategies andalgorithms[J]. Geophysics,75(4): V67-V76.
SCHNEIDER JR W A.1995. Robust and efficient upwind finite-difference traveltimecalculations in three dimensions[J]. Geophysics,60(4):1108-1117.
SEN P N, SCALA C, COHEN M H.1981. A self-similar model for sedimentary rocks withapplication to the dielectric constant of fused glass beads[J]. Geophysics,46(5):781-795.
SETHIAN J A,1999. Leve Sets Methods and Fast Marching Methods, second ed[M].Cambridge Univ. Press.
SETHIAN J A, POPOVICI A M.1999.3-D traveltime computation using the fast marchingmethod[J]. Geophysics,64(2):516-523.
SETHIAN J A.1996. A fast marching level set method for monotonically advancing fronts[J].Proceedings of the National Academy of Sciences,93(4):1591-1595.
SLAWINSKI R A, SLAWINSKI M A.1999. On raytracing in constant velocity-gradientmedia: calculus approach[J]. Canadian journal of exploration geophysics,35:24-27.
SLEIJPEN G L G, FOKKEMA D R.1993. BiCGstab (l) for linear equations involvingunsymmetric matrices with complex spectrum[J]. Electronic Transactions on NumericalAnalysis,1(11):2000.
TARANTOLA A.1984. Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,49(8):1259-1266.
TIKHONOV A N, ARSENIN V I A, JOHN F.1977. Solutions of ill-posed problems[M].Washington, DC: Winston.
TRONICKE J, TWEETON D R, DIETRICH P, ET AL.2001. Improved crosshole radartomography by using direct and reflected arrival times[J]. Journal of AppliedGeophysics,47(2):97-105.
TRONICKE J.2007. The influence of high frequency uncorrelated noise on first-breakarrival times and crosshole traveltime tomography[J]. Journal of Environmental&Engineering Geophysics,12(2):173-184.
TSAI Y H R, CHENG L T, OSHER S, ET AL.2003. Fast sweeping algorithms for a class ofHamilton--Jacobi equations[J]. SIAM journal on numerical analysis,41(2):673-694.
UM J, THURBER C.1987. A fast algorithm for two-point seismic ray tracing[J]. Bulletin ofthe Seismological Society of America,77(3):972-986.
VAN TRIER J, SYMES W W.1991. Upwind finite-difference calculation of traveltimes[J].Geophysics,56(6):812-821.
VIDALE J E.1988. Finite-difference calculation of travel times[J]. Bulletin of theSeismological Society of America,78(6):2062-2076.
VIDALE J E.1990. Finite-difference calculation of traveltimes in three dimensions[J].Geophysics,55(5):521-526.
YATZIV L, BARTESAGHI A, SAPIRO G.2006. O (N) implementation of the fast marchingalgorithm[J]. Journal of computational physics,212(2):393-399.
ZHANG S, ZHOU Z, CHEN L, et al.2007. Seismic ray-tracing method of applyingparabolic Interpolation to travel-time[J]. Progress in Geophysics,1:004.
ZHANG Y T, ZHAO H K, QIAN J.2006. High order fast sweeping methods for staticHamilton–Jacobi equations[J]. Journal of Scientific Computing,29(1):25-56.
ZHAO H.2005. A fast sweeping method for eikonal equations[J]. Mathematics ofcomputation,74(250):603-627.
ZHOU B, FULLAGAR P K.2000. Borehole radar conductivity and permittivitytomography[C]//8th International Conference on Ground Penetrating Radar.International Society for Optics and Photonics,542-547.
ZHOU B, GREENHALGH S A.1998. Crosshole acoustic velocity imaging withfull-waveform spectral data:2.5-D numerical simulations[J]. Exploration Geophysics,29(3/4):680-684.
ZHOU C, LIU L, LANE JR J W.2001. Nonlinear inversion of borehole-radar tomographydata to reconstruct velocity and attenuation distribution in earth materials[J]. Journal ofApplied Geophysics,47(3):271-284.
ZHOU Z, ZHANG S, CHEN L.2004. Seismic ray-tracing calculation based on parabolictravel-time interpolation[J]. Journal of Central South University of Technology,11(2):199-205.
曾昭发,刘四新,王者江,等.2006.探地雷达方法原理及应用[M].北京：科学出版社.
成谷,张宝金,马在田.2007.浅谈反射地震走时层析中的正则化[J].地球物理学进展,22(6):1715-1721.
李灿苹,王南萍,李志宏.2010. VonKarman型自相关函数模拟随机介质[J].物探与化探,34(1):98-102.
李玉喜.2010.钻孔雷达层析成像软件系统的研究与开发[D].吉林：吉林大学.
佩特罗夫斯基.1981.地下无线电波法[M].北京:地质出版社:109-231.
宋文荣,吴仪芳,刘建达,等.1994.电磁波层析成像技术（CT）在地基勘探中的应用[J].地震学刊,2:62-64.
吴俊军.2012.跨孔雷达全波形层析成像方法的研究[D].吉林：吉林大学.
熊金良,李国发,王濮,等.2005.地震波走时计算的逆风差分算法[J].物探化探计算技术,26(4):295-298.
杨文采.1989.地球物理反演和地震层析成像[M].北京：地质出版社.
张东,谢宝莲,杨艳,等.2009.一种改进的线性走时插值射线追踪算法[J].地球物理学报,52(1):200-205.
张旻舳,师学明.2009.电磁波层析成像技术进展[J].工程地球物理学报,6(4):418-425.
张文生,何樵登.1996.有限差分法精确计算地震走时[J].物探化探计算技术,18(4):331-335.
© 2004-2018 中国地质图书馆版权所有 京ICP备05064691号 京公网安备11010802017129号 地址：北京市海淀区学院路29号 邮编：100083 电话：办公室：(+86 10)66554848；文献借阅、咨询服务、科技查新：66554700 |