劈窗演算法最初是為反演海面溫度開發的，具體地說是針對NOAA/AVHRR的4和5通道設計的，後來也被用來反演地表溫度，這種演算法較成熟，精度也高。劈窗演算法以地表熱輻射傳導方程為基礎，利用10~13μm 大氣視窗內，兩個相鄰熱紅外通道(一般為10.5～11.5μm、11.5～12.5μm)對大氣吸收作用的不同，通過兩個通道測量值的各種組合來剔除大氣的影響，進行大氣和地表比輻射率的修正。表示式為：

T _S= T ₄+ A (T ₄- T ₅) + B

其中：T _S為地表真實溫度，T ₄和 T ₅分別為AVHRR的4和5通道，A和B為常量。

      AVHRR 的通道4(10.15~11.13μm) 和通道5 (11.15~12.15μm) 恰與MODIS 的第31 波段 (10.178~ 11.128μm) 和32 波段 (11.177~ 12.127μm ) 的中心波長相對應, 可將MODIS 的31、32 波段資料, 用於劈窗演算法進行地表溫度計算。很多學者對這個演算法進行了推演，得到很多新的演算法，如覃志豪、毛克彪等人。本文就是使用其他學者推演的演算法。

      利用MODIS資料劈窗演算法反演海表溫度技術流程如下圖：

圖：技術流程圖

注：（1）按照本流程反演出來的結果是SST。陸地上的值可以視為無效值，若要得到正確的陸表溫度，需要加入海陸分離的步驟，以及城鎮和自然表面的比輻射率計算。

（2）MODIS資料下載：nasa官網：http://modis.gsfc.nasa.gov/.

下面詳細介紹處理流程，本操作適合ENVI5.3及以上版本

第一步：讀取原始資料

    通過File->Open As->Scientific Formats->HDF4，選擇原始的hdf資料，點選OK，選擇熱紅外原始DN畫素值資料開啟，如下圖所示。

   開啟的是熱紅外原始資料集，第20-36波段，共16個波段，分別是：20、21、22、23、24、25、27、28、29、30、31、32、33、34、35、36波段。

圖：Dataset Brower 

第二步，熱紅外B31、B32輻射亮度定標

     在上一步中的Dataset Brower，找到相應的輻射亮度值定標係數：

     radiance_scales，和radiance_offset這兩項引數代表波段的增益和偏移量，是輻射定標的係數。比如要計算31波段的輻射亮度，讀取到scales為0.00084002，offsets為1577.34，帶入MODIS輻射定標的公式：Radiance=scales*（DN-offsets），即可以得到該波段的輻射亮度。

      開啟/Band Algebra/Band Math工具，輸入公式：0.00084002*(B31-1577.334)，（B31是第31波段的DN值），點選OK，選擇第31波段資料為B31，設定路徑和檔名，點選OK。得到的結果就是31波段的輻射亮度。

     同樣的方法得到32波段的輻射亮度，公式為：0.00072970*(B32-1658.22)。

第三步：幾何校正

      MODIS1b資料是hdf格式，自帶有經緯度座標資訊，可自動進行幾何校正，先對反射率資料集進行自動幾何校正。

反射率資料幾何校正：

（1）直接可選擇MODIS的hdf檔案開啟，ENVI自動開啟為分為三個資料集：熱紅外資料發射率（20-36波段），可見光到短波紅外的輻射率資料（1-19、26波段），可見光到短波紅外的反射率資料（1-19、26波段）；

（2）工具/Geometric Correction/Reproject GLT with Bowtie Correction，選擇反射率資料集檔案，點選Spatial Subset，選擇大亞灣的大概位置，點選Spectral Subet，選擇2和19波段。點選OK；

（3）開啟工具/Geometric Correction/Georeference by Sensor/Georeference MODIS，選擇反射率資料集，點選Spatial Subset，選擇大亞灣的大概位置，點選Spectral Subet，選擇2和19波段。點選OK。

圖：選擇資料同時選擇空間子集和光譜子集

（4）在Georeference MODIS Parameter面板設定為UTM， WGS84，49帶，並輸出GCP檔案。

