陳 露，常龍飛，沈沐傲，張鳴原，李德建*，李英駿，3*

（1.華北科技學院 應急技術與管理學院，河北 三河 065201；2.中國礦業(yè)大學（北京） 深部巖土力學與地下工程國家重點實驗室，北京 100083；3.中國科學院物理研究所，北京 100190）

水對巖石具有明顯的弱化、滲透、水楔和溶蝕作用［1］，巖石表面的水分子可通過吸引力與具有較弱連結強度的礦物分子相結合。水分子進入巖石顆粒間隙，將減弱巖石中礦物成分的粘聚力，使礦物顆粒間起連接作用的膠體礦物與水相互作用產生潤滑效果，進而導致巖石強度降低。因此，明確巖石中的水分含量對于地下工程的安全開展至關重要。

紅外光譜技術能無損快速地獲取物質的光譜信息［2-5］，是動態(tài)監(jiān)測物質含水量的有力手段。Savin等［6］分析了土壤表面光譜反射率與其水分含量的關系，發(fā)現土壤水分含量變化與光譜反射率的變化在短時間內同步。Huang 等［7］探究了受壓階段砂巖的光譜特征，發(fā)現砂巖的光譜特征變化與石英含量有關。張芳等［8］基于最大信息數研究了不同巖性對含水光譜特征的影響。中紅外光譜分析技術 （MIR）具有更廣泛的探測領域，表現出的光譜曲線能夠凸顯物質本身的屬性，在定量分析方面也具有較大潛力，利用MIR分析巖石含水量與其光譜特征的關系仍有待探究。本文利用中紅外光譜技術開展砂巖不同含水量的測試，通過寬波段范圍的光譜特征觀察砂巖對水的敏感響應波段并分析不同含水量砂巖的光譜特征，建立了基于中紅外光譜特征的砂巖含水量反演模型，以期為巖石工程現場開采和勘探提供參考。

1 實驗部分

1.1 樣品制備

砂巖樣品均采自甘孜藏族自治州瀘定縣，將其加工成高度100 mm，直徑50 mm 的標準單軸巖石樣品，選擇質地較為均勻的6塊樣品（分別標記為A-1～A-6）進行光譜采集，樣品如圖1所示。

圖1 砂巖樣品示意圖Fig.1 The diagram of sandstone samples

圖2 不同含水量砂巖樣品的中紅外光譜采集過程Fig.2 The acquisition process of MIR spectra of sandstone samples with different water contents

為避免巖樣內部的初始水含量不同，將樣品平穩(wěn)置于105 ℃的真空干燥箱內烘干3 d，利用精密電子天平測量其初始質量。將烘干后的砂巖樣品分別平鋪在100 ℃的恒溫水浴鍋中浸泡至樣品中自由水基本消失且質量幾乎不變，其質量分布如表1所示。

表1 砂巖樣品質量分布Table 1 The quality distribution of the sandstone samples

1.2 方法與光譜采集

利用Foil 20-Z傅里葉變換中紅外光譜系統(tǒng)進行光譜采集，采集范圍為2.5～25 μm。在 23 ℃且濕度小于26%的恒溫實驗環(huán)境中，將樣品放置在固定位置進行測定。探頭的分辨率為4 cm，探頭與測試樣品之間的距離固定并與樣品成90°夾角。每塊巖樣均經歷從含水量為0到飽水狀態(tài)，將其劃分為11個含水量階段?？紤]到巖石的不均質性，每個巖樣在相同含水量狀態(tài)下分別從樣品側面和底部中心點采集光譜曲線2條。去除差異較大的光譜曲線數據，共獲得不同含水量下較為穩(wěn)定的光譜曲線230條。

