实验用大庆油田松辽盆地某区块采集的岩芯样本共计230个, 在对岩心样本进行分析和数据采集前使用液氮冷冻保存, 按照GB/T19145—2003沉积岩中总有机碳的测定方法测量其总有机碳的含量。

使用傅里叶变换近红外光谱仪TANGO(德国BRUKER公司)采集光谱数据, 波数范围: 11 542~3 940 cm^-1, 分辨率为8 cm^-1, 扫描32次取平均值。 图1所示是全部230个油页岩样品的平均光谱。 从光谱曲线来看, 岩石样本光谱曲线基线漂移较为严重, 在7 300, 5 200和4 500 cm^-1三个波数附近有明显的吸收峰, 同时伴随一定量的噪声。 8 800和4 200 cm^-1附近有吸收峰, 但不是很明显。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 原始光谱图Fig.1 Original spectra

由于受到环境因素、 样本来源多样性以及测量仪器等客观因素的影响, 导致一些样本的光谱偏离样本的总体分布, 这些所谓的异常样本引入将会导致模型的预测精度大幅下降^[10]。 因此, 对异常样本进行剔除是保证定量模型可靠的必要条件。 采用蒙特卡洛交叉验证算法对异常样本进行剔除。 蒙特卡洛随机取样(Monte Carlo sampling, MCS)法每次随机抽取一定比例的样本(占样品量的80%)构成校正集建立偏最小二乘模型, 剩余的20%作验证集对模型进行验证, 计算验证集残差。 经过多次抽样建模后能够得到多个预测残差, 计算出这些预测残差的均值与方差, 将样本预测值误差高于平均残差的样本标记为异常样本。 最后通过校正集相关系数R²、 交叉验证均方差RMSECV、 预测均方差RMSEP对模型进行评价, 验证剔除异常样本是否有利于模型精度的提高。 剔除异常样本前后的建模结果如表1所示, 由表可以看出剔除异常样本后的模型的性能参数有明显改善。

表1

Table 1

表1(Table 1)

表1 剔除异常样本前后的建模结果 Table 1 Modeling results before and after removing outlier samples 剔除方法 样本个数 Rc2 RMSECV Rv2 RMSEP 未剔除 230 0.660 8 0.419 0 0.624 6 0.443 6 MCCV 216 0.718 9 0.366 3 0.688 4 0.386 3

表1 剔除异常样本前后的建模结果 Table 1 Modeling results before and after removing outlier samples

为了消除噪声干扰和基线漂移对模型性能的影响, 一般在建模之前对光谱数据进行预处理。 目前比较常用的近红外光谱预处理方法有Savitzky-Golay(S-G)卷积平滑、 基线校正(baseline correction, BSC)、 标准正态变量变换(sandard normal variate correction, SNV)、 一阶导数(first derivative)、 二阶导数(second derivative)、 去趋势(detrend, DT)等。 S-G能有效提高光谱的平滑性, 减少高频噪声干扰。 SNV主要是减少固体颗粒物大小不均和物体表面散射以及光程变换对光谱数据的影响, 从而达到去除噪声的目的。 DT算法用来处理漫反射光谱基线漂移的问题, 一般和SNV组合使用。 导数方法用来消除背景干扰和基线校正, 提高分辨率和灵敏度。 通过对比不同预处理方法的效果最终确定适合油页岩样本的近红外光谱预处理方法。

近红外光谱筛选波长后所建立的模型相比于全光谱模型, 不仅模型变量数大幅度减少, 而且模型性能也有大幅提升。 进行特征波长选择可以通过简单的模型来提高模型的解释性, 通过减少噪声或者干扰选择出效果更好的变量, 并且能够提高模型的预测能力^[11], 在众多的波长选择算法中, 连续投影(successive projections algorithm, SPA)算法, 无信息变量消除(uninformative variables elimination, UVE)算法和竞争自适应重加权(competitive adaptive reweighted sampling, CARS)算法具有一定的代表性, 波长选择结果较优。

1.5.1 连续投影算法

