基于模型分块交叉移动的学习型模拟退火算法地震反演

张赛民，陈灵君，刘海飞，周竹生

(1.中南大学 地球科学与信息物理学院，湖南 长沙，410083；2.中国石化勘探南方分公司，四川 成都，610041)

地球物理模型与数据之间存在非线性关系[1]。获得模型解的过程是通过最小化目标函数来实现，因存在非线性关系，目标函数存在多极值。模型的非唯一性及数据的不完备也会导致产生局部极值，这导致线性化反演的局部优化技术往往会使全局最小值求解失败，甚至可能获得错误解。在没有充分先验信息情况下，通过全局非线性优化技术反演可以获得更优异的模型解。如Monte Carlo优化技术就是一全局优化技术，该技术在先验约束下，通过产生很多个模型找到最符合数据的模型解。Press[2]应用该方法曾估计地球内部弹性常数；Wiggins[3]应用Monte Carlo优化技术在天然地震中反演上地幔的纵波速度，Jin等[4]用两步法的 Monte Carlo方法进行非线性速度反演。Monte Carlo方法最大的问题在于计算量太大，需要反复的正演计算。为了使算法效率提高，已发展出多种随机优化技术，如遗传算法、蚁群算法、粒子群算法及模拟退火算法等。目前地球物理反演中应用比较多遗传算法是采取进化的观点获得最优解，而模拟退火技术是模拟物质粒子结晶的过程来获得模型解。模拟退火自提出以来，多次改进以适应不同的技术问题。Krikpatrick等[5]应用模拟退火研究晶体；Harold等[6]提出快速模拟退火算法；Ingber[7]提出极快速模拟退火算法。目前模拟退火技术已经广泛应用于地球物理反演中。Ma等[8−14]将模拟退火算法应用于地震资料反演中。人们对粒子群算法在地球物理反演领域的研究较少，易远元[15]利用基本粒子群算法实现了叠后地震波阻抗反演，聂茹等[16]用免疫粒子群算法实现了地震波阻抗反演。在此，本文作者针对地震反演这一具体问题，在讨论模拟退火技术与智能单粒子群算法的基础上，提出一种模型交叉移动学习型模拟退火算法，并通过理论模型及实际资料进行波阻抗反演。

1 模拟退火算法改进

模拟退火技术在地球物理反演过程中是应用较广泛的一种非线性反演技术，对初始模型的依赖性较小，不需要计算目标函数对模型的梯度，理论上，足够的多运算可以获得全局极值。模拟退火要获得质量好的解必须增大计算工作量，因而，该方法的应用受到限制。

纪震等[17]在研究现有粒子群算法的基础上，提出了智能单粒子算法(ISPO)。在 ISPO算法的优化过程中，算法不是对整个速度矢量或位置矢量同时更新，而是先把整个矢量分成若干子矢量，并按顺序循环更新每个子矢量。在子矢量的更新过程中，此算法通过引入一种新的学习策略，使得粒子在更新过程中能够分析之前的更新速度，并决定子矢量在下一次迭代的速度。

地震波阻抗反演的模型量很大，尤其是三维地震数据反演，反演涉及的运算成倍增加，加快收敛速度成为波阻抗反演中的关键问题。

1.1 模型扰动公式改进

模拟退火算法在降温过程中，在每种温度下都需要对模型进行扰动，Ingber[7]提出的算法，新的扰动模型依赖于温度的似 Cauchy分布，该算法称为极速模拟退火算法，模型扰动公式如下：

式中：mi为模型扰动值；为初始模型；和分别为模型范围的最大值、最小值；i为模型序号；Tk为第k步温度；u为随机数，在区间(0,1)内均匀分布。

从式(1)和(2)可以看出：模型的扰动只与当前状态的温度、模型范围及随机数有关，没有任何学习功能，因而，在温度一定时，模型是在似柯西分布概率下随机搜索模型值。

标准粒子群算法粒子飞行的速度及位置变换公式如下：

