基于高光譜特征吸收峰的煤巖識別方法

韋 任, 徐良驥， 孟雪瑩， 吳劍飛， 張 坤

1. 安徽理工大學空間信息與測繪工程學院， 安徽 淮南 232001 2. 深部煤礦采動響應與災害防控國家重點實驗室， 安徽 淮南 232001

引 言

用采煤機切割煤層， 當切割到巖層時， 應用煤巖識別方法可及時調整采煤機滾筒高度來避免造成欠/過切割[1]。 采煤機滾筒調高是通過人工判斷采煤機滾筒截割物料成份實現的[2]。 由于人工識別的效率較為低下， 存在較多不確定因素， 且井下環境復雜， 所以逐步實現自動化無人開采是全世界未來煤炭開采的趨勢。 煤巖識別是實現無人開采的關鍵技術。 如何實現高效準確的煤巖識別一直是國內外的研究熱點。 早在20世紀60年代以來， 美國、 英國這些世界主要產煤國就對煤巖識別技術進行了大量研究[3]。 近些年來， 國內學者在國外學者的理論基礎上提出了很多新的觀點。 李亮等利用探地雷達技術進行煤巖界面探測， 并以郭莊煤礦為試驗地進行試驗， 實現了氣煤與砂巖的區分， 且誤差較小[4]。 孫繼平等提出了基于可見光圖像、 紅外圖像識別、 γ射線、 雷達、 紅外、 有功功率、 震動、 聲音等多參數信息融合煤巖界面識別方法[5]。

縱觀現有的主要的煤巖界面識別方法， 主要是依據煤與巖石的一些物理性質進行區分。 而我國的地下礦井種類較多， 情況復雜， 很難將這些方法普遍用于各地礦井。 高光譜遙感具有波段多而窄、 信息豐富的特點。 HUNT等學者驗證了礦物和有機物在400～2 500 nm間由化學成分和分子結構引起的識別性高光譜吸收帶， 這些吸收帶被廣泛應用于識別礦物和巖石[6]。 以此為基礎， 通過對煤巖的構成成分的差異造成的光譜響應的差別， 結合高光譜每個波段信息含量豐富的特點， 筆者提出了基于光譜的煤巖識別手段。

1 實驗部分

1.1 煤巖樣品的采集

淮南礦區煤種較為豐富， 煤質優良。 煤的品種包括1/3焦煤、 氣煤、 焦煤、 瘦煤等； 巖石種類包括砂巖、 頁巖、 灰質砂巖等。 試驗樣品采集于淮南潘二、 潘四、 謝橋等礦區， 煤樣與巖石均來源于礦區井下， 保持了樣品最初的特征。 為保證煤樣與巖樣的隨機性， 采取均勻采樣的方法， 獲取了總共48組樣本數據。 其中煤樣本23組， 巖樣25組。 ASTER光譜庫的附屬文件指出沉積巖的粗糙表面漫反射光譜反射率等效于其500～1 500 μm粒度粉末試樣的漫反射光譜[7]。 而煤巖均屬于沉積巖。 為了增加模型的泛化能力， 使模型適用于更多場景， 將新鮮的煤樣與巖石碾成塊狀， 而后使用打磨機將其研磨為粉末狀， 采用1 mm篩子對粉末狀的煤與巖石進行篩選， 制成實驗用樣品。

圖1 典型的煤巖樣本

1.2 煤巖光譜曲線的測定與預處理

本實驗所用的是美國ASD公司生產的ASD FieldSpec4高光譜儀， 光譜范圍為350～2 500 nm。 測量時間選在晚上， 保證了暗室的環境。 背景為黑色絨布， 保證視場內僅有白板， 背景反射率為0。 測試光源為1 000 W鹵素燈， 光源照射方向與垂直方向夾角為45°； 光源距離設置為30 cm， 探頭距離5 cm， 探頭與垂直方向保持15°夾角。 測試前進行暗電流校正和白板校正， 測試過程中每隔三十分鐘進行一次白板校正。 每個樣品采集30條光譜曲線， 共取得1 440組數據。

將采集得到的數據進行剔除異常值處理， 而后取算術平均值作為樣品的實際反射率。 因為高光譜數據冗余量大， 煤的整體反射率偏低， 所以煤的光譜曲線部分波段會出現很多“毛刺”， 因此對取得平均值后的數據使用origin軟件進行平滑處理， 達到去除噪聲的效果， 并且能最大程度的保存原始光譜曲線的信息。 部分典型的處理后的煤巖反射圖如圖2所示。

圖2 典型的煤巖光譜曲線圖

1.3 煤巖光譜差異機理分析