连续投影算法是一种前向循环的特征变量选择算法, 即从一个波长开始, 其他波长向这个波长向量的法平面投影, 投影长度最长的波长向量被选择为特征波长, 然后以新选择的波长为基础, 重复上述投影过程, 直到达到指定的波长个数为止。 该算法可以从众多光谱信息中筛选出重要样本变量的波长, 用少数几列光谱数据来概括大部分光谱信息, 降低了模型的复杂度, 有效提高建模的速度和模型的稳定性。

1.5.2 无信息变量消除法

无信息变量消除算法是通过向样本光谱矩阵中人为引入随机噪声, 并在此基础上建立偏最小二乘回归交叉验证模型, 得到的偏最小二乘回归系数均值与标准差的商作为衡量波长重要性的关键指标, 将噪声矩阵的最大值作为阈值, 大于阈值的变量被作为优选的特征向量。 UVE算法可以去除没有贡献的变量, 减少模型的运算量, 增强模型的适应性。

1.5.3 竞争性自适应重加权算法

竞争自适应重加权算法是将每个波长变量看作是一个单位个体, 将适应力弱的个体剔除, 从而保留适应性强的个体。 通过蒙特卡洛采样和PLS模型的测定系数进行特征波长选择, 首先通过蒙特卡罗采样随机选择校正集样本建立PLS模型^[12]。 计算该模型各参数的系数权重, 然后利用CARS算法筛选出PLS模型回归系数绝对值权重大的波长, 去掉权重小的波长, 模型交叉验证的均方根误差最小的波长组合即为选择的特征波长。 该算法可以有效保留特征变量及相关影响变量, 剔除冗余及噪声变量。

1.6.1 支持向量机

支持向量机是建立在统计学习理论的VC维理论和最小化结构风险基础上的一种有监督二分类机器学习算法。 该算法的基本思想是将低维非线性问题转换成高维线性问题来分类, 通过非线性变换将输入变量映射到一个高维的特征空间, 并在新的空间中进行线性回归寻找一个最优的超平面, 使得所有样本到超平面的距离最小, 从而解决常规空间中样本线性不可分的问题^[13]。 支持向量机能够较好地解决高维空间中遇到的维数灾难问题, 具有良好的泛化能力, 在解决小样本, 非线性样本分类以及高维模式识别中表现出很多特有的优势, 并对异常样本及噪声具有很好的鲁棒性。

1.6.2 随机森林

随机森林算法是基于决策树和自助重采样法的一种集成学习算法。 它的基本思想是利用自助重采样法不断生成训练变量集和检验变量集, 由检验变量集随机生成多个分类决策树, 每个决策树节点的分裂变量也是随机产生, 从而形成随机森林。 随机森林通过对产生的决策树进行投票得出最终预测结果, 其结果具有较高的准确性。 该算法优点体现在训练速度快且不易出现过拟合, 能够很好的处理数量大的样本, 并且对噪声具有较强的鲁棒性^[14]。

1.6.3 模型评价指标

通过内部交叉验证和外部验证对模型效果进行检验, 采用模型值和真实值的相关系数R²和均方根误差(root mean square error, RMSE)作为模型的评价指标。 R²越接近1, 模型RMSE值越小, 说明模型的预测效果越好。

针对光谱样本的消噪和基线矫正需求, 分别采用S-G卷积平滑、 SNV、 BSC、 DT、 一阶导数和二阶导数6种方法以及组合处理共8种方法对光谱数据进行预处理, 并以PLS模型的模型参数作为预处理效果评价光谱预处理方法。 不同的光谱预处理方法的结果如表2所示。 其中, 基于DT方法(DT, DT+SNV, DT+BSC等)的预处理方法模型性能优于其他方法。 这是因为, 岩石样本近红外光谱基线漂移较为严重(见图1), DT预处理方法能在一定程度上消除漫反射光谱的基线漂移, 因此效果较好。 所有预处理方法中, 在去趋势基础上进行基线校正的DT+BSC预处理方法的PLS模型性能最优。 经过预处理之后可以看出11 541.94~8 872.71 cm^-1波段共580个波长点没有明显的吸收峰, 且噪声明显。 因此, 在后续的研究中只对DT+BSC处理后的8 864.473~3 946.174 cm^-1波段共1 265个波长点的光谱数据进行分析。

表2

Table 2

表2(Table 2)

