焦云峰, 郭 鵬,高仁鳳，孔令偉

(山東農(nóng)業(yè)大學(xué) 信息科學(xué)與工程學(xué)院，山東 泰安 271018)

干旱自古以來是我國(guó)頻發(fā)的自然災(zāi)害，干旱使水源匱乏，危害作物生長(zhǎng)，造成作物減產(chǎn)，嚴(yán)重時(shí)還會(huì)危害居民生活，影響工業(yè)生產(chǎn)及其他社會(huì)經(jīng)濟(jì)活動(dòng)。山東省是黃淮海地區(qū)的主要產(chǎn)糧大省，受季風(fēng)影響極大，屬于水資源短缺省份[1]，同時(shí)也是受旱災(zāi)影響比較嚴(yán)重的省份之一，有“十年九旱”之說。因此找到一種實(shí)時(shí)、有效、精準(zhǔn)的監(jiān)測(cè)方法來對(duì)其進(jìn)行監(jiān)測(cè)，對(duì)山東省政府進(jìn)行抗旱救援工作有很大的幫助意義。

目前國(guó)內(nèi)外干旱監(jiān)測(cè)方法較多，主要有基于反射率的植被指數(shù)法[2]、基于地表輻射的植被指數(shù)法[3]、熱慣量法[4]、微波遙感法[5]。其中Sandholt[6]等提出的溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index, TVDI)方法被廣泛應(yīng)用，如齊述華[7]等人利用溫度植被干旱指數(shù)(TVDI)進(jìn)行全國(guó)旱情監(jiān)測(cè)研究；閆娜[8]等通過構(gòu)建TVDI特征空間進(jìn)行陜西省干旱監(jiān)測(cè)；Zhang[9]等人基于近紅外和紅色波段數(shù)據(jù)，構(gòu)建了比率干燥監(jiān)測(cè)指數(shù)(RDMI)，發(fā)現(xiàn)該指數(shù)可用于不同尺度的干燥監(jiān)測(cè)；Pandey[10]等人通過地表溫度、蒸散蒸騰量和歸一化植被指數(shù)對(duì)印度本德爾坎德邦進(jìn)行干旱災(zāi)害評(píng)估，研究結(jié)果表明可以在更大尺度上進(jìn)行干旱評(píng)估；Liu[11]等結(jié)合光學(xué)遙感和微波遙感建立微波歸一化植被指數(shù)，經(jīng)實(shí)驗(yàn)發(fā)現(xiàn)該指數(shù)能在不同程度上區(qū)分干旱水平；Park[12]等人研究提出了高分辨率土壤水分干旱指數(shù)，并驗(yàn)證了該指數(shù)特別是在旱季能夠很好地描述氣象干旱；Yang[13]等人利用干旱嚴(yán)重程度指數(shù)、聚類算法和嚴(yán)重-區(qū)域-持續(xù)時(shí)間方法識(shí)別干旱事件，分析中國(guó)各種干旱特征的時(shí)空分布，并對(duì)中國(guó)未來旱情做了預(yù)測(cè)；高雅[14]等人通過構(gòu)建NDVI-LST特征空間計(jì)算得到TVDI，由此對(duì)曲靖市進(jìn)行干旱監(jiān)測(cè)。隨著該方法的廣泛使用，利用溫度植被干旱指數(shù)(TVDI)進(jìn)行干旱監(jiān)測(cè)已成為遙感干旱監(jiān)測(cè)常用方法之一。本文以山東省為研究區(qū)域，利用MODIS的LST(Land Surface Temperature)、NDVI(Normalized Difference Vegetation Index)、EVI(Enhanced Vegetation Index)產(chǎn)品，分別構(gòu)建NDVI-LST、EVI-LST特征空間，計(jì)算TVDI值，基于TVDI值得到山東省干旱監(jiān)測(cè)專題圖，分析山東省干旱情況。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

山東省地處中國(guó)華東沿海地區(qū)，位于北緯34°22.9′～38°24.01′、東經(jīng)114°47.5′～122°42.3′之間，東西長(zhǎng)721.03 km，南北長(zhǎng)437.28 km，全省陸域面積15.58萬km2。山東省內(nèi)最高點(diǎn)位于泰山，海拔高度為1 532.7 m；最低處位于黃河三角洲，一般海拔2～10 m。氣候類型為暖溫帶季風(fēng)氣候，降水集中，雨熱同季，春秋短暫，冬夏較長(zhǎng)，年平均氣溫11～14 ℃，年平均降水量一般在550～950 mm之間，但是降水季節(jié)分布不均衡，極易形成旱澇災(zāi)害，對(duì)農(nóng)業(yè)生產(chǎn)甚至其他社會(huì)經(jīng)濟(jì)活動(dòng)影響較大。因此對(duì)山東省進(jìn)行實(shí)時(shí)的大范圍干旱監(jiān)測(cè)對(duì)社會(huì)經(jīng)濟(jì)發(fā)展有極大的輔助作用。山東省海拔與氣象站點(diǎn)分布如圖1所示。

