[發明專利]一種地表溫度遙感產品不確定度的逐像元估算方法

有效

• 申請號: 202210171560.8
• 申請日: 2022-02-24
• 公開（公告）號: CN114544003B
• 公開（公告）日: 2022-11-08
• 發明（設計）人: 段四波;桂陽;李召良;張霄羽
• 申請（專利權）人: 中國農業科學院農業資源與農業區劃研究所;山西大學
• 主分類號: G01J5/48
• 分類號: G01J5/48;G06F17/15
• 代理公司: 北京恒創益佳知識產權代理事務所(普通合伙) 11556
• 代理人: 付金豹

1.一種地表溫度遙感產品不確定度的逐像元估算方法，其特征在于，包括以下步驟：

步驟1：獲取MODIS數據和數據預處理；

數據預處理過程包括：

(1)利用MODIS地理定位產品MYD03和MODIS再投影工具MRTSwath將MODIS各產品轉換為具有地理坐標的TIFF格式數據；

(2)將MYD11_L2產品中地表發射率和觀測天頂角的計數值分別轉換為地表發射率和觀測天頂角：

VZA＝0.5×DN_VZA (1)

ε₃₁＝0.002×DN_ε31+0.49 (2)

ε₃₂＝0.002×DN_ε32+0.49 (3)

式中，DN_VZA為觀測天頂角的計數值，DN_ε31和DN_ε32分別為MODIS第31和32通道地表發射率的計數值，VZA為觀測天頂角，ε₃₁和ε₃₂分別為MODIS第31和32通道的地表發射率；

(3)將MYD05_L2產品中大氣水汽總量的計數值轉換為大氣水汽總量：

TCW＝0.0010000000474974513×DN_TCW (4)

式中，DN_TCW為大氣水汽總量的計數值，TCW為大氣水汽總量。

(4)將MYD021KM產品中星上亮溫的計數值轉換為星上亮溫：

T₃₁＝1304.413871/ln(1+729.541636/DN_T31) (5)

T₃₂＝1196.978785/ln(1+474.684780/DN_T32) (6)

式中，T₃₁和T₃₂分別為MODIS第31和32通道的星上亮溫，DN_T31和DN_T32分別為MODIS第31和32通道星上亮溫的計數值；

步驟2：建立模擬數據集；

從ASTER地物波譜庫中挑選4條水體和冰雪光譜曲線、3條植被光譜曲線、41條土壤光譜曲線用來表征地表發射率變化；將各個地物光譜曲線與MODIS光譜響應函數卷積得到MODIS第31和32通道的地表發射率；

從TIGR數據庫中挑選946條晴空大氣廓線用來表征大氣狀況變化，每條廓線包括了40層的溫度、濕度、壓強和高程信息；設置6個觀測天頂角：0°、33.56°、44.42°、51.32°、56.25°、60°；將這些大氣廓線輸入到大氣輻射傳輸模型MODTRAN5.2，結合MODIS光譜響應函數，模擬得到每個觀測天頂角下的MODIS第31和32通道的大氣參數：大氣透過率、大氣上行輻射和大氣下行輻射，并計算得到每條大氣廓線的大氣水汽總量；

地表溫度根據每條廓線的最底層空氣溫度T₀進行設置，當T₀大于或等于285K時，地表溫度設置范圍為[T₀-5K,T₀+15K]，步長為5K；當T₀小于285K時，地表溫度設置范圍為[T₀-5K,T₀+5K]，步長為5K；

將模擬的大氣參數、地表溫度和地表發射率輸入到熱紅外輻射傳輸方程得到不同情況下的星上輻亮度，再根據溫度與輻亮度之間建立的查找表，線性插值得到不同大氣狀況和地表條件下的MODIS第31和32通道的星上亮溫；模擬的MODIS第31和32通道星上亮溫、第31和32通道地表發射率、地表溫度、觀測天頂角、大氣水汽總量組成模擬數據集；

步驟3：利用分裂窗地表溫度反演模型反演地表溫度；

分裂窗地表溫度反演模型表示為：

式中，T_s是反演的地表溫度，T₃₁和T₃₂是MODIS第31和32通道的星上亮溫，ε是MODIS第31和32通道地表發射率的均值，即ε＝(ε₃₁+ε₃₂)/2，Δε是MODIS第31和32通道地表發射率的差值，即Δε＝ε₃₁-ε₃₂，A₀、A₁、A₂、A₃、A₄、A₅和A₆為分裂窗模型系數；

利用觀測天頂角、大氣水汽總量、地表溫度作為索引，將模擬數據集劃分成不同的模擬數據子集，通過最小二乘法擬合確定每組模擬數據子集的分裂窗模型系數A₀、A₁、A₂、A₃、A₄、A₅和A₆；其中觀測天頂角劃分6個值：0°、33.56°、44.42°、51.32°、56.25°、60°；大氣水汽總量劃分6個子區間：0-1.5g/cm²、1-2.5g/cm²、2-3.5g/cm²、3-4.5g/cm²、4-5.5g/cm²、5-6.5g/cm²；地表溫度劃分5個子區間：280K、275-295K、290-310K、305-325K、320K；基于分裂窗模型系數，結合獲取的MODIS第31和32通道的地表發射率、第31和32通道的星上亮溫、觀測天頂角、大氣水汽總量，通過分裂窗模型反演得到地表溫度遙感產品；

