基于地质雷达和钻孔数据的三维地层建模

第30卷增刊岩 土 力 学 V ol.30 Supp. 2009年8月 Rock and Soil Mechanics Aug. 2009收稿日期：2009-03-20作者简介：朱发华，男，1983年生，博士生，主要从事岩土工程信息化和可视化技术研究。

文章编号：1000－7598－(2009) 增1―0267―04基于地质雷达和钻孔数据的三维地层建模朱发华1,2，贺怀建1（1. 中国科学院武汉岩土力学研究所，武汉，430071； 2. 中国科学院研究生院，北京，100049）摘 要：为了结合地质雷达探测图像和钻孔数据共同用于地层的三维建模，根据反射波组的同相性和反射波形的相似性，通过与钻孔采样数据的对比，从处理后的地质雷达数据中提取出虚拟钻孔数据用于建模；基于该类数据的一般分布特点，根据曲面光顺的条件构建了网格的插值算法。

采用该方法对某区域进行了三维地层建模。

关 键 词：地质雷达；网格插值；地层建模；三维可视化 中图分类号：P22 文献标识码：A3D stratum modeling based on ground penetrating radar and borehole dataZHU Fa-hua 1,2, HEHuai-jian 1 (1. Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China;2. Graduate School, Chinese Academy of Sciences, Beijing 100049, China)Abstract: In order to model 3D strata using comprehensive data of ground penetrating radar (GPR) graph and borehole, virtual borehole data for modeling are extracting from GPR data by comparing with borehole sampling data according to in-phase of reflection wave group and reflection wave shape analogy; the grid interpolation method based on distribution characteristics of this kind data is given according to surface fairing condition. By using this method, strata of one area is modeling. Key words: GPR; grid interpolation; stratum modeling; 3D visualization1 引 言钻孔数据是三维地层建模的主要数据来源。

通过钻孔采样可以直接获取采样点详细的地层分布状况，但钻孔数据获取成本较高且区域采样数目有限，仅仅依靠钻孔数据是不能准确地构建地层模型，尤其是复杂地层的模型。

地质雷达作为对地下目标有效探测的一种手段，有探测速度快、分辨率高、探测范围广等诸多优点。

利用钻探和地质雷达探测两种方法对地层进行数据采集，可以更加有效地获取和反映地层的情况。

在工程检测和勘查中，就常综合应用地质雷达、钻探和地震勘查等方法探测目标的空间位置、空间形态、埋藏深度和特征。

地质雷达的数据也会被处理渲染为二维或三维的图像来展示探测的成果。

本文借鉴了地质雷达在实际工程中应用和数据处理表现的方法，参照钻孔数据，从处理后的地质雷达数据中提取虚拟钻孔数据，并应用所建立的网格插值算法来构建地层界面格网，从而生成更加准确的地层模型。

2 地质雷达探测的基本原理地质雷达探测是近几十年来发展最快、应用最广的一种地球物理探测方法[1]。

其基本原理是基于电磁波的反射理论，利用发射天线（T ）将宽带短脉冲形式的高频电磁波送入地下，在不均匀介质（主要是介电常数不连续）的交界面上发生反射后返回地面，由接收天线（R ）所接收；将接收到的反射信号加以处理分析，就可得到地表下地质体的特征信息。

地质雷达探测的基本原理如图1所示。

一般在实际工程应用中，发射与接收天线的距离很近（≤2 m ），可近似看作法向入射，故可以由反射时间来近似估算地下目标的深度：/2h tv = （1）岩 土 力 学 2009年（a ）反射原理 （b ）理想GPR 记录图1 地质雷达探测原理示意图 Fig.1 Sketch of detechtion principle of GPR式中：t 为地下目标的雷达波的反射时间（即双程时间）；v 为目标以上介质波速的平均值。

原始的地质雷达数据中存在一些干扰信号，并且所采集的数据也存在不同程度的空间变形，所以在分析解释地质雷达数据之前要对其进行处理。

雷达资料的常规处理方法与地震勘探数据的处理方法基本相同[1–3]。

此外，小波分析也应用到了地质雷达资料的处理中[4–5]。

处理以后的数据图像就能比较真实地反映地下介质的情况，从而就可根据图形的变化对地下的地质状况作出解释了[6–7]。

3 地质雷达数据中地层分界面的提取地质雷达探测数据是电磁波辐射范围内所有地质体散射、反射信号的综合，反映了地下介质的电性分布情况。

不同的地层一般是存在电性差异的，在地层的交界面上会发生反射现象，从而在地质雷达图像上得到反映。

要从地质雷达数据中提取所需要的地层分界面信息，就是要在雷达数据图像中将相应的反射波组找到。

同一地层界面的反射波组的波形、振幅、周期以及其包络线形态等有一定特征，在相邻的记录道上这些特征会保持不变，即所谓的反射波形的相似性。