表2 不同预处理方法下的建模结果 Table 2 Modeling results under different pretreatment methods 光谱预处理方法 Rc2 RMSEC Rv2 RMSEV 无 0.718 9 0.366 3 0.688 4 0.386 3 DT 0.746 0 0.341 4 0.621 6 0.382 0 Derivative-1 0.706 7 0.303 8 0.643 6 0.310 6 Derivative-2 0.719 8 0.365 8 0.663 9 0.342 3 SNV 0.548 7 0.464 2 0.376 0 0.426 5 S-G 0.729 3 0.359 6 0.695 3 0.382 4 BSC 0.658 5 0.403 8 0.570 8 0.449 4 DT+SNV 0.739 3 1.352 7 0.578 0 0.334 8 DT+BSC 0.759 1 0.315 9 0.734 8 0.357 3

表2 不同预处理方法下的建模结果 Table 2 Modeling results under different pretreatment methods

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 DT+BSC预处理的光谱Fig.2 Spectra after DT+BSC pretreatment

预处理后的光谱数据采用SPXY算法按照2:1的比例对剔除异常样本之后的216个油页岩样本进行划分, 得到校正集样品144个, 验证集样品72个。 样本集划分结果见表3, 由表3可知, 校正集和验证集样品的分布比较均匀, 校正集总有机碳样品含量基本涵盖了验证集。

表3

Table 3

表3(Table 3)

表3 样本集划分 Table 3 Sample set division 样本类型 样本数 最小值/% 最大值/% 均值/% 标准差 总体样本 216 0.053 3.922 1.835 0.693 校正集 144 0.053 3.922 1.897 0.747 验证集 72 0.699 2.783 1.713 0.551

表3 样本集划分 Table 3 Sample set division

2.2.1 SPA算法

图3为SPA算法提取不同数量特征波长对应的模型的RMSE, 从图中可以看出RMSE的值随波长数的增加而降低, 当选取波长个数为16时, 模型均方根误差到达稳定且最小。 被选中的波长分别是8 918, 8 790, 7 348, 7 138, 7 068, 7 002, 5 643, 5 412, 5 272, 5 132, 4 526, 4 481, 4 374, 4 353, 4 325和4 213 cm^-1, 其分布情况如图4所示。 这些波长与芳烃和CH₃伸缩振动的倍频峰以及CH₂变形伸缩振动的合频峰一致。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同特征波长数的RMSE值Fig.3 RMSE values of models with different characteristic wavelength numbers

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 基于SPA算法筛选的特征波长Fig.4 Characteristic wavelengths selected based on SPA algorithm

2.2.2 UVE算法

UVE算法筛选的油页岩TOC特征波长结果如图5所示。 其中竖线左侧为1265个光谱变量的稳定性指数分布曲线, 右侧为UVE产生的同光谱变量相同数量的随机变量稳定性指数分布曲线。 以随机变量稳定性指数绝对值的最大值作为筛选变量的阈值, 即稳定性指数分布曲线在两条水平虚线以外的光谱变量被选中, 共选择出253个特征波长, 所有特征波长点位于光谱吸收峰附近, 其分布如图6所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 各波长变量和随机变量下的稳定性指数Fig.5 Stability index of each wavelength variable and random variable

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 基于UVE算法筛选的特征波长Fig.6 Characteristic wavelengths selected based on UVE algorithm

2.2.3 CARS算法

CARS算法筛选的特征波长数、 RMSE以及回归系数随运行次数的变化如图7所示。 从图中可以看出, 当运行次数从1次增加到24次, 特征波长数的下降由快变慢, RMSECV逐渐降低, 表明在1~24次运行过程中剔除了较多的无关光谱变量, 模型精度不断提高。 随着运行次数的继续增加, RMSECV缓慢或者迅速增大, 回归系数不断变大。 在运行次数为24次时, RMSECV值最低, 此时有65个波长被保留下来, 其中大部分波长位于光谱吸收峰附近, 其分布如图8所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 基于CARS算法筛选特征波长的过程Fig.7 Process of selecting characteristic wavelengths based on CARS algorithm

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 基于CARS算法筛选的特征波长Fig.8 Characteristic wavelengths selected based on CARS algorithm

