CN103293552A - 一种叠前地震资料的反演方法及系统 - Google Patents

Abstract

本发明提供一种叠前地震资料的反演方法及系统，所述的方法包括：采集当前地层的叠前地震资料；从所述的叠前地震资料中选取优势频带范围；从所述的优势频带范围中提取角道集地震数据以及地震子波；根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数；根据相关向量机以及所述的角道集反射系数确定超参数；根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。提高了反演结果的分辨率和准确性，提高了识别薄层、薄互层岩性油气藏的地层关系和性质的准确性。

Description

一种叠前地震资料的反演方法及系统

技术领域

本发明关于石油地球物理勘探技术领域，特别是关于石油地球物理勘探中的叠前反演技术，具体的讲是一种叠前地震资料的反演方法及系统。

背景技术

在地震勘探技术领域中，振幅随偏移距的变化（Amplitude Versus Offset，AVO）经常用于地层三参数（纵波速度、横波速度和密度）反射系数的反演，进而获取用于识别含油气地层的纵波速度、横波速度以及密度等基本岩性信息。

现有技术中的基于最小二乘法的线性AVO反演方法，主要是基于纵波反射系数同三参数反射系数的近似关系，联合褶积模型，建立在时间域角道集地震数据同地层三参数反射系数的关系，通过最小二乘法进行求解，进而获得纵波速度、横波速度以及密度三参数反射系数，然后通过道积分获得纵波速度、横波速度以及密度的相对值。在此基础上，可以结合测井数据，通过地质统计的办法，获得纵波速度、横波速度以及密度的绝对值描述地层岩性。上述基于最小二乘法的线性AVO反演方法主要存在如下问题：

1、最小二乘法线性反演本质上无法拓宽地震资料的频带，受地震资料有限带宽的影响，反演结果的频带宽度基本等同于原始地震资料的频带；

2、时间域反演等价于利用了地震资料的全部频率信息。在实际中，受采集条件和噪声的影响，以陆上资料为例，叠前地震资料小于5Hz的低频和高于150Hz的高频成分通常是错误的冗余信息，并不能代表地下真实的地震响应。因此，这部分错误的冗余信息参与反演，会影响反演结果的准确性。而滤除这部分信息，会由于带通滤波窗函数的截断效应而带入新的数据误差。

在地震勘探中，叠前地震资料保留了丰富的波场信息，如何利用这些信息来获取地下地层参数为油气检测和识别提供数据支撑是地震勘探最为重要的工作之一。目前地震勘探中主要以一次反射波的动力学理论为基础，在层状地层假设下实施，诸如上述的AVO技术。然而，受到地震采集技术限制和噪声的污染，地震资料的有效带宽是有限的，低频和高频成分往往是错误的信息，传统的时间域线性AVO反演无法拓宽原始地震资料的频带，反演结果的准确性也较低。

因此，如何根据有限带宽的叠前地震资料来提高AVO解释的精度，进而提高反演结果的准确性是本领域亟待解决的技术难题。

发明内容

为了克服现有技术存在的上述问题，本发明提供了一种叠前地震资料的反演方法及系统，以叠前地震资料中的优势频带范围为数据，在贝叶斯理论框架下，利用相关向量机求解稀疏脉冲反射系数，提高了反演结果的分辨率和准确性，提高了识别薄层、薄互层岩性油气藏的地层关系和性质的准确性。

本发明的目的之一是，提供一种叠前地震资料的反演方法，包括：采集当前地层的叠前地震资料；从所述的叠前地震资料中选取优势频带范围；从所述的优势频带范围中提取角道集地震数据以及地震子波；根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数；根据相关向量机以及所述的角道集反射系数确定超参数；根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

本发明的目的之一是，提供了一种叠前地震资料的反演系统，包括：地震资料采集装置，用于采集当前地层的叠前地震资料；频带范围选取装置，用于从所述的叠前地震资料中选取优势频带范围；角道集数据和子波提取装置，用于从所述的优势频带范围中提取角道集地震数据以及地震子波；角道集反射系数确定装置，用于根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数；超参数确定装置，用于根据相关向量机以及所述的角道集反射系数确定超参数；超级分辨率反射系数确定装置，用于根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；解释装置，用于对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

本发明的有益效果在于，通过以叠前地震资料中的优势频带范围为数据，利用贝叶斯理论框架，考虑了待反演参数的概率分布特征，提高了反演结果的稳定性，采用相关向量机算法得到稀疏的脉冲反射系数，补偿原始记录中没有采集到的频率成分，具有超级分辨率的效果。

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

附图说明

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

图1为本发明实施例提供的一种叠前地震资料反演方法的流程图；

图2为本发明实施例提供的一种叠前地震资料反演方法的实施方式二的流程图；

图3为本发明实施例提供的一种叠前地震资料反演系统的结构框图；

图4为本发明实施例提供的一种叠前地震资料反演系统的实施方式二的结构框图；

图5为纵波入射下的反射纵波和横波以及透射纵波和横波的示意图；

图6为层状薄层模型对应的三参数（纵波速度V_P、横波速度V_S和密度ρ）数据模型图；

图7为层状薄层模型对应的叠前角道集一次反射波反射系数和地震记录图；

图8为层状薄层模型对应的三参数反射系数数据的反演结果（纵波反射系数r_p、横波反射系数r_s和密度反射系数r_ρ）与三参数模型数据的对比图；

图9为推覆体二维模型的三参数数据（纵波速度V_P、横波速度V_S和密度ρ）模型图；

图10为推覆体二维模型对应的三参数反射系数（纵波反射系数r_P、横波反射系数r_S和密度反射系数r_ρ）模型图；

图11为推覆体二维模型的反演结果图（纵波反射系数r_p_inv、横波反射系数r_s_inv和密度反射系数r_ρ_inv）。

具体实施方式

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

图1为本发明实施例提供的一种叠前地震资料反演方法的流程图，由图1可知，该方法具体包括：

S101：采集当前地层的叠前地震资料。

