一种反演方法、系统、装置及存储介质与流程

1.本发明属于油气探勘中地震资料流体识别领域技术，特别涉及一种高精度的流体因子和/或剪切模量和/或密度反演方法、装置系统及存储介质。


背景技术：

2.叠前地震数据中包含了丰富的偏移距信息，利用叠前avo反演可以获取多种弹性参数或弹性参数反射率，从而可以构建对流体因子更为敏感的弹性参数组合，从而实现基于数据的流体识别。smith和goodway提出利用加权叠加技术从叠前道集中获取纵波阻抗和横波阻抗反射率，进而实现流体预测；goodway在此基础上，又提出了利用拉梅参数和密度的组合进行流体识别；george等提出了流体因子角度和交会图角度的概念，通过对这两种角度进行计算可以得到流体因子；mark等提出了泊松阻抗的概念；russell等在考虑多孔饱和流体岩石的前提下，基于biot-gassmann方程对饱和流体条件下的纵波速度方程进行了改写，提出了gassmann流体因子，在此基础上，通过对aki-richards近似进行重组，建立了关于gassmann流体因子、剪切模量与密度反射率表达的russell近似方程，进而可以通过叠前avo反演的方式获取gassmann流体因子。前人的研究表明，gassmann流体因子相对其他弹性参数，对流体具有更高的敏感性。
3.张锐等为了进一步提高反演结果的块化效果，通过结合反射系数奇偶分解理论与信号的稀疏表示理论，实现了稀疏层反演，并将其扩展至叠前反演；稀疏层反演相比稀疏脉冲反演具有更明显的块化效果，可以有效提高岩性或流体分界面的识别精度。
4.通过将叠前ava稀疏层反演和russell近似方程相结合，期望获取具有块化效果的演结果，以提高流体识别的精度，但是在实际应用中发现，由于受限于流体因子对大入射角的敏感性以及稀疏层约束的负面影响，造成常规基于russell近似的叠前稀疏层反演获取的流体因子反演结果的稳定性和精度均较低。


技术实现要素：

