CN111856568A - 基于mwv模型的频率域多道反射系数联合反演方法及系统 - Google Patents

Abstract

本发明提供了一种基于MWV模型的频率域多道反射系数联合反演方法及系统，该方法包括：获取经过预处理的叠后地震资料；根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；从所述叠后地震资料中提取地震子波；基于傅里叶变换提取所述地震子波的频谱信息；根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。本发明有效提高了频率域反射系数反演结果的精度和稳定性。

Description

基于MWV模型的频率域多道反射系数联合反演方法及系统

技术领域

本发明涉及油气田地震勘探领域，尤其涉及一种基于MMV模型的频率域多道反射系数联合反演方法及系统。

背景技术

在油气田勘探开发领域，地震勘探一直是探明地下油藏特征、开展储层精细解释的重要技术手段。地震资料具有覆盖范围广、横向连续性强等特点，其在识别地下薄层，特别是含油气地层在地下空间的三维展布特征等方面具有不可替代的优势。然而在实际生产中，由于地震资料频带宽度有限，制约了地震数据薄层分辨能力。因此基于地震资料开展反射系数反演拓宽其有效频带宽度是提高地震分辨率的关键，为后续薄层识别及储层精细解释奠定基础。

地震反射系数反演可大致分为时间域反演方法和频率域反演方法，其都需要输入地震子波相关信息，并通过剔除带限子波影响从而提高资料分辨率。两类方法不同之处则主要体现在：时间域反射系数反演方法基于时域褶积模型，通过地震子波时间序列构建褶积矩阵，地震波形可视为子波的线性叠加，对应加权系数即为反射系数。频率域反射系数反演方法则是通过构建地震资料频谱与其子波频谱的线性响应关系，视其加权系数为地震道对应的反射系数。对于白噪声，两类反演方法都有较好的压制效果。然而对于有色地震噪声，其有效功率谱通常集中在某一频带范围内，频率域反演方法可通过优选合适频带范围予以处理从而降低有色噪声的影响，而时间域反演方法则难以灵活规避这一问题。此外，频率域反演方法可通过调整频谱采样密度控制正演算子矩阵规模，从而平衡反演精度与运算效率，而时间域反演算法中正演算子矩阵规模受地震道和地震子波的时间长度制约，限制了方法的灵活性。频率域反演方法从频谱角度分析求解问题，更符合拓宽频谱提高分辨率的宗旨。

对于稀疏反射系数反演，其本质是求解稀疏约束即L₀范数约束的最优化问题。由于L₀范数求解属于NP难问题，因此传统基追踪理论将其简化为L₁范数约束的线性规划(LP)问题求解。而稀疏贝叶斯学习理论将L₀范数约束问题转换为一种贝叶斯框架下的稀疏训练问题，其与基追踪的差别在于基追踪理论中L₁范数的权系数为一个预先给定的常数，而稀疏贝叶斯学习理论则对每个基向量设置独立的正则参数，并通过一定的学习机制训练最优正则参数组合。上述理论在频率域反射系数反演中仍有以下问题亟待进一步优化：

1)反演结果的稀疏性。由于L₁范数的稀疏表示能力天然不足，基追踪理论将L₀范数近似转换为L₁范数所求得的解的稀疏性通常不足，其反演结果通常存在结构误差，即所获取的稀疏解并不一定对应于全局最优解。而稀疏贝叶斯学习理论通过贝叶斯学习机制训练开展正则参数训练，较大程度优化了解的稀疏性，但由于当前参数学习准则的精度不足且未考虑反射系数的横向连续性等因素，其求解精度仍有待提高。

2)反演算法的抗噪性。在含噪情况下，基于基追踪理论的稀疏反射系数反演结果的稀疏性会遭到破坏。而对于稀疏贝叶斯学习理论，在含噪情况下其所求的稀疏反射系数与真实情况通常存在一定位置偏差，且偏差程度随噪声的增加而扩大；同时资料噪声也在一定程度上破坏解的稀疏性。

3)反射系数的横向连续性。地震反射系数在地下空间一定范围内表现为垂向稀疏性和横向连续性，传统反演理论，包括基追踪和稀疏贝叶斯学习理论，仅考虑了垂向稀疏性而对其横向连续性则均未加以约束。横向连续性约束的缺失会进一步破坏结果的稀疏性和抗噪性。

发明内容

针对上述问题，本发明提出了一种基于MMV模型的频率域多道反射系数联合反演方法及系统，旨在通过MMV模型定量引入地震反射系数的横向连续性约束，在此基础上联合块稀疏贝叶斯理论在频率域反射系数反演算法中同时加入垂向稀疏性和横向连续性约束，并实现相关参数学习准则的优化，从而提高反演结果的精度和空间连续性，增强地震薄层识别能力，降低油气藏勘探开发风险。

