古浪非均勻近地層地表溫度和感熱通量的衛(wèi)星反演

古浪非均勻近地層地表溫度和感熱通量的衛(wèi)星反演

1地表溫度反演的地表分裂窗算法表面溫度是地氣系統(tǒng)中能量和物質(zhì)交換的重要影響因素。它不僅是衡量全球溫度變化的重要指標，也是區(qū)域和全球規(guī)模中區(qū)域和方向的重要因素。隨著20世紀中期TIROS-Ⅱ衛(wèi)星發(fā)射,利用衛(wèi)星熱紅外資料反演地表溫度開始發(fā)展起來。最初的反演工作僅限于海面溫度,但是隨著大氣輻射傳輸模型的改進及NOAA/AVHRR、MODIS等遙感數(shù)據(jù)的豐富,利用分裂窗算法反演海面溫度的精度提高到1.0K以內(nèi)。海面溫度反演成功之后,利用衛(wèi)星熱紅外數(shù)據(jù)反演陸地表面溫度變得更加具有誘惑力。Price在1984年首先將海面溫度反演的分裂窗法引用到地表溫度的反演中,給出了地表近似黑體時的地表溫度反演算法。隨后,1990年Beckeretal結(jié)合考慮通道比輻射率差異的影響,在對輻射傳輸方程線性近似的基礎(chǔ)上,給出了一個地表溫度反演的局地分裂窗算法(下稱Becker算法),并在1995年考慮水汽的影響對算法做了進一步改進。1992年Kerretal考慮植被的影響,將植被指數(shù)引進算法,得到一個半經(jīng)驗模型的局地分裂窗算法(下稱Kerr算法)。2005年覃志豪等利用MO-DIS數(shù)據(jù)將局地分裂窗算法進行修正,得到一個針對MODIS數(shù)據(jù)反演地表溫度的分裂窗算法(下稱覃志豪算法)。另外,一些研究者也在地表溫度反演上做了不同的嘗試。隨著衛(wèi)星熱紅外通道數(shù)的增加,多通道、多角度方法迅速發(fā)展起來,有的精度能達到1.0K左右。但是由于其計算的復雜性和算法還不夠成熟,所以還未被廣泛使用。而利用熱紅外遙感數(shù)據(jù)反演地表溫度的分裂窗算法是目前地表溫度反演算法中發(fā)展最為成熟的,反演精度一般在3.0K以內(nèi)。在感熱通量的遙感估算中,常用的動力學阻抗法涉及兩個關(guān)鍵因子,即地表溫度和空氣動力學阻抗。如何利用有限的遙感數(shù)據(jù)和觀測資料得到更高精度的感熱通量結(jié)果,研究者們做出了不同的參數(shù)化方案。1993年郭學良等直接利用觀測的感熱通量與地表溫度和氣溫的差值進行線性擬合,得到感熱通量估算的經(jīng)驗公式。1999年Hurtalovaetal根據(jù)近地層風廓線公式,利用不同高度的風速觀測資料得到計算感熱通量的各參量,但這種方法需要至少3個高度的風速資料。Jiaetal利用ATSR衛(wèi)星數(shù)據(jù),使用歸一化植被指數(shù)參數(shù)化動量粗糙度來計算感熱通量,將估算的精度提高到25.5W·m-2。黃妙芬又引入葉面積指數(shù)和土壤調(diào)整植被指數(shù)對粗糙度和零平面位移進行參數(shù)化。在感熱通量的計算中,在已知植被高度的情況下,為簡單起見,可采用植被高度對粗糙度進行粗略的計算。例如,Sellersetal提出的估算動量粗糙度的方法。本文利用野外觀測試驗資料對現(xiàn)有的地表溫度和空氣動力學阻抗算法進行了修正,給出了適用于試驗區(qū)綠洲—沙漠下墊面遙感反演地表溫度和感熱通量的方法,以期為在非均勻地區(qū)利用衛(wèi)星資料反演地表參量提供重要參考依據(jù)。2野外試驗的簡要介紹和數(shù)據(jù)的選擇2.1古浪試驗觀測及儀器布置古浪試驗區(qū)位于甘肅省武威市古浪縣海子鎮(zhèn)東新村,范圍為103°48′-103°51′E、37°37′-37°38′N,統(tǒng)稱為古浪非均勻近地層觀測試驗(GulangHeterogeneousUnderlyingSurfaceLayerExperiment,GHUSLE)。試驗區(qū)綠洲為西南—東北走向,主要種植玉米、小麥和葵花等農(nóng)作物,為典型的植被下墊面。綠洲東西兩側(cè)均為沙漠(下稱東、西沙漠)。這種獨特的自然環(huán)境為研究不同的下墊面陸面參數(shù)特征和綠洲—沙漠相互作用提供了客觀的自然條件,試驗區(qū)站點分布見圖1。本次試驗為陸面過程的研究提供了全面、連續(xù)的觀測資料。GHUSLE在沙漠—綠洲—沙漠東、西剖面的不同下墊面上布設(shè)了3個微氣象觀測站,分別是西沙漠、農(nóng)田和東沙漠站,綠洲是觀測的重點,建有一個32m的觀測塔。試驗主要分為3個部分:(1)基本觀測:東沙漠站觀測項目有2層(高度分別為2m和8m)溫濕觀測、1層(高度為8m)風平均量觀測;6層(深度為5,10,15,20,30,40cm)土壤溫度和含水量以及輻射觀測。農(nóng)田站的32m梯度塔可提供4層(高度為2,8,18,32m)的風儀濕常規(guī)觀測;3層(高度為3,12,28m)渦動相關(guān)觀測;5層(深度為5,10,20,40,80cm)土壤溫度和含水量觀測;3層(深度為5,10,20cm)土壤熱通量觀測及在4m高度上架設(shè)有光徑長度為450m的大孔徑閃爍儀(LargeApertureScintillometer,LAS)通量觀測;植被高度、氣壓和輻射等觀測。西沙漠站的觀測項目有高度為3m渦動相關(guān)儀觀測;5層(深度為5,10,20,40,80cm)土壤溫度和含水量觀測;2層(深度為5,10cm)土壤熱通量觀測;輻射和氣壓等常規(guī)氣象觀測。(2)加強觀測:在2010年6月22日-7月5日加強觀測期間,3個常規(guī)觀測站都增加了蒸發(fā)觀測和系留氣球探空觀測;農(nóng)田站還增加了雷達觀測;在3個常規(guī)觀測站之間還架設(shè)了4個向上長波輻射觀測和蒸發(fā)皿蒸發(fā)量觀測。(3)儀器平行對比觀測:2010年7月15-24日在阿拉善左旗下墊面比較平坦、均勻的沙漠荒灘上(37°45.148′N,103°55.409′E),對古浪試驗主要觀測儀器進行對比觀測。儀器平行架設(shè)在荒灘上,主要包括架設(shè)高度為2.3m溫度、濕度和風速觀測;1.5m輻射各分量觀測;2.9m渦動相關(guān)儀觀測;4.2m高度上架設(shè)有光徑長度為1010m的LAS通量觀測;5層(深度為5,10,20,40,80cm)土壤溫度和含水量觀測;3層(2,5,10cm)土壤熱通量觀測。這些資料主要用于修正使用不同儀器觀測產(chǎn)生的偏差。2.2基于不同數(shù)據(jù)產(chǎn)品的觀測資料MODIS全稱為中分辨率成像光譜儀(Moderate-ResolutionImagingSpectroradiometer),是搭載在Terra和Aqua衛(wèi)星上的一個重要傳感器,是衛(wèi)星上將實時觀測數(shù)據(jù)通過X波段向全世界直接傳播,可以免費接收數(shù)據(jù)并無償使用的星載儀器。全球許多國家和地區(qū)都在接收和使用MODIS數(shù)據(jù)。常用的MODIS數(shù)據(jù)資料格式包括PDS和HDF兩種。MODIS標準數(shù)據(jù)產(chǎn)品根據(jù)內(nèi)容的不同可分為陸地標準數(shù)據(jù)產(chǎn)品、大氣標準數(shù)據(jù)產(chǎn)品和海洋標準數(shù)據(jù)產(chǎn)品等三類主要標準數(shù)據(jù)產(chǎn)品,總計44種標準數(shù)據(jù)產(chǎn)品類型(分別以MOD01,……,MOD44命名),MOD1-3為定標產(chǎn)品。根據(jù)標準數(shù)據(jù)產(chǎn)品分級系統(tǒng)MODIS數(shù)據(jù)又劃分為5個等級。為了消除云對衛(wèi)星反演的影響,選取天氣晴好時的衛(wèi)星觀測資料。具體為2010年5月29日11:55(北京時,下同)、6月10日12:15、12日12:05、14日11:55、16日11:40、17日12:25、21日12:00、23日11:45、7月5日12:10的MODIS一級定標數(shù)據(jù)產(chǎn)品,即MOD021KM和MOD03數(shù)據(jù),時間分辨率為1~2天,空間分辨率為1km×1km。GHUSLE近地層的觀測資料主要包括基本觀測期間的氣壓、植被高度和地表輻射各分量;3m處氣溫、空氣密度、風速和感熱通量等,數(shù)據(jù)全部平均到30min。加強觀測期間每隔2h一次的長波輻射資料。內(nèi)蒙古對比觀測期間30min的通量資料和10min的氣溫、風速資料。3數(shù)據(jù)預處理3.1多通道法分析氣資料本文主要使用的是MOD02的1km數(shù)據(jù)進行計算,基于AVHRR4、5通道發(fā)展起來的分裂窗算法,主要是通過兩個相鄰通道對大氣吸收作用的差異及兩個通道亮溫的組合來剔除大氣影響,避免了衛(wèi)星資料的大氣訂正。在使用衛(wèi)星資料時,首先利用MOD03的經(jīng)緯度將衛(wèi)星資料插值到各站點。本文采用最鄰近內(nèi)插法,即觀測站周圍1km范圍內(nèi)離該測站最近的衛(wèi)星像元值作為該測站的衛(wèi)星資料。其次對衛(wèi)星資料進行云檢測,采用文獻中提到的多通道法對MODIS資料進行云檢測。檢測結(jié)果見表1,當有云存在時MODIS資料不可用。第174天和第186天兩日處在GHUSLE的加強觀測期,期間地面增加了4個移動觀測站(見圖1)。3.2儀器的使用結(jié)果)的資料進行剔除。由于各測站所使用的觀測儀器不同,所獲得的資料可能會存在一定的偏差。為了減小或消除這種偏差,本文利用內(nèi)蒙古的平行觀測數(shù)據(jù)對所使用的感熱通量、氣溫及風速資料進行了修正。在儀器平行對比觀測中,農(nóng)田站和西沙漠站儀器觀測的感熱通量、氣溫和風速變化如圖2所示。從圖2中可看出,農(nóng)田站和西沙漠站儀器的感熱通量、氣溫和風速都吻合得較好,平均偏差分別為3.90W·m-2、0.20℃和0.092m·s-1。所以農(nóng)田站和西沙漠站儀器觀測的感熱、氣溫及風速資料可直接使用。農(nóng)田站和西沙漠站的感熱通量值、空氣密度和風速都使用3m處的觀測值,農(nóng)田站3m處氣溫可用2m處氣溫值通過邊界層位溫廓線公式計算得到,即:其中:θ為位溫;θ0為熱量粗糙度zoh高度處的位溫;θ*為特征位溫;κ為卡曼常數(shù)(取0.35);ψh為穩(wěn)定度修正函數(shù),由Busingeretal給出的普適函數(shù)計算;ζ=z/L,z=Z-d,Z為儀器高度,d為零平面位移,d=32h,h為植被高度,L為莫寧-奧布霍夫長度;熱量粗糙度zoh的計算在下面詳細介紹。各測站的實測地表溫度可根據(jù)地表輻射平衡計算得到,即:其中:Fu、Fd分別為向上和向下長波通量(由于Fd在各測站之間變化很小,故移動點取西沙漠的Fd值),單位為W·m-2;ε為地表發(fā)射率(農(nóng)田站和加強觀測的4個測站取0.99,東、西沙漠站取0.97);σ=5.67×10-8W·m-2·K-4,為Stefan-Boltzmann常數(shù);Ts_o為實際觀測長波輻射計算得到的地表溫度,單位為K。4分裂窗法和各自變量法地表溫度反演算法根據(jù)所選通道數(shù)的不同可分為單通道法、多通道法和分裂窗法,其中分裂窗法是目前發(fā)展最為成熟、運用最為普遍的方法。本文使用三種常見的分裂窗算法估算試驗區(qū)的地表溫度,通過反演結(jié)果的對比分析得到一種較為適合試驗區(qū)的地表溫度估算方法。4.1表溫度和比輻射率該算法的具體形式如下:其中:Ts_b為估算的地表溫度,單位為K;A0為常數(shù);T31、T32分別為MODIS第31、32通道的輻射亮溫;p、m是與比輻射率相關(guān)的計算參數(shù)。詳細算法參見文獻。(2)獨立分裂窗算法組成Kerretal在1992年將植被指數(shù)引入分裂窗算法,這種算法主要由兩個獨立的分裂窗算法組成:其中:Pv為植被覆蓋率;Tveg、Tsoil分別為植被和裸土溫度,Tveg=T31+2.6(T31-T32)-2.4,Tsoil=T31+2.1(T31-T32)-3.1;Ts_k為估算的地表溫度,單位為K。(3)地表溫度變化的計算覃志豪等2005年提出一種適用于MODIS數(shù)據(jù)的分裂窗算法,主要考慮像元為混合像元,即:其中:Ts_t為估算的地表溫度;A0、A1、A2分別為分裂窗計算的中間參數(shù),主要由地表發(fā)射率和大氣透過率計算得到。詳細計算方法參見文獻。4.2反演表面溫度的主要參數(shù)計算其中:r1、r2分別為MODIS第1、2通道的反射率;dNDVI為植被指數(shù)修正項,由衛(wèi)星天頂角θ計算得到:(2)植被覆蓋率Pv為植被覆蓋率,由歸一化植被指數(shù)NDVI計算得到:其中:NDVIS、NDVIV分別為裸土和完全植被覆蓋下的植被指數(shù)值。(3)大氣透過率估算大氣透過率主要受大氣中水汽含量的影響,在計算大氣透過率時,先計算大氣水汽含量,再利用大氣水汽含量與大氣透過率的關(guān)系來估算大氣透過率。而大氣的水汽含量w可由MODIS的第2和第19通道反演得到:其中:a、b為常數(shù);r2、r19分別為第2和第19通道的反射率。(4)b比輻射率測定不同的地表比輻射率計算對地表溫度的反演結(jié)果也有一定的影響。本文采取混合像元的通道比輻射率計算方法:其中:下標v和s分別代表植被和裸土;R為輻射比率;dε為熱輻射相互作用校正項。詳細計算方法參見文獻。5結(jié)果分析5.1地表溫度差的測定圖3a為加強觀測期間(2010年7月5日)各站的地表溫度變化。由于7月4日有降水,所以5日00:00-06:00左右農(nóng)田的地表溫度都低于沙漠的,而在晴天晚上,農(nóng)田的地表溫度高于沙漠的。白天隨著太陽的加熱作用,農(nóng)田和沙漠的地表溫度都迅速升高。由于植被的作用,農(nóng)田升溫到一定程度后,升溫幅度減小。在15:00左右達到峰值,沙漠最高溫度比農(nóng)田的高15.0℃左右。東1點位于東沙漠和農(nóng)田之間,地表為稀疏植被覆蓋,其地表溫度變化始終與沙漠的地表溫度變化相同。東2、東3和西1點都設(shè)在農(nóng)田里,有植被覆蓋,地表溫度變化都與農(nóng)田的地表溫度變化相同?？傮w來說,東1點表現(xiàn)為沙漠特征,東2點為農(nóng)田特征,介于東1點和東2點之間,即為沙漠和農(nóng)田過渡帶;西1點表現(xiàn)出明顯的農(nóng)田特征,介于西1點和西沙漠之間,也為農(nóng)田與沙漠的過渡帶。使用上述試驗區(qū)加強觀測期間7個觀測點的地表溫度,計算相鄰兩個觀測點之間的地表溫度差,得到試驗區(qū)的地表溫度差圖(圖3b)。東1點和東2點之間、西1點和西沙漠站之間的過渡帶地表溫差較大,有明顯的日變化特征;在中午時溫差達到最大,約為17.0℃。而農(nóng)田內(nèi)部觀測點(東2、東3、農(nóng)田、西1)之間的地表溫度差與沙漠內(nèi)部觀測點(東沙漠、東1)之間的地表溫度差都<7.0℃。5.2地表溫度算法修正使用上述三種地表溫度算法估算7個觀測站點的地表溫度,并與實際觀測的地表溫度進行對比。所有7個觀測站點的地表溫度可用于對比的樣本共30個。圖4a1、b1、c1分別是利用上述三種地表溫度估算法估算的地表溫度和觀測的對比結(jié)果?？梢钥吹?兩種地表溫度基本上對稱地分布在1∶1線兩側(cè),進一步分析表明,兩種地表溫度之間的相關(guān)系數(shù)最高達到0.92,均方根誤差都接近于3.0K。由于地表溫度算法中的參數(shù)都是與地表類型及所處的大氣條件有關(guān),所以本文利用實測地表溫度和反演地表溫度之差修正反演算法中與地表發(fā)射率及大氣條件相關(guān)的參數(shù),從而提高地表溫度反演精度。修正后的地表溫度算法如下:Becker算法修正后:Kerr算法修正后:覃志豪算法修正后:圖4a2、b2、c2給出了修正后地表溫度算法估算各站點的地表溫度結(jié)果。修正后的Becker算法相關(guān)系數(shù)最高達到0.93,均方根誤差最小為2.38K,平均相對誤差<4.5%。其他兩種算法修正后相關(guān)系數(shù)也達到0.92左右,均方根誤差都<2.6K,平均相對誤差都<4.5%。所以試驗區(qū)的地表溫度使用修正后的Becker算法估算?？傮w上看,對三種地表溫度算法的修正都能提高地表溫度的估算精度。5.3地表溫度分布選取試驗區(qū)及其周邊區(qū)域為研究對象,選用2010年7月5日衛(wèi)星過境時(12:10)的衛(wèi)星資料研究試驗區(qū)的地表溫度分布特征。根據(jù)觀測站實測地表溫度分析結(jié)果及研究區(qū)NDVI的分布(圖5a),以NDVI為標準區(qū)分綠洲和沙漠下墊面:NDVI<0.2為沙漠下墊面;NDVI≥0.4為綠洲下墊面,介于兩者之間為綠洲和沙漠之間的過渡帶。用修正后的Becker算法計算該區(qū)域的地表溫度。圖5b為2010年7月5日研究區(qū)域地表溫度的平面分布。從圖5b中可看出:(1)綠洲區(qū)的地表溫度主要在32.0~38.0℃之間,地表溫度梯度較小,即綠洲內(nèi)部的地表溫度等值線較稀疏,越往綠洲中心植被指數(shù)越大,地表溫度越低。(2)綠洲和沙漠的過渡帶,地表溫度主要為38.0~46.0℃,介于綠洲和沙漠之間。地表溫度梯度較大,即過渡帶內(nèi)的地表溫度等值線分布密集。(3)沙漠地區(qū)的地表溫度明顯高于過渡帶和綠洲區(qū),為46.0~52.0℃,溫度梯度較小,越往沙漠中心地表溫度梯度越小,等值線分布越稀疏。這種地表溫度的分布特征與5.1節(jié)中討論的不同下墊面地表溫度分布特征相對應。5.4估計ra的方法區(qū)域感熱通量可由下式計算得到:其中:ρ為空氣密度,單位為kg·m-3;Ts_b為遙感反演的地表溫度,單位為℃;Ta為3m處的空氣溫度,單位為℃;ra為空氣動力學阻抗,由Thometal1977年提出的估算ra的方法計算:考慮到在植被較高下墊面zoh和zom相差較大,將算法中的zoh和zom分開計算,公式(16)變?yōu)槠渲?Ux為觀測高度的平均風速,單位為m·s-1;zom為動量輸送的表面粗糙長度,單位為m,按Sellers的經(jīng)驗方法給定:其中:h為植被高度,單位為m;h過渡帶為過渡帶的植被高度,h過渡帶=NDVI/0.4h;zoh為熱量輸送的表面粗糙長度,單位為m,根據(jù)數(shù)量級粗略地估算,zoh=0.1zom;過渡帶的風速、氣溫取沙漠和綠洲風速、氣溫的平均值。5.5測點結(jié)果的誤差感熱能量和空氣動力學阻抗可用于對比的樣本共13個。從整體上看,改進前的空氣動力學阻抗估算結(jié)果總體偏小(圖6),平均偏小20.70%。改進后的空氣動力學阻抗法估算結(jié)果較好(圖7a),均方根誤差僅為17.98s·m-1。誤差較大的測站主要集中在農(nóng)田站,如圖7b所示。西沙漠站的空氣動力學阻抗估算結(jié)果都較好,相關(guān)系數(shù)為0.83,平均相對誤差為11.18%。利用估算的空氣動力學阻抗和估算的地表溫度計算得到的感熱通量如圖8a所示,相關(guān)系數(shù)為0.92,平均相對誤差為0.62%,均方根誤差為23.49W·m-2。估算的感熱通量與實際觀測的感熱通量相比,存在一定的偏差,但估算的感熱通量變化趨勢(農(nóng)田站:y=97.91+3.51x;西沙漠站:y=168.90+4.79x)與實測的感熱通量變化趨勢(農(nóng)田站:y=99.60+3.49x;西沙漠站:y=179.38+4.38x)是一致的(圖8b)。5.6西從t和ra之間的變化感熱通量估算的精度主要取決于地氣溫差ΔT和空氣動力學阻抗ra的估算精度。假定空氣密度為常數(shù)(取本試驗區(qū)的平均值0.97),感熱通量隨ra、ΔT的變化如圖9a所示。當ΔT一定時,ra超過一定的閾值,感熱隨ra變化很小;而在ra小于這個閾值時,感熱通量隨ra迅速減小;閾值隨ΔT的增大而增大。農(nóng)田站的ΔT主要集中在3.0~15.0℃之間(圖9b),ra主要在100~160s·m-1之間(見圖7a),在這個范圍內(nèi)感熱通量隨ra的變化小于隨ΔT的變化:(1)當ΔT一定時,ra從100s·m-1增加到150s·m-1時,感熱通量的增量不超過50W·m-2;(2)當ra一定時,ΔT從5.0℃增大到15.0℃時,感熱通量增加幅度在50~100W·m-2。西沙漠站的ΔT主要集中在10~30℃之間(圖9b),ra在60~100s·m-1之間(見圖7a),所以西沙漠站的感熱通量隨ΔT和ra變化都比較大:(1)當ΔT一定時,ra從70s·m-1變化到100s·m-1時,感熱通量的