圖：Georeference MODIS Parameter面板

（5）點選OK，在結果輸出面板，設定路徑和檔名輸出。

31、32波段的幾何校正：

（1）雙擊/Geometric Correction/Registration/Warp from GCPs: Image to Map Registration，選擇上一步儲存出來的GCP檔案，設定座標系為：UTM WGS84 49帶，設定輸出解析度均為1000米；

（2）點選OK，選擇輻射定標後的31波段資料，點選Spatial Subset，起止行列號設定為和反射率子集一致，點選OK，

（3）輸出面板上，設定輸出的方法為：Triangulation，重取樣方法為：Bilinear，設定檔名輸出；

同樣的方法對32波段輻射定標的結果進行幾何校正。

第四步：海表溫度反演

       本文使用的是《高懋芳，覃志豪等，MODIS資料反演地表溫度的基本引數估計方法》中分裂窗演算法模型進行海表溫度反演，旨在學習ENVI中的操作流程。

演算法為：

T _s = A ₀ + A ₁*T ₃₁ - A ₂* T ₃₂ ( 1)

其中：T s 是 地 表 溫 度 ( K ) , T 31 和 T 32 分 別 是M OD I S 第 31 和 32 波段的亮度溫度； A 0, A 1 和 A 2 是分裂窗演算法的引數，分別定義如下:

A ₀ = [ D ₃₂( 1 - C₃₁ - D ₃₁) / ( D ₃₂ C₃₁ -D ₃₁ C₃₂) ] a₃₁ - [ D ₃₁( 1 - C₃₂ - D ₃₂) /( D ₃₂ C₃₁ - D ₃₁ C₃₂) ] a ₃₂ ( 2)

A ₁ = 1 + D ₃₁/ ( D ₃₂C₃₁ - D ₃₁ C₃₂) + [ D ₃₂( 1 -C₃₁ - D ₃₁) / ( D ₃₂C₃₁ - D ₃₁ C₃₂) ] b₃₁ ( 3)

A ₂ = D ₃₁/ ( D ₃₂ C₃₁ - D ₃₁C₃₂) + [ D ₃₁( 1 - C₃₂ - D ₃₂) / ( D ₃₂ C₃₁ - D ₃₁C₃₂) ] b₃₂ ( 4)

式中， a ₃₁，b₃₁，a ₃₂ 和 b ₃₂ 是常量, 根據 MODIS的波段特徵確定, 在地表溫度 0 ~ 5 e 範圍內, 這些常量 分 別 可 取 a ₃₁ = -64.60363 , b₃₁ = 0.440817, a ₃₂ = -68.72575, b₃₂ = 0.473453。

上述公式的中間引數分別計算如下:

C_i = Ɛ _iτ_i (ɵ) ( 6)

D _i = [ 1 -τ_i (ɵ)] [ 1 + ( 1 - Ɛ _i) τ_i (ɵ)] ( 7)

式中: i 是指 MODIS的第31和32波段, 分別為 i= 31 或 32; τ_i (ɵ)是視角為 ɵ的大氣透過率; Ɛ _i 是波段 i 的地表比輻射率。

      由以上公式可以看出, 該演算法要求衛星遙感器的 31 和 32 波段資料來計算星上亮度溫度, 同時還要求已知大氣透過率和地表比輻射率, 才能進行地表溫度的反演。

下面是詳細操作步驟：

（1）大氣透過率計算

      大氣透過率τ_i (ɵ)是計算地表溫度的基本引數, 通常是通過大氣水汽含量來估計。經過前人研究，可以用MODIS第 2 和 19 波段來反演大氣水分含量，然後再根據大氣水分含量與大氣透過率之間的關係來估計大氣透過率。對於MODIS影象中的任何一個像元，其可能的大氣水分含量用下式估計

（8）

