基于MODIS_TVDI/GNSS_PWV的云南省干旱特征時(shí)空分析

于 維，柯福陽(yáng)，曹云昌

(1.南京信息工程大學(xué)遙感與測(cè)繪工程學(xué)院，南京 210044；2.中國(guó)氣象局氣象探測(cè)中心，北京 100081)

0 引言

我國(guó)是受旱災(zāi)影響最大的國(guó)家之一，近幾十年來(lái)，由于氣溫逐漸升高、降水逐漸減少導(dǎo)致旱災(zāi)頻發(fā)[1]，對(duì)農(nóng)業(yè)經(jīng)濟(jì)造成了巨大損失。傳統(tǒng)的干旱監(jiān)測(cè)多以氣象監(jiān)測(cè)為主，精度高，但效率低、耗費(fèi)人力、應(yīng)用范圍小且受站點(diǎn)分布影響較大[2]，遙感技術(shù)的飛速發(fā)展解決了這一難題。MODIS作為中分辨率遙感衛(wèi)星，是目前干旱監(jiān)測(cè)主要的數(shù)據(jù)來(lái)源之一。溫度植被干旱指數(shù)(temperature vegetation dryness index,TVDI)是遙感干旱監(jiān)測(cè)中常用的方法，它反映了植被覆蓋度與地表溫度變化之間的關(guān)系，進(jìn)而反映土壤濕潤(rùn)狀態(tài)[3]，可以及時(shí)、準(zhǔn)確和有效地監(jiān)測(cè)不同地表面的干旱情況，且數(shù)據(jù)獲取途徑較多，過(guò)程簡(jiǎn)單[4]，因此被諸多學(xué)者進(jìn)行研究。如王美林等利用2000—2015年的長(zhǎng)時(shí)間序列MODIS數(shù)據(jù)提取地表溫度和植被指數(shù)數(shù)據(jù)構(gòu)建TVDI，從年際尺度到季節(jié)尺度反演、解析了瑪曲表層土壤濕度的時(shí)空演變特征，并通過(guò)氣象數(shù)據(jù)及其他多源數(shù)據(jù)進(jìn)行驗(yàn)證，證實(shí)了其結(jié)果的有效性和可靠性[5]；劉英等將TVDI用于陜西省的干旱監(jiān)測(cè)，并探究了其引起干旱的主導(dǎo)因素[6]。雖然TVDI在干旱監(jiān)測(cè)中有一定作用，但MODIS時(shí)間分辨率較低，難以進(jìn)行實(shí)時(shí)旱情監(jiān)測(cè)。

大氣可降水量(precipitable water vapor，PWV)作為大氣中重要的成分之一，獲取手段主要通過(guò)全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)、氣象衛(wèi)星以及探空站等[7]，其中傳統(tǒng)的探空站監(jiān)測(cè)精度較高，但站點(diǎn)數(shù)量少且獲取頻率低，缺乏一定的連續(xù)性[8]，GNSS遙感的出現(xiàn)解決了這一問(wèn)題。由于其全天候、時(shí)空分辨率高、受云雨因素影響較小等諸多優(yōu)點(diǎn)成為 PWV監(jiān)測(cè)的主要方式[9]。PWV是大氣中產(chǎn)生降水的基礎(chǔ)，也是評(píng)估空中水資源含量的重要依據(jù)，與氣象現(xiàn)象密切相關(guān)[10]，因而備受研究人員關(guān)注。F Alshawaf發(fā)現(xiàn)PWV的變化與地表溫度呈一定規(guī)律，即地表溫度每上升1℃，PWV在一定范圍內(nèi)產(chǎn)生波動(dòng)[11]；Wang X M利用PWV監(jiān)測(cè)出了澳大利亞的洪澇災(zāi)害[12]。上述研究表明，PWV在氣象災(zāi)害監(jiān)測(cè)方面具有很大潛力，但在旱災(zāi)監(jiān)測(cè)領(lǐng)域研究較少。

MODIS數(shù)據(jù)廣泛應(yīng)用于干旱監(jiān)測(cè)領(lǐng)域，能夠?qū)崿F(xiàn)較大范圍的干旱監(jiān)測(cè)，但其時(shí)間分辨率較低，GNSS PWV作為一種新的技術(shù)手段，在氣象災(zāi)害監(jiān)測(cè)方面具有很大潛力，且有著高時(shí)空分辨率。為此，本文以云南省為例，利用MODIS的植被指數(shù)產(chǎn)品(normalized difference vegetation index，NDVI)和地表溫度產(chǎn)品(land surface temperature，LST)構(gòu)建TVDI，以驗(yàn)證PWV在干旱監(jiān)測(cè)中的適用性，同時(shí)進(jìn)行干旱特征時(shí)空變化分析。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

