天然气水合物在海洋沉积物孔隙中的赋存规律对储层地震波、电阻率以及渗透率等物理特性具有显著影响^[1-4]，从而在很大程度上决定了水合物地球物理勘探和资源评价的精度。迄今，全球已探明的海洋水合物产状仍以分散型为主，如我国南海神狐海域^[5-6]，水合物主要充填于泥质或砂质沉积物的孔隙中。沉积物孔隙的直径通常在几微米到几百微米之间，且水合物在孔隙中生长具有较大的随机性，水合物赋存形态在不同的饱和度条件下呈现出多样性，目前对此尚未形成统一的认识^[7-8]。

野外常用的X射线岩心扫描仪，可对现场钻探获取的柱状岩心进行整体扫描，空间分辨率一般在厘米级，能够获取岩心内部的沉积特征和较大的孔隙或裂隙空间结构，但对于分散于沉积物孔隙中的水合物微观分布细节缺乏表征能力^[9-10]。微米X射线计算机断层扫描（micro-scale X-ray computed tomography，micro-XCT）能够在不破坏水合物赋存条件下直观呈现沉积物内部微观结构特征及其不同组分的相态变化，在现有含水合物沉积物微观探测研究中应用较为广泛^[11]。然而，在利用CT灰度图像来区分微米级孔隙中水与水合物时仍面临技术挑战，这主要是由水和水合物对X射线的衰减系数相近所造成的，表现为CT图像中二者的灰度区间存在交集。在进行水与水合物相态区分时，由于传统阈值分割方法更多注重图像边缘轮廓的分割，对于水合物与水在图像中处于混合状态的情况下选取的灰度阈值缺乏客观的标准^[12-13]，容易导致图像中部分水和水合物的分布范围发生偏差，从而对准确识别二者的边界造成困难。

现有研究中解决这一相态区分问题的思路主要有两种：一是在水中通过添加原子序数较高的化学物质来增加水与水合物二者之间的衰减系数差，如利用碘化钾或溴化钠溶液代替水^[14-15]，或使用氙气、氩气等来增加水合物的衰减系数^[16-17]，在一定程度上能够增大图像中水与水合物的灰度差，但沉积物孔隙中其他物质的加入也会造成水合物生成和分解动力学性质的改变；另一种思路是在低温高压容器内添加人工合成的块状水合物充当标记物，用以作为沉积物体系中新生成的水合物灰度值提取的参考，实验表明该方法能够取得比较理想的分辨效果^[18-20]，但缺点是转移水合物标记物的过程非常复杂，且容易导致用于充当标记物的水合物样品的分解。

此外，由于甲烷水合物的合成需要较为严格的温度和压力条件，部分学者也有使用氙气^[21-22]、二氧化碳^[23]、四氢呋喃^[24]等作为水合物客体分子进行实验，以此推断含甲烷水合物沉积物的性质。但由于客体分子类型的不同，上述类型水合物与甲烷水合物之间生长动力学及空间赋存规律上的差异尚有诸多不确定性，与野外勘探开发实际也有较大差距。

本研究针对水与水合物相态区分的技术难题，在实验室内获取石英砂颗粒孔隙中甲烷水合物生长过程不同时刻批量CT扫描灰度图像^[19-20]，在此基础上提出一种基于图像灰度直方图归一化的阈值分割方法；通过对不同生长时刻与初始时刻所对应的灰度直方图数据进行对比，即灰度值对应像素点减少的区间分别对应于甲烷气体和水，灰度值对应像素点增加的区间则对应于水合物，石英砂颗粒的灰度保持不变，从而确定图像中4种组分的灰度阈值；基于该灰度阈值，对水合物的空间分布形态进行表征。

1.   CT灰度图像及其直方图归一化方法

1.1   灰度直方图归一化

本文采用0.3～0.8 mm石英砂作为沉积物，将其放进可耐高压的反应容器中，并用去离子水将石英砂饱和。设置反应釜温度为2℃，通入6.2 MP甲烷气体后将反应容器密封，即容器内的甲烷气、水和石英砂为一封闭系统，与外界没有物质交换。在水合物生长的过程中，每间隔2 h进行一次CT扫描，获得反应容器内部的灰度图像。

在水合物生长过程中，除沉积物骨架外，孔隙内部的甲烷气、水合物、水3种组分的空间分布和含量处于动态变化中，同时受CT仪器稳定性和图像噪声的影响，不同阶段CT图像的灰度直方图会有所差异，特定组分的峰值灰度及灰度区间会发生偏移，不利于灰度数值的对比分析。本文提出的灰度直方图归一化方法，其核心是矫正CT图像中特定组分的峰值灰度及灰度区间，将不同时刻CT图像灰度直方图中处于最小值的甲烷气与处于最大值的石英砂二者的峰值分别标定在一个固定的灰度值上，对灰度直方图进行归一化，具体归一化步骤如下：

（1）选取甲烷气与石英砂的峰值灰度基准。

（2）提取甲烷水合物CT图像的有效区域，绘制直方图曲线。

