張琪，李大成,*，歷華，段四波，張霄羽，韓啟金，高海霞

1. 太原理工大學(xué)，太原 030024 2. 中國科學(xué)院 空天信息創(chuàng)新研究院，北京 100094 3. 中國農(nóng)業(yè)科學(xué)研究院，北京 100081 4. 山西大學(xué)，太原 030006 5. 中國資源衛(wèi)星應(yīng)用中心，北京 100094

1 引言

地表溫度(land surface temperature，LST)作為衡量地球表面能量平衡狀態(tài)的重要參數(shù)之一,對地氣循環(huán)系統(tǒng)研究、氣候預(yù)報、生態(tài)環(huán)境、林業(yè)養(yǎng)殖以及地球化學(xué)等領(lǐng)域有深遠意義。因此，對地表溫度的空間和時間上的精確測量對于氣候變化、生態(tài)環(huán)境監(jiān)測等各個方面的應(yīng)用都有重要意義。隨著遙感技術(shù)的快速發(fā)展，熱紅外遙感為高精度地表溫度的獲取提供了有效的途徑[1]。

為了獲得具有較高時空分辨率的地表溫度，國內(nèi)外學(xué)者針對不同的遙感數(shù)據(jù)提出了多種地表溫度反演算法[2-4]。單通道算法是針對一個熱紅外通道提出的地表溫度反演算法[5-6]，輻射傳輸方程法[7]是常用的一種算法。該算法通過輸入大氣廓線數(shù)據(jù)并計算大氣參數(shù)來計算表面溫度。由于實時的大氣參數(shù)一般是利用大氣輻射傳輸模型進行模擬得到的，獲取困難，因此適用性不高。為了解決這一問題，文獻[8]對輻射傳輸方程進行簡化，提出了一種針對Landsat-TM5數(shù)據(jù)的單窗算法(MW算法)，該算法不需要大氣廓線數(shù)據(jù)，只需要根據(jù)大氣透過率、大氣水分含量和平均大氣溫度的經(jīng)驗線性關(guān)系對大氣進行校正，就能計算得到地表溫度。文獻[9]提出了一種只需要輸入大氣水汽含量就可以計算地表溫度的普適性單通道算法(JM&S算法)，該方法通過分析大氣水汽含量與大氣透過率、大氣上行輻射和大氣下行輻射之間的經(jīng)驗關(guān)系來消除大氣影響，計算簡便。

中國于2008年發(fā)射了環(huán)境一號A/B衛(wèi)星(HJ-1A/B)，其中環(huán)境一號B星(HJ-1B)搭載兩臺寬幅多光譜(CCD)相機和一臺紅外相機(IRS)，重訪周期為4 d。CCD相機可提供空間分辨率為30 m的4個波段，IRS相機可提供空間分辨率為300 m的熱紅外波段[10-11]。自衛(wèi)星發(fā)射以來，國內(nèi)學(xué)者針對環(huán)境小衛(wèi)星的單個熱紅外波段的數(shù)據(jù)特點展開了一系列研究工作[12-18]。文獻[12]針對環(huán)境一號衛(wèi)星B星進行數(shù)據(jù)計算，對MW算法和JM&S算法進行對比，結(jié)果表明兩種算法總體接近，JM&S算法精度更高。文獻[13]使用3種單通道算法對HJ-1B數(shù)據(jù)進行地表溫度反演，并對精度進行分析比較，結(jié)果表明JM&S算法精度最高。

在獲得衛(wèi)星影像地表溫度反演結(jié)果后，需要對其真實性進行驗證[19]。對于LST產(chǎn)品的真實性檢驗有4種方法：基于地面的驗證、基于衛(wèi)星LST產(chǎn)品的驗證、基于輻射的驗證和基于時間序列的比較?；诘孛娴尿炞C即將反演結(jié)果與地面測量的溫度進行直接比較，簡單易行，常用于驗證MODIS的地表溫度產(chǎn)品[20]。基于衛(wèi)星LST產(chǎn)品的驗證方法是將兩種衛(wèi)星得到的LST產(chǎn)品進行比較[21]。該方法適用于面積較大的區(qū)域，但是會受到衛(wèi)星傳感器的過境時間、觀測視角以及不同衛(wèi)星數(shù)據(jù)空間分辨率等因素的影響[22]?；谳椛涞尿炞C不需要原位測量，而是需要準(zhǔn)確估算通道的比表面發(fā)射率值以及衛(wèi)星過境時間的大氣水分含量和大氣溫度[23]?；跁r間序列的比較方法通過選取較長時段的遙感地表溫度,分析其自身表現(xiàn)的時序特征來發(fā)現(xiàn)傳感器在運行中可能存在的問題，對地表溫度數(shù)據(jù)進行異常值處理[24]。該方法僅僅能夠得到相對精度,對于真實性的評價存在明顯不足。

