一种基于动态基材料选择的非均质岩心矿物表征方法

1.本发明涉及岩石物理及数字岩心领域，特别涉及一种非均质岩心的数字岩心建模中的矿物表征方法。


背景技术：

2.数字岩石物理(digital rock physics，drp)是石油地质分析的重要手段，包括数字岩心建模和数值模拟计算两个主要步骤。x射线计算机断层成像(x-ray computed tomography，xct)技术因其无损、成像效果好、可重复性高等优点，成为数字岩心建模的热门手段。xct技术依靠线性衰减系数差异来揭示岩心的矿物组分和内部结构。由于线性衰减系数与物质的密度直接相关，因此可以轻易地将岩心分离为骨架和孔隙两相介质，实现岩心的数字化。然而，在非均质岩心中，由于某些矿物的线性衰减系数值非常接近，仅依靠ct图像难以对骨架中不同矿物进行细分，此时的岩心数字模型无法准确表征岩心的性质。因此，从数字岩心建模的角度出发，对非均质岩心进行多矿物组分的无损表征尤其重要。
3.双能ct技术将两个不同能量下测得的线性衰减系数分解为电子密度和有效原子序数，具有更强大的物质表征能力。目前，在利用双能ct进行岩心矿物表征的探索中，alves和martini分别利用基效应分解法和化学计量法完成了岩心的矿物表征。然而，上述两种方法的准确性很大程度上取决于用于拟合的标准材料的数量和纯度，且在致密样品中表现不好。图像域基材料分解(image-domain basis material decomposition，ibmd)方法是双能ct技术中一种典型的图像域材料分解方法，其分解精度依赖于基材料的选择，无需大量标准材料的拟合。在医学应用中，基材料的材质往往通过成分分析方法来确定，选取与样品成分相近的物质作为基材料以提高材料分解的准确度。然而，在非均质岩心表征应用中，由于矿物分布随机、成分未知且种类繁多，成分分析法难以应用，ibmd方法中基材料选择成为难题。此外，对于致密矿物，ibmd方法还存在射束硬化以及电子密度和有效原子序数计算误差大的问题，难以给出高质量的矿物表征效果。上述问题限制了ibmd方法在非均质岩心的矿物表征中的应用。


技术实现要素：

