岩性不同的含水岩石近红外光谱特征对比分析

张 芳, 户佐乐, 王东升，刘雨濛, 谢运鑫, 卓慧慧, 何满潮

1. 中国矿业大学(北京)深部岩土力学与地下工程国家重点实验室，北京 100083 2. 中国矿业大学(北京)力学与建筑工程学院，北京 100083

引 言

近红外光谱分析(near infrared spectroscopy technique, NIRS)是近几十年来发展最快，最引人注目的光谱分析技术之一，被广泛应用于食品、化工、农业等领域[1]。在岩土工程领域应用近红外光谱分析检测、监测含水岩石水分状态，是近些年发展起来的新思路，其原理是根据岩土介质O—H基团的吸收强度，建立含水量与近红外光谱特征之间的定量关系，实现预测岩土介质中含水量的目的。与传统分析方法相比，该方法具有无损、实时、定量的优势。

对于含水量与光谱特征之间建模研究，诸多学者进行了大量的研究工作，尤其在土壤方面研究成果颇多：金慧凝等[2]通过提取反射光谱特征指标，定量分析土壤含水量与反射光谱特征之间关系，建立了土壤水分含量光谱预测模型。包青岭等[3]利用包络线消除法提取反射光谱水分吸收特征，并与土壤含水量进行相关性分析，通过随机森林方法对光谱水分吸收特征参数进行分类，获取各参数对土壤含水量的重要性，运用多元逐步回归方法，建立土壤水分含量反演模型。娄径等[4]通过对光谱数据进行倒数、对数、均方根及一阶导数微分等光谱变换，分析光谱特征，并与土壤含水量进行相关分析，利用多元线性回归分析建立土壤含水量监测模型。

可见，光谱特征选择是构建准确、稳健的定量模型的关键。国内学者对此进行了很多研究，如：孔清清等[5]基于随机森林，结合博弈论进行了近红外光谱特征选择，利用互信息选择出无冗余的特征子集，此法应用于近红外光谱分类中有较高的分类识别率。

利用信息度量法进行波谱的特征选择具有无参、非线性的优势，能有效度量两随机变量之间相关性，能较好地解决特征变量选择问题，该方法在特征选择算法中得到广泛应用。

但是，上述研究成果大多集中于研究土壤含水量与光谱特征之间的相关性，针对岩石的研究工作很少，更没有涉及岩性对含水岩石光谱特征选择的影响，即不同岩性，能否选择相同特征集，或是需要具体问题具体分析，针对不同现场选择对应的特征集合，这个问题对生产实践中构建模型时选择特征变量具有非常重要的指导意义。

故此，利用互信息作为相关程度的度量标准，对比分析不同岩性的含水岩石近红外光谱的特征选择结果，以期评价岩性对含水岩石光谱的影响。

1 基于互信息的特征选择

采用最大信息系数(maximal information coefficient，MIC)[6]进行特征选择，其计算思想是：设任意的行数xi列数yi下的含水量C与特征参数f的散点图网格Gxiyi(j)，计算最大互信息和最大信息系数

Imax(C,f,xi,yi)=maxI((C,f)|Gxiyi(j))，

j=1, 2, …

(1)

(2)

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

最后选择满足预先给定的阈值的特征参数，组成近红外光谱特征集S，完成特征的选择。

2 实验部分

1.1 样品制备

砾岩、粉砂岩采自敦煌莫高窟北区的崖壁，样品信息如表1、表2所示。因为该砾岩呈半胶结状态，不易加工成标准试件，实验时将砾岩加工成尺寸约为80 mm×90 mm×60 mm的不规则形状(图1)，使之能够满足实验仪器的放置要求。

表1 样品基本信息Table 1 Sample basic information

表2 矿物成分信息Table 2 Mineral composition information

图1 岩样及其测试点位置(红点)(a)：砾岩；(b)：粉砂岩；(c)：夯土

夯土试样(表1和表2)采自敦煌莫高窟108洞室，该夯土强度较低，无法直接用钻机取出土样，所以先取土块，然后加工成长方体土样，之后再采用打磨的方法制成φ50mm的标准样品(图1)。

