葉偉林,周自強,石三娥,康麗娟,王得楷

(1.甘肅省科學院地質(zhì)自然災害防治研究所,甘肅 蘭州 730000;2.蘭州大學資源環(huán)境學院,甘肅 蘭州 730000;3.成都理工大學地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室,四川 成都 610059)

甘肅省黑方臺屬于大陸性干旱氣候,降水量少而蒸發(fā)量大,土壤含水量低。但由于20世紀60年代人口遷入,對該區(qū)進行了大量的農(nóng)業(yè)灌溉,導致干旱的區(qū)域土壤含水量增高。同時長期灌溉使得地下水位上升,臺塬邊黃土底部逐漸飽水軟化,引起土體抗剪強度下降,導致黑方臺穩(wěn)定性降低,致使其邊緣成為滑坡的易發(fā)區(qū)[1-2]。據(jù)統(tǒng)計70年代以來黑方臺發(fā)生大小滑坡約150次,造成不同程度的人員傷亡和經(jīng)濟損失,引起了社會的廣泛關注。此外,土壤中的水分是滑坡的主要誘因[3],它是地表層和大氣能量交換過程中的重要參與者,決定著植被的蒸發(fā)以及光合作用,并且對陸地表面蒸發(fā)、水的運移以及碳循環(huán)等有很強的控制作用,對于研究預測區(qū)域干濕情況有著舉足輕重的作用[4-5]。因此,監(jiān)測黑方臺土壤濕度分布及其變化對分析灌溉等因素對研究區(qū)土壤濕度的影響以及灌溉對滑坡的影響具有一定的意義。

土壤濕度的傳統(tǒng)獲取方法主要包括質(zhì)量法、中子儀法在內(nèi)的田間實測法[6],其采用單點進行實地監(jiān)測,雖然可以直接準確地測量土壤剖面的含水量,但是采樣點數(shù)目太多,需要大量的人力物力,不僅費時而且成本高,很難大范圍、高效率地獲取土壤水分。此外,區(qū)域在土壤、地形、植被覆蓋上的空間差異使單點的代表性差,只能夠反映出采樣點局部區(qū)域的土壤濕度。因此,田間實測法在大范圍監(jiān)測土壤濕度的過程中存在很大的局限性[7-9]。遙感手段獲取土壤濕度相對于傳統(tǒng)的手段具有范圍廣、速度快、成本低、可多期監(jiān)測等特點,分為微波遙感和光學遙感,包括熱慣量法[7-9]、蒸散量法[10]、作物缺水指數(shù)法[10-12]、溫度植被干旱指數(shù)(TVDI,temperature vegetation dryness index)法等[13-15]。TVDI通過構建地表溫度(LST,land surface temperature)和歸一化植被指數(shù)(NDVI,normalized difference vegetation index)的對應關系來反映土壤濕度狀況。Price[16]于1990年發(fā)現(xiàn)LST與NDVI呈負相關關系,以它們?yōu)闄M縱軸構成的散點區(qū)域呈三角形,在此基礎上,Sandholt等[17]通過總結前人研究,首次提出了TVDI來監(jiān)測表層土壤水分狀況,獲得了較好的效果;趙杰鵬等[18]通過改進TVDI反演模型,較精確地估算了新疆地區(qū)的土壤水分狀況;王純枝等[19]利用MODIS產(chǎn)品數(shù)據(jù)NDVI和Ts構建了特征空間計算得到TVDI,經(jīng)驗證TVDI能準確的反映黃淮海平原土壤10～20 cm處的土壤水分狀況;王行漢等[20]利用增強溫度植被指數(shù)(ETVDI,enhanced temperature vegetation dryness index)和地表溫度反演TVDI,克服了TVDI在高植被覆蓋區(qū)植被指數(shù)飽和的缺陷,為南部高植被覆蓋區(qū)的土壤水分監(jiān)測提供了一種新的方法。

