基于图像纹理合成的特征可调控的数字岩心重构方法

2021-04-06查文舒刘镇领李道伦

合肥工业大学学报（自然科学版） 2021年3期

胡进, 查文舒, 刘镇领, 李道伦

(1.合肥工业大学数学学院,安徽合肥 230601; 2.渤海钻探油气井测试分公司, 河北廊坊 065007)

0 引言

重构数字岩心的方法[1-4]主要分为物理实验方法和数值重建方法两大类。物理实验方法借助高精度仪器获取岩心的平面图像,然后对平面图像进行重建即可得到数字岩心。物理重建的常用方法有:序列成像法[5]、核磁共振成像法[6]、X射线计算机层析成像(X-CT)扫描法[7]和聚焦离子束电子显微镜(FIB-SEM)扫描法[8]等。这些方法能够对岩心样本的真实孔隙结构进行定性描述, 但视野小,而且花费代价大。为便于分析与获得更大尺度孔隙结构,基于数字岩心重建与重构方法得到了重视。

数值重建方法是利用得到的电镜扫描图片, 结合不同的数学统计方法对岩样进行重构。这类方法统计参数的取得方式比较灵活, 花费的成本比较低，不仅可以大范围地进行重构,而且可以满足重构对精度的要求。基于数值的重构方法主要有过程法[9]、顺序指示模拟技术[10]、模拟退火法[11]、多点地质统计法[12]、马尔科夫链蒙特卡罗方法[13]等。

多点地质统计法一直是数值重建研究的热点,它是相对于两点地质统计法[14]而言的。传统的两点地质统计方法利用变差函数来研究地质变量空间的相关性,主要把握空间上两点之间的相关性,难以表征复杂的空间结构和再现复杂目标的几何形态。针对这种方法的不足,多点地质统计学方法应运而生。在多点地质统计学中, 应用训练图像[15]代替变差函数表达地质变量的空间结构性, 这使得多点地质统计学建模方法与计算机图像纹理合成领域展示出目的的一致性,即都是随机生成与训练图像拥有相似特征的图像[16]。

多点地质统计学方法与计算机图形学算法的融合逐步得到了学者重视。文献[17]首次将计算机图形学算法中的图像缝合算法用于多点地质统计的条件模式模拟,并命名为CIQ算法。CIQ方法与多点地质统计学模式模拟方法直接用新模拟区域的值更新重叠区域值的方法不同,它是将重叠区域沿着重合区域误差最小处进行缝合。缝合方法的使用缓解了原先直接取代方法造成的重叠区域边界突出的缺陷[18],得到了更好的模拟结果。

模拟结果的对比如图1所示。这种方法适用于规则矩形的重合区域。

图1 重叠区域直接覆盖和最小误差拼接的模式模拟效果示意图

文献[19]将计算机图切算法应用于地质岩相反演模拟。通过使用马尔科夫链融合后验信息,图形切割技术能够从训练图像中找到合适的模式加入到模拟结构中。

CIQ算法和计算机图切算法都成功应用在地质模拟中,对于特征信息量较小的河道图像,重构效果好。但是，对于本文中有着大量特征信息的岩心图像,在用上述算法对其重构时,会出现部分特征的缺失和重复等不足。本文在文献[17,19]基础上,针对大样本的页岩岩心数字图像,提出了一种基于多样子图的图像纹理重构算法。这种方法从训练图像中挑选出特征子图,可以避免合成图中的重要信息缺失;同时设置比例,也可避免合成图中的大面积的重复;另外每次都是从挑选的特征子图中选择匹配块,避免了全局搜索,大大提高了合成速度。

1 基于比例调控的页岩图像重构算法

1.1 特征子图挑选及预设比例

页岩气储层孔隙结构特征对页岩气勘探开发具有重要意义[20]。通过扫描电子显微镜可以对这些微裂隙的微观结构进行更直观地观察,获得孔隙的大小、形态及有机质分布等,从而为研究页岩气的成藏机理、运移规律及开发提供基础数据。本文用页岩岩心的电镜扫描数字图像作为原始训练图像,如图2所示。为去除噪声,对图2进行三值化处理[21],通过设定合适的阈值,将图像上像素点的灰度值变为0、128或255,其中0表示有机质,255表示白色矿石结果,如图3所示。

图2 原始训练图像

