储层预测技术

4.1 LPM 储层预测技术
LPM 是斯伦贝谢公司GeoFrame 地震解释系统中最新推出的储层预测软件，利用地震属性体来指导储层参数（如砂岩厚度）在平面的展布，以此来实现储层参数的准确预测。

LPM 预测储层砂体可分两步进行：首先，它是将提取的地震属性特征参数与井孔处的砂岩厚度、有效厚度进行数据分析，将对储层预测起关键作用的地震属性特征参数优选出来，根据线性相关程度的大小，建立线性或非线性方程。

线性方程的建立主要采用多元线性回归方法；非线性方程的建立主要采用神经网络方法；其次，根据建立的方程，利用网格化的地震属性体来指导储层参数（如砂岩厚度）在平面的成图。

4.1.1多元线性回归基本原理
设因变量y 与自变量x 1， x 2 ，…，x m 有线性关系，那么建立y 的ｍ元线性回归模型：
ξβββ++++=m m x x y 110 (4.1)
其中β0，β1，…，βm 为回归系数；ξ是遵从正态分布N(0，σ2)的随机误差。

在实际问题中，对y 与x 1， x 2 ，…，x m 作n 次观测，即x 1t ， x 2t ，…，x mt ，即有：
t m t m t t x x y ξβββ++++= 110 (4.2)
建立多元回归方程的基本方法是：
(1)由观测值确定回归系数β0，β1，…，βm 的估计b 0，b 1， …，b m 得到y t 对x 1t ，x 2t ，…，x mt ；的线性回归方程：
t m t m t t e x x y ++++=βββ 110 (4.3)
其中t y 表示t y 的估计；t e 是误差估计或称为残差。

(2)对回归效果进行统计检验。

(3)利用回归方程进行预报。

