目前，流体因子在储层流体识别中具有非常重要的意义(Smith and Gildow，2000; Quakenbush et al.，2006; 印兴耀等，2013a)，而叠前地震反演是流体识别最有效的手段.传统的作法是首先反演纵波速度、横波速度和密度三个数据，然后用这三个数据间接地计算λρ、μρ、Russell流体因子等流体信息，这种作法的不足 之处是会带来累计误差.针对这一问题，很多学者开展了基于叠前反演的流体因子直接提取方法(宗兆云等，2012; 印兴耀等，2013b)，其优势在于能够减小累计误差，因此结果更加准确可靠.

在叠前反演方面，前人已经做过很多相应的研究.基于波动方程的叠前反演(Mora，1987; Shi et al.，2007)，在理论上比较成熟，但由于其正演模拟过程的复杂性和计算花费大，在实际中并没有广泛应用.基于Zoeppritz近似方程的叠前AVO反演(Hampson，1991;Simmons and Backus，1996; Bul and and Omre，2003)是目前发展最为迅速、深入，也是灵活性和效果最好的叠前反演方法.弹性波阻抗反演(Connolly，1999; Whitcombe，2002; 王保丽等，2005)简洁高效，是目前工业上应用最广泛的叠前反演方法，但是由于该方法使用的是部分叠加的角道集数据，损失了一部分的叠前AVO信息，因此反演结果的准确性方面不如叠前AVO反演.

在叠前流体识别方面，Smith和Gidlow(1987)等提出用加权叠加方法得到纵横波速度相对变化量，并从反演结果估算流体因子和检测气层的方法.Fatti等(1994)用权叠加的方法得到了纵波阻抗的相对变化量和横波阻抗的相对变化量.Goodway等(1997)通过反演得到λρ和μρ等数据体，用于描述储层的岩性和含流体性质，要比用纵波阻抗和横波阻抗更直观.Gray(2002)改进了Goodway方法，从叠前数据中反演出了λ和μ，消去了密度的影响，可以更好地描述岩性和流体.Russell等(2003)在前任研究的基础上，得到了流体饱和条件下的流体因子，称为Russell流体因子，并通过叠前反演的方法间接的获得了该流体因子.在国内方面，也有很多专家学者进行了叠前流体因子提取方面的研究(宁忠华等，2006; 郑静静等，2011; 宗兆云等，2012; 印兴耀等，2013b)，并取得了较好的应用效果.

本文是在前人(Russell et al.，2003; Downton，2005; 宗兆云等，2012)和笔者自己所做研究(Yin and Yang，2008; 杨培杰和印兴耀，2008)的基础上，提出了一种直接提取Gassmann流体因子的新方法.以贝叶斯反演框架为基础，并加入可变数量的点约束信息来稳定提取结果，来构造目标函数，最终实现了Gassmann流体因子的直接提取.该算法的核心是通过可变数量的点来约束提取的过程，可以使用两个点，也可以使用多个点进行约束，可以使用较精确的约束模型进行约束，也可以使用极低频约束模型进行约束，因此约束过程灵活，同时，该方法不需要初始模型的参与，由于没有初始模型的影响，因此该方法的提取结果更加客观准确.

模型试算表明，该方法对于地震数据中的噪声不敏感，即使在信噪比很低的情况下也能较好地提取流体因子.实际应用表明，流体因子提取结果分辨率高、稳定性好，对于不同流体的区分效果较好，并且与已知的流体解释结果能够很好地匹配，同时，该方法既可以使用叠前角道集，又可以使用叠前共成像点道集，适用范围较广，益于进一步推广应用.

将地球物理信息转化为储层信息属于反演的问题，从概率统计的角度，任何的反演问题可以看成是 一种贝叶斯估计问题(Ulrych et al.，2001; Tarantola，2005)，即，通过观测数据不断地更新先验知识，从而得到问题的解，公式为

其中，m 是待估计(反演)的参数空间，c 是归一化的常数，σ_post(m)是后验概率密度函数，p_prio(m)是先验概率密度函数，p_data[ d － f(m)] 是似然函数，d 是观测数据，f(m)表示正演算子.在本方法中 m 表示Gassmann流体因子.