实验中，试样不断吸水，水分自下而上运移，故沿着试样高度方向，选取3个近红外光谱测试点，并尽量避开签字笔标记区域，三种岩性样本的测试点位置如图1中红点位置。

1.2 原理及设备

将岩样放置真空干燥箱内，设置箱内温度105～110 ℃，干燥24 h，取出干燥后岩样冷却12 h称重，利用中国矿业大学(北京)深部岩土力学与地下工程国家重点实验室何满潮[8]自主研发的“深部软岩水理作用智能测试系统”，进行岩样吸水模拟室内实验。该系统主要由主体实验箱、称重系统和数据采集系统三部分组成，如图2所示。

1.3 光谱采集

在试样吸水过程中，采用瑞士万通的XDS SmartProbe近红外光谱分析仪采集不同位置不同时刻的近红外光谱。测试时将光纤探头分别接触试样的3个测量点(图1红点)，自下往上依次测量，测量的频率随试样的吸水速率适时调整。实验参数如表3，实验装置如图3。

图2 岩石吸水过程中的近红外光谱采集实验设备

表3 近红外光谱分析仪的实验参数[9]Table 3 Experimental parameters of the near infrared spectroscopy analyzer[9]

图3 XDS SmartProbe 近红外光谱分析仪

整个实验，在砾岩从干燥到饱和的吸水全过程中，共采集了51条近红外光谱信息，分别为：1号点17条，2号点18条，3号点16条。

在粉砂岩吸水全过程中，共采集了106条近红外光谱信息，分别为：1号点51条，2号点34条，3号点21条。

在夯土吸水全过程中，共采集了149条近红外光谱信息，分别为：1号点24条，2号点59条，3号点66条。

1.4 近红外光谱预处理

利用XDS SmartProbe近红外光谱分析仪配套软件提供的一阶导数法对采集的光谱进行预处理，消除背景的常数平移对近红外光谱的影响，使数据具有更好的连续性，处理前后的光谱如图4—图6所示，限于篇幅，仅列出砾岩、粉砂岩、夯土1号点的光谱图。

图4 砾岩1号点近红外光谱(a)：原始光谱；(b)：一阶导数预处理后光谱

1.5 谱段选择和初始特征变量

分析图4—图6可知，在400～2 500 nm波长范围内有3个明显的吸收峰，分别在1 400，1 900和2 300 nm附近，其光谱反射率随试样含水量变化而变化，依次将其命名为峰R1、峰R2、峰R3。随着含水量的不断增大，R1和R2两个吸收峰的波峰越来越高，峰R1中心点位置最终停留在1 400 nm左右，峰R2中心点位置最终停留在1 900 nm左右，而R3吸收峰的波峰随含水量增加逐渐减小，信号特征逐渐减弱，因2 400 nm之后的噪音干扰强烈，故峰R3不适合作为含水量信息的特征谱段。因此，选择峰R1、峰R2所在的1 400和1 900 nm谱段进行含水试样光谱特征分析，具体提取的特征变量如图7所示，分别为峰面积(Area)、峰高(Height)、半高宽(FWHM)、左肩宽(left half width，LHW)、右肩宽(right half width，RHW)、左右肩宽比(LHW/RHW)共计6个初始特征参数，设定初始特征集F为F={f1,f2,f3,f4,f5,f6}={Area，Height，FWHM，LHW，RHW，LHW/RHW}, 各含水试样近红外光谱的初始特征数值如表4所示。

需要特别强调的是，在整个实验历程中，粉砂岩3号点处采集的21条近红外光谱没有明显的吸收峰，分析其原因。3号点位于粉砂岩顶端，当水分没有达到这个位置并浸润到它时，该点始终处于干燥状态，所以没有吸收峰。3号点采集的21条光谱都没有采集到含水情况下的光谱，为失效光谱。故在后续分析中将该组实验数据去掉。

图5 粉砂岩1号点近红外光谱(a)：原始光谱；(b)：一阶导数预处理后光谱

图6 夯土1号点近红外光谱(a)：原始光谱；(b)：一阶导数预处理后光谱

1.6 特征变量归一化

