地球物理学报  2017, Vol. 60 Issue (9): 3616-3627   PDF    
基于一种反余切紧凑约束的电磁波走时层析研究
胡书凡1, 赵永辉1 , 葛双成2,3, 饶椿风1, 覃谭1, 安聪1     
1. 同济大学海洋与地球科学学院, 上海 200092;
2. 浙江水利水电学院, 杭州 310018;
3. 浙江省水利河口研究院, 杭州 310020
摘要:电磁波走时层析可用于空洞和孤石等离散块体异常的探测中.在传统电磁波走时层析中,通常采用一个低阶差分泛函来稳定反演过程,而这类泛函都具有光滑约束的性质,在成像结果中难以对此类离散块体异常进行准确的解释.本文基于Tikhonov正则化思想,提出了一种反余切泛函,该泛函具有紧凑约束的性质且无需选择一个聚焦因子.结合弯曲胖射线层析理论和重加权正则化共轭梯度反演算法,对两个理论模型进行了成像效果研究.与低阶差分泛函和最小支撑泛函的成像结果相比,该方法能够更好地突出异常的边界,成像结果中的虚假异常也要更少.同时,还分析讨论了激发电磁波的中心频率以及数据噪音对成像结果的影响.此外,针对一个地下连续墙墙体内部缺陷的现场模型,利用该方法取得了理想的成像结果,良好地揭示了缺陷的形态和位置.
关键词: 电磁波      走时层析      离散块体异常      反余切      紧凑约束     
Traveltime tomography of electromagnetic waves based on an anti-cotangent compactness constraints
HU Shu-Fan1, ZHAO Yong-Hui1, GE Shuang-Cheng2,3, RAO Chun-Feng1, QIN Tan1, AN Cong1     
1. School of Ocean & Earth science, Tongji University, Shanghai 200092, China;
2. Zhejiang University of Water Resources and Electronic Power, Hangzhou 310018, China;
3. Zhejiang Institute of Hydraulics and Estuary, Hangzhou 310020, China
Abstract: The traveltime tomography of electromagnetic waves can be used to detect discrete blocky anomalies, such as voids and boulders. The traditional traveltime tomography of electromagnetic waves usually uses the low-order differential operator with smooth constraints to stable the inversion process, finally resulting in a smooth solution, but leading to uncertainties in data interpretation. Here, we propose an anti-cotangent functional based on Tikhonov' regularization theory, which has the property of compactness constraints and no longer requires the user select a focusing parameter. We invert two theoretical models by the bending fat ray theory and Reweighted Regularized Conjugated Gradient (RRCG) algorithm. Compared to the result by low-order differential and minimum support functional, the proposed approach can better highlight the boundary of anomalies and suppress artifacts. In addition, some influence factors to the tomography results, including central frequency of excited electromagnetic waves and the noise contained in the data are also discussed. Furthermore, we apply the proposed method to a diaphragm wall model and obtain an interpreted image, which well reveals the shape and location of the defects distributed in the wall.
Key words: Electromagnetic-wave      Traveltime tomography      Discrete blocky anomalies      Anti-cotangent      Compactness constraints     
1 引言

基于层析成像的电磁波探测,作为一种精准、快速而无损、经济的地球物理方法,广泛应用于环境与工程(Johnson et al., 2007; 陈昌彦等,2012; 李红立等, 2016)、水文与地质(Clement and Barrash, 2006; Balasubramaniam et al., 2013; El-Behiry and Hanafy, 2013)以及矿产资源勘查等领域(Fullagar et al., 2000).该方法主要利用电磁波在不同介质中传播特性的差异来重建地下介质的电性结构.根据反演数据的不同,电磁波层析可以分为走时层析(Irving et al., 2007; 王飞等, 2013)、衰减层析(Johnson et al., 2007; Balasubramaniam et al., 2013; 董一飞等, 2017)以及全波形层析(Cordua et al., 2012; Klotzsche et al., 2013; Busch S et al., 2014; 吴俊军等, 2014; 孟旭等, 2016).全波形层析考虑到了接收信号中所有的波形信息,固然可以获得一个精度更高的结果,但由于每次迭代过程中都需要一个比较耗时的正演过程,使得该方法计算效率比较低(Cordua et al., 2012).走时层析和衰减层析分别利用直达波的到达时和振幅,来重建地下结构的速度和衰减分布,结合两者还可以得到地下的介电常数和电导率分布.而电磁波走时层析更以计算效率高、不受辐射模式和几何扩散的影响等优点(Giroux et al., 2007),在探测实践中得到了更为广泛的应用.