在地震勘探中，叠前地震资料中保留了丰富的波场信息，如何利用这些信息来获取地下地层参数为油气检测和识别提供数据支撑是地震勘探最为重要的工作之一。目前地震勘探中主要以一次反射波的动力学理论为基础，在层状地层假设下实施的，比如AVO技术。

S102：从所述的叠前地震资料中选取优势频带范围。

叠前地震资料由于受到地震采集技术限制和噪声的污染，地震资料的有效带宽是有限的，低频和高频成分往往是错误的信息。传统的时间域线性AVO反演无法提高原始地震资料的分辨率，反演结果的准确性也较低。因此，本发明的步骤S102从叠前地震资料中选取出优势频带范围，从有限带宽的叠前地震资料中选取出部分谱（不包括错误的低频和高频的冗余信息）。在具体的实施方式中，选取的优势频带范围应包含峰值频率附近的频点资料。

因此，如何从有限带宽的叠前地震资料的部分谱信息（不包括错误的低频和高频的冗余信息）中获取宽频带高分辨率的三参数（纵波速度、横波速度和密度）反射系数，补偿原始地震资料中没有的频率成分，是提高AVO解释精度的关键。为了便于进一步理解本发明的技术方案，下面首先介绍平面波弹性动力学AVO理论及模拟。

根据弹性动力学理论，弹性半空间水平分界面上，当一平面纵波震源入射情况下，会产生反射纵波和横波及透射纵波和横波，如图5所示。因此，地震勘探中，叠前地震资料保留了地震波振幅和相位随偏移距变化的动力学信息。在远场情况下，仅仅考虑平面地震波振幅随偏移距变化的动力学信息，又称为地震波的AVO特性。它描述了弹性半空间水平分界面上反射系数和透射系数取决于偏移距对应的入射角和两层介质中的性质（包括纵波速度、横波速度以及密度）。此时，地震波射线路径满足Snell（斯奈尔）定理，不同类型的反射和透射波振幅（反射、透射系数）特性满足Zoeppritz（佐普利兹）关系式，其矩阵形式如下所示：

$\left[\begin{array}{cccc}\stackrel{\↓\↑}{\mathrm{PP}}& \stackrel{\↓\↑}{\mathrm{SP}}& \stackrel{\↑\↑}{\mathrm{PP}}& \stackrel{\↑\↑}{\mathrm{SP}}\\ \stackrel{\↓\↑}{\mathrm{PS}}& \stackrel{\↓\↑}{\mathrm{SS}}& \stackrel{\↑\↑}{\mathrm{PS}}& \stackrel{\↑\↑}{\mathrm{SS}}\\ \stackrel{\↓\↓}{\mathrm{PP}}& \stackrel{\↓\↓}{\mathrm{SP}}& \stackrel{\↑\↓}{\mathrm{PP}}& \stackrel{\↑\↓}{\mathrm{SP}}\\ \stackrel{\↓\↓}{\mathrm{PS}}& \stackrel{\↓\↓}{\mathrm{SS}}& \stackrel{\↑\↓}{\mathrm{PS}}& \stackrel{\↑\↓}{\mathrm{SS}}\end{array}\right]=M^{-1}N---\left(1\right)$

其中，每个矩阵元素代表一种类型的反射、透射系数。第一个字母代表入射波类型，第二个字母代表反射波类型，P代表纵波，S代表横波，箭头↑表示上行传播，↓表示下行传播，其中M和N的矩阵形式如下

$M=\left[\begin{array}{cccc}-\sin {\mathrm{\θ}}_1& -\cos {\mathrm{\φ}}_1& {\sin \mathrm{\θ}}_2& {\cos \mathrm{\φ}}_2\\ \cos {\mathrm{\θ}}_1& -\sin {\mathrm{\φ}}_1& \cos {\mathrm{\θ}}_2& -\sin {\mathrm{\φ}}_2\\ 2{\mathrm{\ρ}}_1{V}_{S1}\sin {\mathrm{\φ}}_1\cos {\mathrm{\θ}}_1& {\mathrm{\ρ}}_1{V}_{S1}(1-2{\sin }^2{\mathrm{\φ}}_1)& 2{\mathrm{\ρ}}_2{V}_{S2}\sin {\mathrm{\φ}}_2\cos {\mathrm{\θ}}_2& {\mathrm{\ρ}}_2{V}_{S2}(1-2{\sin }^2{\mathrm{\φ}}_2)\\ -{\mathrm{\ρ}}_1{V}_{P1}(1-2{\sin }^2{\mathrm{\φ}}_1)& {\mathrm{\ρ}}_1{V}_{S1}\sin 2{\mathrm{\φ}}_1& {\mathrm{\ρ}}_2{V}_{P2}(1-2{\sin }^2{\mathrm{\φ}}_2)& -{\mathrm{\ρ}}_2{V}_{S2}\sin 2{\mathrm{\φ}}_2\end{array}\right]$

$N=\left[\begin{array}{cccc}\sin {\mathrm{\θ}}_1& \cos {\mathrm{\φ}}_1& {-\sin \mathrm{\θ}}_2& {-\cos \mathrm{\φ}}_2\\ \cos {\mathrm{\θ}}_1& -\sin {\mathrm{\φ}}_1& \cos {\mathrm{\θ}}_2& -\sin {\mathrm{\φ}}_2\\ 2{\mathrm{\ρ}}_1{V}_{S1}\sin {\mathrm{\φ}}_1\cos {\mathrm{\θ}}_1& {\mathrm{\ρ}}_1{V}_{S1}(1-2{\sin }^2{\mathrm{\φ}}_1)& 2{\mathrm{\ρ}}_2{V}_{S2}\sin {\mathrm{\φ}}_2\cos {\mathrm{\θ}}_2& {\mathrm{\ρ}}_2{V}_{S2}(1-2{\sin }^2{\mathrm{\φ}}_2)\\ {\mathrm{\ρ}}_1{V}_{P1}(1-2{\sin }^2{\mathrm{\φ}}_1)& -{\mathrm{\ρ}}_1{V}_{S1}\sin 2{\mathrm{\φ}}_1& -{\mathrm{\ρ}}_2{V}_{P2}(1-2{\sin }^2{\mathrm{\φ}}_2)& {\mathrm{\ρ}}_2{V}_{S2}\sin 2{\mathrm{\φ}}_2\end{array}\right]$

