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

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

《ENVI下基于劈窗算法从MODIS数据中反演海表温度.docx》由会员分享，可在线阅读，更多相关《ENVI下基于劈窗算法从MODIS数据中反演海表温度.docx（10页珍藏版）》请在冰豆网上搜索。

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

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

（2013-12-0319:

42:

18）

转载▼

标签：

劈窗算法

modis反演温度

教育

分类：

 遥感解决方案

劈窗算法最初是为反演海面温度开发的，具体地说是针对NOAA/AVHRR的4和5通道设计的，后来也被用来反演地表温度，这种算法较成熟，精度也高。

劈窗算法以地表热辐射传导方程为基础，利用10~13μm 大气窗口内，两个相邻热红外通道（一般为10.5～11.5μm、11.5～12.5μm）对大气吸收作用的不同，通过两个通道测量值的各种组合来剔除大气的影响，进行大气和地表比辐射率的修正。

表达式为：

T S=T 4+A（T 4-T 5）+B

其中：

T S为地表真实温度，T 4和 T 5分别为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    技术流程图

注：

（1）按照本流程反演出来的结果是SST。

陆地上的值可以视为无效值，若要得到正确的陆表温度，需要加入海陆分离的步骤，以及城镇和自然表面的比辐射率计算。

（2）MODIS数据下载：

nasa官网：

http:

//modis.gsfc.nasa.gov/.

    下面详细介绍处理流程。

第一步：

数据打开

Open->Openas->EOS->MODIS，选择MOD021KM.A2013185.0245.005.2013185094144.hdf文件，点击OK打开。

分为三个数据集：

热红外数据（20-36波段），可见光到短波红外的辐射率数据（1-19、26波段），可见光到短波红外的反射率数据（1-19、26波段）。

第二步：

辐射定标

由于ENVI默认的设置是对MODIS数据进行自动的辐射定标，所以打开的数据即是经过了辐射定标的数据。

图2    打开modis02级数据

第三步：

几何校正

MODIS1b数据是hdf格式，自带有经纬度坐标信息，可自动进行几何校正，分别对热红外数据集和反射率数据集进行自动几何校正。

反射率数据集合校正：

（1）           打开工具/GeometricCorrection/GeoreferencebySensor/GeoreferenceMODIS，选择反射率数据集，点击SpatialSubset，选择大亚湾的大概位置，点击SpectralSubet，选择2和19波段。

点击OK。

图3    选择数据同时选择空间子集和光谱子集

（2）           在GeoreferenceMODISParameter面板设置为UTM， WGS84，49带。

图4    GeoreferenceMODISParameter面板

（3）           点击OK，在结果输出面板，设置路径和文件名输出。

热红外数据集几何校正：

方法与反射率数据集校正方法一样，选择31和32两个波段校正。

第四步：

海表温度反演

本文使用的是《高懋芳，覃志豪等，MODIS数据反演地表温度的基本参数估计方法》中分裂窗算法模型进行海表温度反演，旨在学习ENVI中的操作流程。

算法为：

T s =A 0 +A 1*T 31 -A 2*T 32                                        （公式1）

其中：

Ts 是 地 表 温 度 （K）,T31 和 T32 分 别 是MODIS 第 31 和 32 波段的亮度温度； A0,A1 和 A2 是分裂窗算法的参数，分别定义如下:

A 0 =[D 32（1-C31 -D 31）/（D 32 C31 -D 31 C32）]a31 -[D 31（1-C32 -D 32）/（D 32 C31 -D31 C32）]a 32                                                                  （公式2）

A 1 =1+D 31/（D32C31 -D 31 C32）+[D 32（1-C31 -D 31）/（D 32C31 -D31 C32）]b31                                                                                         （公式3）

A 2 =D 31/（D 32 C31 -D 31C32）+[D 31（1-C32 -D 32）/（D 32 C31 -D 31C32）]b32                                                                                              （公式4）

式中， a 31，b31，a 32 和 b 32 是常量, 根据 MODIS的波段特征确定, 在地表温度 0~5e 范围内, 这些常量 分 别 可 取 a 31 = -64.60363,b31 =0.440817,a 32 = -68.72575,b32 =0.473453。

上述公式的中间参数分别计算如下:

         Ci = Ɛ iτi （ɵ）                                      （公式6）

         Di =[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以上，公式为：

τ31=2.89798-1.88366*exp{-[ω/（-21.22704）]}          （公式9）

τ32=-3.59289+4.60414*exp{-[ω/（-32.70639）]}         （公式10）

式中，ω是水汽含量。

      使用bandmath工具计算大气透过率：

31波段大气透过率=2.89798-1.88366*exp（b1/21.22704）

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.992）*b32）

其中 B31：

31波段大气透过率

B32：

32波段大气透过率

（3）          亮度温度的计算

将图像DN值定标维热辐射强度之后, 可用Planck函数求解出星上亮度温度, 计算公式如下：

T i =K i2 /ln（1+K i1 /I i ）

式中,K i1和 K i2 是常量，对于第 i=31 波段, 分别为 K 31,1 =729.541636W•m-2•sr-1•um-1 ，

K 31,2 =1304.413871K; 对于第 i=32 波段, 为 K 32,1 =474.684780W•m-2•sr-1•um-1，,K 31,2=1196.978785K。

     使用bandmath工具计算31和32的亮温。

     31波段亮温=1304.413871/alog（1+729.541636/b31）

     B31:

31波段辐射亮度值

     32波段亮温=1196.978785/alog（1+474.684780/b32）

     B32：

32波段辐射亮度值

（4）       A0,A1 和 A2参数计算

     接下来我们计算A0,A1 和 A2参数，Bandmath表达式分别为：

     A0=b4*（1-b1-b3）/（b4*b1-b3*b2）*（-64.60363）-b3*（1-b2-b4）/（b4*b1-b3*b2）*（-68.72575）

     A1=1+b3/（b4*b1-b3*b2）+b4*（1-b1-b3）/（b4*b1-b3*b2）*0.440817

 A2=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/ComputeStatistics工具统计看到最初的结果会有一些数量很少的异常值，几十个像素，这些异常值大多是处于影像的边缘。

如下图为统计的结果，小于-40的像元有57个，大于32度的34个像元。

可以将这些像元处理，如用以下bandmath处理：

-40>b1<32。

图5    初始结果统计

图6    最终反演结果