CN105388531A - 一种基于支持向量回归机和核fisher分析的岩性识别方法 - Google Patents

Abstract

本发明公开了一种基于支持向量回归机和核FISHER分析的岩性识别方法，解决了现有技术中岩性识别方法识别岩性不准确的缺陷。该方法包括以下步骤：（1）通过实测的含纵横波数据的测井数据，建立利用支持向量回归机预测未知横波模型的步骤；（2）根据预测的横波数据计算岩石模量的步骤；（3）岩石模量结合常规测井属性构建岩性识别的样本的步骤；（4）利用核Fisher判别分析计算投影方向的步骤；（5）将待识别的样本投影到步骤（4）计算得出的投影方向形成交会图，实现岩性的识别的步骤。本发明在利用支持向量回归机获得可靠横波数据的基础上，利用核Fisher判别分析识别岩性，从而提高了岩性识别的准确率。

Description

一种基于支持向量回归机和核FISHER分析的岩性识别方法

技术领域

本发明涉及油气勘探领域，具体的说，是涉及一种基于支持向量回归机和核FISHER分析的岩性识别方法。

背景技术

纵、横波数据是分析岩石物理的基础，在地震模拟、地震资料的AVO分析、流体替换模型的计算、储层预测的流体识别和岩性识别中都需要横波数据。横波预测的方法众多，既有统计经验公式的方法。然而，常规测井只能提供纵波数据，一般的工区由全波测井得到的横波数据较少，要利用测井、地震多属性进行储层预测，就必须开展横波预测工作。很多学者给出了不同的纵横波经验公式，但是，经验公式不具有普遍性，其原因在于，经验公式或方法对于岩性的依赖性很强，且依赖于只适用于当时研究的区域；严格来说，经验公式并不一定适用所有地区。事实上，最可靠和最常用的Vp-Vs关系是实验室或测井数据或者两者的经验拟合。横波速度预测的方法有基于统计理论的方法，如基于最小二乘原理的回归方法，也有非线性的方法。Brian(2004)提出用多重线性回归预测横波速度；Hampson(2001)等人提出利用神经网络建立横波速度数据与测井数据之间的关系进行横波的计算。

岩性识别是储层研究的基础和关键内容之一，是计算储层参数、评价储层的基础。从总体上说，岩性识别可以利用地震资料和测井资料进行识别。利用测井资料进行岩性识别的方法有两种：一是根据各种矿物的测井响应值结合交会图技术(如中子与声波交会，中子与密度交会，M-N值交会等)进行矿物的定性判别；二是在对录井资料和岩心薄片资料的研究分析基础上，运用线性和非线性的方法对岩性识别进行识别。随着地球物理的研究和对致密碎屑岩认识的深入，基于线性理论的方法也不能满足实际需要。因为地层的非均质性和测井响应的非线性特征，采用非线性的方法能够提高碎屑岩储层岩性识别的准确率。基于非线性方法的储层岩性识别包括两个方面：一是特征变量的选取，即选出对岩性识别敏感的测井属性，测井属性的个数应尽可能的少。在不同的研究中，测井属性的选取也不同，如选取自然电位、自然伽玛、声波时差、电阻率与补偿中子作为岩性识别的特征变量；选取声波时差、自然电位、自然伽玛、电阻率、微电位、井径、微梯度来识别岩性。二是岩性识别方法的选取，岩性识别的方法众多，有基于统计理论的方法，也有基于非线性理论的方法，基于传统统计学理论的方法有时因为自身的局限性难于取得理想的效果，但方法大多简单且易于实现；非线性方法如神经网络、支持向量回归机、PSO神经网络、分形等对具有非线性特征的问题一般能够取得好的效果，但方法实现要难于基于传统统计学的方法。识别方法的选取应根据研究对象和方法本身的特点来确定。

发明内容

本发明的目的在于克服上述缺陷，提供一种可提高岩性识别的准确率的基于支持向量回归机和核FISHER分析的岩性识别方法。

为了实现上述目的，本发明采用的技术方案如下：

一种基于支持向量回归机和核FISHER分析的岩性识别方法，包括以下步骤：

(1)通过实测的含纵横波数据的测井数据，建立利用支持向量回归机预测未知横波模型的步骤；

(2)根据预测的横波数据计算岩石模量的步骤；

(3)岩石模量结合常规测井属性构建岩性识别的样本的步骤；

(4)利用核Fisher判别分析计算投影方向的步骤；

