基于天宮二號數(shù)據(jù)的北京市地表溫度反演研究

白沁靈孟 丹*邢 婧

(1.資源環(huán)境與地理信息系統(tǒng)北京市重點實驗室，北京100048；2.首都師范大學資源環(huán)境與旅游學院，北京100048)

1 引言

地表溫度體現(xiàn)地球表面物質(zhì)的能量變化與交換過程，是地球科學研究的一個重要指標。地表溫度的區(qū)域分布在氣候變化、植被生態(tài)、環(huán)境監(jiān)測和城市熱島等研究領(lǐng)域都有著重要的應(yīng)用價值[1]。

隨著遙感衛(wèi)星技術(shù)的發(fā)展，越來越多的衛(wèi)星搭載了波段和分辨率等參數(shù)不一的熱紅外傳感器，國內(nèi)外學者針對不同的傳感器提出多種地表溫度反演方法，這些算法大致可以分為5類：單通道算法、多通道算法、多角度算法、多時相算法和高光譜反演算法[2]。其中多通道算法用2個或2個以上的熱紅外波段反演地表溫度，減少了一些參數(shù)的假設(shè)，因此反演精度較高[3]。多通道算法中，針對2個熱紅外通道數(shù)據(jù)的劈窗算法計算簡單、方法成熟，是目前應(yīng)用最為廣泛的溫度反演算法之一。劈窗算法是由McMillin[4]提出的一種根據(jù)AVHRR數(shù)據(jù)的2個熱紅外波段來反演海面溫度的方法，后被推廣應(yīng)用于陸地表面溫度反演。該算法利用2個相鄰?fù)ǖ缹Υ髿馕兆饔玫牟煌?，通過2個通道輻射亮度的各種線性組合來剔除大氣的影響，反演地表溫度。近年來，國內(nèi)外學者基于 ASTER[5]、MODIS[6]、Landsat-TIRS[7]、Himawari 8 AHI[8]等數(shù)據(jù)開展了一系列的劈窗算法研究和應(yīng)用。

天宮二號寬波段數(shù)據(jù)具有2個相鄰的熱紅外波段(T1∶8.125 ～8.825 μm，T2 ∶8.925 ～9.275 μm)。當前劈窗算法多是針對10～14 μm波長區(qū)間的數(shù)據(jù)，很少涉及8～10 μm區(qū)間內(nèi)的熱紅外通道。孟鵬等[9]證明了將8～10 μm波長區(qū)間的數(shù)據(jù)用于劈窗算法反演地表溫度具有一定的可行性?；谝陨涎芯?，本文在傳統(tǒng)劈窗算法的基礎(chǔ)上，調(diào)整算法參數(shù)，對北京市的地表溫度進行反演，并利用MODIS地表溫度產(chǎn)品對反演結(jié)果進行了精度驗證，最后統(tǒng)計分析北京市不同土地利用類型的地表溫度情況。

2 研究區(qū)與數(shù)據(jù)

2.1 研究區(qū)概況

北京市地處華北平原，位于北緯39°36′到41°03′、東經(jīng) 115°25′到 117°30′之間，以山地和平原地形為主，地勢西北高、東南低。三面環(huán)山，東南部是一片緩緩向渤海傾斜的平原，海拔介于10～2157 m之間，總面積達 16 410.54 km2，屬于暖溫帶半濕潤大陸性季風氣候，夏季高溫多雨，冬季寒冷少雨，春、秋短促，年均降水量 595 mm，其中75%集中在夏季，年平均氣溫11～13℃。

2.2 遙感數(shù)據(jù)源

天宮二號空間實驗室搭載的寬波段成像儀是新一代寬波段、寬視場和圖譜合一的光學遙感器，國內(nèi)外首次在單臺儀器上實現(xiàn)了可見光、短波紅外和熱紅外多光譜大視場全推掃成像的組合集成功能，對地表覆蓋信息、地表溫度等研究具有重要意義[10]。

