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

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

(江蘇師范大學(xué) 地理測繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221116)

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

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

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

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

1.1 研究區(qū)概況

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

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

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

1.2 研究數(shù)據(jù)

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

地表溫度反演精度評定數(shù)據(jù)是來源于地理空間數(shù)據(jù)云MODIS陸地標(biāo)準(zhǔn)產(chǎn)品中的2015年12月17日和2016年4月7日的MOD11A1 1 KM地表溫度/發(fā)射率產(chǎn)品。

1.3 數(shù)據(jù)預(yù)處理

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

2 研究方法

2.1 研究方法及原理

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

(1)

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

2.2 地表參量確定

2.2.1 歸一化植被指數(shù)(NDVI)

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

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

(2)

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

2.2.2 地表溫度(LST)

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

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

(3)

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

黑體輻射亮度公式為：

(4)

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

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

ε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等改進(jìn)及驗證后的單窗算法[11]，其計算公式為：

LST=

(9)

C=ετ

(10)

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

(11)

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

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

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

(12)

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

(13)

T0=t+273.15

(14)

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

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

(15)

(16)

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

3 結(jié)果與分析

3.1 實驗分析

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

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

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

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

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

3.2.1 空間分布特征

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

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

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

3.2.2 時間分布特征

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

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

2014年與2013年相比，潮濕、濕潤、適宜類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從68.272 2 km2變?yōu)?7.918 1 km2，而較旱、干旱類型所占比重則顯著上升，面積分別從13.989 6、3.416 4 km2變?yōu)?6.2032、39.088 8 km2。2016年與2015年相比，潮濕、濕潤、適宜類型所占比重明顯上升，適宜類型比重上升更為顯著，面積從16.213 5 km2變?yōu)?7.483 9 km2，較旱、干旱類型所占比重則明顯下降，干旱類型比重下降更為顯著，面積從54.349 2 km2變?yōu)?.184 3 km2。2018年與2017年相比，潮濕、濕潤、適宜、干旱類型所占比重明顯下降，濕潤類型比重下降更為顯著，面積從12.094 2 km2變?yōu)?.307 7 km2，而較旱類型所占比重則明顯上升，面積從55.027 8 km2變?yōu)?6.551 2 km2。2020年與2019年相比，濕潤、適宜類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從59.208 3 km2變?yōu)?7.027 1 km2，而較旱、干旱類型所占比重則明顯上升，干旱類型比重上升更為顯著，面積從5.507 1 km2變?yōu)?0.470 3 km2。2021年與2020年相比，潮濕、濕潤、適宜、較旱類型所占比重明顯下降，適宜類型比重下降更為顯著，面積從17.027 1 km2變?yōu)?.554 1 km2，而干旱類型所占比重則明顯上升，面積從40.470 3 km2變?yōu)?2.568 8 km2。

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

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

3.2.3 土壤含水量動態(tài)分析

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

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

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

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

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

4 討論

采用植被供水指數(shù)法對土壤含水量進(jìn)行反演，研究結(jié)果與實際調(diào)查結(jié)果基本相同，但由于缺乏實測數(shù)據(jù)，因此VSWI反演的土壤含水量結(jié)果與實際土壤含水量相關(guān)性較弱，且由于影像拍攝時間對應(yīng)的當(dāng)天降雨量均為0 mm，所以并未將降雨量作為影響土壤含水量的因素進(jìn)行分析，后續(xù)研究會將數(shù)據(jù)時空范圍繼續(xù)擴大，并將降雨量、大氣溫度、地形等因素考慮進(jìn)來，便于更精確地監(jiān)測土壤水分的變化。同時，通過長時間序列的反演對比，可清晰地得出該地土壤含水量時空分布特征，對改善當(dāng)?shù)氐纳鷳B(tài)環(huán)境和提高農(nóng)業(yè)產(chǎn)量具有重大的意義。

5 結(jié)論

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

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

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

4)2013—2021年研究區(qū)土壤含水量類型發(fā)生顯著變化，濕潤和適宜類型逐漸變?yōu)檩p旱或干旱類型?？傮w上，9年內(nèi)研究區(qū)土壤含水量變化趨勢較差。

5)煤礦開采過程中排水等一系列操作會使土壤含水量有短暫的升高，但長期開采會導(dǎo)致土壤含水量的降低，所以要及時應(yīng)對煤礦開采對土壤含水量產(chǎn)生的影響，實時監(jiān)測礦區(qū)土壤含水量的變化對改善生態(tài)環(huán)境是必不可少的。