由圖2可知， 煤的整體反射率較低， 而巖石反射率偏高。 在1 450 nm附近處巖石有個較強的吸收谷， 部分煤樣如1/3焦煤、 氣煤存在一個較弱的吸收谷。 這是因為1 450 nm附近的波段主要是水分子的O—H官能基伸縮振動的第一倍頻， 由水分子引起的一個吸收帶[8]。 在1 900 nm附近， 由于巖石中含有的二價鐵離子和煤炭中含有的Al2O3， 巖石與少部分煤樣產生了第二個較強的吸收谷[9]。 在2 130～2 250 nm波段范圍內由于Al元素存在形式的不同， 導致了煤與巖石的光譜產生了較大的差別。 煤的組分中Al元素的存在形式主要是Al2O3， 而巖石中Al元素的存在形式則為Al(OH)3。 Al(OH)3的Al—OH晶格振動使得其在2 209 nm附近具有強吸收峰， 而Al2O3在2 209 nm附近不具有強吸收峰[10]。 根據物質的組分特性的不同進行特征波段選取， 相當于對冗余量較大的高光譜信息進行了“降維”。 提取2 130～2 250 nm的波段作為特征向量進行煤巖識別， 使用隨機森林算法進行分類， 并計算分類精度與Kappa系數， 與其他方法進行比較。

1.4 基于煤巖差異機理的煤巖識別算法構建

基于煤巖在2 130～2 250 nm間的差異性， 采取了連續統去除法、 一階微分法、 二階微分法和SCA-SID四種算法來進行煤巖的識別， 具體的識別流程圖如圖3所示。

圖3 煤巖識別流程圖

1.4.1 連續統去除法

連續統去除法是一種有效增強感興趣吸收特征的光譜分析方法， 它可以有效突出光譜曲線的吸收和反射特征[11]。 連續統去除法選取光譜曲線上的峰值， 用一條相對平滑的曲線將各個峰值連接起來， 這條曲線稱之為包絡線[12]。 連續統去除法的作用是將光譜數值歸一化值0—1之間， 擴大了特征波段的區別， 便于提取差異化特征。 經過連續統去除法后的光譜曲線如圖4所示。

經過包絡線去除后， 利用python編程， 提取幾個典型的吸收特征參數：光譜吸收谷處的位置(P)、 吸收深度(H)、 光譜吸收指數(SAI)、 吸收谷面積之比(S)以及基線的斜率K為特征向量區分煤與巖石。 主要的光譜吸收特征參數如下：

吸收位置P： 連續統去除后的吸收谷中， 反射率值最小的點對應的波長值λmin。

吸收深度H： 吸收位置P點到1的距離， 計算公式為

H=1-Rλmm

(1)

式(1)中，Rλmin表示P點的反射率。

光譜吸收指數SAI計算公式為

(2)

式(2)中，d為對稱度， 計算公式為

(3)

式(3)中，λl為吸收谷左肩對應的波長，λr為吸收谷右肩對應的波長。Rl和Rr分別為λl和λr對應的反射率的值。

圖4 連續統去除后的典型光譜曲線值

吸收谷面積之比S的計算公式為

(4)

式(4)中，Al為吸收谷P點左側的面積，A為吸收谷的整體面積。

基線的斜率K： 基線為連接吸收谷左肩和右肩的直線，K的計算公式為

(5)

提取的五個特征組成了一個五維向量集合{P,H, SAI,S,K}， 將這個五維特征向量輸入到隨機森林算法分類器中， 分類器最大深度為5， 最終得到分類結果混淆矩陣如圖5和表1所示， 0為煤樣本， 1為巖石樣本(后面的混淆矩陣中0均為煤樣， 1均為巖樣)， 有一個煤樣和一個巖樣被誤分， 最終正確率達到83.3%。

圖5 連續統去除法混淆矩陣

表1 煤巖識別算法的結果比較

1.4.2 一階微分變換&二階微分變換

微分是一種被廣泛應用的光譜分析方法， 使用微分變換方法可以去除部分線性或接近線性的噪聲對目標光譜的影響， 從而使得光譜中含有的豐富信息更加真實可靠。 在微分變換中， 一階微分與二階微分是較為常用的方法。 煤與巖石的一階和二階微分如圖6(a,b)所示。

圖6 光譜微分變換

光譜一階微分變換的計算公式為

(6)

光譜二階微分變換的計算公式為

(7)

式(6)和式(7)中，i表示波段號，λi表示第i波段的波長值，Rλi表示第i波段波長λi吩應的反射率的值。 表示λi處光譜的一階微分值， 表示λi處光譜的二階微分值。

圖6(a)是煤巖的一階微分曲線圖， 可以看出圖中存在9個差異較明顯的波段， 其中差異較大的有4個波段， 分別為1 416～1 440， 1 860～1 914， 2 122～2 176和2 209～2 240 nm。 而根據煤巖的物性組分中Al的存在形式的差別， 選取2 122～2 176和2 209～2 240 nm作為煤巖識別的波段， 使得識別精度達到最高。 將得到的特征波段分別采取平均值操作， 構成特征向量， 結果可視化如圖7所示， 煤與大部分巖石的一階微分值差異較大， 比較容易識別。 針對這兩個特征波段， 使用高斯徑向基作為核函數， 采用SVM算法進行分類， 對煤與巖石這兩個類別進行識別， 識別結果如表1和圖8(a)所示。

