基于FY-3D衛(wèi)星的微波與光學(xué)陸表溫度融合研究*

劉勇洪 翁富忠 徐永明 韓秀珍 段四波 唐世浩 葉成志

1 中國氣象局地球系統(tǒng)數(shù)值預(yù)報(bào)中心,北京 100081 2 災(zāi)害天氣國家重點(diǎn)實(shí)驗(yàn)室,北京 100081 3 南京信息工程大學(xué)遙感與測繪工程學(xué)院,南京 210044 4 國家衛(wèi)星氣象中心,北京 100081 5 中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所,北京 100081 6 湖南省氣象局,長沙 410118

提 要：目前還沒有基于國產(chǎn)衛(wèi)星的1 km分辨率的全天候陸表溫度(LST)產(chǎn)品,FY-3D 衛(wèi)星提供了中分辨率成像儀(MERSI)Ⅱ型1 km分辨率晴空LST產(chǎn)品與微波成像儀(MWRI)25 km全天候LST產(chǎn)品,因此可結(jié)合兩者優(yōu)勢開展全天候1 km 分辨率LST的融合研究。基于地理加權(quán)回歸(GWR)方法,選擇海拔、FY-3D 歸一化植被指數(shù)和歸一化建筑指數(shù)等建立GWR模型對FY-3D/MWRI 25 km LST降尺度到1 km,并與MERSI 1 km LST進(jìn)行融合;同時(shí)針對MWRI軌道間隙,利用前后1天融合后的云覆蓋像元1 km LST進(jìn)行補(bǔ)值,可以得到接近全天候下的1 km LST?；谝陨先诤纤惴?選擇了中國區(qū)域多個(gè)典型日期FY-3D/MERSI和MWRI LST官網(wǎng)產(chǎn)品進(jìn)行了融合試驗(yàn),并利用公開發(fā)布的全天候1 km LST產(chǎn)品(TPDC LST)對FY-3D 1 km LST融合結(jié)果進(jìn)行了評估。研究結(jié)果表明,基于GWR法的LST降尺度方法,可以有效避免傳統(tǒng)微波LST降尺度方法中存在的“斑塊”效應(yīng)和局地溫度偏低等問題;LST融合結(jié)果有值率從融合前的22.4%～36.9%可提高到融合后69.3%～80.7%,融合結(jié)果與TPDC LST的空間決定系數(shù)為0.503～0.787,均方根誤差為3.6～5.8 K,其中晴空為2.6～4.9 K,云下為4.1～6.1 K;分析還表明目前官網(wǎng)產(chǎn)品FY-3D/MERSI和MWRI LST均存在缺值較多與精度偏低等問題,顯示其存在較大改進(jìn)潛力,這有利于進(jìn)一步改進(jìn)FY-3D LST融合質(zhì)量。

引 言

陸表溫度(LST)是研究全球與區(qū)域氣候變化的重要參數(shù)(Li et al,2013;Yu et al,2018),特別是獲取1 km空間分辨率的全天候LST已經(jīng)成為相關(guān)研究和應(yīng)用領(lǐng)域的迫切需求(Hijmans et al,2005;Fick and Hijmans,2017;Zhang et al,2020)。基于光學(xué)(可見光與熱紅外)遙感數(shù)據(jù)所獲得的LST,雖然空間分辨率高(一般為1 km或優(yōu)于1 km),反演精度較高(Wan,2008;權(quán)維俊等,2012;Liu et al,2015;Duan et al,2019),但是受限于晴空天氣條件,存在大量的缺失像元或無效像元(Cornette and Shanks,1993),不能進(jìn)行LST的全天候監(jiān)測。不同于光學(xué)波段,微波可以穿透云層,提供了獲取全天候LST的有效手段,但空間分辨率遠(yuǎn)低于光學(xué)反演得到的LST(一般為10～36 km),LST反演精度較低(McFarland et al,1990;Zhou et al,2015),且軌道間存在間隙,反演的LST空間不連續(xù)。因此,結(jié)合兩者優(yōu)勢開展微波與光學(xué)LST的融合目前已成為全天候1 km分辨率LST研究熱點(diǎn)(Duan et al,2017;Xu et al,2019;Yoo et al,2020;Zhang et al,2020;Tang et al,2022)。

目前在微波與光學(xué)LST融合研究中,采用的衛(wèi)星數(shù)據(jù)仍以國外為主,其中微波主要采用AMSR-E(Advanced Microwave Scanning Radiometer for EOS)和AMSR-2(Advanced Microwave Scanning Radiometer 2)數(shù)據(jù),光學(xué)主要為MODIS(Moderate-resolution Imaging Spectroradiometer)數(shù)據(jù),而利用我國衛(wèi)星開展全天候的1 km分辨率LST融合研究尚未見報(bào)道。風(fēng)云三號(hào)氣象衛(wèi)星D星(FY-3D)是我國自主研制的最新一代極軌氣象衛(wèi)星,其上搭載的中分辨率成像光譜儀(MERSI-Ⅱ)是世界上首臺(tái)能夠獲取全球 250 m分辨率紅外分裂窗區(qū)資料的成像儀器,其中熱紅外分裂窗波段第24和第25波段與MODIS分裂窗波段光譜接近(蔣金雄等,2019;Tang et al,2021),已有研究表明FY-3D MERSI-Ⅱ具有與MODIS類似的LST產(chǎn)品反演精度(Wang et al,2019;Du et al,2021;Aveni and Blackett,2022);而FY-3D衛(wèi)星攜帶的微波成像儀(MWRI)提供了5個(gè)觀測頻率波段(10.65、18.7、23.8、36.5和89.0 GHz)垂直極化和水平極化共10個(gè)可接收的亮溫通道(朱瑜馨等,2021;王博等,2022;劉梅等,2022),由于FY-3D MERSI-Ⅱ與MWRI在過境時(shí)間上具有一致性,故利用FY-3D MERSI-Ⅱ反演的陸表溫度(簡稱FY-3D/MERSI LST)與MWRI反演的LST(簡稱FY-3D/MWRI LST)開展全天候1 km分辨率LST融合研究,對于國產(chǎn)氣象衛(wèi)星在全球陸表生態(tài)環(huán)境監(jiān)測與全球氣候變化的應(yīng)用具有重要意義。

