基于混沌天牛群算法的大地电磁反演

谢卓良, 王绪本, 李德伟, 陈先洁, 乃国茹

(成都理工大学 地球物理学院，成都 610059)

0 引言

大地电磁测深法(MT)作为地球物理勘探最主要的手段之一，其反演问题一直是研究热点[1-2]。目前MT一维线性反演方法很多，如梯度法、差商牛顿法、马奎特法和广义逆矩阵法等，此类方法依赖于初始模型的选择，全局寻优能力比较差，容易陷入局部最优解，鲁棒性较差，有时甚至无法收敛[3-4]。鉴于此原因，国内、外学者对非线性反演展开了一系列研究，如模拟退火算法、重力搜索算法、遗传算法、人工神经网络及各种启发式智能仿生算法等。柳建新等[5]将实数编码遗传算法应用于大地电磁二维理论模型的反演；师学明等[6]将阻尼粒子群算法应用于大地电磁的优化反演；胡祖志等[7]将人工鱼群算法应用于大地电磁的最优化约束反演；王鹤等[8-10]基于神经网络对大地电磁的非线性反演做了大量研究；Godio A等[11]再次将粒子群算法引入大地电磁反演，并与Occam反演方法进行了对比；王鹏飞等[12]在大地电磁的理论模型反演中引入改进布谷鸟算法。大量的研究表明，非线性反演方法在地球物理反演中应用前景十分广阔。天牛须搜索算法(Beetle Antennae Search Algorithm,BAS)是理想化甲虫类生物觅食和求偶的活动，并将其建立起数学模型的一种智能算法[13]。算法主要以天牛为例，天牛的头部左右两侧具有两根触须，当天牛附近区域存在目标体(食物或者异性)时，其头部的触须检测目标体所释放的气味，若右触须感知到的气味更强，天牛接下来往右侧运动，反之，就往左侧运动。根据天牛找寻目标体的这一生物活动，建立理想化数学模型，目标体或者说目标体的气味相当于目标函数，然后利用此数学模型进行搜索寻优。

针对标准的BAS算法，一些学者提出天牛须搜索算法虽然能有目标地进行搜索，但是当搜索目标增多，尤其大于4个目标体时，收敛结果容易陷入局部极值，而且在多峰值函数优化时，陷入局部更容易陷入局部最优。赵玉强等[14]提出一类带学习与竞争策略的混沌天牛群搜索算法(Learning and Competing Chaos Beetle Swarm Algorithm,LCCBSA)，这种改进的算法将一个天牛的搜索变成群体天牛的搜索，之后引入参考粒子群算法的指导性学习策略及竞争场策略。引入的策略使得天牛有目的的搜索且避免了陷入局部搜索的问题。这里利用LCCBSA算法对四个测试函数试算并与粒子群算法和遗传算法的结果进行对比，该算法搜索效率和收敛速度相比其他算法有一定的优势，于是进一步应用于大地电磁反演问题中。

1 大地电磁测深反演理论

大地电磁一维层状模型视电阻率可用以下函数表示，即

ρa(f)=F(ρ1,ρ2…,ρn,h1,h2…,hn-1,f)

(1)

其中：ρi(ρ1,ρ2,…,ρn)表示模型电阻率参数;hi(h1,h2,…,hn-1)表示厚度参数;f表示频率。大地电磁的反演的本质是根据实际观测的ρa(f)求解模型参数ρ1、ρ2、…、ρn、h1、h2、…、hn-1。

将混沌天牛群算法应用于一维大地电磁的n层模型反演中，假设有m个天牛，相当于从m个初始模型出发，利用Tent映射将随机向量(ρi,1,ρi,2,…,ρi,n,hi,1,hi,2,hi,n-1)T(i=1,2,…,m)对天牛群进行初始赋值，参数个数D=2n-1,并给出解空间每一个参数的数值大小范围，然后利用混沌天牛群算法在解空间内搜索最优解。选取胡祖志等[15]所提出的大地电磁反演目标函数，如式(2)所示。

(2)

其中：N表示观测频点个数；前一项表示观测视电阻率和计算视电阻率的误差；后一项表示观测相位和计算相位的误差。

2 混沌天牛群搜索算法(LCCBSA)

2.1 标准BAS算法介绍

单个天牛个体朝着触须感知到的适应度优势的方向运动，最终搜索到全局的最优解。假设在解空间为D维，天牛的当前位置为X=(x1,x2,…,xD),则天牛左侧和右侧两只触须的搜索行为可用式(3)定义。

(3)

(4)

2.2 初始种群混沌化

在对群智能算法的研究过程中，很多文献表明算法的收敛性与群体初始值的分布相关，初始群体分布如果具有更好的均匀性和遍历性，算法的收敛性越好。笔者采用Tent映射得到天牛初始种群，Tent映射公式如式(5)。

(5)