圖7 微分變換特征分布

圖6(b)是煤巖的二階微分曲線圖， 在2 203～2 217 nm處煤巖二階微分值差異明顯， 選取此波段為二階微分識別波段， 做平均值處理， 結果可視化如圖7(b)所示。 采用以高斯徑向基為核函數的SVM算法進行分類識別， 最終得到的結果如表1和圖8(b)所示。

圖8 微分變換法混淆矩陣

數據集中包括36個訓練集樣本和12個測試集樣本， 其中含有5個煤樣和7個巖石樣本。 其中， 一階微分算法中有兩個巖石樣本被誤分， 煤樣均被識別， 識別率達到了83.3%。 二階微分算法中僅有一個巖石樣本被誤分， 識別率達到了91.7%。

1.4.3 基于光譜相關角與光譜信息散度模型融合的方法

光譜角匹配算法是一種通過計算光譜間的矢量夾角來比較物質相似性的方法， 常被用來進行礦物之間的識別[13]。 但由于夾角的大小只與光譜矢量方向有關， 與其輻射亮度無關， 因此當兩種礦物的光譜矢量方向相似而輻射亮度大小有差別時， 區分效果較差， 并且光譜角只考慮夾角的絕對值， 而不能識別待匹配光譜的正負相關性[14]。 光譜相關角則可以避免以上問題， 可以反映出光譜相對于其均值的變化。 光譜信息散度則是以信息論為出發點來比較兩光譜的相似性。 它們的計算方法如下：

光譜相關角(SCA)的計算公式為

(8)

式(8)中，μx和μy為兩光譜曲線的均值。

對于光譜曲線上的x和y， SCA的計算公式為

(9)

(10)

光譜信息散度(SID)的計算公式為

ΨSID(x，y)=D(x‖y)+D(y‖x)

(11)

根據信息論， 光譜x和y的自信息為

I(xi)=-logp(xi)

(12)

I(yi)=-logq(yi)

(13)

其中， 式中的p(x)和q(y)分別為

(14)

(15)

信息論中兩光譜曲線的相關熵為

(16)

(17)

光譜信息散度的計算公式則為

(18)

使用單一的SCA模型或者SID模型只能從單一的角度去反映光譜間的相似性。 采用模型融合的方法， 將光譜的形態信息與信息論中的“信息”相融合， 可以顯著提高礦物的識別能力。

SCA模型與SID模型融合后的SCASID(x,y)tan模型計算公式為

SCASID(x，y)tan=ΨSID(x，y)tan(SCA(x，y))

(19)

計算得出的結果如圖9所示。

圖9 SCA-SID特征分布

采用以高斯徑向基為核函數SVM算法對得到的特征數據進行分類， 最終得到的混淆矩陣如圖10所示， 僅有一個巖石樣本被誤分， 正確率達到91.7%。

圖10 SCA-SID法混淆矩陣

2 結果與討論

(1)從表1和混淆矩陣可以看出分類效果最佳的是二階微分算法和SCA-SID算法。 連續統去除法各有一個巖石樣本和煤樣被誤分， 識別率達到了83.3%， Kappa系數為0.66； 一階微分法有兩個巖石樣本被誤分成了煤樣， 識別率為83.3%， Kappa系數為0.68； 二階微分法和SCA-SID模

型法均僅有一個巖石樣本被誤分成了煤樣， 分類的精度均達到了91.7%， Kappa系數均為0.83。

(2)從分類的精度的角度來說一階微分法和連續統去除法分類精度較低， 二階微分法和SCA-SID模型法的精度與Kappa更高。 從模型的復雜度的角度出發， 二階微分法得到特征向量的計算過程和所花費的時間更少， 時間和空間復雜度更低， 所以在實際工程應用中， 二階微分法比SCA-SID模型法更具有高效性和可靠性。

3 結 論

(1)提出了四種基于吸收特征波段的煤巖識別方法， 結合煤巖構成成分Al元素的存在形式的不同產生的特征吸收峰(2 130～2 250 nm)進行了特征向量提取， 并根據實際情況選用了隨機森林算法和SVM算法作為識別的算法。 結果表明， 連續統去除法、 一階微分法、 二階微分法、 SCA-SID模型法的識別準確率達到了83.3%， 83.3%， 91.7%， 91.7%， Kappa系數分別為0.66， 0.68， 0.83， 0.83。 其中二階微分法的模型復雜度較SCA-SID的更低， 更有利于煤巖的識別。

(2)相比于傳統的煤巖識別方法， 使用高光譜技術是建立在對煤巖光譜差異的機理上進行的。 目前， 基于高光譜技術的煤巖識別算法設計相關研究較少， 本工作提出的四種煤巖識別算法較傳統方法更具有可靠性和高效性， 為實際工況下的煤巖識別提供了方法參考。