设初始数据集记作T=(O, F, C), 其中O={o₁, o₂, …, o_n}为近红外光谱样本数据的集合, o_i表示某含水量岩石的近红外光谱, F={f₁, f₂, …, f_m}为初选的近红外光谱特征矩阵, 存储提取的特征变量, f_k为某特征变量, C={c₁, c₂, …, c_l}为岩石含水量级别的集合, c_j表示某含水量级别(某一小区间)。

特征选择过程由产生初始特征子集、 搜索策略、 评价函数、 终止条件四个步骤组成。 本工作采用互信息作为相关程度的评价标准, 对特征变量与含水量级别之间依赖关系的强弱程度进行度量, 并根据其值按降序排序, 然后, 选择前k个特征(或一致性因子等其他终止条件)组成选择子集S。

信息度量法具有无参、 非线性的优势, 且不需要预先知道样本数据的分布, 在特征选择算法中得到广泛应用。

1.2.1 BIF法

BIF(best individual feature)法^[7]直接利用互信息作为评价标准, 是一种简单而直观的信息度量方法, 利用该方法度量含水量级别C和某特征变量f的相关性, 其评价标准可表示为

I(C; F)=∑c∈Cf∈Fp(c, f)logp(c, f)p(c)p(f)(1)

其中, p(c)和p(f)分别为含水量级别c和某特征值f的边缘概率分布; p(c, f)为含水量级别C与某特征值f的联合概率分布。

由此定义式可知, 当含水量级别C和某特征值f完全无关或相互独立时, 互信息I(C; f)=0; 二者间相互依赖程度越高, 互信息I(C; f)的值就越大。 可见, BIF选择算法计算步骤首先按式(1)分别计算特征矩阵F={f₁, f₂, …, f_m}中的特征f与含水量级别C的互信息I(C; f), 然后根据互信息值大小按降序顺序排序, 最后选择前k个特征组成近红外光谱特征集S, 完成特征选择。

1.2.2 MIC法

最大信息系数(maximal information coefficient, MIC)由Resher^[8]于2011年首次提出, 他利用互信息定义了两个变量之间的最大信息系数, 并利用该系数衡量两个变量之间的相关性。 利用该方法度量含水量级别C和特征变量f相关性的具体步骤是: 首先对含水量级别C与特征变量f的散点图进行网格划分, 在固定行数x_i、 列数y_i情况下, 有多种划分方式, 对应着不同的网格 Gxiyi(j)。 然后, 用各散点在网格的子格内的频率来代替变量的概率, 计算含水量级别C与特征变量f之间的互信息, 并取其最大值, 记为最大互信息

Imax(C, f, xi, yi)=maxI((C, f)|Gxiyi(j)), j=1, 2, …(2)

其中, Imax(C, f, xi, yi)表示在固定行列数x_i与y_i情况下, 不同的网格划分方式下的互信息最大值(最大互信息); I((C, f)|Gxiyi(j))表示含水量级别C与特征变量f的散点图在网格 Gxiyi(j)下的互信息。

最后, 对于任意的行数x_i、 列数y_i, i=1, 2, …, 由式(2), 计算最大信息系数

MIC=maxImax(C, f, xi, yi)logmin(xi, yi)(3)

其中x_i, y_i< B(n), B(n)为网格分割细度^[9]。

以最大信息系数为评价标准时, 每个特征变量的MIC值视为其权重, MIC值越大, 表明该特征变量f与含水量级别C之间的相关性越大。 最后选择满足预先给定的阈值的特征变量, 组成近红外光谱特征集S, 完成特征的选择。