圖1 山東省海拔與氣象站點(diǎn)分布圖

1.2 數(shù)據(jù)來源與預(yù)處理

本文使用的MODIS數(shù)據(jù)從NASA網(wǎng)站(https://search.earthdata.nasa.gov/)下載得到，時(shí)間為2018年1月～2018年12月，分辨率為1 km的8天合成的地表溫度產(chǎn)品(MOD11A2)和16天合成的植被指數(shù)產(chǎn)品(MOD13A2)，下載得到的原始數(shù)據(jù)格式均為EOS-HDF。首先利用MRT(MODIS Reprojection Tools)工具對(duì)MODIS原始數(shù)據(jù)進(jìn)行預(yù)處理，包括對(duì)原始影像的鑲嵌，數(shù)據(jù)格式的轉(zhuǎn)換，影像重采樣、重投影、裁剪等批處理。得到研究區(qū)格式為TIFF的23幅16天合成的NDVI/EVI植被指數(shù)產(chǎn)品數(shù)據(jù)和46幅8天合成的陸地表面溫度(LST)產(chǎn)品數(shù)據(jù)，隨后通過最大值合成法(MVC)對(duì)46幅地表溫度產(chǎn)品進(jìn)行處理，最終得到與23幅16天合成的NDVI/EVI植被指數(shù)產(chǎn)品相對(duì)應(yīng)的23幅16天合成的地表溫度(LST)產(chǎn)品，再利用ArcGIS軟件和ENVI軟件進(jìn)行裁剪，得到最終的研究區(qū)域數(shù)據(jù)。

本文使用的三個(gè)國(guó)家級(jí)氣象站點(diǎn)數(shù)據(jù)來源于溫室數(shù)據(jù)共享平臺(tái)(http://data.sheshiyuanyi.com/)。三個(gè)國(guó)家級(jí)氣象站分別為濟(jì)南站、濰坊站、定陶站，數(shù)據(jù)分別是2018年日平均溫度和日平均降水量。數(shù)據(jù)起止時(shí)間為2018年1月1日～2018年12月31日。對(duì)三個(gè)氣象站獲取的日平均溫度和降水量分別按月計(jì)算平均值，得到每月的平均溫度和平均降水量。

1.3 植被指數(shù)-地表溫度特征空間

Price[15]等通過研究發(fā)現(xiàn)在植被覆蓋區(qū)NDVI與LST形成的散點(diǎn)圖大多呈三角形或梯形分布，Sandholt[6]根據(jù)NDVI與LST相對(duì)應(yīng)而構(gòu)建的特征空間，提出了溫度植被干旱指數(shù)(TVDI)，公式如下：

(1)

式中，Ts為在某一分辨率條件下任意像元的地表溫度；Tsmin為某一NDVI對(duì)應(yīng)的最低地表溫度值，稱為濕邊，由地表溫度最小值擬合而成；Tsmax為某一NDVI對(duì)應(yīng)的最高地表溫度值，稱為干邊，由地表溫度最大值擬合而成。Tsmax和Tsmin的計(jì)算公式如下：

Tsmax=a1+b1×NDVI

(2)

Tsmin=a2+b2×NDVI

(3)

式中，a1、b1是干邊擬合方程的系數(shù)；a2、b2是濕邊擬合方程的系數(shù)。

在計(jì)算TVDI時(shí)，NDVI與LST一一對(duì)應(yīng)。由TVDI定義可知，TVDI值位于0～1之間，當(dāng)?shù)乇頊囟鹊扔诘乇碜畹蜏囟葧r(shí)TVDI值為0，地表溫度等于地表最高溫度時(shí)TVDI值為1。因此，若某地區(qū)TVDI值越大、距干邊越近，則該地區(qū)干旱程度越嚴(yán)重，反之TVDI值越小、距離濕邊越近，則干旱程度越輕[16]。

2 結(jié)果與分析

2.1 TVDI計(jì)算

研究表明，當(dāng)利用NDVI與LST計(jì)算TVDI，在NDVI小于0.2時(shí)，意味著該區(qū)域植被覆蓋率小于20%，不能有效地反映該區(qū)域的植被情況；在NDVI大于0.8時(shí)，意味著植被覆蓋率大于80%，此時(shí)植被趨于飽和，NDVI增加不明顯，對(duì)于植被監(jiān)測(cè)有一定的影響。由此確定NDVI的有效范圍是[0.2,0.8][17]。利用ENVI軟件中的擴(kuò)展工具進(jìn)行TVDI的計(jì)算,首先對(duì)NDVI/EVI和LST值分別進(jìn)行波段運(yùn)算，前者分別乘以0.000 1得到實(shí)際的植被指數(shù)影像，后者乘以0.02得到真實(shí)地表溫度；隨后分別以NDVI/EVI作為x軸，LST作為y軸繪制NDVI/EVI-LST的二維散點(diǎn)圖，提取每一步長(zhǎng)的NDVI/EVI相對(duì)應(yīng)LST的最大值與最小值，再將提取出的LST的最大值和最小值與對(duì)應(yīng)的NDVI/EVI建立特征空間得到對(duì)應(yīng)的干濕邊，并對(duì)其進(jìn)行干濕邊的擬合，得到干濕邊的擬合曲線以及干濕邊方程，之后根據(jù)干濕邊方程計(jì)算TVDI值，得到研究區(qū)域的TVDI圖像。

比較NDVI-LST與EVI-LST二者所擬合的干濕邊方程，可發(fā)現(xiàn)在全年大部分時(shí)間內(nèi)EVI-LST擬合干濕邊方程的相關(guān)系數(shù)要大于NDVI-LST擬合干濕邊方程的相關(guān)系數(shù)，雖然有部分日期NDVI-LST擬合方程的相關(guān)系數(shù)要大于EVI-LST所擬合方程，其原因可能是由于EVI-LST干、濕邊擬合線的斜率較小，最高溫度和最低溫度變化區(qū)間小，雖然總體呈現(xiàn)上升或下降的趨勢(shì)，局部的變化趨勢(shì)有反復(fù)，數(shù)據(jù)點(diǎn)在較小的溫度變化區(qū)間內(nèi)分散分布所致[18]。因此選用EVI和LST來進(jìn)行TVDI的計(jì)算。

2.2 研究區(qū)干旱時(shí)空分析

根據(jù)前面的分析，利用EVI-LST計(jì)算出TVDI后，通過ENVI軟件對(duì)TVDI圖像進(jìn)行干旱等級(jí)的劃分，國(guó)家標(biāo)準(zhǔn)中主要以降水量、土壤濕度等5種單項(xiàng)指標(biāo)并綜合其他指標(biāo)將旱情分為5個(gè)等級(jí)：正常或濕澇、輕旱、中旱、重旱和特旱[19]，同時(shí)依據(jù)山東省氣候?qū)嶋H情況，將旱情劃分為五個(gè)等級(jí)，分別是：潮濕(0

圖2 山東省2018年TVDI干旱等級(jí)分布圖

根據(jù)旱情等級(jí)分布圖，2018年山東省的旱情在第65～81天這個(gè)時(shí)間段內(nèi)旱情較為嚴(yán)重；在第145天前后開始出現(xiàn)輕微旱情，到第177天前后旱情加重，隨后旱情逐漸減弱；到第273天前后，出現(xiàn)較為嚴(yán)重的旱情，最后到2018年底都未再次出現(xiàn)大規(guī)模的干旱情況。根據(jù)TVDI在一年內(nèi)對(duì)應(yīng)天數(shù)轉(zhuǎn)換到月份，將以第129天為第一天的16天合成的TVDI值獲得的干旱分布圖定為5月，其余圖像每?jī)煞鶜w為一個(gè)月，從月尺度上來看，山東省2018年干旱主要分布在3月、5～7月、9～10月。2018年3月以來山東省進(jìn)入春旱，至4月份旱情有所好轉(zhuǎn)，這與山東省氣象局資料符合；到6月份進(jìn)入夏季，氣溫迅速升高、植被生長(zhǎng)迅速，出現(xiàn)干旱現(xiàn)象，但隨后到8～9月干旱面積迅速減少，同期氣象資料顯示該時(shí)期受臺(tái)風(fēng)和氣壓影響，降水增多，使得旱情大幅緩解；到2018年冬季全省氣候基本穩(wěn)定，沒有較大旱情發(fā)生。從空間來看，山東省旱情主要發(fā)生在魯南地區(qū)和內(nèi)陸地區(qū)，與山東省氣象資料基本相同。