式中：i=1、2、…、D表示混沌变量；n=1、2、…、N表示种群个体；μ表示混沌形态控制参数，一般情况下，其值在接近2时使得混乱状态较好。为了得到均匀性和遍历性较好的混沌初始种群，首先，种群中第一个天牛个体随机生成，取值范围为[0，1]，然后从第二个天牛个体开始初始值依次迭代产生，直到第N个天牛，最后再映射到解空间中。

2.3 指导性学习策略

在标准BAS算法中，天牛个体的搜索能力取决于触须的长度，从当前位置向相邻位置的移动并不容易，个体之间信息的交流限制较大。为了提升寻优能力，笔者参考粒子群算法个体位置更新的机制，引入指导性学习策略更新天牛个体位置。思路如下：群体中的天牛个体在当前位置会对目标体有一个感知方向和历史群体最优个体的感知方向，当前个体感知目标的同时学习历史群体最优值，向群体最优方向靠拢，到达图中右侧天牛位置(图1)。策略将感知目标这一活动定义为认知部分，学习历史群体最优这一活动定义为社会经验部分。这种兼顾学习和自身判断的策略，由于受到了认知部分和社会经验部分的双重影响，使得算法在多目标寻优中不会陷入局部最优。

图1 天牛寻优示意图

(6)

个体位置更新公式为式(7)。

(7)

(8)

2.4 竞技场策略

在实际实验中，群智能算法寻优在后期搜索时往往会陷入局部最优，难以收敛到全局最优解。Cheng R[16]提出优化大规模问题时可以采用竞技的方法，蒋莹莹等[17]提高算法精度和收敛速度时采取了个体之间相互竞争的方式。笔者提出竞技场策略，思路为将群体中个体的聚集度进行比较，进而竞争产生适应度高的个体，并对适应度小的个体进行位置更新。某代群体中个体聚集度的定义如式(9)所示。

(9)

其中，S表示两只天牛个体的相似度，定义如下：

(10)

从式(9)中发现，天牛群体中任一个体都需要计算与其他个体的欧几里得范数。阈值δ的值本文中取所得欧式范数的中位数，若其大于阈值δ，则两个天牛个体之间无相似性；否则就表示两者相似，相似度增加1，相似度累加值表示个体聚集度。显然，聚集度高的个体越多，个体之间越相似，群体的多样性越少，对算法的全局搜索不利，算法最终容易陷入局部极值。

为了解决群体中上述问题，将某代群体中聚集度的个体降序排序，选取前m个聚集度高个体参加竞争，通过计算它们的适应度，对其中n个适应度较小的个体位置更新。这一部分个体的位置更新采用偏好随机的方式如式(11)所示。

(11)

2.5 算法流程

LCCBSA算法流程如图2所示。

图2 LCCBSA算法简要流程图

3 测试函数及大地电磁应用

3.1 测试函数试算

为了验证本文提出的LCCBSA算法的性能，利用测试函数进行试算是最常用的手段。笔者选取4个多维基准测试函数，函数的基本信息如表1所示，函数表达式中D表示变量的维数。

表1 测试函数

在数值模拟实验中，选择已经较为成熟的粒子群算法(PSO)和遗传算法(GA)与本文的LCCBSA算法进行对比，三种算法种群规模设置N=100，最大迭代次数200次，LCCBSA算法参数设置为：Cmax=2.5，Cmin=1.5，s=1.5，l=2.0，m=N/3，n=m/2；PSO算法参数设置为：惯性权重w=0.9，学习因子C1=C2=1.5；实数GA算法参数设置为：交叉概率Pc=0.8，变异概率Pm=0.1。表2为利用这三种算法对4个基准函数的分别重复进行30次实验最优解的平均值对比。

表2 平均最优解对比

在相同测试函数和测试维度条件下，LCCBSA算法除F2测试函数外，其他三个测试函数较PSO算法和GA算法都有更好的寻优能力，且具有更快的收敛速度。

图3 测试函数目标函数值对比

3.2 MT理论模型反演

3.2.1 三层(H型)理论模型反演

构建H型反演模型，反演过程中，迭代次数设置为2 000次，混沌天牛群其他参数设置与前面的测试函数一致，最小误差ε=1×10-3，表3给出三层H型地电模型参数及LCCBSA算法反演结果，未添加噪声正演响应曲线与反演计算曲线见图4，两者具有很好的拟合效果。

表3 三层H型地电模型LCCBSA算法反演结果

图4 三层地电模型观测数据与反演结果计算值对比

图5 三层地电模型含噪10%观测数据与反演结果计算值对比

3.2.2 六层理论模型反演

设计一个六层地电模型，模型参数、模型值和反演结果见表4。该模型为参考文献[18]中计算大地电磁多尺度反演(MI)时所设计的模型，其主要是为了检验反演方法恢复地下深部低阻层的能力(如表4中9 450 m处的50 Ω·m的低阻地层)。LCCBSA算法反演参数设置同三层理论模型，未添加噪声正演响应曲线与反演计算视电阻率曲线见图6，由表4可以看出，MI反演结果的最大相对误差达到44.6%，远大于LCCBSA的最大相对误差，说明LCCBSA反演结果较MI反演结果更好，而且对深部的低阻层电阻率恢复也更加令人满意。

