基于Landsat 8遙感影像的峁底礦區土壤含水量反演研究

王丹萍，康建榮，練靖文，崔騰飛

(江蘇師范大學 地理測繪與城鄉規劃學院，江蘇 徐州 221116)

煤炭開采引起了一系列的環境問題，其中包括對土壤含水量的影響。煤炭資源的開采往往會造成地面沉陷，從而導致地下含水層遭到破壞，使得地下水位不斷下降，因此造成土壤水分的減少[1]。由于采煤沉陷造成土壤含水量的降低，若使土壤含水量盡可能恢復，則需要一個較長的周期[2]。對于礦區植被的生長及受損后的恢復，土壤含水量是關鍵的環境因子[3]。同時，土壤含水量對農業生產有著不可替代的作用，其直接影響作物產量。所以土壤含水量的監測及預估對于研究區內生態環境、農業等領域的分析和治理具有很大的作用。

峁底煤礦自開采以來，周邊地形發生了較大的變化，改變了土壤中水分含量及分布狀況，影響到植被的生長，使得許多土地類型演變為裸土。由于研究區地形主要為山地，在降水較多時，因為缺少植被的保護，采空沉陷形成的坡面容易引起水土流失，甚至引發泥石流等自然災害的發生。該礦區處于半干旱地區，地表水資源極其寶貴。土壤水分是制約半干旱地區植被生長和生態環境恢復的關鍵性因子[4]。

針對以上問題，為探索研究區內土壤含水量時空分布變化特征，本文采用植被供水指數法對研究區內的土壤含水量進行反演。以峁底煤礦所在鄉鎮為研究區，以長時間序列的Landsat 8遙感影像為數據源，反演研究區2013—2021年的土壤含水量，并對其整體的時空分布特征進行分析，進而揭示研究區土壤含水量的變化情況，最終結合煤礦開采工作對土壤含水量的影響，為礦區內水資源問題研究提供一定理論和數據幫助。

1 研究區與數據

1.1 研究區概況

奧家灣鄉，隸屬于山西省呂梁市興縣，地處興縣東部，位于晉西北丘陵地區，地勢整體東高西低，境內群山蜿蜒，縱橫交錯，丘陵密布，地形支離破碎。該地氣候為溫帶大陸性氣候，四季分明，年平均氣溫較低，一般夏季多雨，冬季少雪，但近幾年春夏季氣候偏旱；境內河道屬黃河流域，主要河道有蔚汾河、嵐漪河；境內有豐富的礦產資源、牧草資源等。

峁底煤礦，隸屬于奧家灣鄉，位于38°25′54″～38°27′39″N，111°09′47″～111°11′23″E。礦區地處蘆芽山東麓，此處地形較復雜，侵蝕沖刷嚴重，形成了接近南北走向的山脊及溝谷，整體地勢南高北低(見圖1)。井田中部為最高處，海拔為1 234.9 m；井田北部為最低處，海拔為1 028.0 m，最大相對高差為206.9 m，屬于中低山區。礦區內的河流屬于黃河流域蔚汾河水系，礦區內無常年性河流，各溝谷一般干涸無水，到雨季時作為泄洪通道。該地區一月份最低平均氣溫為-7.8～11.2 ℃，極端最低為-29.3 ℃；七月份最高平均氣溫為22.3～25 ℃，極端最高為38.4 ℃。年平均降水量為231.4～688.9 mm，60%以上的降水主要集中在夏季及秋初；年平均降水蒸發量為2 090.8 mm，達到降水量的4倍。礦區及周邊地區所覆蓋的植被主要以灌木叢為主，分布較稀疏，同時有少量農田分布。由于地處黃土高原地區，氣候干旱，經氣候變化及煤礦開采的影響，境內山地地表類型逐漸演變為裸土，植被日益減少。

圖1 峁底煤礦西南部及周邊環境的八月份無人機影像數據Fig.1 UAV image data of the southwest and surrounding environment of Maodi coal mine in August

