0 引言

瑞雷面波勘探具有能耗低、效率高、经济便捷、在浅地表分辨率高等优点，在近地表勘探、灾害地质调查、工程无损检测等方面获得了广泛的应用^[1-3]。

非均匀半空间层状介质的横波速度和地层厚度对瑞雷面波的频散曲线的结构和形态影响最大^[4-6]，对瑞雷面波频散曲线反演可获得一维横波速度曲线，能对地下地质结构进行分层。因此频散曲线反演是瑞雷波勘探中最重要的步骤之一^[7-8]。

瑞雷面波频散曲线反演是典型的高度非线性、多参数和多极值的地球物理优化问题。目前反演方法主要分为局部线性优化算法和非线性全局优化算法。局部线性算法包括最小二乘法、最光滑模型反演(OCCAM)算法等。该类方法计算速度快，但其主要缺点是反演结果的可靠性严重依赖初始模型，只有当初始模型接近真实模型时才能获得较好的反演结果，否则很容易陷入局部极值，得出错误的地质解释；同时该类方法还需要计算偏导数信息，雅克比矩阵的求取精度也将直接影响反演结果的可靠性。非线性全局优化算法有蚁群优化(Ant Colony Optimization，ACO)算法、遗传算法、粒子群优化(Particle Swarm Optimization，PSO)算法和模拟退火算法等，该类算法对初始模型要求较低，在全局搜索最优解，避免了局部线性算法存在的问题^[9-11]。

粒子群优化法是一种非线性算法，它通过简单的策略不断地更新位置和速度引导算法搜索到最优解。蔡伟等^[12]将该算法应用于瑞雷波频散曲线反演，结果表明PSO反演具有全局寻优能力强、快速、稳定且适用性强的特点，但是随着搜索次数的增加和反演参数的增加，群体易陷入平衡状态，粒子的进化过程随着时间的推移将会停滞不前。蚁群优化法同样是一种非线性算法，模拟了蚂蚁觅食的行为，通过不断释放信息素缩短觅食路径。王天琦等^[13]将该算法应用于瑞雷波频散曲线反演，结果表明ACO法在求解单峰函数时收敛速度较快、局部寻优能力强、解的精度高，但面对多峰函数却经常发生收敛早熟，无法求出全局最优值。

本文综合考虑了PSO算法和ACO算法各自的优、缺点，在Shelokar等^[14]的研究基础上提出了利用粒子群与蚁群相结合的优化(PSACO)算法进行瑞雷波频散曲线反演。该反演策略充分发挥PSO算法较强的全局寻优能力，在可行解空间进行大范围全局搜索，有效避免了算法陷入局部极小的可能性；又充分利用ACO算法不断更新PSO算法中粒子的位置，粒子的进化不会停滞，迭代后期具有较强的局部寻优能力，能提高解的精度，加快后期收敛速度，从而有效克服了PSO和ACO算法的各自缺点。理论模型的测试结果表明，与单独的PSO算法和ACO算法相比，PSACO算法更优越、更有效；实测资料反演结果验证了PSACO算法的实用性。

1 基本原理及反演流程 1.1 PSO算法

在PSO反演中，种群的候选解称为粒子，它与相邻粒子经验共享，同时共存和进化。每个粒子在问题的搜索空间中飞行时(每个粒子为指定空间的一个多维向量，代表一个可行性解)，通过使用自己的飞行经验(早期飞行中找到的最佳位置的记忆)和邻近粒子的飞行经验修改其速度以找到更好的解决方案(最优解)。粒子的位置和速度更新可表示为

$ \boldsymbol{x}_{i}^{(m+1)}=\boldsymbol{x}_{i}^{(m)}+\boldsymbol{v}_{i}^{(m+1)} $

(1)

$ \begin{gathered} \boldsymbol{v}_{i}^{(m+1)}=w_{m} \boldsymbol{v}_{i}^{(m)}+c_{1} r_{1}\left[\boldsymbol{P}_{i}^{(m)}-\boldsymbol{x}_{i}^{(m)}\right]+ \\ c_{2} r_{2}\left[\boldsymbol{P}_{\mathrm{g}}^{(m)}-\boldsymbol{x}_{i}^{(m)}\right] \end{gathered}$

(2)

