基于梯度法和L-BFGS算法的探地雷达时间域全波形反演

俞海龙， 冯 晅,2， 恩和得力海， 赵建宇， 孙成城

(1.吉林大学 地球探测科学与技术学院，长春 130026;2.地球信息探测仪器教育部重点实验室(吉林大学)，长春 130026)

0 引言

探地雷达(Ground Penetrating Radar)，是一种利用高频无线电波来确定地下物质分布规律的地球物理方法[1]，它具有采样率高、精度高、抗干扰性强、无损探测等优点。由于探地雷达具有能够对浅层地下物质高分辨率成像的特点，因此，它在工程地球物理勘探，环境研究以及基础设施研究中都被广泛地应用[2]。

探地雷达波形反演在最优化理论的框架下，充分利用雷达记录的波形、走时、相位等信息，通过拟合观测数据与模拟数据，使残差达到最小，来反演地下介质的物性参数(如电导率、介电常数、磁导率等)的方法，而电导率σ和介电常数ε最能够反映出地下介质的电磁学性质[3]，因此反演电导率和介电常数具有着重要的意义。

大量的全波形反演方法包括声波、弹性波、粘弹性波、各向异性声波方程等的有限差分或有限元法分别在时间域[4-6]和频率域[7-8]中研究并且应用于地震数据，但是到目前为止，探地雷达反演的发展仍具有一定的局限性，这是因为雷达信号的振幅和相位依赖于天线的方向以及波场的传播路径，所以我们必须了解天线的辐射方向以及雷达数据的矢量特性。对于跨孔雷达，Ernst等[9]、Kuroda等[10]提出了基于二维麦克斯韦方程的时域有限差分全波形反演方法；Meles 等[11]考虑了全波形反演电参数的矢量特性，在反演中允许电导率与介电常数同时更新；Kalogeropoulos等[12]将全波形反演应用于评价混凝土中的湿度和氯化物；王兆磊等[13]利用探地雷达全波形反演方法反演出二维介质电导率和介电常数； Lambot等[14]、Klotzsche等[15]、Streich等[16]分别研究了探地雷达波形反演。

笔者首先对TM模式下的麦克斯韦方程进行时域有限差分正演模拟，并采用单轴各向异性(UPML)吸收边界条件以消除边界反射的影响[17]，分别进行电导率和介电常数的单参数反演，以及电导率和介电常数双参数同时反演，在反演的过程中，由于不同参数之间的相互耦合、相互干涉的影响，继而增加了反演的非线性程度，另外不同参数的物理量纲不同，致使进一步增加了反演的病态性[18-19]，在多参数的全波形反演中，考虑到近似Hessian矩阵的非对角块元素反映不同参数之间的耦合效应，因此在反演中利用近似Hessian矩阵的逆来够消除不同参数之间的影响[20-21]，进而提高目标函数的收敛速度，获得准确的反演结果，采用L-BFGS算法进行探地雷达全波形反演，并和梯度法做了对比，结果显示L-BFGS算法能够更好的消除电导率和介电常数之间的影响。

1 全波形反演方法原理

1.1 探地雷达正演

众所周知，正问题是相对于反问题而言的，正演模拟的精度严重影响着反演的结果，笔者采用TM模式下的麦克斯韦方程组，利用高精度空间10阶进行波场模拟。

TM模式下麦克斯韦旋度方程[17]为式(1)。

(1)

其中：ε为相对介电常数；μ为相对磁导率；σ为电导率；Ez为电场强度z分量；Hx、Hy分别为磁场强度x、y分量。

由于实际计算区域是有限的，当波场传到边界时会发生边界反射，为消除边界反射的影响，采用吸收效果较好的单轴各向异性(UPML)吸收边界条件。

1.2 探地雷达波形反演

1.2.1 目标函数

全波形反演通过拟合正演数据和实测数据，使两者之差的平方达到最小，目标函数可以写为式(2)。

(2)

全波形反演是大规模的无约束非线性反演问题，一般采用局部优化方法来求解，模型参数的迭代更新公式可表示为式(3)。

mk+1=mk+αkpk

(3)

其中：m为待反演的参数，本文中为介电常数ε和电导率σ；αk为迭代步长；pk为迭代方向；k为迭代次数。

目标函数关于模型参数的梯度[11]为式(4)。

(4)

1.2.2 局部优化算法

