基于改进粒子群算法的大地电磁阻抗张量分解方法

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

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

0 引言

在大地电磁测深法中往往存在由于浅地表电性结构不均匀，使得在不均匀体的周围电流密度发生畸变，从而使电场分量发生变化，引起局部的畸变异常[1-4]。为消除此种局部畸变异常带来的影响，Groom和Bailey假设局部三维异常体覆盖在二维区域异常之上，将实际阻抗张量分解为畸变阻抗张量和区域阻抗张量，以此来恢复未受畸变影响的区域阻抗张量(GB分解法)。它的本质是求解由区域阻抗张量和畸变方程组成的非线性超定方程组[6-8]。传统的GB分解法将非线性问题线性化求解，导致求解结果不稳定，易陷入局部极值中，很难得到较为准确的分解结果。李洋等[5]提出来一种混合优化算法，利用模拟退火法得到的全局最优解作为局部非线性最小二乘法的初始值，提高了求解稳定性；尹曜田等[6]研究了基于遗传算法的阻抗张量分解方法，采用遗传算法求解非线性方程组；蔡军涛等[7]提出共轭阻抗变换的概念，不受电场畸变影响，无需地下电性维性信息的假设；谢成良等[8]研究了基于相位张量约束下的阻抗张量分解方法，利用相位张量不受电场畸变影响的优点，分析得到合理的初始模型，提高了计算稳定性和计算效率，但仍然存在对初始值依赖性严重等问题。

笔者在求解GB分解法中的非线性超定方程组时,采用了粒子群优化算法，但传统的粒子群算法极易陷入局部极值当中，无法满足求解要求。为此，笔者在前人研究基础上，在传统粒子群算法中引入振荡环节和随机惯性权重等思想，形成基于随机惯性权重的二阶振荡粒子群算法，并利用该算法实现大地电磁阻抗张量分解，通过理论数据和实测资料的分解处理，验证了改进算法在阻抗张量分解中的适用性和有效性。

1 算法理论

1.1 GB分解原理

GB分解法针对三维/二维电性结构模型，将观测阻抗分解为区域阻抗张量和畸变阻抗张量，以此恢复未受畸变影响的区域阻抗张量，进而求得区域构造走向的角度以及电性主轴的特征。在不考虑磁场的感应型畸变的情况下，观测阻抗张量Z可以分解为：

Z=RCZ2DRT

C=gTSA

(1)

其中：

(2)

上述各矩阵带入式(2)中得到：

(3)

1.2 标准粒子群优化算法

粒子群最早是由Eberhaet等提出[9]，该算法概念源于鸟群 的觅食行为，通过 “鸟群”中 “每只鸟”之间的协作和信息共享来寻找最优解[10]。粒子群算法采用具有速度和位置属性的粒子来模拟鸟群中的鸟，速度代表移动的快慢，位置代表移动的方向。在N维空间中随机生成若干粒子，每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子“共享”，并找到其中最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子，根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解，来调整自己的速度和位置[10-14]。粒子群算法虽然是一个非常简单的算法，但其应用领域十分广泛，主要包括函数优化，约束优化，工程设计，神经网络训练等，具有简单高效且不依赖初始值等特点。

粒子群算法通常的数学描述为设在一个N维空间中，由m个粒子组成的种群X=(x1,…,xi,…,xN)，其中第i个粒子的位置为xi=(xi1,xi2,…,xiN)T, 其中速度为Vi=(vi1,vi2,…,vid,…,viN)T。它的个体极值为pi=(pi1,pi2,…,piN)T,种群的全局极值为pg=(pg1,pg2,…,pgN)T,按照追随当前最优粒子的原理，粒子xi将按照式(4)，式(5)改变自己的速度和位置。

vij(t+1)=ωvij(t)+c1r1(t)(pij(t)-

xij(t))+c2r2(t)(pgj(t)-xij(t))

(4)

xij(t+1)=xij(t)+vij(t+1)

(5)