分析表4可知，由于6个初始特征变量的量纲不同，且特征变量之间的变化幅度不同，可能导致在分析计算过程中，一些数量级较小的特征变量的作用无法体现，因此对上述表中的原始数据进行归一化变换，将所有变量转换成0-1内的数值，消除量纲和变化幅度不同带来的影响。

归一化的方法是将原始数据矩阵的各元素减去该元素所在列的最小值后再除以该列元素的极差，公式如下

图7 近红外光谱的初始特征变量几何意义示意图[9]

Fig.7Schematicdiagramofgeometricmeaningofinitialcharacteristicvariablesofnear-infraredspectroscopy[9]

表4 砾岩在峰R1处的初始特征变量(只列出部分)Table 4 Initial characteristic variables of conglomerate at the peak R1(only some data shown in the table)

续表4

O1-60.858 170.023 0033.442 9915.644 1517.798 840.878 943.097O1-70.707 910.019 5232.949 7315.486 0517.463 680.886 763.408O1-80.757 880.020 1733.442 2416.000 7517.441 490.917 403.733O1-90.985 460.025 9434.016 3716.231 6117.784 760.912 673.984︙︙︙︙︙︙︙︙O3-161.492 560.039 0034.051 8916.336 9917.714 900.922 224.971

归一化结果如表5所示。

表5 砾岩在峰R1处初始特征变量归一化值(只列出部分)

Table5NormalizedvaluesofinitialcharacteristicvariablesofconglomerateatthepeakR1(partiallisting)

近红外光谱特征值Ff1f2f3f4f5f6O1-10.013 460.006 320.190 050.305 060.170 430.669 85O1-20.065 190.029 740.452 660.564 390.409 940.469 75O1-30.423 740.166 490.780 300.845 590.732 680.231 88O1-40.458 730.182 620.780 940.837 510.738 180.218 54O1-50.673 970.251 630.914 860.982 620.853 070.187 38O1-60.525 480.199 530.864 860.869 410.843 450.115 41O1-70.427 240.165 710.829 170.839 610.807 830.134 73O1-80.459 910.172 030.864 800.936 620.805 470.210 43O1-90.608 710.228 110.906 340.980 140.841 950.198 75︙︙︙︙︙︙︙O3-160.940 270.355 040.908 911.000 000.834 530.222 35均值0.443 990.189 430.700 870.749 840.669 980.267 86标准方差0.249 140.148 660.269 650.251 800.257 850.216 78最大值1.000 001.000 001.000 001.000 001.000 001.000 00最小值0.000 000.000 000.000 000.000 000.000 000.000 00

2 结果与讨论

在进行特征选择之前，需要对初始特征集各特征变量之间、特征变量与含水量之间的相关性进行筛选，以便去掉冗余特征，本文参照文献[13]中采用的阈值及结论，取初始特征变量之间的相关系数的阈值为0.95，初始特征变量与含水量之间的相关系数的阈值为0.5。则利用相关系数评价上述变量间的相关程度，选取的特征变量如表6。

表6各试样在峰R1，峰R2处满足相关系数阈值要求的特征变量

Table6CharacteristicvariablessatisfyingthecorrelationcoefficientthresholdatpeaksR1andpeaksR2foreachsample

岩性特征变量R1R2砾岩f1, f5f1, f4粉砂岩f2, f4, f5f1, f5夯土f2, f5f1, f4, f5

2.1 砾岩的特征选择分析

将砾岩特征变量与含水量数据做成散点图，如图8所示。

利用第1节中的公式，分别计算峰R1和R2处的f1，f5,f1，f4与含水量C之间的MIC值，如表7。

表7 砾岩特征变量与含水量间的MIC值Table 7 MIC values between characteristic variables and water content of conglomerate

由表7可知，砾岩在峰R1处有MIC(C,f1)>MIC(C,f5)，说明峰面积与含水量相关关系最强，右肩宽次之；对于峰R2有MIC(C,f1)>MIC(C,f4)，说明峰面积与含水量相关关系最强，左肩宽次之。表7中的MIC值位于0.4～0.55之间，说明特征变量与含水量之间相关性偏弱，这与图8的散点图的规律相一致。究其原因，砾岩形状不规则，因而导致体积计算不准确，含水量计算误差较大，数据规律分散，表现出特征变量之间的相关程度不强，这组数据真正反映的近红外光谱特征与含水量的相关性不具有代表性意义。