云南省位于我國(guó)西南地區(qū)，地處N21°8′～29°15′，E97°31′～106°11′之間。地勢(shì)上，全省地勢(shì)較高，最高海拔高達(dá)6 000多米，呈西北高、東南低態(tài)勢(shì)，為山地高原地形，主要以山地類型為主，約占比全省國(guó)土總面積的84%。氣候上，主要為亞熱帶高原季風(fēng)型氣候，全省氣溫七月達(dá)最高，約為20 ℃，最低在一月份，約為7 ℃，年溫差約為11 ℃；全年干濕分明，有著明顯的季節(jié)性和區(qū)域性降水不均現(xiàn)象，時(shí)間上，降雨主要集中在5—10月，而11月—次年4月降雨較少，空間上，雨量較多地區(qū)可達(dá)2 200～2 700 mm，較少地區(qū)可達(dá)584 mm，表現(xiàn)為夏秋多雨、春冬多旱現(xiàn)象，干旱發(fā)生頻率較高區(qū)域主要在云南東部。圖1為研究區(qū)范圍及CORS站、氣象站分布圖。

圖1 研究區(qū)范圍及CORS站、氣象站分布圖Fig.1 The study area and the distribution map of CORS stations and weather stations

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

1)本文使用的MODIS來(lái)源于LAADS DAAC官網(wǎng)(https://ladsweb.modaps.eosdis.nasa.gov/search/)，分別為8 d合成、1 km分辨率的MOD11A2地表溫度產(chǎn)品和16 d合成、1 km分辨率的MOD13A2植被指數(shù)產(chǎn)品，時(shí)間范圍為2016年—2020年的1—5月(本文將1—5月稱為春季)。利用MRT軟件分別對(duì)目標(biāo)數(shù)據(jù)(NDVI、LST)批量提取、拼接和重投影，格式轉(zhuǎn)換，然后用云南省行政區(qū)劃邊界進(jìn)行掩模裁剪，并同時(shí)對(duì)NDVI和LST進(jìn)行S-G濾波[13]以消除噪音影響，且利用最大值合成[14]法合成為月、春季數(shù)據(jù)。

2)本文使用的GNSS數(shù)據(jù)及氣象數(shù)據(jù)由中國(guó)氣象局氣象探測(cè)中心提供，站點(diǎn)主要包含云南省的35個(gè)連續(xù)運(yùn)行參考站(continuously operating reference stations,CORS)及127個(gè)氣象站，空間分布如圖1。該數(shù)據(jù)主要包含溫度、相對(duì)濕度、PWV、降雨量、氣壓等。

2 研究方法

2.1 溫度植被干旱指數(shù)

TVDI最早由Sandholt[15]提出，該指數(shù)主要考慮了2個(gè)描述土壤表層特征的重要參數(shù)，即地表溫度(LST)和歸一化植被指數(shù)(NDVI)。當(dāng)研究區(qū)植被覆蓋情況為完全裸土到完全覆蓋時(shí)，土壤濕度則由重旱到濕潤(rùn)，理論上，此時(shí)NDVI與LST構(gòu)建的特征空間呈梯形[16]，如圖2。

圖2 Ts-NDVI特征空間Fig.2 The Ts-NDVI feature space

通過(guò)NDVI和LST可構(gòu)建TVDI[17]，其關(guān)系可表示為：

(1)

式中：LST為地表溫度；LSTmax為在NDVI值下的最高溫度值，即特征空間的干邊；LSTmin為在相同的NDVI下的最低溫度值，即特征空間的濕邊；TVDI值域范圍為(0，1)，TVDI越趨向于0，表示土壤濕度越高，植被的蒸散作用增強(qiáng)，使得地表溫度下降，TVDI越趨向于1，表示土壤濕度越低，植被蒸散作用降低，地表溫度升高。通過(guò)線性擬合可得其干、濕邊方程：

(2)

式中：a1，b1，a2，b2分別為干、濕邊方程的擬合系數(shù)。

2.2 GNSS PWV反演計(jì)算