(5)将待识别的样本投影到步骤(4)计算得出的投影方向形成交会图，实现岩性的识别的步骤。

进一步的，所述步骤(1)的具体方式如下：

(11)从实测的含纵横波数据的测井数据中提取训练样本构成训练样本集D＝{(x_i,y_i)|i＝1,2,…,n}，式中，x_i为横波，y_i为纵波；

(12)利用训练样本和支持向量回归机建立横波的预测模型：

$y=\underset{i=1}{\overset{n}{\Σ}}{a}_{i}K(x_i,x)+b$

式中K(x_i,x)))为核函数，b为偏置量，α_i为权向量；

(13)利用支持向量回归机预测的横波与实测横波进行相关分析，若相关较强，则步骤(12)中的预测模型可用于步骤(2)预测横波。步骤(13)是对预测模型可靠性的检验。

进一步的，所述步骤(2)的具体方式如下：

(21)利用所述步骤(2)建立的横波预测模型和实测纵波数据V_P，得到预测横波数据V_S；

(22)利用公式：Vs＝Δts*1000/(12*2.54)得到横波测井时差Δts；

(23)根据横波测井时差Δts和实测的纵波测井时差Δtp可计算得出各模量，分别如下：

扬氏模量E： $E=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2}\×\[\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}\];$

体积模量K： $K=\mathrm{\ρ}\×\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δt}}^2}{3{\mathrm{\Δts}}^2\×{\mathrm{\Δtp}}^2};$

剪切模量μ： $\mathrm{\μ}=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2};$

泊松比v： $v=\frac{0.5{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2};$

拉梅常数λ： $\mathrm{\λ}=\mathrm{\ρ}\left(\frac{{\mathrm{tp}}^2-2{\mathrm{ts}}^2}{{\mathrm{ts}}^2\×{\mathrm{tp}}^2}\right);$

上式中，Δts为横波测井时差，Δtp为纵波测井时差，ρ为密度。

进一步的，所述步骤(3)的具体方式如下：

(31)根据V_S、λ、μ计算纵波横波速度比V_p/v_s、λ·ρ、μ·ρ；

(32)以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL作为特征变量，选取标准层获得训练集和测试集；

(33)选取扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

(34)以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL、扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

其中，步骤(32)、步骤(33)、步骤(34)为并列步骤。

进一步的，所述步骤(4)的具体方式如下：

(41)建立核Fisher判别分析模型：K_bα＝λK_wα；

式中，K_w为核类内散度矩阵，K_b为核类间散度矩阵，α为投影方向

(411)选取特征变量获取数据构成训练样本，并对其进行标准化处理，可供选取的特征变量有常规测井属性、弹性参数；

(412)由式ξ_x＝(K(x,x₁),K(x,x₂),…,K(x,x_N))计算ξ_x；

(413)利用下列公式：

${\mathrm{\μ}}_0={\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K\right(x_1,x_i),\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_2,x_i),...,\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_N,x_i\left)\right)}^T$

${\mathrm{\μ}}_i=\left(\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K\right(x_1,x_j^i)),\mathrm{...},\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K(x_N,x_j^i)$

计算每一类样本的均值μ_i和总体样本的均值μ₀；

(414)利用下列公式：

$K_b=\underset{i=1}{\overset{K}{\Σ}}\frac{N_i}{N}({\mathrm{\μ}}_i-{\mathrm{\μ}}_0){({\mathrm{\μ}}_i-{\mathrm{\μ}}_0)}^T$

