CN103375166A - 一种确定地层中多种矿物组分含量的方法 - Google Patents

Abstract

一种确定地层中的多种矿物组分的相对含量的方法，包括：提供关于地层中的多种矿物组分的信息，所述信息包括多种矿物组分的名称和数目n，其中n为大于等于2的整数；利用至少一个测井仪器对所述地层的m个参数进行测量，其中m为大于等于n的整数；获得地层对所述至少一个测井仪器的m个响应值；基于所述m个响应值和n种矿物组分的信息，建立所述至少一个测井仪器的响应方程组

其中i＝1，2，...，m，A_ij表示关于第i个参数的、第j种矿物组分对所述至少一个测井仪器的响应系数，x_j表示第j种矿物组分的相对含量，和B_i表示关于第i个参数的、地层对所述至少一个测井仪器的第i个响应值；以及求解所述响应方程组以获得n种矿物组份的相对含量。

Description

一种确定地层中多种矿物组分含量的方法

技术领域

本发明大体上涉及石油地质勘探领域中的常规解释处理技术，尤其是涉及一种确定地层中多种矿物组分含量的方法。

背景技术

测井，也被称为地球物理测井或石油测井，是利用岩层的电化学特性、导电特性、声学特性、放射性等地球物理特性来测量地球物理参数的方法。在油田勘探与开发过程中，测井是确定和评价油层、气层的重要手段之一，也是解决一系列地质问题的重要手段。它能直接为石油地质和工程技术人员提供各项资料和数据。

然而使用测井资料和数据对复杂岩性储层进行评价，必须首先识别岩性并且准确确定岩石矿物组分含量。只有得到准确的岩性和各种矿物组分含量，才能有针对性地采用不同的解释模型，从而得到准确的孔隙度和含油(气)饱和度。

特别是碳酸盐地层的矿物组分复杂，所述地层往往含有超过两种以上的矿物组分。现有的砂泥岩地层的处理方法(如POR、SAND等)和以双矿物组分模型为基础的复杂岩性分析方法(如CRA等)均无法有效地确定两种以上矿物组分含量。OPTIMA、ELAN等的其他一些解释处理方法所需的处理参数庞大，主要基于纯粹的数学算法，无法充分利用现有的地质信息并根据地层的实际情况进行处理。

因此，需要一种能够确定地层中多种矿物组分含量的方法。

发明内容

本发明的目的是提供一种基于现有地质信息以及相对简单的处理参数并根据地层的实际情况设置约束、高效准确地估计多种骨架组份含量的解释处理方法。

根据本发明的一方面，提供了一种确定地层中的多种矿物组分的相对含量的方法，包括：提供关于地层中的多种矿物组分的信息，所述信息包括多种矿物组分的名称和数目n，其中n为大于等于2的整数；利用至少一个测井仪器对所述地层的m个参数进行测量，其中m为大于等于n的整数；获得地层对所述至少一个测井仪器的m个响应值；基于所述m个响应值和n种矿物组分的信息，建立所述至少一个测井仪器的响应方程组

其中i＝1，2，...，m，A_ij表示关于第i个参数的第j种矿物组分对所述至少一个测井仪器的响应系数，x_j表示第j种矿物组分的相对含量，和B_i表示关于第i个参数的地层对所述至少一个测井仪器的响应值；以及求解所述响应方程组以获得n种矿物组份的相对含量。

根据本发明的另一方面，所述物理结构信息进一步包括地层中的流体的组分信息。

根据本发明的另一方面，所述m个参数从密度、中子、声波、深电阻率、浅电阻率、自然电位、自然伽马、钍、钾、光电吸收截面指数以及热中子俘获截面中进行选择。

根据本发明的另一方面，所述求解所述响应方程组包括求解通过向所述响应方程组加入约束条件而得到的带约束的线性方程组：

其中c、x_maxj均为常数，x_maxj为第j种组份的最大相对体积，其中c＝1。

根据本发明的另一方面，所述求解所述带约束的线性方程组包括求解目标函数：

根据本发明的另一方面，所述求解目标函数包括标准化处理以便使各个A_ij和B_i值都成为无量纲的数并且处于同一量级上。

根据本发明的另一方面，所述标准化处理包括根据所述m个响应值的质量在相应的方程两边乘上一个权系数以反映相应参数的响应方程对相对含量的最终结果有不同的贡献。

根据本发明的另一方面，所述求解目标函数包括用迭代法进行求解，所述迭代法包括：确定由初始点到下一迭代点的方向；确定由初始点到下一迭代点的最佳步长；以及由确定的方向和最佳步长得到下一迭代点，该点即为下次循环迭代的初始点。

根据本发明的另一方面，所述确定由初始点到下一迭代点的方向采用梯度最小方向。

根据本发明的另一方面，所述确定最佳步长采用一维最优化方法。