M和N包含地层分界面两侧的三（物性）参数和几何参数，这说明平面波在分界面上的反射、透射动力学特征同两侧的物性参数和几何参数（入射角、反射角以及透射角）有关。

对于公式(1)一般不直接求逆，可采用克莱母法则求解。实际应用中，一般采用Aki和Richard（阿基安亿敬一和理查德）给出的解析公式，计算更简单实用：

$R_{\mathrm{PP}}=[(b\frac{\cos {\mathrm{\θ}}_1}{V_{P1}}-c\frac{\cos {\mathrm{\θ}}_2}{V_{P2}})F-(a+d\frac{\cos {\mathrm{\θ}}_1}{V_{P1}}\frac{\cos {\mathrm{\φ}}_2}{V_{S2}}){\mathrm{Hp}}^2]/D---\left(2\right)$

$R_{\mathrm{PS}}=\left[-2\frac{\cos {\mathrm{\θ}}_1}{V_{P_1}}(\mathrm{ab}+\mathrm{cd}\frac{\cos {\mathrm{\θ}}_2}{V_{P2}}\frac{\cos {\mathrm{\φ}}_2}{V_{S2}})p{V}_{P1}\right]/V_{S1}D---\left(3\right)$

其中：a=ρ₂(1-2sin²φ₂)-ρ₁(1-2sin²φ₁)，

b=ρ₂(1-2sin²φ₂)+2ρ₁sin²φ₁，

c=ρ₁(1-2sin²φ₁)+2ρ₂sin²φ₂，

$d=2({\mathrm{\ρ}}_2{V}_{S2}^2-{\mathrm{\ρ}}_1{V}_{S1}^2),$

D=EF+GHp²=(detM)/(V_P1V_P2V_S1V_S2)，

$E=b\frac{\cos {\mathrm{\θ}}_1}{V_{P1}}+c\frac{\cos {\mathrm{\θ}}_2}{V_{P2}},$

$F=b\frac{\cos {\mathrm{\φ}}_1}{V_{S1}}+c\frac{\cos {\mathrm{\φ}}_2}{V_{S2}},$

$G=a-d\frac{\cos {\mathrm{\θ}}_1}{V_{P1}}\frac{\cos {\mathrm{\φ}}_2}{V_{S2}},$

$H=a-d\frac{\cos {\mathrm{\θ}}_2}{V_{P2}}\frac{\cos {\mathrm{\φ}}_1}{V_{S1}},$

为射线参数。

假设层间参数变化微小，Aki和Richards将公式(2)和公式(3)简化成以下近似关系：

$R_{\mathrm{pp}}\≈\frac{1}{2}(1-4p^2{V}_S^2)\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{1}{2{\cos }^2\mathrm{\θ}}\frac{\mathrm{\Δ}{V}_P}{V_P}-4p^2{V}_S^2\frac{\mathrm{\Δ}{V}_S}{V_S}---\left(4\right)$

$R_{\mathrm{PS}}\≈\frac{p{V}_P}{2\cos \mathrm{\φ}}[(1-2V_S^2{p}^2+2V_S^2\frac{\cos \mathrm{\θ}}{V_P}\frac{\cos \mathrm{\φ}}{V_S})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}-(4V_S^2{p}^2-4V_S^2\frac{\cos \mathrm{\θ}}{V_P}\frac{\cos \mathrm{\φ}}{V_S})\frac{\mathrm{\Δ}{V}_S}{V_S}---\left(5\right)$

其中，Δρ=ρ₂-ρ₁，ΔV_P=V_P2-V_P1，ΔV_S=V_S2-V_S1，θ=(θ₁+θ₂)/2，ρ=(ρ₂+ρ₁)/2，V_P=(V_P2+V_P1)/2，V_S=(V_S2+V_S1)/2。公式(2)和公式(3)是地震勘探中通常进行纵波入射AVO正演模拟的公式。实际勘探中，当近似公式(4)和公式(5)同准确公式(2)和公式(3)近似程度较高时，往往采用更为简洁的近似公式(4)和公式(5)，它表征了三参数反射系数同角道集反射系数的线性关系及各自所占的权重，实际应用时非常方便，物理含义也更明确。本发明为了从叠前地震资料中反演出超级分辨率的三参数反射系数，同时考虑到实际勘探中，主要是纵波勘探，因此，以公式(4)为基础推导本发明的实现方法。

S103：从所述的优势频带范围中提取角道集地震数据以及地震子波。假设角道集一次反射波地震数据为s(t,θ)，地震子波为w(t)，角道集反射系数为r(t,θ)。

S104：根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数。根据褶积模型，则有

s(t,θ)=w(t)*r(t,θ)     (6)

根据傅立叶变换转换为频域，则有

S(f,θ)=W(f)R(f,θ)     (7)

其中，S(f,θ)为角道集地震数据，W(f)为地震子波，R(f,θ)为角道集反射系数，f为频率，θ为角度，t为时间。

S105：根据相关向量机以及所述的角道集反射系数确定超参数；

S106：根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；

S107：对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

图2为本发明实施例提供的一种叠前地震资料的反演方法的实施方式二的流程图，由图2可知，在实施方式二中，该方法具体包括：

S201：采集当前地层的叠前地震资料；

S202：从所述的叠前地震资料中选取优势频带范围；

S203：从所述的优势频带范围中提取角道集地震数据以及地震子波；

S204：根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数。上述步骤S201至S204与实施方式一中的步骤S101至S104相同，此处不再赘述。

实施方式一中的步骤S105具体包括：