智能单粒子群算法粒子的飞行速度公式为：

其中：a为多样性因子；p为下降因子；s为收缩因子；b为加速度因子(b≥1)，为常数；t为迭代次数；为粒子更新位置；f( )为计算目标函数值；(a/kp)×r为多样性部分；b×为学习部分。同样地学习部分能“记住”使目标函数减少的方向。

比较粒子群算法与模拟退火算法，粒子群算法里面的算法飞行速度项相当于模拟退火算法中的模型扰动项为避免模拟退火算法中模型扰动太随机、无方向感，本文根据粒子群算法思想在模型扰动项里面加入1项学习项：

式中：α为调节因子，可称为加速因子，控制向最优解飞行的速度；s为收缩因子；ΔE为前、后 2次目标函数变化值。

1.2 分块移动反演

模拟退火算法的性能会随着模型量的增加而变差，标准的模拟退火的模型扰动对所有维进行的，在更新模型里各个数据后，计算目标函数。根据目标函数的减小程度，按一定概率接受产生的新模型。计算的目标函数减少程度与模型的全部数据有关，不能判断部分模型变化方向是否向最优的方向移动，这就导致可能部分模型本以按最优点方向移动，而下次随机取值就可能偏离真值。

图1 模拟退火反演分块移动示意图Fig.1 Diagram of simulated annealing inversion with divided block model cross-move

2 模型分块交叉移动的学习型模拟退火算法流程

地震反演的改进模拟退火算法基本流程可以概括如下。

(1)设置初始状态。确定初始温度T0及模型初值与模型空间范围，设定= 0 ，设定j=k=1。模型空间范围可根据测井信息每个模型的变化范围初始模型可以根据测井曲线建立低频模型，也可在模型范围内随机建立。根据初始模型计算出对应的目标函数值E()。

其中：

j为子块序号；k为循环次数；η为依赖于温度的似Cauchy分布函数；Ti为第i步迭代温度；u为随机数，在区间范围(0,1)内均匀分布。

根据式(10)计算对应的目标函数E，得到计算。

(3)按概率接受解。若ΔE≤0，则接收新模型；若ΔE≤0，则新模型根据概率

与一均匀分布在(0,1)的随机数内进行比较，若概率大于随机数，则接收新模型。式(13)为采取Ingber[4]提出的极速模拟退火算法提到广义Gibbs分布的概率，其中h为实数。

(4)在温度Ti下，若k小于N(其中N为Mapkob链的长度)，在此重复模型扰动与接收过程，即重复步骤(2)和(3)。

(5)逐步降温。根据Ingber[4]提出的极速模拟退火算法，降温方式为：

式中：T0为给定的初始温度；i为迭代次数；c为给定常数；n为待反演的模型参数个数。在实际应用中，将式(14)通常改写成：

选择0.7≤β＜1.0，并采用0.5或1.0代替上式的1/n。

(6)重复步骤(2)～(5)，计算出第j个子块，直至满足收敛条件为止。

(7)j=j+1，设定= 0 ，设定k=1，重复步骤(2)～(6)，获得全部模型最优解。

上述模拟退火算法在2个方面对传统模拟退火算法进行了改进：首先，在模型扰动中加入了根据上次寻优的过程方向的学习来加快收敛速度；另一方面，针对模型非常多，模拟退火算法性能变差的情况，对地震反演采取分块交叉移动的方式反演地震资料所表征的岩性参数。该算法可称为模型分块交叉移动的学习型模拟退火算法。

本文提出的针对地震反演的模型分块交叉移动的学习型模拟退火算法同样可以实现其他多模型复杂的非线性问题的求解。

3 理论模型试算和实际资料地震波阻抗反演

这里应用模型分块交叉移动的学习型模拟退火算法来反演地震波阻抗。首先构造1个理论模型反演波阻抗，然后，将提出的算法用于到某工区的实际地震资料处理。