以上相当宽泛地概括了本发明的特征和技术优势，以便可以更好地理解本发明的以下详细描述。在下文中将描述本发明的附加特征和优势。本领域技术人员应当认识到，所公开的概念和具体实施例可以容易地用作用于修改或设计用于实现本发明的相同目的的其他结构的基础。本领域技术人员也应当认识到，这种等效的构造并不脱离所附权利要求中限定的本发明的精神和范围。

通过结合附图来阅读后面的具体实施方式，可以更好地理解本发明的特征和优点。

附图说明

现在将参照附图来解释本发明的实施例。应当注意，这些实施例用于图示基本原理，使得仅图示为了理解基本原理而必需的那些特征。附图未按比例。另外，相似标号在附图中通篇表示相似特征。

图1是根据本发明的一个实施例的确定地层中多种矿物组分含量的方法的流程图。

图2是根据本发明的一个实施例的由确定地层中多种矿物组分含量的方法对地层所分析的结果图。

图3是根据本发明的另一个实施例的由确定地层中多种矿物组分含量的方法对地层所分析的结果图。

具体实施方式

以下具体描述涉及附图，这些附图通过图示方式示出了可以实施本发明的具体细节和实施例。充分具体描述这些实施例以使本领域的技术人员能够实施本发明。可以利用其它实施例并且可以进行结构、逻辑和电改变而不脱离本发明的范围。各种实施例未必互斥，因为一些实施例可以与一个或者多个其它实施例组合以形成新实施例。

在下面的详细描述中参照形成本详细描述的一部分的附图，在所述附图中通过说明的方式示出了其中可以实践本发明的具体实施例。在这方面，参照所描述的附图的取向使用了诸如“顶”、“底”、“前”、“后”、“首”、“尾”等等的方向术语。由于实施例的组件可以被定位在若干种不同的取向中，因此所述方向术语被用于说明的目的而决不是进行限制。要理解的是，在不背离本发明的范围的情况下可以利用其他实施例并且可以做出结构或逻辑的改变。因此，不要将下面的详细描述视为限制性意义，并且本发明的范围由所附权利要求书限定。

要理解的是，除非具体另行声明，否则这里所描述的各个示例性实施例的特征可以彼此组合。

如在本说明书中所采用的术语“矿物”指的是骨架矿物和/或岩石矿物，除非上下文另外指明。

如在本说明书中所采用的术语“组分”指的是组成地层的矿物和/或流体的名称，除非上下文另外指明。

如在本说明书中所采用的术语“含量”和“相对含量”可互换地使用，一般指的是体积相对含量，除非上下文另外指明。

如在本说明书中所采用的术语“响应值”、“测井值”和“响应系数”可互换地使用，除非上下文另外指明。

下面描述本发明的思想：基于现有地质信息以及相对简单的处理参数并根据地层的实际情况设置约束、高效准确地估计多种矿物组份含量的解释处理方法，也称为地层组份分析解释处理方法。具体而言，在地层组分分析物理模型的基础上，建立各响应方程，考虑约束条件，构建反映实际测井值与理论测井值大小的目标函数并使用最优化算法求最优解，该最优解被认为是符合地层情况的最大几率解，即所求的地层矿物组分含量和流体组分含量。

接下来，具体参考附图来解释本发明的实施例。

参照图1，图1是根据本发明实施例的确定地层中多种矿物组分含量的方法100的流程图。在步骤S1，提供关于地层中的n种矿物组分的信息以建立地层矿物组分的物理模型。在步骤S2，利用至少一个测井仪器对所述地层的m个参数进行测量。在步骤S3，获得地层对所述至少一个测井仪器的m个响应值。在步骤S4，基于所述m个响应值和n种矿物组分的信息，建立所述至少一个测井仪器的响应方程组

其中i＝1，2，...，m，A_ij表示关于第i个参数的第j种矿物组分对所述至少一个测井仪器的响应系数，x_j表示第j种矿物组分的相对含量，和B_i表示关于第i个参数的地层对所述至少一个测井仪器的响应值。在步骤S5，求解所述响应方程组以获得n种矿物组份的相对含量。下面针对每个步骤详细描述本发明的方法。

含油气的储集层可以看成是由不同性质的很多种组分组成的。利用现有的有限测井信息，不能够正确地反演出组成地层的所有组分。因此，需要对储集层做以下简化处理：

把地层中物理性质相近的组分看成是同种组分。例如，绿泥石、伊利石和其他粘土组分可以统称为泥质；以及

把地层中一些含量很小的组分合并到性质与之相近的组分中。例如，当砂岩中含有少量长石时，可将长石合并到石英中，从而认为该砂岩就是纯石英砂岩。

经过以上简化处理后，可以得到地层组分的物理模型，如下表1所示。

表1.地层组分的物理模型

