基于Landsat影像準(zhǔn)東礦區(qū)土壤濕度分布特征研究

趙曉帆,曾 強(qiáng)

(1.新疆大學(xué) 生態(tài)與環(huán)境學(xué)院,新疆 烏魯木齊 830017; 2.新疆大學(xué) 干旱生態(tài)環(huán)境研究所,新疆 烏魯木齊 830017;3.教育部綠洲生態(tài)重點(diǎn)實(shí)驗(yàn)室,新疆 烏魯木齊 830017)

由煤炭開采活動引發(fā)的生態(tài)環(huán)境問題,會對區(qū)域土壤資源造成一定破壞[1]。土壤水分與植被生長、分布、演替過程息息相關(guān),其含量影響區(qū)域性氣候變化。礦區(qū)表層土壤水分的變化還可表征煤炭開采擾動情況。準(zhǔn)東礦區(qū)屬脆弱干旱生態(tài)區(qū),地表植被覆蓋率低,礦區(qū)開發(fā)活動對土壤濕度擾動研究可為評價(jià)其生態(tài)環(huán)境影響提供一定基礎(chǔ)。

土壤含水量的獲取有田間實(shí)測法、土壤水分模型法和遙感方法。相較于傳統(tǒng)方法,遙感方法能對大范圍土壤水分進(jìn)行動態(tài)監(jiān)測,目前應(yīng)用最為廣泛。該方法包含可見光—近紅外遙感,熱紅外遙感與微波遙感[2-3]。礦區(qū)開采活動引發(fā)的地表沉陷積水通過遙感影像進(jìn)行監(jiān)測[4]。ZHAO Tianjie等[5]對微波遙感反演土壤濕度算法進(jìn)行了相關(guān)研究,給出當(dāng)前理解各反演土壤濕度算法的定量依據(jù);LI Zhaoliang等[6]分析多種獲取土壤濕度的理論與算法,指出衛(wèi)星遙感觀測為區(qū)域及全球尺度獲得土壤濕度SM時空數(shù)據(jù)提供了獨(dú)特的手段;吳善玉等[7]基于Sentinel-1、 Sentinel-2 遙感數(shù)據(jù)及地面實(shí)測數(shù)據(jù),構(gòu)建BP神經(jīng)網(wǎng)絡(luò)對西班牙薩拉曼卡地區(qū)土壤濕度進(jìn)行了反演;YI Zhiyu等[8]基于遙感影像和土壤濕度對山西省煤炭開采活動引發(fā)的沉降進(jìn)行監(jiān)測。采用可見光—近紅外技術(shù)表征土壤水分是一種簡便有效的方法:李姝蕊[9]、楊麗萍[10]等基于Landsat8影像數(shù)據(jù)分別反演了黑河流域和居延澤地區(qū)的土壤水分變化;馬保東等[11]利用MODIS數(shù)據(jù)和熱慣量法對神東礦區(qū)地表土壤水分進(jìn)行了反演;SANDHOLT I等[12]基于Ts-NDVI特征空間,提出用于評價(jià)土壤濕度的簡化溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index,TVDI),這是一種基于光學(xué)與熱紅外遙感數(shù)據(jù)進(jìn)行植被覆蓋區(qū)域表層土壤水分反演的方法;趙杰鵬等[13]改進(jìn)了SANDHOLT等提出的TVDI模型,能更好地進(jìn)行大范圍土壤水分的遙感估算;楊玲[14]、王麗霞[15]、徐霞[16]、黃海[17]等基于MODIS或Landsat地表溫度和植被指數(shù)數(shù)據(jù),構(gòu)建Ts-NDVI特征空間,獲取溫度植被干旱指數(shù)TVDI,分別揭示了西遼河流域、香格里拉、布爾臺礦區(qū)和黃土礦區(qū)的土壤濕度時空分布格局和主要影響要素;張喆等[18]研究指出NDVI-LST特征空間更適用于植被生長初期、中期的中低植被覆蓋時反演土壤濕度。

基于Landsat系列遙感衛(wèi)星影像數(shù)據(jù),構(gòu)建溫度植被干旱指數(shù)(TVDI)模型反演準(zhǔn)東礦區(qū)2001—2021年的土壤濕度,并對土壤濕度狀況進(jìn)行分級,分析了準(zhǔn)東礦區(qū)土壤濕度時空分布格局及變化狀況。同時,采用GM(1,1)灰色模型預(yù)測礦區(qū)未來年份土壤濕度變化情況。研究結(jié)果可為評價(jià)該礦區(qū)礦業(yè)開發(fā)活動對生態(tài)環(huán)境的影響提供一定基礎(chǔ)。

1 研究區(qū)概況

