0 引言

地球物理单一方法的局限性、反演的多解性以及地质条件的复杂性使得联合反演成为地球物理领域的必然发展趋势^[1]。联合反演是一种定量的综合地球物理解释技术，可以减少解的非唯一性。随着地球物理勘探向地球深部发展，联合反演的优势也越加突显^[2-4]。基于线性优化算法的联合反演收敛速度快，但容易陷入局部极小。基于非线性优化算法的联合反演具有全局寻优、无需求偏导数矩阵、便于先验信息融入等特点。目前，很多非线性优化算法在地球物理联合反演领域得到了广泛应用，如模拟退火算法^[5-7]、遗传算法^[8-9]、差分进化算法^[10-11]、人工鱼群算法^[12]等。可见基于非线性优化算法的联合反演是地球物理联合反演领域的一个重要发展方向。

人工蜂群(ABC)算法是Karaboga^[13]于2005年提出的一种非线性优化算法。在该算法中，存在引领蜂、跟随蜂和侦察蜂，每种蜜蜂各司其职，并存在独特的角色转换机制。由于该算法是一种较新颖的优化算法，相比传统的非线性优化算法(如模拟退火、遗传算法等)，该算法在地球物理反演领域的应用相对偏少。王猛等^[14]将ABC算法应用于瞬变电磁测深资料反演，并指出ABC算法可以提高瞬变电磁资料的解释精度；候征等^[15]将ABC算法用于一维瑞雷波多阶模式联合反演，研究结果表明该算法可提高反演精度；Wen等^[16]将该算法应用于某煤矿2.5维CSAMT数据反演，揭示出疑似陷落柱的地电结构。需要指出的是，ABC算法虽然具有较好的探索能力，但局部搜索能力较弱，收敛相对较慢^[17-18]。引领蜂在蜜源附近搜索时，只随机地对解的某一个分量进行更新，当待解参数较多时，这样搜索方式效率偏低。此外，优质蜜源对蜂群影响不足。

在非线性优化算法研究领域，双种群架构是增强算法寻优能力的重要思路。赵燕伟等^[19]实现了基于双种群的遗传算法，指出改进算法可以提高算法的全局收敛性能；吴亮红等^[20]提出一种双种群差分进化(DE)算法，不同的种群使用不同的变异策略，仿真实验表明该算法全局搜索能力强、收敛速度快；暴励等^[21]采用双种群架构使ABC算法与DE算法并行进化，每隔一定进化次数就分别比较两个种群的最优解个体，用最优个体替代次优个体，以促进种群间的相互学习；何光杰等^[22]提出了一种双种群的准粒子群算法，两个种群分别使用不同的粒子进化方式，增强了算法的寻优能力；陈亚峰等^[23]借助双种群对萤火虫算法的种群进行了全局和局部种群划分，一定程度上解决了萤火虫算法全局搜索和局部搜索不能兼顾、快速收敛和过早停滞相矛盾的问题。

基于前人的研究成果，本文提出一种基于双种群架构的ABC算法，并以大地电磁测深(MT)和重力联合反演为例，验证该算法的适用性和实用性。

1 基于双种群架构的ABC算法 1.1 方法原理

标准ABC算法的基本流程可以概括为以下四个方面：蜜源的生成、引领蜂寻找更优质蜜源、跟随蜂决定是否跟随及判断是否达到局部搜索限制次数^[24]。在标准ABC算法的基础上，本文提出基于双种群架构的改进ABC算法：将种群一分为二，在种群1中，尝试引入交叉和变异策略以提高ABC算法的搜索效率；在种群2中，尝试引入最优解邻域搜索以提高ABC算法的局部寻优能力。两个种群每经历一定代数后就进行互相交流。

1.1.1 变异和交叉操作

在标准ABC算法生成蜜源阶段，只对模型分量进行随机的一维搜索，即

$ U_{i, d}^{1 \mathrm{a}}=X_{i, d}^{t}+\alpha\left(X_{i, d}^{t}-X_{j, d}^{t}\right) $

(1)

