MODIS數(shù)據(jù)反演三江平原晴空上午陸表溫度比較研究

孔繁艷

（1.南京大學(xué)大氣科學(xué)學(xué)院，南京 210093；2.黑龍江省七臺河市氣象局，黑龍江 七臺河 154600）

三江平原開墾后水熱平衡發(fā)生顯著變化，早在上世紀(jì)80年代，就有對三江平原輻射平衡和小氣候特征的研究[1－2]。后期研究包括王學(xué)雷對三江平原濕地進行的土壤熱學(xué)特性研究[3]、張蕓對三江平原沼澤濕地開墾后的熱量平衡變化的研究[4]等。在應(yīng)用遙感技術(shù)對三江平原研究方面，有張樹清對三江平原濕地的遙感分類研究[5]、靳華安應(yīng)用CBERS－02衛(wèi)星數(shù)據(jù)對三江平原濕地植被葉面積指數(shù)估算[6]等。陸表溫度（Land surface temperature，LST）是指示陸面水熱狀況的重要參數(shù)，大范圍長期的陸面溫度對比觀測較難實現(xiàn)，本文嘗試使用4年MODIS數(shù)據(jù)反演開墾區(qū)和未開墾區(qū)LST，為研究三江平原開墾帶來的水熱狀況變化提供更為精確的定量依據(jù)。

1 數(shù)據(jù)來源及處理

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

本研究使用的MODIS數(shù)據(jù)來源于NASA MODIS網(wǎng)站的Terra衛(wèi)星上午降軌數(shù)據(jù)。從2006～2009年全部L1B數(shù)據(jù)集中篩選出可用軌道55條共67個5 min數(shù)據(jù)塊（Granule）。數(shù)據(jù)分布情況為2006年21個，2007年17個，2008年15個，2009個14，合計67個。

1.2 數(shù)據(jù)采樣點分布

在開墾和未開墾區(qū)域各自選取200個采樣數(shù)據(jù)點。開墾區(qū)域數(shù)據(jù)點從我國黑龍江省撫遠(yuǎn)縣、饒河縣、同江市、富錦市、綏濱縣和蘿北縣境內(nèi)選取，未開墾區(qū)數(shù)據(jù)點從黑龍江北岸烏蘇里江東岸俄羅斯境內(nèi)平坦荒原選取。數(shù)據(jù)點選取采用人工交互方法，參照2007年7月3日上午10:30 Terra MODIS（數(shù)據(jù)集名稱OD021KM.A2007184.0150.005.2010 200125252.hdf）真彩色合成圖。數(shù)據(jù)采樣點空間分布見表1。

開墾區(qū)與未開墾區(qū)數(shù)據(jù)采樣點平均緯度差為0.47°，可以認(rèn)為開墾區(qū)和未開墾區(qū)數(shù)據(jù)樣點緯度差異較小，由于緯度差異對溫度造成的影響可以忽略。

表1 反演LST使用數(shù)據(jù)的空間分布Table 1 Range of data used in land surface temperature retrieval

1.3 數(shù)據(jù)時間

本研究采用的Terra衛(wèi)星探測數(shù)據(jù)時間均為地方時上午，由于衛(wèi)星軌道的準(zhǔn)太陽同步特點，過境時間略有差異。最早9:35，最晚10:30，大部分集中在10:15左右（均為北京時間）。

1.4 數(shù)據(jù)處理

數(shù)據(jù)處理使用了RSD遙感數(shù)據(jù)處理軟件。投影采用阿爾伯斯等積圓錐投影。在數(shù)據(jù)柵格化之前的SWATH數(shù)據(jù)基礎(chǔ)上進行幾何位置精校正，由于數(shù)據(jù)區(qū)域范圍不大，僅適用一次線性幾何校正就已經(jīng)達(dá)到了位置精度要求。

