CN103472482A - 基于基因排序体系的多波地震资料时间域匹配方法及系统 - Google Patents

Abstract

本发明提供一种基于基因排序体系的多波地震资料时间域匹配方法及系统，该方法包括：采集纵波地震道数据；将纵波地震道数据进行反演，得到纵波时间域的横波地震道数据；采集转换波地震道数据；根据横波地震道数据以及所述的转换波地震道数据建立相似性矩阵；根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数；根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。解决多波地震勘探中转换波地震数据与纵波地震数据的时间域匹配问题，为纵波与转换波资料联合处理、解释和反演提供准确的资料，提高综合应用纵横波地震数据属性进行储层流体检测的准确度。

Description

基于基因排序体系的多波地震资料时间域匹配方法及系统

技术领域

本发明是关于地震勘探技术领域，特别是关于石油地球物理勘探领域的转换波地震数据与纵波地震数据时间域的匹配技术，具体的讲是一种基于基因排序体系的多波地震资料时间域匹配方法及系统。

背景技术

随着多波多分量数字地震勘探仪器与技术的发展和采集技术的进步，多分量采集成本不断下降，越来越多的油田进行了多波多分量地震勘探的试验研究。BP、Chevron等大石油公司以及CGG、西方、Veritas等大的地球物理服务公司,在海上和陆上多波多分量勘探方面做了许多工作，在岩性和流体识别、构造成像改善、油藏监测、裂缝与各向异性分析等方面的试验和研究取得了诸多进展。近年来，多波多分量资料的采集技术已经取得了具有突破性的进展。挪威国家石油公司（Statoil）的Berg等人(1994)研究开发了海底电缆多波接收系统，大大改善了多波资料的数据质量，引起了多波勘探技术的又一波研发热潮。本世纪初，I/O和CGG等地球物理服务公司研发了微电子机械系统（micro-electronic–mechanical-system简称MEMS）数字式检波器，代替了传统的数模组合检波器，提高了陆上多波资料的信噪比和分辨率，也极大地降低了陆上多波勘探的成本，带动了陆上多波勘探技术的发展。与此同时，在多波资料处理方面，Tessmer and Behle(1988)和Harrison(1992)等人完成了多波资料常规处理流程的研发工作，其中包括转换波抽道集、动校正、速度分析、倾角时差校正、叠加和叠后偏移。而对于多波资料的后续应用方面，在多波资料用于储层预测和油气预测的综合应用中，虽然投入了同样的人力、物力，效果显著的系统性的综合应用实例并不很多。

目前多波勘探技术与几年前相比有了长足的发展，多波勘探技术已逐步成为地震勘探的有效技术，相应关于衰减各向异性和多波油气预测技术的研究方面也已经受到国内外的高度重视。无论是海上还是陆上，多波勘探已经成为有效的油气勘探开发的手段。海上多波勘探效果显著，出现了许多有标志性和权威性的应用实例。而在陆上，采集成本已有很大降低，由于缺乏很好的技术手段和实用软件，处理和解释成本并未下降，且大部分地区并未得到让勘探家满意的处理、解释成果。同时，陆上多波应用范围基本以岩性和裂缝预测为主，基于全波属性的地质综合应用由于资料处理、联合反演、多波资料综合解释研究等方面还存在许多技术难题。近年实施的多波勘探和试验的项目众多，至今经济效益显著的，但有标志性、权威性和系统性的应用实例仍然不多见。

多波勘探进入21世纪后取得了长足的进步。2003年，中石油在长庆苏里格气田采集了3D多波资料。2004年，中石化分别在胜利和四川都采集了3D3C资料(Tang and Huang2006)。这几次采集都采用了先进MEMS检波器技术，资料品质有了显著提高。为资料处理和解释打下了良好基础。2006-2008年中石化西南油气分公司在新场气田对转换波3D3C勘探技术成果进行了全面的生产应用，取得了良好的勘探效果，深层致密裂缝型气藏的勘探成功率由过去的不到50%提高到了80%以上，并且高产工业气井的比例大大提高，成为近年来多波多分量勘探的标志性项目，在陆上转换波3D3C勘探方面具有明显的示范作用。

多波多分量勘探技术研究，采集技术设计技术已经掌握并能够针对不同的地质需求和地震地质条件设计出满足要求的观测系统。资料采集的生产组织和质量控制已十分成熟；在多波多分量处理方面，尚未形成具有工业化生产水平的多波多分量处理软件系统，国内科研机构开发的软件很不成系统，涉及到多波油气预测的转换波叠前偏移处理、纵横波自动匹配、叠前联合反演、横波分裂分析等核心关键技术的研发力量分散，成效不大。如何有效地利用多波地震资料进行油气预测，以及建立起纵横波地震资料的精确匹配，是当前有效地利用多波地震资料进行储层预测及油气预测的关键所在，是多波勘探资料处理和纵横波联合反演的瓶颈。

现有的纵波与转换波的匹配方法主要包括两类，第一类主要是基于层位控制的，分别在纵波与转换波资料剖面上人工解释相同的地震层位，然后进行二者之间的同相轴时间匹配，该方法需要解释准确的目的层，受人为影响因素较大。另一类是直接利用纵波与转换波的波形相似性求取相关函数谱，然后进行人工拾取速度比进行匹配，该方法的前提假设条件是纵波与转换波具有相同的波形特征，但实际纵波与转换波资料的波形特征是不同的，包括振幅、相位和极性等。