式中: ω是大氣水分含量( g *cm^-2) , α和β常量，分別α= 0.02 和 β=0.6321；ρ19和ρ2分別是MODIS第19和2波段的地面反射率。

    使用bandmath工具計算大氣水分含量：

    表示式：((0.02-alog(b19/b2))/0.6321)^2

    B19：第19波段反射率

    B2：第2波段反射率

     大氣透過率的計算中，水汽是最主要的考慮因素，毛克彪等將MODTRAN等大氣模型模擬出來的兩者的關係，應用到MODIS資料中，提高了地表溫度反演的精度和實時性，本文采用模擬效果較好的指數關係模擬方程，擬合度達到了0.99以上，公式為：

τ ₃₁= 2.89798-1.88366*exp{-[ω/（-21.22704）]} （9）

τ₃₂= -3.59289+4.60414*exp{-[ω/（-32.70639）]} （10）

式中，ω是水汽含量。

使用bandmath工具計算大氣透過率：

31波段大氣透過率表示式：2.89798-1.88366*exp(b1/21.22704)

B1：大氣水分含量。

32波段大氣透過率表示式：-3.59289+4.60414*exp(b1/32.70639)

B1：大氣水分含量。

（2）地表比輻射率的估算

      地表比輻射率主要取決於地表的物質結構，對MODIS來說，大致分為水面、城鎮和自然表面。對於反演來說，利用混合像元分解的方法，根據植被覆蓋率來計算自然表面和城鎮的比輻射率，水體的可以用常量：

Ɛ _31水體=0.996_，Ɛ _32水體=0.992。

有了這些引數，我們就可以計算C31、C32、D31、D32中間引數，BandMath表示式分別為：

C31=0.996*b31

C32=0.992*b32

D31=(1-b31)*(1+(1-0.996)*b31)

D32=(1-b32)*(1+(1-0.996)*b32)

其中 B31：31波段大氣透過率

B32：32波段大氣透過率

（3）亮度溫度的計算

      將影象DN值定標維熱輻射強度之後, 可用Planck函式求解出星上亮度溫度, 計算公式如下：

T _i = K _{i 2} / l n ( 1 + K _{i 1} /I _i )

式中, K _{i 1}和 K _{i 2} 是常量，對於第 i = 31 波段, 分別為 K _{31 , 1} = 729 .541636 W•m^-2•sr^-1•um^-1 ，

K _{31 , 2} = 1304.413871K ; 對於第 i = 32 波段, 為 K _{32 , 1} = 474 . 6 84780 W•m^-2•sr^-1•um^-1，, K _{31 , 2} = 1196 . 978785 K。

    使用bandmath工具計算31和32的亮溫。

   31波段亮溫：1304.413871/alog(1+729.541636/b31)

    B31: 31波段輻射亮度值

    32波段亮溫：1196.978785 /alog(1+474.684780/b32)

    B32：32波段輻射亮度值

（4）A 0, A 1 和 A 2引數計算

     接下來我們計算A 0, A 1 和 A 2引數，Bandmath表示式分別為：

      A 0=b4*(1-b1-b3)/(b4*b1-b3*b2)*(-64.60363)-b3*(1-b2-b4)/(b4*b1-b3*b2)*(-68.72575)

      A 1=1+b3/( b4*b1-b3*b2)+b4*(1-b1-b3)/( b4*b1-b3*b2)* 0.440817

     A 2=b3/(b4*b1-b3*b2)+b3*(1-b2-b4)/(b4*b1-b3*b2)* 0.473453

其中，b1：C31

B2：C32

B3：D31

B4：D32

（5）溫度計算

   把這些引數帶入公式1中計算溫度值，Bandmath表示式為：

T _s=b0+b1*b31-b2*b32-273

其中：b0：A0引數

B1：A1引數

B2：A2引數

B31：B31亮溫值

B32：B32亮溫值

      如下為得到的最終地表溫度反演結果，單位為攝氏度。通過/Statistics/Compute Statistics工具統計看到最初的結果會有一些數量很少的異常值，幾十個畫素，這些異常值大多是處於影像的邊緣。如下圖為統計的結果，小於-40的像元佔0.6%，大於32度像元佔0.1%。可以將這些像元處理，如用以下bandmath處理：-40>b1<32。

圖：初始結果統計

圖：最終反演結果