0 引言

AVO反演是地球物理油气勘探领域进行流体识别的有效方法，AVO以Zoeppritz(1919)方程为理论基础来计算反射界面上下的弹性参数(Chopra and Castagna, 2014).但是受限于Zoeppritz方程复杂的形式，人们很难从振幅变化中直接获得弹性参数信息，也无法直接描述地层参数的变化对反射系数的影响，于是出现了多种Zoepprtitz的近似方程，这些近似方程给出了不同弹性参数对反射系数变化的直观表示(Aki and Richards, 1980; Shuey, 1985; Smith and Gidlow, 1987; Fatti et al., 1994)，且计算简单，促进了在地下弹性参数估计方面的发展.

从Ostrander (1984)首次将AVO反演技术成功运用到含气砂岩的检测中开始，许多研究者基于纵波反射地震数据开始了AVO反演的研究(Smith and Gidlow, 1987; Buland and More, 2003a; Downton, 2005; 彭真明等，2008；Ursenbach and Stewart, 2008; Alemie and Sacchi, 2011；张丰麒等，2013)，并且取得了不同程度的进展.然而，纵波反射系数对横波速度和密度是不够敏感的，在有噪声的情况下，单纯利用纵波反射数据进行反演并不能得到可靠的结果，如果没有准确的横波速度和密度信息，可能会导致后续储层的错误解释.Jin等(2000)指出，转换横波地震数据可以提供更准确地横波速度和密度估计.随后，Kurt (2007)验证了基于纵波和转换横波地震数据的联合反演可以得到更好的弹性参数估计.随着纵波激发、三分量接收技术的发展，大大降低了利用转换横波勘探的成本，也促进了PP-PS联合反演的发展(傅戈平, 2014).Buland和Omre(2003b)通过马尔可夫蒙特卡洛方法进行了联合反演；Rabben等(2008)对Zoeppritz方程进行了二阶近似，引进Metropolis-Hastings算法进行了联合反演；Kato(2010)通过合成数据和Hangingstone油田实际数据论证了联合反演的优越性.Lu等(2015)开展了基于精确Zoeppritz方程的最小二乘联合反演.由于地震资料质量不高、子波提取不准等问题，地震反演本身就是一个不适定的问题，对误差极为敏感，联合反演在一定程度上降低了解的不确定性，文中我们还通过引入正则化技术来解决不适定性问题.不适定性是球物理反演问题内在的基本属性，除非附加的先验信息或某种知识，比如单调性、光滑性、边界条件或数据的误差界等进入反演模式中，否则反演的困难并不会降低(王彦飞等，2011).在地震反演中，许多参数本身是具有稀疏性质的，如地震反射系数，用传统的l₂范数极小化模型误差求解是不合适的，因此，引入求解最优解的新观点—求解l₁极小化问题.基于l₁范数的极小化问题能够获得问题的稀疏解，在求解该问题时，一类常用的算法为匹配追踪法(Davis et al., 1997；Donoho and Tsaig, 2006)，当求解的问题十分稀疏时，该方法优势明显，但是存在计算量大的问题；另外可以利用内点法(Lustig et al., 1991; 姜志侠, 2008; 王彦飞等, 2011; 孙晨, 2012)进行求解.该方法的代表性工作是Karmarkar (1984)提出的Karmarkar线性规划算法，该算法在多项式时间内通过遍历可行域内部来寻找最佳解.Nemirovskii和Nesterov(1994)将内点法推广到凸规划问题中，通过在可行域边界设置自适应障碍函数，不断缩小凸集范围来获得问题的解.Forsgren等(1996)、Forsgren和Gill(1998)研究了不同内点法中矩阵分解的稳定性问题.由于内点法通常需要求解对偶问题，而求解对偶问题时，原问题和对偶问题的间隙(Dual gap)又很难控制.因此，本文转而考虑l₁范数极小化问题的梯度投影算法.

