CN113378469B - 一种基于卡尔曼滤波与支持向量机的测井曲线预测方法 - Google Patents

Abstract

本发明涉及一种基于卡尔曼滤波与支持向量机的测井曲线预测方法,包括以下步骤：首先选择足够数量的完备井作为标准井，并对测井曲线进行标准化采样处理，建立状态空间模型，其次，基于熵离散化测井曲线，建立状态迁移变换矩阵和服从高斯分布的概率密度函数，接着，基于支持向量回归方法建立三次多项式非线性观测函数，并建立服从高斯分布的概率密度函数，最后，基于状态空间模型和卡尔曼滤波算法预测测井曲线。本发明方法具有快速有效、预测测井曲线准确性更高的优点，有效解决了常规利用回归方法无法有效预测测井曲线空间序列所蕴含的特征和结构等问题。

Description

一种基于卡尔曼滤波与支持向量机的测井曲线预测方法

技术领域：

本发明涉及油气井的测井曲线数据预测方法，尤其是一种快速有效、预测准确性更高的的基于卡尔曼滤波与支持向量机的测井曲线预测方法。

背景技术：

随着油气开采技术的不断进步，在钻井过程中一般采用专业设备测量采集包括声波曲线AC、自然伽马曲线GR、中子曲线CNL、电阻率曲线RT、井径曲线CAL等相关测井曲线数据，以便后期对油气田所处区域地质结构、油气储层以及油气开采研究提供基础数据。然而，一个油气田的老井，可能没有获取收集测井曲线数据，或者由于技术条件限制，获取收集的测井曲线数据不全，又由于老井在钻井过程中，都基本采用三层套管，套管之间采用砼浇筑，在钻井过程中没有采集的测井曲线数据，后期无法再进行采集。面对如此情况，现在一般基于测井曲线完备的井使用回归方法建立预测模型，实现对缺失曲线的预测，然而这些预测模型没有太多地关注测井曲线作为空间序列所蕴含的特征与结构，预测效果不是特别的理想，因此需要一种更加有效的方法提高预测缺失测井曲线的准确性。

发明内容：

为了解决上述技术问题，本发明提供了一种基于卡尔曼滤波与支持向量机的测井曲线预测方法，包括以下步骤：

(1)选择足够数量的完备井为标准井，并对测井曲线进行标准化采样处理。标准井的曲线被划分为两组，一组对应目标井中缺失的，需要进行预测的曲线列表，设为X，另一组对应目标井中已知的，作为预测参数的曲线列表，设为Y。X与Y是2维数组，每列代表一条曲线，每行代表每条曲线的一个采样点，X与Y的行数必须相同，但列数可以不同，Y的列数要大于或等于X的列数。X与Y的每行代表一个特征向量，特征向量的值为在一定值域范围内的连续变量。

所述标准井为测井曲线较全的完备井，这些完备井的测井曲线为常规测井曲线，常规测井曲线是指声波曲线AC、自然伽马曲线GR、中子曲线CNL、电阻率曲线RT、井径曲线CAL等具体测井曲线；缺失井与标准井相比，缺少我们需要的曲线，设此曲线为X。

(2)建立状态空间模型。状态空间模型是动态时域模型，以隐含着的时间为自变量。使用标准井的两组曲线X、Y建立状态空间模型，如式(1)和式(2)所示：

X^t＝f^t(X^t-1)+Q^t (1)

Y^t＝o^t(X^t)+R^t (2)

其中，t代表时序指标，也就是测井曲线的采样位置指标，X^t代表X曲线列表的第t个采样点特征向量，一个特征向量就代表一种状态；f^t为状态转移函数，代表从第t-1个采样点的状态向地第t个采样点的状态转移的规律；o^t为观测函数，代表从状态X^t观测到Y^t的规律；Y^t代表Y曲线列表第t个采样点的特征向量；Q^t，R^t是服从于零均值正态分布的扰动向量，Q^t代表状态转移规律的随机性，R^t代表观测规律的随机性。建立状态空间模型的目标是利用标准井的两组曲线X和Y作为样本，通过某种方法推断或生成函数f^t、o^t、以及随机变量Q^t、R^t分布函数。