在全天候1 km分辨率光學(xué)與微波LST融合研究中,最為關(guān)鍵的技術(shù)為LST降尺度技術(shù)(Zhan et al,2013;Chen et al,2014)。目前LST降尺度方法主要有三大類:時(shí)空相鄰融合法、機(jī)器學(xué)習(xí)法和統(tǒng)計(jì)回歸法。其中,時(shí)空相鄰融合法是對具有不同時(shí)間尺度和不同空間尺度的同一遙感參數(shù)數(shù)據(jù)(如LST、輻亮度、反射率等),利用粗、高分辨率影像光譜相似性原理和混合像元分解法,并結(jié)合時(shí)空相鄰差異與地表類型差異構(gòu)建尺度轉(zhuǎn)換因子從而達(dá)到LST降尺度,主要代表性模型有時(shí)空自適應(yīng)反射率融合模型(STARFM)(Gao et al,2006)、自適應(yīng)時(shí)空融合變化檢測模型(STAARCH)(Hilker et al,2009)、增強(qiáng)時(shí)空自適應(yīng)反射率融合模型(ESTARFM)(Zhu et al,2010;Long et al,2020)、自適應(yīng)數(shù)據(jù)融合模型(FSDAF)(Zhu et al,2016),這些模型主要針對具有相似熱紅外光譜特征的粗、高分辨率LST融合,而微波與熱紅外光譜特征存在明顯差異,因此這些方法難以應(yīng)用于微波LST降尺度中。而機(jī)器學(xué)習(xí)法通常采用神經(jīng)網(wǎng)絡(luò)(ANN)和隨機(jī)森林(RF)模型進(jìn)行降尺度,如Shwetha and Kumar(2016)采用神經(jīng)網(wǎng)絡(luò)模型,考慮地形、土地覆蓋類型、微波植被指數(shù)等對AMSR-E/2 25 km 亮溫降尺度到1 km來模擬全天候LST;Yoo et al(2020)針對韓國地區(qū)采用RF模型進(jìn)行AMSR-2 LST降尺度,除了考慮地形因子和不透水蓋度等外,還加入了該地區(qū)與LST有關(guān)的3個(gè)年尺度參數(shù)(年均LST、LST振幅和相對于春分的LST時(shí)相差),而這些變量需要利用長時(shí)間序列LST資料估算獲取;Zhong et al(2021)針對中國青藏高原地區(qū)利用AMSR-2衛(wèi)星數(shù)據(jù)LST降尺度中,在RF模型中考慮了地形因子、歸一化植被指數(shù)(NDVI)和3個(gè)微波通道垂直極化亮溫等。由于機(jī)器學(xué)習(xí)模型需要構(gòu)建大量的代表性訓(xùn)練樣本,且輸入的變量較多不宜獲取,建模過程較為復(fù)雜,建立的模型在其他地區(qū)并不一定適用。

基于統(tǒng)計(jì)回歸的LST降尺度方法假定在低空間分辨率像元尺度上建立的LST與其他相關(guān)地表參數(shù)(或回歸因子)之間的統(tǒng)計(jì)關(guān)系具有空間尺度不變性,可將其應(yīng)用到高空間分辨率像元尺度上,因此,該方法不受光譜特征的限制。最初廣為應(yīng)用的LST統(tǒng)計(jì)降尺度方法是DisTrad(disaggregation procedure for radiometric surface temperature)模型(Kustas et al,2003),該方法將LST與NDVI的負(fù)相關(guān)關(guān)系應(yīng)用于LST的降尺度中;考慮到NDVI在高密度植被區(qū)的飽和性,Agam et al(2007)利用植被覆蓋度替代NDVI發(fā)展了TsHARP(thermal data sharpening)模型進(jìn)行LST的降尺度。進(jìn)一步,考慮到城市地表的特性,Dominguez et al (2011)在NDVI的基礎(chǔ)上加入地表反照率,可以有效提高城市地區(qū)LST降尺度精度;Wang et al (2020)研究顯示,在城市地區(qū)采用歸一化建筑指數(shù)(NDBI)較NDVI具有更好的LST降尺度效果。以上研究主要利用單變量或者多變量與LST之間回歸關(guān)系來預(yù)測LST,這些回歸模型均假定LST和回歸變量之間的關(guān)系具有空間不變性。事實(shí)上,LST與這些回歸變量(例如NDVI、地表高程DEM)之間的關(guān)系是隨地理位置變化的,例如,Jeganathan et al (2011)基于TsHARP模型對比了全局與局部并結(jié)合土地利用的分層LST降尺度模型,研究發(fā)現(xiàn),局部的LST模型最為可靠。因此,由于地形地理因素,各回歸變量如LST、NDVI、NDBI存在空間異質(zhì)性,隨著樣本的地理空間位置改變,變量之間的關(guān)系結(jié)構(gòu)也發(fā)生相應(yīng)變化,即LST與回歸變量之間的關(guān)系具有空間局地性特征。由此,本文引入考慮局地空間異質(zhì)性的地理加權(quán)回歸(GWR)法(Duan and Li,2016)來開展FY-3D/MWRI LST從25 km 分辨率降尺度到1 km 分辨率研究,并與過境時(shí)間相同的FY-3D/MERSI LST進(jìn)行融合,從而開展基于FY-3D的全天候1 km LST融合研究,這對于發(fā)展基于國產(chǎn)氣象衛(wèi)星的長時(shí)間序列衛(wèi)星氣候數(shù)據(jù)集具有積極意義。

