砂岩铀矿岩心颜色的成像光谱自动识别与定名方法研究

2023-12-09潘蔚张元涛师俞晨陈雪娇李子颖

铀矿地质 2023年6期

潘蔚，张元涛，师俞晨，陈雪娇，李子颖

（核工业北京地质研究院遥感信息与图像分析技术国家级重点实验室，北京 100029）

沉积岩石的颜色不仅可以用来推断其形成的古地理环境，还可以用来分析其后生作用过程和铀成矿条件，在砂岩铀矿勘查中具有十分重要的作用。目前岩心编录中颜色的辨识和描述依赖编录人员的感知与经验指示来完成，同种颜色可能由于编录人员的先验知识差异和环境的变化被编录为不同的颜色，从而影响后期的使用和分析。因此，十分必要引入自动化手段，实现颜色识别的定量化和描述的规范化。

遥感技术，尤其是成像光谱遥感技术作为一种数字化获取地面景物色彩与纹理空间分布信息的新手段，具有定量化表述颜色的天然优势。除了SeaStar 和我国海洋水色卫星等专用平台外，ENVISAT 和ADEOS 等陆地卫星也搭载有水色扫描仪，特别是MODIS 和MERISD等中等分辨率卫星的水色传感器，其为海洋水色和陆地目标的颜色识别提供了可靠的卫星遥感数据源。目前，国内外利用地物反射率数据识别目标颜色的方法主要有真彩色三通道法［1-3］和全谱段识别法［4-6］。利用上述两种方法可获取特定反射率下地物相应的颜色数值，但是建立颜色RGB 数值与颜色名称之间的联系仍然是关键。

本文尝试利用成像光谱数据识别岩心颜色，结合石油地质行业的SY/T5751 规范［7］中颜色名称与RGB 数值关系，开展砂岩铀矿相关岩性的颜色识别与编码方法研究，以期为岩心数字化自动编录提供颜色识别与表述的技术支撑。

1 颜色的基本理论与表述方法

1.1 物体颜色

现实世界的物体可以分为自发光和非自发光两种，相应的颜色测量为光源色测量及物体色测量。光源色测量涉及光源的光谱功率及视觉，受光源和观察者影响。光源的相对光谱功率分布通常是将波长为555 nm 处功率标定为100，并将其他光谱波段功率与之进行比较而得到的［8］。光谱功率分布既是光源自身显色的根本原因，又是照射其他物体使其呈色的关键因素。物体色测量涉及光源、物体及视觉，即其颜色除了与自身的光谱反射率（或透射率）有关外还受观察者及光源的影响。

1.2 颜色表示方法

颜色可以用颜色空间（彩色空间或彩色系统）来表示，通常用三维模型来表示，任何一种颜色可用色彩空间的一个点来标识。常见的色彩模型有3 种：CIE 1931 RGB、CIE 1931 XYZ和 CIE 1976 L*a*b*。

CIE 1931 RGB 颜色空间实现了对所有可见光谱色的匹配。其标准色度观察者三刺激值的曲线见图1a。CIE 1931 XYZ 颜色空间消除了红色刺激值的负值区域，其标准色度观察者三刺激值曲线见图1b。该颜色空间的提出奠定了颜色测量的基础［9］，并可用光谱仪测量结果进行计算。CIE 色度系统颜色三刺激值X，Y，Z 定义为：

图1 CIE 1931 RGB 和CIE 1931 XYZ 刺激值曲线与CIE L*a*b*颜色空间Fig.1 The CIE 1931 RGB and CIE 1931 XYZ color matching functions and CIE L*a*b* color space model

式中：φ（λ）为进入人眼并能产生颜色感觉的光能量，称为颜色刺激函数；xˉ（λ）、yˉ（λ）、zˉ（λ）为标准观察者不同波长的三刺激值，其中λ为波长；dλ为波长间隔，k 为归一化常数，目的是使照明体（或光源）的Y值为100。针对非自发光不透明物体，φ（λ）=S（λ）×R（λ），S（λ）为照明光源的相对光谱功率分布，通常取D65 光源的相对光谱功率分布，R（λ）为被观察物体的反射光谱［10］。