由于噪音和观测空间限制的影响,电磁波层析成像通常是一个不适定的问题,解的不稳定性和不唯一性是影响成像结果的最主要因素.对于这种不适定问题,通常采用Tikhonov正则化方法来稳定反演过程(Tikhonov and Arsenin, 1977).传统的Tikhonov正则化理论中通常采用一个作用于模型参数的低阶差分算子(零阶、一阶或二阶)来稳定反演过程,从而得到一个最小、平整或光滑的成像结果.当地下异常呈块状分布时,如空洞和孤石等离散异常体,这种最小、平整或者光滑的反演结果并不能良好的描述地下异常体的真实分布,甚至会产生一些虚假异常以至于得到一些错误的解释.特别是在电磁波层析问题中,由于观测空间的限制,其横向分辨率往往不足,而利用这种平滑约束更加会降低成像结果的横向分辨率.

许多学者提出了一些提高块体结构分辨率的方法.Becht等(2004)利用不同射线覆盖角度间隔的数据进行反演并解释得到一个非均匀的初始模型,并利用此模型对全射线覆盖角度数据进行反演得到了一个高块体分辨率的成像结果.Irving等(2007)提出了一种结合高射线覆盖角度数据的方法来提高电磁波走时层析的分辨率.这些方法通过一些反演策略来提高电磁波走时层析的块体分辨率,在实际应用中却是有些费时和繁琐的.另一些方法可以归结为使用不同的稳定泛函来提高块体的分辨率.Rudin等(1992)提出了一种全变差(TV)稳定泛函,该泛函使用一范数将模型参数的一阶差分最小化从而提高块体的分辨率,然而该泛函在零处不可导.Acar和Vogel (1994)在TV泛函中增加了一个小值从而解决TV泛函在零处不可导的问题.然而,TV泛函对模型参数变化的边界有惩罚作用,从这层意义上来说它仍然有光滑作用(Zhdanov, 2002).Last和Kubik (1983)提出了最小支撑(MS)泛函,该泛函能够将模型参数压缩到一个更加紧凑的区域,他们将其应用于位场数据反演中,得到的反演结果具有良好的块体分辨率.Ajo-Franklin等(2007)将MS泛函应用于地震跨孔走时层析中,重建的P波速度结构具有非常尖锐的边界.Portniaguine和Zhdanov (1999)在MS泛函的基础上提出了最小梯度支撑(MGS)泛函,该泛函将紧凑约束应用于模型的空间梯度,在位场数据反演中得到了一个相比于TV泛函和低阶差分泛函具有更高块体分辨率的反演结果.Youzwishen和Sacchi (2006)利用MGS泛函对理论合成的地震数据进行了声波速度重建.

在MS泛函中,重建图像的块体分辨率与聚焦因子的选取有关.聚焦因子太小会引起重建图像“过聚”,导致反演图像失真;而当聚焦因子过大时,重建图像会失去聚焦效果.Zhdanov等(2004)提出了一种类似于L曲线的方法来选取一个折中的聚焦因子,而Ajo-Franklin等(2007)表明当聚焦因子接近仪器的精度时,该方法会变得难以实施,反演也会变得不稳定.同时,也应该认识到在MGS泛函中会存在同样的问题.本文基于Tikhonov正则化理论提出了一种反余切(AC)紧凑约束,该方法无需选取一个聚焦因子,从而提供了一种简便的方法来提高块体的分辨率;结合弯曲胖射线理论和重加权正则化共轭梯度(RRCG)算法对电磁波走时数据进行了反演;分析和对比了低阶差分算子、MS泛函以及我们提出的AC泛函的块体分辨效果.在此基础上,讨论了激发电磁波的中心频率和数据噪音对成像结果的影响.此外,利用该方法,对一个地下连续墙缺陷的现场物理模型进行了电磁波走时层析分析,并对所有重建图像结果进行了讨论.

2 方法 2.1 电磁波走时层析理论

电磁波走时层析通常在一个孔中发射高频电磁波,电磁波在介质中以速度v传播并在另一个孔中被接收.接收的信号中包含了直达波的到达时t,通常又称作初至.可以利用该初至时来重建地下介质的速度分布,该关系可以表示为

(1)

其中G(x)描述的是初至对每个模型网格单元中慢度的敏感性程度,称为敏感核.该敏感核的计算可以基于不同的正演假设,包括射线理论(Vidale, 1988; Giroux et al., 2007)、胖射线理论(Jensen et al., 2000; Spetzler and Snieder, 2004)以及Born近似(Buursink et al., 2008; Liu et al., 2009).相比于其他两种方法,胖射线理论是对波动方程的有限频率近似,且具有较高的精度和较快的计算速度,因此在这里选择弯曲胖射线理论来对走时进行正演.

考虑到所有的观测数据,电磁波走时层析可以简化为求解如下方程,公式为

(2)