步驟4：構建地表溫度遙感產品的不確定度估算模型

(1)地表溫度反演模型誤差引起的地表溫度不確定度S_a

用分裂窗模型反演的地表溫度和模擬數據中的真實地表溫度的均方根誤差表示地表溫度反演模型誤差引起的地表溫度不確定度S_a：

式中，E[]代表求算術平均值，T_s,ret為分裂窗模型反演的地表溫度，T_s,act為模擬數據中的真實地表溫度；

(2)傳感器噪聲引起的地表溫度不確定度S_n

通過分裂窗模型中T_s對T₃₁和T₃₂進行求導，估算傳感器噪聲引起的地表溫度不確定度S_n：

式中，和分別代表分裂窗模型中T_s對T₃₁和T₃₂求導，σ_T31和σ_T32分別表示MODIS第31和32通道的傳感器噪聲；

(3)地表發射率誤差引起的地表溫度不確定度S_e

通過分裂窗模型中T_s對ε₃₁和ε₃₂進行求導，估算地表發射率誤差引起的地表溫度不確定度S_e：

式中，和分別代表分裂窗模型中T_s對ε₃₁和ε₃₂求導，σ_ε31和σ_ε32分別表示MODIS第31和32通道的地表發射率誤差；

(4)大氣水汽總量誤差引起的地表溫度不確定度S_w

將使用正確的大氣水汽總量分組對應的模型系數反演的地表溫度與使用大氣水汽總量分組發生變化后的模型系數反演的地表溫度的均方根誤差表示大氣水汽總量誤差引起的地表溫度不確定度S_w：

式中，E[]代表求算術平均值，T_s,corr為使用正確的大氣水汽總量分組對應的模型系數反演的地表溫度，T_s,incorr為使用大氣水汽總量分組發生變化后的模型系數反演的地表溫度；

(5)地表溫度遙感產品總不確定度S_GSW

根據誤差傳遞理論，將地表溫度遙感產品不確定度的各分量合成總不確定度S_GSW：

步驟5：參數化地表溫度遙感產品的不確定度分量；

對地表溫度遙感產品的不確定度分量進行參數化，建立地表溫度遙感產品的不確定度參數化模型，便于逐像元估算地表溫度遙感產品的總不確定度；

(1)從式(8)和(11)可以看出，基于模擬數據集通過統計分析可以估算得到S_a和S_w的值，其值是常量；

(2)MODIS第31和32通道的傳感器噪聲相等，即σ_T31＝σ_T32＝σ_T，將式(9)展開、合并及優化后，確定S_n表達式為：

式中，B₀、B₁和B₂是S_n 模型系數。

(3)MODIS第31和32通道的地表發射率誤差相等，即σ_ε31＝σ_ε32＝σ_ε，將式(10)展開、合并及優化后，確定S_e表達式為：

式中，T是MODIS第31和32通道的星上亮溫均值，即T＝(T₃₁+T₃₂)/2；ΔT是MODIS第31和32通道的星上亮溫差值，即ΔT＝T₃₁-T₃₂；C₀、C₁、C₂和C₃是S_e 模型系數；

(4)基于式(12)，結合S_a、S_n、S_e和S_w的表達式，確定地表溫度遙感產品總不確定度S_GSW表達式為：

式中，D₀、D₁、D₂、D₃、D₄和D₅是地表溫度遙感產品總不確定度S_GSW模型系數；

步驟6：逐像元估算地表溫度遙感產品總不確定度；利用步驟2構建的模擬數據集，將觀測天頂角、大氣水汽總量、地表溫度作為索引，通過最小二乘法擬合確定每組模擬數據子集的地表溫度遙感產品總不確定度模型系數D₀、D₁、D₂、D₃、D₄和D₅；任意觀測天頂角下的地表溫度遙感產品總不確定度估算模型系數可以通過6個預設的觀測天頂角下的模型系數線性插值得到；

根據獲取的MODIS觀測天頂角和大氣水汽總量以及反演的地表溫度，逐像元確定地表溫度遙感產品總不確定度估算模型系數；將獲取的MODIS第31和32通道的地表發射率、第31和32通道的星上亮溫輸入到地表溫度遙感產品總不確定度估算模型，逐像元估算地表溫度遙感產品總不確定度。

2.根據權利要求1所述的方法，其特征在于，所述步驟1中，從美國宇航局NASAEARTHDATA網站下載搭載在Aqua衛星上的MODIS系列產品，包括MYD11_L2、MYD05_L2、MYD021KM產品；其中，MYD11_L2產品用于獲取地表發射率和觀測天頂角，MYD05_L2用于獲取大氣水汽總量，MYD021KM用于獲取星上亮溫。