CIE 1931XYZ颜色空间虽然克服了RGB颜色空间的缺点，但他的色品图不均匀，致使XYZ颜色空间中，相同的距离不能表示相同的颜色感知差异，不利于实际应用。为此，国际照明委员会CIE（Commission Internationale de l’Eclairage，CIE）1976 年提出CIE1976 L*a*b*模型（图1c）。

图1 中，L*代表明度，表示黑/白。a*和b*代表色度，a*表示红绿，+a*表示红色，−a*表示绿色。b*表示黄蓝，+b*表示黄色，−b*表示蓝色。样品颜色L*，a*，b*的值可通过对X，Y，Z三刺激值的非线性变换得到。相应的公式为：

式中：X、Y、Z为颜色样品的三刺激值。Xn、Yn、Zn为CIE 标准照明体照射到完全漫反射体后再反射到观测者眼中的三刺激值。对于D65 光源（日光，色温6 500 K），Xn，Yn和Zn分为95.047、100 以及108.883［1］。实际应用中，除了由样品的三刺激值X，Y，Z求解L*，a*，b*之外，也常用到相反的过程，即由样品的L*，a*，b*坐标反算样品的三刺激值X，Y，Z。

除了以上三种色彩模型外，惠普和微软公司于1996 年联合提出一种专门用于计算机显示器、打印机和因特网的标准颜色空间——sRGB 颜色空间。其作为硬件设备本身的颜色空间，优点在于与传统颜色管理相比可不必进行独立设备颜色空间和依赖设备颜色空间的相互转换。

1.3 色差的计算与评价

为了测量、比较样本与参考目标颜色的差异，前人相继提出了多种色差公式，如CIE 1976 L*a*b*色差公式和CIEDE2000 色差公式等。

CIE 1976 L*a*b*色差公式广泛应用于色差计算，其通过计算两种颜色样本的欧几里得距离评价二者的差异。在CIEL*a*b*均匀色空间中，两种颜色的总色差可用下列公式计算。

式中：△L*为明度差，△a*和△b*为色度差。当总色差很小接近于0 时反映二者颜色相同。

由于CIE 1976 L*a*b*色差公式缺乏感知一致性，因而产生了其他色差评价指标。CIEDE2000 是CIE 于2000 年确定的一个新的色差评价公式［11］，旨在产生更加接近人类视觉系统的感知结果，色差可表达为：

式中：L为亮度，△L’为明度差；C为色度，△C’为彩度差；H为色调，△H’为色相差；RT为旋转函数，用来校正椭圆轴方向蓝色区域的偏转；kL，kC以及kH为权重因子；SL，SC以及SH为对LCH 的补偿因子［12］。为保持与传统地质编录颜色描述依靠人眼感知的一致性，本次色差评价采用CIEDE2000 公式进行计算。

2 砂岩铀矿岩心颜色的成像光谱识别

2.1 砂岩铀矿相关岩性的颜色特征

砂岩铀矿是一种产于沉积盆地中的矿床，与砂岩铀矿相关的岩性是碎屑岩，其中主要有砾岩、砂岩、泥岩、少量煤层和页岩。碎屑岩的颜色取决于岩石的碎屑矿物成分，特别是染色成分［13］。由长石、石英和白云母等组成的岩石主要为浅色调，如灰褐色、灰色、灰白色；当岩石中染色物质，如铁化合物（绿色、红色、褐色、黄色）和游离碳（灰色、黑色）等物质占据一定量时，岩石会呈现丰富的颜色。染色成分的存在与碎屑岩的形成环境和后生改造相关，因此岩石的颜色可以反映其形成环境和后生改造条件（表1）。

表1 岩石颜色与其形成环境（后生改造条件）关系Table 1 The relationship between rock color and its formation environment or epigenetic transformation conditions