因此，由于地下介质的复杂性及地震子波的带限特性，使得纵波与转换波的时间匹配目标函数成强非线性特性，而且，纵波与转换波资料的波形特征存在很大的差异，需要使用稳定、有效的非线性全局最优求解方法。

发明内容

为了克服现有技术存在的上述问题，本发明提供了一种基于基因排序体系的多波地震资料时间域匹配方法及系统，解决多波地震勘探中转换波地震数据与纵波地震数据的时间域匹配问题，为纵波与转换波资料联合处理、解释和反演提供准确的资料，提高综合应用纵横波地震数据属性进行储层流体检测的准确度。

本发明的目的之一是，提供一种基于基因排序体系的多波地震资料时间域匹配方法，包括：采集纵波地震道数据；将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据；采集转换波地震道数据；根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵；根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数；根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

本发明的目的之一是，提供一种基于基因排序体系的多波地震资料时间域匹配系统，所述的系统包括：纵波地震道数据采集装置，用于采集纵波地震道数据；反演装置，用于将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据；转换波地震道数据采集装置，用于采集转换波地震道数据；相似性矩阵建立装置，用于根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵；最优化路径确定装置，用于根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数；匹配装置，用于根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

本发明的有益效果在于，提供了一种基于基因排序体系的多波地震资料时间域匹配方法及系统，可以从叠前纵波道集数据中反演得到横波信息数据，利用反演得到横波数据与转换波地震数据，计算最优的时间匹配路径，建立起时间转换目标函数，可以实现准确的高精度转换波地震资料时间匹配，同时还可以得到纵横波速度比模型，为多波地震资料联合反演与解释提供准确的资料，并提高储层油气检测的准确性。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为纵波与转换波传播路径示意图；

图2(a)为PS时间域的转换波地震数据示意图;

图2(b)为PP时间域的转换波地震数据示意图;

图2(c)为PP时间域的纵波地震数据示意图;

图2(d)为PP时间域的反演横波信息数据示意图；

图2(e)为纵波数据与转换波数据的时移相关谱示意图；

图2(f)为反演横波信息与转换波数据的时移相关谱示意图；

图3(a)为基因自相似性形成的一矩阵示意图，自相似的矩阵单元为1；

图3(b)、3（c）、3（d）为自相似性矩阵计算方法及过程示意图；

图3(e)、3（f）为基因对自相似性矩阵寻优结果示意图；

图4(a)为按地震子波时窗转换波地震道与纵波地震道相互的相似性矩阵示意图；

图4(b)为通过基因排序方法寻优求解后得到的时间匹配路径示意图；

图5(a)为测井数据计算得到纵波与横波的速度比示意图；

图5(b)从上到下依次为合成的纵波地震道(PP)、转换波地震数据(PS)、反演得到横波信息数据(PPS)和匹配到纵波旅行时后的转换波地震数据示意图；

图6(a)为时间匹配后的转换波地震剖面（左）与纵波地震剖面（右）在纵波旅行时间域的对比图；

图6(b)为纵横波速度比模型示意图；

图7为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配方法的流程图；

图8为图7中的步骤S105的具体流程图；

图9为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配系统的结构框图；

图10为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配系统中最优化路径确定装置的结构框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提出了一种新的基于基因排序体系的多波地震资料时间域匹配方法及系统，可以从叠前纵波道集数据中反演得到横波信息数据，利用反演得到横波数据与转换波地震数据，计算最优的时间匹配路径，建立起时间转换目标函数，解决多波地震勘探中转换波地震数据与纵波地震数据的时间域匹配问题，为纵波与转换波资料联合处理、解释和反演提供准确的资料，提高综合应用纵横波地震数据属性进行储层流体检测的准确度。

图7为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配方法的流程图，由图7可知，该方法具体包括：

S101：采集纵波地震道数据，假设输入的纵波地震道数据为：

pp=[w_pp(1),w_pp(2),…,w_pp(i),…,w_pp(n)]

其中，i表示为行方向，j是列方向。

S102：将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据。

现有的陆上多波地震资料，主要是以纵波作为震源，利用多分量数字检波器作为接收，得到纵波地震资料及转换波地震资料。因此，转换波地震资料与纵波地震资料之间的差别主要是反射波的波型、路径和传播时间不一样。假设地下一水平层（如图1所示），可以得到如图所示的纵波与转换波传播路径及其相互关系图。图中S为震源点，P为纵波接收点，V为转换波接收点，C为转换波形成的转换点，Vp为地层的纵波速度，Vs为地层的横波速度，h为地层的深度，x为纵波炮检点距，x_o为震源点到中心点（转换点）的距离，x_s为中心点O以转换波接收点的距离，x_ps为转换波接收点到纵波接收点的距离。

从图1中，可以看到纵波由于入射为纵波，反射也为纵波，其传播速度是一样的。而对于转换波地震资料，入射为纵波，反射为横波，其传播速度和路径都不一样。因此，反射纵波地震资料与转换波地震资料对于同一反射界面的同相轴在时间域内不同，二者的之间的时差关系可以写成：

$t_{\mathrm{pp}}=\frac{2h}{V_p}---\left(1\right)$

$t_{\mathrm{ps}}=\frac{h}{V_P}+\frac{h}{V_s}---\left(2\right)$

$\frac{t_{\mathrm{pp}}}{t_{\mathrm{ps}}}=\frac{\frac{2h}{V_p}}{\frac{h}{V_p}+\frac{h}{V_s}}=\frac{2}{1+\frac{V_p}{V_s}}---\left(3\right)$

从式（3）可以看出，纵波与转换波的时间匹配关系与纵波与横波速度比相关。