对于式(1)的求解方法一般有两种(Bosch et al.，2010)，第一种是通过最优化目标函数(Yang et al.，2008)来求解，属于确定性反演方法；第二种是通过随机模拟(Haas and Dubrule，1994)的方法来实现，称为随机反演方法.本文采用第一种求解的思路，其优点是计算速度快、多解性小.

Russell等人基于Biot-Gassmann理论对饱含 流体多孔介质的AVO理论进行了研究，并推导了 包含Gassmann流体因子的反射系数近似公式(Russell et al.，2011)，公式为

其中，R_pp(θ)表示不同角度的地震反射系数，μ 表示剪切模量，ρ 表示密度，γ_dry2和γ_sat2 分别表示干燥岩石和饱和流体岩石的纵横波速度平方比平均值的平方，F_g 表示介质的Gassmann流体因子，其公式为(Russell et al.，2003; 张世鑫，2012)

其中，I_p、I_s 分别表示纵波阻抗和横波阻抗，c=γ_dry²，c 值可通过实验室测量的方法获得，如先计算干岩的泊松比，然后再得到c值(Russell et al.，2003)，或是在实际应用过程中通过反复试算来确定最佳c值，c的取值一般在2.25~3.0之间.

Gassmann流体因子可直接作为流体因子来识别储层流体类型，目前已有学者对不同流体因子的敏感性进行了定量的分析(Chi and Chi，2006; 张世鑫，2012)，认为Gassmann流体因子对于不同流体的识别最为敏感，无论是对气、油、水的区分度最高.式(2)和式(3)即为Gassmann流体因子直接提取的理论基础.需要指出的是，通过式(2)可提取出三个不同的参数，即Gassmann流体因子、剪切模量和密度，本文只考虑Gassmann流体因子的提取.

将公式(2)按偏移距写成矩阵形式为

其中，θ₁，θ₂，…，θ_K 表示分界面处的平均入射角度，R_f、R_μ、R_ρ 分别表示Gassmann流体因子反射系数、剪切模量反射系数和密度反射系数.

将上式扩展为

其中，F =[ R _f ]^T，N表示叠前地震数据的时间点数，K表示地震数据的覆盖次数，当使用地震角道集进行流体因子提取时，K一般取3.

文献(Downton，2005; 杨培杰，2008)指出，由于待反演的三个参数之间一般是统计相关的，因此需要应用三者之间的协方差矩阵对参数进行去相关处理，目的是提高解的稳定性，经过处理后的待反演参数之间是相互独立的.

公式(5)进一步变成：

其中，G ′ _NK×3N= G _NK×3N× V _3N×3N，V _3N×3N 为通过测井数据所获得的协方差矩阵(杨培杰，2008).公式(6)就是本方法中所使用的正演公式.

假定地震数据的噪声服从正态分布，且独立，F 服从柯西分布(杨培杰，2008)，用公式(1)将似然函数、先验分布以及正演公式结合起来，就得到了初步的目标函数为

其中，σ_n² 表示地震数据噪声的方差，σ² 为柯西分布中的一个参数，类似于高斯分布中的方差.

如果对式(7)的目标函数直接进行求解，由于该目标函数中没有考虑模型的约束作用，所以会导致反演结果的稳定性不好，特别是三维的反演结果，横向上的连续性会很差，因此，我们提出了一种新的约束思路，通过在提取过程中加入可变数量的点约束来达到稳定反演结果、提高反演结果的客观性、准确性的目的，公式推导过程如下.

在反射系数较小时有：

其中，r_g 为Gassmann流体因子反射系数.上式相对于时间积分，可以得到相对流体因子为

其中，F_g(t₀)为初始流体因子的数值.将上式简记为

其中，.

进一步地将式(10)进行扩展，就可以得到可变点约束的矩阵方程为

其中，n是约束点的数量，是变化的，m≥n>1，m是待反演的地震道的总点数，C 是1行m列的向量.将点约束矩阵方程简记为

那么根据贝叶斯公式，定义新的目标函数为

当点约束的矩阵中的n选不同的值时，就可以 实现待反演数据不同数量和不同时间点处的点约束.

式(13)即为最终的目标函数，下面将目标函数的每一项对参数求导，并令导数为零，最终可以得到式(14)的矩阵方程组为