本发明一方面提供了一种基于MWV模型的频率域多道反射系数联合反演方法，所述方法包括：

步骤S1、获取经过预处理的叠后地震资料；

步骤S2、根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；

步骤S3、从所述叠后地震资料中提取地震子波；

步骤S4、基于傅里叶变换提取所述地震子波的频谱信息；

步骤S5、根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；

步骤S6、通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；

步骤S7、基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

在一实施例中，所述预处理包括静校正、去噪、振幅补偿、速度分析、动校正、偏移中的至少一种。

在一实施例中，所述频率域多道反射系数联合反演参数先验信息至少包括：有效频谱范围、时间采样密度、频率采样密度、噪声方差、收敛阈值及最大迭代次数。

在一实施例中，在步骤S3中，如所述叠后地震资料中包含测井资料，则基于所述测井资料中的测井曲线及井旁地震道进行地震子波提取；如所述叠后地震资料中不包含测井资料，则对所述叠后地震资料进行统计性子波提取以获取地震子波。

在一实施例中，在步骤S5中，根据所述叠后地震资料的时间序列长度及时间采样密度确定所述正演算子矩阵的列数，根据所述叠后地震资料的有效频谱范围及频率采样密度确定所述正演算子矩阵的行数。

在一实施例中，所述步骤S7包括：

根据所述地震子波的频谱信息以及所选地震道的频谱信息，基于频谱响应关系计算反射系数频谱；

结合所述正演算子矩阵构建频率域多道反射系数联合反演的目标函数；

通过MMV模型量化叠后地震资料的横向连续性，联合块稀疏贝叶斯学习理论迭代求解所述目标函数，以获取反射系数模型参数。

另一方面，本发明实施例还提供了一种基于MWV模型的频率域多道反射系数联合反演系统，所述系统包括：

叠后地震资料获取单元，用于获取经过预处理的叠后地震资料；

信息获取单元，用于根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；

地震子波提取单元，用于从所述叠后地震资料中提取地震子波；

第一频谱信息获取单元，用于基于傅里叶变换提取所述地震子波的频谱信息；

正演算子矩阵构建单元，用于根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；

第二频谱信息获取单元，用于通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；

反演单元，用于基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

在一实施例中，所述预处理包括静校正、去噪、振幅补偿、速度分析、动校正、偏移中的至少一种。

在一实施例中，所述频率域多道反射系数联合反演参数先验信息至少包括：有效频谱范围、时间采样密度、频率采样密度、噪声方差、收敛阈值及最大迭代次数。

在一实施例中，如所述叠后地震资料中包含测井资料，则所述地震子波提取单元基于所述测井资料中的测井曲线及井旁地震道进行地震子波提取；如所述叠后地震资料中不包含测井资料，则所述地震子波提取单元对所述叠后地震资料进行统计性子波提取以获取地震子波。

在一实施例中，所述正演算子矩阵构建单元根据所述地震资料的时间序列长度及时间采样密度确定所述正演算子矩阵的列数，根据所述地震资料的有效频谱范围及频率采样密度确定所述正演算子矩阵的行数。

在一实施例中，所述反演单元包括：

反射系数频谱获取模块，用于根据所述地震子波的频谱信息以及所选地震道的频谱信息，基于频谱响应关系计算反射系数频谱；

目标函数构建模块，用于结合所述正演算子矩阵构建频率域多道反射系数联合反演的目标函数；

反演计算模块，用于通过MMV模型量化叠后地震资料的横向连续性，联合块稀疏贝叶斯学习理论迭代求解所述目标函数，以获取反射系数模型参数。

本发明由于采取以上技术方案，其具有以下优点：1)本技术方案所述的基于MMV模型的频率域多道反射系数联合反演方法，与当前现有技术相比，方案通过MMV模型量化引入反射系数的横向连续性特征，可大幅度提高反演结果的空间连续性；2)本技术方案联合块稀疏贝叶斯学习理论求解多个模型参数，包括噪声方差、正则参数、协方差矩阵等，实现多道参数的联合优化并提高了参数优化精度；3)本技术方案求解基于MMV模型的频率域多道反射系数联合反演目标函数，反演结果较传统反演方法结果稀疏性更好，更接近全局最优解；4)本技术方案的抗噪能力优于传统理论；5)本技术方案同时考虑了反射系数的垂向稀疏性和横向连续性，能够有效提高频率域反射系数反演结果的精度和稳定性。

附图说明

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

图1是本发明实施例提供的基于MMV模型的频率域多道反射系数联合反演方法的流程示意图；

图2a、图2b分别是本发明实施例中输入的二维楔状地层模型的反射系数及其频谱信息；

图2c、图2d分别是本发明实施例中输入的二维楔状地层模型在30Hz处的零相位雷克子波及其频谱信息；

