CN112731541A - 核磁回波数据优化反演方法及孔隙结构参数确定方法 - Google Patents

Abstract

本发明提供了核磁回波数据优化反演方法及孔隙结构参数确定方法，建立球缝模型，根据球缝模型计算得到新的核磁回波数据反演的布点；扫描确定的第一特定参数，第一特定参数为球缝模型中裂缝的半径宽度与球缝模型中球形孔的半径的比值；根据每个第一特定参数路径下利用布点确定的T₂谱时间核磁回波数据反演的布点值，反演核磁回波数据，得到不同的T2谱分布，计算T₂谱分布与回波信号之间拟合协方差，其中，协方差最小的T2谱分布，就是优化的T₂谱分布。本发明降低了核磁测井回波数据的反演误差；本发明提供的孔隙结构参数确定方法充分提取了核磁数据中蕴含的关于岩石孔隙结构的信息，实现了孔隙结构参数的定量计算。

Description

核磁回波数据优化反演方法及孔隙结构参数确定方法

技术领域

本发明属于常规油气藏地球物理勘查技术，具体涉及页岩气核磁共振测井资料回波数据反演及页岩气储层渗透率计算算法设计。

背景技术

就已经发表的研究成果而言，核磁共振测井和岩芯分析数据是研究人员认识复杂储层岩芯及流体类型的有利资料，这些数据在研究储层孔隙结构及识别孔隙流体类型方面具有独特的优势。但是，核磁共振测井和岩芯核磁实验的直接观测数据是核磁回波，核磁回波数据本身与岩石孔隙结构和孔隙流体有关，但是，核磁回波数据并不能以直观的方式展示这种关系，所以，在实际应用中，很少直接采用回波数据研究岩石的孔隙结构和流体性质。常用的方式是首先采用特定的数学方法，例如奇异值分解(SVD)、稳固反复迭代算法(SIRT)将核磁回波数据反演成T₂谱分布，这种T₂谱分布反映了岩石孔隙结构的数量、大小形态与横向弛豫时间之间的关系。因此，T₂谱分布才是研究人员用于储层岩石孔隙结构识别、流体识别直接资料。但是，比较麻烦的问题是，将核磁回波反演成T₂谱分布的问题是高度病态的，反演系数矩阵的条件数可以到10¹⁶数量级。高度病态问题对于反演结果的影响是，回波数据中微小的扰动，都可能导致反演T₂谱分布的极大变化。

为了避免核磁反演的多解问题，有学者提出“先滤波、再反演了”策略，先对回波做滤波平滑处理，消除部分随机噪声干扰，再对滤波后的回波数据做反演。这个策略的确可以压制部分随机噪声对于反演结果的干扰，但是，其缺点是不可避免也会压制部分有用信号，使得反演结果对于孔隙结构和流体类型的刻画能力降低。另外，现有的核磁数据回波反演，一般只反演一次，就输出反演结果。这样的反演结果很可能陷于局部极值，而无法获得全局极值。局部极值可能无法全面刻画岩石的真实孔隙结构和流体类型，甚至是虚假的信息。

为了避免以上的问题，充分挖掘核磁共振测井和岩芯分析数据蕴藏的信息，有必要研究新的、能更有效客观反演回波数据的新方法。

发明内容

本发明针对目前核磁回波数据反演方法存在的技术问题，提供一种核磁回波数据优化反演方法。

本发明采用以下技术方案。

本发明提供一种核磁回波数据优化反演方法，包括以下步骤：

获取核磁回波数据，所述核磁回波数据与样本的孔隙度线性相关；

建立球缝模型确定球缝模型的几何关系，包括球冠高、球冠半径、球冠面积和球冠体积；采用球缝模型优化反演程序得到优化的T₂谱分布，所述T₂谱分布与特定的球缝模型相关；所述球缝模型优化反演程序包括以下步骤：

根据球缝模型的几何关系确定初始核磁回波数据反演布点及等效球的半径，根据球缝模型计算得到新的核磁回波数据反演的布点；扫描确定的第一特定参数，所述第一特定参数为球缝模型中裂缝的半径宽度与球缝模型中球形孔的半径的比值；根据每个第一特定参数路径下利用所述布点确定的T₂谱时间核磁回波数据反演的布点值，反演核磁回波数据，得到不同的T2谱分布，计算T₂谱分布与回波信号之间拟合协方差，其中，协方差最小的T2谱分布，就是优化的T₂谱分布。