其中，G ′^T G ′和G ′^T d 描述实际叠前道集与合成叠前道集间的拟合程度； θ Q 用来约束解的稀疏程度，即先验项，，λ₁ 为稀疏约束的权系数，ε 为实际叠前道集与反演合成的叠前道集间的误差，σ₁² 为协方差，Q 是一个斜对角加权矩阵； α C ′^T C ′和α C ′^T ξ 即为参数点约束项，α为点约束的权系数.通过求解 该方程组，就可以直接地提取出Gassmann流体因子.

首先构造砂泥岩互层模型(图 1)，然后基于Gassman理论，分别用气、油、水对砂岩进行流体替代，并计算砂岩含气、含油和含水后的纵波速度、横波速度和密度，如图 2所示，其中，气砂孔隙度 30%，含气20%；油砂孔隙度30%，含油90%；水砂孔隙度30%，含水90%.

图1

Fig. 1

图 1 砂泥岩薄互层模型 Fig. 1 Thin sand-shale interbed model

图2

Fig. 2

图 2 流体替代后砂泥岩薄互层模型 Fig. 2 Thin sand-shale interbed model after fluid substitution

用公式(3)计算流体替代后模型的Gassmann流体因子，如图 3所示，该模型数据被加入了一定的噪声，目的是使之更接近真实的地层参数.

图3

Fig. 3

图 3 模型的Gassmann流体因子 Fig. 3 Gassmann fluid factor of the model

以图 2中的纵横波速度、密度为输入，应用30 Hz 的混合相位子波，进行基于Zoeppritz方程的叠前地震正演，为了说明本文提出的方法优秀的抗噪性，在道集中加入信噪比2∶1的随机噪声，如图 4所示.

图4

Fig. 4

图 4 叠前合成记录(S/N=2) Fig. 4 Prestack synthetic record (S/N=2)

图 5为没有可变点约束的流体因子直接提取结果，由于叠前合成记录的信噪比很低，并且没有模型的约束，提取结果的稳定性和分辨率都比较差，虽然气层有较好的提取结果，但是对于油层，特别是水层的结果不理想.

图5

Fig. 5

图 5 无约束提取结果(S/N=2) Fig. 5 Unconstrained extraction results (S/N=2)

下面以该流体替代后的模型为例，来说明点约束的三种不同模式及其应用效果.

约束模式1：约束点的数量n=2，此时，认为待反演地震道的第一个点和最后一个点的流体因子值是已知的，其他所有的点都是未知的.α的取值应大一些，一般100>α>50.此时的可变点约束的矩阵方程如式(15)所示，其点约束的效果如图 6所示，公式(15)为

图6

Fig. 6

图 6 模式1提取结果(S/N=2) Fig. 6 Extraction result with pattern 1

图7

Fig. 7

图 7 模式2提取结果(S/N=2) Fig. 7 Extraction result with pattern 2

约束模式3：n=m，即使用和待反演的地震道相同的点数来进行约束，这时不使用高频的初始模型，而是往往使用一个极低频的模型(极端情况下是使用一个常数模型进行约束)，此时α的取值应最小，一般1>α，则反演结果不会和约束模型很像，而是往模型靠拢，即能够起到稳定反演结果的目的.此时的可变点约束的矩阵方程如式(17)所示，其点约束的效果如图 8所示，公式(17)为

图8

Fig. 8

图 8 模式3提取结果(S/N=2) Fig. 8 Extraction result with pattern 3

为了便于理解和应用，对不同的可变点约束模式与原则进行了总结，如表 1所示.

表1(Table 1)

表 1 约束模式与原则 Table 1 Constraint patterns and principles 模式 约束点数量(n) 约束点权重α) 约束原则 模式1 n=2，即待反演地 震道的第一个点 和最后一个点 1000>α>50 使用精确的模型，α的 取值应大一些，则可 以将结果紧紧地约束 在待约束点上 模式2 m>n>2，即待反 演地震道的第一 个点和最后一个 点，再根据实际情 况选择若干个点 50>α>1 使用较精确的模型，α 的取值应较大，则可 以将结果约束在待约 束点上 模式3 n=m，即使用和 待反演的地震 道相同的点数来进 行约束 1>α 使用极低频模型或是 常数模型，α的取值应 小，能够起到稳定反 演结果的目的

表 1 约束模式与原则 Table 1 Constraint patterns and principles