其中：ω是惯性权重系数，其值较大时，全局寻优能力强，局部寻优能力弱；其值较小时，全局寻优能力弱，局部寻优能力强;j=1、2、…、N，i=1、2、…、m；m为种群规模；t为当前进化代数；r1、r2为分布于[0,1]之间的随机数；c1、c2为加速常数。从式(4)中可知，每个粒子由三部分组成：第一部分为“记忆”部分，表示上次速度大小和方向的影响；第二部分为“认知”部分，是从当前点指向粒子自身最好点的一个矢量，表示粒子自身的“思考”；第三部分为“社会”部分，是一个从当前点指向种群最好点的矢量，表示粒子之间的” 信息共享与协同合作”[14]。在使用粒子群算法时，需要重点设置待定参数变化范围和参数的每步迭代最大允许值，尤其是参数的每步迭代最大允许值，一般为变化范围的10%～20%，该值越小，收敛分辨率越高，但收敛速度慢，反之，该值越大，收敛速度越快，但易跳出全局最优值，故需要反复试验，以达到最佳优化效果。图1为粒子群算法流程图。

图1 粒子群算法流程图

1.3 改进的粒子群优化算法

标准的粒子群算法的性能很大程度上依赖于其参数(惯性权重、随机因子、加速常数、最大迭代步长等)，这就导致在应用过程中出现容易丧失种群多样性，过早陷入局部极值，精度不高等问题。为此国内、外学者提出了许多种改进方法，其中包括将其他优化算法思想引入到粒子群算法中(如基于模拟退火的粒子群优化算法[15])，另外通过改进标准粒子群中的速度和位置公式(式(4)和式(5))，提高算法性能也是行之有效的方法(如在速度公式中增加二次项，粒子位置采用高斯分布[16]等)。不管从什么角度对粒子群算法进行改进，都能一定程度上提高粒子群算法的性能。笔者在二阶振荡粒子群算法的基础上，将随机惯性权重引入该算法中，以此改善算法的全局收敛性，提高算法精度。

二阶振荡粒子算法由国内学者胡建秀等[17]提出，在速度方程中增加了二次项和振荡环节，使算法在前期具有较强全局收敛能力，呈振荡收敛，而在后期加强了局部搜索能力，呈渐进收敛。蒋丽等[18]通过改进在振荡环节中的参数及其取值范围，进而加强全局搜索能力，改进后的速度和位置公式如下：

vij(t+1)=ωvij(t)+c1r1(t)(pij(t)-

(1+ξ1)xij(t)+ξ2xij(t-1))+

c2r2(t)(pgj(t)-(1+ξ3)xij(t)+

ξ4xij(t-1))

(6)

xij(t+1)=xij(t)+vij(t+1)

(7)

惯性权重对粒子群算法来说是一个非常重要的参数，它的物理意义是对当前粒子速度的继承，标准粒子群算法的惯性权重，通常是一个常数或者随着迭代次数的增加成线性递减，但在实际应用中，有其不适用之处，容易导致收敛过慢或者陷入局部极值中。若将惯性权重设定为服从某种随机分布的随机数，这样就可以在一定程度上克服传统惯性权重的不足。如设ω'=ω-c1r1-c2r2，因为r1、r2是(0,1)之间的随机数，所以ω'实际上也是一个随机数。因此，式(6)和式(7)就可以改写为式(8)和式(9)。

vij(t+1)=(ω-c1r1(t)-c2r2(t))vij(t)+

c1r1(t)(pij(t)-(1+ξ1)xij(t)+

ξ2xij(t-1))+c2r2(t)(pgj(t)-

(1+ξ3)xij(t)+ξ4xij(t-1))

(8)

xij(t+1)=xij(t)+vij(t+1)

(9)

基于此，笔者将随机惯性权重引入到改进后的二阶振荡粒子群中，使算法的全局搜索能力和收敛速度得到一定程度地提高。

2 算法验证

2.1 仿真实验

利用两个测试函数来对上述改进方法进行试验比较，测试函数分别是Matyas函数和Sphere函数。为体现仿真实验的有效性，利用标准粒子群算法(算法1)，恒定惯性权重二阶振荡粒子算法(算法2)，线性惯性权重二阶振荡粒子群算法(算法3)，随机惯性权重二阶振荡粒子群算法(算法4)，对测试函数分别进行100次仿真计算实验，并统计分析计算后的数据，从最优值，平均值，标准差，平均耗时等对算法进行统计评估。

Matyas函数：

Sphere函数：

测试函数及算法参数如表1所示，Matyas函数维数为2，Sphere函数维度为50，除惯性权重，其余参数均一致。通过高低两个维度的测试函数来测试在各算法的性能表现，实验结果如表2所示，算法1在低维度函数中勉强达到求解要求，但在高维度函数中计算结果出现较大偏差；算法2的计算表现稍强于算法1，但不及算法3和算法4，平均耗时却少于算法3；算法3是计算结果最好的，但平均耗时不如算法2和算法4；算法4计算表现虽略逊于算法3，但平均耗时却强于算法3。总的来说，除算法1以外，其余算法均能达到测试计算要求。