根据上述的地层组分的物理模型，通过各种测井仪器对地层进行测量，可以得到地层对各种测井仪器的响应值。具体而言，各种测井仪器对地层的测量包括但不限于如下测井。

自然电位测井

井下自然电场是由钻开岩层时井内钻井液的矿化度与地层水矿化度不同，井壁附近出现电化学活动产生的。自然电场的分布特点取决于井孔剖面岩层的性质。自然电位主要是测量沿井轴自然电位的变化。

声波测井

声波测井测量地层声波速度。根据声波在地层中的传播速度，可以确定地层孔隙度、岩性以及孔隙流体性质。

电阻率测井

电阻率测井测量地层电阻率。通过研究地层电阻率差异，可对岩性、储层物性、孔隙流体性质等进行分析。

自然伽马测井

自然伽马测井通过测量地层中天然放射性的高低来研究地层性质，主要是基于岩石中所含的放射性核的种类和数量不同，放射性强度也不同。岩石的自然伽马放射性水平主要取决于铀、钍、钾的相对含量。测量过程如下：给下井仪供电，探测器工作；提升下井仪经过不同地层；当伽马射线照射探测器，探测器输出相应数目的电脉冲；脉冲信号放大，传至地面，单位时间的脉冲数被转化成相应电位差值，即测量值。

中子测井

测量过程如下：通过中子源(同位素中子源)发射快中子，照射地层(与地层介质作用)减速形成热中子或超热中子，中子探测器探测热中子或超热中子的密度，不同地层，减速能力不同，计数率不同，一次寻找储集层、确定孔隙度。

密度测井

测量岩石的体积密度以及光电吸收截面指数。测井过程是由源发射0.661Mev的射线，照射地层发生康普顿效应，散射射线到达探测器，获得测井值。

因此，通过上面提及的测井(为简化起见有些测井在这里没有展开描述)并且利用各种测井仪器的响应值，可以建立线性响应方程和非线性响应方程(后者将在下面进一步描述)。所建立的线性响应方程例如分别为密度、中子、声波、自然电位、自然伽马、钍、钾、光电吸收截面指数、热中子俘获截面等等的响应方程。首先仅举几个线性响应方程的例子，其中在下面的公式中x_or，x_om，x_fw，x_bw，x_gas，x_sh，x_ma1，x_ma2，...x_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的相对含量。

密度响应方程

ρ_b＝ρ_orx_or+ρ_omx_om+ρ_fwx_fw+ρ_bwx_bw+ρ_gasx_gas+ρ_shx_sh+

                                                                (1-1)

  ρ_ma1x_ma1+ρ_ma2x_ma2+…+ρ_makx_mak

其中：ρ_b表示密度测井值(g/cm³)，即关于密度的地层对测井仪器的响应值；ρ_or，ρ_om，ρ_fw，ρ_bw，ρ_gas，ρ_sh，ρ_ma1，ρ_ma2，...ρ_mak分别表示地层中不动油(例如束缚烃)、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的体积密度值(g/cm³)。

中子响应方程

H＝H_orx_or+H_omx_om+H_fwx_fw+H_bwx_bw+H_gasx_gas+H_shx_sh+

                                                                (1-2)

 H_ma1x_ma1+H_ma2x_ma2+…+H_makx_mak

其中：H表示中子测井值，即关于中子的地层对测井仪器的响应值；H_or，H_om，H_fw，H_bw，H_gas，H_sh，H_ma1，H_ma2，...H_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的中子测井响应系数。

声波响应方程(线性)

Δt＝cp*(Δt_orx_or+Δt_omx_om+Δt_fwx_fw+Δt_bwx_bw+Δt_gasx_gas)+Δt_shx_sh+

                                                                    (1-3)

   Δt_ma1x_ma1+Δt_ma2x_ma2+…+Δt_makx_mak

其中：Δt表示声波时差测井值(us/ft)，即关于声波时差的地层对测井仪器的响应值；cp为压实校正系数；

Δt_or，Δt_om，Δt_fw，Δt_bw，Δt_gas，Δt_sh，Δt_ma1，Δt_ma2，...Δt_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的声波时差系数。

自然伽马响应方程(线性)

GR＝GR_orx_or+GR_omx_om+GR_fwx_fw+GR_bwx_bw+GR_gasx_gas+GR_shx_sh+

                                                                    (1-4)

  GR_ma1x_ma1+GR_ma2x_ma2+…+GR_makx_mak

其中：GR表示自然伽马测井值(gAPI)，即关于自然伽马的地层对测井仪器的响应值；

GR_or，GR_om，GR_fw，GR_bw，GR_gas，GR_sh，GR_ma1，GR_ma2，...GR_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的自然伽马测井响应系数(gAPI)。

钾响应方程