图2e、图2f分别是本发明实施例中输入的二维楔状地层模型的合成地震记录及其频谱信息；

图3a、图3b分别是以图2e中二维楔状地层模型的合成地震记录为输入得到的基追踪反演结果及对应误差；

图3c、图3d分别是以图2e中二维楔状地层模型的合成地震记录为输入得到的基于稀疏贝叶斯学习的反演结果及对应误差；

图3e、图3f分别是以图2e中二维楔状地层模型的合成地震记录为输入得到的基于MMV模型的频率域多道反射系数联合反演结果及对应误差，其中，相关道数L＝3；

图4a、图4b分别是本发明实施例输入的二维水平地层模型的反射系数与合成地震数据(信噪比为3:1)；

图5a、图5b分别是以图4b中二维水平地层模型合成地震记录(信噪比为3:1)为输入得到的基追踪反演结果及对应误差；

图5c、图5d分别是以图4b中二维水平地层模型合成地震记录(信噪比为3:1)为输入得到的稀疏贝叶斯学习反演结果及对应误差；

图5e、图5f分别是以图4b中二维水平地层模型合成地震记录(信噪比为3:1)为输入得到的基于MMV模型的频率域多道反射系数联合反演结果及对应误差(L＝3)；

图6a、图6b分别是本发明实施例输入的二维起伏地层模型的反射系数及合成地震数据(信噪比为3:1)；

图7a、图7b分别是以图6b中二维起伏地层模型合成地震记录(信噪比为3:1)为输入得到的基追踪反演结果及对应误差；

图7c、图7d分别是以图6b中二维起伏地层模型合成地震记录(信噪比为3:1)为输入得到的稀疏贝叶斯学习反演结果及对应误差；

图7e、图7f分别是以图6b中二维起伏地层模型合成地震记录(信噪比为3:1)为输入得到的基于MMV模型的频率域多道反射系数联合反演结果及对应误差(L＝1)；

图8是二维楔状地层模型、水平地层模型及起伏地震模型的反演结果的均方误差(MSE)统计；

图9是本发明实施例提供的一种基于MWV模型的频率域多道反射系数联合反演系统的结构示意图；

图10为图9所示系统中反演单元70的组成结构示意图。

具体实施方式

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

针对传统频率域地震反射系数反演所存在的抗噪性差、求解精度不高以及未考虑横向连续性等问题，本发明的目的是提供一种基于MMV模型的频率域多道反射系数联合反演方法。本发明是在研究了以下问题的基础之上提出的：1)当前频率域地震反射系数反演均未考虑反射系数的横向连续性，反演结果横向稳定性通常较差；2)在频率域反射系数反演中通过MMV模型量化引入地震邻道横向连续性并将其作为约束条件加入反演流程，不仅能够有效提高反演结果的横向稳定性，而且有益于反演结果整体精度的提升；3)在含噪情况下，传统理论求解目标函数通常不能有效收敛于其全局最优解，结果稀疏性不足反演精度较低；4)在基于MMV模型的频率域多道反射系数联合反演方法中，通过联合MMV模型和块稀疏贝叶斯学习理论，实现了相关参数学习准则的优化，提高了反演结果的稀疏性和精确性。在输入单道地震资料情况下，其反演精度依然优于传统反演结果。

图1是本发明实施例提供的基于MMV模型的频率域多道反射系数联合反演方法的大致流程图，具体包括以下步骤：

步骤S1、获取经过预处理的叠后地震资料。

步骤S2、根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息。

步骤S3、从所述叠后地震资料中提取地震子波。

步骤S4、基于傅里叶变换提取所述地震子波的频谱信息。

步骤S5、根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵。

步骤S6、通过傅里叶变换逐道提取所述叠后地震资料的频谱信息。

步骤S7、基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

在一实施例中，步骤S1中所获取的经过预处理的叠后地震资料是指经过一系列地震资料处理环节的叠后地震资料，这些处理环节包括但不限于：静校正、去噪、振幅补偿、速度分析、动校正、偏移等。前期各环节的处理效果对后续频率域多道反射系数联合反演具有较大影响，因此高质量高保真的前期数据处理是确保后续反演质量的必要条件。最终处理的叠后地震资料等同于自激自收条件下的采集数据，可记为s(t,x)，其中，t表示地震波双程旅行时，即地震道时间序列，x为对应的地震道号。

获取经过预处理的叠后地震信息后对其开展分析，确定多个重要的反演参数，例如，通过频率域分析确定有效频谱范围f_span、频率采样间隔df、频率域采样点数N等参数；通过时间域分析可确定诸如时间采样间隔dt、时间采样点数M等参数。最终联合叠后地震资料时频特征共同确定算法收敛阈值tol、最大迭代次数MaxIter、噪声方差λ、相关道数L等重要参数。