1.2 研究數據

隨著遙感領域研究的不斷發展，目前已經實現了大范圍區域內土壤含水量的實時監測和動態預測，通過遙感影像監測土壤含水量能夠較好地反映出土壤含水量在一定時間和空間上土壤水分的變化特點[5]。Landsat 8衛星攜帶的陸地成像儀(OLI)和熱紅外傳感器(TIRS)的波段劃分精細，對植被、溫度以及土壤含水量都有極佳的監測效果[6]。Landsat 8影像數據量大且種類豐富，可以完全覆蓋監測區域，因此，本文所采用的影像資料是來自于地理空間數據云(https：//www.gscloud.cn/)2013—2021年的Landsat 8 OLI_TIRS影像數據，云量均在10%以下，軌道號為126/33。

地表溫度反演精度評定數據是來源于地理空間數據云MODIS陸地標準產品中的2015年12月17日和2016年4月7日的MOD11A1 1 KM地表溫度/發射率產品。

1.3 數據預處理

在ENVI5.3中對Landsat 8遙感影像數據進行預處理，包括裁剪、輻射定標、大氣校正等操作，以改善遙感影像的質量，為后續實驗作好準備。同樣在ENVI5.3中對MODIS地表溫度產品數據進行預處理，包括格式轉換、重采樣、重投影、裁剪、溫度換算等，為后續精度評定做準備。

2 研究方法

2.1 研究方法及原理

植被供水指數(VSWI)是一種綜合植被狀態和地表溫度信息，反演有植被覆蓋地區表層土壤含水量的監測方法。在植物供水狀態正常的條件下，一定生長期內的植被指數和植被冠層溫度都能維持在一定范圍內，但在植被處于缺水狀態時，植被的生長情況會遭受影響，植被指數會隨之降低，而植被冠層溫度則會隨之上升[7]。VSWI值可表示植被相對干旱程度，VSWI值越小，植被的干旱程度越嚴重，土壤含水量越低，反之土壤含水量越高。公式為：

(1)

式中，NDVI為歸一化植被指數；LST為地表溫度。

2.2 地表參量確定

2.2.1 歸一化植被指數(NDVI)

歸一化植被指數一般與植被的繁盛程度有關，因此，可以利用歸一化植被指數對植被覆蓋度等進行研究分析。一般情況下，當作物缺水時，其生長情況會受到影響，植被指數也會隨之降低。公式為：

NDVI=(ρNIR-ρRED)/(ρNIR+ρRED)

(2)

式中，ρNIR為近紅外波段反射率；ρRED為紅光波段反射率。本研究使用Landsat 8 OLI傳感器的第四波段紅光波段和第五波段近紅外波段來計算NDVI值。

2.2.2 地表溫度(LST)

地表溫度是一個反映地表能量平衡的客觀指標。對于裸地，地表溫度被稱為土壤溫度；對于植被，地表溫度被稱為冠層溫度[8]。由于土壤溫度的變化能夠間接反映土壤水分的變化，當冠層溫度升高時，植被的缺水程度也逐漸升高[9]。本文中采用了大氣校正法、單窗算法、基于影像的反演算法(Image-based Method，簡稱IB算法)來反演地表溫度，并將三種方法的反演結果與MODIS地表溫度產品進行精度驗證，選出相關性最高的反演結果進行VSWI值的提取。

1)大氣校正法。利用大氣校正法反演地表溫度的計算公式為：

(3)

式中，B(Ts)為黑體輻射亮度；K1、K2為指定波段的熱轉換參數，分別為774.885 3 W·(m2·sr·μm)-1和1 321.078 9 K。

黑體輻射亮度公式為：

(4)

