探地雷达宽角反射图形拟合方法

doi:10.11720/wtyht.2021.2506

探地雷达宽角反射图形拟合方法

蔡连初¹, 缪念有²

1. 浙江华东建设工程有限公司,浙江杭州 310014

2. 华东勘测设计研究院,浙江杭州 310014

The graph fitting method of wide angle reflection of GPR

CAI Lian-Chu¹, MIAO Nian-You²

1. Zhejiang Huadong Construction Engineering Co., Ltd.,Hangzhou 310014, China

2. Huadong Engineering Co., Ltd., Hangzhou 310014, China

责任编辑: 沈效群

收稿日期: 2019-10-30 修回日期: 2020-08-29 网络出版日期: 2021-02-20

Received: 2019-10-30 Revised: 2020-08-29 Online: 2021-02-20

作者简介 About authors

蔡连初(1969-),男,高级工程师,主要从事工程物探及地质超前预报工作。

摘要

探地雷达共中心点(CMP)法仅用于局部速度测试,为使用宽角反射法进行剖面探测,可固定一个雷达天线而将另一个天线沿测线移动,从而获取整个地质界面的反射信号。本文提出宽角反射剖面探测方法,基于平直岩层界面导出其反射信号的关系式,并开发出解释软件,可精确计算岩层几何参数和层速度,从而进行地质分层。

关键词： 探地雷达 ; 宽角反射 ; 图形拟合 ; 时距曲线 ; 层速度

Abstract

The Geological Radar Common Middle Point (CMP) method is only used for local velocity testing. For the purpose of using a wide-angle reflection method for profiling, one radar antenna can be fixed and another antenna moves along the line of measurement to obtain the reflection signal of the entire geological interface. In this paper, the method of wide angle reflection profile detection is proposed. The relations of reflection signals are derived based on the interface of straight rock formations, and the interpretation software is developed to accurately calculate the geometric parameters and layer velocity of rock formations, so as to carry out geological stratification.

Keywords： GPR ; wide angle reflection ; graph fitting ; time-distance curve ; layer velocity

PDF (2348KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

蔡连初, 缪念有. 探地雷达宽角反射图形拟合方法. 物探与化探[J], 2021, 45(1): 239-244 doi:10.11720/wtyht.2021.2506

CAI Lian-Chu, MIAO Nian-You. The graph fitting method of wide angle reflection of GPR. Geophysical and Geochemical Exploration[J], 2021, 45(1): 239-244 doi:10.11720/wtyht.2021.2506

0 引言

探地雷达探测技术的应用非常普遍,但仍存在诸多硬件和软件上的技术问题亟须解决,例如提高探测深度与精度,研制集束、宽频带、高发射率、体积小重量轻的发射天线,提高发射机发射能量,消除噪声的数据处理方法等,在野外数据采集方法、测线布置方式及反演计算等方面也应进行深入的研究。探地雷达的野外工作必须根据所要研究的地质、岩土工程问题,采取合适的观测方式,正确选择测量参数,才能得到高质量的、满足解释要求的原始记录。

野外观测一般使用收发一体天线或将发射天线和接收天线靠近同时沿测线移动,反射界面的反射波沿入射波原路返回,称为窄角反射法,如此获取整条测线的连续振幅—时间数据序列,可以接收到最强的反射信号;如果将收发天线分离,使反射界面的反射波和入射波呈一定夹角,称为宽角反射法^[1]。相较于窄角反射法,宽角反射的信号强度降低,但当可观测连续反射波同相轴均来自同一个平直界面时,可以计算测线和反射界面之间的电磁波速度和介电常数,并借此进行地质分层。

1 测线布置方法

宽角反射法测线布置如图1所示,固定发射天线(或接收天线)于测线上一点,另一个天线沿测线均匀移动,即可得到宽角反射雷达图像。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 宽角反射原理

Fig.1 Schematic diagram of wide angle reflection

典型宽角反射雷达图像如图2所示,图中信号主要包括头部的空气直达波、岩体直达波和岩体反射波,当在地下洞室中测试时,还能采集到信号很强的空气反射波,甚至形成多次反射。在空气中传播的直达波和反射波速度明显高于岩体界面反射波。因此根据雷达反射信息准确计算雷达图像中各种同相轴的传播速度,不仅是地质分层的需要,也可以借此将岩体界面反射波和其他干扰波区分开来。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 典型的宽角反射图像

Fig.2 Typical wide angle reflection image

2 倾斜界面正演

如图1所示,根据反射原理,倾斜界面上的反射波按照法向镜像路径回到测线,建立射线路径总长度与天线收发间距的函数关系,则反射波的时距曲线方程为:

(1)

v^{2} t^{2} - {(x - 2 h \sin α)}^{2} = 4 h^{2} {(\cos α)}^{2}

式(1)为双曲线方程,式中:v为介质电磁波速度,α为倾斜界面与测线夹角,h为固定天线与倾斜界面的距离;如果给定此3个初始条件,就可得到宽角反射正演曲线。建立三层地质模型(表1),代入各项参数,按式(1)计算,正演时距曲线如图3所示。

表1 三层宽角反射界面模型

Table 1 Three layer wide angle reflection interface model

	v/(m·ns^-1)	α/(°)	h/m
第一层	0.1	0	6
第二层	0.09	20	10
第三层	0.11	-30	20

注:v为测线至界面之间的平均速度或叠加速度。

新窗口打开| 下载CSV

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 三层模型正演时距曲线

Fig.3 Time-distance curve of forward modeling of three-layer model

3 反演计算

从雷达图像中拾取同相轴读数,根据式(1)进行数值拟合,可计算其参数,但读数及计算工作量大^[2]。本文从正演图形入手,对照实测雷达同相轴,采取直观图形对比方法,实时修正正演曲线,使所绘曲线与雷达图像同相轴吻合达到最佳,可得到式(1)中参数的最优解。据此思路开发了图形拟合软件,主要步骤如下。

1)雷达图像预处理。原始雷达数据受地质、地形、现场干扰等多种因素的制约,形成雷达图像之前,应对数据进行预处理以达到最利于进行地质解释的目的。数据处理手段众多,宽角反射法得到的双曲线形态的同相轴反演时因干扰造成的误差会被放大,因此应更加精细处理数据^[3]。数据处理主要包括以下几方面。