砾岩岩样采自敦煌莫高窟北区崖壁, 基本信息如表1所示, 其中干燥后质量为砾岩样品在设定温度为105~110 ℃真空干燥箱里烘干24 h后冷却至室温所测得质量。 为了进一步了解砾岩所含矿物成分, 采用日本理学电机公司(Rigaku)生产的D/MAX2500 X射线衍射仪对砾岩样品进行X射线衍射实验, 其矿物成分分析结果见表2所示。 因为砾岩结构疏松, 呈半胶结状态, 遇水易崩解, 不易加工成标准试件, 将其加工成尺寸约为80 mm×90 mm×60 mm的不规则形状[图2(a)], 以供实验所需。

表1

Table 1

表1(Table 1)

表1 试样基本信息表 Table 1 Basic information of the conglomerate samples 岩性 外貌描述 颜色 干燥后质量/g 尺寸/mm 砾岩 疏松多孔, 块状结构 青灰色 922.60 约为80×90×60

表1 试样基本信息表 Table 1 Basic information of the conglomerate samples

表2

Table 2

表2(Table 2)

表2 砾岩矿物成分分析结果 Table 2 Mineral composition of the conglomerate samples 岩性 石英/% 钾长石/% 钠长石/% 方解石/% 白云石/% 角闪石/% 石盐/% 黏土矿物总量/% 砾岩 30.4 2.2 23.5 11.8 9.7 2.0 6.9 13.5

表2 砾岩矿物成分分析结果 Table 2 Mineral composition of the conglomerate samples

利用中国矿业大学(北京)深部岩土力学与地下工程国家重点实验室何满潮^[10]自主研发的“ 深部软岩水理作用智能测试系统” 实现岩样吸水模拟。 该系统主要由主体实验箱、 称重系统和数据采集三部分组成, 同时采用由瑞士万通(Metrohm)生产的NIRS XDS OptiProbe Analyzer(近红外光谱分析仪)采集岩样的近红外光谱信息(实验参数见表3), 该近红外仪器由光谱分析仪、 电控平移操作台以及计算机组成, 其中光谱分析仪可精确快速地捕捉波长在400~2 500 nm范围的光谱信号, 光纤探头可采集样品固定位置处的反射光谱信号; 电控平移操作台可实现0~1 500 mm行程内的精确定位, 确保了不同岩石样品之间或同一试样不同时刻取样点位置的一致性。 实验测试装置和测试原理如图1所示。

表3

Table 3

表3(Table 3)

表3 近红外光谱仪采集系统的实验参数 Table 3 Experimental parameters of near infrared spectrometer acquisition system 探头类型 反射探头 采样方式 漫反射-固体 光谱范围 400~2 500 nm 测样角度 90° 测试方式 直接与样品接触

表3 近红外光谱仪采集系统的实验参数 Table 3 Experimental parameters of near infrared spectrometer acquisition system

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 (a)吸水实验系统及近红外测试装置; (b)实验原理示意图Fig.1 (a) Water sorption test system and near-infrared testing device; (b) schematic diagram of the experimental setup

考虑到岩石材料的不均匀性, 沿着砾岩高度方向分别采集了样品3个不同位置(图2所示红点位置), 共51条近红外光谱信息, 这些光谱涵盖了砾岩3个测试点从初始干燥状态到饱和状态的全过程, 表征了砾岩材料中O—H基团数量变化与光谱特征参数之间的相关程度。 将干燥后的砾岩岩样放置于图1所示位置, 岩样开始吸水。 在岩样吸水过程中, 采用近红外光谱分析仪采集岩样的近红外光谱。 测试时将光纤探头分别接触岩样的3个测量点, 自下而上依次测量, 测量的频率随岩石的吸水速度适时调整。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 砾岩吸水过程可见光图 (a): 409 s; (b): 558 s; (c): 704 s; (d): 775 s; (e): 1 224 sFig.2 Visible light map of the conglomerate water absorption process (a): 409 s; (b): 558 s; (c): 704 s; (d): 775 s; (e): 1 224 s

