0 引言

瑞雷面波勘探可获得地下空间横波速度结构，已广泛应用于城市地下空间开发、地质灾害调查等领域^[1-6]。瑞雷面波频散曲线反演是面波勘探的重要步骤，是多极值、多参数、非线性、多模态的迭代优化过程^[7]，主要包括线性局部优化算法和非线性全局优化算法两大类。常见的线性局部优化算法有阻尼最小二乘法^[8-9]、奇异值分解算法^[10]及Occam算法^[11-12]等。线性局部优化算法对初始模型的选择以及对偏导数的计算精度要求均较高，限制了它的应用与发展。非线性全局优化算法有遗传算法^[13-16]、模拟退火算法^[17-18]、粒子群算法^[19-21]、麻雀搜索算法^[22]、蚁狮优化算法^[23]、海洋捕食者优化算法^[24]等。非线性全局优化算法因结构简单、计算效率高、不依赖建立初始模型，能更好地获得全局最优解而备受重视。粒子群算法具有收敛速度快、涉及参数少、编程易于实现的优点，广泛应用于瑞雷面波频散曲线反演中。Song等^[25]2012年首次将粒子群算法应用于基阶瑞雷面波频散曲线反演；其后两年，杨钊^[26]采用频散函数的模作为目标函数，将粒子群算法首次应用到多模态瑞雷面波频散曲线反演，并验证了该算法具有一定的抗噪性；张晓阳等^[27]结合细化分层理论与粒子群算法反演瑞雷面波频散曲线，取得了较好效果；刘隼等^[28]提出了基于混沌状态下的粒子群算法反演瑞雷面波频散曲线的方法；蔡伟等^[29]利用压缩因子和自组织的惯性权重优化策略对粒子群算法进行了改进，实现了快速、稳定的瑞雷面波频散曲线反演；符健^[30]提出基于正切函数的惯性权重和加入遗传算法的交叉变异操作两种改进策略对粒子群算法反演瑞雷面波频散曲线进行了完善，取得了较好的效果；杨博等^[31]采用自适应惯性权重、设置粒子的节速度、引入遗传算法的交叉和变异操作、单维全分量的混沌局部搜索等四种优化策略对粒子群算法进行改进，取得了较好的瑞雷面波频散曲线反演结果；王一鸣等^[32]提出将粒子群优化算法与蚁群优化算法相结合的优化算法(PSO-ACO)反演瑞雷面波频散曲线，充分发挥了粒子群全局寻优能力和蚁群算法的局部搜索能力。

传统粒子群算法容易早熟、收敛，陷入局部极值。针对粒子群算法存在的问题，学者们提出了诸多优化策略对其进行改进，其中有非线性的自适应惯性权重方法^[33]、添加压缩因子^[34]、设置粒子边界条件约束^[35]等优化策略。非线性的自适应惯性权重策略在前期具有较大的惯性权重w，具备良好的全局搜索能力；在算法后期较小的w具备良好的局部搜索能力，可以提高反演精度。添加压缩因子有助于提高算法的收敛速度和精度。设置粒子的边界条件有助于提高算法的多样性，提高寻优解的准确性。本文将这三种优化策略进行融合，并将其应用于瑞雷面波频散曲线反演，提高了反演的精度。

为分析杂交粒子群优化(Hybrid Particle Swarm Optimization, HPSO)算法反演瑞雷面波频散曲线的效果与优势，建立了典型地质模型，进行无噪声和含噪声的基阶理论频散曲线反演、多模态理论频散曲线反演，并将其与边界约束粒子群算法的反演结果进行对比，结果表明本文方法具有更高的反演精度和抗噪性能。最后，本文方法在实测面波勘探数据的应用证明了其实用性。

1 杂交粒子群优化算法 1.1 粒子群算法的基本原理

Kennedy等^[36]在观察鸟群觅食的社会行为基础上，提出一种进化算法——粒子群算法。该算法将鸟类捕食过程看作是粒子在一定范围内的寻优过程。粒子具备速度和位置两大属性，优化问题的解就是粒子本身的位置。假设粒子i在m维度的寻优空间表示为x_i=(x_{i, 1}, x_{i, 2}, …，x_{i, m})，每个粒子位置的优劣由目标函数表示的适应度函数进行评估，粒子的飞行速度可表示为v_i=(v_{i, 1}, v_{i, 2}, …，v_{i, m})，每个粒子当前的最优位置可表示为p_i=(p_{i, 1}, p_{i, 2}, …，p_{i, m})，当前群体中粒子g的最优位置可表示为p_g=(p_{g, 1}, p_{g, 2}, …，p_{g, m})。粒子的速度和位置可用下式迭代更新

