基于遗传退火进化算法的瞬变电磁非线性反演

司元雷 郑龙金

1(中国矿业大学资源与地球科学学院 江苏 徐州 221116) 2(江苏建筑职业技术学院信电工程学院 江苏 徐州 221116) 3(江西有色地质勘查五队 江西 九江 332000)

0 引 言

瞬变电磁法(简称TEM)属于时间域电磁法，它采用不接地的回线或接地线源向地层发射一次脉冲磁场[1](也称一次场)，在发射波形关断的期间，通过线圈或者感应探头来接收地下目标体被激发而产生的二次磁场，以此来解决与之相关的地球物理探测问题[2]。TEM自20世纪被提出后，经过几十年的发展，目前已在多个勘测领域得到了广泛的应用，例如能源开发、矿产勘探、水文地质、工程物探、考古调查、军事和环境保护等领域均有所涉及[3]。

目前，瞬变电磁数据的高维精准反演是摆在广大研究人员面前的一道大难题。尽管TEM的高维反演理论研究工作早已展开，但由于高维正反演的数值计算量非常之巨大，因此在这一方面一直未曾有较大的突破。正演是反演的基础，瞬变电磁高维反演还是更多地停留在理论研究阶段，大规模应用于工业生产还有一段很长的路要走，瞬变电磁法在目前阶段的反演解释主要还是使用较为实用的一维反演方法[4]。实测资料的反演一直是TEM的核心问题之一。为此，许多地球物理学家也提出了非常多的反演方法。例如烟圈反演法、浮动薄板解释法、人机对话自动反演法、成像类反演算法、马奎特法和广义逆反演法等，而且这些方法在国内外也都得到了不同程度的应用。然而，传统的反演方法大多基于最小方差原理[5]，其反演结果依赖于初始模型的选择，容易陷于局部极值，并且巨大的矩阵运算量在反演过程中需要耗费很多时间，但是反演结果却仍然不稳定。从这些反演方法的实质上来讲，它们都是将非线性的问题进行线性化来处理[5]，并非真正意义上的非线性方法。因此，采用完全非线性方法反演瞬变电磁数据势在必行，对于提升反演结果的稳定性和精度也同样有着非常重要的实际意义[6]。

鉴于大多数地球物理反演问题具有高度的非线性，因此，自非线性反演提出以来，便引起了相关研究学者的高度重视。它们通过各种方式非间接地求解非线性问题，将数据空间映射到模型空间，进而反演求解问题被归化为求一个泛函极值问题的解，并采用各种优化方法来解决对应的问题，从而在求解非线性地球物理反演问题的过程中能够快速地得到与之对应的模型参数[1]，最终更加精准快速地解译野外实测物探资料。凡事都具有两面性，非线性反演方法也有不足之处。它虽然可以在全局中自动寻找最优化的解，但是对算法参数的要求却较为严格，一旦相关参数选取不合适或将会促使反演出现早熟，或出现无法收敛的情况[1]。目前，现有的非线性反演算法和线性反演算法都依赖于初始模型的设置，若设置不得当，很有可能得到的并不是最优解。除此之外，也存在着容易陷入局部极值、收敛速度缓慢等一系列问题[7]。例如：遗传算法(Genetic Algorithms，GA)可以从概率的角度随机找到最优解，且把握全局搜索过程的能力很强[8]，但在实际运用中也存在一些问题，例如早熟现象和较差的局部优化性能；模拟退火(Simulated Annealing，SA)算法在局部搜索方面具有优势，可以避免在搜索过程当中落入局部最优解，但它在把握全局搜索空间方面力不从心，搜索过程不能进入最有希望的寻优区域，运行效率低下[6]。本文首先通过测试函数对遗传退火进化算法(Genetic Annealing Evolutionary Algorithm，GAEA)进行数值实验，然后将其应用于瞬变电磁地电模型理论数据的反演，从而验证反演算法的有效性，以期提升瞬变电磁反演解释工作的准确性。

1 遗传退火进化算法

GA是基于自然遗传学进化思想的启发式搜索算法，它基于对随机搜索的智能化开发，将其转移到获得问题最优解的方向。遗传算法的概念是约翰·霍兰德(John Holland)提出的，它是一种解决优化问题并克服传统优化方法所带来的局限性的方法。搜索过程从一组称为种群的染色体开始，每个染色体表示为问题解的编码。每一代种群通过使用遗传算子的重组过程随着时间的推移而进化，并朝着更好的个体进化以产生新种群。为了量化个体的适应度，使用了适应度函数，基于适者生存原理的自然选择和遗传操作生成新一代个体。GA是一个迭代过程，需要按周期进行操作，直到满足收敛或停止标准为止[7]。SA算法被认为是一种模拟金属退火过程的迭代随机搜索算法。SA算法的基本思想是搜索当前个体的邻居以找到新的个体，它具有避免陷入局部最优状态的能力，这是因为SA算法在每次迭代中都会将新个体与当前个体进行比较。如果新的个体更好(基于其适应性值判断)，那么它将被选作下一次迭代的基础；否则，根据收受概率，SA算法可能会采用此个体[9-10]。