（3）分别拟合步骤（2）直方图曲线中的甲烷气与石英砂的高斯曲线，得到甲烷水合物CT图像直方图曲线中甲烷气与石英砂的峰值灰度。

（4）以步骤（3）中拟合得到的高斯曲线峰值灰度作为当前灰度直方图曲线中甲烷气与石英砂的峰值灰度。

（5）通过拟合步骤（1）中选取的峰值灰度基准与步骤（4）中得到的当前直方图峰值灰度的一次函数得到函数系数$a，b$，将灰度坐标$g\Big( {g \in [0,1,2, \cdots ,255]} \Big)$作为自变量带入系数为$a，b$的一次函数中，计算得到新的灰度坐标${g^\prime }\Big( {{g^\prime } = \{ {{x_0},{x_1},{x_2}, \cdots ,{x_n},\;\,n = 0,1,2, \cdots ,255} \}} \Big)$。

（6）以灰度坐标$g$作为归一化直方图曲线横坐标，以灰度坐标${g^\prime }$在原始灰度直方图曲线上对应灰度坐标的纵坐标作为归一化灰度直方图曲线纵坐标，得到归一化灰度直方图曲线。图1为直方图归一化流程图。

图  1  CT灰度图像直方图归一化流程图

Figure  1.  Histogram normalization flow chart of CT gray image

下载: 全尺寸图片 幻灯片

针对水合物生长过程中不同时刻的CT图像，挑选反应时间依次增长的5个时刻的数据进行分析，反应时间从短到长依次为0、30、34、48和72 h。每个时刻的CT切片图像数量为790张，从样品底部到上部的切片编号依次为1～790，每一张CT切片图像尺寸均为530×530像素。图2为选取的5个不同反应时间下，切片编号均为150位置处的CT图像直方图归一化结果，根据图像组分不超出灰度范围的原则，对比分析直方图曲线的峰值和峰宽，选取的甲烷气与石英砂峰值灰度基准分别为35和220。其中图2（a）为归一化之前的CT图像直方图曲线，图2（b）为归一化之后的CT图像直方图曲线。

图  2  灰度直方图归一化对比（切片编号150）

Figure  2.  Comparison of histogram curves before and after normalization（slice 150）

下载: 全尺寸图片 幻灯片

图3为图2（b）中灰度坐标59～135范围内水合物与水混合区域灰度直方图曲线，图中随着反应时间的增长，水合物和水的混合峰的峰值灰度逐渐向低灰度值（水合物）方向偏移，表明了水合物含量的增加和水含量的减少。根据灰度直方图的这一变化趋势，通过定量统计的方法获得增加和减少的灰度区间，可为水合物与水的灰度阈值划分提供依据。

图  3  水合物与水混合峰区域

Figure  3.  Mixing peak region of hydrate and water

下载: 全尺寸图片 幻灯片

1.2   CT灰度图像归一化

在CT图像灰度直方图归一化的基础上，利用直方图归一化的中间参数信息，对图像每个像素点重新赋值，将原始的CT图像与归一化的直方图曲线匹配，获得归一化的灰度图像，有助于识别水合物的分布特征。用公式（1）对原始甲烷水合物CT图像中每个像素点的灰度重新赋值，式中I(x,y)为输入图像不同位置处的灰度值，O(x,y)为输出的对应位置处的灰度值，$ {x_0} $与$ {x_n} $为直方图归一化步骤（5）中灰度坐标${g^\prime }$的最小值与最大值。通过遍历整个输入图像，得到归一化输出图像。

图4为选取的5个不同反应时间下的CT灰度图像归一化前后对比（以切片编号150位置处图像为例），其中图4左列为归一化之前的CT切片图像，图4右列为归一化之后的CT切片图像。图像归一化结果表明，不同时刻CT图像的亮度、对比度等更加一致，且并未影响图像中各组分的空间分布及含量信息。

图  4  图像归一化对比（切片编号150）及三维数字岩心图像

Figure  4.  Image normalization contrast (slice 150) and 3D digital core

下载: 全尺寸图片 幻灯片

2.   CT图像灰度阈值分割

2.1   确定各组分灰度阈值

灰度图像阈值分割是一种基于区域的图像分割技术，核心是按照灰度等级对像素集合进行划分，得到的每个子集形成一个与样品中不同组分相对应的灰度区间。水合物CT灰度图像中甲烷气体、水合物、水和砂各自占据不同的灰度级范围，阈值分割能够更直观地从图像上观察沉积物孔隙中水合物的空间分布特征。根据CT成像原理，在水合物与水的密度相近的情况下，CT图像中水合物与水的灰度分布范围相互重叠，因此在灰度直方图中几乎不会出现水合物与水相互独立的两个峰，从而难以直接从单次CT扫描图像或灰度直方图中对二者进行准确区分。

通过对比归一化后的水合物生长不同时刻CT图像灰度直方图，在水合物与水混合灰度区间内筛选出强度持续增加的灰度范围，即代表水合物的灰度区间，其中，该灰度区间的最小灰度为甲烷气体与水合物的分割阈值，最大灰度为水合物与水的分割阈值。对于水与石英砂这两种组分边界，取甲烷水合物与水的整体曲线的峰值灰度与石英砂峰值灰度的中点作为其分割阈值。