眾多學者采用TVDI方法對不同區(qū)域土壤濕度進行了大量的研究,得到了很好的效果。黑方臺為低植被區(qū)域,用TVDI來研究土壤濕度效果更好。基于此,利用2019年11月遙感數(shù)據(jù)反演TVDI分布狀況,之后通過相同日期的采樣數(shù)據(jù)建立TVDI與不同深度處土壤濕度的關系,分析TVDI與不同土壤深度關系的可靠性和準確性。在此基礎上,分析黑方臺土壤濕度空間分布以及動態(tài)變化,探討灌溉等因素對黑方臺土壤濕度的影響以及灌溉對滑坡的影響。

1 研究地區(qū)與研究方法

1.1 研究區(qū)概況

甘肅黑方臺地處西北黃土高原,為黃河Ⅳ級階地臺塬,海拔1 700 m左右,由黑臺和方臺組成,屬于溫帶大陸性氣候,全年降水量少,晝夜溫差大,日照時間長,四季分明。多年平均降水量為287.6 mm,年蒸發(fā)量達1 600 mm,最大降水量為431.9 mm,最小降水量為178.8 mm,年蒸發(fā)量約為年降水量的5.4倍[21-22]。20世紀60年代,由于劉家峽與鹽鍋峽水庫的修建,大量的人員移居在此,為解決農(nóng)業(yè)缺水問題,自1960年開始,大量黃河水被抽至臺塬進行大水漫灌,長年灌溉使黃土底部形成了數(shù)十米厚的地下水位。灌溉水的入滲使得土體含水量增加,飽和度增大,引起土體抗剪強度下降,從而導致黑方臺穩(wěn)定性降低,容易發(fā)生滑坡[22]。

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