GNSS衛(wèi)星發(fā)射的信號(hào)在穿過(guò)大氣層時(shí)，由于電離層和對(duì)流層的影響，產(chǎn)生信號(hào)延遲，記為大氣總延遲(zenith total delay，ZTD)。ZTD由電離層延遲和對(duì)流層延遲組成，其中電離層折射引起的延遲可通過(guò)雙頻接收機(jī)消除99%的影響[9]；對(duì)流層延遲主要有靜力延遲(zenith hydrostatic delay，ZHD)和濕延遲(zenith wet delay，ZWD)，即ZWD=ZTD-ZHD，ZTD計(jì)算公式為：

(3)

式中：k1，k2，k3分別為大氣折射常數(shù)；Rd為干空氣氣體常數(shù)；ρ為干空氣總質(zhì)量密度；Pw和Zw分別為水汽局部氣壓和可壓縮系數(shù)；H和H0分別為對(duì)流層頂層高度及站點(diǎn)高度。式(3)等式右側(cè)兩項(xiàng)分別表示ZHD和ZWD，其中ZHD可通過(guò)站點(diǎn)信息及地面氣壓求得[18]，公式為：

(4)

f(φ,H)=1-0.002 66cos(2φ)-0.000 28H。

(5)

式(4)為解算ZHD的常用模型，Saastamonien模型。式中：P0為測(cè)站地面氣壓，單位為hPa；φ為測(cè)站地理緯度；H為測(cè)站點(diǎn)海拔，單位為km。由于氣壓測(cè)量精度較高，該模型所估算的ZHD可達(dá)mm級(jí)精度[19]。計(jì)算公式為：

(6)

式中：ρω為水密度;R=461(J·kg-1·K-1);κ=(3.776±0.014)×105(K2·hPa-1);κ'=16.48(K2·hPa-1);Tm為大氣加權(quán)平均溫度;Π為無(wú)量綱水汽轉(zhuǎn)換系數(shù)，僅與大氣加權(quán)平均溫度Tm有關(guān)，常用取值范圍為6.0～6.5[20]。

2.3 相關(guān)性分析

由中國(guó)氣象局氣象探測(cè)中心提供的云南省的35個(gè)CORS站數(shù)據(jù)和128個(gè)氣象站點(diǎn)數(shù)據(jù)，包括降雨量、溫度、PWV、相對(duì)濕度。借助Pearson相關(guān)系數(shù)分析TVDI與PWV、降雨量、溫度及相對(duì)濕度的之間的相關(guān)關(guān)系。公式可表示為：

(7)

式中：R表示Pearson相關(guān)系數(shù)；N代表樣本個(gè)數(shù)；xi，yi分別表示TVDI和氣象因子；R取值范圍為[-1,1]，|R|值越大相關(guān)性越強(qiáng)。當(dāng)R<0時(shí)，表示兩變量呈負(fù)相關(guān)關(guān)系，反之則表示為正相關(guān)，R=0時(shí)表示兩變量之間不相關(guān)。

3 結(jié)果與分析

3.1 溫度植被干旱指數(shù)反演

由于植被覆蓋度過(guò)低或者過(guò)高，均會(huì)對(duì)監(jiān)測(cè)結(jié)果產(chǎn)生影響。過(guò)低時(shí)，NDVI難以顯示區(qū)域的真實(shí)植被覆蓋情況；過(guò)高時(shí)，植被達(dá)到過(guò)飽和，生長(zhǎng)緩慢。因此，本文NDVI取0.2～0.8，以0.01為步長(zhǎng)，計(jì)算每個(gè)NDVI像元對(duì)應(yīng)的最大、最小地表溫度，然后對(duì)每期特征空間的干、濕邊進(jìn)行擬合，便可得到Ts-NDVI特征空間，如圖3，擬合結(jié)果見表1，最后將擬合方程代入式(1)即可求得每個(gè)像元的TVDI。

根據(jù)圖3可知，2016—2020年每1—5月特征空間都具有相似的梯形形狀，干、濕邊都具有相同的變化趨勢(shì)，其中地表溫度最大/最小值隨著NDVI的增大而減小/增大；由表1可知，LSTmax與NDVI呈負(fù)相關(guān)，且干邊擬合較好，擬合系數(shù)基本均高達(dá)0.8，但是每年1月，擬合系數(shù)較小，呈現(xiàn)弱相關(guān)性，根據(jù)孫麗等的研究[21]可知，這主要是由于植物1月份生長(zhǎng)期緩慢，影響了NDVI的真值。

(a)2016年1月 (b)2016年2月 (c)2016年3月 (d)2016年4月 (e)2016年5月