对于垂直入射的SS波，其零偏移距反射率具有如下形式：

$R_{\mathrm{SS}}\left(0\right)=\frac{{\mathrm{\ρ}}_1{\mathrm{\β}}_1-{\mathrm{\ρ}}_2{\mathrm{\β}}_2}{{\mathrm{\ρ}}_1{\mathrm{\β}}_1+{\mathrm{\ρ}}_2{\mathrm{\β}}_2}---\left(4\right)$

根据波阻抗的定义I_S=ρβ，上述关系式可以写成：

$R_{\mathrm{SS}}\left(0\right)=\frac{I_{S1}-I_{S2}}{I_{S1}+I_{S2}}=-\frac{1}{2}\frac{\mathrm{\Δ}{I}_S}{I_S}---\left(5\right)$

其中，I_S1，I_S2，I_S和ΔI_S分别表示反射界面两侧的波阻抗、平均波阻抗和波阻抗差。

进一步假设

ΔT_S≈Δρβ+ρΔβ   (6)

则有

$R_{\mathrm{SS}}\left(0\right)\≈-\frac{1}{2}(\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{\mathrm{\Δ\β}}{\mathrm{\β}})---\left(7\right)$

上式即是横波入射，反射得到的零偏移距情况下的横波的反射系数公式。当然，上式的反射系数的关系式，也可以从Aki和Richards等(1980)的SS波的反射系数近似公式的一阶近似得到。

Aki等(1980)给出在弱反差条件下P-SV波的反射系数表达式为

其中，为射线参数。

Δα=(α₂-α₁)  Δβ=(β₂-β₁)  Δρ=(ρ₂-ρ₁)

α=(α₂+α₁)/2  β=(β₂+β₁)/2  ρ=(ρ₂+ρ₁)/2

θ=(i₁+i₂)/2

上式中，下脚标1表示入射波和反射波所在的介质为1；下脚标2表示透射波所在的介质为2；α₁和α₂、β₁和β₂、ρ₁和ρ₂分别表示介质1和介质2中的纵波速度、横波速度和密度；i₁和j₁、i₂和j₂分别表示纵波的入射角和横波的反射角及纵波和横波的透射角。

如果对上式近似公式进一步假设，忽略p2项有

当θ，

都比较小，且β/α→0.5时，则有：

因此，与SS波零偏移距反射系数对比有

$\left\{\begin{array}{c}{R}_{\mathrm{PS}}\left(\mathrm{\θ}\right)\≈2\frac{\mathrm{\β}}{\mathrm{\α}}{R}_{\mathrm{SS}}\left(0\right)\sin 2\mathrm{\θ}\≈R_{\mathrm{SS}}\left(0\right)\sin 2\mathrm{\θ}\\ R_{\mathrm{SS}}\left(0\right)\≈\frac{\mathrm{\α}{R}_{\mathrm{PS}}\left(\mathrm{\θ}\right)}{2\mathrm{\β}\sin 2\mathrm{\θ}}\≈\frac{R_{\mathrm{PS}}\left(\mathrm{\θ}\right)}{\sin 2\mathrm{\θ}}\end{array}\right.---\left(11\right)$

为了能够建立反映P－SV转换波反射系数与纯横波反射率之间的高精度近似公式，考虑垂直入射的SS波反射率具有如下的表达式

$R_{\mathrm{SS}}\left(0\right)=\frac{{\mathrm{\ρ}}_1{\mathrm{\β}}_1-{\mathrm{\ρ}}_2{\mathrm{\β}}_2}{{\mathrm{\ρ}}_2{\mathrm{\β}}_2+{\mathrm{\ρ}}_1{\mathrm{\β}}_1}\≈\frac{1}{2}(\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{\mathrm{\Δ\β}}{\mathrm{\β}})---\left(12\right)$

对Aki ＆ Richards给出的P-SV波反射系数公式进行泰勒展开，并忽略射线参数的高阶项，可以将P-SV转换波反射系数表示为横波反射率的函数：

R_PS(θ)≈2γR_SS(0)sin2θ   (13)

从上面不同的近似公式可以看出，转换波反射系数与纯横波的反射系数之间只存在与角度相关的系数，反映的只是波的振幅大小的变化，与波的相位没有关系。因此利用相似系数作为目标函数不会影响其大小，完全可以改进相似系数的大小，从而达到更好的匹配效果。

因此，本发明建立起利用叠前纵波资料反演得到横波信息（PPS）剖面，改进纵波与转换波资料之间的波形相似征。其反演矩阵可以写成：

$\left(\begin{array}{cc}\underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i^2& \underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{b}_i\\ \underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{b}_i& \underset{i}{\overset{N}{\mathrm{\Σ}}}{b}_i^2\end{array}\right)\left(\begin{array}{c}{I}_{\mathrm{PP}}\\ I_{\mathrm{PPS}}\end{array}\right)=\left(\begin{array}{c}\underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{X}_i\\ \underset{i}{\overset{N}{\mathrm{\Σ}}}{b}_i{X}_i\end{array}\right),$ 即AI=d   （14）

式中，A为已知参数的矩阵，d为纵波地震道数据向量，I为反演未知参数，具体的，i为叠前纵波数据偏移距道数，N为叠前数据偏移距总道数，a_i为叠前反演公式已知系数向量，b_i为叠前反演公式已知系数向量，I_PP为反演的未知参数纵波阻抗，I_PPS反演的未知参数横波信息，X_i为叠前反演资料的输入的纵波数据

图2(a)至图2(f)为纵波地震数据（PP）、转换波地震数据（PS）及其与反演横波信息（PPS）之间的波形相似性及其相关谱比较。其中，图2(a)为PS时间域的转换波地震数据;图2(b)为PP时间域的转换波地震数据;图2(c)为PP时间域的纵波地震数据;图2(d)为PP时间域的反演横波信息数据；图2(e)为纵波数据与转换波数据的时移相关谱;图2(f)为反演横波信息与转换波数据的时移相关谱。纵波数据与转换波数据之间的波形特征差，其相似性系数只有0.473，而反演横波信息与转换波数据之间的波形特征好，其相似性系数达到0.933。也即通过线性化反演上述的线性方程组，可以得到横波信息剖面（I_PPS）。对比反演的横波信息和原始的纵波资料与转换波资料的相似性，可以看到横波信息具有与转换波资料更高的相似性（图2(a)至图2(f)）。图2(a)至图2(f)中反演横波信息与转换波数据的时移相关谱能量集中，结果更稳定，同时波形特征相关性，具有更高的波形相似性。

S103：采集转换波地震道数据，假设输入的转换波地震道数据为：

ps=[w_ps(1),w_ps(2),…,w_ps(j),…,w_ps(m)]

如图4(a)、图4(b)所示，式中i表示为行方向，j是列方向。

S104：根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵。相似性矩阵如下所示：

其中， $c_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{k}{\mathrm{\Σ}}}{P}_i\left(t\right)C_j\left(t\right)}{\sqrt{\underset{t=1}{\overset{k}{\mathrm{\Σ}}}{P}_i^2\left(t\right)\underset{t=1}{\overset{k}{\mathrm{\Σ}}}{C}_j^2\left(t\right)}}$

