envi反演水质参数_ENVI5.2中基于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)按照本流程反演出来的结果是SST。陆地上的值可以视为无效值,若要得到正确的陆表温度,需要加入海陆分离的步骤,以及城镇和自然表面的比辐射率计算。
下面详细介绍处理流程。
第一步:读取原始数据
打开File->Open As->Generic
formats->HDF4,选择hdf文件,选择发射率数据集,点击OK,选择数据存储方式BSQ,点击OK。
打开的是热红外原始数据集,第20-36波段,共16个波段,分别是:20、21、22、23、24、25、27、28、29、30、31、32、33、34、35、36波段。
第二步,辐射亮度定标
打开/Raster Management/Data-Specific Utilities/View HDF Dataset
Attributes,选择原始的hdf数据,点击OK,选择相应的热红外数据集:Earth View 1KM Emissive
Bands Scaled Integers,点击OK,得到一个文件说明面板:
radiance_scales,和radiance_offset这两项参数代表波段的增益和偏移量,是辐射定标的系数。比如要计算31波段的辐射亮度,读取到scales为0.00084002,offsets为1577.33972168,带入MODIS辐射定标的公式:Radiance=scales*(DN-offsets),即可以得到该波段的辐射亮度。
打开/Band
Algebra/Band
Math工具,输入公式:0.00084002*(B31-1577.33972168),(B31是第31波段的DN值),点击OK,选择第31波段数据为B31,设置路径和文件名,点击OK。得到的结果就是31波段的辐射亮度。
同样的方法得到32波段的辐射亮度,公式为:0.00072970*(B32-1658.22131348)
第三步:几何校正
MODIS1b数据是hdf格式,自带有经纬度坐标信息,可自动进行几何校正,先对反射率数据集进行自动几何校正。
反射率数据几何校正:
(1)打开反射率数据集:File->Open As->EOS->MODIS,选择hdf数据文件打开;
(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 0 + A 1*T
31 - A 2* T 32 ( 1)
其中:T s 是 地 表 温 度 ( K ) , T 31 和 T 32 分 别 是M OD I S 第 31 和 32
波段的亮度温度; A 0, A 1 和 A 2 是分裂窗算法的参数,分别定义如下:
A 0 = [ D 32( 1 - C31 - D
31) / ( D 32 C31 -D 31
C32) ] a31 - [ D 31( 1 -
C32 - D 32) /( D 32 C31
- D 31 C32) ] a 32 ( 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 ~ 5 e 范围内, 这些常量 分 别 可 取
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 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。
图:初始结果统计
图:最终反演结果
envi反演水质参数_ENVI5.2中基于MODIS数据的海表温度反演相关推荐
- envi反演水质参数_科技前沿基于GOCI静止水色卫星数据的长江口及邻近海域Kd(490)遥感反演及其在机载激光测深预评估中的应用...
点击上方"溪流之海洋人生"即可订阅哦 长江是我国的第一大河流,每年输入东海的悬浮泥沙含量可达1.81亿t(据大通站2000-2011年输沙量资料),巨量的悬浮泥沙在口门堆积,造成长 ...
- 基于Modis数据的地表温度反演
基于Modis数据的北京市地表温度反演 评论区有下载原文和相关资料的链接,自己翻找即可. 操作平台 ENVI 5.5 ArcGIS 10.2 数据源 MODIS B1产品(包含1km 热红外波段) 数 ...
- matlab导入s2p,如何将S参数导入matlab中可用的数据文件
如何将S参数导入matlab中可用的数据文件 导出成*.txt,就可以 Q: The export function of CST MWS support data format as Magnitu ...
- oracle自动释放表空间,Oracle中关于清除数据和释放表空间
一.表的重命名 flashback table test2 to before drop rename to test3;--[to test3]将表重命名 drop table test3 purg ...
- 【Envi】基于单窗算法的地表温度反演实验操作记录
文章目录 比前言还前 前言 参考博客 技术流程 1.数据预处理 1.1数据获取 1.2辐射亮度温度 1.2.1热红外波段辐射定标与亮度 1.2.2辐射亮度温度计算 1.3地表比辐射率 1.3.1多光谱 ...
- Google earth engine(GEE):基于MODIS的LST(地表温度数据)计算一定时间序列的城市热岛强度(UHI),并绘制直方图
2023.01.06更新: 完整流程是这样的,首先第一个代码可以得到城市斑块,第二个代码可以计算这个地区的城市热岛效应强度.有朋友私信问全部的代码,同时自己也修改了一下计算城市热岛强度代码的一些比较繁 ...
- 基于MODIS数据的NDVI与LST相关性分析(IDL代码实现)
1 数据预处理 (1)数据提取 我们可以选取2018年5月初华东地区MODIS中的MOD11A2和MOD13A2的16天合成LST和NDVI产品数据,下载地址:MODIS数据下载 网站下载数据需要注册 ...
- 基于MODIS数据的滁州市冬小麦长势遥感监测研究
目录 摘要 1 1 引言 2 1.1 选题背景与意义 2 1.2 研究现状 2 1.3 研究目标.研究内容及技术路线 3 1.3.1 研究目标 3 1.3.2 研究内容 3 1.3.3 技术路线图 4 ...
- python脚本:向表中插入新数据,删除表中最旧的数据
一张表存储历史数据,最多存储HISTORY_TABLE_MAX_ROWS条数据,当表中数据未达到HISTORY_TABLE_MAX_ROWS,直接插入:如果达到的话需要保证插入新数据的时候将最旧的数据 ...
最新文章
- ES shrink ——一般是结合rollover一起使用的,一开始没有看懂官方shrink文档,当看了这个之后就明白了...
- 第10章 嵌入式linux的调试技术
- JS闭包的理解及常见应用场景
- Java - Jackson JSON Java Parser API
- mysql接口测试_用python实现接口测试(四、操作MySQL)
- Lua 脚本获取 EVAL EVALSHA 命令的参数
- cp分解实现_如何用贝叶斯高斯张量分解修复缺失数据?(Jupyter notebook - Python)
- 兰州中考计算机考试,宜昌、兰州发布中考新政新消息:增加口语考试,采取人机对话形式...
- MNIST数据集格式ubyte转png
- windows 10 Tera Term显示乱码
- 康涅狄格大学计算机科学排名,2019美国硕士研究生cs专业前100排名出炉,四校共坐榜首!...
- 反问疑问_反问、疑问还是设问?
- 经典再现,看到就是赚到。尚硅谷雷神 - SpringBoot 2.x 学习笔记 - 核心功能篇
- 为什么用于开关电源的开关管一般用MOS管而不是三极管
- 用C语言编辑一光年相当于多少米,一光年到底有多远?是光速跑了365天的距离,这样说你就少算了...
- 京东到家订单派发的技术实战
- leetcode1646. 获取生成数组中的最大值
- 天龙架设一条龙教程_新手福利,天龙一条龙优先级顺序分享
- STC89C52RC特点及引脚介绍
- 专家称谷歌收购摩托罗拉意在专利
热门文章
- python作业——随机生成不重复的20以内加法算式
- QTableView右键显示菜单
- 净推荐值软件市场现状、行业发展机遇、主要驱动因素及未来发展趋势
- 如何保证原有特效不变让EDIUS替换素材
- C语言while循环模板,完整版)1《while循环》教学设计模板
- P1099 树的直径 DFS + 二分 / 尺取法
- Android-扫描更新媒体库(图库相册)
- 大数据24小时:“中译语通”完成3.34亿元C轮融资,姚劲波和左凌烨成猎豹移动独立董事
- android打不开ios分享微博,iOS 新浪微博方法[WeiboSDK registerApp:K_SINA_APP_KEY]; 崩溃解决办法...
- 计算机硬盘显示打印机名,当连接打印机提示错误代码0X00000709怎么回事?