格拉姆相关性约束叠前地震反演方法

1.本发明属于石油地球物理勘探领域，涉及利用叠前地震反演结果以进行油气储层预测与流体识别等多个环节，具体技术为格拉姆相关性约束叠前地震反演方法，解决现有叠前地震反演采用的相关性正则化约束方法不能够准确刻画非平稳相关性的问题，提高叠前地震反演结果的稳定性。
技术背景
2.随着油气勘探不断深入和勘探技术不断进步，基于zoeppritz方程的叠前地震反演一直是学术界和工业界研究与应用的前沿热点。叠前地震反演利用地震振幅随偏移距或入射角、甚至方位角的变化信息进行反演，获取地下储层弹性参数、裂缝参数、各向异性参数、流体敏感参数、脆性指数以及地应力参数等，为油气储层预测与流体识别等提供可靠的数据支撑。
3.叠前地震反演是严重病态的，而造成叠前地震反演问题的不适定性具有多方面因素。首先，实际观测到的叠前地震数据或多或少都包含诸如随机噪声、相干噪声等各种噪声。其次，用于叠前地震反演求解的反问题算子几乎都是病态的，微小的观测数据扰动都会造成反推出的模型参数发生剧烈变化。再者，叠前地震反演同时反演多个不同的模型参数，而不同模型参数之间存在着相关性，进一步导致反演结果不稳定。
4.求解不适定反演问题的有效方法就是对反演进行正则化约束。正则化约束可将不适定的病态反演问题转换成有条件的适定问题，从而提高反演稳定性。在叠前地震反演中，目前最常用的正则化约束有tikhonov正则化、稀疏性正则化，以及模型参数相关性正则化等。
5.在众多正则化方法中，tikhonov正则化在数学物理反问题中应用较为广泛。在叠前地震反演的最初历程中，学者们通常采用吉洪诺夫正则化对模型参数进行约束。但是，tikhonov正则化采用平滑方式处理模型中的不连续性特征，得到的模型参数具有明显的平滑性，不能清晰表征地下地层边界与地质体边缘特征，反演结果的分辨率比较低。稀疏性约束作为一种正则化方法，近些年来，已越来越得到地球物理界特别是地震反演工作者的青睐。在稀疏性约束叠前地震反演中，认为地下地层的弹性参数相对变化率是稀疏分布的，并对其进行稀疏性正则化约束。稀疏性正则化不仅能够实现稳定反演求解的目的，还能更好地刻画地震数据以及模型参数的稀疏性特征，提高叠前地震反演结果的分辨率。但无论是tikhonov正则化还是稀疏性正则化，本质上都是直接对单个模型参数进行约束，并没有考虑叠前地震反演中模型参数间的相关性。
6.目前，关于叠前地震反演中模型参数相关性的研究，包括不同模型参数之间的相关性，以及同一模型参数在相邻地层间的相关性，即同一模型参数在不同采样点的相关性。对于同一模型参数在不同采样点的相关性，可基于tikhonov正则化，在先验低频模型参数约束中引入同一模型参数不同采样点间的协方差进行解决。引入同一模型协方差矩阵的先验低频模型约束相当于在贝叶斯反演理论中，假设同一模型参数同其先验低频模型之间的
差异服从多元高斯分布，因此协方差矩阵刻画了同一模型参数在不同采样点之间的相关性。
7.相对于同一模型参数在相邻地层间的相关性，在叠前地震反演中，不同模型参数间的相关性对反演结果影响更大，是造成反演结果不稳定的一个重要因素。本发明解决的问题是不同模型参数间的相关性引起的反演不稳定性。对于该问题，其解决办法是对不同模型参数进行相关性正则化约束。
8.目前，常用的不同模型参数间相关性正则化约束方法不外乎两种。第一种为硬约束方法。将确定的岩石物理模型、或经验岩石物理关系等相关关系作为硬约束等式条件，引入到叠前地震反演中。第二种即为软约束方法。基于确定的岩石物理模型、经验岩石物理关系等相关关系，构建先验协方差表征不同模型参数间的相关性。无论是硬约束采用的等式约束条件，还是在软约束方法构建的先验协方差，均需要确定的模型参数相关关系，其形式与系数都是确定的，即认为相关性是平稳的。例如，常用的相关关系有纵波速度和密度之间的gardner经验关系式，纵波速度与横波速度之间的castagna经验关系式，各向异性参数岩石物理模型关系等。但是，对于不同岩性或深度的储层，经验关系式的系数通常是不一样的，所采用的岩石物理模型也不一样，即地下介质不同模型参数之间的相关性本质上是非平稳的。现有相关性正则化约束方法均不能准确刻画非平稳相关性，是限制叠前地震反演实际应用效果的一个重要因素。更为合理地进行不同模型参数间相关性正则化约束，提高叠前地震反演稳定性，稳健地提取地下储层的各种参数，仍然是地震反演工作者面临的关键问题。


