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

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

(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；

2. 中国石化勘探南方分公司，四川 成都，610041)

摘要：针对地震非线性反演问题，提出一种基于模型分块交叉移动的学习型模拟退火的全局优化地震反演方法。其步骤为：首先，在模拟退火算法及粒子群算法基础上，在算法模型扰动项里面加入1个向目标优化的方向移动的学习项；其次，针对地震反演模型数量多及地震记录为褶积形式的特点，采用模型分块交叉移动的方法来实施模拟退火反演，给出模型分块交叉移动的学习型模拟退火算法流程。研究结果表明：该方法具有收敛速度快、精度高、实现简单、高效的特点，可以用于其他多维多极值的目标函数反演。

关键词：全局优化；模拟退火；地震反演；模型分块交叉移动

中图分类号：P631.4           文献标志码：A         文章编号：1672-7207(2012)03-1040-07

Seismic inversion based on simulated annealing of learning and cross-moving with divided block model

ZHANG Sai-min¹, CHEN Lin-jun², LIU Hai-fei¹, ZHOU Zhu-sheng¹

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. SINOPEC Southern Exploration Company, Chengdu 610041, China)

Abstract: To solve the seismic nonlinear inverse problem, a seismic inversion method was presented based on simulated annealing of learning and cross-moving with divided block model. The algorithm was proposed in two aspects. Firstly, based on the algorithm of simulated annealing and particle swarm optimization, simulated annealing formula of perturbed model was added with a learning item which moved to the direction of objective optimization. In addition, because the seismic inversion model’s number are large and form of seismic records is convolution, simulated annealing algorithm was implemented by the way of divided block model cross-moving. The results show that the proposed method has fast convergence and high accuracy. The new simulated annealing is simple and efficient, and it can be used to other multidimensional object function inversion.

Key words: global optimization; simulated annealing; seismic inversion; divided block model cross-move

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

1  模拟退火算法改进

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

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

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

1.1  模型扰动公式改进

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

 [U4]           (1)

      (2)

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

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

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

    (3)

              (4) 

式中：为飞行速度(粒子移动的距离)；r₁和r₂为均匀分布在(0,1)区间随机数；c₁和c₂为学习因子，也称加速因子；为个体粒子最优位置；为整个粒子群全局最优位置；t为迭代次数；为粒子更新位置，相当于模拟退火算法的新的模型扰动；i为粒子序号。式(3)后两项因利用了粒子群的“历史”最优位置，粒子将沿着使目标函数减少的方向位置移动，能“记住”使目标函数减少的方向。[U7] 

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

            (5)

          (6)

        (7)

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

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

   (8)

                (9)

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

从式(8)可以看出：[U11] 依然保留在似柯西概率分布下寻找随机特性的解，加入的能使模型扰动向具有使目标函数减少的优化方向寻找解。该模型扰动公式修改后具有向“记忆”较优的方向靠拢的功能，从而加快了模拟退火收敛，减少反演运算时间。

1.2  分块移动反演

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

如记为地震数据、为子波数据、为反射系数数据。由形成地震数据的褶积模型公式可知，地震数据与子波数据和[U12] 1个子波长度l的反射系数有关，也就是说只有这个子波长度的模型(反射系数)改变才会对该点的地震数据起作用。进一步，对于旁瓣趋近于0的零相位子波[U13] ，旁瓣位置的模型改变对地震数据起作用小，作用最大的是对应子波主瓣位置[U14] ，主瓣位置为地震数据对应的模型(反射系数)位置。反演过程中，计算地震反演目标函数时，因上部分时窗的地震数据只与上部分的岩性有关，下一部分地震数据只与下部分的岩性有关，这就使分块反演提供可能。其基本思想是：固定一部分初始模型，时间域从低到[U15] 高依次反演。因每个地震数据取决于子波和[U16] 1个子波长度的模型信息，每个子块可按1个子波长度划分，每次只更新该子块的模型值，其他子块保持不变。地震记录最后的形成取决于所有地震模型量，若模型模型扰动按每个子块依次更新，则目标函数最优化过程可能陷入局部极值。为解决这一问题，在地震反演可设置按1个子波长度划分子块，然后，按时间每半个子块从低到高移动依次更新每个子块，这样，[U17] 每个子块的模型更新都与前后子块有重复的半个子块。更新模型过程可以用图1表示，图1[U18] 中l为子波长度，目标寻优过程实现了模型分块交叉移动模拟退火算法来进行地震反演，其实现方式是每次模型的扰动在子波长度范围内，并以半个子波长度滑动，最后实现模型空间所有维数据的反演。

图1  模拟退火反演分块移动示意图

Fig.1  Diagram of simulated annealing inversion with divided block model cross-move

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

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

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

(2) 新模型扰动。对当前模型采用带学习功能扰动模型，有

       (10)

其中：

      (11)

           (12)

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

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

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

         (13)

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

若按概率接收了模型，则设置，。

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

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

           (14)

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

              (15)

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

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

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

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

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

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

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

在反演过程中，选取初始温度T₀=10.0，终止温度为T_end=10.0，温度衰减系数，接收概率中实数h=5，加速因子，收缩因子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]的反演理论，反演目标函数为：

  (16)

式中：L=lnZ，Z为波阻抗列向量；W为子波矩阵；D为系数矩阵；S为理论模型合成地震记录列向量；为势函数；在反演过程中取，。式(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.

(编辑 陈灿华)

收稿日期：2011-05-10；修回日期：2011-07-28

基金项目：国家科技支撑计划项目(2011BAB04B08)；国家自然科学基金资助项目(41074085，40804027)

通信作者：张赛民(1977-)，男，湖南桃江人，博士，讲师，从事地震数据处理、资料解释及反演研究；电话：13107417858；E-mail: saiminzhang@126.com