a.滤波,主要除去与优势频率相差较大的频率成分,突出目的体的有用信息,降低背景噪声和余振影响,提高图像品质^[4,5,6]。

b.调整增益,由于有耗介质对电磁波的吸收作用,随着探测深度的增加,信号振幅会很快衰减,另外由于宽角反射法随着入射角的增大,反射信号振幅也会减小,因此应调整增益,使深部的信号达到可明显分辨的状态。调整增益时一般采用指数增益方法,增益曲线应平滑连续,不可突变,增益强度一般以背景噪声尽量弱,有用信号明显可分辨为原则。

c.对雷达图像进行x轴插值校正,用以修正因雷达天线不均匀移动而造成的x轴方向伸缩变形,校正的关键在于采集数据时必须在天线经过某些特征点或按一定间距做上适当的标记。

d.如果测线经过地面或地下洞室的开挖表面不足够平直,则需进行地形测量工作,再根据地形对雷达图像做地形校正。

因解释软件不对雷达数据直接成图,需对雷达图像截图保存为图形文件后再供解释软件使用,并记录截取的测线长度和时窗长度。需要说明的是图形拟合本身是对结构面产状的反演计算,因此不需要作其他偏移归位处理^[7]。

2)计算双曲线参数。在雷达图像上拾取双曲线同相轴顶点(x₁,t₁),其顶点坐标为(2hsin α,2hcos α)。在同相轴上拾取另一点(x₂,t₂),联立方程求解:

(2)

\{\begin{array}{l} x_{1} = 2 h \sin α, \\ t_{1} = 2 h \cos α / v, \\ v^{2} t_{2}^{2} - (x_{2} {- 2 h \sin α)}^{2} = 4 h^{2} \cos α^{2} \end{array}

解得曲线方程的参数为:

(3)

\{\begin{array}{l} α = ta n^{- 1} (\frac{x_{1} \cdot \sqrt[]{t_{2}^{2} - t_{1}^{2}}}{(x_{2} - x_{1}) \cdot t_{1}}), \\ h = \sqrt[]{x_{1}^{2} + (x_{2} - x_{1})^{2} \cdot \frac{t_{1}^{2}}{t_{2}^{2} - t_{1}^{2}}} / 2, \\ v = \frac{x_{2} - x_{1}}{\sqrt[]{t_{2}^{2} - t_{1}^{2}}} \end{array}

3)重画时距曲线。根据拾取的两点和计算得到的3个参数即可构建双曲线方程,并在雷达图像的对应图层中重画如图3的正演时距曲线。

4)对比重画的曲线和实测雷达同相轴,修正曲线顶点和曲率以达到最佳拟合效果。

