基于MODIS數據的山東省2014-2016年干旱監測分析

王正東, 郭 鵬, 萬 紅, 楊 綱

(山東農業大學 信息科學與工程學院, 山東 泰安 271018)

干旱是我國乃至世界上許多國家主要的自然災害之一，不僅會對生態系統和環境造成極大破壞，還會嚴重影響社會經濟活動以及居民生活。在全球變暖的大趨勢下，干旱的發生變得越來越頻繁，干旱監測也成為全球關注的一個科學熱點問題。山東省地處黃淮海地區，是主要的產糧大省，受季風影響明顯，屬于水資源短缺省份，近幾年，出現了持續的少雨及異常氣候事件頻率增多的現象，受旱災影響嚴重。因此，如何及時、有效、準確地監測黃淮海地區干旱的發生、旱情的發展已成為有關部門關注的焦點。

常規的土壤水分監測手段通過人工測墑完成，耗時、耗力。利用遙感衛星監測干旱則具有周期短、范圍廣等優點，通過對研究區的遙感影像進行分析，可以快速獲得干旱數據，確定干旱等級，從而為防災減災工作提供依據[1]。溫度植被干旱指數(Temperature Vegetation Dryness Index，TVDI)是遙感監測干旱的常用方法之一，利用該方法可以及時、準確和有效地監測干旱情況，能較好地監測不同地表面，且具有計算較簡單，容易實現等優點[2]，許多學者已進行過相關研究。如齊述華等[3]利用NOAA-AVHRR數據，采用溫度植被干旱指數法反映表層土壤水分變化情況。王婷婷等[4]利用趨勢線法對2002—2009年每年8月份的TVDI值進行回歸分析，來研究2002—2009年松遼平原的干旱變化趨勢。

本文以山東省為研究區，利用MODIS的NDVI，EVI，LST產品，在使用Savitzky-Golay(S-G)濾波方法進行重構，填補缺失數據的基礎上，分別構建NDVI-LST和EVI-LST特征空間，建立基于時間序列的溫度植被干旱指數對山東省2014—2016年干旱變化的時空演變特征進行分析，并探討其與氣象因子之間的關系。

1 數據與方法

1.1 研究區概況

山東省位于中國東部沿海，地處114°20′—122°43′E， 34°23′—38°23′N，地幅遼闊，陸地南北最長距離約420 km，東西最寬約700 km，全省總面積為15.8萬km2，境內最高處位于泰山，海拔為1 524 m，最低處位于東北部的黃河三角洲，海拔為2 m左右。山東的氣候類型屬暖溫帶季風氣候，降水集中，四季分明，冬夏較長，全年平均溫度為10～14℃，年平均降水量為550～950 mm，但是降水季節分布不均衡，對農業生產影響較大，因此尋求一種大范圍、實時的旱情監測手段對于山東省的社會經濟發展具有現實推動意義。

1.2 數據來源與預處理

本研究所使用的MODIS遙感數據從MODIS官網(https:∥ladsweb.modaps.eosdis.nasa.gov/)下載得到，主要有8 d合成的地表溫度產品(MOD11A2)和16 d合成的植被指數產品(MOD13A2)，分辨率為1 km，時間范圍為2014年1月—2016年12月，MOD11A2和MOD13A2數據的原始格式都為EOS-HDF。在對數據預處理時，首先利用MRT(MODIS Reprojection Tools)對兩種產品數據進行格式轉換、拼接、重投影等批處理，然后根據山東省的行政區劃邊界裁剪得到研究區影像，最后得到研究區3年的歸一化植被指數產品(NDVI)數據、增強型植被指數產品(16 d合成數據)和陸地表面溫度數據(8 d合成數據)以及相應的質量控制文件。

溫度和降水等氣象數據來源于中國氣象科學數據共享服務網(http:∥cdc.nmic.cn/)提供的中國地名氣候資料月值數據集，本文選取惠民、定陶、濟南、兗州、濰坊5個國家級氣象站的每月月平均降水量和月平均溫度數據，起止時間為2014—2016年，研究區氣象站點的空間分布如圖1所示。

圖1 研究區范圍與氣象站點位置

1.3 數據重建

