張 芳, 戶佐樂, 王東升，劉雨濛, 謝運鑫, 卓慧慧, 何滿潮

1. 中國礦業(yè)大學(北京)深部巖土力學與地下工程國家重點實驗室，北京 100083 2. 中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083

引 言

近紅外光譜分析(near infrared spectroscopy technique, NIRS)是近幾十年來發(fā)展最快，最引人注目的光譜分析技術之一，被廣泛應用于食品、化工、農(nóng)業(yè)等領域[1]。在巖土工程領域應用近紅外光譜分析檢測、監(jiān)測含水巖石水分狀態(tài)，是近些年發(fā)展起來的新思路，其原理是根據(jù)巖土介質(zhì)O—H基團的吸收強度，建立含水量與近紅外光譜特征之間的定量關系，實現(xiàn)預測巖土介質(zhì)中含水量的目的。與傳統(tǒng)分析方法相比，該方法具有無損、實時、定量的優(yōu)勢。

對于含水量與光譜特征之間建模研究，諸多學者進行了大量的研究工作，尤其在土壤方面研究成果頗多：金慧凝等[2]通過提取反射光譜特征指標，定量分析土壤含水量與反射光譜特征之間關系，建立了土壤水分含量光譜預測模型。包青嶺等[3]利用包絡線消除法提取反射光譜水分吸收特征，并與土壤含水量進行相關性分析，通過隨機森林方法對光譜水分吸收特征參數(shù)進行分類，獲取各參數(shù)對土壤含水量的重要性，運用多元逐步回歸方法，建立土壤水分含量反演模型。婁徑等[4]通過對光譜數(shù)據(jù)進行倒數(shù)、對數(shù)、均方根及一階導數(shù)微分等光譜變換，分析光譜特征，并與土壤含水量進行相關分析，利用多元線性回歸分析建立土壤含水量監(jiān)測模型。

可見，光譜特征選擇是構建準確、穩(wěn)健的定量模型的關鍵。國內(nèi)學者對此進行了很多研究，如：孔清清等[5]基于隨機森林，結合博弈論進行了近紅外光譜特征選擇，利用互信息選擇出無冗余的特征子集，此法應用于近紅外光譜分類中有較高的分類識別率。

利用信息度量法進行波譜的特征選擇具有無參、非線性的優(yōu)勢，能有效度量兩隨機變量之間相關性，能較好地解決特征變量選擇問題，該方法在特征選擇算法中得到廣泛應用。

但是，上述研究成果大多集中于研究土壤含水量與光譜特征之間的相關性，針對巖石的研究工作很少，更沒有涉及巖性對含水巖石光譜特征選擇的影響，即不同巖性，能否選擇相同特征集，或是需要具體問題具體分析，針對不同現(xiàn)場選擇對應的特征集合，這個問題對生產(chǎn)實踐中構建模型時選擇特征變量具有非常重要的指導意義。

故此，利用互信息作為相關程度的度量標準，對比分析不同巖性的含水巖石近紅外光譜的特征選擇結果，以期評價巖性對含水巖石光譜的影響。

1 基于互信息的特征選擇

采用最大信息系數(shù)(maximal information coefficient，MIC)[6]進行特征選擇，其計算思想是：設任意的行數(shù)xi列數(shù)yi下的含水量C與特征參數(shù)f的散點圖網(wǎng)格Gxiyi(j)，計算最大互信息和最大信息系數(shù)

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

j=1, 2, …

(1)

(2)

其中，Imax(C,f,xi,yi)表示在固定行列數(shù)xi與yi情況下，不同的網(wǎng)格劃分方式下的互信息最大值(最大互信息)；I((C,f)|Gxiyi(j))表示含水量C與特征參數(shù)f的散點圖在網(wǎng)格Gxiyi(j)下的互信息，xi,yi

最后選擇滿足預先給定的閾值的特征參數(shù)，組成近紅外光譜特征集S，完成特征的選擇。

2 實驗部分

1.1 樣品制備