整个砾岩样品历时1.15 h从干燥状态达到饱和状态。 在1, 2和3号点分别采集近红外光谱曲线17条、 18条、 16条, 分别记作O1-1, O1-2, …, O1-17; O2-1, O2-2, …, O2-18; O3-1, O3-2, …, O1-16, 共51条近红外光谱。 不同采集点位置达到吸水饱和状态的时间是不一样的, 整个实验过程3个采集点位置从干燥到饱和状态所用的时间分别为2 998, 3 048和3 099 s。

3.1.1 光谱预处理

通过对比常用预处理方法的处理效果, 选择多点平滑与基线校正相结合的方法(NPS+B-corr)对原始近红外光谱进行预处理, 既去除了无效噪音, 又校正了基线漂移; 限于篇幅, 文中仅列出砾岩1号点处的近红外光谱图, 如图3所示, 其中(a)为砾岩原始近红外光谱图, (b)为经过NPS+B-corr法预处理得到的砾岩近红外光谱图。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 砾岩1号点处近红外光谱 (a): 原始光谱; (b): 经NPS+B-corr预处理后光谱Fig.3 Near-infrared spectrum of conglomerate point 1 (a): Original spectrum; (b): NPS+B-corr pre-processed spectrum

3.1.2 谱段选择

由图3可知, 在400~2 500 nm波长范围内分别在1 450, 1 930和2 300 nm附近有明显的吸收峰, 依次命名为峰R₁、 峰R₂、 峰R₃; 光谱的吸光度随岩石含水量变化明显。

分析1号点光谱曲线(图3)可知, R₁和R₂两个吸收峰的波峰随含水量的增大而增高, 峰顶中心位置逐渐右移, 当含水量达2.467%之后峰形变化越来越小, 峰值波动减弱, 最后峰R₁的中心点位于1 450 nm左右, 峰R₂中心点位于1 930 nm左右。 随着含水量的增加, 峰R₁、 峰R₂的光谱信号越来越强, 峰R₃的信号则不断减弱, 加上临近噪音波段信号干扰, 峰R₃的特征信号越来越难提取。 综上分析, 本文主要分析含水岩石在峰R₁(1 450 nm)、 峰R₂(1 930 nm)谱段处的光谱特征。

3.2.1 初始特征集F的建立

构建特征集是建立岩石含水量与近红外光谱之间的数学模型中至关重要的一环, 它是后续特征选择的基础, 决定了特征选择的范围和质量。

分析经NPS+B-corr法预处理后的近红外光谱图的峰R₁、 峰R₂特点, 提取下列特征变量: 峰高(Height), 半高宽(FWHM), 峰面积(Area), 左肩宽度(Left Half Width, LHW), 右肩宽度(Right Half Width, RHW), 左右肩比(LHW/RHW)共6个初始特征变量, 其几何示意图见图4, 因此设定的初始特征集F={f₁, f₂, f₃, f₄, f₅, f₆}={Height, Area, FWHM, LHW, RHW, LHW/RHW}, 其具体数值如表4所示。 限于篇幅, 此处只列出1 450 nm谱段处部分数据。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 初始特征变量示意图Fig.4 Schematic diagram of initial characteristic variables

3.2.2 特征变量归一化

由表4可以看出6个初始特征变量的量纲不同, 且变量之间的变化幅度也不同, 因此需要对原始数据进行归一化处理, 来消除量纲和变化幅度不同带来的影响。

表4

Table 4

表4(Table 4)

表4 1 450 nm谱段近红外光谱的初始特征变量(只列出部分) Table 4 Initial characteristic variables of near-infrared spectra at 1 450 nm (only part of data listed) 特征值F f1 /nm f2 /nm f3 /nm f4 /nm f5 /nm f6 O1-1 0.006 36.702 0.294 17.641 19.061 0.925 O1-2 0.014 103.242 1.628 21.971 81.271 0.270 O1-3 0.071 118.620 9.1463 52.161 66.460 0.784 8 ︙ ︙ ︙ ︙ ︙ ︙ ︙ O3-16 0.157 119.694 20.562 53.409 66.285 0.806 均值 0.079 110.146 10.439 47.652 62.494 0.770 标准方差 0.041 26.540 5.512 12.990 14.850 0.120 最大值 0.177 123.649 23.952 54.434 81.271 0.938 最小值 0.005 36.032 0.279 16.305 18.589 0.270