回归系数的最小二乘法估计
根据最小二乘法，要选择这样的回归系数b 0，b 1， …，b m 使
∑∑∑===----=-==n
t n t mt m t t t t n t t
x b x b b y y y e Q 11211012
)()( (4.4) 达到极小。

为此，将Q 分别对b 0，b 1， …，b m 求偏导数，并令
0=∂∂b
Q ，经化简整理可以得到b 0，b 1， …，b m ，必须满足下列正规方程组：
⎪⎪⎩
⎪⎪⎨⎧=+++=+++=+++my m mm m m y m m y m m S b S b S b S S b S b S b S S b S b S b S
22112222212111212111 (4.5) m m x b x b x b y b ----= 22110 (4.6)
其中
∑==n
t t y n y 1
1 (4.7) m i x n x n
t it i ,,2,111
==∑= (4.8) ),,2,1())((1))((11
11m i x x n x x x x x x S S n
t n t jt it jt n t it j jt i n t it ji ij =-=--==∑∑∑∑==== (4.9) ),,2,1())((1))((1111m i y x n y x y y x x S n
t n t t it n t t it t i n t it iy =-=--=∑∑∑∑==== (4.1
0)
解线性方程组(4.5)，即可求得回归系数i b ，将i b 代入(4.6)式可求出常数项0b 。

4.1.2 BP 网络网络基本原理
多层感知器具有独特的学习算法，该学习算法就是著名的BP 算法，所以多层感知器常常被称为BP 网络。

BP 网络是一种层状结构的前馈神经网络，它是由输入层、隐含层（一层或者多层）和输出层构成（图4-3）；输入层神经元的个数为输入信号的维数，隐含层个数视具体情况而定，输出层神经元个数为输入信号维数。

BP 神经网络输入层中的每个源节点的激励模式（输入向量）单元组成了应用于第二层（如第一隐层）中神经元（计算节点）的输入信号，第二层输出信号称为第三层的输入，其余层类似。

网络每一层神经元只含有作为他们输入前一层的输出信号，网络输出层（终止层）神经元的输出信号组成了对网络中输入信号（起始层）源节点产生激励模式的全部响应。

即信号从输入层输入，经隐层传给输入层，由输出层得到输出信号。

BP 学习过程可以描述如下：
工作信号正向传播：输入信号从输入层经隐单元传向输出层，在输出端产生输出信号，这是工作信号的正向传播。

在信号向前传递过程中网络的权值是固定不变的，每一层神经元的状态只影响下一层神经元的状态。

如果输出层不能得到期望的输出，则转入误差信号的反向传播。

误差信号的反向传播：网络实际输出与期望输出之间的差值即为误差信号，误差信号由输出端开始逐层向前传播，这是误差信号的反向传播。

在误差信号反向传播过程中，网络权值由误差反馈进行调节。

经过权值的不断修正使网络实际输出更接近期望输出。

误差反传播算法（BP 算法）利用梯度下降技术使实际输出y (t)与期望输出d (t)的误差能量最小。

∑-=2/)(2)()()(t t t y d e （4.11）
网络学习时，开始取一小的随机数作为网络权值和内部阈值的初值，然后反复输入训练样本，计算实际输出与期望输出的差值，据此调整权值，直至权值收敛，并使代价函数降至可接受值。

对训练样本集中第P 个输入及其期望输出，网络的第j 个节点与第i 个节点的联接权修正量p ji W ∆可用下式计算
ij pj pj W O ηδ∆=∙∙ （4.12）
式中 η—学习率
pi δ—误差项
图4-1 BP 人工神经元模型
pi O —节点i 的输出
其中误差信号δ对输出层和隐含层分别为：
)()1()2()()()2()()2()(2t t t t t O d O O -∙-∙=）（δ （4.13）
∑-=++∙-∙=1
0)1()1()()()()()()1(k n i k ji
k t i k t j k t j k t W O O δδ）（ （4.14） 调整后的权值与阈值分别为：
)()()()()()1(k t ij k t ij k t ij W W W ∆+=+ （4.15）
)(
)()(
)()(
)1(k t j k t j k t j δηθθ∙-=+ （4.16）
节点输出O 由前向传播算得：
)()0(
)(t j t j x O = （4.17）
()⎪⎪⎭
⎫ ⎝⎛-∙=∑-=-10)(1)()
1()(k n i k j k i k ij t j O W f O θ （4.18） (2)
()()t j t y O = （4.19）
以上各式中k ＝0，1，2代表输入层、第一隐层和输出层；N k 为第k 层节点个数；j=0，1，…，（N k -1）。

BP 算法的步骤可归纳为
第一步 设置变量和参量：
X k =[x k1，x k2，…，x kM ]，(k =1，2， …，N )为输入向量，或称训练样本，N 为训练样本个数。

Y k (n)=[y k1(n)， y k2(n)，…，y kM (n)]，(k =1，2， …，N )为第n 次迭代时网络的输出。

η为学习率
n 为迭代次数
第二步 初始化，赋给初始权值和初始阈值较小的随机非零值。

第三步 随机输入样本X k ，n ＝0。

第四步 对输入样本前向计算BP 网络每层神经元的输入信号u 和输出信号o 。

第五步 由期望输出d k 和上一步求得的实际输出Y k (n )计算误差E (n )，判断是否满足要求，若满足转至第八步；不满足转至第六步。

第六步 判断n +1是否大于最大迭代次数，若大于转至第八步，若不大于，对输入样本X k ，反向计算每层神经元的局部梯度δ。

第七步 计算权值修正量Δw ，并修正权值；n =n +1，转至第四步。

第八步判断是否学完所有的训练样本，是则结束，否则转至第三步。

上述BP学习过程中要注意几点：
（1）BP学习时权值的初始值是很重要的。

初始值过大，过小都会影响学习速度，因此权值的初始值应选为均匀分布的小数经验值，大概为（-2.4/F，2.4/F）之间（也有人建议在（
F为所连单元的输入端的个
数，另外，为避免每一步的权值的调整方向是同向的（即权值同时增加或同时减小），应将初始权值设为随机数。

（2）神经元的激励函数是Sigmoid函数，如果Sigmoid函数的渐近值为+α和-α，则期望值只能趋近于+α和-α，而不能达到+α和-α。

为避免学习算法不收敛，提高学习速度，应设期望输出为相应的小数，若逻辑函数渐进值为1和0，此时设相应的期望输出为0.99和0.01等小数，而不应设为1和0。

（3）用BP算法训练网络时有两种方式：一种是顺序方式，即每输入一个训练样本修改依次权值；另一种是批处理方式，即待组成训练周期的全部样本都一次输入网络后，以总的平均误差能量为学习目标函数的修正值的训练方式。

顺序方式所需的临时存储空间较批处理方式小，但顺序方式的误差收敛条件难以建立，而批处理方式能够精确的计算出梯度向量，收敛条件非常简单，易于并行处理。

（4）BP学习中，学习步长η的选择比较重要。

η值大权值变换就大，则BP学习的收敛速度就快，但是η值过大引起振荡即网络不稳定；η值小可避免网络不稳定，但是收敛速度就慢了。

要解决这一矛盾最简单的方法就是加入“动量项”。

（5）要计算多层感知器局部梯度δ，需要知道神经元的激励函数f(·)的导数。

（6）在BP算法第五步需要判断误差E(n)是否满足要求，这里的要求是：对顺序方式，误差小于我们的设定值，即E(n)<ε；批处理方式，每个周期的平均误差其变E av化量在0.1%到1%之间，我们就认为满足误差要求了。

（7）在分类问题中，我们会碰到同一类的的训练样本有几组，在第一步设置变量时，一般使用同一类的训练样本其期望输出相同。

4.1.3 LPM储层预测步骤
1、油层标定
为了分析储层砂体在地震剖面上的反射特征，必须将油层在地震剖面上进行准确的标定。

油层的标定是岩性解释的基础，标定是否准确直接影响到岩性预测的精度。

本次资料解释工作共预测4个油层组，即扶Ｉ油层、扶Ｉ油层上、中、下部油层的砂岩厚度。

这四个油层分别对应于地震时间剖面上的T2层至T2层+45ms时窗内的反射
波。

以T2为基点，将13口井的合成地震记录与时间剖面进行对比，从而确定四个油层组在时间剖面上的反射位置（见图4-2）。

最后，我们根据各油层在时间剖面上占的范围确定各油层岩性预测的时窗长度，使时窗长度既包括了各油层的完整波形，又最大限度地减少了相邻油层波形的进入，保证了砂岩预测的准确性。

2、地震属性参数的提取
地震属性参数是从GeoFrame地震解释系统中提取出来的，它是下一步进行储层预测的关键。

根据井孔的目的层的厚度开时窗，读取时窗内的地震数据，进行地震特征参数的提取。

可供提取的地震参数有能量、频谱、波形等三大类32种参数。

沿层地震属性是以解释层位为基础，在地震数据体(剖面)中提取的属性，它的数值对应一个层位或一套地层，每个属性值对应一个x、y坐标。

提取方式有两类：沿一个解释层开一个常数时窗，在此时窗内提取地震属性，提取方式有4种(图4-3)。

用两个解释层提取某一段地层对应的地震属性，提取方式也有4种(图4-4)。

图4-3 单层位地震属性提取方式
图4-4层间地震属性提取方式
(1)、均方根振幅(RMS Amplitude)：均方根振幅是将振幅平方的平均值开平方，由于振幅值在平均前平方，因此，它对特别大的振幅非常敏感。

(2)、平均绝对值振幅(Average Absolute Amplitude)：平均绝对值振幅没有均方根振幅那样对特别大的振幅敏感。

(3)、最大波峰振幅(Maximum Peak Amplitude)：最大波峰振幅的求取方法是，对于每一道，在分析时窗里做一抛物线，恰好通过最大正的振幅值和它两边的两个采样点，沿着这曲线内插可得到最大波峰值振幅值。

(4)、平均波峰振幅(Average Peak Amplitude)：平均峰值振幅是对每一道在分析时窗里的所有正振幅值相加，得到总数除以时窗里的正振幅值采样数得到的。

(5)、最大波谷振幅(Maximum Trough Amplitude)：最大波谷振幅的求取方法是，对于每一道，在分析时窗里做一抛物线，恰好通过最大负的振幅值和它两边的两个采样点，沿着这曲线内插可得到最大波谷振幅值。

(6)、平均波谷振幅(Average Trough Amplitude)：平均波谷振幅是对每一道在分析时窗里的所有负振幅值相加，得到总数除以时窗里的负振幅值采样数得到的。

(7)、最大绝对值振幅(Maximum Absolute Amplitude)：计算每道的最大绝对值振幅的求取方法是，首先在分析时窗内计算出波峰和波谷的值，得出最大的波峰或波谷值，然后，画一抛物线，恰好通过最大波峰或波谷振幅值和它两边的两个采样点，沿着这曲线内插可得到最大绝对值振幅值。

(8)、总绝对值振幅(Total Absolute Amplitude)：总绝对值振幅是计算确定时窗内的所有道的绝对值振幅值。

(9)、总振幅(Total Amplitude)：每一道的总振幅是在层内对采样点求取总
的振幅值。

(10)、平均能量(Average Energy)：对于每一道的平均能量的求取方法是对分析时窗内的振幅值平方相加，对总数除以时窗内的采样数求得。

(11)、总能量(Total Energy)：对于每一道总能量的求取方法是对分析时窗内的振幅值平方相加求和得到的。

(12)、平均振幅(Mean Amplitude)：对于每一道的平均振幅的求取方法是对分析时窗内的振幅值相加，总数除以非零采样点数得到的。

(13)、振幅的平方差(Variance in Amplitude)：对于每一道的振幅的平方差的求取方法是，对分析时窗内的每个振幅值减去平均值累加，总数除以非零采样点数得到的。

(14)、振幅的立方差 (Skew in Amplitude)：对于每一道的振幅的立方差的求取方法是，对分析时窗内的所有采样点求取平均值，然后减去每道的平均值，计算差值的立方，求出这些值的总和，除以采样点数就可得到。

(15)、振幅的峰态 (Kurtosis in Amplitude)：对于每一道的振幅的峰态的求取方法是，对分析时窗内的所有采样点求取平均值，然后减去每道的平均值，计算差值的四次方，求出这些值的总和，除以采样点数就可得到。

(16)、有效带宽(Effective Bandwidth)：数据体时窗的有效带宽是由数据体的零延时的自相关除以采样周期与道两边所有自相关的总和之积而求得的。

(17)、弧长(Arc Length)：弧长是作为地震道的波形长度来定义的，它是在时窗内对所有地震道的变化范围的比例测量。

(18)、过零值平均频数(Average Zero Crossings Frequency)：过零值平均频率的计算方法是通过数据体时窗中的过零点的个数，和求出第一个通过零值的反射时间和最后一个通过零值的反射时间。

(19)、主频系列F1、F2、F3(Dominant Frequency Series F1、F2、F3)：对于所确定时窗的每一个输入道的估算值是由能量谱中的三个最主要频率分量组成。

F1是低频段中的峰值，F2是中间频段中的峰值，F3是高频段中的峰值。

(20)、峰值谱频率(Peak Spectral Frequency)：对于所确定时窗内的每一输入道，峰值谱频率的估算值是由能量谱中单一的最主要的频率组分组成。

峰值谱频率相似于主频系列，主频系列估算值是由能量谱中的三个最主要的频率段组成。

(21)、从谱的峰值到最高频率的斜率(Spectral Slope from Peak to Maximum Frequency)：这个属性表明了在分析时窗内高频成分被吸收的特点。

确定了一个感兴趣的最大值，就计算出在谱中的峰值频率到你设定的门槛值衰减比率。

如果斜率是一个高值，高频成分很快被吸收；如果斜率是一个低值，就没有信息被吸收。

(22)、大于门槛值百分比(Percent Greater than Threshold)：对于每一道来说，在分析时窗中，大于设定的门槛值的采样个数除以总采样个数，乘以100。

(23)、小于门槛值百分比(Percent Less than Threshold)：对于每一道来说，在分析时窗中，小于设定的门槛值的采样个数除以总采样个数，结果乘以100。

(24)、能量半衰时(Energy Half-Time)：在研究的时窗内，从上到下根据样点数求能量累加之和。

当能量之和达到计算时窗内总能量的一半时，到这点的样点个数除以总的样点个数为这点的能量半衰时。

(25)、能量半衰时斜率(Slope at Energy Half-Time)：能量半衰时斜率所计算的是当所累计的能量是总能量的一半时所需时间的能量曲线的斜率。

(26)、正采样点数与负采样点数的比率(Ratio of Positive to Negative Samples)：在分析时窗内对于每一道正采样数到负采样数的比率是由正采样数除以负采样数得到的。

(27)、波峰数(Number of Peaks)：波峰数计算的是分析时窗内的正波峰数，这个结果总是整数。

因为波峰在这里被认为是任意相对的最大值。

(28)、波谷数(Number of Troughs)：波谷数计算时窗内负波谷数，这个结果总是整数，因为波谷在这里被认为是任意相对的负最小值。

(29)、协变系数(Covariance Coefficient to Next CDP)：协变性系数是由两个相邻道之间所求的标准互相关计算得到的。

(30)、相关时移(Correlation Window Time Shift to Next CDP)：选择这个属性是用于计算时移的，是在一道和它的相邻道之间互相关。

(31)、平均信噪比(Average Signal- to-Noise Ratio)：平均信噪比是在层间计算时窗中多道的中心道平均信噪比。

(32)、相关长度(Correlation Length)：相关长度计算时窗里中心道和相邻道之间的相关系数减小的快慢的属性，长度的计算是中心道两边相关值变化，当相关值到达0.5时，中心道到这两点的平均距离。

这个平均距离是通过线性内插所估算的。

3、储层预测结果
建立在目的层时窗内提取的地震属性参数的基础上，利用地震资料，测井资料，地质资料，采用LPM储层预测技术，通过对地震属性参数的分析，根据井孔处的砂岩值与地震属性参数的线性相关程度，分别对扶余油层建立线性和非线性方程，预测了扶Ｉ油层、扶Ｉ油层上、中、下部油层的砂岩厚度（图4-5-图4-8）。