其中，P(t)为纵波的地震记录，C(t)为转换波的地震记录，i为转换波的时窗数，j为纵波的时间窗数，0<i<m，0<j<n，m为转换波的时窗长度，n为纵波的时窗长度，通常m>n，k为时窗大小，E_cc为相似性矩阵。

S105：根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。本发明是基于基因排序体系高精度多波地震资料匹配方法，因此，首先介绍基因排序的非线性求解算法。

根据纵波与转换波资料的波传播路径原理(图1)，对于同一地层的反射同相轴，转换波资料与纵波资料时差关系为：

$\frac{t_{\mathrm{pp}}}{t_{\mathrm{ps}}}=\frac{2}{1+{\mathrm{\&gamma;}}_0}---\left(15\right)$

其中，γ₀为纵横波速度比。对于转换波地震资料的高精度自动匹配，需要寻找非线性的全局最优算法求解方法。因为，由于地下介质的复杂性及地震子波的带限特性，转换波与纵波的匹配过程中的速度比γ₀求解具强非线性特征。

(1)基因排序的非线性求解算法

在人类基因工程中，对于当染色体的配对分析中，也需要非线性的全局最优算法。例如：对于下列两个基因序列：ABCNYRQCLCRPM和AYCYNRCKCRBP，对上述基因进行排序的非线性最优算法可以用图3所示的流程来表示。图3(a)至图3(f)为基因序列对排序算法流程示意图。图3(a)为基因自相似性形成的一矩阵，自相似的矩阵单元为1；图3(b)、图3（c）、图3（d)为自相似性矩阵计算方法及过程；图3(e)、图3（f)为基因对自相似性矩阵寻优结果。

在图3(a)至图3(f)所示的基因排序中，两个不同长度的字符串形成矩阵，沿横向和纵向，首先在矩阵单元中分配1，如果列字母和行字母相同时，否则为空，这样代表基因对的字符排序问题就转化为数值问题，这样就形成相似性矩阵，如图3（a）。从相似性矩阵，从右下角开始，形成分值排序矩阵。具体计算方法是当前行列的值加上右下行列中的最大值(如图3（b）、图3（c）)，通过遍历相似性矩阵中的所有单元，最后形成如图3（d）所示的矩阵。最终的分值矩阵就展示了有多少个字符匹配及其最优的排序，这也是一全部相似性的数值计算。最后一步是追踪，从最左上角的最大值开始，沿矩阵对角线方向寻找最大值路径，这样可以得到整个路径上的最大值，当然可能有几个最大值路径（如图3（e）、图3（f）），在这里的例子中具有两个最大值路径。该方法的最大优点是能找到全局的最大值，也可以分成几个部分进行局部最大值路径寻优。

(2)转换波与纵波地震数据时间匹配

本发明中，建立基于基因排序体系的转换波与纵波时间匹配方法，达到高精度自动时间匹配。图8为图7中的步骤S105的具体流程图，由图8可知，步骤S105具体包括：

S201：将所述的转换波地震道数据按横向方向排列；

S202：将所述的纵波地震道数据按纵向方向排列；

S203：获取预先设定的时间间隔，此处的时间间隔即为选择的合适的时间间隔。

S204：根据所述的时间间隔以及所述的相似性矩阵确定所述转换波地震道数据与纵波地震道数据之间的相似系数矩阵；

S205：根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。

如图4(a)至图4(b)所示，将转换波地震道数据、纵波地震道数据按横向与纵向方向排列，选择合适的时间间隔计算转换波地震道与纵波地震道数据之间的相似系数矩阵（图4a）,最后通过基因排序方法进行非线性最优求解得到时间匹配路径，通过时间匹配路径进行自动匹配。图4(a)至图4(b)为基于基因排序体系的转换波与纵波资料匹配目标函数求解示意图。图4(a)为按地震子波时窗进行转换波地震道与纵波地震道相互的相似性矩阵生成；图4(b)为通过基因排序方法寻优求解后得到的时间匹配路径。