由于观测系统的限制，地球物理反演问题多是离散、稀疏的.因此，如果一个算法不具有强适应性，先验约束的l₁空间极小化解可能比传统的正则化效果好(王彦飞等，2011).考虑到反射系数的稀疏性，引入l₁范数正则化模型，同时引入梯度投影解法(Wang and Ma, 2007)来对模型进行求解.

为了加快梯度投影算法收敛的速度，同时避免陷入局部极小值的情况，本文引入了粒子滤波算法来为梯度投影算法提供初始值.随着计算机计算性能的不断增强和计算成本的降低，粒子滤波(Particle Filter, PF)(Gordon et al., 1993)逐渐成为解决当前非线性，非高斯系统的重要方法.粒子滤波直接采用样本(粒子)形式而非函数的形式对系统概率分布进行描述，不需要对状态变量的概率分布做过多假设.由于粒子滤波在解决非线性，非高斯问题上面的优势，国内外学者已经将其运用到了各个领域，如计算机视觉、目标追踪、通讯与信号处理等领域(王法胜等，2014).据我们调研的文献所知，Baziw(2005)可能是最早将粒子滤波运用到地球物理领域的，他利用Rao-Blackwellized粒子滤波来预测低信噪比情况下的地震事件.Peng等(2008)将粒子滤波引入到AVO反演中，并且证明了该方法在反演中的有效性.Tang和Wang(2017)将粒子滤波引入贝叶斯反演框架中，实现了AVO参数联合反演.

本文工作的新颖之处在于将粒子滤波引入到稀疏反演问题的梯度投影下降法中，由于粒子滤波属于递归贝叶斯滤波方法的一种，它结合先验模型和观测数据的影响产生的后验分布能够给稀疏反演算法提供一个较好的可行解，加快收敛速度的同时也提高了解的精度.

1 方法理论 1.1 正演理论

Aki和Richards(1980)基于相邻两侧介质参数变化较小的假设，在前人的基础上得出了以纵横波速度以及密度相对变化作为变量的三参数模型，公式为

(1)

其中：

(2)

此外，V_P、V_S和ρ分别为上下层纵波速度、横波速度和密度的平均值；ΔV_P、ΔV_S和Δρ为上下层介质之差；θ为上下介质纵波的平均入射角；φ为上下介质转换横波的平均入射角；γ为V_S/V_P的比值.

Stolt和Weglein(1985)从适用于单一均匀界面的Aki-Richards近似公式出发，推导出时间连续方程为

(3)

其中，a_p(t, θ)、a_S(t, θ)、a_ρ(t, θ)、b_S(t, θ)、b_ρ(t, θ)是公式(2)随时间变化的概括表达.通过地震褶积模型，可得到含噪声的AVO联合反演形式为

(4)

其中，关于观测数据d和矩阵算子G的具体表达式见(Tang and Wang, 2017).在正演形式已知的情况，利用叠前地震数据可直接反演出纵波速度、横波速度以及密度项.

1.2 稀疏反演理论

我们考虑AVO反演中模型参数的稀疏性和非负性，用l₁范数对问题进行求解时，该问题具有形式为(Wang et al., 2009)：

(5)

其中S为待求参数的解的集合.求解问题(5)的一种方式是建立一个对偶问题.但求解对偶问题时，原问题和对偶问题的间隙很难控制，同时用内点法进行求解时，要求初始值尽量靠近中心路径并接近真实的解，这在实际应用中很难达到.为此，本文研究一种非线性规划方法，利用Lagrange乘子法，可以将上面(5)式转化为有界约束的二次规划问题进行求解，公式为

(6)