本研究使用天宮二號寬波段成像光譜儀的可見光近紅外圖像和熱紅外圖像對地表溫度進行反演，數(shù)據(jù)來源于載人航天空間應(yīng)用數(shù)據(jù)推廣平臺(http：//www.msadc.cn/)。本次研究選取2018年4月27日經(jīng)過系統(tǒng)輻射和幾何校正后的兩景天宮二號寬波段L2級產(chǎn)品數(shù)據(jù)，圖像質(zhì)量良好，無云。地表溫度反演結(jié)果驗證借助MODIS地表溫度數(shù)據(jù)產(chǎn)品 MOD11 A(https：//modis.gsfc.nasa.gov/)。土地利用數(shù)據(jù)來源于中國科學院資源環(huán)境科學數(shù)據(jù)中心(http：//www.resdc.cn/)。

3 研究方法

本文在傳統(tǒng)劈窗算法的基礎(chǔ)上，針對天宮二號熱紅外通道的特點，重新定義算法參數(shù)。

首先對輻射傳輸方程中的Planck函數(shù)進行線性簡化，建立2個熱紅外通道的亮度溫度與輻射亮度之間的關(guān)系；接著利用天宮二號數(shù)據(jù)的第2和第4波段反演大氣水汽數(shù)據(jù)，并計算其大氣透過率；然后根據(jù)ASTER光譜數(shù)據(jù)庫提供的地物波譜數(shù)據(jù)以及利用可見光數(shù)據(jù)計算的歸一化植被指數(shù)NDVI(Normalized Difference Vegetation Index)計算地表比輻射率；最后聯(lián)立2個熱紅外波段的輻射傳輸方程，即可求得地表溫度。技術(shù)流程圖如圖1所示。

圖1 技術(shù)流程圖Fig.1 Technology flowchart

3.1 圖像預(yù)處理

獲取的天宮二號數(shù)據(jù)是經(jīng)過系統(tǒng)輻射校正、幾何校正后的(入瞳處)輻亮度產(chǎn)品，DN(Digital Number)值是經(jīng)過放大的輻亮度，經(jīng)過縮小還原后得到天宮二號數(shù)據(jù)的輻射亮度。采用Landsat 8影像作為基準參考影像，對天宮二號影像數(shù)據(jù)進行幾何精糾正，校正精度控制在1個像元內(nèi)。采用FLAASH大氣校正模塊對天宮二號數(shù)據(jù)進行大氣校正，設(shè)置天宮二號的軌道參數(shù)，并導(dǎo)入其波譜響應(yīng)函數(shù)。最后利用行政邊界數(shù)據(jù)對影像進行裁切，得到天宮二號寬波段成像儀的北京市遙感影像(圖2)。

圖2 2018年4月27日北京市天宮二號影像圖(紅：V9，綠：V10，藍：V12)Fig.2 TG-2 image of Beijing on April 27， 2018(R：V9， G：V10， B：V12)

3.2 劈窗算法的推導(dǎo)

地表溫度反演通過建立輻射傳輸方程來實現(xiàn)，熱輻射傳輸方程是熱紅外遙感和地表溫度反演的基礎(chǔ)[5]，地表輻射經(jīng)過大氣到達傳感器的過程可以簡單的描述為式(1)：

式中Bi(Ti)是遙感器所接收的輻射亮度；Ts是地表溫度；τi(θ)是通道i在遙感器視角為θ下從地面到遙感器的大氣透過率；εi是比輻射率；Bi(Ts)是在地表溫度為Ts時的黑體輻射亮度；Ti是星上亮度溫度；是大氣下行輻射；是大氣上行輻射。覃志豪等[11]對有詳細的推導(dǎo)，有近似解見式(2)：

式中Ta為大氣向上平均溫度，為向下平均作用溫度。

熱輻射傳輸方程涉及Planck函數(shù)的計算，由于Planck函數(shù)十分復(fù)雜，所以在進行地表溫度反演研究時，通常對其進行線性簡化。覃志豪[12]研究AVHRR數(shù)據(jù)反演地表溫度和 Offer[7]研究Landsat 8數(shù)據(jù)反演地表溫度研究時，均對Planck函數(shù)進行泰勒展開，見式(3)：

將線性簡化的 Planck函數(shù)式(3)帶入式(1)，對于天宮二號的2個熱紅外波段輻射傳輸方程簡化為式(4)：

為了方便計算，定義系數(shù)為式(5)：

這樣，式(4)可簡寫為式(6)：

聯(lián)立2個方程組，求解出地表溫度為式(7)：

ai、aj、bi、bj為線性簡化的 Planck 函數(shù)中的系數(shù)。這樣，只需確定Planck函數(shù)參數(shù)、亮度溫度、大氣透過率、地表比輻射率即可得到地表溫度Ts。

3.3 Planck函數(shù)線性簡化

本文利用Planck函數(shù)分別對天宮二號寬波段數(shù)據(jù)的2個熱紅外波段的輻射亮度在溫度273～322 K區(qū)間[13]內(nèi)的變化關(guān)系進行計算，如圖3所示，輻射亮度隨溫度的變化近似線性關(guān)系。進行線性回歸得到方程式(8)。

圖3 天宮二號2個熱紅外波段的輻射亮度與溫度關(guān)系圖Fig.3 Relationship between radiation intensity and temperature in two thermal infrared bands of TG-2

3.4 亮度溫度的計算

天宮二號寬波段數(shù)據(jù)為經(jīng)輻射校正和系統(tǒng)幾何校正后的產(chǎn)品，利用Planck函數(shù)求解亮度溫度[14]，計算公式如式(9)所示：

式中Ii為輻射亮度，Ki，1＝2hc2/λ5i，Ki，2＝hc/(kλi)，h為普朗克常數(shù)，約為 6.62606896×10-34J˙s，c為光速，為 2.99792×108m˙s-2，k 為玻爾茲曼常數(shù)，為 1.3806505×10-23J˙K-1，λi為第i通道內(nèi)的中心波長。

天宮二號寬波段數(shù)據(jù)在熱紅外數(shù)據(jù)T1波段的中心波長λ1＝8.475 μm，在熱紅外數(shù)據(jù)T2波段的中心波長λ2＝9.1 μm，根據(jù)Planck公式得出天宮二號寬波段數(shù)據(jù)的Ki，1和Ki，2。 對于熱紅外波段T1，K1，1＝2724.1 W˙m-2˙sr-1˙μm-1，K1，2＝1697.3 K，熱紅外波段T2，K2，1＝1908.6 W˙m-2˙sr-1˙μm-1，K2，2＝1580.7 K。

3.5 大氣水汽含量與大氣透過率估算

大氣水汽含量與大氣透過率是地表溫度反演過程中的關(guān)鍵參數(shù)。在MODIS的大氣水汽反演中利用MODIS數(shù)據(jù)的第2波段(0.841～0.876 μm)和第 19 波段(0.915～0.965 μm)計算其大氣水汽含量[15]。天宮二號寬波段數(shù)據(jù)的第4波段(0.845～0.885 μm)和第 2波段(0.930～0.950 μm)的波長范圍與MODIS數(shù)據(jù)的第2波段和第19波段波長范圍相近，故天宮二號寬波段數(shù)據(jù)反演大氣水汽含量的估算公式如式(10)所示：

式中ω是大氣水汽含量(g/cm-2)，取α＝0.02和β＝0.6321[15]；ρ2和ρ4是天宮二號數(shù)據(jù)第2和第4波段的地面反射率。

張曉等[16]模擬出的8～9.3 μm波長范圍內(nèi)2個熱紅外波段的夏季大氣透過率與水汽含量的關(guān)系，利用該關(guān)系計算大氣透過率。計算公式如式(11)所示：

式中τ1、τ2為天宮數(shù)據(jù)2個熱紅外波段的大氣透過率。

3.6 比輻射率的估算

地表比輻射率主要由地物表面的結(jié)構(gòu)和波長范圍決定，天宮二號熱紅外數(shù)據(jù)的波段范圍為8～9.3 μm。在劈窗算法地表溫度反演中，通常假定地表輻射率已知，估算方法主要有經(jīng)驗值法和混合像元法。本文采用覃志豪等[17]提出的混合像元法對地表比輻射進行估計?；旌舷裨ㄖ?，將遙感影像分為水體、自然表面和城鎮(zhèn)表面，自然表面看作植被和裸土的混合，城鎮(zhèn)表面看作城鎮(zhèn)和植被的混合，其計算公式如式(12)所示：

式中εsurface為自然表面的比輻射率，εbuilding為城鎮(zhèn)表面的比輻射率，Pv為植被覆蓋度，εiv、εis、εim為植被、裸土、建筑在熱紅外波段內(nèi)的比輻射率，Rv、Rs、Rm為植被、裸土、建筑的溫度比率；dε為比輻射率的修正項，在地表相對平整的情況下，這一項的值就很小[18]，研究區(qū)域?qū)儆谄皆貐^(qū)，因此取 dε＝0；

