地球物理学报  2011, Vol. 54 Issue (7): 1900-1911   PDF    
基于离散余弦变换的集合Kalman滤波方法对河流相油藏渗透率场的反演
王玉斗1 , 何健1 , 李高明2     
1. 中国石油大学(华东)物理科学与技术学院,东营 257061;
2. The University of Tulsa, Tulsa Oklahoma, USA 74104
摘要: 针对河流相油藏参数场的非高斯性及油藏生产数据与地质模型之间的非线性关系,提出了耦合离散余弦变换的半迭代集合Kalman滤波方法对河流相油藏进行生产历史拟合及反演.将河流相油藏非高斯渗透率场利用离散余弦变换变为类似高斯分布的DCT(Discrete cosine transform)系数,利用能量压缩特性对DCT系数进行截断,使得反演渗透率场的分布变得光滑,较好地解决了河流相油藏的非高斯分布问题.半迭代集合Kalman滤波法很好地解决了油藏生产数据与地质模型之间的非线性问题.理论证明了半迭代集合Kalman滤波与EnRML方法的等效性.比较了标准EnKF、耦合DCT的EnKF及耦合DCT的半迭代EnKF法对河流相油藏渗透率场的反演.讨论了DCT系数保留模式及系数保留数量对耦合DCT的半迭代集合Kalman滤波方法反演结果的影响,结果表明三角形系数保留模式可以用更少的系数保留参数场的地质特征,当保留2.5%的DCT系数时,反演结果已经充分反映了河流相油藏的地质特征.
关键词: 离散余弦变换      集合Kalman滤波      河流相      历史拟合      反演     
Inverse of the fluvial channel reservoir permeability field using ensemble Kalman filter based on discrete cosine transform
WANG Yu-Dou1, HE Jian1, LI Gao-Ming2     
1. College of Physics Science & Technology, China University of Petroleum, Dongying 257061, China;
2. The University of Tulsa, Tulsa Oklahoma 74104, USA
Abstract: Focusing on the non-Gaussian distribution of the field and the nonlinearity between the production data and the reservoir model of the fluvial channel reservoir, the HIEnKF method coupling discrete cosine transform was used to invert for the fluvial channel reservoir model by automatic history matching method. The non-Gaussian permeability field is transformed to nearly Gaussian distributed DCT coefficient by discrete cosine transform. DCT exhibits excellent energy compaction for highly correlated images. The estimated permeability field is smoother if the low frequency DCT coefficients are retained. The non-Gaussian problem of the fluvial channel reservoir can be solved by the discrete cosine transform. The HIEnKF method can solve the nonlinearity of the production data and the reservoir model. The equivalence of the HIEnKF to EnRML is proved theoretically. The results of standard EnKF, standard EnKF coupling with DCT and HIEnKF coupling with DCT were compared. The effects of retained mode and the retained number of the DCT coefficients were discussed. The results show that less DCT coefficients need to be retained to capture the geological structure when the triangular area mode was used to truncate the coefficients. The estimated results can get the main geological structure of the fluvial channel reservoir when only 2.5% DCT coefficients were retained..
Key words: Discrete cosine transform      Ensemble Kalman filter      Fluvial channel      History matching      Inverse     
1 引言

油藏地质模型的测量手段包括地震勘探、测井、试井和岩心分析等手段.但是由于经济或技术的原因,这些测量手段获得的地质信息有很大的分散性和随机性,使得地质模型具有多解性和不确定性.利用油藏生产历史拟合方法,根据大量的、带有相关地质及流体信息的生产数据(井底流压、产油量、含水率、生产气油比等)来反演油藏地质模型,是油藏描述的重要手段.但是传统的自动历史拟合方法是基于梯度的方法,需要利用伴随矩阵方法(Adjoint method)求解敏感系数矩阵(Sensitivity matrix),即目标函数对模型参数的导数,但求解步骤繁琐且费时较长,算法稳定性差,对硬件要求高,其应用一直受到限制[12].集合Kalman 滤波(EnKF)方法是新近兴起的一种数据拟合技术,1994 年第一次由Evensen[3]在海洋动力学中提出.2003 年该理论被Naevdal[4]第一次引入石油领域后,在生产数据、地震数据拟合及油藏地质参数描述等方面已引起越来越多研究者的注意,现已被成功地应用于一些实际油藏[5~7].国内也开始将该方法应用于地表参数反演和岩石力学参数估计等方面[8~11],但是在石油地质和油藏工程领域应用很少[12].

