本发明提供基于网格点的随机耦合四维地震反演油藏监测方法及装置，该方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；根据多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；利用波阻抗变化量处于任意区间的概率体，监测油藏变化。本发明直接使用两次采集的叠后地震资料作为输入数据，减少了其它操作带来的误差；只反演波阻抗的变化量，将计算量限定在较小的范围内；给出波阻抗变化量在任意区间的概率体，定量描述了波阻抗的变化量及其不确定性，可以有效得反映生产过程中油藏的变化。
1.一种基于网格点的随机耦合四维地震反演油藏监测方法，其特征在于，所述基于网格点的随机耦合四维地震反演油藏监测方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；
根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；
利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化；
其中，所述基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；
步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm，三维网格中各道的网格点数设为n,步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；
步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点的原始岩性值；
步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；
步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是 新的波阻抗值是 分别在该点使用原始的波阻抗值和新的波阻抗值计算该点所在道的合成地震记录，记作syn和 设该道时间1采集的实际地震记录是s1；
按照Metropolis-Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始岩性值和波阻抗值；具体做法是计算因子H1，如果 大于等于1，则用新的岩性和波阻抗值代替原始值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值；H1因子计算公式如下：其中，m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率，Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率；
步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；
步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；
步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；
步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2；再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H2,如果 大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于 则仍用新值代替原始值，否则保留原始值，放弃新值；H2按如下公式计算：其中，n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始波阻抗变化量和新波阻抗变化量在概率分布PΔimp中的概率；
步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。
2.如权利要求1所述基于网格点的随机耦合四维地震反演油藏监测方法，其特征在于，重复上述步骤，获得多个波阻抗变化量数据体，包括：整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于
3.一种基于网格点的随机耦合四维地震反演油藏监测装置，其特征在于，所述基于网格点的随机耦合四维地震反演油藏监测装置包括：地震反演单元，用于基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；
获取单元，用于根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；
监测单元，用于利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化；
其中，所述地震反演单元是基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；
步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm，三维网格中各道的网格点数设为n,步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；
步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点的原始岩性值；
步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；
步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新岩性值是 新波阻抗值是 分别在该点使用原始波阻抗值和新波阻抗值计算该点所在道的合成地震记录，记作syn和 设该道时间1采集的实际地震记录是s1；
按照Metropolis-Hasting方法判断是否用该点新岩性值和波阻抗值代替原始岩性值和波阻抗值；具体做法是计算因子H1，如果 大于等于1，则用新的岩性和波阻抗值代替原始的值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值；H1因子计算公式如下：其中，m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率，Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率；
步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；
步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；
步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；
步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2；再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H2,如果 大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于 则仍用新值代替原始值，否则保留原始值，放弃新值；H2按如下公式计算：其中，n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始的波阻抗变化量和新的波阻抗变化量在概率分布PΔimp中的概率；
步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。
4.如权利要求3所述基于网格点的随机耦合四维地震反演油藏监测装置，其特征在于：所述地震反演单元，进一步用于整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于30的自然数。
基于网格点的随机耦合四维地震反演油藏监测方法及装置
技术领域
[0001] 本发明涉及油藏地球物理技术，具体地说涉及基于网格点的随机耦合四维地震反演油藏监测方法及装置。
背景技术
[0002] 随机地震反演适用于井资料丰富的地区，可以给出高分辨率的反演结果以及反演不确定性的评估，是一种重要的油藏地球物理技术。随机地震反演可以分为基于地震道的反演方法和基于网格点的反演方法。四维地震反演技术是利用两次或多次采集的地震资料，在互均化处理的基础上进行地震反演获得油藏变化信息的技术，是老油田增加新储量，进行开发方案调整，提高采收率的重要方法。四维地震反演可以分为非耦合反演方法、耦合反演方法和差异反演方法。非耦合反演就是对多次地震资料单独进行反演得到不同时期的油藏模型，然后比较它们的差异得到油藏的变化量；耦合反演是先对第一次地震资料进行反演得到油藏模型，然后让这个模型参与到其它时期地震资料的反演中，最后比较多个反演结果获得油藏的变化量或者直接获得油藏的变化量；差异反演是直接对两次地震资料的差异进行反演，获得油藏的变化量。随机四维地震反演结合了随机反演和四维反演的优势，不少学者对其开展了研究。
[0003] Leon Barens和Aline Franche(SEG，2005年会论文集)进行了基于道的随机非耦合四维地震反演数值实验。他们分别对两次地震资料完成了基于道的随机反演，获得对应两个时期的两套反演结果，通过对比这两套反演结果纵横波阻抗的交汇图，指出在反演使用子波正确的情况下该方法可以发现生产引起的油藏变化。
[0004] Arild Buland等人在2006年(GEOPHYSICS,71卷3期)提出基于地震道的叠前随机差异四维反演方法。该方法认为不同时期油藏弹性参数比值的对数值满足高斯分布，并利用不同时期地震资料的差异与弹性参数差异的线性关系构建似然概率分布，进而获得不同时期油藏弹性参数比值的对数值满足的后验概率分布，最终解析地给出上述后验概率分布的期望值和标准差表征反演结果和反演结果的不确定性。
[0005] Youli Quan等人(SEG2008年年会论文集)提出另一种基于地震道的随机耦合四维地震反演方法，该方法将第一次反演的结果作为第二次反演的初始模型，利用集合卡尔曼滤波方法获得不同时刻的油藏速度模型。
[0006] 后两种基于道的反演均没有考虑反演结果中各道之间的相关性。与之相比，基于网格点的方法则考虑三维空间中各网格点与周围网格点取值的相关性。
[0007] Helene Hafslund Veire等人在2006年(GEOPHYSICS,71卷5期)提出一种基于网格点的随机四维地震反演方法。该方法属于差异反演，利用不同时期地震资料的AVO截距和梯度的差异构建似然概率分布，使用马尔可夫链蒙特卡洛方法（MCMC）对后验概率分布采样获得油藏压力、饱和度等参数的变化量。由于该方法使用地震资料的AVO截距和梯度作为反演的输入而不是直接使用地震资料，AVO截距和梯度的提取误差可能会影响后续反演的精度。
[0008] T.Hong(SEG2007年年会论文集)等提出一种基于网格点的随机耦合四维反演方法。该方法构建了一个包含油藏的渗透率、含油饱和度、压力以及岩石物理模型参数等的后验概率分布，使用马尔可夫链蒙特卡罗方法采样这个后验概率分布获得反演结果。由于后验概率分布涉及的参数较多，该方法在实际应用时可能面临巨大的计算量。
[0009] Long Jin(SEG2007年年会论文集)等提出另一种基于网格点的随机四维地震反演方法。该方法结合油藏模拟和岩石物理方法，使用非常快速模拟退火算法获得满足四维地震资料和井生产曲线的孔隙度模型。与马尔可夫链蒙特卡洛算法相比，非常快速模拟退火算法不能很好得反映反演结果的不确定性。
[0010] 在实现本发明过程中，发明人发现现有技术中至少存在如下问题：现有技术亟待需要一种较准确监测油藏变化的技术方案。
发明内容
[0011] 本发明实施例提供一种基于网格点的随机耦合四维地震反演油藏监测方法及装置，以提供一种较准确监测油藏变化的技术方案。
[0012] 一方面，本发明实施例提供了一种基于网格点的随机耦合四维地震反演油藏监测方法，所述基于网格点的随机耦合四维地震反演油藏监测方法包括：
[0013] 基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；
[0014] 重复上述步骤，获得多个波阻抗变化量数据体；
[0015] 根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；
[0016] 利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化。
[0017] 优选的，在本发明一实施例中，所述基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：
[0018] 步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；
[0019] 步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm，三维网格中各道的网格点数设为n,
[0020] 步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；
[0021] 步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点原始的岩性值；
[0022] 步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；
[0023] 步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是 ，新的波阻抗值是 ；分别在该点使用原始波阻抗值和新波阻抗值计算该点所在道的合成地震记录，记作syn和 ；设该道时间1采集的实际地震记录是s1；
[0024] 按照Metropolis-Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始的岩性值和波阻抗值；具体做法是计算因子H1，如果 大于等于1，则用新的岩性和波阻抗值代替原始的值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值；H1因子计算公式如下：
[0025]
[0026] 其中，m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率，Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率；
[0027] 步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；
[0028] 步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；
[0029] 步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；
[0030] 步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2；再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H2,如果 大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于 则仍用新值代替原始值，否则保留原始值，放弃新值；H2按如下公式计算：
[0031]
[0032] 其中，n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始波阻抗变化量和新波阻抗变化量在概率分布PΔimp中的概率；
[0033] 步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。
[0034] 优选的，在本发明一实施例中，所述重复上述步骤，获得多个波阻抗变化量数据体，包括：整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于30的自然数。
[0035] 另一方面，本发明实施例提供了一种基于网格点的随机耦合四维地震反演油藏监测装置，所述基于网格点的随机耦合四维地震反演油藏监测装置包括：
[0036] 地震反演单元，用于基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；
[0037] 获取单元，用于根据所述多个波阻抗变化量数据体，获得处波阻抗变化量处于任意区间的概率体；
[0038] 监测单元，用于利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化。
[0039] 优选的，在本发明一实施例中，所述地震反演单元基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：
[0040] 步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；
[0041] 步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm，三维网格中各道的网格点数设为n,
[0042] 步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；
[0043] 步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点的原始岩性值；
[0044] 步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；
[0045] 步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是 新的波阻抗值是 分别在该点使用原始的波阻抗值和新的波阻抗值计算该点所在道的合成地震记录，记作syn和 设该道时间1采集的实际地震记录是s1；
[0046] 按照Metropolis-Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始的岩性值和波阻抗值；具体做法是计算因子H1，如果 大于等于1，则用新的岩性和波阻抗值代替原始的值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值；H1因子计算公式如下：
[0047]
[0048] 其中，m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率，Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率；
[0049] 步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；
[0050] 步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；
[0051] 步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；
[0052] 步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2；再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H2,如果 大于等于1，则用新的波阻抗变化量代替原始的波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于 ，则仍用新值代替原始值，否则保留原始值，放弃新值；H2按如下公式计算：
[0053]
[0054] 其中，n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始波阻抗变化量和新波阻抗变化量在概率分布PΔimp中的概率；
[0055] 步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。
[0056] 优选的，在本发明一实施例中，所述地震反演单元，进一步用于整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于30的自然数。
[0057] 上述技术方案具有如下有益效果：因为采用所述基于网格点的随机耦合四维地震反演油藏监测方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化的技术手段，所以达到了如下的技术效果：与上述现有方法相比，本发明是一种基于网格点的反演方法，考虑了油藏模型中各网格点与周围点取值的相关性。其次，本发明直接使用两次采集的叠后地震资料作为输入数据，减少了其它操作带来的误差。再次，本发明只反演波阻抗的变化量，将计算量限定在较小的范围内。最后，本发明使用贝叶斯框架和马尔可夫链蒙特卡洛采样方法，给出波阻抗变化量在任意区间的概率体，定量描述了波阻抗的变化量及其不确定性，可以有效得反映生产过程中油藏的变化。
附图说明
[0058] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0059] 图1为本发明实施例一种基于网格点的随机耦合四维地震反演油藏监测方法流程图；
[0060] 图2为本发明应用实例时间1的实际地震记录与合成地震记录示意图；
[0061] 图3为本发明应用实例时间2的实际地震记录与合成地震记录示意图；
[0062] 图4为本发明应用实例波阻抗变化量大于500g/cc*m/s的概率体的平面图；
[0063] 图5为本发明应用实例波阻抗变化量大于500g/cc*m/s的概率体的剖面图；
[0064] 图6为本发明实施例一种基于网格点的随机耦合四维地震反演油藏监测装置结构示意图。
具体实施方式
[0065] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0066] 如图1所示，为本发明实施例一种基于网格点的随机耦合四维地震反演油藏监测方法流程图，所述基于网格点的随机耦合四维地震反演油藏监测方法包括：
[0067] 101、基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；
[0068] 102、重复上述步骤，获得多个波阻抗变化量数据体；
[0069] 103、根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；
[0070] 104、利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化。
[0071] 优选的，所述基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：
[0072] 步骤1：构建时间域三维网格，空间采样率与地震数据相同。三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值。
[0073] 步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm。三维网格中各道的网格点数为n, 这样做保证三维网格各道的时间分辨率均大于预设时间分辨率。目标层的顶底时间差在三维网格范围内不应变化太大，否则，应适当调整目标层的顶底时间。
[0074] 步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点。
[0075] 步骤4：通过对井资料的统计，获得目标层岩性的分布直方图和变差函数，然后使用序贯指示模拟将三维网格中井曲线以外各点的岩性值填充，作为这些点的原始岩性值。
[0076] 步骤5：对井资料分岩性统计波阻抗的直方图和变差函数。通过序贯高斯模拟填充三维网格中井曲线以外各点的波阻抗值，作为这些点的原始波阻抗值。在模拟每一点的波阻抗值时，要参考该点的岩性值，使用该点附近相同岩性点的波阻抗值、该岩性的波阻抗直方图和变差函数完成模拟。
[0077] 步骤6：对于井曲线以外的每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值。设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是 新的波阻抗值是 再分别在该点使用原始波阻抗值和新波阻抗值计算该点所在道的合成地震记录，记作syn和 设该道时间1采集的实际地震记录是s1。按照Metropolis-Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始的岩性值和波阻抗值。具体做法是计算因子H1，如果 大于等于1，则用新的岩性和波阻抗值代替原始值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值。H1因子计算公式如下：
[0078]
[0079] 其中m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率。Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率。
[0080] 步骤7：重复45次步骤6，最终获得三维岩性数据litho和波阻抗数据imp。步骤4到7实现了基于网格点的三维随机地震反演方法，这是用马尔可夫链蒙特卡罗算法对波阻抗和岩性的后验概率分布采样的过程。如图2所示，为本发明应用实例时间1的实际地震记录与合成地震记录示意图，其中，中间部分为目标层，深黑色为地震记录，浅黑灰色为用波阻抗数据imp正演获得的合成地震记录，二者在目标层差异很小，验证了该三维反演的可靠性。
[0081] 步骤8：使用油藏数模结果，统计在两次地震采集时刻之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布，并将这个概率分布拟合为一个高斯分布。
[0082] 步骤9：对步骤1和2构建的三维网格，在已知步骤8获得的波阻抗变化量的概率分布后，使用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值。
[0083] 步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量。设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn。再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2。再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量。具体做法是计算因子H2，如果 大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b。如果b小于 则仍用新值代替原始值，否则保留原始值，放弃新值。H2因子计算公式如下：
[0084]
[0085] 其中n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始波阻抗变化量和新波阻抗变化量在概率分布PΔimp中的概率；
[0086] 步骤11：重复45次步骤10，获得一个三维波阻抗变化量数据体Δimp，从而通过两次随机三维地震反演获得波阻抗的变化量。上述步骤9到11完成对波阻抗变化量的后验概率分布采样的过程。用步骤7获得的波阻抗模型与步骤11获得的波阻抗变化量模型相加可获得时间2对应的波阻抗模型，即imp+Δimp。如图3所示，为本发明应用实例时间2的实际地震记录与合成地震记录示意图，其中，中间部分为目标层，深黑色为地震记录，浅黑灰色为用波阻抗模型imp+Δimp正演获得的合成地震记录，二者在目标层差异很小，进一步验证了该四维反演的可靠性。
[0087] 优选的，所述重复上述步骤，获得多个波阻抗变化量数据体，包括：整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于30的自然数。例如，重复32次上述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，获得32个波阻抗变化量体，然后利用这32个波阻抗变化量体统计各点波阻抗变化量大于500g/cc*m/s（此处可以取波阻抗值的任意区间，该实施例中以大于500g/cc*m/s为例）的概率。所得概率体的平面图和剖面图分别如图4和图5所示；图4为本发明应用实例波阻抗变化量大于500g/cc*m/s的概率体的平面图；图5为本发明应用实例波阻抗变化量大于500g/cc*m/s的概率体的剖面图（中间部分为目标层）。
[0088] 对应于上述方法实施例，如图6所示，为本发明实施例一种基于网格点的随机耦合四维地震反演油藏监测装置结构示意图，所述基于网格点的随机耦合四维地震反演油藏监测装置包括：
[0089] 地震反演单元61，用于基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；
[0090] 获取单元62，用于根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；
[0091] 监测单元63，用于利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化。
[0092] 优选的，所述地震反演单元61基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：
[0093] 步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；
[0094] 步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值tm，三维网格中各道的网格点数设为n,
[0095] 步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；
[0096] 步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点的原始岩性值；
[0097] 步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；
[0098] 步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布Plitho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布Pimp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是 新的波阻抗值是 分别在该点使用原始的波阻抗值和新的波阻抗值计算该点所在道的合成地震记录，记作syn和 设该道时间1采集的实际地震记录是s1；
[0099] 按照Metropolis-Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始的岩性值和波阻抗值；具体做法是计算因子H1：如果 大于等于1，则用新的岩性和波阻抗值代替原始值；如果 小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于 则仍用新值代替原始值，否则保留原始值，抛弃新值；H1因子计算公式如下：
[0100]
[0101] 其中，m表示合成地震记录的采样点数，n1表示时间1采集的地震记录的噪音程度，Plitho(litho)和 表示原始岩性值litho和新岩性值 在概率分布Plitho中的概率，Pimp(imp)和 表示原始波阻抗值imp和新波阻抗值 在概率分布Pimp中的概率；
[0102] 步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；
[0103] 步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；
[0104] 步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；
[0105] 步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布PΔimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为 用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和imp+ 作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s2；再次使用Metropolis-Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H2,如果 大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果 小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于 则仍用新值代替原始值，否则保留原始值，放弃新值；H2按如下公式计算：
[0106]
[0107] 其中，n2表示时间2采集的地震数据的噪音程度，PΔimp(Δimp)和 分别表示原始波阻抗变化量和新波阻抗变化量在概率分布PΔimp中的概率；
[0108] 步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。
[0109] 优选的，所述地震反演单元61，进一步用于整体重复N次所述步骤4、步骤5、步骤6、步骤7、步骤9、步骤10、步骤11，相当于对波阻抗变化量的后验概率分布采样N次，获得N个波阻抗变化量数据体，其中N为大于或等于30的自然数。
[0110] 本发明实施例上述技术方案具有如下有益效果：因为采用所述基于网格点的随机耦合四维地震反演油藏监测方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化的技术手段，所以达到了如下的技术效果：与上述现有方法相比，本发明是一种基于网格点的反演方法，考虑了油藏模型中各网格点与周围点取值的相关性。其次，本发明直接使用两次采集的叠后地震资料作为输入数据，减少了其它操作带来的误差。再次，本发明只反演波阻抗的变化量，将计算量限定在较小的范围内。最后，本发明使用贝叶斯框架和马尔可夫链蒙特卡洛采样方法，给出波阻抗变化量在任意区间的概率体，定量描述了波阻抗的变化量及其不确定性，可以有效得反映生产过程中油藏的变化。
[0111] 本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块（illustrative logical block），单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性（interchangeability），上述的各种说明性部件（illustrative components），单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
[0112] 本发明实施例中所描述的各种说明性的逻辑块，或单元都可以通过通用处理器，数字信号处理器，专用集成电路（ASIC），现场可编程门阵列或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。
[0113] 本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。
[0114] 在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线（DSL）或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片（disk）和磁盘（disc）包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。
[0115] 以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