$K_w=\frac{1}{N}\underset{i=1}{\overset{K}{\Σ}}\underset{j=1}{\overset{N_i}{\Σ}}({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i){({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i)}^T$

计算核类内散度矩阵K_w和核类间散度矩阵K_b，得到K_bα＝λK_wα；

(42)正则参数优化：采用正则化方法将K_bα＝λK_wα转换为K_bα＝λ(K_w+μI)α；

(43)用核Fisher判别分析方法进行岩性识别：

(431)利用正则化参数优化方法计算最优的分类参数值μ_opt；

(432)计算K_bα＝λ(K_w+μ_optI)α的特征值和特征向量，取最大及次大的特征值所对应的特征向量作单位化得到投影方向α_first和α_second；

(433)将待判岩性样本x_i(i＝1,2,…,m)投影到α_first和α_second得到投影值和

(434)根据投影值和作交会图，根据交会图对储层岩性进行判别。

进一步的，所述步骤(431)的具体方式：

(4311)初始化μ₀＝20，迭代次数n＝1；

(4312)计算μ＝μ₀-0.01n；

(4313)计算当f(μ)＝|K_w+μI|，以f(μ)纵坐标，以μ横坐标，作f(μ)随μ的变化图；

(4314)当f(μ)＝|K_w+μI|趋于0且存在明显的拐点时，μ的值为最优的分类参数值；否则迭代次数增加一次，并转步骤(4312)。

与现有技术相比，本发明具有以下有益效果：

本发明利用支持向量回归机和核Fisher判别分析来识别岩性；首先实测的测井数据，利用支持向量回归机预测横波数据，并根据预测的横波数据计算各种岩石模量；结合常规测井属性，构建岩性识别的样本集；利用核Fisher判别分析计算投影方向，将待识别的样本投影到上述方向形成交会图实现岩性的识别；本发明综合了支持向量回归机和核Fisher判别分析两种方法的优点，取长补短，在利用支持向量回归机获得可靠横波数据的基础上，利用核Fisher判别分析识别岩性，从而提高了岩性识别的准确率。

附图说明

图1为本发明的原理示意图。

图2为本发明中支持向量机预测横波的原理示意图。

图3为本发明中正规化参数优化的原理示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。

实施例

如图1所示，本实施例提供了一种基于支持向量回归机和核Fisher分析的岩性识别方法，该方法综合了支持向量回归机和核Fisher判别分析来识别岩性，在利用支持向量回归机获得可靠横波数据的基础上，利用核Fisher判别分析识别岩性。该方法包括以下步骤：(1)通过实测的含纵横波数据的测井数据，建立利用支持向量回归机预测未知横波模型的步骤；(2)根据预测的横波数据计算岩石模量的步骤；(3)岩石模量结合常规测井属性构建岩性识别的样本的步骤；(4)利用核Fisher判别分析计算投影方向的步骤；(5)将待识别的样本投影到步骤(4)计算得出的投影方向形成交会图，实现岩性的识别的步骤。

一、

支持向量回归机预测横波

1.1支持向量回归机模型

最小二乘支持向量机是Suykens等(1998)提出的用于解决分类问题和回归问题的新型支持向量机方法。最小二乘支持向量回归机(LeastSquaresSupportVectorRegression,LS-SVR)采用最小二乘线性系统作为损失函数，代替传统的支持向量机采用的二次规划方法。最小二乘支持向量机的训练只要求解一个线性方程组，因而在模式识别和函数估计时取得了较好的效果。

设有n个样本的训练集S＝{(x_i,y_i)|i＝1,2,…,n},x_i∈R^d,y_i∈R，样本集中x_i是输入数据，y_i是输出数据。最小二乘支持向量机就是在权w空间中求解下面的优化问题(式3-5)：

式(3-5)中是核空间映射函数，e_i∈R为误差变量，b是偏差量，γ为可调的常数。通过采用核函数，从输入空间中抽取特征，从而将输入空间中的样本映射为高维特征空间中的一个向量，解决输入空间中线性不可分的问题。式(3-5)的对偶问题的Lagrange函数为(式3-6)

对函数L(w,b,e,a)进行优化，即对w,b,e,a求偏导并令其等于零，得到

消除方程给(式3-7)中的变量w,e可得到等价的方程，其矩阵形式如下：

式(3-8)中，l_v＝[1,1,…,1]，a＝[a₁,a₂,…,a_n]，y＝[y₁,y₂,…,y_n]，k,l＝1,2,…,n。矩阵Ω+γ^-1I为一个对称半正定矩阵，解线性方程组(3-8)得到解(式3-9)：

根据Mercer条件，存在映射函数和核函数K(x_i,x)

最小二乘支持向量机的函数估计为

1.2利用支持向量回归机预测横波的步骤

如图2所示，首先，从实测的含纵横波数据的测井数据中提取训练样本构成训练样本集D＝{(x_i,y_i)|i＝1,2,…,n}，式中，x_i为横波，y_i为纵波；

其次，利用训练样本和支持向量回归机建立横波的预测模型：

式中，K(x_i,x)))为核函数，b为偏置量，α_i为权向量；

最后，利用支持向量回归机预测的横波与实测横波进行相关分析，若相关较强，则步骤(12)中的预测模型可用于步骤(2)预测横波。

二、根据预测的横波数据计算岩石模量的步骤