由于河流相油藏的复杂性,EnKF 方法对此类油藏的历史拟合及反演效果并不好[13],主要原因有以下两点:一是河流相油藏渗透率、孔隙度等参数的分布属于典型的非高斯分布;二是由于河流相储层高度的非均质性及储层构型的复杂性,导致地下油水运动十分复杂,油藏生产数据与地质模型之间具有很强的非线性关系.对于非高斯场问题常常通过数学变换将非高斯分布场转换为高斯分布场进行处理.截断的pluri-Gaussian方法[1415]、N-score变换[16]等常被用来研究沉积相和饱和度场的非高斯分布问题,但是变换后的高斯场与真实的非高斯场之间无法建立直接的有效联系.并且虽然模型参数转换为多参数高斯分布,但是油藏生产数据与实际模型参数的关系仍是强非线性关系.Jafarpour等[17]2009年将离散余弦变换(DCT)技术引入到油藏反演中,非高斯场的DCT 系数接近于高斯分布,并且通过保留低频系数重构渗透率场,保证了渗透率场的连续性.Jafarpour 将DCT 技术与标准EnKF 耦合进行了耦合,非线性问题没有很好地得到解决,但DCT 技术使得拟合和反演结果得到了改进.为解决油藏非线性问题,常常引入迭代算法.Li等[18]2007 年提出一种迭代方法,该方法采用与标准EnKF相同的方法计算协方差矩阵,但仍然利用伴随矩阵技术计算梯度,工作量太大.Gu等[19]在随机极大似然法[20](RML)的基础上提出一种集合RML 方法(EnRML),每次迭代后,所有的油藏模型都需要从时间零运行油藏数值模拟软件,计算量仍然很大.并且由于平均敏感系数矩阵的应用,在后续的迭代中不能保证搜索沿目标函数下降的方向进行,多次迭代的效果并不好.Wang 等[21]2010 年提出一种半迭代EnKF(HIEnKF)法,对油藏流体界面位置及相对渗透率进行了估计,并与EnRML 方法进行了比较.该方法充分利用了标准EnKF 法的优点,既避免了复杂的敏感系数计算,又较好地解决了非线性问题.相对于其他迭代方法,该方法计算量大幅度降低.

本文以HIEnKF 技术为基础,理论上证明了HIEnKF与EnRML 方法的等效性,将DCT 技术与HIEnKF方法进行了耦合,解决了河流相油藏储层参数生产历史反演的非线性和非高斯分布问题,充分体现了河流相油藏的复杂结构及非均质性,实现了河流相油藏的历史拟合及反演.

2 油藏渗流模型

油藏中流体的非线性渗流由下式的初边值问题表示:

(1)

其中,m表示油藏静态模型,包括油藏孔隙度场、渗透率场等;p表示油藏动态场(如压力、饱和度、溶解油气比等);p0 表示油藏初始动态场.该问题为非线性初边值问题,必须进行时间、空间离散化进行数值求解.离散模型表示如下:

(2)

其中上标n表示时间tnFd表示离散后的隐式数学模型,通过逐步求解可以得到动态场随时间、空间的变化规律.时间、空间离散后,油藏静态模型mNm×1的列向量;油藏动态场pNp×1 的列向量.式(2)也可写为

(3)

其中Fdn为离散后n时刻动态场数学模型,表示n时刻动态场完全由油藏模型与初始条件决定.通过方程(2)或(3)的求解还可以得到不同时间的油藏生产预测数据dn,其与动态场和油藏模型之间有以下非线性关系:

(4)

其中gn为油藏数值模型,表示tn时刻预测数据只与油藏地质模型和该时刻的动态场有关,dnNd×1的列向量,包括地震数据、生产油气比、含水率、井底流压等.实际上,方程(2)、(3)或(4)的求解就是运行油藏数值模拟软件的过程.

假设油田可在不同时刻t1t2,…,tn得到测量数据,tn时刻的生产测量数据dobsn 可看作是随机向量,由下式表示:

(5)

εdn 表示测量误差,假设不同时间的测量误差没有相关性,且测量误差属于高斯分布,其平均值为0,协方差矩阵为CDn.

上述油藏渗流模型为油藏渗流正问题.油藏渗流的反问题是在得到地面测量数据的条件下,反演油藏模型m或油藏初始条件p0,使得式(5)成立,也称之为油藏历史拟合.

3 集合Kalman滤波反演方法 3.1 标准集合Kalman滤波方法

集合Kalman 滤波(EnKF)是一种免梯度(Derivative-free)反演算法,相对于传统的基于梯度的算法,其计算效率大大提高;并且EnKF 方法是一种按时间递推的实时反演方法,随着生产数据的不断获得,可以实时进行反演.

定义状态向量

(6)

向量yn的维数为(Nm+Np)×1.集合Kalman滤波法的目的是得到一系列可以拟合油田生产历史数据的油藏模型,对于高斯分布模型,通过对下面目标函数进行极小值求解得到可以拟合生产测量数据的油藏模型集合:

(7)