由于碎屑颗粒经过搬运和分选，矿物成分相对简单，沉积成岩时有限范围内的环境相对统一，因此碎屑岩的颜色相对均匀。特别是在钻孔岩心的孔径内，颜色的变化主要表现在深度方向。因此本文钻孔岩心的识别与编码方法研究主要针对钻孔深度方向。

2.2 岩心成像光谱数据的特点

本文使用的岩心成像光谱数据由挪威NEO 公式生产的HySpex 地面成像光谱仪获得，该仪器由独立成像的可见—近红外成像光谱仪和短波红外成像光谱仪组成。其中，可见—近红外成像光谱仪VNIR-1600 探测波长范围为400～1 000 nm，光谱采样间隔为3.7 nm，所获数据包含160 个波段（表2）。

表2 VNIR-1600 主要技术指标Table 2 The main specifications of VNIR-1600

虽然VNIR-1600 获得的可见—近红外数据的光谱范围与颜色三刺激函数光谱范围略有不同，但覆盖了肉眼最为敏感的400～700 nm 颜色光谱，光谱响应范围占三刺激函数的93.3%以上，有效响应波长在99.9%以上。

2.3 岩心颜色的全谱段识别方法

全谱段识别是指利用成像光谱仪采集的可见—近红外（VNIR）波长范围内全部反射率数据，根据色度学原理和颜色表示方法，对岩心颜色进行分类和识别。本文设计的HySpex VNIR 数据岩石颜色识别流程主要步骤包括X、Y、Z 三刺激值计算，CIE 1976 L*a*b*色度值计算、色差计算及颜色名称确定等（图2）。

图2 岩心颜色的成像光谱可见—近红外全谱段识别技术流程Fig.2 The technical process of full VNIR imaging spectrum identification for core color

1）X、Y、Z 三刺激值计算

1931 年国际照明委员会（CIE）规定了标准观察者的光谱三刺激值，从而为颜色测量奠定了基础［9］。由地物光谱反射率分布到颜色三刺激值X、Y、Z，可利用式1 实现。实际计算中，照明光源的相对光谱功率分布S（λ）采用D65 光源的值，且由于所用光源光谱功率分布的光谱范围为380～700 nm，光谱分辨率为10 nm，而HySpex VNIR-1600 获取影像波谱范围为400～1 000 nm，光谱分辨率为3.7 nm，为此将影像光谱进行插值使其与前者匹配。

2）CIE 1976 L*a*b*色度值计算

L*a*b*色度值通过对颜色的量化描述，反映了颜色的亮度、色调以及饱和度信息。利用公式（2）和（3）可将颜色值从CIE XYZ 空间转到CIELAB 均匀色空间。

3）标准岩石颜色代码转换

石油地质岩石名称及颜色代码（SY/T 5751—2012）对112 种常用岩石标准颜色进行了编码，包括岩石名称、孟塞尔颜色及RGB 颜色［7］。部分颜色名称与RGB 参数如表3 所示。

表3 石油行业沉积岩部分颜色名称与RGB 参数（SY/T5751—2012）Table 3 Color name and RGB parameters of some sedimentary rock color in the petroleum industry (SY/T5751—2012)

为了在CIE 1976 L*a*b*均匀色空间进行颜色对比分析，需将包括表3 在内的112 种标准颜色的RGB 值转换为CIE 1976 L*、a*、b*值。主要包括RGB 值转CIE XYZ 值以及CIE XYZ值转CIELAB 值两个步骤。

RGB 值转CIE XYZ 值计算公式：

式中：R、G、B 数据在转换前需要缩放到0 和1之间［1］。

CIE XYZ 转CIELAB 的转换关系如公式（2）和公式（3）所示，通过非线性变换实现。

4）颜色比对与名称确定

利用CIEDE2000 色差公式，分别计算待定颜色与112 种参考颜色之间的差异，将色差最小的参考色名称及RGB 值作为待定岩石的颜色名称及RGB 值，以实现岩心颜色识别与定名。

2.4 岩心颜色的真彩色三通道识别方法

