ENVI下基于劈窗算法从MODIS数据中反演海表温度

劈窗算法最初是为反演海面温度开发的，具体地说是针对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。

图：初始结果统计

图：最终反演结果