表1 测试函数及算法参数

表2 实验结果

2.2 理论模型

从计算精度和计算效率两方面考虑，算法4相对更加适用于大地电磁阻抗张量分解。为验证其有效性，利用单个频点数据和三维模型数据进行分解试验。

2.2.1 单频点数据

利用算法4，经过反复试验，其中主要几个控制参数设置如下：种群中粒子数为24个，加速常数c1、c2均为2，最大迭代次数为2 000次，阈值为1e-25，最大迭代步长为取值范围的15%。 对上述畸变数据计算100次，计算出的统计结果为：扭曲因子为2.101°，剪切因子为24.68°，构造走向角为0.070°。所得结果如图2所示，其中横轴为计算次数，纵轴为计算结果。

从图2可看出， 采用改进算法计算出100次结果非常稳定，并且每次的结果都十分接近真实值，误差在允许范围之内，说明改进算法是能够准确地计算出模型的各个参数，进而可以达到理想的分解效果。

图2 单频点理论模型GB分解结果

2.2.2 三维正演模型

为进一步验证本文算法的分解效果，依据GB分解法要求的假设条件，建立一个三维/二维电性结构模型。其中背景区域电阻率为500 Ω·m，二维电性结构为一个无限延伸的低阻体，电阻率为100 Ω·m，局部三维块体电阻率为1 Ω·m。计算所用网格为44×45×60，计算频段为0.01 Hz～100 Hz，共20个频点。这里选取典型测点进行分析，测点与异常体相对位置如图3所示。图4是畸变前、后的视电阻率曲线，畸变前曲线表现出明显的二维特征，受到局部畸变影响后，视电阻率曲线出现平移，表现出强烈的三维性。

图3 三维/二维构造模型示意图

图4 畸变前后视电阻率曲线

利用改进算法进行分解处理，得到如图5所示的分解结果。从结果上看，在正演频段内，剪切角分布在0°到20°之间，扭曲角分布在0°到10°之间，电性主轴角分布在5°到33°之间。这表明利用基于粒子群算法的GB分解法所计算出的模型各参数均十分稳定且可靠。

图5 三维/二维模型GB分解结果

2.3 实际资料分解

BC87数据集(数据来源：https://www.mtnet.info/data/bc87/bc87.html)位于大不列颠哥伦比亚省东南部，大地电磁剖面长度为150 km，总计27个测点,观测周期为0.002 s～1 820 s.数据质量较好，但地质构造复杂，受到的畸变效应严重，为我们研究阻抗张量分解方法，恢复区域阻抗响应研究提供了有效的实测数据集。该地区的地质略图如图6所示，众多学者对该数据集进行了分析研究[17-18]。笔者选择其中位于Valhalla(瓦尔哈拉片麻岩)附近的lit000号测点(图6中红色测点)进行阻抗张量分解研究，并通过商业软件(MT-Pioneer)处理，得到分解结果真值，对比结果如图7所示，其中走向角(Strike)集中分布在0°到60°，扭曲因子(Twist)集中分布在-60°到30°，剪切因子(Shear)集中分布在0°到30°。利用本文算法得到的计算结果与商业软件的分解结果基本一致，但由于实测数据存在噪声干扰，对分解结果造成一定干扰，也就导致某些频点的分解结果与真实值存在差异。

图6 工区地质略图[20]

图7 实际资料分解结果

3 结论

在前人的研究基础上，将粒子群算法引入到大地电磁阻抗张量分解方法中，并对算法进行改进，提高其算法性能。利用测试函数对算法进行仿真实验，并经过对理论数据和实测数据的阻抗张量分解试验，主要得到以下结论：

1)单纯使用标准粒子算法无法满足计算要求，而该进后的随机惯性权重二阶振荡粒子群算法，无论在计算精度还是效率上都能满足需求。

2)GB分解算法对初始模型依赖性较强，初始模型选择对算法结果影响较大，如选择不当，会造成计算结果偏差较大，而粒子群算法随机生成初始值，在一定程度上降低模型初始值对分解效果的影响。再者，若初始值并不理想，改进算法中的随机惯性权重也有可能产生相对较大的值，加强全局搜索能力，避免陷入局部极值。

3)改进的粒子群算法，具有较强的随机性，控制参数较多，为达到满意的计算结果，可能需要反复试验，多次计算，但过程简单方便，易于实现。