1 試驗(yàn)數(shù)據(jù)

1.1 研究數(shù)據(jù)

選取中國區(qū)域一年代表月份(1、4、7、10月)典型日(2019年7月14—16日、2019年10月14—16日、2020年1月12—14日、2020年4月16—18日)白天FY-3D/MERSI LST和25 km FY-3D/MWRI LST白天(升軌)產(chǎn)品以及所在旬FY-3D/MERSI NDVI產(chǎn)品(王圓圓和李貴才,2022),這些產(chǎn)品均來源于中國氣象局風(fēng)云衛(wèi)星產(chǎn)品官方網(wǎng)站“風(fēng)云衛(wèi)星遙感數(shù)據(jù)服務(wù)網(wǎng)”(http:∥satellite.nsmc.org.cn/portalsite/default.aspx)。其中FY-3D/MERSI LST產(chǎn)品空間分辨率為1 km,采用Hammer投影,LST主要采用Becker局地分裂窗改進(jìn)算法(董立新等,2012;楊軍和董超華,2011)反演得到;FY-3D/MWRI LST產(chǎn)品空間分辨率為25 km,采用全球圓柱等積投影(EASE-GRID),LST反演算法采用基于土地覆蓋分類的多通道回歸模型(楊軍和董超華,2011),并假設(shè)每種類型地表在各個(gè)被動(dòng)微波通道具有一致的比輻射率;引入土地覆蓋分類保證了每種類型的發(fā)射率的相對穩(wěn)定,可提高LST的反演精度(武勝利和楊虎,2007)。

文中所涉及的中國地圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號(hào)GS(2019)1822號(hào)的標(biāo)準(zhǔn)地圖制作,底圖無修改。

1.2 產(chǎn)品驗(yàn)證數(shù)據(jù)

由于沒有足夠的均一性站點(diǎn)LST觀測數(shù)據(jù)進(jìn)行1 km尺度的LST驗(yàn)證,這里選取由國家青藏高原數(shù)據(jù)中心(TPDC)網(wǎng)站(https:∥data.tpdc.ac.cn/zh-hans/)公開發(fā)布的“中國陸域及周邊逐日1 km 全天候地表溫度數(shù)據(jù)集”(簡稱TPDC LST)(周紀(jì)等,2021)對FY-3D 1 km LST融合結(jié)果進(jìn)行評估。TPDC LST基于1 km MODIS/Aqua LST產(chǎn)品和0.25°空間分辨率GLDAS(Global Land Data Assimilation System)以及0.0625°空間分辨率CLDAS(China Land Data Assimilation System)等再分析數(shù)據(jù)生成(Zhang et al,2021),其中再分析數(shù)據(jù)主要用于云下LST的融合,且CLDAS主要針對中國區(qū)域,并融合了地面氣象站點(diǎn)LST觀測數(shù)據(jù)(孫帥等,2017;2022;師春香等,2018;王莉莉和龔建東,2018);利用分布于黑河流域、東北、華北和華南地區(qū)的15個(gè)站點(diǎn)實(shí)測數(shù)據(jù)對TPDC LST產(chǎn)品評估表明,其平均偏差(MBE)為-1.17～-0.06 K,均方根誤差(RMSE)為1.52～3.71 K,且在晴空與非晴空條件下無顯著區(qū)別。

此外,為了評估FY-3D/MERSI LST質(zhì)量,這里還選取了與FY-3D/MERSI LST日期相同的MODIS Aqua LST產(chǎn)品V6版本(MYD11A1)進(jìn)行分析,空間分辨率為1 km,采用正弦曲線投影,數(shù)據(jù)來源于美國航天局網(wǎng)站(https:∥earthdata.nasa.gov/)。

2 研究方法

2.1 地理加權(quán)回歸模型

地理加權(quán)回歸(GWR)方法是傳統(tǒng)標(biāo)準(zhǔn)回歸方法(如普通最小二乘方法)的拓展,該方法是進(jìn)行局部而不是全局參數(shù)估算(Fotheringham et al,2002)。通過利用地理空間變化的回歸系數(shù),來研究自變量與因變量之間的動(dòng)態(tài)回歸關(guān)系。GWR模型可表示為:

(1)

式中:μi和υi是第i個(gè)位置的地理坐標(biāo),yi和χik分別是第i個(gè)位置的因變量和第k個(gè)自變量,β0和βk分別為第i個(gè)位置的回歸模型的截距和系數(shù)(斜率),εi為第i個(gè)位置的回歸殘差。

β(μi,υi)=[XTW(μi,υi)X]-1XTW(μi,υi)Y

(2)

式中:β(μi,υi)是各回歸系數(shù)β的無偏估計(jì),W(μi,υi)是權(quán)重矩陣,其確保離目標(biāo)點(diǎn)更近的相鄰觀測點(diǎn)有更大權(quán)重,X和Y分別是自變量和因變量矩陣。