利用SPA, UVE和CARS进行特征波长筛选后, 特征波长数明显少于全光谱波长的1 265个波长, 分别是全谱波长的1.26%, 20%和5.14%。 说明特征波长筛选对于简化模型、 提高模型效率能表现出较好的效果。 此外从筛选出的波长来看, 有明显吸收峰或者吸收峰附近的波长被保留了下来。

为比较不同建模方法和不同特征波长提取方法对油页岩中有机碳含量的预测效果, 分别采用PLS、 SVM和RF建模方法对全光谱波段、 CARS、 SPA、 UVE筛选出的特征波长建立油页岩TOC含量的预测模型(表3)。 SVM进行建模, 以径向基函数作为模型核函数, 根据网格搜索法优选惩罚因子和核函数参数。 建立RF模型时, 以不同特征波长提取方法下验证集相关系数最高时的决策树数量作为RF模型的最优决策树数目。

由表4可知, 基于CARS, UVE和SPA三种特征波长的模型精度均高于全谱波长模型的精度。 这是因为利用全光谱进行建模时, 包含的变量较多, 变量间存在有冗余信息的干扰, 而特征波长提取可以有效去除冗余信息, 提取后的波长能充分代表原始光谱的有效信息, 从而提高模型质量。 在三种特征波长选择方法中, 基于CARS算法提取特征波长之后所建立的模型效果最好, 尤其是CARS-SVM模型, 其验证集测定系数由未进行特征波长选择时的0.793 0提升到0.906 6, 均方根误差由0.286 8降低到0.222 0, 是所有模型中最优的, 该模型可以应用于油页岩总有机碳含量预测。

表4

Table 4

表4(Table 4)

表4 不同方法建模结果 Table 4 Modeling results of PLS, RF and SVM 建模 方法 特征波段 筛选方法 波段 数 校正集 验证集 Rc2 RMSEC Rv2 RMSEV SVM 全谱波段 1 265 0.867 8 0.260 5 0.793 0 0.286 8 CARS 65 0.912 4 0.204 6 0.906 6 0.222 0 UVE 253 0.882 8 0.222 6 0.806 5 0.282 0 SPA 16 0.883 5 0.244 0 0.815 3 0.278 0 RF 全谱波段 1 265 0.785 9 0.274 6 0.717 4 0.342 7 CARS 65 0.812 5 0.297 6 0.765 9 0.316 9 UVE 253 0.857 2 0.285 1 0.799 3 0.278 0 SPA 16 0.842 5 0.293 7 0.812 7 0.323 5 PLS 全谱波段 1 265 0.765 5 0.314 8 0.742 1 0.342 5 CARS 65 0.803 8 0.322 4 0.761 6 0.280 9 UVE 253 0.834 5 0.284 3 0.754 1 0.310 0 SPA 16 0.774 8 0.316 2 0.762 7 0.320 6

表4 不同方法建模结果 Table 4 Modeling results of PLS, RF and SVM

三种模型中, SVM模型的效果优于RF模型和PLS模型, 说明基于SVM法建立的预测模型能较好地应用于近红外光谱检测模型, 其预测精度高、 实用性强, 对快速准确检测油页岩TOC含量有实际价值^[15]。 此外, 相较于近红外光谱领域常用的线性的PLS模型而言, SVM和RF两种非线性的建模方法得到模型在校正集上的表现均有不同程度的提高, 其中SVM模型无论在校正集还是在验证集均有明显提升。 这是因为油页岩样本中的碳存在于各类烃的中, 由于不同类别含烃基团的吸收峰之间相互影响, 使得油页岩TOC含量和近红外光谱数据之间存在着复杂的非线性关系。 作为线性分析的PLS模型, 其表达光谱数据和浓度数据之间的非线性关系能力不及SVM和RF, 后者能够有效地对自变量和因变量之间的非线性关系进行描述, 从而增强了光谱数据和样本理化值之间的相关性, 使得所建立的非线性模型效果要略优于线性的PLS模型。 该结果与陈华舟^[4]等基于近红外光谱针对鱼粉蛋白进行定量分析所建立的非线性模型SVM和RF效果优于线性模型PLS的效果结论一致。 建立非线性模型比线性模型的R²更好, 准确度更高^[16]。 因此, 可以通过非线性建模近红外光谱法实现油页岩总有机碳含量快速检测。