表4 1 450 nm谱段近红外光谱的初始特征变量(只列出部分) Table 4 Initial characteristic variables of near-infrared spectra at 1 450 nm (only part of data listed)

归一化的方法是将原式数据矩阵的各元素减去该元素所在列的最小值后再除以该列元素的极差(详见本刊40卷3期971页), 有

x'ij=xij-min(xj)max(xj)-min(xj)(4)

归一化结果如表5所示。

表5

Table 5

表5(Table 5)

表5 1 450 nm谱段近红外光谱的初始特征变量归一化值(只列出部分) Table 5 Normalized initial characteristic variables of near-infrared spectra at 1 450 nm (only part of data listed) 特征值F f1/nm f2/nm f3/nm f4/nm f5/nm f6 O1-1 0.005 8 0.007 6 0.000 6 0.035 0 0.007 5 0.980 5 O1-2 0.052 3 0.767 1 0.057 0 0.148 6 1.000 0 0.000 0 O1-3 0.383 7 0.942 6 0.374 6 0.940 4 0.763 7 0.770 7 ︙ ︙ ︙ ︙ ︙ ︙ ︙ O3-16 0.883 7 0.954 9 0.856 8 0.973 1 0.760 9 0.802 4 均值 0.432 6 0.845 9 0.429 2 0.822 1 0.700 4 0.747 9 标准方差 0.040 6 26.540 2 5.511 8 12.989 6 14.850 3 0.120 3 最大值 1.000 0 1.000 0 1.000 0 1.000 0 1.000 0 1.000 0 最小值 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0

表5 1 450 nm谱段近红外光谱的初始特征变量归一化值(只列出部分) Table 5 Normalized initial characteristic variables of near-infrared spectra at 1 450 nm (only part of data listed)

3.2.3 含水量级别C的确定

含水量级别C={c₁, c₂, …, c_l}为岩石含水量的级别集合, c_j表示某含水量级别。 实验中砾岩不断动态吸水, 根据总吸水量与测量时得到的近红外光谱对应关系, 绘制各特征变量随含水量变化曲线, 这里仅列出峰高随含水量的变化曲线, 如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 含水量-峰高曲线及含水量级别分区 (a): 1号点; (b): 2号点; (c): 3号点Fig.5 Peak height of water content and water content level partition (a): Point 1; (b): Point 2; (c): Point 3

分析图5可知, 含水量与各光谱特征参数的关系曲线有相似的变化趋势, 曲线大致呈“ S” 型, 可将其划分为三个阶段: 低含水量时特征参数平稳变化, 中间段特征参数迅速增长, 高含水量时特征参数基本稳定。 当高含水量达到一定程度时(图5中各点含水量超过5%), 特征参数有所增加, 这一现象可以解释为: 在岩石自下而上动态吸水的过程中, 虽然监测点的含水量已经达到饱和, 但水分仍然在持续不断地路经监测点, 向上吸附至试件上端未饱和区域, 由此, 最终在测试点表面形成一层水分薄膜, 使得该处岩石表面具有较高反射率, 对应的光谱特征参数出现进一步增加的现象。

结合吸水过程分析可知, 这三个阶段水分尚未到达测试点, 水分接触到测试点并开始浸润该位置, 最终水分使测试点区域内岩石吸水饱和。 测试点吸水过程是短暂迅速的, 利用近红外光谱分析技术可以很精确的对岩石上各点的含水量大小进行监测识别。 可见, 测试点的光谱特征随含水量变化规律, 有着非常明显的三个阶段, 对应着各测试点的三个含水量级别, 据此可设置含水量级别C: C={c₁, c₂, …, c_l}={干燥状态, 吸水状态, 饱和状态}, 该级别属于文本类别, 为了计算最大相关系数MIC评价变量间相关程度, 将其定义为: C={c₁, c₂, c₃}={干燥状态, 吸水状态, 饱和状态}={1, 2, 3}。