混合像元法中提出了植被、裸土和建筑表面的溫度比率Rv、Rs、Rm的估算公式如式(13)所示：

其中Pv為植被覆蓋度，由公式Pv＝NDVINDVIs/NDVIv-NDVIs求出，NDVIv和NDVIs分別表示純植被像元和純裸土像元的歸一化植被指數(shù)NDVI值。當NDVI≤NDVIs時，植被覆蓋度Pv＝0；當NDVI≥NDVIv時，植被覆蓋度Pv＝1。 由于受到季節(jié)、地表濕度、大氣條件等因素的影響，NDVIv和NDVIs并不是一個固定值[19]，本文通過統(tǒng)計天宮二號數(shù)據(jù)影像NDVI直方圖，取5%和95%的累積百分比作為置信區(qū)間，取NDVI的上下限閾值分別為研究區(qū)的NDVIv和NDVIs來進行植被覆蓋度的近似估計。

根據(jù)ASTER光譜數(shù)據(jù)庫提供的典型地物波譜數(shù)據(jù)計算典型地物的比輻射率，水體、植被、裸地和建筑在8～9.3 μm波段范圍內(nèi)的反射率變化如圖4所示。由于該波譜庫的數(shù)據(jù)為反射率數(shù)據(jù)，根據(jù)基爾霍夫定律，處于某一溫度的物質(zhì)比輻射率與吸收率值相同。除水體外，對于一般地物都可不考慮透射，即比輻射率＝吸收率＝1-反射率，水體的比輻射率在熱紅外波段內(nèi)變化較小，取0.995[17]。選擇幾種典型的地物計算其平均比輻射率，典型地物水體、植被、裸地和建筑在天宮二號數(shù)據(jù)的2個熱紅外波段T1和T2范圍內(nèi)的比輻射率分別為：0.995、0.9802、0.9457、0.9378 和0.995、0.9761、0.9432、0.9408。

