随着国民经济建设的不断加速发展，我国对铜、金、银、铀等金属矿产的需求不断增加.探地雷达通过处理和分析接收到的电磁波反射信号，根据波形、强度和双程走时等参数探测地下隐蔽目标物(杨峰，2010)，是一种高分辨率、高效率的无损探测技术，已广泛应用到地质勘察工程中(杨峰，2010；李华等，2010).金属矿区的地质结构比油气等能源更加复杂，除少量的层状金属矿床外，大多金属矿地层倾角较大，与隐伏岩体有关(程久龙等，2000；王妙月，2010).目前应用于金属矿区探测的地球物理手段主要包括重力勘探、磁力勘探、电法和探地雷达等(刘光鼎和郝天珧，1995；严加永等，2008；徐耀鉴等，2009).其中，探地雷达对浅层地段不仅有很高的探测精度和高分辨率，能达到对浅部采空区的精确定位，但其反射波勘探方法的横向分辨率较低，难以解决金属矿勘查中的陡倾斜地层、复杂构造和不规则形态的隐伏岩体等地质问题.因此，迫切需要寻找新的方法和技术，提高探地雷达的分辨率.

从上世纪90年代开始，国外对探地雷达层析成像技术(Computerized Tomography，CT)进行了广泛研究和实践应用，如地层结构探测(Knight，2001)，桥梁溶洞的探测(Tan et al., 2012)，污染物示踪物的探测(Clement and Barrash， 2006；Chang and Alumbaugh， 2011)，矿体圈定(Bellefleur et al., 2000)等.有学者通过探地雷达层析成像技术同其反射法勘探(Gr and jean et al., 2000)与电阻率层析成像技术的对比研究(Pellicer and Gibson， 2011；Haarder et al., 2012)，证实了该技术在探测复杂地质构造方面的优越性.由于国内用于层析成像探测的雷达设备都处于研制阶段(Liu et al., 2012)，只能依赖于进口国外的钻孔雷达，因此相关的实践应用较少.有学者利用采集到的信号中的所有信息进行全波形反演(Cai et al., 1996；Alumbaugh and Newman， 2000；吴俊军，2012)，但由于这种技术需要大量的计算资源，目前应用最广泛的还是基于射线的层析成像技术.

根据反演采用的数据类型不同，层析成像分为波速层析和衰减层析两类.当探测介质是低损耗介质时，只利用采集到的电磁波走时数据进行波速反演就可以计算介质的相对介电常数(杨峰，2010).由于绝大多数岩土介质均属于低损耗介质(杨峰，2010)，本文将波速层析用于探地雷达对金属矿区的勘查，通过建立速度模型，对模型数据进行正反演计算，并从模型整体和具体切片两个角度将计算结果与模型数据进行对比，分析该方法的有效性.该项技术可提高探地雷达的横向分辨率，实现对金属矿区倾角较大的地层和复杂的地质构造及隐伏岩体等的高精度探测. 1 雷达波速层析的原理和方法

1.1 岩土的电磁特性

探地雷达采用高频电磁波进行探测，根据电磁波传播理论，高频电磁波在介质中的传播满足麦克斯韦方程组.在低损耗和非色散介质中，雷达电磁波波速可用如下方程(Balanis，1989)描述：

相对介电常数是反映地下介质电性的一个重要参数.各类岩石、土的电磁学性质有了很多的研究和测定.空气是自然界中介电常数最小的介质，电磁波速最高，而水介电常数最大，电磁波速最小.干燥的岩石和土的电磁参数差异不大，但由于各类岩土的孔隙率和含水饱和程度不同，相对介电常数差异较大，使得不同岩性对应不同的波速.

1.2 层析探测方法

在进行探地雷达波速层析测量时，在待测剖面两侧布置均匀、对称的测点，将雷达的发射天线与接收天线分离，分别置于待测剖面的两侧.将发射天线固定在一个测点，接收天线在另一侧沿各测点移动，采集相应的透射波的走时数据，将发射天线置于各个测点，均重复以上过程，由此得到整个剖面所有发射点与接收点的组合的走时数据.

2 正反演算法

层析成像技术分为波速层析和衰减层析两类.其中，波速层析是利用接收天线接收波形的初至时间来反演计算对探测剖面介质的波速数据.首先将探测剖面划分为均匀的网格，每个网格中的介质的波速视为一个均值，将连续分布的波速离散化.介质的速度差异小于15%时，可将电磁波近似视为直线传播，则每一道射线的走时即为射线路径中穿过的网格的慢度线性积分，即这些网格中射线的长度与对应网格的慢度值之积的求和，正演方法如图 1所示.

图 1

Fig. 1

图 1 雷达层析正演方法 Fig. 1 GPR tomography forward method

基于上述分析，层析成像数据处理的实质是解一个大型的矩阵方程组

