一种地震数据的叠前同步反演方法与流程

本发明涉及油气勘探中地震资料解释与反演技术领域，尤其是涉及一种地震数据的叠前同步反演方法。

背景技术：

叠前地震数据中包含了丰富的偏移距信息，振幅随偏移距、入射角的变化揭示了地下介质的岩性变化和孔隙内流体成分变化。因此利用叠前ava(amplitudeversusangle)同步反演可以从角度部分叠加地震数据中提取多种岩石弹性参数，进一步降低了储层预测和流体识别中的不确定性。

然而常规叠前地震反演存在病态性和多解性，其原因主要有三方面：(1)子波的带限特性。地震子波频带中缺失低频信息和高频信息，只保留了地震的中频信息，因此反演中的子波褶积矩阵是不满秩的，这是所有地震反演多解性的根本原因；(2)介质弹性参数反射率向角度反射系数正演的不精确性以及地震噪声的存在。通常正演矩阵的广义逆矩阵的条件数较高，当地震数据中不含噪声时，可以较为精确地获取地下介质参数，但当地震数据中含有随机噪声时，较高的条件数会急剧降低反演的稳定性；(3)通常在有地质背景的条件下，三参数反射率(纵波阻抗/速度反射率、横波阻抗/速度反射率和密度反射率)并不是独立不相关的，这种不相关会造成反演的不适定性，尤其是密度的反演结果无法准确获取，然而密度结果通常对岩性和流体识别具有重要的意义。

针对上述问题，前人对叠前地震反演问题进行了深入的研究。buland和more在2002年针对叠前地震反演的强不适定性，将贝叶斯理论引入地震反演，即以通过假设地震数据的似然函数和模型参数的先验分布均服从多变量高斯分布，给出了模型参数后验分布的均值和方差的解析解，并指出了确定性反演的解为模型参数后验分布的期望；而随机反演的解，可以通过mcmc等技术从后验概率中抽样实现。downton在2004年在buland和more的基础上，提出“长尾巴”分布相对高斯分布可以进一步提高反演的垂向分辨率，并借助于参数协方差矩阵的去相关技术去除三参数的相关性，提高了叠前反演的不适定性，针对长尾巴分布中的lp范数分布、单变量柯西分布和huber分布分别实现了叠前二项和三项反演。hampson和rusell在2005年通过引入岩石物理经验公式作为约束，进一步提高了叠前反演的稳定性，实现了基于角道集的叠前三参数同步反演算法，该算法为hrs软件的叠前反演模块的核心算法。

然而上述叠前反演算法仅仅考虑了反射系数的稀疏约束，在实际资料应用中单项约束并不能获取稳定的反演结果；另外地震资料中缺失0-10hz的低频成分，因此直接从地震数据中获取的反演结果的低频信息是不准确，反演剖面通常出现不连续的“门帘”效应，需要后续进行低频补偿等二次处理流程，这势必会增加额外的计算量。

技术实现要素：

针对现有技术中所存在的上述技术问题，本发明提出了一种地震数据的叠前同步反演方法，提高了三参数反演结果的稀疏性和垂向分辨率，进一步降低了三参数反演的不适定性，尤其是提高了密度结果的稳定性。该发明同时涉及到模型趋势约束项，实现了反演结果的低频成分同模型的低频成分相一致，保证了反演的低频成分的正确性。

本申请提出了一种地震数据的叠前同步反演方法，包括以下步骤：

s1、构建地震数据拟合差项；

s2、构建反射系数稀疏约束项；

s3、构建模型趋势约束项；

s4、根据所述地震数据拟合差项、反射系数稀疏约束项和模型趋势约束项获取三参数反演结果。

在一个实施方式中，所述步骤s1具体包括：

s11、获取角度子波和角度部分叠加地震数据；

1)读入m个角度子波w(θ1)，w(θ2)，…w(θm)，其中每个角度子波包含了子波起始时间、子波长度、子波采样率和子波数据本身够长，并统计出每个角度子波的绝对值最大值w(θi)max；

2)读入m个角度部分叠加数据d(θ1)，d(θ2)，…d(θm)，其中每个角度部分叠加数据包含了n个采样点；

s12、考虑到反演求逆过程的稳定性，对角度子波进行归一化处理，即：