雖然MODIS L1B數(shù)據(jù)集是經(jīng)過刪選后獲得，但仍無法避免研究區(qū)域有部分的云遮擋。從400個數(shù)據(jù)樣點回取數(shù)據(jù)時難免取到云數(shù)據(jù)，給研究帶來不利影響，因此必須進行云去除處理。使用RSD的剪切工具，交互式從每幅圖面上對云以及相應(yīng)的陰影區(qū)域剪切。

對67個MODIS L1B數(shù)據(jù)集投影重采樣形成柵格數(shù)據(jù)后，應(yīng)用RSD軟件樣本數(shù)據(jù)分析功能，依據(jù)樣本點的經(jīng)緯度位置，從67個柵格數(shù)據(jù)文件中提取樣本數(shù)據(jù)供分析使用。

基于LST季節(jié)變化大于同日序年際變化的假設(shè)，將所有MODIS數(shù)據(jù)按一年中日序進行排序，并用滑動平均的方法消除年際之間差異，以彌補由于云遮擋等原因造成的MODIS逐日數(shù)據(jù)缺失。

2 LST反演算法

常用的MODIS數(shù)據(jù)反演LST算法有Price算法[7]、Becker－Li算法[8]、推廣的分裂窗算法[9]和 Kerr算法、Prata－Platt算法、Coll算法等。本研究使用了推廣的分裂窗算法。

2.1 LST推廣的分裂窗算法

推廣的分裂窗算法是在Becker－Li算法基礎(chǔ)上改進而來的，是MODIS溫度產(chǎn)品的官方算法之一，其公式為

式中，Ts為陸表溫度，T31是MODIS 31通道亮溫，T32是MODIS 32通道亮溫，ε為兩個熱紅外通道的平均比輻射率。A0、A1、A2、A3、B1、B2、B3均為參數(shù)，由衛(wèi)星數(shù)據(jù)和地面實測數(shù)據(jù)回歸求得，本文使用適用于AVHRR數(shù)據(jù)的參數(shù)代替，具體數(shù)值為 A01.274，A11.0，A20.15616，A30.482，B16.26，B23.98，B338.33。

2.2 亮溫計算

亮溫的求算公式：

式中，c2=1.438833×10－2mK，c1=1.1910659×10－22W·m－2，L為定標(biāo)后探測值（W·m－2·μm－1·s－1）。λ是波長（μm）。

2.3 熱紅外通道平均比輻射率

平均比輻射率的計算使用了Sobrino方法[10]。該方法引自NOAA AVHRR比輻射率計算，分三種情況：

①INDV＜0.2時，地表為裸土或稀疏植被，植被覆蓋率幾乎為零，ε=0.980－0.042ρ1，Δε=－0.03－0.029ρ2。ε是兩個長波波段的平均比輻射率，Δε為兩個長波波段的比輻射率之差，INDV為植被指數(shù)，ρ1和ρ2分別表示波段1和波段2的反射率；

② 0.2≤INDV≤0.5，植被與裸土混合，植被覆蓋率不為零，則ε=0.971+0.018Pv，Δε=－0.006表示植被覆蓋程度。

③INDV＞0.5，地表完全被植被覆蓋，植被覆蓋率近似為1，則ε=0.985，Δε=0。

3 結(jié)果與分析

3.1 LST年變化

首先，將開墾區(qū)和未開墾區(qū)上午LST反演結(jié)果按照一年對照，分析年變化。

如圖1所示，為開墾區(qū)與未開墾區(qū)LST的年變化圖，圖中紫線部分為開墾區(qū)數(shù)據(jù)，藍(lán)線部分為未開墾區(qū)數(shù)據(jù)。從圖中可以看出，無論是開墾區(qū)還是未開墾區(qū)，LST變化趨勢相同，地表溫度表現(xiàn)出的隨年季變化曲線較為一致。日序為100（約為3月份）～300（約為10月份）之間MODIS數(shù)據(jù)反演出的上午晴空LST在絕對零度以上。

圖1 開墾區(qū)/未開墾區(qū)LST年變化Fig.1 LST annual variation in reclaimed and unreclaimed area

