② 中国地质大学(武汉)湖北省地球内部多尺度成像重点实验室, 湖北武汉 430074
② Hubei Subsurface Multi-scale Imaging Laboratory(SMIL), China University of Geosciences(Wuhan), Wuhan, Hubei 430074, China
瑞雷波由纵波与横波的干涉形成,沿着自由表面传播。英国学者瑞雷首先发现并证明均匀半空间中瑞雷面波的存在[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]。本文针对传统算法中出现的问题,引入一种新型智能算法——蚱蜢算法,对瑞雷波频散曲线进行反演以获取近地表横波速度剖面。对不同地质模型的理论频散曲线进行反演试验,并进行多次、多角度、深层次的探索,最后利用实测数据检验其实用性。
1 蚱蜢算法基本原理及验证 1.1 基本原理受蚱蜢不同时期觅食行为的启发,澳洲学者Saremi等[32]提出了一种新的仿生学算法——蚱蜢算法(Grasshopper Optimization Algorithm,GOA),算法模拟了蚱蜢在幼虫和成虫两个阶段搜索食物时所表现出的不同行为。蚱蜢是自然界中一种常见的昆虫,它们成群结队地出现,在幼虫时依靠强有力的双腿进行跳跃和移动,虽移动缓慢,但能够较为详尽地对所经区域是否存在食物进行判断;当其进入成虫期后长出翅膀,可在空中飞翔,移动变得十分迅速,且移动范围迅速扩大,搜索范围得到大幅度的扩大,但同时对食物的搜索精确度会得到一定程度的削弱。
蚱蜢算法将蚱蜢个体模拟为搜索算子,将成年蚱蜢大范围且迅速移动和幼虫期局部移动模拟为算法寻优过程中随迭代次数的增加所进行的全局探查与局部搜索两个阶段。在全局探查阶段,算子在模型内进行大范围的移动,获取全局信息,并最终固定为某局部区域;在局部搜索阶段,算子的搜索范围缩减至某局部区域,在该区域内不断迭代优化所得解,不断提高解的精度。算法主要包括两个要素:自身所处位置以及与其他算子的位置关系。其中,自身位置体现了算子当前的优劣,蚱蜢与其他算子的位置关系可表示为:搜索算子的移动受群体中其他算子的影响,当算子相隔较近(距离小于2.079个单位)时会出现彼此排斥现象,此时称蚱蜢处于排斥域范围内;当算子相隔大于该范围时又会出现吸引现象,称其位于吸引域内;当距离恰好相差2.079个单位时,彼此之间不会发生吸引或排斥现象,蚱蜢位于舒适域。算法通过不断调整自身位置及与其他算子的相对位置实现优化。需要说明的是,这两个阶段之间没有明显的界限,随着搜素区间逐步缩小,算法由全局阶段逐渐转至局部阶段。蚱蜢算法的优化原理描述如下。
蚱蜢算子xi受蚱蜢算子xj的影响因子定义为
$ F({d_{i, j}}) = \left[ {h{{\rm{e}}^{ - \frac{{{d_{i, j}}}}{g}}} - {{\rm{e}}^{ - {d_{i, j}}}}} \right]\frac{{{x_{j, k}} - {x_{i, k}}}}{{{d_{i, j}}}} $ | (1) |
假设i个蚱蜢算子当前位置为xi, k,移动时会受到所有其他算子的影响,因此移动后的位置为
$ {\hat x_{i, k}} = {c_1}\{ \sum\limits_{j = 1, j \ne i}^{{N_{\rm{L}}}} {\left[ {{c_2}\frac{{u\left( k \right) - l\left( k \right)}}{2}F\left( {{d_{i, j}}} \right)} \right]\} + {x_{i, 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虽在算法中所起作用不尽相同,但均为随迭代进行自适应衰减的变量,可统一表示如下
$ c = {c_{{\rm{max}}}} - {I_{{\rm{iter}}}}\frac{{{c_{{\rm{max}}}} - {c_{{\rm{min}}}}}}{{{N_{\rm{I}}}}} $ | (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.2 算法测试不同算法在处理问题上的能力存在差异,同一算法处理不同问题的能力也有高低。故在全局优化算法领域,通常使用一些特定的函数对算法进行测试,基于生成解的精度对算法的优化能力进行评价。本文主要研究将蚱蜢算法应用于多极值的频散曲线反演,多模态复合测试函数与反演中的目标函数存在一定的相似性,即具有多个局部最优解,因此选取两个含有多峰值的函数进行全局寻优能力测试。
(1) Griewank函数。如图 1a所示(为方便显示,仅显示局部结果),该函数存在许多局部极小值,其数目与问题的维数有关,全局最小值0在(x1, x2, …, xD)=(0, 0, …, 0)处,因此该函数是典型的非线性多模态函数,具有广泛的搜索空间,通常被认为是优化算法很难处理的复杂多模态问题,在程序中搜索算子个数NL取500。其函数表达式为
$ \begin{array}{l} f(x_i) = \sum\limits_{i = 1}^D {\frac{{x_i^2}}{{4000}} - \mathop \Pi \limits_{i = 1}^D {\rm{cos}}\frac{{{x_i}}}{{\sqrt i }} + 1} \\ \;\;\;\;\;\;\;\;\;\;\;{x_i} \in \left[ { - 600, 600} \right] \end{array} $ | (4) |
(2) Rastrigrin函数。如图 1b所示,该函数是一个多峰值函数,在(x1, x2, …, xD)=(0, 0, …, 0)处取得全局最小值0,在{xi∈(-5.12, 5.12), i=1, 2, …, D}范围内大约有10D个局部极小值,此函数与Griewank函数类似,也是一种典型的非线性多模态函数,峰形呈高低起伏不定跳跃性地出现,所以很难优化查找到全局最优解。其函数表达式为
$ \begin{array}{l} f({x_i}) = \sum\limits_{i = 1}^D {[x_i^2 - 10{\rm{cos}}(2{\rm{ \mathsf{ π} }}{x_i}) + 10]} \\ \;\;\;\;\;\;\;\;\;\;\;{x_i} \in \left[ { - 5.12, {\rm{ }}5.12} \right] \end{array} $ | (5) |
测试中,种群数量NL取100,迭代次数NI取700。表 1列出了蚱蜢算法对两个测试函数的优化结果,该结果为5次试验的平均值,可看出算法的所得解均收敛到函数最小值附近,收敛误差接近于零,说明算法具备较强的全局寻优能力。图 2a为Griewank函数的算法拟合误差随迭代次数的变化,可见迭代逐渐收敛并最终趋近于零,说明GOA法通过多次迭代能够有效提高解的精度。另外,图中也显示出拟合误差发生了不确定性变化,一方面是因为算法中引入了贪婪准则,保留当前最优解,算法在迭代之后所得解的精度未得到提高时会出现误差不变的情况;另一方面,高排斥率导致蚱蜢算子在全局范围内发生大范围的移动,使迭代误差发生较大的改变。图 2b所示为算法在多峰值、多模态测试函数Rastrigrin函数中的表现,图中每个点均代表一个蚱蜢算子在某次迭代中的搜索位置。可以看出:搜索过程中算子已基本覆盖所有目标区域,可见算法具有较高的搜索能力;图中红点表示全局最优解,算子在该点附近大量聚集,说明算法并未陷入局部最优解,而是逐步逼近最优解。另外,图中有“抱团”现象,即算子在某些区域聚集,原因可能是该点为局部最优点,蚱蜢算子在该处附近迭代次数较其他位置稍大。
在瑞雷波勘探的应用中,依据面波频散资料推算工程场地的剪切波速度剖面时,通常认为地下为层状介质。同时,基阶波在面波频散资料中能量最高、最易观测,也应用最为广泛,故本文主要采用基阶模式频散曲线反演获取地下层状地质结构。夏江海等[6]的研究表明,与瑞雷波频散曲线特征关系最密切的是地层的横波速度(vS)和厚度(H),这两个变量可通过反演求取,地层密度、泊松比为已知参数(对频散曲线影响微弱,可通过工区信息大致估算)。另外,为更符合实际情况(通常没有足够的先验信息估计浅地表的横波速度和厚度准确范围),文中使用较大的搜索范围,在下文理论模型测试中,搜索区域的下限和上限设定为与真实值相差50%或更多。如前文所述,蚱蜢算法主要受到三个参数控制:NL、NI和c。反演程序中NL取值为10D;NI取值200;cmax和cmin分别取1和1.0×10-5。
频散曲线反演是一个求解目标函数最小值的优化过程,算法的好坏取决于能否在给定范围内充分搜索以找到反演参数的最优解。同时,由于反演中存在局部极值的干扰,算法能否跳出局部极值、找到全局最优解也是评价其质量的标准之一。目标函数是通过算法反演得到的模型能否精确解释观测资料所确定的,应用最为广泛的均方差函数是
$ F = \sqrt {\frac{1}{M}\sum\limits_{i = 1}^M {{{\left[ {\mathit{\boldsymbol{v}}_{\rm{R}}^{{\rm{obs}}}\left( i \right) - \mathit{\boldsymbol{v}}_{\rm{R}}^{{\rm{cal}}}\left( i \right)} \right]}^2}} } $ | (6) |
式中:vRobs是实测面波的相速度;vRcal是快速矢量传递算法正演模拟得到的相速度[33]; M是频点数。一般情况下,vRobs和vRcal均为一维向量,代表基阶波频散曲线,向量长度等于M。借助式(6)反演频散曲线, 可得到准确的地层参数,但各频率下的相速度较大且采样的频点数较多,使用该公式计算量极大,影响算法的反演效率。因此引入一种新型的目标函数
$ \left\{ \begin{array}{l} F = \frac{{100}}{M}\sum\limits_{i = 1}^M {{F_i}} \\ {F_i} = \frac{{\left| {v_{\rm{R}}^{{\rm{obs}}}(i) - v_{\rm{R}}^{{\rm{cal}}}\left( i \right)} \right|}}{{v_{\rm{R}}^{{\rm{obs}}}(i)}} \end{array} \right. $ | (7) |
目标函数具有以下几个特点:
① 目标函数为M个频点(即频散曲线上的采样点)上的相速度共同作用。值得注意的是,采样数目过大并不会对反演过程产生显著影响,反而会影响反演效率[26]。由于所研究面波的频散曲线频率多在5~100Hz范围,试验发现频点数取25~30能平衡算法的高效性与采样的保真度。
② 目标函数为多参数、多极值函数。式中vRcal和vRobs两个向量中所有参数共同作用于函数,使函数产生大量的局部极值解。
③ 目标函数与反演模型的参数并不是直接相关的。算法通过调节模型参数改变相速度vRcal间接影响函数值,可更直观地研究反演过程中迭代误差的变化,在设定函数时将迭代误差结果扩大100倍。
针对目标函数的特点,蚱蜢算法在处理反演问题上具有独特优势,体现在算法能够保证目标函数不困于局部极值,而是进行全局范围的搜索;同时不断优化所得解,提高反演精度。算法引入排斥域、吸引域和舒适域机制:当某位置上多个算子之间距离过小时(在反演中表现为迭代误差相近)会相互排斥,算子在下次迭代中远离该位置,在算法初期需深入分析变量空间,避免受困于局部极值;当距离过大时认定进入吸引域,算子之间会产生吸引力,体现在中后期的“开发工作”中,随着算法不断运行,大量的非可行区间被排除,搜索范围会不断被压缩至最优解附近,算法对此区域不断开发,对所得解不断优化以求提高反演精度;当距离适当时,算子彼此之间不会吸引或排斥,此时位于舒适域。为更好地平衡“开发”(即全局探查)与“改进”(即局部搜索)两个过程,本算法配备了一个适应性缩减三个作用域面积的系数,随着迭代的进行,算法会逐渐由全局探查转为针对最优解附近区域的局部搜索。由此可知,通过群体获得的最佳解决方案是蚱蜢开发和改进的结果。因此,在解决一系列反演问题时,蚱蜢算法具有超越当前几种算法的潜力。
2.1 理论模型反演为检验蚱蜢算法在瑞雷波频散曲线反演中效果,使用两个理论模型进行测试,二者均为实际浅层工程勘察中经常遇到的地质模型。如表 2所示,模型A代表横波速度递增的四层地质模型;模型B代表在两个较硬地层之间夹杂低速软弱地层,该地层剪切波速度低于上覆地层速度,常见的有路面结构等。
对理论模型进行蚱蜢算法反演试验,模型A的数值模拟及反演结果见图 3和图 4,模型B的数值模拟及反演结果见5和图 6。图 3a、图 5a分别为两个模型波场正演得到的地震记录,以此模拟实测资料并作为频散曲线反演的依据。图 4a、图 6a为拟合频散曲线,图 3b、图 5b是通过高分辨率线性拉东变换提取的模拟地震记录得到的频散能量图。值得注意的是,由于含低速软弱夹层模型中频散曲线的特殊性,即基阶模式在高频范围内(该模型中为频率大于60Hz时)携带的能量并不是最高的,因此文中拾取的基阶波频散曲线也有一定的范围。
可以看出,在没有足够先验信息的情况下,频散曲线吻合度较高,且模型中的参数均可被精确地反演和重建。就解的精确性而言,两个模型中参数的最大相对误差相对分别为2.50%、2.89%。由此可见,通过蚱蜢算法获得的地下横波速度与真实情况相差不大,说明本算法可以对理论瑞雷波频散曲线进行很好地反演,对频散曲线影响较大的参数(各层横波速度以及层厚度)均能被精确地反演和重建。
2.1.1 频散曲线联合反演测试瑞雷波是频散且多阶的,某频率下的最小传播速度被称为基阶模式,其余按相速度大小依次称为一阶高阶模式、二阶高阶模式……。在瑞雷波频散曲线的应用中,由于基阶波能量较强且易观测等原因,通常使用基阶波频散曲线进行非线性反演。但在某些地形条件下,如文中的模型B,高阶模式的频散曲线在高频范围内携带了更多的能量,此时若只使用基阶波进行反演可能存在有效信息不足的问题。为进一步验证算法的反演能力,本次测试针对模型B,利用蚱蜢算法对基阶和高阶模式频散曲线进行联合反演。联合反演的频散曲线共N条,故目标函数相应变为
$ \left\{ \begin{array}{l} F = \frac{{100}}{N}\sum\limits_{i = 1}^N {\frac{{\sum\limits_{j = 1}^{{M_i}} {F\left( {i, j} \right)} }}{{{M_i}}}} \\ F\left( {i, j} \right) = \frac{{\left| {\mathit{\boldsymbol{v}}_{\rm{R}}^{{\rm{obs}}}(i, j) - \mathit{\boldsymbol{v}}_{\rm{R}}^{{\rm{cal}}}\left( {i, j} \right)} \right|}}{{\mathit{\boldsymbol{v}}_{\rm{R}}^{{\rm{obs}}}(i, j)}} \end{array} \right. $ | (8) |
式中vRobs和vRcal均为二维矩阵,矩阵中每列向量都代表一种模式频散曲线下对应的相速度。
图 7b所示为横波速度剖面的对比。表 3为基阶波反演与联合反演的对比,可以看出联合反演结果较基阶波反演得到进一步优化,各参数基本接近真实值,反演误差趋于零,此结果进一步证实蚱蜢算法反演精度较高。
由于野外环境复杂,采集到的地震记录中不可避免含有噪声,含噪数据可能降低算法精度及稳定性,并影响反演结果。对一种新的反演算法需检验其抗噪能力。本文对理论频散曲线进行加噪处理(加入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所示。由图可见,在先验信息有限、搜索空间更为广泛的前提下,反演过程耗时增加,反演结果仍十分接近真实值,并未发生显著变化,说明算法的全局搜索能力较强,反演结果具备较高的可信度。
在实地勘测中获得的瑞雷波频散曲线具有先天的不完整性,存在部分频段受干扰或缺失等情况。实际上,在松散的沉积物上施加低频源(如大锤)所产生的频谱往往会缺乏高频率数据。针对丢失频段的瑞雷波频散曲线反演时,首先应考虑目标层的敏感频率范围。本文针对模型A首先分析基阶频散曲线对各层横波速度的敏感度。从图 10a可以看出,基阶波对各层横波速度变化敏感的频段随深度变化而变化,表层的敏感频段为高频,由于穿透深度的影响,深部层敏感频段的频率变得更低更窄,即面波的高频成分对浅部地层的横波速度较敏感,而低频部分对深度层和半空间的横波速度更敏感。基于此,截取0~50Hz范围的频散曲线进行基阶波频散曲线反演,通过对比横波速度剖面分析蚱蜢算法反演能力。从图 10b可见,在缺失部分频段的情况下,多次试验结果虽然出现一定程度的波动,误差较之前也增大,但均未实质性地偏离真实值,且多次反演的均值与真实剖面吻合度较高,说明算法的反演能力值得信赖。
通过理论模型运算及四种测试全方位地分析了蚱蜢算法的反演能力,其中稳定性、收敛性和有效性是评价反演算法的重要指标。为避免单次反演的偶然性,对上述测试均进行20次独立反演运算,并取均值讨论。受篇幅限制,只列出速度递增模型A的反演结果统计表(表 4)。可以看出:对于无噪数据,模型反演结果与真实值相符,相对误差和均方根误差均趋近于零,说明结果准确可靠;对于含噪测试中含噪数据(10%随机噪声),模型的7个参数中相对误差最大不超过9.5%,均方差最大不超过14.32,这说明噪声会对反演精度造成一定的影响,但所得解与真实解相差并不大,精度能满足实际需求。原因在于蚱蜢算法引入舒适域、排斥域和吸引域三个作用域,算子会受到群体中其余所有算子的共同作用,全局寻优能力更强,不易出现早熟收敛,能够有效避免其他算法中由于随机化个体单一、受当前最优解影响过大导致的全局寻优能力不足等问题。
另外,由搜索区间测试中的收敛曲线可以看出,无论搜索区间的上、下边界是否靠近真实值,算法的最优拟合值在初期的迭代过程中均迅速下降,说明算法在此阶段迅速排除搜索空间内大量非可行区间,类比于蚱蜢在成虫期的觅食中能够迅速排除非食物源区域。在随后迭代过程中误差下降趋势变缓,且在后期趋近于零,说明算法的反演过程在此之前已基本完成,各个参数已被精准地推算。值得注意的是,任何一种群智能算法都会出现多次迭代中解的精度没有得到提高的情况,这是因为在某几次或多次迭代过程中,种群中没有找到更优个体,根据贪婪准则保留当前最优解,解的精度没有得到提高,故而误差不变。
同时,频散曲线反演是一个非线性、多极值问题,不同的模型参数组合也可能获得相似的拟合值,因此还需要分析每次反演所得解中模型参数的分布概率,以验证反演结果的有效性。图 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。可见反演所得模型参数均大致分布在真实值附近,进一步说明蚱蜢算法是一种可靠、准确、收敛、有效的反演算法,能够有效地应用于瑞雷波频散曲线的反演。
为验证蚱蜢算法在处理地球物理反演问题中的适用性,将该算法与粒子群算法(Particle Swarm Optimization,PSO)进行比较。粒子群算法是由Kennedy等[34]提出,后得到众多学者的改进和优化[35],已成为当今最流行、最成熟的自然启发优化算法之一,并在地震勘探中已获得广泛应用[36-38]。前人已成功将其应用于频散曲线反演[27],本文就基于新型算法与粒子群算法的频散曲线反演结果进行对比,证明前者具有较高的可信度。
这两种算法具备如下几个相同点:均为全局优化算法,算法在全局解空间内进行搜索,并将搜索重点集中在性能高的部分;都是随机搜索算法,且具备记忆性能,最优解都被保留下来;具有隐含并行性搜索特性,搜索过程均从解空间内的一个集合开始。但蚱蜢算法较粒子群算法而言,在处理地球物理反演问题上具有不易陷入局部极值的优势。后者根据粒子本身历史最优位置和当前群体内的最优算子的位置决定下一步的移动方向,因此可能会出现如下现象:由于目标函数存在大量的局部极小值,若最优算子进入某个极值位置,则算法接下来会在该位置附近多次移动,反演效率不高。而蚱蜢算法中算子会受到群体内其他所有算子的共同作用,体现在程序中算子每一步的移动方向不仅受当前最优算子的影响,也会受到群体内其余算子的共同作用,一定程度上降低了上述现象产生的可能性。同时,当算法运行至中后期时,算法通过缩小与迭代次数成比例的排斥面积,会减少排斥现象的产生,故不会对算法的收敛性造成过大的影响。
以模型B为例,分别采用上述两种算法对模型B不含噪声的理论数据进行10次独立反演测试。为排除其他因素的干扰,采用相同的群体规模、搜索空间和迭代次数,M取值30。测试分析一次反演的结果而不是10次反演迭代误差的均值,目的是更准确直观地分析算法的性能。反演结果如图 12所示,由图可以看出,如果迭代截至50~100次(图 12a),迭代误差的整体变化情形与算法优化所得反演解的能力基本一致,两种算法效果不相上下,说明这两种算法均具备较强的收敛能力和收缩搜索区间性能;但迭代进入中后期(图 12b),即两种算法运行100~200次,蚱蜢算法比粒子群算法收敛速度更快,且收敛精度也有一定提高,说明蚱蜢算法在处理地球物理反演问题时具备较高的寻优性能。
前文充分验证了蚱蜢算法在理论瑞雷波频散曲线反演中可行且有效,所得解精度高。为进一步检验算法的适用性,利用蚱蜢算法对美国怀俄明某地获得的地震数据进行了分析。图 13a所示为采集的地震记录,数据采集采用了48个8Hz垂直分量检波器,道间距为0.9m,最小炮检距为0.9m,震源采用锤击震源。图 13b所示为对应的频散能量图,通过采集能量最大点提取基阶波频散曲线。与反演理论模型类似,固定地层的密度及泊松比,只反演横波速度vS与厚度H两个参量。根据测井数据将地层大致划分为5个层位,并粗略估计搜索范围。反演模型的搜索空间、泊松比以及密度等参数设置见表 5。
图 14a所示为算法所得解与实测数据的基阶波频散曲线拟合情况,通过拾取频散能量图中能量最大值(图 13b中圆点)得到不同频率下的相速度,可以看出频散曲线拟合情况较好。文中采用基阶波反演的原因在于能量图中出现较为严重的“模式接吻”以及“模式跳跃”现象,在拾取高阶频散曲线的过程中很容易出现模式误判,从而影响反演结果,若采用基阶波反演精度会更高。图 14b为两种算法反演时迭代拟合误差随迭代次数的变化曲线,为保证反演的精准性、降低因单次反演中出现的偶然误差对结果造成影响,针对野外数据进行了多次反演,舍弃异常的反演解,并对剩余解进行均值化处理,因此迭代误差曲线较图 12更为平滑。图 14c为横波速度剖面对比,可见通过有限频段的频散曲线反演所得速度模型与测井资料吻合较好,尤其是浅地表部分。相较于PSO算法,GOA算法能够准确地分辨厚度为3~5m的软弱夹层。值得注意的是,较前文测试,此次实测资料反演的信息量有限,与模型试算结果相比反演效果并不突出,但借助有限的频散信息所得的反演结果能基本满足实际工作需求。另外,较深地层反演速度剖面与测井资料出现了一定偏离,这是由于采用锤击震源,瑞雷波频散能量在低频段分辨率低,导致提取误差增大;另外一个原因是算法只设定5个层位,较深的地层被认定为均匀介质,即横波速度没有变化。
本文介绍了一种全新的仿生学算法——蚱蜢算法,并通过两个测试函数验证其全局搜索能力,证明了其优异的收敛能力以及全局寻优能力。将蚱蜢算法应用于瑞雷波频散曲线反演,通过设置较大的搜索范围以模拟实际工作中缺少足够先验信息的情况,并进行多角度、深层次的探索试验,深入分析了算法的反演能力。最后通过实测数据检测了本算法的适用性。理论数据和实测资料反演的结果表明:
(1) 蚱蜢算法具有较高的全局寻优能力以及抗局部最优解能力,所得解与函数最小解十分接近,说明这是一种有效的全局优化算法;
(2) 蚱蜢算法可应用于瑞雷波频散曲线反演,与传统算法相比,反演过程中目标函数曲线收敛更加稳定,最终迭代误差更接近于零。证明基于蚱蜢算法的瑞雷波频散曲线反演具有较高的精度,是一种简单、稳健的反演算法。
同时,蚱蜢算法也存在一定的缺陷,比如算法的收敛过程相对较长。在今后的研究中,可与其他算法如遗传算法、粒子群算法或人工神经网络算法相结合,探讨如何进一步提高算法的收敛效率。
[1] |
张立, 刘争平. 水平层状介质中基阶瑞利面波椭圆极化特征数值分析与研究[J]. 地球物理学报, 2013, 56(5): 1686-1695. ZHANG Li, LIU Zhengping. A study of the elliptic polarization characteristics of fundamental mode Rayleigh wave based on numerical simulation[J]. Chinese Journal of Geophysics, 2013, 56(5): 1686-1695. |
[2] |
李庆春, 邵广周, 刘金兰, 等. 瑞雷面波勘探的过去、现在和未来[J]. 地球科学与环境学报, 2006, 28(3): 74-77. LI Qingchun, SHAO Guangzhou, LIU Jinlan, et al. Past, present and future of Rayleigh surface wave exploration[J]. Journal of Earth Sciences & Environment, 2006, 28(3): 74-77. DOI:10.3969/j.issn.1672-6561.2006.03.016 |
[3] |
潘佳铁, 李永华, 吴庆举, 等. 中国东北地区地壳上地幔三维S波速度结构[J]. 地球物理学报, 2014, 57(7): 2077-2087. PAN Jiatie, LI Yonghua, WU Qingju, et al. 3-D S-wave velocity structure of crust and upper-mantle beneath the northeast China[J]. Chinese Journal of Geophysics, 2014, 57(7): 2077-2087. |
[4] |
滕吉文, 胡家富, 张中杰, 等. 喜马拉雅碰撞造山带和羌塘盆地的瑞雷波频散与三维速度结构[J]. 石油地球物理勘探, 1998, 33(5): 632-648. TENG Jiwen, HU Jiafu, ZHANG Zhongjie, et al. The frequency dispersion of Rayleigh waves and 3-D velocity structures in Himalayan collision orogenic zone and Qingtang basin[J]. Oil Geophysical Prospecting, 1998, 33(5): 632-648. |
[5] |
黄忠贤, 李红谊, 胥颐. 中国西部及邻区岩石圈S波速度结构面波层析成像[J]. 地球物理学报, 2014, 57(12): 3994-4004. HUANG Zhongxian, LI Hongyi, XU Yi. Lithospheric S-wave velocity structure of west China and neighboring areas from surface wave tomography[J]. Chinese Journal of Geophysics, 2014, 57(12): 3994-4004. DOI:10.6038/cjg20141212 |
[6] |
夏江海, 高玲利, 潘雨迪, 等. 高频面波方法的若干新进展[J]. 地球物理学报, 2015, 58(8): 2591-2605. XIA Jianghai, GAO Lingli, PAN Yudi, et al. New findings in high-frequency surface wave method[J]. Chinese Journal of Geophysics, 2015, 58(8): 2591-2605. |
[7] |
刘江平, 罗银河, 何伟兵. 相邻道瞬态瑞雷面波法与压实度检测[J]. 岩土工程学报, 2009, 31(11): 1652-1659. LIU Jiangping, LUO Yinhe, HE Weibing. Method of neighboring trace transient Rayleigh wave and its application in compactness inspection[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(11): 1652-1659. DOI:10.3321/j.issn:1000-4548.2009.11.002 |
[8] |
Miller R D, Xia J, Park C B, et al. Multichannel analysis of surface wave to map bedrock[J]. The Leading Edge, 1999, 18(12): 1392-1396. DOI:10.1190/1.1438226 |
[9] |
祁生文, 孙进忠, 何华. 瑞雷波勘探的研究现状及展望[J]. 地球物理学进展, 2002, 17(4): 630-635. QI Shengwen, SUN Jinzhong, HE Hua. Review of Rayleigh wave exploration[J]. Progress in Geophy-sics, 2002, 17(4): 630-635. DOI:10.3969/j.issn.1004-2903.2002.04.011 |
[10] |
鲁来玉, 张碧星, 汪承灏. 基于瑞利波高阶模式反演的实验研究[J]. 地球物理学报, 2006, 49(4): 1082-1091. LU Laiyu, ZHANG Bixing, WANG Chenghao. Experi-ment and inversion studies on Rayleigh wave consi-dering higher modes[J]. Chinese Journal of Geophy-sics, 2006, 49(4): 1082-1091. DOI:10.3321/j.issn:0001-5733.2006.04.021 |
[11] |
陈祥, 孙进忠. 改进的等效半空间法及瑞雷波频散曲线反演[J]. 地球物理学报, 2006, 49(2): 569-576. CHEN Xiang, SUN Jinzhong. An improved equivalent homogenous half-space method and reverse fitting analysis of Rayleigh wave dispersion curve[J]. Chinese Journal of Geophysics, 2006, 49(2): 569-576. DOI:10.3321/j.issn:0001-5733.2006.02.033 |
[12] |
李子伟, 刘学伟. 空间假频对频率-波数域瑞利面波频散曲线反演的影响[J]. 石油地球物理勘探, 2013, 48(3): 395-402. LI Ziwei, LIU Xuewei. Effects of spatial aliasing on frequency-wavenumber domain inversion of Rayleigh-wave dispersion curve[J]. Oil Geophysical Prospecting, 2013, 48(3): 395-402. |
[13] |
马见青, 李庆春, 樊金生, 等. 基于S变换的瑞利面波频散分析[J]. 地球科学与环境学报, 2010, 32(3): 319-323. MA Jianqing, LI Qingchun, FAN Jinsheng, et al. Rayleigh wave dispersion analysis based on S transform[J]. Journal of Earth Sciences and Environment, 2010, 32(3): 319-323. DOI:10.3969/j.issn.1672-6561.2010.03.019 |
[14] |
邵广周, 李庆春, 梁志强. 液体表层层状介质导波频散曲线研究[J]. 地球物理学报, 2007, 50(3): 915-920. SHAO Guangzhou, LI Qingchun, LIANG Zhiqiang. A study on dispersion curves of guided wave in layered media with overlying liquid surface[J]. Chinese Journal of Geophysics, 2007, 50(3): 915-920. DOI:10.3321/j.issn:0001-5733.2007.03.033 |
[15] |
熊章强, 张大洲, 秦臻, 等. 瑞雷波数值模拟中的边界条件及模拟实例分析[J]. 中南大学学报(自然科学版), 2008, 39(4): 824-830. XIONG Zhangqiang, ZHANG Dazhou, QIN Zhen, et al. Boundary conditions and case analysis of numerical modeling of Rayleigh wave[J]. Journal of Central South University(Science and Technology), 2008, 39(4): 824-830. |
[16] |
邵广周, 李庆春, 吴华. 基于波场数值模拟的瑞利波频散曲线特征及各模式能量分布[J]. 石油地球物理勘探, 2015, 50(2): 306-315. SHAO Guangzhou, LI Qingchun, WU Hua. Characteristics of Rayleigh wave dispersion curve based on wave field numerical simulation and energy distribution of each mode[J]. Oil Geophysical Prospecting, 2015, 50(2): 306-315. |
[17] |
袁士川, 宋先海, 蔡伟, 等. 瑞雷波有限差分数值模拟中不同自由表面边界条件的对比[J]. 石油地球物理勘探, 2017, 52(6): 1156-1169. YUAN Shichuan, SONG Xianhai, CAI Wei, et al. Comparison of different free-surface boundary conditions for Rayleigh waves finite-difference modeling[J]. Oil Geophysical Prospecting, 2017, 52(6): 1156-1169. |
[18] |
周红, 陈晓非. 凹陷地形对Rayleigh面波传播影响的研究[J]. 地球物理学报, 2007, 50(4): 1182-1189. ZHOU Hong, CHEN Xiaofei. A study on the effect of depressed topography on Rayleigh surface wave[J]. Chinese Journal of Geophysics, 2007, 50(4): 1182-1189. DOI:10.3321/j.issn:0001-5733.2007.04.027 |
[19] |
张碧星, 肖柏勋, 杨文杰, 等. 瑞利波勘探中"之"字型频散曲线的形成机理及反演研究[J]. 地球物理学报, 2000, 43(4): 557-567. ZHANG Bixing, XIAO Boxun, YANG Wenjie, et al. Mechanism of zigzag dispersion curves in Rayleigh exploration and its inversion study[J]. Chinese Journal of Geophysics, 2000, 43(4): 557-567. DOI:10.3321/j.issn:0001-5733.2000.04.017 |
[20] |
Lu L, Wang C, Zhang B. Inversion of multimode Rayleigh waves in the presence of a low-velocity layer:numerical and laboratory study[J]. Geophysical Journal International, 2007, 168(3): 1235-1246. DOI:10.1111/gji.2007.168.issue-3 |
[21] |
Zhang S X, Chan L S, Chen C Y, et al. Apparent phase velocities and fundamental-mode phase velocities of Rayleigh waves[J]. Soil Dynamics and Earthquake Engineering, 2003, 23(7): 563-569. DOI:10.1016/S0267-7261(03)00069-1 |
[22] |
高静怀, 何洋洋, 马逸尘. 黏弹性与弹性介质中Rayleigh面波特性对比研究[J]. 地球物理学报, 2012, 55(1): 207-218. GAO Jinghuai, HE Yangyang, Ma Yichen. Comparison of the Rayleigh wave in elastic and viscoelastic media[J]. Chinese Journal of Geophysics, 2012, 55(1): 207-218. DOI:10.6038/j.issn.0001-5733.2012.01.020 |
[23] |
Yamanaka H, Ishida H. Application of genetic algorithms to an inversion of surface-wave dispersion data[J]. Bulletin of the Seismological Society of America, 2008, 86(2): 436-444. |
[24] |
Beaty K S, Schmitt D R, Sacchi M. Simulated annealing inversion of multimode Rayleigh wave dispersion curves for geological structure[J]. Geophysical Journal International, 2002, 151(2): 622-631. DOI:10.1046/j.1365-246X.2002.01809.x |
[25] |
Shirazi H, Abdallah I, Nazarian S. Developing artificial neural network models to automate spectral analysis of surface wave method in pavements[J]. Journal of Materials in Civil Engineering, 2009, 21(12): 722-729. DOI:10.1061/(ASCE)0899-1561(2009)21:12(722) |
[26] |
Song X, Gu H, Zhang X, et al. Pattern search algo-rithms for nonlinear inversion of high-frequency Rayleigh-wave dispersion curves[J]. Computers & Geosciences, 2008, 34(6): 611-624. |
[27] |
Song X, Tang L, Lv X, et al. Application of particle swarm optimization to interpret Rayleigh wave dispersion curves[J]. Journal of Applied Geophysics, 2012, 84(9): 1-13. |
[28] |
印兴耀, 孔栓栓, 张繁昌, 等. 基于差分进化算法的叠前AVO反演[J]. 石油地球物理勘探, 2013, 48(4): 591-596. YIN Xingyao, KONG Shuanshuan, ZHANG Fanchang, et al. Pre-stack AVO inversion based on the differential evolution algorithm[J]. Oil Geophysical Prospecting, 2013, 48(4): 591-596. |
[29] |
刘晓晶, 印兴耀, 吴国忱, 等. 基于正交匹配追踪算法的叠前地震反演方法[J]. 石油地球物理勘探, 2015, 50(5): 925-935. LIU Xiaojing, YIN Xingyao, WU Guochen, et al. Pre-stack seismic inversion based on orthogonal matching pursuit algorithm[J]. Oil Geophysical Prospecting, 2015, 50(5): 925-935. |
[30] |
孙彩堂, 李玲, 黄维宁, 等. 基于自适应遗传算法的CSAMT一维反演[J]. 石油地球物理勘探, 2017, 52(2): 392-397. SUN Caitang, LI Ling, HUANG Weining, et al. 1D inversion of CSAMT data based on adaptive genetic algorithm[J]. Oil Geophysical Prospecting, 2017, 52(2): 392-397. |
[31] |
张广智, 赵晨, 涂奇催, 等. 基于量子退火Metropolis-Hastings算法的叠前随机反演[J]. 石油地球物理勘探, 2018, 53(1): 153-160. ZHANG Guangzhi, ZHAO Chen, TU Qicui, et al. Prestack stochastic inversion based on the quantum annealing Metropolis-Hastings algorithm[J]. Oil Geophysical Prospecting, 2018, 53(1): 153-160. |
[32] |
Saremi S, Mirjalili S, Lewis A. Grasshopper optimization algorithm:theory and application[J]. Advances in Engineering Software, 2017, 105: 30-47. DOI:10.1016/j.advengsoft.2017.01.004 |
[33] |
凡友华, 刘家琦, 肖柏勋. 计算瑞利波频散曲线的快速矢量传递算法[J]. 湖南大学学报:自然科学版, 2002, 29(5): 25-30. FAN Youhua, LIU Jiaqi, XIAO Boxun. Fast vector-transfer algorithm for computation of Rayleigh wave dispersion curves[J]. Journal of Hunan University (Natural Sciences Edition), 2002, 29(5): 25-30. |
[34] |
Kennedy J, Eberhart R C.Particle swarm optimization[C]. IEEE International Conference on Neural Networks, 1995, 1942-1948.
|
[35] |
Shi Y, Eberhart R C.A modified particle swarm optimizer[C]. IEEE World Congress on Computational Intelligence, 1998, 69-73. https://www.researchgate.net/publication/3755900_A_Modified_Particle_Swarm_Optimizer
|
[36] |
李擎, 张超, 陈鹏, 等. 一种基于粒子群参数优化的改进蚁群算法[J]. 控制与决策, 2013, 28(6): 873-878. LI Qing, ZHANG Chao, CHEN Peng, et al. Improved ant colony optimization algorithm based on particle swarm optimization[J]. Control and Decision, 2013, 28(6): 873-878. |
[37] |
蔡涵鹏, 贺振华, 黄德济. 基于粒子群优化算法波阻抗反演的研究与应用[J]. 石油地球物理勘探, 2008, 43(5): 535-539. CAI Hanpeng, HE Zhenhua, HUANG Deji. Study and application of wave impedance inversion based on particle swarm optimized algorithm[J]. Oil Geophysical Prospecting, 2008, 43(5): 535-539. DOI:10.3321/j.issn:1000-7210.2008.05.009 |
[38] |
谢玮, 王彦春, 刘建军, 等. 基于粒子群优化最小二乘支持向量机的非线性AVO反演[J]. 石油地球物理勘探, 2016, 51(6): 1187-1194. XIE Wei, WANG Yanchun, LIU Jianjun, et al. Nonlinear AVO inversion based on particle swarm optimization least squares support vector machines[J]. Oil Geophysical Prospecting, 2016, 51(6): 1187-1194. |