真彩色三通道法是利用成像光谱数据对应于3 个真彩色通道：蓝色通道（450～520 nm）、绿色通道（520～600 nm）和红色通道（600～690 nm）光谱反射率平均值乘以255 得到颜色的亮度值。主要步骤包括平均反射率计算，RGB 值计算，RGB转CIELAB以及颜色比对与名称确定等步骤（图3）。

图3 岩心颜色的成像光谱真彩色三通道识别技术流程Fig.3 The technique flow of true RGB color bands recognition for core color

1）平均反射率计算

分别读取岩心光谱在600～690 nm（红色通道），520～600 nm（绿色通道）和450～520 nm（蓝色通道）的光谱反射率数据，并计算平均值。

2）RGB 值计算

对红、绿、蓝3 个通道获取的反射率平均值分别乘以255 进行RGB 颜色编码。

3）RGB 转CIELAB

由RGB（sRGB）到CIELAB 可分两步完成，主要包括RGB（sRGB）转CIE XYZ 以及CIE XYZ转CIELAB。前者可借助公式（6）完成，后者可结合公式（2）和公式（3）实现。

4）颜色比对与名称确定

三通道法色彩计算与颜色定名与全谱段法采用相同的计算公式和步骤完成。

2.5 编程实现

为实现上述两种颜色自动识别算法与定名，本文基于Python 编程语言进行了程序开发，涉及的库有pandas、pectral、colormath 及scipy 等。主要技术思路是先将岩心箱扫描图像经过单根裁切后的岩心图像平均切分为50小段，以小段内所有像元的平均光谱计算段内岩心所呈现出的总体颜色信息。

岩心小段的切分采用滑动窗口进行，每次以固定的方框获取相应的岩心成像光谱影像。方框在岩心深度方向以固定的步长进行移动，水平方向以岩心影像中心为基准左右扩展1/3图像的宽度进行统计，而不是整个岩心大小。这是因为岩心裁切时边缘可能存在岩心箱残留和阴影，如果计算在内会干扰岩心反射率计算的准确性。统计2/3 岩心宽度是考虑到沉积岩在水平方向上物质成分的稳定性。

为建立颜色参数与颜色名称之间的联系，先将石油行业中的112 种颜色名称和其对应的RGB值存储为数据字典，然后通过循环，在CIELAB 均匀色空间分别计算待判定颜色与112 种颜色之间的CIEDE2000 色差，最后将色差最小的颜色名称赋予待定颜色，并与对应的RGB 值一同输出。

颜色识别结果最终以CSV 文件进行输出。文件中每一行记录每个小段的起止深度、RGB计算值、颜色名称和相应的参考RGB 值。

2.6 干扰去除与结果优化

岩心提取和放入岩心箱时，通常会产生折断甚至破碎，致使岩心段之间存在一定间隙（图4）。由图4 可见，岩心箱中某格岩心存在6处间隙（红色矩形标注），影响了统计单元内的平均光谱反射率计算值，造成本属于同一种颜色的岩心段内出现孤立异常值。此外，岩心回次标签等杂物也会影响识别结果。

图4 岩心段之间的空隙Fig.4 The gap between core segments

本文采用了比较当前分段与前后相邻分段的颜色异同的方法对这些干扰带来的异常进行校正，即当当前统计小段与前后相邻分段颜色相同时，则将当前分段的颜色判断为与前后段相同。若不同，则不作处理。图5显示灰红色异常颜色被修正为浅棕色，表明该干扰去除方法能有效去除异常颜色。

图5 异常色与命名结果优化示意Fig.5 The schematic diagram of abnormal color removal and nomination optimization

3 典型钻孔岩心颜色识别实验

3.1 实验对象

本次实验采集ZKH3 的钻孔位于鄂尔多斯盆地西南部的庆阳地区，选取岩性复杂、颜色丰富的6 段岩心作为实验对象。钻孔钻遇的地层主要有下白垩统罗汉洞组（K1lh）、径川组（K1j）、新近系（N）及第四系（Q）。罗汉洞组主体为风成沙漠相夹河流相沉积，岩性为浅黄色、黄褐色、浅红色中细砂岩，夹薄层灰色、青灰色的粉砂岩、粉砂质泥岩；泾川组主要为河湖相沉积，主要岩性为褐黄色、灰绿色中细砂岩、浅红色泥岩、青灰色泥灰岩，泥质成分增加，泥灰岩沉积较多。

