煤礦開采中SOM的遙感估算和時空動態分析

高文龍， 張圣微,2,3， 林 汐， 雒 萌， 任照怡

(1.內蒙古農業大學水利與土木建筑工程學院，呼和浩特 010018； 2.內蒙古自治區水資源保護與利用重點實驗室，呼和浩特 010018； 3.內蒙古自治區農牧業大數據研究與應用重點實驗室，呼和浩特 010018)

0 引言

全球和局部碳循環研究一直是近年生態環境領域的研究熱點[1]。土壤是最大的陸地碳儲層，其儲存的碳量是大氣的兩倍以上[2]。土壤有機質(soil organic matter, SOM)包括各種動植物殘體、微生物以及分解和合成的各種有機物質，是土壤中碳的主要表現形式，被稱作植物的“養分銀行”[3-5]。因此，測量SOM是判斷土壤肥力的重要途徑之一[6]。常規實驗測定方法費事且成本較高，無法滿足大范圍應用、快速有效監測的需求，而多波段、信息豐富的遙感技術為解決這一問題提供了一種新途徑[7]。但衛星遙感數據也存在空間分辨率較低、對地觀測受天氣和大氣條件影響較多等問題[8]。以往的研究主要集中在提高農業生產、環境污染監測、土地利用變化等方面的SOM預測[9-11]，而對開采中的煤礦區的研究卻很少見。煤礦開采過程中，土壤理化性質會發生很大程度的改變[12]，因此及時有效地監測礦區土壤質量狀況、快速全面地得到區域性量化指標、準確估算礦區SOM對于其生態環境監測和修復具有十分重要的意義[13]。

近地高光譜遙感相較于衛星遙感來說，能準確地得到實時地表狀態、快速獲取土壤連續光譜數據，可獲得成百甚至上千個窄波段光譜信息，進而得到完整連續的地物光譜曲線[14]。近幾年學者通過高光譜技術對SOM與反射光譜之間的關系、土壤反射光譜特性、光譜數據變換方法、有機質敏感波段和定量反演建模方法等進行了深入研究，取得了很多重要成果[15-19]。同時也發現，土壤有機質的含量變化引起的光譜差異，會受到土壤水分和土壤質地等因素影響[20-21]。在利用高光譜進行SOM反演過程中要綜合考慮土壤水分和質地等影響因素，同時由于近地高光譜監測范圍有限，只有利用其獲取的精確土壤光譜數據與衛星遙感數據相結合，才能更好的對區域SOM進行準確的估算[22-23]。

因此，本文以位于陜西省榆林市和內蒙古鄂爾多斯市交界處的中煤集團榆林礦區為研究區，通過全區域樣點土壤取樣，測得土壤有機質含量； 利用成像高光譜儀獲得樣品光譜數據，采用光譜分析、相關性分析和逐步回歸分析等方法，構建基于高光譜和Landsat8 OLI衛星多光譜數據的SOM反演模型。使用所建立的反演模型，對研究區內的大海則、營盤壕、納林河二號、巴拉素等重點煤礦進行SOM反演，并對礦區開采前后SOM的空間變化狀況進行探討。這些結果將對礦區土壤資源和生態環境定量研究、管理以及可持續發展提供依據。

1 研究區概況與數據源

1.1 研究區概況

研究區位于陜西省榆林市與內蒙古鄂爾多斯市境內，地理位置在E108°38′～109°26′、N37°59′～38°42′之間，北連鄂爾多斯市市區，南接榆林市橫山縣，總面積約為3 800 km2； 地勢自西北向東南降低，海拔范圍在989～1 409 m，橫跨烏審旗、榆陽區和橫山縣； 年降水量200～400 mm，屬干旱半干旱地區。研究區內遍布煤礦，以井工煤礦為主，包括： 中煤集團的大海則煤礦、納林河二號煤礦，兗礦集團營盤壕煤礦，延長石油巴拉素礦等15個采樣點； 流域水系主要為納林河、海流兔河和硬地梁河并全部匯入無定河流域(圖1)。

圖1 研究區概況及采樣點位置Fig.1 Overview of study area and sampling point location

1.2 遙感數據獲取