wuni(θi)＝w(θi)/w(θi)max(2)

然后获取角度子波褶积矩阵；

s13、对角度部分地叠加地震数据进行归一化处理，即：

s(θi)＝d(θi)/w(θi)max(4)

s14、利用角度子波褶积矩阵和归一化处理后的角度部分叠加地震数据构建地震数据拟合差项。

在一个实施方式中，所述步骤s2具体包括：

s21、从地震数据的共深度点道集中获取纵波阻抗模型p_mod、横波阻抗模型s_mod和密度模型ρ_mod；

s22、获取纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率；

s23、获取纵波阻抗模型的反射率的t分布的尺度矩阵、横波阻抗模型的反射率的t分布的尺度矩阵和密度模型的反射率的t分布的尺度矩阵，构造反射系数稀疏约束项。

所述步骤s23中，考虑到纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率在有地质背景的情况下并非独立不相关，假设其满足t分布，在n个采样点的情况下，其联合概率密度函数为：

其中，mp，ms，mρ分别代表纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数，rp，rs，rρ分别为纵波阻抗模型的弹性参数反射率、横波阻抗模型的弹性参数反射率和密度模型的弹性参数反射率；

rp，rs，rρ同反演参数m之间存在如下关系：

其中，d为一阶差分矩阵：

结合公式(5)和(6)，则弹性参数的自然对数服从如下概率分布：

其中ψ为t分布的尺度矩阵，其中γi为3×3n维矩阵，其表达式为：

在一个实施方式中，所述步骤s3具体包括：

s31、获取纵波阻抗模型、横波阻抗模型和密度模型的自然对数，相当于对纵波阻抗、横波阻抗和密度反射率做道积分。

s32、分别对所述纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数进行0-10hz的低通滤波，获取纵波阻抗自然对数的低频趋势、横波阻抗自然对数的低频趋势和密度自然对数的低频趋势，如公式(10)、(11)和公式(12)所示；

ptrend＝filter_lowpass[log(p_mod)](10)

strend＝filter_lowpass[log(s_mod)](11)

ρtrend＝filter_lowpass[log(ρ_mod)](12)

s33、根据所述纵波阻抗模型、横波阻抗模型和密度模型获取低通滤波矩阵l；

其中θ＝f^-1ωf，f表示离散傅立叶变换矩阵其表达式如下所示：

其中ω表示由汉明窗函数低通滤波器组成的对角阵。

以纵波阻抗为例，fmp表示将反演结果变换到频率域，ωfmp表示在频率域取反演结果的低频成分，f^-1ωfmp表示将低频成分变换到时间域，然而由于该计算过程涉及到离散傅立叶变换，当反演的计算时窗较大时，计算效率较低，为了提高计算效率，θ的计算可以变换到时间域中进行，根据傅立叶变换的性质，频率域相乘等价在时间域中褶积。在频率域中，低通滤波可以看作为：

diag(ω)·fft(mp)，fft表示傅立叶变换，diag表示取矩阵ω的对角线元素；该过程在时间域中对应于ifft表示傅立叶反变换，表示圆卷积，而圆卷积又可以表示矩阵相乘的形式，因此θ的快速构建形式中避免了复系数矩阵的相乘操作，可以大大提高计算效率，并且θ可以直接以显式的方式给出矩阵中每个元素赋值，其计算步骤如下：

首先取矩阵ω的对角线元素diag(ω)，然后计算对角线元素的反傅里叶变换为xx＝ifft[diag(ω)]，最后将圆卷积表达成矩阵形式就是θ：

xxi表示xx向量中第i个元素；可以看出θ的快速计算过程中只涉及到了一次反快速傅里叶变换，因此计算效率大大提高了。

s34、根据所述低通滤波矩阵和纵波阻抗模型的自然对数的低频趋势、横波阻抗模型的自然对数的低频趋势和密度模型的自然对数的低频趋势构建模型趋势约束项。

在一个实施方式中，所述步骤s4具体包括：

s41、根据纵波阻抗模型的先验分布、横波阻抗模型的先验分布和密度模型的先验分布以及地震噪声的似然函数，获取纵波阻抗模型的后验分布、横波阻抗模型的后验分布和密度模型的后验分布；