在一实施例中，当利用步骤S3中基于叠后地震资料开展地震子波提取时，在步骤S1中获得的叠后地震资料中包含测井资料的条件下，可通过该测井资料中的测井曲线及井旁地震道进行地震子波提取，所提取的地震子波应符合测井资料和地震数据约束；当所述叠后地震资料不包含测井资料时，可在一定假设条件下通过地震资料开展统计性子波提取，子波提取的精度取决于实际子波与相关假设条件的匹配程度。提取的地震子波可记为w(t)，其中t为子波时间序列。

利用步骤S4基于傅里叶变换提取地震子波的频谱信息时，频谱信息可以通过复数表征，复数的模对应于地震子波的频率-能量分布规律。地震子波频谱可记为F_w(f)，其中f表示频率，单位为Hz。地震子波w(t)与其频谱F_w(f)的响应关系为：

F_w(f)＝∫w(t)e^-j2πftdt (1)

其中，F_w(f)代表复数频谱，f表示频率，单位为Hz；j为虚数，满足j×j＝-1；复频谱的模对应于地震子波的频率-能量分布规律。在实际生产中，地震子波频谱信息F_w(f)通常由离散傅里叶变换关系计算获得：

其中，M表示地震子波时间序列长度，k为离散频谱F_w(k)的对应下标。

在一实施例中，利用步骤S5构建频率域多道反射系数联合反演的正演算子矩阵时，正演算子矩阵的列数取决于上述地震资料的时间序列长度t及时间采样密度dt，行数取决于叠后地震资料的有效谱范围f_span以及频率采样密度df。正演算子矩阵可记为D(f,t)，其中，t表示时间，其中f表示频率，单位为Hz。

设叠后地震数据频率域采样点数为N，时间域采样点数为M，则正演算子矩阵D(f,t)的规模为N×M。

时间域地震褶积模型可表示为：

其中，s(t)、w(t)和r(t)分别表示时间域地震道、地震子波以及反射系数序列，

代表褶积运算。对应于频率域，地震资料频谱F_s(f)、地震子波频谱F_w(f)以及反射系数频谱F_r(f)三者响应关系可表示为：

F_s(f)＝F_w(f)F_r(f) (4)

在上式基础上，可通过地震道与地震子波的频谱计算对应的反射系数频谱：

而反射系数频谱F_r(f)与反射系数序列r(t)的响应关系亦可表征为：

F_r(f)＝D(f,t)r(t) (6)

其中，D(f,t)为正演算子矩阵，写为：

其中，t₁,t₂,...,t_M为反射系数对应的时间序列，长度为M；f₁,f₂,...,f_N为有效频谱内的频率采样成分，长度为N。

在一实施例中，在利用步骤S6获取叠后地震资料s(t,x)的频谱信息时，通过傅里叶变换逐道提取上述叠后地震资料的频谱信息，频谱信息可通过复数表征，复数的模对应于所选地震道的频率-能量分布规律。地震资料频谱可记为F_s(f)，其中f表示频率，单位为Hz。

设s(t)＝s(t,x_m)表示第m个叠后地震道，满足1≤m≤M。s(t)与其频谱F_s(f)的响应关系为：

F_s(f)＝∫s(t)e^-j2πftdt (8)

形如式(2)，通过离散傅里叶变换计算s(t)的频谱：

其中，M表示s(t)的序列长度，k为离散频谱F_s(k)的对应下标。

传统反演方法通过求解式(6)获得单道反射系数信息，没有考虑地震邻道之间的横向相关性特征，反演结果的精度和稳定性通常较差。本发明引入MMV模型量化表征地震邻道横向相关性，并结合块稀疏贝叶斯学习理论和最大期望算法优化相关参数学习准则，在算法中联合引入垂向稀疏性和横向连续性双重约束，从而综合提高反演结果质量。

在一实施例中，利用步骤S7进行反演时，通过上述地震子波频谱信息F_w(f)以及所选地震道的有效频谱信息F_s(f)，基于频谱响应关系F_s(f)＝F_w(f)F_r(f)计算反射系数频谱F_r(f)，并在此基础上结合正演算子矩阵构建频率域多道反射系数联合反演的目标函；最后，通过MMV模型量化地震资料横向连续性，并联合块稀疏贝叶斯学习理论迭代求解目标函数，以获取反射系数模型参数。

MMV模型对具有相同稀疏结构的多个信号进行观测，综合考虑多个信息的重构问题。在地震勘探中，由于相邻地震道之间的波形及其频谱特征通常具有一定的相似性，因此其所对应的反射系数也存在一定的横向连续性，即相邻道之间非零反射系数位置基本相同，符合MMV模型应用的基本假设条件。通过MMV模型可构建多道反射系数频谱与其反射系数时间序列的响应关系：

F_{r_matrix}(f)＝D(f,t)r_matrix(t) (10)