隨著中國遙感技術(shù)的發(fā)展，針對國產(chǎn)衛(wèi)星的數(shù)據(jù)反演及真實性檢驗工作也不斷深入。大多數(shù)研究對于反演結(jié)果的驗證較少，通常只是利用同一時刻的衛(wèi)星LST產(chǎn)品進行驗證，缺少對反演結(jié)果的定量評價，無法精確應(yīng)用于其他研究。因此，在HJ-1B衛(wèi)星方面地表溫度的反演及產(chǎn)品的真實性檢驗尚需要完善。本文以黑河流域多個類型的下墊面為研究區(qū)，采用普適性單通道算法對覆蓋研究區(qū)的2012-2018年清晰無云的HJ-1B遙感影像進行地表溫度反演，并利用對應(yīng)時刻的MOD11A1地表溫度和地面站點數(shù)據(jù)對反演結(jié)果進行驗證。由于驗證數(shù)據(jù)在時序上比反演數(shù)據(jù)多，并不完全匹配，為了充分利用驗證數(shù)據(jù)，引入時序相似性檢測思想，提出利用動態(tài)時間規(guī)整(dynamic time warping，DTW)將時間序列進行拉伸或收縮匹配來計算兩個時序之間的相似性。

2 研究區(qū)及數(shù)據(jù)預(yù)處理

研究區(qū)位于甘肅省黑河流域中上游，東經(jīng)97.1°～102.0°，北緯37.7°～42.7°，該地區(qū)氣候干旱，降水較少。本研究區(qū)內(nèi)用于驗證的地面站點有阿柔超級站(ArouCJZ)、黃藏寺站(HZS)、埡口站(YK)、花寨子荒漠站(HZZ)、神沙窩沙漠站(SSW)、張掖濕地站(SD)和黑河遙感站(HHYG)，站點的地表覆蓋類型主要包括高寒草地、戈壁、濕地、荒漠、沙漠、農(nóng)田等，具體情況如表1所示。

表1 研究區(qū)站點情況

本文所選用的遙感數(shù)據(jù)是HJ-1B數(shù)據(jù)與MODIS數(shù)據(jù)，其中HJ-1B數(shù)據(jù)包括多光譜數(shù)據(jù)和熱紅外數(shù)據(jù)，時間2012-2018年，共60個時相，均可覆蓋研究區(qū)并清晰無云。MODIS數(shù)據(jù)包括用于對反演結(jié)果進行驗證的MOD11A1地表溫度產(chǎn)品、用于計算大氣水汽含量的MOD021KM以及用于對MOD021KM進行幾何校正的MOD03數(shù)據(jù)。

在進行地表溫度反演之前首先對HJ-1B數(shù)據(jù)進行幾何校正，以消除影像因傳感器的高度、姿態(tài)等因素影響造成的幾何畸變，并以Landsat-8遙感影像對其進行影像配準(zhǔn)。隨后利用定標(biāo)系數(shù)對HJ-1B衛(wèi)星的CCD影像進行輻射定標(biāo)，定標(biāo)公式如下：

式中：L為定標(biāo)后的輻亮度，單位為W·m-2·sr-1·μm-1；a為絕對定標(biāo)系數(shù)增益；L0為絕對定標(biāo)系數(shù)偏移量，具體值可以在影像的頭文件中獲取。

衛(wèi)星傳感器在成像過程中會受到大氣、光照等因素的影響，因此需要對其進行大氣校正，本文采用基于MORTRAN4的FLAASH(fast line-of-sight atmospheric analysis of spectral hypercubes)大氣校正模型對定標(biāo)后的HJ-1B衛(wèi)星的CCD數(shù)據(jù)進行大氣校正。

