基于非常快速模拟退火算法的地震谱反演方法①

张广智，李 宁，李 超，陈怀震，陈 雷

（中国石油大学地球科学与技术学院，山东 青岛 266555）

0 引言

随着石油勘探开发的不断深入，地震勘探已经从简单的构造识别转向复杂地表、复杂构造、裂缝、薄储层和老油田剩余油的勘探。勘探难度的加大也对能够得到具有更高信噪比和分辨率的地震资料提出了挑战。地震子波大约在30Hz左右，对于25m以下的薄层很难区分它的顶部和底部［1］。但是工业的目标地质体的厚度一般在10m甚至更小。反演薄层在储层识别和油藏描述中起到了重要的作用［2－3］。

Charles I．Puryear和 John P．Castagna［4］在2008年提出了谱反演的方法。谱反演是一种基于谱分解技术，在时频域建立目标函数的一种地震反演技术，它能够反演出小于调谐厚度的薄层。在过去的几十年里人们一直利用 Widess模型［5］理论得到的调谐厚度分析方法指导薄层的识别。Partyka等［6］以及 Marfurt和 Krilin［7］采用基于离散傅里叶变换得到的频谱分解结果来求取薄层厚度的方法。但是如果地震频带宽度不能够清晰识别谱峰和陷波变化规律时，谱分解对于分辨薄层还是存在困难。这 也 推 动 了 Partyka、Portniaguine、Puryear 和Castagna等人对谱反演的探索和发现。

由于地球物理问题中普遍存在的非线性和不适定性。本文采用全局寻优的非常快速的模拟退火算法对上述的时频域目标函数求最小值。该方法不依赖初始模型，能够有效的跳出局部极值，并且对先验信息获知较少的情况下也能得到较好的结果。

1 谱反演原理

谱反演研究思路如图1所示。

图1 谱反演流程图Fig．1 The flowchart of spectral inversion．

谱反演方法是一种不同于传统地震反演的新方法，可以用于薄层成像。谱反演是利用时频分析和谱分解技术获取地震资料的局部频谱信息，首先在时域将反射系数分为奇偶分量，然后在频域建立目标函数，再利用改进的模拟退火算法求解。其中关键技术是目标函数的建立和求解。在理想的情况下可以识别任意薄层，并且精确的反演出反射系数的值。该方法不需要任何先验模型、任何反射系数假设、任何层位约束、也不需要井强制约束等优点。下面重点介绍目标函数的建立和模拟退火的寻优算法。

1．1 目标函数的建立

目标函数是由频率域的地震数据、地震子波及时域的反射系数的奇、偶分量组成。如果反演过程中反射系数正确，目标函数能够达到最小值，作为判断反演结果正确与否的依据。

基于图2的多层模型建立目标函数，经过复杂推导，最后整理得目标函数表达式为［4］

式中S（t，f）、W（t，f）分别为短时傅里叶变换地震记录的振幅谱和子波的振幅谱；re（t），ro（t）为图2模型中反射系数对的偶分量和奇分量；ae，ao为权衡系数，依据信噪比的高低来确定二者的比值；fl，fh分别为反演频带的低截频和高截频；tw为反演的半个时窗长度。Ti为反射系数对之间的时间厚度。红色三角表示分析点的位置。

图2 多层模型Fig．2 Multi－layers model．

构建目标函数理论是基于任何一个反射系数序列都可以分解成奇偶分量的和。奇分量具有相等的大小和不同的符号，偶分量具有相等的大小和相同的符号；奇分量不利于检测薄层，而偶分量可以提高分辨薄层的能力。谱反演的实质就是利用偶分量在厚度趋于零时的有效干涉来提高地震资料的分辨率，同时偶分量也有利于抗噪。

目标函数是通过结合小时窗内地震资料时频分析的结果和反射系数的奇偶分量，类似于频率域褶积原理建立起来的。谱反演的目标函数与时间域的褶积残差目标函数相比具有更好的收敛性和约束能力，能够有效的减少地震反演中的多解性问题。在反演的过程中因为地震资料是带限的，所以不能在全频带进行反演，只能选取信噪比较高的优势频带进行反演。因此不能得到类似于脉冲信号的反射系数，只能得到一个提高分辨率的剖面，但它代表一定的反射系数的信息。

1．2 模拟退火算法

模拟退火法（Simulated Annealing Method，简称SA），最早是由 Metropolis等人［8］在1953年提出，1983年 Kirkpatrick等［9］成功将它应用于组合优化问题中。Rothmann最早将模拟退火算法引入到地球物理问题中，解决地震资料的剩余静校正问题。随后人们开始利用其进行一维地震波形反演，一维电测深反演以及叠前速度分析［10］。由于模拟退火的计算过程简单、通用，鲁棒性强，适用于并行运算，不过多的依赖初始模型，能够有效的跳出局部极小，所以它成为人们解决非线性优化问题常用的算法［11］。