3.2.4 特征选择计算与结果分析

在进行特征选择之前, 首先要进行初始特征集各特征变量间相关性筛选, 以便去掉冗余特征, 根据文献[11]中计算方法, 进行特征变量筛选。 将初始特征集F简化为S={f₁, f₄, f₅}, 即由峰高(Height), 左肩宽度(LHW), 右肩宽度(RHW)三个特征变量构成特征集。

(1)BIF法的特征选择

利用式(1)计算特征S={f₁, f₄, f₅}与含水量级别C={c₁, c₂, c₃}之间的互信息I(C, f)。 首先根据3.2.2小节中光谱特征变量归一化之后的计算结果, 将各特征变量值按照[max, min]进行分组, 并计算每组所得的概率, 然后再根据3.2.3小节含水量级别C的分类结果分别计算每一级别所占的概率, 最后根据各特征变量和含水量级别的分组结果计算两者的联合分布概率, 再代入到1.2.1小节中的公式, 具体计算结果见表6。

表6

Table 6

表6(Table 6)

表6 特征变量与含水量级别的互信息值 Table 6 Mutual information values of characteristic variables and water content levels 吸收峰 I(C, Height) I(C, LHW) I(C, RHW) 峰R1 0.618 8 0.553 8 0.675 3 峰R2 0.584 0 0.529 8 0.530 6

表6 特征变量与含水量级别的互信息值 Table 6 Mutual information values of characteristic variables and water content levels

特征变量与含水量级别之间的互信息值越大, 表示二者之间的相关性越大, 分析表6可知, 按互信息I(C, f)由大到小排序, 则峰R1的特征变量排序结果为I(C, RHW)> I(C, Height)> I(C, LHW), 可知相关性由高到低为右肩宽度, 峰高, 左肩宽度。 同理峰R₂的特征变量排序结果为I(C, Height)> I(C, RHW)> I(C, LHW), 相关性由高到低为峰高, 右肩宽度, 左肩宽度。

根据BIF法得出特征变量的重要性排序, 据此确定最终的近红外光谱特征选择, 可以选择相关性最大的特征, 也可以选择相关性较大的多个特征的组合作为最优特征子集。

(2)MIC法的特征选择

基于第1.2.2节中的计算方法, 首先将特征f与含水量级别C={干燥状态, 吸水状态, 饱和状态}={1, 2, 3}组成散点图D, 以峰R₁为例, 各特征变量与含水量级别C之间的散点图D_i如图6所示, 然后依次计算不同网格划分方式 Gxiyi(j)下最大互信息值I_max(C, f, x_i, y_i), 接着组成最大互信息标准化矩阵 Imax(C, f, xi, yi)logmin(xi, yi), 再取这个矩阵的最大值, 确定最大相关系数MIC=max Imax(C, f, xi, yi)logmin(xi, yi), 计算结果如表7所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 含水量级别C-特征f散点图 (a): Height-C散点图; (b): LHW-C散点图; (c): RHW-C散点图Fig.6 Water content level C-feature f scatter plot (a): Height-C; (b): LHW-C; (c): RHW-C

表7

Table 7

表7(Table 7)

表7 特征变量与含水量级别间的MIC值 Table 7 MIC values for correlation between the characteristic variables and water content levels 吸收峰 MIC(C, Height) MIC(C, LHW) MIC(C, RHW) 峰R1 0.692 3 0.787 1 0.609 5 峰R2 0.716 3 0.787 1 0.370 6

表7 特征变量与含水量级别间的MIC值 Table 7 MIC values for correlation between the characteristic variables and water content levels

