臧金龍，國巧真，吳歡歡，喬悅，付盈

(天津城建大學(xué) 地質(zhì)與測繪學(xué)院，天津 300384)

0 引言

衛(wèi)星遙感在地表溫度(land surface temperature，LST)反演方面展示出巨大潛力，利用熱紅外波段數(shù)據(jù)進(jìn)行地表溫度反演，可快速且方便地獲得大范圍、較準(zhǔn)確、連續(xù)變化的地表溫度[1]。然而，城市環(huán)境的復(fù)雜性要求盡可能高的空間分辨率LST來進(jìn)行城市不同地表覆蓋類型熱特性的準(zhǔn)確描述[2]。由于現(xiàn)有星載熱紅外傳感器技術(shù)的限制，這些星載熱紅外影像的空間分辨率普遍偏低，限制了高空間分辨率的地表溫度信息在城市熱環(huán)境研究、城市能效監(jiān)測評估、建筑熱能消耗檢測等研究中的應(yīng)用[3-4]。為獲取滿足應(yīng)用需求的LST產(chǎn)品，利用高空間分辨率影像的波段信息進(jìn)行地表溫度降尺度的技術(shù)應(yīng)運(yùn)而生。

目前LST的空間降尺度方法主要分為：數(shù)理統(tǒng)計回歸法、光譜混合模型法和調(diào)制分配法，其中較為常用的是數(shù)理統(tǒng)計回歸法。數(shù)理統(tǒng)計回歸法的基礎(chǔ)假設(shè)是關(guān)系尺度不變，即LST與回歸核的統(tǒng)計關(guān)系在各個尺度上保持不變，回歸核包括反照率、組分權(quán)重和光譜指數(shù)等[5]。Kustas等[6]首次提出利用LST與NDVI的關(guān)系來實(shí)現(xiàn) LST降尺度的DisTrad算法，通過建立二者間的二次回歸模型，成功實(shí)現(xiàn)對千米級到百米級LST的降尺度。Agam等[5]對DisTrad算法進(jìn)行了改進(jìn)，通過在對 LST與4種植被指數(shù)進(jìn)行回歸分析的基礎(chǔ)上提出了TsHARP算法，該算法認(rèn)為利用LST與植被覆蓋度之間的一元線性關(guān)系能夠達(dá)到更好的降尺度效果。Essa等[7]指出對于下墊面復(fù)雜的城市地表，利用不透水面積指數(shù)建立模型比NDVI的效果更好。近幾年也有研究表明，對于不同的地表覆蓋特征，綜合利用NDVI、EVI、NDBI及UI等光譜指數(shù)對地表溫度降尺度的效果優(yōu)于單一指數(shù)的效果[8-10]。

上述研究已經(jīng)取得突破性的進(jìn)展，但仍存在一些不足。在分辨率方面，目前許多衛(wèi)星攜帶的多光譜成像儀獲得的影像分辨率都優(yōu)于30 m，將地表溫度表現(xiàn)在更小尺度內(nèi)的方法有待研究。在精度方面，許多研究先將高分辨率的LST數(shù)據(jù)聚合至低分辨率后進(jìn)行降尺度處理，原始LST數(shù)據(jù)在重采樣或升尺度過程中會產(chǎn)生新的誤差，對結(jié)果產(chǎn)生較大干擾。在結(jié)果驗證方面，多數(shù)研究用反演后的相對較高分辨率的LST數(shù)據(jù)對降尺度結(jié)果進(jìn)行驗證，由于溫度反演過程本身存在誤差，驗證結(jié)果的可靠性不強(qiáng)。針對以上問題，本文利用Sentinel-2衛(wèi)星的10 m分辨率數(shù)據(jù)，直接對反演的30 m分辨率LST進(jìn)行降尺度研究。對于因城市下墊面復(fù)雜導(dǎo)致的降尺度精度問題，基于對溫度圖像梯度的分析，通過利用空間濾波方式改善溫度梯度的分布來降低結(jié)果誤差。對于最終降尺度到10 m分辨率的LST數(shù)據(jù)，利用實(shí)地測量的溫度進(jìn)行驗證。

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

