基于TVDI 的旱情時空動態變化監測
——以神東礦區為例

邵天意，包斯琴*，王 楠，韓阿茹汗

（1.內蒙古農業大學 沙漠治理學院，呼和浩特 010020；2.中國地質大學（北京）土地科學技術學院，北京 100089）

0 引 言

【研究意義】干旱是自然生態系統產量的重要輸入因子，同時也是植被發揮蒸騰作用及光合作用的限制因子，其在一定程度上能夠決定一個地區的植被類型和植被生長結構[1]，更是衡量荒漠化程度以及指導生態恢復措施實施的重要影響因子，對準確且快速的掌握旱情狀況對土地復墾、生態修復和農業生產具有指導意義[2]。【研究進展】傳統長時間序列的旱情監測方法需要消耗巨大的人力、物力和財力。相比于傳統的監測方法，遙感監測具有監測范圍廣、時效性長等特點[3]。諸多學者[4-6]針對遙感監測領域進行的大量研究表明，表觀熱慣量模型進行干旱遙感監測需要地表反照率和太陽輻射作為參數，其更適合于裸土的研究區，同時還需考慮土壤類型變化對監測結果的影響。基于冠層溫度的遙感監測方法在實際監測中易受自然氣候條件，衛星傳感器等影響會產生不確定的誤差[7-8]。綜合作物長勢進行干旱監測以植被長勢狀態表征干旱程度存在相對大的滯后性，預警時效較差。微波遙感監測旱情對衛星影像數據的精度要求極高[9-10]。綜上所述，利用研究區歸一化植被指數（Normalized Difference Vegetation Index，NDVI）與地表溫度（Land Surface Temperature，LST）建立溫度干旱植被指數（Temperature Drought Vegetation Index，TVDI）干旱監測模型能夠消除了單一因素對監測結果的影響，更具時效性，且該模型監測適用于中高植被覆蓋地區[11]。TVDI模型是一種采用光學和熱紅外遙感通道數據進行地表水分的遙感反演方法[12]。對于植被覆蓋區而言，地表水分一定程度上決定著植被的冠層溫度，在一定條件下，植被的冠層溫度能間接反映植被的供水狀況[13-15]。據此，Sandholt等[16]基于NDVI與LST建立了Ts-NDVI特征空間，計算得到的TVDI能夠很好地反演出地表水分，其反演精度也得到了驗證[17]。諸多學者[18-20]也基于不同數據源利用TVDI模型進行地表水分反演并進行精度驗證，發現TVDI在0～10 cm 土層深度具有良好的反演效果。【切入點】國內基于TVDI模型進行地表水分動態監測的研究很多，但針對地勢復雜溝壑縱橫的礦區進行的研究尚不多見。神東礦區受地理位置的影響干旱程度高，加之礦區內常年進行煤炭開采作業，導致地表結構被破壞，土壤水分流失嚴重，干旱程度加劇。【擬解決的關鍵問題】本文基于TVDI模型，利用神東礦區1991—2018 年植被生長季遙感數據進行干旱程度動態監測，以期為研究區的生態恢復和土地復墾提供指導建議。

1 材料與方法

1.1 研究區概況

神東礦區位于陜西省榆林市北部和內蒙古自治區鄂爾多斯市南部（109°51′—110°46′E，38°52′—39°41′N），區內大部分為典型的風成沙丘及沙灘地地貌。海拔800～1 385 m，溝壑縱橫，地表支離破碎，水土流失十分嚴重。神東礦區地處毛烏素沙地與黃土高原的過渡地帶，屬半干旱大陸性季風氣候，溫度年較差較大，年平均氣溫6.6 ℃，1 月平均氣溫-10.1 ℃，7 月平均氣溫20.5 ℃。區內干旱少雨，降水集中，汛期為6—9 月，占全年降水量的76%，降水年際變化較大，最大年達919.10 mm，最小量年僅為108.60 mm。區內以沙生植物為主，植被每年4 月返青，10月葉落，郁閉較差。區內風沙土（Aeolian sandy soil）占礦區總面積的50%，黃土性土（Loess soil）占總面積的30%，紅土性土（Red native soil）占總面積的10%，土壤具有質地較粗，結構不良，肥力低，抗侵蝕性差的特點。研究區位置如圖1 所示