第二方面，为了获得更加精确的孔隙结构参数，本发明提供了孔隙结构参数确定方法，包括采用第一方面任意一种可能的实施例提供的核磁回波数据优化反演方法获得核磁回波数据的优化的T₂谱分布；得到的最优T₂谱进行球形孔和裂缝孔的分解；在完成最优谱的分解之后，分别得到球形孔和裂缝孔分布；分别做球形孔和裂缝谱的累计曲线，再计算不同频率处的孔隙半径值，或裂缝宽度值，并使用分选系数、半径均值或峰态参数的定义计算孔隙结构参数。

本发明所取得的有益技术效果：

本发明提供的核磁回波数据优化反演方法降低了核磁测井回波数据的反演误差；

本发明利用提供的核磁回波数据优化反演方法获得反演数据，并基于反演数据提供了孔隙结构参数确定方法充分提取了核磁数据中蕴含的关于岩石孔隙结构的信息，实现了孔隙结构参数的定量计算。

附图说明

图1为本发明具体实施例中回波数据反演误差对比分析结果示意图；

图2为本发明具体实施例中计算储层岩石孔隙结构参数结果示意图；图3为本发明具体实施例核磁共振测井球缝模型优化反演处理结果示意图；

图4为本发明具体实施例不同孔隙结构导致的岩石电阻率差异；

图5为本发明具体实施例球缝模型示意图；

图6为本发明具体实施例扫描确定C_d值路径图；

图7为本发明具体实施例实现流程图。

具体实施方式

以下结合说明书附图和具体实施例对本发明做进一步说明。

实施例一、核磁回波数据优化反演方法，包括以下步骤：

使用孔隙度刻度模块，刻度核磁测井仪器，或核磁共振岩芯分析仪器，采集储层的核磁共振测井数据，或实验室岩芯的核磁共振测井数据，并确保仪器测量的核磁测井回波数据与储层岩石，或实验室核磁仪器的回波数据与岩石样品的孔隙度线性相关；

采用球缝模型优化反演程序，反演核磁共振测井或核磁共振岩芯分析数据，得到T₂谱分布，且这个T₂谱分布与特定的球缝模型相关。

本实施例球缝模型优化反演程序包括如下步骤(如图7所示)：

第一步，建立球缝模型(附图5)。即一个球形孔与一个裂缝组合，构成球缝模型。

第二步，推导球缝模型的几何关系与参数的计算方程(方程1-方程28)。

球冠高：h＝R_s-R_f＝(1-C_d)R_s (1)

其中，h,球冠的高度，μm；R_s，R_f，分别为球形孔的半径和裂缝的半宽度，单位均为μm。C_d的含义由方程(2)定义。

球冠半径：

r,球冠半径，μm。

球冠表面积：S_s＝2π×r×h (4)S_S，球冠表面积，μm²。

球冠体积：

初始核磁回波数据反演布点(T_2i)的计算方法及Re的计算，采用方程(6)计算。

其中，T_2max,T_2min,分别是反演布点的最大值和最小值，ms。例如，可以分别取为3000ms和0.3ms；

T_2i,反演布点的当前分组弛豫时间，ms；

N，i，分别是反演布点的个数和反演布点的当前编号。

ρ₂，横向弛豫率，cm/ms。

R_e，是等效球半径，μm。

b＝16R_e·C_d·R_s (9)

a,b,c,分别由方程(8)，(9)，(10)定义的中间参数，

S_S，S_f，分别是球形孔和裂缝的表面积，μm²；

V_s，V_f，分别是球形孔和裂缝的体积，μm³；

C_d＝1时，球缝模型退化为纯裂缝，模型方程由方程(16)定义，

R_s＝0 (18)

V_s＝S_s＝0 (19)

C_d＝0时，球缝模型退化为纯球形孔，据此可以得到模型方程(方程23)，

R_f＝0 (21)

R_s＝R_e (24)

V_f＝S_f＝0 (25)

其中，ρ₂为横向弛豫率，取为0.022μm/ms。