式中，Lλ為影像熱紅外波段進行輻射校正后得到的星上輻射亮度；ε為地表比輻射率；τ、Lup、Ldown分別為大氣透過率、大氣向上輻射亮度、大氣向下輻射亮度，三者可從NASA官網(https：//atmcorr.gsfc.nasa.gov/)獲取。

地表比輻射率的計算可通過基于NDVI閾值的方法[10]進行估算，公式為：

εwater=0.995 (NDVI<0)

(5)

(0≤NDVI≤0.7)

(6)

(NDVI>0.7)

(7)

式中，PV為植被覆蓋度，其計算公式為：

(8)

式中，NDVIsoil和NDVIveg分別為裸地和植被完全覆蓋的像元NDVI值，本研究中選取NDVIsoil=0.05，NDVIveg=0.7。

2)單窗算法。本文中的單窗算法使用的是WANG和QIN等改進及驗證后的單窗算法[11]，其計算公式為：

LST=

(9)

C=ετ

(10)

D=(1-τ)[1+(1-ε)τ]

(11)

式中，ε為地表比輻射率；T6為地面亮溫值；a、b為常數，分別為-67.355 351和0.458 606；Ta為大氣平均作用溫度；τ為大氣透過率，可通過NASA官網獲取。

大氣平均作用溫度公式為：

Ta=16.0110+0.92621T0(中緯度夏季)

(12)

Tb=19.2704+0.91118T0(中緯度冬季)

(13)

T0=t+273.15

(14)

式中，t為大氣平均氣溫，本文通過天氣網(http：//lishi.tianqi.com/)查詢到當天最高和最低氣溫，根據遙感影像拍攝時間進行估計。

3)IB算法。利用IB算法反演地表溫度的計算公式為：

(15)

(16)

式中，T為輻射亮溫，單位為K；λ為熱紅外波段的中心波長值，對于Landsat 8影像，其值為10.9 μm；δ為Boltzmann常數，值為1.38×10-23J/K；h為Plank常數，值為6.626×10-34Js；c代表光速，值為2.998×108ms；ε為地表比輻射率。

3 結果與分析

3.1 實驗分析

為了獲取最佳的地表溫度數據，將三種方法反演出的地表溫度結果分別與MODIS溫度產品進行相關性分析(圖2)。大氣校正法、單窗算法、IB算法與MODIS溫度產品分析結果從整體上來看，都具有良好的線性關系，相關系數分別為0.677 7、0.677 6、0.677 6。

圖2 三種不同算法反演的地表溫度與MODIS溫度產品線性擬合結果Fig.2 Linear fitting results of land surface temperature and MODIS temperature products by inversion of three different algorithms

根據相關性比較，本研究采用了大氣校正法反演地表溫度，由于無法精確地收集到研究區的大氣平均氣溫，所以反演得到的地表溫度精度會略低，但與MODIS溫度產品數據仍具有較強的相關性。

通過采用植被供水指數法反演得出研究區土壤含水量的時空分布情況。從分布情況來看，顏色由淺到深，土壤含水量則由低到高，且土壤含水量程度主要以較旱為主。研究區春夏季氣候較干旱，所以大部分地區土壤含水量較低，得到的結果與實際調查結果基本一致，此方法對于反演研究區土壤含水量具有一定的可行性。

3.2 土壤含水量時空分布特征

3.2.1 空間分布特征

根據研究區土壤含水量分布情況及相關資料，本文將土壤含水量等級分為5級：干旱(<20%)、輕旱(20%～40%)、適宜(40%～60%)、濕潤(60%～80%)、潮濕(>80%)。從2013—2021年研究區土壤含水量的分布情況(圖3)來看，土壤含水量等級主要以輕旱和干旱為主。研究區北部植被覆蓋度高，主要以潮濕和濕潤為主；中部主要以適宜和輕旱為主，境內有省道貫通，所以省道所處位置表現為干旱；南部植被覆蓋度低，有大量裸土分布，主要以輕旱和干旱為主。