$ \left\{\begin{aligned} v_{i, m}^{(t+1)}= & {wv}_{i, m}^{(t)}+c_1 \operatorname{rand}\left(p_{i, m}-x_{i, m}^{(t)}\right)+ \\ & c_2 \operatorname{rand}\left(p_{g, m}-x_{i, m}^{(t)}\right) \\ x_{i, m}^{(t+1)}= & x_{i, m}^{(t)}+v_{i, m}^{(t)} \end{aligned}\right. $

(1)

式中：粒子i=1, 2, …，N，其中N为粒子的总个数；t为迭代次数；c₁、c₂分别表示个体学习因子与群体学习因子，分别调节向个体最好粒子和全局最好粒子方向飞行的最大步长；rand(·)表示取0到1之间的随机数。

为进一步优化算法，设置粒子的边界条件，当粒子超出搜索范围，需将逃逸出的粒子在一定的搜索范围随机分布，以保证种群的多样性^[35]

$ \left\{ {\begin{array}{*{20}{l}} {{x_{i, m}} = {x_{{\rm{ub}}}} - B{\rm{rand}}\left( {{x_{{\rm{ub}}}} - {x_{{\rm{lb}}}}} \right)}&{{x_{i, m}} > {x_{{\rm{ub}}}}}\\ {{x_{i, m}} = {x_{{\rm{lb}}}} + B{\rm{rand}}\left( {{x_{{\rm{ub}}}} - {x_{{\rm{lb}}}}} \right)}&{{x_{i, m}} < {x_{{\rm{lb}}}}} \end{array}} \right. $

(2)

式中：B∈[0, 1]；x_ub、x_lb分别为搜索范围的上、下边界。

1.2 HPSO的改进策略

粒子群优化算法是一种传统的自然启发式优化算法，涉及的参数简单，不依赖于初始模型，编程易于实现，且粒子具有记忆功能，在迭代优化的过程中，粒子之间相互传递信息，并且利用粒子的速度迅速寻优，收敛速度较快，全局寻优能力强。尽管粒子群算法存在诸多优点，但其缺点同样不容小觑。虽然粒子群算法在寻优早期，粒子都朝着最优位置方向前进，收敛速度较快，但后期粒子群趋于稳定，收敛速度慢，种群缺乏多样性，这容易导致算法收敛、早熟，陷入局部极值，收敛精度不高。此外，粒子群算法最终的反演结果质量易受种群大小及其算法参数大小的影响，反演结果不稳定^[26-31]。

目前对粒子群算法的改进主要表现在几个方面：①惯性权重的改进，主要包括线性递减的权重方法、自适应调整的惯性权重方法等，惯性权重的改进使粒子群算法在算法前期对全局进行充分探索，后期也能够满足局部搜索的需求^[33]；②引入压缩因子，改进粒子更新速度更新方式，提高算法的收敛速度与收敛精度^[34]；③采用变学习因子，学习因子的变化以采用同步变化或者异步变化进行调整，学习因子的调整同样有助平衡算法的全局搜索和局部搜索能力^[37]；④边界条件的处理，可将粒子的位置限制在搜索范围，提高种群的多样性^[34]；⑤将粒子群算法与遗传算法、模拟退火等其他智能优化算法相结合，取长补短，形成混合粒子群算法^[38-40]。

本文将边界条件约束的粒子群改进方法应用于面波频散曲线反演研究。此外，将其与前人研究的惯性权重系数和引入压缩因子两种改进策略融合，进一步优化粒子群算法反演面波频散曲线的过程。

具体优化过程如下：

首先，在粒子群算法的改进中采用了动态惯性权重系数^[30]，其表达式为

$ w= \begin{cases}w_{\min }+\frac{\left(w_{\max }-w_{\min }\right)\left(f-f_{\min }\right)}{f_{\mathrm{avg}}-f_{\min }} & f \leq f_{\mathrm{avg}} \\ w_{\max } & f>f_{\mathrm{avg}}\end{cases} $

(3)

式中：w_min、w_max分别表示最小和最大惯性权重，本文取值分别为0.6和0.9^[28]；f、f_avg、f_min分别表示当前迭代下的目标函数值、平均目标函数值及最小目标函数值。较大的惯性权重有助于进行全局搜索，对目标函数大于平均目标函数值的粒子，其惯性权重增大，有利于跳出局部极值。粒子的目标函数一旦小于平均目标函数值，当前惯性权重动态调整为较小的惯性权重，从而增强局部搜索能力，有助于算法收敛到最优值。动态的惯性权重调整在一定程度上平衡了全局和局部的搜索能力。

其次，在采用动态惯性权重的基础上引入压缩因子，对粒子的速度向量进行更新，在有效控制粒子速度的同时，增强算法在局部空间搜索的能力^[34]。则粒子更新表示为

$ \left\{\begin{aligned} & v_{i, m}^{(t+1)}= \gamma\left[{wv}_{i, m}^{(t)}+c_1 \operatorname{rand}\left(p_{i, m}-x_{i, m}^{(t)}\right)+\right. \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\left.c_2 \operatorname{rand}\left(p_{g, m}-x_{i, m}^{(t)}\right)\right] \\ & x_{i, m}^{(t+1)}= x_{i, m}^{(t)}+v_{i, m}^{(t)} \\ & \gamma=\frac{2}{\left|2-\left(c_1+c_2\right)-\sqrt{\left(c_1+c_2\right)^2-4\left(c_1+c_2\right)}\right|} \end{aligned}\right. $

(4)

式中γ为压缩因子，可以加快收敛速度，使粒子不同时期在有效的搜索范围内进行搜索，得到高质量的解。

最后，在前两种优化策略基础上，设置粒子的边界条件。当粒子的位置范围超过设定的搜索范围时，利用式(2)重新将逃逸出的粒子随机分布在设定的边界范围内，提高种群的多样性。

2 理论模型试算

建立三层速度递增的地质模型(模型A)和三层含低速夹层的地质模型(模型B)，通过基阶理论频散曲线反演各模型的横波速度V_S、纵波速度V_P和层厚h，验证本文算法的有效性。以V_S和h反演参数赋值的35%作为搜索的上、下限。V_S、V_P以及h的具体参数见表 1。为防止单独一次反演算法结果出现的随机性，取10次独立反演的平均值作为最终反演结果。反演迭代次数为20次，种群大小为40。采用快速矢量传递算法计算模型的理论频散曲线，频散曲线反演的目标函数为

$ \phi=\frac{\left\|V_{\mathrm{R}}^{\text {obs }}-V_{\mathrm{R}}^{\text {theo }}\right\|}{\sqrt{M}} $

(5)

式中：V_R^obs、V_R^theo分别为观测频散曲线数据及理论频散曲线数据；M为瑞雷面波的频散点数。

2.1 无噪声基阶频散曲线反演

瑞雷面波具有多模态特征，但大量实测数据表明基阶面波占据主要能量，且分布较为广泛，最容易观测得到。通过模型A、模型B的基阶频散曲线的反演验证本文方法的有效性。频率范围选择5~70 Hz，反演结果如图 1所示，可以发现模型B的每一次独立反演结果相对模型A波动较大，这可能是地层速度倒转，横波速度剧烈变化，导致反演结果不稳定。尽管反演难度加大，但模型B的反演频散曲线相对均匀分布在理论频散曲线的上、下两侧，反演模型与已知模型基本一致，模型A、模型B的反演参数平均误差分别为1.88%、3.10%(表 2)，表明HPSO可以有效进行无噪声基阶频散曲线反演。

2.2 含噪声频散曲线反演

为进一步说明本文算法的适应性，采用加入噪声的频散曲线反演测试。加噪公式可以表示为

$ d_n=d+2(0.5-B) \mathrm{d} n $

(6)

式中：d_n表示加噪后的频散曲线数据；d为无噪声的理论频散曲线数据；n表征噪声的程度。

在模型A、模型B的基阶频散曲线中加入10%的随机噪声，HPSO反演结果如图 2和表 2所示。模型A、B不含噪声数据横波速度反演的平均相对误差分别为0.61%、1.67%，层厚反演平均相对误差分别为3.80%、5.25%；模型含10%噪声数据横波速度反演的平均相对误差分别为1.21%、1.69%，层厚反演平均相对误差分别为2.25%、7.63%。通过对比可以发现，频散曲线对横波速度的反演结果明显优于对层厚的反演结果。此外，加入噪声数据后，横波速度和层厚反演整体反演结果误差有所增加，但变化不大，最终反演结果依旧与理论频散曲线拟合较好，表明HPSO具有很强的抗噪能力。

2.3 多阶频散曲线的联合反演

在常见的速度递增地质模型中，实测瑞雷面波频散曲线常以基阶为主，但当地质模型为含低速夹层等复杂模型，实测面波频散曲线往往包含高阶模式频散曲线。为了反演得到更加准确的地层信息，多模式频散曲线反演是非常有必要的。利用HPSO进行模型A、模型B的基阶与二阶频散曲线联合反演，验证本文方法多阶面波频散曲线的反演效果。反演得到的频散曲线和横波速度剖面如图 3所示，反演的模型参数值如表 3所示。模型A多模态频散曲线联合反演的横波速度以及层厚的平均相对误差分别为0.36%、2.45%，而对应基阶反演的横波速度、层厚平均相对误差分别为0.61%、3.80%；模型B多模态频散曲线联合反演的横波速度以及层厚的平均相对误差分别为0.92%、3.75%，对应基阶反演的横波速度、层厚平均相对误差分别为1.67%、5.25%。反演结果表明，无论是三层速度递增模型还是含低速夹层模型，多模态频散曲线联合反演的横波速度和层厚平均相对误差均明显低于仅基阶反演的平均相对误差，多模态频散曲线联合反演结果与已知地质模型高度一致，各阶面波频散曲线拟合度都较高(图 3)。上述结果表明HPSO同样适用于多模态频散曲线联合反演。

2.4 与边界约束粒子群算法的对比分析

为证明HPSO的优越性，设置相同反演种群大小(N=40)、迭代次数(20次)及独立反演次数(10次), 将HPSO与边界约束粒子群算法反演结果进行对比分析，边界约束粒子群算法反演结果如图 4和表 4所示。反演结果表明：边界约束的粒子群算法每次独立反演结果的稳定性较差，在模型A、模型B中对横波速度反演平均标准差分别为9.07%、14.95%，而平均相对误差分别为1.03%、1.83%；HPSO对模型A、模型B的横波速度反演平均标准差分别为7.89%、12.88%，平均相对误差分别为0.61%、1.67%。从平均标准差以及平均相对误差的结果来看，HPSO稳定性相对较好，反演结果精度也相对更高。

3 实测资料分析

为了检验HPSO的实用性，利用在冰岛南部Ölfusá河岸的Arnarbæli周边试验场地A2采集到的多道面波数据(图 5a)进行分析。共24道接收，主频为4.5 Hz的检波器，道间距为1 m，炮检距为10 m，采样频率为1000 Hz，震源采用6.3 kg的大锤^[41]。利用倾斜叠加方法得到频谱能量图(图 5b)拾取频散曲线。从拾取的频散曲线可以看出，在5 Hz前及42 Hz后拾取的值偏离基阶模式频散曲线较远，误差较大，因此本文选取5~42 Hz段内的频散曲线作为实测基阶频散曲线进行反演。参照文献[30]频散曲线反演，地层层数设置为10，纵波速度均设置为1440 m/s，密度均设置为1.85 g/cm³，横波速度和厚度的搜索范围如表 5所示。反演结果表明，HPSO反演结果的频散曲线与实测频散曲线拟合更好，其均方差误差为2.01 m/s，远小于参考文献[41]中的5.18 m/s。本文方法的反演结果与文献[41]反演结果在5 m以内基本一致，5 m以下横波速度存在一定偏差，但横波速度结构的整体趋势基本保持一致，证明了HPSO的实用性。鉴于本文方法反演结果的频散拟合程度更高，可以认为本文反演的横波速度结果具有更高的精度。

4 结束语

因为粒子群算法反演过程容易陷入局部极值，导致算法早熟、收敛，本文综合非线性的自适应惯性权重、引入压缩因子以及边界条件约束的三种优化策略对粒子群算法进行了改进，建立了杂交粒子群优化算法(HPSO)。利用HPSO对速度递增模型与含低速夹层模型的无噪声、含噪声和多阶模式的理论频散曲线进行了反演，验证了该方法的有效性、抗噪能力和多模态频散曲线反演能力。通过与边界约束的粒子群算法反演结果进行对比，证明了HPSO比边界约束的粒子群算法具有更高的反演精度，凸显了HPSO的优越性。最后，通过冰岛Arnarbæli地区的实测数据的反演，得到与前人研究基本一致的反演结果和更高的频散曲线拟合度，检验了该算法的实用性。