FY-3D/MWRI L1B亮溫LST反演與降尺度研究

朱瑜馨，吳門新，鮑艷松，李鑫川，張錦宗

(1.淮陰師范學院城市與環境學院，淮安 223300；2.國家氣象中心，北京 100081；3.南京信息工程大學大氣物理學院，南京 210044)

0 引言

陸表溫度(land surface temperature,LST)是指陸地表面最上層的熱力學溫度，是研究地表水熱平衡和全球氣候變化的重要參數。精確的陸表溫度產品是陸面水文模型、數字天氣及氣候預報模式的輸入參數，也是利用微波遙感技術進行土壤水分反演的一個關鍵輸入參量[1]。目前，陸表溫度數據的獲取方式主要包括：地面氣象站點監測與遙感反演。地面氣象站點監測數據精度高，但時空不連續；基于遙感數據的LST反演可以獲取大范圍時空連續的數據，隨著遙感技術及計算機技術的發展，成為主要的LST獲取方式?；诳梢姽饧t外遙感數據的LST反演，空間分辨率高，精度相對較高[2]，常用的方法有單通道算法、劈窗算法、溫度與發射率分離算法等[3]，如MODIS(Moderate Resolution Imaging Spectroradiometer)，FY-2/3 VIRR(FengYun-2/3 Visible Infrared Scanning Radiometer)，SEVIRI(Spinning Enhanced Visible and InfraRed Imager)產品等。但基于可見光紅外遙感數據的LST反演需要晴空條件以及儀器的無誤差條件[4]，因此，受云、大氣等的影響較大，存在大量的缺失像元或無效像元，不能進行陸表溫度的全天候監測?；谖⒉ㄝ椛溆嫷牧翜財祿囱軱ST，具有多極化及高時間分辨率等特點[1]，受大氣、云等環境影響相對較小，且能全天候、全天時監測，但空間分辨率較低。二者各有優缺點，又相互補充，為充分利用二者間的互補性進行LST反演并降尺度，提供了可行性。

基于被動微波遙感的地表溫度反演研究起始于20世紀80年代末期[5]。如利用SSM/I數據，AMSR數據建模進行陸表溫度反演[6-7]，模型精度在1.5～3 K之間；以及以AMSR-E數據為基礎的陸表溫度反演等[8-11]。與可見光紅外數據相比，極少受到云或水蒸汽的影響，但空間分辨率遠低于可見光紅外反演LST，且傳感器缺陷及軌道掃描間隙也導致了空間的不連續性。

根據Planck黑體輻射原理和Jeans近似，一般地物的微波輻射亮溫與其真實溫度存在簡單的線性關系[12]，因此，本文利用FY-3D MWRI日亮溫數據與FY-3C VIRR LST日數據，建立單通道、多極化回歸模型，選擇相關系數最高的頻率的水平和垂直極化數據建立二元一次回歸模型，進行MWRI LST反演，并將反演后的LST數據通過層次貝葉斯模型與VIRR LST融合進行降尺度，最終得到1 km空間連續的LST日數據。

目前，我國FY-3系列微波輻射計亮溫反演LST數據較少，僅僅有FY-3D MWRI 25km LST產品，但空間不連續，軌道間存在大量的缺失像元。本文基于FY-3C 1 km VIRR LST產品和FY-3D 10km MWRI L1B亮溫數據，充分利用二者在空間完整性及時空分辨率等方面的互補性，進行LST反演與降尺度，具有較強的應用價值。本文以18°～54°N，73°～135°E區域為樣例區，進行FY-3D MWRI L1B亮溫LST反演與降尺度研究，技術路線如圖1。

圖1 技術路線Fig.1 Technology roadmap

1 數據源及其預處理

1.1 數據源

本文用到的遙感數據為2020年2月1日FY-3D MWRI L1B降軌和升軌日數據、FY-3C VIRR LST日數據、MODIS LST日數據。

FY-3D MWRI的L1B降軌和升軌全球日數據，空間分辨率10 km，MWRI在10.65～89 GHz頻段內設置5個頻率：10.65 GHz，18.70 GHz，23.80 GHz，36.50 GHz和89.00 GHz，每個頻率有垂直和水平極化模式。降軌數據自東向西通過樣例區的時間分別為：15∶53，17∶35，19∶16和20∶58；升軌數據經過樣例區的軌道數據時間分別為：03∶12，04∶54，06∶35和08∶16。數據存儲格式為HDF(hierarchical data format)。