将纵波地震道pp和转换波地震道ps分别表示成在不同时间（i或j）处的地震子波形式。根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，即对于时间匹配的最优路径求解，就是求解在不同的时间处地震子波之间的相似系数。因此，对于相似性矩阵上每一点，我们最终是寻找下一点(E)，右边点(F)，或是沿对角点(G)中的最大值。为了稳定转换波地震道与纵波地震道的时间匹配路径的求解过程，我们在计算过程中加入惩罚因子来调节由于地震子波交叉或噪声的影响，包括两个因子p和q分别进行沿对角线附近的惩罚。最终的最优化过程可以写成：

V(i,j)=max[G(i,j),F(i,j),E(i,j)]   (17)

G(i,j)=σ(i,j)+V(i+1,j+1)   (18)

F(i,j)=-(p+q)+max[V(i+1,j),F(i+1,j)+p]   (19)

E(i,j)=-(p+q)+max[V(i,j+1),E(i,j+1)+q]   (20)

上述式中，E为所定义的矩阵，V是矩阵中当前点的下一行点，G是矩阵当前点沿对角线上的点，F为矩阵当前点右边一列点，p为惩罚因子参数，q为惩罚因子参数。并且，

$\mathrm{\&sigma;}(i,j)=\left\{\begin{array}{cc}1& i=j\\ 0& i\&NotEqual;j\end{array}\right.$

通过递归的方法，利用上述方程式可以得到最优的时间匹配路径，根据时间路径建立起将转换波地震道转换波到纵波旅行时间域的地震道。

假设纵波和转换波反射系数的横向位置在偏移后是正确的，转换波匹配仅处理垂向的时间或是深度的时差。同时，纵波和转换波反资料中只含有一次反射波有效信号，地表相关的多波和层间多次波等已经消除。因此，转换波与纵波匹配的目标函数函数就是最大化纵波与转换波成像剖面之间的相似性函数，以补偿纵波与转换波反射系数之间的差异。

$\max \left(\frac{\mathrm{\&Sigma;P}\left(t\right)C\left(w\left(t\right)\right)}{\sqrt{\mathrm{\&Sigma;}{P}^2\left(t\right)\mathrm{\&Sigma;}{C}^2\left(w\left(t\right)\right)}}\right)---\left(16\right)$

其中：且w(t)是转换波时间(T_ps)到纵波时间(T_pp)的转换函数，由上步求取得到，P(t)为纵波的地震记录。

S106：根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

如此，本发明提供的基于基因排序体系的多波地震资料时间域匹配方法可以实现准确的高精度转换波地震资料时间匹配，同时还可以得到纵横波速度比模型，为多波地震资料联合反演与解释提供准确的资料，并提高储层油气检测的准确性。

图9为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配系统的结构框图，由图9可知，该系统具体包括：

纵波地震道数据采集装置100，用于采集纵波地震道数据，假设输入的纵波地震道数据为：

pp=[w_pp(1),w_pp(2),…,w_pp(i),…,w_pp(n)]

其中，i表示为行方向，j是列方向。

反演装置200，用于将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据。

现有的陆上多波地震资料，主要是以纵波作为震源，利用多分量数字检波器作为接收，得到纵波地震资料及转换波地震资料。因此，转换波地震资料与纵波地震资料之间的差别主要是反射波的波型、路径和传播时间不一样。假设地下一水平层（如图1所示），可以得到如图所示的纵波与转换波传播路径及其相互关系图。图中S为震源点，P为纵波接收点，V为转换波接收点，C为转换波形成的转换点，Vp为地层的纵波速度，Vs为地层的横波速度，h为地层的深度，x为纵波炮检点距，x_o为震源点到中心点（转换点）的距离，x_s为中心点O以转换波接收点的距离，x_ps为转换波接收点到纵波接收点的距离。

从图1中，可以看到纵波由于入射为纵波，反射也为纵波，其传播速度是一样的。而对于转换波地震资料，入射为纵波，反射为横波，其传播速度和路径都不一样。因此，反射纵波地震资料与转换波地震资料对于同一反射界面的同相轴在时间域内不同，二者的之间的时差关系可以写成公式（1）、（2）、（3）。从式（3）可以看出，纵波与转换波的时间匹配关系与纵波与横波速度比相关。

对于垂直入射的SS波，其零偏移距反射率具有公式（4）所示的形式。根据波阻抗的定义I_S=ρβ，上述关系式可以写成公式(5)，其中，I_S1，I_S2，I_S和ΔI_S分别表示反射界面两侧的波阻抗、平均波阻抗和波阻抗差。进一步假设ΔI_S≈Δρβ+ρΔβ，则有 $R_{\mathrm{SS}}\left(0\right)=-\frac{1}{2}(\frac{\mathrm{\&Delta;\&rho;}}{\mathrm{\&rho;}}+\frac{\mathrm{\&Delta;\&beta;}}{\mathrm{\&beta;}}).$

上式即是横波入射，反射得到的零偏移距情况下的横波的反射系数公式。当然，上式的反射系数的关系式，也可以从Aki和Richards等(1980)的SS波的反射系数近似公式的一阶近似得到。Aki等(1980)给出在弱反差条件下P-SV波的反射系数表达式如公式(8)所示。其中，

为射线参数。

Δα=(α₂-α₁)  Δβ=(β₂-β₁)  Δρ=(ρ₂-ρ₁)

α=(α₂+α₁)/2  β=(β₂+β₁)/2  ρ=(ρ₂+ρ₁)/2

θ=(i₁+i₂)/2

上式中，下脚标1表示入射波和反射波所在的介质为1；下脚标2表示透射波所在的介质为2；α₁和α₂、β₁和β₂、ρ₁和ρ₂分别表示介质1和介质2中的纵波速度、横波速度和密度；i₁和j₁、i₂和j₂分别表示纵波的入射角和横波的反射角及纵波和横波的透射角。

如果对上式近似公式进一步假设，忽略p²项有公式(9)，当θ，都比较小，且β/α→0.5时，则有公式(10)。因此，与SS波零偏移距反射系数对比有公式(11)。

为了能够建立反映P－SV转换波反射系数与纯横波反射率之间的高精度近似公式，考虑垂直入射的SS波反射率具有公式(12)的表达式。对Aki ＆Richards给出的P-SV波反射系数公式进行泰勒展开，并忽略射线参数的高阶项，可以将P-SV转换波反射系数表示为横波反射率的函数，如公式(13)所示。

从上面不同的近似公式可以看出，转换波反射系数与纯横波的反射系数之间只存在与角度相关的系数，反映的只是波的振幅大小的变化，与波的相位没有关系。因此利用相似系数作为目标函数不会影响其大小，完全可以改进相似系数的大小，从而达到更好的匹配效果。

因此，本发明建立起利用叠前纵波资料反演得到横波信息（PPS）剖面，改进纵波与转换波资料之间的波形相似征。其反演矩阵可以写成公式（14）。

式中，A为已知参数的矩阵，d为纵波地震道数据向量，I为反演未知参数，具体的，i为叠前纵波数据偏移距道数，N为叠前数据偏移距总道数，a_i为叠前反演公式已知系数向量，b_i为叠前反演公式已知系数向量，I_PP为反演的未知参数纵波阻抗，I_PPS反演的未知参数横波信息，X_i为叠前反演资料的输入的纵波数据

图2(a)至图2(f)为纵波地震数据（PP）、转换波地震数据（PS）及其与反演横波信息（PPS）之间的波形相似性及其相关谱比较。图2(a)为PS时间域的转换波地震数据;图2(b)为PP时间域的转换波地震数据;图2(c)为PP时间域的纵波地震数据;图2(d)为PP时间域的反演横波信息数据；图2(e)为纵波数据与转换波数据的时移相关谱;图2(f)为反演横波信息与转换波数据的时移相关谱。纵波数据与转换波数据之间的波形特征差，其相似性系数只有0.473，而反演横波信息与转换波数据之间的波形特征好，其相似性系数达到0.933。也即通过线性化反演上述的线性方程组，可以得到横波信息剖面（I_PPS）。对比反演的横波信息和原始的纵波资料与转换波资料的相似性，可以看到横波信息具有与转换波资料更高的相似性（图2(a)至图2(f)所示）。图2(a)至图2(f)中反演横波信息与转换波数据的时移相关谱能量集中，结果更稳定，同时波形特征相关性，具有更高的波形相似性。

转换波地震道数据采集装置300，用于采集转换波地震道数据，假设输入的转换波地震道数据为：

ps=[w_ps(1),w_ps(2),…,w_ps(j),…,w_ps(m)]

如图4(a)、图4(b)所示，式中i表示为行方向，j是列方向。

相似性矩阵建立装置400，用于根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵。

最优化路径确定装置500，用于根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。本发明是基于基因排序体系高精度多波地震资料匹配系统，因此，首先介绍基因排序的非线性求解算法。

根据纵波与转换波资料的波传播路径原理(图1)，对于同一地层的反射同相轴，转换波资料与纵波资料时差关系如公式（15）所示。其中，γ₀为纵横波速度比。对于转换波地震资料的高精度自动匹配，需要寻找非线性的全局最优算法求解方法。因为，由于地下介质的复杂性及地震子波的带限特性，转换波与纵波的匹配过程中的速度比γ₀求解具强非线性特征。

(1)基因排序的非线性求解算法

在人类基因工程中，对于当染色体的配对分析中，也需要非线性的全局最优算法。例如：对于下列两个基因序列：ABCNYRQCLCRPM和AYCYNRCKCRBP，对上述基因进行排序的非线性最优算法可以用图3所示的流程来表示。图3(a)至图3(f)为基因序列对排序算法流程示意图。图3(a)为基因自相似性形成的一矩阵，自相似的矩阵单元为1；图3(b)、图3（c）、图3（d)为自相似性矩阵计算方法及过程；图3(e)、图3（f)为基因对自相似性矩阵寻优结果。