图8 砾岩特征变量与含水量散点图(归一化)(a)：峰R1；(b)：峰R2

2.2 粉砂岩的特征选择分析

将粉砂岩特征变量与含水量数据做成散点图，如图9所示。

图9 粉砂岩特征变量与含水量散点图(归一化)(a): 峰R1; (b): 峰R2

利用第1节中的公式，分别计算峰R1和R2处的f2,f4,f5,f1,f5与含水量之间的MIC值，如表8。

由表8可知，粉砂岩在峰R1处有MIC(C, f2)>MIC(C, f5)>MIC(C, f4)，说明峰高与含水量相关性最强，其次是右肩宽，最后是左肩宽。对于峰R2有MIC(C, f5)>MIC(C, f1)，说明右肩宽与含水量相关关系最强，峰面积次之。表8中的MIC值位于0.48～0.90之间，说明特征变量与含水量之间相关程度较强，这与图9的散点图表现出来的规律相一致。

表8 粉砂岩特征变量与含水量间的MIC值Table 8 MIC values between characteristic variables and water content of siltstone

2.3 夯土的特征选择分析

将夯土特征变量与含水量数据做成散点图，如图10所示。

利用第1节中的公式，分别计算峰R1和R2处的f2，f5,f1，f4，f5与含水量之间的MIC值，如表9。

图10 夯土特征变量与含水量散点图(归一化)(a)：峰R1；(b)：峰R2

表9 夯土特征变量与含水量间的MIC值Table 9 MIC values between characteristic variables and water content of rammed soil

由表9可知，夯土在峰R1处有MIC(C,f5)>MIC(C,f2)，说明右肩宽与含水量相关性最强，峰高次之。对于峰R2有MIC(C,f5)>MIC(C,f1)>MIC(C,f4)，说明右肩宽与含水量相关性最强，峰面积次之，最后是左肩宽。表9中的MIC值位于0.62～0.95之间，说明特征变量与含水量之间相关程度较强，这与图10的散点图表现出来的规律相一致。

2.4 不同岩性近红外光谱特征选择对比分析

以最大相关系数MIC值作为指标，评价岩石近红外光谱的特征变量与其含水量之间的相关性，对于砾岩、粉砂岩、夯土各特征变量按照相关性大小排序结果具体见表10。

由表10中可知，对于粉砂岩和夯土，若只取两个特征变量，峰R2处只需要关注右肩宽和峰面积即可。峰R1处只需要关注峰高和右肩宽即可，二者的特征变量选择一致。因砾岩形状不规则，计算体积时产生误差较大，不予以考虑。

综上所述，对于含水岩土介质的近红外光谱的特征选择，在1 400 nm附近可选择峰高、右肩宽作为特征变量，参与到含水量预测模型的构建中；在1 900 nm附近可选择峰面积和右肩宽作为特征变量，参与含水量预测模型的构建。

表10 近红外光谱特征变量相关程度排序Table 10 Near-infrared spectral feature variables correlation degree ranking

3 结 论

(1)砾岩、粉砂岩、夯土三种含水岩石的近红外光谱在1 400和1 900 nm附近都有着明显的吸收峰，且随着含水量的变化，吸收强度越来越强，与含水量大小有明显的相关性，故波段1 400和1 900 nm附近可作为分析光谱特征的基本谱段。

(2)对比砾岩、粉砂岩、夯土近红外光谱特征变量与其含水量的最大相关系数MIC值，表明夯土的近红外光谱与含水量之间的相关性最强。

(3)不同岩性的近红外光谱各个特征值与含水量的相关程度不同，具体表现为在1 400 nm附近峰高、右肩宽与含水量都具有较高的相关系数，只是相关性大小会因岩性不同而不同；1 900 nm附近的右肩宽和峰面积与含水量都具有较高相关系数，且右肩宽的相关程度高于峰面积。

(4)不同岩性的含水岩土介质近红外光谱的特征变量与含水量相关性具有相似规律：峰高、右肩宽、峰面积是相关程度最高的三个特性。