3.2 识别过程

对采集的3 段岩心的可见—近红外成像光谱数据，分别采用全谱段法和真彩色三通段法开展颜色识别实验。两种方法的技术流程是对每段岩心进行颜色参数计算、色差计算和颜色定名，并利用干扰去除和优化流程，对每段岩心的识别结果进行检查和结果优化。表4 展示了两种方法对孔深1 183.6～1 185.2 m 岩心段的识别结果。表中给出了每个颜色段的起止深度、计算得到的RGB 参数值和对应的颜色名称。

表4 ZKH3 孔某岩心段颜色识别结果对比Table 4 Comparison of color recognition results for a core section of ZKH3

3.3 识别结果分析

基于优化后的RGB 值及起止深度，绘制了部分岩心段岩心颜色分布柱状图，并将其与地质资料、彩色数字相机照片以及高光谱真彩色通道合成影像进行对比分析。

图6是243.9～245.2 m岩心段颜色识别效果，该段岩心地质编录为浅灰白粉砂质泥岩，上段颜色深、色杂，下段岩心整体呈灰白色。全谱段方法将岩心分为12 个色段，颜色类型和分段位置与彩色相机照片的基本一致，色彩分段没有受岩心段之间空缺形成的阴影影响。真彩色三通道法将此段岩心分为15 个色段，颜色类型和分带位置与高光谱真彩色波段合成影像基本相同，但与彩色相机照片相比，颜色总体色调偏暗分段位置相差大。与全谱段方法比较，孔深245.0 m 处的深灰色段、244.22～244.32 m 和244.62～244.68 m 处的灰色段错误识别主要受到了岩心段之间空缺的影响，说明三通道法不仅颜色复原效果不理想，而且在去除岩心段间空隙方面效果不好。

图6 ZKH3 孔孔深243.9～245.2 m 两种方法岩心颜色识别效果对比Fig.6 Comparison of the two core color recognition methods in ZKH3 at depth 243.9～245.2 m

图7 展示了孔深976.6～977.9 m 岩心段颜色识别效果，该段岩心为灰色泥岩、浅灰细砂岩及浅灰色泥岩。全谱段法将此段岩心识别为13 个色段，颜色类型与分段位置和相机照片显示的颜色几乎一致。唯一的不同是在977.8 m以下靠近岩心箱的位置出现了2.6 cm 深灰色段，是受残留的岩心箱背景的干扰所致。真彩色三通道法将此段岩心识别为8 个色段，与真彩色相机照片和高光谱真彩色通道合成影像的颜色都有明显的差别，特别是在孔深977.4 m以后，不同程度的灰色没能识别出来。

图7 ZKH3 孔孔深976.6～977.9 m 两种方法岩心颜色识别效果对比Fig.7 Comparison of two core color recognition methods in ZKH3 at depth 976.6～977.9 m

图8 展示了1 183.6～1 185.2 m 岩心段颜色识别效果。该段岩心上段为砖红色中粗砂岩，下段是黄色粗砂岩。全谱段法将此段岩心识别为11 个色段，颜色类型和分段位置与彩色相机照片基本一致，仅1 183.60～1 183.63 m 处深褐色编码段明显受岩心箱干扰，出现了偏差。真彩色三通道法将此段岩心识别为8 个颜色段，总体与高光谱真彩色波段合成影像的颜色类型相符合，但在分段位置上有差别，主要在1 184.83 m 以后。在真彩色图像上，1 184.83 m以后出现了两次灰白相间的颜色变化，但颜色识别编码结果上没有体现，表明真彩色三通道法对颜色细微差别的识别存在不足。同时，三通道法识别效果与全谱段法比，不仅色调偏暗，而且颜色类型与分段位置差别明显。