在图3(a)至图3(f)所示的基因排序中，两个不同长度的字符串形成矩阵，沿横向和纵向，首先在矩阵单元中分配1，如果列字母和行字母相同时，否则为空，这样代表基因对的字符排序问题就转化为数值问题，这样就形成相似性矩阵，如图3（a）。从相似性矩阵，从右下角开始，形成分值排序矩阵。具体计算方法是当前行列的值加上右下行列中的最大值(如图3（b）、图3（c）)，通过遍历相似性矩阵中的所有单元，最后形成如图3（d）所示的矩阵。最终的分值矩阵就展示了有多少个字符匹配及其最优的排序，这也是一全部相似性的数值计算。最后一步是追踪，从最左上角的最大值开始，沿矩阵对角线方向寻找最大值路径，这样可以得到整个路径上的最大值，当然可能有几个最大值路径（如图3（e）、图3（f）），在这里的例子中具有两个最大值路径。该方法的最大优点是能找到全局的最大值，也可以分成几个部分进行局部最大值路径寻优。

(2)转换波与纵波地震数据时间匹配

本发明中，建立基于基因排序体系的转换波与纵波时间匹配方法，达到高精度自动时间匹配。图10为本发明实施例提供的一种基于基因排序体系的多波地震资料时间域匹配系统中最优化路径确定装置的结构框图，由图10可知，最优化路径确定装置500具体包括：