式中：x_i^(m)为第i个粒子、第m次搜索时在解空间中的位置；P_i^(m)为第i个粒子、前m次搜索时的最佳位置；v_i^(m)表示第i个粒子、第m次搜索时的速度，其值控制在[-v_max, v_max]，避免粒子在搜索空间外的过度漫游；P_g^(m)为第m次搜索时所有粒子中的全局最优解；r₁、r₂为在(0，1)内均匀分布的随机数；c₁、c₂为加速系数，用于调整算法的收敛策略，加速对全局最优值的搜索；w_m为惯性权重系数。为加速对全局最优值的搜索，应尽量满足3.0≤c₁+c₂≤4.0；如果c₁与c₂设置不当，则容易使算法早熟收敛^[12]。w_m动态变化，前期较大，让算法对全局最优解敏感，也让算法不易陷入局部最优；在后期较小，以减少搜索面积，加强局部搜索能力，以提高解的精度和缩短搜索时间。w_m定义为

$ w_{m}= \begin{cases}-\left(\frac{m}{M}\right)^{2}+w_{\mathrm{s}} & 1 \leqslant m \leqslant \frac{M}{2} \\ \left(\frac{m}{M}-1\right)^{2}+w_{\mathrm{e}} & \frac{M}{2}<m \leqslant M\end{cases} $

(3)

式中：M为最大迭代次数；w_s、w_e分别为初始和终止惯性权重系数。

1.2 ACO算法

ACO算法模拟蚂蚁的觅食行为。当蚂蚁从食物源到巢穴的往返过程中，它们会在地上释放信息素，并形成信息素轨迹；蚂蚁能感知信息素的气味和浓度，当它们在选择路径时，会选择信息素浓度强(较短)的路径。该算法可以解决复杂的组合优化问题，例如旅行商(TSP)问题和二次分配问题。以经典的TSP问题为例，蚂蚁算子的移动搜索策略遵循