T₂(R_ei,C_di),就是根据球缝模型计算得到的新弛豫时间布点，ms。

在根据方程(11)计算球半径时，唯一的未知数是C_d。确定C_d后，R_s，R_f，Vs，V_f，S_s，S_f均可计算得到。最后，根据方程(28)计算T_2i，也就是核磁回波数据反演的布点公式。

第三步，扫描确定C_d。

参见球缝模型设计示意图(附图4，即由球形孔和裂缝组合构成的储层岩石孔隙模型)，可以知道，C_d的变化区间是[0,1]。0对应于纯球形孔，1对应于纯裂缝孔。当C_d位于0,1之间时，就是缝孔隙与球孔隙以某种形式叠加形成的球缝模型。确定孔隙结构的关键问题是确定C_d的数值，理论上，每个组分内的球管模型理论上都有无数种，为了方便计算，提出了计算孔隙结构参数的扫描方法(附图6)。

在每个球管模型的分组内部，设定缝宽半径与球半径的比值。对以上图形进行扫描，得到描述岩石孔隙特征的路径。每一个路径上，任何一点的横坐标就是核磁反演的布点序号，每一个布点序号对应于一个弛豫时间布点。任何一点的纵坐标就是C_d值。

第四步，确定最优反演T₂谱及孔隙结构参数计算

根据每个C_d路径下计算的T₂时间反演布点值，反演核磁回波数据，得到不同的T₂谱分布，计算T₂谱分布与回波信号之间拟合协方差。其中，协方差最小的T₂谱，就是优化的T₂谱。

实施例二、本实施例以核磁共振测井和核磁共振岩芯分析数据为基础，采用本次研究提出的球缝模型优化反演算法，优化反演核磁共振测井，或核磁共振岩芯分析的回波数据。根据优化反演结果，得到页岩油气储层中孔隙结构与球形孔、板形孔最近似的结构参数。然后，根据这个最近似的结构参数，计算储层岩石中球形孔均值半径、裂缝均值宽度、球形孔隙度、裂缝孔隙度、球形孔分选系数及裂缝孔隙的分选系数等孔隙结构参数。

就球形孔、裂缝孔的力学特性而言，球形孔具有最大的抗压强度，在裂缝孔的垂直方向压裂才会取得最好的压裂效果。基于这样的力学原理，在了解岩石孔隙结构的基础上，有利于制定合理高效的压裂方案。

本实施例提供了孔隙结构参数确定方法，本实施例在实施例一的基础上利用实施例一提供的核磁回波数据优化反演方法，获得核磁回波数据的优化的T₂谱分布，根据最优T₂谱对应的C_d路径，即可计算储层岩石的孔隙结构参数，即球半径均值，缝宽均值，球形孔隙度，缝孔隙度，球形孔分选系数，裂缝孔隙度等孔隙结构参数。

对于已经得到的最优T₂谱进行球形孔和裂缝孔的分解(方程29)。

S_ti＝S_si+S_fi

其中，A_i，A_si，A_fi，分别是最优谱分布幅度，球形孔谱分布幅度，裂缝孔谱分布幅度，标准化数据无量纲；

S_ti，S_si，S_fi，分别是最优谱的面积，球形孔面积，裂缝孔面积，μm²；

在完成最优谱的分解之后，可以分别得到球形孔和裂缝孔分布(附图3)。分别作球形孔和裂缝谱的累计曲线，再计算不同频率处的孔隙半径值，或裂缝宽度值，并使用分选系数、半径均值、峰态等参数的定义计算孔隙结构参数(方程30-方程34)。

例如，从连续曲线上读出当频率等于0.5时对应的球形孔半径，写成R_s50。

R_s50就是球形孔半径中值，μm。

依次类推，读出频率分别为0.84，0.16，0.95，0.05对应的球形孔半径，分别记为R_s84，R_s16，R_s95，R_s5，μm。也就是累计频率对应于0.84,0.16,0.95,0.05处，球形孔具有半径。

球形孔的分选系数计算公式(S_p)

S_p,球形孔分选系数，μm。

球形孔的歪度为(SK_p)

SK_p,球形孔歪度，无量纲。

球形孔的均值半径D_M(μm)