图6 六层地电模型观测数据与反演结果计算值对比

3.2.3 含噪理论数据分析

由于MT野外采集数据受到各种噪声的影响，造成反演算法的不稳定性，为了检验本文算法的抗噪能力，对上述三层和六层理论正演计算结果数据添加10%的高斯噪声，之后使用LCCBSA算法对含噪数据进行反演计算。三层模型和六层模型含噪数据反演结果分别见表3和表4，正演响应曲线与反演计算曲线分别见图5和图7。

从表3和表4中可以看出，理论模型含噪数据相对于无噪数据的反演结果，与真实模型值的差异增大，可见噪声会对反演结果的准确性产生一定的影响，三层模型和六层模型参数的最大相对误差分别为4.72%和14.85%，说明随着模型的复杂度增加，添加噪声后的数据反演结果差异也会变大。虽然噪声一定程度上会影响反演结果的准确性，但是反演结果的相对误差较小，与真实值相差不大，且图5和图7可以看出反演计算视电阻率曲线和相位曲线能与正演理论曲线拟合效果较好，说明含噪反演结果准确可靠，LCCBSA算法的抗噪性较好。

表4 六层地电模型LCCBSA算法反演结果对比表

图7 六层地电模型含噪10%观测数据与反演结果计算值对比

3.3 MT实测数据反演

应用新疆塔里木盆地布设的某条MT实测数据，验证LCCBSA的反演效果。工区内第四系和新近系沉积大套砾石层，但其横纵向分布范围变化剧烈，造成地震资料反射凌乱，影响圈闭的准确落实。在实际生产工作中发现这套砾石层在电测井资料显示为高阻特征，利用大地电磁刻画可以取得较好的效果，且Bostick反演结果与地质规律吻合率较高。本文中选取的测线长为24.075 km，点距为500 m，频点38个。对该测线进行分别进行LCCBSA反演和Bostick反演。

LCCBSA算法反演实测数据时，建立了首层厚度为50 m，递增因子为1.1，总计37层的初始层厚模型，反演的初始电阻率模型为经过预处理后的电阻率，反演时保持层厚度不变，只对38个电阻率参数进行反演，电阻率搜索空间为经预处理后电阻率的50%～150%。LCCBSA算法计算所需要的天牛群数100个，迭代次数2 000次，认知因子Cmax=2.5，Cmin=1.5，学习因子s=1.5，天牛须长度l=2.0，迭代终止条件ε=1×10-3。反演时对测线上的每个测点独立反演20次，求取20次结果的平均值作为每个测点反演结果。

图8为bz1井旁测点的反演结果，图8(a)显示实际观测数据与LCCBSA反演结果计算的视电阻率值变化趋势一致，且拟合效果较好;图8(b)为LCCBSA反演结果与Bostick反演结果和测井曲线的对比，可以发现LCCBSA和Bostick反演的电阻率与测井曲线具有很好地吻合，但在深部6 000 m处LCCBSA相对于Bostick反演电阻率偏大。

图8 bz1井旁测点实测数据反演结果

图9为实测电阻率和LCCBSA反演计算视电阻率拟断面图的对比，可以看出反演计算的视电阻率和实测电阻率的一致性较好，说明反演的电阻率较为可靠。图10是LCCBSA和Botick拟二维反演解释剖面图的对比，两种方法反演剖面具有一定的相似性，推断解释存在三套地层，第一层为高电阻率的第四系地层(Q)，电阻率一般在1 000 Ω·m左右，主要分布未成岩砾石；第二层为次高阻的新近系地层(N)，一般电阻率为几十至几百欧姆米，主要分布准成岩和成岩砾石；第三层低电阻率的古近系地层(E)和白垩系地层(K)，主要岩性为砂岩、泥岩。同时，LCCBSA算法反演剖面深部低阻范围较Bostick反演剖面低阻范围大，是由于初始层厚模型递增变化导致深部层厚大，可以利用测井和地震解释层位约束，反演出更为准确的低阻层。

图9 实测视电阻率和反演计算视电阻率拟断面图对比

图10 Bostick和LCCBSA算法反演剖面对比

4 结论

1)测试函数的试算表明，本文介绍的混沌天牛群搜索算法(LCCBSA)，在大部分情况下相对于粒子群算法和遗传算法收敛速度较快，具有一定的优势性。

2)MT理论模型和实测数据的试算结果表明，LCCBSA算法能够较好地恢复地下介质的电阻率，是一种比较可行的非线性反演方法，避免了传统线性反演方法初始模型的选取、陷入局部最优和计算雅可比偏导数矩阵的问题，并且可以用于航空电磁、瞬变电磁、地震等地球物理方法的反演。

3)LCCBSA算法跟其他非线性反演方法一样，模型空间的选取、参数的改变都对反演效率有较大的影响，当前LCCBSA算法尚不成熟，应用于大地电磁的实验是一种尝试，并且其计算时间长的缺陷，极大地影响其向二、三维反演方向发展，但是将其和神经网络结合是一个具有研究意义的方向。