将方程(1)称为“状态”或“转移”方程(2)称为“信号”或“测量”方程。扰动向量，Q^t，R^t，f^t，o^t被称为系统矩阵。系统矩阵Q^t、R^t、f^t、o^t可以依赖于一个未知参数的集合。状态空间模型的一个主要的任务就是估计这些参数。为了便于表述，设θ＝(Q^t、R^t、f^t、o^t)为模型的参数。基于状态空间模型，任何测井曲线预测工作可表述为：在给定模型参数θ及批量观测数据y＝{y^t-n+1，y^t-n+2，...，y^t-1，y^t}的条件下，求模型的状态序列x＝{x^t-n+1，x^t-n+2，...，x^t-1，x^t}，其中n为批量观测数据y的长度。基于贝叶斯框架，可建立随机变量X^t的后验概率密度函数p(X^t|y，θ)＝c*p(y|X^t，θ)*p(X^t|θ)，在此后验概率密度函数p(X^t|y，θ)确定之后，建立随机序列X＝{X^t-n+1，X^t-n+2，...，X^t-1，X^t}的联合后验概率密度函数：p(X|Y，θ)＝p(X^t-n+1X^{t -n+2}，...，X^t-1，X^t|Y，θ)＝p(X^t-n+1|Y，θ)p(X^t-n+2|X^t-n+1，Y，θ)p(X^t-n+3|X_t-n+2，Y，θ)...p(X^t|X^t-1，Y，θ)。

(3)基于熵离散化测井曲线。

由于X^t表示的特征值是连续变量，代表了无限的状态空间，使用有限的样本推断或生成连续状态转移函数f^t特别困难，使用有限的离散化状态空间可以极大地降低推断或生成状态转移函数的难度，基于离散状态的状态空间模型可调整为如下式(3)至式(6)：

S^t＝S(X^t) (3)

S^t＝F^t(S^t-1) (4)

X^t＝F^t(S^-1(s^t-1))+Q^t (5)

Y^t＝o^t(X^t)+R^t (6)

其中，S^t为第t个连续型特征值(连续状态)对应的离散型特征值(离散状态)，S为基于熵的离散化函数，代表从连续状态到离散状态的转换，S^-1为逆离散化函数，代表从离散状态到连续状态的转换，F^t为离散状态转移函数，称为状态转移矩阵。任何离散化都会有信息的损失，基于熵的离散化可以最大地保留原始数据的有效信息，损失无效信息。基于熵的离散化是一种监督的、自顶向下的分裂技术。它在计算和确定分裂点时利用数据的分布信息。例如，为了离散化属性A，该方法选择A中最小熵的值作为分裂点，并递归地划分结果区间，得到分层离散化。

为了度量某一划分之后得到完全的分类还需要信息，引入期望信息需求的概念，期望信息需求由下式给出：

其中D₁和D₂分别对应于D中满足条件A≤split_point和A≥split_point的元组，|D|是D中的元组的个数，如此等等。集合中的熵函数根据下式来计算，假设集合D₁中的元素分别属于m个类，它们分别为C₁，C₂，...，C_m，D₁的熵是

其中，p_i是D₁中元组属于C_i的概率，由D₁中的C_i类元组数除以D₁中的元组总数|D₁|确定。这样在选择属性A的分裂点时，我们希望产生使得期望信息需求最小的属性值split_point作为分裂点，使得用A≤split_point和A＞split_point划分之后，对元组完全分类还需要的信息量最小。确定分裂点的过程递归地作用于所得到的每个划分，直到满足某个终止标准，如当所有候选点上的最小信息需求小于一个阈值，或者当区间的个数大于阈值max_interval时终止。本发明方法把上述基于熵离散化样本数据过程记为S()。设X^t＝f^t(X^t-1)+Q^t，其中X^t∈Rⁿ为模型状态，R为实数，n为状态X^t的维度，如我们需要同时预测声波曲线AC、自然伽马曲线GR、中子曲线CNL，可令n＝3，X^t∈R³，

其中

为声波曲线AC，

为自然伽马曲线GR，

为中子曲线CNL。

(4)建立状态迁移变换矩阵。

根据式(3)对连续状态X^t进行离散化后，得到状态空间{1、2、3...k}，S^t∈{1、2、3...k}为状态空间的一个状态，k为大于1的整数，根据式(3)的逆变换可得

其中S为基于熵离散化样本数据的函数，称为离散化变换，S^-1为离散化逆变换。根据式(3)、式(4)和式(5)可得式(7)