Kirkpatrick提出的算法是以Gaussian概率分布的准平衡Boltzmann－Gibbs统计为基础的，为经典的模拟退火算法（MSA）；Szu和 Hartley用Cauchy－Lorentz分布代替Gaussian分布，形成快速模拟退火算法（FSA）；Stariolo和Tsallis泛化FSA方法，提出广义模拟退火算法（GSA），包含了CSA和FSA，计算速度比FSA 快；Ingber［12］则采用了依赖于温度的似Cauchy分布产生新的扰动模型，提出非常快速的模拟退火算法（VFSA）；张霖斌等［13］以广义Boltzmann－Gibbs统计理论为基础，采用依赖于温度的似Cauchy分布产生新的扰动模型，建立一种新的快速模拟退火算法。

陈华根等［14］深入分析VFSA的核心技术，发现VFSA存在的缺陷，将VFSA的退火计划与模型扰动相匹配，提出了改进的模拟退火算法（MVFSA）。思路是：算法开始在降温计划的控制下，模型作全局扰动以搜索并锁定最优解区间；在锁定最优解空间后，在新的降温计划下，模型在已被接受模型周围做局部扰动，以近似局部搜索的模拟退火方式，逐步逼近最优解；新的退火计划作适当的回火升温，给予已寻找到的最优解空间再一次跳出局部极小值的机会，使得找到的最优解更可靠。同等条件下，MVFSA所花费的时间要比VFSA少20%～30%，算法显得更稳健。

模拟退火算法之所以能够称为全局寻优算法，是因为它在搜索过程中以一定的概率接受“恶化解”，温度较高时，接受的概率近似为1，随着温度的降低，接受“恶化解”的概率逐渐降低，最后变为局部寻优。

图3 全局搜索与局部搜索Fig．3 Global search and loacal search．

随着模拟退火算法应用的广泛，人们对其进行了多种改进，现在应用较多为非常快速的模拟退火算法。本文采用了VFSA算法并将其和经典的模拟退火算法在机理上进行了比较研究［13－14］。

经典的模拟退火算法多采用均匀的模型扰动方式即：mj＝mi＋step×（rand－0．5）。这里mj表示新解；mi表示前一个解；step表示步长（通常与温度无关）；rand表示0～1之间的随机数。而VFSA的模型扰动方式采用依赖于温度的似Cauchy分布法［12］：

由Boltzmann－Gibbs［14］分布，可得出新的接受概率计算公式为

本文中h＝－5；ΔE扰动前后的能量差；T为当前的温度。T的计算公式为

式中T0为初始温度，本文中设为T0＝100；α＝0．96；c＝1；N＝2。

图4 模型扰动量和接受概率随迭代次数的变化Fig．4 Variation of the model turbulent and acceptance probability with iteration．

经典的模拟退火算法与VFSA的模型扰动量、接受概率随着迭代次数的变化关系如图4所示。可以看出VFSA在高温时搜索范围比MSA大，在低温时搜索范围比MSA小，这样更有利于跳出局部极值，加快了收敛的速度。VFSA接受概率随迭代次数的变化较MSA要快，并且曲线的尾端较MSA略有上翘，这是由于在低温时VFSA的扰动量变小，使得反演解更靠近全局最小解，反演精度更高。

利用好模拟退火这个全局寻优的优化算法需要模型扰动、接受概率、降温函数3者之间的有效配合。

2 模型试算

2．1 单道多层模型试算

主频为30Hz，4ms采样的雷克子波与反射系数模型（图5（b）1）褶积合成地震记录，加上5%的高斯随机噪声，得到观测值（图5（a）1）。利用非常快速模拟退火算法对依式（1）建立的目标函数进行求解。由于噪声的存在，且从噪声分布的频带上看，噪音使得高频带的地震记录信噪比变低，使用5～90 Hz频带信号参与反演。

图5 观测与反演的地震记录和反射系数比较Fig．5 The seismic waves and reflectivities from record and the inversion．

反演地震记录（图5（a）2）与观测地震记录的相对误差为0．246%，反演反射系数（图5（b）2）与真实反射系数的相对误差为3．814%。通过图5可以看出利用谱反演方法能够准确的反演出反射系数的大小和位置。不能够在地震记录中分辨的薄层，能够在反射系数中得到分辨，提高了分辨率。

2．2 楔形模型试算

主频为30Hz，4ms采样的雷克子波与反射系数模型褶积，加上5%的高斯随机噪声，得到图6（a）的楔形层状地震记录，作为反演输入。

图中黑色实竖线的位置即利用λ／4确定的调谐厚度的位置。从反演结果（图6（b））与输入地震记录的比较，可以看出谱反演可以识别出小于调谐厚度的薄层。在250～350ms（矩形框）的位置存在3个反射面，在地震剖面上不能清楚的分辨，但是从反演的结果能够清楚的识别。理论上谱反演可以反演出任意厚度的薄层，实际由于噪声的存在只能得到类似反射系数的剖面。

图6 模型数据的地震记录与谱反演的结果Fig．6 The model of seismic record and spectral inversion result．