S205：根据所述的褶积模型建立角道集反射系数的反演矩阵方程。公式(7)中的角道集反射系数满足：

$r(t,\mathrm{\θ})\≈a(t,\mathrm{\θ})\frac{\mathrm{\Δ\ρ}\left(t\right)}{\mathrm{\ρ}\left(t\right)}+b(t,\mathrm{\θ})\frac{\mathrm{\Δ}{V}_P\left(t\right)}{V_P\left(t\right)}+c(t,\mathrm{\θ})\frac{\mathrm{\Δ}{V}_S\left(t\right)}{V_S\left(t\right)}---\left(8\right)$

其中，三参数的权系数为

$\left\{\begin{array}{c}a(t,\mathrm{\θ})=\frac{1}{2}(1-4p^2{V}_S^2\left(t\right))\\ b(t,\mathrm{\θ})=\frac{1}{2{\cos }^2\mathrm{\θ}}\\ c(t,\mathrm{\θ})=-4p^2{V}_S^2\left(t\right)\end{array}\right.---\left(9\right)$

实际中，由于层速度通常都是未知的，当速度变化不大时，利用平均速度代替层速度

，x=P,S，通常可以对式(9)进一步取近似关系

$\left\{\begin{array}{c}a(t,\mathrm{\θ})\≈\frac{1}{2}(1-4p^2{\stackrel{\‾}{V}}_S^2)\\ b(t,\mathrm{\θ})\≈\frac{1}{2{\cos }^2\mathrm{\θ}}\\ c(t,\mathrm{\θ})\≈-4p^2{\stackrel{\‾}{V}}_S^2\end{array}\right.---\left(10\right)$

考虑到频率域角道集反射系数满足：

$R(f,\mathrm{\θ})=\underset{t}{\mathrm{\Σ}}\exp (-i2\mathrm{\πf})r(t,\mathrm{\θ})---\left(11\right)$

由此可建立基于角道集反射系数R(f,θ)=S(f,θ)/W(f)，f∈[f₁,f_M]与三参数反射系数的矩阵关系为：

$\left[\begin{array}{c}R(f_1,\mathrm{\θ})\\ R(f_2,\mathrm{\θ})\\ .\\ .\\ .\\ R(f_M,\mathrm{\θ})\end{array}\right]=F\left[\begin{array}{ccc}A& B& C\end{array}\right]\left[\begin{array}{c}{r}_{\mathrm{\ρ}}\\ r_P\\ r_S\end{array}\right]---\left(12\right)$

其中，傅立叶变换阵 $F=\left[\begin{array}{cccc}{e}^{-i2\mathrm{\π}{f}_1{t}_1}& e^{-i2\mathrm{\π}{f}_1{t}_2}& ...& e^{-i2\mathrm{\π}{f}_1{t}_N}\\ e^{-i2\mathrm{\π}{f}_2{t}_1}& e^{-i2\mathrm{\π}{f}_2{t}_2}& ...& e^{-i2\mathrm{\π}{f}_2{t}_N}\\ .& .& ...& .\\ .& .& & .\\ .& .& & .\\ e^{-i2\mathrm{\π}{f}_M{t}_1}& e^{-i2\mathrm{\π}{f}_M{t}_2}& ...& e^{-i2\mathrm{\π}{f}_M{t}_N}\end{array}\right]$ ，三参数反射系数向量及其相应的权系数对角矩阵为：

，A=diag[a(t₁,θ)a(t₂,θ)...a(t_N,θ)]；

，B=diag[b(t₁,θ)b(t₂,θ)...b(t_N,θ)]；

，C=diag[c(t₁,θ)c(t₂,θ)...c(t_N,θ)]。

因此，令R(θ_i)=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，若选取角度道集θ_i∈[θ₁,θ_K]，K≥3可以建立三参数反射系数的反演矩阵方程为：

令

由此，可将公式(13)用矩阵符号简化为d_K*M×1=G_K*M×3Nm_3N×1，其中，m=[r_ρ r_P r_S]^T，d=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，K为选取的角度总数，M为选取的离散频率总数，N为反射系数采样总数，F为傅立叶变换算子，A为密度反射系数的权系数构成的对角矩阵，B为纵波速度反射系数的权系数构成的对角矩阵，C为横波速度反射系数的权系数构成的对角矩阵，r_ρ为密度反射系数，r_P为纵波速度反射系数，r_S为横波速度反射系数。

S206：根据似然函数、以零值为中心的标准正态分布以及所述的反演矩阵方程确定超参数。

令观测数据误差为n，则公式(13)的反演矩阵的简化方程为：

d=Gm+n=[G₁ G₂...G_k...G_3N][m₁ m₂...m_k...m_3N]^T+n  (14)

其中， $\left[\begin{array}{c}{r}_{\mathrm{\ρ}}={\left[\begin{array}{cccc}{m}_1& m_2& ...& m_N\end{array}\right]}^T\\ r_P={\left[\begin{array}{cccc}{m}_{N+1}& m_{N+2}& ...& m_{2N}\end{array}\right]}^T\\ r_S={\left[\begin{array}{cccc}{m}_{2N+1}& m_{2N+2}& ...& m_{3N}\end{array}\right]}^T\end{array}\right].$

相关向量机的核心思想是利用贝叶斯框架，在最大后验概率下顺序寻找相关向量，并对应求解权值。此处的G_k表示相关向量，m_k表示相关向量的权重。本质地，获得相关向量的过程就是求解反射系数时间位置或确定地层边界的过程，获得权值的过程实质就是求取反射系数大小的过程。

假设n服从均值为0、方差为σ²的高斯分布，即n～N(0,σ²)，在给定模型参数m和σ²的条件下，d的条件概率，也称为似然函数可表示为：

p(d|m,σ²)=(2πσ²)^-K*Mexp[-(d-Gm)^T(d-Gm)/(2σ²)]  (15)

为了获得受地质假设驱动的解，稀疏贝叶斯学习方法需要对m施加先决概率条件，根据贝叶斯理论，限制m的概率分布是以零值为中心的标准正态分布，即

$p\left(m\right|h)={\left(2\mathrm{\π}\right)}^{-3N/2}\underset{k=1}{\overset{3N}{\mathrm{\Π}}}{h}_k^{1/2}\exp ({-h}_k{m}_k^2/2)---\left(16\right)$

其中，h=[h₁,h₂,...,h_3N]^T代表3N个相互独立的超参数，每一个超参数分别控制其对应反射系数大小的先验信息。

表示反射系数m_k服从均值为0，方差为

的高斯分布，其中，k=1，2，…，3N。当

时，m_k→0，表示位置k处没有脉冲反射系数，如此最终会获得稀疏的超级分辨率脉冲反射系数序列。

为了得到反射系数，首先要确定超参数h和方差σ²的最佳值。根据贝叶斯框架，超参数的边缘分布通过下式计算

p(d|h,σ²)=(2π)^3N|Q|exp(d^TQ^-1d)    (17)

其中， $Q={\mathrm{\σ}}^{2}I+{\mathrm{GH}}^{-1}{G}^T={\mathrm{\σ}}^{2}I+\underset{k\≠j}{\mathrm{\Σ}}{h}_k^{-1}{G}_k{G}_k^T+h_j^{-1}{G}_j{G}_j^T$ ，对公式(17)求最小即可求出超参数h，此处的N为反射系数采样总数，p为条件概率。

S207：根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数。该步骤与实施方式一中的步骤S106相同。

由贝叶斯准则，在d，h和σ²已知的条件下，联立公式(15)和公式(16)，则有超级分辨率稀疏脉冲反射系数m的条件概率或后验概率密度分布：

$p\left(m\right|d,h,{\mathrm{\σ}}^2)=C{\left|\mathrm{\Σ}\right|}^{-\frac{1}{2}}\exp [-\frac{1}{2}{(m-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(m-\mathrm{\μ})]---\left(18\right)$

其中， $\left\{\begin{array}{c}\mathrm{\Σ}={(H+{\mathrm{\σ}}^{-2}{G}^{T}G)}^{-1}\\ \mathrm{\μ}={\mathrm{\σ}}^{-2}\mathrm{\Σ}{G}^{T}d\end{array}\right.---\left(19\right)$

其中，C为一常数，Σ为协方差，μ为均值，即超级分辨率稀疏脉冲反射系数，H=diag(h₁,h₂,...,h_3N)是一个对角矩阵。超级分辨率稀疏脉冲反射系数m的估计由反射系数后验分布的均值μ给出。

在具体的实施方式中，叠前超级分辨率反演方法利用相关向量机算法对公式(17)和公式(19）进行循环迭代求解。每次迭代只更新一个反射系数脉冲m_k对应的超参数h_k和一个相关向量G_k，但每次更新后目标函数或公式(17)都会减小。反复迭代，利用公式(17)计算h，公式(19)计算μ，直到前后两次迭代的目标函数之差达到容忍误差。最后，输出的μ即是估计的稀疏反射系数。实质上，通过公式(17)每次迭代求出超参数h，相当于获得稀疏脉冲的位置，而通过公式(19)求出μ，相当于获得稀疏脉冲的大小。由于相关向量机在迭代运算时，只是通过添加或删除基向量G_k来实时更新脉冲的位置和大小，不需要进行大矩阵的求逆运算，所以该方法具有很快的收敛速度。

S208：对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。该步骤与实施方式一中的步骤S107相同。超级分辨率稀疏脉冲反射系数包括超级分辨率稀疏脉冲纵波速度的反射系数、超级分辨率稀疏脉冲横波速度的反射系数以及超级分辨率稀疏脉冲密度的反射系数。在实际的使用过程中，本发明对反射系数进行综合解释，特别能识别薄层。

为了从有限带宽的低分辨率叠前地震资料获取能够识别薄层、薄互层岩性油气藏的高分辨率的地层三参数（纵波速度、横波速度以及密度）信息，本发明提供了一种叠前地震资料的反演方法。该方法以叠前地震资料相对准确的部分谱信息为数据，在贝叶斯理论框架下，利用相关向量机求解三参数的稀疏脉冲反射系数，而非传统线性AVO反演的带限版本，获得了超级分辨率的效果，提高了反演结果的分辨率和准确性，其计算效率也相对较高，为识别薄层、薄互层岩性油气藏的地层关系和性质提供了技术支持。

图3为本发明实施例提供的一种叠前地震资料反演系统的结构框图，由图3可知，反演系统具体包括：

地震资料采集装置100，用于采集当前地层的叠前地震资料。

在地震勘探中，叠前地震资料中保留了丰富的波场信息，如何利用这些信息来获取地下地层参数为油气检测和识别提供数据支撑是地震勘探最为重要的工作之一。目前地震勘探中主要以一次反射波的动力学理论为基础，在层状地层假设下实施的，比如AVO技术。

频带范围选取装置200，用于从所述的叠前地震资料中选取优势频带范围。

叠前地震资料由于受到地震采集技术限制和噪声的污染，地震资料的有效带宽是有限的，低频和高频成分往往是错误的信息。传统的时间域线性AVO反演无法提高原始地震资料的分辨率，反演结果的准确性也较低。因此，本发明的步骤S102从叠前地震资料中选取出优势频带范围，从有限带宽的叠前地震资料中选取出部分谱（不包括错误的低频和高频的冗余信息）。在具体的实施方式中，选取的优势频带范围应包含峰值频率附近的频点资料。

角道集数据和子波提取装置300，用于从所述的优势频带范围中提取角道集地震数据以及地震子波。假设角道集一次反射波地震数据为s(t,θ)，地震子波为w(t)，角道集反射系数为r(t,θ)。

角道集反射系数确定装置400，用于根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数。根据褶积模型，则有

s(t,θ)=w(t)*r(t,θ)     (6)

根据傅立叶变换转换为频域，则有

S(f,θ)=W(f)R(f,θ)     (7)

其中，S(f,θ)为角道集地震数据，W(f)为地震子波，R(f,θ)为角道集反射系数，f为频率，θ为角度，t为时间。

超参数确定装置500，用于根据相关向量机以及所述的角道集反射系数确定超参数；

超级分辨率反射系数确定装置600，用于根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；

解释装置700，用于对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

图4为本发明实施例提供的一种叠前地震资料反演系统的实施方式二的结构框图，由图4可知，超参数确定装置500具体包括：

矩阵方程建立单元501，用于根据所述的褶积模型建立角道集反射系数的反演矩阵方程。令R(θ_i)=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，若选取角度道集θ_i∈[θ₁,θ_K]，K≥3可以建立三参数反射系数的反演矩阵方程为：

令

，由此，可将公式(13)用矩阵符号简化为d_K*M×1=G_K*M×3Nm_3N×1，其中，m=[r_ρ r_P r_S]^T，d=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，K为选取的角度总数，M为选取的离散频率总数，N为反射系数采样总数，F为傅立叶变换算子，A为密度反射系数的权系数构成的对角矩阵，B为纵波速度反射系数的权系数构成的对角矩阵，C为横波速度反射系数的权系数构成的对角矩阵，r_ρ为密度反射系数，r_p为纵波速度反射系数，r_S为横波速度反射系数。

超参数确定单元502，用于根据似然函数、以零值为中心的标准正态分布以及所述的反演矩阵方程确定超参数。令观测数据误差为n，则公式（13）的反演矩阵的简化方程为：

d=Gm+n=[G₁ G₂...G_k...G_3N][m₁ m₂...m_k...m_3N]^T+n   (14)

其中， $\left[\begin{array}{c}{r}_{\mathrm{\ρ}}={\left[\begin{array}{cccc}{m}_1& m_2& ...& m_N\end{array}\right]}^T\\ r_P={\left[\begin{array}{cccc}{m}_{N+1}& m_{N+2}& ...& m_{2N}\end{array}\right]}^T\\ r_S={\left[\begin{array}{cccc}{m}_{2N+1}& m_{2N+2}& ...& m_{3N}\end{array}\right]}^T\end{array}\right].$

相关向量机的核心思想是利用贝叶斯框架，在最大后验概率下顺序寻找相关向量，并对应求解权值。此处的G_k表示相关向量，m_k表示相关向量的权重。本质地，获得相关向量的过程就是求解反射系数时间位置或确定地层边界的过程，获得权值的过程实质就是求取反射系数大小的过程。

假设n服从均值为0、方差为σ²的高斯分布，即n～N(0,σ²)，在给定模型参数m和σ²的条件下，d的条件概率，也称为似然函数可表示为：

p(d|m,σ²)=(2πσ²)^-K*Mexp[-(d-Gm)^T(d-Gm)/(2σ²)]   (15)

为了获得受地质假设驱动的解，稀疏贝叶斯学习方法需要对m施加先决概率条件，根据贝叶斯理论，限制m的概率分布是以零值为中心的标准正态分布，即

$p\left(m\right|h)={\left(2\mathrm{\π}\right)}^{-3N/2}\underset{k=1}{\overset{3N}{\mathrm{\Π}}}{h}_k^{1/2}\exp ({-h}_k{m}_k^2/2)---\left(16\right)$

其中，h=[h₁,h₂,...,h_3N]^T代表3N个相互独立的超参数，每一个超参数分别控制其对应反射系数大小的先验信息。

表示反射系数m_k服从均值为0，方差为

的高斯分布，其中，k=1，2，…，3N。当

时，m_k→0，表示位置k处没有脉冲反射系数，如此最终会获得稀疏的超级分辨率脉冲反射系数序列。

为了得到反射系数，首先要确定超参数h和方差σ²的最佳值。根据贝叶斯框架，超参数的边缘分布通过下式计算

p(d|h,σ²)=(2π)^3N|Q|exp(d^TQ^-1d)     (17)

其中， $Q={\mathrm{\σ}}^{2}I+{\mathrm{GH}}^{-1}{G}^T={\mathrm{\σ}}^{2}I+\underset{k\≠j}{\mathrm{\Σ}}{h}_k^{-1}{G}_k{G}_k^T+h_j^{-1}{G}_j{G}_j^T$ ，对公式(17)求最小即可求出超参数h，此处的N为反射系数采样总数，p为条件概率。

本发明的反演系统中，超级分辨率反射系数确定装置600，用于根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数。由贝叶斯准则，在d，h和σ²已知的条件下，联立公式(15)和公式(16)，则有超级分辨率稀疏脉冲反射系数m的条件概率或后验概率密度分布：

$p\left(m\right|d,h,{\mathrm{\σ}}^2)=C{\left|\mathrm{\Σ}\right|}^{-\frac{1}{2}}\exp [-\frac{1}{2}{(m-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(m-\mathrm{\μ})]---\left(18\right)$

其中， $\left\{\begin{array}{c}\mathrm{\Σ}={(H+{\mathrm{\σ}}^{-2}{G}^{T}G)}^{-1}\\ \mathrm{\μ}={\mathrm{\σ}}^{-2}\mathrm{\Σ}{G}^{T}d\end{array}\right.---\left(19\right)$

其中，C为一常数，Σ为协方差，μ为均值，即超级分辨率稀疏脉冲反射系数，H=diag(h₁,h₂,...,h_3N)是一个对角矩阵。超级分辨率稀疏脉冲反射系数m的估计由反射系数后验分布的均值μ给出。

在具体的实施方式中，叠前超级分辨率反演方法利用相关向量机算法对公式(17）和公式(19)进行循环迭代求解。每次迭代只更新一个反射系数脉冲m_k对应的超参数h_k和一个相关向量G_k，但每次更新后目标函数或公式(17)都会减小。反复迭代，利用公式(17)计算h，公式(19)计算μ，直到前后两次迭代的目标函数之差达到容忍误差。最后，输出的μ即是估计的稀疏反射系数。实质上，通过公式(17)每次迭代求出超参数h，相当于获得稀疏脉冲的位置，而通过公式(19)求出μ，相当于获得稀疏脉冲的大小。由于相关向量机在迭代运算时，只是通过添加或删除基向量G_k来实时更新脉冲的位置和大小，不需要进行大矩阵的求逆运算，所以该方法具有很快的收敛速度。

反演装置700，用于对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。超级分辨率稀疏脉冲反射系数包括超级分辨率稀疏脉冲纵波速度的反射系数、超级分辨率稀疏脉冲横波速度的反射系数以及超级分辨率稀疏脉冲密度的反射系数。在实际的使用过程中，本发明对反射系数进行综合解释，特别能识别薄层。

本发明只需要叠前地震资料和地震子波的部分谱信息，所需数据满足目前勘探的采集条件。利用贝叶斯理论框架，考虑到了待反演参数的概率分布特征，提高了反演结果的稳定性。采用的相关向量机算法可以得到稀疏的脉冲反射系数，补偿原始记录中没有采集到的频率成分，具有超级分辨率的效果，与传统的线性反演结果的带限性有着本质不同。此外，二维反演结果的横向特征保持较好，迭代收敛速度较快，各中心点位置的参数独立反演，因此可以方便地进行并行计算，为实际应用提供了广阔前景。

为了测试本发明提供的叠前地震资料的反演方法以及系统的可行性，下面结合具体的模型进行试验。

1、以一层状薄层模型为例进行试验，其对应的三参数数据模型（纵波速度V_P、横波速度V_S和密度ρ）如图6所示。由该模型得到的三参数反射系数数据模型（纵波反射系数r_p、横波反射系数rs和密度反射系数r_ρ）如图8中的连续线脉冲。叠前角道集一次反射波地震记录模拟和相应的三参数反射系数反演均选取主频为30Hz的雷克子波。根据瑞利（Rayleigh）准则，可以计算该子波相应的薄层调谐厚度约16ms。换言之，模型中厚度小于16ms的地层均可认为是薄层，如果连续2层以上的层厚小于16ms，则可看成薄互层。

通过公式(3)和公式(6)可以合成试验用的叠前角道集一次反射波反射系数和地震记录，如图7所示。从图6中可以看出，该层状模型是独立薄层和薄互层的组合模型，有10个地层界面，地层速度变化不是特别大，因此模拟的叠前角道集是合理的。其相应的三参数反射系数振幅强弱相间，比较符合实际地层的变化规律。由合成的角道集叠前反射系数和地震记录对比可知，带限的地震记录无法准确的识别出薄层的界面和位置，如果能够反演出如图8所示的结果（连续线脉冲），就能清楚地识别薄层、薄互层厚度和界面位置以及三参数的准确变化规律。

图8为利用叠前地震资料的反演方法反演的三参数反射系数结果（图中断线为反演的结果，连续线脉冲是原始模型），反演带宽选择5-120Hz。从图8中可见，选取三组角道集叠前有限带宽地震数据的优势频带范围，均能正确地反演出三参数脉冲反射系数，达到了识别薄层的超级分辨率效果，而且幅值也同原始模型匹配，保持了原始模型的全部动力学性质。

2、为了进一步验证该反演技术的可行性，以一个推覆体二维模型进行试验。该模型的三参数数据（纵波速度V_P、横波速度V_S和密度ρ）如图9所示，该模型中有薄层、薄互层、独立的透镜体和断层结构。最薄地层厚度为5ms，且存在过渡层（反射系数为较弱的连续脉冲），通常这种近似的层状模型反演难度较高。角道集地震数据的正演和三参数反射系数的反演均采用主频为30Hz的雷克子波。该模型对应的三参数反射系数（纵波反射系数r_p、横波反射系数r_s和密度反射系数r_ρ）如图10所示。相应的反演结果（纵波反射系数r_p_inv、横波反射系数r_s_inv和密度反射系数r_ρ_inv）如图11所示。反演带宽为5-100Hz，输入的叠前角道集地震数据的角度θ=[10263650]。从图11的反演结果可以看出，二维模型中的薄层、薄互层、断层界面以及透镜体边界都刻画的非常清晰，幅值也基本同原始模型匹配，保持了原始模型的动力学特性。此外，对于该二维模型，反演结果保持了很好的横向稳定性，反演出了非带限的脉冲反射系数，实现了超级分辨率的效果。因此，本发明提供的叠前地震资料的反演方法及系统是行之有效的。

综上所述，本发明的有益成果是：提供了一种叠前地震资料的反演方法及系统，通过以叠前地震资料中的优势频带范围为数据，利用贝叶斯理论框架，考虑了待反演参数的概率分布特征，提高了反演结果的稳定性，采用相关向量机算法得到稀疏的脉冲反射系数，补偿原始记录中没有采集到的频率成分，具有超级分辨率的效果。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一般计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-OnlyMemory，ROM）或随机存储记忆体（Random Access Memory，RAM）等。

本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

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

Claims (12)

1.一种叠前地震资料的反演方法，其特征是，所述的方法包括：

采集当前地层的叠前地震资料；

从所述的叠前地震资料中选取优势频带范围；

从所述的优势频带范围中提取角道集地震数据以及地震子波；

根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数；

根据相关向量机以及所述的角道集反射系数确定超参数；

根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；

对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

2.根据权利要求1所述的方法，其特征是，根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数通过如下公式进行：

S(f,θ)=W(f)R(f,θ)

其中，S(f,θ)为角道集地震数据，W(f)为地震子波，R(f,θ)为角道集反射系数，f为频率，θ为角度。

3.根据权利要求2所述的方法，其特征是，根据相关向量机以及所述的角道集反射系数确定超参数具体包括：

根据所述的褶积模型建立角道集反射系数的反演矩阵方程；

根据似然函数、以零值为中心的标准正态分布以及所述的反演矩阵方程确定超参数。

4.根据权利要求3所述的方法，其特征是，根据所述的褶积模型建立的角道集反射系数的反演矩阵方程为：

上述公式对应的矩阵简化方程为d_K*M×1=G_K*M×3Nm_3N×1；

其中，R(θ_i)=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，θ_i∈[θ₁,θ_K]，m=[r_ρ r_P r_S]^T，K为选取的角度总数，M为选取的离散频率总数，N为反射系数采样总数，F为傅立叶变换算子，A为密度反射系数的权系数构成的对角矩阵，B为纵波速度反射系数的权系数构成的对角矩阵，C为横波速度反射系数的权系数构成的对角矩阵，r_ρ为密度反射系数，r_p为纵波速度反射系数，r_S为横波速度反射系数。

5.根据权利要求4所述的方法，其特征是，根据似然函数、以零值为中心的标准正态分布以及所述的反演矩阵方程确定超参数通过如下公式进行：

p(d|h,σ²)=(2π)^3N|Q|exp(d^TQ^-1d)

其中，p为条件概率，h为超参数，N为反射系数采样总数，σ²为方差， $Q={\mathrm{\σ}}^{2}I+{\mathrm{GH}}^{-1}{G}^T={\mathrm{\σ}}^{2}I+\underset{k\≠j}{\mathrm{\Σ}}{h}_k^{-1}{G}_k{G}_k^T+h_j^{-1}{G}_j{G}_j^T.$

6.根据权利要求5所述的方法，其特征是，根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数通过如下公式进行：

$p\left(m\right|d,h,{\mathrm{\σ}}^2)=C{\left|\mathrm{\Σ}\right|}^{-\frac{1}{2}}\exp [-\frac{1}{2}{(m-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(m-\mathrm{\μ})]$

$\left\{\begin{array}{c}\mathrm{\Σ}={(H+{\mathrm{\σ}}^{-2}{G}^{T}G)}^{-1}\\ \mathrm{\μ}={\mathrm{\σ}}^{-2}\mathrm{\Σ}{G}^{T}d\end{array}\right.$

其中，C为一常数，Σ为协方差，μ为均值，即超级分辨率稀疏脉冲反射系数，H=diag(h₁,h₂,...,h_3N)为一个对角矩阵。

7.一种叠前地震资料的反演系统，其特征是，所述的系统包括：

地震资料采集装置，用于采集当前地层的叠前地震资料；

频带范围选取装置，用于从所述的叠前地震资料中选取优势频带范围；

角道集数据和子波提取装置，用于从所述的优势频带范围中提取角道集地震数据以及地震子波；

角道集反射系数确定装置，用于根据褶积模型、所述的角道集地震数据、所述的地震子波确定角道集反射系数；

超参数确定装置，用于根据相关向量机以及所述的角道集反射系数确定超参数；

超级分辨率反射系数确定装置，用于根据所述的超参数、所述的角道集反射系数以及贝叶斯准则确定超级分辨率稀疏脉冲反射系数；

解释装置，用于对所述的超级分辨率稀疏脉冲反射系数进行综合解释，得到当前地层的层位分布。

8.根据权利要求7所述的系统，其特征是，所述的角道集反射系数确定装置通过如下公式进行：

S(f,θ)=W(f)R(f,θ)

其中，S(f,θ)为角道集地震数据，W(f)为地震子波，R(f,θ)为角道集反射系数，f为频率，θ为角度。

9.根据权利要求8所述的系统，其特征是，所述的超参数确定装置具体包括：

矩阵方程建立单元，用于根据所述的褶积模型建立角道集反射系数的反演矩阵方程；

超参数确定单元，用于根据似然函数、以零值为中心的标准正态分布以及所述的反演矩阵方程确定超参数。

10.根据权利要求9所述的系统，其特征是，所述的矩阵方程建立单元通过如下公式进行：

上述公式对应的矩阵简化方程为d_K*M×1=G_K*M×3Nm_3N×1；

其中，R(θ_i)=[R(f₁,θ_i)R(f₂,θ_i)...R(f_M,θ_i)]^T，θ_i∈[θ₁,θ_K]，m=[r_ρ r_P r_S]^T，K为选取的角度总数，M为选取的离散频率总数，N为反射系数采样总数，F为傅立叶变换算子，A为密度反射系数的权系数构成的对角矩阵，B为纵波速度反射系数的权系数构成的对角矩阵，C为横波速度反射系数的权系数构成的对角矩阵，r_ρ为密度反射系数，r_p为纵波速度反射系数，r_S为横波速度反射系数。

11.根据权利要求10所述的系统，其特征是，所述的超参数确定单元通过如下公式进行：

p(d|h,σ²)=(2π)^3N|Q|exp(d^TQ^-1d)

其中，p为条件概率，h为超参数，σ²为方差，N为反射系数采样总数， $Q={\mathrm{\σ}}^{2}I+{\mathrm{GH}}^{-1}{G}^T={\mathrm{\σ}}^{2}I+\underset{k\≠j}{\mathrm{\Σ}}{h}_k^{-1}{G}_k{G}_k^T+h_j^{-1}{G}_j{G}_j^T.$

12.根据权利要求11所述的系统，其特征是，所述的超级分辨率反射系数确定装置通过如下公式进行：

$p\left(m\right|d,h,{\mathrm{\σ}}^2)=C{\left|\mathrm{\Σ}\right|}^{-\frac{1}{2}}\exp [-\frac{1}{2}{(m-\mathrm{\μ})}^T{\mathrm{\Σ}}^{-1}(m-\mathrm{\μ})]$

$\left\{\begin{array}{c}\mathrm{\Σ}={(H+{\mathrm{\σ}}^{-2}{G}^{T}G)}^{-1}\\ \mathrm{\μ}={\mathrm{\σ}}^{-2}\mathrm{\Σ}{G}^{T}d\end{array}\right.$

其中，C为一常数，Σ为协方差，μ为均值，即超级分辨率稀疏脉冲反射系数，H=diag(h₁,h₂,...,h_3N)为一个对角矩阵。

Images