為了進行下一步計算及驗證，對于MOD021KM數(shù)據(jù)通過MOD03數(shù)據(jù)進行幾何校正，隨后將MOD021KM和MOD11A1數(shù)據(jù)的投影坐標(biāo)轉(zhuǎn)換成研究區(qū)所在的投影，并將像元尺度重采樣到30 m。

本文所研究的地面站點數(shù)據(jù)來源于黑河流域遙感觀測自動氣象站，站點每30 s觀測一次，每10 min輸出一次結(jié)果，該結(jié)果為這10 min內(nèi)的平均值。這些站點的觀測設(shè)備不完全相同，主要有CNR1四分量輻射儀和SI-111紅外輻射計兩種觀測設(shè)備，兩種設(shè)備的計算公式分別為：

(1)

B(TS)=[B(Tr)-(1-ε)Lsky]/ε

(2)

式中：TS為地表溫度；B為普朗克公式；F↑為上行短波輻射強度；F↓為下行長波輻射強度；σ為斯蒂芬玻爾茲曼常數(shù)，σ=5.67×10-8W/(m2·K4)；εb為寬波段發(fā)射率，根據(jù)ASTER GED的窄波段發(fā)射率或波普庫數(shù)據(jù)可以計算得到；Tr為SI-111紅外輻射計觀測得到的站點實測溫度；Lsky為大氣下行輻射強度；ε為SI-111紅外輻射計實測的發(fā)射率數(shù)據(jù)。

3 研究方法

本文采用普適性單通道算法(JM&S算法)對HJ-1B衛(wèi)星數(shù)據(jù)進行地表溫度反演，并利用黑河流域遙感觀測站實測數(shù)據(jù)和MODIS地表溫度產(chǎn)品對反演結(jié)果進行驗證，基于動態(tài)時間規(guī)整(dynamic time warping，DTW)對反演結(jié)果和驗證數(shù)據(jù)的時序相似性進行評價，具體流程如圖1所示。

圖1 研究技術(shù)路線Fig.1 Research technology roadmap

3.1 基于JM&S的HJ-1B LST產(chǎn)品反演

普適性單通道算法(generalized single-channel method，SC)又稱JM&S算法，是由Jiménez-Muňoz和Sobrino在2003年針對只有一個熱紅外通道的遙感數(shù)據(jù)提出的地表溫度反演算法，該算法是通過對普朗克函數(shù)做一階泰勒級數(shù)展開所得到的，具體計算公式如下：

TS=γ[ε-1(ψ1Lsensor+ψ2)+ψ3]+δ

δ=-γLsensor+Tsensor

式中：Lsensor為衛(wèi)星傳感器測得的輻射強度，Tsensor為亮溫，單位為K；λ為有效波長；C1和C2為普朗克函數(shù)常量，C1=1.191 04×108W·μm4·m-2·sr-1，C2=1.438 77×104μm·K；大氣函數(shù)ψ1、ψ2和ψ3是關(guān)于大氣水分含量ω的函數(shù)。對于大氣函數(shù)的取值，文獻[12]根據(jù)HJ-1B衛(wèi)星的特征對函數(shù)進行重新擬合，得到的計算公式如下：

ψ1=0.041 2ω2+0.093 6ω+0.985 6

ψ2=-0.717 4ω2-0.881 2ω+0.394 1

ψ3=0.263 9ω2+0.649 9ω+0.470 3

地物比輻射率ε和大氣水分含量ω是JM&S算法的關(guān)鍵參數(shù)。大氣水分含量反映了大氣水分在空間上的不同分布，在地表溫度反演中有重要作用。考慮到數(shù)據(jù)獲取的便捷，MODIS數(shù)據(jù)的通道比值法在計算大氣水分含量方面應(yīng)用廣泛。因此本文采用與HJ-1B衛(wèi)星同時相的MODIS數(shù)據(jù)的第2波段(大氣窗口波段)和第19波段(水汽強烈吸收波段)計算大氣水分含量，計算公式如下：