用{S^t}表示建立在状态空间{1、2、3...k}上的平稳马尔可夫链随机过程，F^t为随机过程{S^t}的状态迁移变换矩阵，是一个与时间无关的变换方阵，因此F^t可写为F，F是一个k阶方阵，是一个马尔可夫链状态迁移概率矩阵，F_ij矩阵F中的一个元素，代表从状态i向状态j转移的概率。建立状态迁移变换矩阵F的过程就是使用标准井的曲线列表X分组，首先用离散化变换S对样本进行离散化处理，即{s^t}＝S(X)，其中{s^t}为样本数据的离散状态序列，作为马尔可夫链状态迁移概率矩阵F的统计数据样本，统计出F的每个元素F_ii＝P(s^t＝j|s^t-1＝i)，其中P(s^t＝j|s^t-1＝i)表示在第t-1个状态为i的条件下第t个状态为j的概率。

(5)为状态残差随机变量建立服从高斯分布的概率密度函数。

根据式(7)可得

其中Q^t为状态残差随机变量，由于S^t-1为离散随机变量，S^t＝F^t(S^t-1)亦为离散随机变量，

作为状态S^t的对应连续随机变量X^t的期望值，把X^t看作k维高斯分布，即

其中ε²为X^t的协方差阵，是一个k阶方阵，因为

所以Q^t也为k维高斯分布，即得到式(8)：

Q^t～N(0，ε²) (8)

其中0为其均值，ε²为其协方差阵，与X^t的协方差阵相同。根据充分的样本数据{x^t}与步骤(1)和步骤(2)中建立的状态迁移变换矩阵F、离散化变换S以及离散化逆变换S^-1，计算得到数据残差序列{q^t}，其中q^t＝x^t-S^-1(s^t)，s^t＝S(x^t)，即q^t＝x^t-S^-1(S(x^t))，即{q^t}＝{x^t-S^-1(S(x^t))}。因为{q^t}为充分的样本数据，又因Q^t～N(0，ε²)，可根据{q^t}方便地统计推断出ε²。

(6)基于支持向量回归方法建立三次多项式非线性观测函数。

支持向量回归(support vector regression)简称SVR，是SVM(支持向量机support vector machine)对回归问题的一种运用。根据式(6)状态空间模型观测函数Y^t＝o^t(X^t)+R^t，其中o^t，代表状态观测函数，令

即得

{x^t}为样本中得状态数据序列，{y^t}为对应得观测数据序列，{x^t}与{y^t}组成状态观测对序列{(x^ty^t)}，使用SVR可以从{(x^t，y^t)}学习出o^t。学习得过程如下：1)对X^t进行升维处理，即把原有得n元1次的随机变量升维为n元3次随机变量；假如令X^t为2维随机变量，即

设升级后的随机变量为

即

即

为9维向量；设Y^t为4维的随机变量，显然

R^t亦同为4为向量，

本发明方法把其中的o^t看作线性变换，即o^t为4行9列的矩阵，记为W，即

称W为观测矩阵。因此样本数据{(x^t，y^t)}可表示为

即可根据样本数据可建立SVR的目标函数

其中

即目标函数为式(9)：

其中|W|为矩阵W行向量的模的和，c为正则化因子。通过SVR学习方法，求得W。

(7)为观测误差随机变量建立服从高斯分布的概率密度函数。

根据

可得

其中

R^t为观测残差随机变量，

作为状态

对应的观测随机变量Y^t的期望值，把Y^t看作m维高斯分布，即

其中μ²为Y^t的协方差阵，是一个m阶方阵，因为

所以R^t也为m维高斯分布，即得式(10)：

R^t～N(0，μ²) (10)

其中，0为其均值，ε²为其协方差阵，与Y^t的协方差阵相同。根据充分的样本数据{y^t}与步骤(2)～步骤(5)中建立的状态迁移变换矩阵F、离散化变换S、离散化逆变换S^-1以及观测矩阵W，计算得到观测数据残差序列{r^t}，其中

为从样本状态数据x^t升维得到的数据，即

即

因为{r^t}为充分的样本数据，又因R^t～N(0，μ²)，可根据{r^t}方便地统计推断出μ²。(8)基于状态空间模型利用卡尔曼滤波算法进行测井曲线预测。