綜上，自然表面和城鎮(zhèn)表面在天宮二號的熱紅外波段范圍內(nèi)的比輻射率由公式(14)進行估算。

4 結(jié)果分析

4.1 地表溫度反演結(jié)果驗證

圖4 典型地物在8～9.4 μm波段范圍內(nèi)的反射率Fig.4 Reflectivity of typical ground objects in the range of 8～ 9.4 μm

采用上述劈窗算法，得到天宮二號地表溫度反演結(jié)果(圖5)，選取2018年4月27日的MODIS地表溫度產(chǎn)品(圖6)與天宮二號寬波段數(shù)據(jù)反演結(jié)果進行一致性分析。天宮二號熱紅外波段的分辨率為400 m，而MODIS地表溫度產(chǎn)品空間分辨率為1 km，將MODIS溫度產(chǎn)品空間分辨率重采樣為400 m，進而與同分辨率的反演結(jié)果進行對比。從圖5、圖6可以看出，天宮二號寬波段數(shù)據(jù)反演得出的北京市地表溫度最低溫度為292.2 K，低溫部分主要出現(xiàn)在密云水庫，最高溫度為310.4 K，高溫區(qū)出現(xiàn)在中心城區(qū)和延慶地區(qū)；MODIS溫度產(chǎn)品的最低溫度為286.8 K，最高溫度為309.0 K，高溫區(qū)和低溫區(qū)與天宮二號寬波段數(shù)據(jù)反演的結(jié)果一致。從空間上看，天宮二號寬波段數(shù)據(jù)反演得出的地表溫度結(jié)果與MODIS地表溫度產(chǎn)品的空間分布情況基本一致，北京地區(qū)的地表溫度空間分布趨勢表現(xiàn)出明顯的城市熱島效應(yīng)，城鎮(zhèn)溫度高于自然地表的溫度，水體的地表溫度最低，從圖5可以看出密云水庫明顯邊界。

圖5 基于天宮二號數(shù)據(jù)反演的北京市地表溫度Fig.5 Beijing LST retrieval from TG-2 data

圖6 基于MODIS數(shù)據(jù)反演的北京市地表溫度Fig.6 Beijing LST retrieval from MODIS data

圖7 西南-東北方向地表溫度剖面圖Fig.7 LST profile in SW-NE direction

選取研究區(qū)地表溫度起伏較多的西南至東北方向?qū)μ鞂m二號和MODIS反演的地表溫度做剖面分析，以剖面線的西南端起點為原點，剖面線的位置見圖5、圖6，分析結(jié)果如圖7所示。剖面線上，天宮二號數(shù)據(jù)反演的地表溫度整體比MODIS地表溫度產(chǎn)品數(shù)據(jù)高；地表溫度剖面起伏情況基本一致；因為天宮數(shù)據(jù)的原始分辨率較高，其地表溫度剖面的起伏較多；高溫中心均出現(xiàn)在市區(qū)和城鎮(zhèn)地區(qū)，從市區(qū)和城鎮(zhèn)地區(qū)到郊區(qū)地表溫度逐漸降低；人類活動頻繁地帶溫度較高，自然表面的溫度較低，水體的溫度最低。