式中：ρ2和ρ19分別為MODIS第2波段和第19波段的表觀反射率；α和β為常量，其值依據(jù)地表不同而有所不同。本研究區(qū)內(nèi)地表覆蓋類型包含植被、裸土及其他混合地表，因此本文采用混合地表類型，對于混合型地表，α=0.02，β=0.651。

按照地物類型的區(qū)別，影像分類一般包括植被、水體、裸土及建筑物表面這4類。在估算地物比輻射率方面，常用NDVI閾值法，該方法對地物在多光譜影像不同波段的光譜信息結(jié)合，通過NDVI的大小來區(qū)分不同地物。針對環(huán)境衛(wèi)星，根據(jù)地物的波譜特征，在裸土地區(qū)，地物比輻射率取值0.972，純植被地區(qū)的地物比輻射率取為0.99，由于水體接近于黑體，水體的比輻射率取為0.995[25]。混合像元區(qū)域地物類型復(fù)雜，可由如下公式計算地物比輻射率：

ε=εvPv+εs(1-Pv)+dε

dε=(1-εs)(1-Pv)Fεv

式中：Pv為植被覆蓋度，dε為地表幾何像元內(nèi)部的影響；NDVI為歸一化植被指數(shù)。對于環(huán)境衛(wèi)星，ρ3和ρ4分別為環(huán)境一號衛(wèi)星第3和第4波段的反射率；εv為完全植被覆蓋時的比輻射率；εs為完全裸土覆蓋時的比輻射率；F為地形因子，F(xiàn)=0.55。

3.2 基于DTW的星地時序非對應(yīng)LST產(chǎn)品檢驗

現(xiàn)有的對地表溫度反演結(jié)果的驗證都是基于對應(yīng)實點的驗證，無論是直接驗證還是間接驗證都無法對時序進行評定。用于驗證的數(shù)據(jù)往往在時序上連續(xù)且比反演數(shù)據(jù)的時序更多，當(dāng)驗證數(shù)據(jù)與反演數(shù)據(jù)時間不匹配時，為了充分利用地面站點數(shù)據(jù)，引入時序相似性思想，通過DTW對兩種數(shù)據(jù)的時序相似性進行質(zhì)量評價。

DTW通過將時間序列拉伸或收縮匹配來計算兩個時間序列之間的相似性，該方法相較于傳統(tǒng)歐式距離測度檢驗方法，克服了時間扭曲問題，在兩個時間序列之間尋找一個靈活的時間匹配路徑使兩個序列對齊，具有更加廣泛的實用性。

計算DTW首先需要構(gòu)建時間序列。分別構(gòu)建反演結(jié)果時序{H}={H1,H2,…,Hm}和檢驗數(shù)據(jù)時序{L}={L1,L2,…,Ln}，m和n分別為兩個時序的長度。其次，計算兩組時序之間任意兩個元素之間的差值，生成距離矩陣D(Li,Hj)=(Li-Hj)2，其中1≤i，j≤m。然后根據(jù)距離矩陣尋找最短路徑，即最佳規(guī)整路徑。該過程是基于動態(tài)規(guī)劃的思想，假如當(dāng)前節(jié)點是D(Li,Hj)，那么下一個節(jié)點必須是在D(Li+1,Hj)，D(Li,Hj+1)，D(Li+1,Hj+1)三者之間選擇，并且路徑必須是最短的。最后，通過如下遞歸式計算DTW累計距離：

DTW(i,j)=D(Li,Hj)+

4 結(jié)果與分析

本研究選取覆蓋黑河流域研究區(qū)的2012-2018年清晰無云的HJ-1B衛(wèi)星影像，采用JM&S的普適性單通道算法對其熱紅外波段進行地表溫度反演，并采用基于地面和基于衛(wèi)星地表溫度產(chǎn)品的兩種方法對結(jié)果進行精度驗證和分析。分別統(tǒng)計各個站點處的實測溫度與MODIS地表溫度，將反演結(jié)果與這兩者進行對比，計算平均偏差Bias、標(biāo)準(zhǔn)差STD、均方根誤差RMSE、動態(tài)時間規(guī)整DTW和決定系數(shù)R2。