GAEA是一种混合遗传模拟退火算法[11]，它通过组合GA和SA算法的优势来求解最优化问题[12-22]。GA是此混合算法的主要框架，而SA算法被用作局部搜索策略，以帮助GA跳出局部最优值。本文中的遗传退火进化算法主要参考邢文训等[10]及何则干等[16]给出的方法，这里仅给出该算法的流程图，如图1所示，其余部分将不再赘述。对于瞬变电磁测深的反演建模，首先选择一组(或总体)参数模型，并定义了模型的电阻率和层厚度等参数，然后设置每个模型参数的搜索范围[23]：

(1)

式中：mk,j是参数向量；n是参数个数。在本文使用的遗传算法中，模型的所有参数都以二进制形式进行编码。

2 具体实现

2.1 使用经典测试函数进行对比

在进行瞬变电磁非线性反演之前，为了检验这种GAEA的有效性，这里首先选用两个经典的测试函数：Schaffer函数及Alpine函数。其中Schaffer函数是一种含有无数个极小值点的多维函数，由于其分布呈现强烈震荡形态，因而寻找全局最优值的难度非常大；而Alpine函数是一种多峰最小化测试函数。当函数趋于无穷大时，该函数将沿着自变量的方向生成大量不同的局部极值，这给优化带来了很大的困难。因此，可将GA、SA和GAEA三种算法进行比较。两个测试函数分别为：

-100≤xi≤100i=1,2

(2)

-10≤xi≤10i=1,2,…,n

(3)

式中：Schaffer函数在其定义域内只有一个全局最小点f(0,0)=0；x*为当Alpine函数取得全局最优值时xi的值，当x*=0时，Alpine函数取得全局最优值f(x*)=0。函数图形如图2所示。

(a) Schaffer函数

(b) Alpine函数图2 两个测试函数的立体图

2.1.1测试函数一

将种群设置为60，精度设置为1e-20，使用GA、SA和GAEA三种算法分别求解Schaffer函数的最小值点f(0,0)=0，图3是通过三种算法计算所得到的每一代种群中最小目标值的寻优过程。可以看出GAEA精度最高，SA算法精度最低。另外，相对于GA，GAEA不仅寻优能力更强且收敛速度更快。

图3 三种算法在Schaffer函数寻优过程的目标值图

鉴于上述三种算法的搜索方式具有随机性，为了更客观地评价算法的性能，这里将算法的容许误差设置为1e-20，三种算法独立地运行60次，得到的寻优结果如表1所示。

表1 三种算法的计算结果对比

通过对比表1中三种算法的计算结果，GAEA相对GA而言精度更高，GAEA得到的最优解与Schaffer函数的全局最小值点最相近，远比SA算法优越。

2.1.2测试函数二

将种群设置为60，精度设置为1e-20，使用GA、SA和GAEA三种算法分别求解Alpine函数的最小值点f(x*)=0，图4是通过三种算法计算所得到的每一代种群中最小目标值的寻优过程。可以看出GAEA精度最高，SA算法精度最低。另外，相对于GA，GAEA不仅寻优能力更强且收敛速度更快。

图4 三种算法在Alpine函数寻优过程的目标值图

同样地，这里将算法的容许误差设置为1e-20，三种算法独立地运行60次，得到的寻优结果如表2所示。

表2 三种算法的计算结果对比

从运算结果中可知，GAEA有12次得到全局最优解0；远远好于SA算法和GA。综合两个函数的测试结果可知，GA和SA算法均以一定的概率逼近全局最优解。从比较结果不难发现，GAEA逼近全局最优解的概率较大，落入局部极值的概率较小。

2.2 层状地电模型反演

2.2.1反演目标函数的确定

通过反演算法对瞬变电磁数据进行的分析包括使用感应电动势或视电阻率的值来估计被分析介质的真实电阻率、厚度或深度值。对于一维情况，水平分层介质被广泛使用。该模型的自由参数是每层的电阻率ρn和厚度hn，由向量m表示，同时地电模型参数可被记为x。在这种情况下，反演问题即可被转化为一个求解多目标非线性优化问题。因此，构建目标函数S(x)，层状地电模型的最优解xopt需满足S(xopt)=min(S(x))，即：

(4)

2.2.2遗传退火进化算法反演参数的选择