其中e为所有组份都为1的向量，λ为Lagrange乘子.我们采用梯度投影方法对(6)式进行求解(Figueiredo et al., 2007；Wang and Ma, 2007).设一个迭代算法的第k次迭代点为m^(k)，目标函数的梯度值为g^(k)=∇J(m^(k))，则迭代公式可以写成：m^(k+1)=m^(k)－τ^(k)g^(k)，其中τ^(k)为迭代步长.由于(6)的解集合为S，则梯度投影公式可以定义为(Wang and Ma, 2007)：m^(k+1)=P_S(m^(k)－τ^(k)g^(k))，其中P_S(·)为投影算子，定义作：P_S(·)=(·)₊=max{0, ·}.

当我们进行迭代求解时，沿着负梯度方向进行搜索，并执行回溯线性搜索来确保迭代点的可行性，其基本思想是从一列步长点列中选择满足非精确线搜索条件的一个最优步长.初始步长由极小化问题τ₀=argmin_τJ(m^(k)－τg^(k))获得，即根据(7)式来却确定初始步长，公式为

(7)

总结上面的推导过程，可以得到该极小化问题的梯度投影算法步骤如算法1所示.

算法1  AVO参数反演的稀疏解算法

1) 初始化：

(a) 选择最大迭代步数k_max，β∈(0, 1), μ∈(0, )，令k=0;

(b) 令m₀为方程Gm=d的一个初始解(该解由1.3节粒子滤波方法得到).

2) 根据公式(7)计算τ₀的值.

3) 回溯线性搜索，在序列{τ₀, βτ₀, β²τ₀, …}中第一个满足以下不等式的值赋为τ^(k)：

其中：

4) 进行终止条件判断，达到停机要求则输出，否则令k←k+1，返回第2)步继续进行迭代计算.

(1) 在算法1中，要求初始值必须满足非负约束条件.为了保证每一步迭代过程中迭代点列的可行性，引入了投影技巧(由公式(·)₊=max{0, ·}保证)，这样算法的收敛性就有了保障.可以证明，当dist(m, S)≤κ‖min(m, g(m))‖时(κ>0为一常数)，上述梯度投影算法是收敛的.停机准则我们采用的是：‖min(m, g(m))‖≤ε，其中ε>0是一个小的正数.‖min(m, g(m))‖≤ε意味着dist(m, S)足够小，这里dist(·, ·)为距离算子.dist(m^(k), S)表示的是迭代点列距解的可行集S的距离.(2)为了加快收敛速度，需要为方程Gm=d提供一个合适的初始解m₀，本文提出获得m₀的粒子滤波算法.

1.3 粒子滤波

粒子滤波是一种递归贝叶斯蒙特卡罗方法，它的基本思想是对状态向量的经验条件分布进行采样，并利用所观察到的数据不断更新样本的位置和权重，最终修正后的条件分布将收敛到真实状态向量的条件分布.根据我们的具体问题，我们考虑式(8)的离散时间状态空间模型，公式为

(8)

其中，m_k∈R^3n×1为k时刻的状态向量，y_k∈R^2n×1相应的测量值.h(·, ·):R^3n×1×R^2n×1→R^2n×1为状态向量的更新方程，h(·, ·):R^3n×1×R^2n×1→R^2n×1描述了状态向量和观测值之间的关系.μ_k和v_k分别为k时刻的状态扰动和观测噪声.我们假设入射角θ随时间变化，则可以得到相应的时间维度.其中，离散时间序列可由垂直入射角到临界入射度代替.

1.3.1 蒙特卡罗积分

当选择最小方差作为损失函数时，状态估计问题退化为计算后验分布的期望值，可以写为

(9)

如果这个问题可以描述为线性高斯状态空间模型，采用卡尔曼滤波(Kalman, 1960)可以直接得到该问题的解析解.由于本文特定的反演问题，观测值和入射角呈非线性关系.因此，对于公式(9)并没有明确的表达形式.蒙特卡罗方法将积分看成某种随机变量的数学期望，并通过采样求和的方法加以估计，其实质是用一组采样值去近似目标分布.中央极限定理确保了该方法的收敛性，收敛的误差阶为O(1/)(Särkkä, 2013).根据后验概率密度分布进行采样，即m⁽ⁱ⁾~p(m|d_1:T), i=1, …N，则上面的式子可以写为