在反演过程中，选取初始温度T0=10.0，终止温度为Tend=10.0，温度衰减系数α= 0.98，接收概率中实数h=5，加速因子α=0.8，收缩因子s=2。

3.1 理论模型反演

本文设计了1个7个地质层和8个地震道的二维理论模型，时窗长为500 ms，采样间隔为2 ms。表1所示为各个层的波阻抗。为了解薄层反演效果，设计的模型对应的5层厚度逐渐变小。图2(a)所示为理论波阻抗模型。模型合成记录使用的地震子波主频为35 Hz，最小相位的理论雷克子波，长度为 100 ms。图2(b)所示为用雷克子波与反射系数褶积后的合成记录。

表1 地质理论模型波阻抗值Table 1 Impedance value of theoretic model

图3(a)所示为新模拟退火算法使用的初始波阻抗模型。初始模型设计为3层，模拟退火模型扰动的上下边界为此初始模型值−25%～25%范围内。为避免无用计算，在程序中设计了：若扰动超过此范围，则重新产生在此模型范围的值。

对初始波阻抗模型合成记录与理论模型合成记录进行相关运算，其相关系数仅为0.51。根据张赛民[18]的反演理论，反演目标函数为：

式中：L=lnZ，Z为波阻抗列向量；W为子波矩阵；D为系数矩阵；S为理论模型合成地震记录列向量；φ为势函数；在反演过程中取δ=0.2，μ=0.02。式(16)具有边界保持块约束的特点。图3(b)所示为新模拟退火算法波阻抗反演结果，与图1(a)所示理论模型比较可以看出：反演结果与理论模型基本一致，算法达到了预定目标，精度得到满足。根据该反演结果合成地震记录与真模型合成的记录进行相关运算，其相关系数达到0.992，改进的模拟退火算法收敛性很好。

图2 波阻抗理论模型与地震合成记录Fig.2 Theoretic model of impedance and seismic synthetic seismogram

图3 初始波阻抗与反演波阻抗Fig.3 Initial impedance and invert impedance

本文提出的模型分块交叉移动学习型模拟退火算法与常用的快速模拟退火算法两者收敛速度效果的对比结果如图4所示。从图4可以看出：本文提出的模拟退火算法较常用的模拟退火算法在相同条件下收敛速度更快，提高了地震反演效率。

3.2 实际资料反演

对某工区抽取1个地震剖面，应用本文提出的新模拟退火算法进行波阻抗反演。算法目标函数为边界保持的块约束目标函数，算法选择的参数与前面的相同。图5所示为原始地震剖面，图6所示为反演波阻抗震剖面，时窗为780～1 175 ms。

在反演过程中，利用测井资料反演子波，同时利用测井资料建立初始模型，并设定模拟退火算法模型扰动范围在初始模型−25%～25%范围内，反演后的合成记录与原始地震记录相关系数达0.977。从图6所示反演剖面可以看出：在测井位置反演的波阻抗与测井资料计算的波阻抗形态基本一致，从波阻抗整个剖面可以看出纵向岩层界面清晰，横向岩性变化延续性也显示清楚。该反演结果提供了详细的波阻抗信息，可为储层预测提供可靠的基础资料。

图4 模拟退火算法效果比较图Fig.4 Comparison of simulated annealing

图5 原始地震剖面Fig.5 Section of seismic data

图6 反演波阻抗剖面Fig.6 Section of invert impedance

4 结论

(1)对模型扰动公式进行了改进，传统的模型的扰动只与当前状态的温度、模型范围及随机数有关，没有任何学习功能。为避免模拟退火算法中模型扰动太随机，无方向感，根据粒子群思想在模型扰动项里面加入1个学习项，该学习项具有使扰动向优化方向移动的功能。

(2)针对模拟退火算法的性能随着模型量的增大而变差的情况，提出分块交叉移动的方法来加快模拟退火收敛速度，提高了算法性能。

(3)理论模型及实际资料反演结果表明了该方法的有效性、可行性和可靠性。该方法还可用于其他多维多极值的非线性反演。

