吳迪， 陳健， 石滿， 覃幫勇， 李盛陽

(1.南京信息工程大學(xué)遙感與測繪工程學(xué)院，南京 210044； 2.中國科學(xué)院空間應(yīng)用工程與技術(shù)中心太空應(yīng)用重點實驗室，北京 100094)

0 引言

地表溫度(land surface temperature, LST)在地—氣間物質(zhì)與能量交換過程中起著重要作用，是全球及區(qū)域尺度地表物理過程研究的關(guān)鍵參數(shù)，對生態(tài)、環(huán)境、水文、氣象氣候、生物地球化學(xué)以及農(nóng)業(yè)生產(chǎn)等研究具有重要意義[1-2]。傳統(tǒng)氣象站點可實時觀測獲取LST，但受站點分布、地形等因素影響，具有一定局限性且空間不連續(xù)[3-4]。隨著遙感技術(shù)的發(fā)展，目前已有許多傳感器可提供空間連續(xù)的LST產(chǎn)品，如EOS/MODIS，NOAA/AVHRR和FY/VISSR等，其中MODIS LST產(chǎn)品在天氣晴朗條件下精度可達1 K[5]。然而，受到云、氣溶膠、觀測角度和太陽光照角度等影響，各種LST產(chǎn)品在時間和空間上均存在不同程度的缺失[6]，這限制了LST遙感產(chǎn)品的應(yīng)用。時間序列數(shù)據(jù)重建技術(shù)可以充分利用多種統(tǒng)計或數(shù)值分析方法，模擬LST的時間變化規(guī)律，從而插補缺失的觀測值[7-8]，因此利用時間序列重建技術(shù)重建LST數(shù)據(jù)，進而得到較高時空分辨率的LST產(chǎn)品具有重要意義。

針對遙感反演的LST存在大量缺失值的情況，已有許多學(xué)者基于MODIS等數(shù)據(jù)開展了LST的重建研究。Nguyen等[9]利用歸一化植被指數(shù)(normalized difference vegetation index,NDVI)時間序列諧波分析法(harmonic analysis of NDVI time-series，HANTS)重建越南紅河三角洲MODIS LST數(shù)據(jù)，并分析了該地區(qū)LST及土地利用時空變化規(guī)律； Neteler[10]提出了一種基于LST直方圖識別云污染像元的方法，并根據(jù)溫度梯度重建了MODIS LST數(shù)據(jù)； 臧琳等[11]基于統(tǒng)計模型和濾波算法聯(lián)合的LST重建方法建立LST與地表亮溫的統(tǒng)計模型，然后利用濾波算法對基于統(tǒng)計模型的結(jié)果進行插值，基于LST時間序列的周期性進一步控制反演誤差。以上研究均有效重建了MODIS的LST數(shù)據(jù)，并進行了更加深入的應(yīng)用。但是，這些研究大多是針對MODIS等極軌衛(wèi)星數(shù)據(jù)開展的研究，目前，針對靜止衛(wèi)星LST產(chǎn)品的重建研究仍比較缺乏。相較于極軌衛(wèi)星，靜止衛(wèi)星能夠提供1 d內(nèi)多個時相的觀測數(shù)據(jù)，可更準確地獲取LST的時間變化規(guī)律，為準確進行LST重建提供了更多的信息。

風云2號衛(wèi)星是我國自主研發(fā)的第一代地球同步軌道靜止氣象衛(wèi)星，由2顆試驗星(FY-2A，B)和5顆業(yè)務(wù)星(FY-2C，D，E，F(xiàn)，G)組成，并根據(jù)觀測數(shù)據(jù)制作了多種定量產(chǎn)品數(shù)據(jù)。風云衛(wèi)星可以提供包括LST在內(nèi)的多種遙感反演產(chǎn)品，其中LST分為逐時、日平均、候平均、旬平均及月平均產(chǎn)品。本文基于FY-2F LST日數(shù)據(jù)，利用LST隨時間變化規(guī)律，基于Savitzky-Golay(S-G)濾波算法，對2013年長江三角洲(下文簡稱長三角)地區(qū)LST進行數(shù)據(jù)重建，以填補LST的缺失值，提高LST產(chǎn)品的可用性。基于該方法重建的FY-2F LST時間序列數(shù)據(jù)，可為長三角地區(qū)熱環(huán)境研究提供較好的數(shù)據(jù)基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)概況