研究區(qū)包括天津的中心城區(qū)及附近郊區(qū)，如圖1所示。其地理范圍為116° 54′ 53″E～117° 22′ 48″ E、38° 55′ 27″N～39° 17′ 6″N，研究區(qū)屬于暖溫帶半濕潤大陸性季風(fēng)氣候，地勢平坦，地表覆蓋類型豐富。街道、建筑物、植被和水域等地物交錯分布，空間異質(zhì)性高，作為地表溫度降尺度的研究對象具有代表性。天津市春季降雨少、天氣晴朗，衛(wèi)星影像質(zhì)量高，因此選取的研究時段為4月。

圖1 研究區(qū)與測量方法

本文用于地表溫度反演的數(shù)據(jù)來自Landsat-8的OLI和TIRS傳感器，Landsat-8 OLI的1～7波段分辨率為30 m，TIRS的分辨率為100 m。成像時間為2017年4月，軌道號為122-33，影像下載于美國地質(zhì)勘探局(http：//glovis.usgs.gov/)。本文使用的高分辨率影像來自Sentinel-2 MSI傳感器。Sentinel-2 MSI共有13個波段的影像，其中，10 m高分辨率的波段有藍(lán)波段、綠波段、紅波段和近紅外波段。影像成像時間為2017年4月，軌道號為N0205-R032-T50SMJ和N0205-R032-T50SNJ，影像下載于科學(xué)網(wǎng)(https：//scihub.copernicus.eu/)。實(shí)地測量數(shù)據(jù)在2017年4月Landsat-8衛(wèi)星過境時采集，測量位置與方法如圖1所示。使用溫度計為主要測量工具，測量每一個點(diǎn)位時，將溫度計的感溫泡與地面充分接觸，并用手持GPS進(jìn)行定位。所有實(shí)測點(diǎn)均在衛(wèi)星過境時刻附近完成測量，因受影像成像時間限制，共得到33個實(shí)測點(diǎn)數(shù)據(jù)(表1)。

表1 測量地點(diǎn)名稱和LST

2 研究方法

2.1 地表溫度反演

有研究表明，單通道算法(single channel，SC)基于Landsat-8 TIRS數(shù)據(jù)反演地表溫度具有很好的精度和敏感性[11]，故本文使用SC進(jìn)行LST的反演。SC的計算方法如式(1)所示。

K2DTa}/(K2C)

C=τε

D=(1-τ)[1+(1-ε)τ]

Ta=16.011 0+0.926 21t0

T10=K2/ln(1+K1/L10)

L10=(gain)·DN+(offset)

(1)

式中：Ts為反演的地表溫度；T10為TIRS10的亮度溫度，單位為K；Ta為大氣平均作用溫度，單位為K；C和D為中間變量；τ為大氣透過率；ε為地表比輻射率；t0為地面附近的氣溫，單位為K；對于TIRS10，K1=774.89 W·m-2·sr-1·μm-1，K2=1 321.08 K；L10為TIRS10的熱輻射強(qiáng)度值，單位為W·(m2·sr·μm)-1；DN為影像的灰度值；gain為增益值；offset為偏移量，可以從數(shù)據(jù)頭文件中直接讀取。

對于ε的計算，城鎮(zhèn)和自然表面的地表比輻射率εbuilding和εsurface可由公式(2)求得。

εbuilding=PvRvεv+(1-Pv)Rmεm+dε

εsurface=PvRvεv+(1-Pv)Rsεs+dε

Rv=0.933 2+0.058 5Pv

Rm=0.988 6+0.128 7Pv

Rs=0.990 2+0.106 8Pv

Pv=(NDVI-NDVIs)/(NDVIv-NDVIs)

(2)

式中：Rv、Rm和Rs分別為植被、建筑和裸土的溫度比率；εv、εm和εs分別為植被、建筑和裸土純凈像元的地表發(fā)射率，其值分別取0.986、0.980和0.972；dε為比輻射率修正項，當(dāng)Pv≤0.5時，dε為0.003 8Pv，當(dāng)Pv>0.5時，dε為0.003 8(1-Pv)；Pv表示植被覆蓋度情況，可由NDVI來計算；NDVIv和NDVIs分別表示純植被像元與純裸土像元的NDVI值。本文通過統(tǒng)計影像NDVI直方圖，取5%和95%的累積百分比作為置信區(qū)間，截取NDVI的上下限閾值分別作為研究區(qū)的NDVIs和NDVIv值。