礫巖、粉砂巖采自敦煌莫高窟北區(qū)的崖壁，樣品信息如表1、表2所示。因為該礫巖呈半膠結狀態(tài)，不易加工成標準試件，實驗時將礫巖加工成尺寸約為80 mm×90 mm×60 mm的不規(guī)則形狀(圖1)，使之能夠滿足實驗儀器的放置要求。

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

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

圖1 巖樣及其測試點位置(紅點)(a)：礫巖；(b)：粉砂巖；(c)：夯土

夯土試樣(表1和表2)采自敦煌莫高窟108洞室，該夯土強度較低，無法直接用鉆機取出土樣，所以先取土塊，然后加工成長方體土樣，之后再采用打磨的方法制成φ50mm的標準樣品(圖1)。

實驗中，試樣不斷吸水，水分自下而上運移，故沿著試樣高度方向，選取3個近紅外光譜測試點，并盡量避開簽字筆標記區(qū)域，三種巖性樣本的測試點位置如圖1中紅點位置。

1.2 原理及設備

將巖樣放置真空干燥箱內(nèi)，設置箱內(nèi)溫度105～110 ℃，干燥24 h，取出干燥后巖樣冷卻12 h稱重，利用中國礦業(yè)大學(北京)深部巖土力學與地下工程國家重點實驗室何滿潮[8]自主研發(fā)的“深部軟巖水理作用智能測試系統(tǒng)”，進行巖樣吸水模擬室內(nèi)實驗。該系統(tǒng)主要由主體實驗箱、稱重系統(tǒng)和數(shù)據(jù)采集系統(tǒng)三部分組成，如圖2所示。

1.3 光譜采集

在試樣吸水過程中，采用瑞士萬通的XDS SmartProbe近紅外光譜分析儀采集不同位置不同時刻的近紅外光譜。測試時將光纖探頭分別接觸試樣的3個測量點(圖1紅點)，自下往上依次測量，測量的頻率隨試樣的吸水速率適時調(diào)整。實驗參數(shù)如表3，實驗裝置如圖3。

圖2 巖石吸水過程中的近紅外光譜采集實驗設備

表3 近紅外光譜分析儀的實驗參數(shù)[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近紅外光譜分析儀配套軟件提供的一階導數(shù)法對采集的光譜進行預處理，消除背景的常數(shù)平移對近紅外光譜的影響，使數(shù)據(jù)具有更好的連續(xù)性，處理前后的光譜如圖4—圖6所示，限于篇幅，僅列出礫巖、粉砂巖、夯土1號點的光譜圖。

圖4 礫巖1號點近紅外光譜(a)：原始光譜；(b)：一階導數(shù)預處理后光譜

1.5 譜段選擇和初始特征變量

分析圖4—圖6可知，在400～2 500 nm波長范圍內(nèi)有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個初始特征參數(shù)，設定初始特征集F為F={f1,f2,f3,f4,f5,f6}={Area，Height，F(xiàn)WHM，LHW，RHW，LHW/RHW}, 各含水試樣近紅外光譜的初始特征數(shù)值如表4所示。

需要特別強調(diào)的是，在整個實驗歷程中，粉砂巖3號點處采集的21條近紅外光譜沒有明顯的吸收峰，分析其原因。3號點位于粉砂巖頂端，當水分沒有達到這個位置并浸潤到它時，該點始終處于干燥狀態(tài)，所以沒有吸收峰。3號點采集的21條光譜都沒有采集到含水情況下的光譜，為失效光譜。故在后續(xù)分析中將該組實驗數(shù)據(jù)去掉。

圖5 粉砂巖1號點近紅外光譜(a)：原始光譜；(b)：一階導數(shù)預處理后光譜

圖6 夯土1號點近紅外光譜(a)：原始光譜；(b)：一階導數(shù)預處理后光譜

1.6 特征變量歸一化

分析表4可知，由于6個初始特征變量的量綱不同，且特征變量之間的變化幅度不同，可能導致在分析計算過程中，一些數(shù)量級較小的特征變量的作用無法體現(xiàn)，因此對上述表中的原始數(shù)據(jù)進行歸一化變換，將所有變量轉(zhuǎn)換成0-1內(nèi)的數(shù)值，消除量綱和變化幅度不同帶來的影響。

歸一化的方法是將原始數(shù)據(jù)矩陣的各元素減去該元素所在列的最小值后再除以該列元素的極差，公式如下

圖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)