其中t为一个m×1的矩阵,包含了所有的走时数据;mn×1的矩阵,描述地下的慢度分布;Gm×n的敏感核矩阵.注意到敏感核的计算是依赖于慢度场的,因此电磁波走时层析是一个强非线性问题.

对于电磁波走时层析这种不适定问题,通常采用正则化方法来稳定反演结果.基于Tikhonov正则化思想下的反演目标函数可以定义为

(3)

其中φ(m)为数据拟合泛函,表征合成数据与观测数据的匹配程度;s(m)为稳定泛函,用来稳定反演过程;λ为正则化参数.

在二范数下,目标函数可以表示为

(4)

其中,Wd为数据加权矩阵;Wm为模型加权矩阵;mpriori为先验初始模型.在这里选择单位阵作为数据加权矩阵,意味着所有数据的权重一致.

2.2 稳定泛函及RRCG算法

惩罚泛函的引入,可以简便的将约束问题转化为无约束问题.在电磁波走时层析中,许多不同种类的惩罚泛函可以作为稳定泛函,而这些泛函又对应着不同类型的模型约束,从而得到不同类型的反演结果.在大多数实际应用中,通常选择最小范数(MN)泛函作为稳定泛函,可以表示为

(5)

可以看出此时的模型加权矩阵是一个单位阵.

同样,也可以将最小范数准则应用于模型参数的梯度Δm (Tikhonov et al., 1977, Constable et al., 1987),称作最光滑泛函(maxsm),公式为

(6)

此时的模型加权矩阵为一个一阶差分算子.在很多情况下,梯度算子也可以用拉普拉斯算子代替(Constable et al., 1987; Smith et al., 1991),公式为

(7)

相应的模型加权矩阵为一个二阶差分算子.

然而,这些低阶差分算子得到的最小、平整或光滑的地质模型,在很多情况下并不能恰当的描述地下介质的真实分布,例如空洞和孤石等离散块体异常.特别是在电磁波层析问题中,光滑约束往往会使得成像结果的横向分辨率变得更为糟糕.Last和Kubik (1993)提出了一种能够提高块体分辨率的稳定泛函,称作最小支撑(MS)泛函,公式为

(8)

其中,β称作聚焦因子,可以避免当m=mpriori时引起的奇异性.然而,利用MS泛函重建得到的图像好坏很容易受到β的影响,如果β太小,那么重建图像会“过聚”而失真;如果太大则会失去聚焦效果,无法提高块体的分辨率.Zhdanov和Tolstaya (2004)提出了一种类似于L曲线的方法来选取一个折中的聚焦因子,而Ajo-Franklin等(2007)表明当聚焦因子接近仪器的精度时,该方法会变得难以实施,反演也会变得不稳定.

在这里,基于反余切(AC)函数的性质,提出了一种新的稳定泛函,可以表示为

(9)

其中,β为一个很小的常数,用来确保该泛函在零处是可微的,在这里只需让β→0即可.

根据Tikhonov正则化理论,稳定泛函可以表示为模型参数的伪二次泛函形式,公式为

(10)

其中We可以认为是一个与模型加权矩阵Wm类似的重加权矩阵,该加权矩阵中的权重与当前的模型参数有关,可以表示为

(11)

其中,e是一个很小的常数,并且应该让e→0.可以看出,将稳定泛函表示为伪二次泛函形式后,不同稳定泛函的区别就主要在于重加权矩阵We,因此可以根据这个重加权矩阵We来分析稳定泛函的性质.

图 1为MS泛函与AC泛函的重加权矩阵随模型参数的变化关系.可以看出这两个泛函在总体上有相同的性质:与先验模型的偏差非常小的模型参数会得到一个较大的权重,使得这部分的模型参数更依赖于先验模型;而与先验模型偏差较大的模型参数权重则较小,使得这部分的模型参数更依赖于数据.这样在反演迭代过程中,与先验模型偏差较大的模型参数会得以突出,而与先验模型偏差较小的模型参数会得到压制,并且最终使得异常的结构呈块体分布.同时也可以看出,当聚焦因子较小时,MS泛函的权重主要集中于模型参数与先验模型偏差较小的区域内,但难以压制一些与先验模型偏差稍大的虚假异常;当聚焦因子较大时,权重的范围会变的更广,但过大又会使得重建图像失去聚焦效果.AC泛函的权重在模型参数与先验模型偏差较小的区域内比较接近于聚焦因子较小的MS泛函,在模型参数与先验模型偏差较大的区域比较接近聚焦因子较大的MS泛函.因此,AC泛函能够弥补MS泛函在聚焦因子过小时难以压制与先验模型偏差较大虚假异常的缺点,也不会因为聚焦因子过大而导致不聚,从而能够更好地提高块体的分辨率.

