于东凯 宋先海*② 张学强 赵素涛 蔡 伟
(①中国地质大学(武汉)地球物理与空间信息学院,湖北武汉 430074;②中国地质大学(武汉)湖北省地球内部多尺度成像重点实验室,湖北武汉 430074)
瑞雷波由纵波与横波的干涉形成,沿着自由表面传播。英国学者瑞雷首先发现并证明均匀半空间中瑞雷面波的存在[1]。瑞雷波勘探具有无损、高效、经济等特点,被视为未来最有希望的技术之一[2],现已被广泛用于浅地表、岩石圈横波速度结构信息[3-5]及地层品质因子[6]的计算、路基路面检测[7]、基岩面刻画[8]等。与常规勘探方法相比,瑞雷波勘探具有分辨率高、应用范围广、检测设备简单、检测速度快等优点[9]。近年来,诸多学者进行了相关研究,如瑞雷波反演[10-11]、瑞雷波频散研究[12-14]、瑞雷波数值模拟研究[15-17]、地形对瑞雷波传播的影响[18]等。
瑞雷波勘探可分为三个过程:面波数据采集、频散曲线提取和频散曲线反演[19-20]。拾取瑞雷波频散曲线后,利用瑞雷波频散特性可求得近地表横波速度,是许多工程的关键参数之一[21-22]。多数局部线性化优化方法被用来反演频散曲线,并在各种商业软件中占主导地位,如最速下降法、L-M算法等。然而,与大多数其他地球物理问题一样,瑞雷波频散曲线反演是高度非线性、多参数、多极值的,局部线性化方法易收敛到局部极小值,反演结果在很大程度上取决于初始模型的选择以及偏导数计算的准确性。全局最优化算法由于没有上述限制,被越来越多地应用于频散曲线反演,如Yamanaka等[23]借助遗传算法反演频散曲线研究地下结构信息; Beaty等[24]利用模拟退火算法反演多模式瑞雷波频散曲线; Shirazi等[25]通过人工神经网络算法反演频散曲线以获得路面结构; Song等[26-27]先后将模式识别算法、粒子群算法应用于频散曲线反演。传统算法,如遗传算法,在解决反演问题中具有较高的应用前景,然而往往存在计算量偏大、效率不高并且容易早熟收敛等缺陷,在反演中受到限制。
随着智能算法的发展,新型算法不断被应用于地球物理反演[28-31]。本文针对传统算法中出现的问题,引入一种新型智能算法——蚱蜢算法,对瑞雷波频散曲线进行反演以获取近地表横波速度剖面。对不同地质模型的理论频散曲线进行反演试验,并进行多次、多角度、深层次的探索,最后利用实测数据检验其实用性。
受蚱蜢不同时期觅食行为的启发,澳洲学者Saremi等[32]提出了一种新的仿生学算法——蚱蜢算法(Grasshopper Optimization Algorithm,GOA),算法模拟了蚱蜢在幼虫和成虫两个阶段搜索食物时所表现出的不同行为。蚱蜢是自然界中一种常见的昆虫,它们成群结队地出现,在幼虫时依靠强有力的双腿进行跳跃和移动,虽移动缓慢,但能够较为详尽地对所经区域是否存在食物进行判断;当其进入成虫期后长出翅膀,可在空中飞翔,移动变得十分迅速,且移动范围迅速扩大,搜索范围得到大幅度的扩大,但同时对食物的搜索精确度会得到一定程度的削弱。
蚱蜢算法将蚱蜢个体模拟为搜索算子,将成年蚱蜢大范围且迅速移动和幼虫期局部移动模拟为算法寻优过程中随迭代次数的增加所进行的全局探查与局部搜索两个阶段。在全局探查阶段,算子在模型内进行大范围的移动,获取全局信息,并最终固定为某局部区域;在局部搜索阶段,算子的搜索范围缩减至某局部区域,在该区域内不断迭代优化所得解,不断提高解的精度。算法主要包括两个要素:自身所处位置以及与其他算子的位置关系。其中,自身位置体现了算子当前的优劣,蚱蜢与其他算子的位置关系可表示为:搜索算子的移动受群体中其他算子的影响,当算子相隔较近(距离小于2.079个单位)时会出现彼此排斥现象,此时称蚱蜢处于排斥域范围内;当算子相隔大于该范围时又会出现吸引现象,称其位于吸引域内;当距离恰好相差2.079个单位时,彼此之间不会发生吸引或排斥现象,蚱蜢位于舒适域。算法通过不断调整自身位置及与其他算子的相对位置实现优化。需要说明的是,这两个阶段之间没有明显的界限,随着搜素区间逐步缩小,算法由全局阶段逐渐转至局部阶段。蚱蜢算法的优化原理描述如下。
蚱蜢算子xi受蚱蜢算子xj的影响因子定义为
(1)
假设i个蚱蜢算子当前位置为xi,k, 移动时会受到所有其他算子的影响,因此移动后的位置为
]}+xi,k
(2)
式中:i和j表示不同蚱蜢算子;NL表示算子总数;di,j表示第i个与第j个算子的空间距离;h和g均为评价算子受其他算子影响的参数,前者表示吸引力强度,后者是衡量空间距离的尺度,二者与影响因子F(di,j)紧密相关,不同组合会导致算子不同的行为习惯,算法中g通常取1.5,h取0.5;u(k)和l(k)分别表示算子的第k个分量在区间内的上、下限;xi,k表示当前蚱蜢算子所处位置;c1和c2为自适应的递减变量,c1与粒子群算法中惯性权重相似,作用是减少蚱蜢算子围绕目标值的移动,c2作用是收缩算子群体的舒适域、排斥域和吸引域,通过对自身的调节可在有限迭代次数中达到全局探查以及局部搜索的相对平衡。
总体来说,式(2)的前半部分实现了其他蚱蜢算子的约束作用,后半部分则对当前算子移动至全局最优解的行为趋势起到一定的控制作用。c1和c2虽在算法中所起作用不尽相同,但均为随迭代进行自适应衰减的变量,可统一表示如下
(3)
式中:Iiter代表当前迭代次数;NI为最大迭代次数;cmax为自适应参数最大值,本文cmax取值为1;cmin为最小值,文中取1.0×10-5。可根据解决问题的不同设定变量c1或c2的起始值。
本算法主要受三个参数控制:种群数量NL、迭代次数NI以及自适应参数c1和c2。种群数量与迭代次数对算法的影响相似:随着种群数量取值的增大,算法在变量空间中搜索会变得详尽,但同时也可能出现过多的冗余搜索,影响算法的效率;若种群数量取值过小,则会出现搜索不充分,无法找到全局最优解。因此,取值要根据所解决的具体问题来设定。自适应参数的取值与算法的寻优功能紧密相关,具体来说,为避免所得解陷入局部最优,设置如下机制:蚱蜢群体对某算子会产生排斥或吸引现象,其中对算子的排斥作用能有效避免局部最优解的产生,阻止算法陷入局部最优以及停滞不前。同时,由于自适应参数的引入,算法会缩小与迭代次数成比例的排斥面积,并不会对算法的收敛性造成过大影响,因此蚱蜢算子在开始阶段能有效避免局部最优解,在最后的局部搜索阶段也能改进优化所得解。
蚱蜢优化算法的基本流程如下。
(1)参数初始化。在模型搜索空间内随机生成蚱蜢算子的初始位置xi(i=1,2,…,NL),每个算子的位置可表示为xi=(xi,1,xi,2,…,xi,k,…,xi,D),这里D表示算子中包含的变量个数。
(2)根据蚱蜢算子所处位置,计算各个算子的目标函数值,并找出最优解及相应的算子位置。
(3)更新变量c1和c2,计算蚱蜢之间空间距离,并对所有蚱蜢算子遍历,根据式(2)更新空间位置; 同时,若新位置跃出上下边界,则给予约束。
(4)当达到最大搜索次数或解满足给定条件时,转至步骤(5); 否则,搜索次数增加1,转至步骤(2),进行下一次搜索。
(5)输出全局最优解及其对应的最优拟合值。
不同算法在处理问题上的能力存在差异,同一算法处理不同问题的能力也有高低。故在全局优化算法领域,通常使用一些特定的函数对算法进行测试,基于生成解的精度对算法的优化能力进行评价。本文主要研究将蚱蜢算法应用于多极值的频散曲线反演,多模态复合测试函数与反演中的目标函数存在一定的相似性,即具有多个局部最优解,因此选取两个含有多峰值的函数进行全局寻优能力测试。
(1)Griewank函数。如图1a所示(为方便显示,仅显示局部结果),该函数存在许多局部极小值,其数目与问题的维数有关,全局最小值0在(x1,x2,…,xD)=(0,0,…,0)处,因此该函数是典型的非线性多模态函数,具有广泛的搜索空间,通常被认为是优化算法很难处理的复杂多模态问题,在程序中搜索算子个数NL取500。其函数表达式为
xi∈[-600,600]
(4)
(2)Rastrigrin函数。如图1b所示,该函数是一个多峰值函数,在(x1,x2,…,xD)=(0,0,…,0)处取得全局最小值0, 在{xi∈(-5.12, 5.12),i=1,2,…,D}范围内大约有10D个局部极小值,此函数与Griewank函数类似,也是一种典型的非线性多模态函数,峰形呈高低起伏不定跳跃性地出现,所以很难优化查找到全局最优解。其函数表达式为
xi∈[-5.12, 5.12]
(5)
测试中,种群数量NL取100,迭代次数NI取700。表1列出了蚱蜢算法对两个测试函数的优化结果,该结果为5次试验的平均值,可看出算法的所得解均收敛到函数最小值附近,收敛误差接近于零,说明算法具备较强的全局寻优能力。图2a为Griewank函数的算法拟合误差随迭代次数的变化,可见迭代逐渐收敛并最终趋近于零,说明GOA法通过多次迭代能够有效提高解的精度。另外,图中也显示出拟合误差发生了不确定性变化,一方面是因为算法中引入了贪婪准则,保留当前最优解,算法在迭代之后所得解的精度未得到提高时会出现误差不变的情况; 另一方面,高排斥率导致蚱蜢算子在全局范围内发生大范围的移动,使迭代误差发生较大的改变。图2b所示为算法在多峰值、多模态测试函数Rastrigrin函数中的表现,图中每个点均代表一个蚱蜢算子在某次迭代中的搜索位置。可以看出:搜索过程中算子已基本覆盖所有目标区域,可见算法具有较高的搜索能力;图中红点表示全局最优解,算子在该点附近大量聚集,说明算法并未陷入局部最优解,而是逐步逼近最优解。另外,图中有“抱团”现象,即算子在某些区域聚集,原因可能是该点为局部最优点,蚱蜢算子在该处附近迭代次数较其他位置稍大。
图1 测试函数图
图2 蚱蜢算法中Griewank函数的迭代误差曲线(a)和Rastrigrin函数的算子位置图(b)
测试函数Griewank函数Rastrigrin函数理论极值00实际极值0.07880.0037算法所得解(4.17×10-4,7.01×10-3)(2.30×10-6,-1.35×10-6)
在瑞雷波勘探的应用中,依据面波频散资料推算工程场地的剪切波速度剖面时,通常认为地下为层状介质。同时,基阶波在面波频散资料中能量最高、最易观测,也应用最为广泛,故本文主要采用基阶模式频散曲线反演获取地下层状地质结构。夏江海等[6]的研究表明,与瑞雷波频散曲线特征关系最密切的是地层的横波速度(vS)和厚度(H),这两个变量可通过反演求取,地层密度、泊松比为已知参数(对频散曲线影响微弱,可通过工区信息大致估算)。另外,为更符合实际情况(通常没有足够的先验信息估计浅地表的横波速度和厚度准确范围),文中使用较大的搜索范围,在下文理论模型测试中,搜索区域的下限和上限设定为与真实值相差50%或更多。如前文所述,蚱蜢算法主要受到三个参数控制:NL、NI和c。反演程序中NL取值为10D;NI取值200;cmax和cmin分别取1和1.0×10-5。
频散曲线反演是一个求解目标函数最小值的优化过程,算法的好坏取决于能否在给定范围内充分搜索以找到反演参数的最优解。同时,由于反演中存在局部极值的干扰,算法能否跳出局部极值、找到全局最优解也是评价其质量的标准之一。目标函数是通过算法反演得到的模型能否精确解释观测资料所确定的,应用最为广泛的均方差函数是
(6)
(7)
目标函数具有以下几个特点:
①目标函数为M个频点(即频散曲线上的采样点)上的相速度共同作用。值得注意的是,采样数目过大并不会对反演过程产生显著影响,反而会影响反演效率[26]。由于所研究面波的频散曲线频率多在5~100Hz范围,试验发现频点数取25~30能平衡算法的高效性与采样的保真度。
针对目标函数的特点,蚱蜢算法在处理反演问题上具有独特优势,体现在算法能够保证目标函数不困于局部极值,而是进行全局范围的搜索;同时不断优化所得解,提高反演精度。算法引入排斥域、吸引域和舒适域机制:当某位置上多个算子之间距离过小时(在反演中表现为迭代误差相近)会相互排斥,算子在下次迭代中远离该位置,在算法初期需深入分析变量空间,避免受困于局部极值;当距离过大时认定进入吸引域,算子之间会产生吸引力,体现在中后期的“开发工作”中,随着算法不断运行,大量的非可行区间被排除,搜索范围会不断被压缩至最优解附近,算法对此区域不断开发,对所得解不断优化以求提高反演精度;当距离适当时,算子彼此之间不会吸引或排斥,此时位于舒适域。为更好地平衡“开发”(即全局探查)与“改进”(即局部搜索)两个过程,本算法配备了一个适应性缩减三个作用域面积的系数,随着迭代的进行,算法会逐渐由全局探查转为针对最优解附近区域的局部搜索。由此可知,通过群体获得的最佳解决方案是蚱蜢开发和改进的结果。因此,在解决一系列反演问题时,蚱蜢算法具有超越当前几种算法的潜力。
为检验蚱蜢算法在瑞雷波频散曲线反演中效果,使用两个理论模型进行测试,二者均为实际浅层工程勘察中经常遇到的地质模型。如表2所示,模型A代表横波速度递增的四层地质模型;模型B代表在两个较硬地层之间夹杂低速软弱地层,该地层剪切波速度低于上覆地层速度,常见的有路面结构等。
对理论模型进行蚱蜢算法反演试验,模型A的数值模拟及反演结果见图3和图4,模型B的数值模拟及反演结果见5和图6。图3a、图5a分别为两个模型波场正演得到的地震记录,以此模拟实测资料并作为频散曲线反演的依据。图4a、图6a为拟合频散曲线,图3b、图5b是通过高分辨率线性拉东变换提取的模拟地震记录得到的频散能量图。值得注意的是,由于含低速软弱夹层模型中频散曲线的特殊性,即基阶模式在高频范围内(该模型中为频率大于60Hz时)携带的能量并不是最高的,因此文中拾取的基阶波频散曲线也有一定的范围。
表2 两种地质模型参数及反演搜索范围
图3 模型A模拟地震记录(a)及频散能量图(b)
图4 基于蚱蜢算法的模型A反演结果
图5 模型B模拟地震记录(a)以及频散能量图(b)
图6 基于蚱蜢算法的模型B反演结果
可以看出,在没有足够先验信息的情况下,频散曲线吻合度较高,且模型中的参数均可被精确地反演和重建。就解的精确性而言,两个模型中参数的最大相对误差相对分别为2.50%、2.89%。由此可见,通过蚱蜢算法获得的地下横波速度与真实情况相差不大,说明本算法可以对理论瑞雷波频散曲线进行很好地反演,对频散曲线影响较大的参数(各层横波速度以及层厚度)均能被精确地反演和重建。
2.1.1 频散曲线联合反演测试
瑞雷波是频散且多阶的,某频率下的最小传播速度被称为基阶模式,其余按相速度大小依次称为一阶高阶模式、二阶高阶模式……。在瑞雷波频散曲线的应用中,由于基阶波能量较强且易观测等原因,通常使用基阶波频散曲线进行非线性反演。但在某些地形条件下,如文中的模型B,高阶模式的频散曲线在高频范围内携带了更多的能量,此时若只使用基阶波进行反演可能存在有效信息不足的问题。为进一步验证算法的反演能力,本次测试针对模型B,利用蚱蜢算法对基阶和高阶模式频散曲线进行联合反演。联合反演的频散曲线共N条,故目标函数相应变为
(8)
图7b所示为横波速度剖面的对比。表3为基阶波反演与联合反演的对比,可以看出联合反演结果较基阶波反演得到进一步优化,各参数基本接近真实值,反演误差趋于零,此结果进一步证实蚱蜢算法反演精度较高。
图7 基于蚱蜢算法的模型B联合反演结果
反演方式vS1m·s-1vS2m·s-1vS3m·s-1vS4m·s-1vS5m·s-1H1mH2mH3mH4m基阶波反演195.48147.14204.30304.01398.572.012.013.996.01联合反演 199.27149.58201.29301.74399.502.002.014.006.01
2.1.2 含噪测试
由于野外环境复杂,采集到的地震记录中不可避免含有噪声,含噪数据可能降低算法精度及稳定性,并影响反演结果。对一种新的反演算法需检验其抗噪能力。本文对理论频散曲线进行加噪处理(加入5%、10%及20%噪声),即对各频率下的相速度以理论值为基准在一定区间内进行随机程度的增减,以模拟实际情况中混杂了环境噪声的相速度,含噪频散曲线处理结果见图8。
图8a中杂乱不规则的蓝色、蓝绿色以及黑色实线所示为加入不同程度随机噪声的频散曲线,图8b为算法对含噪数据的反演结果。由图8可见,横波速度与真实模型之间有一定的偏离,三种不同噪声水平的数据反演结果最大相对误差分别增加到3.5%、9.5%和13.0%,由此可见噪声对反演解的准确性有一定的影响,但反演速度剖面均未实质性地偏离真实模型,因此由算法得到的反演结果仍然具有较高的可信度,说明蚱蜢算法对随机噪声具有一定的容忍度,用于反演瑞雷波频散曲线具有较高的稳定性。
2.1.3 搜索区间测试
前文设定的搜索范围是根据前期测井等数据对整个工区层位进行推断,并以此估算地震记录下对应各层横波速度及厚度的大致范围。倘若地质环境复杂,地层可能发生突变,则上述层位信息在当前地震记录下是不适用的,有必要设定更多数量的薄层以求准确模拟近地表地质情形。这是因为应用快速矢量传递算法计算频散曲线时需要输入各层横波速度和厚度,比如n层模型需输入(2n-1)个参数,故在反演之前有必要对模型层数进行设定。本文为进一步检测算法的反演能力,以模型A为例,将模型层数设定为7,各层层厚均设定为2m,且横波速度搜索范围均设定为100~800m/s,频散曲线反演是在给定的范围内找出最优模型,因此要预先设定反演中横波速度的搜索范围,反演结果如图9所示。由图可见,在先验信息有限、搜索空间更为广泛的前提下,反演过程耗时增加,反演结果仍十分接近真实值,并未发生显著变化,说明算法的全局搜索能力较强,反演结果具备较高的可信度。
图8 含不同噪声水平的模型A理论数据反演结果
图9 扩大搜索空间时基于蚱蜢算法的模型A反演结果
2.1.4 频散曲线缺失频段测试
在实地勘测中获得的瑞雷波频散曲线具有先天的不完整性,存在部分频段受干扰或缺失等情况。实际上,在松散的沉积物上施加低频源(如大锤)所产生的频谱往往会缺乏高频率数据。针对丢失频段的瑞雷波频散曲线反演时,首先应考虑目标层的敏感频率范围。本文针对模型A首先分析基阶频散曲线对各层横波速度的敏感度。从图10a可以看出,基阶波对各层横波速度变化敏感的频段随深度变化而变化,表层的敏感频段为高频,由于穿透深度的影响,深部层敏感频段的频率变得更低更窄,即面波的高频成分对浅部地层的横波速度较敏感,而低频部分对深度层和半空间的横波速度更敏感。基于此,截取0~50Hz范围的频散曲线进行基阶波频散曲线反演,通过对比横波速度剖面分析蚱蜢算法反演能力。从图10b可见,在缺失部分频段的情况下,多次试验结果虽然出现一定程度的波动,误差较之前也增大,但均未实质性地偏离真实值,且多次反演的均值与真实剖面吻合度较高,说明算法的反演能力值得信赖。
通过理论模型运算及四种测试全方位地分析了蚱蜢算法的反演能力,其中稳定性、收敛性和有效性是评价反演算法的重要指标。为避免单次反演的偶然性,对上述测试均进行20次独立反演运算,并取均值讨论。受篇幅限制,只列出速度递增模型A的反演结果统计表(表4)。可以看出:对于无噪数据,模型反演结果与真实值相符,相对误差和均方根误差均趋近于零,说明结果准确可靠; 对于含噪测试中含噪数据(10%随机噪声),模型的7个参数中相对误差最大不超过9.5%,均方差最大不超过14.32,这说明噪声会对反演精度造成一定的影响,但所得解与真实解相差并不大,精度能满足实际需求。原因在于蚱蜢算法引入舒适域、排斥域和吸引域三个作用域,算子会受到群体中其余所有算子的共同作用,全局寻优能力更强,不易出现早熟收敛,能够有效避免其他算法中由于随机化个体单一、受当前最优解影响过大导致的全局寻优能力不足等问题。
图10 基于蚱蜢算法的模型A反演结果(缺失频段)
参数真实值不含噪声含10%随机噪声反演均值相对误差/%标准差反演均值相对误差/%标准差vS1/(m·s-1)200.00203.841.921.57208.114.0614.32vS2/(m·s-1)250.00252.400.963.02257.402.965.60vS3/(m·s-1)300.00303.241.082.14290.243.2510.50vS4/(m·s-1)400.00396.810.803.26394.011.507.86H1/(m)2.001.952.501.202.042.001.02H2/(m)2.002.031.500.372.199.500.79H3/(m)4.004.061.500.414.082.000.85
另外,由搜索区间测试中的收敛曲线可以看出,无论搜索区间的上、下边界是否靠近真实值,算法的最优拟合值在初期的迭代过程中均迅速下降,说明算法在此阶段迅速排除搜索空间内大量非可行区间,类比于蚱蜢在成虫期的觅食中能够迅速排除非食物源区域。在随后迭代过程中误差下降趋势变缓,且在后期趋近于零,说明算法的反演过程在此之前已基本完成,各个参数已被精准地推算。值得注意的是,任何一种群智能算法都会出现多次迭代中解的精度没有得到提高的情况,这是因为在某几次或多次迭代过程中,种群中没有找到更优个体,根据贪婪准则保留当前最优解,解的精度没有得到提高,故而误差不变。
同时,频散曲线反演是一个非线性、多极值问题,不同的模型参数组合也可能获得相似的拟合值,因此还需要分析每次反演所得解中模型参数的分布概率,以验证反演结果的有效性。图11所示为模型A中各个参数在真实值附近的分布概率。图11a所示为第一层横波速度vS1在20次反演结果中的分布概率,真实值为200m/s,反演所得解在真实值±10m/s范围内的概率接近0.80。由图11可见,其余模型参数在其各自真实值附近(速度±10m/s,厚度±0.1m)的分布概率分别为0.75、0.60、0.70、0.65、0.70和0.65。可见反演所得模型参数均大致分布在真实值附近,进一步说明蚱蜢算法是一种可靠、准确、收敛、有效的反演算法,能够有效地应用于瑞雷波频散曲线的反演。
图11 模型A反演参数的分布概率直方图
为验证蚱蜢算法在处理地球物理反演问题中的适用性,将该算法与粒子群算法(Particle Swarm Optimization,PSO)进行比较。粒子群算法是由Kennedy等[34]提出,后得到众多学者的改进和优化[35],已成为当今最流行、最成熟的自然启发优化算法之一,并在地震勘探中已获得广泛应用[36-38]。前人已成功将其应用于频散曲线反演[27],本文就基于新型算法与粒子群算法的频散曲线反演结果进行对比,证明前者具有较高的可信度。
这两种算法具备如下几个相同点:均为全局优化算法,算法在全局解空间内进行搜索,并将搜索重点集中在性能高的部分;都是随机搜索算法,且具备记忆性能,最优解都被保留下来;具有隐含并行性搜索特性,搜索过程均从解空间内的一个集合开始。但蚱蜢算法较粒子群算法而言,在处理地球物理反演问题上具有不易陷入局部极值的优势。后者根据粒子本身历史最优位置和当前群体内的最优算子的位置决定下一步的移动方向,因此可能会出现如下现象:由于目标函数存在大量的局部极小值,若最优算子进入某个极值位置,则算法接下来会在该位置附近多次移动,反演效率不高。而蚱蜢算法中算子会受到群体内其他所有算子的共同作用,体现在程序中算子每一步的移动方向不仅受当前最优算子的影响,也会受到群体内其余算子的共同作用,一定程度上降低了上述现象产生的可能性。同时,当算法运行至中后期时,算法通过缩小与迭代次数成比例的排斥面积,会减少排斥现象的产生,故不会对算法的收敛性造成过大的影响。
以模型B为例,分别采用上述两种算法对模型B不含噪声的理论数据进行10次独立反演测试。为排除其他因素的干扰,采用相同的群体规模、搜索空间和迭代次数,M取值30。测试分析一次反演的结果而不是10次反演迭代误差的均值,目的是更准确直观地分析算法的性能。反演结果如图12所示,由图可以看出,如果迭代截至50~100次(图12a),迭代误差的整体变化情形与算法优化所得反演解的能力基本一致,两种算法效果不相上下,说明这两种算法均具备较强的收敛能力和收缩搜索区间性能;但迭代进入中后期(图12b),即两种算法运行100~200次,蚱蜢算法比粒子群算法收敛速度更快,且收敛精度也有一定提高,说明蚱蜢算法在处理地球物理反演问题时具备较高的寻优性能。
图12 模型B不含噪声时的GOA与PSO迭代误差对比
前文充分验证了蚱蜢算法在理论瑞雷波频散曲线反演中可行且有效,所得解精度高。为进一步检验算法的适用性,利用蚱蜢算法对美国怀俄明某地获得的地震数据进行了分析。图13a所示为采集的地震记录,数据采集采用了48个8Hz垂直分量检波器,道间距为0.9m,最小炮检距为0.9m,震源采用锤击震源。图13b所示为对应的频散能量图,通过采集能量最大点提取基阶波频散曲线。与反演理论模型类似,固定地层的密度及泊松比,只反演横波速度vS与厚度H两个参量。根据测井数据将地层大致划分为5个层位,并粗略估计搜索范围。反演模型的搜索空间、泊松比以及密度等参数设置见表5。
图14a所示为算法所得解与实测数据的基阶波频散曲线拟合情况,通过拾取频散能量图中能量最大值(图13b中圆点)得到不同频率下的相速度,可以看出频散曲线拟合情况较好。文中采用基阶波反演的原因在于能量图中出现较为严重的“模式接吻”以及“模式跳跃”现象,在拾取高阶频散曲线的过程中很容易出现模式误判,从而影响反演结果,若采用基阶波反演精度会更高。图14b为两种算法反演时迭代拟合误差随迭代次数的变化曲线,为保证反演的精准性、降低因单次反演中出现的偶然误差对结果造成影响,针对野外数据进行了多次反演,舍弃异常的反演解,并对剩余解进行均值化处理,因此迭代误差曲线较图12更为平滑。图14c为横波速度剖面对比,可见通过有限频段的频散曲线反演所得速度模型与测井资料吻合较好,尤其是浅地表部分。相较于PSO算法,GOA算法能够准确地分辨厚度为3~5m的软弱夹层。值得注意的是,较前文测试,此次实测资料反演的信息量有限,与模型试算结果相比反演效果并不突出,但借助有限的频散信息所得的反演结果能基本满足实际工作需求。另外,较深地层反演速度剖面与测井资料出现了一定偏离,这是由于采用锤击震源,瑞雷波频散能量在低频段分辨率低,导致提取误差增大;另外一个原因是算法只设定5个层位,较深的地层被认定为均匀介质,即横波速度没有变化。
表5 GOA反演中模型搜索范围及模型参数
图13 怀俄明某地区地震记录(a)与频散能量图(b)
图14 怀俄明某地区地震数据蚱蜢算法反演结果
本文介绍了一种全新的仿生学算法——蚱蜢算法,并通过两个测试函数验证其全局搜索能力,证明了其优异的收敛能力以及全局寻优能力。将蚱蜢算法应用于瑞雷波频散曲线反演,通过设置较大的搜索范围以模拟实际工作中缺少足够先验信息的情况,并进行多角度、深层次的探索试验,深入分析了算法的反演能力。最后通过实测数据检测了本算法的适用性。理论数据和实测资料反演的结果表明:
(1)蚱蜢算法具有较高的全局寻优能力以及抗局部最优解能力,所得解与函数最小解十分接近,说明这是一种有效的全局优化算法;
(2)蚱蜢算法可应用于瑞雷波频散曲线反演,与传统算法相比,反演过程中目标函数曲线收敛更加稳定,最终迭代误差更接近于零。证明基于蚱蜢算法的瑞雷波频散曲线反演具有较高的精度,是一种简单、稳健的反演算法。
同时,蚱蜢算法也存在一定的缺陷,比如算法的收敛过程相对较长。在今后的研究中,可与其他算法如遗传算法、粒子群算法或人工神经网络算法相结合,探讨如何进一步提高算法的收敛效率。