将位于雷达测线上的钻孔取样与该测点的波形数据相互参照，从波形数据中取得与地层界面对应的反射波特征。

根据所取得的反射波特征，找出所有记录道上具有该特征的一系列反射波组，连接其同相位处构成一与地层界面对应的反射层的同相轴，这样就基本上将地层分界面提取出来了，但要获得地层面的深度的具体数值，还涉及到地层电磁波速的计算以及时深的换算。

由于地质雷达探测所获得的电磁波双程旅行时间是相当准确的，在对原始的雷达数据进行处理之后，雷达图像上波的反射时间与地下介质深度的关系就更加明确，可以由公式/2h tv =直接来换算。

所以，对目标深度定位的问题最终转化为对电磁波在地下传播速度的确定。

在高阻介质中，电磁波的传播速度可以表示为[1]v c = （2）式中：c 为真空中电磁波的传播速度，8310m/s c =×；ε为介质的介电常数（在地层深度的测定中，ε应为反射界面以上的介质的介电常数平均值，可以由钻孔取样来测定）。

4 地层分界面格网的插值与构建在从地质雷达数据中提取了地层分界面（一系列地层剖面上的分界线）后，为了将其应用于分界面格网的构建，需要将它们按一定的取样率转化为虚拟钻孔。

由于这种基于钻孔和地质雷达对地层探测的方法所获得的数据分布不均匀——数据点都分布在初始网格线上（如图2所示），所以有必要设计出一种加密格网的插值方法来适应该采样数据的特点，从而使生成的新格网更贴合实际分界面的情 况。

图2 初始网格及数据点示意图Fig.2 Sketch of initial grids and data points268增刊 朱发华等：基于地质雷达和钻孔数据的三维地层建模本文采用规则格网建模。

为了使生成的格网所渲染出的面更加平滑，基于曲面光顺的条件[8]——二阶导数平方和最小，由曲面上纵横方向曲线的平均斜率改变连续线性变化为条件获得网格点初始值，再迭代求取网格点终值。

具体的插值算法如下（见图3）：图3 插值网格示意图Fig.3 Sketch of interpolation grids（1）根据初始网格纵横边界上数据点（钻孔和虚拟钻孔点）的个数确定新网格。

以图3为例，初始网格横向边界的数据点有m +1个，纵向边界的数据点有n +1个，可以将初始网格划分为mn 个新网格，相应地在初始网格内部要增加(m －1)(n －1)个插值点。

（2）曲面上插值点在纵向和横向曲线上的微分用格网中相应的插值点与其相邻点的差分来近似。

首先计算插值点的初始值，令该点纵横方向上的二阶差分平均值等于原始网格纵横边界相应数据点的二阶差分的反距离加权和。

对于图中的P ij 点，有：l l clc 0022c im i nj j ijijj m j i n i m m n n −−∆+∆+∆+∆∆+∆=（3）式中：l ij ∆，cij ∆分别为格网中P ij 点的纵横二阶差分，可分别表示为l (1)(1)c(1)(1)2;2ij i j i j ij ij i j i j ij z z z z z z +−+−⎫∆=+−⎪⎬∆=+−⎪⎭（4） 其他类似的符号可做相似的解释。

对格网中的每一个插值点都列出这样一个关系式，将这些关系式联立组成方程组求解便得到所有插值点的初始值。

（3）用格网中所有插值点的纵横二阶差分的平方和最小来取代插值点二阶导数平方和最小，作为曲面光顺的条件，即()()22lc ,min []ij ij i j∆+∆∑。

根据极值条件，可以求取出一系列方程，它们表示了插值点与围绕它的相邻点的数值关系——本文根据插值点所处的不同位置给出了如图4所示的3种系数模板图，其中黑色点代表了插值点，空心点是与插值点数值相关的相邻点。

插值点的数值就等于所有相邻点的数值与其上的系数乘积之和。

将(2)计算出来的初始值代入到这一系列方程中进行迭代即可求出最终的网格插值点的值。

图4 系数模板Fig.4 Coefficient stencits对所有的地层界面，用上述插值算法计算出网格插值点的数值以构建出地层界面格网，然后将格网进行渲染以及其他一些几何图形处理便可生成地层的多层DEM 模型[9]。

5 实 例对某200 m×150 m 区域进行地层三维建模，该区域上的钻孔和地质雷达测线分布如图5所示，区内有10个钻孔点、4条雷达测线，均位于边界上。

对钻孔和地质雷达数据分析对比后，在每条纵横边界上各提取出9个和8个虚拟钻孔点，从而将该区域以雷达测线为边界的初始格网进行插值，总共分11行9列进行插值。

图5 区域钻孔和地质雷达测线分布Fig.5 Distributions of boreholes and GPR surveying linesin the area本文根据上述方法，在Windows 环境下采用基于C++语言的可视化工具包VTK [10–11]进行地层的三维建模，建模的效果如图6所示。