光譜數據包括衛星遙感影像和近地土壤高光譜數據。2019年分別在15個采樣點位置，在不同月份，對近地土壤高光譜和土壤樣本數據進行3次采樣，采集時間分別為： 6月1日、7月4日和9月21日。衛星遙感影像取自美國地質調查局(https: //earthexplorer.usgs.gov/)空間分辨率30 m的Landsat8 OLI影像的3景少云數據，并且遙感數據是在實地采樣前后5 d內(近同步)獲取的。近地高光譜和土樣采集分布在大海則礦、營盤壕礦、納林河二號、巴拉素礦及流域周邊(海流兔河匯流口、納林河上煤礦首采區等)(圖1)。

近地高光譜數據在10：00—14：00之間，日照充足、無云的條件下完成采樣[24]。土樣光譜反射率測定通過便攜式高光譜成像系統完成，該系統主要由成像高光譜相機Resonon 的Pika L、水平旋轉云臺(Isuzu Optics Corp，China)、筆記本電腦、數據采集軟件(SpectrononPro)、三腳架和電源裝置等組成。波長觀測范圍約為400～1 000 nm。為保證SOM含量與實際情況高度一致，高光譜波段曲線提取直接在野外布點的0.5 m×0. 5 m樣方內完成，光譜測量幾何條件為： 光源方位角70°，光源照射方向與垂直夾角15°，傳感器垂直放置在土壤表層中心上方50 cm處，光源到土壤表層距離100 cm。每次測量前進行一次標準白板標定，每個樣點重復測量5次，取5條高光譜曲線反射率的算術平均值作為該土樣的實際光譜反射率。在實際測定過程中發現高光譜成像儀獲取的波譜范圍為380～1 020 nm； 因此，為了降低噪聲和水汽吸收產生的影響，需要去除信噪比較低的邊緣波段380～399 nm和1 001～1 020 nm[25]。

土壤有機質樣本采集方式為十字法，在光譜采集結束后，取同一位置深度為0～10 cm的土壤樣本5個，并利用全球定位系統記錄采樣點位置。土樣帶回實驗室后均勻混合，3次合計得到樣本45個，經自然風干、研磨過100目孔篩后，采用重鉻酸鉀氧化外加熱法測定SOM[26]，計算公式為：

(1)

式中：C為0.5 mol·L-1(1/6 K2Cr2O7)標液濃度,mol·L-1；V0為空白滴定用去的FeSO4體積,mL；V為樣品滴定用去的FeSO4體積,mL。

1.3 高光譜數據數學變換和擬合評價

1.3.1 數學變換

為了建立SOM含量與高光譜反射率的敏感關系，采用光譜曲線數學變換方程、相關分析等方法對土壤光譜反射率(R)以及其變換形式—反射率的倒數(1/R)、反射率的對數(lgR)、反射率的一階微分(R′)進行分析，數學變換方程為：

(2)

lgR=lgR(λi)

(3)

(4)

式中：λi-1，λi，λi+1為光譜相鄰波長；R(λi)為光譜反射率；R為λi波長的倒數反射率； 1/R為λi波長的倒數反射率； lgR為λi波長的倒數反射率；R′為λi波長的一階微分反射率。

1.3.2 擬合評價

擬合評價包括決定系數(R2)和均方根誤差(RMSE)，計算公式為：

(5)

(6)

此外，本文還通過比較實測值和預測值的偏離程度，得到特征提取波段組中的最優組合，進行區域性SOM反演。

2 結果與分析

2.1 高光譜土壤特征波段提取

SOM是土壤重要的賦色成分，土壤光譜反射率對其含量的高低會產生一定響應。但在獲取土壤樣本和測量樣本理化性質時，檢測各條光譜反射率是否存在異常值是必不可少的。本文通過人工逐條檢驗法，對45個反射率光譜曲線進行異常樣本排除。同時利用Savitzky-Golay濾波對全部土壤樣本在高光譜400～1 000 nm譜段的光譜反射率數據進行平滑去噪處理[27]。

對土壤光譜反射率進行數學變換，既能降低背景噪聲對目標光譜的影響，又能將非線性關系變為線性關系[28]。本文將原始反射率和3種數學變換(倒數、對數和一階微分)結果用于光譜曲線特征波段的選取和統計分析，以確定最佳光譜變換形式，從而得到相應的特征波段，如圖2所示。

(a) 原始光譜反射率 (b) 倒數反射率 (c) 對數反射率 (d) 一階微分反射率

(e) 原始光譜相關性 (f) 倒數相關性 (g) 對數相關性 (h) 一階微分相關性圖2 土壤有機質實測值與光譜相關性關系Fig.2 Correlation between SOM and Spectrum

將45個土壤光譜反射率及反射率變換處理(倒數、對數和一階微分)結果分別與SOM含量進行相關性分析(圖2(b)—(h))，并基于0.05水平上做相關系數(r)的顯著相關檢驗。結果表明，SOM與原始高光譜反射率在400～642 nm波段范圍內相關性較好； 而反射率的倒數(1/R)和反射率的對數(lgR)與SOM含量相關性絕對值呈弱相關(相關系數最優分別為-0.21和0.28)，并未達到顯著水平； 但在反射率一階微分(R′)中： 波長范圍400～404 nm，410～414 nm，420 nm，432～451 nm，457～459 nm，465～469 nm，490～492 nm，541～544 nm，592 nm，602 nm，623～628 nm，634～636 nm，640～644 nm，666～668 nm，674～679 nm，767～769 nm，773～778 nm，806 nm，868 nm和885 nm處均處于0.05顯著相關，相關順序為R′>R>lgR>1/R，并且一階微分與SOM含量之間正相關性最高(相關系數為0.453)，負相關性最高(相關系數為-0.477)。

結合土壤反射率及其與3種類型數學變換結果的相關系數，不難發現使用一階微分變換的反射率反演區域SOM效果最佳。為保證冗余度低，凸顯有效信息，采用了400 nm，412 nm，420 nm，436 nm，457 nm，467 nm，490 nm，546 nm，592 nm，602 nm，625 nm，636 nm，644 nm，666 nm，676 nm，767 nm，773 nm，806 nm，868 nm和885nm單波段處的變換反射率，其中包括9個正相關波段、11個負相關波段。

2.2 SOM空間分布格局分析

將2019年6月1日、7月4日、9月16日獲取的3景Landsat8 OLI多光譜影像的7個可見光波段進行預處理(輻射定標、大氣校正、幾何糾正、研究區提取等)后，再與每一個高光譜特征波段做相關分析。結果表明，相關性大于0.7的波段分別為546 nm(綠波段)、666 nm(紅波段)和868 nm(近紅外波段)，R2分別為0.70，0.79和0.82。將3次取樣的45個樣本按照2∶1的比例，隨機抽取30個樣本用于模型建立(表1)。變異系數(CV)反映特性參數的空間變異程度，揭示區域化變量的離散程度； 變異系數越大，說明數據的差異和離散程度越大。一般認為，CV<0.1為弱變異性，0.1≤CV≤1為中等變異性，CV>1為高等變異性[29]。從變異系數來看，SOM指標空間變異不大，為中等變異。將綠波段、紅波段、近紅外波段及其波段組合與SOM做線性回歸分析，其余15個樣本用于精度驗證(表2)，得到SOM與Landsat8 OLI多波段數據的最優線性回歸方程為：

表1 土壤有機質含量基本統計特征Tab.1 Basic statistic characteristic values of SOM content

Y=0.08X綠+0.19X紅+0.18X近紅-1.9 。

(7)

表2 線性回歸比較Tab.2 Linear regression comparison

根據線性回歸得出的結論，對6月1日、7月4日和9月21日獲取的3景Landsat8圖像反演研究區的SOM。為消除植被對SOM反演的影響，對研究區內NDVI>0.15的區域進行剔除，得到不同月份礦區的SOM(圖3)。

(a) 2019年6月1日數據反演SOM(b) 2019年7月4日數據反演SOM

(c) 2019年9月21日數據反演SOM(d) 主要礦區及流域不同位置6—9月SOM均值圖3 不同時間SOM的空間分布情況Fig.3 Spatial distribution of SOM in different times

從圖3(a)—(c)可以看出，研究區內土壤有機質含量分布不均勻，空間差異較大，水系流域周邊的SOM明顯高于無水系區域； 東部、南部地區普遍高，SOM平均含量達3.98 g/kg； 西部、西北部地區普遍較低，SOM平均含量僅1.57 g/kg。整個研究區SOM含量處于較低水平，7月SOM總體均值低于6月和9月，這是因為研究區位于毛烏素沙地腹地，土壤表層為栗鈣土，加之常年處于煤礦開采作業的影響下，導致地表植被較少，SOM較低。

對2019年6—9月份研究區內的幾個主要煤礦周邊及納林河上、中、下游及流域出口位置的SOM進行統計分析(圖3(d))?？梢钥闯?，營盤壕礦、納林河二號礦相較于大海則礦、巴拉素礦SOM含量高出32%，這是因為大海則和巴拉素礦目前正在建設中，還沒有正式開始開采。此外，對于納林河流域而言，上、中、下游到流域出口SOM逐漸減少，總體下降0.57 g/kg。

3 問題與討論

3.1 水分對有機質反演的影響

利用遙感技術研究表層SOM空間分布具有省時省力且可信度較高的特點[30]，為土壤環境監測與可持續利用開辟了一條新路[31]。然而土壤含水量(soil water content，SWC)被普遍認為是SOM光譜特征提取的最主要干擾因子[32]，在氣溫氣候、地形地貌、土壤類型等影響因素基本相同的條件下，水與有機質對土壤原始光譜的相互影響是客觀存在的，并且對土壤光譜的作用程度從大到小依次為： 水、有機質和二者相互作用[33-34]。當SWC>10%時，與SOM含量呈現顯著的負相關關系[35-36]； 但當SWC<10%時，原始光譜能較好反映有機質的作用。

本研究中同時獲取了SWC與SOM數據， 圖4為6月1日、7月4日和9月21日實測的SOM和SWC的關系以及反演樣點處的SOM含量。

(a) 7月4日SOM和SWC影響關系(b) 6月1日SOM和SWC影響關系(c) 9月21日SOM和SWC影響關系圖4 土壤含水量對土壤有機質的影響Fig.4 Influence of SWC to SOM

從圖4(a)—(c)可以看出，土壤樣品SWC值都處于10%以下，其中84%的土樣SWC值處于5%以下，除部分樣點部分時期SOM與SWC呈相反關系外，其余大部分沒有明顯關聯性，證明此區域SOM并未受到SWC的影響。

3.2 采礦前后有機質含量的變化

圖5為采礦前SOM反演圖。煤礦的開采會對礦區及其周邊生態環境造成重要影響，SOM也會發生變化。利用2003年9月9日研究區煤礦建設開采前的Landsat5多光譜數據，使用本文構建的線型模型反演了采礦前研究區SOM的空間分布狀況(圖5)。通過對比發現，開礦前SOM明顯高于開采后，SOM范圍為0～10 g/kg，說明采礦過程對于SOM影響極大[37-38]。

圖5 采礦前SOM反演Fig.5 SOM inversion before mining

3.3 SOM回歸反演模型的精度評價

一階微分光譜的線性回歸模型是本研究內的最優估算模型，取得良好的預測結果，其實測SOM和預測SOM量化數據的相關系數達到0.90，如圖6所示。采用該方法反演干旱半干旱區的井工礦區域具有很高的可行性，也為更加快速、準確地估測SOM含量提供了一個新思路。

圖6 SOM反演精度評價Fig.6 Accuracy evaluation of SOM inversion

干旱半干旱地區SOM含量被長期忽視，可能有大量的碳循環過程的存在，其碳循環內部以及煤礦開采過程、生態恢復過程中動態含量變換仍無法在我國西北地區完全精準提取SOM含量[39]，主要由于干旱半干旱地區地貌廣闊，本次研究區的范圍對于整個西北干旱半干旱地區仍然面積較小，因此，基于本次SOM空間分布格局研究仍為將來實時監測SOM奠定了理論基礎。

4 結論

1)對土壤高光譜反射率及數學變換結果與實測SOM作相關分析分析，得到相關性為R′>R>lgR>1/R； 并通過0.05顯著相關性檢驗，提取關系顯著的特征波段為： 400 nm，412 nm，420 nm，436 nm，457 nm，467 nm，490 nm，546 nm，592 nm，602 nm，625 nm，636 nm，644 nm，666 nm，676 nm，767 nm，773 nm，806 nm，868 nm和885 nm。

2)對高光譜、多光譜和SOM進行相關性分析，得到了最優線性回歸方程，并反演預測SOM，得到2019年6月1日、7月4日、9月21日的區域SOM量化數值，量化數據的相關性達到0.90。

3)對研究區SOM空間分布格局及流域不同位置處SOM的分析表明： 研究區總體SOM普遍處于較低水平，流域上、中、下游到流域出口的SOM逐漸減少，且采礦前SOM含量明顯高于采礦后，說明采礦過程對于SOM影響極大。