2.3 TVDI與地面氣象站實(shí)測(cè)數(shù)據(jù)的相關(guān)性分析

由于衛(wèi)星數(shù)據(jù)受天氣、云層等影響過大，本文還通過對(duì)濟(jì)南、濰坊、定陶三個(gè)國(guó)家級(jí)地面氣象站收集到的每日的平均溫度和平均降水量按月合成，分別與TVDI值進(jìn)行相關(guān)性分析，結(jié)果如表1所示。通過分析表1的結(jié)果可發(fā)現(xiàn)，三個(gè)氣象站的月平均溫度與TVDI值的相關(guān)系數(shù)均在0.6以上，平均降水與TVDI值的相關(guān)系數(shù)在0.31以上，這說明通過對(duì)從衛(wèi)星獲取的地表溫度與植被指數(shù)產(chǎn)品進(jìn)行反演，進(jìn)而獲取的TVDI指數(shù)可以作為判斷是否干旱的指標(biāo)。同時(shí)還可以通過對(duì)比溫度和降水量二者分別與TVDI的相關(guān)系數(shù)，發(fā)現(xiàn)溫度與TVDI的相關(guān)系數(shù)大于降水量與TVDI的相關(guān)系數(shù)，由此可發(fā)現(xiàn)TVDI受溫度影響多大于受降水影響。因?yàn)門VDI值受溫度與降水的影響不同，所以在之后的研究中可以對(duì)溫度和降水按照不同的權(quán)重進(jìn)行相關(guān)性分析，使通過TVDI值進(jìn)行的干旱劃分監(jiān)測(cè)更加準(zhǔn)確可靠，相關(guān)性更高。此外，由于氣象站點(diǎn)收集的數(shù)據(jù)代表的范圍較小，即使使用ArcGIS軟件對(duì)氣象站點(diǎn)進(jìn)行緩沖區(qū)分析，擴(kuò)大站點(diǎn)的覆蓋范圍，仍不能較好地代表大面積地區(qū)的實(shí)際情況，從而使得部分站點(diǎn)的氣溫與降水量和TVDI的相關(guān)性系數(shù)較低。因此在以后進(jìn)行相關(guān)性分析時(shí)，可增加氣象站點(diǎn)的個(gè)數(shù)，使之盡量均勻地分布在研究區(qū)域內(nèi)，以此來更好地對(duì)衛(wèi)星數(shù)據(jù)反演進(jìn)行驗(yàn)證。