圖2為反演溫度在7個站點處與地面溫度及站點處MODIS地表溫度的散點圖，圖2(a)(c)(e)(g)(i)(k)(m)分別為7個站點的反演溫度和地面溫度的相關(guān)性散點圖，圖2(b)(d)(f)(h)(j)(l)(n)分別為7個站點的反演溫度和MOD11A1地表溫度的相關(guān)性散點圖，圖上決定系數(shù)R2是在對反演溫度和驗證溫度做線性回歸模擬時計算得到的評價系數(shù)，無單位。圖3為7個站點的反演溫度與2種驗證溫度之間偏差值的折線圖。表2為兩種驗證方式在各個站點處的評價指標(biāo)值。

對于下墊面為沙漠、荒漠的神沙窩沙漠站和花寨子荒漠站，反演精度最高。圖3(a)(b)分別為花寨子荒漠站和神沙窩沙漠站的偏差值曲線，可以看出這兩個站點的反演溫度和地面溫度之間的差值較小，與MODIS地表溫度的偏差值稍大些但比其他站點小。從時序上來看沒有明顯受到時間的影響，可能是由于兩個站點的下墊面為荒漠和沙漠，空間異質(zhì)性較小，地表輻射率隨時間變化較小，因此反演精度較高。

對于下墊面為濕地、雪地的張掖濕地站和啞口站，反演精度相對其他站點較低。埡口站的Bias相對其他站點較高，與地面站點相比平均低估3.33 K，與MODIS地表溫度相比平均低估4.13 K。張掖濕地站的反演結(jié)果與地面站點相比平均低估2.24 K，與MODIS地表溫度相比平均低估3.12 K。但這兩個站點與驗證數(shù)據(jù)的決定系數(shù)R2均在0.9以上，表明反演結(jié)果與驗證數(shù)據(jù)的相關(guān)性并不差。造成這一結(jié)果的原因可能是對于濕地和雪地區(qū)域，地表覆蓋類型復(fù)雜，NDVI閾值法不能準(zhǔn)確地計算出該區(qū)域的地物比輻射率。

對于下墊面為植被或農(nóng)作物的3個站點，從圖3(c)(d)(e)可看出，阿柔超級站在7～10月反演結(jié)果與驗證數(shù)據(jù)差別較大，可能是因為阿柔超級站地表覆蓋類型為高寒草地，在夏季溫度上升植被生長，對傳感器的觀測造成影響。黃藏寺站點在3～5月的反演結(jié)果與驗證數(shù)據(jù)相差較大，在其他日期相差較小。這可能是因為該站點的地表覆蓋類型為小麥，冬小麥的生長季節(jié)在春季，成長季地表輻射率變化較大，因此反演結(jié)果的精度相對較低。黑河遙感站在夏季的反演結(jié)果與驗證數(shù)據(jù)相差比其他季節(jié)大，該站點的地表覆蓋類型為玉米，生長季節(jié)在6～10月，植物生長對地物比輻射率的計算造成影響，因此在植被生長季時反演效果不如其他季節(jié)好。

對于7個不同下墊面的站點研究區(qū)，反演得到的LST與溫度的偏差值均小于其與MODIS地表溫度的偏差值，RMSE和STD的表現(xiàn)與Bias一致。圖3的曲線圖也可看出偏差值大多為負值，表明反演結(jié)果大多低于驗證數(shù)據(jù)?；诘孛鏈囟鹊尿炞C顯示，除了埡口站的平均偏差最大為-3.33 K外，其他6個站點的與實測溫度的平均偏差均在2 K以內(nèi)。反演結(jié)果與地面溫度的RMSE最大為埡口站的3.37 K，除埡口站之外其他站點的RMSE最大都不超過3 K。所有站點的反演結(jié)果與地面溫度的STD都在2 K以內(nèi)。

