崔岩,王彦飞

1 山东科技大学地球科学与工程学院，青岛 266590 2 河北地质大学，京津冀城市群地下空间智能探测与装备重点实验室，石家庄 050031 3 中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室，北京 100029 4 中国科学院大学,北京 100049 5 中国科学院地球科学研究院，北京 100029

0 引言

近年来，与环境、工程、能源密切相关的近地表地球物理探测技术及方法成为研究热点，横波速度是近地表最关键的地震参数之一，目前主要用面波来提取横波速度.面波速度建模有三个步骤：(1)数据的采集和处理；(2)频散曲线拾取；(3)频散曲线反演.频散曲线反演由频散曲线获取地下横波速度，是其中的关键步骤.目前人们最常使用的面波是瑞雷波，它有很强的刻画浅地表/浅水介质(沉积物)的能力(夏江海，2015；Wang et al.,2020).瑞雷波在层状介质中发生频散，并且会发育多阶模态，即一个频率对应多个不同的相速度值(牛滨华和孙春岩，2002).对某给定频率，它所对应的最小相速度值称为该频率的基阶模式相速度(一阶模式相速度)，其余相速度值可统称为该频率的高阶模式相速度.高阶模式相速度按照从低到高的顺序，又分别称为一阶高阶模式相速度(二阶模式相速度)、二阶高阶模式相速度(三阶模式相速度)等.在对瑞雷波频散曲线进行反演时，既可以使用单阶模式，也可以使用多阶模式.一般情况下，地层刚度随深度逐层增加，此时瑞雷波能量大部分集中在基阶模式，使用基阶模式就能够取得好的效果.但由于实际地层结构是复杂多变的，瑞雷波能量并不总是集中在基阶模式，这就使得仅用基阶模式频散曲线反演出的横波速度偏离真实值.我们可以从数学角度把这看作一个加权问题，各阶模式的加权系数是其能量占总能量的百分比，它决定了各阶模式频散曲线在反演横波速度模型时影响的大小.前人通过分析指出，在某些特殊地质结构中(例如含低速夹层)，高阶模式较基阶模式在某些频段可能更加发育，仅用基阶模式频散曲线反演出的横波速度存在失真(徐佩芬等，2020).而多阶模式频散曲线联合反演不仅能使反演过程稳定、提高横波速度的反演精度和计算结果的可靠性，还能够提高模型纵向分辨率(宋先海，2008；潘冬明等，2010；夏江海等，2015).在实际工作中，频散曲线提取的误差会随着阶数的增加而增大，为了保证结果的可靠性，本文研究的瑞雷波多阶模式频散曲线反演使用基阶、一阶高阶和二阶高阶模式频散曲线.

瑞雷波多阶模式频散曲线反演是一个典型的地球物理反问题，任何反问题的研究都涉及正演模拟、反演建模和数值实现三个方面(王彦飞等，2011).利用频散曲线反演地下速度结构，频散曲线的计算是首先要解决的正问题，对此前人已进行了大量研究，研究基于水平层状介质模型，主要算法有Haskell算法(Haskell,1953)、Schwab-Knopoff算法(Schwab and Knopoff,1970)、广义反射-透射系数算法(Chen,1993,1999)等.其中广义反射-透射系数算法在高频更为稳定、精确，在求取频散曲线的同时还能给出瑞雷波本征位移、应力及理论地震图的求解公式.在反演模型的建立上，现有文献大多仅考虑数据的拟合，将频散曲线理论计算值与实测值之间误差的L2范数作为目标函数(蔡伟等，2017).但从连续介质物理学的概念出发，地球物理模型最好采用“逐块”光滑的间断函数集来描述，允许空间存在间断面，间断面之间的介质参数是连续且光滑的(杨文采，1997)，这种物理上的间断可以通过某些模型参数在数学上的稀疏性进行描述.L1范数正则化将待反演参数的L1范数或其稀疏表示后的L1范数最小作为约束条件，目前已应用于反射系数反演(Wang et al.，2016)、波阻抗反演(Liu et al.，2018)、叠前弹性参数反演(刘晓晶等，2016)和基阶模式频散曲线反演(Cui and Wang,2021)等，研究表明L1范数正则化可以提高分辨率，增强对非高斯噪声的鲁棒性(王彦飞等，2021).在反问题的数值实现上，现有方法可以分为线性和非线性两大类(Wang et al.，2011b)，在瑞雷波频散曲线反演中使用的线性方法主要有最小二乘法(Gabriels et al.，1987)、奇异值分解和Levenberg-Marquardt法(Xia et al.，1999)、Occam算法(Lai et al.，2002)等；非线性方法有遗传算法(Yamanaka and Ishida，1996)、小波分析(Tillmann，2005)、粒子群算法(Song et al.，2012)，洗牌蛙跳算法(Sun et al.，2017)、蜂群算法(于东凯等，2018)等.线性方法结构简单、计算效率高、所需内存少，但对初始模型依赖性较强，容易陷入局部极值；非线性方法对初始模型的要求不高，具有较强的全局寻优能力，但收敛速度较慢，运算量大，目前线性方法的应用更为广泛.