式中：$ U_{i, d}^{1{\rm{a}}} $表示种群1标准ABC算法搜索更新的蜜源，其中i表示蜜源个体编号，d表示蜜源的第d个分量；α表示[-1, 1]的随机数；$ X_{i, d}^t $和$ X_{j, d}^t $分别表示第t次迭代第i个和第j个蜜源的第d个分量。该搜索策略效率偏低，本文借鉴DE和遗传算法里的搜索策略，将变异和交叉操作引入种群1。

在标准ABC算法搜索策略更新的蜜源U^1a劣于上代蜜源时，按照下式进行搜索以提高寻优效率

$ U_{i, d}^{1 \mathrm{~b}}= \begin{cases}X_{i, d}^{t}+F\left(X_{j, d}^{t}-X_{i, d}^{t}\right) & \operatorname{rand}(0, 1)<\mathrm{CR} \\ X_{i, d}^{t} & \text { 其他 }\end{cases} $

(2)

式中：U^1b是根据变异和交叉操作搜索出的新蜜源；F和CR分别表示变异因子和交叉因子，取值范围均为[0, 1];rand (·)表示取随机数。

1.1.2 最优解邻域搜索

上节介绍的变异和交叉操作加快了对个体的开发，但是优质蜜源在蜂群中仍然无法传播。故此，本文在种群2中加入最优解邻域搜索策略，利用下式对优质蜜源进行充分开采

$ U_{i, d}^{2}=X_{\text {best }}+\left(X_{\text {best }}-X_{i, d}\right) $

(3)

式中：U²是种群2更新的蜜源；X_best表示种群2中的最优个体，该参数的引入旨在提高ABC算法局部寻优能力并加快收敛速度。具体而言，式(3)表示将上次迭代中的最优个体作为被继承的父代，在最优解邻域进行探索。

两个种群每隔一定的代数就分别随机选取若干个个体并交换其在种群中的位置，进而促进优良蜜源信息在种群间的传播。该操作可在平衡全局寻优和局部寻优能力的同时维护种群的多样性，进一步提高解的搜索效率。

1.1.3 基本流程

基于双种群架构的ABC算法流程见图 1，具体如下。

(1) 生成初始蜜源：设置最大蜜源数、局部限制次数limit、最大迭代次数T。将蜜源平均分为两部分，即种群1和种群2。

(2) 种群1中的引领蜂搜寻蜜源：在种群1中，引领峰按照标准ABC算法搜索策略(式(1))搜索新蜜源，对比新蜜源适应度f(U^1a)与上一次迭代发现蜜源的适应度f(X_i)并择优选择，若发现新蜜源的适应度优于上次迭代发现蜜源，则保留新蜜源，直接进入步骤(3)；反之，按照式(2)借鉴交叉和变异的思想搜索新蜜源，再对比新蜜源的适应度f(U^1b)与上一次迭代发现蜜源的适应度，择优选择新蜜源。

(3) 种群2中的引领峰搜索蜜源：在种群2中，按照式(3)在最优解附近搜索蜜源，同样对比新蜜源适应度f(U²)与上一次迭代发现蜜源的适应度，择优选择新蜜源。每迭代一定次数，两个种群随机交换种群中若干个蜜源的位置, 生成新蜜源。

(4) 跟随蜂的局部寻优：完成上述过程后，引领蜂飞回，与跟随蜂交流，与标准ABC算法一样，利用轮盘赌的方式计算跟随概率，然后选择是否更新。

(5) 判断蜜源X_i经历的局部搜索次数k是否大于限制次数limit：若是，则该引领蜂转变为侦察蜂，同时按照初始蜜源的生成方式产生新蜜源; 若否，保留该蜜源。

(6) 转到步骤(2)进行下一次循环, 直到满足最大迭代次数。

1.2 算法测试

本文选择三种典型的测试函数^[25]检验改进ABC算法的性能。

Griewank函数是典型的非线性多模态函数，其全局最小值0在(x₁, x₂, …, x_n)=(0, 0, …, 0)处取得。此函数存在许多局部极小值点，并随着问题维数的增加而增加。此函数在优化过程中，具有广阔的搜索空间，通常被认为是优化算法很难处理的复杂多模态问题。其函数形式为