横向排列单元501，用于将所述的转换波地震道数据按横向方向排列；

纵向排列单元502，用于将所述的纵波地震道数据按纵向方向排列；

时间间隔获取单元503，用于获取预先设定的时间间隔；

相似系数矩阵确定单元504，用于根据所述的时间间隔以及所述的相似性矩阵确定所述转换波地震道数据与纵波地震道数据之间的相似系数矩阵；

最优化路径确定单元505，用于根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。

如图4(a)至图4(b)所示，将转换波地震道数据、纵波地震道数据按横向与纵向方向排列，选择合适的时间间隔计算转换波地震道与纵波地震道数据之间的相似系数矩阵（图4a）,最后通过基因排序方法进行非线性最优求解得到时间匹配路径，通过时间匹配路径进行自动匹配。图4(a)至图4(b)为基于基因排序体系的转换波与纵波资料匹配目标函数求解示意图。图4(a)为按地震子波时窗进行转换波地震道与纵波地震道相互的相似性矩阵生成；图4(b)为通过基因排序方法寻优求解后得到的时间匹配路径。

将纵波地震道pp和转换波地震道ps分别表示成在不同时间（i或j）处的地震子波形式。根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，即对于时间匹配的最优路径求解，就是求解在不同的时间处地震子波之间的相似系数。因此，对于相似性矩阵上每一点，我们最终是寻找下一点(E)，右边点(F)，或是沿对角点(G)中的最大值。为了稳定转换波地震道与纵波地震道的时间匹配路径的求解过程，我们在计算过程中加入惩罚因子来调节由于地震子波交叉或噪声的影响，包括两个因子p和q分别进行沿对角线附近的惩罚。最终的最优化过程可以写成公式（17）至（20）。

上述式中，E为所定义的矩阵，V是矩阵中当前点的下一行点，G是矩阵当前点沿对角线上的点，F为矩阵当前点右边一列点，p为惩罚因子参数，q为惩罚因子参数。并且，

$\mathrm{\&sigma;}(i,j)=\left\{\begin{array}{cc}1& i=j\\ 0& i\&NotEqual;j\end{array}\right.$

通过递归的方法，利用上述方程式可以得到最优的时间匹配路径，根据时间路径建立起将转换波地震道转换波到纵波旅行时间域的地震道。

假设纵波和转换波反射系数的横向位置在偏移后是正确的，转换波匹配仅处理垂向的时间或是深度的时差。同时，纵波和转换波反资料中只含有一次反射波有效信号，地表相关的多波和层间多次波等已经消除。因此，转换波与纵波匹配的目标函数函数就是最大化纵波与转换波成像剖面之间的相似性函数，以补偿纵波与转换波反射系数之间的差异。

公式（16）中的

且w(t)是转换波时间(T_ps)到纵波时间(T_pp)的转换函数，由上步求取得到，P(t)为纵波的地震记录。

匹配装置600，用于根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

如此，本发明提供的基于基因排序体系的多波地震资料时间域匹配系统可以实现准确的高精度转换波地震资料时间匹配，同时还可以得到纵横波速度比模型，为多波地震资料联合反演与解释提供准确的资料，并提高储层油气检测的准确性。

下面集合具体的实施例，详细介绍本发明提供的基于基因排序体系的多波地震资料时间域匹配方法。

在本实例中，我们选取一过井点处的多波地震资料，且该井有横波测井数据。利用测井得到纵波与横波数据，我们建立了合成地震数据来展示反演横波信息提高波形相似。图5（a）展示了测井数据计算得到纵波与横波的速度比，图中点线为测井数据计算得到的初始速度比，实线为时间匹配后得到最终速度比。

图5（b）的地震道数据从上到下依次为：合成的纵波地震道(PP)、转换波地震数据(PS)、反演得到横波信息数据(PPS)和匹配到纵波旅行时后的转换波地震数据。对比分析图中的地震波形，可以看到反演得到的横波信息PPS与转换波地震道PS之间具有更好的波形相似性，相比于原始的纵波地震数据PP。证明了本发明中的反演方法的正确性，同时验证了直接基于纵波与转换波资料波形相似性匹配的不准确性。

采用本发明中的技术流程，我们实现了实际三维多波地震数据的匹配，图6（a）、图6（b）中展示了最终的匹配效果。图6（a）中采用变密度的方式显示了时间匹配后的转换波剖面（图中左边）和纵波地震剖面（图中右边），对比可以发现，经过匹配后的转换波剖面与纵波地震剖面具有很好的同相轴排列。同时，利用本技术流程，我们还可建立准确的纵波与横波速度比模型，图6（b）显示实例中，三维数据体时间匹配后的速度比模型，利用纵横波速度比可以进行多波地震数据的联合反演与解释，提高油气检测能力。

综上所述，本发明提供了一种基于基因排序体系的多波地震资料时间域匹配方法及系统，主要应用于多波地震资料联合处理与联合反演，解决多波地震资料的深入应用研究。该项技术的特点主要有三点：全局非线性寻优方法；高精度波形相似性匹配；全自动匹配技术。可以从叠前纵波道集数据中反演得到横波信息数据，利用反演得到横波数据与转换波地震数据，计算最优的时间匹配路径，建立起时间转换目标函数，可以实现准确的高精度转换波地震资料时间匹配，同时还可以得到纵横波速度比模型，为多波地震资料联合反演与解释提供准确的资料，并提高储层油气检测的准确性。

本发明是参照根据本发明实施例的方法、设备（系统）、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于基因排序体系的多波地震资料时间域匹配方法，其特征是，所述的方法包括：

采集纵波地震道数据；

将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据；

采集转换波地震道数据；

根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵；

根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数；

根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

2.根据权利要求1所述的方法，其特征是，将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据通过如下公式进行：

AI=d

其中，A为已知参数的矩阵，d为纵波地震道数据向量，I为反演的未知参数。

3.根据权利要求2所述的方法，其特征是，根据所述的纵波地震道数据以及所述的转换波地震道数据建立的相似性矩阵如下所示：

其中， $c_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{P}_i\left(t\right)C_j\left(t\right)}{\sqrt{\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{P}_i^2\left(t\right)\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{C}_j^2\left(t\right)}}$