状态空间模型为：X^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^t，本发明方法把其中的X^t、Q^t、Y^t、R^t看作服从于高斯分布的随机向量，把o^t()看作非线性变换矩阵，又根据步骤(6)，低维的非线性变换转变为高维的线性变换，因此模型为升维的线性高斯状态空间模型，通过前面的步骤与计算，f^t、Q^t、R^t，o^t为已知量，把缺失井的相应的源曲线作为Y^t的输入参数，调用卡尔曼滤波算法，可求解状态序列

即为目标预测曲线。

本发明的基于卡尔曼滤波与支持向量机的测井曲线预测方法特点是：本方法运用状态空间模型与支持向量机等机器学习技术建立预测模型，并把低维的非线性高斯状态空间模型通过多项式升维，变为线性高斯状态空间模型，最后调用卡尔曼滤波算法对缺失的测井曲线进行预测。

本发明与上述背景技术相比较可具有如下有益效果：

(1)常规做法是基于测井曲线完备的井使用回归方法建立预测模型，实现对缺失曲线的预测这些预测模型，但没有太多地关注测井曲线作为空间序列所蕴含的特征与结构，预测效果不理想。

(2)将本发明“一种基于卡尔曼滤波与支持向量机的测井曲线预测方法”应用于测井曲线预测中，包含所要预测的缺失曲线的井(称为缺失井)和与这些井所处同一油气田的至少一口测井曲线完备的井(称为完备井)，己知完备井的多条测井曲线，运用线性高斯模型与支持向量机等机器学习技术建立预测模型，并把低维的非线性高斯状态空间模型通过多项式升维，变为线性高斯状态空间模型，然后利用该预测模型和己知完备井的测井曲线调用卡尔曼滤波算法对缺失井进行测井曲线预测。

附图说明：

图1为本发明的完备井实测的声波曲线和预测的声波曲线示意图。

图2为本发明的缺失井的待预测的声波曲线示意图。

具体实施方式：

下面结合附图及实施例对本发明作进一步说明：

为使本发明的目的、技术方案和优点更加清楚，下面将以大庆探区松辽盆地古龙凹陷为例，结合附图对本发明实施方式作进一步地详细描述。

一种基于卡尔曼滤波与支持向量机的测井曲线预测方法，包括以下步骤：

(1)工区内选择足够数量的完备井为标准井，本例中选择213口，并对测井曲线进行标准化采样处理。标准井的曲线被划分为两组，一组对应目标井中缺失的，需要进行预测的曲线列表，设为X，另一组对应目标井中现存的，作为预测参数的曲线列表，设为Y。X与Y是2维数组，每列代表一条曲线，每行代表每条曲线的一个采样点，X与Y的行数必须相同，但列数可以不同，Y的列数要大于或等于X的列数。X与Y的每行代表一个特征向量，特征向量的值为在一定值域范围内的连续变量。

所述标准井为测井曲线较全的完备井，这些完备井的测井曲线为常规测井曲线，常规测井曲线是指声波曲线AC、自然伽马曲线GR、中子曲线CNL、电阻率曲线RT、井径曲线CAL等具体测井曲线；缺失井与标准井相比，缺少我们需要的曲线，设此曲线为X。

本例中，所用的输入Y＝[RHOB、CALI、GR、LLD]；输出为X＝ttt(DT)，ttt为声波测井曲线DT，用有声波曲线的井建立了基于卡尔曼滤波与支持向量机的预测模型。

(2)建立状态空间模型。状态空间模型是动态时域模型，以隐含着的时间为自变量。使用标准井的两组曲线X、Y建立状态空间模型的数学表达式(1)、(2)所示。

X^t＝f^t(X^t-1)+Q^t (1)

Y^t＝o^t(X^t)+R^t (2)

其中，t代表时序指标，也就是测井曲线的采样位置指标，X^t代表X曲线列表的第t个采样点特征向量，一个特征向量就代表一种状态；f^t为状态转移函数，代表从第t-1个采样点的状态向地第t个采样点的状态转移的规律；o^t为观测函数，代表从状态X^t观测到Y^t的规律；Y^t代表Y曲线列表第t个采样点的特征向量；Q^t，R^t是服从于零均值正态分布的扰动向量，Q^t代表状态转移规律的随机性，R^t代表观测规律的随机性。建立状态空间模型的目标是利用标准井的两组曲线X和Y作为样本，通过某种方法推断或生成函数f^t、o^t、以及随机变量Q^t、R^t分布函数。

