【作者】 王飞；

【导师】 刘四新；

【作者基本信息】 吉林大学 ， 地球探测与信息技术， 2014， 博士

【摘要】 电磁波层析成像技术（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算法的反演结果更好一些。减小反演网格距可以提高反演精度，但是反演精度的提高是以牺牲计算时间为代价的。但是，在基于二维弯曲射线追踪的走时层析成像算法中，这个计算时间的提高还是可以接受的。因此，在跨孔雷达走时层析成像中，我们可以通过减小反演网格距的尺寸来提高反演的精度。通过减小天线移动步长能获得更多的射线数（观测数据），从而能提高反演的精度；但是同时也会引入更多的大角度射线，从而降低了模型重建的效果。当我们减小天线移动步长并且只使用小角度射线进行反演时，我们即提高了反演精度，又能获得更好的反演结果。更多还原

【Abstract】 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.更多还原