長三角位于中國大陸東部沿海，是中國第一大經(jīng)濟區(qū)。根據(jù)國務(wù)院2010年批準的《長江三角洲地區(qū)區(qū)域規(guī)劃》，長三角包括江蘇省、浙江省和上海市2省1市，區(qū)域面積為21.07萬km2。其中上海市為平原； 江蘇省以平原為主，間或有一些山地丘陵； 浙江省多山地。整個長三角地區(qū)水網(wǎng)密布，自然資源豐富。屬亞熱帶濕潤季風氣候區(qū)，夏季高溫多雨，年均氣溫為15～17 ℃，年降水量約為1 000～1 800 mm。

長三角地區(qū)地理位置優(yōu)越，并且具有良好的經(jīng)濟基礎(chǔ)，在中心城市帶動下經(jīng)濟快速蓬勃發(fā)展，城市化進程也隨之加速[12]。伴隨著經(jīng)濟活動、城市化進程，長三角地區(qū)城市空間熱環(huán)境日益惡化，對環(huán)境、生態(tài)及人體健康造成極大的危害[13-14]。而LST是研究區(qū)域熱環(huán)境的重要因子，因此對長三角地區(qū)LST進行重建，從而獲得較高質(zhì)量的產(chǎn)品數(shù)據(jù)，對研究區(qū)域熱環(huán)境變化情況具有重要的意義。

1.2 數(shù)據(jù)源

FY-2F是靜止氣象衛(wèi)星風云二號03批的首發(fā)星，于2012年1月成功發(fā)射，搭載了2個主要載荷(掃描輻射計和空間環(huán)境監(jiān)測器)。FY-2F提供包括LST在內(nèi)的多種遙感反演產(chǎn)品，LST產(chǎn)品數(shù)據(jù)可從國家衛(wèi)星氣象中心的風云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)(http: //satellite.nsmc.org.cn/portalsite/default.aspx)免費訂購下載。

FY-2F衛(wèi)星能夠提供時間分辨率高達1 h的LST產(chǎn)品，為探究不同時間分辨率的數(shù)據(jù)質(zhì)量，分別統(tǒng)計了2013年長三角地區(qū)逐時產(chǎn)品及日產(chǎn)品的空間缺值率(即一景圖像空間上的缺值率)，結(jié)果如表1所示。

表1 FY-2F LST產(chǎn)品空間缺值率統(tǒng)計

從表1中可以清楚地看到，F(xiàn)Y-2F LST逐時產(chǎn)品空間缺值率較大，全年空間缺值率超過90%的影像高達21.78%； 相較于逐時產(chǎn)品，日LST數(shù)據(jù)質(zhì)量顯著提高，全年有57.80%的數(shù)據(jù)空間缺值率低于10%，僅有1.73%的影像空間缺值率大于90%?？紤]到逐時產(chǎn)品的缺值率過高，重建難度較高，且在研究長時間尺度地表熱環(huán)境變化時，應(yīng)用日產(chǎn)品效果更佳，因此將對FY-2F LST日產(chǎn)品進行重建以便更好地進行后續(xù)研究。本研究數(shù)據(jù)采用2013年FY-2F LST日均值產(chǎn)品，該數(shù)據(jù)存儲為HDF5格式，星下點空間分辨率為5 km，每幅影像可獲取覆蓋地球表面約1/3的圓盤標稱影像。由于數(shù)據(jù)不具備投影和經(jīng)緯度信息，因此根據(jù)風云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)提供的經(jīng)緯度數(shù)據(jù)進行了幾何糾正。