局部优化算法分为梯度类算法和牛顿类算法，梯度法具有工作量少，所需存储变量少等优点，但是收敛速度慢，计算效率低，易陷入局部极小点，而牛顿法具有精度高，收敛速度快等优点，但是所需存储变量多，Hessian矩阵以及Hessian矩阵的逆求取困难，为了避免Hessian矩阵的直接求取，进而提出了拟牛顿算法，其中L-BFGS算法是一种公认的最有效的拟牛顿类算法[22]。

梯度法中模型迭代公式为式(5)。

εk+1=εk-αkε▽Sε

σk+1=σk-αkσ▽Sσ

(5)

拟牛顿法中模型迭代公式为式(6)。

(6)

(7)

对式(7)进行m次递推可得式(8)。

8)

因此，只需记录m个向量{pi,qi}，i=k-m、…、k-1，就可以按式(8)构造出近似矩阵，m的取值一般为3～20，本文取为5。

1.2.3 搜索步长

搜索步长的选取对反演结果有着重要的影响，搜索步长太小，需要增加迭代次数以使得目标函数收敛，进而增加了计算量；搜索步长太大，可能会破坏迭代过程的稳定性，且极其容易陷入局部极小值。当前常用的步长搜索方法有线性搜索方法以及抛物线拟合法等，线性搜索方法虽然简单，但是对于高维数问题的求解耗时较多；抛物线法在每一次迭代中都需进行多次正演以得到最优步长，而全波形反演需要多次迭代才能够得到最终结果，因此抛物线法的计算量也是相当大的。

笔者采用A. Pica*等[23]提出的步长求解方法，根据式(2)全波形反演目标函数为式(9)。

S(σk+1,εk+1)=S(σk+αk▽Skσ,εk+αk▽Skε)=

(9)

对目标函数S(σk+αk▽Skσ,εk+αk▽Skε)求导，令导数为零，以求取使得目标函数达到最小时所对应的步长为式(10)。

(10)

最终得到步长的表达式为式(11)。

(11)

(12)

式(11)中介电常数和电导率采用相同的步长，Meles等[11]通过试验说明式(11)求取的步长主要取决于介电常数的梯度，而电导率的梯度贡献很小，几乎不起作用，对此，分别求取电导率和介电常数的步长，电导率和介电常数的步长求取如下

(13)

(14)

其中：稳定因子δσ、δε的选取需满足式(12)，从式(13)、式(14)可以看出，用此种方法求取步长只需做两次正演模拟即可，相比较抛物线法可以大大减少正演次数，进而提高计算效率。

1.2.4 反演策略

反演主要分三步：①对介电常数单独反演，即假定电导率已知，介电常数为未知参数；②对电导率单独反演，即假定介电常数已知，电导率为未知参数；③对介电常数和电导率同时反演，即电导率和介电常数都为未知参数。在双参数同时反演时，分别采用梯度法和L-BFGS法，并就两者反演的精度和计算效率进行了对比。图1是全波形反演计算流程，该计算流程是对介电常数和电导率同时反演，对于单参数反演，只需计算相应参数的梯度、步长以及模型更新量即可。

图1 探地雷达全波形反演计算流程Fig.1 Flow chat of GPR full waveform inversion

1.2.5 时间复杂度分析

模型的空间采样间隔为0.01 m，网格数为50*50，时间采样间隔为0.02 ns，FDTD时间迭代次数为1 000次，共有9个场源位置，采用单精度(4 bit)保存，电场值保存需要内存为50*50*1000*9/1024/1024=21.457 7 M，计算反传波场同样需要内存21.457 7 M，共计42.915 4 M内存。在一次迭代过程中每一个源都需要4次正演计算，第一次得到正演数据，第二次波场反传计算梯度，另外两次得到迭代步长，因为共有9个场源，因此每一次迭代共需要36次正演计算。

2 数值试算

笔者设计的试算模型为方格模型(图(2))，图2(a)中，背景模型的相对介电常数为6，左右两个方形异常体的相对介电常数分别为5和7。图2(b)中，背景模型的电导率为0.05 S·m-1，左右两个方形异常体的电导率分别为0.01 S·m-1和0.1 S·m-1。模型大小为0.5 m×0.5 m，网格大小为50×50，网格间距为0.01 m，采用一点激发多点接收的采集方式，发射天线位于垂直位置为0 m处，发射天线初始点位于第5个网格点，发射天线之间间隔为5个网格点，共有9个发射天线，接收天线位于垂直位置为0 m以及水平位置为0 m和0.5 m处，每个位置均有50个接收点，时间采样间隔为0.02 ns，记录时间为12 ns，场源函数选用主频为800 MHz的雷克子波。