3.2 ΔLST年變化

為比較開墾區(qū)和未開墾區(qū)之間的LST差異情況，將反演所得的開墾區(qū)和未開墾區(qū)LST值相減求得逐日溫差。將這55個時間序列的開墾區(qū)和未開墾區(qū)上午LST取平均值，開墾區(qū)和未開墾區(qū)上午LST平均值分別為16.28和15.97，相差僅0.31。考慮到開墾區(qū)和未開墾區(qū)數(shù)據(jù)樣點的平均緯度差異，平均LST差異不明顯，但在季節(jié)分布上差異明顯。開墾后區(qū)域與未開墾區(qū)域上午LST差ΔLST按一年中日序排列結(jié)果見表2。

表2 開墾后區(qū)域與未開墾區(qū)域上午LST差Table 2 LST Difference(ΔLST)between Reclaimed and unreclaimed area(°C)

可見ΔLST最大負(fù)值出現(xiàn)在第130天，為－7.07，最大正值出現(xiàn)在第184天，為7.59。滑動平均后的逐日ΔLST分布（圖1）表現(xiàn)為表現(xiàn)為兩谷一峰。具體特征為：

①自3月中下旬開始至5月中旬表現(xiàn)為第一個溫度谷。這個時期陸表冰雪以及凍土融化，開墾區(qū)和未開墾區(qū)均少有植物生長。開墾區(qū)以裸露土壤為主，上午升溫較慢，未開墾區(qū)有大量枯萎植物覆蓋上午升溫較快，ΔLST為負(fù)值。

②5月中旬至7月下旬表現(xiàn)為一個溫度峰。這一時期是農(nóng)田播種和作物旺盛生長的前期，太陽輻射增加。未開墾區(qū)植物大量生長，上午升溫較慢。開墾區(qū)仍然以裸露土壤為主，上午晴天陸表急劇升溫，表現(xiàn)為一個很強的溫度峰。

③7月下旬至8月末表現(xiàn)為第二個溫度谷，是農(nóng)田作物生長茂盛的時期，作物的植被效應(yīng)強于未開墾荒原的植被，是一個較弱的溫度谷。

④9月以后農(nóng)田作物和荒原植被開始枯黃，入冬至次年3月大部分時間無論開墾區(qū)或未開墾區(qū)均為凍土或被冰雪覆蓋，ΔLST變化不明顯。

圖2 開墾后區(qū)域與未開墾區(qū)域上午ΔLSTFig.2 Morning ΔLST of Reclaimed and unreclaimed area

4 討論

a.MODIS的高光譜特點能夠反演除LST外的多種陸表參數(shù)，本研究為應(yīng)用衛(wèi)星遙感數(shù)據(jù)綜合研究三江平原陸表參數(shù)變化的部分做基礎(chǔ)工作。由于時間以及技術(shù)水平的不足，研究中存在不足：①本文反演出的LST并未涉及下午和夜間LST變化，有待結(jié)合Aqua數(shù)據(jù)進一步研究；②采用光學(xué)遙感衛(wèi)星數(shù)據(jù)反演LST只能獲取晴天數(shù)據(jù)，對于有云或陰影遮擋處的數(shù)據(jù)本研究中采取直接剔除的方式，這使得數(shù)據(jù)的代表性受到一定限制；③研究中對于LST的反演方法，直接選擇應(yīng)用較為廣泛的Becker－Li分裂窗算法，對于其它方法的適用性以及反演效果并未探討；④對于反演精度問題沒有涉及，直接將反演LST的結(jié)果進行對比分析。

b.本文的創(chuàng)新之處在于系統(tǒng)的大規(guī)模使用遙感數(shù)據(jù)來研究三江平原開墾區(qū)與未開墾區(qū)地溫變化，首次進行大范圍的取點對比研究，彌補了利用小范圍未開墾區(qū)進行此類工作因受周邊環(huán)境影響而導(dǎo)致對未開墾環(huán)境的失真的不足，并填補了在俄羅斯境內(nèi)地溫實況數(shù)據(jù)的空白，使研究更具代表性。