圖1 研究區位置示意圖Fig.1 Schematic diagram of study area location

1.2 數據來源

遙感影像數據來源于地理空間數據云（http://www.gscloud.cn/）中的美國國家航空航天局（National Aeronautics and Space Administration,NASA）發射的陸地系列衛星（Landsat）所搭載的30 m 空間分辨率數據，包括 Landsat5-TM 和Landsat8-OLI，時間分辨率為16 d 的2 種傳感器。數據時相選取神東礦區1991、2002、2007、2010、2014、2018 年植被生長季（5—10 月）共36 景影像，在1991—2018 年數據質量上均選擇云量小于10%反演效果好的數據單元，借助ENVI 5.3 進行輻射定標、大氣校正、鑲嵌等預處理，然后利用數據估算得到空間分辨率30 m 的NDVI及LST數據集。

1.3 研究方法

1.3.1TVDI計算方法

以提取得到研究區的歸一化植被指數（NDVI）為橫坐標，地表溫度最大值和最小值為縱坐標進行線性擬合，形成關于Ts-NDVI的三角形特征空間，將該特征空間的斜率作為TVDI的等值線，TVDI自下而上升高，TVDI的絕對值越大，干旱程度越大。因此，在研究區擬合出其特征空間的干濕邊方程，即可得到每個像元的干旱指數，其計算式為：

式中：TVDI為溫度干旱植被指數；Ts為地表真實溫度；Ts_max為地表最高溫度；Ts_min為地表最低溫度。TVDI值越大，代表地表水分越低，干旱程度越大。TVDI的干濕邊擬合方程為：

式中：a1、b1、a2、b2分別為干濕邊線性擬合方程的系數。

在干濕邊方程中，其斜率在生態學意義上代表著土壤含水率在飽和與不足時的地表溫度值，將各像元的NDVI值帶入所對應的干濕邊方程，計算出最高地表溫度（Ts_max）與最低地表溫度（Ts_min），利用式（1）計算各像元對應的TVDI值。本文采取適用干旱半干旱區的TVDI分級標準進行分級，如表1 所示[21]。

表1 TVDI 分等定級Table 1 TVDI classification and classification

1.3.2 偏差分析法

本文采用偏差分析法[22]分析研究區1991—2018年植被生長季的TVDI時間分布特征。偏差分析法可以在時間尺度上通過某年TVDI與多年平均的TVDI距離，表示某一時期TVDI偏離多年平均TVDI的程度。TVDI偏離值為正值，表示年均TVDI高于多月平均水平，TVDI在時間序列上呈上升趨勢，TVDI偏離值為負值，表示年均TVDI低于多年平均水平，TVDI在時間序列上呈下降趨勢。

1.3.3 Mann-Kendall 檢驗法[23]

Mann-Kendall 檢驗法（M-K 檢驗法）是世界氣象組織推薦并廣泛使用的非參數檢驗方法。其不需要服從特定的分布，亦不受樣本值的影響，因此廣泛應用于分析一般數據在時間序列上的趨勢檢驗和突變點檢驗。對于時間變量（x1,x2, …,xn），n為時間序列的長度，M-K 法定義了其統計量S：

式中：Sgn（）為函數符號，規則如下：

S為正態分布，其均值為 0 ，方差Var(s)=(n-1)(2n+5)/18，當n＞10 時，正態分布統計量計算為：

若Z＞0，表明TVDI在時間序列上呈上升趨勢，若Z＜0，反之，絕對值越大，趨勢越明顯；且Z≥1.28、1.96、2.32 時，分別通過了90%、95%、99%水平的信度檢驗，在該水平上顯著。UF曲線與UB曲線的交點，為TVDI在時間序列上的突變點。

1.3.4 趨勢分析和F檢驗

線性傾向率（Bslope）能夠反映TVDI在空間上隨時間變化的的上升或下降趨勢。因此，本文逐像元計算1991、2002、2007、2010、2014、2018 年的TVDI均值，得到相應年份的TVDI空間分布，并利用線性傾向率計算研究區1991—2018 年植被生長季的TVDI在空間上的線性變化趨勢。其計算式為：