$ \left\{\begin{array}{l} A_{i j}^{k}=\frac{\left[\tau_{i j}(t)\right]^{\alpha}\left[\eta_{i j}(t)\right]^{\beta}}{\sum\limits_{s}\left[\tau_{i s}(t)\right]^{\alpha}\left[\eta_{i s}(t)\right]^{\beta}} \\ \eta_{i j}(t)=\frac{1}{D_{i j}} \end{array}\right. $

(4)

式中：i、j表示两个目标城市；s表示城市i与城市j间的蚂蚁k未到达的地点；τ_ij(t)表示t时刻位于路径(i，j)上的信息素；η_ij(t)为启发函数，为距离的倒数，两个城市之间的距离D_ij越短启发函数的值越大；指数α以及β分别表示影响因子，影响着算法的对解的寻优能力；A_ij^k表示第k只蚂蚁根据信息素浓度和启发函数的大小由城市i选择到城市j的概率。当所有蚂蚁到达节点后，信息素将进行更新

$ \tau_{i j}(t+1)=(1-\rho) \tau_{i j}(t)+\Delta \tau_{i j}(t) $

(5)

式中：ρ为信息素挥发系数，经验取值一般约为0.9；Δτ_ij(t)为本次循环路径(i，j)上信息素的增量，表示每只蚂蚁算子在本次循环留在路径上信息素含量之和。

1.3 PSACO算法的原理

在一次总迭代过程中交替运行PSO算法和ACO算法实现PSACO算法。ACO算法是一种局部搜索算法，蚂蚁利用信息素引导机制更新粒子在早期发现的位置。在PSACO中，一种简单的信息素引导的ACO机制应用于局部最优值的搜索。ACO算法中生成的蚂蚁个数等于PSO算法中粒子的个数。每个蚂蚁根据下式围绕全局最佳位置生成对应的粒子

$ \boldsymbol{z}_{i}^{(m)} \sim N\left[\boldsymbol{P}_{\mathrm{g}}^{(m)}, \sigma_{m}\right] $

(6)

式(6)表示生成的解向量z_i^(m)满足以P_g^(m)为期望、以σ_m为标准差的高斯分布。当m=1时，σ₁=1，并在每次搜索结束时更新σ_m+1=σ_m×d，其中参数d满足0.250≤d≤0.997。令10^-4≤σ_min≤10^-2，如果σ_m < σ_min，则σ_m=σ_min。再计算粒子z_i^(m)的适应度，通过比较蚁群和粒子群适应度f的大小进行位置更换^[14]。若f [z_i^(m)] < f [x_i^(m)]，则用z_i^(m)替换x_i^(m)，f [z_i^(m)]替换f[x_i^(m)]。这种简单的信息素引导机制能够快速地找到并锁定全局最优值。起初蚂蚁在标准差σ_m较大时，搜索范围更广，可以避免搜索陷入局部最优值，随着搜索次数的增加和标准差的逐渐减小，让所有蚂蚁都更加接近全局最优值，加快最优值的搜索。因此，ACO不仅有助于PSO算法有效地搜索全局以快速获得可行的解空间，而且还有助于有效接近或达到最优解。

1.4 混合优化算法瑞雷面波频散曲线反演步骤

瑞雷面波频散曲线反演实质是一个对目标函数求全局最优解的过程。衡量反演结果优劣最直接的因素是反演的目标函数大小。本文瑞雷面波频散曲线反演的目标函数(适应度)定义为

$ f=\sqrt{\frac{\sum\limits_{n=1}^{J}\left(V_{\mathrm{R}, n}^{\mathrm{obs}}-V_{\mathrm{R}, n}^{\mathrm{al}}\right)^{2}}{J}} $

(7)

式中：J为待拟合频散曲线的频点数；V_R^obs、V_R^cal分别为待拟合频散数据和反演的地质模型在对应频率下正演的相速度。

混合算法中的粒子个数应与蚁群个数保持相同，一般解向量的维度越低，设置的种群数较少，反之则更多。针对不同维度的问题一般设置种群的数量为2L(2L-1)，其中L表示地层层数。通过该经验公式可更好地平衡反演精度与反演用时。

PSACO算法实现流程如下：

(1) 初始化。设置粒子个数I、最大迭代次数M等参数；

(2) 拾取待拟合频散曲线信息并确定搜索边界，在搜索范围内随机产生I个粒子；

(3) 计算每个解的目标函数，并令最小目标函数值对应的粒子为全局最优解，单个粒子为局部最优解；

(4) 利用粒子群优化算法更新粒子并计算适应度，若优于局部最优解和全局最优解，则进行粒子替换；

(5) 利用蚁群算法在全局最优解附近产生I个蚂蚁(可能解)，并计算每个解的适应度；

(6) 比较蚂蚁与粒子的适应度，若蚂蚁解更优，则用蚂蚁替代粒子；

(7) 若全局最优解的适应度大于精度要求或没有达到最大迭代次数，则返回步骤(4)；

(8) 输出全局最优解。

2 理论模型试算

应用PSACO算法对瑞雷面波勘探中常遇到的地质结构进行反演，分析不同模型的反演结果，评估算法的应用效果。

首先，设置每个模型的层数、每层的纵横波速度(V_P、V_S)、密度以及层厚度(H)；然后，用快速矢量算法正演模拟瑞雷波频散曲线(搜索步长为2，频带范围为0~100Hz)；最后，用PSACO算法对频散数据进行反演。通过分析目标函数的大小、各层横波速度和厚度随搜索次数的变化过程等，探究PSACO算法的实用性和稳定性。

2.1 三层地质模型

为研究算法的适用性，建立三个典型的三层地质模型：速度递增型、含低速软夹层型、含高速硬夹层型^[15]。三个模型都采用相同的算法参数：初始惯性权重系数w_s=0.7，终止惯性权重系数w_e=0.4，加速系数c₁=2、c₂=2，初始标准差为σ₁=1，控制标准差变化的参数d=0.9。

2.1.1 模型1

模型1为三层速度递增型，相关参数以及搜索范围如表 1所示。利用PSACO算法对模型1的频散曲线进行反演，共计反演20次，每次反演迭代100次(每个粒子都进行100次位置更新)，去除两次与平均值相差最大的反演结果后，剩余结果取平均值(图 1)。

图 1a为目标函数随迭代次数变化曲线，可以看出，前15次目标函数下降较快，到第70次时几乎下降为0。由图 1b~图 1f可知：当搜索到100次时，反演的每个参数与真值基本一致；收敛曲线在60次后仅有极小幅度的波动。

由图 2可知，反演的频散曲线和理论频散线拟合很好，反演的横波速度和层厚度与真实模型高度吻合。表 2为反演20次的结果统计。

由表 2可知，模型1经过多次反演，PSACO算法所得结果的相对误差都接近于0，标准差同样很小，验证了本文算法精度高、稳定性好。

2.1.2 模型2

用PSACO算法对三层含低速软夹层地质模型模拟的频散曲线进行反演，其模型的相关参数以及搜索范围如表 3所示。共计反演20次，每次反演迭代100次，去除两次与平均值相差最大的反演结果，剩余结果取平均如图 3所示。

图 3a为目标函数随迭代次数的变化曲线，可见：前15次目标函数下降较快，到第80次时几乎下降到零；与模型1相比，搜索难度有所增加。由图 3b~图 3f可以看出，V_S1、V_S3收敛结果与真值非常接近，但V_S2、H₁、H₂收敛精度有所下降，稳定性也有所降低，但仍与真值接近。

模型2的最终反演结果如图 4所示，反演的横波速度剖面和理论模型几乎重合。PSACO算法20次反演的统计结果如表 4所示。

由表 4可知，经过多次反演，模型2的PSACO算法各参数反演的相对误差以及标准差均接近于0，说明PSACO算法对三层含低速软夹层模型的反演结果精度高、稳定性好。

2.1.3 模型3

模型3为三层含高速硬夹层型，相关参数和参数的搜素范围如表 5所示。利用PSACO算法对频散曲线反演共计20次，每次反演迭代100次，去除两次与平均值相差最大的反演结果后，剩余结果取平均，如图 5所示。

由图 5a目标函数随搜索次数变化曲线可以看出，前15次目标函数下降较快，到50次左右时趋近于0。由图 5b~图 5f可见，每个参数反演都能在搜索100次时与真值重合，收敛曲线在后期未发生明显波动。

模型3反演的最终结果如图 6所示，反演的横波速度剖面和理论模型几乎重合。统计PSACO算法模型3反演20次的结果，如表 6所示。

由表 6可知，经过多次反演，PSACO算法的各参数反演均值与理论值十分接近，相对误差较低，标准差控制在合理范围内。

三个典型模型数据反演结果均表明，PSACO算法精度高、稳定性好。

2.2 PSACO与单独的PSO、ACO算法对比

四层含低速软夹层地质模型(表 7)参数较多，适合用来对比不同算法的稳定性、准确性。为有效对比，PSACO、ACO、PSO算法采用的反演参数保持一致：初始惯性权重系数w_s=0.7，终止惯性权重系数w_e=0.4；加速系数c₁=2、c₂=2；初始标准差为1，控制σ变化的参数d=0.9；影响因子α=0.4，β=0.4；挥发系数ρ=0.9。

图 7a为三种方法目标函数随迭代次数的变化曲线对比，可见：PSACO算法能在目标函数快速下降的同时也能保证算法不陷入局部最优；PSO算法虽然在整个迭代过程中保持下降的趋势，但是后期的搜索能力不强，解的精度不高；ACO算法在前期收敛速度快，但很明显在后期出现收敛停滞，解的误差偏大。由图 7b~图 7h可知，与PSO和ACO算法相比，PSACO算法搜索到100次时的结果更接近真值，在搜索的过程中稳定性更高，表明PSACO算法具有更高的稳定性和更强适用性。

图 8为四层模型三种方法反演结果的对比，正演频散曲线和反演剖面基本相同，都接近真值。表 8为四层模型三种方法反演结果的统计，可见，PSACO算法的7个反演参数的相对误差和标准差都低于PSO和ACO算法。相对误差低，说明该算法的全局寻优能力强；标准差小，说明该算法的稳定性更高。

3 实例分析

通过对意大利北部某垃圾填埋场的实例反演试算，检验PSACO算法对实测数据的反演能力，验证该算法的实用性。

3.1 实例背景介绍与参数取值

实例数据来自意大利东北边的垃圾填埋场^[16]，该区域的基岩层是灰岩，上层由未胶结的堆积物构成，约18m厚。由该地点的钻孔资料可知每层的厚度及其分层情况。

同前文对理论模型反演的方法类似，泊松比和密度视为已知(通过已掌握的实测资料进行估算^[15])，横波速度和厚度为待反演参数，初始惯性权重系数w_s=0.7，终止惯性权重系数w_e=0.4；加速系数c₁=2、c₂=2；初始标准差为1，控制σ变化的参数d=0.9。模型参数及参数搜索范围如表 9所示。

3.2 反演效果分析

由PSACO、ACO、PSO反演算法目标函数随迭代次数变化曲线(图 9a)可见：PSACO算法能在目标函数快速下降的同时也能保证算法不陷入局部最优；PSO算法在后期的搜索能力不强，解的精度不高；ACO算法在前期收敛速度快，但很在后期出现收敛停滞，解的误差偏大。由于地层厚度和介质速度是未知的，无法直观对比三种算法反演参数V_S、H与实测值的近似程度，由图 9b~图 9j曲线的稳定性可知，PSACO算法稳定性最高。

由频散曲线对比(图 10a)可知，反演的频散曲线与实测频散点拟合度高。通过将不同算法反演的横波速度剖面图(图 10b)和钻孔资料(图 10c)对比可以发现，PSACO算法与钻孔资料的匹配度更高，反演的横波速度剖面更接近实际地层。

4 结束语

本文通过理论模型试算、不同算法之间的对比、实测资料分析验证了PSACO算法应用于瑞雷面波频散曲线非线性反演的可行性、适用性、稳定性。PSACO算法具有优秀的全局寻优能力，在迭代后期局部寻优能力也较强；该算法的成功率、精度、稳定性在总体上都优于PSO算法和ACO算法；PSACO算法应用至实测资料的瑞雷面波频散曲线反演结果较其他算法更接近真值。