由表7可知, 峰R₁各特征变量f与含水量级别C之间的MIC值, 由大到小排序为MIC(C, LHW)> MIC(C, Height)> MIC(C, RHW), 其中任一特征值对应的MIC值均> 0.5, 说明通过峰分析提取的特征变量, 随着岩石含水量的变化, 均会发生一定程度的变化, 特征值f与岩样吸水状态C之间存在着一定的相关性, 特别地LHW的MIC值接近于0.8, 说明其与岩石含水量存在较强的相关性, 可作为有效表征岩石含水量的特征变量。 综上, 根据MIC法得到特征变量排序依次为左肩宽度, 峰高, 右肩宽度。

同理, 对峰R₂的特征值f进行最优选择, MIC值由大到小排序为MIC(C, LHW)> MIC(C, Height)> MIC(C, RHW), 除特征变量右肩宽RHW外, 任一特征变量对应的MIC值均> 0.5, 其中, LHW的MIC值接近于0.8, 表明其与岩石含水量有较强的相关性。 综上, 依MIC值为权重进行特征值重要性排序由高到低依次为左肩宽度, 峰高, 右肩宽度。

根据MIC法得出特征变量的重要性排序及其相关权重, 据此可以选择相关性最大的特征变量, 也可以选择相关性较强的多个特征变量的组合作为最有特征子集。

3.2.5 BIF与MIC方法性能评估

为了评估BIF和MIC法的性能, 我们使用决策树来评估所选特征的有效性。 根据所采集到的3个点共51条近红外光谱特征, 利用BIF法和MIC法分别计算峰R₁(1 450 nm)和峰R₂(1 930 nm)处6个特征变量的平均值, 为了计算方便, 我们对含水量取插值平均得到决策树中的k值, 通过软件matlab导入决策树数据流, 计算各级别的平均分类准确度。 从图7看出, MIC的准确度较高。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 峰R1 (a)和峰R2 (b)在决策树分类器下的准确率Fig.7 The accuracies of the decision tree classifier for peak R1 (a) and peak R2 (b)

基于互信息的特征选择算法的评价标准是属于信息论的相关度量, 主要用于描述两个随机变量之间相互依存关系的强弱。 采用BIF方法和MIC方法对含水岩石近红外光谱进行特征选择, 均得到了近红外光谱特征参数与岩石含水量级别之间的相关程度强弱, 且各个特征参数具有明显几何意义, 通过BIF法与MIC法的特征选择结果, 后期在利用近红外光谱特征建立砾岩含水量检测模型中, 可以将峰高、 左肩宽、 右肩宽三个特征参数直接参与模型的构建中。 总体来说互信息方法思路简单, 概念清晰, 计算过程简单, 利用此方法获得的与含水岩石水分信息最相关的最优光谱特征, 可以为建立更稳定、 精确、 高效的预测模型做基础。

此外两种方法在评价含水岩石近红外光谱特征变量与含水量级别相关程度强弱的结果并不完全一致, 而且MIC方法比BIF方法对含水量级别的识别精度更高些。 其实两种方法都是通过计算变量间的概率来获得最优相关子集, 所不同的是BIF法直接计算两个变量的边缘概率和联合概率分布, 而MIC方法先对变量间通过不同网格划分的形式, 用变量在网格内出现的频率来代替概率, 然后用不同网格划分方式下的互信息最大值计算最大相关系数, MIC方法所获得的样本数量更为精确, 进一步考虑了所选变量间的相关性, 因此计算结果更为准确。 本文目前只针对砾岩一种岩性进行研究, 对于MIC方法的准确性还需要在更多的样本集上来进行验证, 同时基于上述两种特征选择算法, 得出特征变量重要性排序, 关于最优特征子集中变量类型与个数的确定, 需要后期构建预测模型的精度来佐证, 是一个反复验证与选择的过程。