其中,Ne表示集合中的模型数,yjnp表示在第n个拟合时间步的第j个模型的预测状态向量,由tn-1时刻的更新状态向量运行油藏模拟器至tn时刻得到,上标T 表示矩阵的转置.duc, jn由下式表示

(8)

其中测量误差协方差矩阵CDn=LDLDTZDNd×1维的正态分布的无相关的随机数,其平均值为0,标准偏差为1.相应的协方差矩阵由下式计算:

(9)

其中Ynp的第j列表示第j个模型的预测状态向量,Ynp的每一列都表示集合预测状态向量的平均值ynpΔYnp的第j列表示第j个模型与集合平均模型的偏差.假设G 表示集合平均模型ynp处的敏感系数矩阵,进行泰勒展开有

(10)

对式(7)利用高斯-牛顿法求极值,取第一次迭代,并将式(9)、(10)代入,得到标准集合Kalman 滤波方法的分析方程

(11)

其中,上标up分别表示更新变量和预测变量;CDnDn表示预测生产数据dnp的自协方差矩阵;CYnDn表示状态向量ynp和预测生产数据dnp的交互协方差矩阵,分别由下式估算:

(12)

(13)

其中Dnp的第j列表示第j个模型的预测生产数据,Dnp的每一列都表示集合预测生产数据的平均值dnpΔDnp的第j列表示第j个模型的预测值与平均集合预测值的偏差.CYnDn是(Nm+NpNd的大矩阵,但是在集合Kalman 滤波中该矩阵不用直接计算和存储,这也是集合Kalman 滤波方法计算效率高的一个反映.式(11)可以分为静态场、动态场分析方程:

(14)

(15)

CMnDnCPnDn分别表示模型静态场、动态场与预测生产数据dnp的交互协方差矩阵.应用集合Kalman滤波分析方程对生产数据进行拟合,就可以得到更新的状态向量,即更新的油藏模型以及更新的油藏动态参数场,然后利用这些更新的参数场运行油藏数值模拟软件到下一个需要拟合数据的时间步,进行下一次的数据拟合,标准集合Kalman 滤波过程示意图如图 1 所示.所以集合Kalman 滤波是一种按时间循序推进的数据拟合方法,集合中每一个模型只需要运行一次数值模拟软件,相对于基于梯度的历史拟合方法,极大地节省了计算量和计算要求.

图 1 标准集合Kalman滤波过程 Fig. 1 Process of the standard ensemble Kalman filter

但是标准集合Kalman滤波方法作为按时间循序推进的数据拟合方法的前提是在每一数据拟合时间步得到的更新模型与更新动态场是一致的,即利用更新的油藏模型(式(14))从时间零运行油藏数值模拟软件得到的动态场(式(3))与更新动态场(式(15))一致,即下式成立

(16)

Thulin等[22]已经证明在线性模型中该前提是正确的,但对于强非线性问题标准集合Kalman 滤波常常会得出错误的结果[21].

3.2 半迭代集合Kalman滤波

为了克服非线性油藏问题中更新模型与更新动态场的非一致性,对标准集合Kalman 滤波方法进行改进,得到半迭代集合Kalman滤波[21].半迭代集合Kalman滤波中,状态向量只包括油藏模型参数,分析方程变为

(17)

与式(14)不同的是预测数据djnp是利用油藏模型mjnp从初始状态运行油藏数值模拟软件到时间tn得到的,所有更新参数对生产数据的影响可以准确地反映出来,从而正确地调整油藏模型.拟合完成后利用所有的更新模型从时间零重新运行数值模拟软件到下一个数据拟合时间tn+1 得到djn+1,p,进行下一次数据拟合,半迭代集合Kalman 滤波方法过程示意图如图 2所示.因为只对油藏模型进行估计更新,并且每次完成数据拟合后,都要利用更新的油藏模型从时间零运行油藏数值模拟软件进行下一时刻生产预测,这样,既可准确地反映出模型对渗流的非线性影响,也保证了式(16)的成立,避免了非线性引起的更新模型与更新动态场的不一致性问题.相对于其他迭代集合Kalman 滤波法,该方法计算量要小得多,只需要在每一个数据拟合时间点对所有的油藏模型从时间零开始运行一次油藏数值模拟软件即可.

图 2 半迭代集合Kalman滤波过程 Fig. 2 Process of the haff iterative ensemble Kalman filter
3.3 半迭代集合Kalman 滤波法与EnRML 方法的等效性

可以从数学上证明半迭代集合Kalman滤波法与EnRML 方法是等效的.如果选取迭代步长为1,且第一次迭代的初始模型是先验模型,则EnRML[16]方法的第一次迭代方程如下:

(18)

其中Gnn时刻平均敏感系数矩阵的近似,由下式计算:

(19)

其中ΔMnp的第j列表示第j个模型与集合平均模型的偏差,是Nm×Ne的矩阵,ΔDnp的第j列表示第j个模型的预测值与集合平均预测值的偏差,是Nd×Ne的矩阵,与式(14)、(15)不同的是,所有预测数据都从时间零运行数值模拟软件得到.一般情况下ΔMnp是奇异阵,不可逆,进行奇异值分解:

(20)

其中UV都是正交阵,满足以下条件:

(21)

(22)

ΛNm×Ne的矩阵.一般情况下NmNe,则

(23)

λiΔMnp的第i个非零奇异值,ΛNe是对角线奇异值不为零的Ne×Ne的对角阵.所以

(24)

UNeNm×Ne的矩阵,其第j列是对应第j个非零奇异值λj的左奇异向量,VNeNe×Ne的矩阵,其第j列是对应第j个非零奇异值λj的右奇异向量.所以平均敏感系数矩阵可近似为

(25)

因为

(26)

所以

(27)

同理

(28)

将式(27)和(28)代入EnRML 第一次迭代时的分析方程(18),方程(18)变为半迭代EnKF 的分析方程(17),即半迭代EnKF 相当于EnRML 方法第一次迭代.由于平均敏感系数矩阵的应用,EnRML 在后续的迭代中不能保证搜索沿目标函数下降的方向进行,实际应用中发现,EnRML 第二次迭代的效果已经较差[21].因此半迭代EnKF 的拟合及反演结果与EnRML 结果区别不大,但是半迭代EnKF 的工作量要远远低于EnRML.

4 基于离散余弦变换的集合Kalman滤波

离散余弦变换是将一系列的离散点表示为频率不同的余弦函数的叠加.由于良好的能量压缩性能,离散余弦变换已经成功地应用于信号及图像处理领域,增加计算效率.非高斯参数场经过离散余弦变换后,DCT系数更接近于高斯分布,这符合集合Kalman滤波的基本前提.离散余弦变换的另外一个优势是当对非高斯参数场进行反演时,由于良好的能量压缩性能,可以只保留部分低频系数来代表参数场,保证参数场分布的光滑性.

设二维油藏离散渗透率场,k(ij)表示在(ij)网格处的渗透率值,沿x方向有M个网格,沿y方向有N个网格.对该渗透率场进行离散余弦变换:

(29)

其中,

D是渗透率场离散余弦变换后的DCT 系数.由式(29)可以看出,随着mn增加,离散余弦变换基函数的频率增加,即在原点处频率为0,随着mn的增加,频率增加,如图 3所示.定义DCT 变换矩阵WMWN,其元素分别定义如下

图 3 DCT 系数保留模式示意图 (a)矩形区域;(a)三角形区域. Fig. 3 The sketch map of the retained mode of the DCT coefficients (a) Rectangular area; (b) Triangular area.

式(29)可以写为矩阵的形式:

(30)

其中K表示二维渗透率场.因为DCT 变换矩阵WMWN是正交矩阵,

(31)

(32)

所以离散余弦逆变换可以写为

(33)

因此一个M×N的渗透率场K经过离散余弦变换后变成了一个M×N的DCT 系数矩阵D,经过离散余弦逆变换可以由DCT 系数D重新得到渗透率场K.在图像处理领域,由于图像是已知的,DCT 系数可以由式(30)计算得到,利用离散余弦变换良好的能量压缩性能,通过保留大系数来保存大部分的图像信息.对于油藏生产历史拟合,实际的油藏模型是未知的,无法判断较大DCT 系数的位置.但是对于一般的油藏参数场,根据离散余弦变换的性质,大部分大系数都位于低频部分,通过保留部分低频系数进行逆变换,可以得到近似的原始参数场[17].如果只是保留M1×N1 区域的DCT 系数,则进行正、逆变换时无需计算所有的系数,变换矩阵WM变为M1×M矩阵,WN变为N1×N矩阵.系数保留区域可以为矩形或三角形的情况,如图 3所示.

经过离散余弦变换后,标准集合Kalman 滤波方程变为

(34)

其中下标DCT 表示参数场离散余弦变换后的DCT系数.对于标准集合Kalman滤波,油藏静态参数场和动态场要同时进行变换,以保证模型与动态场的一致性.对于改进的集合Kalman滤波,其分析方程如下

(35)

反演过程中只对静态参数场进行离散余弦变换,然后利用所有的更新模型的DCT 系数逆变换得到更新模型,从时间零重新运行数值模拟软件到下一个数据拟合时间tn+1 得到djn+1,p,进行下一次数据拟合.

5 算例 5.1 已知渗透率场的离散余弦变换

图 4a是一假定河流相油藏渗透率场,该油藏网格数为45×45×1,每个网格大小为10m×10m×10m.其中河流相砂体渗透率为5000×10-15m2(黑色所示),泥岩渗透率为500×10-15 m2(浅灰色所示).图 4b是经过对数变换后的渗透率分布柱状图,很明显该渗透率场分布具有两个峰值,是非高斯场.而高斯分布场是EnKF 算法的前提[23],对非高斯场,EnKF 常常会得到不合理的反演结果[13],通常的做法是利用数学变换将非高斯分布场转换为高斯分布场进行处理[14~16].图 4c是该渗透率场变换后的DCT 系数分布,系数分布范围较宽,大数值的系数(黑色数据点)个数相对较少.图 4d是DCT 系数分布柱状图,经过DCT 变换后,非高斯渗透率场变成了近似高斯分布的DCT 系数,这对于基于高斯分布假设的反演理论是非常有利的.图 4c中大部分数值大的DCT 系数都集中在低频部分(原点附近),因为DCT 变换具有很好的能量压缩特性,大数值的系数包含了渗透率场的大部分信息,所以可以通过图 3所示的系数保留模式进行逆变换.

图 4 渗透率场的离散余弦变换 (a)原始渗透率场;(b)对数渗透率场分布柱状图;(c)DCT 系数分布;(d)DCT 系数柱状图. Fig. 4 Discrete cosine transform of the permeability field (a) Initial permeability field; (b) Histogram of the log permeability field; (c) Distribution of the DCT coefficients;(d) Histogram of the DCT coefficients.

图 5比较了三种不同系数保留模式对渗透率场逆变换的影响.图 5a所示是保留300 个最大DCT系数逆变换得到的渗透率场,可以看出逆变换渗透率场与原始渗透率场已非常接近.由于这300 个系数包含一些高频基(见图 4c),反演得到的渗透率场中同一沉积相内渗透率分布非常不均匀,网格间渗透率变化较大.图 5b是保留M1×N1=20×20的矩形区域内的系数(图 4c所示实线包围的矩形区域)逆变换得到的渗透率场,虽然有部分数值较大的系数在该区域之外,但是大部分大数值的系数在该区域之内,所以逆变换得到的结果仍然较好.因为该区域包含了部分数值较低的系数,虽然保留的系数数目比图 5a所示情况要多,但由于丢失了部分信息,逆变换效果明显比图 5a差,砂岩边界变得有些模糊.但是由于只保留低频系数,在同一沉积相内渗透率的分布变得相对均匀.图 5c是保留M1 ×N1 =20×20的三角形区域内的系数(图 4c所示虚线包围的三角形区域)反演得到的渗透率场,与图 5b相比,保留系数减少一半,逆变换渗透率场进一步变差,砂岩边界越来越模糊,但是反演渗透率场仍然反映了原始渗透率场的主要结构与特征.因为频率进一步降低,在同一沉积相内渗透率变得更为均匀.

图 5 DCT 系数保留模式对离散余弦拟变换的影响 (a)最大系数(300);(b)M1×N1=20×20(矩形);(c)M1×N1=20×20(三角形) Fig. 5 Effect of the DCT coefficients retaining mode on discrete cosine transform (a) Retaining largest coefficients; (b) M1X N1 = 20 X 20(rectangle area);(c) M1X N1 = 20 X 20 (triangle area).

在油藏历史拟合过程中,真实的渗透率场是未知的,因此无法准确知道DCT 系数的分布,通过保留一些大数值系数来表示渗透率场的方式是不可行的.根据离散余弦变换的性质[24]及上面讨论结果,可以通过保留低频区域的DCT 系数来表示原渗透率场,无论保留低频部分矩形区域或者三角形区域的DCT 系数,都能反映原始渗透率场的主要结构与特征.并且忽略了高频基,使得同一沉积相中的渗透率分布更光滑.

5.2 耦合离散余弦变换的EnKF 法和HIEnKF 法对河流相油藏的反演

利用标准EnKF、耦合DCT 的EnKF 及耦合DCT 的HIEnKF方法分别对一假定河流相油藏进行生产历史拟合来反演其渗透率场.该油藏网格数为45×45×1,每个网格大小为10m×10m×10m, 真实渗透率场如图 6a所示.该油藏共4 口生产井,4口注入井.利用多点地质统计软件Snesim[25]产生200个先验模型,在构造先验模型时使用了井点网格渗透率硬数据,其先验平均模型如图 6b所示,训练图像如图 6c所示.图 7表示了3个先验模型的渗透率场,除井点符合渗透率硬数据外,其余网格的渗透率具有很大的差异.图 8 表示了200 个先验模型对P-2井的产水量(qw)的生产预测曲线,可以看出,先验模型不能拟合生产测量数据.并且生产预测曲线具有很大的展布,这既说明了先验模型的差异,也说明了先验模型集合具有较大的自由度用以保证取得较好的历史拟合及反演结果.该河流相油藏水驱开采1080天,注入井的测量数据为井底流动压力,生产井的测量数据为油、水产量.通过对这些生产数据的历史拟合来反演该河流相油藏渗透率场分布.

图 6 油藏基本信息 (a)真实模型;(b)先验平均模型;(c)训练模型. Fig. 6 Basic information of the reservoir (a) Real permeability field ; (b) Prior mean field of the permeability; (c) Training image.
图 7 先验油藏渗透率场 (a)模型1;(b)模型20;(c)模型60. Fig. 7 Prior reservoir permeability fields (a) Realization 1 ; (b) Realization 20 ; (c) Realization 60.
图 8 先验油藏模型集合预测的P-2井产水量 Fig. 8 Water production prediction of P-2 from prior reservoir ensemble

将离散余弦变换分别与标准EnKF法和HIEnKF法进行耦合,通过拟合各井的生产历史数据反演该河流相油藏的渗透率场分布.图 9a为未耦合DCT的标准EnKF反演结果,可以看出经过历史拟合得到的渗透率场反映了实际渗透率场的主要结构.但是砂岩与泥岩的边界却比较模糊,位置也与真实渗透率场不同,并且在砂岩与泥岩内部,渗透率分布不均匀.主要原因有以下两点,一是河流相油藏非线性渗流问题的影响,标准EnKF 在反演过程中的更新渗透率场与更新动态场不一致;第二点是EnKF 处理非高斯场会产生较大的误差.图 9b 为耦合DCT的标准EnKF反演结果,虽然反演得到的渗透率场的主要结构与图 9a相同,但是由于耦合DCT 后,只保留了M1 ×N1 =20×20 矩形区域低频部分的DCT系数,所以无论在高渗透区域还是在低渗透区域,其连续性得到了很大的改进.图 9c是耦合DCT的HIEnKF的反演结果,HIEnKF 方法由于每次拟合数据时都利用更新模型从时间零运行油藏数值模拟软件,保证了反演过程中估计渗透率场与估计动态场的一致性,其反演结果明显得到了改善,砂岩位置接近于真实渗透率场,由于DCT 的作用,各沉积相分布相对比较均匀,砂岩与泥岩边界也相对清晰.

图 9 渗透率场的反演 (a)EnKF;(b)EnKF+DCT;(c)HIEnKF + DCT. Fig. 9 Estimation of the permeability field(a)EnKF;b)EnKF+DCT;c)IIIEnKF + DCT. 也可以从生产历史的拟合曲线看地质模型的反演结果.图 10表示了利用最终的反演渗透率场从时间零重新运行油藏数值模拟软件得到的P-2井产水量的拟合曲线.可以看出,由于标准EnKF 对非高斯渗透率场的反演结果不是很好,其历史拟合效果较差,后期预测曲线明显偏离了实际生产曲线与测量数据.与DCT 耦合后,只保留M1×N1=20×20矩形区域低频部分的DCT 系数,非高斯渗透率场反演结果得到改进,历史拟合曲线(图 10b)相对于未耦合DCT 的标准EnKF 结果也得到了改进,拟合偏差明显减小.耦合DCT 的HIEnKF 拟合结果如图 10c所示,该方法既解决了油藏渗流的非线性,又改进了非高斯场的反演结果,拟合结果进一步改善,P-2井产水量可以得到很好的拟合.
图 10 中文P-2井产水量的拟合标题 (a)EnKF; (b)EnKF + DCT; (c) IIIEnKF + DCT. Fig. 10 History matching of the water production rate of P-2 (a)EnKF;b)EnKF + DCT;(c)IIIEnKF + DCT.
5.3 DCT系数保留模式对反演结果的影响