砂巖樣品在從干燥達到飽和狀態(tài)的過程中，分別在樣品底部（118 條）和側面（112 條）采集不同含水量的光譜曲線。巖石是一種非均質多組分的晶體結構，其內部含有大量的孔隙和節(jié)理，表面平整度較差且相對粗糙，導致所得光譜整體反射率偏低，并且砂巖樣品的光譜曲線在部分波段存在多而雜的“毛刺”。為更直觀地獲取光譜信息，利用FTIRDas系統(tǒng)對校正后的光譜數據進行SG（Savitzky-Golay）濾波平滑［9］。經篩選，確定窗口長度為5 的二階SG 濾波平滑效果最佳。以相同含水量砂巖光譜樣本的算術平均值作為樣品光譜，不同含水量砂巖樣品底部測試點處理后的光譜曲線如圖3所示。

圖3 不同含水量砂巖底部測試點的中紅外光譜曲線Fig.3 The MIR spectral curves of test points at sandstone bottom with different water contents

2 結果與討論

2.1 光譜差異分析

由圖3 可知，不同含水量樣品的光譜曲線分布基本一致，隨著含水量的增加，光譜強度明顯增加。圖4 為不同含水量砂巖的二階微分中紅外光譜曲線，從圖中可以看出在7.23 μm 和10.41 μm 附近存在明顯的吸收峰（依次命名為峰R1和峰R2），且這兩個吸收峰的波峰隨含水量的增加而增強，峰頂位置逐漸向右偏移，說明砂巖樣品分子在相應光譜波段范圍內與水分子產生了相同頻率的振動，導致能量躍遷。通常，水在近紅外光譜中的吸收峰主要集中在O—H 伸縮振動的一級倍頻1 440 nm及水分子的組合頻1 940 nm 附近［10］，是其中紅外吸收基頻的倍頻及合頻。相較而言，中紅外波段范圍內的光譜曲線更能說明水分變量的特征，水在整個中紅外光譜波段范圍內的振動較為活躍，完全干燥與初始含水狀態(tài)下的砂巖樣品光譜差異最為明顯，這是因為初始條件下水分覆蓋在樣品表面，樣品表面與光波產生同頻振動的分子數量明顯增加，吸收的入射光增多，光譜曲線明顯上升；隨著水分的增多，水分子開始向巖石內部顆粒間隙擴散并逐漸趨于飽和。此時樣品表面水分子含量變化不明顯，不同含水量之間的光譜差異逐漸縮小，當波長超過11 μm時，不同含水量砂巖樣品的中紅外光譜的吸收能力逐漸下降。

圖4 不同含水量砂巖的中紅外光譜二階微分曲線Fig.4 The second order differential curves of MIR spectral of sandstone with different water contents

2.2 基于不同含水量砂巖光譜的預測模型構建

為消除散射、樣本不均勻等誤差影響，提高被分析信號中重要信息的權重，采用多元散射校正、歸一化、一階微分方法對光譜數據進行預處理。

2.2.1 結合多元散射校正與偏最小二乘法（PLS）的預測模型多元散射校正（MSC）［11］是將所有光譜數據取平均值作為“理想光譜”，并將其余光譜與平均光譜數據進行一元線性回歸從而校正光譜。計算公式為：

式（1）中，n為樣本數量，xi，j為由樣本數和波長點數構建的光譜矩陣，為所有樣本經簡單處理后在相應波點處的平均值。式（2）中，xi表示單一樣本的光譜數據，ki和bi分別代表光譜和平均光譜的相對偏移系數和平移系數。利用Matlab 對簡單處理后的光譜數據進行多元散射校正，能夠有效消除原始光譜數據由于噪聲造成的基線平移問題。

將處理后的光譜曲線以3∶1的比例劃分為訓練集和測試集進行主成分分析［12］，通過比較主成分分析中前15個特征值的貢獻率，發(fā)現前5個主成分的累計貢獻率超過90%。根據這5個主成分得分計算各樣本到訓練集中心的馬氏距離［13］，篩選異常樣品并剔除離群點，在此基礎上結合PLS［14］建立砂巖樣本含水量的預測模型。馬氏距離的計算公式為：