空間權(quán)重W(μi,υi)是GWR模型的核心,常用的有距離閾值函數(shù)法、距離反比法、指數(shù)函數(shù)法、高斯函數(shù)法和雙平方函數(shù)法,本研究采用自適應(yīng)雙平方函數(shù)(bi-square)法:

(3)

式中:Wi,j是當(dāng)前回歸點(diǎn)位置i與鄰近觀測點(diǎn)位置j之間的權(quán)重,dij是兩點(diǎn)之間的歐式距離,hi是自適應(yīng)核的波寬(或帶寬)。在本研究中利用改進(jìn)的AIC(Akaike information criterion)準(zhǔn)則來確定寬帶:

(4)

式中:n是觀測樣本數(shù),S是帽子矩陣,tr(S)為矩陣的秩,RSS是剩余平方和。

2.2 LST降尺度變量選擇

在LST降尺度中,選擇合適的輔助變量進(jìn)行LST降尺度是非常關(guān)鍵的。許多研究表明,NDVI和DEM與LST關(guān)系密切,這是由于作為與地表覆蓋相關(guān)的植被活動(dòng)和生物量的指示器,NDVI可用來考慮不同地表覆蓋類型對地表熱過程的影響,而DEM被認(rèn)為是表征LST變化的另一個(gè)重要因素(吳芳營等,2022)。此外,地理經(jīng)度(Lon)反映了從海岸到陸地區(qū)域土壤濕度對LST的影響,地理緯度(Lat)反映了不同緯度帶太陽輻射對LST的影響(Ke et al,2013)。另外,一些研究還表明在LST降尺度過程中增強(qiáng)型植被指數(shù)(EVI)較NDVI具有更好的表現(xiàn)(Qiu et al,2018),歸一化沙塵指數(shù)(NDSI)在荒漠地區(qū)有較好表現(xiàn)(Pan et al,2018),在半干旱地區(qū)及城市地區(qū)歸一化建筑指數(shù)NDBI較NDVI具有更好表現(xiàn)Wang et al(2020)。

本研究針對FY-3D/MWRI LST建模,選擇變量的一個(gè)重要原則是考慮數(shù)據(jù)業(yè)務(wù)生產(chǎn)的便利性,盡量避免選擇難以獲取的變量,而是從FY-3D產(chǎn)品中直接獲取相關(guān)變量。由于受云雨影響,1 km分辨率的 FY-3D MERSI日NDVI產(chǎn)品存在大量缺值,而旬NDVI產(chǎn)品則基本上克服了這種缺值問題。假定旬內(nèi)NDVI變化很小,由此可選用FY-3D旬NDVI合成產(chǎn)品作為輔助變量;此外旬NDVI合成產(chǎn)品中包含了短波紅外通道(中心波長為1.380 μm)和近紅外通道(中心波長為0.865 μm)反射率值,由此,可以基于短波紅外通道反射率RSWIR和近紅外通道反射率RNIR,估算得到歸一化建筑指數(shù)(NDBI):

(5)

由此,本研究選擇FY-3D旬歸一化植被指數(shù)NDVI、NDBI產(chǎn)品和DEM作為LST降尺度變量;由于地理經(jīng)度Lon和緯度Lat在GWR模型中己經(jīng)是輸入?yún)?shù),因此不再于模型中顯示。

2.3 FY-3D/MWRI LST與MERSI LST融合

FY-3D/MWRI 25 km LST與FY-3D/MERSI 1 km LST融合技術(shù)流程如圖1所示,主要包括FY-3D/MWRI LST訂正、基于GWR模型的FY-3D/MWRI LST降尺度、軌道間隙(Gap)補(bǔ)值等步驟。

圖1 逐日FY-3D 25 km MWRI LST與MERSI 1 km LST融合技術(shù)流程圖

2.3.1 FY-3D/MWRI LST訂正

對2019年7月14—15日、2019年10月15—16日和2020年4月17—18日FY-3D/MWRI LST分析顯示(圖2),晴空天氣條件下白天FY-3D/MWRI 25 km LST較FY-3D/MERSI 1 km LST明顯偏低,而且是系統(tǒng)性的偏低,平均偏差MBE(MWRI LST與MERSI LST之差平均值)小于-5 K,部分在-10 K以下,這與高浩等(2018)對FY-3C/MWRI LST評估顯示其LST明顯偏低結(jié)果相似,表明FY-3/MWRI LST存在明顯的低估。由于MWRI LST與MERSI LST存在明顯空間相關(guān),決定系數(shù)R2一般在0.75以上,因此可以利用線性回歸模型基于FY-3D/MERSI LST對FY-3D/MWRI LST進(jìn)行偏差訂正,在此基礎(chǔ)上開展兩者融合,主要步驟為:

注:深藍(lán)色線為LST 1∶1線,紅色線為線性回歸方程線。

(1)25 km分辨率晴空率判別:根據(jù)1 km分辨率MERSI LST文件中的有效值,估算25 km分辨率的晴空像元百分比;同時(shí)對MERSI 1 km分辨率LST文件,利用均值法重采樣生成MERSI 25 km分辨率LST。

