煤礦開采中SOM的遙感估算和時(shí)空動(dòng)態(tài)分析

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

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，呼和浩特 010018； 2.內(nèi)蒙古自治區(qū)水資源保護(hù)與利用重點(diǎn)實(shí)驗(yàn)室，呼和浩特 010018； 3.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)大數(shù)據(jù)研究與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，呼和浩特 010018)

0 引言

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

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

因此，本文以位于陜西省榆林市和內(nèi)蒙古鄂爾多斯市交界處的中煤集團(tuán)榆林礦區(qū)為研究區(qū)，通過全區(qū)域樣點(diǎn)土壤取樣，測(cè)得土壤有機(jī)質(zhì)含量； 利用成像高光譜儀獲得樣品光譜數(shù)據(jù)，采用光譜分析、相關(guān)性分析和逐步回歸分析等方法，構(gòu)建基于高光譜和Landsat8 OLI衛(wèi)星多光譜數(shù)據(jù)的SOM反演模型。使用所建立的反演模型，對(duì)研究區(qū)內(nèi)的大海則、營(yíng)盤壕、納林河二號(hào)、巴拉素等重點(diǎn)煤礦進(jìn)行SOM反演，并對(duì)礦區(qū)開采前后SOM的空間變化狀況進(jìn)行探討。這些結(jié)果將對(duì)礦區(qū)土壤資源和生態(tài)環(huán)境定量研究、管理以及可持續(xù)發(fā)展提供依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

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

圖1 研究區(qū)概況及采樣點(diǎn)位置Fig.1 Overview of study area and sampling point location

1.2 遙感數(shù)據(jù)獲取

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

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

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

(1)

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

1.3 高光譜數(shù)據(jù)數(shù)學(xué)變換和擬合評(píng)價(jià)

1.3.1 數(shù)學(xué)變換

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

(2)

lgR=lgR(λi)

(3)

(4)

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

1.3.2 擬合評(píng)價(jià)

擬合評(píng)價(jià)包括決定系數(shù)(R2)和均方根誤差(RMSE)，計(jì)算公式為：

(5)

(6)

此外，本文還通過比較實(shí)測(cè)值和預(yù)測(cè)值的偏離程度，得到特征提取波段組中的最優(yōu)組合，進(jìn)行區(qū)域性SOM反演。

2 結(jié)果與分析

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

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

對(duì)土壤光譜反射率進(jìn)行數(shù)學(xué)變換，既能降低背景噪聲對(duì)目標(biāo)光譜的影響，又能將非線性關(guān)系變?yōu)榫€性關(guān)系[28]。本文將原始反射率和3種數(shù)學(xué)變換(倒數(shù)、對(duì)數(shù)和一階微分)結(jié)果用于光譜曲線特征波段的選取和統(tǒng)計(jì)分析，以確定最佳光譜變換形式，從而得到相應(yīng)的特征波段，如圖2所示。

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

(e) 原始光譜相關(guān)性 (f) 倒數(shù)相關(guān)性 (g) 對(duì)數(shù)相關(guān)性 (h) 一階微分相關(guān)性圖2 土壤有機(jī)質(zhì)實(shí)測(cè)值與光譜相關(guān)性關(guān)系Fig.2 Correlation between SOM and Spectrum

將45個(gè)土壤光譜反射率及反射率變換處理(倒數(shù)、對(duì)數(shù)和一階微分)結(jié)果分別與SOM含量進(jìn)行相關(guān)性分析(圖2(b)—(h))，并基于0.05水平上做相關(guān)系數(shù)(r)的顯著相關(guān)檢驗(yàn)。結(jié)果表明，SOM與原始高光譜反射率在400～642 nm波段范圍內(nèi)相關(guān)性較好； 而反射率的倒數(shù)(1/R)和反射率的對(duì)數(shù)(lgR)與SOM含量相關(guān)性絕對(duì)值呈弱相關(guān)(相關(guān)系數(shù)最優(yōu)分別為-0.21和0.28)，并未達(dá)到顯著水平； 但在反射率一階微分(R′)中： 波長(zhǎng)范圍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顯著相關(guān)，相關(guān)順序?yàn)镽′>R>lgR>1/R，并且一階微分與SOM含量之間正相關(guān)性最高(相關(guān)系數(shù)為0.453)，負(fù)相關(guān)性最高(相關(guān)系數(shù)為-0.477)。

結(jié)合土壤反射率及其與3種類型數(shù)學(xué)變換結(jié)果的相關(guān)系數(shù)，不難發(fā)現(xiàn)使用一階微分變換的反射率反演區(qū)域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個(gè)正相關(guān)波段、11個(gè)負(fù)相關(guān)波段。

2.2 SOM空間分布格局分析

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

表1 土壤有機(jī)質(zhì)含量基本統(tǒng)計(jì)特征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

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

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

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

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

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

3 問題與討論

3.1 水分對(duì)有機(jī)質(zhì)反演的影響

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

本研究中同時(shí)獲取了SWC與SOM數(shù)據(jù)， 圖4為6月1日、7月4日和9月21日實(shí)測(cè)的SOM和SWC的關(guān)系以及反演樣點(diǎn)處的SOM含量。

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

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

3.2 采礦前后有機(jī)質(zhì)含量的變化

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

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

3.3 SOM回歸反演模型的精度評(píng)價(jià)

一階微分光譜的線性回歸模型是本研究?jī)?nèi)的最優(yōu)估算模型，取得良好的預(yù)測(cè)結(jié)果，其實(shí)測(cè)SOM和預(yù)測(cè)SOM量化數(shù)據(jù)的相關(guān)系數(shù)達(dá)到0.90，如圖6所示。采用該方法反演干旱半干旱區(qū)的井工礦區(qū)域具有很高的可行性，也為更加快速、準(zhǔn)確地估測(cè)SOM含量提供了一個(gè)新思路。

圖6 SOM反演精度評(píng)價(jià)Fig.6 Accuracy evaluation of SOM inversion

干旱半干旱地區(qū)SOM含量被長(zhǎng)期忽視，可能有大量的碳循環(huán)過程的存在，其碳循環(huán)內(nèi)部以及煤礦開采過程、生態(tài)恢復(fù)過程中動(dòng)態(tài)含量變換仍無(wú)法在我國(guó)西北地區(qū)完全精準(zhǔn)提取SOM含量[39]，主要由于干旱半干旱地區(qū)地貌廣闊，本次研究區(qū)的范圍對(duì)于整個(gè)西北干旱半干旱地區(qū)仍然面積較小，因此，基于本次SOM空間分布格局研究仍為將來實(shí)時(shí)監(jiān)測(cè)SOM奠定了理論基礎(chǔ)。

4 結(jié)論

1)對(duì)土壤高光譜反射率及數(shù)學(xué)變換結(jié)果與實(shí)測(cè)SOM作相關(guān)分析分析，得到相關(guān)性為R′>R>lgR>1/R； 并通過0.05顯著相關(guān)性檢驗(yàn)，提取關(guān)系顯著的特征波段為： 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)對(duì)高光譜、多光譜和SOM進(jìn)行相關(guān)性分析，得到了最優(yōu)線性回歸方程，并反演預(yù)測(cè)SOM，得到2019年6月1日、7月4日、9月21日的區(qū)域SOM量化數(shù)值，量化數(shù)據(jù)的相關(guān)性達(dá)到0.90。

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