1)介电常数单独反演，介电常数真实模型为图2(a)，电导率为0.05 S·m-1，观测数据是由真实模型正演得到的，初始模型为相对介电常数为6，电导率为0.05 S·m-1的均匀半空间，反演采用L-BFGS算法，迭代次数为73次，图2(c)是介电常数反演结果，图3(a)是在深度为0.25 m处穿过两异常体中心的横向剖面，从图2(c)和图3(a)中可以看出，异常体的形态能够非常准确地刻画出来，并且在数值上，反演结果和真实模型也非常接近。图4(a)为目标函数随迭代次数变化的收敛曲线，从图4(a)中可以看出，采用L-BFGS算法，目标函数收敛速度很快，迭代35次即可收敛。

2)电导率单独反演，电导率真实模型为图2(b)，相对介电常数为6，观测数据是由真实模型正演得到的，初始模型为电导率为0.05 S·m-1，相对介电常数为6的均匀半空间，反演采用L-BFGS算法，迭代次数为100次，图2(d)是电导率反演结果，图3(b)是在深度为0.25 m处穿过两异常体中心的横向剖面，从图2(d)和图3(b)可以看出，无论在异常体的形态上还是在数值上，反演结果都接近于真实模型。图4(b)为目标函数随迭代次数变化的收敛曲线，从图中可以看出，采用L-BFGS算法，目标函数收敛速度很快，迭代40次即可收敛。

图2 梯度法单参数反演结果Fig.2 Radient method of single parameter of the full waveform inversion results(a)真实相对介电常数；(b)真实电导率；(c)单参数介电常数反演结果；(d)单参数电导率反演结果

图3 穿过两异常体中心的单参数反演结果横向剖面Fig.3 The transverse profile of through the two abnormal body center of single parameter inversion results(a)相对介电常数；(b)电导率

图4 目标函数Fig.4 Objective function(a)相对介电常数；(b)电导率

图5 梯度法双参数全波形反演结果Fig.5 Gradient method of double parameters of the full waveform inversion results(a)相对介电常数；(b)电导率,L-BFGS双参数反演结果；(c)相对介电常数；(d)电导率

图6 穿过两异常体中心的双参数反演结果横向剖面Fig.6 Through the center of the two abnormal body double parameters inversion results transverse section(a)相对介电常数；(b)电导率

3)介电常数与电导率同时反演，真实模型分别为图2(a)、图2(b)，初始模型为相对介电常数为6、电导率为0.05 S·m-1的均匀半空间，分别采用梯度法和L-BFGS算法进行反演，图5(a)、图5(b)分别是梯度法的介电常数和电导率的反演结果，图5(c)、图5(d)分别是L-BFGS算法的介电常数和电导率反演结果，图6分别是介电常数(图6(a))和电导率(图6(b))在深度为0.25 m处穿过两异常体中心的横向剖面。

从图6中可以看出，L-BFGS算法的反演结果比梯度法的反演的结果，在异常体的形态上相差不大，但是在数值上，L-BFGS算法更加接近于真实模型，并且L-BFGS算法的收敛速度和精度都高于梯度法，这是因为在多参数全波形反演中，Hessian矩阵的非对角块元素反应不同参数之间的耦合效应，而L-BFGS算法是一种拟牛顿算法，即构造出一个不含二阶导数的矩阵来近似Hessian矩阵，可以有效地解决多参数反演中不同参数之间的影响。图7为目标函数随迭代次数变化的收敛曲线，从图7中可以看出，L-BFGS算法比梯度法的目标函数收敛速度更快，目标函数下降的更小。

图7 目标函数Fig.7 Objective function

3 结论与展望

从TM模式下麦克斯韦方程组出发，采用时域有限差分的方法进行雷达波场的数值模拟，并采用单轴各向异性(UPML)吸收边界条件以消除边界反射的影响，对于探地雷达波形反演，给出了电导率和介电常数梯度方向的计算方法，并以步长为自变量通过求取目标函数为极值的方式来确定最优步长，分别对介电常数、电导率进行单参数反演，另外也对介电常数与电导率双参数同时反演，并利用拟牛顿算法L-BFGS法来消除不同参数之间的耦合影响，从反演的结果可以看出，对于单参数反演，无论是电导率还是介电常数，反演结果都十分接近于真实模型，对于双参数同时反演，反演结果的异常体形态接近于真实模型，但是由于电导率和介电常数之间的耦合影响，致使在数值上不能像单参数反演那样得到准确的反演结果，但是和真实模型相差不大。本文波形反演是在时间域进行的，今后将采用并行算法以解决巨量的正演计算，进而提高计算效率。