2.2 地表溫度降尺度

LST從高尺度向低尺度轉(zhuǎn)換時，通常需要補(bǔ)充額外的信息來構(gòu)造趨勢面。最早提出的DisTrad算法就是通過引入NDVI或其他指數(shù)構(gòu)造趨勢面，其基本假設(shè)是LST和植被指數(shù)的關(guān)系在各個尺度上基本一致[6]。以Landsat-8數(shù)據(jù)為例，首先建立30 m尺度上LST與趨勢面因子(NDVI)的關(guān)系，通過回歸分析構(gòu)造LST與NDVI之間的函數(shù)關(guān)系。假設(shè)此函數(shù)關(guān)系仍然適用于10 m尺度上的LST與NDVI，此不變性是降尺度技術(shù)的基礎(chǔ)，也是關(guān)鍵所在。最后，將10 m尺度上Sentinel-2的NDVI值代入此函數(shù)關(guān)系中，計算得到10 m尺度上的LST。在構(gòu)造30 m尺度上二者的函數(shù)關(guān)系時，應(yīng)選擇盡可能多的隨機(jī)點(diǎn)，以保證關(guān)系的可靠性。比較多種回歸方式，選擇最佳關(guān)系式可以確保得到最小殘差和最大相關(guān)系數(shù)。

2.3 溫度梯度評價與誤差評價

梯度反映了相鄰像元的亮度變化率，圖像中如果存在地物邊緣或者像元亮度差較大時，會出現(xiàn)較高的梯度值[12]；在地表溫度的圖像中，溫度像元差值較大時，亦是如此。由于光譜指數(shù)法降尺度的原則是，相同光譜指數(shù)值具有相同的溫度值，這就導(dǎo)致降尺度后的地表溫度在地物交界處易出現(xiàn)較大的波動。運(yùn)用梯度分析的方法可以有效檢測出這種波動，常用的梯度分析方法有羅伯特梯度、索伯爾梯度和拉普拉斯算法。索伯爾梯度與其他2種算法相比，既考慮了更多鄰域點(diǎn)的關(guān)系，也較好地保留了原圖像的特征[13]，故本文采用索伯爾梯度對溫度圖像的梯度值進(jìn)行計算，其計算如(3)所示。

|gradf|≌|t1|+|t2|

(3)

式中：gradf為梯度計算的結(jié)果；t1，t2為計算模板。

溫度圖像中的突變梯度值分布越多，溫度的高低變化越不連續(xù)。為了評價圖像中突變梯度值的分布，本文提出了梯度異常頻率(gradient anomaly frequency，GAF)，即圖像中突變梯度值占全部梯度值的比例，表達(dá)如式(4)所示。

(4)

式中：n為圖像像元的總個數(shù)；nm為突變梯度值的個數(shù)。對于突變梯度值的界定，本文選取了降尺度之前圖像中的100個地物邊界附近的梯度值并求其平均值，將大于或等于這個平均值的梯度值界定為梯度突變值。同時，本文也選擇了平均絕對誤差(MAE)和均方根誤差(RMSE)2個指標(biāo)進(jìn)行誤差評價，MAE和RMSE的計算分別如式(5)和式(6)所示。

(5)

(6)

2.4 濾波處理

熱傳導(dǎo)在二維平面的傳播可簡化為式(7)[14]。

(7)

式中：?u/?t為平面中一點(diǎn)的溫度對時間的變化率；?2u/?x2與?2u/?y2為溫度對2個坐標(biāo)軸的二次導(dǎo)數(shù)；k為熱擴(kuò)散率。此公式表明，任何一個不同的初始溫度，經(jīng)過一定時間都會與周圍形成熱平衡。因此，熱方程的結(jié)果具有將初始溫度平滑化的特質(zhì)，這與遙感圖像處理中的均值平滑模板作用十分相似。利用NDVI與LST的關(guān)系進(jìn)行降尺度后，由于沒有考慮到相鄰像元的關(guān)系，圖像上存在很多梯度異常值。根據(jù)溫度的熱傳導(dǎo)特性，本文利用均值濾波方法處理降尺度后的溫度圖像，以提高圖像中溫度的連續(xù)性。均值濾波的計算如式(8)所示。