[1]王家映.地球物理反演理论[M].北京: 高等教育出版社,2002: 1−2.WANG Jia-yin.Geophysical inversion theory[M].Beijing:Higher Education Press,2002: 1−2.

[2]Press F.Earth models obtained by Monte-Carlo inversion[J].Journal of Geophysical Research,1968,73(1): 5223−5234.

[3]Wiggins R A.Minimum entropy deconvolution[J].Geoexploration,1978,16(1): 21−35.

[4]Jin S,Madiriaga R.Nonlinear velocity inversion by a two-step Monte Carlo method[J].Geophysics,1994,59(4): 577−590.

[5]Kirkpatrick S,Gelatt C P,Vecchi M P.Optimization by simulated annealing[J].Science,1983,220(4598): 671−680.

[6]Harold S,Ralph H.Fast simulated annealing[J].Physics Letters,1987,122(3/4): 157−162.

[7]Ingber L.Very fast simulated re-annealing[J].Mathl Comput Modelling,1989,12(8): 967−973.

[8]MA Xin-quan.Global joint inversion for the estimation of acoustic and shear impedances from AVO derived P- and S-wave reflectivity data[J].First Break,2001,19(10): 557−566.

[9]MA Xin-quan.Simultaneous inversion of prestack seismic data for rock properties using simulated annealing[J].Geophysics,2002,67(6): 1877−1885.

[10]Misra S,Sacchi M D.Nonlinear one dimensional prestack seismic inversion with edge preserving smoothing filter[C]//Extended Abstract,EAGE 69th Conference and Exhibition.London,United Kingdom,2007: 26−30.

[11]Misra S,Sacchi M D.Global optimization with model space preconditioning: Application to AVO inversion[J].Geophysics,2008,73(5): 71−82.

[12]Varela O J,Tomes-Verdin C,Sen M K.Enforcing smoothness and assessing uncertainty in non-linear one dimensional prestack seismic inversion[J].Geophysical Prospecting,2006,54(3):239−259.

[13]杨午阳.模拟退火叠前AVA同步反演方法[J].地球物理学进展 ,2010,25(1): 219−224.YANG Wu-yang.Pre-stack simultaneous AVA inversion algorithm with simulated annealing[J].Progress in Geophysics,2010,25(1): 219−224.

[14]刘全新,方光建.利用模拟退火算法实现地震叠前反演[J].岩性油气藏,2010,22(1): 87−92.LIU Quan-xin,FANG Guang-jian.Using simulated annealing algorithm to pre-stack seismic inversion[J].Lithologic Reservoirs,2010,22(1): 87−92.

[15]易远元.地震波阻抗反演的粒子群算法实现[J].石油天然气学报,2007,29(3): 79−81.YI Yuan-yuan.Implementation of particle swarm algorithm in seismic wave impedance inversion[J].Journal of Oil and Gas Technology,2007,29(3): 79−81.

[16]聂茹,岳建华,邓帅奇.地震波阻抗反演的免疫粒子群算法[J].中国矿业大学学报,2010,39(5): 733−739.NIE Ru,YUE Jian-hua,DENG Shuai-qi.Immune particle swarm optimization for seismic wave impedance inversion[J].Journal of China University of Mining & Technology,2010,39(5):733−739.

[17]纪震,周家锐,廖惠连.智能单粒子优化算法[J].计算机学报,2010,33(3)3: 556−561.JI zhen,ZHOU Jia-rui,LIAO Hui-lian.A novel intelligent single particle optimizer[J].Chinese Journal of Computers,2010,33(3):556−561.

[18]张赛民.基于边界保持块约束叠后波阻抗及叠前AVA三参数同步反演研究[D].长沙: 中南大学地球科学与信息物理学院,2011: 70−73.ZHANG Sai-min.Research of post-stack impedance inversion and AVA three-term simultaneous inversion based on edge preserving blocky constrain[D].Changsha: Central South University.School of Geosciences and Info-Physics,2011:70−73.