图5展示了水合物生长过程中反应时间30、34、48和72 h下CT切片图像（以切片编号150位置处图像为例）灰度直方图与反应时间0 h对应的灰度直方图差值的变化规律，其纵坐标负值表示在这个过程中强度在减少的灰度范围，正值代表强度增加的灰度范围，从而确定了各种组分的分割阈值（表1）。

图  5  不同时刻CT图像直方图曲线差值（切片编号150）

Figure  5.  Histogram curves of CT images subtracted at different times（slice 150）

下载: 全尺寸图片 幻灯片

表  1  CT图像中各种组分的灰度区间

Table  1.  Gray intervals of various components of CT images

反应时间/h	灰度区间（0～255）
	甲烷气	甲烷水合物	水	砂
0	0～58	－	59～156	157～255
30	0～58	59～96	97～156	157～255
34	0～58	58～96	97～156	157～255
48	0～58	59～96	97～156	157～255
72	0～61	62～105	106～156	157～255

下载: 导出CSV 

| 显示表格

2.2   图像四值化

图6是选取的5个不同反应时间下，CT图像的阈值分割后对应的伪彩色图像以及局部放大图像，其中灰黑色代表甲烷气、黄色代表水合物、青色代表水、白色代表石英砂。结果表明，水合物在孔隙空间中的生长是一个具有空间非均匀性的替代过程，水合物容易在较大的孔隙空间中赋存；生成的水合物固体颗粒内部，存在较多粒内孔，填充了部分甲烷气体或水，其直径主要介于20～100 μm。

图  6  基于本文阈值分割方法得到的伪彩色三维数字岩心及图像的孔隙局部放大图像

Figure  6.  Pseudo-color digital core and local pore magnification image obtained based on threshold segmentation method used in this study

下载: 全尺寸图片 幻灯片

2.3   有效性验证

在CT图像中各组分灰度阈值确定后，基于像素统计的方法即可获得含水合物沉积物各个组分的体积百分含量，因此可采用水合物饱和度参数这一定量指标来对阈值分割的有效性进行验证。沉积物孔隙部分包括甲烷气、水合物和水，骨架部分为石英砂，水合物饱和度定义为水合物体积占据孔隙体积的百分比。

文献[19]提出了一种基于物质界面处灰度值随距离变化曲线识别方法，通过不同物质之间的灰度值梯度变化来确定其灰度分割阈值；文献[20]则采用外加人工水合物颗粒作为标定的方法，确定水合物的灰度值区间，从而判别水合物与气体、水之间的灰度阈值。为验证本文阈值分割结果的准确性，采用以上文献中的两种方法分别对本文CT图像进行阈值分割并计算水合物饱和度参数，用于对比分析，如图7所示。

图  7  饱和度参数对比图

Figure  7.  Comparison of hydrate saturation

下载: 全尺寸图片 幻灯片

结果表明：本文阈值分割方法所获得水合物饱和度数值介于文献[19]和文献[20]方法计算结果之间，与二者平均绝对误差分别为7.8％ 和4.2％，且水合物饱和度参数随着水合物生长过程的演变具有较好的一致性。另外，相比文献[19]中的阈值分割方法，其需要挑选感兴趣区域进行局部灰度值变化趋势分析，不同部位因水合物含量的多少会对阈值选取造成影响，而本研究是基于每次CT扫描的整体三维数据体进行的总体分析，从而避免了因局部水合物分布的空间各向异性造成的差异。

同时，本研究相比文献[20]所采用的外加水合物颗粒作为标定物方法，实验操作更加简便，也能够避免因先后两次合成的水合物笼占有率不同而对CT图像灰度产生的误差。最后，本研究方法的分割阈值是通过一定的模型从数据中合理预测得到的，相比于其他通过直方图曲线的波峰/波谷确定分割阈值的方法，本文中的方法能在一定程度上避免主观因素的干扰。

3.   结论

为准确区分天然气水合物CT图像中水合物与水的边界，提高阈值分割精确度，本文围绕沉积物孔隙中水合物生长动态过程的CT扫描图像创新性提出了基于水合物生长时间演化的阈值分割方法，表征了水合物微观分布形态特征。实验结果表明：

（1）提出的阈值分割方法能够在一定程度上降低天然气水合物CT图像灰度偏移，更为准确、便捷地获得封闭体系中甲烷气、水合物、水、石英砂4种组分的灰度阈值；若实验样品在开放的环境中，气、水或沉积物与外界发生物质迁移，则本阈值分割方法不再适用。

（2）通过该阈值计算的水合物饱和度参数与文献中常用的饱和度计算方法相比，计算结果具有较高的一致性。

（3）在阈值分割的基础上，可对甲烷气、水合物、水以及石英砂骨架四种组分的空间分布进行表征，水合物在孔隙空间中的生长行为具有较大的空间非均匀性，且生成的水合物固体颗粒内部存在较多粒内孔，这些空间分布特征对天然气水合物储层地震、电阻等现场勘探数据的解释具有显著的影响，值得进一步深入研究。