对于矿物成分非磁性、低电导率类岩体,在高频天线工作模式下,符合低耗介质极限情况,电磁波速度仅与相对介电常数相关,而水的相对介电常数比一般岩石的大得多,因此岩体中的含水量对电磁波速度的影响非常大。例如:在图2中,因岩体表面含水量低于深部,岩体直达波的速度比岩体界面反射波速度可提高约10%;在地下洞室中,常温下当空气湿度达到饱和时,空气中的直达波和反射波速度可降低到真空电磁波速度的75%。可见通过速度测定方法探测地下水是行之有效的,对岩体破碎带的分析也应首先考虑含水量的变化^[8]。

反演得到的速度为测线与对应反射界面之间的加权平均速度或叠加速度,仍需计算各界面之间的层速度才能用于地质解释。通过图形拟合方法得到的叠加速度可视为均方根速度v_R。

一般利用Dix公式计算层速度^[9]:

(4)

v_{i}^{2} = \frac{v_{R_{i}}^{2} \cdot t_{0 i} - v_{R_{i - 1}}^{2} \cdot t_{0 i - 1}}{t_{0 i} - t_{0 i - 1}}

式中:t₀为固定天线至反射界面垂直点的单程走时,得到的层速度无需做倾角校正。

4 工程实例

某水电站3#施工支洞岩性为泥质白云岩和硅质白云岩,断层发育,沿断层和裂隙有强烈溶蚀现象,多见闪长玢岩岩脉,围岩类别Ⅲ~Ⅳ类。开挖期间开展了探地雷达超前预报工作,桩号3支0+333~343 m段左侧边墙的雷达图像如图4所示,采用收发一体天线的窄角反射方法,从图中可明显分辨出4条信号强烈的同相轴,进行偏移归位处理后其解释成果为4条结构面,其中J1、J2、J3相互平行,与测线的夹角为-32.1°,J4与测线夹角为17.2°。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 窄角反射雷达图像及解释成果

Fig.4 Narrow angle reflection radar image and interpretation results

为进行对比分析,在3#施工支洞相同部位进行了宽角反射探地雷达试验。在3支0+336放置一个固定天线,测线为3支0+325~351,雷达图像反演计算如图5所示。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 宽角反射反演

Fig.5 Wide angle reflection inversion

图5中可分辨出4支较明显的双曲线同相轴对应图4中的4条结构面,分别计算其各项参数,并和窄角反射法所测得的成果进行对照,得到解释成果见表2。

表2 宽角反射与窄角反射成果对照

Table 2 Comparison of wide angle reflection and narrow angle reflection results

结构面	测试方法	α/(°)	h/m	v_R/(m·ns^-1)
J1	窄角	-32.1	2.2	/
	宽角	-30.5	2.7	0.103
J2	窄角	-32.1	5.6	/
	宽角	-30.9	5.8	0.098
J3	窄角	-32.1	7.8	/
	宽角	-34.8	7.5	0.093
J4	窄角	17.2	11.2	/
	宽角	18.5	10.2	0.095

新窗口打开| 下载CSV

将反演结果代入式(4),分别计算各层电磁波速度v_i和相对介电常数ε_r_i,结果见表3。

表3 层速度计算结果

Table 3 Interval velocity calculation table

i	v_R/(m·ns^-1)	t₀_i/ns	V_i/(m·ns^-1)	ε_r_i
1	0.103	24.2	0.103	8.5
2	0.098	51.9	0.093	10.3
3	0.093	73.4	0.08	14.1
4	0.095	103.4	0.1	9.0

新窗口打开| 下载CSV

地质解释成果见图6。图中结构面J1~J4之间速度分层较明显,第2、3层速度低于第1、4层的速度,推测JI~J3之间大致平行发育3条结构面,岩体相对破碎,孔隙度增大,致使含水量增加或夹有泥质,因而介电常数变大,速度降低^[8]。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 宽角反射成果

Fig.6 Results of wide angle reflection

窄角反射和宽角反射的雷达测线均布置在4条结构面出露范围之外,按照反演结果将4条结构面按平直方向延伸,J1、J2、J3、J4与隧洞相交部位分别为3支0+331、3支0+325、3支0+323、3支0+369。经开挖验证,4条结构面均在隧洞壁出露,且相交部位误差不超过2 m,J1、J2、J3为3条大致平行的溶蚀裂隙,见泥质填充物。