(8)

式中：r(i，j)為計算結(jié)果；I(m，n)為降尺度后圖像的活動窗口；t(m，n)為均值濾波模板。

3 結(jié)果分析

3.1 地表溫度反演及降尺度

通過SC算法反演的30 m分辨率LST如圖2(a)所示。圖中LST變化比較連續(xù)，但由于分辨率低，許多細(xì)節(jié)處的溫度無法表現(xiàn)。在研究區(qū)內(nèi)隨機(jī)選取5 000個點(diǎn)，提取對應(yīng)點(diǎn)位的LST和NDVI數(shù)值并進(jìn)行擬合。經(jīng)過多次實(shí)驗發(fā)現(xiàn)，二次擬合模型(式(9))具有較好的相關(guān)性，其R2為0.615，擬合結(jié)果如圖2(c)所示。

圖2 反演的30 m分辨率LST與NDVI的擬合結(jié)果

LST=-213.19NDVI2+68.131NDVI+27.538

(9)

對Sentinel-2 MSI的波段數(shù)據(jù)進(jìn)行輻射校正和大氣校正后，計算得到10 m分辨率NDVI。利用NDVI與LST的二次擬合模型(式(9))，計算得到10 m分辨率LST，結(jié)果如圖3所示。

圖3 降尺度后的10 m分辨率LST

由圖3可以看出，光譜指數(shù)法降尺度的溫度大體分布特點(diǎn)和反演后的結(jié)果相同。隨著圖像分辨率的提高，更多地物溫度的細(xì)節(jié)被表現(xiàn)出來，對比反演的細(xì)節(jié)部分結(jié)果(圖3(b)、圖3(d))，很大程度上提高了圖像信息量(圖3(c)、圖3(e))，在地物覆蓋類型復(fù)雜的區(qū)域，能夠準(zhǔn)確地區(qū)分不同地物之間溫度的差別。然而此結(jié)果也存在較多誤差，由于未考慮相鄰像元的關(guān)系，圖像連續(xù)性較差，許多位置出現(xiàn)了較大的溫度差。尤其是河流與陸地之間，其溫度界限過于“明顯”，在實(shí)際環(huán)境中，這些溫差較大的區(qū)域并不能達(dá)到熱平衡狀態(tài)。

3.2 溫度梯度與誤差評價

利用索伯爾梯度算法對降尺度前后的溫度圖像梯度進(jìn)行檢測(圖4(a)、圖4(c))，并選取降尺度之前溫度圖像中地物邊界附近的100個梯度值的平均值作為突變梯度值nm，對GAF進(jìn)行計算，結(jié)果如圖4(b)、圖4(d)所示。由于反演的LST是根據(jù)熱紅外波段影像計算的結(jié)果，其保留了實(shí)際溫度的宏觀分布特點(diǎn)，所以將反演后LST的nm作為參考值計算其他溫度梯度圖像的GAF。

圖4 降尺度前后溫度圖像梯度與GAF結(jié)果

由圖4(a)可知，降尺度之前的溫度圖像梯度整體差別較小，異常梯度分布比較稀疏，異常梯度值只占了總數(shù)的3.1%。一般情況下，在地物屬性差別較大的區(qū)域(如水體和地面分界處)，不可避免會出現(xiàn)少量梯度異常值(圖4(b))，所以可以推斷，降尺度之前的溫度圖像基本符合熱平衡規(guī)律。而圖4(c)中，溫度梯度值高低交錯、變化頻繁，異常梯度分布十分密集，異常梯度值幾乎占總數(shù)的一半，尤其是中心區(qū)域，遍布了大量的高異常梯度值(圖4(d))。由此可見，光譜指數(shù)法降尺度后，像元存在嚴(yán)重的獨(dú)立性，幾乎和相鄰像元脫離了關(guān)系，造成溫度變化非常不均勻。通過圖5的誤差評價發(fā)現(xiàn)，溫度的準(zhǔn)確度也大幅降低。