利用耦合DCT 的HIEnKF 方法讨论DCT 系数保留模式对渗透率场反演结果的影响.图 11讨论了保留不同矩形区域的DCT 系数时渗透率场的反演结果.可以看出当M1×N1=10×10时,虽然只是保留了4.9%的DCT 系数,但渗透率场反演结果已相当不错,砂体位置、边界比较清晰,同一沉积相内渗透率的分布也比较均匀,但部分砂体缺失比较厉害,比如P-2、P-3 井之间的砂体不明显.当M1×N1=20×20 时,随着系数保留的数量越来越多(占总系数的19.75%),渗透率场反演结果越来越好,砂体位置、边界愈加清晰.随着系数保留的数量越来越多,频率也越来越高,但同一沉积相内渗透率的分布仍然比较均匀.当M1 ×N1 =30×30,已经保留了44.44%的系数,虽然河道结构与保留系数较少时相同,但由于频率越来越高,渗透率的分布开始变得粗糙.如果保留全部的DCT 系数,则反演结果与HIEnKF反演结果一致.图 12是各种DCT系数保留模式对P-2井产水量的拟合曲线,可以看出各种系数保留模式对生产曲线的拟合基本一致,都能对生产测量数据进行较好的拟合.

图 11 DCT 系数保留数量对渗透率场反演的影响(矩形区域)(a)M1×N1= 10×10;b)M1×N1 = 20×20;(c)M1×N1=30×30. Fig. 11 Effect of the number of the retaining DCT coefficients on inversion of the permeability field (rectangular area)
图 12 DCT 系数保留数量对P-2井产水量拟合的影响(矩形区域) (a)MiXNi = 10X10;b)MiXNi = 20X20;c)MiXNi=30X30. Fig. 12 Effect of the number of the retaining DCT coefficients on matching of the water production rate of P-2 (rectangular area)

