基于网格点的随机耦合四维地震反演油藏监测方法及装置

摘要

本发明提供基于网格点的随机耦合四维地震反演油藏监测方法及装置，该方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；根据多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；利用波阻抗变化量处于任意区间的概率体，监测油藏变化。本发明直接使用两次采集的叠后地震资料作为输入数据，减少了其它操作带来的误差；只反演波阻抗的变化量，将计算量限定在较小的范围内；给出波阻抗变化量在任意区间的概率体，定量描述了波阻抗的变化量及其不确定性，可以有效得反映生产过程中油藏的变化。

主权项

一种基于网格点的随机耦合四维地震反演油藏监测方法，其特征在于，所述基于网格点的随机耦合四维地震反演油藏监测方法包括：基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量；重复上述步骤，获得多个波阻抗变化量数据体；根据所述多个波阻抗变化量数据体，获得波阻抗变化量处于任意区间的概率体；利用所述波阻抗变化量处于任意区间的概率体，监测油藏变化；其中，所述基于网格点的随机耦合四维地震反演方法，通过耦合的两次随机三维地震反演获得波阻抗的变化量，包括：步骤1：构建时间域三维网格，空间采样率与地震资料相同，三维网格最上层各网格点的时间按照目标层顶界面的时间赋值，三维网格最下层各网格点的时间按照目标层底界面的时间赋值；步骤2：预设时间分辨率Δt，搜索目标层中顶底时间差的最大值t_m，三维网格中各道的网格点数设为n,步骤3：将测井获得的波阻抗曲线和岩性曲线转换到时间域，并将它们的值赋给距离最近的网格点；步骤4：通过序贯指示模拟将三维网格中井曲线以外各点的岩性值充填，作为这些点的原始岩性值；步骤5：通过序贯高斯模拟将三维网格中井曲线以外各点的波阻抗值充填，作为这些点的原始波阻抗值；步骤6：对于井曲线以外每一个网格点，利用指示克里金建立该点的岩性概率分布P_litho，利用序贯指示模拟获得该点的新岩性值，再利用简单克里金构建该点的波阻抗概率分布P_imp，利用序贯高斯模拟获得该点的新波阻抗值；设该点的原始岩性值为litho，原始波阻抗值是imp，新的岩性值是新的波阻抗值是分别在该点使用原始的波阻抗值和新的波阻抗值计算该点所在道的合成地震记录，记作syn和设该道时间1采集的实际地震记录是s₁；按照Metropolis‑Hasting方法判断是否用该点新的岩性值和波阻抗值代替原始岩性值和波阻抗值；具体做法是计算因子H₁，如果大于等于1，则用新的岩性和波阻抗值代替原始值；如果小于1，则用满足0到1之间概率均匀分布的随机数产生器生成一个随机数b，如果b小于则仍用新值代替原始值，否则保留原始值，抛弃新值；H₁因子计算公式如下：$\begin{array}{c}{H}_1=\left[-\frac{1}{2n_1^n}\underset{i=1}{\overset{m}{\Sigma }}{\left({\overline{syn}}^i-s_1^i\right)}^2+\ln \left(P_{litho}\left(litho\right)\right)+\ln \left(P_{imp}\left(imp\right)\right)\right]\\ -\left[-\frac{1}{2n_1^2}\underset{i=1}{\overset{m}{\Sigma }}{\left({\mathrm{syn}}^i-s_1^i\right)}^2+\ln \left(P_{litho}\left(\overline{litho}\right)\right)+\ln \left(P_{imp}\left(\overline{imp}\right)\right)\right]\end{array},$其中，m表示合成地震记录的采样点数，n₁表示时间1采集的地震记录的噪音程度，P_litho(litho)和表示原始岩性值litho和新岩性值在概率分布P_litho中的概率，P_imp(imp)和表示原始波阻抗值imp和新波阻抗值在概率分布P_imp中的概率；步骤7：多次重复步骤6，重复次数大于30次，最终获得三维岩性数据litho和波阻抗数据imp；步骤8：使用油藏数模结果，统计在两次地震采集时期之间油藏含油饱和度和孔隙度变化量的概率分布，用岩石物理方法将其转换为波阻抗变化量的概率分布；步骤9：使用步骤8获得的波阻抗变化量的概率分布，用序贯高斯模拟对每个网格点计算其波阻抗的变化量，作为这些点的原始波阻抗变化量值；步骤10：对三维网格的每一点，利用简单克里金计算各点波阻抗变化量的概率分布P_Δimp，用序贯高斯模拟获得该点新的波阻抗变化量；设原始波阻抗变化量为Δimp，新的值为用步骤7获得的该点的波阻抗值与原始波阻抗变化量的和imp+Δimp作为该点的波阻抗值，做该点所在道的合成地震记录Δsyn；再用步骤7获得的该点的波阻抗值与新的波阻抗变化量的和作为该点的波阻抗值，再做该点所在道的合成地震记录设时间2采集的地震记录为s₂；再次使用Metropolis‑Hasting方法判断是否保留新的波阻抗变化量；具体做法是计算因子H₂,如果大于等于1，则用新的波阻抗变化量代替原始波阻抗变化量；如果小于1，则用满足0到1之间均匀概率分布的随机数产生器生成随机数b；如果b小于则仍用新值代替原始值，否则保留原始值，放弃新值；H₂按如下公式计算：$\begin{array}{c}{H}_2=\left[-\frac{1}{2n_2^2}\underset{i=1}{\overset{m}{\Sigma }}{\left({\overline{\Delta syn}}^i-s_2^i\right)}^2+\ln \left(P_{\Delta imp}\left(\Delta imp\right)\right)\right]\\ -\left[-\frac{1}{2n_2^2}\underset{i=1}{\overset{m}{\Sigma }}{\left({\mathrm{\Delta syn}}^i-s_2^i\right)}^2+\ln \left(P_{\Delta imp}\left(\overline{\Delta imp}\right)\right)\right]\end{array},$其中，n₂表示时间2采集的地震数据的噪音程度，P_Δimp(Δimp)和分别表示原始波阻抗变化量和新波阻抗变化量在概率分布P_Δimp中的概率；步骤11：多次重复步骤10，重复次数大于30次，获得一个三维波阻抗变化量数据体，从而通过两次随机三维地震反演获得波阻抗的变化量。