将方程(1)称为“状态”或“转移”方程(2)称为“信号”或“测量”方程。扰动向量，Q^t、R^t、f^t、o^t被称为系统矩阵。系统矩阵Q^t、R^t、f^t、o^t可以依赖于一个未知参数的集合。状态空间模型的一个主要的任务就是估计这些参数。为了便于表述，设θ＝(Q^t、R^t、f^t、o^t)为模型的参数。基于状态空间模型，任何测井曲线预测工作可表述为：在给定模型参数θ及批量观测数据y＝{y^t-n+1，y^t-n+2，...，y^t-1，y^t}的条件下，求模型的状态序列x＝{x^t-n+1，x^t-n+2，...，x^t-1，x^t}，其中n为批量观测数据y的长度。基于贝叶斯框架，可建立随机变量X^t的后验概率密度函数p(X^t|y，θ)＝c*p(y|X^t，θ)*p(X^t|θ)，在此后验概率密度函数p(X^t|y，θ)确定之后，建立随机序列X＝{X^t-n+1，X^t-n+2，...，X^t-1，X^t}的联合后验概率密度函数：p(X|Y，θ)＝p(X^t-n+1，X^{t -n+2}，...，X^t-1，X^t|Y，θ)＝p(X^t-n+1|Y，θ)p(X^t-n+2|X^t-n+1，Y，θ)p(X^t-n+3|X_t-n+2，Y，θ)...p(X^t|X^t-1，Y，θ)。

(3)基于熵离散化测井曲线。

由于X^t表示的特征值是连续变量，代表了无限的状态空间，使用有限的样本推断或生成连续状态转移函数f^t特别困难，使用有限的离散化状态空间可以极大地降低推断或生成状态转移函数的难度，基于离散状态的状态空间模型可调整为如下式(3)、式(4)、式(5)、式(6)：

S^t＝S(X^t) (3)

S^t＝F^t(S^t-1) (4)

X^t＝F^t(S^-1(s^t-1))+Q^t (5)

Y^t＝o^t(X^t)+R^t (6)

其中，S^t为第t个连续型特征值(连续状态)对应的离散型特征值(离散状态)，S为基于熵的离散化函数，代表从连续状态到离散状态的转换，S^-1为逆离散化函数，代表从离散状态到连续状态的转换，F^t为离散状态转移函数，称为状态转移矩阵。任何离散化都会有信息的损失，基于熵的离散化可以最大地保留原始数据的有效信息，损失无效信息。基于熵的离散化是一种监督的、自顶向下的分裂技术。它在计算和确定分裂点时利用数据的分布信息。例如，为了离散化属性A，该方法选择A的具有最小熵的值作为分裂点，并递归地划分结果区间，得到分层离散化。

为了度量某一划分之后得到完全的分类还需要信息，引入期望信息需求的概念，期望信息需求由下式给出：

其中D₁和D₂分别对应于D中满足条件A≤split_point和A≥split_point的元组，|D|是D中的元组的个数，如此等等。集合中的熵函数根据下式来计算，假设集合D₁中的元素分别属于m个类，它们分别为C₁，C₂，...，C_m，D₁的熵是

其中，p_i是D₁中元组属于C_i的概率，由D₁中的C_i类元组数除以D₁中的元组总数|D₁|确定。这样在选择属性A的分裂点时，我们希望产生使得期望信息需求最小的属性值split_point作为分裂点，使得用A≤splitpoint和A＞split_point划分之后，对元组完全分类还需要的信息量最小。确定分裂点的过程递归地作用于所得到的每个划分，直到满足某个终止标准，如当所有候选点上的最小信息需求小于一个阈值，或者当区间的个数大于阈值max_interval时终止。本发明方法把上述基于熵离散化样本数据过程记为S()。设X^t＝f^t(X^t-1)+Q^t，其中X^t∈Rⁿ为模型状态，R为实数，n为状态X^t的维度，如我们需要同时预测声波曲线AC、自然伽马曲线GR、中子曲线CNL，可令n＝3，X^t∈R³，

其中

为声波曲线AC，

为自然伽马曲线GR，

为中子曲线CNL。

