基于FY-3D衛星的微波與光學陸表溫度融合研究*

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

1 中國氣象局地球系統數值預報中心,北京 100081 2 災害天氣國家重點實驗室,北京 100081 3 南京信息工程大學遙感與測繪工程學院,南京 210044 4 國家衛星氣象中心,北京 100081 5 中國農業科學院農業資源與農業區劃研究所,北京 100081 6 湖南省氣象局,長沙 410118

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

引 言

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

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

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

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

1 試驗數據

1.1 研究數據

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

文中所涉及的中國地圖基于國家測繪地理信息局標準地圖服務網站下載的審圖號GS(2019)1822號的標準地圖制作,底圖無修改。

1.2 產品驗證數據

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

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

2 研究方法

2.1 地理加權回歸模型

地理加權回歸(GWR)方法是傳統標準回歸方法(如普通最小二乘方法)的拓展,該方法是進行局部而不是全局參數估算(Fotheringham et al,2002)。通過利用地理空間變化的回歸系數,來研究自變量與因變量之間的動態回歸關系。GWR模型可表示為:

(1)

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

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

(2)

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

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

(3)

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

(4)

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

2.2 LST降尺度變量選擇

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

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

(5)

由此,本研究選擇FY-3D旬歸一化植被指數NDVI、NDBI產品和DEM作為LST降尺度變量;由于地理經度Lon和緯度Lat在GWR模型中己經是輸入參數,因此不再于模型中顯示。

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

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

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

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明顯偏低,而且是系統性的偏低,平均偏差MBE(MWRI LST與MERSI LST之差平均值)小于-5 K,部分在-10 K以下,這與高浩等(2018)對FY-3C/MWRI LST評估顯示其LST明顯偏低結果相似,表明FY-3/MWRI LST存在明顯的低估。由于MWRI LST與MERSI LST存在明顯空間相關,決定系數R2一般在0.75以上,因此可以利用線性回歸模型基于FY-3D/MERSI LST對FY-3D/MWRI LST進行偏差訂正,在此基礎上開展兩者融合,主要步驟為:

注:深藍色線為LST 1∶1線,紅色線為線性回歸方程線。

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

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

經過上述步驟后,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降尺度

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

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

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

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

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

(6)

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

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

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

(7)

2.3.3 當日1 km LST融合

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

2.3.4 Gap補值

由于MWRI全球各掃描軌道之間尤其是中低緯度存在明顯空隙(Gap),且每天均不固定,如圖4所示,因此MWRI和 MERSI融合后的LST產品仍在Gap區域存在大量缺值,在這里采用Duan et al(2017)的補值法進行Gap處理,即假設當天云區下的LST值與前1天和后1天云區下的LST值相當,因此根據1 km云區Gap位置,需讀取前后各1天融合后的1 km LST值進行補值。如果前后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進行LST降尺度的方法(以下簡稱Duan法),其核心思想是把AMSR-E LST 25 km 分辨率的像元分為三類:完全晴空、完全云覆蓋、部分云覆蓋進行融合,針對每類像元分別進行處理:

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

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

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

3 結果分析

3.1 GWR模型效果分析

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

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

對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;與傳統的TsHARP模型(Agam et al,2007)相比較,GWR模型的決定系數R2顯著提高(平均值從0.36提高到0.85),而RMSE明顯降低(平均值從4.3 K降低到2.2 K)。該結果表明GWR模型比TsHARP模型可更好地預測大范圍區域LST。

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

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

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

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具有良好一致性,空間細節增加明顯,還能凸顯出北京中心城市地區的高溫現象,而在原始25 km LST影像中這種高溫現象不明顯。

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

3.3 當日LST融合效果

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

注:圓圈為降尺度中出現的局部低值現象。

3.4 最終LST融合結果

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

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

3.5 融合LST質量分析

3.5.1 LST融合產品缺值分析

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

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

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

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

3.5.2 LST融合產品精度分析

在這里利用全天候1 km TPDC LST產品對1 km FY-3D Fusion LST(融合產品)分為晴空、云天和全天候條件下進行產品交叉驗證評估。為了減少產品之間的配準誤差以及不同點擴散函數導致的地理定位的不確定性,產品間交叉驗證時使用最小一致空間支持法(MCSS)(Baret et al,2013),即像元值采用3×3像元中有效像元比例不低于5/9的情況下有效像元平均值?；贛CSS法,對2019年7月15日、2019年10月15日、2020年4月17日1 km LST融合產品評估結果如圖12所示,無論是晴空、云天還是全天候下白天1 km FY-3D/MERSI LST與MODIS/Aqua LST存在明顯空間相關性,其中晴空天氣下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);而全天候則是晴天和云天綜合的結果,R2為0.503～0.787(平均為0.660),MBE為-1.5～0.8 K,RMSE為3.6～5.8 K(平均為4.8 K)。

注:深藍色線為LST 1∶1線,紅色線為線性回歸方程線。

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

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

4 結論與討論

4.1 結 論

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

與傳統的全局參數回歸方法TsHARP相比,基于GWR建立的回歸模型預測的MWRI LST的R2從0.363增加到0.85,均方根誤差從4.26 K降至2.17 K?；谠摶貧w模型對FY-3D/MWRI 25 km LST降尺度到1 km LST,可以有效解決傳統的基于海拔與混合像元分解的LST降尺度方法存在明顯“斑塊”效應和局地溫度偏低等問題。對多個典型日期中國FY-3D MERSI和MWRI 1 km LST融合效果表明,原始1 km晴空的LST有值率為22.4%～36.9%,融合后的1 km LST有值率可達69.3%～80.7%;利用基于MODIS LST和GLDAS/CLDAS的全天候1 km LST產品對本文1 km 融合LST評估表明,基于FY-3D的融合LST的空間決定系數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的融合產品精度,表明基于GWR法,可以實現FY-3D/MWRI與MERSI LST的有效融合。

4.2 討 論

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

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

(1)FY-3D/MWRI LST誤差:從FY-3氣象衛星數據服務官網獲取的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對其進行偏差訂正從而降低誤差,但如果MWRI LST與MERSI LST空間相關性并不好時,基于線性回歸模型對25 km MWRI LST進行訂正仍會造成較大誤差,從而造成融合后LST云下誤差較大,因此,還需改進FY-3D/MWRI LST反演算法,使其平均誤差降至5 K 以下。

(2)FY-3D/MERSI LST誤差:從FY-3氣象衛星數據服務官網獲取的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,這個產品精度均明顯低于Du et al(2021)基于改進的分裂窗算法反演的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)。因此,還需改進FY-3氣象衛星數據服務官網的FY-3D/MERSI LST反演算法,使其接近AQUA LST產品精度。

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

致謝：感謝國家衛星氣象中心董立新正研級高工給予的技術指導。