CN111767650B - 一种数字岩芯等效弹性参数估算方法及装置 - Google Patents
一种数字岩芯等效弹性参数估算方法及装置 Download PDFInfo
- Publication number
- CN111767650B CN111767650B CN202010595422.3A CN202010595422A CN111767650B CN 111767650 B CN111767650 B CN 111767650B CN 202010595422 A CN202010595422 A CN 202010595422A CN 111767650 B CN111767650 B CN 111767650B
- Authority
- CN
- China
- Prior art keywords
- digital core
- equivalent
- pressure
- elastic wave
- elastic
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 239000011435 rock Substances 0.000 title claims abstract description 42
- 238000004088 simulation Methods 0.000 claims abstract description 74
- 229910052500 inorganic mineral Inorganic materials 0.000 claims abstract description 51
- 239000011707 mineral Substances 0.000 claims abstract description 51
- 238000004458 analytical method Methods 0.000 claims abstract description 47
- 238000012937 correction Methods 0.000 claims abstract description 30
- 239000002245 particle Substances 0.000 claims abstract description 27
- 238000004364 calculation method Methods 0.000 claims abstract description 25
- 238000012545 processing Methods 0.000 claims abstract description 18
- 230000011218 segmentation Effects 0.000 claims abstract description 18
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 12
- 239000013078 crystal Substances 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 12
- 238000010586 diagram Methods 0.000 claims description 8
- 239000000126 substance Substances 0.000 claims description 7
- 238000003709 image segmentation Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 2
- 230000003068 static effect Effects 0.000 abstract description 15
- 238000005315 distribution function Methods 0.000 abstract description 3
- 230000006870 function Effects 0.000 description 22
- 238000012512 characterization method Methods 0.000 description 4
- 101001093143 Homo sapiens Protein transport protein Sec61 subunit gamma Proteins 0.000 description 3
- 101000694017 Homo sapiens Sodium channel protein type 5 subunit alpha Proteins 0.000 description 3
- 102100027198 Sodium channel protein type 5 subunit alpha Human genes 0.000 description 3
- 230000014509 gene expression Effects 0.000 description 3
- 239000011148 porous material Substances 0.000 description 3
- 101100120905 Saccharomyces cerevisiae (strain ATCC 204508 / S288c) TDH1 gene Proteins 0.000 description 2
- 230000000295 complement effect Effects 0.000 description 2
- 239000007771 core particle Substances 0.000 description 2
- 238000002050 diffraction method Methods 0.000 description 2
- 241001270131 Agaricus moelleri Species 0.000 description 1
- 229910021532 Calcite Inorganic materials 0.000 description 1
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 235000015076 Shorea robusta Nutrition 0.000 description 1
- 244000166071 Shorea robusta Species 0.000 description 1
- TZCXTZWJZNENPQ-UHFFFAOYSA-L barium sulfate Chemical compound [Ba+2].[O-]S([O-])(=O)=O TZCXTZWJZNENPQ-UHFFFAOYSA-L 0.000 description 1
- 229910052601 baryte Inorganic materials 0.000 description 1
- 239000010428 baryte Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000004927 clay Substances 0.000 description 1
- 239000002734 clay mineral Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000010433 feldspar Substances 0.000 description 1
- 238000011545 laboratory measurement Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000011028 pyrite Substances 0.000 description 1
- 229910052683 pyrite Inorganic materials 0.000 description 1
- NIFIFKQPDTWWGU-UHFFFAOYSA-N pyrite Chemical compound [Fe+2].[S-][S-] NIFIFKQPDTWWGU-UHFFFAOYSA-N 0.000 description 1
- 239000010453 quartz Substances 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 229910021646 siderite Inorganic materials 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N silicon dioxide Inorganic materials O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
- G06T5/92—Dynamic range modification of images or parts thereof based on global image properties
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Evolutionary Computation (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Operations Research (AREA)
- Evolutionary Biology (AREA)
- Artificial Intelligence (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Computing Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种数字岩芯等效弹性参数估算方法及装置，方法包括：从扫描装置中获得数字岩芯图像，并对数字岩芯图像进行阈值分割得到多个矿物类型，从而构建无压力数字岩芯模型，对其进行取向分布函数分析及矿物组分等效模量计算得到等效模拟的参数；然后基于颗粒边缘与压力的关系采用分水岭算法对数字岩芯模型处理得到多个有压力数字岩芯模型；根据各向同性弹性波一阶速度‑应变差分方程分别在多个压力数字岩芯模型的模拟传播得到等效速度平均值；对等效速度平均值的校正分析得到校正结果。本发明有效地改善了弹性参数估算精度，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
Description
技术领域
本发明主要涉及岩石弹性参数估算技术领域，具体涉及一种数字岩芯等效弹性参数估算方法及装置。
背景技术
数字岩芯是估算岩石弹性参数的一类常用方法。岩石的弹性参数对认识和了解其性质具有重要的意义。目前基于数字岩芯的弹性参数等效计算方法常见的可分为基于岩石物理等效表征和静态模量数值模拟。常见的岩石物理弹性等效解析模型主要有(MavkoG.2009)：微分等效介质模型(DEM)、自相容近似模型(SCA)等。邓继新等(2015)基于龙马溪组页岩数字岩芯，对比研究了DEM、SCA、Backus平均三种模型的弹性等效建模精度。赵峦啸等(2016)综合利用Xu&White、SCA、DEM等模型，提出了表征不同成熟度页岩的弹性等效参数。这些岩石物理解析等效弹性参数估算因假设前提条件的理论局限性，极大限制了其使用范围，且精度有限。静态模量数值模拟方法是近十年发展起来的，例如Arns等(2002)采用有限元法模拟计算了不同孔隙度下Fontainebleau砂岩数字岩芯的弹性参数，并与实验测量结果进行对比验证。Knackstedt等(2009)讨论了复杂碎屑岩和碳酸盐岩数字岩芯的弹性等效数值建模计算问题。孙建孟等(2014)利用扩展有限元法求取了低渗透砂岩数字岩芯的弹性模量。张文辉等(2016)基于龙马溪组页岩数字岩芯实现弹性等效有限元数值模拟，讨论了不同孔隙度和干酪根对岩石弹性模量的影响。静态法虽然适用性广，但是由于线性等效计算模型的局限性，对于复杂介质弹性等效参数的计算精度有限，很难精确模拟。另外，静态法对数字岩芯分辨率有很高要求，且模拟估算过程效率不高。
发明内容
本发明所要解决的技术问题是针对现有技术的不足，提供一种数字岩芯等效弹性参数估算方法及装置。
本发明解决上述技术问题的技术方案如下：一种数字岩芯等效弹性参数估算方法，包括如下步骤：
从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
对所述等效速度平均值进行校正分析，得到校正结果。
本发明解决上述技术问题的另一技术方案如下：一种数字岩芯等效弹性参数估算装置，包括：
数字岩芯图像分割模块，用于从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
弹性波分析模块，用于对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
有效压力模型建立模块，用于基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
压力数字岩芯模型模拟模块，用于根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
校正分析模块，用于对所述等效速度平均值进行校正分析，得到校正结果。
本发明的有益效果是：通过对数字岩芯图像的阈值分割分析得到多个矿物类型，并根据多个矿物类型构建无压力数字岩芯模型，对无压力数字岩芯模型的弹性波模拟分析得到各向同性弹性波一阶速度-应变差分方程，基于颗粒边缘与压力的关系采用分水岭算法对无压力数字岩芯模型处理得到多个有压力数字岩芯模型，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，根据各向同性弹性波一阶速度-应变差分方程分别在多个有压力数字岩芯模型的模拟传播得到等效速度平均值，对等效速度平均值的校正分析得到校正结果，有效地改善了弹性参数估算精度，且推广性和适用性强对等效速度平均值的校正分析得到校正结果，提供了一种适用于强非均质岩芯模型、优势互补的等效弹性参数估算方法，为认识和了解岩石性质提供了可靠依据。
附图说明
图1为本发明实施例提供的数字岩芯等效弹性参数估算方法的流程示意图；
图2为本发明实施例提供的400×400页岩数字岩芯原始数据实验图；
图3为本发明实施例提供的阈值分割后的数字岩芯图；
图4为本发明实施例提供的取向分布函数坐标转换Euler空间示意图；
图5为本发明实施例提供的各向同性弹性波方程模拟时间0.07us的波场快照图；
图6为本发明实施例提供的各向同性弹性波方程模拟时间0.1us的波场快照图；
图7为本发明实施例提供的数字岩芯颗粒接触表征过程的阈值分割图；
图8为本发明实施例提供的数字岩芯颗粒接触表征过程的二值化图；
图9为本发明实施例提供的数字岩芯颗粒接触表征过程的岩石颗粒边缘结构图；
图10为本发明实施例提供的速度随压力变化的数值模拟结果与实验结果的对比图；
图11为本发明实施例提供的根据校正策略得到的结果对比图；
图12为本发明实施例提供的数字岩芯等效弹性参数估算装置的模块框图。
具体实施方式
以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
图1为本发明实施例提供的数字岩芯等效弹性参数估算方法的流程示意图。
如图1所示，一种数字岩芯等效弹性参数估算方法，包括如下步骤：
S1：从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
S2：对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
S3：基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
S4：根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
S5：对所述等效速度平均值进行校正分析，得到校正结果。
上述实施例中，通过对数字岩芯图像的阈值分割分析得到多个矿物类型，并根据多个矿物类型构建无压力数字岩芯模型，对无压力数字岩芯模型的弹性波模拟分析得到各向同性弹性波一阶速度-应变差分方程，基于颗粒边缘与压力的关系采用分水岭算法对无压力数字岩芯模型处理得到多个有压力数字岩芯模型，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，根据各向同性弹性波一阶速度-应变差分方程分别在多个有压力数字岩芯模型的模拟传播得到等效速度平均值，对等效速度平均值的校正分析得到校正结果，有效地改善了弹性参数估算精度，且推广性和适用性强对等效速度平均值的校正分析得到校正结果，提供了一种适用于强非均质岩芯模型、优势互补的等效弹性参数估算方法，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述得到多个矿物类型的过程包括：
对所述数字岩芯图像进行灰度图像处理，得到灰度图像，所述灰度图像包括多个像素点；
对各个所述像素点进行概率分布计算，得到多个相对概率分布直方图；
利用阈值分割法对多个所述相对概率分布直方图进行阈值分割，得到多个最优阈值；
根据所述多个最优阈值对所述数字岩芯图像进行矿物分类，得到多个矿物类型。
具体地，如图2至3所示，图2中，不同灰度代表不同的矿物组分，阈值分割法具体步骤如下：
SS1、所述数字岩芯图像可以看作是一个个像素组成，选取其分辨率为400×400像素，大小为400×400μm2。
SS2、假设图像灰度级为S，图像中灰度级为i的像素点有mi(i＝0,1,2,…,S)个，则总像素有
某一像素呈现的概率为Pi＝mi/M，通过计算各像素点灰度的概率分布，获得相对概率分布直方图，图中出现一些峰值及拐点，根据这些点可以将原始灰度图像大致划分成N(N>3)类，需要用N–1个分割阈值，设阈值集合t＝{tkk＝1,2,3,...,N-1}，同时令t0＝0,tN＝S，并通过第五式得到多阈值的类内方差
所述第五式为：
其中，
其中，ωk为第k类的概率，
为第k类的类方差。
SS3、当类内方差最小，通过第六式得到最优阈值，所述第六式为：
其中，
分别为多个最优阈值，下标0＜t1＜t2＜...＜tN-1为初始阈值，即根据步骤SS2、SS3获取的。
SS4、根据上述步骤SS1、SS2、SS3的描述，根据不同矿物的X射线吸收率，令N＝6将所述数字岩芯图像分割为六组矿物类型：孔隙；TOC；粘土；石英；方解石、长石、重晶石等；黄铁矿、菱铁矿等。
上述实施例中，对所述数字岩芯图像的灰度图像处理得到灰度图像，所述灰度图像包括多个像素点；对各个所述像素点的概率分布计算得到多个相对概率分布直方图；利用阈值分割法对多个所述相对概率分布直方图的阈值分割得到多个最优阈值；根据所述多个最优阈值对所述数字岩芯图像的矿物分类得到多个矿物类型，为后续处理提供模型，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述得到各向同性弹性波一阶速度-应变差分方程的过程包括：
对所述无压力数字岩芯模型进行取向分布分析，得到取向密度函数；
利用矿物组分等效模量法对所述取向密度函数进行弹性波模拟参数分析，得到等效弹性模量；
将所述等效弹性模量作为各向同性弹性波模拟参数，对所述各向同性弹性波模拟参数进行弹性波模拟计算，得到各向同性弹性波一阶速度-应变差分方程。
上述实施例中，对所述无压力数字岩芯模型的取向分布分析得到取向密度函数；利用矿物组分等效模量法对所述取向密度函数的弹性波模拟参数分析得到等效弹性模量；将所述等效弹性模量作为各向同性弹性波模拟参数，对所述各向同性弹性波模拟参数的弹性波模拟计算得到各向同性弹性波一阶速度-应变差分方程，为后续处理提供基础，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述对所述数字岩芯模型进行取向分布分析，得到取向密度函数的过程包括：
建立一个以
为坐标系的Euler空间；
通过第一式对所述无压力数字岩芯模型进行取向分布计算，得到取向密度函数，所述第一式为：
其中，
应理解地，利用晶体学中的取向分布函数(ODF)来分析所述无压力数字岩芯模型中矿物颗粒取向分布方位特征，从而证明其各向同性性质，并得到取向密度函数。
具体地，依据晶体学基本原理，不同的矿物具有不同的优势排列方向。特别是富含黏土矿物的页岩，矿物颗粒、层理及微裂缝存在明显的分布方向特征，因此进行矿物颗粒特征分析。如图4所示，首先通过任意三次旋转来建立一个以
为坐标系的Euler空间，常选择取向密度函数
来表达多晶集合体的取向分布情况。定义处于取向
周围的取向单元
中的晶粒的体积百分比如第七式，所述第七式为：
其中，dV为取向落在该取向元的晶粒体积，V为样品的体积，Kw为比例系数，常取值为1。
由于ODF无法直接测量，只得通过傅里叶分析中的级数展开原理，用调和函数分析法，利用第八式进行计算，所述第八式为：
其中，Wlmn为级数的第lmn项系数，l为级数展开项。同时根据第九式、第十式和第十一式求取Wlmn之后，反代入所述第八式中，从而计算出
所述第九式、所述第十式和所述第十一式分别为：
其中，qi(α,β)为晶面i的极图函数，α为晶面法线的倾角，β为晶面法线的旋转角，ei和e–i为单独一项自然对数函数，i2＝-1，
与Φi分别为晶面i的法线在晶体坐标系的倾角及辐角，Zlmn(cosθ)为归一化Jacobi多项式，Wlmn为级数的第lmn项系数，Pl m(cosα)为归一化的连带Legendre多项式，
为级数的第lm项系数。
上述实施例中，通过第一式对所述无压力数字岩芯模型的取向分布计算得到取向密度函数，为后续参数处理提供数据支持，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述利用矿物组分等效模量法对所述取向密度函数进行弹性波模拟参数分析，得到等效弹性模量的过程包括：
导入参数信息，所述参数信息从岩石物理手册中获得，并依据所述取向密度函数在导入的参数信息中筛选矿物的基本弹性参数，得到多个基本弹性参数；
利用Hill平均法对多个所述基本弹性参数进行等效值计算，得到等效弹性模量。
应理解地，将所述取向密度函数作为参考，选择利用矿物组分等效模量法得到各类矿物组分的所述基本弹性参数，从而获取所述等效弹性模量。
具体地，在求取所述取向密度函数的基础上，将其作为参考依次选择岩石物理手册(Mavko等，2009)可知各种矿物的所述基本弹性参数(如拉梅系数、密度等)，通过数字岩芯阈值分割数据体的矿物含量分析，得到该样本各矿物组分的百分含量及其基本弹性参数，同时利用Hill平均法计算得到阈值分割后六种矿物组分类型的等效弹性模量。表1为该岩芯样本各矿物组分的百分含量及其基本弹性参数。
表1
上述实施例中，导入参数信息，所述参数信息从岩石物理手册中获得，并依据所述取向密度函数在导入的参数信息中筛选矿物的基本弹性参数，得到多个基本弹性参数；利用Hill平均法对多个所述基本弹性参数的等效值计算得到等效弹性模量，为后续计算提供数据支持，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述对所述各向同性弹性波模拟参数进行弹性波计算，得到各向同性弹性波一阶速度-应变差分方程的过程包括：
所述各向同性弹性波模拟参数包括介质密度、拉梅系数和剪切模量；
通过各向同性弹性波差分方程组对所述介质密度、所述拉梅系数和剪切模量进行弹性波模拟计算，得到各向同性弹性波一阶速度-应变差分方程，所述各向同性弹性波差分方程组为：
其中，ρ为介质密度，λ为拉梅系数，μ为剪切模量，U为X方向速度的离散值，V为Z方向速度的离散值，R为XX方向应力的离散值，T为ZZ方向应力的离散值，H为XZ方向应力的离散值，k为时间节点，i为X方向网格节点，j为Z方向网格节点，Cn为差分系数,Δt为时间差分增量，Δx为空间差分增量，Δz为空间差分增量。
具体地，如图5至6所示，本发明利用高精度的交错网格有限差分法求解各向同性弹性波方程来模拟超声波，各向同性弹性波一阶速度-应变差分结构如下：
二维情况下(XOZ平面)，弹性波方程为第一方程组，所述第一方程组为：
其中，ρ为介质密度，λ为拉梅系数，μ为剪切模量，Vx,Vz,σxx,σzz,σxz分别为X方向速度、Z方向速度、XX方向应力、ZZ方向应力、XZ方向应力。
将高阶差分近似代入到所述第一方程组中得到第二方程组，就得到所述各向同性弹性波一阶速度-应变差分方程，所述第二方程组为：
其中，U,V,R,T,H分别表示Vx,Vz,σxx,σzz,σxz的离散值，ρ为介质密度，λ为拉梅系数，μ为剪切模量，k为时间节点，i为X方向网格节点，j为Z方向网格节点，Cn为差分系数,Δt为时间差分增量，Δx为空间差分增量，Δz为空间差分增量。
表2为各向同性弹性波方程模拟参数表。
表2
上述实施例中，通过各向同性弹性波差分方程组对所述介质密度、所述拉梅系数和剪切模量的弹性波模拟计算得到各向同性弹性波一阶速度-应变差分方程，为后续处理提供超声波，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述得到多个有压力数字岩芯模型的过程包括：
利用分水岭算法对所述无压力数字岩芯模型进行颗粒边缘的接触关系处理，得到多个颗粒间接触关系；
通过第二式分别对所述多个颗粒间接触关系进行计算，得到多个有效压力，所述第二式为：
Pf＝c+wMgtg，
其中，Mgtg∈[0,1]为颗粒间接触关系，Pf为有效压力，c为预设的拟合系数，w为预设的拟合系数；
根据多个所述有效压力和所述无压力数字岩芯模型进行构建，得到多个有压力数字岩芯模型。
优选地，所述预设的拟合系数c＝0.9743Mpa，所述预设的拟合系数w＝200.86Mpa。
应理解地，为后续计算方便，本发明在得到的多个有压力数字岩芯模型上下各加入一层厚度为100像素的均匀介质，大小为600×400μm2，即将模型网格大小定为600×400。
具体地，如图7至9所示，得到多个有压力数字岩芯模型其步骤如下：
SSS1:在所述无压力数字岩芯模型的基础上，获得二值化图像；
SSS2:在SSS1的基础上采用倒角距离变换，对原始灰度图像进行梯度计算，按照每个像素的灰度级进行从低到高排序，同时对每个局部极小值的影响域慢慢向外扩展至边缘进行标注，获取数字岩芯中的岩石颗粒(含孔隙)之间的接触关系图；
SSS3:颗粒间边缘变化情况与压力之间存在第十五式的关系，所述第十五式为：
Pf＝c+wMgtg，
其中，Mgtg∈[0,1]为颗粒间接触关系，Pf为有效压力，c为预设的拟合系数，w为预设的拟合系数，前人研究确定c＝0.9743MPa和w＝200.86MPa。因此根据该关系，改变颗粒接触关系Mgtg，重复SSS1，SSS2，构建0～50MPa下的有压力数字岩芯模型。
上述实施例中，利用分水岭算法对所述无压力数字岩芯模型的颗粒边缘的接触关系处理得到多个颗粒间接触关系；通过第二式分别对所述多个颗粒间接触关系的计算得到多个有效压力，根据多个所述有效压力和所述无压力数字岩芯模型的构建得到多个有压力数字岩芯模型，为之后的计算提供数据基础，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述得到等效速度平均值的过程包括：
通过第三式分别对所述多个有压力数字岩芯模型进行计算，得到多个有压力数字岩芯模型等效速度，所述第三式为：
其中，
其中，
其中，Vp2为有压力数字岩芯模型等效速度，H2为有压力数字岩芯模型高度，Δt：为振幅峰值时间差，Vp1为预设无压力数字岩芯模型等效速度，t1为参考模型时间，t2：为有压力数字岩芯模型时间，H1：为预设无压力数字岩芯模型高度；
对多个所述有压力数字岩芯模型等效速度进行平均值计算，得到等效速度平均值。
具体地，为了质控不同压力下的数字岩芯模拟的条件，采用相同的初始参数、相同的震源位置和接收条件进行模拟。在模型顶部激发震源，底部放置检波器，模拟过程中不同于常规的利用均匀介质作为背景参考计算岩石等效速度，本发明采用无压力时的数字岩芯作为背景参考来计算压力对等效速度的影响。
应理解地，假设所述第三式中Vp2为待求结果，H1＝H2为激发点到接收点的距离，先求出参考模型时间t1，再分别利用多个有压力数字岩芯模型与0MPa数字岩芯振幅峰值时间差来计算Vp2；获取数字岩芯模型底部第检波器的接收结果，计算每道结果Vn p2，取每个检波器计算结果的平均值Vnum作为最后结果。
上述实施例中，通过第三式分别对所述多个有压力数字岩芯模型的计算得到多个有压力数字岩芯模型等效速度，并对多个所述有压力数字岩芯模型等效速度的平均值计算得到等效速度平均值，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
可选地，作为本发明的一个实施例，所述对所述等效速度平均值进行校正分析，得到校正结果的过程包括：
通过第四式对所述等效速度平均值进行校正计算，得到关系因子s和关系因子g，所述第四式为：
Vexp＝sVnum+g，
其中，Vexp为预设实验室测量速度，Vnum为等效速度平均值；
利用最小二乘法对所述关系因子s和所述关系因子g进行拟合，得到校正结果。
应理解地，图10至11中，圆圈为数值模拟结果，方框为实验结果，所述等效速度平均值与超声实验二者保持较好的相对变化关系，但绝对值有较大差别。原因可能是数字岩芯的分辨能力有限，大量的微孔、裂纹、节理和颗粒接触特征等细节未能刻画并参与模拟计算。另外，页岩矿物组分复杂，以六类矿物分类较粗，没有考虑到所有的矿物。鉴于此情况，本发明提出以实验室测量速度为准对所述等效速度平均值进行校正，并用最小二乘法拟合得到：s＝0.4486，g＝2385.3。由此可见，在数字岩芯弹性参数模拟中对其可进行校准以获得最佳结果。
上述实施例中，通过第四式对所述等效速度平均值的校正计算得到关系因子和关系因子，并利用最小二乘法对所述关系因子和所述关系因子的拟合得到校正结果，有效地改善了弹性参数估算精度，且推广性和适用性强，有效地解决常规岩石物理等效弹性模拟及静态数值模拟的局限性问题，为认识和了解岩石性质提供了可靠依据。
图12为本发明实施例提供的数字岩芯等效弹性参数估算装置的模块框图。
可选地，作为本发明的另一个实施例，如图12所示，一种数字岩芯等效弹性参数估算装置，包括：
数字岩芯图像分割模块，用于从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
弹性波分析模块，用于对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
有效压力模型建立模块，用于基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
压力数字岩芯模型模拟模块，用于根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
校正分析模块，用于对所述等效速度平均值进行校正分析，得到校正结果。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。
作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本发明实施例方案的目的。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以是两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。用于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。
Claims (10)
1.一种数字岩芯等效弹性参数估算方法，其特征在于，包括如下步骤：
从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
对所述等效速度平均值进行校正分析，得到校正结果。
2.根据权利要求1所述的数字岩芯等效弹性参数估算方法，其特征在于，所述得到多个矿物类型的过程包括：
对所述数字岩芯图像进行灰度图像处理，得到灰度图像，所述灰度图像包括多个像素点；
对各个所述像素点进行概率分布计算，得到多个相对概率分布直方图；
利用阈值分割法对多个所述相对概率分布直方图进行阈值分割，得到多个最优阈值；
根据所述多个最优阈值对所述数字岩芯图像进行矿物分类，得到多个矿物类型。
3.根据权利要求1所述的数字岩芯等效弹性参数估算方法，其特征在于，所述得到各向同性弹性波一阶速度-应变差分方程的过程包括：
对所述无压力数字岩芯模型进行取向分布分析，得到取向密度函数；
利用矿物组分等效模量法对所述取向密度函数进行弹性波模拟参数分析，得到等效弹性模量；
将所述等效弹性模量作为各向同性弹性波模拟参数，对所述各向同性弹性波模拟参数进行弹性波模拟计算，得到各向同性弹性波一阶速度-应变差分方程。
4.根据权利要求3所述的岩芯等效弹性参数估算方法，其特征在于，所述得到取向密度函数的过程包括：
建立一个以
为坐标系的Euler空间；
通过第一式对所述无压力数字岩芯模型进行取向分布计算，得到取向密度函数，所述第一式为：
其中，级数的第lmn项系数Wlmn是通过第十一式计算得到，并将所述级数的第lmn项系数Wlmn代入所述第一式中，所述第十一式为：
其中，级数的第lm项系数
是通过第十式计算得到，并将所述级数的第lm项系数
代入所述第十一式中，所述第十式为：
5.根据权利要求3所述的数字岩芯等效弹性参数估算方法，其特征在于，所述利用矿物组分等效模量法对所述取向密度函数进行弹性波模拟参数分析，得到等效弹性模量的过程包括：
导入参数信息，所述参数信息从岩石物理手册中获得，并依据所述取向密度函数在导入的参数信息中筛选矿物的基本弹性参数，得到多个基本弹性参数；
利用Hill平均法对多个所述基本弹性参数进行等效值计算，得到等效弹性模量。
6.根据权利要求3或5所述的岩芯等效弹性参数估算方法，其特征在于，所述得到各向同性弹性波一阶速度-应变差分方程的过程包括：
所述各向同性弹性波模拟参数包括介质密度、拉梅系数和剪切模量；
通过各向同性弹性波差分方程组对所述介质密度、所述拉梅系数和剪切模量进行弹性波模拟计算，得到各向同性弹性波一阶速度-应变差分方程，所述各向同性弹性波差分方程组为：
其中，ρ为介质密度，λ为拉梅系数，μ为剪切模量，U为X方向速度的离散值，V为Z方向速度的离散值，R为XX方向应力的离散值，T为ZZ方向应力的离散值，H为XZ方向应力的离散值，k为时间节点，i为X方向网格节点，j为Z方向网格节点，Cn为差分系数,Δt为时间差分增量，Δx为空间差分增量，Δz为空间差分增量。
7.根据权利要求1所述的数字岩芯等效弹性参数估算方法，其特征在于，所述得到多个有压力数字岩芯模型的过程包括：
利用分水岭算法对所述无压力数字岩芯模型进行颗粒边缘的接触关系处理，得到多个颗粒间接触关系；
通过第二式分别对所述多个颗粒间接触关系进行计算，得到多个有效压力，所述第二式为：
Pf＝c+wMgtg，
其中，Mgtg∈[0,1]为颗粒间接触关系，Pf为有效压力，c为预设的拟合系数，w为预设的拟合系数；
根据多个所述有效压力和所述无压力数字岩芯模型进行构建，得到多个有压力数字岩芯模型。
8.根据权利要求1所述的数字岩芯等效弹性参数估算方法，其特征在于，所述得到等效速度平均值的过程包括：
通过第三式分别对所述多个有压力数字岩芯模型进行计算，得到多个有压力数字岩芯模型等效速度，所述第三式为：
其中，
其中，
其中，Vp2为有压力数字岩芯模型等效速度，H2为有压力数字岩芯模型高度，Δt为振幅峰值时间差，Vp1为预设无压力数字岩芯模型等效速度，t1为参考模型时间，t2为有压力数字岩芯模型时间，H1为预设无压力数字岩芯模型高度；
对多个所述有压力数字岩芯模型等效速度进行平均值计算，得到等效速度平均值。
9.根据权利要求8所述的数字岩芯等效弹性参数估算方法，其特征在于，所述对所述等效速度平均值进行校正分析，得到校正结果的过程包括：
通过第四式对所述等效速度平均值进行校正计算，得到关系因子s和关系因子g，所述第四式为：
Vexp＝sVnum+g，
其中，Vexp为预设实验室测量速度，Vnum为等效速度平均值；
利用最小二乘法对所述关系因子s和所述关系因子g进行拟合，得到校正结果。
10.一种数字岩芯等效弹性参数估算装置，其特征在于，包括：
数字岩芯图像分割模块，用于从扫描装置中获得数字岩芯图像，并对所述数字岩芯图像进行阈值分割分析，得到多个矿物类型，并根据多个所述矿物类型构建无压力数字岩芯模型；
弹性波分析模块，用于对所述无压力数字岩芯模型进行弹性波模拟分析，得到各向同性弹性波一阶速度-应变差分方程；
有效压力模型建立模块，用于基于颗粒边缘与压力的关系采用分水岭算法对所述无压力数字岩芯模型处理得到多个有压力数字岩芯模型；
压力数字岩芯模型模拟模块，用于根据所述各向同性弹性波一阶速度-应变差分方程分别在所述多个有压力数字岩芯模型进行模拟传播，得到等效速度平均值；
校正分析模块，用于对所述等效速度平均值进行校正分析，得到校正结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010595422.3A CN111767650B (zh) | 2020-06-28 | 2020-06-28 | 一种数字岩芯等效弹性参数估算方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010595422.3A CN111767650B (zh) | 2020-06-28 | 2020-06-28 | 一种数字岩芯等效弹性参数估算方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111767650A CN111767650A (zh) | 2020-10-13 |
| CN111767650B true CN111767650B (zh) | 2022-05-03 |
Family
ID=72722370
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010595422.3A Active CN111767650B (zh) | 2020-06-28 | 2020-06-28 | 一种数字岩芯等效弹性参数估算方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111767650B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112394072B (zh) * | 2020-11-26 | 2021-10-22 | 西安石油大学 | 一种基于Micro-CT岩芯宽频介电常数表征方法和装置 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105095631A (zh) * | 2014-05-21 | 2015-11-25 | 中国石油化工股份有限公司 | 一种页岩各向异性岩石物理建模方法 |
| CN107045580A (zh) * | 2017-04-27 | 2017-08-15 | 中国石油大学(华东) | 一种基于数字岩心的页岩力学参数快速计算方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105393110B (zh) * | 2013-08-06 | 2019-03-29 | Bp北美公司 | 在模拟应力和应变条件下的岩石物理属性的基于图像的直接数值模拟 |
| CN109581492A (zh) * | 2017-09-29 | 2019-04-05 | 中国石油化工股份有限公司 | 基于地震波模拟的岩石物理参数计算方法及系统 |
-
2020
- 2020-06-28 CN CN202010595422.3A patent/CN111767650B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105095631A (zh) * | 2014-05-21 | 2015-11-25 | 中国石油化工股份有限公司 | 一种页岩各向异性岩石物理建模方法 |
| CN107045580A (zh) * | 2017-04-27 | 2017-08-15 | 中国石油大学(华东) | 一种基于数字岩心的页岩力学参数快速计算方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111767650A (zh) | 2020-10-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Costanza‐Robinson et al. | Representative elementary volume estimation for porosity, moisture saturation, and air‐water interfacial areas in unsaturated porous media: Data quality implications | |
| Saxena et al. | Estimating elastic moduli of rocks from thin sections: Digital rock study of 3D properties from 2D images | |
| CA2700666C (en) | Method for determining the properties of hydrocarbon reservoirs from geophysical data | |
| US9070049B2 (en) | Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties | |
| US10223782B2 (en) | Digital rock physics-based trend determination and usage for upscaling | |
| Buscombe et al. | Grain‐size information from the statistical properties of digital images of sediment | |
| Jouini et al. | Numerical estimation of carbonate rock properties using multiscale images | |
| Zeghadi et al. | Ensemble averaging stress–strain fields in polycrystalline aggregates with a constrained surface microstructure–Part 1: Anisotropic elastic behaviour | |
| Karimpouli et al. | Estimating 3D elastic moduli of rock from 2D thin-section images using differential effective medium theory | |
| Smal et al. | An automatic segmentation algorithm for retrieving sub-resolution porosity from X-ray tomography images | |
| CN106525680B (zh) | 岩心孔隙度参数场的获取方法 | |
| US20140270393A1 (en) | Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties | |
| EP2090907A1 (en) | Method for determining the properties of hydrocarbon reservoirs from geophysical data | |
| Callow et al. | Optimal X-ray micro-CT image based methods for porosity and permeability quantification in heterogeneous sandstones | |
| Claes et al. | A three-dimensional classification for mathematical pore shape description in complex carbonate reservoir rocks | |
| Taguas et al. | Simulation and testing of self-similar structures for soil particle-size distributions using iterated function systems | |
| WO2014142976A1 (en) | Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties | |
| CN111767650B (zh) | 一种数字岩芯等效弹性参数估算方法及装置 | |
| Scott et al. | Multi-scale image-based pore space characterisation and pore network generation: Case study of a north sea sandstone reservoir | |
| CN111426616B (zh) | 碳酸盐岩弹性性质与孔隙结构获取方法、装置及存储介质 | |
| AU2014290779A1 (en) | System and method for estimating porosity distribution in subterranean reservoirs | |
| Barker et al. | A parametrization of 3‐D subgrid‐scale clouds for conventional GCMs: Assessment using A‐Train satellite data and solar radiative transfer characteristics | |
| Sharma et al. | Sensitivity of digital rock method for pore-space estimation to heterogeneity in carbonate formations | |
| Datta et al. | Determination of porosity of rock samples from photomicrographs using image analysis | |
| Naraghi et al. | 3D reconstruction of porous media from a 2D section and comparisons of transport and elastic properties |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | 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 |