K＝K_orx_or+K_omx_om+K_fwx_fw+K_bwx_bw+K_gasx_gas+K_shx_sh+

                                                                    (1-5)

 K_ma1x_ma1+K_ma2x_ma2+…+K_makx_mak

其中：K为钾的测井值(％)，即关于钾的地层对测井仪器的响应值；K_or，K_om，K_fw，K_bw，K_gas，K_sh，K_ma1，K_ma2，...K_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的钾的测井响应系数(％)。

光电吸收截面指数方程

U＝U_orx_or+U_omx_om+U_fwx_fw+U_bwx_bw+U_gasx_gas+U_shx_sh+

                                                                    (1-6)

 U_ma1x_ma1+U_ma2x_ma2+…+U_makx_mak

其中：U为光电吸收截面指数测井值，即关于光电吸收截面指数的地层对测井仪器的响应值；U_or，U_om，U_fw，U_bw，U_gas，U_sh，U_ma1，U_ma2，...U_mak分别表示地层中不动油、可动油、自由水、束缚水、气、泥质、矿物组分(1～k种)的光电吸收截面指数的测井响应系数(barn/cc)。

此外根据本发明的其他实施例，还可以包括其他参数的响应方程，比如自然电位、钍、热中子俘获截面等等。

此外需要提及的是，上面是以线性响应方程为例展开的，但是根据本发明的其他实施例，也可以结合深电阻率和浅电阻率等非线性响应方程。然而，比如深电阻率响应方程和浅电阻率响应方程的非线性响应方程不同于线性方程直接参与处理。作为实例，非线性响应方程需要通过两次迭代进行修正后参与处理，具体解释如下：

深电阻率响应方程

S_wx_or+S_wx_om+(S_w-1)x_fw+(S_w-1)x_bw+S_wx_gas＝0         (1-7)

其中：S_w为含水饱和度(是未知的)，初始设置第一个采样点为1，将其转换为线性响应方程后参与解释处理得出第一个点的各组份含量后利用饱和度公式计算出该采样点的S_w，将此次计算的S_w重新代入深电阻率响应方程S_wx_or+S_wx_om+(S_w-1)x_fw+(S_w-1)x_bw+S_wx_gas＝0再参与一次处理后得到的组份含量为准确的地层组份含量。其他采样点第一次参与处理使用的S_w为上一个采样点计算出的值，第二次参与处理使用第一次处理后通过饱和度公式计算出的S_w，其中饱和度公式可以采用现有技术中的多种形式，比如Archie、Poupon、Simandoux、Indonesian、Nigerian、Fertl、Borai、Shell等等。

类似地，浅电阻率响应方程

S_xox_or+(S_xo-FI)x_om+(S_xo-1)x_fw+(S_xo-1)x_bw+(S_xo-FI)x_gas＝0  (1-8)

其中：S_xo为冲洗带含水饱和度；FI为侵入因子(数值范围为0～1)。

因此可以将上述提及的这些线性响应方程以及经过修正处理后的非线性响应方程进行组合以形成方程组，由以下式子表示：

$\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{x}_j=B_i$ (i＝1，2，...，m)           (2)

其中n表示组成地层的组分个数，

A_ij表示关于第i个参数的第j种矿物组分对第i个测井仪器(例如当一个测井仪器对应于一个测井参数时)的响应系数，

x_j表示第j种矿物组分的相对含量，

m表示由测井仪器得到的响应方程的个数(例如为测井仪器的个数，当一个测井仪器对应于一个测井参数时)，即测井参数的数目，以及

B_i表示关于第i个参数的地层对第i个测井仪器(例如当一个测井仪器对应于一个测井参数时)的响应值。

需要注意的是，在上面式子(2)中等式左边为理论测井值而右边为实际测井值(即响应值)。注意，由于它们二者都随着地层深度而变化，所以实际上可能得到的是相应的理论测井曲线和实际测井曲线。

另外，在上面式子(2)中每个测井仪器可以对应于一个测井参数，但是实际情况可以不限于此，即单个测井仪器可以对多个测井参数进行测量并且因此对应于多个测井参数。

通过求解以上由m个方程组成的方程组，就可以求得x_j，即地层中的矿物组分的相对含量，这就是待解决的反演问题。

对于上面式子(2)而言：当m＜n时，方程组有多个解，无实际意义；当m＝n时，以上方程组有唯一解；但是为了充分利用测井信息并且提高测井解释的可靠性，一般情况下m＞n，此时方程组为超定线性方程组，它具有一最优解。这一最优解

有可能出现x_j ^*＜0或x_j ^*＞1的现象，这种结果在地质上是不存在的或无意义的。为了使求解的结果合乎地质意义并符合地层实际情况，需要在式子(2)中加入有关的约束条件，即：

写成更一般的形式

其中：c、x_maxj均为常数，x_maxj为第j种组份的最大相对体积。例如在地层组份分析方法中可假设c＝1，当然也可以取范围为0到1的其他数值。