其中p(d)为一常数项不影响后验概率密度函数的形状，因此：

p(m|duni)∞p(duni|m)p(m)(17)

假设地震数据中的噪声服从多变量的高斯分布，则似然函数可以表示为：

其中cn表示噪声的协方差矩阵；g代表ava正演算子，其表达式为：

其中a(θi)，b(θi)，c(θi)表示fatti近似式中的纵波阻抗、横波阻抗和密度项的系数，其表达式为：

其中，γ表示背景横纵波速度比。

s42、整合上述公式(8)和(18)，获取纵波阻抗模型、横波阻抗模型和密度模型的后验概率密度函数；

对上述公式(23)两边取对数并乘-1得到反演的目标函数：

但该反演结果中只引入了反射系数稀疏约束项，缺失低频信息，因此还需引入模型趋势约束项。通过低通滤波矩阵l，在l2范数(向量各元素的平方和然后求平方根)准则下使得反演结果的低频成分同模型的低频成分只差最小，其表达式为：

consttrend＝α(lm-mtrend)^t(lm-mtrend)(25)

其中，

ptrend，strend，ρtrend分别为纵波阻抗模型的自然对数的低频趋势，横波阻抗模型的自然对数的低频趋势和密度模型的自然对数的低频趋势，mp，ms，mρ分别代表纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数。

s43、构建反演目标函数，所述反演目标函数表达式为：

γi为3×3n维矩阵，ψ为尺度矩阵，μ为稀疏约束项权重，α为模型趋势约束项权重，duni为归一化处理后的角度部分叠加数据，g为ava正演算子，l为低通滤波矩阵；

d为一阶差分矩阵。

s44、根据反演目标函数获取纵波阻抗模型、横波阻抗模型和密度模型的反演结果。具体的，对公式(1)两边求关于m的一阶偏导数，得出：

m＝(g^tg+μq+αl^tl)^-1(g^tduni+αl^tmtrend)(26)

其中μ为稀疏约束项权重，α为模型趋势约束项权重，q为非一致性加权矩阵，其表达式为：

由于非一致性加权矩阵q中含有m，因此公式(27)不能直接求取，需给定初始值后通过迭代算法求取。

对获取的反演参数m求其指数，便得到最终的反演结果。

p_inv＝exp[m(1：n)](28)

s_inv＝exp[m(1+n：2n)](29)

ρ_inv＝exp[m(1+2n：3n)](30)

在一个实施方式中，尺度矩阵ψ的获取步骤包括：

其中，rp，rs，rρ为纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率，上标j代表第j次迭代，下标i代表第i个采样点。

在一个实施方式中，所述纵波阻抗模型的先验分布、横波阻抗模型的先验分布和密度模型的先验分布以及地震噪声的似然函数均服从多变量高斯分布。

在一个实施方式中，利用最大期望算法获取纵波阻抗模型、横波阻抗模型和密度模型的反射率的尺度矩阵。

与现有技术相比，本发明的优点在于，本发明提高了三参数反演结果的稀疏性和垂向分辨率，进一步降低了三参数反演的不适定性，尤其是提高了密度结果的稳定性。该发明同时涉及到模型趋势约束项，实现了反演结果的低频成分同模型的低频成分相一致，保证了反演的低频成分的正确性。

附图说明

下面将结合附图来对本发明的优选实施例进行详细地描述。在图中：

图1显示了根据本发明的实施例所述的地震数据的叠前同步反演方法的流程图。

图2(a)显示了实际测井数据示意图。

图2(b)显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图2(a)所述的实际测井数据进行正演获取的角道集示意图。

图3显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图2(a)所述的测井数据进行反演的结果对比示意图。

图4显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图2(a)所述的测井数据进行反演的反射率对比示意图。

图5(a)显示了小角度地震叠加数据示意图。

图5(b)显示了中角度地震叠加数据示意图。

图5(c)显示了大角度地震叠加数据示意图。

图6(a)显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图5(a)所述的小角度地震叠加数据进行反演获取的纵波阻抗、横波阻抗和密度的剖面示意图。

图6(b)显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图5(b)所述的中角度地震叠加数据进行反演获取的纵波阻抗、横波阻抗和密度的剖面示意图。