图 1 MS泛函和AC泛函的重加权矩阵We随模型参数的变化关系 Fig. 1 Reweighted matrix of MS and AC functional versus model parameters

当稳定泛函表示为模型参数的伪二次泛函形式后,使用RRCG算法(Portniaguine and Zhdanov, 1999)便可以统一不同泛函的反演迭代流程.其迭代流程为:

(12)

其中,为线性搜索步长;为目标函数的共轭梯度方向.在这里采用动态的自适应正则化参数(Zhdanov, 2002).

反演结束后,采用均方根误差(RMS)来表征数据的拟合程度,均方根误差的计算方式为

(13)

其中,M为走时数据的个数,tobstcal分别为观测走时数据和计算走时数据.

3 理论模型计算

为了验证AC泛函对块体的分辨能力,对两个理论模型进行了试验,并与低阶差分泛函和MS泛函的成像结果进行了对比.首先使用一个2D时间域有限差分算法(Irving and Knight, 2006)对模型区域的电磁场传播进行模拟,发射的子波是一个中心频率为100 MHz的Blackman-Harris脉冲.在得到了理论电磁波的波形记录后,采用改进的Coppens方法(Sabbione and Velis, 2010)来提取电磁波的初至时.最后利用RRCG算法对模型区域进行速度成像.其中MS泛函中聚焦因子的选取采用Zhdanov和Tolstaya (2004)提出的方法,首先利用MN泛函得到的反演解m0来计算MS泛函随聚焦因子β的变化曲线,公式为

(14)

然后将曲线上的最大曲率点作为一个折中值应用于反演过程中.

3.1 双椭圆模型

图 2a所示,在相对介电常数εr=9的均匀背景中存在两个相对介电常数εr=5的椭圆异常,模型中所有的电导率σ=1 mS·m-1为一个常数,磁导率等于真空中的磁导率μ0.该模型与Becht等(2004)中的模型类似.两孔之间的距离为10 m,左边有21个发射装置,间距为0.5 m,用圆圈表示,接收装置的排布与发射装置一致,用叉号表示.在正演过程中采用0.025×0.025 m的网格剖分.图 2b为发射装置位于7.5 m深处接收得到的电场垂向分量(Ez),红线为提取的初至走时.总共提取了441个初至走时数据,反演过程中采用0.2×0.2 m的网格对射线覆盖区域进行剖分.

图 2 双圆模型及初至走时拾取 (a)双椭圆模型(Becht et al., 2004),○和×分别表示发射装置和接收装置;(b)发射装置位于7.5 m深处接收到的记录,红线代表提取的初至走时. Fig. 2 Double ellipsoid model and first arrival picking (a) Double ellipsoid model (Becht et al., 2004), ○ and × represent the transmitters and receivers, respectively; (b) Received record for source located at 7.5 m depth. Red line indicates picked up first arrival.

图 3bcd为使用低阶差分算子对双椭圆模型进行图像重建的结果.可以看出,这些成像结果的均方根误差都比较大,模型参数与背景的差异都非常小,并且异常在水平方向上都非常发散.利用零阶差分算子重建的图像中异常的最大幅值约为0.117 m·ns-1;一阶差分算子重建的图像中约为0.1125 m·ns-1;二阶差分算子重建的图像中约为0.117 m·ns-1.这与理论值v=0.134 m·ns-1差异非常大.同时,在图像中难以区分两个椭圆异常,在成像结果中很容易把异常解释为单个连通的高速体.

图 3 双椭圆模型成像结果 (a)初始模型;(b)—(f)使用零阶差分算子、一阶差分算子、二阶差分算子、MS泛函(β=0.48) 和AC泛函得到的成像结果. Fig. 3 Tomographic images of double ellipsoid model (a) Initial model; (b)—(f) Reconstructed images using zero-order differential, first-order differential, and second-differential operators, and MS (β=0.48) and AC functional.

图 3ef为利用MS泛函和AC泛函得到的成像结果.成像结果明显地揭示了模型中的两个高速异常,相比于低阶差分泛函,MS泛函和AC泛函得到的结果中异常的幅值更大,成像结果要更加接近于真实模型,利用MS泛函得到的成像结果中异常的最大幅值约为0.123 m·ns-1;AC泛函得到的成像结果中约为0.126 m·ns-1.同时,虽然在MS泛函的成像结果能够将两个异常体在水平方向区分开,但是背景中含有许多虚假异常,而利用AC泛函得到的成像结果中异常边界要更为明显,且背景中的虚假异常更少,均方根误差也最低.

3.2 多块体模型