式子(4)为一个带约束的超定线性方程组。

由线性最小二乘原理，解式子(4)这一带约束的线性方程组的问题可转换成以下求极值问题：

由于不同的测井值的量纲不一样，而且它们的测井值大小的差别也很大。因此在实际计算中可以将式子(5)的目标函数的系数A及B进行标准化处理，以便使各种仪器的A和B值都成为无量纲的数并在同一数量级上，这样可使得各种测井方法对最终结果具有相同的贡献。标准化处理的方法是将方程的两边同时除以一个系数P(式子6)。

$P=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{R}_i^2}---\left(6\right)$

其中：n为地层组份个数，Ri为第i种组份对本测仪器的响应系数。

实际中，根据测井数据的质量好坏，在方程两边乘了一个权系数，本权系数隐含为1，其大小可由用户修改，是测井解释参数。

当n≤m时，线性方程组即式子(4)的最小二乘解是唯一的。式子(5)构成了地层组份分析方法的数学模型，为目标函数。

下面详细解释根据本发明实施例的对目标函数求解过程。

用迭代法确定n维空间函数最小值一般有三大步：第一步，确定由初始点到下一迭代点的方向；第二步，确定由初始点到下一迭代点的最佳步长；第三步，由确定的方向和最佳步长得到下一迭代点，该点即为下次循环迭代的初始点。不同的最优化算法确定下一迭代点方向的方法是不同的，本算法用梯度最小方向，但确定最佳步长的方法一般相同，本算法用一维最优化方法，即最佳步长处目标函数值最小。

现假定方程组即式子(4)已进行标准化处理(注意，这不是必需的)，下面讨论极小化问题即式子(5)的求解方法。

第一步：

任取解空间R中的一个点

在

处

的线性逼近函数为

$f_L\left(\stackrel{\→}{x}\right)=f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)+{[\▿f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)]}^T(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(0\right)})---\left(7\right)$

其中： $\▿f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)=[\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_1},\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_2},...,\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_n}]$

显然因均与无关，求线性规划问题

的最优解，等价于求线性规划问题

的最优解。令

为式(5)的最优解，

$\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_i}=[\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_1},\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_2},...,\frac{\∂f\left({\stackrel{\→}{x}}^{\left(0\right)}\right)}{{\∂x}_n}],$ 由线性规划的性质

可知，

必为R的一个顶点，因此可得

$x_{\mathrm{FLj}}^{\left(0\right)}=\left\{\begin{array}{cc}0& (j\≠i)\\ x_{\max j}& (j=i)\end{array}\right.,(j=\mathrm{1,2},...,n)---\left(8\right)$

与

的连线即为搜索方向。

第二步：

下面分两种情况讨论：

(1)当 ${[\▿f\left({\stackrel{\→}{x}}^{\left(0\right)}\right]}^T({\stackrel{\→}{x}}_{\mathrm{FL}}^{\left(0\right)}-{\stackrel{\→}{x}}^{\left(0\right)})=0$ 时(实际中，

(与

两点目标函数值的相对误差)来判断，若D足够小，小于ε，则条件成立)，就是线性规划问题的解，迭代停止。

(2)当 ${[\▿f\left({\stackrel{\→}{x}}^{\left(0\right)}\right]}^T({\stackrel{\→}{x}}_{\mathrm{FL}}^{\left(0\right)}-{\stackrel{\→}{x}}^{\left(0\right)})\≠0$ 时，则问题变为极值问题：

$\underset{\mathrm{\λ}\∈\left[\mathrm{0,1}\right]}{\min }f[{\stackrel{\→}{x}}^{\left(0\right)}+\mathrm{\λ}({\stackrel{\→}{x}}_{\mathrm{FL}}^{\left(0\right)}-{\stackrel{\→}{x}}^0)]---\left(9\right)$

的最优解λ₀，这时必有0≤λ₀≤1。

即求当λ为多大时目标函数在处的值最小，该值即为最佳步长λ₀，最佳步长λ₀对应的点

就是下次迭代初始点。这是个一维最优化问题，不同于式子(5)的多维最优化问题。

本一维极值问题用高等数学中的求导法解决。今目标函数对λ的导数为0可得到

${\mathrm{\λ}}_0=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[B_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}({\stackrel{\→}{x}}_{\mathrm{FLj}}^{\left(0\right)}-{\stackrel{\→}{x}}_j^{\left(0\right)})-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\left(A_{\mathrm{ij}}{\stackrel{\→}{x}}_j^{\left(0\right)}\right)\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}({\stackrel{\→}{x}}_{\mathrm{FLj}}^{\left(0\right)}-{\stackrel{\→}{x}}_j^{\left(0\right)})}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}({\stackrel{\→}{x}}_{\mathrm{FLj}}^{\left(0\right)}-{\stackrel{\→}{x}}_j^{\left(0\right)})\right)}^2}---\left(10\right)$