FY-3C VIRR LST日數據，空間分辨率1 km，投影類型為Hammer，10°×10°分塊，經過樣例區的經緯度編號如表1。數據存儲格式為HDF。為了盡量降低FY-3D MWRI亮溫數據與MODIS LST 日數據之間的時間差異，本文采用MYD11A1的day LST來驗證FY-3D MWRI反演LST降軌和升軌降尺度數據，空間分辨率為1 km，投影類型為正弦。經過樣例區的行列號如表2。數據存儲格式為HDF。

表1 FY-3C VIRR LST 樣例區經緯度編號及相對位置Tab.1 Longitude and latitude number and relative position of sample area of FY-3C VIRR LST

表2 MYD11A1 LST 樣例區行列號及相對位置Tab.2 Tile and relative position of sample area of MYD11A1 LST

1.2 數據預處理

1)投影與轉換。對于風云3D MWRI亮溫數據，利用其自帶的經緯度圖層構建GLT(geometric lookup file)，進行幾何校正，輸出WGS-84地理坐標系統。對FY-3C VIRR LST，首先進行Hammar投影定義，然后轉換為WGS-84地理坐標系統。對于MODIS LST，通過MRT工具直接進行tile數據的拼接與投影轉換。

2)軌道拼接和裁剪并進行區域裁剪，得到樣例區范圍的數據。

3)根據數據提供的有效范圍及質量控制層進行數據的質量控制及異常值剔除。

根據MODIS LST效值范圍(7 500,65 535)剔除異常值，再根據質量控制層選擇質量標識為0的質量好的像元。VIRR LST是經過了云處理后的數據。根據MWRI BT數據的有效范圍(-32 766,32 767)剔除異常值。

4)亮溫和陸表溫度值的計算。根據公式(1)—(3)，分別計算FY-3D MWRI亮溫、FY-3C VIRR LST和MODIS LST，即

BTMWRI=0.01·DN+327.68，

(1)

LSTVIRR=0.1·DN，

(2)

LSTMODIS=0.02·DN，

(3)

在進行LST反演時，由于FY-3D MWRI亮溫與FY-3C VIRR LST存在空間分辨率差異，因此本文采用最近鄰插值算法進行重采樣，將FY-3C VIRR LST升尺度至10 km。

2 偏差校正

因為統計回歸模型和層次貝葉斯融合降尺度模型均不能降低數據的系統誤差和隨機誤差，因此，文章首先以MODIS LST為參考值，進行VIRR LST偏差校正。偏差校正前VIRR LST與MODIS LST散點圖，如圖2。

圖2 偏差校正前VIRR LST與MODIS LST散點圖Fig.2 Scatter plot of VIRR LST against MODIS LST before bias correction

VIRR LST與MODIS LST之間的相關系數為0.94，平均偏差為-10.82 K；誤差標準差為5.25 K，均方根誤差為12.03 K。

以VIRR LST為自變量，MODIS LST為因變量，建立回歸模型(公式(4))，R2=0883 2，P=0，<0.5，說明回歸模型成立，即

MODISLST=-1.490 3+1.045 1·VIRRLST。

(4)

將VIRR LST帶入回歸模型，完成偏差校正，校正后的VIRR LST與MODIS LST之間的散點圖如圖3。二者之間的平均偏差為4.67e-13 K；誤差標準差為5.21 K；均方根誤差為5.21 K。偏差校正后，平均偏差和均方根誤差明顯降低，散點圖分布更加對稱，說明校正效果明顯。校正后的VIRR LST如圖4。

圖3 偏差校正后VIRR LST與MODIS LST散點圖Fig.3 Scatter plot of VIRR LST against MODIS LST after bias correction

圖4 校正后的VIRR LSTFig.4 Spatial distribution of bias corrected VIRR LST

3 陸表溫度反演

反演算法采用單頻率水平極化和垂直極化的二元線性統計回歸模型：直接建立FY-3D MWRI L1B亮溫為自變量，FY-3C VIRR LST為因變量的統計回歸關系，通過求解方程系數得到最小二乘解。FY-3C VIRR LST，分別與FY-3D MWRI L1B亮溫數據5個頻率的水平和垂直極化數據進行回歸分析，選擇相關系數最大的頻率的水平和垂直極化數據建立地表溫度反演模型，公式為：

LSTVIRR=ao+a1·BTH+a2·BTV，

(5)

式中：LSTVIRR為FY-3C VIRR LST；BTH為FY-3D MWRI L1B水平極化亮溫；BTV為FY-3D MWRI L1B垂直極化亮溫。

3.1 模型構建

FY-3D MWRI降軌亮溫數據與FY-3C VIRR LST的相關系數如表3。根據相關系數，選擇36.50 GHz頻率的水平和垂直極化數據與FY-3C VIRR LST建立回歸模型，公式為：

表3 FY-3D MWRI降軌數據與FY-3C VIRR LST的相關系數Tab.3 Correlation coefficient between FY-3D MWRI orbit dropping data and FY-3C VIRR LST