從研究區(qū)隨機選取1000個采樣點，將天宮二號寬波段數(shù)據(jù)反演的地表溫度結(jié)果與MODIS地表溫度產(chǎn)品進行對比分析，二者的誤差均值為3.5 K，誤差的標準差為1.4 K；對其進行線性回歸分析，可以看出，二者有著較高的相關(guān)性，其相關(guān)系數(shù)達到0.79，均方根誤差為1.28 K(圖8)。由于天宮二號數(shù)據(jù)與MODIS數(shù)據(jù)過境時間、空間分別率等存在差別，導(dǎo)致兩者反演的地表溫度存在差異，具體原因還有待下一步工作的研究。

圖8 基于MODIS與天宮二號數(shù)據(jù)反演的地表溫度散點圖Fig.8 LST scatter plot based on MODIS and TG-2 data inversion

4.2 北京市不同土地利用類型地表溫度分析

結(jié)合北京市2018年的土地利用圖(圖9)，統(tǒng)計各土地利用類型地表溫度的均值、中位數(shù)、最大值、最小值、標準差，如表1所示。

不同土地的利用類型存在差異。其中，平均地表溫度較高的是建設(shè)用地，為305.40 K，這是由于地表溫度的高低主要取決于地面物體的熱學性質(zhì)，如熱容量、熱慣量、熱傳導(dǎo)率以及地表潛熱通量[20]。城市人工構(gòu)筑物吸熱快而熱容量小，升溫快，同時，城市人為的熱排放也是造成建設(shè)用地的地表溫度較高的原因之一。水域的平均地表溫度最低為301.1 K，原因是水體的熱容量大，升溫慢，且受熱后會以水分蒸發(fā)的形式把熱量散發(fā)出來，降低自身的溫度。耕地的平均地表溫度僅次于建設(shè)用地，原因是其植被覆蓋度較低，植被蒸騰作用弱。不同土地利用類型的地表溫度的平均值從高到低的順序為：建設(shè)用地＞耕地＞草地＞林地＞未利用土地＞水域。

利用不同土地利用類型所對應(yīng)的地表溫度繪制箱線圖，如圖10所示。箱線圖可以看出反演的溫度值是否穩(wěn)定，而中位數(shù)可以對比不同類型數(shù)據(jù)的趨勢。從箱線圖可以看出：水域地表溫度數(shù)據(jù)值的箱體較長，說明該地表溫度數(shù)據(jù)變化較大，較不穩(wěn)定，可能是因為天宮二號寬波段數(shù)據(jù)的熱紅外波段的400 m的分辨率容易形成混合像元，且由于河流兩岸多為耕地、草地、建設(shè)用地等，吸收更多的太陽輻射，混合像元會表現(xiàn)出更高的表面溫度，因此會導(dǎo)致河流表面溫度過高；而其他土地利用類型的箱體相對較短，說明其地表溫度較為穩(wěn)定。

圖9 2018年北京市土地利用圖Fig.9 Land use map of Beijing in 2018

表1 不同土地利用類型地表溫度情況統(tǒng)計(單位：K)Table 1 Statistics of LST for different Land use(Unit：K)

5 結(jié)論

1)天宮二號寬波段數(shù)據(jù)波段設(shè)置可以滿足地表溫度反演的參數(shù)需求。其中的第2波段和第4波段，可以反演出大氣水汽含量，進而計算出大氣透過率。

圖10 北京市不同土地利用類型地表溫度箱線圖Fig.10 LST Box Map of different Land use types in Beijing

2)基于天宮二號數(shù)據(jù)反演出的結(jié)果與MODIS地表溫度產(chǎn)品具有較高的一致性，該參數(shù)可推廣應(yīng)用于熱紅外波段在8～10 μm波長的傳感器。驗證了天宮二號寬波段數(shù)據(jù)可以利用改進參數(shù)后的劈窗算法進行地表溫度反演，天宮二號為地表溫度反演提供了新的數(shù)據(jù)源。

3)北京地區(qū)不同土地利用類型地表溫度不同，地表溫度情況從高到低的順序為：建設(shè)用地＞耕地＞草地＞林地＞未利用土地＞水域。

致謝：感謝載人航天工程提供天宮二號寬波段成像儀數(shù)據(jù)產(chǎn)品。