第三步：

令： ${\stackrel{\→}{x}}^{\left(1\right)}={\stackrel{\→}{x}}^{\left(0\right)}+{\mathrm{\λ}}_0({\stackrel{\→}{x}}_{\mathrm{FL}}^{\left(0\right)}-{\stackrel{\→}{x}}^{\left(0\right)}),$ 把

作为

继续用上述方法线性逼近目标函数

并重复以上步骤，直到满足精度为止，就可求得带约束线性超定方程组(4)的解。

在极小值问题(5)中，很容易写出

的一阶导数，因此，寻优计算工作量小，并且目标函数

不存在多个局部极小点，所以对迭代初始值的要求不严，不管迭代初始值怎样选择，在有限步内，总可收敛到同一极小点。

实验结果

实例1

参照图2，图2是根据本发明的方法所分析的结果图。该结果图是使用声波响应方程、自然伽马响应方程、钾响应方程、中子响应方程、密度响应方程、深电阻率响应方程以及浅电阻率响应方程而获得的，其中生成该结果图使用的响应系数，如下表I所示。另外，根据实际需要可以选择上面提及的和/或其他未具体描述的各种响应方程。

表I.实例1中使用的各种组分的响应系数

如图2所示，该图总共有7列，分别代表7个道：深度，DTC(声波曲线道)，GR(自然伽马曲线道)，POTA(钾曲线道)，NPHI(中子曲线道)，RHOB(密度曲线道)以及组分分析(示出确定的组分含量)，这将在下面进行进一步的解释。注意，为清楚起见，可能没有示出所有方程的曲线道，仅选择性地示出其中上面提及的7个道。

如图2中的最右边的组分分析道所示，本发明的方法确定了多种矿物组分的相对含量，具体地是泥质、石英、钾长石、斜长石、方解石、白云石和角闪石的相对含量，如下表2所示，其中特殊矿物1为钾长石，特殊矿物2为斜长石，特殊岩石1为角闪石。然后根据上表2，进而可以得出矿物组分含量的平均值，如下表3所示。

表2.通过本发明的方法确定的矿物组分含量(单位：％)

表3.通过本发明的方法确定的矿物组分含量的平均值(单位：％)

从表3可以看出，方解石和白云石的相对含量基本上为零，因此为清楚起见在图2中的组分分析道中没有示出方解石和白云石的相对含量曲线。

下表4是3100m-3110m的岩心分析数据，代表的是该井段的实际矿物组分含量的平均值。

表4.实际矿物组分含量的平均值(单位：％)

对比表3与表4，可以看出使用本发明的方法确定的结果与岩心分析结果接近，即符合该地层的地质特征。

然后通过将确定的各组分含量代入上述提及的各种响应方程，可以得出测井理论曲线，即如图2中的第2列到第6列中的相应虚线DTCT、GRT、POTAT、NPHIT和RHOBT所示。在图2中的第2列到第6列中的实线分别代表原始测井曲线DTC、GR、POTA、NPHI和RHOB，而虚线分别代表在将确定的各组分含量代入相应的响应方程后得出的测井理论曲线DTCT、GRT、POTAT、NPHIT、RHOBT。

从图2可以看出，将确定的各组分含量代入相应的响应方程后得出的测井理论曲线DTCT、GRT、POTAT、NPHIT、RHOBT与原始测井曲线DTC、GR、POTA、NPHI、RHOB的吻合度非常好，进一步说明由本发明的方法确定的组分含量是准确的，符合该地层的地质特征。

下面再次参照图2，具体解释图2中的各列中的相关内容：

第一列深度道示出处理深度(单位为m)，例如3100m到3130m；

第二列DTC(声波道)示出声波原始测井曲线DTC和将确定的各组分含量代入相应的响应方程后得出的测井理论曲线DTCT，单位为US/FT(微妙/英尺)，示出刻度为左刻度140，右刻度40，[DP0]为系统标识，无实际物理含义；

第三列GR(自然伽马道)示出自然伽马原始测井曲线GR和将确定的各组分含量代入相应的响应方程后得出的测井理论曲线GRT，单位为GAPI(自然伽马的单位)，示出刻度为左刻度0，右刻度200，[DP0]为系统标识，无实际物理含义；

第四列POTA(钾道)示出钾原始测井曲线POTA和将确定的各组分含量代入相应的响应方程后得出的测井理论曲线POTAT，单位为％，示出刻度为左刻度0，右刻度10，[DP0]为系统标识，无实际物理含义；

第五列NPHI(中子道)示出中子原始测井曲线NPHI和将确定的各组分含量代入相应的响应方程后得出的测井理论曲线NPHIT，单位为none(空)，示出刻度为左刻度0.45，右刻度-0.15，[DP0]为系统标识，无实际物理含义；