(2)MWRI LST偏差訂正:根據(jù)25 km分辨率的晴空像元百分比,取值為1.0時(shí)(完全晴空),獲取相應(yīng)位置的FY-3D/MWRI 25 km分辨率LST值和重采樣后的25 km分辨率的MERSI LST值,對這二者LST進(jìn)行最小二乘法線性擬合,生成自變量為MWRI LST、因變量為MERSI LST的關(guān)系式,依據(jù)該關(guān)系式對原始FY-3D/MWRI 25 km分辨率LST進(jìn)行偏差訂正。

經(jīng)過上述步驟后,MWRI 25 km LST偏低問題可得到明顯改善所示,如圖3所示,2019年7月15日、2019年10月15日和2020年4月17日MWRI LST值得到明顯提升,與FY-3D LST的平均絕對偏差分別降至2.2、2.3和2.7 K。

圖3 典型日期下的FY-3D/MWRI 25 km LST訂正前后對比

2.3.2 基于GWR模型的FY-3D/MWRI LST降尺度

對經(jīng)過系統(tǒng)偏差訂正后的FY-3D/MWRI 25 km LST,利用GWR模型降尺度到1 km LST,主要步驟包括:

(1) 空間聚合:根據(jù)有效像元平均法將當(dāng)日所在旬合成1 km NDVI、NDBI影像以及DEM進(jìn)行空間聚合,生成25 km分辨率的NDVI、NDBI及DEM影像。

(2) GWR建模:利用GWR法建立訂正后25 km 分辨率MWRI LST與25 km分辨率NDVI、NDBI和DEM之間的動(dòng)態(tài)回歸關(guān)系,其可表示為:

LSTi=β0(μi,υi)+β1(μi,υi)NDVIi+

β2(μi,υi)NDBIi+β3(μi,υi)DEMi+εi

(6)

各參數(shù)含義見式(1)。由于覆蓋中國區(qū)域的FY-3D/MWRI LST產(chǎn)品需要3條軌道,而不同軌道相隔較遠(yuǎn),中間存在大量無值區(qū)域,如圖4所示,為避免樣本在空間不均一帶來的誤差,在建立該景LST地理加權(quán)回歸模型時(shí),需對每條軌道覆蓋中國區(qū)域LST樣本單獨(dú)進(jìn)行建模,即每景中國區(qū)域MWRI 25 km LST影像需要建立3個(gè)GWR回歸模型。

圖4 典型日期下的FY-3D/MWRI 25 km LST有效覆蓋范圍

(4)1 km LST模擬:利用前面得到的1 km回歸系數(shù)和1 km 分辨率NDVI、NDBI和DEM等變量,進(jìn)行1 km LST模擬預(yù)測:

(7)

2.3.3 當(dāng)日1 km LST融合

基于當(dāng)日FY-3D/MERSI 1 km LST與降尺度后的MWRI 1 km LST進(jìn)行融合,即:MERSI LST有值部分仍保留,缺值部分用降尺度后的1 km LST模擬值代替。

2.3.4 Gap補(bǔ)值

由于MWRI全球各掃描軌道之間尤其是中低緯度存在明顯空隙(Gap),且每天均不固定,如圖4所示,因此MWRI和 MERSI融合后的LST產(chǎn)品仍在Gap區(qū)域存在大量缺值,在這里采用Duan et al(2017)的補(bǔ)值法進(jìn)行Gap處理,即假設(shè)當(dāng)天云區(qū)下的LST值與前1天和后1天云區(qū)下的LST值相當(dāng),因此根據(jù)1 km云區(qū)Gap位置,需讀取前后各1天融合后的1 km LST值進(jìn)行補(bǔ)值。如果前后2天均無LST值,仍作為缺值處理。

2.4 Duan et al(2017)LST降尺度法

為了對比本文25 km FY-3D/MWRI LST降尺度效果,這里引入Duan et al(2017)基于MODIS/Aqua 1 km LST對AMSR-E 25 km LST進(jìn)行LST降尺度的方法(以下簡稱Duan法),其核心思想是把AMSR-E LST 25 km 分辨率的像元分為三類:完全晴空、完全云覆蓋、部分云覆蓋進(jìn)行融合,針對每類像元分別進(jìn)行處理:

(1)完全晴空像元,融合后的該像元在1 km分辨率的LST仍為原始MODIS LST值。

(2)完全云覆蓋像元:利用該25 km像元內(nèi)1 km 分辨率的DEM值進(jìn)行海拔高度差訂正(海拔每下降1 km,LST降低6.5 K)。

(3)部分云覆蓋混合像元:基于混合像元分解法原理,把25 km混合像元LST近似看成晴空子像元LST及面積比例和云覆蓋子像元LST及面積比例的線性加權(quán)和,根據(jù)該混合像元LST、晴空子像元平均LST及面積比例、云覆蓋子像元面積比例,即可求取所有云覆蓋子像元平均LST,然后根據(jù)子像元的海拔高度進(jìn)行海拔訂正,即可得到每個(gè)云覆蓋子像元LST。

3 結(jié)果分析

3.1 GWR模型效果分析

對MWRI 25 km LST 模擬分別選用NDVI、NDBI、DEM和NDVI、DEM的GWR回歸模型分析結(jié)果如圖5所示:針對中國區(qū)域,在NDVI和DEM基礎(chǔ)上,加入NDBI后,25 km LST模擬效果總體決定系數(shù)R2得到提高;而所有局地回歸模型平均決定系數(shù)R2增加。

圖5 基于FY-3D/MWRI LST的不同回歸變量建立的GWR模型LST模擬結(jié)果比較