为了避免合成结果中出现大面积的特征信息重复或缺失,本文方法不直接从图3中挑选匹配块,而是首先人工选取m幅的特征子图(本文中m取11),如图4所示。这些特征子图包含着原始页岩岩心图像中的矿石、裂缝、有机质等重要信息。然后给特征子图预设定在合成中的比例,依次为a1,a2,…,am。通过设置不同的ai值(i=1,2,…,m),可以控制各组分在合成图像中的比例,达到比例可调控的目的。

图4 11幅特征子图

1.2 匹配过程中特征子图的选取

本文算法的合成过程与文献[22-23]类似, 每次从训练图像中挑选ε个最相似的匹配块随机选择一个进行匹配,然后用Graph Cut算法决定哪部分将输出到合成图中, 这样不断地扩充合成区域直至完成整幅图像。

不同的是,本文中的训练图像是从图3中挑选出的多幅特征子图,即有多个训练图像,因此在匹配的过程中第1步要解决的问题就是训练图像的选择。

(1)

bm表示在合成过程中每幅子图所占当前合成区域的比例与其预设比例的差值。在合成过程中尽量要求每幅子图的Nm/N接近它的预设比例am,因此每次选择最小的bm值对应的特征子图作为下一个训练图像,这样就能保证在最后的合成图像中每幅特征子图所占比例与预设比例大致相同。

1.3 Graph Cut 算法

重叠区域的拼接也是影响重构结果的关键因素。本文沿用计算机图切技术中的Graph Cut 能量最小化算法完成重叠区域的拼接,它最早是由文献[24]提出的。

本文通过复制训练图像中的矩形像素块到输出图像中合成纹理,再通过图切算法在矩形纹理块中选出最优的不规则区域用于合成。

在Graph Cut算法执行过程中,最佳路径由许多相关的像素点组合而成,像素点的生成是从每一对像素点中选择而来,最简单的匹配标准是比较相邻像素点的颜色误差。设t、u为重叠区域中的2个相邻像素点,B(t)、C(t)分别为像素点t在旧纹理块和新纹理块上的像素颜色值,B(u)、C(u)分别为像素点u在旧纹理块和新纹理块上的像素颜色值。连接t、u2个顶点的边的值定义为：

M(t,u,B,C)=|B(t)-C(t)|+

|B(u)-C(u)|

(2)

这样的边值对应着重叠区域的模式差异值大小,使用基于最小能量的图像剪切技术可以将2个模式以最小误差的方式拼接起来,过程如图5所示。图5中的左边为重叠区域示意图，右边红线代表接缝。

图5 2个模式以最小误差方式的拼接过程

1.4 算法的主要流程

(1) 根据实验需要,对原始的数字岩心图像做三值化处理,提取图中的矿石、有机质、裂缝等重要信息。

(2) 从处理后的三值化图像中选取m幅特征子图作为训练图像。

(3) 输入m幅训练图像、初始参数(匹配块大小p、重叠区域大小O和候选匹配块数量ε等)。

(4) 计算bm的值,从对应最小的bm值的训练图像中随机选取一个匹配块开始重构。

(5) 迭代直到整个图像重构完成。步骤如下:① 根据重叠区域的相似性大小,在对应bm值最小的训练图像中搜索找到下一个匹配块;② 从ε个最相似的匹配块中随机选择一个输入到待合成区域中;③ 在待合成区域中使用Graph Cut技术将待合成图像与匹配块的重叠区域进行误差最小化拼接。

2 实验与结果分析

2.1 图像重构

三值化后的页岩岩心数字图像(图3)大小为1 279×887,特征子图为图像(图4)大小依次为209×206、210×202、212×210、221×211、222×211、220×209、200×205、222×212、618×497、562×494、586×468。初始参数p、O和ε的设定对于本文的重构结果也具有一定的影响[25]。从经验上讲,重构过程中匹配块大小为原始训练图像的1/4～1/6之间为最好,故本文匹配块大小p=200×200。对于重叠区域,如果使用较小的重叠区域,无法保持空间结构的连续性,但是使用过大的重叠区域,会导致重构时间的大大增加,故重叠区域大小设置为匹配块的1/6为宜(O=35×200)。而对于候选匹配块数量ε,取值太小容易出现大面积的训练图像复制,增大候选匹配块数目减少了对训练图像的复制,增加了匹配的多样性,故本文中ε=10。

本文算法实现的部分重构结果如图6所示,大小为1 000×1 000。CIQ算法在相同的初始参数条件下实现的部分重构结果如图7所示,大小为1 000×1 000。