c.這種通過衛(wèi)星數(shù)據(jù)獲取大范圍陸表參數(shù)的方法為區(qū)域內(nèi)氣象因子或其他因子研究提供一種新手段，相較于傳統(tǒng)數(shù)據(jù)采集方式具有數(shù)據(jù)獲取及時、數(shù)據(jù)量大以及數(shù)據(jù)覆蓋范圍廣泛等優(yōu)點。

d.本文只應(yīng)用MODIS傳感器進行研究，另有多種衛(wèi)星傳感器也可進行LST反演的計算。在以后的研究工作中可著眼于應(yīng)用我國自己的衛(wèi)星數(shù)據(jù)進行LST等陸表參數(shù)的反演。由此可見，這種通過衛(wèi)星數(shù)據(jù)進行參數(shù)反演對區(qū)域內(nèi)氣象條件變化進行研究、在一定程度上指導(dǎo)實踐的方式必將成為未來發(fā)展的前沿方向。

5 結(jié) 論

三江平原開墾后農(nóng)田晴天上午LST變化發(fā)生明顯改變，表現(xiàn)為兩谷一峰型，晴天上午LST年均值變化不明顯。

a.無論是開墾區(qū)還是未開墾區(qū)，LST變化趨勢相同，表現(xiàn)出的隨年季變化曲線較為一致，反演結(jié)果與實際情況大致表現(xiàn)出相同的勢。

b.年中開墾區(qū)與未開墾區(qū)晴空上午LST差值，即ΔLST的變化表現(xiàn)出負(fù)－正－負(fù)的形式，ΔLST存在按照年季間變化的特點。具體為：LST差年變化表現(xiàn)為兩谷一峰：春季解凍后未開墾區(qū)有大量枯萎植物覆蓋，開墾后區(qū)域裸土為主，表現(xiàn)為溫度谷；5月下旬～7月中下旬，未開墾區(qū)植物大量生長，開墾區(qū)農(nóng)田前期多為裸土或植株矮小，表現(xiàn)為一個很強的溫度峰；7月末～8月末農(nóng)田作物茂密生長，表現(xiàn)為另一個溫度谷；9月以后～次年4月份大部分時間為凍土或被冰雪覆蓋，開墾區(qū)和未開墾區(qū)LST差異不明顯。

[1]劉興土.三江平原的輻射氣候特征[J].地理科學(xué),1983(1):27－32.

[2]劉興土.三江平原沼澤輻射平衡與小氣候基本特征[J].地理科學(xué),1988(2):126－127.

[3]王學(xué)雷.沼澤土壤熱學(xué)特性研究[J].地理科學(xué),1993(1):85－86.

[4]張蕓,呂憲國,朱誠.三江平原沼澤濕地開墾后的熱量平衡變化[J].南京大學(xué)學(xué)報:自然科學(xué)版,2002,38(6):11.

[5]張樹清,陳春,萬思璞.三江平原溫地遙感受分類模式研究[J].遙感技術(shù)與應(yīng)用,1999,14(1):54－58.

[6]靳華安,等.三江平原濕地植被葉面積指數(shù)遙感估算模型[J].生態(tài)學(xué)雜志,2008(5):803－808.

[7]Price J C.Land surface temperature measurements from the split window channels of the NOAA7 advanced very high resolution radiometer[J].Journal of Geophysical Research,1984,89:7231－7237.

[8]Becker F,Li Z L.Towards a local split window method over land surface[J].International Journal of Remote Sensing,1990(3):369－393.

[9]Wan Z,Dozier J.A generalized split－window algorithm for retrieving land surface temperature from space[J].IEEE Transctions on Geosciences and Remote Sensing,1996,34:892－905.

[10]Sobrino J A,Raissouni N,Li Z L.A comparative study of land surface emissivity retrieved from NOAA data[J].Remote Sensing of Environment,2001,75(1):256－266.