技术领域
[0001] 本
发明涉及
地震勘探技术领域，尤其涉及一种基于希尔伯特谱和BP神经网络的储层孔隙度预测方法。
背景技术
[0002] 储层孔隙度是衡量油气储层
岩石中所含孔隙体积多少的一种参数，孔隙度反映岩石储存
流体的能
力，是油气勘探中的一项重要指标。由于由地震记录获取的弹性属性与不同储层孔隙度间复杂的非线性映射关系，传统的模型驱动的孔隙度预测方法往往难以取得较好的效果。
[0003] 同时，地震记录的希尔伯特谱特征与包括孔隙度在内的储层特性有密切联系，但目前在包括
人工智能方法在内的地震
数据处理方法中没有得到充分应用。
发明内容
[0004] 本发明针对
现有技术的以上不足，提出一种基于希尔伯特谱和BP神经网络的
地震勘探记录处理方法，该方法可用于预测储层孔隙度。同
时针对勘探区域钻井/
测井普遍较少的情况，提出一种训练、测试神经网络的方法。
[0005] 本发明的目的可以通过以下技术方案来实现：
[0006] S1：读取勘探区域的叠后或自激自收单道地震记录并应用以下公式进行预处理：
[0007]
[0008] 其中os(t)是叠后或自激自收单道地震记录，N是
叠加次数。
[0009] S2：对经S1处理后的单道记录进行经验模态分解(EMD)，得到本征模态函数(IMF)。具体包括以下步骤：
[0010] S2-1：找出s(t)的所有极大值，用三阶样条函数进行插值，得到其上包络序列max(t)。
[0011] 用同样的方法得到其下包络序列min(t)。对上下包络取平均得到平均值序列：
[0012]
[0013] S2-2：从s(t)减去m(t)：
[0014] h(t)＝s(t)-m(t)
[0015] 考察h(t)是否满足IMF的两个条件：(1)极值个数和零点个数相差不超过1；(2)在序列的任一点，上下包络均值为0.
[0016] 若h(t)不能同时满足IMF的两个条件，则将h(t)作为原序列，重复S2-1、S2-2，直至h(t)同时满足IMF的两个条件。
[0017] 当h(t)同时满足IMF的两个条件，则h(t)是第一个IMF，记作c1(t)。
[0018] S2-3：从原序列减去所得到的IMF，得到余值序列：
[0019] r1(t)＝s(t)-c1(t)
[0020] 以r1(t)作为原序列，重复S2-1到S2-3，得到c2(t)和r2(t)，…，ci(t)和ri(t)，…，直到rn(t)成为单道序列或常数列。此时EMD完成，原序列s(t)被分解为n个IMF和一个余值序列。
[0021] S3：对各IMF作卷积运算，得到其希尔伯特变换：
[0022]
[0024]
[0025] S4：求取s(t)的希尔伯特瞬时
能量、希尔伯特边际谱：
[0026] 按下式求取s(t)的希尔伯特谱：
[0027]
[0028] 计算希尔伯特谱幅值平方对
频率的积分，得到希尔伯特瞬时能量：
[0029]
[0030] 计算希尔伯特谱幅值对时间的积分，得到希尔伯特边际谱：
[0031]
[0032] S5：获取神经网络的输入：
[0033] 假设地震记录
采样数为N，则S4所得希尔伯特瞬时能量和希尔伯特边际谱均为包含N个值的序列，可分别表示为：
[0034] A1＝[A1(1)，A1(2)，…，A1(N)]
[0035] A2＝[A2(1)，A2(2)，…，A2(N)]
[0036] 将2N维向量
[0037] X＝[x1，x2，…，x2N]＝[A1，A2]＝[A1(1)，A1(2)，…，A1(N)，A2(1)，A2(2)，…，A2(N)][0038] 作为神经网络的输入。
[0039] S6：选取训练-测试样本并获取训练-测试样本的期望输出：
[0040] 将钻井/测井
位置的
地震道记录作为训练-测试样本，根据由地震数据速度模型进行时深转换，得到N个采样时间对应的N个深度，由钻井/测井数据获取各个深度处的孔隙度组成N维向量
[0041] Y＝[y1，y2，…，yN]
[0042] 作为训练-测试样本的期望输出。
[0043] S7：建立并初始化BP神经网络。所用BP神经网络包括
输入层、隐层、
输出层。输入层有2N个神经元，输出层有N个神经元，隐层的层数、各隐层神经元个数以及激励函数可根据实际情况确定和调整。相邻层神经元间全连接，同层神经元间无连接。对神经网络各层间的权值矩阵和各神经元的
阈值随机初始化。
[0044] S8：训练、测试BP神经网络。假设共有M个样本用于训练和测试，将样本从1到M编号，形成含M个样本的有序数据集。对此数据集执行训练-测试流程。
[0045] 所述训练-测试流程，是针对地震勘探中钻井/测井数量教少，因而可用于训练和测试的地震道较少的情况，所提出的一种充分利用由钻井/测井获取的储层参数，将训练和测试统一于同一过程的方法。具体执行过程通过如下递归方法进行：
[0046] 假设使用含有N个从1到N编号的样本的有序数据集执行训练-测试流程，i号样本的输入和期望输出分别用Xi和Yi表示，过程如下：
[0047] (1)若N＝1：
[0048] 使用N号样本对S7建立的神经网络使用BP
算法进行
迭代训练，直到输出与YN的误差在预设的范围之内；
[0049] (2)若N＞1：
[0050] 先使用由1号到N-1号样本组成的有序数据集对S7建立的神经网络执行训练-测试流程，完成后再使用N号样本进行测试，若输出与YN的误差在预设的范围之内，流程结束；
[0051] 否则将原来的i号(1≤i≤N-1)样本序号变为i+1，N号样本序号变为1，重新执行上段所述操作。
[0052] S9：使用BP神经网络对其余检波点孔隙度随深度的分布进行预测。
[0053] 将非钻井/测井位置的由S5获得的2N维向量作为经过S8训练、测试的BP神经网络的输入，得到的N维向量表示N个采样时间对应的孔隙度。对该道地震记录进行时深转换得到N个采样时间对应的深度，即可得到该检波点处孔隙度随深度的分布。
附图说明
[0054] 图1为本发明流程示意图；
[0055] 图2为本发明
实施例所用BP神经网络示意图；
[0056] 图3为本发明提出的训练-测试流程示意图。
具体实施方式
[0057] 下面结合具体实施例对本发明进行详细说明。该实施例描述较为具体和详细，但并不能理解为对本发明
专利保护范围的限制。在本发明所述的技术范围内，可轻易想到的各种等效的
修改都应涵盖在本发明的保护范围之内。
[0058] 该实施例包含以下步骤：
[0059] S1：选取在部分检波点位置有钻井/测井数据的勘探工区，由多次
覆盖技术获得的地震数据作为本方法的应用对象。
[0060] S2：读取勘探工区的叠后单道地震记录，应用以下公式进行预处理：
[0061]
[0062] 其中os(t)是叠后单道地震记录，N是叠加次数。
[0063] S3：对经S2处理后的单道记录进行经验模态分解(EMD)，得到本征模态函数(IMF)。具体包括以下步骤：
[0064] S3-1：找出s(t)的所有极大值，用三阶样条函数进行插值，得到其上包络序列max(t)。
[0065] 用同样的方法得到其下包络序列min(t)。对上下包络取平均得到平均值序列：
[0066]
[0067] S3-2：从s(t)减去m(t)：
[0068] h(t)＝s(t)-m(t)
[0069] 考察h(t)是否满足IMF的两个条件：
[0070] (1)极值个数和零点个数相差不超过1；
[0071] (2)在序列的任一点，上下包络均值为0。
[0072] 若h(t)不能同时满足IMF的两个条件，则将h(t)作为原序列，重复S3-1、S3-2，直至h(t)同时满足IMF的两个条件。
[0073] 当h(t)同时满足IMF的两个条件，则h(t)是第一个IMF，记作c1(t)。
[0074] 在实际应用中，假设m(t)为h(t)上下包络的均值序列，计算
[0075]
[0076] 当e的值小于0.3时，可认为h(t)满足IMF的第二个条件。
[0077] S3-3：从原序列减去所得到的IMF，得到余值序列：
[0078] r1(t)＝s(t)-c1(t)
[0079] 以r1(t)作为原序列，重复S3-1到S3-3，得到c2(t)和r2(t)，…，ci(t)和ri(t)，…，直到rn(t)成为单道序列或常数列。此时EMD完成，原序列s(t)被分解为n个IMF和一个余值序列。
[0080] S4：对各IMF作卷积运算，得到其希尔伯特变换：
[0081]
[0082] 则第k个IMF对应的解析信号为：
[0083]
[0084] S5：求取s(t)的希尔伯特瞬时振幅谱、希尔伯特边际谱：
[0085] 按下式计算s(t)的希尔伯特谱：
[0086]
[0087] 计算希尔伯特谱平方对频率的积分，得到希尔伯特瞬时能量：
[0088]
[0089] 计算希尔伯特谱对时间的积分，得到希尔伯特边际谱：
[0090]
[0091] S6：获取神经网络的输入：
[0092] 假设地震记录采样数为N，则S5所得希尔伯特瞬时能量和希尔伯特边际谱均为包含N个值的序列，可分别表示为：
[0093] A1＝[A1(1)，A1(2)，…，A1(N)]
[0094] A2＝[A2(1)，A2(2)，…，A2(N)]
[0095] 将2N维向量
[0096] X＝[x1，x2，…，x2N]＝[A1，A2]＝[A1(1)，A1(2)，…，A1(N)，A2(1)，A2(2)，…，A2(N)][0097] 作为神经网络的输入。
[0098] S7：选取训练-测试样本并获取训练-测试样本的期望输出：
[0099] 将钻井/测井位置的地震道记录作为训练-测试样本，根据地震数据速度模型进行时深转换，得到N个采样时间对应的N个深度，由钻井/测井数据获取各个深度处的孔隙度组成N维向量
[0100] Y＝[y1，y2，…，yN]
[0101] 作为训练-测试样本的期望输出。
[0102] S8：建立并初始化BP神经网络。所用BP神经网络包括输入层、两个隐层、输出层。输入层和两个隐层各有2N个神经元，输出层有N个神经元。相邻层神经元间全连接，同层神经元间无连接。隐层和输出层均使用线性激励函数。对神经网络各层间的连接权值矩阵和各神经元的阈值随机初始化。
[0103] S9：训练、测试BP神经网络。假设S7共获取了M个样本用于训练和测试，构造由这M个样本的M个神经网络输入向量组成的数组X，和由这M个样本的M个期望输出向量组成的数组Y。这两个数组的第i个元素X(i)和Y(i)分别是2N维向量和N维向量。以M、X、Y为参数执行递归函数LT对S8建立的神经网络执行训练-测试流程。
[0104] 对LT函数作如下描述：
[0105] LT函数的参数是正整数W，含有W个元素的数组A、B，A和B的元素都是向量。
[0106] LT函数的执行过程如下：
[0107] 若W＝1：
[0108] 以A(W)为输入、B(W)为标记对神经网络使用BP算法进行迭代训练，更新各层间的连接权值矩阵和各神经元的阈值直到输出与B(W)的误差在预设的范围之内；
[0109] 若W＞1：
[0110] 先以W-1和A的前W-1个元素构成的数组、B的前W-1个元素构成的数组为参数执行LT，再以A(W)为输入、B(W)为期望输出对神经网络进行测试。若输出与B(W)的误差在预设的范围之内，函数执行结束；
[0111] 否则用原来的A(i)、B(i)(1≤i≤W-1)替换A(i+1)、B(i+1)，用A(W)、B(W)替换原来的A(1)、B(1)，重新执行上段所述操作。
[0112] 误差计算公式：
[0113]
[0114] Y为由钻井/测井技术得到的标记向量/期望输出向量，Z为实际输出向量。
[0115] 使用某个样本对神经网络进行训练时，当e＜δ，训练结束。使用某个样本对神经网络进行测试时，若e＜δ，测试合格。δ的值参考其他孔隙度预测技术的精确度确定。
[0116] S10：使用BP神经网络对其余检波点孔隙度随深度的分布进行预测。
[0117] 将非钻井/测井位置的由S6获得的2N维向量作为经过S9训练、测试的BP神经网络的输入，得到的N维向量表示N个采样时间对应的孔隙度。对该道地震记录进行时深转换，即可得到该检波点处孔隙度随深度的分布。