本文针对目前瑞雷波多阶模式频散曲线反演中仅考虑数据的拟合，缺乏对模型约束的问题，研究了瑞雷波多阶模式频散曲线的稀疏正则化反演方法.首先讨论瑞雷波多阶模式频散曲线正演的原理及数值计算，原理采用经典的广义反射-透射系数法，数值计算上采用具有三阶收敛速度的快速求根方法.然后是反演模型的建立及反问题的数值实现，通过在目标函数中引入模型的L1范数施加稀疏正则化约束，使反演结果更加符合地质实际，并针对该模型提出一种隐式迭代正则化算法.最后将新的反演方案应用于理论模型，并与传统方法进行对比.

1 正问题

1.1 瑞雷波频散曲线正演原理

所有反演方法都涉及理论数据的模拟，利用频散曲线反演地下速度结构，频散曲线的计算是首先要解决的正问题，本文选用广义反射-透射系数法.简正振型(Seismic Normal Modes)是无源弹性动力学方程在特定边界条件下的非零解，从这个基本原理出发，计算任意给定频率和任意固体介质参数下水平层状介质中的简正振型，从而得到频散曲线.

图1 广义反射-透射系数法示意图Fig.1 Schematic diagram of generalized reflection-transmission coefficient method

(1)

当j=1,2,…,N时

(2a)

当j=0时(自由界面处)

(2b)

(3)

由此可得

(4)

线性齐次代数方程组存在非零解的充要条件是系数行列式为零，因此

上式左端是瑞雷波频散函数(久期函数)，它是一个关于模型参数(拉梅常数、密度、层厚)、频率和瑞雷波相速度的函数.给定模型参数，选取某一频率，只有有限个相速度能够满足方程(5)，它们就是频散函数的根.这一组根由小到大分别对应于该频率成分的基阶模式相速度(一阶模式相速度)、一阶高阶模式相速度(二阶模式相速度)、二阶高阶模式相速度(三阶模式相速度)，以此类推，所有频率成分和它们所对应的某模式相速度之间的关系就是该模式的频散曲线.

1.2 瑞雷波多阶模式频散曲线的数值计算

频散函数的快速求解直接影响反演过程的稳定性和结果的准确性.传统的二等分方法从某个初始值(通常是第一层S波速度的0.8倍)开始，并沿相速度方向以一定步长执行根搜索.当满足给定的误差精度时，认为获得了该频率下的相速度.然后以相同的方式滑动频率以获取整个频散曲线.该方法不具有超线性收敛性，计算效率低.当复杂的地层模型引起“模式接吻”或“模式跳变”时，可能会发生严重的缺陷，例如丢根和模式误判(夏江海，2015).因此有必要在正演模拟中引入先进的寻根方法，以提高计算的效率和准确性，故本文考虑快速根搜索方法，该方法具有三阶收敛速度(Wang and Xiao，2001).

对一个地质模型，给定频率ω，瑞雷波相速度VR要满足方程(5).为了简化公式，将方程(5)的左侧定义为非线性函数φ(x)，其中x指相速度VR.方程(5)转换为寻找非线性方程φ(x)=0的根.众所周知，φ(x)相对于x是无限可微的，在项(x-xn)3处截断泰勒展开式，可以得到

(6)

其中ξn的取值范围在x和xn之间.

通过式(6)可以获得迭代公式

xn+1=xn-

(7)

上述频散曲线求根的迭代过程可以通过图2所示的频散曲线计算流程实现.

图2 瑞雷波多阶模式频散曲线计算流程图Fig.2 Flowchart of Rayleigh wave multiple mode dispersion curve extraction

通过上述算法得出Rayleigh波多阶模式频散曲线，该曲线作为理论模拟数据用于反演地下横波速度.

2 反问题

2.1 瑞雷波多阶模式频散曲线反演建模

通过上面的讨论可以看出，瑞雷波频散曲线与模型参数有关.正演是根据模型计算频散曲线，反演则是根据频散曲线反推模型.要讨论瑞雷波多阶模式频散曲线联合反演的建模问题，首先需要对单阶模式频散曲线反演建模.

2.1.1 单阶模式

根据式(5)，对于瑞雷波单阶模式的频散曲线，频率成分ωi和它所对应的相速度VRi之间的关系可表示为

Fi(ωi,VRi,λ,μ,ρ,h)=0i=1,2,…,M

(8)