选用的多块体模型与Irving等(2007)中的模型类似,模型中包含多个块体异常.图 4a给出了各块体的相对介电常数分布,模型大小为5×12 m,正演网格剖分为0.025×0.025 m,均匀背景的相对介电常数εr=25,红色块体异常的相对介电常数εr=28,蓝色块体异常的相对介电常数εr=22,模型中的电导率σ=1 mS·m-1为一个常数,磁导率等于真空中的磁导率μ0.两孔之间的水平距离为4 m,分别放置了发射装置和接收装置,发射装置从0.5 m深移动到11.5 m深,间隔0.25 m,接收装置的排布与发射装置一致.图 4b为发射装置位于5.25 m深处接收得到的电场垂向分量(Ez),红线为提取的初至走时.总共提取了2025个初至走时数据,采用0.1×0.1 m的网格剖分对射线覆盖区域进行反演.

图 4 多块体模型及初至走时拾取 (a)多块体模型(Irving et al., 2007),○和×分别表示发射装置和接收装置;(b)发射装置位于5.25 m深处接收到的记录,红线代表提取的初至走时. Fig. 4 Multi-block model and first arrival picking (a) Multi-block model (Irving et al., 2007), ○ and × represent the transmitters and receivers, respectively; (b) Received record for source located at 5.25 m depth. Red line indicates picked up first arrival.

图 5bcd为使用低阶差分算子得到的成像结果,可以看出这些重建图像可以反映出不同差分算子的性质.在零阶差分算子的重建图像中,异常的模型参数与背景的差异非常小,在3.5 m深的小块体的速度约为0.05925 m·ns-1,这与初始模型中的背景值0.05996 m·ns-1非常接近.在一阶差分算子的重建图像中,异常体的边界非常光滑,浅部的两个低速体容易被当作是连通的,而中间的两个高速体像是一个阶梯状异常,并且图像中存在很多虚假异常,特别是在发射装置和接收装置处.相比于一阶差分泛函,在二阶差分算子的重建图像中的虚假异常有所改善,同时图像也要更为光滑.可以看出,当异常呈块体分布时,使用低阶差分算子难以准确的描述异常的位置和形态.

图 5 多块体模型成像结果 (a)初始模型;(b)零阶差分算子; (c)一阶差分算子; (d)二阶差分算子; (e) MS泛函(β=0.39); (f) AC泛函. Fig. 5 Tomographic images of multi-block model (a) Starting model; (b) Zero-order differential operator; (c) First-order differential operator; (d) Second-differential operator; (e) MS functional (β=0.39); (f) AC functional.

图 5ef为MS泛函和AC泛函的成像结果.利用这两个泛函重建图像的均方差比低阶差分泛函更小,并且能够良好的揭示中间两个高速块体异常,异常的幅值约为0.0642 m·ns-1,这与理论值0.064 m·ns-1非常接近.同时可以看出,AC泛函重建的图像中块体分辨能力比MS泛函更高,3.5 m处的小块体异常的形态更接近于真实模型,背景中的虚假异常也较MS泛函更少.

在实际应用中,数据通常是含有噪音的.因此,在多块体模型的初至数据中加入了均值为0方差为1 ns的高斯噪音,图 6为对加噪数据利用不同泛函得到的成像结果.与无噪音的数据类似,相比于低阶差分泛函,具有紧凑约束性质的泛函能够压制由于观测空间限制引起的虚假异常.同时,紧凑约束还能将噪音引起的虚假异常也压缩到最小的区域内(图 6de),从而得到一个更易于解释的成像结果,但是异常的形态和位置会发生一定的变化.此外,可以看出AC泛函得到的成像结果中虚假异常要比MS泛函中的更少,且均方根误差也要更低,为1.00 ns,这与加入高斯噪音的量级非常吻合,这说明AC泛函比MS泛函更能够处理含有噪音的数据.

图 6 含噪音多块体模型的成像结果 (a)零阶差分算子; (b)一阶差分算子; (c)二阶差分算子; (d) MS泛函(β=0.39); (e) AC泛函. Fig. 6 Tomographic images of noise-bearing multi-block model (a) Zero-order differential operator; (b) First-order differential operator; (c) Second-differential operator; (d) MS functional (β=0.39); (e) AC functional.

为了分析激发电磁波的中心频率对成像结果的影响,采用250 MHz的Blackman-Harris激发脉冲对多块体模型中的电磁波传播进行模拟并提取了相应的初至.图 7a为利用AC泛函得到的成像结果,相比于利用低中心频率重建的图像(图 5f),高中心频率重建的图像中块体分辨率更高,异常的形态也更接近于真实模型,5 m左右深处的两个高速异常体更容易区分开.但是,利用高中心频率重建的图像也更容易受到观测空间限制的影响,在发射装置和接收装置处出现了很多虚假异常.在利用高中心频率生成的数据中加入了图 6中同样的高斯噪音,成像结果如图 7b所示.可以看出当激发脉冲的中心频率变高时,异常的形态更容易受到噪音的影响而发生变化,同时重建图像中的虚假异常也会变多.因此,在实际工作中选择发射电磁波的中心频率时,不仅要考虑不同中心频率电磁波具有不同的传播距离和分辨率这一特点,还需考虑高中心频率的成像结果更易受到噪音的影响这一因素.