图6(c)显示了利用本发明的实施例所述的地震数据的叠前同步反演方法对图5(c)所述的大角度地震叠加数据进行反演获取的纵波阻抗、横波阻抗和密度的剖面示意图。

图7显示了根据本发明的实施例所述的地震数据的叠前同步反演方法的井旁道反演结果与实际井曲线对比示意图。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例绘制。

具体实施方式

以下结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以互相结合，所形成的技术方案均在本发明的保护范围之内。

如图1所示，本发明公开了提出了一种地震数据的叠前同步反演方法，包括以下步骤：

s1、获取角度子波和角度部分叠加地震数据，并构建地震数据拟合差项；

s2、获取纵波阻抗模型、横波阻抗模型和密度模型，构建反射系数稀疏约束项；

s3、获取纵波阻抗模型的自然对数的低频趋势、横波阻抗模型的自然对数的低频趋势和密度模型的自然对数的低频趋势，并结合低通滤波矩阵构建模型趋势约束项；

s4、根据所述地震数据拟合差项、反射系数稀疏约束项和模型趋势约束项获取三参数反演结果。

在一个实施方式中，所述步骤s1具体包括：

s11、获取角度子波和角度部分叠加地震数据；

1)读入m个角度子波w(θ1)，w(θ2)，…w(θm)，其中每个角度子波包含了子波起始时间、子波长度、子波采样率和子波数据本身够长，并统计出每个角度子波的绝对值最大值w(θi)max；

2)读入m个角度部分叠加数据d(θ1)，d(θ2)，…d(θm)，其中每个角度部分叠加数据包含了n个采样点；

s12、考虑到反演求逆过程的稳定性，对角度子波进行归一化处理，即：

wuni(θi)＝w(θi)/w(θi)max(2)

然后获取角度子波褶积矩阵；

s13、对角度部分地叠加地震数据进行归一化处理，即：

s(θi)＝d(θi)/w(θi)max(4)

将归一化处理后的m个角度部分地叠加数据组成一个mn行1列的矩阵duni。

s14、利用角度子波褶积矩阵和归一化处理后的角度部分叠加地震数据构建地震数据拟合差项。

在一个实施方式中，所述步骤s2具体包括：

s21、从地震数据的共深度点道集中获取纵波阻抗模型p_mod、横波阻抗模型s_mod和密度模型ρ_mod；

s22、获取纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率；反射率的计算公式为：

其中，rx表示纵波阻抗反射率、横波阻抗反射率或密度反射率，x表示纵波阻抗模型、横波阻抗模型或密度模型。

s23、获取纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率的t分布的尺度矩阵，构造反射系数稀疏约束项。

所述步骤s23中，考虑到纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率在有地质背景的情况下并非独立不相关，假设其满足t分布，在n个采样点的情况下，其联合概率密度函数为：

其中，mp，ms，mρ分别代表纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数，rp，rs，rρ分别为纵波阻抗模型的弹性参数反射率、横波阻抗模型的弹性参数反射率和密度模型的弹性参数反射率；

rp，rs，rp同反演参数m之间存在如下关系：

其中，d为一阶差分矩阵：

结合公式(5)和(6)，则弹性参数的自然对数服从如下概率分布：

其中ψ为t分布的尺度矩阵，其中γi为3×3n维矩阵，其表达式为：

在一个实施方式中，所述步骤s3具体包括：

s31、获取纵波阻抗模型、横波阻抗模型和密度模型的自然对数，相当于对纵波阻抗、横波阻抗和密度反射率做道积分。

s32、分别对所述纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数进行0-10hz的低通滤波，获取纵波阻抗自然对数的低频趋势、横波阻抗自然对数的低频趋势和密度自然对数的低频趋势，如公式(10)、(11)和公式(12)所示；

ptrend＝filter_lowpass[log(p_mod)](10)

strend＝filter_lowpass[log(s_mod)](11)

ρtrend＝filter_lowpass[log(ρ_mod)](12)

s33、根据所述纵波阻抗模型、横波阻抗模型和密度模型获取低通滤波矩阵l；

其中θ＝f^-1ωf，f表示离散傅立叶变换矩阵其表达式如下所示：