降尺度后的MAE和RMSE達(dá)到了降尺度之前的2倍以上，大部分溫度與實(shí)際測量溫度差別較大，最大誤差達(dá)到了11.7 ℃，比降尺度之前的最大誤差高出6.9 ℃。由圖5同時可以發(fā)現(xiàn)，降尺度之前25～33號點(diǎn)位的溫度普遍低于實(shí)測溫度，而這部分點(diǎn)位大多處于高溫段。其原因可能是由于反演后溫度圖像分辨率低，溫度像元覆蓋的地面范圍比較大，不能準(zhǔn)確體現(xiàn)小尺度區(qū)域的高溫點(diǎn)信息所致。

圖5 降尺度前后LST誤差評價

3.3 LST降尺度結(jié)果優(yōu)化

對降尺度后的溫度圖像進(jìn)行均值濾波處理的結(jié)果如圖6所示。由圖6可知，圖像在經(jīng)過濾波處理后，溫度的整體分布趨勢無較大變化，而圖像中高溫點(diǎn)和低溫點(diǎn)明顯有所減少。這說明由于均值濾波的作用，圖像中的極端溫度值被削弱，極高或極低的溫度值減少。由溫度圖像的細(xì)節(jié)可以看出，溫度變化的連續(xù)性有大幅提高，溫度高低交錯的情況明顯減少(圖6(b)、圖6(c))。可見均值濾波處理改善了像元之間的關(guān)系，使大部分溫度呈均勻變化，均值濾波后的梯度和誤差如圖7所示。

圖6 均值濾波處理后的溫度圖像

由圖7(a)和圖4(c)對比可知，均值濾波處理后，大量高梯度值變成了低梯度值，部分梯度值降低到了0，整體梯度水平下降。GAF降低到了32.1%(圖7(b))，異常梯度值的分布相比處理之前更加稀疏，說明溫度突變值有相應(yīng)減少。對比圖7(c)與圖5(b)，MAE和RMSE均有小幅度下降，最大誤差值也下降了2.6 ℃，說明均值濾波處理后溫度圖像的精度有一定提升。然而，均值濾波后的GAF與反演后溫度圖像的GAF仍相差較大，若以反演后溫度圖像梯度作為參考的熱平衡狀態(tài)下的梯度，均值濾波后溫度像元之間的關(guān)系并不理想，GAF仍有很大的改善空間。

圖7 均值濾波處理后梯度檢測和誤差

4 討論

本文將梯度分析的方法引入到溫度圖像的分析中，獲得的溫度梯度圖像直觀地展示了溫度像元的連續(xù)性，是LST降尺度結(jié)果分析的一個重要補(bǔ)充。本文提出的GAF是對溫度圖像梯度的量化評價，同時包含了溫度的連續(xù)性信息。GAF可作為溫度圖像評價的一個重要指標(biāo)，對溫度反演及計算結(jié)果的合理性有一定指示作用。對于LST降尺度結(jié)果的優(yōu)化和降尺度方法的適用性，本文進(jìn)行了以下討論。

由于用3×3均值濾波模板處理后的溫度圖像梯度還有一定的改善空間，為了使GAF進(jìn)一步減小，本文通過擴(kuò)大模板選區(qū)來提升濾波強(qiáng)度，討論了用5×5和7×7像元范圍的濾波模板對降尺度后溫度圖像進(jìn)行處理的效果(圖8)。

由圖8可以看出，隨著濾波強(qiáng)度的增加，圖像梯度進(jìn)一步改善，溫度精度也再一次提高。5×5和7×7模板對圖像梯度的處理效果依次增強(qiáng)(圖8(a)、圖8(b))。尤其是圖8(b)中，GAF已經(jīng)非常接近反演后的水平，而且，異常梯度分布非常稀疏，主要分布在地物邊界附近。然而，在用5×5模板以上的濾波強(qiáng)度時，溫度圖像的精度幾乎無變化(圖8(c)、圖8(d))，說明圖像梯度改善到一定程度時，其對溫度精度的提高作用不明顯，而此時的溫度誤差已不再是由像元不連續(xù)造成的。7×7模板處理后的圖像不僅精度高，且具有最佳的GAF值，因此7×7模板達(dá)到了最優(yōu)的均值濾波效果。圖9展示了運(yùn)用7×7模板處理后溫度圖像更多的細(xì)節(jié)。