式中：Bslope為線性傾向趨勢；i為36 個月變量（i=1,2, 3, …, 36），n=36，TVDIi為第i月的TVDI值。若Bslope＞0，表示干旱程度在空間上隨時間變化呈上升趨勢；Bslope＜0，表示研究區干旱程度在空間上隨時間變化呈下降趨勢。用F檢驗對TVDI的空間變化趨勢進行顯著性分析，根據F檢驗顯著性臨界值（α=0.05），將F值劃分為顯著下降、顯著上升、不顯著下降、不顯著上升4 種趨勢。

1.3.5 空間轉移矩陣

空間轉移矩陣來源于系統分析中對系統狀態與狀態轉移的定量描述，可以定量識別TVDI不同等級在某一時間間隔的空間格局變化，不僅可以反映TVDI不同等級的面積變化，還可以直觀反映TVDI各等級面積轉入轉出情況[25]。

2 結果與分析

2.1 神東礦區植被生長季TVDI 時間變化特征

2.1.1TVDI月際變化特征

利用ArcGis10.7 像元統計工具，統計TVDI值，得到神東礦區1991—2018年植被生長季TVDI月際變化如圖2 所示。由圖2 可知，神東礦區1991—2018年TVDI月際變化非常明顯，月均TVDI在0.35～0.85之間波動。2007 年6 月TVDI均值最高為0.85，說明地表干旱程度最大。2018年9月TVDI均值最低為0.39，說明地表干旱程度最小。月均TVDI在每年的6—9 月出現最低值，這是由于神東礦區降水集中導致的，與神東礦區6—9 月為汛期的自然條件相吻合。

圖2 神東礦區1991—2018 年植被生長季TVDI 月際變化Fig.2 Intermonthly TVDI variation of vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.2TVDI突變分析

采用M-K 突變檢驗法對神東礦區1991—2018 年植被生長季的月均TVDI進行檢驗，如圖3 所示。圖3 檢驗結果表明，UF～UB曲線幾乎處于顯著性水平α=0.05 置信區間內，TVDI值呈下降趨勢（UF＜0），與神東礦區月際變化規律相吻合。1991 年7—10 月存在4 個突變點，2014 年9 月存在1 個上升趨勢突變點，2018 年5 月、7 月左右存在一個上升趨勢、一個下降趨勢突變點。

圖3 1991—2018 年神東礦區植被生長季TVDI M-K 檢驗Fig.3 TVDI M-K test of vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.3TVDI年際變化特征

計算得到神東礦區各年植被生長季TVDI均值如圖4 所示。從年際變化來開，年均值波動于0.56～0.69之間，按照地表干旱等級類型劃分，神東礦區干旱程度基本處于正常、較干旱等級，這與TVDI的空間分布相吻合。在研究時段內，2002 年TVDI達到峰值為0.686，自2002 年后TVDI逐年下降，2018 年達到最低的0.561，下降速率為0.02/30 a。該結果在實地調研中得到了驗證，這是由于神東礦區在煤炭開采后進行了大量的土地復墾工作，水土流失強度減弱，種植大量的防風固沙植被，神東礦區地表涵養水源的能力得到提升，因此地表水分上升，地表干旱程度呈下降趨勢。

圖4 神東礦區1991—2018 年植被生長季TVDI 年際變化Fig.4 Interannual TVDI variation during the vegetation growing season in Shendong Mining Area from 1991 to 2018

2.1.4TVDI偏差分析

利用偏差分析法分析神東礦區1991—2018 年年際TVDI變化趨勢，如圖5 所示。神東礦區1991、2002、2007、2010 年TVDI偏離值＞0，說明1991—2010 年TVDI均值呈上升趨勢，干旱程度具有增大趨勢。其中1991—2002 年TVDI偏離程度最大，上升趨勢最明顯，高于TVDI均值0.05。2014、2018 年TVDI偏離值＜0，說明2010—2018 年TVDI呈下降趨勢，干旱程度減小。2014—2018 年偏離程度最大為0.56，低于TVDI均值（0.08）。該結果與神東礦區TVDI年際變化趨勢一致。

圖5 神東礦區1991—2018 年TVDI 變化趨勢Fig.5 TVDI trend of Shendong Mining Area from 1991 to 2018