LSTFY-D=118.631 3+0.866 1·BTMWRIH-0.211 3·BTMWRIV，

(6)

FY-3D MWRI升軌亮溫數據與FY-3C VIRR LST的相關系數如表4。

表4 FY-3D MWRI升軌數據與FY-3C VIRR LST的相關系數Tab.4 Correlation coefficient between FY-3D MWRI orbit lifting data and FY-3C VIRR LST

根據相關系數，選擇36.50 GHz頻率的水平和垂直極化數據與FY-3C VIRR LST建立回歸模型，即

LSTFY_R=124.091 9+0.818 6·BTMWRIH-0.217 7·BTMWRIV。

(7)

3.2 反演與驗證

根據公式(6)和(7)反演得到FY-3D降軌和升軌2020年2月1日10 km 空間分辨率LST。以MYD11A1 Day LST為參考值，分別計算FY LST與MYD LST之間的平均偏差、誤差標準差、均方根誤差及相關系數，對反演得到的升軌和降軌數據進行驗證，驗證結果如表5。從驗證結果看，升軌數據與VIRR LST吻合度稍高，但升軌和降軌數據與VIRR LST的吻合度差異不大。反演結果見圖5。

表5 FY LST與MYD11A1 Day LST驗證結果Tab.5 Validation results between FY-LST and MYD11A1 LST (K)

(a)FY-3D降軌反演LST(b)FY-3D升軌反演LST

圖6為MODIS LST與FY LST之間的散點圖。升軌數據中，絕對偏差在1 K內的點對數為2 350；絕對偏差在1～2之間的點對數為2 096；絕對偏差在2～3之間的點對數為1 986；絕對偏差在3～4之間的點對數為2 046；絕對偏差大于4 K的點對數為12 270；最小絕對偏差為5.11e-05 K；最大絕對偏差為29.81 K。整體來看，匹配點基本對稱分布在1∶1線兩側。降軌數據中，絕對偏差在1 K內的點對數為2 301；絕對偏差在1～2之間的點對數為2 265；絕對偏差在2～3之間的點對數為2 155；絕對偏差在3～4之間的點對數為1 999；絕對偏差大于4 K的點對數為16 222；最小絕對偏差為0.0 015 K；最大絕對偏差為39.50 K。升軌數據與降軌數據相比，升軌數據更貼近1∶1線，說明升軌數據偏差較降軌數據小。

(a)FY降軌LST與MODISLST散點圖(b)FY升軌LST與MODISLST散點圖

4 陸表溫度降尺度

在層次貝葉斯框架下，通過建立參數的層次嵌套結構進行偏差校正后的1 km空間分辨率VIRR LST與10 km空間分辨率FY-3D MWRI 反演LST的融合，得到1 km空間分辨率、較高精度的、時空完整的、局部細節信息豐富的融合LST，實現降尺度。

層次貝葉斯將隨機變量的聯合分布分解為一系列條件概率乘積的模式，核心算法為：通過構建一個馬爾科夫鏈，迭代抽樣，每一步迭代的結果都依賴于前一次迭代抽樣的結果，最后收斂至平穩的靜態分布——后驗分布。此方法可以概述為：某一時空過程為Y，在該時空過程支配下的觀測數據為Z，與時空過程Y相關的一系列參數及與觀測數據模型相關的一系列參數組成的參數集為θ，在層次貝葉斯框架下，將Z，Y，θ的聯合概率分布表示為一組條件概率分布的乘積：[Z,Y,θ]=[Z|Y,θ][Y|θ][θ]。根據貝葉斯定理，結合觀測數據，更新、修正未知變量的先驗，得到其后驗分布[13-14]，具體公式為：

[Y,θ|Z]∝[Z|Y,θ][Y|θ][θ]。

(8)

層次的嵌套結構表現為利用最近鄰距離插值得到空間完整的10 km LST作為層次貝葉斯過程參數的先驗均值。

4.1 模型構建

層次貝葉斯分3個階段建模：數據模型、過程模型和參數模型。

1)數據模型。FY-3C VIRR LST與FY-3D 反演LST的觀測誤差模型(公式(9)和(10))。

3)尺度無縫轉換。1 km VIRR LST與10 km MWRI LST間存在尺度差異，通過建立嵌套的參數正態分布，實現無縫尺度轉換(公式(11)和(12))[13-14]。

(9)

(10)

(11)

(12)

4.2 降尺度結果

丟棄前1 000次不穩定迭代，經過500次迭代抽樣，單鏈收斂。表6和表7分別為部分節點抽樣統計。

