遥感实验7教材.docx

遥感实验7教材.docx

《遥感实验7教材.docx》由会员分享，可在线阅读，更多相关《遥感实验7教材.docx（20页珍藏版）》请在冰豆网上搜索。

遥感实验7教材

实验7利用ETM+数据反演地表温度

1实验目的与要求

实习目的：

利用成都市Landsat7/8遥感影像中的热红外波段数据，计算两个时期内地表亮度温度，并分析地表温度的空间特征。

实习要求：

⑴掌握利用遥感影像进行地表温度反演的原理与方法。

⑵几何校正。

由于热红外波段影像丢失空间参考信息，故需要进行几何校正（配准）处理，实习中建议使用多项式法。

⑶地表亮温空间特征分析。

对计算结果进行分析，重点分析地表温度的空间分异特征，并分析其与地表覆被变化之间的关系。

⑷采用ArcGIS、ENVI或Erdas软件进行地表温度。

⑸详细记录操作过程，比较各种算法的优缺点并整理成实验报告上交。

2实验数据

⑴成都市Landsat7/8遥感影像中的热红外波段数据。

⑵矢量数据：

成都市行政边界。

3实验分组与时间

⑴本次实验采取分组方式进行，每组4-5人，设组长一名，组长负责协调并记录每人的操作内容（工作量）。

⑵每组交一份实验报告。

⑶时间为两次实验课，如果可能完成不了，每组利用课后时间完成。

4实验性质

本实验属于软件操作练习。

5参考材料

5.1计算公式

⑴辐射定标运算。

即由像元DN值计算辐射亮度L：

⑵亮度温度计算。

即由辐射亮度L计算亮度温度：

其中，K1的取值为666.09W/（m2srμm），K2的取值为1282.71K。

5.2ENVI下利用ETM+数据反演地表温度

地表温度作为地球环境分析的重要指标，而遥感技术作为现代重要的对地观测手段，使得基于遥感图像的地表温度反演的研究越来越多。

主要的地表温度反演方法有：

大气校正法，单窗算法，单通道法等等。

本文介绍用辐射传输方程法对地表温度进行反演。

技术流程：

例子数据为2002年9月2日的襄樊市LandsatETM+数据。

根据数据的特点以及地表温度反演研究的技术要求，采用的技术路线为：

先对LandsatETM+数据进行预处理：

数据读取、辐射定标、大气校正、襄樊区域裁剪，利用大气校正，即：

辐射传输方程法对其影像热红外波段数据进行操作反演，实现襄樊市地区的地表真实温度的反演研究。

具体的处理流程如下：

具体的实现步骤如下：

第一步：

准备数据

热红外数据使用的是Landsat的第六波段，已经做了传感器定标、几何校正、工程区裁剪，详细流程参考上面的流程图。

文件为TM6-rad-subset-jz-xiangfan.img。

由TM影像（已经过大气校正）生成的NDVI数据，已经利用主菜单->BasicTools->ResizeData（SFatial/SFectral）重采样为60米分辨率，与TMi6数据保持一致，文件名为：

TM-NDVI-60m.img。

第二步：

地表比辐射率计算

物体的比辐射率是物体向外辐射电磁波的能力表征。

它不仅依赖于地表物体的组成，而且与物体的表面状态（表面粗糙度等）及物理性质（介电常数、含水量等）有关，并随着所测定的波长和观测角度等因素有关。

在大尺度上对比辐射率精确测量的难度很大，目前只是基于某些假设获得比辐射率的相对值，本文主要根据可见光和近红外光谱信息来估计比辐射率。

（一）植被覆盖度计算

计算植被覆盖度Fv采用的是混合像元分解法，将整景影像的地类大致分为水体、植被和建筑，具体的计算公式如下：

FV =（NDVI-NDVIS）／（NDVIV -NDVIS） 

（2）

其中，NDVI为归一化差异植被指数，取NDVIV =0.70和NDVIS =0.00，且有，当某个像元的NDVI大于0.70时，FV取值为1；当NDVI小于0.00，FV取值为0。

利用ENVI主菜单->BasicTools->BandMath，在公式输入栏中输入：

（b1gt0.7）*1+（b1lt0.）*0+（b1ge0andb1le0.7）*（（b1-0.0）/（0.7-0.0））

b1：

选择NDVI图像

得到植被盖度图像。

（二）地表比辐射率计算

根据前人的研究，将遥感影像分为水体、城镇和自然表面3种类型。

本专题采取以下方法计算研究区地表比辐射率：

水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式（3）（4）进行计算：

εsurface =0.9625+0.0614FV -0.0461FV2  （3）

εbuilding =0.9589+0.086FV -0.0671FV2   （4）

式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率。

利用ENVI主菜单->BasicTools->BandMath，在公式输入栏中输入：

（b1le0）*0.995+（b1gt0andb1lt0.7）*（0.9589+0.086*b2-0.0671*b2^2）+（b1ge0.7）*（0.9625+0.0614*b2-0.0461*b2^2）

b1：

NDVI值；

b2：

植被覆盖度值。

得到地表比辐射率数据。

第三步：

计算相同温度下黑体的辐射亮度值

辐射传输方程法，又称大气校正法，其基本思路为：

首先利用与卫星过空时间同步的大气数据来估计大气对地表热辐射的影响。

然后把这部分大气影响从卫星高度上传感器所观测到的热辐射总量中减去。

从而得到地表热辐射强度．再把这一热辐射强度转化为相应的地表温度.

卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：

大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。

卫星传感器接收到的热红外辐射亮度值的表达式可写为（辐射传输方程）：