表1 特征空間干濕邊擬合方程及相關(guān)系數(shù)(2016—2020)Tab.1 The fitting equation and correlation coefficient of dry and wet edge in characteristic space (2016—2020)

(續(xù)表)

3.2 TVDI與PWV相關(guān)性分析

在關(guān)于TVDI和PWV的相關(guān)性分析方面，主要從時(shí)間和空間上做具體分析。其中在時(shí)間上表現(xiàn)在季尺度、月尺度和日尺度。在季尺度上，TVDI與PWV相關(guān)性較高，相關(guān)系數(shù)基本均大于0.5；在月尺度上，對(duì)研究區(qū)內(nèi)各個(gè)時(shí)期的TVDI和PWV的月均值做統(tǒng)計(jì)分析。結(jié)果如圖4所示，TVDI與PWV月值具有一致變化規(guī)律，每年特征表現(xiàn)呈先上升后下降趨勢(shì)，分別于3，4月達(dá)到峰值，造成這一現(xiàn)象的原因可能在于該時(shí)間段內(nèi)高溫少雨，空氣中PWV含量較低導(dǎo)致PWV變化較小，此外在2020年2—5月出現(xiàn)了步調(diào)相反的變化趨勢(shì)，這與統(tǒng)計(jì)資料顯示的2020年特大春旱情況相符合；在日尺度上，由于PWV與降雨量具有較強(qiáng)的瞬時(shí)性，月均值難以表現(xiàn)PWV變化的細(xì)節(jié)特征，因此選取并計(jì)算了云南墨江(YNMJ)站PWV與降雨量的逐日值進(jìn)行分析，如圖5，在逐日的時(shí)間序列中，每次降雨必然會(huì)伴隨著PWV的陡然上升和下降的趨勢(shì)，此時(shí)TVDI也有一定程度的降低，當(dāng)降雨較少或者不降雨時(shí)，TVDI又逐漸增大，表明PWV的變化具有一定的干旱特征信號(hào)[22]。在空間上，通過(guò)對(duì)兩者進(jìn)行Pearson相關(guān)性分析，結(jié)果見表2，PWV與TVDI之間具有較強(qiáng)的相關(guān)性，表現(xiàn)出相似的空間分布特征。

表2 TVDI與PWV相關(guān)性系數(shù)(2016—2020)Tab.2 Correlation coefficient between TVDI and PWV (2016—2020)

圖4 PWV與TVDI變化趨勢(shì)圖(2016—2020)Fig.4 Change trend chart of PWV and TVDI (2016—2020)

圖5 云南墨江(YNMJ)站TVDI、降雨量、PWV變化趨勢(shì)圖Fig.5 Change trend chart of TVDI,rainfall and PWV of YNMJ Station

3.3 干旱時(shí)空特征分析

根據(jù)前文計(jì)算的TVDI，采用齊述華等[23]對(duì)全國(guó)旱情監(jiān)測(cè)的等級(jí)劃分法對(duì)其進(jìn)行等級(jí)劃分，分為濕潤(rùn)、正常、輕旱、中旱及重旱5個(gè)等級(jí)，由此便獲得云南省干旱等級(jí)分布圖。如圖6，在時(shí)間上，2016—2020年5 a間，旱情變化趨勢(shì)一致，旱情程度逐月上升，在每年1，2月相對(duì)緩和，但在3，4月均較為嚴(yán)重，且均主要分布在滇南、滇西南地區(qū)，在5月有緩和趨勢(shì)，主要是由于由春入夏，雨量相對(duì)增多，因此旱情有所緩解。年際間，每年旱情均表現(xiàn)出一定程度的春旱，這與曹影等[24]的研究結(jié)果具有一致的特征。

(a)2016年1月 (b)2016年2月 (c)2016年3月 (d)2016年4月 (e)2016年5月

在空間上，根據(jù)云南省區(qū)域特征將云南行政區(qū)域劃分為滇中、滇東北、滇東、滇東南、滇南、滇西南、滇西、滇西北，并利用Python，以上述為矢量范圍掩模，對(duì)TVDI進(jìn)行分幅裁剪，以獲取各區(qū)域TVDI分布情況，然后對(duì)每個(gè)區(qū)域輕旱、中旱及重旱面積進(jìn)行統(tǒng)計(jì)。從圖7可以看出云南省在2016—2020年各區(qū)域干旱面積占比及空間分布特點(diǎn)，大部分地區(qū)均以中旱為主，尤其是滇中、滇東地區(qū)，中旱面積分別高達(dá)68%和80%，滇東北主要以輕旱為主，且輕旱面積逐年上升，滇西南、滇東南地區(qū)主要以重旱為主，其中滇西南重旱面積有逐年上升趨勢(shì)，2020年高達(dá)63.64%，滇東南呈現(xiàn)逐漸下降趨勢(shì)，由2016年的65.43%降低到2020年的26.13%。