图 13为保留不同三角形区域的DCT 系数时渗透率场的反演结果.可以看出当M1×N1=10×10时,虽然保留DCT 系数数量只是相同情况下矩形区域的一半,但渗透率场反演结果仍然不错,甚至比矩形区域的情况(图 11a)还要好.当M1×N1=20×20和30×30,与矩形区域的情况比较可以看出,由于保留较少的DCT 系数,且系数的频率更低,同一沉积相内渗透率的分布更加均匀.图 14是各种三角形DCT 系数保留模式对P-2井产水量的拟合曲线,可以看出各种系数保留模式对生产曲线的拟合基本一致,都能对生产测量数据进行较好的拟合.与图 12所示的矩形系数保留模式的拟合曲线相比较,拟合结果基本相似.

图 13 DCT 系数保留数量对渗透率场反演的影响(三角形区域) M1×N1= 10×10;b)M1×N1 = 20×20;(c)M1×N1=30×30. Fig. 13 Effect of the number of the retaining DCT coefficients on inversion of the permeability field (triangular area)
图 14 DCT 系数保留数量对P-2井产水量拟合的影响(三角形区域) M1×N1= 10×10;b)M1×N1 = 20×20;(c)M1×N1=30×30. Fig. 14 Effect of the number of the retaining DCT coefficients on matching of the water production rate of P-2 (triangular area)