岩石的弹性参数受深度、地层的各向异性、地壳运动和岩性等地质因数的影响。岩石的结构和地质发展史在决定岩石的弹性性质上，比矿物组分上更为重要，通常结晶的岩石比破碎的岩石呈现出更大的弹性模量值，较大的地壳运动和出现胶结物质脱水时，也会增加弹性模量值。从预测岩石力学性质的实际观点出发，根据测井所做的动力学测量，是最有效的方法，连续的测井资料突出地表现了有关地层的这些性质，这样可以用来对力学性质不同的目的层的预测结果进行比较，因此利用测井资料计算各种岩石力学参数。

首先，利用建立的横波预测模型和实测纵波数据V_P，得到预测横波数据V_S；

其次，利用公式：Vs＝Δts*1000/(12*2.54)得到横波测井时差Δts；

最后，根据横波测井时差Δts和实测的纵波测井时差Δtp可计算得出各模量，分别如下：

扬氏模量E(Gpa)：它的物理意义是单位截面积的杆件伸长一倍所需要的应用大小，E值越大说明固体介质对拉伸力的阻力越大，弹性越好，其计算公式如(式4-1)。

$E=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2}\×\[\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}\]---(4-1)$

体积模量K(Gpa)：其代表物体的抗压缩性质，体现了岩石的耐压程度(式4-2)。

$K=\mathrm{\ρ}\×\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δt}}^2}{3{\mathrm{\Δts}}^2\×{\mathrm{\Δtp}}^2}---(4-2)$

剪切模量μ(Gpa)：物理意义是切应力与切应变的比，它是阻止剪切应变的一个度量(式4-3)。

$\mathrm{\μ}=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2}---(4-3)$

泊松比v：物理意义为杆件受力的作用相对缩短量(或伸长量)与其截面尺寸相对增大量(缩小量)之比(式4-4)。岩石中孔隙岩性的性质对泊松比有明显的影响，因此它是反映岩性和含油气性的重要参数。

$v=\frac{0.5{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}---(4-4)$

拉梅常数λ(Gpa)

$\mathrm{\λ}=\mathrm{\ρ}\left(\frac{{\mathrm{tp}}^2-2{\mathrm{ts}}^2}{{\mathrm{ts}}^2\×{\mathrm{tp}}^2}\right)---(4-5)$

式中，Δts为横波测井时差，Δtp为纵波测井时差，ρ为密度。

三、岩石模量结合常规测井属性构建岩性识别的样本

首先，根据V_S、λ、μ计算纵波横波速度比V_p/v_s、λ·ρ、μ·ρ；

其次，以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL作为特征变量，选取标准层获得训练集和测试集；

再其次，选取扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

最后，以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL、扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

其中，步骤(32)、步骤(33)、步骤(34)为并列步骤。

四、核Fisher判别分析预测岩性

首先，建立核Fisher判别分析模型：K_bα＝λK_wα；式中，K_w为核类内散度矩阵，K_b为核类间散度矩阵，α为投影方向；

具体的建立方式如下：

假设一个包含N个d维的K类岩性样本集S＝{x₁,x₂,…,x_N}，其中N为总岩性样本数，d为观测变量数，Ni个样本属于第i类，N＝N₁+N₂+…+N_k。通过映射把样本投影到一个高维的特征空间F中，从而实现样本的最优分类。在特征空间F中，表示第j类的第i个样本的投影值，设总体的平均值向量为ω_i类的总体样本平均值向量为m^φ(式5-2)

$m_j^{\mathrm{\φ}}=\frac{1}{N_j}\underset{i=1}{\overset{N_j}{\Σ}}\mathrm{\φ}\left(x_i^j\right)---(5-1)$

$m^{\mathrm{\φ}}=\underset{j=1}{\overset{K}{\Σ}}\underset{i=1}{\overset{N_j}{\Σ}}\mathrm{\φ}\left(x_i^j\right)---(5-2)$

在特征空间F中定义总体散度矩阵、类间散度矩阵、类内散度矩阵。总体散度矩阵S_t为(式5-3)

$S_t^{\mathrm{\φ}}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}(\mathrm{\φ}\left(x_i\right)-m^{\mathrm{\φ}}){(\mathrm{\φ}\left(x_i\right)-m^{\mathrm{\φ}})}^T---\left(5.3\right)$

类内离散度矩阵S_w为(式5-4)

$S_w^{\mathrm{\φ}}=\frac{1}{N}\underset{j=1}{\overset{K}{\Σ}}\underset{i=1}{\overset{N_j}{\Σ}}(\mathrm{\φ}\left(x_i^j\right)-m_j^{\mathrm{\φ}}){(\mathrm{\φ}\left(x_i^j\right)-m_j^{\mathrm{\φ}})}^T---\left(5-4\right)$