表1 溫度和降水量與TVDI相關(guān)系數(shù)表

3 結(jié)論與討論

(1)本研究利用MODIS植被指數(shù)產(chǎn)品(MOD13A2)和地表溫度產(chǎn)品(MOD11A2)分別構(gòu)建了NDVI-LST特征空間和EVI-LST特征空間。通過對(duì)比兩種特征空間，發(fā)現(xiàn)NDVI-LST特征空間更為平緩，像元點(diǎn)范圍更大；而EVI-LST特征空間則更為緊湊，像元更為集中。因此在進(jìn)行干濕邊曲線擬合時(shí)，EVI-LST特征空間干濕邊擬合程度較好，相關(guān)性更高。

(2)從時(shí)間上來看，山東省2018年旱情主要出現(xiàn)在春季，夏季和秋季也曾出現(xiàn)過短暫的旱情；從空間上來看，山東省2018年的干旱地區(qū)主要分布在魯南、魯西等內(nèi)陸地區(qū)和魯中與半島交界處，呈現(xiàn)出一種圍繞山東省內(nèi)陸邊界依次干旱的現(xiàn)象。濕潤(rùn)及正常范圍主要分布在魯中以及膠東半島等地區(qū)。

(3)通過三個(gè)國(guó)家級(jí)氣象站的溫度和降水與TVDI之間的相關(guān)系數(shù)可以發(fā)現(xiàn)，TVDI受溫度影響明顯大于受降水影響。此外降水對(duì)于TVDI值的影響與降水量有關(guān)，降水越多TVDI值越大。

由于本研究所需的數(shù)據(jù)產(chǎn)品為地表溫度產(chǎn)品和植被指數(shù)產(chǎn)品，受外界因素如云層等影響較大，容易產(chǎn)生噪聲使LST、NDVI、EVI產(chǎn)生異常值，在TVDI值計(jì)算時(shí)會(huì)產(chǎn)生相應(yīng)的影響，因此對(duì)衛(wèi)星數(shù)據(jù)進(jìn)行了Hants濾波，使其光譜曲線變平滑，但是效果提升有限。在查閱部分文獻(xiàn)后發(fā)現(xiàn)S-G濾波相較于Hants濾波處理效果更優(yōu)，因此在以后進(jìn)行相關(guān)研究時(shí)，可利用相應(yīng)的濾波技術(shù)如S-G濾波對(duì)地表溫度和植被指數(shù)進(jìn)行重建，以此來減小噪聲去除異常值。