下面用一个二维Gassmann流体因子模型来说明不同约束模式的适用条件和应用效果，如图 9所示，并设置一口虚拟井Well A，用于提供约束模型.图 10是首先根据二维Gassmann模型的纵波速度、横波速度、密度进行叠前正演，得到叠前道集后，再对叠前道集进行叠加的剖面，在该剖面上解释了5个层位，用来提供约束点信息，从深到浅分别是h₀、h₁、h₂、h₃、h₄.此次模型反演使用的是三个角度的角道集数据，信噪比为4∶1.

图9

Fig. 9

图 9 二维Gassmann流体因子模型 Fig. 9 2D Gassmann fluid factor model

图10

Fig. 10

图 10 二维Gassmann模型叠加剖面(S/N=4) Fig. 10 Stacked profile of 2D Gassmann model(S/N=4)

图 11为模式1的Gassmann流体因子直接提取结果.在该模式中，分别选取层位h₁和h₂作为两个约束点(在二维应用中，约束点就变成了约束线)，两个约束点(线)的流体因子值从模型中得到，通过反复试验，设定约束点权重α=200.可以看出，在该模式中，由于没有初始模型的作用，因此流体因子的提取结果客观准确.

图11

Fig. 11

图 11 二维模式1提取结果 Fig. 11 2D Extraction result with pattern 1

图 12为模式2的Gassmann流体因子直接提取结果.在该模式中，分别选取层位h₀、 h₁、 h₂、 h₃、 h₄作为约束点(线)，约束点的流体因子值从模型中得到，通过反复试验，设定约束点权重α=20.该模式由于也没有初始模型的作用，因此流体因子的提取结果较为客观准确.

图12

Fig. 12

图 12 二维模式2提取结果 Fig. 12 2D Extraction result with pattern 2

图 13b为模式3的Gassmann流体因子直接提取结果.在该模式中，使用一个极低频模型的所有点 进行约束，该极低频模型来自于对原始模型的400次平滑处理，如图 13a所示，通过反复试验，设定约束点权重α=0.02.可以看出，相对于模式1和模式2，模式3的流体因子提取结果分辨率和准确度略有降低，但是横向上的稳定性却得到了有效的提高.同时需要指出的是，在模式3中，低频约束模型的作用只是为了稳定反演结果，并不参与反演解的迭代过程，因此，不会影响到结果的客观准确性.

图13

Fig. 13

图 13 二维模式3提取结果 (a) 约束所用的极低频模型；(b)提取结果. Fig. 13 2D Extraction result with pattern 3 (a) Extremely low frequency model；(b) Extraction result.

从提取效果可以看出，该方法类似于基于模型的方法，但是与模型反演又是不相同的，基于模型的反演方法受模型的影响比较大，且分辨率会比较低，而该方法的灵活度高，且受模型影响的程度小，其优点是可以控制模型对反演结果的影响，在稳定性和分辨率之间寻求一个平衡.

一般来说，模式1适合于约束层位间时间段比较小的情况，即属于一种针对目标体的提取模式.实际应用中总结发现，当反演时间段大于一个子波的长度，而小于两个子波的长度时，应用模式1会有较好的效果.

模式2可以看成是模式1的不同组合形式，即模式2可以用多个不同的层位作为约束点，进行约束，该方法适用于解释的层位比较多的研究工区.

模式3使用约束模型的所有点进行约束，模式3更适合于约束层位间时间段比较大的、信噪比低的地震资料的流体因子提取.同时，如果反演的数据段不是很长，可以使用一个常数模型进行点约束，也会有很好的反演结果，即认为在这个时间段内地层的流体因子的变化趋势是个常数.

在实际应用中，对于点约束的模式怎么选取并没有什么规则可言，可尝试不同的约束模式进行约束，并分析不同模式的提取结果和效果.最后需要指出的是，这三种约束模式是最常用的，但并不是唯一的，也可以根据工区的实际情况来制定不同的约束模式，此时只需要改变点约束的矩阵方程即可，或是将这三种模式进行组合使用，同样会有很好的效果.

流体因子直接提取实际应用来自胜利油田埕岛地区，工区主力含油层系是新近系，是目前勘探开发的热点地区.