圖3 2013—2021年研究區土壤含水量空間分布情況圖Fig.3 Spatial distribution of soil moisture content in the study area from 2013 to 2021

2013年，土壤含水量整體偏高，總體上呈“北部和中部高，南部較低”的特點，北部土壤水分含量最高，以潮濕為主，而省道及附近村莊建筑土壤含水量則最低。2014年，土壤含水量呈現出“北部和東南部高，西部低”的特點，北部土壤水分含量最高，以潮濕為主，省道、附近鄉村建筑及西部地區土壤水分含量較少。2015年、2017年和2020年土壤含水量分布特點均與2014年基本一致，土壤含水量總體上與2014年大致相同。2016年和2019年土壤含水量的分布特點與2013年比較接近，2016年土壤含水量總體上與2013年大致相同，但2019年整體土壤含水量低于二者。2018年，土壤含水量呈“北部和東部高，西部低”的分布特點，東北部土壤水分含量最高，以適宜為主，而省道、附近鄉村建筑及西部地區則相對較少。2021年，土壤含水量總體偏低，呈現出“北高南低”的特點，北部土壤含水量偏高，以輕旱為主，省道、附近村莊建筑及中、南部大部分地區土壤水分含量都比較低。

3.2.2 時間分布特征

2013—2021年土壤含水量情況表現出整體波動降低的態勢?？傮w上，2013—2015年，土壤水分含量呈遞減趨勢；2015—2018年，土壤水分含量也呈現出明顯的降低趨勢；而2018—2021年，土壤含水量出現了先增加后降低的變化。2013—2021年土壤含水量各類型(表1)所占面積及比例變化明顯。

表1 2013—2021年研究區土壤含水量時間變化情況Table 1 Temporal change of soil moisture content in the study area from 2013 to 2021

2014年與2013年相比，潮濕、濕潤、適宜類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從68.272 2 km2變為17.918 1 km2，而較旱、干旱類型所占比重則顯著上升，面積分別從13.989 6、3.416 4 km2變為56.2032、39.088 8 km2。2016年與2015年相比，潮濕、濕潤、適宜類型所占比重明顯上升，適宜類型比重上升更為顯著，面積從16.213 5 km2變為57.483 9 km2，較旱、干旱類型所占比重則明顯下降，干旱類型比重下降更為顯著，面積從54.349 2 km2變為2.184 3 km2。2018年與2017年相比，潮濕、濕潤、適宜、干旱類型所占比重明顯下降，濕潤類型比重下降更為顯著，面積從12.094 2 km2變為1.307 7 km2，而較旱類型所占比重則明顯上升，面積從55.027 8 km2變為86.551 2 km2。2020年與2019年相比，濕潤、適宜類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從59.208 3 km2變為17.027 1 km2，而較旱、干旱類型所占比重則明顯上升，干旱類型比重上升更為顯著，面積從5.507 1 km2變為40.470 3 km2。2021年與2020年相比，潮濕、濕潤、適宜、較旱類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從17.027 1 km2變為3.554 1 km2，而干旱類型所占比重則明顯上升，面積從40.470 3 km2變為72.568 8 km2。

從整體上來看，2014年與2015年、2017年、2020年相比，土壤含水量各類型面積差別較小，土壤含水量較低，集中分布在20%～40%，面積分別占研究區總面積的43.68%、36.96%、42.77%、43.40%；2013年與2016年相比，土壤含水量各類型面積差別也較小，而土壤含水量則相對較高，集中分布在40%～60%，面積分別占研究區總面積的53.06%、44.68%?？傮w上，2013年、2016年及2019年的土壤含水量分布狀況要優于其他年份，從土壤含水量程度占比(圖4)中可明顯看出，這三個年份的土壤濕度相對較高。

圖4 2013—2021年土壤含水量程度占比柱狀圖Fig.4 Bar chart of the percentage of soil moisture content from 2013 to 2021

3.2.3 土壤含水量動態分析