由此可以看出,DCT 变换后,三角形区域的DCT 系数保留模式对参数场的反演比较有利,保留2.5%~10%的DCT 系数,既可以降低EnKF 方法对硬件的要求,又可以保留足够多的地质信息,保证反演结果的正确,并且由于只保留低频系数部分,同一沉积相内渗透率的分布变得比较均匀,较好地反映了非高斯场的特征.

6 结语

理论证明了半迭代集合Kalman 滤波方法与EnRML 方法的等效性.将离散余弦变换与半迭代集合Kalman滤波相耦合对河流相油藏进行历史拟合与反演.离散余弦变换将非高斯分布的渗透率场变换为高斯分布的DCT 系数;利用离散余弦的能量压缩特性对DCT 系数进行截断,既降低了算法对硬件的要求,也降低了在求解过程中产生病态矩阵的几率;保留频率较低的系数,使得反演渗透率场的分布变得均匀,较好地解决了河流相油藏的非高斯分布问题;半迭代集合Kalman 滤波很好地解决了油藏生产数据与地质模型之间的非线性关系.油藏反演中可以通过矩形及三角形模式对DCT 参数进行保留,其中三角形模式效果更好,保留系数数量不宜太多,2.5%~10%的比例比较合适,既降低了EnKF 方法对硬件的要求,又可以保留足够多的地质信息,保证反演结果的正确,较好反映非高斯场的特征.