叠前数据使用的是部分叠加角道集数据，以Nm41和T0作为控制层位，上下各开50 ms进行反演.通过分析，认为该区的地震资料的信噪比较低并且连续性较差，且反演的时间段较长(约240 ms)，所以采用约束模式3进行约束，由于考虑到待反演的时间段只有200 ms左右，可以近似的将该时间段内的流体因子值的变化趋势看成是一个直线，根据 工区内的测井资料计算可知，工区的流体因子平均 值在1.05×10⁷左右，即提取过程使用一个值为1.05×10⁷的常数模型进行约束.同时，根据工区测井资料和岩石物理的分析结果，并不断的反复试 算，认为当 γ_dry²=2.25，γ_sat²=4.0，λ₁=0.65，α=0.1 时的流体因子提取结果最为理想.图 14为Gassmann流体因子提取结果过A井和B井的一条连井剖面.

图14

Fig. 14

图 14 Gassmann流体因子提取剖面 (a) 叠后地震剖面；(b) 直接提取结果；(c) 间接提取结果. Fig. 14 Gassmann fluid factor extraction profiles (a) Poststack seismic section；(b) Direct extraction；(c) Indirect extraction.

进一步将流体因子直接提取的结果和间接求取的结果进行了对比分析.图 14a为叠后地震剖面，黑色曲线为伽玛测井曲线，粉红色曲线为油水指示曲线，地震反射属于中强度振幅，信噪比较低，很难通过叠后亮点或是分频技术有效的识别含油气性；图 14b为Gassmann流体因子直接提取的结果，图 14c为间接提取的结果，其中白色曲线为伽玛曲线.对于A井，井旁道1.06 s和1.077 s处分别有一个6 m、10 m的油层，1.11 s处有一水层，图 14b的流体因子直接提取结果在这三处都表现为低值，并且油层的流体因子相对于水层的流体因子更低一些.而图 14c的流体因子间接提取结果对于6 m油层的反映并不明显，对于10 m的水层没有反映.

对于B井，井旁道1.055 s处、1.075 s处、1.105 s处分别有一个11 m、9 m、7 m的水层，图 14b的流体因子直接提取结果对于这三处都能较好的识别，而图 14c的流体因子间接提取结果只是在1.06 s处有个很低的值，对于上面的两个水层，也没有分开，分辨率明显低于流体因子直接提取的结果.

通过以上的对比可以看出，流体因子在油层、水层处都表现为低值，但是油层的流体因子相对于水层的流体因子会更低一些.

通过提取两个数据体的沿层结果进行比较，来进一步对比和分析本文所提出的方法的有效性和先进性.如图 15所示，图 15a为本文方法的流体因子直接提取结果，图 15b为流体因子间接提取的结果，沿层结果为拾取的Nm41层位往下40 ms的流体因子均方根振幅值.本工区目的层系为河流相储层，流体因子直接提取的结果(图a)很好地展现了河流相砂体储层的形态，准确度高.然而，相比于流体因子直接提取的结果，流体因子间接提取结果(图b)并没有较好地展现河道砂体储层的形态，特别是B井在沿层切片上也表现出很低的流体因子值，与实钻结果不吻合，准确度不高.

图15

Fig. 15

图 15 Gassmann流体因子提取结果的沿层切片 (a) 直接提取结果；(b) 间接提取结果. Fig. 15 Layer-along slices of extracted Gassmann fluid factor (a) Direct extraction；(b) Indirect extraction.

从实际应用结果可以看出，相对于间接提取的方法，本文提出的Gassmann流体因子直接提取结果的分辨率更高、准确性更好，从而实现了不同流体的有效区分，并且能够与已知的流体解释结果很好地匹配，说明本文提出的算法是准确可靠的.

(1)本文提出的可变点约束Gassmann流体因子直接提取算法约束过程灵活，可以使用两个点，也可以使用多个点进行约束，可以使用较精确的模型进行约束，也可以使用极低频模型进行约束，并可以通过调节α的大小来控制点约束的权重.

(2)由于最大程度的弱化了对于反演初始模型的依赖，因此提取过程更加客观准确，并且提取结果的稳定性好，实际应用范围也更加广泛，益于进一步的推广应用.

(3)该方法不仅可以提取Gassmann流体因子，也可以提取其它类型的流体因子，此时只需得到该流体因子的近似方程，而可变点约束的模式和原则不用改变.