图 7 利用AC泛函对多块体模型的成像结果,激发脉冲的中心频率为250 MHz. (a)使用无噪音数据;(b)使用含噪音数据. Fig. 7 Tomographic images of multi-block model by AC functional with central frequency of the exciting pulse 250 MHz (a) Using noise-free data; (b) Using noise-bearing data.
4 物理模型试验

为测试本文所提出的基于AC泛函的电磁波走时层析成像方法应对实际数据的有效性,将该方法应用于一个地下连续墙物理模型的试验中.该模型位于上海市新江湾城的一个实验基地,地下连续墙模型长、宽、深分别为5 m、3 m、4 m.墙体材料为钢筋混凝土,混凝土强度C30,在混凝土内布置了参照地下连续墙实体的钢筋笼.模型照片如图 8a所示.在墙1、2、3内布置了缺陷,缺陷位置及尺寸如图 8bcde所示.缺陷1、缺陷2、缺陷3、缺陷4、缺陷5分别为夹泥量为30%、纯夹泥(100%泥)、夹泥量为30%、离析(100%砂)和离析(100%骨料)的长方体,平行设置在墙身内部,不贯通连续墙.每个孔中都放置了一个PVC测管,管中充填了水.使用中心频率为500 MHz,带宽为100~800 MHz的步进频率式钻孔雷达进行数据测量,发射装置和接收装置的排布如图 8cde所示.

图 8 地下连续墙物理模型及缺陷分布 (a)地下连续墙物理模型;(b)模型平面图,图中单位为毫米;(c)—(e)墙1、2、3的缺陷布置图,实心圆圈和叉号分别代表发射装置和接收装置. Fig. 8 Diaphragm showing wall physical model and distribution of defects (a) Wall model; (b) Plane graph of model, unit: mm; (c)—(e) Location of defects in walls 1, 2 and 3. Solid circles and crosses denote transmitter and receiver, respectively.

在开始测量前,把发射和接收装置放置在空气中测得了理论时延.在测量完成后,对接收到的数据采用了一些数据处理方法来提高数据的质量,包括零时校正和一些频率滤波方法.然后采用改进的Coppens方法提取了直达波的初至,分别提取到了810、899和870个走时数据.由于长发射-接收距离的电磁波信号较弱,容易被噪音掩盖,使用这类数据进行成像会严重影响成像质量,因此去除了一些信噪比非常低的数据,分别剩下605、593和637个走时数据.最后采用AC泛函和RRCG算法对三个墙内的速度分布进行了成像.

图 9abc分别为墙1、2、3的电磁波走时层析成像结果,从成像的结果中可以良好地判断缺陷的位置和边界.墙1的成像结果中良好地揭示了两个低速异常缺陷,缺陷的位置与真实模型非常接近,重建图像中的虚假异常也非常少;由于去除了一些大射线覆盖角度的数据,墙2的成像结果在水平方向有些发散,但在垂向方向仍然能够良好地确定缺陷的边界,缺陷的位置与实际模型较吻合,同时可以看到,由于没有足够的射线覆盖缺陷3(缺陷2则是无射线覆盖),其引起的异常发散到了发射装置和接收装置处,但仍然有低速异常的反映;墙3中两孔之间的距离最短,相应的成像结果中分辨率也最高,有效地揭示了缺陷的位置和形态,从图像中异常的速度值v≈0.034 m·ns-1可以推断该空洞中充填了水.

图 9 基于AC泛函的地下连续墙内部缺陷成像结果 (a)墙1(缺陷1及空洞1); (b)墙2(缺陷2~5及空洞2); (c)墙3(空洞3). Fig. 9 Tomographic images for the diaphragm wall model based on AC functional (a) Wall 1 (defect 1 and void 1); (b) Wall 2 (defects 2~5 and void 2); (c) Wall 3 (void 3).
5 结论

本文基于Tikhonov正则化思想,提出了一种带有紧凑约束的反余切稳定泛函,该泛函无需选择聚焦因子,因此可以更加方便快速的应用于电磁波走时层析中.利用重加权正则化共轭梯度算法对2个理论模型和1个物理模型进行了成像,得到如下结论:

(1) 当地下介质异常呈块体分布时,使用传统的低阶差分泛函难以准确的描述异常的真实形态,不能对真实的地下异常做出合理的解释,而使用具有紧凑约束性质的泛函可以提高块体的分辨率,重建图像能够良好的突出块体的边界,异常的形态更加接近于真实模型.