式中，Si為各訓練集樣本主成分的得分向量，Sˉ為總訓練集樣本的平均得分矩陣，M為訓練集主成分得分的協方差矩陣，Di為各訓練集樣本的馬氏距離，m為訓練集的樣本個數。

訓練集異常樣本閾值（Dth）的計算公式為：

式中，Dm為訓練集樣本馬氏距離的平均值，σd為相應的標準差，e為權重系數。

通過計算得到預處理后訓練集光譜樣本的馬氏距離如圖5 所示。若數據的標準差較大，說明數據分布較為分散，可以適當放寬閾值，避免漏檢異常樣本；若數據的標準差較小，則說明數據分布比較集中，可適當縮小閾值，以減少誤檢正常樣本。使用最優(yōu)閾值調整權重系數（e=2.35）篩選離群光譜數據點，共篩選出馬氏距離大于閾值的異常離群點38個（圖5）。

圖5 基于多元散射校正光譜的訓練集樣本的馬氏距離和最佳閾值Fig.5 Mahalanobis distance and optimal threshold of training set based on MSC spectra

對經多元散射校正、主成分分析及馬氏距離異常值篩選的訓練集光譜樣本進行PLS建模，分析不同含水量砂巖樣本特征值與含水量之間的相關關系。圖6 分別為PLS 模型下訓練集和測試集砂巖含水量與真實含水量的對比結果。

圖6 PLS模型預測的含水量與真實含水量對比結果Fig.6 Comparison of water content between predicted by PLS model and measured

將真實含水量與預測含水量的相關系數（R2)、訓練集的均方根誤差（RMSEC）和測試集的均方根誤差（RMSEP）作為模型評價指標。相關系數越高、均方根誤差越低且主成分數越小，PLS模型的預測效果越好。

結果表明，基于多元散射校正、主成分分析與馬氏距離構建的PLS 模型能夠很好地預測砂巖含水量，訓練集（173個樣本）的R2為0.985，RMSEC為0.068；測試集（57個樣本）的R2為0.985，RMSEP為0.074。此時的主成分數為11。

2.2.2 結合歸一化特征向量提取與隨機森林法的預測模型為更直觀地獲取不同含水量砂巖的有效信息，從光譜數據中選取R1和R2兩個峰附近相對變化較為明顯的波段標記為A1（7.023～7.612 μm）和A2（9.895～10.989 μm），以峰高H、半高寬FWHM、峰面積A、肩寬比S、基線的斜率K為特征變量。峰高H定義為：以橫軸為零基線，經濾波和平滑校正后吸收峰所對應的光譜強度；半高寬FWHM 定義為：吸收譜帶為最大高度一半時譜帶的寬度；峰面積A定義為：以橫坐標為基線所得的曲線積分面積；肩寬比S定義為：右肩寬（R）與左肩寬（L）長度之比?；€（連接吸收峰左右兩肩的直線）斜率K按下式計算：

式中，λ1和λ2分別為左肩及右肩對應的波長，R1和R2分別為λ1和λ2對應的光譜吸光度。

上述5個特征變量構成了初始特征集N=｛H，FWHM，A，S，K｝。由于各個特征變量之間的量綱不同，導致提取相關性較強的特征變量變得困難，因此對原始特征變量進行歸一化處理以快速篩選相關性較高的初始特征，避免因量綱不同導致的信息誤差，歸一化公式如下：