基于上述的目标函数，针对算法参数进行了详细的试算与筛选，初始群体由随机概率产生，种群的大小设置为60，最大迭代次数设置为30 000，种群中个体的交叉概率设置为0.7，每一位个体的变异概率设置为0.1，个体的选择方法使用轮盘赌法，采用浮点法进行编码，允许误差设置为1e-04。模拟退火初始温度T0为1 500℃，降温方式T=T0×0.98k-1，模拟退火的终止步数k设置为20。另外，在模拟退火算法部分，随机生成小的扰动。

2.2.3不同地电理论模型试算

为验证上述反演算法在瞬变电磁资料解释中的可行性，使用不同瞬变电磁理论模型进行了反演试算。其中瞬变电磁正演计算部分采用文献[24-27]中的方法。反演过程中先建立n层介质模型，其参数为λ=(ρ1,ρ2,…,ρn,h1,h2,…,hn-1)，ρ和h分别为各层介质的电阻率和厚度；使用式(4)构建反演算法的目标函数，设定初始模型后进行迭代计算直至达到拟合终止条件。计算中，发射线圈和接收线圈的边长分别设置为100 m和1 m，匝数均为1，供电电流为10 A。为了与传统线性方法进行比较，除采用GA、SA算法、GAEA方法外，还采用传统线性化反演方法的代表马奎特方法(Marquardt)进行了反演计算[28]，每种算法独立运行20次，然后取20次反演结果的均值作为模型参数的估计值。其中：GA的种群的大小设置为60，最大迭代次数为30 000，个体的交叉概率为0.7，变异概率为0.1，选择方法使用轮盘赌法，采用浮点法编码，允许误差设置为1e-04。在交叉与变异运算中，进行交叉运算的个体和交叉的位置及变异位，均通过产生的随机数与交叉概率与变异概率的大小比较来判定。SA算法的模型群体数量为60，最大迭代次数为30 000，初始温度T0为1 500 ℃，降温方式T=T0×0.98k-1，模拟退火的终止步数k设置为20，随机产生小的扰动。下面将分别针对选取的三层低阻模型及五层地电模型进行反演试算，并分析相应的反演结果。

(1) 三层H型地电模型。表3为三层H型地电模型的参数及多种算法的反演结果，同时，GAEA的反演结果经过正演计算得到感应电位之后，与真实模型的正演响应进行了对比(如图5所示)，二者的感应电动势曲线基本重合，证明反演结果是有效的。另外，将表3中的反演结果绘制成柱状图可以更加直观地看出每种算法反演的效果，如图6所示，可以看出GAEA反演效果最好，其各参数反演结果的误差均控制在2%以内，可以满足实际需要。

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

图5 原始模型响应和GAEA反演结果响应的对比

图6 反演结果柱状图对比

(2) 五层HKH型地电模型。表4为五层HKH型地电模型的参数及多种算法的反演结果，GAEA的反演结果经过正演计算得到感应电位之后，与真实模型的正演响应进行了对比，如图7所示，二者的感应电动势曲线基本重合，证明反演结果是有效的。从图8中的直观结果易知，与Marquardt法、GA、SA算法相对比，GAEA在反演TEM数据方面更为准确和有效。因为五层模型参数相较于三层更多，所以对应的反演难度加大，从而导致了结果的误差也相应增大，但是该模型的关键信息仍可以由GAEA获取，其大部分参数的反演结果的误差控制在6%以内，仅个别参数误差超过10%(第三层电阻率参数ρ3的误差为12.9%)，这在瞬变电磁法勘探的实际生产中，仍能满足实际需要。

表4 五层HKH型地电模型反演结果

图7 原始模型响应和GAEA反演结果响应的对比

图8 反演结果柱状图对比

3 结 语

本文将一种新的GAEA应用于瞬变电磁法的反演中，与传统线性反演方法的对比实验验证了GAEA的有效性和可行性。主要结论如下：

(1) 经融合后得到的遗传退火进化算法，结合了GA的优秀全局搜索性能和SA算法的局部搜索能力强的特点。与此同时，种群在模拟退火操作中还能充分利用融合算法所带来的全局信息。

(2) GAEA方法在1D合成数据中的实现和应用对于更好地理解这些方法的优缺点具有重要意义。合成数据的反演结果表明，基于遗传退火进化算法的瞬变电磁非线性反演比线性反演更有效，为2D反演算法的发展奠定了基础。

(3) 将GAEA应用于两个测试函数，即Schaffer函数及Alpine函数，该技术在寻找测试函数的最优极值方面具有非常好的效果；将GAEA应用于瞬变电磁地电模型的合成数据，模型的关键信息仍可以由GAEA获取，但是需要大量的处理时间。这验证了该算法的适用性，其缺点可以在将来的实施中加以改进。