在獲取遙感影像時，因為存在傳感器設計、大氣傳輸、太陽光照角度、觀測視角等隨機干擾因素，影像中會存在很多噪聲，當有云覆蓋時，會造成數據的缺失，所以獲取的影像數據有時無法準確表達地表情況，因此有必要進行數據重建。由Savitzky等[5]提出的移動窗口的加權平均算法S-G濾波重建能夠較好地完成數據平滑，內插，反映真實的植被生長、氣象變化等地表情況，適用于NDVI，EVI等數據的重建[6]。S-G濾波是一種特殊的低通濾波器，用于數據的噪聲平滑，具有實現簡單、先驗知識少等優點，可以較真實地體現植被的生長變化情況[7]。

S-G濾波是以平滑時間序列數據為移動窗口的加權平均算法[8]，該方法中的加權系數取決于一個濾波窗口內給定高階多項式的最小二乘擬合次數。若某個數據集內包含n個點，半窗口的寬度為m，濾波器的長度等于滑動數組的寬度等于N。Ci為平滑窗口內對應每個點的卷積系數，則對于第j個點來說，平滑后的值即為：

(1)

式中：Yj*為擬合后的結果。

根據像元可信度值(pixel reliability)和權重選擇質量控制(quality control)文件分別對影像的每個像元賦予權重，通過編程分別對NDVI，EVI和LST進行S-G時序濾波處理。除因時間而異的觀測條件會使時間序列呈現類似于鋸齒狀的不規則波動，在其中的某一次觀測時由于不同設備間的電磁干擾，或者觀測設備的位置偏移等原因會產生高頻噪聲與干擾，所以還需要進行空間濾波，目的是提高單次觀測獲取的影像的質量，因此本研究使用中值濾波消除空間上的異常值。

1.4 數據重建效果分析

圖2 S-G濾波重構值、站點實測值與原始值的對比

1.5 植被指數-地表溫度特征空間

Lambin等[9]指出了地表溫度與歸一化植被指數兩者存在的關系，建立TVDI模型對于消除土壤背景的影響效果較好。如果研究區是從裸土到植被完全覆蓋，土壤濕度從極濕潤變化到極干旱時，以LST和VI為縱橫坐標構建的散點圖呈三角形。植被覆蓋區的TVDI理論取值應處于(0，1)之間[10]。TVDI值越低則表明干旱程度越輕，反之TVDI值越高，則干旱程度越嚴重。TVDI由植被指數(NDVI/EVI)和地表溫度(LST)計算得到，其定義為公式:

(2)

式中：TS為在一定的分辨率條件下任意像元的地表溫度；TSmin表示某一NDVI對應的最低地表溫度值，即濕邊；TSmax為某一NDVI對應的最高地表溫度值，即干邊。

TSmax=a1+b1×NDVI

TSmin=a2+b2×NDVI

(3)

式中：a1，b1是干邊擬合方程的系數；a2，b2是濕邊擬合方程的系數[10]。由地表溫度最大值擬合而成的稱為干邊，最小值擬合而成的則稱為濕邊。由TVDI的原理可知，地表溫度最大值與植被指數呈負線性關系。若植被覆蓋度小于20%，低植被覆蓋度對應的NDVI不能準確反映出區域內植被生長情況；若植被覆蓋度大于80%時，其對應的NDVI增長速度會逐漸減慢并趨于平緩，出現一種飽和狀態，并且對于植被監測的靈敏度也下降。由此得出，NDVI更適用于中等植被覆蓋率的情況。因此，在擬合特征空間中每一步長NDVI對應LST最大值與最小值時，選擇處于0.2～0.8范圍內的VI值[11]。

2 結果與分析

2.1 TVDI計算