氣象數(shù)據(jù)為2013年長三角地區(qū)氣象站點日平均LST觀測數(shù)據(jù)。長三角地區(qū)共有46個氣象站點(江蘇省23個，浙江省22個，上海市1個)，空間分布如圖1所示，該數(shù)據(jù)可從中國氣象數(shù)據(jù)網(wǎng)免費下載(http: //data.cma.cn/site/index.html)。

圖1 研究區(qū)氣象站點分布示意圖

1.3 基于S-G濾波的LST重建方法

受研究區(qū)內(nèi)云量、傳感器等因素影響，LST遙感數(shù)據(jù)存在數(shù)值異?；蛉敝惮F(xiàn)象，F(xiàn)Y-2F LST逐時數(shù)據(jù)即存在大量缺值，即使是每日合成產(chǎn)品仍存在很多無效像元，嚴重影響了數(shù)據(jù)的可用性，因此能否精確地重建LST數(shù)據(jù)就尤為重要。本文基于LST長時序變化特征，并結(jié)合濾波算法重建LST缺值。S-G濾波由Savitzky和Golay[15]在1964年首次提出，是一種在時域內(nèi)基于局部多項式最小二乘擬合的濾波方法。通過取點xi相鄰時相(前后)固定個數(shù)點擬合多項式，多項式在xi的值即為S-G濾波的模擬值。該濾波器最大的特點為在濾除噪聲的同時可以保持曲線的變化趨勢。在對FY-2F LST數(shù)據(jù)重建的過程中，S-G濾波方法能夠最大程度地保留原始數(shù)據(jù)，并且根據(jù)相鄰數(shù)據(jù)模擬出缺值。S-G濾波過程可表示為

，

(1)

當像元長時間連續(xù)缺失較多數(shù)據(jù)時，若單純使用S-G濾波進行數(shù)據(jù)重建則可能與真實值相差較大。針對這一問題本文篩選出某像元一段不存在缺值的數(shù)據(jù)作為研究對象，從該段數(shù)據(jù)中分別提取1～5個相鄰時相數(shù)據(jù)不參加S-G濾波，模擬出像元不同連續(xù)缺值情況，并將重建前、后LST值進行對比。圖2中繪制出原始LST和5條不同連續(xù)缺值重建后曲線，方框代表模擬缺值，從圖中可以清楚地看到連續(xù)缺失數(shù)據(jù)越多，模擬結(jié)果與原始LST差值越大，當連續(xù)缺失達到5個數(shù)據(jù)時誤差大于3 K。

圖2 不同連續(xù)缺值重建效果

因此，本文首先遍歷每一像元所有時相LST數(shù)據(jù)，并根據(jù)像元連續(xù)缺值情況將全年數(shù)據(jù)拆分為若干組，當連續(xù)缺失5個及以上時相數(shù)據(jù)時該段LST不進行重建。隨后，根據(jù)實測數(shù)據(jù)統(tǒng)計出LST最大、最小值作為閾值，判斷LST異常點，利用接口描述語言(interface description language，IDL)函數(shù)庫提供的S-G濾波函數(shù)重建異常點及缺失數(shù)據(jù)，并且保留原始有效數(shù)據(jù)不變，通過反復(fù)迭代達到最佳擬合效果。S-G濾波函數(shù)需設(shè)定窗口寬度m及多項式擬合階數(shù)Degree，通常情況下m值越大濾波結(jié)果越平滑；Degree一般設(shè)定在2～4范圍內(nèi)，階數(shù)越低濾波結(jié)果平滑，但會帶來誤差，較高階數(shù)可降低此項誤差但“過擬合”會帶來噪聲影響濾波結(jié)果[16-17]。為了確定S-G濾波LST重建最優(yōu)參數(shù)，本文在連續(xù)30 d數(shù)據(jù)中隨機選擇5 d數(shù)據(jù)并賦為缺值，再按照窗口寬度m=3,5,7及階數(shù)Degree=2,3,4組成9組參數(shù)組合進行S-G濾波，經(jīng)實驗平滑窗口寬度m取5，階數(shù)Degree為4時模擬LST與原始值最為接近。

2 結(jié)果與分析

2.1 LST重建結(jié)果

本研究對2013年FY-2F LST時間序列數(shù)據(jù)進行了重建，通過IDL函數(shù)庫提供的S-G濾波函數(shù)迭代擬合，最終達到最優(yōu)擬合效果。S-G濾波能夠很好地將FY-2F LST數(shù)據(jù)缺失值補充完整，圖3展示了長三角地區(qū)每一像元點重建前、后時間序列缺值率(即表示每個像元點長時間序列上的缺值率，下文簡稱時相缺失率)。經(jīng)統(tǒng)計，長三角地區(qū)重建前FY-2F LST平均時相缺失率為19.43%，其中單一像元最少缺失9.54%數(shù)據(jù)，最高時相缺失率則達到39.31%，整個長三角地區(qū)南部數(shù)據(jù)時相缺失率高于北部。經(jīng)過重建后，F(xiàn)Y-2F LST時相缺失率明顯降低，長三角大部分區(qū)域LST補充完整，平均時相缺失率下降至1.69%，可見S-G濾波方法具有良好的重建效果。

(a) 重建前時相缺失率 (b) 重建后時相缺失率

為了更直觀地展現(xiàn)LST重建效果，圖4選取不同空間缺值率的FY-2F LST原始數(shù)據(jù)與重建后數(shù)據(jù)進行對比。從圖中可以看到，S-G濾波利用LST時序變化特征進行數(shù)據(jù)重建，重建效率較高，即便研究區(qū)內(nèi)空間缺值率超過80%也能很好地完成重建工作，且能夠較好地保證LST空間一致性。

(a) 20130204重建前(b) 20130316重建前(c) 20131221重建前(d) 20130716重建前

(e) 20130204重建后(f) 20130316重建后(g) 20131221重建后(h) 20130716重建后

為了從LST時間序列角度評價重建效果，本文進一步選取了南京站點作為驗證點，對比站點實測值與重建前、后LST隨時間變化曲線，如圖5所示。

(a) 2013年

由南京站點2013年LST隨時間變化總體曲線(圖5(a))可以看出，本文重建方法可以有效地保留原始LST，并達到對缺值進行補充的目的。實測數(shù)據(jù)與重建前、后LST相比變化趨勢一致，全年整體隨時間變化先升高后降低。從圖中也能看到，重建后的LST時間序列與實測數(shù)據(jù)仍然存在一定的誤差，其中8月份的差值達到了10 K左右。為了更清楚地觀察重建效果，分別提取4月及8月份數(shù)據(jù)如圖5(b)和(c)。S-G濾波方法對LST時間序列進行濾波獲取其長時間變化趨勢，結(jié)合相鄰時相點對缺值進行補充，使其符合LST變化規(guī)律。圖5(b)中，實測數(shù)據(jù)與重建后FY-2F LST變化趨勢較為一致，相鄰2個時相FY-2F數(shù)據(jù)波動更大； 8月份(圖5(c))存在連續(xù)缺失數(shù)據(jù)較多的情況，為了保證重建數(shù)據(jù)準確性該段缺值未進行補充，其余部分缺值經(jīng)S-G濾波擬合效果較好，但與實測LST差值偏大，這種情況受數(shù)據(jù)本身精度的影響較大。

2.2 精度驗證

為了進一步定量檢驗重建后LST結(jié)果的精度，根據(jù)經(jīng)緯度提取了長三角地區(qū)46個氣象站點重建前、后LST每日數(shù)據(jù)，并與相應(yīng)日期實測LST數(shù)據(jù)進行對照。圖6(a)和(b)分別展示了氣象站點實測數(shù)據(jù)與LST重建前、后對比情況。從圖中可以看到FY-2F LST原始產(chǎn)品與實測數(shù)據(jù)具有較好的擬合精度，判定系數(shù)R2為0.75，而平均絕對誤差(mean absolute error,MAE)為4.88 K，原始LST產(chǎn)品精度較低。重建后樣本數(shù)由11 201個增加至13 621個，R2及MAE與重建前基本保持一致，因此利用S-G濾波對FY-2F LST進行重建能夠保證數(shù)據(jù)精度，重建結(jié)果較為可信。

(a) 實測值與重建前LST (b) 實測值與重建后LST

經(jīng)過對比發(fā)現(xiàn)，氣象站點實測值與重建后LST值誤差較大，主要有以下3點原因： ①FY-2F原始LST產(chǎn)品在某些地區(qū)誤差較大，甚至超過10 K,由圖5可以看出，高誤差像元主要出現(xiàn)在夏季7—8月份，這可能受FY-2F LST產(chǎn)品本身反演算法的影響，也可能是該地區(qū)夏季云量較多，導(dǎo)致逐時數(shù)據(jù)缺失，進而影響了日平均LST數(shù)據(jù)的計算； ②氣象站點實測值為“點數(shù)據(jù)”，而FY-2F LST數(shù)據(jù)為“面數(shù)據(jù)”，且空間分辨率為5 km，導(dǎo)致像元多為混合像元，由于兩者尺度不同，在匹配時必然帶來一定誤差[18-19]； ③本文是在FY-2F LST產(chǎn)品的基礎(chǔ)上進行重建的，受原始產(chǎn)品質(zhì)量影響，無法解決一些時相的高誤差問題。因此，重建后產(chǎn)品與重建前產(chǎn)品相比更加合理。

綜合以上分析，為了更加明確S-G濾波方法重建效果，本文隨后采取人工模擬缺值的方法排除由于數(shù)據(jù)問題導(dǎo)致的影響，對比重建前、后數(shù)據(jù)驗證S-G濾波重建數(shù)據(jù)的準確性(圖7)。

圖7 重建前、后LST散點圖

通過隨機選取一些非缺值區(qū)域不參與S-G濾波，然后跟濾波之后的LST結(jié)果進行對比。選取驗證點的時候，要求在長三角地區(qū)內(nèi)均勻分布，且影像按照不同季節(jié)隨機抽取，最大程度保證了驗證結(jié)果的準確性。從圖7中可見，重建前、后LST的MAE為1.35 K，R2為0.95，大部分樣本都位于1∶ 1線附近，重建后LST與原始值吻合度較好。因此，除去數(shù)據(jù)質(zhì)量等其他因素影響，利用S-G濾波方法進行數(shù)據(jù)重建效果良好。

為了進一步評價本文重建方法的魯棒性，本文使用相同方法對2014年和2015年FY-2F LST日均值數(shù)據(jù)進行了重建，評價其R2和MAE，結(jié)果如表2所示。2014—2015年間站點實測值與重建前、后LST的R2和MAE基本持平，能夠保證重建后LST數(shù)據(jù)精度不變； 人工模擬缺值點驗證結(jié)果良好，與2013年結(jié)果大體相同，模型魯棒性較好。

表2 2014—2015年LST重建精度評價

3 結(jié)論

本文利用S-G濾波方法，結(jié)合LST時間序列變化特征，對長三角地區(qū)FY-2F LST產(chǎn)品進行了數(shù)據(jù)重建，并對重建結(jié)果進行了評價。結(jié)果表明：

1)長三角地區(qū)FY-2F LST原始產(chǎn)品平均時相缺失率為19.43%，S-G濾波方法能夠有效地重建LST，重建后平均時相缺失率降低為1.69%。并且對于不同空間缺值率影像，該方法均能在保持LST空間一致性的基礎(chǔ)上對其進行有效重建。