致谢

本研究部分内容是笔者在美国Tulsa大学石油工程系做博士后期间在Reynolds博士的指导下完成的,在此深表谢意!

参考文献
[1] 刘福平, 杨长春. 渗透率场敏感系数的数值计算. 地球物理学报 , 2003, 46(1): 131–137. Liu F P, Yang C C. The numerical calculation of sensitivity coefficients of permeability field. Chinese J. Geophys. (in Chinese) , 2003, 46(1): 131-137.
[2] Gao G H. Data integration and uncertainty evaluation for large scale automatic history matching. . USA: Petroleum Engineering Department of Tulsa University, 2005
[3] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys Res , 1994, 99(C5): 10143-10162.
[4] Naevdal G, Johnsen L M, Aanonsen S I, et al. Reservoir monitoring and continuous model updating using ensemble Kalman filter. SPE Annual Technical Conference and Exhibition, Transactions, SPE Paper 84372, October 2003
[5] Haugen V, Natvik L J, Evensen G, et al. History matching using the ensemble Kalman filter on a North Sea field case. SPE Annual Technical Conference and Exhibition, Transactions, SPE Paper 102430, September 2006
[6] Evensen G, Hove J, Meisingset H C, et al. Using the EnKF for assisted history matching of a North Sea reservoir model. SPE Reservoir Simulation Symposium, SPE Paper 106184, February 2007
[7] Seiler A, Evensen G, Skjervheim J A, et al. Advanced reservoir management workflow using an EnKF based assisted history matching method. SPE Reservoir Simulation Symposium, SPE Paper 118906, February 2009
[8] 秦军, 阎广建, 刘绍民, 等. 集合卡曼滤波在遥感反演地表参数中的应用—以核驱动模型反演BRDF为例. 中国科学D辑(地球科学) , 2005, 35(8): 790–798. Qin J, Yan G J, Liu S M, et al. Application of EnKF in the inversion of land surface parameter by the Remote Sensing-Kernal-Driven Model BRDF study. Sci. China Ser. D (Earth Sci.) (in Chinese) , 2005, 35(8): 790-798.
[9] 韩旭军, 李新. 非线性滤波方法与陆面数据同化. 地球科学进展 , 2008, 23(8): 813–820. Han X J, Li X. Review of the nonlinear filters in the land data assimilation. Adv. Earth Sci. (in Chinese) , 2008, 23(8): 813-820.
[10] 赵红亮, 冯夏庭, 张东晓, 等. 岩土力学参数空间变异性的集合Kalman滤波估值. 岩土力学 , 2007, 28(10): 2219–2223. Zhao H L, Feng X T, Zhang D X, et al. Spatial variability of geomechanical parameter estimation via ensemble Kalman filter. Rock and Soil Mechanics (in Chinese) , 2007, 28(10): 2219-2223.
[11] 贾炳浩, 谢正辉, 田向军, 等. 基于微波亮温及集合Kalman滤波的土壤湿度同化方案. 中国科学D辑(地球科学) , 2010, 40(2): 239–251. Jia B H, Xie Z H, Tian X J, et al. A soil moisture assimilation scheme based on the ensemble Kalman filter using microwave brightness temperature. Sci. China Ser. D (Earth Sci.) (in Chinese) , 2010, 40(2): 239-251.
[12] 李晓静, 程时清, 李彦尊, 等. EnKF在油藏描述及开发优化中的应用研究与进展. 内蒙古石油化工 , 2009, 35(4): 1–5. Li X J, Cheng S Q, Li Y Z, et al. Review on reservoir description and optimization problems with ensemble Kalman filter. Inner Mongolia Petrochemical Industry (in Chinese) , 2009, 35(4): 1-5.
[13] Aanonsen S I, Nvdal G, Oliver D S, et al. The ensemble Kalman filter in reservoir engineering-a review. SPE Journal , 2009, 14(3): 393-412. DOI:10.2118/117274-PA
[14] Agbalaka C C, Oliver D S. Application of the EnKF and localization to automatic history matching of facies distribution and production data. Mathematical Geostatistics , 2008, 40(4): 353-374.
[15] Zhao Y. Ensemble Kalman Filter Method for Gaussian and Non-Gaussian Priors. USA: Petroleum Engineering Department of Tulsa University, 2008 .
[16] Gu Y Q, Oliver D S. The ensemble Klaman filter for continuous updating of reservoir simulation models. Journal of Energy Resources Technology , 2006, 128(1): 79-87. DOI:10.1115/1.2134735
[17] Jafarpour B, Mclaughlin D B. Reservoir characterization with the discrete cosine transform. SPE Journal , 2009, 14(1): 182-201. DOI:10.2118/106453-PA
[18] Li G M, Reynolds A C. Iterative ensemble Kalman filter for data assimilation. SPE Journal , 2009, 14(3): 496-505. DOI:10.2118/109808-PA
[19] Gu Y Q, Oliver D S. An iterative ensemble Kalman filter for multiphase fluid flow data assimilation. SPE Journal , 2007, 12(4): 438-446. DOI:10.2118/108438-PA
[20] Oliver D S, Reynolds A C, Liu N. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge: Cambridge University Press, 2008 .
[21] Wang Y D, Li G M, Reynolds A C. Estimation of depths of fluid contacts and relative permeability curves by history matching using iterative ensemble Kalman smoothers. SPE Journal , 2010, 15(2): 509-525. DOI:10.2118/119056-PA
[22] Thulin K, Li G M, Aanonsen S I, et al. Estimation of initial fluid contacts by assimilation of production data with EnKF. SPE Annual Technical Conference and Exhibition, Transactions, SPE Paper 109975, November 2007
[23] Evensen G. Data assimilation: the Ensemble Kalman Filter. Berlin: Springer, 2007 .
[24] Gilbert S. The discrete cosine transform. SIAM Review , 1999, 41(1): 135-147. DOI:10.1137/S0036144598336745
[25] Strebelle S B, Journel A G. Reservoir modeling using multiple-point statistics. SPE Annual Technical Conference and Exhibition, Transactions, SPE Paper 71324, September 2001