圖8 不同強(qiáng)度的均值濾波算子對GAF和誤差的影響

從圖9的三維溫度細(xì)節(jié)中可以看出，7×7均值濾波模板處理后的圖像效果很好，像元高低變化連續(xù)，無溫度突變的情況。從圖9中一處水域附近的溫度分布來看，極高溫度與極低溫度之間溫度變化自然、均勻，已基本符合熱平衡狀態(tài)的特點(diǎn)。

圖9 7×7均值濾波模板處理后的溫度圖像空間立體展示

均值濾波處理雖然有效改善了降尺度后圖像的梯度，也提高了溫度圖像的精度，但提高后的精度距離反演后的精度還有一定差距。通過對比圖10(a)和圖10(b)發(fā)現(xiàn)，均值濾波處理后的圖像在高溫區(qū)域和低溫區(qū)域還存在比較明顯的誤差，而這些誤差主要是在降尺度過程中出現(xiàn)的。對于光譜指數(shù)降尺度方法本身存在的缺陷，未來研究可以將反演的LST作為參考值，和降尺度后對應(yīng)位置的溫度像元建立聯(lián)系，從而控制降尺度后的溫度分布。若能發(fā)揮反演后LST的監(jiān)督作用，再利用本文的處理方法，應(yīng)能使降尺度后溫度圖像的精度進(jìn)一步提高。

圖10 高、低溫區(qū)域存在的誤差

本研究同時發(fā)現(xiàn)，GAF的降低雖然能夠改善光譜指數(shù)法降尺度后溫度圖像的精度，但它的值并不是越小越好。在本文所選的研究區(qū)中，GAF在10%左右變動時，溫度精度幾乎無變化，此時的溫度像元達(dá)到一個模擬的熱平衡狀態(tài)，其他區(qū)域是否也滿足這個規(guī)律還有待實(shí)驗驗證。此外，運(yùn)用SC算法反演的LST本身就存在誤差，若能利用真實(shí)測量溫度與反演的LST建立函數(shù)關(guān)系，并對反演的LST做初步優(yōu)化，再進(jìn)行后續(xù)操作，可能會取得更理想的效果。

5 結(jié)束語

本文以天津市中心城區(qū)及附近郊區(qū)為研究區(qū)，利用SC算法對LST進(jìn)行反演，并用光譜指數(shù)法對反演的LST進(jìn)行了空間降尺度，成功將30 m分辨率LST降低到了10 m。降尺度后的LST通過均值濾波處理進(jìn)行了改善，并基于梯度分析和真實(shí)測量溫度對結(jié)果進(jìn)行了評價。

運(yùn)用光譜指數(shù)法對反演的LST降尺度到10 m分辨率，可表現(xiàn)出更多地物溫度的細(xì)節(jié)。但由于該方法未考慮相鄰像元的關(guān)系，導(dǎo)致溫度圖像連續(xù)性差、GAF劇增、溫度精度降低。通過3×3均值濾波模板處理后，圖像連續(xù)性大幅提高；異常溫度梯度減少，GAF降低；溫度精度也有小幅度提高。加大均值濾波的強(qiáng)度，溫度梯度與溫度精度會進(jìn)一步改善。而濾波強(qiáng)度增大到7×7像元范圍時，溫度精度便很難再提升。運(yùn)用7×7模板處理后溫度圖像的效果最理想，像元高低變化連續(xù)，基本符合熱平衡狀態(tài)的特點(diǎn)。

均值濾波處理雖然有效改善了降尺度后溫度圖像的梯度和精度，但對于光譜指數(shù)法本身存在的缺陷無法彌補(bǔ)。未來的研究還需要綜合反演LST與真實(shí)測量溫度的特點(diǎn)建立相關(guān)關(guān)系，進(jìn)一步優(yōu)化LST降尺度的結(jié)果。