5 图形拟合与数值拟合结果对比

和数值拟合方法相比,图形拟合方法不需要从雷达图像逐点读取数值,大大减少了计算量,且可以实时将实测曲线和正演曲线进行对比,即时修正,自动剔除了突变干扰点,最大限度减少读数误差,因而更加准确、快速,且操作简单。但对于拟合结果的准确性必须以实例计算结果进行对比。

将式(1)移项得到:

(5)

t^{2} = \frac{1}{v^{2}} \cdot x^{2} - \frac{4 h \cdot \sin α}{v^{2}} \cdot x + \frac{4 h^{2}}{v^{2}},

令:

$τ = t^{2}, a_{0} = \frac{4 h^{2}}{v^{2}}, a_{1} = - \frac{4 h \cdot \sin α}{v^{2}}, a_{2} = \frac{1}{v^{2}};$

将时距曲线构造为(x、τ)的标准二次多项式:

(6)

τ = a_{0} + a_{1} \cdot x + a_{2} \cdot x^{2}

对式(6)可采用最小二乘法进行数值拟合计算时距曲线的参数。以上述工程实例中结构面J4为例进行数值拟合对比,等距选取J4同相轴上10组读数(表4),计算出正规方程组中的各项数值(表5),由此得到正规方程组:

$[\begin{array}{l} 10 & 10 & 340 \\ 10 & 340 & 1000 \\ 340 & 1000 & 21328 \end{array}] [\begin{array}{l} a_{0} \\ a_{1} \\ a_{2} \end{array}] = [\begin{array}{l} 509017 \\ 54152 \\ 17463105 \end{array}],$

解得:a₀=48 241.22,a₁=-1 630.81,a₂=126.21,α=19.3,h=9.78,v=0.089。

表4 J4同相轴读数

Table 4 J4 cophase axis reading table

i	x_i	t_i	τ_i
1	-8	262.34	68 821.92
2	-6	251.54	63 273.59
3	-4	238.30	56 788.48
4	-2	228.15	52 052.01
5	0	220.44	48 592.48
6	2	211.70	44 815.48
7	4	208.87	43 626.22
8	6	207.20	42 930.99
9	8	208.87	43 626.22
10	10	210.93	44 489.55

新窗口打开| 下载CSV

表5 J4同相轴读数计算结果

Table 5 J4 Cophase axis reading calculation table

i	x_i	τ_i	$x_{i}^{2}$	$x_{i}^{3}$	$x_{i}^{4}$	x_iτ_i	$x_{i}^{2}$ τ_i
1	-8	68 822	64	-512	4 096	-550 575	4 404 603
2	-6	63 274	36	-216	1 296	-379 642	2 277 849
3	-4	56 788	16	-64	256	-227 154	908 616
4	-2	52 052	4	-8	16	-104 104	208 208
5	0	48 592	0	0	0	0	0
6	2	44 815	4	8	16	89 631	179 262
7	4	43 626	16	64	256	174 505	698 019
8	6	42 931	36	216	1296	257 586	1 545 516
9	8	43 626	64	512	4096	349 010	2 792 078
10	10	44 490	100	1 000	10 000	444 895	4 448 955
$\sum$	10	509 017	340	1 000	21 328	54 152	17 463 105

新窗口打开| 下载CSV

将最小二乘法数值拟合的结果和表4中图形拟合的结果对比,二者数值接近,可相互印证。可见采用图形拟合,只需在图中拾取2个点,使正演曲线与同相轴直观最佳匹配,就可达到与数值拟合相同的准确度。

6 结语

在笔者多年的探地雷达探测工作实践中,宽角反射探测方法取得了一定的成效,相比于窄角反射,从雷达信号中提取出速度参数,更能综合评价结构面对岩体完整程度的影响。为了计算更加准确,探测应注重细节,包括测线测量、采样点数、测点步距以及数据处理等,尽量做到精细化操作。

高频电磁波在低阻介质中具有高衰减性^[10],宽角反射法因雷达射线入射角与反射角变小,信号幅度衰减更快,在某些岩性中信噪比可能降低到难以分辨有效信号^[11],因此需注意岩层的岩性和完整程度。一般而言,在灰岩、白云岩、大理岩、砂岩等完整程度较好的沉积岩或变质岩中探测效果好,而在完整程度较差的火成岩中探测效果差。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

雷宛

. 工程与环境物探教程[M]. 北京: 地质出版社, 2006.