與MODIS地表溫度相比，除了埡口站和張掖濕地站的平均偏差稍高外，其他站點的偏差值都在3 K以內(nèi)，RMSE指數(shù)最小可達1.26 K(花寨子荒漠站)。STD最大為阿柔超級站的2.04 K，最小為花寨子荒漠站的0.67 K。反演結(jié)果與MODIS地表溫度的差異可能是由于HJ-1B衛(wèi)星在本研究區(qū)的過境時間并不固定，而MODIS衛(wèi)星傳感器的一景影像的覆蓋區(qū)域較大，兩種衛(wèi)星的過境時間在不同日期相差不同，因此HJ-1B地表溫度的反演結(jié)果與MODIS的地表溫度數(shù)據(jù)差別并不穩(wěn)定。另外，HJ-1B CCD衛(wèi)星的原始空間分辨率為300 m，MODIS地表溫度產(chǎn)品的原始分辨率為1 000 m，由于最終使用的是重采樣得到的結(jié)果，因此插值計算也可能對結(jié)果造成影響。

反演結(jié)果與地面溫度相比，花寨子荒漠站和神沙窩沙漠站的DTW最小，為0.001 7，表明這兩個站點的反演溫度與地面溫度的相似性最高，這與Bias、RMSE和STD表現(xiàn)一致。從表2其他站點的驗證指標(biāo)也可看出，當(dāng)Bias、RMSE和STD較小時，DTW時序檢驗結(jié)果也相對較小，表明反演溫度與驗證溫度的相似性較高；反之當(dāng)Bias、RMSE和STD較大時，DTW時序檢驗結(jié)果也相對較大。在評價時序預(yù)測精度上，DTW與已有的評價指標(biāo)(如Bias、RMSE和STD)具有較高的一致性。

圖2 HJ-1B反演溫度與地面站點、MODIS溫度產(chǎn)品散點圖Fig.2 Scatter plots of HJ-1B retrieval LST and ground temperature,and of HJ-1B retrieval LST and MODIS LST

續(xù)圖2Fig.2 Continued

續(xù)圖2Fig.2 Continued

圖3 各站點反演溫度與驗證溫度偏差值Fig.3 The bias between retrieved temperature and verified temperature at each site

續(xù)圖3 Fig.3 Continued

表2 站點處反演溫度評價指標(biāo)

5 結(jié)論

本文使用Jiménez-Muňoz和Sobrino在2003年提出的普適性單通道算法對覆蓋黑河流域研究區(qū)的2012-2018年的HJ-1B衛(wèi)星數(shù)據(jù)進行地表溫度反演，并基于地面實測站點數(shù)據(jù)和MODIS地表溫度進行了精度分析和驗證，利用DTW對時序不匹配的數(shù)據(jù)進行了時序相似性評價。通過對結(jié)果的分析和討論，可以得出以下結(jié)論:

1)在本研究區(qū)內(nèi)，對于地表覆蓋類型為沙漠、荒漠的站點，由于地表輻射率變化不大，反演結(jié)果與地面站點相差最小(最小可達0.015 K)，精度最高；對于地表覆蓋類型為植被或農(nóng)作物的區(qū)域，地表輻射率隨時間發(fā)生變化，反演結(jié)果精度較好，平均偏差均在2 K以內(nèi)，相關(guān)性較高；對于地表覆蓋類型復(fù)雜的區(qū)域，可能是由于地物比輻射率不能準(zhǔn)確計算，反演結(jié)果與站點數(shù)據(jù)偏差較大。

2)HJ-1B反演結(jié)果與MODIS地表溫度產(chǎn)品對比時，對于所有站點，可能是由于環(huán)境衛(wèi)星和MODIS的過境時間不同，反演結(jié)果與驗證數(shù)據(jù)的差值均大于其與地面觀測溫度的差值，但相差不足1 K，可以認為兩種驗證方法結(jié)果一致。

3)從反演結(jié)果來看，在本研究區(qū)內(nèi)，下墊面為植被和農(nóng)作物的站點在植被生長季時反演精度較低，在其他季節(jié)反演精度較高。因此在對植被覆蓋區(qū)域進行地表溫度反演預(yù)測時，應(yīng)該選擇植被非生長季時相的影像。

4)對于時序不匹配的驗證數(shù)據(jù)和反演數(shù)據(jù)，利用DTW對其進行動態(tài)規(guī)整以評價兩組數(shù)據(jù)的時序相似性，DTW的評價結(jié)果與現(xiàn)有的評價指標(biāo)表現(xiàn)一致，為遙感影像的時序檢驗提供一種新的途徑。