2)提取長三角地區(qū)46個氣象站重建前、后LST值，結(jié)合站點LST進行精度驗證與誤差分析。重建后結(jié)果精度與重建前FY-2F LST產(chǎn)品基本持平，但與實測值間誤差較大，其原因為： 產(chǎn)品本身質(zhì)量影響了重建結(jié)果精度； 進行驗證時，實測值為“點數(shù)據(jù)”，與遙感“面數(shù)據(jù)”因尺度問題會進一步擴大誤差。

3)通過人為模擬驗證區(qū)缺值，將重建后結(jié)果與重建前對比，發(fā)現(xiàn)重建精度較高，兩者R2為0.95，MAE為1.35 K。通過重建2014—2015年間FY-2F LST日均值產(chǎn)品驗證S-G濾波模型適用性，其結(jié)果與2013年大體相同，表明S-G濾波法可以用來進行LST重建。

本文在保持FY-2F LST原始產(chǎn)品精度的基礎(chǔ)上進行數(shù)據(jù)重建，能夠有效提高產(chǎn)品質(zhì)量，為地表熱環(huán)境研究提供數(shù)據(jù)基礎(chǔ)。然而文中重建采用的S-G濾波方法僅考慮了LST時序變化規(guī)律，并未利用到其空間信息。在今后的研究中可結(jié)合多源遙感數(shù)據(jù)和下墊面特性等因子進行空間上的重建，進一步提高產(chǎn)品數(shù)據(jù)質(zhì)量。