表6 部分節點抽樣統計(升軌數據)Tab.6 Sample statistics of partial nodes (orbit lifting)

表7 部分節點抽樣統計(降軌數據)Tab.7 Sample statistics of part nodes (orbit dropping)

圖7為降軌和升軌LST融合降尺度結果。與10 km空間分辨率反演LST(圖4)對比分析，降尺度數據空間完整，局部細節信息增多。由于10 km空間分辨率原始亮溫數據海洋像元標識因軌道縫隙部分缺失，導致圖6中融合降尺度數據的東南方向未能做海洋像元屏蔽。

(a)降軌降尺度LST(b)升軌降尺度LST

圖8為降尺度LST與MODISLST散點圖。降尺度為1 km的降軌數據中，絕對偏差在1 K內的點對數為602 484；絕對偏差在1～2 K之間的點對數為575 318；絕對偏差在2～3 K之間的點對數為502 309；絕對偏差在3～4 K之間的點對數為425 718；絕對偏差大于4 K的點對數為1 535 756；最小絕對偏差為0 K；最大絕對偏差為39.47 K。降尺度為1 km的升軌數據中，絕對偏差在1 K內的點對數為582 063；絕對偏差在1～2之間的點對數為563 119；絕對偏差在2～3之間的點對數為501 844；絕對偏差在3～4之間的點對數為433 734；絕對偏差大于4 K的點對數為1 589 183；最小絕對偏差為0 K；最大絕對偏差為35.42 K。整體來看，匹配點基本對稱分布在1:1線兩側。升軌數據與降軌數據相比，升軌數據更貼近1:1線，說明升軌數據偏差較降軌數據小。

(a)降軌數據與MODIS散點圖(b)升軌數據與MODIS散點圖

表7為降尺度數據與MODISLST驗證結果，由表可知，降軌和升軌融合降尺度LST精度相似。平均偏差，升軌降尺度LST更小。對比分析VIRRLST與MODISLST驗證結果，可以看出，融合降尺度數據精度高于原始VIRRLST，接近偏差校正后的VIRRLST，說明融合降尺度不能降低數據的誤差，融合降尺度前的偏差校正是必要的。

表7 FY LST降尺度與MYD11A1 day LST驗證結果Tab.7 Validation results between FY downscaling LST and MYD11A1 day LST (K)

層次貝葉斯的優勢除了能保留被融合數據的局部細節信息外，還可以通過融合后驗標準差反映融合降尺度結果的不確定性。圖9顯示了融合降尺度結果的融合后驗標準差。對比分析VIRRLST和10 km空間分辨率FYLST，可以看出在VIRRLST和FYLST均有有效值的區域，融合后驗標準差較小，而在VIRRLST和FYLST缺值區域融合后驗標準差偏大。

(a)降軌數據融合后驗標準差(b)升軌數據融合后驗標準差

5 結論與討論

根據LST與微波亮溫之間存在的線性關系，建立FY-3D MWRI 單通道水平和垂直極化二元線性回歸模型進行LST反演，并對反演LST通過層次貝葉斯融合模型進行降尺度，得到了空間局部細節信息豐富的高空間分辨率LST數據，以MYD11A1 day LST數據為參考數據進行驗證。反演LST平均偏差分別為-1.28 K(降軌)和-0.81 K(升軌)；降尺度數據平均偏差分別為0.50 K(降軌)和0.25 K(升軌)。本研究結果表明，單通道水平和垂直極化二元線性模型反演LST簡單可行，層次貝葉斯融合降尺度模型保持了原有細尺度數據的局部細節信息，通過融合后驗標準差可以清楚表達融合結果的不確定性，為空間降尺度提供了新的思路與方法。

在LST反演與降尺度過程中，考慮到直接建立回歸模型得到1 km反演LST需要其他地表參數的遙感產品作為輔助，會引入更多的不確定性，同時受FY-3D MWRI L1B軌道間隙及輔助參數空間完整性的影響，不能得到空間完整的降尺度LST。因此，文章沒有直接建立回歸模型得到1 km反演LST，采用的策略為先將FY-3C VIRR LST升尺度至10 km來匹配FY-3D MWRI L1B數據，建立回歸模型得到粗尺度反演LST，再通過層次貝葉斯融合降尺度模型實現降尺度，這樣既可以得到空間完整的降尺度LST，且降尺度LST又同時保持了細尺度原始數據的細節信息，信息量大。本文在融合降尺度過程中僅僅采用了參數鏈接的先驗均值方法，還有進一步改進的空間，在今后的研究中將側重研究LST空間變異過程模擬，并將此過程模型代替參數先驗均值嵌入層次貝葉斯融合降尺度模型中，以解決在兩種數據均缺值區域融合不確定性偏大的問題。