利用樣本數據（xij）減去其相同列的最小值（min（xj）），再除以所在列的極差（max（xj）-min（xj）），即得到歸一化后不同含水量砂巖的光譜（xij*）。A1波段和A2波段中紅外光譜的初始特征變量數值分別如表2和表3所示，將全部230條中紅外光譜曲線按照樣品含水量不同進行標記，含水量為0 的測試樣品的光譜數據記作0-1，0-2，0-3，...，含水量為10%的樣品光譜數據記作1-1，1-2，1-3，...，以此類推。計算5 種特征變量與砂巖樣品含水量之間的線性相關系數，如表4 所示。由表4 可知，A2波段中峰高H、半峰寬FWHM、峰面積A和肩寬比S，4個特征變量與含水量之間的相關性明顯優(yōu)于A1波段，由于該波段特征變量峰高H和峰面積A與含水量之間的相關系數均達到0.9 以上，其余特征變量與含水量的相關性較差，故將初始特征集N簡化為F=｛H，A｝，即由峰高H和峰面積A兩個特征變量構成特征集。

表2 A1波段歸一化中紅外光譜的初始特征變量Table 2 Initial characteristic variables of normalized MIR during A1 waveband

表3 A2波段歸一化中紅外光譜的初始特征變量Table 3 Initial characteristic variables of normalized MIR during A2 waveband

表4 兩個敏感波段的特征變量與水分含量相關系數的絕對值Table 4 Two sensitive wavebands and absolute correlation coefficients between different spectral features and moisture contents

隨機森林［15］具有可處理離散型、高維度無規(guī)則數據的優(yōu)勢，將上述兩種特征變量結合隨機森林進行建模。決策樹的個數取100，特征變量為2。圖7 為隨機森林法得到的訓練集與測試集砂巖預測含水量與真實含水量的對比結果?？梢钥闯觯柧毤蜏y試集的預測含水量與真實含水量的R2均為0.995，RMSEC和RMSEP分別為0.021、0.022，說明該模型預測能力較好。

圖7 隨機森林法預測的含水量與真實含水量對比結果Fig.7 Comparison of water content between predicted by random forest method and measured

2.2.3 基于一階微分法的支持向量機預測模型微分處理［16］是光譜基線校正的方法之一，可以去除背景噪聲或漂移對光譜信號產生的影響，其中一階微分法可有效消除曲線背景噪聲，提高光譜分辨率。基本公式如下：

式中，x'n為一階微分處理后的光譜數據，xn為第n個樣品的光譜，xn+g為第n個樣品經過g窗口寬度變換的光譜數據，g為窗口寬度。

對一階微分變換后的光譜數據進行支持向量機（SVM）分類，確保學習過程中得到的極值解為最優(yōu)解。利用高斯核函數（RBF kernel）選擇最優(yōu)相關參數，SVM模型的核參數γ為RBF函數自帶的參數，懲罰系數C為1.0?；谝浑A微分變換采用SVM 算法［17］進行定量分析，結果如圖8 所示。數據集中包括173 個訓練集和57 個測試集樣本，訓練集預測的含水量與真實含水量的R2為0.997，RMSEC 為0.034，測試集的R2為0.951，RMSEP為0.137。訓練集與測試集的均方根誤差相差較大，說明模型的預測能力一般。

圖8 基于SVM的預測含水量與真實含水量對比結果Fig.8 Comparison of water content between predicted by SVM and measured

由上述3 種建模結果可知，不同預處理方法下的建模效果不同?？梢钥闯?，結合SG 平滑的3 種處理方法中，由于特征變量篩選可對原始光譜數據進行降維，得到強相關性的參數（峰高和峰面積），避免了全波段多特征的復雜原始數據難以準確分辨特征變量的問題，得到了最佳的建模效果，即基于歸一化特征變量提取結合隨機森林的建模效果最佳。

3 結 論

本文通過分析砂巖的中紅外光譜隨樣品含水量的變化，選取與含水量高度相關、位于7.23 μm 和10.41 μm附近的兩個特征峰的峰高和峰面積作為特征變量，采用隨機森林法建立了砂巖含水量的預測模型，獲得了比基于MSC 光譜建立的PLS模型及基于一階微分光譜建立的SVM 模型更為精準的預測結果，為基于紅外光譜技術實現砂巖含水量的無損、快速分析提供了解決方案。