技术实现要素：

9.本发明的目的在于提供格拉姆相关性约束叠前地震反演方法。该方法利用能够准确刻画非平稳相关性的格拉姆行列式，作为不同模型参数间的相关性测度，不需要确定的模型参数相关关系。再将其作为叠前地震反演的相关性正则化约束，解决现有相关性正则化约束方法不能够准确刻画非平稳相关性的问题。本发明综合了叠前地震反演理论与方法、最优化方法、泛函分析、变分分析、参数估计等多个学科，将叠前地震反演、格拉姆相关性测度、梯度下降算法、无偏风险估计有机结合，有效提高叠前地震反演的稳定性，适用于油气储层预测与流体识别等油气工业多个重要环节。
10.为达到以上技术目的，本发明提供以下技术方案。
11.格拉姆相关性约束叠前地震反演方法，其关键核心技术包括构建格拉姆相关性约束叠前地震反演目标泛函、梯度下降法目标泛函优化求解、格拉姆行列式展开特性分析、无偏风险估计自适应求取格拉姆相关性约束系数。该方法依次包括以下步骤：
12.(1)输入保幅叠前地震数据；
13.(2)输入地震子波并建立叠前地震反演正问题方程；
14.(3)构建初始模型参数与先验低频模型参数；
15.(4)基于正则化反演理论，构建格拉姆相关性约束叠前地震反演目标泛函；
16.(5)基于格拉姆行列式的展开性质，对格拉姆相关性约束进行展开特性分析，计算格拉姆相关性约束的梯度方向；
17.(6)基于格拉姆相关性约束的梯度方向，采用梯度下降算法对反演目标函数进行
迭代求解；
18.(7)判断反演目标函数是否满足迭代停止容限条件；
19.(8)若不满足迭代停止容限条件，通过无偏风险估计自适应求取格拉姆相关性约束系数，并转下一次迭代；
20.(9)若满足迭代停止容限条件，则停止迭代，输出反演得到的模型参数，并换算成最终的反演结果，即各种弹性参数、裂缝参数、各向异性参数、流体敏感参数等。
21.具体步骤为：
22.步骤(1)：输入保幅叠前地震数据，是经过动校正的未叠加的叠前cmp道集数据，或者是已经过部分叠加的多个角度叠加地震数据。叠前地震反演所采用的地震数据经过前期精细地震处理的保幅数据。
23.步骤(2)：输入地震子波并建立叠前地震ava反演正问题方程。
24.基于zoeppritz方程的叠前地震反演；正问题包括描述平面弹性波在弹性分界面上反射与透射现象的zoeppritz方程以及地震数据的褶积模型；
25.对于叠前地震ava反演，其正问题方程的建立过程如下。
26.弹性各向同性介质zoeppritz方程的aki-richards近似式为：
27.r
pp
(θ)＝a(θ)r
p
+b(θ)rs+c(θ)r
ρ
ꢀꢀꢀꢀꢀꢀ
(1)
28.其中，θ为入射角，r
pp
(θ)为不同入射角的纵波反射系数，r
p
、rs与r
ρ
分别为纵波速度、横波速度与密度的相对变化率，即弹性参数相对变化率，a(θ)＝(1+tan2θ)、b(θ)＝-8γ2sin2θ与c(θ)＝(1-4γ2sin2θ)分别纵波速度、横波速度与密度的系数因子，γ为背景纵横波速度比。
29.假设在叠前地震数据中包含m个入射角，且每道地震数据有n个采样点，可将式(1)扩展成矩阵形式：
30.r＝ar
ꢀꢀꢀꢀꢀꢀ
(2)
31.其中，r为角度反射系数矢量，由不同入射角与不同采样点的纵波反射系数组成：
[0032][0033]
r为弹性参数相对变化率矢量，由不同采样点弹性参数相对变化率组成：
[0034]
r＝[r
p
(1),...,r
p
(n),rs(1),...,rs(n),r
ρ
(1),...,r
ρ
(n)]
t
ꢀꢀꢀꢀꢀꢀ
(4)
[0035]
a为由不同入射角对应的a(θ)、b(θ)和c(θ)构成的系数矩阵：
[0036][0037]
由地震数据的褶积模型可知，地震数据为反射系数与地震子波的褶积结果：
[0038]
d＝w*r
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0039]
其中，d是地震数据，w是地震子波，r是纵波反射系数。
[0040]
地震反演的实际数值求解常采用离散形式，因此，式(6)被离散化为：
[0041]
d＝wr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0042]
这时，w是子波褶积矩阵。对于叠前地震反演，d则为叠前地震数据，r则为角度反射系数。
[0043]
合并式(2)与式(7)可得叠前地震ava反演的正问题方程为：
[0044]
d＝wr＝war＝fr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0045]
其中，f＝wa是子波褶积矩阵和系数矩阵a组合成的联合矩阵，即ava反演的正问题算子。
[0046]
步骤(3)：构建初始模型参数与先验低频模型参数。
[0047]
对实际测井资料或伪测井数据进行基线校正、环境校正、刻度均衡等预处理后，在地质解释层位约束下，由克里金插值法建立测井插值模型，作为反演求解的初始模型参数。再对测井插值模型进行低通滤波，得到先验低频模型参数。
[0048]
步骤(4)：基于正则化反演理论，构建格拉姆相关性约束叠前地震反演目标泛函。
[0049]
在叠前地震反演中，目标泛函为：
[0050]
minf(m)＝||d-f(m)||2+αq(m)+μ||m-m
prior
||2,
ꢀꢀꢀꢀꢀ
(9)
[0051]
其中，f则为叠前地震反演的正问题算子，m为叠前地震反演的所有模型参数，对于叠前地震ava反演，m即为式(8)中的r，q(m)为模型参数修正柯西约束，m
prior
为先验低频模型参数，最后一项为先验低频模型参数约束，α与μ为相应的约束系数。
[0052]
将格拉姆相关性测度作为叠前地震反演目标泛函的约束条件，通过求解下面的反演目标泛函：
[0053][0054]
上式中，σ为格拉姆相关性测度约束的容限条件参数，则为模型参数格拉姆相关性测度，其表达式为：
[0055][0056]
其中，m1、m2、
…
、mn为m中包含的多个不同模型参数。对于叠前地震ava反演，则n为3，而m1、m2、m3则为纵波速度、横波速度与密度的相对变化率等，即r
p
、rs与r
ρ
，(m1,m2)为模型参数m1与m2的内积。
[0057]
目标泛函(10)等价地写为如下形式：
[0058]
minf(m)＝||d-f(m)||2+αq(m)+μ||m-m
prior
||2+λgm,
ꢀꢀꢀꢀꢀꢀ
(12)
[0059]
上式(12)即为格拉姆相关性约束叠前地震反演的目标泛函，其中，λ为格拉姆相关性约束系数。
[0060]
步骤(5)：基于格拉姆行列式的展开性质，对格拉姆相关性约束进行展开特性分析，计算格拉姆相关性约束的梯度方向。
[0061]
采用梯度下降算法对目标泛函(12)进行迭代求解。由目标泛函(12)可知，在采用
梯度下降算法进行迭代求解时，除了需要计算观测数据拟合项即目标泛函第一项、修正柯西约束项与先验低频模型参数约束项的梯度方向外，还需要计算格拉姆相关性约束项gm关于模型参数m的梯度方向。
[0062]
对于格拉姆相关性约束项的梯度方向，通过变分法进行求解。格拉姆相关性约束项的一阶变分为：
[0063][0064]
其中，(δmi,mi)g为格拉姆相关性测度关于模型参数mi的变分：
[0065][0066]
由格拉姆矩阵行列式的展开特性，将上式(14)写为：
[0067][0068]
其中，g
ij
是消除式(11)中格拉姆矩阵第i列和第j行形成的子格拉姆矩阵行列式。
[0069]
格拉姆相关性约束项的一阶变分简化为：
[0070][0071]
其中，g(mi)即为格拉姆相关性约束项关于模型参数mi的梯度方向：
[0072][0073]
格拉姆相关性约束项的梯度方向是梯度下降算法的主要组成部分。
[0074]
步骤(6)：基于格拉姆相关性约束的梯度方向，采用梯度下降算法对反演目标函数进行迭代求解。
[0075]
输入当前模型参数，采用梯度下降算法对目标函数(12)进行迭代求解，其迭代公式为：
[0076][0077]
上式(18)中，τk是第k次迭代的迭代步长，由一维搜索法得到，mk为第k次迭代的当前模型参数，第1次迭代的当前模型参数为初始模型参数，m
k+1
为更新后的模型参数，为第k次迭代的梯度方向：
[0078][0079]
其中q为修正柯西约束关于mk的偏导数矩阵，上标t为矩阵转置，为格拉姆相关性约束项在mk的梯度方向。
[0080]
步骤(7)：判断反演目标函数是否满足迭代停止容限条件。
[0081]
若满足下面的迭代停止容限条件，则转步骤(9)：
[0082][0083]
其中ε是迭代容限参数。否则，转步骤(8)并进行下一次迭代。
[0084]
步骤(8)：若不满足迭代停止容限条件，通过无偏风险估计自适应求取格拉姆相关性约束系数，并进行下一次迭代。
[0085]
采用无偏风险估计法，在反演的迭代求解过程中，自适应求取格拉姆相关性约束系数。对于一个给定的格拉姆相关性约束系数，通过梯度下降算法求解反演目标泛函(12)，得到对应于该约束系数的模型参数m
*
。模型参数m
*
是格拉姆相关性约束系数λ的函数：
[0086]m*
＝m
*
(λ).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0087]
无偏风险估计法通过求解模型参数m
*
的无偏风险估计泛函的极小化，在迭代求解过程中，自适应求取格拉姆相关性约束系数。无偏风险估计泛函为：
[0088][0089]
无偏风险估计泛函采用蒙特卡洛法对其进行求解。求出格拉姆相关性约束系数后λ，转步骤(6)进行下一次迭代。
[0090]
步骤(9)：若满足迭代停止容限条件，则停止迭代，输出反演得到的模型参数，并换算成最终的反演结果，即各种弹性参数、裂缝参数、各向异性参数、流体敏感参数。
[0091]
通过弹性参数相对变化率与弹性参数间的积分关系，将弹性参数相对变化率换算成最终反演结果，即纵波速度、横波速度和密度等弹性参数。换算公式为：
[0092][0093][0094][0095]
式(23)至式(25)中，t为时间采样点，v
p
(0)、vs(0)和ρ(0)分别是初始纵波速度、横波速度和密度。
附图说明
[0096]
图1是本发明的流程图。
具体实施方式
[0097]
叠前地震反演按所采用的正问题方程，包括多种具体反演方法，例如ava反演、弹性阻抗反演、方位ava即avaz反演、以及方位弹性阻抗反演等。下面以叠前地震ava反演为例，对本发明的具体实施流程作进一步详细描述。
[0098]
以下步骤对应的流程图见附图1所示。
[0099]
步骤(1)：输入保幅叠前地震数据，可以是经过动校正的未叠加的叠前cmp道集数据，也可以是已经过部分叠加的多个角度叠加地震数据。叠前地震反演所采用的地震数据
必须是经过前期精细地震处理的保幅数据，特别是待反演目的层段地震数据的振幅，必须至少是相对保幅的。
[0100]
步骤(2)：输入地震子波并建立叠前地震ava反演正问题方程。
[0101]
本发明涉及的具体反演方法为基于zoeppritz方程的叠前地震反演。因此，其正问题包括描述平面弹性波在弹性分界面上反射与透射现象的zoeppritz方程以及地震数据的褶积模型。例如，对于ava反演，正问题方程表示弹性各向同性介质zoeppritz方程或其近似式与褶积模型；对于avaz反演，正问题方程则为弹性各向异性介质的zoeppritz方程或其r
ü
ger近似式与褶积模型；对于弹性阻抗反演，则正问题方程为弹性阻抗方程以及褶积模型。对于叠前地震ava反演，其正问题方程的建立过程如下。
[0102]
弹性各向同性介质zoeppritz方程的aki-richards近似式为：
[0103]rpp
(θ)＝a(θ)r
p
+b(θ)rs+c(θ)r
ρ
ꢀꢀꢀꢀꢀ
(1)
[0104]
其中，θ为入射角，r
pp
(θ)为不同入射角的纵波反射系数，r
p
、rs与r
ρ
分别为纵波速度、横波速度与密度的相对变化率，即弹性参数相对变化率，a(θ)＝(1+tan2θ)、b(θ)＝-8γ2sin2θ与c(θ)＝(1-4γ2sin2θ)分别纵波速度、横波速度与密度的系数因子，γ为背景纵横波速度比。
[0105]
假设在叠前地震数据中包含m个入射角，且每道地震数据有n个采样点，可将式(1)扩展成矩阵形式：
[0106]
r＝ar
ꢀꢀꢀꢀꢀꢀ
(2)
[0107]
其中，r为角度反射系数矢量，由不同入射角与不同采样点的纵波反射系数组成：
[0108][0109]
r为弹性参数相对变化率矢量，由不同采样点弹性参数相对变化率组成：
[0110]
r＝[r
p
(1),...,r
p
(n),rs(1),...,rs(n),r
ρ
(1),...,r
ρ
(n)]
t
ꢀꢀꢀꢀꢀ
(4)
[0111]
a为由不同入射角对应的a(θ)、b(θ)和c(θ)构成的系数矩阵：
[0112][0113]
由地震数据的褶积模型可知，地震数据为反射系数与地震子波的褶积结果：
[0114]
d＝w*r
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0115]
其中，d是地震数据，w是地震子波，r是纵波反射系数。
[0116]
地震反演的实际数值求解常采用离散形式，因此，式(6)被离散化为：
[0117]
d＝wr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0118]
这时，w是子波褶积矩阵。对于叠前地震反演，d则为叠前地震数据，r则为角度反射系数。
[0119]
因此，合并式(2)与式(7)可得叠前地震ava反演的正问题方程为：
[0120]
d＝wr＝war＝fr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0121]
其中，f＝wa是子波褶积矩阵和系数矩阵a组合成的联合矩阵，即ava反演的正问题算子。
[0122]
步骤(3)：构建初始模型参数与先验低频模型参数。对实际测井资料或伪测井数据进行基线校正、环境校正、刻度均衡等预处理后，在地质解释层位约束下，由克里金插值法建立测井插值模型，作为反演求解的初始模型参数。再对测井插值模型进行低通滤波，得到先验低频模型参数。
[0123]
步骤(4)：基于正则化反演理论，构建格拉姆相关性约束叠前地震反演目标泛函。
[0124]
在叠前地震反演中，目前常用的目标泛函为：
[0125]
minf(m)＝||d-f(m)||2+αq(m)+μ||m-m
prior
||2,
ꢀꢀꢀꢀ
(9)
[0126]
其中，f则为叠前地震反演的正问题算子，m为叠前地震反演的所有模型参数，对于叠前地震ava反演，m即为式(8)中的r，q(m)为模型参数修正柯西约束，m
prior
为先验低频模型参数，最后一项为先验低频模型参数约束，α与μ为相应的约束系数。
[0127]
格拉姆行列式可作为不同模型参数间的非平稳相关性测度，其值越小则相关性越强。因此，将格拉姆相关性测度作为叠前地震反演目标泛函的约束条件，通过求解下面的反演目标泛函：
[0128][0129]
即为要求不同模型参数间具有非平稳相关关系。上式中，σ为格拉姆相关性测度约束的容限条件参数，则为模型参数格拉姆相关性测度，其表达式为：
[0130][0131]
其中，m1、m2、
…
、mn为m中包含的多个不同模型参数。对于叠前地震ava反演，则n为3，而m1、m2、m3则为纵波速度、横波速度与密度的相对变化率等，即r
p
、rs与r
ρ
，(m1,m2)为模型参数m1与m2的内积。因此，格拉姆相关性测度即为多个模型参数之间的内积组成的格拉姆矩阵的行列式，仅需要不同模型参数间的内积，不需要确定的模型参数相关关系，能够准确度量不同模型参数间的非平稳相关性。
[0132]
通过求解目标泛函(10)，可实现不同模型参数间的非平稳相关性约束，解决现有的相关性约束方法不能准确刻画非平稳相关性的问题。
[0133]
由正则化反演理论可知，目标泛函(10)中的格拉姆相关性测度约束本质上是一个正则化稳定性泛函。因此，目标泛函(10)可等价地写为如下形式：
[0134]
minf(m)＝||d-f(m)||2+αq(m)+μ||m-m
prior
||2+λgm,
ꢀꢀꢀꢀꢀ
(12)
[0135]
上式(12)即为格拉姆相关性约束叠前地震反演的目标泛函，其中，λ为格拉姆相关性约束系数。
[0136]
步骤(5)：基于格拉姆行列式的展开性质，对格拉姆相关性约束进行展开特性分析，计算格拉姆相关性约束的梯度方向。
[0137]
本项目采用梯度下降算法对目标泛函(12)进行迭代求解。由目标泛函(12)可知，在采用梯度下降算法进行迭代求解时，除了需要计算观测数据拟合项即目标泛函第一项、修正柯西约束项与先验低频模型参数约束项的梯度方向外，还需要计算格拉姆相关性约束项gm关于模型参数m的梯度方向。
[0138]
对于格拉姆相关性约束项的梯度方向，通过变分法进行求解。格拉姆相关性约束项的一阶变分为：
[0139][0140]
其中，(δmi,mi)g为格拉姆相关性测度关于模型参数mi的变分：
[0141][0142]
由格拉姆矩阵行列式的展开特性，将上式(14)写为：
[0143][0144]
其中，g
ij
是消除式(11)中格拉姆矩阵第i列和第j行形成的子格拉姆矩阵行列式。
[0145]
因此，格拉姆相关性约束项的一阶变分简化为：
[0146][0147]
其中，g(mi)即为格拉姆相关性约束项关于模型参数mi的梯度方向：
[0148][0149]
格拉姆相关性约束项的梯度方向是梯度下降算法的主要组成部分。
[0150]
步骤(6)：基于格拉姆相关性约束的梯度方向，采用梯度下降算法对反演目标函数进行迭代求解。
[0151]
输入当前模型参数，采用梯度下降算法对目标函数(12)进行迭代求解，其迭代公式为：
[0152][0153]
上式(18)中，τk是第k次迭代的迭代步长，由一维搜索法得到，mk为第k次迭代的当前模型参数，第1次迭代的当前模型参数为初始模型参数，m
k+1
为更新后的模型参数，为第k次迭代的梯度方向：
[0154][0155]
其中q为修正柯西约束关于mk的偏导数矩阵，上标t为矩阵转置，为格拉姆相关性约束项在mk的梯度方向。
[0156]
步骤(7)：判断反演目标函数是否满足迭代停止容限条件。
[0157]
若满足下面的迭代停止容限条件，则转步骤(9)：
[0158][0159]
其中ε是迭代容限参数。否则，转步骤(8)并进行下一次迭代。
[0160]
步骤(8)：若不满足迭代停止容限条件，通过无偏风险估计自适应求取格拉姆相关性约束系数，并进行下一次迭代。
[0161]
目标泛函(12)中，修正柯西约束与先验低频模型参数约束的正则化约束系数α与μ的求取方法有多种，例如l曲线法等，在反演求解前明确地给出。然而，由于格拉姆相关性测度的取值范围大，在反演目标泛函的迭代求解过程中，每一次迭代得到的模型参数，其格拉姆相关性测度的数值变化很大。不能采用目前成熟的方法，在反演求解前明确给出格拉姆相关性约束系数λ。
[0162]
本发明采用无偏风险估计法，在反演的迭代求解过程中，自适应求取格拉姆相关性约束系数。对于一个给定的格拉姆相关性约束系数，可通过梯度下降算法求解反演目标泛函(12)，得到对应于该约束系数的模型参数m
*
。对于不同的格拉姆相关性约束系数，反演得到的模型参数m
*
是不同的，因此模型参数m
*
是格拉姆相关性约束系数λ的函数：
[0163]m*
＝m
*
(λ).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0164]
无偏风险估计法通过求解模型参数m
*
的无偏风险估计泛函的极小化，在迭代求解过程中，自适应求取格拉姆相关性约束系数。无偏风险估计泛函为：
[0165][0166]
极小化无偏风险估计泛函的非线性虽较为严重，但求解它只是关于格拉姆相关性约束系数的单变量优化问题，本发明采用蒙特卡洛法对其进行求解。求出格拉姆相关性约束系数后λ，转步骤(6)进行下一次迭代。
[0167]
步骤(9)：若满足迭代停止容限条件，则停止迭代，输出反演得到的模型参数，并换算成最终的反演结果，即各种弹性参数、裂缝参数、各向异性参数、流体敏感参数等。
[0168]
对于叠前地震ava反演，模型参数为弹性参数的相对变化率。因此，需要通过弹性参数相对变化率与弹性参数间的积分关系，将弹性参数相对变化率换算成最终反演结果，即纵波速度、横波速度和密度等弹性参数。换算公式为：
[0169][0170][0171][0172]
式(23)至式(25)中，t为时间采样点，v
p
(0)、vs(0)和ρ(0)分别是初始纵波速度、横波速度和密度。