Lλ=[ε·B（TS）+（1-ε）L↓]·τ+L↑  （4）

这里，ε为地表辐射率，TS为地表真实温度，B（TS）为普朗克定律推到得到的黑体在TS的热辐射亮度，τ为大气在热红外波段的透过率。

则温度为T的黑体在热红外波段的辐射亮度B（TS）为：

B（TS）=[Lλ-L↑-τ·（1-ε）L↓]/τ·ε   （5）

在NASA官网（http:

//atmcorr.gsfc.nasa.gov/）中输入成影时间以及中心经纬度，则会提供上式中所需要的参数。

本专题输入的数据是襄樊市地区2002年9月2日北京时间10:

30成像的Landsat7ETM+影像，影像中心的经纬度为：

32.51N，111.81E。

得到下图参数图：

大气在热红外波段的透过率τ为0.6，大气向上辐射亮度L↑为3.39W/（m2·sr·μm），大气向下辐射亮辐射亮度L↓为5.12W/（m2·sr·μm）。

图4.412002年9月2日LandsatETM+数据的大气辅助参数

利用ENVI主菜单->BasicTools->BandMath，在公式输入栏中输入：

（b2-3.39-0.6*（1-b1）*5.12）/（0.6*b1）

b1：

60m分辨率的地表比辐射率值；

b2：

表示热红外波段的辐射定标值。

得到了温度为T的黑体在热红外波段的辐射亮度值。

第四步：

反演地表温度

在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS：

TS =K2/ln（K1/B（TS）+1）

对于ETM+，K1=666.09W/（m2·sr·μm），K2=1282.71K。

利用ENVI主菜单->BasicTools->BandMath，在公式输入栏中输入：

（1282.71）/alog（666.09/b1+1）-273

b1：

温度为T的黑体在热红外波段的辐射亮度值。

得到真实的地表温度值，单位是摄氏度。

第五步：

结果浏览与输出

在DisFlay中显示温度值，是一个灰度的单波段图像。

（1）选择Tools->ColorMaFFing->DensitySlice，单击ClearRange按钮清除默认区间。

（2）选择OFions->AddNewRanges，增加以下四个区间：

   39℃以上，红色

   35℃至39℃，黄色

   30℃至35℃，绿色

   低于30℃，蓝色

（3）单击Apply。

（4）选择File->OutFutRangetoClassImage，可以将反演结果输出。

图4.42地表温度反演结果

地表温度反演结果可以用于很多地方，如城市热岛监测、计算土壤湿度指数（NDVI/T）等。

5.3基于MODIS数据的海表温度反演

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

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

表达式为：

TS=T4+A（T4-T5）+B

其中：

TS为地表真实温度，T4和T5分别为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/.

下面详细介绍处理流程。

第一步：

读取原始数据

打开File->OpenAs->Genericformats->HDF4，选择hdf文件，选择发射率数据集，点击OK，选择数据存储方式BSQ，点击OK。

打开的是热红外原始数据集，第20-36波段，共16个波段，分别是：

20、21、22、23、24、25、27、28、29、30、31、32、33、34、35、36波段。

第二步，辐射亮度定标

打开/RasterManagement/Data-SpecificUtilities/ViewHDFDatasetAttributes，选择原始的hdf数据，点击OK，选择相应的热红外数据集：

EarthView1KMEmissiveBandsScaledIntegers，点击OK，得到一个文件说明面板：

radiance_scales，和radiance_offset这两项参数代表波段的增益和偏移量，是辐射定标的系数。

比如要计算31波段的辐射亮度，读取到scales为0.00084002，offsets为1577.33972168，带入MODIS辐射定标的公式：

Radiance=scales*（DN-offsets），即可以得到该波段的辐射亮度。

打开/BandAlgebra/BandMath工具，输入公式：

0.00084002*（B31-1577.33972168），（B31是第31波段的DN值），点击OK，选择第31波段数据为B31，设置路径和文件名，点击OK。

得到的结果就是31波段的辐射亮度。

同样的方法得到32波段的辐射亮度，公式为：

0.00072970*（B32-1658.22131348）

第三步：

几何校正

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

反射率数据几何校正：

（1）打开反射率数据集：

File->OpenAs->EOS->MODIS，选择hdf数据文件打开；

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

点击OK；

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

点击OK。

图：

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

（4）在GeoreferenceMODISParameter面板设置为UTM，WGS84，49带，并输出GCP文件。

图：

GeoreferenceMODISParameter面板

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

31、32波段的几何校正：

（1）双击/GeometricCorrection/Registration/WarpfromGCPs:

ImagetoMapRegistration，选择上一步保存出来的GCP文件，设置坐标系为：

UTMWGS8449带，设置输出分辨率均为1000米；

（2）点击OK，选择辐射定标后的31波段数据，点击SpatialSubset，起止行列号设置为和反射率子集一致，点击OK，

（3）输出面板上，设置输出的方法为：

Triangulation，重采样方法为：

Bilinear，设置文件名输出；

同样的方法对32波段辐射定标的结果进行几何校正。

第四步：

海表温度反演

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

算法为：

Ts =A0 +A1*T31 -A2*T32 

（1）

其中：

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

A0 =[D32（1-C31-D31）/（D 32 C31-D31C32）]a31 -[D31（1-C32-D 32）/（D32 C31 -D31 C32）]a32 

（2）

A 1 =1+D 31/（D 32C31 -D 31 C32）+[D 32（1-C31 -D 31）/（D 32C31 -D 31 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）

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

τ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）

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 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的像元占0.6%，大于32度像元占0.1%。

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

-40>b1<32。

图：

初始结果统计

图：

最终反演结果