(1) 遙感數(shù)據(jù) 選用2017—2019年間26景LandsatOLI/TIRS遙感影像(影像獲取于https://www.usgs.gov/)。遙感影像預處理具體如下:首先利用Envi5.3所提供的輻射定標工具(radiometric calibration)對Band10、Band4、Band5進行輻射定標,將原始的DN值轉化為輻射亮度值,后利用FLAASH模型進行大氣校正,減少大氣散射、吸收、反射引起的誤差,得到真實的地表反照率,再對數(shù)據(jù)進行裁剪等處理。

(2) 野外采樣數(shù)據(jù) 為了精確提取土壤濕度,在遙感影像獲取當天進行了土壤數(shù)據(jù)采集,土壤濕度在該時間段內(nèi)保持在穩(wěn)定狀態(tài)。采樣點分布如圖1所示,12個采樣點分別分布在不同的土地利用類型和典型的滑坡體上。每個采樣點分別在10 cm、20 cm、30 cm、40 cm、50 cm的埋深處采樣,最終采集到5個深度的土壤質(zhì)量水。

圖1 研究區(qū)概況及采樣點分布Fig.1 General data and distribution of sampling points of study area

(3)氣溫降水數(shù)據(jù)和其他數(shù)據(jù) 從中國氣象數(shù)據(jù)網(wǎng)(https://data.cma.cn/)獲取了永靖站點2017—2018年月平均氣溫、月降水、月日照時數(shù)、月相對濕度等氣象數(shù)據(jù);由于永靖站缺少月蒸發(fā)量數(shù)據(jù),其經(jīng)緯度為(103°18′E,35°58′N),而榆中縣、民和縣經(jīng)緯度分別為(104°09′E,35°52′N)、(102°5′E,36°2′N),發(fā)現(xiàn)永靖縣恰處于中間且離它們最近,因此利用榆中縣和民和縣的月蒸發(fā)量均值來代替永靖縣;每日的灌溉量數(shù)據(jù)來自于永靖縣鹽鍋峽鎮(zhèn)水管所。

1.3 研究方法

(1) 歸一化植被指數(shù) 基于地表植被對紅光波段和近紅外波段不同的吸收和反射能力,通過比值運算可以有效地反映土壤背景信息和地表植被覆蓋狀況。生長較茂盛的植被在近紅外和紅外波段內(nèi)表現(xiàn)出強度的反射差異,而生長不好的植被反射差異較少[23]。植被的生長狀況以及分布范圍能被NDVI有效提取,因此NDVI被廣泛應用于自然環(huán)境中植被信息的提取,其反演公式為

(1)

其中:ρ1、ρ2分別為Landsat8近紅外和紅外波段的反射率。

(2) 地表溫度 應用Landsat熱紅外波段反演地表溫度主要有3種方法:大氣校正法、單窗算法、單通道算法?；诖髿庑Ｕ?利用Landsat8 TIRS反演地表溫度[24-25],反演公式為

Lλ=[εB(TS)+(1-ε)L↓]τ+L↑,

(2)

(3)

(4)

其中:Lλ是衛(wèi)星傳感器接收到的熱紅外輻射亮度值;L↑是大氣向上的輻射亮度;L↓是大氣向下的輻射亮度;ε是地表比輻射率;B(TS)是黑體的輻射亮度;TS是地表真實的溫度;對于TIRS Band10,K1=774.89 W/(m2×μm×sr),K2=1 321.08 W/(m2×μm×sr)。

(3) 溫度植被干旱指數(shù) Price[16]發(fā)現(xiàn),當研究區(qū)植被覆蓋滿足從裸地到完全植被覆蓋時,植被指數(shù)和LST構成的散點圖呈現(xiàn)出三角形。Sandholt 等[17]基于NDVI和LST建立了LSTi-NDVI三角形特征空間,提出了溫度植被干旱指數(shù)的概念TVDI,并證明地表溫度LSTi與NDVI存在明顯的負相關關系,任意一個NDVI值對應唯一一組LSTmax(干邊)和LSTmin(濕邊),因而LSTi-NDVI特征空間的擬合線斜率可以反映區(qū)域土壤濕度的情況(見圖2)。TVDI計算公式為

圖2 LSTi-NDVI的特征空間Fig.2 The feature space of LSTi-NDVI

(5)

LSTmin=a1-b1NDVI,
LSTmax=a2-b2NDVI,

(6)

(7)

其中:LSTmax(干邊)和LSTmin(濕邊)分別表示當NDVI等于某一特定值時,地表溫度的最大值和最小值,即特征空間的干邊和濕邊;a1、b1、a2、b2是干邊、濕邊擬合系數(shù);LSTi表示任一像元地表溫度。

(4) 偏差分析 偏差為每月TVDI值距多月平均TVDI的距離,可表示某一時段內(nèi)TVDI偏離多月TVDI均值的程度,其中正值表示高于多年平均水平,反之則表示低于多年平均水平[26],計算公式為

Di=TVDIi-TVDIm,

(8)

其中:Di為第i月偏差值;TVDIi為第i月TVDI值;TVDIm為多月平均TVDI;i為月份。

(5) 線性傾向率 利用線性傾向率分析2017—2019年TVDI的月際變化趨勢,若為正值則表示TVDI呈現(xiàn)增加趨勢,否則表示下降趨勢?；谙裨木€性傾向率計算公式為

(9)

其中:Bslope為線性傾向值;TVDIi為第i月的TVDI值;i為月份變量(i=1,2,3,…,26),n=26。當Bslope>0時,表示隨時間增加TVDI呈上升趨勢,反之,TVDI呈下降趨勢,Bslope的大小表示TVDI上升或下降的速率。

(6) 相關分析 為了研究黑方臺土壤濕度狀況與自然和人為影響因素的線性關系,采用卡爾·皮爾森(KarlPearson)在1880年提出的Pearson相關系數(shù),該系數(shù)能檢驗出變量因子間相關性的強弱和方向,已廣泛應用于各學科領域,其計算公式為

(10)

其中:rxy為相關系數(shù),值域范圍是(-1,1),值大于零表示正相關,值小于零表示負相關,其絕對值越大表示相關性越顯著;xi為第i月的平均TVDI值;yi為當月相關變量的均值。

2 結果與分析

2.1 土壤濕度監(jiān)測結果檢驗

國內(nèi)學者比較了各溫度植被干旱指數(shù)與土壤濕度之間的關系,結果表明TVDI相對于其他指數(shù)能更有效反映土壤濕度變化[27-28]。為了驗證研究區(qū)TVDI與表層土壤濕度的相關性,選取12個野外采樣點,得到土壤不同深度60個實測土壤濕度數(shù)據(jù),將其分別與樣點處TVDI進行線性回歸(見圖3)。發(fā)現(xiàn)10 cm土壤深度處(見圖3(a)),P=0.002,R2=0.63;20 cm土壤深度處(見圖3(b)),P=0,R2=0.79,相關性最好;30 cm土壤深度處(見圖3(c)),P=0.001,R2=0.68,相關性次之;在40 cm(見圖3(d))、50 cm(見圖3(e))處隨著土壤深度的加深,相關性呈減小趨勢,P值分別為0.007、0.028。結果表明研究區(qū)TVDI和實測土壤濕度在20 cm、30 cm土壤深度處呈較強的負相關關系,而在40 cm、50 cm處隨土壤深度逐漸加深,TVDI指數(shù)和土壤濕度的相關性減小。

圖3 土壤濕度和TVDI線性擬合Fig.3 Linear fit of soil moisture and TVDI

2.2 土壤濕度空間分布特征

為了分析黑方臺土壤濕度的空間分布狀況,對2017—2019年26個月的TVDI均值進行統(tǒng)計并按自然斷點法(natural breaks)將其分為5個等級(見圖4(a)),1～5等級月均TVDI的范圍值分別為0.05～0.32、0.32～0.50、0.50～0.61、0.61～0.74、0.74～0.99,并對各等級的面積比進行了統(tǒng)計分析(見圖4(b)),其中等級越高代表土壤越干旱,土壤濕度越低。黑方臺3等級占研究區(qū)面積最大,為32.44%,說明對研究區(qū)而言,整體土壤濕度較高,干旱程度較低。等級4和等級5的區(qū)域主要分布在黑方臺臺塬邊周圍,占總區(qū)域的35.79%,這些區(qū)域多為裸地,沒有農(nóng)業(yè)灌溉,加之干旱的氣候導致該區(qū)域在20 cm、30 cm處土壤濕度非常低。臺塬內(nèi)部土壤濕度以2、3等級為主,占總區(qū)域的58.5%,且4等級呈零星狀分布其中。由于靠近黃河,離水源近,臺塬內(nèi)部以農(nóng)用水澆地為主,包括耕地、林地等,因此土壤濕度較臺塬邊緣高。綜上說明黑方臺土壤含水量呈現(xiàn)臺塬內(nèi)高,臺塬邊低的空間分布格局,這主要是受人類活動農(nóng)業(yè)灌溉的影響。

圖4 月均TVDI的空間分布和各等級面積占比Fig.4 The spatial distribution of monthly TVDI mean and the area ratio of each grade

2.3 土壤濕度空間變化特征

(1) 時間變化特征 從黑方臺各月TVDI均值來看(見圖5),干旱月份較少,土壤濕度較高。TVDI值呈現(xiàn)不規(guī)則的波動變化趨勢,峰值出現(xiàn)在2018年7月和8月,分別為0.67和0.65,超過平均值0.57;低值區(qū)出現(xiàn)在2017年11月、12月,2019年1月為0.47,低于平均值,因為此時研究區(qū)有一定量的灌溉且蒸發(fā)量較少。各月TVDI均值在時間尺度上呈微小的下降趨勢,表明土壤濕度在該時間尺度上呈微小的上升趨勢。應用偏差分析得到2017年3月—2019年11月各月TVDI偏離多月平均水平的程度,TVDI的偏差值表現(xiàn)的特征為在每年3—8月值為正,表現(xiàn)為增加的趨勢,而在每年的10月至次年1月值為負,呈減小的趨勢。據(jù)黑方臺灌溉調(diào)查資料顯示,每年春末、秋初及整個夏季黑方臺灌溉量增大,土壤濕度會隨著灌溉量的增大而增大,但由于研究區(qū)隨溫度升高,蒸發(fā)量急劇增大,導致這些月份土壤濕度反而較小。而在冬季等氣溫較冷的月份,蒸發(fā)量比夏季小,加之有少量的灌溉,導致研究區(qū)臺塬土壤濕度較高。

圖5 2017年3月—2019年11月 TVDI 月變化Fig.5 Monthly variations of TVDI during March 2017 to November 2019

(2) 空間變化特征 為了進一步探討研究區(qū)土壤濕度的空間分異特征,采用最小二乘法,逐像元對研究區(qū)年均TVDI進行回歸分析,結果見圖6。由圖6看出,整個黑方臺土壤濕度表現(xiàn)出明顯的空間分異格局,其中2017年3月—2019年11月TVDI上升的區(qū)域主要分布在黑方臺的西北區(qū)域,說明該時間段內(nèi)黑方臺的土壤濕度呈減小的趨勢。在臺塬內(nèi)部發(fā)現(xiàn),黑臺中部的小片區(qū)域內(nèi),變化趨勢小于零,說明在該時間段內(nèi)TVDI越來越小,土壤濕度在變大。此外,臺塬上的其他區(qū)域變化率不太顯著,說明臺塬上大部分區(qū)域在2017年3月—2019年11月土壤濕度沒有發(fā)生明顯變化。這主要是由于黑方臺區(qū)域土壤濕度在長時間尺度上受人類定時灌溉活動的影響,其臺塬上TVDI基本保持不變。

圖6 2017—2019年TVDI線性傾向率Fig.6 The tendency rate of TVDI during 2017-2019

2.4 土壤濕度影響因素

黑方臺年降水少而蒸發(fā)大,農(nóng)業(yè)用地均為水澆地,鑒于黑方臺以上環(huán)境條件,從自然因素和人為因素2個方面選取月平均氣溫(℃)、降水(mm)、日照時數(shù)(h)、蒸發(fā)量(mm)、相對濕度(%)和月灌溉量(m3)6個因子來分析其對土壤濕度的影響情況。

(1) 自然因素 通過對月平均氣溫、月降水量、月日照時數(shù)、月蒸發(fā)量、月相對濕度與月均TVDI進行Pearson統(tǒng)計(見表1),發(fā)現(xiàn)除相對濕度與TVDI呈不顯著負相關外,其余因子均與TVDI呈正相關。從Pearson相關系數(shù)|R|來看,各影響因子與TVDI關系密切程度依次為月蒸發(fā)量>月平均氣溫>月降水量>月日照時數(shù)>月相對濕度。結果表明土壤濕度與月平均氣溫、月蒸發(fā)量、月降水量、月日照時數(shù)均在0.01水平上呈負相關,相關系數(shù)依次為-0.789、-0.733、-0.621、-0.557。黑方臺2017—2019年氣象數(shù)據(jù)見圖7。由于研究區(qū)年降水量稀少且集中在夏季,與此同時研究區(qū)溫度也急劇升高,導致夏季蒸發(fā)量增高而土壤濕度呈降低趨勢。

表1 月均TVDI與各因素相關系數(shù)分析

圖7 研究區(qū)2017年1月—2019年10月氣象數(shù)據(jù)Fig.7 The meteorological data during January 2017 to October 2019

(2) 人為因素 黑方臺屬于大陸性干旱氣候,年降水量較少,而蒸發(fā)量巨大,為了深入分析黑方臺農(nóng)用灌溉量對土壤含水量變化趨勢的影響情況,獲取2017—2019年該區(qū)域月灌溉量(見圖8),通過計算每月的累計量來研究月均TVDI與灌溉量的Pearson關系,其相關系數(shù)為0.409,呈不顯著正相關。研究表明在年際內(nèi)隨著月份的增大,農(nóng)業(yè)需水量增多,其灌溉量也增多,灌溉水不僅滿足地表水,隨著淺層土壤水飽和而入滲成為地下水,造成該區(qū)域成為滑坡的高發(fā)區(qū);與此同時大氣溫度逐漸升高,土壤溫度也在升高,研究區(qū)蒸發(fā)量快速上升,導致在溫度較高的季節(jié)隨灌溉量的增多,TVDI值增大,而土壤中的含水率減小。

圖8 研究區(qū)2017年3月—2019年11月灌溉量Fig.8 The irrigation volume during March 2017 to November 2019

3 討論與結論

3.1 討論

研究黑方臺土壤濕度空間分布和時間變化發(fā)現(xiàn),該區(qū)域年蒸發(fā)總量達到1 000 mm以上,成為繼農(nóng)業(yè)灌溉量后土壤濕度的第二大影響因素。由于研究區(qū)長期的農(nóng)業(yè)灌溉引起地下水位上升,造成土壤抗剪強度下降,加之西高東低的地勢導致研究區(qū)東南方向成為滑坡的高發(fā)區(qū)。許多學者對研究區(qū)灌溉引起的滑坡機理、過程等進行了大量的研究,張茂省等[29-30]研究了黑方臺地表水的入滲、灌溉與地下水的演化關系;谷天峰等[31]研究了黑方臺地下水位上升引起的黃土斜坡內(nèi)部滲流場的變化特征,探討了地下水位變化對斜坡穩(wěn)定性的影響;亓星等[1]通過研究黑方臺地下水位與臺塬灌溉的響應規(guī)律,發(fā)現(xiàn)地下水具有133 d的滯后期;嚴冬冬[32]研究了灌溉引起的黑方臺區(qū)域地下水的演化過程以及水位上升速度;朱立峰[33]對黑方臺滑坡群密集發(fā)育的內(nèi)在控制條件和外動力誘發(fā)因素進行了研究。

綜上,這些研究都忽略了研究區(qū)氣候因素——高蒸發(fā)量與滑坡的關系。上述研究表明,農(nóng)業(yè)灌溉是黑方臺土壤濕度上升的主要影響因素,而高蒸發(fā)量是土壤濕度降低的重要因素,說明用于農(nóng)業(yè)灌溉的一部分水被蒸發(fā)損失。因此在考慮灌溉量對滑坡影響的研究中,蒸發(fā)量應作為不可忽視的一部分水分被考慮在黑方臺滑坡影響因素內(nèi)。在研究其他區(qū)域滑坡影響以及機理的過程中,研究區(qū)諸如蒸發(fā)量等氣候因子是不容忽視的重要因素。

3.2 結論

研究利用多源遙感數(shù)據(jù),構建指示生態(tài)干旱強弱的溫度植被干旱指數(shù),從而研究整個流域的土壤濕度變化特征,探討其對氣候變化的響應。主要結論如下:

(1) 從2017年3月—2019年11月TVDI月均值的空間分布來看,黑方臺土壤濕度表現(xiàn)出臺塬內(nèi)高,邊緣低的明顯空間差異,這主要受研究區(qū)人為干擾的影響,臺塬內(nèi)多為可利用的水澆地,而臺塬邊緣由于經(jīng)常發(fā)生滑坡,多為廢棄的未利用地。TVDI月均值為0.57,說明對降水稀少的研究區(qū)而言,灌溉對土壤濕度影響很大。

(2) 從TVDI月均值的時間變化特征來看,年際土壤濕度存在減—增—減的變化態(tài)勢,主要受蒸發(fā)量及灌溉量的影響,總體時間尺度上呈現(xiàn)土壤濕度增加的態(tài)勢。

(3) 從黑方臺土壤濕度的空間變化特征來看,臺塬內(nèi)部有小部分土壤濕度增大的區(qū)域,而大部分區(qū)域保持土壤濕度不變,此外臺塬邊緣土壤濕度也基本保持不變,主要受到人類活動的影響。土壤濕度減小的區(qū)域多分布在黑方臺西北方向。

(4) 從土壤濕度的多個影響因子來看,月平均氣溫與其負相關關系最為顯著,也進一步驗證了研究區(qū)土壤濕度的大小主要受蒸發(fā)量和灌溉量影響,而空間分布格局受灌溉影響。