(4)建立状态迁移变换矩阵。

根据式(3)对连续状态X^t进行离散化后，得到状态空间{1、2、3...k}，S^t∈{1、2、3...k}为状态空间的一个状态，k为大于1的整数，根据式(3)的逆变换可得

其中S为基于熵离散化样本数据的函数，称为离散化变换，S^-1为离散化逆变换。根据式(3)、(4)、(5)可得式(7)：

用{S^t}表示建立在状态空间{1、2、3...k}上的平稳马尔可夫链随机过程，F^t为随机过程{S^t}的状态迁移变换矩阵，是一个与时间无关的变换方阵，因此F^t可写为F，F是一个k阶方阵，是一个马尔可夫链状态迁移概率矩阵，F_ij矩阵F中的一个元素，代表从状态i向状态j转移的概率。建立状态迁移变换矩阵F的过程就是使用标准井的曲线列表X分组，首先用离散化变换S对样本进行离散化处理，即{s^t}＝S(X)，其中{s^t}为样本数据的离散状态序列，作为马尔可夫链状态迁移概率矩阵F的统计数据样本，统计出F的每个元素F_ij＝P(s^t=j|s^t-1＝i)，其中P(s^t＝j|s^t-1＝i)表示在第t-1个状态为i的条件下第t个状态为j的概率。

(5)为状态残差随机变量建立服从高斯分布的概率密度函数。

根据式(7)可得

其中Q^t为状态残差随机变量，由于S^t-1为离散随机变量，S^t＝F^t(S^t-1)亦为离散随机变量，

作为状态S^t的对应连续随机变量X^t的期望值，把X^t看作k维高斯分布，即

其中ε²为X^t的协方差阵，是一个k阶方阵，因为

所以Q^t也为k维高斯分布，即得到式(8)：

Q^t～N(0，ε²) (8)

其中0为其均值，ε²为其协方差阵，与X^t的协方差阵相同。根据充分的样本数据{x^t}与步骤(1)和步骤(2)中建立的状态迁移变换矩阵F、离散化变换S以及离散化逆变换S^-1，计算得到数据残差序列{q^t}，其中q^t＝x^t-S^-1(s^t)，s^t＝S(x^t)，即q^t＝x^t-S^-1(S(x^t))，即{q^t}＝{x^t-S^-1(S(x^t))}。因为{q^t}为充分的样本数据，又因Q^t～N(0，ε²)，可根据{q^t}方便地统计推断出ε²。

(6)基于支持向量回归方法建立三次多项式非线性观测函数。

支持向量回归(support vector regression)简称SVR，是SVM(支持向量机support vector machine)对回归问题的一种运用。根据式(6)状态空间模型观测函数Y^t＝o^t(X^t)+R^t，其中o^t代表状态观测函数，令

即得

{x^t}为样本中得状态数据序列，{y^t}为对应得观测数据序列，{x^t}与{y^t}组成状态观测对序列{(x^t，y^t))}，使用SVR可以从{(x^t，y^t)}学习出o^t。学习得过程如下：1)对X^t进行升维处理，即把原有得n元1次的随机变量升维为n元3次随机变量；假如令X^t为2维随机变量，即

设升级后的随机变量为

即

即

为9维向量；设Y^t为4维的随机变量，显然

R^t亦同为4为向量，

本发明方法把其中的o^t看作线性变换，即o^t为4行9列的矩阵，记为W，即