其中，F_{r_matrix}(f)和r_matrix(t)分别为多道反射系数频谱矩阵和时间域反射系数矩阵，满足F_{r_matrix}(f)＝[F_r1(f),F_r2(f),...,F_rL(f)]和r_matrix(t)＝[r₁(t),r₂(t),...,r_L(t)]。F_ri(f)、r_i(t)分别为第i道地震数据对应的反射系数频谱与反射系数时间序列。

在上式基础上，将MMV模型构建的多地震道频谱信息与多道反射系数序列的对应关系转换为单道表征关系：

F_{r_mmv}(f)＝D_mmv(f,t)r_mmv(t) (11)

其中，F_{r_mmv}(f)＝vec(F_{r_matrix}(f))，r_mmv(t)＝vec(r_matrix(t))，序列长度分别为L_{f_mmv}＝N×L和L_{t_mmv}＝M×L，L为多道地震道数且满足1≤L≤M；vec(.)表示矩阵按列所组成的向量；D_mmv(f,t)为对应算子矩阵，满足

表示克罗内克积。

对于反射系数时间序列r(t)的垂向稀疏性，稀疏贝叶斯学习理论假设任意时间位置处反射系数均满足零均值且方差为ε_i(i＝1,...,M)的高斯分布：

p(r(t_i)|ε_i)＝N(0,ε_i) (12)

当ε_i＝0时，反射系数r(t_i)＝0。因此对于稀疏信号，通过参数训练使大部分ε_i(i＝1,...,M)等于零，从而满足稀疏性约束。但上式未考虑反射系数的横向连续性特征。

为此，对于多道反射系数序列r_mmv(t)，联合MMV模型和块稀疏贝叶斯学习理论将多道反射系数的横向连续性特征在算法中量化表征，设r_mmv(t)亦满足零均值高斯分布：

p(r_mmv|ε,Γ)＝N(0,Σ_r) (13)

其中，Γ表示多道反射系数的协方差矩阵，矩阵规模为L×L；Σ_r为高斯分布的协方差矩阵，矩阵规模为L_{t_mmv}×L_{t_mmv}且满足：

其中，Ψ＝diag(ε₁,ε₂,...,ε_M)。

在此先验信息基础上，设频率域多道反射系数联合反演误差亦满足零均值高斯分布，则其似然函数可表示为：

其中，λ为噪声方差，I(L_{f_mmv})表示行\列数为L_{f_mmv}的单位矩阵。

联合式(13)和式(15)可知，频率域多道反射系数联合反演中多道反射系数序列r_mmv(t)的后验概率分布可写为：

其中，p(F_{r_mmv}|ε,Γ,λ)为正则化项；μ和Σ分别为高斯分布的均值向量与协方差矩阵，满足：

由上式可知，频率域多道反射系数联合反演的目标参数μ是地震频谱F_{r_mmv}、正则参数ε＝[ε₁,...,ε_M]、多道协方差矩阵Γ以及噪声方差λ的函数。因此，在频谱信息F_{r_mmv}已知的情况下，计算上述三类参数ε、Γ和λ信息即可获取反射系数模型参数。

设上述待求参数集合为Υ＝[ε,Γ,λ]＝[ε₁,...,ε_M,Γ,λ]，当Υ等于真实参数时，F_{r_mmv}出现的概率应为最大，即p(F_{r_mmv}|Υ)最大。因此，可认为令p(F_{r_mmv}|Υ)取得最大概率值的Υ即为最优参数组合。p(F_{r_mmv}|Υ)可表示为：

其中，V＝λI+D_mmvΣ_rD_mmv ^T。

设L(Υ)＝-logp(F_{r_mmv}|Υ)，上式中概率最大化问题即可转换为L(Υ)的极小化问题。L(Υ)为：

由于式中未包含变量r_mmv，因此不能直接通过对L(Υ)分别求偏导数获取最优参数组合Υ＝[ε₁,...,ε_M,Γ,λ]。为此，引入最大期望算法开展上述参数估计。本发明中，最大期望算法以r_mmv为隐变量，通过使得如下目标函数最大化获取最优参数组合：

其中，E(.)表示期望；Υ^(old)为前一次迭代所估计的参数组合。

上式中，Q(Υ)的第一项仅与λ有关而第二项仅与ε,Γ有关。因此上式可写为：

Q(Υ)＝Q(λ)+Q(ε,Γ) (21)

在此基础上，可通过Q(ε,Γ)对ε,Γ分别求导，Q(λ)对λ求导，获取最优参数组合Υ＝[ε₁,...,ε_M,Γ,λ]。其中Q(ε,Γ)对ε的偏导数为：

式中，Tr(.)表示矩阵的迹；μⁱ表示μ的第(i-1)L+至i×L个元素所组成的向量；Σⁱ表示矩阵Σ第(i-1)*L+1至i×L行列重叠元素所组成的方阵。由上式可知，正则参数ε_i的更新关系可写为：