其中，M为频散曲线数据点的个数，ωi为频率，VRi为频率ωi对应的瑞雷波相速度，λ=[λ(1),λ(2),…,λ(N+1)]、μ=[μ(1),μ(2),…,μ(N+1)]、ρ=[ρ(1),ρ(2),…,ρ(N+1)]、h=[h(1),h(2),…,h(N)]分别为每一层的拉梅常数、密度、层厚组成的向量，它们共同描述了该水平层状模型.

Fi(ωi,VRi,VP,VS,ρ,h)=0i=1,2,…,M

(9)

由上式可知，横波速度、纵波速度、密度和层厚都会对瑞雷波频散曲线产生影响.但如果同时反演横波速度、纵波速度、密度和层厚，反演过程将极不稳定，存在严重的多解性.为此，我们希望减少方程中的未知量，专注于主要参数.

研究表明，横波速度是瑞雷波频散曲线最主要的影响因素.因此在反演过程中可以假定密度和纵波速度已知(夏江海，2015)，将层厚固定为合理的薄层(常数).这样只有横波速度是未知的，式(9)可以简化为

Fi(ωi,VRi,VS)=0,i=1,2,…,M

(10)

式中，VRi为该频率对应的瑞雷波某阶模式相速度，VS为横波速度向量.

2.1.2 多阶模式

本文研究基阶模式、一阶高阶模式和二阶高阶模式频散曲线的联合反演.由于高阶模式是在频率大于一定数值(截止频率)后才逐渐出现的，因此式(10)改写为

Fi0(ωi0,VR0,i0,VS)=0 (i0=1,2,…,M0),

Fi1(ωi1,VR1,i1,VS)=0 (i1=1,2,…,M1),

Fi2(ωi2,VR2,i2,VS)=0 (i2=1,2,…,M2),

(11)

2.2 非线性反演问题的正则化

2.2.1 Tikhonov正则化

广义Tikhonov正则化是解下面的最优化问题

(12)

方程右侧第一项表示数据的拟合，第二项Ω(·)是稳定因子，α>0是正则化因子，m是模型参数，d是观测数据，F(m)是由模型到数据的正演算子.

2.2.2 稀疏约束正则化

选择不同的稳定因子Ω(·)会导致不同的正则化方案，考虑解的稀疏约束，这里的“稀疏”一词用于描述具有相对较少非零值的时间序列，在其他量为常数的情况下，稀疏度越大，反演结果中的层数越少.由于反演的横波速度模型是一维的，因此考虑具有一维稀疏性的L1范数正则化.L1范数正则化将待反演参数的L1范数或其稀疏表示后的L1范数最小作为约束条件，可以提高分辨率，增强对非高斯噪声的鲁棒性(王彦飞等,2011；Wang et al.，2011a).此处选取离散模型的l1范数作为Ω(·)，即

=J1(m)+J2(m).

(13)

依据式(13)，我们构造瑞雷波基阶、一阶高阶和二阶高阶模式频散曲线联合反演的目标函数

(14)

注：(1)瑞雷波多阶模式频散曲线反演在方法上与基阶模式频散曲线反演区别不大，只需依据数据精度给予高阶模式数据相同或不同的权重，然后将其作为数据集加入反演过程(Xia et al.，2003).在数值实验中，数据通常由正演模拟产生，不同模式频散曲线数据的精度一致，故选择相同权重具有合理性；而在实际资料中，由于数据的采集、处理和频散曲线提取都会对数据精度产生影响，不同模式频散曲线数据的精度可能存在较大差异，故需要根据具体情况合理设置权重.(2)正则化因子影响反演结果的稳定性，根据反问题的正则化理论，正则化因子应当是噪声水平δ的等价无穷小(王彦飞，2007).在数值实验中，噪声水平δ是已知的，当信噪比远远大于1时，可以直接通过δ2-ε(ε是一个小的正数)计算正则化因子；在实际资料中，噪声水平通常不确定，但可依据上述参数选择理论对正则化因子进行几何级数选取.

2.2.3 隐式迭代正则化算法

编程最简单，最容易实现的梯度方法是最速下降法

mk+1=mk+ωksk,(15)

在第k步，梯度可由下式计算

(16)

mk+1=mk+Δm,(17)

其中Δm是更新项.

Fl(m)的线性化格式可以写作

Fl(mk+1)=Fl(mk+Δm)=Fl(mk)+HlkΔm

+O‖(Δm)2‖,(18)

(19)

J2可近似为

(20)

(21)

其中γ(·)是向量函数，表示为

(22)

因此目标函数Jα(mk+1)的梯度可以近似为

+αγ(Δmp),(23)

从而得到隐式迭代公式

(24)

即

(25)

给定Δmp的初始猜测值，下一步迭代将满足

(26)

这个线性系统可以通过共轭梯度法求解，且‖Δmp+1-Δmp‖≤δ时迭代终止，其中δ>0是一个很小的正数.