第六列RHOB(密度道)示出密度原始测井曲线RHOB和将确定的各组分含量代入相应的响应方程后得出的测井理论曲线RHOBT，单位为g/cm³(克/立方厘米)，示出刻度为左刻度2，右刻度3.5，[DP0]为系统标识，无实际物理含义；以及

第七列示出确定的各组分含量，从左到右依次为泥质、特殊岩石1(在本实例中为角闪石)、石英、特殊矿物1(在本实例中为钾长石)、特殊矿物2(在本实例中为斜长石)、可动油、束缚烃以及水的相对含量。

实例2

参照图3，图3是根据本发明的方法所分析的另一结果图。该结果图是使用声波响应方程、自然伽马响应方程、中子响应方程、密度响应方程、深电阻率响应方程，浅电阻率响应方程以及光电吸收截面指数响应方程而获得的，其中生成该结果图使用的响应系数，如下表II所示。另外，根据实际需要可以选择上面提及的和/或其他未具体描述的各种响应方程。

表II.实例2中使用的各种组分的响应系数

如图3所示，该图共有6列，分别代表6个道：深度，RT(深电阻率曲线道)，DTC(声波曲线道)，NPHI(中子曲线道)，RHOB(密度曲线道)，组份分析(示出确定的组份含量)，这将在下面进行进一步的解释。注意，为清楚起见，可能没有示出所有方程的曲线道，仅选择性地示出其中上面提及的6个道。

如图3中的最右边的组分分析道所示，本发明的方法确定了多种矿物组分的相对含量，具体地是泥质，石灰岩，砂岩，硬石膏，白云岩和斜长石的相对含量，如下表5所示。然后根据表5，进而可以得出矿物组分含量的平均值，如下表6所示。

表5.通过本发明的方法确定的矿物组分含量(单位：％)

表6.通过本发明的方法确定的矿物组分含量的平均值(单位：％)

从表6可以看出，斜长石的相对含量基本上为零，因此为清楚起见在图3中的组分分析道中没有示出斜长石的相对含量曲线。

下表7是2190m-2194m的岩心分析数据，代表的是该井段的实际矿物组分含量的平均值。

表7.实际矿物组分含量的平均值(单位：％)

对比表6与表7，可以看出使用本发明的方法确定的结果与岩心分析结果接近，即符合该地层的地质特征。

然后通过将确定的各组分含量代入上述提及的各种响应方程，可以得出测井理论曲线，即如图3中的第2列到第5列中的相应虚线RTT、DTCT、NPHIT和RHOBT所示。在图3中的第2列到第5列中的实线分别代表原始测井曲线RT、DTC、NPHI和RHOB，而虚线分别代表在将确定的各组分含量代入相应的响应方程后得出的测井理论曲线RTT、DTCT、NPHIT和RHOBT。

从图3可以看出，将确定的各组份含量代入相应的响应方程后得出的测井理论曲线RTT、DTCT、NPHIT、RHOBT与原始测井曲线RT、DTC、NPHI、RHOB的吻合度非常好，进一步说明由本发明的方法确定的组分含量是准确的，符合该地层的地质特征。

下面再次参照图3，具体解释图3中的各列中的相关内容：

第一列深度道示出处理深度(单位为m)，例如2190m到2210m；

第二列RT(深电阻率曲线道)示出深电阻率原始测井曲线RT和通过计算得出的测井理论曲线RTT，单位为ohm.m(欧姆.米)，示出刻度为左刻度0.2，右刻度2000，[DP0]为系统标识，无实际物理含义；

第三列DTC(声波曲线道)示出声波原始测井曲线DTC和将确定的各组份含量代入相应的测井响应方程后得出的测井理论曲线DTCT，单位为us/ft(微妙/英尺)，示出刻度为左刻度140，右刻度40，[DP0]为系统标识，无实际物理含义；

第四列NPHI(中子道)示出中子原始测井曲线NPHI和将确定的各组份含量代入相应的测井响应方程后得出的测井理论曲线NPHIT，单位为空，示出刻度为左刻度0.45，右刻度-0.15，[DP0]为系统标识，无实际物理含义；

第五列RHOB(密度道)示出密度原始测井曲线RHOB和将确定的各组份含量代入相应的测井响应方程后得出的测井理论曲线RHOBT，单位为g/cm³(克/立方厘米)，示出刻度为左刻度1.95，右刻度2.95，[DP0]为系统标识，无实际物理含义。

第六道示出确定的各组份含量，从左到右依次为泥质，石灰岩，砂岩，硬石膏，白云岩，可动油，束缚烃以及水的相对含量。