在VI為0.2～0.8的范圍內，以0.01為步長，提取每一步長NDVI，EVI像元對應的所有LST像元中的最大值與最小值，并采用最小二乘擬合方法對每期特征空間的干濕邊進行線性擬合，得到干濕邊的擬合方程及其相關系數，結果見表1—2。從干濕邊的擬合結果來看，NDVI與LST的最大值呈負相關關系，與最小值呈正相關關系。EVI-LST構成的特征空間干濕邊斜率與NDVI-LST一致，但EVI與LST構建的特征空間中干濕邊擬合的相關系數更高，主要因為EVI加入了藍色波段來增強植被信號，消除土壤背景和氣溶膠散射的影響，因此能夠克服NDVI在高植被覆蓋地區增長緩慢趨于平緩，容易飽和，在低植被覆蓋地區容易受到土壤植被影響等缺點。進入7月份后，相關系數均降低，原因可能是進入夏季后，山東大部分地區地表溫度普遍較高，因此隨著NDVI，EVI的變化，地表溫度變化較小。

圖3為2015年1月濟南市NDVI，EVI與LST構建的散點圖，由圖3可知，相比NDVI，EVI總體偏低，適用于高植被覆蓋區域，在構建特征空間時，擬合的干濕邊更穩定，趨勢更明顯，干濕邊的相關系數分別達到0.84，0.86，而NDVI的干濕邊相關性僅為0.79，0.74，因此，本文選擇EVI和LST構建TVDI模型進行下一步的研究。

表1 NDVI-LST干濕邊擬合方程及其相關系數

表2 EVI-LST干濕邊擬合方程及其相關系數

2.2 研究區干旱狀況時空分析

根據EVI-LST計算得到山東省2014—2016年每個月TVDI圖像，再以3—5月份為春季，6—8月份為夏季，9—11月份為秋季，12月—翌年2月份為冬季，制作每個季度平均TVDI圖像。然后按照整個研究區域內像元TVDI值的分布直方圖和中國土壤濕度界定干旱的標準對研究區干旱情況進行分級[12]，將干旱劃分為5個等級，分別是潮濕(0≤TVDI≤0.2)、濕潤(0.2≤TVDI≤0.4)、正常(0.4≤TVDI≤0.6)、干旱(0.6≤TVDI≤0.8)、重旱(0.8≤TVDI≤1.0)，分級結果如圖4所示。

為了分析山東省2014—2016年土壤水分時空分布變化趨勢，統計山東省每年各干旱等級土壤面積以及百分比(表3)。從分析結果看，2014年上半年旱情較嚴重，進入秋季之后，重旱與干旱區域的面積顯著下降，表明干旱得到一定程度的緩解；2015年旱情為3年中最嚴重的一年，重旱與干旱區域的面積比例分別高達7.46%，34.24%，表明在一年中山東省有近一半面積處于干旱狀態，全年各干旱等級土壤面積比例與2014年基本相一致，但是2015年是在進入冬季之后，干旱才得到減弱，東北部地區表現為相對濕潤，其他地區基本屬于正常狀態，只有局部地區存在干旱情況；2016年正常和濕潤區域面積占比分別達到44.52%，42.54%，相對于前兩年，干旱面積顯著下降。

圖3NDVI，EVI-LST特征空間干濕邊方程

圖4研究區3年旱情等級時空變化

表3 不同干旱等級面積及比例

2.3 TVDI與氣象數據的相關性分析

濕潤指數中涉及多種氣象因子，為了進一步探究氣象因子與TVDI變化的關系，本文研究了主要氣象因子氣溫和降水的年際變化及其與TVDI的關系。利用山東省的濟南、兗州、濰坊、定陶和惠民5個國家級氣象站收集的每月平均溫度和月均降水量與TVDI進行相關性分析，根據溫度、降水量與TVDI的相關系數(表4)可以得出，TVDI與溫度均通過0.01水平的顯著性檢驗，表明相比于降水，TVDI對溫度更為敏感，在溫度升高時，作物生長需水量大，地表水分蒸發速度加快，蒸發量大于降雨量，土壤含水量降低，使得地表較為干旱，導致TVDI值升高；TVDI與降水也均通過0.05水平的顯著性檢驗，且TVDI與降水量存在顯著負相關關系，當存在連續降水時，TVDI的值會呈降低趨勢，但兩者并不同步，TVDI的變化存在一定的滯后性，當降水量變化穩定時，TVDI也隨之趨于穩定。以濰坊市為例，2014—2016年3年中的降水量明顯少于其他區域，從氣象站統計3年年均溫度與年均降水量中可以看出3年年平均降水量僅為378，435，478 mm，相比于同時期其他地區，明顯偏少，所以濰坊在3年中一直處于TVDI偏高的情況。