5.针对现有技术的不足，本发明的目的在于提供一种高精度的流体因子和/或剪切模量和/或密度反演方法。该方法有效解决了常规的基于russell近似的叠前稀疏层反演稳定性较差的问题，能够获取高精度、高鲁棒性的流体因子和/或剪切模量和/或密度反演结果。
6.为了实现上述目的，本发明提供了一种反演方法，其中，该方法包括：
7.利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果；
8.基于所述纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量；
9.基于所述非零反射率的垂向位置向量，构建分界面约束矩阵；
10.以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
11.在上述反演方法中，优选地，所述利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果包括：
12.基于fatti近似方程和角度子波褶积矩阵确定第一avo正演矩阵；
13.基于所述第一avo正演矩阵并引入反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的第一地震数据拟合差项；
14.基于反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的稀疏层约束项；
15.基于所述第一地震数据拟合差项、所述稀疏层约束项构建基于fatti近似的叠前稀疏层反演目标函数即基于fatti近似的关于弹性参数反射率对应的奇偶分量系数的叠前反演目标函数；
16.基于构建的所述目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数；
17.基于所述纵波阻抗反射率对应的奇偶分量系数，得到所述纵波阻抗反射率数据体即为所述目标区尖脉冲化的纵波阻抗反射率反演结果。
18.更优选地，所述第一avo正演矩阵为：
[0019][0020]
其中，g
fatti
为基于fatti近似方程和角度子波褶积矩阵确定的第一avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个所述角度叠加地震数据)；a(θj)、b(θj)和c(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由fatti近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0021][0022][0023][0024]
更优选地，所述第一地震数据拟合差项为：||d-g
fatti
γx||2；
[0025]
其中，g
fatti
为由fatti近似和子波褶积矩阵确定的第一avo正演矩阵；d为角度叠加地震数据；γ为由奇偶分解矩阵确定的块化矩阵；x为弹性参数反射率对应的奇偶分量系数；其中，
[0026]
k指k个所述角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0027]
其中d为奇偶分解矩阵；其中，
[0028]
d＝[de,do]；其中，
[0029][0030][0031]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0032]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0033]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0034]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0035]
更优选地，所述稀疏层约束项为：λ|x|；
[0036]
其中，λ为稀疏层约束权重；x为弹性参数反射率对应的奇偶分量系数；其中，
[0037]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0038]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0039]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0040]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0041]
更优选地，所述基于fatti近似的叠前稀疏层反演目标函数为：
[0042]
f(x)＝||d-g
fatti
γx||2 λ|x|；
[0043]
其中，||d-g
fatti
γx||2为第一地震数据拟合差项；λ|x|为稀疏层约束项；x为弹性参数反射率对应的奇偶分量系数。
[0044]
更优选地，所述基于构建的所述目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数利用原始-对偶对数障碍法进行。
[0045]
更优选地，所述基于所述纵波阻抗反射率对应的奇偶分量系数，得到所述纵波阻抗反射率数据体通过下述公式进行：
[0046]ri
＝dxi[0047]
式中，ri为纵波阻抗反射率；d为奇偶分解矩阵；xi为纵波阻抗反射率对应的奇偶分量系数。在上述反演方法中，优选地，在所述基于所述纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量过程中，使用如下准则确定非零反射率的垂向位置向量：
[0048]
pos(j)＝i,if[δri(i)＜0&δri(i 1)＞0&ri(i)≤-ε]||[δri(i)＞0&δri(i 1)＜0&ri(i)≥ε]；i＝1,2,
…
n；j＝1,2,
…m[0049]
其中，δri(i)＝ri(i 1)-r
i(i)[0050]
即，如果满足：δri(i)＜0并且δri(i 1)＞0并且ri(i)≤-ε，
[0051]
或者满足：δri(i)＞0并且δri(i 1)＜0并且ri(i)≥ε
[0052]
那么i垂向位置即为第j个非零反射率对应的垂向位置(即第j个非零反射率的垂向位置为第i个采样点)；
[0053]
其中，pos(j)为第j个非零反射率对应的垂向位置；n为纵波阻抗反射率的采样个数；m为非零反射率的总个数；ri(i)为i垂向位置处的纵波阻抗反射率；ε为阀值；ε为一个非常小的正数，其数量级通常为10-6
。
[0054]
在上述反演方法中，优选地，构建的所述分界面约束矩阵为系数矩阵；更优选地，
构建的所述分界面约束矩阵为：
[0055]
p
pos(j),j
＝1；j＝1,2,
…m[0056]
其中，p
pos(j),j
为所述分界面约束矩阵的第pos(j)行、第j列元素；pos(j)为第j个非零反射率对应的垂向位置；m为非零反射率的总个数。
[0057]
在上述反演方法中，优选地，所述以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取反演结果包括：
[0058]
基于russell近似方程和角度子波褶积矩阵确定第二avo正演矩阵；
[0059]
基于所述第二avo正演矩阵并引入分界面约束矩阵进行约束，构建关于地层分界面处的反射率的第二地震数据拟合差项；
[0060]
构建引入分界面约束矩阵进行约束的关于地层分界面处的反射率的趋势软约束项；
[0061]
基于所述地震数据拟合差项、所述趋势软约束项构建基于russell近似的分界面约束最小二乘反演目标函数即基于russell近似的关于地层分界面处的反射率的最小二乘反演目标函数；
[0062]
基于构建的所述目标函数进行反演得到地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率；
[0063]
基于所述地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率，计算所有采样点的流体因子反射率和/或剪切模量反射率和/或密度反射率并利用道积分获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0064]
更优选地，所述第二avo正演矩阵为：
[0065][0066]
其中，g'
russell
为基于russell近似方程和所述角度子波褶积矩阵确定的第二avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个所述角度叠加地震数据)；a'(θj)、b'(θj)和c'(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由russell近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0067]
[0068][0069][0070]
更优选地，所述第二地震数据拟合差项为：
[0071]
其中，g'
russell
为基于russell近似方程和所述角度子波褶积矩阵确定的第二avo正演矩阵；d为角度叠加地震数据为地层分界面处的反射率；其中，
[0072]
k指k个所述角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0073]
p为分界面约束矩阵；
[0074]
其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率。
[0075]
更优选地，所述趋势软约束项为：
[0076]
其中，α为趋势软约束权重；为地层分界面处的反射率；m
trend
为地层分界面处的反射率自然对数的低频趋势；ω为由积分矩阵、低通滤波矩阵以及分界面约束矩阵确定的块化矩阵；其中，
[0077]
其中，σ为由积分矩阵；l为低通滤波矩阵；p为分界面约束矩阵；
[0078]
其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率；
[0079]
其中，lnf
trend
为gassmann流体因子自然对数的低频趋势；lnμ
trend
为剪切模量自然对数的低频趋势；lnρ
trend
为密度自然对数的低频趋势。
[0080]
更优选地，所述基于russell近似的分界面约束最小二乘反演目标函数为：
[0081][0082]
其中，为第二地震数据拟合差项；为趋势软约束项。
[0083]
更优选地，计算所有采样点的流体因子反射率通过下述公式进行：
[0084][0085]
其中，p为分界面约束矩阵；为地层分界面处的流体因子反射率；rf为采样点的流体因子反射率；
[0086]
更优选地，计算所有采样点的剪切模量反射率通过下述公式进行：
[0087][0088]
其中，p为分界面约束矩阵；为地层分界面处的剪切模量反射率；r
μ
为采样点的剪切模量反射率；
[0089]
更优选地，计算所有采样点的密度反射率通过下述公式进行：
[0090][0091]
其中，p为分界面约束矩阵；为地层分界面处的密度反射率；r
ρ
为采样点的密度反射率；
[0092]
更优选地，利用道积分获取流体因子反演结果通过如下公式进行：
[0093]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0094]
其中，f0表示第一个采样点对应的流体因子；rf(j)为第j个采样点的流体因子反射率；
[0095]
更优选地，利用道积分获取流体因子反演结果通过如下公式进行：
[0096]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0097]
其中，μ0表示第一个采样点对应的剪切模量；r
μ
(j)为第j个采样点的剪切模量反射率；
[0098]
更优选地，利用道积分获取流体因子反演结果通过如下公式进行：
[0099]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0100]
其中，ρ0表示第一个采样点对应的密度；r
ρ
(j)为第j个采样点的密度反射率。
[0101]
其中，角度子波褶积矩阵可以通过常规方式获取，例如通过下述方式获取：读入目标区角度叠加地震数据以及对应的角度子波，并形成角度子波褶积矩阵。
[0102]
本发明还提供了一种反演系统，其中，该系统包括：
[0103]
纵波阻抗反射率反演模块：用于利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果；
[0104]
非零反射率垂向位置确定模块：用于基于所述纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量；
[0105]
分界面约束矩阵构建模块：用于基于所述非零反射率的垂向位置向量，构建分界面约束矩阵；
[0106]
分界面约束反演模块：用于以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0107]
在上述反演系统中，优选地，所述纵波阻抗反射率反演模块包括：
[0108]
第一确定子模块：用于基于fatti近似方程和角度子波褶积矩阵确定第一avo正演矩阵；
[0109]
第一构建子模块：用于基于所述第一avo正演矩阵并引入反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的第一地震数据拟合差项；
[0110]
第二构建子模块：用于基于反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的稀疏层约束项；
[0111]
第三构建子模块：用于基于所述第一地震数据拟合差项、所述稀疏层约束项构建基于fatti近似的叠前稀疏层反演目标函数即基于fatti近似的关于弹性参数反射率对应的奇偶分量系数的叠前反演目标函数；
[0112]
第一反演子模块：用于基于构建的所述目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数；
[0113]
第二确定子模块：用于基于所述纵波阻抗反射率对应的奇偶分量系数，得到所述纵波阻抗反射率数据体即为所述目标区尖脉冲化的纵波阻抗反射率反演结果。
[0114]
在上述反演系统中，优选地，所述非零反射率垂向位置确定模块在进行所述基于所述纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量过程中，使用如下准则确定非零反射率的垂向位置向量：
[0115]
pos(j)＝i,if[δri(i)＜0&δri(i 1)＞0&ri(i)≤-ε]||
[0116]
[δri(i)＞0&δri(i 1)＜0&ri(i)≥ε]；i＝1,2,
…
n；j＝1,2,
…m[0117]
其中，δri(i)＝ri(i 1)-r
i(i)[0118]
其中，pos(j)为第j个非零反射率对应的垂向位置；n为纵波阻抗反射率的采样个数；m为非零反射率的总个数；ri(i)为i垂向位置处的纵波阻抗反射率；ε为阀值；ε为一个非常小的正数，其数量级通常为10-6
。
[0119]
在上述反演系统中，优选地，所述分界面约束矩阵构建模块构建的所述分界面约束矩阵为系数矩阵；更优选地，构建的所述分界面约束矩阵为：
[0120]
p
pos(j),j
＝1；j＝1,2,
…m[0121]
其中，p
pos(j),j
为所述分界面约束矩阵的第pos(j)行、第j列元素；pos(j)为第j个非零反射率对应的垂向位置；m为非零反射率的总个数。
[0122]
在上述反演系统中，优选地，所述分界面约束反演模块包括：
[0123]
第三确定子模块：用于基于russell近似方程和角度子波褶积矩阵确定第二avo正演矩阵；
[0124]
第四构建子模块：用于基于所述第二avo正演矩阵并引入分界面约束矩阵进行约束，构建关于地层分界面处的反射率的第二地震数据拟合差项；
[0125]
第五构建子模块：构建引入分界面约束矩阵进行约束的关于地层分界面处的反射率的趋势软约束项；
[0126]
第六构建子模块：用于基于所述地震数据拟合差项、所述趋势软约束项构建基于russell近似的分界面约束最小二乘反演目标函数即基于russell近似的关于地层分界面处的反射率的最小二乘反演目标函数；
[0127]
第二反演子模块：用于基于构建的所述目标函数进行反演得到地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率；
[0128]
第四确定子模块：用于基于所述地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率，计算所有采样点的流体因子反射率和/或剪切模量反射率和/或密度反射率并利用道积分获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0129]
本发明还提供了一种反演装置，包括处理器及存储器；其中，
[0130]
存储器，用于存放计算机程序；
[0131]
处理器，用于执行存储器上所存放的程序时，实现上述反演方法的步骤。
[0132]
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现上述反演方法的步骤。
[0133]
现有技术中的反演方法例如基于russell近似的叠前稀疏层反演的稳定性较差。本发明技术方案提出了一种两步法的高精度反演方法流程，首先利用基于fatti近似的叠前稀疏层反演获取尖脉冲化纵波阻抗反射率(在所有弹性参数中，纵波阻抗反射率反演结果是最稳定的)，在此基础上，利用一定的准则确定纵波阻抗反射率中非零反射率的垂向位置向量，并建立分界面约束矩阵，最后利用建立的分界面约束矩阵约束基于russell近似的最小二乘反演获取得到高精度、高鲁棒性的流体因子和/或剪切模量和/或密度反演。总而言之，本发明提供的技术方案通过引入分界面约束，提高了流体因子和/或剪切模量和/或密度反演的稳定性。
附图说明
[0134]
图1为本发明一实施例提供的反演方法的流程示意图。
[0135]
图2为本发明一实施例提供的反演方法中步骤s1的优化示意图。
[0136]
图3为本发明一实施例提供的反演方法中步骤s4的优化示意图。
[0137]
图4为本发明一实施例提供的反演方法的流程示意图。
[0138]
图5为本发明一实施例提供的反演系统的结构示意图。
[0139]
图6为本发明一实施例提供的反演系统中纵波阻抗反射率反演模块：的优化示意图。
[0140]
图7为本发明一实施例提供的反演系统中分界面约束反演模块的优化示意图。
[0141]
图8为本发明一实施例提供的反演装置的结构示意图。
[0142]
图9a为本发明一实施例中入射角为13度的角度叠加地震数据图。
[0143]
图9b为本发明一实施例中入射角为17度的角度叠加地震数据图。
[0144]
图9c为本发明一实施例中入射角为21度的角度叠加地震数据图。
[0145]
图10为本发明一实施例中所使用的角度子波图。
[0146]
图11a为本发明一实施例中流体因子反演结果剖面图。
[0147]
图11b为本发明一实施例中剪切模量反演结果剖面图。
[0148]
图11c为本发明一实施例中密度反演结果剖面图。
[0149]
图12a为本发明一实施例中流体因子井旁道反演结果与井曲线叠合图。
[0150]
图12b为本发明一实施例中剪切模量井旁道反演结果与井曲线叠合图。
[0151]
图12c为本发明一实施例中密度井旁道反演结果与井曲线叠合图。
具体实施方式
[0152]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明的保护范围。
[0153]
下面参考本发明的若干代表性实施方式，详细阐述本发明的原理和精神。
[0154]
参见图1，为了实现上述目的，本发明提供了一种反演方法，其中，该方法包括：
[0155]
步骤s1：利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果；
[0156]
步骤s2：基于纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量(这些位置对应着地层的分界面)；
[0157]
步骤s3：基于非零反射率的垂向位置向量，构建分界面约束矩阵；
[0158]
步骤s4：以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0159]
其中，基于fatti近似的叠前稀疏层反演是指将fatti近似作为叠前反演的正演公式，通过引入例如反射系数奇偶分解理论进行稀疏层反演，将常规的反射系数反演问题转化为薄层厚度和薄层位置的反演问题，相比传统的稀疏脉冲反演，可以获取较高的垂向分辨率。由于角度反射系数对纵波阻抗的变化最敏感，因此在所有弹性参数中纵波阻抗的反
演精度是最高的，且鲁棒性和抗噪性也是最高的，由此本发明利用鲁棒性和抗噪性较高的纵波阻抗反射率反演结果建立分界面约束矩阵，进行后续获取流体因子和/或剪切模量和/或密度反演，更能够保证流体因子和/或剪切模量和/或密度反演结果的高精度以及高鲁棒性。
[0160]
在一实施方式中，所述反演方法在步骤s1之前还包括：获取目标区地震数据；其中，步骤s1以及步骤s2的反演过程中需要基于获取的目标区地震数据进行；
[0161]
所述获取目标区地震数据可以包括获取目标区角度叠加地震数据以及对应的角度子波。
[0162]
参见图2，在一实施方式中，步骤s1利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果进一步包括：
[0163]
步骤s11：基于fatti近似方程和角度子波褶积矩阵确定第一avo正演矩阵；
[0164]
步骤s12：基于第一avo正演矩阵并引入反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的第一地震数据拟合差项；
[0165]
步骤s13：基于反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的稀疏层约束项；
[0166]
步骤s14：基于第一地震数据拟合差项、稀疏层约束项构建基于fatti近似的叠前稀疏层反演目标函数即基于fatti近似的关于弹性参数反射率对应的奇偶分量系数的叠前反演目标函数；
[0167]
步骤s15：基于构建的目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数；
[0168]
步骤s16：基于纵波阻抗反射率对应的奇偶分量系数，得到纵波阻抗反射率数据体即为目标区尖脉冲化的纵波阻抗反射率反演结果。
[0169]
进一步地，所述反演方法在步骤s11之前还包括：
[0170]
基于目标区地震数据形成角度子波褶积矩阵；例如，读入目标区角度叠加地震数据以及对应的角度子波，从而形成角度子波褶积矩阵。
[0171]
进一步地，第一avo正演矩阵为：
[0172][0173]
其中，g
fatti
为基于fatti近似方程和角度子波褶积矩阵确定的第一avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个角度叠加地震数据)；a(θj)、b(θj)和c(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由fatti近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0174]
[0175][0176][0177]
进一步地，第一地震数据拟合差项为：||d-g
fatti
γx||2；
[0178]
其中，g
fatti
为由fatti近似和子波褶积矩阵确定的第一avo正演矩阵；d为角度叠加地震数据；γ为由奇偶分解矩阵确定的块化矩阵，该矩阵起到将奇偶分量系数转化为弹性参数反射率的作用；x为弹性参数反射率对应的奇偶分量系数；其中，k指k个角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0179]
其中d为奇偶分解矩阵；其中，
[0180]
d＝[de,do]；其中，
[0181][0182]
[0183]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0184]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0185]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0186]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0187]
进一步地，稀疏层约束项为：λ|x|；
[0188]
其中，λ为稀疏层约束权重(λ可以结合地震资料品质进行确定，例如可以为0.01)；x为弹性参数反射率对应的奇偶分量系数；其中，
[0189]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0190]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0191]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0192]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0193]
进一步地，基于fatti近似的叠前稀疏层反演目标函数为：
[0194]
f(x)＝||d-g
fatti
γx||2 λ|x|；
[0195]
其中，||d-g
fatti
γx||2为第一地震数据拟合差项；λ|x|为稀疏层约束项；x为弹性参数反射率对应的奇偶分量系数。
[0196]
进一步地，基于构建的目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数利用原始-对偶对数障碍法进行；
[0197]
反演基于fatti近似的叠前稀疏层反演目标函数等价于基追踪优化问题，因此可
以利用原始-对偶对数障碍法进行最优化求解得到包括纵波阻抗反射率对应的奇偶分量系数在内的弹性参数反射率对应的奇偶分量系数。
[0198]
进一步地，基于纵波阻抗反射率对应的奇偶分量系数，得到纵波阻抗反射率数据体通过下述公式进行：
[0199]ri
＝dxi[0200]
式中，ri为纵波阻抗反射率；d为奇偶分解矩阵；xi为纵波阻抗反射率对应的奇偶分量系数。
[0201]
在一实施方式中，在基于纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量过程中，使用如下准则确定非零反射率的垂向位置向量：
[0202]
pos(j)＝i,if[δri(i)＜0&δri(i 1)＞0&ri(i)≤-ε]||[δri(i)＞0&δri(i 1)＜0&ri(i)≥ε]；i＝1,2,
…
n；j＝1,2,
…m[0203]
其中，δri(i)＝ri(i 1)-r
i(i)[0204]
其中，pos(j)为第j个非零反射率对应的垂向位置；n为纵波阻抗反射率的采样个数；m为非零反射率的总个数；ri(i)为i垂向位置处的纵波阻抗反射率；ε为阀值；ε为一个非常小的正数，其数量级通常为10-6
。
[0205]
基于fatti近似的叠前稀疏层反演获取的纵波阻抗反射率为一系列的尖脉冲响应，在这些尖脉冲响应中，有些是实际地层引起的有效反射，有些是噪声引起的弱小反射，通过引入一个非常小的阈值ε，可以确定有效非零反射率对应的垂向位置。
[0206]
在一实施方式中，构建的分界面约束矩阵为系数矩阵；
[0207]
进一步地，构建的分界面约束矩阵为：
[0208]
p
pos(j),j
＝1；j＝1,2,
…m[0209]
其中，p
pos(j),j
为分界面约束矩阵的第pos(j)行、第j列元素；pos(j)为第j个非零反射率对应的垂向位置；m为非零反射率的总个数。
[0210]
参见图3，在一实施方式中，步骤s4以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合进一步包括：
[0211]
步骤s41：基于russell近似方程和角度子波褶积矩阵确定第二avo正演矩阵；
[0212]
步骤s42：基于第二avo正演矩阵并引入分界面约束矩阵进行约束，构建关于地层分界面处的反射率的第二地震数据拟合差项；
[0213]
步骤s43：构建引入分界面约束矩阵进行约束的关于地层分界面处的反射率的趋势软约束项；
[0214]
步骤s44：基于地震数据拟合差项、趋势软约束项构建基于russell近似的分界面约束最小二乘反演目标函数即基于russell近似的关于地层分界面处的反射率的最小二乘反演目标函数；
[0215]
步骤s45：基于构建的目标函数进行反演得到地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率；
[0216]
步骤s46：基于地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率，计算所有采样点的流体因子反射率和/或剪切模量反射率和/或密度反射率并利用道积分获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0217]
进一步地，所述反演方法在步骤s41之前还包括：
[0218]
基于目标区地震数据形成角度子波褶积矩阵；例如，读入目标区角度叠加地震数据以及对应的角度子波，从而形成角度子波褶积矩阵。
[0219]
进一步地，第二avo正演矩阵为：
[0220][0221]
其中，g'
russell
为基于russell近似方程和角度子波褶积矩阵确定的第二avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个角度叠加地震数据)；a'(θj)、b'(θj)和c'(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由russell近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0222][0223][0224][0225]
进一步地，第二地震数据拟合差项为：
[0226]
其中，g'
russell
为基于russell近似方程和角度子波褶积矩阵确定的第二avo正演矩阵；d为角度叠加地震数据；为地层分界面处的反射率；其中，
[0227]
k指k个角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0228]
p为分界面约束矩阵；
[0229]
其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率。
[0230]
进一步地，趋势软约束项为：
[0231]
其中，α为趋势软约束权重(α的值可以根据地震资料品质进行确定的，例如可以为1)；为地层分界面处的反射率；m
trend
为地层分界面处的反射率自然对数的低频趋势；ω为由积分矩阵、低通滤波矩阵以及分界面约束矩阵确定的块化矩阵；其中，其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率；
[0232]
其中，lnf
trend
为gassmann流体因子自然对数的低频趋势；lnμ
trend
为剪切模量自然对数的低频趋势；lnρ
trend
为密度自然对数的低频趋势；其中，地层分界面处的反射率自然对数的低频趋势优选借助测井数据沿层位横向插值后并进行低通滤波获取。
[0233]
进一步地，基于russell近似的分界面约束最小二乘反演目标函数为：
[0234][0235]
其中，为第二地震数据拟合差项；为趋势软约束项。
[0236]
进一步地，计算所有采样点的流体因子反射率通过下述公式进行：
[0237][0238]
其中，p为分界面约束矩阵；为地层分界面处的流体因子反射率；rf为采样点的流体因子反射率。
[0239]
进一步地，计算所有采样点的剪切模量反射率通过下述公式进行：
[0240][0241]
其中，p为分界面约束矩阵；为地层分界面处的剪切模量反射率；r
μ
为采样点的剪切模量反射率。
[0242]
进一步地，计算所有采样点的密度反射率通过下述公式进行：
[0243][0244]
其中，p为分界面约束矩阵；为地层分界面处的密度反射率；r
ρ
为采样点的密度反射率。
[0245]
进一步地，利用道积分获取流体因子反演结果通过如下公式进行：
[0246]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0247]
其中，f0表示第一个采样点对应的流体因子；rf(j)为第j个采样点的流体因子反射率。
[0248]
进一步地，利用道积分获取流体因子反演结果通过如下公式进行：
[0249]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0250]
其中，μ0表示第一个采样点对应的剪切模量；r
μ
(j)为第j个采样点的剪切模量反射率。
[0251]
进一步地，利用道积分获取流体因子反演结果通过如下公式进行：
[0252]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0253]
其中，ρ0表示第一个采样点对应的密度；r
ρ
(j)为第j个采样点的密度反射率。
[0254]
本发明的又一实施例提供了一种反演方法，该方法用于获取某目标区的流体因子，其中，该方法流程如图4所示，该方法包括：
[0255]
步骤一：读入目标区角度叠加地震数据以及对应的角度子波，并形成角度子波褶积矩阵；
[0256]
读入k个角度叠加地震数据d(θ1),d(θ2),
…
d(θk)，读入k个对应的角度子波w(θ1),w(θ2),
…
w(θk)，并形成角度子波褶积矩阵w(θ1),w(θ2),
…
w(θk)；
[0257]
步骤二：利用基于fatti近似的叠前稀疏层反演，获取纵波阻抗反射率ri、横波阻抗反射率rj和密度反射率r
ρ
数据体；其中，纵波阻抗反射率ri即为目标区尖脉冲化的纵波阻抗反射率反演结果；
[0258]
图9a-图9c为本实施读入的实际地震数据；共有3个角度部分叠加数据，分别是入射角为13度的角度部分叠加数据(见图9a)，入射角为17度的角度部分叠加数据(见图9b)以及入射角为21度的角度部分叠加数据(见图9c)；
[0259]
图10为本实施读入的角度子波；其中黑色实线为入射角13度的角度子波，灰色实线为入射角17度的角度子波，灰色虚线为入射角21度的角度子波；
[0260]
2.1：基于fatti近似方程和角度子波褶积矩阵确定第一avo正演矩阵；基于第一avo正演矩阵并引入反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的第一地震数据拟合差项；基于反射系数奇偶分解理论，构建关于弹性参数反射率对应
的奇偶分量系数的稀疏层约束项；基于第一地震数据拟合差项、稀疏层约束项构建基于fatti近似的叠前稀疏层反演目标函数即基于fatti近似的关于弹性参数反射率对应的奇偶分量系数的叠前反演目标函数；
[0261]
构建的基于fatti近似的叠前稀疏层反演目标函数为：f(x)＝||d-g
fatti
γx||2 λ|x|；
[0262]
其中，||d-g
fatti
γx||2为第一地震数据拟合差项；λ|x|为稀疏层约束项；x为弹性参数反射率对应的奇偶分量系数；
[0263]
其中，g
fatti
为基于fatti近似方程和角度子波褶积矩阵确定的第一avo正演矩阵；
[0264]
w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵；a(θj)、b(θj)和c(θj)均为对角矩阵(其中j＝1,2,
…
k)，其对角线元素由fatti近似方程的系数确定；θj为第j个入射角；其中，
[0265][0266][0267][0268]
其中，d为角度叠加地震数据；
[0269]
其中，γ为由奇偶分解矩阵确定的块化矩阵，该矩阵起到将奇偶分量系数转化为
弹性参数反射率的作用；其中，d为奇偶分解矩阵，d＝[de,do]；其中，
[0270][0271]
其中，x为弹性参数反射率对应的奇偶分量系数；其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数；
[0272]
其中，λ为稀疏层约束权重，此实施例中取0.01；
[0273]
2.2：利用原始-对偶对数障碍法，对构建的目标函数进行最优化求解，得到弹性参数反射率对应的奇偶分量系数x；
[0274]
2.3：基于弹性参数反射率对应的奇偶分量系数x，利用如下公式计算高分辨率的弹性参数反射率：
[0275]ri
＝dxi,rj＝dxj,r
ρ
＝dx
ρ
[0276]
其中，获得的纵波阻抗反射率数据体即为目标区尖脉冲化的纵波阻抗反射率反演结果；
[0277]
步骤三：基于纵波阻抗反射率反演结果，利用如下准则确定非零反射率的垂向位置向量(这些位置对应着地层的分界面)：
[0278]
pos(j)＝i,if[δri(i)＜0&δri(i 1)＞0&ri(i)≤-ε]||[δri(i)＞0&δri(i 1)＜0&ri(i)≥ε]；i＝1,2,
…
n；j＝1,2,
…m[0279]
其中，δri(i)＝ri(i 1)-r
i(i)[0280]
其中，pos(j)为第j个非零反射率对应的垂向位置；n为纵波阻抗反射率的采样个数；m为非零反射率的总个数；ri(i)为i垂向位置处的纵波阻抗反射率；ε为一个非常小的正数，取值10-6
；
[0281]
步骤四：基于非零反射率的垂向位置向量，构建分界面约束矩阵p(p为一个稀疏矩阵)；构建的分界面约束矩阵p为：
[0282]
p
pos(j),j
＝1；j＝1,2,
…m[0283]
其中，p
pos(j),j
为分界面约束矩阵的第pos(j)行、第j列元素；pos(j)为第j个非零反射率对应的垂向位置；m为非零反射率的总个数；
[0284]
步骤五：以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取地层分界面处流体因子反射率、剪切模量反射率和密度反射率；
[0285]
5.1：基于russell近似方程和角度子波褶积矩阵确定第二avo正演矩阵；基于第二avo正演矩阵并引入分界面约束矩阵进行约束，构建关于地层分界面处的反射率的第二地震数据拟合差项；构建引入分界面约束矩阵进行约束的关于地层分界面处的反射率的趋势软约束项；基于地震数据拟合差项、趋势软约束项构建基于russell近似的分界面约束最小二乘反演目标函数即基于russell近似的关于地层分界面处的反射率的最小二乘反演目标函数；
[0286]
构建的基于russell近似的分界面约束最小二乘反演目标函数为：
[0287][0288]
其中，为第二地震数据拟合差项；为趋势软约束项；
[0289]
其中，g'
russell
为基于russell近似方程和角度子波褶积矩阵确定的第二avo正演矩阵；其中，
[0290]
w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵；a'(θj)、b'(θj)和c'(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由russell近似方程的系数确定；其中，
[0291]
[0292][0293][0294]
其中，d为角度叠加地震数据；
[0295]
其中，p为分界面约束矩阵；
[0296]
其中，为地层分界面处的反射率；其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率；
[0297]
其中，其中，σ为由积分矩阵；l为低通滤波矩阵；p为分界面约束矩阵；ω为由积分矩阵、低通滤波矩阵以及分界面约束矩阵确定的块化矩阵；
[0298]
其中，α为趋势软约束权重，此实施例取值1；
[0299]
其中，m
trend
为地层分界面处的反射率自然对数的低频趋势；lnf
trend
为gassmann流体因子自然对数的低频趋势；lnμ
trend
为剪切模量自然对数的低频趋势；lnρ
trend
为密度自然对数的低频趋势；其中，地层分界面处的反射率自然对数的低频趋势优选借助测井数据沿层位横向插值后并进行低通滤波获取；
[0300]
5.2：基于构建的目标函数进行反演得到地层分界面处流体因子反射率、剪切模量反射率和密度反射率；
[0301]
步骤六：基于地层分界面处流体因子反射率，计算所有采样点的流体因子反射率
并利用道积分获取流体因子反演结果；
[0302]
基于地层分界面处剪切模量反射率，计算所有采样点的剪切模量反射率并利用道积分获取剪切模量反演结果；
[0303]
基于地层分界面处密度反射率，计算所有采样点的密度反射率并利用道积分获取密度反演结果；
[0304]
计算所有采样点的流体因子反射率通过下述公式进行：
[0305][0306]
其中，p为分界面约束矩阵；为地层分界面处的流体因子反射率；rf为采样点的流体因子反射率；
[0307]
计算所有采样点的剪切模量反射率通过下述公式进行：
[0308][0309]
其中，p为分界面约束矩阵；为地层分界面处的剪切模量反射率；r
μ
为采样点的剪切模量反射率；
[0310]
计算所有采样点的密度反射率通过下述公式进行：
[0311][0312]
其中，p为分界面约束矩阵；为地层分界面处的密度反射率；r
ρ
为采样点的密度反射率；用道积分获取流体因子反演结果通过如下公式进行：
[0313]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0314]
其中，f0表示第一个采样点对应的流体因子；rf(j)为第j个采样点的流体因子反射率；
[0315]
利用道积分获取流体因子反演结果通过如下公式进行：
[0316]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0317]
其中，μ0表示第一个采样点对应的剪切模量；r
μ
(j)为第j个采样点的剪切模量反射率；
[0318]
利用道积分获取流体因子反演结果通过如下公式进行：
[0319]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0320]
其中，ρ0表示第一个采样点对应的密度；r
ρ
(j)为第j个采样点的密度反射率。
[0321]
图11a为本实施例获取的高精度流体因子反演剖面、图11b为本实施例获取的剪切模量反演剖面、图11c为本实施例获取的密度反演剖面，可以看出，相比角度叠加地震数据，反演剖面展示出了丰富的储层及流体的细节信息，并且由于本发明引入了地层分界面约
束，反演剖面呈现明显的块化效果，可有效降低流体及岩性识别过程中的不确定性。验证本实施例所选取的角度叠加地震数据最大入射角仅为21度，本实施例获取的密度反演剖面也是非常稳定的，这是由于该发明通过引入地层分界面约束，有效降低了密度反演的不确定性。
[0322]
图12a-图12c为基于本实施例所获取的井旁道反演结果与井曲线叠合显示(该井位于inline 4003，xline 2135附近)，其中黑线为实际井曲线，灰线为井旁道反演结果。为了进一步验证本实施例所获取流体因子反演结果的精度，图12a为实测含水饱和度曲线与流体因子反演结果的叠合显示，可以看出流体因子反演结果与实测的含水饱和度曲线具有较高的吻合度，并且反演结果呈现明显的块化效果，因此可提高流体界面的识别精度；图12b、图12c分别是剪切模量反演结果与剪切模量实测曲线的叠合显示以及密度反演结果与密度实测曲线的叠合显示，同样其吻合精度也较高。
[0323]
本发明实施例还提供了一种流体因子反演系统，优选地，该系统用于实现上述的方法实施例。
[0324]
图5是根据本发明实施例的流体因子反演系统的结构框图，如图5所示，该系统包括：
[0325]
纵波阻抗反射率反演模块51：用于利用基于fatti近似的叠前稀疏层反演获取目标区尖脉冲化的纵波阻抗反射率反演结果；
[0326]
非零反射率垂向位置确定模块52：用于基于纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量(这些位置对应着地层的分界面)；
[0327]
分界面约束矩阵构建模块53：用于基于非零反射率的垂向位置向量，构建分界面约束矩阵；
[0328]
分界面约束反演模块54：用于以构建的分界面约束矩阵作为约束进行基于russell近似的最小二乘反演获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0329]
在一实施方式中，流体因子反演系统还包括：
[0330]
第一读取模块：用于读入目标区角度叠加地震数据以及对应的角度子波，并形成角度子波褶积矩阵。
[0331]
参见图6，在一实施方式中，纵波阻抗反射率反演模块51进一步包括：
[0332]
第一确定子模块511：用于基于fatti近似方程和角度子波褶积矩阵确定第一avo正演矩阵；
[0333]
第一构建子模块512：用于基于第一avo正演矩阵并引入反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的第一地震数据拟合差项；
[0334]
第二构建子模块513：用于基于反射系数奇偶分解理论，构建关于弹性参数反射率对应的奇偶分量系数的稀疏层约束项；
[0335]
第三构建子模块514：用于基于第一地震数据拟合差项、稀疏层约束项构建基于fatti近似的叠前稀疏层反演目标函数即基于fatti近似的关于弹性参数反射率对应的奇偶分量系数的叠前反演目标函数；
[0336]
第一反演子模块515：用于基于构建的目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数；
[0337]
第二确定子模块516：用于基于纵波阻抗反射率对应的奇偶分量系数，得到纵波阻抗反射率数据体即为目标区尖脉冲化的纵波阻抗反射率反演结果。
[0338]
在一实施方式中，纵波阻抗反射率反演模块51还包括：
[0339]
第一读取子模块：用于读入目标区地震数据形成角度子波褶积矩阵；例如，用于读入目标区角度叠加地震数据以及对应的角度子波，并形成角度子波褶积矩阵。
[0340]
在一实施方式中，第一确定子模块511确定的第一avo正演矩阵为：
[0341][0342]
其中，g
fatti
为基于fatti近似方程和角度子波褶积矩阵确定的第一avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个角度叠加地震数据)；a(θj)、b(θj)和c(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由fatti近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0343][0344][0345][0346]
在一实施方式中，第一构建子模块512构建的第一地震数据拟合差项为：||d-g
fatti
γx||2；
[0347]
其中，g
fatti
为由fatti近似和子波褶积矩阵确定的第一avo正演矩阵；d为角度叠加地震数据；γ为由奇偶分解矩阵确定的块化矩阵，该矩阵起到将奇偶分量系数转化为弹性参数反射率的作用；x为弹性参数反射率对应的奇偶分量系数；其中，
[0348]
k指k个角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0349]
其中d为奇偶分解矩阵；其中，
[0350]
d＝[de,do]；其中，
[0351][0352][0353]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0354]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0355]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0356]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0357]
更优选地，第二构建子模块513构建的稀疏层约束项为：λ|x|；
[0358]
其中，λ为稀疏层约束权重(λ可以结合地震资料品质进行确定，例如可以为0.01)；x为弹性参数反射率对应的奇偶分量系数；其中，
[0359]
其中，xi为纵波阻抗反射率对应的奇偶分量系数，xj为横波阻抗反射率对应的奇偶分量系数，x
ρ
为密度反射率对应的奇偶分量系数；其中，
[0360]
其中，x
i_e
为纵波阻抗反射率对应的偶分量系数，x
i_o
为纵波阻抗反射率对应的奇分量系数；
[0361]
其中，x
j_e
为横波阻抗反射率对应的偶分量系数，x
j_o
为横波阻抗反射率对应的奇分量系数；
[0362]
其中，x
ρ_e
为密度反射率对应的偶分量系数，x
ρ_o
为密度反射率对应的奇分量系数。
[0363]
在一实施方式中，第三构建子模块514构建的基于fatti近似的叠前稀疏层反演目标函数为：
[0364]
f(x)＝||d-g
fatti
γx||2 λ|x|；
[0365]
其中，||d-g
fatti
γx||2为第一地震数据拟合差项；λ|x|为稀疏层约束项；x为弹性参数反射率对应的奇偶分量系数。
[0366]
在一实施方式中，第一反演子模块515在基于构建的目标函数进行反演得到纵波阻抗反射率对应的奇偶分量系数过程中利用原始-对偶对数障碍法进行；
[0367]
反演基于fatti近似的叠前稀疏层反演目标函数等价于基追踪优化问题，因此可以利用原始-对偶对数障碍法进行最优化求解得到包括纵波阻抗反射率对应的奇偶分量系数在内的弹性参数反射率对应的奇偶分量系数。
[0368]
在一实施方式中，第二确定子模块516在基于纵波阻抗反射率对应的奇偶分量系数，得到纵波阻抗反射率数据体过程中通过下述公式进行：
[0369]ri
＝dxi[0370]
式中，ri为纵波阻抗反射率；d为奇偶分解矩阵；xi为纵波阻抗反射率对应的奇偶分量系数。
[0371]
在一实施方式中，非零反射率垂向位置确定模块52在进行基于纵波阻抗反射率反演结果，确定非零反射率的垂向位置向量过程中，使用如下准则确定非零反射率的垂向位置向量：
[0372]
pos(j)＝i,if[δri(i)＜0&δri(i 1)＞0&ri(i)≤-ε]||[δri(i)＞0&δri(i 1)＜0&ri(i)≥ε]；i＝1,2,
…
n；j＝1,2,
…m[0373]
其中，δri(i)＝ri(i 1)-r
i(i)[0374]
其中，pos(j)为第j个非零反射率对应的垂向位置；n为纵波阻抗反射率的采样个数；m为非零反射率的总个数；ri(i)为i垂向位置处的纵波阻抗反射率；ε为阀值；ε为一个非常小的正数，其数量级通常为10-6
。
[0375]
基于fatti近似的叠前稀疏层反演获取的纵波阻抗反射率为一系列的尖脉冲响应，在这些尖脉冲响应中，有些是实际地层引起的有效反射，有些是噪声引起的弱小反射，通过引入一个非常小的阈值ε，可以确定有效非零反射率对应的垂向位置。
[0376]
在一实施方式中，分界面约束矩阵构建模块53构建的分界面约束矩阵为系数矩阵；
[0377]
进一步地，构建的分界面约束矩阵为：
[0378]
p
pos(j),j
＝1；j＝1,2,
…m[0379]
其中，p
pos(j),j
为分界面约束矩阵的第pos(j)行、第j列元素；pos(j)为第j个非零反射率对应的垂向位置；m为非零反射率的总个数。
[0380]
参见图7，在一实施方式中，分界面约束反演模块54进一步包括：
[0381]
第三确定子模块541：用于基于russell近似方程和角度子波褶积矩阵确定第二avo正演矩阵；
[0382]
第四构建子模块542：用于基于第二avo正演矩阵并引入分界面约束矩阵进行约束，构建关于地层分界面处的反射率的第二地震数据拟合差项；
[0383]
第五构建子模块543：构建引入分界面约束矩阵进行约束的关于地层分界面处的反射率的趋势软约束项；
[0384]
第六构建子模块544：用于基于地震数据拟合差项、趋势软约束项构建基于russell近似的分界面约束最小二乘反演目标函数即基于russell近似的关于地层分界面处的反射率的最小二乘反演目标函数；
[0385]
第二反演子模块545：用于基于构建的目标函数进行反演得到地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率；
[0386]
第四确定子模块546：用于基于地层分界面处流体因子反射率和/或剪切模量反射率和/或密度反射率，计算所有采样点的流体因子反射率和/或剪切模量反射率和/或密度反射率并利用道积分获取流体因子反演结果、剪切模量反演结果和密度反演结果中的一种或两种的组合。
[0387]
在一实施方式中，分界面约束反演模块54还包括：
[0388]
第二读取子模块：用于读入目标区地震数据形成角度子波褶积矩阵；例如，用于读入目标区角度叠加地震数据以及对应的角度子波，并形成角度子波褶积矩阵。
[0389]
在一实施方式中，第三确定子模块541确定的第二avo正演矩阵为：
[0390][0391]
其中，g'
russell
为基于russell近似方程和角度子波褶积矩阵确定的第二avo正演矩阵；w(θ1),w(θ2),
…
w(θk)为角度子波褶积矩阵(k指k个角度叠加地震数据)；a'(θj)、b'(θj)和c'(θj)(其中j＝1,2,
…
k)均为对角矩阵，其对角线元素由russell近似方程的系数确定；θj为第j个入射角，其中j＝1,2,
…
k；其中，
[0392][0393][0394][0395]
在一实施方式中，第四构建子模块542构建的第二地震数据拟合差项为：
[0396][0397]
其中，g'
russell
为基于russell近似方程和角度子波褶积矩阵确定的第二avo正演矩阵；d为角度叠加地震数据；为地层分界面处的反射率；其中，
[0398]
k指k个角度叠加地震数据；θj为第j个入射角，其中j＝1,2,
…
k；
[0399]
其中，p为分界面约束矩阵；
[0400]
其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率。
[0401]
在一实施方式中，第五构建子模块543构建的趋势软约束项为：
[0402]
其中，α为趋势软约束权重(α的值可以根据地震资料品质进行确定的，例如可以为1)；为地层分界面处的反射率；m
trend
为地层分界面处的反射率自然对数的低频趋势；ω为
由积分矩阵、低通滤波矩阵以及分界面约束矩阵确定的块化矩阵；其中，
[0403]
其中，σ为由积分矩阵；l为低通滤波矩阵；p为分界面约束矩阵；
[0404]
其中，为地层分界面处的流体因子反射率；为地层分界面处的剪切模量反射率；为地层分界面处的密度反射率；
[0405]
其中，lnf
trend
为gassmann流体因子自然对数的低频趋势；lnμ
trend
为剪切模量自然对数的低频趋势；lnρ
trend
为密度自然对数的低频趋势；其中，地层分界面处的反射率自然对数的低频趋势优选借助测井数据沿层位横向插值后并进行低通滤波获取。
[0406]
在一实施方式中，第六构建子模块544构建的基于russell近似的分界面约束最小二乘反演目标函数为：
[0407][0408]
其中，为第二地震数据拟合差项；为趋势软约束项。
[0409]
在一实施方式中，第二反演子模块545在计算所有采样点的流体因子反射率过程中通过下述公式进行：
[0410][0411]
其中，p为分界面约束矩阵；为地层分界面处的流体因子反射率；rf为采样点的流体因子反射率。
[0412]
在一实施方式中，第二反演子模块545在计算所有采样点的剪切模量反射率过程中通过下述公式进行：
[0413][0414]
其中，p为分界面约束矩阵；为地层分界面处的剪切模量反射率；r
μ
为采样点的剪切模量反射率。
[0415]
在一实施方式中，第二反演子模块545在计算所有采样点的密度反射率过程中通过下述公式进行：
[0416][0417]
其中，p为分界面约束矩阵；为地层分界面处的密度反射率；r
ρ
为采样点的密度反射率。
[0418]
在一实施方式中，第四确定子模块546在利用道积分获取流体因子反演结果过程中通过如下公式进行：
[0419]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0420]
其中，f0表示第一个采样点对应的流体因子；rf(j)为第j个采样点的流体因子反射率。
[0421]
在一实施方式中，第四确定子模块546在利用道积分获取流体因子反演结果过程中通过如下公式进行：
[0422]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0423]
其中，μ0表示第一个采样点对应的剪切模量；r
μ
(j)为第j个采样点的剪切模量反射率。
[0424]
在一实施方式中，第四确定子模块546在利用道积分获取流体因子反演结果过程中通过如下公式进行：
[0425]
其中，i＝1,2,
……
n，n为反演时窗确定的采样点数总个数；
[0426]
其中，ρ0表示第一个采样点对应的密度；r
ρ
(j)为第j个采样点的密度反射率。
[0427]
图8是根据本发明实施例的反演装置的示意图。图8所示的反演装置为通用数据处理装置，其包含通用的计算机硬件结构，其至少包含处理器1000、存储器1111；所述处理器1000用于执行所述存储器中存储的反演程序，以实现各方法实施例所述的反演方法(具体方法参见上述方法实施例的描述，在此不再赘述)。
[0428]
本发明实施例还提供了一种计算机可读存储介质，所述存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现各方法实施例所述的反演方法(具体方法参见上述方法实施例的描述，在此不再赘述)。
[0429]
以上参照附图描述了本发明的优选实施方式。这些实施方式的许多特征和优点根据该详细的说明书是清楚的，因此权利要求旨在覆盖这些实施方式的落入其真实精神和范围内的所有这些特征和优点。此外，由于本领域的技术人员容易想到很多修改和改变，因此不是要将本发明的实施方式限于所例示和描述的精确结构和操作，而是可以涵盖落入其范围内的所有合适修改和等同物。
[0430]
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0431]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程
图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0432]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0433]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0434]
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。