$ \begin{aligned} f_{1}\left(x_{n}\right) &=\sum\limits_{n=1}^{N} \frac{x_{n}^{2}}{4000}-\prod\limits_{n=1}^{N} \cos \frac{x_{n}}{\sqrt{n}}+1 \\ &-600 \leqslant x_{n} \leqslant 600 \end{aligned} $

(4)

式中N表示函数的维度。

Resenbrock函数是一个测试最优化算法性能的非凸函数，全局最小值0在(x₁, x₂, …, x_n)=(1, 1, …, 1)处取得，可用来检测算法的精度。其函数形式为

$ \begin{gathered} f_{2}\left(x_{n}\right)=\sum\limits_{n=1}^{N-1}\left[100\left(x_{n+1}-x_{n}^{2}\right)^{2}+\left(x_{n}-1\right)^{2}\right] \\ -15 \leqslant x_{n} \leqslant 15 \end{gathered} $

(5)

Rastrigin函数存在多个波峰与波谷，易陷入局部极小，全局最小值0在(x₁, x₂, …x_n)=(0, 0, …, 0)处取得，可用来检测算法跳出局部极小的能力。其函数形式为

$ \begin{array}{c} {f_3}\left( {{x_n}} \right) = \sum\limits_{n = 1}^N {\left[ {x_n^2 - 10\cos \left( {2{\rm{ \mathsf{ π} }}{x_n}} \right) + 10} \right]} \\ - 15 \le {x_n} \le 15 \end{array} $

(6)

分别利用标准和改进的ABC算法对上述三种函数进行寻优测试，迭代次数设为500，每种算法独立重复试验30次，试验结果见表 1和图 2。可看出，与标准ABC算法相比，改进后的ABC算法不仅有较快的收敛速度，收敛精度也有较大提高，具有较强的全局寻优能力。

2 MT和重力数据联合反演目标函数

根据正则化反演理论^[26]，MT数据和重力数据正则化联合反演的目标函数为

$ \begin{gathered} P^{\beta}(\boldsymbol{m})=w^{\mathrm{MT}}\left(\left\|\boldsymbol{d}_{\mathrm{obs}}^{\mathrm{MT}}-\boldsymbol{d}_{\mathrm{cal}}^{\mathrm{MT}}\right\|^{2}+\beta^{\mathrm{MT}}\left\|\nabla \boldsymbol{m}^{\mathrm{MT}}\right\|^{2}\right)+ \\ w^{\mathrm{G}}\left(\left\|\boldsymbol{d}_{\mathrm{obs}}^{\mathrm{G}}-\boldsymbol{d}_{\mathrm{cal}}^{\mathrm{G}}\right\|^{2}+\beta^{\mathrm{G}}\left\|\nabla \boldsymbol{m}^{\mathrm{G}}\right\|^{2}\right) \end{gathered} $

(7)

式中：P^β(m)表示联合反演目标函数，m表示模型向量，β是正则化因子；$ \boldsymbol{d}_{\mathrm{obs}}^{\mathrm{MT}}、\boldsymbol{d}_{\mathrm{obs}}^{\mathrm{G}} $和$ \boldsymbol{d}_{\mathrm{cal}}^{\mathrm{MT}} 、\boldsymbol{d}_{\mathrm{ca}}^{\mathrm{G}} $中，d代表地球物理数据，上标“G”和“MT”分别代表重力和MT地球物理方法，下标“obs”和“cal”分别代表实际观测值与理论计算值；w为加权因子。

本文采取双正则化因子自适应调整方案^[27-28]。MT数据正演采用有限差分法，重力数据正演采用网格累加法。联合反演建模采用于鹏等^[29]提出的物性参数随机分布的共网格模型建模技术，该技术已被应用于多种地球物理方法的二维联合反演^{[5-6, 30-33]}。

3 模型测试 3.1 模型试验一