續(xù)表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]中采用的閾值及結論，取初始特征變量之間的相關系數(shù)的閾值為0.95，初始特征變量與含水量之間的相關系數(shù)的閾值為0.5。則利用相關系數(shù)評價上述變量間的相關程度，選取的特征變量如表6。

表6各試樣在峰R1，峰R2處滿足相關系數(shù)閾值要求的特征變量

Table6CharacteristicvariablessatisfyingthecorrelationcoefficientthresholdatpeaksR1andpeaksR2foreachsample

巖性特征變量R1R2礫巖f1, f5f1, f4粉砂巖f2, f4, f5f1, f5夯土f2, f5f1, f4, f5

2.1 礫巖的特征選擇分析

將礫巖特征變量與含水量數(shù)據(jù)做成散點圖，如圖8所示。

利用第1節(jié)中的公式，分別計算峰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的散點圖的規(guī)律相一致。究其原因，礫巖形狀不規(guī)則，因而導致體積計算不準確，含水量計算誤差較大，數(shù)據(jù)規(guī)律分散，表現(xiàn)出特征變量之間的相關程度不強，這組數(shù)據(jù)真正反映的近紅外光譜特征與含水量的相關性不具有代表性意義。

圖8 礫巖特征變量與含水量散點圖(歸一化)(a)：峰R1；(b)：峰R2

2.2 粉砂巖的特征選擇分析

將粉砂巖特征變量與含水量數(shù)據(jù)做成散點圖，如圖9所示。

圖9 粉砂巖特征變量與含水量散點圖(歸一化)(a): 峰R1; (b): 峰R2

利用第1節(jié)中的公式，分別計算峰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的散點圖表現(xiàn)出來的規(guī)律相一致。

表8 粉砂巖特征變量與含水量間的MIC值Table 8 MIC values between characteristic variables and water content of siltstone

2.3 夯土的特征選擇分析

將夯土特征變量與含水量數(shù)據(jù)做成散點圖，如圖10所示。

利用第1節(jié)中的公式，分別計算峰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的散點圖表現(xiàn)出來的規(guī)律相一致。

2.4 不同巖性近紅外光譜特征選擇對比分析

以最大相關系數(shù)MIC值作為指標，評價巖石近紅外光譜的特征變量與其含水量之間的相關性，對于礫巖、粉砂巖、夯土各特征變量按照相關性大小排序結果具體見表10。

由表10中可知，對于粉砂巖和夯土，若只取兩個特征變量，峰R2處只需要關注右肩寬和峰面積即可。峰R1處只需要關注峰高和右肩寬即可，二者的特征變量選擇一致。因礫巖形狀不規(guī)則，計算體積時產(chǎn)生誤差較大，不予以考慮。

綜上所述，對于含水巖土介質(zhì)的近紅外光譜的特征選擇，在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)對比礫巖、粉砂巖、夯土近紅外光譜特征變量與其含水量的最大相關系數(shù)MIC值，表明夯土的近紅外光譜與含水量之間的相關性最強。

(3)不同巖性的近紅外光譜各個特征值與含水量的相關程度不同，具體表現(xiàn)為在1 400 nm附近峰高、右肩寬與含水量都具有較高的相關系數(shù)，只是相關性大小會因巖性不同而不同；1 900 nm附近的右肩寬和峰面積與含水量都具有較高相關系數(shù)，且右肩寬的相關程度高于峰面積。

(4)不同巖性的含水巖土介質(zhì)近紅外光譜的特征變量與含水量相關性具有相似規(guī)律：峰高、右肩寬、峰面積是相關程度最高的三個特性。