CN106503456B - 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法 - Google Patents

基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法 Download PDF

Info

Publication number
CN106503456B
CN106503456B CN201610943462.6A CN201610943462A CN106503456B CN 106503456 B CN106503456 B CN 106503456B CN 201610943462 A CN201610943462 A CN 201610943462A CN 106503456 B CN106503456 B CN 106503456B
Authority
CN
China
Prior art keywords
state
prediction
matrix
state vector
production data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610943462.6A
Other languages
English (en)
Other versions
CN106503456A (zh
Inventor
梁晓东
李太福
周伟
易军
张元涛
刘媛媛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing Juxin Intelligent Technology Research Institute Co ltd
Original Assignee
Chongqing University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing University of Science and Technology filed Critical Chongqing University of Science and Technology
Priority to CN201610943462.6A priority Critical patent/CN106503456B/zh
Publication of CN106503456A publication Critical patent/CN106503456A/zh
Application granted granted Critical
Publication of CN106503456B publication Critical patent/CN106503456B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16ZINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
    • G16Z99/00Subject matter not provided for in other main groups of this subclass

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,包括:步骤1:初始化油藏模型的集合,集合包括静态参数、动态参数以及油井生产数据;步骤2:将静态参数中的渗透率进行超球体变换,并构造新的状态向量集合;步骤3:将新的状态向量集合中的每个状态向量输入到油藏模拟器中进行预测获得状态预测值,每个状态向量及其状态预测值构成预测集合;步骤4:根据预测集合计算卡尔曼增益矩阵;步骤5:根据预测集合、卡尔曼增益矩阵以及测量的油井生产数据对预测集合进行更新,获得更新后的静态参数、动态参数以及油井生产数据。利用本发明能够提高油藏历历史拟合精度的准确性和减少人工历史拟合的盲目性。

Description

基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法
技术领域
本发明涉及油田开发技术领域,具体涉及一种基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法。
背景技术
油藏历史拟合是油田生产过程中不可缺少的环节之一。在石油工业中,闭环油藏管理的概念在智能油田的概念上得到了油田工作者的极大关注,对生产制度的制定、最大限度的开发现有油气资源有着重要的意义。传统的油藏历史拟合方法为油藏工程师根据其经验来调整模型参数,通过油藏模拟器运算的结果来判定设定的参数能否与生产历史数据相吻合,最后用能够最佳拟合历史生产数据的一组模型参数对未来的生产数据进行预测。
在油藏历史拟合过程中需要调整的参数多,模型的自由度较大,传统的人工历史拟合方法非常耗时耗力,基于梯度的历史拟合方法需要计算灵敏度系数,难以对不确定性进行评估。
发明内容
本发明通过提供一种基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,以解决上述背景技术中所提出的问题。
本发明提供的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,包括:
步骤S1:初始化油藏模型的集合;其中,集合由状态向量的矩阵组成,状态向量由状态变量组成,状态变量包括静态参数、动态参数和油井生产数据;静态参数包括油藏模型每个网格的渗透率和孔隙度,动态参数包括油藏模型每个网格的压力和饱和度,油井生产数据包括井底压力、油井产油量和油井产水量;
步骤S2:对集合中静态参数的渗透率进行超球体变换,获得Sigma点集与孔隙度组成新的静态参数;其中,为nK×(nK+2)维矩阵,表示集合中第i个静态参数的渗透率经超球体变换之后第j维的第l个粒子,第l个粒子的第一列为Ki,nK为每组渗透率的维数;
步骤S3:基于新的静态参数、动态参数和油井生产数据构造k时刻的状态向量集合
步骤S4:将状态向量集合中的每个状态向量输入到油藏模拟器,以对每个状态向量k+1时刻的值进行预测,获得每个状态向量在k+1时刻的状态预测值,所有的状态向量及其状态预测值形成预测集合,同时计算预测集合中每个状态向量的协方差矩阵;
步骤S5:根据预测集合和预测集合中每个状态向量的协方差矩阵计算卡尔曼增益矩阵Gk+1
步骤S6:根据卡尔曼增益矩阵Gk+1和测量的油井生产数据dk+1,对预测集合中每个状态向量的状态预测值进行更新,获得更新后的静态参数、动态参数以及油井生产数据。
与现有技术相比,本发明提供的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法的优点是:通过对静态参数中的渗透率进行超球体转换,以提高油藏历历史拟合精度的准确性,以及,集合卡尔曼滤波算法具有计算效率高、鲁棒性强的特点,使油井生产数据可以持续、实时、快速地被吸收,可缩短计算周期,减少人工历史拟合的盲目性。
具体实施方式
本发明的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,包括:
步骤S1:初始化油藏模型的集合。
初始化形成的油藏模型的集合为:
其中,xn,j表示在时刻tn的状态矢量的第j个集合元素。ms和md分别为静态参数和动态参数;其中,静态参数包括油藏模型每个网格的渗透率和孔隙度,动态参数包括油藏模型每个网格的含水饱和度和压力;d为油井生产数据,包括井底压力、油井产油量和油井产水量。
静态参数、动态参数和油井生产数据构成状态变量,状态变量组成的向量为状态向量,所有状态向量组成的矩阵为油藏模型的集合。
步骤S2:对集合中静态参数的渗透率进行超球体变换,获得Sigma点集与孔隙度组成新的静态参数;其中,为nK×(nK+2)维矩阵,表示集合中第i个静态参数的渗透率经超球体变换之后第j维的第l个粒子,第l个粒子的第一列为Ki,nK为每组渗透率的维数。
对集合中静态参数的渗透率进行超球体变换的过程,包括:
步骤S21:随机选择0≤W0≤1。
W0为计算Sigma点的权值时的初始权重。
步骤S22:计算第i个静态变量第j维对应的Sigma点的权值。
计算Sigma点权值的公式为:
步骤S23:初始化向量序列。
其中,初始化公式为:
步骤S24:当输入维数j=2,3,L,n时,对初始化的向量序列进行迭代生成Sigma点集。
其中,迭代公式为:
其中,为第j维的第i个粒子点;
步骤S25:对Sigma点集中的每个Sigma点加入状态变量的均值和协方差矩阵Pxx
Sigma点集中的每个Sigma点在加入状态变量的均值和协方差矩阵Pxx后变为
从以上采样算法可以看出,除了原点以外其他采样点具有相同的权值,而且都位于半径为的超球体上。
步骤S3:基于新的静态参数、动态参数和油井生产数据构造k时刻的状态向量集合
构造的k时刻的状态向量集合为:
其中,φi为每个静态参数的孔隙度;0和φi分别为nK×(nK+1)维矩阵,a第k时刻第i个状态向量的更新值。
步骤S4:将状态向量集合中的每个状态向量输入到油藏模拟器,以对每个状态向量k+1时刻的值进行预测,获得每个状态向量在k+1时刻的状态预测值,所有的状态向量及其状态预测值形成预测集合,同时计算预测集合中每个状态向量的协方差矩阵。
将状态向量集合中的每个状态向量输入到油藏模拟器,以对每个状态向量k+1时刻的值进行预测的公式为:
其中,为状态预测值,F为油藏模拟器,对中的进行加权得到:
每个状态向量在k+1时刻的状态预测值为:
其中,为加权后的状态预测值。
预测集合中状态预测值的均值为:h为预测集合中状态向量的数量;
预测集合中每个状态向量的协方差矩阵为:
步骤S5:根据预测集合和预测集合中每个状态向量的协方差矩阵计算卡尔曼增益矩阵Gk+1
根据预测集合和所述预测集合中每个状态向量的协方差矩阵计算卡尔曼增益矩阵Gk+1的公式为:
其中,H为状态变量与状态变量中的生产数据相关联的矩阵算子,表示为的形式,0为全部元素都为0的nd,k+1×(nx,k+1-nd,k+1)矩阵,为nd×nd的单位矩阵,Cd,k+1为油井生产数据误差的协方差矩阵,其维数为nd,k+1×nd,k+1
步骤S6:根据卡尔曼增益矩阵Gk+1和预测的油井生产数据dk+1,对预测集合中每个状态向量的状态预测值进行更新,获得更新后的静态参数、动态参数以及油井生产数据。
根据卡尔曼增益矩阵Gk+1和预测的生产数据dk+1,对预测集合中每个状态向量的状态预测值进行更新的公式为:
预测集合中更新后的状态预测值的均值为:
预测集合中每个更新后状态预测值所对应的状态向量的协方差矩阵为:

Claims (6)

1.一种基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,包括如下步骤:
步骤S1:初始化油藏模型的集合;其中,所述集合由状态向量的矩阵组成,所述状态向量由状态变量组成,所述状态变量包括静态参数、动态参数和油井生产数据;所述静态参数包括所述油藏模型每个网格的渗透率和孔隙度,所述动态参数包括所述油藏模型每个网格的压力和饱和度,所述油井生产数据包括井底压力、油井产油量和油井产水量;
步骤S2:对所述集合中静态参数的渗透率进行超球体变换,获得Sigma点集与所述孔隙度组成新的静态参数;其中,为nK×(nK+2)维矩阵,表示所述集合中第i个静态参数的渗透率经超球体变换之后第j维的第l个粒子,第l个粒子的第一列为Ki,nK为每组渗透率的维数;
步骤S3:基于新的静态参数、所述动态参数和所述油井生产数据构造k时刻的状态向量集合
步骤S4:将所述状态向量集合中的每个状态向量输入到油藏模拟器,以对每个状态向量k+1时刻的值进行预测,获得每个状态向量在k+1时刻的状态预测值,所有的状态向量及其状态预测值形成预测集合,同时计算所述预测集合中每个状态向量的协方差矩阵;
步骤S5:根据所述预测集合和所述预测集合中每个状态向量的协方差矩阵计算卡尔曼增益矩阵Gk+1
步骤S6:根据所述卡尔曼增益矩阵Gk+1和测量的k+1时刻的油井生产数据dk+1,对所述预测集合中每个状态向量的状态预测值进行更新,获得更新后的静态参数、动态参数以及油井生产数据。
2.根据权利要求1所述的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,对所述集合中静态参数的渗透率进行超球体变换的过程,包括:
步骤S21:随机选择0≤W0≤1;其中,W0为计算Sigma点的权值时的初始权重;
步骤S22:计算第i个静态变量第j维对应的Sigma点的权值;其中,
步骤S23:初始化向量序列;其中,初始化公式为:
步骤S24:当输入维数j=2,3,L,n时,对初始化的向量序列进行迭代生成所述Sigma点集;其中,迭代公式为:
其中,为第j维的第i个粒子点;
步骤S25:对所述Sigma点集中的每个Sigma点加入状态变量的均值和协方差矩阵Pxx后变为:
3.根据权利要求2所述的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,构造的k时刻的状态向量集合为:
其中,0和φi分别为nK×(nK+1)维矩阵,其中,0表示一个矩阵,矩阵里面的元素全部是0,φi为nK×(nK+1)维矩阵,nK×(nK+1)为φi矩阵的维度,而φi为每个静态参数的孔隙度;ms为所述静态参数;md为所述动态参数;d为k时刻的油井生产数据;di表示第i个状态的油井生产数据;md,i表示第i个状态向量的动态参数;a为第k时刻第i个状态向量的更新值;表示渗透率Ki经过超球体变换产生的渗透率值,为nK×(nK+1)的矩阵。
4.根据权利要求3所述的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,将所述状态向量集合中的每个状态向量输入到油藏模拟器,以对每个状态向量k+1时刻的值进行预测的公式为:
其中,为状态预测值,F为油藏模拟器,对中的进行加权得到:
其中,表示k+1时刻对状态向量进行变换后第i个静态参数第j维的第l个粒子;Wj表示第i个静态变量第j维对应的Sigma点的权值;
每个状态向量在k+1时刻的状态预测值为:
其中,为加权后的状态预测值;Ki'表示k+1时刻第i个静态参数的渗透率;
所述预测集合中状态预测值的均值为:h为所述预测集合中状态向量的数量;
所述预测集合中每个状态向量的协方差矩阵为:
Ck+1 f表示k+1时刻状态的预测协方差矩阵。
5.根据权利要求4所述的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,根据所述预测集合和所述预测集合中每个状态向量的协方差矩阵计算卡尔曼增益矩阵Gk+1的公式为:
其中,H为状态变量与状态变量中的生产数据相关联的矩阵算子,表示为的形式,0为全部元素都为0的nd,k+1×(nx,k+1-nd,k+1)矩阵,为nd×nd的单位矩阵,Cd,k+1为油井生产数据误差的协方差矩阵,其维数为nd,k+1×nd,k+1;nd,k+1表示k+1时刻生产数据d的维数;nx,k+1表示k+1时刻状态向量x的维数。
6.根据权利要求5所述的基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法,其特征在于,根据所述卡尔曼增益矩阵Gk+1和测量的K+1时刻的油井生产数据dk+1,对所述预测集合中每个状态向量的状态预测值进行更新的公式为:
其中,dk+1,i表示第i个状态向量k+1时刻的生产数据;
所述预测集合中更新后的状态预测值的均值为:
所述预测集合中每个更新后状态预测值所对应的状态向量的协方差矩阵为:
CN201610943462.6A 2016-10-26 2016-10-26 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法 Expired - Fee Related CN106503456B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610943462.6A CN106503456B (zh) 2016-10-26 2016-10-26 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610943462.6A CN106503456B (zh) 2016-10-26 2016-10-26 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法

Publications (2)

Publication Number Publication Date
CN106503456A CN106503456A (zh) 2017-03-15
CN106503456B true CN106503456B (zh) 2019-02-22

Family

ID=58321085

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610943462.6A Expired - Fee Related CN106503456B (zh) 2016-10-26 2016-10-26 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法

Country Status (1)

Country Link
CN (1) CN106503456B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109214013A (zh) * 2017-06-29 2019-01-15 中国石油化工股份有限公司 一种集合卡尔曼滤波方法及装置
CN108038330B (zh) * 2017-12-26 2022-02-08 重庆科技学院 基于sukfnn算法的铝电解工耗模型构建方法
CN113033862B (zh) * 2020-01-10 2023-09-22 重庆科技学院 稠油油藏超临界多元热流体吞吐产能预测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101356537A (zh) * 2005-11-21 2009-01-28 切夫里昂美国公司 使用集合卡尔曼滤波进行实时油藏模型更新的方法、系统和装置
CN104750896A (zh) * 2013-12-31 2015-07-01 中国石油化工股份有限公司 一种缝洞型碳酸盐岩油藏数值模拟方法
CN105808311A (zh) * 2014-12-29 2016-07-27 中国石油化工股份有限公司 一种基于降维策略的油藏模拟快速拟合方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101356537A (zh) * 2005-11-21 2009-01-28 切夫里昂美国公司 使用集合卡尔曼滤波进行实时油藏模型更新的方法、系统和装置
CN104750896A (zh) * 2013-12-31 2015-07-01 中国石油化工股份有限公司 一种缝洞型碳酸盐岩油藏数值模拟方法
CN105808311A (zh) * 2014-12-29 2016-07-27 中国石油化工股份有限公司 一种基于降维策略的油藏模拟快速拟合方法

Also Published As

Publication number Publication date
CN106503456A (zh) 2017-03-15

Similar Documents

Publication Publication Date Title
CN110705029B (zh) 一种基于迁移学习的振荡扑翼能量采集系统流场预测方法
CN113378939B (zh) 基于物理驱动神经网络的结构数字孪生建模与参数识别法
CN105808311B (zh) 一种基于降维策略的油藏模拟快速拟合方法
CN102831269A (zh) 一种流程工业过程工艺参数的确定方法
CN103730006A (zh) 一种短时交通流量的组合预测方法
CN110224862B (zh) 基于多层感知器的多智能体系统网络容侵能力评估方法
CN104537033A (zh) 基于贝叶斯网络和极限学习机的区间型指标预报方法
CN104732303A (zh) 一种基于动态径向基函数神经网络的油田产量预测方法
CN102982250B (zh) 基于随机响应面估计参数不确定性的随机模型修正方法
CN106503456B (zh) 基于超球体变换的集合卡尔曼滤波油藏动态历史拟合方法
CN110262245A (zh) 基于事件触发机制的多智能体系统的控制器设计方法
CN106355003B (zh) 基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及系统
CN105425583B (zh) 基于协同训练lwpls的青霉素生产过程的控制方法
CN109829217A (zh) 压裂性裂缝油藏产能模拟方法及装置
CN111308896A (zh) 基于可变误差的非线性系统自适应最优控制方法
CN110807544A (zh) 一种基于机器学习的油田剩余油饱和度分布的预测方法
CN112016253A (zh) 一种适用于cfd不确定度量化的高保真度混沌多项式修正方法
Grema et al. Optimal feedback control of oil reservoir waterflooding processes
CN105678015B (zh) 一种高超声速三维机翼的非概率可靠性气动结构耦合优化设计方法
CN115375031A (zh) 一种产油量预测模型建立方法、产能预测方法及存储介质
CN109736720B (zh) 一种基于改进Kriging模型的深海连接器密封结构优化方法
CN112113146B (zh) 供水管网管道粗糙系数和节点需水量同步自适应校核方法
CN103258131A (zh) 基于正交学习粒子群的功率电路元件优化方法
CN118551290B (zh) 一种基于深度神经算子的水驱油藏生产动态预测方法
CN118734457A (zh) 一种基于多保真流形神经算子的复合材料构件固化状态场预测方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20210113

Address after: Room 5-2, No.1-5, Tieshan Road, Biquan street, Bishan District, Chongqing 402760

Patentee after: Chongqing Juxin Intelligent Technology Research Institute Co.,Ltd.

Address before: No. 20, East Road, University City, Chongqing, Shapingba District, Chongqing

Patentee before: Chongqing University of Science & Technology

TR01 Transfer of patent right
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190222

CF01 Termination of patent right due to non-payment of annual fee