类间离散度矩阵S_b为(式5-5)

$S_b^{\mathrm{\φ}}=\underset{j=1}{\overset{K}{\Σ}}\frac{N_j}{N}(m_j^{\mathrm{\φ}}-m^{\mathrm{\φ}}){(m_j^{\mathrm{\φ}}-m^{\mathrm{\φ}})}^T---\left(5-5\right)$

S_t、S_w和S_b均为非负定矩阵，三者之间有如下的关系：

$S_b^{\mathrm{\φ}}+S_w^{\mathrm{\φ}}=S_t^i$

通过寻找最佳的投影向量使得同类岩性的样本有最大的类间散度，而不同类岩性的样本有最小类内散度，即通过最优化核Fisher准则函数(式5-6)找到一个最优的投影矩阵V_opt

$J\left(v\right)=\frac{v^T{S}_b^{\mathrm{\φ}}v}{v^T{S}_w^{\mathrm{\φ}}v}---(5-6)$

把输入空间中的样本数据通过非线性映射φ投影到高维的特征空间F中。但由于特征空间H的维数很高乃至于是无穷维的，因此不可能直接计算出最优Fisher判别向量，解决这一问题的方法就是核技巧，无须对数据进行非线性映射，使它只包含投影后数据的内积运算(式5-7)

K(x_i,x_j)＝(φ(x_i),φ(x_j))(5-7)

上述的高维空间可以认为是一个再生核希尔伯特空间，根据再生核的理论(etal.,1999)，任何一个最优化准则函数的解向量v一定位于由特征空间F中的数据φ(x₁),φ(x),…,φ(x_N)所张成的空间内(式5-8)

$v=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\α}}_i\mathrm{\φ}\left(x_i\right)---(5-8)$

式(5-8)中的α称为对应于特征空间F中的最优特征向量v的最优判别向量。在特征空间F中把任一测试样本x投影到v上有(5-9)

$\begin{array}{c}{v}^T\mathrm{\φ}\left(x\right)=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\α}}_i\left(\mathrm{\φ}\left(x\right),\mathrm{\φ}\left(x_i\right)\right)\\ ={\mathrm{\α}}^T\left(\left(\mathrm{\φ}\left(x\right)\mathrm{\φ}\left(x_1\right)\right),\left(\mathrm{\φ}\left(x\right),\mathrm{\φ}\left(x_2\right)\right),\mathrm{...},\left(\mathrm{\φ}\left(x\right),\mathrm{\φ}\left(x_N\right)\right)\right)\\ ={\mathrm{\α}}^T\left(K\left(x,x_1\right),K\left(x,x_2\right),\mathrm{...},K\left(x,x_N\right)\right)\end{array}---\left(5-9\right)$

记ξ_x＝(K(x,x₁),K(x,x₂),…,K(x,x_N))，在特征空间F中把类内均值(式5-10)和总体均值m^φ投影到v上有(式5-11)

${\mathrm{vm}}_i^{\mathrm{\φ}}={\mathrm{\α}}^T{\mathrm{\φ}}^T\frac{1}{N}\underset{j=1}{\overset{N_j}{\Σ}}\mathrm{\φ}\left(x_j^i\right)---(5-10)$

${\mathrm{vm}}^{\mathrm{\φ}}={\mathrm{\α}}^T{\mathrm{\φ}}^T\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}\mathrm{\φ}\left(x_j\right)---(5-11)$

可以得到投影后的类内均值μ_i(式5-12)和总体均值μ₀(式5-13)，再定义核类内散度矩阵K_w(式5-14)和核类间散度矩阵K_b(式5-15)后，在高维特征空间F中的判别准则函数(式5-6)等价于(式5-16)

${\mathrm{\μ}}_0={\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K\right(x_1,x_i),\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_2,x_i),...,\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_N,x_i\left)\right)}^T---\left(5-12\right)$

${\mathrm{\μ}}_i=\left(\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K\left(x_1,x_j^i\right)\right),\mathrm{...},\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K\left(x_N,x_j^i\right)---\left(5-13\right)$

$K_b=\underset{i=1}{\overset{K}{\Σ}}\frac{N_i}{N}({\mathrm{\μ}}_i-{\mathrm{\μ}}_0){({\mathrm{\μ}}_i-{\mathrm{\μ}}_0)}^T---\left(5-14\right)$