由于射线分布不均，且每条射线只经过整个网格剖面的少数网格，L 矩阵通常是稀疏的.该反演方程组一般是病态的，无法用直接法求解，在诸多迭代算法中，最小平方正交二乘算法(Least Square QR，LSQR)可利用矩阵的稀疏性简化计算(Paige and Saunders， 1982; Sanchez-Avila and Garcia-Moreno， 1990)，具有计算量小、收敛速度快的优点，在地球物理反演问题中得到广泛应用(Bunse-Gerstner et al.，2006; Jiang et al.，2009a; Jiang et al.，2009b).

3 建模与仿真

如图 2a所示，假定待测剖面的空间范围为200 m×100 m，在剖面两侧分别布置了40个发射点和40个接收点.在使用探地雷达进行波速层析测量时，一共可采集到1600道射线的走时数据.

图 2

Fig. 2

图 2 金属矿区速度模型 (a)布点示意图；(b)模型一：断层模型；(c)模型二：复杂构造模型. Fig. 2 Velocity models of metallic mine

针对金属矿区的地质特点，设计了两个具有不同结构的模型.考虑到随着深度的增加，地下介质的含水量增大，模型 中介质随深度的增加电磁波波速降低.如图 2b所示，模型一包含断层界面，三层介质的 波速分别为0.12 m/ns，0.11 m/ns，0.105 m/ns.如图 2c所示，模型二为包含陡峭、倾斜界面和多种岩石介质的金属矿区复杂速度模型.为了观察波速层析技术对均匀介质、纵向变化、横向变化的分辨率，在模型一、二中分别选取了两个切片进行观察，如图 2b和图 2c中直线所示.

将剖面划分为40行×20列的网格模型，网格为边长5 m的正方体，得到模型的速度分布 V₁.根据模型中的速度分布计算走时，将计算结果作为实测走时T .介质的速度差异小于15%时，可将电磁波近似视为直线传播，计算出每条射线在每个网格中的长度，得到系数矩阵 L，由此构建出公式(2)给出的反演方程组.采用Lsqr算法求解方程组，迭代初始值为零向量，解得模型中的慢度分布 S，代入公式有

即求得反演计算得到的模型速度分布V₂.误差值的计算方法为

4.1 模型整体结果分析

图 3a为模型一反演结果的等值线图，成像结果很好的反演出界面的起伏形态，可以辨认出断层的断面.图 3b为模型二反演结果的等值线图，模型的速度结构基本得到恢复，各个不同速度区域的轮廓都比较清晰.

图 3

Fig. 3

图 3 反演结果 Fig. 3 Inversion result

网格编号规则为以网格的左下角为起点，自下而上、自左向右逐行进行编号，即左下角为1号，右上角为800号.从图 4可以看出，模型一的误差分布范围为[-0.02%，0.02%]，并且大部分集中在[-0.01%，0.01%]的范围内.模型二的误差分布范围约为[-0.1%，0.1%]，误差分不均匀，1~500号网格的误差离散程度大，分布比较松散，而500~800号网格大部分集中在[-0.05%，0.05%]的范围内.从绝对的角度看，模型一和模型二的反演结果误差都是1‰的数量级，结果的准确度很高，与模型的速度分布相吻合.从相对的角度看，模型二比模型一的速度结构更复杂，而误差也更严重，并且模型二的速度不均匀分布也体现了速度结构复杂性对反演结果的干扰.

图 4

Fig. 4

图 4 反演结果误差分布 Fig. 4 Error analysis of inversion result

基于对模型整体的结果分析，观察四个切片的反演结果，其中切片B的20个波速值为从40个网格中抽取偶数号组成.对比反演值与模型的波速值(图 5)，模型一中的切片A与切片B的反演结果与模型值比较吻合.模型二中的切片C与切片D的反演结果与模型值的偏离程度较为严重，但整体的升降变化趋势与模型值相同.

图 5

Fig. 5

图 5 切片反演结果 Fig. 5 Inversion result of sections

如图 6所示，进一步分析四个切片的误差分布.模型一中的切片A与切片误差很小，与真值吻合程度高.模型二中的切片C与切片D变化趋势相同，切片D的误差更大.在模型值发生突变的位置，即不同速度值的区域交界处，正误差与负误差发生了转换.切片D中的速度差异性大于切片C，其引起的结果偏离也更明显.

图 6

Fig. 6

图 6 切片误差对比 Fig. 6 Error comparison of sections

本文提出将探地雷达与层析成像技术相结合，用于金属矿区复杂地质结构的勘探，运用LSQR算法迭代求解反演方程组.对理论模型的处理达到了预期的效果，反演数据与模型数据相吻合，证明探地雷达层析成像技术用于金属矿区勘探是可行的.相对于其他物探技术包括探地雷达反射波探测技术，它具有探测速度快、分辨率高等优点，为金属矿区勘探技术的研究提供了新的思路.

探地雷达层析成像技术的纵向分辨率很高，并且获得了较高的横向分辨率，但反演误差随速度结构复杂程度而增大，尤其是陡峭倾斜界面处对反演结果的干扰较大.今后工作可以改进反演算法，使该技术得到进一步完善.