2015年降水量相比于2014年有所增加，但是受到“歷史最強”厄爾尼諾的影響，干旱情況更加嚴重，2015年TVDI與溫度的相關系數較高，除定陶外均達到0.76以上，而同年降水與TVDI的相關系數僅有0.5左右，最高僅達到0.69，因此溫度成為2015年TVDI值升高、干旱加重的主要影響因素，所以即使在降水量增加的情況下，干旱程度并沒有減弱反而增強。2016年年平均降水量明顯升高，濰坊、定陶、兗州等地降水量均達到3年最高水平，TVDI反演的季度干旱情況也表明全省干旱情況有所減弱。

表4 溫度、降水量與TVDI的相關系數

注：*表示達到0.01的顯著性水平；**表示達到0.05的顯著性水平；剩余表示未通過顯著性檢驗。

3 結論與討論

(1) 利用S-G濾波與空間濾波對植被指數與地表溫度進行重建，可以有效消除衛星數據在時間和空間上存在的缺失值和噪聲的影響，對濾波后的溫度值與實測值做相關性分析，相關系數均在0.91以上，說明經過重構后的數據與站點實測數據在趨勢上基本一致，因此可以利用重構后的數據進行TVDI模型的構建。

(2) NDVI/EVI-LST特征空間的干邊隨著植被指數的增大而呈遞減趨勢，而濕邊在氣溫較低月份與植被指數關系明顯，另外NDVI與EVI在構建特征空間上也有所不同，NDVI與LST構成的特征空間更平緩；而EVI與LST構成的特征空間范圍更集中，在干濕邊擬合時，趨勢更穩定，相關性更高，因此運用EVI-LST構建TVDI更為合理。

(3) 從時間上看，山東省2014年干旱主要集中在春季和夏季，全年平均干旱面積占到全省面積的37.62%，在進入春季之后，溫度迅速回升，植被生長迅速，蒸發旺盛，干旱面積顯著提高，進入8月之后，除局部發生重旱外，大部分地區的旱情得到緩解。2015年干旱分布情況基本與2014年一致，但持續時間更長，全年平均干旱面積比例增加，占全省面積的41.7%。只有當進入12月份之后，干旱才減弱。2016年由于降水量顯著增多，只有在2,3月份有一次范圍較大的旱情，其余月份除局部地區外全省大部分地區相對濕潤。

(4) 從空間上看，2014年和2015年年均重旱面積比例分別達到4.93%，7.46%，其中魯中與半島交界處、魯西南地區是干旱較為嚴重的地區，根據3年的變化情況分析得出，干旱分布體現為由東南向西北轉移的趨勢。濕潤以及正常的區域主要分布在中部的山區、最北部的平原地區以及膠東地區。

(5) 對TVDI與降水量、溫度數據做相關性分析，均達到顯著相關，并發現TVDI與溫度的相關性更強，且通過0.01水平的顯著性檢驗，表明改進的TVDI指數能夠較好地反映山東省2014—2016年的旱情變化情況，對研究區旱情的快速準確檢測和干旱演變過程的研究具有一定的參考價值，有助于制定相應的防災和減災決策。

受限于MODIS數據產品本身的問題，本文仍存在一些不足之處：溫度是TVDI的影響因子之一，而高程的變化是影響地表溫度的重要條件。山東省地處丘陵地帶，境內高程變化較大，因此還應該考慮高程對溫度變化的影響，利用高程對地表溫度進行校正從而計算TVDI。另外TVDI與溫度、降水數據與有些氣象站點的相關關系較弱，其原因可能是：① 由于范圍不同，整個行政區觀測數據與3 km×3 km緩沖區數據的相關性不可能達到很高的水平；② 雖然氣象資料是準確的，但參與TVDI建模的地表溫度和植被指數都是按固定期數合成的數據，盡管假設兩者在時相上是相一致的，但實際上可能并不完全一致，存在部分偏差；③ TVDI與氣象因子對干旱的反映側重點不同，TVDI是地面植被、土壤等受干旱影響的響應程度，氣象因子則比較直觀地表現地面濕潤情況。