$K_w=\frac{1}{N}\underset{i=1}{\overset{K}{\Σ}}\underset{j=1}{\overset{N_i}{\Σ}}({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i){({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i)}^T---\left(5-15\right)$

$J\left(v\right)=\frac{v^T{S}_b^{\mathrm{\φ}}v}{v^T{S}_w^{\mathrm{\φ}}v}=\frac{{\mathrm{\α}}^T{K}_b\mathrm{\α}}{{\mathrm{\α}}^T{K}_w\mathrm{\α}}=J\left(\mathrm{\α}\right)---(5-16)$

由广义Rayleigh商的极值性质，使准则函数(5-6)极大化的最优解向量等价于求解广义特征方程(式5-17)

K_bα＝λK_wα(5-17)

正则参数优化

在K_w非奇异的情况下，最大化式(5-13)的最优解向量α为广义特征方程(5-17)的前m个最大特征值所对应的特征向量，此时式(5-17)可以等价为式(5-18)。

$K_w^{-1}{K}_b\mathrm{\α}=\mathrm{\λ}\mathrm{\α}---(5-18)$

在使用本发明解决实际问题时，此条件一般不成立：从n个训练样本估计n维的方差结构，所以K_w是奇异的。从数值稳定性或分类器的推广性角度考虑，我们可以采用正则化的方法，即在矩阵的对角元素加一个很小的正数(式5-19)

(K_w)_μ＝K_w+μI(5-19)

使得(K_w)_μ为非奇异矩阵。此时式(5-17)可表示为式(5-20)

K_bα＝λ(K_w+μI)α(5-20)

在一般的应用中，常根据经验或是实验的结果来确定参数μ。本实施例采用行列式分析的方法来解决参数μ的选取：当f(μ)＝|K_w+μI|趋于0且存在明显的拐点时，μ的值为最优的分类参数值，如图3所示：

首先，初始化μ₀＝20，迭代次数n＝1；

其次，计算μ＝μ₀-0.01n；

再其次，计算当f(μ)＝|K_w+μI|，以f(μ)纵坐标，以μ横坐标，作f(μ)随μ的变化图；

最后，当f(μ)＝|K_w+μI|趋于0且存在明显的拐点时，μ的值为最优的分类参数值；否则迭代次数增加一次，即n＝n+1，并转步骤：计算μ＝μ₀-0.01n。

四、用核Fisher判别分析方法进行岩性识别

首先，利用正则化参数优化方法计算最优的分类参数值μ_opt；

其次，计算K_bα＝λ(K_w+μ_optI)α的特征值和特征向量，取最大及次大的特征值所对应的特征向量作单位化得到投影方向α_first和α_second；

再其次，将待判岩性样本x_i(i＝1,2,…,m)投影到α_first和α_second得到投影值和

最好，根据投影值和作交会图，根据交会图对储层岩性进行判别。

预测实例

XC气田位于四川盆地川西坳陷中段，川西及邻区上三叠统-侏罗系地层是充填川西前陆盆地巨厚实体的一部分，总厚达6500m。XC气田位于坳陷中段XQ-FG北东东向的大型隆起带西段，须家河组须二构造为XQ与XC地区的多个局部高点构成的复式背斜(吴朝容，2005；刘爱疆，2010)。本次的基于本发明的弹性参数岩性识别研究的目的层为须二段。须二段地层主要分布在川西成都、德阳、绵阳、江油一带，厚度一般为400-700m，在龙门山前安县附近最厚可达1000m。

(一)岩石类型

须家河组为一套砂岩、粉砂岩、泥岩、页岩以及煤系地层组成的互层，其岩石类型较为复杂。通过对须二段薄片鉴定成果的统计和分析表明，岩石类型有岩屑石英砂岩、岩屑砂岩、长石岩屑石英砂岩、长石岩屑砂岩、石英砂岩、粉砂岩、泥岩等。须二段与须家河其余段相比较，长石含量较高的砂岩较为发育是其另外一个显著特征。须二段中高长石含量砂岩主要分布在中亚段中，上亚段、下亚段砂体中长石含量较低。须二段部分层段的高含量长石加上较好的岩石结构和可溶的基性岩岩屑，因此，形成的溶蚀孔隙对储层储集性有改善作用。由于须家河组各段中均发育较多的烃源岩，这些烃源岩在生排烃的多个时期也生成了较为充分的酸性岩性，因此只要砂岩中有可溶矿物，就可能被溶蚀形成发育的溶蚀孔隙。砂岩粒度组成从粗砂岩、中砂岩、细砂岩到粉砂岩均发育，主要是中砂岩、细砂岩。岩石普遍致密。

(二)岩石成分特征

据岩石薄片鉴定资料统计，碎屑物以石英为主，其含量为60％-90％，平均含量为73.9％；其次岩屑含量为10％-40％，平均含量21.1％；长石大多分布于中亚段，其含量较低为2％-5％，平均值为8。2％。岩屑成份包括凝灰岩、石英岩、花岗岩、板岩、泥岩、灰岩岩屑等。储层砂岩具有很低的杂基含量，平均值不到2％。胶结物以硅质为主，含量为1.37％；方解石含量普遍较少，分布不均；含少量白云石，泥质含量平均为1.05％。胶结类型以压结式为主，孔隙式－压结式次之，分选好，磨园度以次棱角状为主。

在应用本发明识别碎屑岩岩性时，采用三种实验方案，即第一种选取声波时差x₁、自然伽玛x₂、密度x₃、补偿中子x₄作为特征变量，获得样本集；第二种选取纵波速度V_p(x₁)、横波速度V_s(x₂)、弹性参数泊松比(x₃)、以及λ·ρ(x₄)、μ·ρ(x₅)作为特征变量获得样本集；第三种选取声波时差(x₁)、自然伽玛(x₂)、密度(x₃)、补偿中子(x₄)纵波速度V_p(x₅)、横波速度V_s(x₆)、弹性参数泊松比(x₇)、以及λ·ρ(x₈)、μ·ρ(x₉)作为特征变量，获得样本集。将三种样本集分为两部分：训练集与测试集。通过对四种核函数的对比实验，本发明采用高斯核函数时的效果最佳。

选取三种特征变量即测井属性、弹性参数、以及二者的组合，从采用高斯核函数的本发明预测结果来看，能够很好的区分泥岩；利用测井属性与或弹性参数识别岩性时，砂岩与粉砂岩的部分样本是重叠或混杂在一起的，但采用弹性参数识别时，砂岩与粉砂岩的样本聚合程度比采用测井属性时更高，识别的效果更好一些；利用测井属性与弹性参数的组合识别岩性时，三类岩性的样本在图中能够分开，且聚合程度较好，反映出利用综合的属性识别岩性，效果是最佳的，这也充分说明了本发明在致密碎屑岩岩性识别中的有效性，展示了该方法的应用前景。

上述实施例仅为本发明的部分实施例而已，并不用以限制本发明，而且本发明中各项参数可以根据实际使用效果随意选取组合。基于上述结构设计的前提下，为解决同样的技术问题，即使在本发明上做出的一些无实质性的改动或润色，所采用的技术方案的实质仍然与本发明一致的，也应当在本发明的保护范围内。

Claims (6)

1.一种基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，包括以下步骤：

(1)通过实测的含纵横波数据的测井数据，建立利用支持向量回归机预测未知横波模型的步骤；

(2)根据预测的横波数据计算岩石模量的步骤；

(3)岩石模量结合常规测井属性构建岩性识别的样本的步骤；

(4)利用核Fisher判别分析计算投影方向的步骤；

(5)将待识别的样本投影到步骤(4)计算得出的投影方向形成交会图，实现岩性的识别的步骤。

2.根据权利要求1所述的基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，所述步骤(1)的具体方式如下：

(11)从实测的含纵横波数据的测井数据中提取训练样本构成训练样本集D＝{(x_i,y_i)|i＝1,2,…,n}，式中，x_i为横波，y_i为纵波；

(12)利用训练样本和支持向量回归机建立横波的预测模型：

$y=\underset{i=1}{\overset{n}{\Σ}}{a}_{i}K(x_i,x)+b$

式中K(x_i,x)))为核函数，b为偏置量，α_i为权向量；