其中ω表示由汉明窗函数低通滤波器组成的对角阵。

以纵波阻抗为例，fmp表示将反演结果变换到频率域，ωfmp表示在频率域取反演结果的低频成分，f^-1ωfmp表示将低频成分变换到时间域，然而由于该计算过程涉及到离散傅立叶变换，当反演的计算时窗较大时，计算效率较低，为了提高计算效率，θ的计算可以变换到时间域中进行，根据傅立叶变换的性质，频率域相乘等价在时间域中褶积。在频率域中，低通滤波可以看作为：

diag(ω)·fft(mp)，fft表示傅立叶变换，diag表示取矩阵ω的对角线元素；该过程在时间域中对应于ifft表示傅立叶反变换，表示圆卷积，而圆卷积又可以表示矩阵相乘的形式，因此θ的快速构建形式中避免了复系数矩阵的相乘操作，可以大大提高计算效率，并且θ可以直接以显式的方式给出矩阵中每个元素赋值，其计算步骤如下：

首先取矩阵ω的对角线元素diag(ω)，然后计算对角线元素的反傅里叶变换为xx＝ifft[diag(ω)]，最后将圆卷积表达成矩阵形式就是θ：

xxi表示xx向量中第i个元素；可以看出θ的快速计算过程中只涉及到了一次反快速傅里叶变换，因此计算效率大大提高了。

s34、根据所述低通滤波矩阵和纵波阻抗模型的自然对数的低频趋势、横波阻抗模型的自然对数的低频趋势和密度模型的自然对数的低频趋势构建模型趋势约束项。

在一个实施方式中，所述步骤s4具体包括：

s41、根据纵波阻抗模型的先验分布、横波阻抗模型的先验分布和密度模型的先验分布以及地震噪声的似然函数，获取纵波阻抗模型的后验分布、横波阻抗模型的后验分布和密度模型的后验分布；

其中p(d)为一常数项不影响后验概率密度函数的形状，因此：

p(m|duni)∞p(duni|m)p(m)(17)

假设地震数据中的噪声服从多变量的高斯分布，则似然函数可以表示为：

其中cn表示噪声的协方差矩阵；g代表ava正演算子，其表达式为：

其中a(θi)，b(θi)，c(θi)表示fatti近似式中的纵波阻抗、横波阻抗和密度项的系数，其表达式为：

其中，γ表示背景横纵波速度比。

s42、整合上述公式(8)和(18)，获取纵波阻抗模型、横波阻抗模型和密度模型的后验概率密度函数；

对上述公式(23)两边取对数并乘-1得到反演的目标函数：

但该反演结果中只引入了反射系数稀疏约束项，缺失低频信息，因此还需引入模型趋势约束项。通过低通滤波矩阵l，在l2范数准则下使得反演结果的低频成分同模型的低频成分只差最小，其表达式为：

consttrend＝α(lm-mtrend)^t(lm-mtrend)(25)

其中，

ptrend，strend，ρtrend分别为纵波阻抗模型的自然对数的低频趋势，横波阻抗模型的自然对数的低频趋势和密度模型的自然对数的低频趋势，mp，ms，mρ分别代表纵波阻抗模型的自然对数、横波阻抗模型的自然对数和密度模型的自然对数。

s43、构建反演目标函数，所述反演目标函数表达式为：

γi为3×3n维矩阵，ψ为尺度矩阵，μ为稀疏约束项权重，α为模型趋势约束项权重，duni为归一化处理后的角度部分叠加数据，g为ava正演算子，l为低通滤波矩阵；

d为一阶差分矩阵。

s44、根据反演目标函数获取纵波阻抗模型、横波阻抗模型和密度模型的反演结果。具体的，对公式(1)两边求关于m的一阶偏导数，得出：

m＝(g^tg+μq+αl^tl)^-1(g^tduni+αl^tmtrend)(26)

其中μ为稀疏约束项权重，α为模型趋势约束项权重，q为非一致性加权矩阵，其表达式为：

由于非一致性加权矩阵q中含有m，因此公式(27)不能直接求取，需给定初始值后通过迭代算法求取。

对获取的反演参数m求其指数，便得到最终的反演结果。

p_inv＝exp[m(1：n)](28)