4.本发明的目的是提供一种基于动态基材料选择的非均质岩心矿物表征方法，以至少部分地解决以上所提出的技术问题。技术方案如下：
5.一种基于动态基材料选择的非均质岩心矿物表征方法，包括如下步骤：
6.s1：利用岩心先验信息建立基材料库；
7.s2：对基材料模体和岩心样品进行双能ct成像，获取基材料库能量图点云数据和样品能量图数据，随后进行动态基材料选择，为岩心样品逐像素选取最佳基材料材质和直径；
8.s3：根据基材料选择结果利用基于最佳虚拟单色图的电子密度和有效原子序数计算方法计算样品电子密度和有效原子序数；
9.s4：依据电子密度和有效原子序数信息通过阈值分割完成矿物表征。
10.进一步，上述步骤s1中所述用于岩心矿物表征的基材料库由以下5种标准单质材料组成：石墨、硅、铝、钛、铁，将每种标准单质材料制作成一组直径不同的圆柱形模体。
11.进一步，上述步骤s2具体包含以下子步骤：
12.s21：根据下式计算基材料能量图点云库f和基材料能量图点云质心库o：
[0013][0014]
式中，e代表双能ct成像中的低能能量l或高能能量h，表示能量e下直径为t的基材料k(k＝i，j且i，j∈n
*
；i，j≤5；i≠j)的线性衰减系数，均为常数，通过拟合不同直径的基材料模体重建线性衰减系数获得，为能量e下直径为t的基材料i与基材料j的线性衰减系数线性组合结果，简写为a1和a2为线性组合系数；
[0015]
为获得基材料能量图点云库f和基材料能量图点云质心库o，设定基材料直径t和线性组合系数a1，直径t取值如下：t＝t
min
，t
min
+δt，t
min
+2δt，...，t
max
，a1取值如下：a1＝0，δa1，2δa1，...，1；随后以低能衰减值μ
l
为x轴，高能衰减值μh为y轴，建立能量图；则基材料能量图点云库基材料能量图点云质心库o＝{(x
ij
，y
ij
)|i，j∈n
*
；i，j≤5；i≠j}，其中，(x
ij
，y
ij
)为基材料i与基材料j的能量图点云f
ij
的质心；
[0016]
s22：动态基材料材质选择：根据样品双能ct重建结果，在基材料能量图点云质心库o中逐像素搜索与样品能量图点欧氏距离最近的质心o
ij
，并将o
ij
对应的一组基材料材质确定为样品进行材料分解的最佳基材料材质(j
opt
，j
opt
)；
[0017]
s23：动态基材料直径选择：根据样品双能ct重建结果，在最佳基材料材质(i
opt
，j
opt
)的能量图点云中逐像素搜索与样品能量图点欧氏距离最近的点并将对应的基材料直径确定为样品进行材料分解的最佳基材料直径t
opt
。
[0018]
进一步，上述步骤s3所述电子密度和有效原子序数计算，具体包含以下步骤：
[0019]
s31：根据动态基材料选择结果，联立下式求解基材料分解系数：
[0020][0021]
式中，μ
l
和μh分别表示低能l和高能h下样品的重建线性衰减系数，分别表示低能l和高能h下样品的重建线性衰减系数，分别表示低能l和高能h下直径为t
opt
的基材料i
opt
和基材料j
opt
的线性衰减系数，由式(1)中基材料重建线性衰减系数及直径关系公式的拟合结果给出，b1和b2即为所求材料分解系数；
[0022]
s32：按照下式完成基于最佳虚拟单色图的电子密度和有效原子序数计算：
[0023][0024]
式中，和为最佳虚拟单色能量e1和e2下的样品线性衰减系数，根据基材料分解模型获得，和分别为基材料i
opt
和基材料j
opt
在能量e下的理论线性衰
减系数，由nist xcom数据库根据基材料成分及密度给出，b1和b2分别为两种基材料的分解系数，ρe和z
eff
分别为所求电子密度和有效原子序数，m(e1)、k(e2，e1)、g(e2，e1)和n均为常数，通过拟合给出；需要注意的是，上述拟合过程一般在准备阶段完成，拟合数据由nistxcom数据库根据常见成岩矿物的化学式和密度给出，无需进行真实实验。
[0025]
进一步，所述步骤s32中最佳虚拟单色能量e1和e2分别为1mev和10kev，通过分析成岩矿物的衰减特性确定。
[0026]
进一步，所述步骤s4具体包含以下步骤：
[0027]
s41：依据岩心样品的先验信息，获得该类岩心中常见矿物的化学式并计算该类岩心中常见矿物的电子密度和有效原子序数真值，将不同矿物之间的电子密度阈值和有效原子序数阈值设定为真值的
±
5％；
[0028]
s42：结合电子密度和有效原子序数计算结果进行阈值分割，将分割结果映射为矿物，获得岩心矿物分布。
[0029]
进一步，所述步骤s41中岩心样品的先验信息包括岩心名称、种类，根据先验信息可获得该岩心样品中常见矿物；若无法获得岩心的任何先验信息，则可将常见成岩矿物作为先验信息来确定分割阈值，常见成岩矿物为本领域悉知内容，在此不做赘述。
[0030]
基于上述技术方案，本发明的有益效果是：
[0031]
1.本发明基于基材料分解模型计算电子密度和有效原子序数，规避了大量标准材料的拟合，操作简便；
[0032]
2.本发明通过动态基材料选择为非均质岩心样品逐像素选择最佳基材料材质，实现了材料分解中基材料的自动选择，使基材料分解算法适用于岩心矿物表征；
[0033]
3.本发明通过动态基材料选择为非均质岩心样品逐像素选择最佳基材料直径，从结果上减弱了射束硬化对材料分解的影响；
[0034]
4.本发明提出了基于最佳虚拟单色图的电子密度和有效原子序数计算方法，以代替传统的经验计算公式，解决了传统公式在致密材料上电子密度和有效原子序数计算误差大的问题，提高了计算精度，保证了表征的准确性。
附图说明
[0035]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单介绍，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0036]
图1为本发明一实施例中基于动态基材料选择的非均质岩心矿物表征方法的流程示意图；
[0037]
图2为本发明一实施例中部分基材料能量图点云分布、基材料能量图点云质心分布和火成岩样品切片能量图点分布示意图；
[0038]
图3为本发明一实施例中基于最佳虚拟单色图的电子密度和有效原子序数计算公式拟合结果示意图，其中图3(a)为基于最佳虚拟单色图的电子密度计算公式拟合结果，图3(b)为基于最佳虚拟单色图的有效原子序数计算公式拟合结果；
[0039]
图4为本发明一实施例中高、低能下火成岩线性衰减系数重建切片示意图，其中图4(a)为火成岩低能线性衰减系数重建图，图4(b)为火成岩高能线性衰减系数重建图；
[0040]
图5为本发明一实施例中火成岩样品最佳基材料材质及直径选择结果切片示意图，其中图5(a)为火成岩样品最佳基材料材质选择结果，图5(b)为火成岩样品最佳基材料直径选择结果；
[0041]
图6为本发明一实施例中7种标准矿物电子密度及有效原子序数计算值及真值对比结果示意图，其中图6(a)为电子密度计算结果及真值对比结果示意图，图6(b)为有效原子序数计算结果及真值对比结果示意图；
[0042]
图7为本发明一实施例中传统基材料分解方法和基于动态基材料选择的非均质岩心矿物表征方法在致密菱铁矿样品电子密度和有效原子序数计算中射束硬化减弱程度对比结果图，其中图7(a)为致密菱铁矿样品重建切片灰度图，图7(b)为图7(a)中第87行的电子密度计算结果，图7(c)为图7(a)中第87行的有效原子序数计算结果；
[0043]
图8为本发明一实施例中火成岩样品的矿物表征结果三维示意图。
[0044]
附图标记：
[0045]
a-火成岩；b-正长石；c-拉长石；d-中长石；e-石英；f-辉石；g-铁矿石；h-其它矿物；i-孔隙及裂痕。
具体实施方式
[0046]
为使本发明的目的、技术方案和优点更加清晰明白，以下结合具体实施例，并参照附图，对本发明作进一步详细说明。以下实施例或者附图用于说明本发明，但并不用来限制本发明的范围。
[0047]
本发明的目的在于提出一种基于动态基材料选择的非均质岩心矿物表征方法，该方法首先利用岩心先验信息建立基材料库；然后对基材料模体和岩心样品进行双能ct成像，获取基材料库能量图点云数据和样品能量图数据，随后进行动态基材料选择，为岩心样品逐像素选取最佳基材料材质和直径；然后利用基于最佳虚拟单色图的电子密度和有效原子序数计算方法根据基材料选择结果计算样品电子密度和有效原子序数；最后依据电子密度和有效原子序数信息通过阈值分割完成矿物表征。本发明基于图像域基材料分解模型进行非均质岩心的矿物表征，避免了大量标准材料的拟合，不仅操作简便，还确保了表征的准确性。
[0048]
参照图1，在一种具体的实施例中，本发明提供的一种基于动态基材料选择的非均质岩心矿物表征方法包括以下内容：
[0049]
一、基材料库数据及公式拟合结果获取
[0050]
1.执行步骤s01：分析常见成岩矿物的先验信息，确定基材料库成员及基材料模体直径。
[0051]
在一个具体的实施例中，首先分析常见成岩矿物的成分信息，确定基材料库由以下5种标准单质材料组成：石墨、硅、铝、钛、铁，所有标准单质材料均被制成一组直径不同的圆柱形模体，直径分别为：4mm，6mm，8mm，10mm，12mm，15mm。
[0052]
2.执行步骤s02：对基材料进行双能ct成像，获取基材料能量图点云库及基材料能量图点云质心库。
[0053]
在一个具体的实施例中，利用基于光子计数探测器的能谱ct对基材料模体进行双能ct成像。其中，低能能谱电压e
l
为：10kv～27kv，高能能谱电压eh为：27kv～55kv。需要说明
square均大于0.98，这证实了本发明中提出的基于最佳虚拟单色图的电子密度和有效原子序数计算公式在矿物表征中的适用性。
[0062]
表1：基于最佳虚拟单色图的电子密度和有效原子序数计算公式拟合所用数据
[0063][0064]
特别指出，上述步骤s01，s02，s03均属于非均质岩心矿物表征方法中的预先准备工作，无需重复进行；上述步骤s01，s02，s03一旦完成即可应用于类似岩心的表征工作，但需要保证岩心样品的双能ct扫描参数与基材料双能ct扫描参数一致；其中，类似岩心指的是岩心中基本矿物所含元素相近。
[0065]
二、岩心样品数据获取
[0066]
执行步骤s11：对非均质岩心样品进行双能ct成像，获取样品能量图数据。
[0067]
在一个具体的实施例中，非均质岩心样品为火成岩样品，首先利用基于光子计数探测器的能谱ct在低能能谱电压e
l
：10kv～27kv和高能能谱电压eh：27kv～55kv下对火成岩进行双能ct成像，其低能和高能线性衰减系数重建结果切片如图4所示；随后将火成岩的低、高能线性衰减系数映射于能量图上，结果见图2。
[0068]
三、算法处理
[0069]
1.执行步骤s21：利用动态基材料选择算法确定最佳基材料材质和直径。
[0070]
在一个具体的实施例中，如图2所示，逐像素在基材料能量图点云质心库o中搜索与样品能量图点欧氏距离最近的基材料能量图点云质心点o
ij
，并将o
ij
对应的一组
基材料材质作为样品进行材料分解的最佳基材料材质(i
opt
，j
opt
)；然后，利用k最近邻搜索算法，逐像素在最佳基材料材质(i
opt
，j
opt
)的能量图点云中搜索与样品能量图点欧氏距离最近的一个点并将该点对应的基材料直径作为样品进行材料分解的最佳基材料直径t
opt
；火成岩样品基材料材质选择结果(i
opt
，j
opt
)如图5(a)所示，基材料直径选择结果t
opt
如图5(b)所示。
[0071]
2.执行步骤s22：根据基材料选择结果进行材料分解，将材料分解系数代入基于最佳虚拟单色图的电子密度和有效原子序数公式获取电子密度和有效原子序数。
[0072]
在一个具体的实施例中，步骤s22包括以下子步骤：
[0073]
(1)材料分解：根据动态基材料选择结果，联立下式求解基材料分解系数b1和b2：
[0074][0075]
式中，和分别表示低能能谱电压e
l
和高能能谱电压eh下样品的重建线性衰减系数，系数，和分别表示低能能谱电压e
l
和高能能谱电压eh下直径为t
opt
的基材料i
opt
和基材料k
opt
的线性衰减系数，由式(4)给出，b1和b2即为所求材料分解系数。
[0076]
(2)电子密度和有效原子序数计算：根据步骤s03的拟合结果及b1，b2的计算结果，利用下式计算火成岩样品的电子密度ρe及有效原子序数z
eff
分布：
[0077][0078]
式中，和为最佳虚拟单色能量e1和e2下的样品线性衰减系数；和分别为基材料i
opt
和基材料j
opt
在最佳虚拟单色能量e1和e2下的理论线性衰减系数，由nist xcom数据库根据基材料成分及密度给出，b1和b2分别为两种基材料的分解系数，ρe和z
eff
分别为所求电子密度和有效原子序数。
[0079]
在一个具体的实施例中，为了验证本发明方法中电子密度和有效原子序数的计算准确性，特利用本发明的方法计算了7种标准矿物的电子密度和有效原子序数，计算结果如图6所示。可以看出，7种矿物的电子密度和有效原子序数计算结果相对误差均在5％以内，说明本发明的方法可以为矿物表征提供足够精度的电子密度和有效原子序数计算结果。
[0080]
在一个具体的实施例中，为验证本发明方法对致密矿物的射束硬化效应的减弱程度，利用传统图像域基材料分解方法计算标准菱铁矿样品的电子密度和有效原子序数，并与本发明方法的计算结果对比，结果见图7。可以看出，利用传统方法，射束硬化导致的误差被传递到电子密度和有效原子序数计算结果中；而本发明方法受射束硬化的影响较传统方法小，计算结果精度较传统方法高。
[0081]
3.执行步骤s23：根据样品的电子密度和有效原子序数分布进行阈值分割，完成非均质岩心的矿物表征。
[0082]
在一个具体的实施例中，首先查阅资料获取火成岩样品中常见矿物信息，并计算
常见矿物电子密度和有效原子序数真值，结果如下表2所示；随后将电子密度和有效原子序数在真值的
±
5％以内的像素点表征为对应矿物，表征结果如图8所示。
[0083]
表2：火成岩中常见矿物及其电子密度和有效原子序数真值
[0084]
矿物电子密度有效原子序数石英2.64611.75正长石2.54913.34中长石2.64813.25拉长石2.67816.61辉石3.36014.85铁矿石4.55021.32
[0085]
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明。应理解的是，以上所述仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，可利用上述揭示的技术内容做出更动或修饰为等同变化的等效实施例。但凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。