對2019年7月14—16日、2019年10月14—16日、2020年4月16—18日FY-3D/MWRI LST建立的GWR模型擬合效果(圖6)可以看出,GWR模型對MWRI LST模擬的R2為0.77～0.93,均方根誤差(RMSE)為1.2～3.2 K;與傳統(tǒng)的TsHARP模型(Agam et al,2007)相比較,GWR模型的決定系數(shù)R2顯著提高(平均值從0.36提高到0.85),而RMSE明顯降低(平均值從4.3 K降低到2.2 K)。該結(jié)果表明GWR模型比TsHARP模型可更好地預(yù)測大范圍區(qū)域LST。

圖6 GWR模型與TsHARP模型對FY-3D/MWRI LST模擬結(jié)果的比較

如圖7所示為2020年4月17日FY-3D 25 km MWRI LST GWR模型回歸系數(shù)1 km分辨率空間插值結(jié)果,各回歸系數(shù)存在明顯空間差異,而局地回歸模型模擬效果也存在較大空間差異,但局地決定系數(shù)R2一般在0.40以上,平均為0.568。

圖7 2020年4月17日FY-3D/MWRI 25 km LST GWR模型回歸系數(shù)空間插值1 km結(jié)果

3.2 GWR模型降尺度效果

如圖8所示為2019年10月15日FY-3D/MWRI 25 km LST降尺度到1 km前后LST對比,可以看出,基于GWR模型,25 km LST降尺度到1 km后,1 km LST空間分布與原始25 km LST具有良好一致性,空間細(xì)節(jié)增加明顯,還能凸顯出北京中心城市地區(qū)的高溫現(xiàn)象,而在原始25 km LST影像中這種高溫現(xiàn)象不明顯。

圖8 2019年10月15日FY-3D/MWRI 25 km LST降尺度到1 km前后LST對比

3.3 當(dāng)日LST融合效果

如圖9所示為2019年7月15日FY-3D/MWRI 25 km LST降尺度到1 km前后LST對比,可以看出,與原始MWRI LST相比,采用Duan法融合后的LST與本算法融合的LST空間分辨率明顯增加,且均具有較好的空間分布一致性,但Duan法出現(xiàn)明顯“斑塊”效應(yīng),且易出現(xiàn)局部低值現(xiàn)象(圖中藍(lán)色圓圈所示)。斑塊化的原因是由于采用海拔差進(jìn)行25 km LST降尺度訂正時(shí),由于不同尺度海拔存在明顯差異造成降尺度后的相鄰像元LST不連續(xù);局部低值則是采用的混合像元(包括晴空和云覆蓋區(qū))法所導(dǎo)致,例如25 km分辨率像元LST為295 K,晴空覆蓋率為0.90,如果晴空子像元LST為300 K,則云覆蓋區(qū)域LST僅為250 K,明顯偏低;而本文算法(GWR)則有效消除了“斑塊”效應(yīng)和局部低值問題。

注:圓圈為降尺度中出現(xiàn)的局部低值現(xiàn)象。

3.4 最終LST融合結(jié)果

在當(dāng)日FY-3D 1 km MERSI LST與25 km MWRI LST融合基礎(chǔ)上,進(jìn)行前后1天融合LST軌道間隙(Gap)區(qū)域補(bǔ)值,生成最終當(dāng)日LST融合結(jié)果,如圖10所示為2019年7月15日、2019年10月15日、2020年4月17日生成的當(dāng)日FY-3D/MERSI 1 km LST(圖10a)、與FY-3D/MWRI 25 km LST融合后生成的當(dāng)日1 km LST(圖10b)以及前后1天Gap LST補(bǔ)值結(jié)果(圖10c),可見當(dāng)日FY-3D 1 km MERSI LST有值率(有效值面積百分比)僅為22.4%～36.9%,與25 km MWRI LST融合后,LST有值率增加到59.3%～65.3%,經(jīng)過前后1天LST補(bǔ)值后,LST有值率增加到69.3%～80.7%,但仍遠(yuǎn)未達(dá)到全天候LST的融合比例(100%)。

圖10 典型日期下(a)FY-3D/MERSI 1 km LST,(b)當(dāng)日1 km融合結(jié)果,(c)軌道間隙前后1天補(bǔ)值結(jié)果及LST有效值占中國區(qū)域面積百分比

3.5 融合LST質(zhì)量分析

3.5.1 LST融合產(chǎn)品缺值分析

為了分析LST融合產(chǎn)品有值率明顯小于全天候比例(100%)原因,這里統(tǒng)計(jì)了FY-3D/MWRI衛(wèi)星軌道(FY-3D/MWRI Orbit)覆蓋中國區(qū)域百分比與25 km 分辨率MWRI LST有值率,以及1 km FY-3D/MERSI LST和時(shí)間相近的MODIS/Aqua LST有值率,如圖11所示。可以看出,FY-3D/MWRI衛(wèi)星軌道一般覆蓋中國區(qū)域64.0%～73.5%,而MWRI LST有值率則明顯偏低,尤其是冬季如2020年1月12—14日MWRI LST有值率僅為7.0%～11.1%,這也是本文未對冬季FY-3D/MWRI LST與MERSI LST進(jìn)行融合的主要原因;而其他季節(jié)MWRI LST有值率為38.2%～56.9%,較相應(yīng)MWRI軌道覆蓋率降低15.0%～27.6%;結(jié)合圖4、圖8和圖9 MWRI過境軌道LST有效值分析,缺值主要出現(xiàn)在青藏高原地區(qū)與部分低溫地區(qū),這與FY-3D/MWRI LST反演之前需要做冰雪和降雨等像元剔除原因密切相關(guān)(楊軍和董超華,2011),而且MWRI LST反演算法在冬天容易把大量低溫像元識(shí)別為冰雪,因此造成冬季MWRI LST大量缺值,而在其他季節(jié)也容易把青藏高原低溫地區(qū)識(shí)別為冰雪,從而造成青藏高原地區(qū)與低溫地區(qū)MWRI LST缺值明顯。