DM，球形孔均值半径，μm。

球形孔的峰态(Kp)

Kp，球形孔峰态，无量纲。

球形孔的结构系数(F_y)

F，孔隙结构系数，无量纲；

平均孔隙半径，μm；

孔隙度，小数，K，渗透率，μm²；

本实施例根据特定的球缝模型，计算储层岩石最近似的孔隙结构参数，包括球形孔、裂缝孔的均值半径和均值宽度，球形孔、裂缝孔孔隙度，球形孔、裂缝孔的分选系数、峰态等孔隙结构参数。

在以上实施例的基础上，根据球缝模型的优化反演结果，参考岩芯分析的渗透率数据，确定束缚球形孔、束缚裂缝孔与自由球形孔、自由裂缝孔之间的弛豫时间划分界限，把优化反演得到的球形孔谱分布和裂缝谱分布分别划分成束缚球形孔、束缚裂缝孔和自由球形孔、自由裂缝孔2个部分。

应用本发明得到的岩石孔隙结构参数，可以计算孔隙结构导致的岩石电阻率变化。在此基础上研发新的饱和度模型，提高储层饱和度的计算精度(附图4)。

首先，本发明降低了核磁测井回波数据的反演误差。以Hechuan5井的核磁测井回波数据为例，对比常规反演算法和优化反演算法的误差，并进行误差分析(附表1，附图1)。在该井的11个测井回波数据点的计算结果表明，优化反演的拟合均方差比常规反演的拟合均方差低12％左右。对于更多数据计算结果，结论依然成立。

表1回波反演的误差分析

其次，充分提取了岩芯核磁数据和测井核磁数据中蕴含的关于岩石孔隙结构的信息，实现了孔隙结构参数的定量计算(附图2)。

基于本发明算法编写了测井数据计算机处理程序，并完成了实际测井数据的数值处理，得到了深度剖面上连续的核磁测井数据处理结果。这些结果参数包括总T2谱分布，球形孔T2谱，均值球形孔半径，均值裂缝宽度，球形孔分选系数，裂缝分选系数等孔隙结构参数，即实现了本发明优化反演算法的实际测井数据处理的应用(附图3)。

应用本发明计算得到的岩石孔隙结构参数，可以计算分析孔隙结构导致的岩石电阻率差异，在此基础上，研究岩石的饱和度模型，提高储层饱和度的计算精度(附图4)。

本实施例以核磁共振测井和核磁共振岩芯分析数据为基础，采用本次研究提出的球缝模型优化反演算法，优化反演核磁共振测井，或核磁共振岩芯分析的回波数据。根据优化反演结果，得到页岩油气储层中孔隙结构与球形孔、板形孔最近似的结构参数。然后，根据这个最近似的结构参数，计算储层岩石中球形孔均值半径、裂缝均值宽度、球形孔隙度、裂缝孔隙度、球形孔分选系数及裂缝孔隙的分选系数等孔隙结构参数。

采用优化反演的方法得到岩石的孔隙结构参数以后，可以用于这些孔隙结构参数，研究孔隙结构对于岩石导电特征的影响，并研究新的饱和度模型，提高饱和度模型的数值计算精度。另外，还可以基于这些孔隙结构参数，研究岩石的力学参数，对于致密储层岩石和页油岩岩石，有利于制定合理高效的压裂方案。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (9)

1.核磁回波数据优化反演方法，其特征在于，包括以下步骤：

获取核磁回波数据，所述核磁回波数据与样本的孔隙度线性相关；

建立球缝模型确定球缝模型的几何关系，包括球冠高、球冠半径、球冠面积和球冠体积；采用球缝模型优化反演程序得到优化的T₂谱分布，所述T₂谱分布与特定的球缝模型相关；所述球缝模型优化反演程序包括以下步骤：

根据球缝模型的几何关系确定初始核磁回波数据反演布点及等效球的半径，根据球缝模型计算得到新的核磁回波数据反演的布点；扫描确定的第一特定参数，所述第一特定参数为球缝模型中裂缝的半径宽度与球缝模型中球形孔的半径的比值；根据每个第一特定参数路径下利用所述布点确定的T₂谱时间核磁回波数据反演的布点值，反演核磁回波数据，得到不同的T2谱分布，计算T₂谱分布与回波信号之间拟合协方差，其中，协方差最小的T2谱分布，就是优化的T₂谱分布。