(13)利用支持向量回归机预测的横波与实测横波进行相关分析，若相关较强，则步骤(12)中的预测模型可用于步骤(2)预测横波。

3.根据权利要求2所述的基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，所述步骤(2)的具体方式如下：

(21)利用所述步骤(2)建立的横波预测模型和实测纵波数据V_P，得到预测横波数据V_S；

(22)利用公式：Vs＝Δts*1000/(12*2.54)得到横波测井时差Δts；

(23)根据横波测井时差Δts和实测的纵波测井时差Δtp可计算得出各模量，分别如下：

扬氏模量E： $E=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2}\×\[\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}\];$

体积模量K： $K=\mathrm{\ρ}\×\frac{3{\mathrm{\Δts}}^2-4{\mathrm{\Δt}}^2}{3{\mathrm{\Δts}}^2\×{\mathrm{\Δtp}}^2};$

剪切模量μ： $\mathrm{\μ}=\frac{\mathrm{\ρ}}{{\mathrm{\Δts}}^2};$

泊松比v： $v=\frac{0.5{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2}{{\mathrm{\Δts}}^2-{\mathrm{\Δtp}}^2};$

拉梅常数λ： $\mathrm{\λ}=\mathrm{\ρ}\left(\frac{{\mathrm{tp}}^2-2{\mathrm{ts}}^2}{{\mathrm{ts}}^2\×{\mathrm{tp}}^2}\right);$

上式中，Δts为横波测井时差，Δtp为纵波测井时差，ρ为密度。

4.根据权利要求3所述的基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，所述步骤(3)的具体方式如下：

(31)根据V_S、λ、μ计算纵波横波速度比V_p/v_s、λ·ρ、μ·ρ；

(32)以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL作为特征变量，选取标准层获得训练集和测试集；

(33)选取扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

(34)以声波时差AC、自然伽玛GR、密度DEN、补偿中子CNL、扬氏模量E、体积模量K、纵波横波速度比V_p/v_s、泊松比、以及λ·ρ、μ·ρ作为特征变量，选取与步骤(32)相同的层位获得训练集和测试集；

其中，步骤(32)、步骤(33)、步骤(34)为并列步骤。

5.根据权利要求4所述的基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，所述步骤(4)的具体方式如下：

(41)建立核Fisher判别分析模型：K_bα＝λK_wα；

式中，K_w为核类内散度矩阵，K_b为核类间散度矩阵，α为投影方向

(411)选取特征变量获取数据构成训练样本，并对其进行标准化处理，可供选取的特征变量有常规测井属性、弹性参数；

(412)由式ξ_x＝(K(x,x₁),K(x,x₂),…,K(x,x_N))计算ξ_x；

(413)利用下列公式：

${\mathrm{\μ}}_0={\left(\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K\right(x_1,x_i),\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_2,x_i),...,\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}K(x_N,x_i\left)\right)}^T$