其中，P(t)为纵波的地震记录，C(t)为转换波的地震记录，i为转换波的时窗数，j为纵波的时间窗数，0<i<m，0<j<n，m为转换波的时窗长度，n为纵波的时窗长度，m>n，k为时窗大小，E_cc为相似性矩阵。

4.根据权利要求3所述的方法，其特征是，根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径具体包括：

将所述的转换波地震道数据按横向方向排列；

将所述的纵波地震道数据按纵向方向排列；

获取预先设定的时间间隔；

根据所述的时间间隔以及所述的相似性矩阵确定所述转换波地震道数据与纵波地震道数据之间的相似系数矩阵；

根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。

5.根据权利要求4所述的方法，其特征是，根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径通过如下公式进行：

V(i,j)=max[G(i,j),F(i,j),E(i,j)]

G(i,j)=σ(i,j)+V(i+1,j+1)

F(i,j)=-(p+q)+max[V(i+1,j),F(i+1,j)+p]

E(i,j)=-(p+q)+max[V(i,j+1),E(i,j+1)+q]

其中，E为所定义的矩阵，V是矩阵中当前点的下一行点，G是矩阵当前点沿对角线上的点，F为矩阵当前点右边一列点，p为惩罚因子参数，q为惩罚因子参数，且

$\mathrm{\&sigma;}(i,j)=\left\{\begin{array}{cc}1& i=j\\ 0& i\&NotEqual;j\end{array}\right..$

6.一种基于基因排序体系的多波地震资料时间域匹配系统，其特征是，所述的系统包括：

纵波地震道数据采集装置，用于采集纵波地震道数据；

反演装置，用于将所述的纵波地震道数据进行反演，得到纵波时间域的横波地震道数据；

转换波地震道数据采集装置，用于采集转换波地震道数据；

相似性矩阵建立装置，用于根据所述的横波地震道数据以及所述的转换波地震道数据建立相似性矩阵；

最优化路径确定装置，用于根据所述的相似性矩阵以及基因排序的非线性求解算法确定时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数；

匹配装置，用于根据所述的时间转换函数将所述的转换波地震道数据匹配到纵波时间域。

7.根据权利要求6所述的系统，其特征是，所述的反演装置通过如下公式进行：

AI=d

其中，A为已知参数的矩阵，d为纵波地震道数据向量，I为反演的未知参数。

8.根据权利要求7所述的系统，其特征是，所述的相似性矩阵建立装置建立的相似性矩阵如下所示：

其中， $c_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{P}_i\left(t\right)C_j\left(t\right)}{\sqrt{\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{P}_i^2\left(t\right)\underset{t=1}{\overset{k}{\mathrm{\&Sigma;}}}{C}_j^2\left(t\right)}}$

其中，P(t)为纵波的地震记录，C(t)为转换波的地震记录，i为转换波的时窗数，j为纵波的时间窗数，0<i<m，0<j<n，m为转换波的时窗长度，n为纵波的时窗长度，m>n，k为时窗大小，E_cc为相似性矩阵。

9.根据权利要求8所述的系统，其特征是，所述的最优化路径确定装置具体包括：

横向排列单元，用于将所述的转换波地震道数据按横向方向排列；

纵向排列单元，用于将所述的纵波地震道数据按纵向方向排列；

时间间隔获取单元，用于获取预先设定的时间间隔；

相似系数矩阵确定单元，用于根据所述的时间间隔以及所述的相似性矩阵确定所述转换波地震道数据与纵波地震道数据之间的相似系数矩阵；

最优化路径确定单元，用于根据基因排序的非线性求解算法进行非线性最优求解，得到时间匹配的最优化路径，并建立起纵波与转换波之间的时间转换函数。

10.根据权利要求9所述的系统，其特征是，所述的最优化路径求解单元通过如下公式进行：

V(i,j)=max[G(i,j),F(i,j),E(i,j)]

G(i,j)=σ(i,j)+V(i+1,j+1)

F(i,j)=-(p+q)+max[V(i+1,j),F(i+1,j)+p]

E(i,j)=-(p+q)+max[V(i,j+1),E(i,j+1)+q]

其中，E为所定义的矩阵，V是矩阵中当前点的下一行点，G是矩阵当前点沿对角线上的点，F为矩阵当前点右边一列点，p为惩罚因子参数，q为惩罚因子参数，且

$\mathrm{\&sigma;}(i,j)=\left\{\begin{array}{cc}1& i=j\\ 0& i\&NotEqual;j\end{array}\right..$

Images