圖11 (a)FY-3D/MWRI衛(wèi)星軌道覆蓋面積百分比與25 km 分辨率MWRI LST有值率,(b)1 km FY-3D/MERSI LST和同時(shí)間MODIS/Aqua LST有值率對比

另外,FY-3D/MERSI LST的有值率僅為15.1%～38.9%,較相應(yīng)MODSI/Aqua LST有值率偏低1.8%～22.8%,其中冬季2020年1月12—14日偏低均在19%以上,而其他季節(jié)也明顯偏低,表明FY-3D/MERSI LST很可能過高地估計(jì)了云覆蓋率或冰雪覆蓋率(目前FY-3對冰雪地區(qū)未做LST反演),而在冬季尤為明顯。

綜上所述,FY-3D/MWRI 和FY-3D/MERSI LST均存在缺值率較為明顯問題,這是兩者LST融合產(chǎn)品有值率未能達(dá)到全天候比例(100%)的一個(gè)重要原因。

3.5.2 LST融合產(chǎn)品精度分析

在這里利用全天候1 km TPDC LST產(chǎn)品對1 km FY-3D Fusion LST(融合產(chǎn)品)分為晴空、云天和全天候條件下進(jìn)行產(chǎn)品交叉驗(yàn)證評估。為了減少產(chǎn)品之間的配準(zhǔn)誤差以及不同點(diǎn)擴(kuò)散函數(shù)導(dǎo)致的地理定位的不確定性,產(chǎn)品間交叉驗(yàn)證時(shí)使用最小一致空間支持法(MCSS)(Baret et al,2013),即像元值采用3×3像元中有效像元比例不低于5/9的情況下有效像元平均值?；贛CSS法,對2019年7月15日、2019年10月15日、2020年4月17日1 km LST融合產(chǎn)品評估結(jié)果如圖12所示,無論是晴空、云天還是全天候下白天1 km FY-3D/MERSI LST與MODIS/Aqua LST存在明顯空間相關(guān)性,其中晴空天氣下R2為0.718～0.914(平均為0.813),平均偏差(MBE)為-0.9～0.5 K,RMSE為2.6～4.9 K(平均為3.5 K);云天下R2較晴空有明顯下降,為0.457～0.628(平均為0.530),MBE為-1.6～0.9 K,RMSE為4.1～6.1 K(平均為5.1 K);而全天候則是晴天和云天綜合的結(jié)果,R2為0.503～0.787(平均為0.660),MBE為-1.5～0.8 K,RMSE為3.6～5.8 K(平均為4.8 K)。

注:深藍(lán)色線為LST 1∶1線,紅色線為線性回歸方程線。

除了利用站點(diǎn)和已有LST產(chǎn)品進(jìn)行評估外,與其他類似相關(guān)LST融合產(chǎn)品的精度相互比較也能一定程度看出該產(chǎn)品的質(zhì)量總體狀況(Duan et al,2017;Yoo et al,2020;Zhong et al,2021)。表1列出了近些年已有研究對基于微波傳感器(主要為AMSR-E/2)與光學(xué)傳感器(主要為MODIS Aqua)1 km LST白天融合產(chǎn)品精度評估結(jié)果。由于采用的不同的融合算法和評估方法,各融合LST產(chǎn)品精度變化較大,如決定系數(shù)R2變動(dòng)較大(0.51～0.98),全天候條件下Yoo et al(2020)的RMSE最小(1.9～3.3 K),而Xu et al(2019)的RMSE可達(dá)4.27～8.29 K;此外,晴空天氣下的LST精度一般好于云下,如Duan et al(2017)、Zhang et al (2020)和Zhong et al(2021),但晴空LST的RMSE一般小于4.2 K,云下的RMSE可達(dá)6 K以上。與這些已有研究結(jié)果相比,本文基于FY-3D的LST融合產(chǎn)品精度雖稍差于基于AMSR-E/2和MODIS/Aqua LST的融合產(chǎn)品,但無論是平均R2,還是平均RMSE,均在上述研究結(jié)果最大最小精度范圍之內(nèi),表明基于本文的融合方法,仍可以實(shí)現(xiàn)FY-3D MWRI與MERSI LST的有效融合。

表1 本文1 km LST 融合產(chǎn)品精度與近些年已有研究結(jié)果(白天)比較

4 結(jié)論與討論

4.1 結(jié) 論

1 km分辨率的全天候陸表溫度(LST)已經(jīng)成為全球與區(qū)域氣候變化領(lǐng)域研究的迫切需求,目前還沒有基于國產(chǎn)衛(wèi)星的1 km分辨率的全天候LST產(chǎn)品。FY-3D 中分辨率成像儀Ⅱ型(MERSI-Ⅱ)提供了1 km分辨率的晴空LST產(chǎn)品,而FY-3D微波成像儀(MWRI)能獲取不受云影響的LST產(chǎn)品,但空間分辨率較粗(25 km),因此可結(jié)合兩者優(yōu)勢開展全天候1 km分辨率LST的融合研究。本文基于局部參數(shù)估計(jì)的地理加權(quán)回歸(GWR)方法,選擇海拔、FY-3D 歸一化植被指數(shù)(NDVI)和歸一化建筑指數(shù)(NDBI)等變量開展25 km MWRI LST降尺度到1 km,并針對MWRI過境軌道之間間隙(Gap)無LST情況,利用前后1天融合后的1 km LST云覆蓋像元LST進(jìn)行補(bǔ)值,從而開展了基于FY-3D的全天候LST融合產(chǎn)品研究。研究結(jié)果表明:

與傳統(tǒng)的全局參數(shù)回歸方法TsHARP相比,基于GWR建立的回歸模型預(yù)測的MWRI LST的R2從0.363增加到0.85,均方根誤差從4.26 K降至2.17 K。基于該回歸模型對FY-3D/MWRI 25 km LST降尺度到1 km LST,可以有效解決傳統(tǒng)的基于海拔與混合像元分解的LST降尺度方法存在明顯“斑塊”效應(yīng)和局地溫度偏低等問題。對多個(gè)典型日期中國FY-3D MERSI和MWRI 1 km LST融合效果表明,原始1 km晴空的LST有值率為22.4%～36.9%,融合后的1 km LST有值率可達(dá)69.3%～80.7%;利用基于MODIS LST和GLDAS/CLDAS的全天候1 km LST產(chǎn)品對本文1 km 融合LST評估表明,基于FY-3D的融合LST的空間決定系數(shù)R2為0.457～0.914,RMSE為3.6～5.8 K,其中晴空2.6～4.9 K,云下4.1～6.1 K,接近或稍差于AMSR-E/2和MODIS/Aqua LST的融合產(chǎn)品精度,表明基于GWR法,可以實(shí)現(xiàn)FY-3D/MWRI與MERSI LST的有效融合。

4.2 討 論

本文基于GWR所獲得的融合LST產(chǎn)品最高有值率僅為80%左右,這遠(yuǎn)未達(dá)到100%的全天候LST要求。從前面分析可知,FY-3D 25 km MWRI LST在青藏高原地區(qū)和低溫地區(qū)缺值率較高,而FY-3D/MERSI LST在四季均存在缺值偏高情況,尤其是在冬季。因此,需要一方面改進(jìn)官網(wǎng)產(chǎn)品FY-3D 1 km MERSI LST的云檢測算法;另一方面還需大力改進(jìn)官網(wǎng)產(chǎn)品FY-3 LST在青藏高原地區(qū)和低溫地區(qū)缺值問題,這是由于這些地區(qū)可能判識(shí)為冰雪,而FY-3 MWRI LST和MERSI LST目前的算法均不支持對冰雪地區(qū)的LST反演,因此開展冰雪地區(qū)LST反演是未來提高FY-3 LST有值率的重要發(fā)展方向。

另外,基于FY-3D 的LST融合產(chǎn)品誤差較大的原因,主要有:

(1)FY-3D/MWRI LST誤差:從FY-3氣象衛(wèi)星數(shù)據(jù)服務(wù)官網(wǎng)獲取的FY-3D/MWRI LST誤差較大。雖然目前基于微波反演得到LST誤差較大,一般在3～6 K(Fily et al,2003;Holmes et al,2009;Zhou et al,2015;Shwetha and Kumar,2016),但本文的FY-3D/MWRI LST誤差一般在5 K以上,雖然可以利用FY-3D/MERSI LST對其進(jìn)行偏差訂正從而降低誤差,但如果MWRI LST與MERSI LST空間相關(guān)性并不好時(shí),基于線性回歸模型對25 km MWRI LST進(jìn)行訂正仍會(huì)造成較大誤差,從而造成融合后LST云下誤差較大,因此,還需改進(jìn)FY-3D/MWRI LST反演算法,使其平均誤差降至5 K 以下。

(2)FY-3D/MERSI LST誤差:從FY-3氣象衛(wèi)星數(shù)據(jù)服務(wù)官網(wǎng)獲取的FY-3D/MERSI LST誤差較大。從圖12晴空天氣下的分析可知,FY-3D/MERSI LST與MODIS/Aqua LST(全天候TPDC的晴空像元即為MODIS/Aqua LST)之間的R2為0.717～0.914,RMSE為2.6～4.9 K,這個(gè)產(chǎn)品精度均明顯低于Du et al(2021)基于改進(jìn)的分裂窗算法反演的FY-3D/MERSI LST 精度(R2為0.84～0.96,RMSE為1.57～2.52 K)以及Aveni and Blackett(2022)反演的FY-3D/MERSI LST精度(R2為0.92,RMSE為0.75 K)。因此,還需改進(jìn)FY-3氣象衛(wèi)星數(shù)據(jù)服務(wù)官網(wǎng)的FY-3D/MERSI LST反演算法,使其接近AQUA LST產(chǎn)品精度。

(3)融合算法誤差:本文基于GWR建立的局地回歸模型一個(gè)前提假設(shè)是MWRI LST與NDVI、NDBI、DEM、Lat、Lon存在明顯線性關(guān)系,但受其他因素影響,可能這種關(guān)系并不顯著,例如圖7中建立的局地回歸模型平均R2不到0.60,因此勢必會(huì)存在許多像元LST與回歸變量之間的線性相關(guān)性比較低,從而會(huì)造成LST預(yù)測以及降尺度產(chǎn)生明顯誤差。因此,還需不斷改進(jìn)融合算法來減小LST預(yù)測誤差。

致謝：感謝國家衛(wèi)星氣象中心董立新正研級(jí)高工給予的技術(shù)指導(dǎo)。