2.根据权利要求1所述的核磁回波数据优化反演方法，其特征在于，确定球缝模型的几何关系的计算公式如下：

h＝R_s-R_f＝(1-C_d)R_s (1)

S_s＝2π×r×h (4)

其中h是球冠高；R_s是球缝模型中球形孔的半径，C_d是第一定义参数，无量纲；R_f是球缝模型中裂缝的半宽度；r是球冠底圆的半径；S_s是球冠面积，V_s是球冠体积，V_s是球冠体积；

根据球缝模型的几何关系确定等效球的半径及初始核磁回波数据反演布点的计算公式如下：

其中，Re是等效球的半径，T_2max，T_2min，分别是反演布点的最大值和最小值，T_2i为核磁回波数据反演的第i个布点的当前分组弛豫时间，N为反演布点的个数，i为反演布点的当前编号，ρ₂为横向弛豫率。

3.根据权利要求2所述的核磁回波数据优化反演方法，其特征在于，根据球缝模型计算得到新的核磁回波数据反演的布点的方法如下：

b＝16R_e·Cd·R_s (9)

a，b，c，分别由方程(8)，(9)，(10)定义的中间参数，

S_S，S_f，分别是球形孔和裂缝的表面积，μm²；

V_s，V_f，分别是球形孔和裂缝的体积，μm³；

C_d＝1时，球缝模型退化为纯裂缝，模型方程由方程(16)定义，

R_s＝0 (18)

V_s＝S_s＝0 (19)

C_d＝0时，球缝模型退化为纯球形孔，据此可以得到模型方程：

R_f＝0 (23)

R_s＝R_e (24)

V_f＝S_f＝0 (25)

核磁回波数据反演的布点计算公式表示为：

其中，Re是等效球的半径，C_d是第一定义参数，R_ei为核磁回波数据反演的第i个布点的等效球的半径，C_di是核磁回波数据反演的第i个布点的第一定义参数；ρ₂为横向弛豫率，T₂(R_ei，C_di)为根据球缝模型计算的横向弛豫时间，i是核磁回波数据反演的布点编号。

4.根据权利要求1所述的核磁回波数据优化反演方法，其特征在于，扫描确定第一定义参数的方法为：在每个球缝模型的分组内部确定所述第一定义参数，扫面确定根据该第一定义参数确定的路径图；对所述路径图形进行扫描，得到描述岩石孔隙特征的路径，每一个路径上，任何一点的横坐标就是核磁反演的布点序号，每一个布点序号对应于一个弛豫时间布点，任何一点的纵坐标就是第一定义参数值。

5.孔隙结构参数确定方法，其特征在于，采用如权利要求1～4任意一项权利要求所述的核磁回波数据优化反演方法获得核磁回波数据的优化的T₂谱分布；对于得到的最优T₂谱进行球形孔和裂缝孔的分解。

6.根据权利要求5所述的孔隙结构参数确定方法，其特征在于，对于得到的最优T₂谱进行球形孔和裂缝孔的分解的方程如下：

其中，A_i，A_si，A_fi，分别是最优谱分布幅度，球形孔谱分布幅度，裂缝孔谱分布幅度，标准化数据无量纲；S_ti，S_si，S_fi，分别是最优谱的面积，球形孔面积，裂缝孔面积。

7.根据权利要求5所述的孔隙结构参数确定方法，其特征在于，在完成最优谱的分解之后，分别得到球形孔和裂缝孔分布。

8.根据权利要求5所述的孔隙结构参数确定方法，其特征在于，分别做球形孔和裂缝谱的累计曲线，再计算不同频率处的孔隙半径值，或裂缝宽度值，并使用分选系数、半径均值或峰态参数的定义计算孔隙结构参数。

9.根据权利要求8所述的孔隙结构参数确定方法，其特征在于，定义孔结构参数包括：球形孔的分选系数S_p计算公式

球形孔的歪度为SK_p，

球形孔的均值半径为D_M

球形孔的峰态为Kp

球形孔的结构系数为F

为平均孔隙半径；

为孔隙度，K为渗透率。

Images