(10)

1.3.2 重要性采样

在实际模型中，我们无法直接从p(m|d_1:T)中进行采样，于是我们采用重要性采样的方法，对易于采样的近似分布π(x)进行采样，并且每一个粒子{m_k¹, …, m_k^N}都具有相应的权值{w_k¹, …, w_k^N}.权值是对每个采样值跟真值之间近似情况的一种度量.上述所求的积分值可以通过以下加权求和的形式获得，公式为

(11)

其中：

(12)

w_kⁱ为均一化的重要性权值.

总结粒子滤波的基本思想，我们得出了以下标准粒子滤波的算法步骤如算法2.在应用算法2时，我们给出一些说明：定义为重采样后的粒子和粒子权值；定义N_eff为有效粒子容量.迭代公式(8)的具体形式为

(13)

其中G(k)为第k个入射角时的正演算子.

算法2  粒子滤波算法

1) 初始化.由先验概率p(m₀)产生N个初始粒子，且权值均为1/N；给定有效粒子数阈值N_thr.给定模型的初始分布为

2) 对于k=1, …K，选取合适的重要性密度函数，抽取N个样本：

3) 计算权值，假设重要性分布具有马尔科夫性，公式为

4) 均一化权值，公式为

5) 当N_eff有效粒子容量较低时，当N_eff < N_thr进行重采样.其中：

6) 计算粒子滤波估计值，公式为

7) k:=k+1, 返回第1步.

2 数值实验

为了检验算法的正确性，同时更好的模拟实际数据，本文建立了不同信噪比情况下的地震响应模型.正演计算采用Stolt-Weglein连续时间方程：d=Gm=WADm+e, 其中W为卷积算子，e为噪声，W、A和D分别表示为下面的形式，公式为

(14)

其中，dt为时间导数算子.有关上述角道集生成的具体实现见(Buland and More, 2003a).

Buland和More(2003a)阐述了在界面弹性参数变化较小情况下，采用该近似方程的合理性.在进行PP-PS联合反演前，首先需要将转换横波时间剖面与纵波剖面上来自同一反射层的同相轴进行匹配.这就需要将PS道集压缩至PP波时间后再提取PS地震子波.在我们的数值模拟中，我们分别利用主频为50 Hz和35 Hz的理论雷克子波合成入射角在0~30°内的PP和PS角道集，采样率为2 ms.理论模型的P波速度，转换S波速度和密度如图 1所示，其中红色曲线代表了真实的模型值.该模型是不光滑的，对于传统的l₂范数下的反演方法会无法反演出真实的结果.通过地震子波和反射系数的褶积我们获得PP和PS的合成记录，为了探讨不同信噪比情况下反演算法的可行性，我们模拟了三种随机噪声水平的地震数据，即：S/N=100；S/N=10；S/N=5(其中，S/N表示信噪比).具有不同噪声水平的合成数据如图 2所示.

在稀疏反演算法中，其初始值采用粒子滤波的结果，反演结果如图 3所示.其中红色曲线为稀疏反演方法的反演结果，蓝色曲线为传统贝叶斯反演结果，黑色为真实模型.从图 3的对比结果，我们可以看出，新方法在不同的信噪比资料中都能取得较好的效果，并且在低信噪比情况下，相比传统反演方法能获得更多参数的细节变化.

3 结论与认识

针对油气勘探领域流体识别的AVO反演问题，本文建立了基于l₁范数优化的稀疏反演模型，给出了求解反演问题的稀疏优化算法，并将粒子滤波方法与稀疏优化算法结合起来，通过模型试验证明了该方法的有效性，同时通过与传统反演方法相对比，发现在信噪比较低情况下，稀疏反演方法具有更好地刻化模型的能力.

致谢

感谢两位审稿人和编辑提出的宝贵意见，使得论文内容更加充实.