此外，计算Q(ε,Γ)对Γ的求导：

由此推导出Γ的更新关系：

对于Q(λ)，可近似简化为：

由此推导出λ新的更新关系为：

通过式(23)、式(25)、式(27)所示学习准则迭代计算参数Υ＝[ε₁,...,ε_M,Γ,λ]，然后带入式(17)可获取频率域多道反射系数联合反演的多道反射系数序列r_mmv。

然而，在实际生产中由于D_mmv的矩阵规模较大，为L_{f_mmv}×L_{t_mmv}，极大制约了算法运行效率和稳定性。因此，可以利用如下近似以提高参数学习准则的运行效率：

当λ＝0或Γ为单位矩阵时，上式左右两侧完全相等。

在此基础上，可对ε_i的学习准则进行简化，令：

其中，D_i表示矩阵D第i列所组成的向量。对于μ在此基础上也可进行简化：

令

则ε_i的学习准则可简化为：

其中，X_i为矩阵X＝ΨD^T(λI+DΨD^T)^-1F_{r_mmv}第i列元素所组成的向量。

为简化Γ的学习准则，可令：

基于定点学习准则(fixed-point learning rule)，Γ的更新式可写为：

其中，||.||_f表示矩阵的弗罗贝尼乌斯范数；η为大于零的常数，调节反演稳定性：当噪声水平较低时，可令η≈0；随着噪声强度的增强，η适当增大。

对于λ的学习准则亦可简化为：

由上述学习准则(式(31、33、34))可快速获得最优参数组合Υ＝[ε₁,...,ε_M,Γ,λ]，并由此获取频率域多道反射系数联合反演的反射系数模型参数。

总体算法流程可简单归纳为：

(1)输入多道频谱信息F_{f_mmv}，子波频谱信息F_w以及相关先验参数f_span、df、λ、MaxIter、tol等；

(2)通过式(31)、式(33)、式(34)的学习准则计算最优参数组合Υ＝[ε₁,...,ε_M,Γ,λ]；

(3)计算当前多道反射系数模型参数μ^k；

(4)计算当前多道反射系数模型参数μ^k与上一次迭代模型参数μ^k-1的最大绝对误差err^k＝max(|μ^k-μ^k-1|)；

(5)判断是否满足err^k≤tol或k≥MaxIter；

若满足，输入下一组数据，直至结束；若不满足，重复步骤(2)～步骤(5)重新迭代。

在一实施例中，可以采取以下工作步骤来实现上述技术方案：1)采集叠后地震资料；2)根据所述的叠后地震资料获取频率域多道反射系数联合反演参数先验信息；3)根据所述的叠后地震资料开展地震子波提取；4)根据所述的地震子波开展傅里叶变换获取地震子波频谱信息；5)构建频率域多道反射系数联合反演的正演算子矩阵；6)根据所述的叠后地震资料开展傅里叶变换获取地震资料频谱信息；7)根据所述的频谱信息构建频率域多道反射系数联合反演目标函数，结合MMV模型和块稀疏贝叶斯学习理论求解目标函数获取反射系数模型参数；8)对每一组地震道组合重复步骤7)直至反演残差小于给定阈值或迭代次数大于给定的最大迭代次数；9)对所有地震道组合重复步骤6)～8)直至结束；10)最终输出反射系数模型参数。

图2a是本发明一实施例中提供的输入的二维楔状地层模型的反射系数，图2b是本发明实施例中输入的二维楔状地层模型的反射系数及其频谱信息。图2c、图2d分别是本发明实施例中输入的二维楔状地层模型在30Hz处的零相位雷克子波及其频谱信息。图2e、图2f分别是本发明实施例中输入的二维楔状地层模型的合成地震记录及其频谱信息。

图3a～图3f分别是以图2e中二维楔状地层模型的合成地震记录为输入，在不含噪情况下，分别基于基追踪理论、稀疏贝叶斯学习理论以及基于MMV模型的频率域多道反射系数联合方法(相关道数L＝3)开展频率域反射系数反演得到的反演结果及对应的误差。黑色虚线表示实际非零反射系数所在位置，如反演结果所示，基于MMV模型的频率域多道反射系数联合反演结果(图3e、图3f)的薄层反射系数识别能力明显优于传统基追踪理论(图3a、图3b)和稀疏贝叶斯学习理论(图3c、图3d)，体现了本发明的理论优势。