s_inv＝exp[m(1+n：2n)](29)

ρ_inv＝exp[m(1+2n：3n)](30)

在一个实施方式中，尺度矩阵ψ的获取步骤包括：

其中，rp，rs，rρ为纵波阻抗模型的反射率、横波阻抗模型的反射率和密度模型的反射率，上标j代表第j次迭代，下标i代表第i个采样点。

在一个实施方式中，所述纵波阻抗模型的先验分布、横波阻抗模型的先验分布和密度模型的先验分布以及地震噪声的似然函数均服从多变量高斯分布。

在一个实施方式中，利用最大期望算法获取纵波阻抗模型、横波阻抗模型和密度模型的反射率的尺度矩阵。

图2(a)和图2(b)给出了理论数据，图2(a)给出了用于验证该反演方法的弹性参数模型数据，实线1，实线3和实线5分别是纵波阻抗、横波阻抗和密度，虚线2，虚线4和虚线6分别是纵波阻抗的低频成分、横波阻抗的低频趋势和密度的低频趋势。在给定入射角分别为10度、20度、30度和40度，利用zoeppritz方程计算角度反射系数，在忽略投射损失和多次波的情况，利用40hz的雷克子波与角度反射系数褶积形成图2(b)中的角道集。

图3为基于图2(a)和图2(b)中的理论数据，给出了反射系数稀疏约束项权重为1，模型趋势约束项权重为1，迭代次数为5次的反演结果，从左到右依次是纵波阻抗、横波阻抗、密度以及纵横波速度比的理论数据与反演数据叠合对比。图中实线a，实线c和实线e分别为纵波阻抗、横波阻抗、密度的反演结果，虚线b，虚线d和虚线f分别为纵波阻抗、横波阻抗、密度反演结果的低频成分，线g为纵横波速度比的理论数据，线h为纵横波速度比的反演结果。从图中的对比可以看出，反演结果与模型数据吻合度较高，并且反演结果的低频成分同模型数据的低频成分吻合度也高，反映了模型趋势约束项在反演中起到了重要的作用。

图4为基于图2(a)和图2(b)中的理论数据，给出了纵波阻抗反射率、横波阻抗反射率和密度反射率反演结果对比，其中实线a，实线c和实线e为模型数据的反射率，虚线b，虚线d和虚线f为反演结果得到的反射率，可以看出，反射率反演结果呈现明显的尖脉冲效果，说明该反演方法可以有效压缩子波，提高反演的垂向分辨率。

图5(a)、5(b)和5(c)为验证本发明所选用的实际地震数据，图5(a)、5(b)和5(c)分别表示小角度叠加数据、中角度叠加数据和大角度叠加数据，其中心角度分别为6度、18度和28度，反演时窗为图中显示的两个层位之间。

图6(a)、6(b)、6(c)给出了对图5(a)、5(b)和5(c)中实际数据反演获取的纵波阻抗、横波阻抗和密度剖面，相对地震剖面，反演剖面物理意义更加明确丰富。图6(a)、6(b)、6(c)分别为小角度地震叠加数据、中角度地震叠加数据和大角度地震叠加数据的纵波阻抗、横波阻抗、密度反演结果，从图中可以看出，反演剖面相对地震剖面能更加清晰的展示储层岩性、物性的横向变化。

图7给出了井旁道反演结果与实际测井数据的对比显示，该井位于inline1200，xline1388处，图中线4，线6和线8分别为经过250hzbackus滤波后的纵波阻抗，横波阻抗和密度曲线，线5，线7和线8分别为反演得到的纵波阻抗，横波阻抗和密度曲线结果，从图中对比可以看出，反演结果与实测井曲线吻合度较高。

与现有技术相比，本发明的优点在于，本发明提高了三参数反演结果的稀疏性和垂向分辨率，进一步降低了三参数反演的不适定性，尤其是提高了密度结果的稳定性。该发明同时涉及到模型趋势约束项，实现了反演结果的低频成分同模型的低频成分相一致，保证了反演的低频成分的正确性。

以上所述仅为本发明的优选实施方式，但本发明保护范围并不局限于此，任何本领域的技术人员在本发明公开的技术范围内，可容易地进行改变或变化，而这种改变或变化都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。