称W为观测矩阵。因此样本数据{(x^t，y^t）}可表示为

即可根据样本数据可建立SVR的目标函数

其中

即目标函数为式(9)：

其中|W|为矩阵W行向量的模的和，c为正则化因子。通过SVR学习方法，求得W。

(7)为观测误差随机变量建立服从高斯分布的概率密度函数。

根据

可得

其中

R^t为观测残差随机变量，

作为状态

对应的观测随机变量Y^t的期望值，把Y^t看作m维高斯分布，即

其中μ²为Y^t的协方差阵，是一个m阶方阵，因为

所以R^t也为m维高斯分布，即得式(10)：

R^t～N(0，μ²) (10)

其中，0为其均值，ε²为其协方差阵，与Y^t的协方差阵相同。根据充分的样本数据{y^t}与步骤(2)～步骤(5)中建立的状态迁移变换矩阵F、离散化变换S、离散化逆变换S^-1以及观测矩阵W，计算得到观测数据残差序列{r^t}，其中

为从样本状态数据x^t升维得到的数据，即

即

因为{r^t}为充分的样本数据，又因R^t～N(0，μ²)，可根据{r^t}方便地统计推断出μ²。(8)基于状态空间模型利用卡尔曼滤波算法进行测井曲线预测。

状态空间模型为：X^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^t，本发明方法把其中的X^t、Q^t、Y^t、R^t看作服从于高斯分布的随机向量，把o^t()看作非线性变换矩阵，又根据步骤(6)，低维的非线性变换转变为高维的线性变换，因此模型为升维的线性高斯状态空间模型，通过前面的步骤与计算，f^t、Q^t、R^t，o^t为已知量，把缺失井的相应的源曲线作为Y^t的输入参数，调用卡尔曼滤波算法，可求解状态序列

即为目标预测曲线。

图1中波峰较低(即数字2表示的曲线)的曲线是实测的声波曲线，波峰较高的曲线(即数字1表示的曲线)是预测的声波曲线，图2为预测的目标井的声波曲线(即目标井缺失的声波曲线)。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于卡尔曼滤波与支持向量机的测井曲线预测方法，包括以下步骤：

(1)选择足够数量的完备井为标准井，并对测井曲线进行标准化采样处理，标准井的曲线被划分为两组，一组对应目标井中缺失的，需要进行预测的曲线列表，设为X，另一组对应目标井中已知的，作为预测参数的曲线列表，设为Y，X与Y是2维数组，每列代表一条曲线，每行代表每条曲线的一个采样点，X与Y的行数必须相同，但列数不同，Y的列数要大于或等于X的列数，X与Y的每行代表一个特征向量，特征向量的值为在一定值域范围内的连续变量；

(2)建立状态空间模型：状态空间模型是动态时域模型，以隐含着的时间为自变量，使用标准井的两组曲线X、Y建立状态空间模型的数学表达式(1)、(2)所示，

X^t＝f^t(X^t-1)+Q^t (1)

Y^t＝o^t(X^t)+R^t (2)

其中，t代表时序指标，也就是测井曲线的采样位置指标，X^t代表X曲线列表的第t个采样点特征向量，一个特征向量就代表一种状态；f^t为状态转移函数，代表从第t-1个采样点的状态向地第t个采样点的状态转移的规律；o^t为观测函数，代表从状态X^t观测到Y^t的规律；Y^t代表Y曲线列表第t个采样点的特征向量；Q^t，R^t是服从于零均值正态分布的扰动向量，Q^t代表状态转移规律的随机性，R^t代表观测规律的随机性；

(3)基于熵离散化测井曲线：

由于X^t表示的特征值是连续变量，代表了无限的状态空间，使用有限的样本推断或生成连续状态转移函数f^t特别困难，使用有限的离散化状态空间极大地降低推断或生成状态转移函数的难度，基于离散状态的状态空间模型调整为如下式(3)至式(6)：

S^t＝S(X^t) (3)

S^t＝F^t(S^t-1) (4)

X^t＝F^t(S^-1(s^t-1))+Q^t (5)

Y^t＝o^t(X^t)+R^t (6)

其中，S^t为第t个连续型特征值对应的离散型特征值，S为基于熵的离散化函数，代表从连续状态到离散状态的转换，S^-1为逆离散化函数，代表从离散状态到连续状态的转换，F^t为离散状态转移函数，称为状态转移矩阵；

(4)建立状态迁移变换矩阵：

根据式(3)对连续状态X^t进行离散化后，得到状态空间{1、2、3...k}，S^t∈{1、2、3…k}为状态空间的一个状态，k为大于1的整数，根据式(3)的逆变换得

其中S为基于熵离散化样本数据的函数，称为离散化变换,S^-1为离散化逆变换，根据式(3)、(4)、(5)得式(7)：

用{S^t}表示建立在状态空间{1、2、3...k}上的平稳马尔可夫链随机过程，F^t为随机过程{S^t}的状态迁移变换矩阵，是一个与时间无关的变换方阵，因此F^t写为F，F是一个k阶方阵，是一个马尔可夫链状态迁移概率矩阵，F_ij矩阵F中的一个元素，代表从状态i向状态j转移的概率，建立状态迁移变换矩阵F的过程就是使用标准井的曲线列表X分组，首先，用离散化变换S对样本进行离散化处理，即{s^t}＝S(X)，其中{s^t}为样本数据的离散状态序列，作为马尔可夫链状态迁移概率矩阵F的统计数据样本，统计出F的每个元素F_ij＝P(s^t＝j|s^t-1＝i)，其中P(s^t＝j|s^t-1＝i)表示在第t-1个状态为i的条件下第t个状态为j的概率；

(5)为状态残差随机变量建立服从高斯分布的概率密度函数：

根据式(7)得

其中Q^t为状态残差随机变量，由于S^t-1为离散随机变量，S^t＝F^t(S^t-1)亦为离散随机变量，

作为状态S^t的对应连续随机变量X^t的期望值，把X^t看作k维高斯分布，即

其中ε²为X^t的协方差阵，是一个k阶方阵，因为

所以Q^t也为k维高斯分布，即得到式(8)：

Q^t～N(0，ε²) (8)

其中0为其均值，ε²为其协方差阵，与X^t的协方差阵相同，根据充分的样本数据{x^t}与步骤(1)和步骤(2)中建立的状态迁移变换矩阵F、离散化变换S以及离散化逆变换S^-1，计算得到数据残差序列{q^t}，其中q^t＝x^t-S^-1(s^t)，s^t＝S(x^t)，即q^t＝x^t-S^-1(S(x^t))，即{q^t}＝{x^t-S^-1(S(x^t))}，因为{q^t}为充分的样本数据，又因Q^t～N(0，ε²)，根据{q^t}统计推断出ε²；

(6)基于支持向量回归方法建立三次多项式非线性观测函数：

根据式(6)状态空间模型观测函数Y^t＝o^t(X^t)+R^t，其中o^t,代表状态观测函数，令

即得

{x^t}为样本中得状态数据序列，{y^t}为对应得观测数据序列，{x^t}与{y^t}组成状态观测对序列{(x^t，y^t)}，使用SVR从{(x^t，y^t)}学习出o^t；

(7)为观测误差随机变量建立服从高斯分布的概率密度函数：

根据

可得

其中

R^t为观测残差随机变量，

作为状态

对应的观测随机变量Y^t的期望值，把Y^t看作m维高斯分布，即

其中μ²为Y^t的协方差阵，是一个m阶方阵，因为

所以R^t也为m维高斯分布，即得式(10)：

R^t～N(0，μ²) (10)

其中，0为其均值，ε²为其协方差阵，与Y^t的协方差阵相同；根据充分的样本数据{y^t}与步骤(2)～步骤(5)中建立的状态迁移变换矩阵F、离散化变换S、离散化逆变换S^-1以及观测矩阵W，计算得到观测数据残差序列{r^t}，其中

为从样本状态数据x^t升维得到的数据，即

即

因为{r^t}为充分的样本数据，又因R^t～N(0，μ²)，可根据{r^t}方便地统计推断出μ²；

(8)基于状态空间模型利用卡尔曼滤波算法进行测井曲线预测：

状态空间模型为：X^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^t，其中的X^t、Q^t、Y^t、R^t看作服从于高斯分布的随机向量，把o^t()看作非线性变换矩阵，又根据步骤(6)，低维的非线性变换转变为高维的线性变换，因此，是升维的线性高斯状态空间模型，通过前面的步骤与计算，f^t、Q^t、R^t，o^t为已知量，把缺失井的相应的源曲线作为Y^t的输入参数，调用卡尔曼滤波算法，可求解状态序列

得到目标预测曲线。

2.根据权利要求1所述的测井曲线预测方法，其特征在于，步骤(6)中的学习过程如下：1)对X^t进行升维处理，即把原有得n元1次的随机变量升维为n元3次随机变量；假如令X^t为2维随机变量，即

设升级后的随机变量为

即

即

为9维向量；设Y^t为4维的随机变量，

R^t亦同为4为向量，

其中的o^t看作线性变换，即o^t为4行9列的矩阵，记为W，即

称W为观测矩阵，因此样本数据{(x^t，y^t)}表示为

根据样本数据可建立SVR的目标函数

其中

即目标函数为式(9)：

其中|W|为矩阵W行向量的模的和，c为正则化因子，通过SVR学习方法，求得W。

Images