图4a、图4b分别是本发明实施例输入的二维水平地层模型的反射系数与合成地震数据，其中信噪比为3:1，地震子波为30Hz零相位雷克子波，如图所示，反射系数振幅值向右向下逐渐衰减。图5a～图5f分别是以图4b中二维水平地层模型合成地震记录(信噪比为3:1)为输入，在SNR＝3:1情况下，分别基于基追踪理论、稀疏贝叶斯学习理论以及基于MMV模型的频率域多道反射系数联合方法(L＝3)开展频率域反射系数反演得到的反演结果及对应误差。如反演结果所示，在数据含噪情况下，基于MMV模型的频率域多道反射系数联合方法的反演结果(图5e、图5f)精度优于传统的基追踪理论(图5a、图5b)和稀疏贝叶斯学习理论(图5c、图5d)，体现了本发明理论优势。

图6a、图6b分别是本发明实施例输入的二维起伏地层模型的反射系数及合成地震数据，信噪比为3:1，地震子波为30Hz零相位雷克子波。图7a～图7f分别是以图4b中二维起伏地层模型合成地震记录为输入得到的反演结果及对应误差，在SNR＝3:1情况下，分别基于基追踪理论、稀疏贝叶斯学习理论以及基于MMV模型的频率域多道反射系数联合方法(L＝1)开展频率域反射系数反演得到的反演结果及对应误差。如反演结果所示，基于MMV模型的频率域多道反射系数联合反演在未考虑反射系数横向连续性的情况下，反演结果(图7e、图7f)精度依然优于传统基追踪理论(图7a、图7b)和稀疏贝叶斯学习理论(图7c、图7d)，体现了本发明的理论优势。

图8是二维楔状地层模型、水平地层模型及起伏地震模型的反演结果的均方误差(MSE)统计。如图8所示，对所有实施例输入的三个地层模型的反演结果开展均方误差(MSE)统计分析。统计结果表明：1)基于MMV模型的频率域多道反射系数联合反演方法的反演精度高于传统理论的反演结果，提高了反演精度；2)MMV模型的引入提高了频率域反射系数反演结果的横向连续性和稳定性；3)即使在不考虑地震横向相关性的情况下，基于MMV模型的频率域多道反射系数联合反演结果精度依然优于传统理论的反演精度。

基于与图1所示基于MWV模型的频率域多道反射系数联合反演方法相同的发明构思，本申请实施例还提供了一种系统，如下面实施例所述。由于该系统解决问题的原理与图1中基于MWV模型的频率域多道反射系数联合反演方法相似，因此该系统的实施可以参见图1的基于MWV模型的频率域多道反射系数联合反演方法的实施，重复之处不再赘述。

在另一实施例中，本发明还提供了一种基于MWV模型的频率域多道反射系数联合反演系统，其结构如图9所示，该系统包括：叠后地震资料获取单元10、信息获取单元20、地震子波提取单元30、第一频谱信息获取单元40、正演算子矩阵构建单元50、第二频谱信息获取单元60及反演单元70。

其中，叠后地震资料获取单元10用于获取经过预处理的叠后地震资料；信息获取单元20用于根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；地震子波提取单元30用于从所述叠后地震资料中提取地震子波；第一频谱信息获取单元40用于基于傅里叶变换提取所述地震子波的频谱信息；正演算子矩阵构建单元50用于根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；第二频谱信息获取单元60用于通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；反演单元70用于基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

在一实施例中，所述预处理包括静校正、去噪、振幅补偿、速度分析、动校正、偏移中的至少一种。

在一实施例中，所述频率域多道反射系数联合反演参数先验信息至少包括：有效频谱范围、时间采样密度、频率采样密度、噪声方差、收敛阈值及最大迭代次数。

在一实施例中，如所述叠后地震资料中包含测井资料，则地震子波提取单元30基于所述测井资料中的测井曲线及井旁地震道进行地震子波提取；如所述叠后地震资料中不包含测井资料，则地震子波提取单元30对所述叠后地震资料进行统计性子波提取以获取地震子波。

在一实施例中，所述正演算子矩阵构建单元50根据所述地震资料的时间序列长度及时间采样密度确定所述正演算子矩阵的列数，根据所述地震资料的有效频谱范围及频率采样密度确定所述正演算子矩阵的行数。

在一实施例中，反演单元70的结构如图10所示，其主要包括：反射系数频谱获取模块71用于根据所述地震子波的频谱信息以及所选地震道的频谱信息，基于频谱响应关系计算反射系数频谱；目标函数构建模块72用于结合所述正演算子矩阵构建频率域多道反射系数联合反演的目标函数；反演计算模块73用于通过MMV模型量化叠后地震资料的横向连续性，联合块稀疏贝叶斯学习理论迭代求解所述目标函数，以获取反射系数模型参数。