研究區各類型土壤水分的面積變化有明顯的差異(表2)，其中變化最大的是干旱類型，土壤含水量<20%，減少了69.1524 km2，變化率為53.75%；面積變化最小的是潮濕類型，土壤含水量>80%，增加了11.3724 km2，變化率為8.84%。

表2 2013—2021年研究區土壤含水量面積變化差異統計Table 2 Difference statistics of soil moisture content area change in the study area from 2013 to 2021

將土壤含水量面積變化率劃分為7個等級：劇烈減少區(變化率<-30%)，減少區(-30%≤變化率<-10%)，輕度減少區(-10%≤變化率<-5%)，穩定區(-5%≤變化率<5%)，輕度增加區(5%≤變化率<10%)，增加區(10%≤變化率<30%)和劇烈增加區(變化率≥30%)[5]。從土壤含水量類型面積變化情況來看，干旱類型面積變化最大，可歸為劇烈增加區；適宜類型面積變化程度僅次于干旱類型，可歸為劇烈減少區?？傮w來看，2013—2021年研究區土壤含水量呈下降趨勢。

3.3 煤礦開采對土壤含水量的影響

2013—2021年，峁底煤礦均處于開采狀態。采礦區與省際干線公路相鄰，位于公路南側，從圖3中可以看出采礦區的土壤含水量屬于較高的區域，由于礦井會進行一系列的排水工作使礦區的開采能正常進行，所以此處的土壤含水量會較其他地區略高，而變化趨勢與整體變化趨勢基本一致。礦區與其他區域的土壤含水量在研究年份前幾年差異比較顯著，但后面幾年差異逐漸縮小，考慮到煤礦開采會導致土壤受損，且土壤的自主修復周期較長，會產生許多無法閉合的裂縫，而這些裂縫會加速土壤水分的流失和蒸發[4]，因此礦區土壤含水量的變化趨勢是整體降低的。針對礦區土壤含水量減少的問題，應及時對受損土壤進行人工修復，及時填補裂縫，減少土壤水分的流失。

4 討論

采用植被供水指數法對土壤含水量進行反演，研究結果與實際調查結果基本相同，但由于缺乏實測數據，因此VSWI反演的土壤含水量結果與實際土壤含水量相關性較弱，且由于影像拍攝時間對應的當天降雨量均為0 mm，所以并未將降雨量作為影響土壤含水量的因素進行分析，后續研究會將數據時空范圍繼續擴大，并將降雨量、大氣溫度、地形等因素考慮進來，便于更精確地監測土壤水分的變化。同時，通過長時間序列的反演對比，可清晰地得出該地土壤含水量時空分布特征，對改善當地的生態環境和提高農業產量具有重大的意義。

5 結論

1)基于Landsat 8遙感影像反演的土壤含水量情況與實際調查結果基本一致，所以基于Landsat 8數據和植被供水指數法反演峁底礦區及周邊環境的土壤含水量具有一定的可行性。

2)大氣校正法作為研究區反演地表溫度的方法效果更佳，與MODIS溫度產品數據進行相關性分析，通過了雙尾顯著性檢驗，相關性顯著。

3)在2013—2021年共9年的時間跨度下，研究區土壤含水量的整體變化呈波動降低的趨勢，空間整體分布上以“北部高，東部較高，西部低”趨勢為主，所以利用遙感監測研究土壤含水量，能較好地分析土壤含水量的時空變化特征。

4)2013—2021年研究區土壤含水量類型發生顯著變化，濕潤和適宜類型逐漸變為輕旱或干旱類型?？傮w上，9年內研究區土壤含水量變化趨勢較差。

5)煤礦開采過程中排水等一系列操作會使土壤含水量有短暫的升高，但長期開采會導致土壤含水量的降低，所以要及時應對煤礦開采對土壤含水量產生的影響，實時監測礦區土壤含水量的變化對改善生態環境是必不可少的。