为了验证ABC算法在二维MT和重力数据联合反演中的应用效果，设计了图 3所示的电阻率和密度模型。模型的背景电阻率为1000Ω·m，密度为2.68g/cm³。模型中包含四个异常体。MT数据的模拟频率范围为0.001~320Hz，按对数等间距取38个频点，测点间距为1km。

设计了三种反演方案：方案一为基于标准ABC算法的联合反演；方案二为基于双种群架构的ABC算法单独反演；方案三为基于双种群架构的ABC算法联合反演。这三种方案的蜜源数量最大值均设为300，最大迭代次数为50。种群中初始蜜源位置按照真实模型上下扰动一定范围随机生成：异常体的上下界面深度扰动范围为真实值±30%，电阻率扰动范围为真实值±30%，密度扰动范围为真实值±0.02g/cm³。

初始模型及反演结果见图 4。图 4a、图 4b是随机生成的电阻率和密度的初始扰动模型。首先，对比基于改进ABC算法的联合反演(方案三)和单独反演(方案二)结果，可以看到，虽然重力异常曲线(图 4j)的拟合效果都较好，但是单独反演结果(图 4f)的物性及界面还原效果较差。不同方案的反演结果与真实模型的残差(图 5)也表明，联合反演优于单独反演。

对比基于标准ABC算法的联合反演(方案一)与基于改进ABC算法的联合反演(方案三)结果可见，在标准ABC算法的反演结果(图 4c、图 4d)中，界面起伏严重，物性差异大；而采用基于双种群架构的联合反演(图 4g、图 4h)，物性界面起伏小且物性分布连续。由图 5的残差图也可以看出，标准ABC算法的反演结果偏离真实模型较大。此外，从迭代误差曲线(图 4i)也可以看出，改进ABC算法收敛快、寻优精度高。

3.2 模型试验二

为了检验本文反演算法对复杂模型的适应性，本文设计了图 6所示的电阻率界面与密度界面不完全一致的地球物理模型。该模型的背景电阻率为500Ω·m，背景密度为2.65g/cm³。反演过程中，假设电阻率和密度都是完全共界面进行扰动，其余参数设置与前述模型试验相同。初始扰动模型见图 6c和图 6d，基于改进ABC算法的联合反演结果见图 6e和图 6f。可以看出，反演的电阻率和密度异常分布与理论模型基本相符。从图 6i所示的迭代误差曲线及图 6j所示的重力异常拟合曲线也可以看出，反演得到的异常体电阻率和密度分布基本可靠。

4 实测数据应用

为进一步验证改进ABC算法的实用性，选择中国下扬子地区某测线MT和重磁实测数据进行联合反演。

依据测线实测数据的OCCAM反演结果、物性统计及文献[33]所揭示的构造格架建立初始模型，并对待解参数进行随机扰动(图 7)。反演过程对Tg印支面以上的界面不做反演，对其以下界面进行开放反演。物性约束范围根据该测区的岩石物性资料统计给定，见表 2。

对比电阻率和密度的初始扰动结果，可见电阻率和密度的联合反演剖面(图 8)具有明显的成层性。此外，在反演剖面横向185~200km、纵向-5~-8km范围内，存在物性分布变化剧烈的地层，结合两侧的物性分布特征，推测此处可能存在断裂。反演的密度和电阻率模型的正演响应与实测数据的对比分别见图 9和图 10。上述认识与文献[33]基于模拟退火算法的联合反演结果也基本相符，验证了本文提出的双种群ABC算法的实用性。

5 结束语

(1) 本文通过双种群架构将蜂群一分为二。种群1结合标准ABC算法和变异交叉操作策略，保证了算法的全局寻优能力，又提高了解的搜索效率；种群2则利用最优解邻域的搜索，增强了算法的局部寻优能力。算法测试和模型试验结果证明双种群架构可以提高ABC算法的寻优能力。

(2) MT和重力数据联合反演模型试验及实际数据处理结果表明，MT与重力数据的联合反演优于单独MT或重力数据反演，改进ABC算法可以应用于实际数据处理，具有一定的实用性。

需要指出的是，ABC算法具有很强的开发潜力，多种群的设置方案也是不唯一的，如何进一步提升该算法的寻优能力值得进一步研究。