本发明由于采取以上技术方案，其具有以下优点：1)本技术方案所述的基于MMV模型的频率域多道反射系数联合反演方法，与当前现有技术相比，方案通过MMV模型量化引入反射系数的横向连续性特征，大幅度提高了反演结果的横向稳定性与反演精度；2)本技术方案通过块稀疏贝叶斯学习理论，并借助最大期望算法构建参数学习准则，实现了包括噪声方差、正则参数、协方差矩阵等多参数联合训练，提高了参数训练的精度；3)本技术方案求解频率域多道反射系数反演目标函数，反演结果较传统理论更接近全局最优解，也保证了反演结果的稀疏性；4)本技术方案的抗噪性明显优于传统理论；5)本技术方案同时考虑了反射系数的垂向稀疏性和横向连续性，并在算法中联合加入上述约束条件，能够有效提高频率域反射系数反演的精度和空间连续性。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

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

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

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

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

Claims (12)

1.一种基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，所述方法包括：

步骤S1、获取经过预处理的叠后地震资料；

步骤S2、根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；

步骤S3、从所述叠后地震资料中提取地震子波；

步骤S4、基于傅里叶变换提取所述地震子波的频谱信息；

步骤S5、根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；

步骤S6、通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；

步骤S7、基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

2.根据权利要求1所述的基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，所述预处理包括静校正、去噪、振幅补偿、速度分析、动校正、偏移中的至少一种。

3.根据权利要求1所述的基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，所述频率域多道反射系数联合反演参数先验信息至少包括：有效频谱范围、时间采样密度、频率采样密度、噪声方差、收敛阈值及最大迭代次数。

4.根据权利要求1所述的基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，在步骤S3中，如所述叠后地震资料中包含测井资料，则基于所述测井资料中的测井曲线及井旁地震道进行地震子波提取；如所述叠后地震资料中不包含测井资料，则对所述叠后地震资料进行统计性子波提取以获取地震子波。

5.根据权利要求1所述的基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，在步骤S5中，根据所述叠后地震资料的时间序列长度及时间采样密度确定所述正演算子矩阵的列数，根据所述叠后地震资料的有效频谱范围及频率采样密度确定所述正演算子矩阵的行数。

6.根据权利要求1所述的基于MWV模型的频率域多道反射系数联合反演方法，其特征在于，所述步骤S7包括：

根据所述地震子波的频谱信息以及所选地震道的频谱信息，基于频谱响应关系计算反射系数频谱；

结合所述正演算子矩阵构建频率域多道反射系数联合反演的目标函数；

通过MMV模型量化叠后地震资料的横向连续性，联合块稀疏贝叶斯学习理论迭代求解所述目标函数，以获取反射系数模型参数。

7.一种基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，所述系统包括：

叠后地震资料获取单元，用于获取经过预处理的叠后地震资料；

信息获取单元，用于根据所述叠后地震资料获取频率域多道反射系数联合反演参数先验信息；

地震子波提取单元，用于从所述叠后地震资料中提取地震子波；

第一频谱信息获取单元，用于基于傅里叶变换提取所述地震子波的频谱信息；

正演算子矩阵构建单元，用于根据所述频率域多道反射系数联合反演参数先验信息构建频率域多道反射系数联合反演的正演算子矩阵；

第二频谱信息获取单元，用于通过傅里叶变换逐道提取所述叠后地震资料的频谱信息；

反演单元，用于基于所述频谱信息构建频率域多道反射系数联合反演的目标函数，并通过MMV模型和块稀疏贝叶斯学习理论求解所述目标函数，以获取反射系数模型参数。

8.根据权利要求7所述的基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，所述预处理包括静校正、去噪、振幅补偿、速度分析、动校正、偏移中的至少一种。

9.根据权利要求7所述的基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，所述频率域多道反射系数联合反演参数先验信息至少包括：有效频谱范围、时间采样密度、频率采样密度、噪声方差、收敛阈值及最大迭代次数。

10.根据权利要求7所述的基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，如所述叠后地震资料中包含测井资料，则所述地震子波提取单元基于所述测井资料中的测井曲线及井旁地震道进行地震子波提取；如所述叠后地震资料中不包含测井资料，则所述地震子波提取单元对所述叠后地震资料进行统计性子波提取以获取地震子波。

11.根据权利要求7所述的基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，所述正演算子矩阵构建单元根据所述地震资料的时间序列长度及时间采样密度确定所述正演算子矩阵的列数，根据所述地震资料的有效频谱范围及频率采样密度确定所述正演算子矩阵的行数。

12.根据权利要求7所述的基于MWV模型的频率域多道反射系数联合反演系统，其特征在于，所述反演单元包括：

反射系数频谱获取模块，用于根据所述地震子波的频谱信息以及所选地震道的频谱信息，基于频谱响应关系计算反射系数频谱；

目标函数构建模块，用于结合所述正演算子矩阵构建频率域多道反射系数联合反演的目标函数；

反演计算模块，用于通过MMV模型量化叠后地震资料的横向连续性，联合块稀疏贝叶斯学习理论迭代求解所述目标函数，以获取反射系数模型参数。

Images