(2) 最小支撑泛函可以用来提高块体的分辨率,但需要人为的选择一个聚焦因子,本文提出的反余切泛函无需选择聚焦因子.通过理论模型对比表明,相比于最小支撑泛函,反余切泛函可以更好压制成像结果中的伪像异常,成像结果更接近于实际地下模型.

(3) 紧凑约束不仅能够压制由于观测空间限制引起的虚假异常,同时还能将数据噪音引起的虚假异常压缩到最小的区域内.此外,相比于最小支撑泛函,反余切泛函能够更好地压制数据噪音引起的虚假异常,从而得到一个更易于解释的成像结果.

(4) 当激发电磁波的中心频率更高时,重建图像的块体分辨率也将提高,但成像结果中的虚假异常也会变得更多.因此,在实际选取激发电磁波的中心频率时,不仅要考虑不同中心频率电磁波具有不同的传播距离和分辨率这一特点,还需考虑高中心频率的成像结果更易受到噪音的影响这一因素.

参考文献
Acar R, Vogel C R. 1994. Analysis of bounded variation penalty methods for ill-posed problems. Inverse Problems, 10(6): 1217-1229. DOI:10.1088/0266-5611/10/6/003
Ajo-Franklin J B, Minsley B J, Daley T M. 2007. Applying compactness constraints to differential traveltime tomography. Geophysics, 72(4): R67-R75. DOI:10.1190/1.2742496
Balasubramaniam V R, Jha P C, Chandrasekhar E, et al. 2013. Imaging weak zones in the foundation using frequency domain attenuation tomography. Journal of Applied Geophysics, 97: 97-106. DOI:10.1016/j.jappgeo.2013.06.003
Becht A, Tronicke J, Appel E, et al. 2004. Inversion strategy in crosshole radar tomography using information of data subsets. Geophysics, 69(1): 222-230. DOI:10.1190/1.1649390
Busch S, Van Der Kruk J, Vereecken H. 2014. Improved characterization of fine-texture soils using on-ground GPR full-waveform inversion. IEEE Transactions on Geoscience and Remote Sensing, 52(7): 3947-3958. DOI:10.1109/TGRS.2013.2278297
Buursink M L, Johnson T C, Routh P S, et al. 2008. Crosshole radar velocity tomography with finite-frequency Fresnel volume sensitivities. Geophysical Journal International, 172(1): 1-17. DOI:10.1111/j.1365-246X.2007.03589.x
Chen C Y, Shen X K, Su Z F, et al. 2012. Application of electromagnetic wave tomography on geotechnical investigation of the complex geological slope engineering. Progress in Geophysics, 27(2): 796-803. DOI:10.6038/j.issn.1004-2903.2012.02.048
Clement W P, Barrash W. 2006. Crosshole radar tomography in a fluvial aquifer near boise, idaho. Journal of Environmental and Engineering Geophysics, 11(3): 171-184. DOI:10.2113/JEEG11.3.171
Constable S C, Parker R L, Constable C G. 1987. Occam's inversion:A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. DOI:10.1190/1.1442303
Cordua K S, Hansen T M, Mosegaard K. 2012. Monte Carlo full-waveform inversion of crosshole GPR data using multiple-point geostatistical a priori information. Geophysics, 77(2): H19-H31. DOI:10.1190/geo2011-0170.1
Dong Y F, Zhang Z F. 2017. Regularization method of radio wave tomography based on Fresnel zone theory. Progress in Geophysics, 32(1): 18-25. DOI:10.6038/pg20170103
El-Behiry M G, Hanafy S M. 2013. Imaging soil moisture using GPR tomography and reflection field experiments. Arabian Journal of Geosciences, 6(9): 3493-3503. DOI:10.1007/s12517-012-0577-7
Fullagar P K, Livelybrooks D W, Zhang P, et al. 2000. Radio tomography and borehole radar delineation of the McConnell nickel sulfide deposit, Sudbury, Ontario, Canada. Geophysics, 65(6): 1920-1930. DOI:10.1190/1.1444876
Giroux B, Gloaguen E, Chouteau M. 2007. Bh_tomo-a Matlab borehole georadar 2D tomography package. Computers & Geosciences, 33(1): 126-137. DOI:10.1016/j.cageo.2006.05.014
Irving J, Knight R. 2006. Numerical modeling of ground-penetrating radar in 2-D using MATLAB. Computers & Geosciences, 32(9): 1247-1258. DOI:10.1016/j.cageo.2005.11.006
Irving J D, Knoll M D, Knight R J. 2007. Improving crosshole radar velocity tomograms:A new approach to incorporating high-angle traveltime data. Geophysics, 72(4): J31-J41. DOI:10.1190/1.2742813
Jensen J M, Jacobsen B H, Christensen-Dalsgaard J. 2000. Sensitivity kernels for time-distance inversion. Solar Physics, 192(1-2): 231-239.
Johnson T C, Routh P S, Barrash W, et al. 2007. A field comparison of Fresnel zone and ray-based GPR attenuation-difference tomography for time-lapse imaging of electrically anomalous tracer or contaminant plumes. Geophysics, 72(2): G21-G29. DOI:10.1190/1.2431638
Klotzsche A, Van Der Kruk J, Linde N, et al. 2013. 3-D characterization of high-permeability zones in a gravel aquifer using 2-D crosshole GPR full-waveform inversion and waveguide detection. Geophysical Journal International, 195(2): 932-944. DOI:10.1093/gji/ggt275
Last B J, Kubik K. 1983. Compact gravity inversion. Geophysics, 48(6): 713-721. DOI:10.1190/1.1441501
Li H L, Pan D M, Zhang H, et al. 2016. Comparison study on tomography technique based on elastic wave velocity and electromagnetic wave attenuation properties of spherical weathered granite. Geotechnical Investigation & Surveying, 44(10): 72-78.
Liu Y Z, Dong L G, Wang Y W, et al. 2009. Sensitivity kernels for seismic Fresnel volume tomography. Geophysics, 74(5): U35-U46. DOI:10.1190/1.3169600
Meng X, Liu S X, Fu L, et al. 2016. Frequency domain waveform inversion of cross-hole GPR data based on a logarithmic objective function. Chinese Journal of Geophysics, 59(5): 1875-1887. DOI:10.6038/cjg20160530
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596
Rudin L I, Osher S, Fatemi E. 1992. Nonlinear total variation based noise removal algorithms. Physica D:Nonlinear Phenomena, 60(1-4): 259-268. DOI:10.1016/0167-2789(92)90242-F
Sabbione J I, Velis D. 2010. Automatic first-breaks picking:New strategies and algorithms. Geophysics, 75(4): V67-V76. DOI:10.1190/1.3463703
Smith R T, Zoltani C K, Klem G J, et al. 1991. Reconstruction of tomographic images from sparse data sets by a new finite element maximum entropy approach. Applied Optics, 30(5): 573-582. DOI:10.1364/AO.30.000573
Spetzler J, Snieder R. 2004. The Fresnel volume and transmitted waves. Geophysics, 69(3): 653-663. DOI:10.1190/1.1759451
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. New York:Winston.
Vidale J. 1988. Finite-difference calculation of travel times. Bulletin of the Seismological Society of America, 78(6): 2062-2076.
Wang F, Liu S X, Qu X X, et al. 2013. Crosshole radar traveltime tomography without ray tracing using the high accuracy fast marching method. Chinese Journal of Geophysics, 56(1): 3896-3907. DOI:10.6038/cjg20131131
Wu J J, Liu S X, Li Y P, et al. 2014. Study of cross-hole radar tomography using full-waveform inversion. Chinese Journal of Geophysics, 57(5): 1623-1635. DOI:10.6038/cjg20140525
Youzwishen C F, Sacchi M D. 2006. Edge preserving imaging. Journal of Seismic Exploration, 15(1): 45-57.
Zhdanov M, Tolstaya E. 2004. Minimum support nonlinear parametrization in the solution of a 3D magnetotelluric inverse problem. Inverse Problems, 20(3): 937-952. DOI:10.1088/0266-5611/20/3/017
Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. New York:Elsevier.
陈昌彦, 沈小克, 苏兆锋, 等. 2012. 电磁波层析成像技术在复杂地质边坡工程勘察中的应用研究. 地球物理学进展, 27(2): 796–803. DOI:10.6038/j.issn.1004-2903.2012.02.048
董一飞, 张致付. 2017. 基于菲涅尔带理论的无线电波层析正则化方法研究. 地球物理学进展, 32(1): 18–25. DOI:10.6038/pg20170103
李红立, 潘冬明, 张华, 等. 2016. 基于孤石的弹性波波速和电磁波衰减特性层析成像技术对比研究. 工程勘察, 44(10): 72–78.
孟旭, 刘四新, 傅磊, 等. 2016. 基于对数目标函数的跨孔雷达频域波形反演. 地球物理学报, 59(5): 1875–1887. DOI:10.6038/cjg20160530
王飞, 刘四新, 曲昕馨, 等. 2013. 基于HAFMM的无射线追踪跨孔雷达走时层析成像. 地球物理学报, 56(11): 3896–3907. DOI:10.6038/cjg20131131
吴俊军, 刘四新, 李彦鹏, 等. 2014. 跨孔雷达全波形反演成像方法的研究. 地球物理学报, 57(5): 1623–1635. DOI:10.6038/cjg20140525