获得最佳解Δmp后，可以通过以下方式更新外部迭代点

mk+1=mk+Δmp.

(27)

由于极小化问题的凸性，采用上述共轭梯度算法可以收敛到极小化问题的解.上述迭代求解横波速度的方法定义为瑞雷波多阶模式频散曲线稀疏正则化反演算法(SIRWDC，Sparse inversion of Rayleigh wave dispersion curve)，算法描述如下：

第一步，给定最大外部迭代步数kmax，最大内部迭代步数pmax，Δmp的初始猜测值，终止误差δ，p=0，k=0；

第二步，使用共轭梯度法根据式(26)进行迭代，直到‖Δmp+1-Δmp‖≤δ或p=pmax，得到Δmp；

第三步，根据式(27)更新模型；

第四步，如果‖mk+1-mk‖≤δ，迭代终止；否则令k:=k+1，回到第二步.

3 数值实验

为了测试新算法的有效性和可靠性，设计半空间上覆四套地层的理论模型，模型参数如表1所示.基于广义反射-透射系数法对上述模型的瑞雷波多阶模式频散曲线(基阶、一阶高阶和二阶高阶)进行正演模拟，总频带范围取实际生产中常用的1～100 Hz，在数值计算时采用快速根搜索方法，结果如图3所示.

表1 地层模型参数Table 1 Layer model′s parameters

图3 正演模拟得到的瑞雷波多阶模式(基阶、一阶高阶和二阶高阶)频散曲线Fig.3 Forward simulation of Rayleigh wave multiple mode dispersion curves (fundamental mode,first-order mode,second-order mode)

3.1 理想数据

将正演模拟得到的多阶模式频散曲线作为观测数据，分别采用SIRWDC方法和传统方法反演横波速度模型，3次迭代后的结果如图4所示.从图中可以看出，相较于传统方法，SIRWDC方法对速度变化的刻画更为精细，反演得到的速度模型更逼近真实.但由于反问题存在多解性，在缺乏先验信息(例如地质认识)的情况下，反演结果不会随着迭代次数的增加而得到显著改善.

图4 SIRWDC与传统方法反演结果的对比(理想数据)Fig.4 Comparison of SIRWDC with traditional algorithm

3.2 含噪声的数据

理想情况下，频散分析应该完整地识别和提取每个模式的频散曲线.而实际中，利用频散能量成像提取频散曲线会遇到各种挑战，因此我们对理想频散曲线数据进行删减和加噪，使其更加接近真实情况.令提取的基阶频散曲线的频率范围在5～50 Hz，一阶高阶频散曲线的频率范围在40～100 Hz，二阶高阶频散曲线的频率范围在70～100 Hz，并且在频率范围10～30 Hz处添加噪声，噪声用来模拟基阶模式数据被高阶模式波或体波污染的情况，结果如图5所示.

图5 含噪声的瑞雷波多阶模式(基阶、一阶高阶和二阶高阶)频散曲线Fig.5 Multiple mode dispersion curves in noisy case (fundamental mode,first-order mode and second-order mode)

图6 SIRWDC与传统方法反演结果的对比(含噪声数据)Fig.6 Comparison of SIRWDC with the traditional method (noisy case)

4 结论

本文针对目前瑞雷波多阶模式频散曲线反演中仅考虑数据的拟合，缺乏对模型的约束，不能很好地刻画地层间断面的问题，研究了瑞雷波多阶模式频散曲线稀疏正则化反演方法，从正演模拟、反演建模和数值实现三个方面展开工作，提出了一套完整的瑞雷波多阶模式频散曲线反演方案：

(1)频散曲线正演基于经典的广义反射-透射系数法，采用一种快速求根方法进行瑞雷波相速度数值计算，与文献中常用的二分类方法相比，能够在很短的时间内达到最优的收敛效果.

(2)反演建模时考虑了地球物理模型的“逐块”光滑特性，由于这种物理上的间断可以通过模型参数在数学上的稀疏性进行描述，通过在目标函数中引入模型的L1范数对模型进行稀疏性刻画，使反演出的横波速度模型更加符合地质实际，并且能够增强反演算法对非高斯噪声的鲁棒性.

(3)在反问题的数值实现上，针对上述稀疏正化模型提出一种隐式迭代正则化算法，其迭代算子具有非膨胀特性，可以收敛到极小化问题的解，是一种结构简单、计算效率较高的反演计算方法.

数值实验结果表明，新的反演方案具有计算效率高，模型“逐块”光滑的特性刻画好、对非高斯噪声鲁棒性强的特点.在下一步的研究工作中，我们将针对实际资料进一步完善方案(包括地震记录中噪声影响的分析及压制，不同模态频散曲线的分离和提取，以及各阶模态加权系数的确定方法)，使其能够更好地应用于实际资料.

致谢感谢两位审稿人提出的非常具体的修改意见.