2.2 神東礦區TVDI 空間變化特征

2.2.1 神東礦區植被生長季TVDI干旱等級空間分布

神東礦區1991—2018 年植被生長季干旱等級的空間分布如圖6 所示。1991—2002 年研究區西部干旱等級下降，東部區域的干旱等級上升。2002—2007年研究區中部、南部干旱等級上升。2007—2010 年研究區中部、西部干旱等級基本達到正常級別，東部地區干旱等級上升。2010—2014 年東部區域干旱等級下降，西部區域干旱等級略有提升。2014—2018年研究區西南部干旱等級增大，中部及東北部的干旱等級有小范圍的下降。從干旱等級劃分來看，神東礦區干旱等級始終表現為西南高于東北，其主要原因為神東礦區的降水量從西南到東北呈遞減趨勢，空間上基本處于正常、較干旱和干旱3 種等級。

圖6 神東礦區1991—2018 年植被生長季干旱等級的空間分布Fig.6 Spatial distribution of drought grades in the vegetation growing season from 1991 to 2018 in Shendong Mining Area

2.2.2 神東礦區線性傾向趨勢分析

利用ArcGIS10.7，基于神東礦區1991—2018 年植被生長季TVDI的空間分布，運用式（7）計算其線性傾向率（Bslope），分析其線性傾向趨勢，并對其線性傾向趨勢空間分布在α=0.05 顯著性水平上進行F檢驗，如圖7 所示。由圖7 可知，神東礦區1991—2018年TVDI線性變化趨勢下降趨勢面積大于上升趨勢面積，Bslope＜0 的區域集中在研究區的大部分地區，該區域TVDI呈下降趨勢，地表干旱程度變小。Bslope＞0的區域零散分布于研究區的北部、東部、南部，這些區域TVDI呈上升趨勢。地表干旱程度變大。根據F檢驗的空間分布得知，在TVDI呈下降趨勢的區域中，大部分地區下降趨勢顯著，下降趨勢不顯著的地帶則零星分布于其中。在TVDI呈上升趨勢的區域中，上升顯著的地區大部分集中在神東礦區的東北部，該區域內分布著神東礦區的礦井群，這是由于采煤后塌陷，造成植被破壞嚴重，地表水分下滲速度加快，地表水分減小，所以干旱程度增大。神東礦區TVDI空間變化結果與時間變化的趨勢一致

圖7 1991—2018 年神東礦區植被生長季TVDI 線性變化趨勢空間分布及F 檢驗Fig.7 Spatial distribution and F-test of TVDI linear trend in vegetation growing season in Shendong Mining Area from 1991 to 2018

通過對神東礦區TVDI線性傾向趨勢面進行計算發現：神東礦區1991—2018 年植被生長季TVDI下降趨勢面積為3 788.43 km2，上升趨勢面積為1 215.33 km2，下降趨勢面積大于上升趨勢面積。從F顯著性檢驗結果來看，神東礦區顯著下降面積最大為2 394.74 km2，顯著上升面積最小為425.91 km2，說明研究區的地表干旱程度得到緩解，該地區生態恢復效果顯著。

2.2.3 神東礦區TVDI轉移矩陣

利用ArcGIS10.7 空間分析工具，選擇1991、2002、2010、2018 年的TVDI干旱等級分布進行空間轉移矩陣的計算，得到神東礦區1991—2018 年植被生長季干旱等級轉移情況，如圖8 所示。由圖8 可知，1991—2002 年，整個研究區的西北、西南、東南大范圍區域由干旱變化為其他4 種干旱等級。較濕潤、正常、較干旱變為干旱零星分布于中東部，該部分區域旱情有所加重，其他區域基本保持原有的干旱等級。2002—2010 年，研究區大部分區域干旱等級由較干旱向正常、較濕潤等級轉化，只有西南部小范圍由干旱變為其他等級，東部由較濕潤、正常、較干旱轉變為干旱。2010—2018 年，研究區中部、東部干旱等級由干旱轉為其他等級，旱情程度下降，中部及北方大部分干旱等級則由較干旱變為正常，西南小部分區域由較干旱轉換為干旱，這一時段研究區旱情轉移基本表現為由干旱變為較干旱、較干旱變為正常，研究區大部分旱情得到一定程度的緩解。