注意，上面列举的矿物组分以及各列参数都仅仅是作为实例给出的，它们可以根据实际应用以及需要加以改变和调整。例如，各列参数可以例如包括但不限于密度、中子、声波、深电阻率、浅电阻率、自然电位、自然伽马(线性、非线性)、钍、钾、光电吸收截面指数、热中子俘获截面中的一个或多个。

综上所述，本发明例如具有如下优点：可以利用现有地质信息，参数处理相对简单，根据地层的实际情况设置约束，高效准确地估计多种骨架组份含量。

虽然以上参考附图示出并描述了本发明的各种示例性实施例，但是本领域技术人员将清楚可以做出将实现本发明一些优点的各种改变和修改而未脱离本发明的精神和范围。因此，只要本发明的这些修改、变型和替换属于本发明的权利要求书及其等同技术方案的范围之内，则本发明也意图包含这些修改、变型和替换。此外，就在详细描述或权利要求书中使用术语“包括”、“具有”、“带有”或它们的其他变体来说，这样的术语意图以与术语“包含”类似的方式是包含性的。

应当理解的是，本发明的实施例可以被实施在分立电路、部分集成电路或完全集成电路或者编程装置中。此外，术语“示例性”仅仅意味着作为一个实例，而不是最佳的或最优的。还要认识到，这里所描绘的特征和/或元素是出于简单及易于理解的目的而以相对于彼此的特定尺寸示出的，实际的尺寸可能与这里所示出的尺寸显著不同。

本领域合理技术人员将明白，可以用执行相同功能的其它部件适当地进行替换。应当提到，即使在尚未明确提到这一点的那些情况下，参照具体图所解释的特征可以与其它图的特征组合。另外，可以在使用适当处理器指令的全软件实施方式中或者在利用硬件逻辑与软件逻辑的组合以实现相同结果的混合实施方式中实现本发明的方法。这样对发明概念的修改旨在于由所附权利要求覆盖。

此外，本申请的范围不旨在限于说明书中描述的过程、机器、制造、物质组成、装置、方法和步骤的具体实施例。如本领域普通技术人员从本发明的公开中容易认识到的，根据本发明可以利用与这里描述的对应实施例执行实质上相同的功能或实现实质上相同的结果的、当前存在或以后待开发的过程、机器、制造、物质组成、装置、方法或步骤。据此，所附权利要求旨在在其范围内包括这些过程、机器、制造、物质组成、装置、方法或步骤。

Claims (10)

1.一种确定地层中的多种矿物组分的相对含量的方法，包括：

提供关于地层中的多种矿物组分的信息，所述信息包括多种矿物组分的名称和数目n，其中n为大于等于2的整数；

利用至少一个测井仪器对所述地层的m个参数进行测量，其中m为大于等于n的整数；

获得地层对所述至少一个测井仪器的m个响应值；

基于所述m个响应值和n种矿物组分的信息，建立所述至少一个测井仪器的响应方程组

$\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{x}_j=B_i$

其中每个参数对应于一个响应方程，i＝1，2，...，m，

A_ij表示关于第i个参数的、第j种矿物组分对所述至少一个测井仪器的响应系数，

x_j表示第j种矿物组分的相对含量，和

B_i表示关于第i个参数的、地层对所述至少一个测井仪器的第i个响应值；以及

求解所述响应方程组以获得n种矿物组份的相对含量。

2.根据权利要求1的方法，其中所述信息进一步包括地层中的流体组分的信息。

3.根据权利要求1的方法，其中所述m个参数从密度、中子、声波、深电阻率、浅电阻率、自然电位、自然伽马、钍、钾、光电吸收截面指数以及热中子俘获截面中进行选择。

4.根据权利要求1的方法，其中所述求解所述响应方程组包括求解通过向所述响应方程组加入约束条件而得到的带约束的线性方程组：

其中c、x_maxj均为常数，x_maxj为第j种组份的最大相对体积，其中c＝1。

5.根据权利要求4的方法，其中所述求解所述带约束的线性方程组包括求解目标函数：

6.根据权利要求5的方法，其中所述求解目标函数包括标准化处理以便使各个A_ij和B_i值都成为无量纲的数并且处于同一量级上。

7.根据权利要求6的方法，其中所述标准化处理包括根据所述m个响应值的质量在相应的响应方程两边乘上一个权系数以反映相应参数的响应方程对所述相对含量的最终结果有不同的贡献。

8.根据权利要求5的方法，其中所述求解目标函数包括用迭代法进行求解，所述迭代法包括：确定由初始点到下一迭代点的方向；确定由初始点到下一迭代点的最佳步长；以及由确定的方向和最佳步长得到下一迭代点，该点即为下次循环迭代的初始点。

9.根据权利要求8的方法，其中所述确定由初始点到下一迭代点的方向采用梯度最小方向。

10.根据权利要求8的方法，其中所述确定最佳步长采用一维最优化方法。

Images