圖7 云南省各區(qū)輕旱、中旱及重旱面積占比Fig.7 The proportion of light drought,medium drought and severe drought in Yunnan Province

3.4 TVDI，PWV與氣象因子的關(guān)系

為了進(jìn)一步研究溫度、相對(duì)濕度等氣象因子與植被干旱指數(shù)和PWV之間的關(guān)系，本文計(jì)算并提取了云南省2016—2020年1—5月各CORS站點(diǎn)的溫度、相對(duì)濕度及PWV，觀察其時(shí)間序列變化特點(diǎn)。結(jié)果顯示：如圖8所示，1，2月溫度約以0.48℃/年、0.32℃/年的速度上升。同時(shí)以CORS站為中心，計(jì)算其周圍3像元×3像元范圍內(nèi)的TVDI值。最后將各TVDI值、PWV值分別與溫度、相對(duì)濕度進(jìn)行Pearson相關(guān)性分析，分析結(jié)果如圖9，TVDI，PWV均與溫度呈現(xiàn)正相關(guān)關(guān)系，相關(guān)系數(shù)高達(dá)0.85，且通過(guò)了P<0.01顯著性檢驗(yàn)，但與相對(duì)濕度呈弱相關(guān)或不相關(guān)，這表明溫度與TVDI和PWV密切相關(guān)，可作為干旱的重要評(píng)價(jià)因子。

表3 干旱等級(jí)分級(jí)Tab.3 Drought grade

圖8 溫度、相對(duì)濕度變化趨勢(shì)Fig.8 Trend of temperature and relative humidity

(a)TVDI與溫度相關(guān)性 (b)TVDI與相對(duì)濕度相關(guān)性

4 結(jié)論

針對(duì)近幾年嚴(yán)峻的干旱形勢(shì)，本文利用TVDI分析了云南省2016—2020年春季的干旱特征時(shí)空變化，同時(shí)通過(guò)均勻分布在云南省的35個(gè)CORS站解算各站點(diǎn)PWV，并利用TVDI驗(yàn)證了PWV在干旱監(jiān)測(cè)中的適用性，得出以下結(jié)論：

1)TVDI能較好地監(jiān)測(cè)干旱情況，在云南省具有較好的適用性。由于區(qū)域性不均勻降水導(dǎo)致云南省干旱常年呈滇西北向滇東南增強(qiáng)趨勢(shì)，主要集中在滇中、滇東、滇東南、滇南地區(qū)；由于時(shí)間性不均勻降水導(dǎo)致云南省干旱主要以冬旱、春旱為主，季內(nèi)旱情呈先遞增后減緩趨勢(shì)，尤其3—4月份旱情變化特征最為明顯，年間旱情呈波動(dòng)變化特點(diǎn),沒(méi)有出現(xiàn)明顯減少或增加的趨勢(shì)。

2)GNSS PWV在干旱監(jiān)測(cè)領(lǐng)域中具有一定潛力。基于Pearson相關(guān)分析發(fā)現(xiàn)PWV和TVDI存在較強(qiáng)的相關(guān)性，在季尺度上，相關(guān)系數(shù)基本均大于0.5；在月尺度上，PWV變化趨勢(shì)與TVDI變化趨勢(shì)基本一致，但TVDI變化有一定的時(shí)間延遲；在日尺度上，尤其時(shí)降雨時(shí)期，PWV變化和TVDI變化幅度契合度更高，表現(xiàn)出了一定的干旱特征信號(hào)，因此PWV為旱災(zāi)監(jiān)測(cè)提供了一種新的技術(shù)手段。

本研究結(jié)果較好地反映了云南省近5 a的干旱演變特征，為防災(zāi)減災(zāi)提供了理論依據(jù)，有一定的參考價(jià)值。同時(shí)驗(yàn)證了GNSS PWV在干旱響應(yīng)上具有一定潛力，但本文僅做了一些定性分析，對(duì)于GNSS PWV在干旱監(jiān)測(cè)中的定量化分析還需深入研究。