準(zhǔn)東礦區(qū)位于準(zhǔn)噶爾盆地東部阜康市至木壘哈薩克自治縣的狹長地帶,北緯44°00′～ 45°15′,東經(jīng)88°42′～90°36′,面積15 334 km2。準(zhǔn)東礦區(qū)包括五彩灣、將軍廟、大井、西黑山和老君廟5個礦區(qū),目前開采方式多為露天開采。準(zhǔn)東礦區(qū)地理位置如圖1所示。

圖1 研究區(qū)地理位置

準(zhǔn)東礦區(qū)生態(tài)環(huán)境脆弱,屬溫帶大陸性氣候,降水量少,蒸發(fā)量大,氣候干燥。地表植被較為稀疏,荒漠廣布,生物種群少,其自我修復(fù)能力較弱,具有明顯的干旱區(qū)特征。研究區(qū)土壤類型主要為灰棕模土、石膏灰棕模土、荒漠風(fēng)沙土和荒漠堿土,表層土壤含水率在10%以下[19]。

2 數(shù)據(jù)來源及處理

2.1 數(shù)據(jù)來源

2.1.1 遙感數(shù)據(jù)

考慮植被生長情況和云量等氣象因素的影響,研究數(shù)據(jù)需覆蓋整個研究區(qū)。從美國地質(zhì)勘探局(USGS)網(wǎng)站(https://glovis.usgs.gov/)和地理數(shù)據(jù)空間云網(wǎng)站(http://www.gscloud.cn/)下載2001、2003、2005、2007、2009、2011年Landsat5 TM和2013、2015、2017、2019、2021年Landsat8 OIL_TIRS,每年7—9月的遙感影像,共11年22幅作為計(jì)算準(zhǔn)東礦區(qū)TVDI的遙感數(shù)據(jù)。

2.1.2 氣象數(shù)據(jù)

從氣象大數(shù)據(jù)系統(tǒng)獲取距離礦區(qū)最近氣象站點(diǎn)(奇臺縣站)2001—2021年逐年氣象數(shù)據(jù),包括降水量、氣溫等。

2.2 數(shù)據(jù)處理

運(yùn)用 ENVI軟件對相應(yīng)年份的 Landsat 系列影像,依次進(jìn)行輻射定標(biāo)、大氣校正、影像拼接和裁剪,完成圖像預(yù)處理工作。

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

歸一化植被指數(shù)INDV(Normalized Difference Vegetation Index,NDVI),是由近紅外波段和可見光紅光波段計(jì)算得到表征植被的一種指數(shù)形式。其計(jì)算方式如下:

(1)

式中:RNI為近紅外波段光譜值;R為可見光紅光波段光譜值。

2.2.2 地表溫度LST反演

利用Landsat數(shù)據(jù)反演地表溫度(Land Surface Temperature,LST)的方法有多種,包括基于輻射傳輸方程的算法、單窗算法、單通道算法和分裂窗算法等。考慮到TIRS(Thermal Infrared Sensor,紅外熱傳感器)波段的定標(biāo)偏差[20]和反演精度、普適性問題,采用覃志豪等[21]改進(jìn)的單窗算法,反演地表溫度TLS:

(2)

式中:a和b為普朗克方程相關(guān)經(jīng)驗(yàn)系數(shù);Ta為大氣平均作用溫度,K;T6為衛(wèi)星上遙感器所觀測到的亮度溫度,K;C和D為中間變量。

其中:

(3)

C=ετ

(4)

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

(5)

式中:K1、K2為衛(wèi)星發(fā)射前預(yù)設(shè)輻射常量;E為中間變量;ε為地表比輻射率;τ為當(dāng)天的大氣透射率。

2.3 構(gòu)造TVDI 模型反演礦區(qū)土壤濕度

準(zhǔn)東礦區(qū)植被覆蓋稀疏,屬于低植被覆蓋區(qū),以NDVI-LST特征空間為基礎(chǔ),構(gòu)造溫度植被干旱指數(shù)ITVD(Temperature Vegetation Dryness Index,TVDI)模型反演準(zhǔn)東礦區(qū)土壤濕度的時空變化情況:

(6)

TLSmax=a1+b1INDV

(7)

TLSmin=a2+b2INDV

(8)

式中:TLSmax為給定NDVI值條件下的地表溫度最高值,對應(yīng)NDVI-LST特征空間的“干邊”(ITVD=1);TLSmin為地表溫度最低值,對應(yīng)NDVI-LST特征空間的“濕邊”(ITVD=0);a和b為系數(shù),可由NDVI-LST特征區(qū)間“干邊”“濕邊”擬合得到。

TLS越接近“濕邊”,ITVD越小,土壤濕度越高;TLS越接近“干邊”,ITVD越大,土壤濕度越低。

利用ITVD值作為不同土壤濕度分級指標(biāo),結(jié)合準(zhǔn)東礦區(qū)實(shí)際情況,將土壤濕度劃分為5個等級,分別為:一級(0≤ITVD<0.76)、二級(0.76≤ITVD<0.82)、三級(0.82≤ITVD<0.88)、四級(0.88≤ITVD<0.94)和五級(0.94≤ITVD<1.00)。等級越低,土壤濕度越高;等級越高,土壤濕度越低,即五級土壤濕度最低。

2.4 GM(1,1)灰色模型預(yù)測

灰色預(yù)測模型(Grey Forecast Model)是基于一階常微分方程建立的,故稱為一階一元灰色模型,記為GM(1,1)。該方法利用已知少量、不完全的數(shù)據(jù),通過構(gòu)建數(shù)學(xué)模型,對未來的趨勢做出預(yù)測。其原理如下:

設(shè)初始非負(fù)數(shù)據(jù)序列為:

X(0)={X(0)(i),i=1,2,…,n}

(9)

式中n為序列長度。

首先對X(0)進(jìn)行一次累加,生成一次累加序列X(1):

X(1)={X(1)(k),k=1,2,…,n}

(10)

(11)

(12)

灰色序列矩陣B:

(13)

Yn=(X(0)(2),X(0)(3),…,X(0)n)T

(14)

(15)

(16)

3 結(jié)果分析與討論

3.1 礦區(qū)土壤濕度時空變化

利用ENVI軟件對2001—2021年、每2年1期的遙感影像進(jìn)行處理,計(jì)算各年土壤濕度情況并進(jìn)行分級,得到礦區(qū)2001—2021年不同等級土壤濕度分布情況,如圖2所示。

圖2 準(zhǔn)東礦區(qū)2001—2021年不同等級土壤濕度分布圖

由圖2可知,2001—2021年,準(zhǔn)東礦區(qū)內(nèi)各年土壤濕度ITVD值均為 0～1。中部和西南部地區(qū)土壤濕度較低,多為五級土壤濕度狀態(tài),呈現(xiàn)淺藍(lán)色;北部和東部區(qū)域土壤濕度較高,呈現(xiàn)深藍(lán)色;研究區(qū)內(nèi)一至三級土壤濕度面積總體變化不大,五級土壤濕度面積有縮小趨勢,四級土壤濕度面積逐年有所增加。可知,礦區(qū)內(nèi)土壤濕度總體較低,目前稍有改善,這與礦區(qū)煤炭開采過程相應(yīng)的生態(tài)環(huán)境保護(hù)措施有關(guān)。

對準(zhǔn)東礦區(qū)5種等級(一至五級)的土壤濕度面積(S1,S2,S3,S4,S5)進(jìn)行統(tǒng)計(jì),其年降水量與不同等級土壤濕度面積見表1。

表1 準(zhǔn)東礦區(qū)年降水量與不同等級土壤濕度面積

由表1可知,2001年和2003年的一級、二級土壤濕度面積均為0;2001—2021年,一級土壤濕度面積占總研究區(qū)面積<5%,其中2013年面積最大,達(dá)到350.906 km2;二級土壤濕度面積占比略小于一級土壤濕度面積,最大占比為2021年的2.18%,土壤濕度面積為264.546 km2;三級土壤濕度面積占研究區(qū)總面積0～15%,2007年三級土壤濕度面積最大,為1 702.101 km2;四級和五級土壤濕度面積在研究區(qū)內(nèi)面積最大,其面積總和占比達(dá)85%以上,2019年面積總和最大,占總面積99.76%;除2001年和2013年外,研究區(qū)內(nèi)四級土壤濕度面積始終大于五級土壤濕度面積。

準(zhǔn)東礦區(qū)不同等級土壤濕度面積與年降水量變化情況如圖3所示。

圖3 不同等級土壤濕度面積與年降水量變化情況

由圖3可知,2001—2021年準(zhǔn)東礦區(qū)年降水量總體呈波動上升趨勢,與區(qū)域氣候變化有關(guān),2019年降水量突增至1 319.78 mm,2021年繼續(xù)增加至1 496.06 mm。2001—2007年,降水量增加了164.34 mm,五級土壤濕度面積縮小,三級、四級土壤濕度面積增大,總體土壤濕度增大。其中,2007—2011年,降水量下降117.10 mm,三級、四級土壤濕度面積縮小,五級土壤濕度面積增大,總體土壤濕度減小,土壤濕度變化一定程度上受到降水因素影響;2011—2015年,降水量增加,五級土壤濕度面積先增大后縮小,四級土壤濕度面積先縮小后增大;2017—2021年,降水量有較大增長,但土壤濕度略有降低。

礦區(qū)內(nèi)四級、五級土壤濕度面積變化較為明顯,呈現(xiàn)出互為相反的變化趨勢(見圖3),且兩者面積變動始終保持在1 000～11 000 km2;一級、二級、三級土壤濕度面積處于0～2 000 km2,年際波動平緩,無明顯變化。

總體上,準(zhǔn)東礦區(qū)無明顯地表水補(bǔ)給,植被覆蓋率低,土壤濕度變化受降水等氣象因素影響較大,礦區(qū)煤炭開發(fā)活動對土壤水分含量分布具有一定影響。研究區(qū)內(nèi)四級、五級土壤濕度面積變化較大程度地表征了區(qū)域土壤濕度整體狀況,四級土壤濕度面積增大,五級土壤濕度面積減小,土壤濕度年際間變化呈現(xiàn)良好態(tài)勢。

3.2 GM(1,1)灰色模型預(yù)測分析

以2001—2021年準(zhǔn)東礦區(qū)各級土壤濕度面積為基礎(chǔ)數(shù)據(jù),采用GM(1,1)灰色預(yù)測模型進(jìn)行預(yù)測分析。設(shè)置一級土壤濕度面積為M(0)(k),二級土壤濕度面積為N(0)(k),三級土壤濕度面積為O(0)(k),四級土壤濕度面積為P(0)(k),五級土壤濕度面積為Q(0)(k),得到準(zhǔn)東礦區(qū)各級土壤濕度面積序列:

1 321.534 416

(17)

(18)

4 312.719 038

(19)

1 528 021.607 119

(20)

1 390 657.819 37

(21)

由式(17)～(21)計(jì)算可得研究區(qū)內(nèi)2001—2039年一至五級土壤濕度面積預(yù)測值。2001—2039年,準(zhǔn)東礦區(qū)一級至五級土壤濕度面積預(yù)測結(jié)果如表2～3所示。

表2 準(zhǔn)東礦區(qū)2001—2021年各級土壤濕度面積基礎(chǔ)數(shù)據(jù)與預(yù)測結(jié)果

表3 準(zhǔn)東礦區(qū)2023—2039年各級土壤濕度面積預(yù)測結(jié)果

由表2可知,四級、五級土壤濕度面積預(yù)測結(jié)果較好,最大相對誤差為1.24%;其次是一級、三級土壤濕度面積預(yù)測,一級土壤濕度面積預(yù)測相對誤差最大為27.33%,三級土壤濕度面積預(yù)測相對誤差最大為14.90%;二級土壤濕度面積預(yù)測結(jié)果存在較大誤差,與遙感影像反演中云量產(chǎn)生的影響有關(guān)。

由表2～3可得準(zhǔn)東礦區(qū)一級至五級土壤濕度面積變化趨勢,如圖4所示。

由圖4可知,2021—2039年期間,準(zhǔn)東礦區(qū)一級土壤濕度面積在0～1 000 km2內(nèi),略有增大;二級土壤濕度面積逐漸減小并趨于0;三級土壤濕度面積在0～1 000 km2,并逐年下降;四級和五級土壤濕度面積仍占據(jù)研究區(qū)主要部分,四級土壤濕度面積呈現(xiàn)逐年上升趨勢,五級土壤濕度面積穩(wěn)定在3 000 km2左右,逐年有所下降。綜上可知,準(zhǔn)東礦區(qū)未來17年總體土壤濕度呈降低趨勢。

圖4 準(zhǔn)東礦區(qū)土壤濕度面積變化趨勢

4 結(jié)論

1)2001—2021年,準(zhǔn)東礦區(qū)各年土壤濕度ITVD值均為0～1,總體土壤濕度較低,以四級、五級土壤濕度為主,占礦區(qū)總面積的85%以上;北部和東部區(qū)域土壤濕度較高,中部和西南部地區(qū)土壤濕度較低,多為五級土壤濕度狀態(tài)。

2)2001—2021年,一級、二級、三級土壤濕度面積年際變化不大,在0～2 000 km2內(nèi)波動;受降水等氣象因素影響,四級和五級土壤濕度的面積變動始終保持在1 000～11 000 km2,兩者表現(xiàn)出相反的變化趨勢。

3)基于灰色預(yù)測模型的預(yù)測結(jié)果表明:在預(yù)測期內(nèi),四級、五級土壤濕度面積仍占據(jù)研究區(qū)主要部分,四級土壤濕度面積呈現(xiàn)逐年上升趨勢,五級土壤濕度面積穩(wěn)定在3 000 km2左右,并逐年下降;在0～1 000 km2內(nèi),一級土壤濕度面積略有增大,二級土壤濕度面積逐漸減小并趨于0,三級土壤濕度面積逐年下降;準(zhǔn)東礦區(qū)未來17年土壤濕度總體呈下降趨勢。