针对上述楔形层状模型的地震记录、反演结果、真实反射系数做出频谱图（图7）。曲线1为真实地震记录频谱，可以看出原始地震记录带限，缺失高频信息。曲线3为真实反射系数的频谱频带最宽，是最理想的反演结果。由于噪声的存在和地震记录的带限性，谱反演没有恢复出全频带的反射系数，但是对于地震记录的拓频作用明显，极大的提升了60～90Hz频带范围内的有效信号。

图7 地震记录的频谱、反演反射系数的频谱与真实反射系数的频谱的比较Fig．7 Comparison of spectrums of seismic record、inversion reflectivity and true reflectivity．

2．3 marmous模型

下面截取部分marmous模型（纵波速度）求取反射系数与4ms采样的雷克子波褶积合成地震记录。如（图8（a））所示，在150道至350道间有一透镜状含气砂体，但是由于子波的影响，其边界不是很明显。在已知地震记录和准确子波的情况下，对地震道进行谱反演试算。

通过谱反演试算，得到结果（图8（b））。可以看出谱反演能够有效的去除子波的影响，清楚地识别出砂体及其他反射体的分界面；能够准确的反演出反射系数的位置和大小；极大的提高了数据体的分辨率，有利于精细的地震解释。

3 实际数据处理

该工区为加拿大Alberta白垩系浅层的部分地震数据，2ms采样，在第60～80道0．02～0．03处存在一个气层。应用测井方法从实际数据中提取较为准确的子波，然后对其进行谱反演（图9）。

图8 原始叠后地震剖面和反演结果Fig．8 The original seismic record and spectral inversion result．

图9 加拿大某工区原始叠后地震剖面和谱反演结果Fig．9 The original seismic record and spectral inversion result at one site in Canada．

通过反演前后数据的比较，可以看出数据的分辨率得到了显著的提高，同相轴较反演前变细，边界更加明显。地震解释结果比原始数据更加可靠。

4 结论

通过对谱反演方法的推导以及模型数据和实际资料的处理，可以得出以下几点结论：（1）全局搜索的模拟退火算法可以有效的跳出局部极值，收敛到全局极小，但是收敛速度慢。本文采用的改进的模拟退火算法VFSA提高了搜索过程中的收敛速度，并且保证了解的正确性和精度。（2）谱反演作为一种精细识别薄层、反演反射系数的方法，可以准确的反演出反射系数的大小和位置，其目标函数与时间域的目标函数相比具有更好的收敛性。（3）虽然使用部分频带信号参与反演，但是从反演的效果看反演有效地去除了子波的影响，明显的拓宽了原始资料的频带，极大地提高了资料的分辨率，对于精细的地震解释起到了作用。

在反演过程中，要反复的对地震数据和子波进行傅里叶变换，因此谱反演相对时域的其它反演方法更费时；因为要变换到频率域处理，所以该方法对于噪音较敏感。谱反演方法在国内虽然还不十分成熟，但是对于薄层的识别，地震资料分辨率的提高具有深刻意义。

［1］ Chopra S，J P Castagna，O Portniaguine．Seismic resolution and thin－bed reflectivity inversion［J］．Canadian Society of Exploration Geophysicists Recorder，2006，31：19－25．

［2］ 秦德文．基于谱反演的薄层预测与反演方法研究［D］．青岛：中国石油大学（华东），2009．

［3］ 陈科．基于模拟退火的谱反演方法研究［D］．青岛：中国石油大学（华东），2010．

［4］ Charles，Puryear，Castagna John P．Layer－thickness determination and stratigraphic interpretation using spectral inversion：Theory and application［J］．Geophysics，2008，73（2）：37－48．

［5］ Widess M．How thin is a thin bed？［J］．Geophysics，1973，38：1176－1180．

［6］ Partyka G A，Gridley J A，Lopez J A．Interpretational aspects of spectral decomposition in reservoir characterization［J］．The Leading Edge，1999，18（3）：353－360．

［7］ Marfurt K J，Kirlin R L．Narrow－band spectral analysis and thin－bed tuning［J］．Geophysics，2001，66：1274－1283．

［8］ Metropolis N，Rosenbluth A，Rosenbluth M，et al．Equation of state calculations by fast computing machines［J］．J．Chem．Phys，1953，21：1087－1092．

［9］ Rothman D H．Nonlinear inversion statistical mechanics and residual statics estimation［J］．Geophysi－cs，1985，50（12）：2784－2796．

［10］ 师学明，王家映．模拟退火法［J］．工程地球物理学报，2007，4（3）：165－174．

［11］ 康立山，谢云，尤矢勇，等．非数值并行算法－模拟退火算法［M］．科学出版社，2003．

［12］ Inger L．Very fast simulated annealing［J］．Math Conput Modeling，1989，12（8）：967－973．

［13］ 张霖斌，姚振兴，纪晨，等．快速模拟退火算法及应用［J］．石油地球物理勘探，1997，32（5）：654－661．

［14］ 陈华根，吴健生，王家林，等．模拟退火算法机理研究［J］．同济大学学报（自然科学版），2004，32（6）：802－805．