${\mathrm{\μ}}_i=\left(\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K\right(x_1,x_j^i\left)\right),...,\frac{1}{N_i}\underset{j=1}{\overset{N_i}{\Σ}}K(x_N,x_j^i)$

计算每一类样本的均值μ_i和总体样本的均值μ₀；

(414)利用下列公式：

$K_b=\underset{i=1}{\overset{K}{\Σ}}\frac{N_i}{N}({\mathrm{\μ}}_i-{\mathrm{\μ}}_0){({\mathrm{\μ}}_i-{\mathrm{\μ}}_0)}^T$

$K_w=\frac{1}{N}\underset{i=1}{\overset{K}{\Σ}}\underset{j=1}{\overset{N_i}{\Σ}}({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i){({\mathrm{\ξ}}_{x_j}-{\mathrm{\μ}}_i)}^T$

计算核类内散度矩阵K_w和核类间散度矩阵K_b，得到K_bα＝λK_wα；

(42)正则参数优化：采用正则化方法将K_bα＝λK_wα转换为K_bα＝λ(K_w+μI)α；

(43)用核Fisher判别分析方法进行岩性识别：

(431)利用正则化参数优化方法计算最优的分类参数值μ_opt；

(432)计算K_bα＝λ(K_w+μ_optI)α的特征值和特征向量，取最大及次大的特征值所对应的特征向量作单位化得到投影方向α_first和α_second；

(433)将待判岩性样本x_i(i＝1,2,…,m)投影到α_first和α_second得到投影值和

(434)根据投影值和作交会图，根据交会图对储层岩性进行判别。

6.根据权利要求5所述的基于支持向量回归机和核FISHER分析的岩性识别方法，其特征在于，所述步骤(431)的具体方式：

(4311)初始化μ₀＝20，迭代次数n＝1；

(4312)计算μ＝μ₀-0.01n；

(4313)计算当f(μ)＝|K_w+μI|，以f(μ)纵坐标，以μ横坐标，作f(μ)随μ的变化图；

(4314)当f(μ)＝|K_w+μI|趋于0且存在明显的拐点时，μ的值为最优的分类参数值；否则迭代次数增加一次，并转步骤(4312)。