圖8 神東礦區1991—2018 年植被生長季TVDI 轉移矩陣Fig.8 TVDI transfer matrix of vegetation growing season from 1991 to 2018 in Shendong Mining Area

綜上所述，神東礦區在研究區間內，干旱等級變化均有不同程度的波動，但在空間上基本表現為由干旱等級向較干旱、正常等級轉變，較干旱等級向正常、較濕潤等級進行轉變，說明神東礦區進行的一系列生態恢復及水土保持措施達到了一定成效。神東礦區地處我國西北內陸，空間跨度較大，因此水熱狀況差異較大，且隨著地表與地下的采煤擾動，地表覆蓋變化復雜，加之自然因素及人類活動的影響導致了地表水分在時間上和空間上變化明顯[26]。

3 討 論

本文基于TVDI模型對神東礦區1991—2018年植被生長季的旱情進行監測。從TVDI的月際變化來看，神東礦區每年植被生長季內6—9 月TVDI出現最低值，該結果與神東礦區降水集中的特點相吻合[27]。TVDI的年際變化顯示，神東礦區自1991—2018 年，TVDI整體呈下降趨勢，下降速率為0.02/10 a，地表水分逐年上升，旱情逐步得到緩解，該結果與劉英等[28]基于梯度結構相似度的神東礦區土壤濕度空間分析得到的神東礦區土壤濕度在空間上呈增加趨勢，旱情緩解的結果一致。對TVDI時間變化進行趨勢分析發現，神東礦區2018 年TVDI下降趨勢最為顯著，這是由于神東礦區在2018 年前后實施了大量的水土保持措施，在不考慮自然因素對各年TVDI變化影響的情況下，研究結果與實地調查結果一致。神東礦區1991—2018 年干旱等級的空間分布表明，研究區干旱程度，始終為西南部高于東北部，這是由于該地區西南部地處毛烏素沙地邊緣[29]。干旱等級面積減少，正常及較濕潤面積增加，該結果與TVDI的時間變化特征相吻合。從線性傾向趨勢的空間分布來看，研究區整體TVDI下降趨勢面積遠遠大于TVDI上升趨勢面積，且TVDI呈上升趨勢的區域集中在神東礦區的礦井群附近，這是由于采煤后地表塌陷，地表發生變形，地表水下滲速度加快，加快該區域的水土流失[30]。利用轉移矩陣分析發現，整個研究區在1991—2018年間，干旱等級向較干旱、正常轉化，而較干旱等級大多向正常、較濕潤轉化，旱情逐步得到緩解。

本文是在保證TVDI模型對神東礦區旱情監測具有較高的精度下進行研究，研究發現的TVDI在時間序列中存在多個突變點問題，還需要在后續研究中對突變點的變化情況進行相關驗證。下一步還需要對造成其發生時空演變的自然及人為因素進行更深一步的研究，以期為神東礦區的水土保持措施及生態恢復提供參考。

4 結 論

1）1991—2018 年神東礦區植被生長季TVDI在0.35～0.85 之間波動，每年6—9 月出現低值。TVDI月際變化整體呈下降趨勢。

2）年際變化表明，1991、2002、2010 年TVDI偏離值＜0，TVDI呈下降趨勢，2014、2018 年TVDI偏離值＞0，TVDI呈上升趨勢。TVDI在時間序列上整體呈下降趨勢，下降速率為0.02/10 a。

3）神東礦區1991—2018 年植被生長季TVDI干旱等級空間分布，始終表現為西南高于東北。研究區干旱等級面積逐年減少，正常、較濕潤等級面積逐年增加。

4）研究區TVDI下降趨勢面積遠大于TVDI上升趨勢面積。顯著下降趨勢面積最大為2 394.74 km2，顯著上升面積最小為425.91 km2。

5）1991—2018 年，TVDI的干旱等級逐漸向較干旱、正常等級轉變。TVDI的較干旱等級逐漸向正常和較濕潤等級轉變。向干旱及較干旱等級轉變的范圍較小。

（作者聲明本文無實際或潛在的利益沖突）