第一次用 ENVI5.3 把 Landsat 8 的单窗算法跑完出图那一刻我盯着屏幕愣了半分钟研究区东北角那片水库的温度比旁边的城区还高出 6K而水体的比热容决定了它白天升温不可能比水泥地快。后来顺着参数一条一条往回查才发现是大气平均作用温度那一项单位搞错了摄氏度和开尔文混用白送进去将近 2K 的偏差再经过公式里系数放大最后就变成了水库发烧。温度反演这件事难的地方从来不在公式本身。覃志豪单窗算法的表达式统共不超过两行任何一个会按计算器的人都能算。真正拉开差距的是那七八个输入参数——亮温、比辐射率、大气透射率、大气平均作用温度、植被覆盖度——每一个都有自己的取值范围、获取途径和单位约定任何一个环节的疏忽都会在最终结果里被放大而且是那种看起来正常、细看不对的错。这篇东西适合三类人看刚接触地表温度反演、手里有 ENVI5.3 和 Landsat 数据但不知道从哪下手的人跑出结果但心里没底、想搞清楚每个参数来龙去脉的人以及需要批量处理时间序列、要求结果可复现、打算把流程固化下来的人。下面我把整个链条按实际操作顺序拆开包括每一步的 ENVI 表达式、参数取值依据、我自己踩过的坑以及一个用来判断参数精度值多少钱的灵敏度小实验。1. 开工前把参数清单钉在墙上1.1 单窗算法的输入门槛一个热红外波段就够单窗算法这个名字里的单窗指的就是只用一个热红外通道。这一点决定了它的适用面凡是搭载了单热红外波段、并且能得到定标系数的传感器理论上都能用。Landsat 5 的 TM 第 6 波段、Landsat 7 的 ETM 第 6 波段、Landsat 8/9 的 TIRS 第 10 波段是最常见的选择国产的高分五号、资源系列热红外数据也有对应的系数但系数来源要自己查文献不能拿 Landsat 的往上报。我这里主要讲 Landsat 8 第 10 波段因为它数据最容易拿、系数最稳定、社区里讨论也最多。为什么不用第 11 波段从 2014 年前后开始TIRS 第 11 波段被证实存在明显的杂散光污染定量应用里基本被弃用做温度反演就老老实实用第 10 波段。这一点很多刚入门的人不知道看到 TIRS 有两个热红外波段就想两个都用上做劈窗结果越做越乱。数据产品版本也要盯住。Landsat 8 的 Level-1 产品分 Collection 1 和 Collection 2两者的定标系数、尤其是热红外的辐射偏移量不一样。如果你从网上抄了一份参数笔记先看它是哪一年的、对应哪个 Collection。我的习惯是把每一景影像 MTL.txt 里的 RADIANCE_MULT_BAND_10、RADIANCE_ADD_BAND_10、K1_CONSTANT_BAND_10、K2_CONSTANT_BAND_10 四个数抄进一个表按景核对不靠记忆。1.2 K1、K2 和定标系数参考值可以背但要用MTL核对亮温反算要用的 K1、K2 是普朗克函数的简写形式常数不同传感器的值不一样。下面这张表里的数是参考值实际跑的时候还是以 MTL 里的为准传感器热红外波段K1 (W/(m²·sr·µm))K2 (K)Landsat 5 TMBand 6607.761260.56Landsat 7 ETMBand 6666.091282.71Landsat 8 TIRSBand 10774.88531321.0789Landsat 8 TIRSBand 11480.88831201.1442注意一个细节K1 的单位里带 µm也就是说这些系数是基于波段有效波长处、按波段响应函数积分之后的结果不是中心波长的单色值。有些老教程里给的是单色普朗克常数用它算出来的亮温会系统性偏那么零点几开尔文。零点几开尔文听着不多但前面说过亮温误差在公式里是被放大的放大倍数后面会算给你看。Landsat 8 的辐射定标系数MULT 和 ADD我是强烈建议直接从 MTL 里读不要背。原因很简单USGS 对 TIRS 做过重新定标不同版本的偏移量不一样背错了就是整幅图系统性偏差而且这个偏差在图上完全看不出来——图很漂亮色调很合理就是数值不对。1.3 那个最容易被敷衍过去的输入近地表气温单窗算法的参数表里大气平均作用温度 Ta这一项不是直接从影像上算出来的它需要用成像时刻的近地表气温 T0 去估算。这就是整个流程里唯一一个影像之外的输入也是最容易被敷衍的地方。常见的敷衍做法有两种。一种是拿研究区所在城市当天的日最高气温顶上去另一种是随手填个 25。这两种做法带来的误差量级不一样但都在一度以上。正经一点的做法是找成像时刻Landsat 8 过境大约是当地太阳时 10:30 前后最近整点的气象站实测气温多个站做空间插值插到与影像相同的网格上。站点少的时候可以用反距离权重站点够多、地形起伏大的时候用带高程协变量的方法更稳一些。如果研究区附近压根没有气象站或者你需要处理的是一二十年的历史影像那就只能靠再分析资料了——ERA5 的 2m 气温产品是小时分辨率用成像时刻那一层去取就行。再分析资料的空间分辨率比较粗做区域尺度研究够用做城市街区尺度的话要意识到它能提供的只是背景场。提示气象站测的是 1.5 米高度的气温不是地表温度。这一点在后面的精度验证环节会再强调一次因为拿气温去验证反演出来的地表温度是最常见的错误对比方式。2. 亮温反算两行算式三种典型错法2.1 辐射定标在ENVI 5.3里有两条路把 DN 值变成辐亮度ENVI 5.3 提供了工具箱路径Toolbox 里找 Radiometric Correction 分类下的 Radiometric Calibration选好输入文件、传感器类型、输出辐亮度单位W/(m²·sr·µm)勾上要处理的波段直接出结果。这条路适合只有一两景数据的时候界面友好不容易出错。但如果要处理时间序列我更推荐直接读 MTL 里的系数自己用 Band Math 算。理由有三个一是批量的时候工具箱得一个个点效率太低二是自己算的时候系数是你亲眼看到、亲手填进去的出了问题能回溯三是自己算的时候可以在同一个表达式里把无效值处理掉工具箱给你的结果里 DN0 的填充区域会变成一个很小的辐亮度值混进后续计算里。Archive自己算的表达式就是最普通的线性定标Lλ RADIANCE_MULT_BAND_10 * b1 RADIANCE_ADD_BAND_10b1 就是定标前的 DN 波段。填系数的时候把 MTL 里的实际数值代进去别用变量名占位ENVI 的 Band Math 不认 MTL 变量。2.2 浮点类型ENVI Band Math 最阴的一个坑Band Math 的输出类型默认是取所有输入波段里最高的那个数据类型。如果你输入的是一个 16 位整型的原始波段输出很可能会是整型那么这个表达式774.8853 / alog(1321.0789 / float(b1) 1)里面的除法就会变成整数除法结果直接被截断成整数。亮温被截断到整数开尔文等于你主动放弃了所有小数位最后反演出来的地表温度精度也就跟着崩了。我的做法是养成习惯所有表达式里第一个出现的波段变量都用float()包一层输入文件也尽量先用 ENVI 的 Convert Data 转成浮点型。输出类型那一栏手动选 Float32 或者 Float64不要用默认。ENVI 里对数的函数名是alog()自然对数别写成log()——IDL 里log()是常用对数用错了结果差得离谱而且不会报错。另外提一句b1 在 Band Math 里代表你映射的第一个波段这个映射关系是手动指定的。一个表达式里用了 b1、b2、b3 三个变量就要在下方列表里把三个波段都配上少配一个 ENVI 会直接报错配错顺序它可不报错——这是我认为 ENVI 最该改的一个交互设计。所以每次点 OK 之前我会把映射列表念一遍。2.3 亮温量级自查出图之前先看直方图亮温算完别急着往下走。有一个三秒钟的自查打开统计看最小值、最大值、均值。中纬度地区白天 Landsat 过境时刻地表亮温一般落在 285K 到 320K 之间均值大概在 295K 到 305K。如果你看到均值只有 250K那八成是辐亮度单位错了比如把 W/(m²·sr·µm) 换成了 W/(m²·sr·nm) 没换算回来如果看到 2000 这种数那是把开尔文和摄氏度搞混了或者 K1、K2 填反了。还有一个更隐蔽的情况结果看起来正常但研究区边缘有一圈异常高值。这通常是影像边缘的填充像元DN0参与计算造成的辐亮度算出来是一个很小的值套进对数公式之后亮温被抬得很高。解决办法是在定标那一步就加掩膜或者在整个流程开始前用 QA_PIXEL 波段把无效区域、云、云影都标出来做成一个掩膜文件最后一步统一应用。3. 地表比辐射率动 0.01 值多少度3.1 给整幅图赋一个 0.95 是最省事也最不负责的做法比辐射率是单窗算法里对结果影响最直接、同时最容易被简化处理的参数。我见过太多人直接给整幅影像赋一个 0.95 或者 0.97理由听起来也挺合理研究区大部分是植被和建筑取个平均值差不多。问题在于比辐射率的空间差异没那么小。干燥裸土的比辐射率在 0.92 上下全植被像元可以到 0.986水体能到 0.995。这三类地物在同一个城市研究区里往往同时存在跨度接近 0.08。按后面灵敏度实验算出来的比例0.08 的差异足以让不同地类的温度偏差超过 3K——而且这个偏差是系统性的城区偏冷、水体偏热恰好会把你想要分析的城市热岛信号搅乱甚至得出相反结论。所以比辐射率这一步必须按像元算不能赋常数。常用的是 NDVI 阈值法思路朴素但有效植被越多比辐射率越高裸土越多比辐射率越低中间状态按植被覆盖度加权再补一个混合像元的腔体效应修正项。3.2 NDVI、植被覆盖度 Pv 的 ENVI 表达式先算 NDVI。Landsat 8 的近红外是第 5 波段红光是第 4 波段(float(b1) - b2) / (float(b1) b2)这里 b1 是近红外b2 是红光。同样float()不能省尤其分母里如果有整型除法会出问题。然后是植被覆盖度 Pv用分段的方式把 NDVI 映射到 0 到 1 之间(b1 ge 0.157) * (b1 le 0.727) * ((b1 - 0.157) / 0.57)^2 (b1 gt 0.727) * 1b1 这里是 NDVI。0.157 和 0.727 这两个阈值来自覃志豪等人的研究是国内文献里用得最多的一组。为什么用乘法而不是用 AND因为 ENVI 不同版本对逻辑运算符的解析偶尔会出问题用(条件1) * (条件2) * 表达式这种乘法写法条件成立返回 1不成立返回 0乘法天然起到与的作用兼容性最好。这一招我在处理十几年的 Landsat 时间序列时救过很多次场。这段表达式其实同时干了三件事NDVI 小于 0.157 的像元纯裸土或近似裸土Pv 得到 0NDVI 大于 0.727 的像元全植被Pv 得到 1中间的像元得到按平方关系插值的覆盖度。平方关系来自像元的几何分布假设不是随手加的替换成线性插值会让中段像元的 Pv 系统性偏高。3.3 水体、混合像元修正和全植被端的取值细节拿到 Pv 之后陆地部分的比辐射率可以合成一个线性表达式0.923 0.0668 * b1b1 是 Pv。推导过程是裸土端 ε 取 0.923植被端取 0.986混合像元再加一个腔体效应修正 dε 0.0038 × Pv平地情况合起来就是 0.923 × (1-Pv) 0.986 × Pv 0.0038 × Pv 0.923 0.0668 × Pv。这里有个值得注意的细节也是我看过好几篇论文都没讲清楚的地方这个线性式在 Pv 1 的时候得到的是 0.9898而不是 0.986。因为修正项在端点没有减掉。0.004 的差值对应大约 0.2K 的温度偏差。要不要在意我的做法是分两段处理——NDVI 大于 0.727 的像元直接赋 0.986其余像元用线性式多一次 Band Math但心里踏实。如果你对精度要求没那么苛刻用统一线性式也说得过去只是要在方法描述里写明白。水体单独处理NDVI 小于 0 的像元赋 0.995。合成表达式可以写成(b1 lt 0) * 0.995 (b1 ge 0) * (0.923 0.0668 * b2)b1 是 NDVIb2 是 Pv。这个写法逻辑上没错但要提醒一句NDVI 小于 0 的不一定是水体云、云影、深色屋顶、以及部分不透水面都可能出现负值。在城区研究里用 NDVI 负值判水体一定要配合掩膜检查否则你会在城市中心给一堆屋顶赋上 0.995 的水体比辐射率这些像元最后会莫名其妙偏热。3.4 一个灵敏度实验参数精度值多少钱讲到这里很多人心里会有一个疑问这么多参数到底哪个的精度最重要我在自己的数据上做过一个简单的敏感性测试参数基准设为比辐射率 0.98、大气透射率 0.85、亮温 27.0℃、大气平均作用温度 22.0℃跑出来地表温度约 26.96℃。然后每次只动一个参数参数变化变化幅度反演结果偏差比辐射率 0.98 → 0.990.0127.43℃0.47K大气透射率 0.85 → 0.80-0.0527.40℃0.44K大气平均作用温度 22 → 23℃1K26.78℃-0.18K亮温 27 → 28℃1K28.15℃1.19K这张表信息量很大。第一亮温的误差在公式里被放大了 1.19 倍所以前面定标那一步的精度比什么都重要第二比辐射率每偏差 0.01 就带来约 0.5K 的误差按这个比例给整幅图赋常数 0.95 而真实值在 0.92 到 0.99 之间浮动误差能到两度以上第三反而是大气平均作用温度的敏感性最低1K 的气温误差只带来 0.18K 的结果偏差——这不代表它不重要而是说这个参数允许你在数据源上稍微宽容一点。要注意这张表是特定参数组合下的结果换一组 τ 和 ε具体数值会变但量级关系是稳的。你也可以用同样的方法给自己的研究区做一遍把基准值换成你研究区的典型值这样你就知道该在哪一步多花时间。4. 大气透射率和大气平均作用温度参数表怎么看4.1 水汽含量 w 的三条路各有各的适用场景大气透射率 τ 不能直接测通常先拿到大气水汽含量 w再用经验关系换算。水汽含量的获取路径有三条。第一条是用在线的大气校正参数计算器。输入影像的成像日期、时间、中心经纬度它会返回水汽含量和各个波段的透射率估计值。这条路最省事适合零星几景影像。缺点也明显它给的是一个点的值整幅影像用同一个 τ空间上是均一的。对于一景 185 公里幅宽的影像如果研究区里既有平原又有山地这个均一假设会有点勉强。第二条是用中分辨率成像光谱仪的水汽产品空间分辨率 1 公里能反映出区域尺度的水汽差异。用的时候要按成像日期选最接近的一天然后重采样到 Landsat 的网格上。这条路的精度够用而且能保留空间变化是我做区域研究时最常用的。第三条是用再分析资料的总柱水汽含量时间分辨率高空间尺度粗适合长时间序列批量处理。缺点是分辨率低做小区域研究时要接受它提供的是背景值。4.2 τ 的经验拟合式系数不能跨传感器借拿到 w 之后很多文献给了 τ 和 w 的线性拟合式。比如针对 TM 和 ETM 热红外波段就有在特定水汽区间上的线性关系形式是 τ k1 - k2 × w两个系数在不同的水汽区间里取不同的值。这里我要非常明确地提醒一件事**这类拟合式的系数是跟传感器、波段、以及拟合时使用的大气廓线数据集绑定在一起的不能跨传感器借。**TM 第 6 波段的系数套到 Landsat 8 第 10 波段上有多少误差取决于两者的波段响应函数差异没人能给你一个通用的换算系数。我见过有人直接拿 TM 的系数跑 Landsat 8结果整体偏冷两三度怎么调都调不回来。所以对 Landsat 8 第 10 波段我的建议是优先用在线计算器给出的 τ 值或者用 MODTRAN 之类的辐射传输模型按实际大气廓线算一遍。经验拟合式可以拿来做交叉验证——如果两种途径给出的 τ 差了 0.02 以上就要回头查是不是水汽含量取错了或者拟合式的适用水汽区间不对。顺便说一个容易忽略的点τ 的取值一般在 0.5 到 0.95 之间。如果你算出来一个大于 1 的值那一定是系数用错了或者水汽单位搞错了g/cm² 和 kg/m² 差 10 倍这个坑很深。看到超出范围的值先别往下算。4.3 Ta 的季节和纬度带选择大气平均作用温度 Ta 的估算用的是近地表气温 T0 的线性关系适用区域估算式热带Ta 17.9769 0.91715 × T0中纬度夏季Ta 16.0110 0.92621 × T0中纬度冬季Ta 19.2704 0.91118 × T0三个式子的形式一模一样差别在截距和斜率。选哪一个看研究区的纬度和成像季节。国内大部分地区属于中纬度夏季和冬季各有一套春秋季怎么选是个让人纠结的事。我自己的处理原则是看成像月份和当地物候四月到九月按夏季式十月到次年三月按冬季式三月和十月这种过渡月份两个式子各算一遍取平均。T0 的单位是摄氏度这一点必须记住因为 Ta 也是摄氏度最后代入主公式时也要保持摄氏度。前面开头讲的那个水库发烧的故事就是某一步把 T0 从摄氏转成了开尔文再套式子Ta 被抬高了 273 度——不对是抬高了 273 倍的 0.926反正是彻底跑飞了。提示Ta 这个参数的敏感性虽然不高但它的错误来源往往是单位而不是数值精度。用同一套换算函数处理所有输入不要在中途手工转换。5. 单窗算法公式落地C、D 先算主式后算5.1 公式结构拆解覃志豪单窗算法的主式Ts [a(1-C-D) (b(1-C-D) C D) × T6 - D × Ta] / C其中 C ε × τD (1-τ) × [1 (1-ε)τ]。a 和 b 是拟合系数跟温度区间有关0 到 30℃ 区间取 a -60.3263、b 0.434360 到 70℃ 区间取 a -67.355351、b 0.458606。选哪一组看研究区的温度范围一般区域研究用第一组地表温度可能很高的干旱区用第二组。这个选择不是随便的两个系数是从普朗克函数在不同温度区间做线性近似拟合出来的用错了区间在边缘温度上会偏。在 ENVI 里我的做法是先算 C 和 D 两个中间量各存一个文件再把主式写成一步。为什么不分步算更多中间量因为 1-C-D 这个量在分子里出现了两次算一次存下来能省一次运算而且万一分母的 C 出了负数或者零你从中间文件能立刻看出问题——比如 ε 被误赋成了 0C 就是 0主式会直接除零ENVI 给出的结果是一堆极大值图上是一片刺眼的白。C 的表达式b1 * b2b1 是比辐射率b2 是大气透射率。D 的表达式(1 - b2) * (1 (1 - b1) * b2)5.2 单位一致性一个能白白吃掉 2K 的坑主式里的 T6 和 Ta必须和 a、b 的单位区间一致也就是都用摄氏度。但亮温算出来天然是开尔文所以要先减 273.15。我用基准参数验证过这件事T6 300.15K、Ta 295.15K如果直接代进公式而不是先转成 27℃ 和 22℃得到的结果是 302.14K也就是 28.99℃而正确的答案是 26.96℃。差了整整 2.03K。这个误差不会让图看起来有任何异常只会让所有数值系统性偏高两度非常危险。反过来如果你用的是某个只给开尔文版本的公式有些文献会把公式和系数一起改成开尔文形式那就全程用开尔文别再减 273.15。关键是公式和系数、还有输入单位三者必须是一套的混搭就是找死。5.3 越界值、无效值和掩膜处理主式跑完先做两件事。第一件是看统计量地表温度应该在什么范围中纬度夏季白天植被覆盖区大约 20 到 35℃裸露地表和建筑屋顶到 45℃ 甚至更高都有可能水体在 20 到 28℃。如果你看到负值、或者看到 100 以上的数值别犹豫回头查参数。第二件是处理无效值。主式计算过程里会有一些像元的输入本身就是填充值或者被云污染的极端值这些像元算出来的结果会落在物理上不可能的范围。处理方式是先建掩膜再做统计而不是直接把异常像元删掉——保留原始结果用掩膜控制哪些像元参与统计和出图这样你随时能回溯。掩膜的来源可以是 QA_PIXEL 波段解析出来的云、云影、水体标志位也可以自己用亮温阈值加 NDVI 阈值组合一个。我的习惯是两条路都做取交集宁可保守一点也不要把云边缘的半透明像元放进结果里。5.4 结果转摄氏度、配色和归档单窗算法的输出本身就是摄氏度前提是你按 5.2 说的做了如果你用的是开尔文版本的公式最后减 273.15 得到摄氏度。ENVI 里转一次很方便b1 - 273.15出图配色不要用彩虹色带。温度是连续变量用彩虹色带会在数值梯度平缓的地方造出假的边界看图的读者容易误判。用单色系的渐变色带或者分位数断点都更稳妥。同时记得在图上标出像元值的取值范围不要只给一个色带。归档这件事我特别想强调。ENVI 的 Band Math 表达式在项目切换、软件升级之后经常找不回来我踩过这个坑。现在我的做法是所有表达式写在一个 txt 文件里按步骤编号所有中间结果存成 ENVI 格式文件名里带上关键参数比如LST_swa_tau0.85_Ta22_eps_ndvi.dat同一个项目建一个 README写清楚数据版本、MTL 系数、气象数据来源。半年后要复现十分钟就能跑通。6. 结果对不对三条校验路径和一张异常速查表6.1 和官方温度产品对拍最快的一道体检Landsat 的 Level-2 产品里直接带了地表温度波段 ST_B10量化方式是 DN × 0.00341802 149.0结果单位是开尔文另外还有一个 ST_QA 波段给出每个像元的精度估计。这个产品是用辐射传输加再分析资料大气廓线算出来的跟单窗算法的原理不完全一样但都是同一个传感器、同一时刻的数据拿来对拍非常方便。我的做法是把官方 ST_B10 转到摄氏度和我的单窗结果做逐像元差值看差值的均值和标准差。正常情况下两者应该差不多差值均值在 1K 以内标准差 1.5K 左右散点分布在植被、裸土、城区比辐射率估算的差异会体现出来。如果差值均值超过 3K那基本可以确定是我的某个参数出了问题而不是两种方法本来就有差异。要注意 ST_QA 这一步不能省。官方产品里质量差的像元本身误差就大拿这些像元去对拍等于用坏尺子校准好尺子。6.2 和 MODIS LST、地面观测对比的正确姿势和 MODIS 地表温度产品对比有两个必须处理的问题。第一是时间匹配Terra 星过境时间和 Landsat 接近都是上午十点半左右用 Terra 的 MOD11A1 日产品比用 Aqua 的 MYD11A1 更合理。第二是空间尺度MODIS 是 1 公里Landsat 是 30 米直接逐像元比毫无意义。正确做法是先把 30 米的结果按 1 公里网格求平均升尺度再和 MODIS 像元比。同时用 MODIS 的 QC 层筛掉质量等级差的像元。和地面观测对比是最容易出错的一环。常规气象站测的是 1.5 米高度的气温和地表温度是两个物理量白天晴朗天气下两者差 10K 以上很正常。拿气温验证地表温度得到的误差里绝大部分是物理量本身的不同不是反演误差。真正可用的地面数据是红外测温仪测的表面辐射温度、或者长波辐射表反演的地表温度而且要注意观测的视场角、表面类型和像元尺度是否匹配。如果手上只有常规气象站数据那就老老实实只做相对分析不要报绝对精度。6.3 异常现象速查表下面这张表是我自己攒的按现象反查原因比从头查参数快得多现象最可能的原因排查方向整体偏高 2K 左右单位混用开尔文和摄氏同时进入主式检查 T6、Ta 是否都已减 273.15整体偏高 5K 以上比辐射率被赋成常数且偏低检查 ε 图层的取值范围整体偏低定标系数用了旧版本 Collection 的参数核对 MTL 里的 MULT/ADD水体异常偏热水体像元 ε 赋成了陆地值检查 NDVI 负值区的处理逻辑城区异常偏冷混合像元 dε 修正项漏掉检查线性式里的 0.0668 系数边缘一圈高值填充像元 DN0 参与计算加无效值掩膜结果图上出现条带输入波段本身有条带或掩膜未对齐检查原始数据质量和各图层网格一致性研究区一半是负值NDVI 计算时波段映射错位核对红光和近红外的波段编号倒数第二条值得展开说一句。我在做多景时间序列的时候遇到过一次某一年的结果上有规则条带查了半天最后发现是当年那景影像本身在传感器异常期采集的数据质量就有问题。所以在做趋势分析之前一定要按景看数据质量别把传感器问题当成地表变化。6.4 什么情况下不该硬用单窗算法单窗算法有它明确的假设前提地表近似朗伯体、大气在垂直方向均匀、天气晴朗无云。这几个假设在某些场景下会明显失效。地形起伏剧烈的山区坡度和坡向会改变局部入射角和大气路径长度单窗算法里没有对应的修正项结果会有系统性偏差通常表现为阳坡偏冷、阴坡偏暖。如果研究区是山区要么换用带地形修正的方法要么至少用 DEM 做一个分区分析别把整片山当成平地处理。水汽含量很高的地区比如夏季的南方τ 会比较低而 τ 越低公式对 τ 误差的敏感性越高前面那张灵敏度表里 0.05 的 τ 误差对应 0.44K如果 τ 本身只有 0.6这个放大效应会更明显。这种情况下要么把水汽估计做得更细比如用逐像元的水汽产品而不是单一值要么考虑换用对水汽不那么敏感的方法。还有一种情况是数据本身的限制。如果你手上只有单热红外波段那单窗算法几乎是唯一选择前提是把参数做实。如果你手上有两个热红外波段并且质量可靠劈窗算法在理论上更有优势因为它用两个通道的差异直接消掉了一部分大气影响减少了对辅助大气参数的依赖。但要注意前提是质量可靠——Landsat 8 的第 11 波段就不满足这个前提硬做劈窗反而更差。我自己在实际操作中的体会是单窗算法的门槛不在数学在数据。把 MTL 系数核对清楚、把比辐射率按像元算出来、把水汽含量找对来源、把单位统一好这四件事做到位结果基本就能用。剩下那部分误差更多来自数据本身的时空代表性不是靠调参数能解决的。每次跑完做一遍敏感性测试知道自己的结果对哪个参数最敏感比反复调参数追求看起来合理的图要靠谱得多。
ENVI5.3 Landsat 8 单窗算法地表温度反演参数与单位避坑
第一次用 ENVI5.3 把 Landsat 8 的单窗算法跑完出图那一刻我盯着屏幕愣了半分钟研究区东北角那片水库的温度比旁边的城区还高出 6K而水体的比热容决定了它白天升温不可能比水泥地快。后来顺着参数一条一条往回查才发现是大气平均作用温度那一项单位搞错了摄氏度和开尔文混用白送进去将近 2K 的偏差再经过公式里系数放大最后就变成了水库发烧。温度反演这件事难的地方从来不在公式本身。覃志豪单窗算法的表达式统共不超过两行任何一个会按计算器的人都能算。真正拉开差距的是那七八个输入参数——亮温、比辐射率、大气透射率、大气平均作用温度、植被覆盖度——每一个都有自己的取值范围、获取途径和单位约定任何一个环节的疏忽都会在最终结果里被放大而且是那种看起来正常、细看不对的错。这篇东西适合三类人看刚接触地表温度反演、手里有 ENVI5.3 和 Landsat 数据但不知道从哪下手的人跑出结果但心里没底、想搞清楚每个参数来龙去脉的人以及需要批量处理时间序列、要求结果可复现、打算把流程固化下来的人。下面我把整个链条按实际操作顺序拆开包括每一步的 ENVI 表达式、参数取值依据、我自己踩过的坑以及一个用来判断参数精度值多少钱的灵敏度小实验。1. 开工前把参数清单钉在墙上1.1 单窗算法的输入门槛一个热红外波段就够单窗算法这个名字里的单窗指的就是只用一个热红外通道。这一点决定了它的适用面凡是搭载了单热红外波段、并且能得到定标系数的传感器理论上都能用。Landsat 5 的 TM 第 6 波段、Landsat 7 的 ETM 第 6 波段、Landsat 8/9 的 TIRS 第 10 波段是最常见的选择国产的高分五号、资源系列热红外数据也有对应的系数但系数来源要自己查文献不能拿 Landsat 的往上报。我这里主要讲 Landsat 8 第 10 波段因为它数据最容易拿、系数最稳定、社区里讨论也最多。为什么不用第 11 波段从 2014 年前后开始TIRS 第 11 波段被证实存在明显的杂散光污染定量应用里基本被弃用做温度反演就老老实实用第 10 波段。这一点很多刚入门的人不知道看到 TIRS 有两个热红外波段就想两个都用上做劈窗结果越做越乱。数据产品版本也要盯住。Landsat 8 的 Level-1 产品分 Collection 1 和 Collection 2两者的定标系数、尤其是热红外的辐射偏移量不一样。如果你从网上抄了一份参数笔记先看它是哪一年的、对应哪个 Collection。我的习惯是把每一景影像 MTL.txt 里的 RADIANCE_MULT_BAND_10、RADIANCE_ADD_BAND_10、K1_CONSTANT_BAND_10、K2_CONSTANT_BAND_10 四个数抄进一个表按景核对不靠记忆。1.2 K1、K2 和定标系数参考值可以背但要用MTL核对亮温反算要用的 K1、K2 是普朗克函数的简写形式常数不同传感器的值不一样。下面这张表里的数是参考值实际跑的时候还是以 MTL 里的为准传感器热红外波段K1 (W/(m²·sr·µm))K2 (K)Landsat 5 TMBand 6607.761260.56Landsat 7 ETMBand 6666.091282.71Landsat 8 TIRSBand 10774.88531321.0789Landsat 8 TIRSBand 11480.88831201.1442注意一个细节K1 的单位里带 µm也就是说这些系数是基于波段有效波长处、按波段响应函数积分之后的结果不是中心波长的单色值。有些老教程里给的是单色普朗克常数用它算出来的亮温会系统性偏那么零点几开尔文。零点几开尔文听着不多但前面说过亮温误差在公式里是被放大的放大倍数后面会算给你看。Landsat 8 的辐射定标系数MULT 和 ADD我是强烈建议直接从 MTL 里读不要背。原因很简单USGS 对 TIRS 做过重新定标不同版本的偏移量不一样背错了就是整幅图系统性偏差而且这个偏差在图上完全看不出来——图很漂亮色调很合理就是数值不对。1.3 那个最容易被敷衍过去的输入近地表气温单窗算法的参数表里大气平均作用温度 Ta这一项不是直接从影像上算出来的它需要用成像时刻的近地表气温 T0 去估算。这就是整个流程里唯一一个影像之外的输入也是最容易被敷衍的地方。常见的敷衍做法有两种。一种是拿研究区所在城市当天的日最高气温顶上去另一种是随手填个 25。这两种做法带来的误差量级不一样但都在一度以上。正经一点的做法是找成像时刻Landsat 8 过境大约是当地太阳时 10:30 前后最近整点的气象站实测气温多个站做空间插值插到与影像相同的网格上。站点少的时候可以用反距离权重站点够多、地形起伏大的时候用带高程协变量的方法更稳一些。如果研究区附近压根没有气象站或者你需要处理的是一二十年的历史影像那就只能靠再分析资料了——ERA5 的 2m 气温产品是小时分辨率用成像时刻那一层去取就行。再分析资料的空间分辨率比较粗做区域尺度研究够用做城市街区尺度的话要意识到它能提供的只是背景场。提示气象站测的是 1.5 米高度的气温不是地表温度。这一点在后面的精度验证环节会再强调一次因为拿气温去验证反演出来的地表温度是最常见的错误对比方式。2. 亮温反算两行算式三种典型错法2.1 辐射定标在ENVI 5.3里有两条路把 DN 值变成辐亮度ENVI 5.3 提供了工具箱路径Toolbox 里找 Radiometric Correction 分类下的 Radiometric Calibration选好输入文件、传感器类型、输出辐亮度单位W/(m²·sr·µm)勾上要处理的波段直接出结果。这条路适合只有一两景数据的时候界面友好不容易出错。但如果要处理时间序列我更推荐直接读 MTL 里的系数自己用 Band Math 算。理由有三个一是批量的时候工具箱得一个个点效率太低二是自己算的时候系数是你亲眼看到、亲手填进去的出了问题能回溯三是自己算的时候可以在同一个表达式里把无效值处理掉工具箱给你的结果里 DN0 的填充区域会变成一个很小的辐亮度值混进后续计算里。Archive自己算的表达式就是最普通的线性定标Lλ RADIANCE_MULT_BAND_10 * b1 RADIANCE_ADD_BAND_10b1 就是定标前的 DN 波段。填系数的时候把 MTL 里的实际数值代进去别用变量名占位ENVI 的 Band Math 不认 MTL 变量。2.2 浮点类型ENVI Band Math 最阴的一个坑Band Math 的输出类型默认是取所有输入波段里最高的那个数据类型。如果你输入的是一个 16 位整型的原始波段输出很可能会是整型那么这个表达式774.8853 / alog(1321.0789 / float(b1) 1)里面的除法就会变成整数除法结果直接被截断成整数。亮温被截断到整数开尔文等于你主动放弃了所有小数位最后反演出来的地表温度精度也就跟着崩了。我的做法是养成习惯所有表达式里第一个出现的波段变量都用float()包一层输入文件也尽量先用 ENVI 的 Convert Data 转成浮点型。输出类型那一栏手动选 Float32 或者 Float64不要用默认。ENVI 里对数的函数名是alog()自然对数别写成log()——IDL 里log()是常用对数用错了结果差得离谱而且不会报错。另外提一句b1 在 Band Math 里代表你映射的第一个波段这个映射关系是手动指定的。一个表达式里用了 b1、b2、b3 三个变量就要在下方列表里把三个波段都配上少配一个 ENVI 会直接报错配错顺序它可不报错——这是我认为 ENVI 最该改的一个交互设计。所以每次点 OK 之前我会把映射列表念一遍。2.3 亮温量级自查出图之前先看直方图亮温算完别急着往下走。有一个三秒钟的自查打开统计看最小值、最大值、均值。中纬度地区白天 Landsat 过境时刻地表亮温一般落在 285K 到 320K 之间均值大概在 295K 到 305K。如果你看到均值只有 250K那八成是辐亮度单位错了比如把 W/(m²·sr·µm) 换成了 W/(m²·sr·nm) 没换算回来如果看到 2000 这种数那是把开尔文和摄氏度搞混了或者 K1、K2 填反了。还有一个更隐蔽的情况结果看起来正常但研究区边缘有一圈异常高值。这通常是影像边缘的填充像元DN0参与计算造成的辐亮度算出来是一个很小的值套进对数公式之后亮温被抬得很高。解决办法是在定标那一步就加掩膜或者在整个流程开始前用 QA_PIXEL 波段把无效区域、云、云影都标出来做成一个掩膜文件最后一步统一应用。3. 地表比辐射率动 0.01 值多少度3.1 给整幅图赋一个 0.95 是最省事也最不负责的做法比辐射率是单窗算法里对结果影响最直接、同时最容易被简化处理的参数。我见过太多人直接给整幅影像赋一个 0.95 或者 0.97理由听起来也挺合理研究区大部分是植被和建筑取个平均值差不多。问题在于比辐射率的空间差异没那么小。干燥裸土的比辐射率在 0.92 上下全植被像元可以到 0.986水体能到 0.995。这三类地物在同一个城市研究区里往往同时存在跨度接近 0.08。按后面灵敏度实验算出来的比例0.08 的差异足以让不同地类的温度偏差超过 3K——而且这个偏差是系统性的城区偏冷、水体偏热恰好会把你想要分析的城市热岛信号搅乱甚至得出相反结论。所以比辐射率这一步必须按像元算不能赋常数。常用的是 NDVI 阈值法思路朴素但有效植被越多比辐射率越高裸土越多比辐射率越低中间状态按植被覆盖度加权再补一个混合像元的腔体效应修正项。3.2 NDVI、植被覆盖度 Pv 的 ENVI 表达式先算 NDVI。Landsat 8 的近红外是第 5 波段红光是第 4 波段(float(b1) - b2) / (float(b1) b2)这里 b1 是近红外b2 是红光。同样float()不能省尤其分母里如果有整型除法会出问题。然后是植被覆盖度 Pv用分段的方式把 NDVI 映射到 0 到 1 之间(b1 ge 0.157) * (b1 le 0.727) * ((b1 - 0.157) / 0.57)^2 (b1 gt 0.727) * 1b1 这里是 NDVI。0.157 和 0.727 这两个阈值来自覃志豪等人的研究是国内文献里用得最多的一组。为什么用乘法而不是用 AND因为 ENVI 不同版本对逻辑运算符的解析偶尔会出问题用(条件1) * (条件2) * 表达式这种乘法写法条件成立返回 1不成立返回 0乘法天然起到与的作用兼容性最好。这一招我在处理十几年的 Landsat 时间序列时救过很多次场。这段表达式其实同时干了三件事NDVI 小于 0.157 的像元纯裸土或近似裸土Pv 得到 0NDVI 大于 0.727 的像元全植被Pv 得到 1中间的像元得到按平方关系插值的覆盖度。平方关系来自像元的几何分布假设不是随手加的替换成线性插值会让中段像元的 Pv 系统性偏高。3.3 水体、混合像元修正和全植被端的取值细节拿到 Pv 之后陆地部分的比辐射率可以合成一个线性表达式0.923 0.0668 * b1b1 是 Pv。推导过程是裸土端 ε 取 0.923植被端取 0.986混合像元再加一个腔体效应修正 dε 0.0038 × Pv平地情况合起来就是 0.923 × (1-Pv) 0.986 × Pv 0.0038 × Pv 0.923 0.0668 × Pv。这里有个值得注意的细节也是我看过好几篇论文都没讲清楚的地方这个线性式在 Pv 1 的时候得到的是 0.9898而不是 0.986。因为修正项在端点没有减掉。0.004 的差值对应大约 0.2K 的温度偏差。要不要在意我的做法是分两段处理——NDVI 大于 0.727 的像元直接赋 0.986其余像元用线性式多一次 Band Math但心里踏实。如果你对精度要求没那么苛刻用统一线性式也说得过去只是要在方法描述里写明白。水体单独处理NDVI 小于 0 的像元赋 0.995。合成表达式可以写成(b1 lt 0) * 0.995 (b1 ge 0) * (0.923 0.0668 * b2)b1 是 NDVIb2 是 Pv。这个写法逻辑上没错但要提醒一句NDVI 小于 0 的不一定是水体云、云影、深色屋顶、以及部分不透水面都可能出现负值。在城区研究里用 NDVI 负值判水体一定要配合掩膜检查否则你会在城市中心给一堆屋顶赋上 0.995 的水体比辐射率这些像元最后会莫名其妙偏热。3.4 一个灵敏度实验参数精度值多少钱讲到这里很多人心里会有一个疑问这么多参数到底哪个的精度最重要我在自己的数据上做过一个简单的敏感性测试参数基准设为比辐射率 0.98、大气透射率 0.85、亮温 27.0℃、大气平均作用温度 22.0℃跑出来地表温度约 26.96℃。然后每次只动一个参数参数变化变化幅度反演结果偏差比辐射率 0.98 → 0.990.0127.43℃0.47K大气透射率 0.85 → 0.80-0.0527.40℃0.44K大气平均作用温度 22 → 23℃1K26.78℃-0.18K亮温 27 → 28℃1K28.15℃1.19K这张表信息量很大。第一亮温的误差在公式里被放大了 1.19 倍所以前面定标那一步的精度比什么都重要第二比辐射率每偏差 0.01 就带来约 0.5K 的误差按这个比例给整幅图赋常数 0.95 而真实值在 0.92 到 0.99 之间浮动误差能到两度以上第三反而是大气平均作用温度的敏感性最低1K 的气温误差只带来 0.18K 的结果偏差——这不代表它不重要而是说这个参数允许你在数据源上稍微宽容一点。要注意这张表是特定参数组合下的结果换一组 τ 和 ε具体数值会变但量级关系是稳的。你也可以用同样的方法给自己的研究区做一遍把基准值换成你研究区的典型值这样你就知道该在哪一步多花时间。4. 大气透射率和大气平均作用温度参数表怎么看4.1 水汽含量 w 的三条路各有各的适用场景大气透射率 τ 不能直接测通常先拿到大气水汽含量 w再用经验关系换算。水汽含量的获取路径有三条。第一条是用在线的大气校正参数计算器。输入影像的成像日期、时间、中心经纬度它会返回水汽含量和各个波段的透射率估计值。这条路最省事适合零星几景影像。缺点也明显它给的是一个点的值整幅影像用同一个 τ空间上是均一的。对于一景 185 公里幅宽的影像如果研究区里既有平原又有山地这个均一假设会有点勉强。第二条是用中分辨率成像光谱仪的水汽产品空间分辨率 1 公里能反映出区域尺度的水汽差异。用的时候要按成像日期选最接近的一天然后重采样到 Landsat 的网格上。这条路的精度够用而且能保留空间变化是我做区域研究时最常用的。第三条是用再分析资料的总柱水汽含量时间分辨率高空间尺度粗适合长时间序列批量处理。缺点是分辨率低做小区域研究时要接受它提供的是背景值。4.2 τ 的经验拟合式系数不能跨传感器借拿到 w 之后很多文献给了 τ 和 w 的线性拟合式。比如针对 TM 和 ETM 热红外波段就有在特定水汽区间上的线性关系形式是 τ k1 - k2 × w两个系数在不同的水汽区间里取不同的值。这里我要非常明确地提醒一件事**这类拟合式的系数是跟传感器、波段、以及拟合时使用的大气廓线数据集绑定在一起的不能跨传感器借。**TM 第 6 波段的系数套到 Landsat 8 第 10 波段上有多少误差取决于两者的波段响应函数差异没人能给你一个通用的换算系数。我见过有人直接拿 TM 的系数跑 Landsat 8结果整体偏冷两三度怎么调都调不回来。所以对 Landsat 8 第 10 波段我的建议是优先用在线计算器给出的 τ 值或者用 MODTRAN 之类的辐射传输模型按实际大气廓线算一遍。经验拟合式可以拿来做交叉验证——如果两种途径给出的 τ 差了 0.02 以上就要回头查是不是水汽含量取错了或者拟合式的适用水汽区间不对。顺便说一个容易忽略的点τ 的取值一般在 0.5 到 0.95 之间。如果你算出来一个大于 1 的值那一定是系数用错了或者水汽单位搞错了g/cm² 和 kg/m² 差 10 倍这个坑很深。看到超出范围的值先别往下算。4.3 Ta 的季节和纬度带选择大气平均作用温度 Ta 的估算用的是近地表气温 T0 的线性关系适用区域估算式热带Ta 17.9769 0.91715 × T0中纬度夏季Ta 16.0110 0.92621 × T0中纬度冬季Ta 19.2704 0.91118 × T0三个式子的形式一模一样差别在截距和斜率。选哪一个看研究区的纬度和成像季节。国内大部分地区属于中纬度夏季和冬季各有一套春秋季怎么选是个让人纠结的事。我自己的处理原则是看成像月份和当地物候四月到九月按夏季式十月到次年三月按冬季式三月和十月这种过渡月份两个式子各算一遍取平均。T0 的单位是摄氏度这一点必须记住因为 Ta 也是摄氏度最后代入主公式时也要保持摄氏度。前面开头讲的那个水库发烧的故事就是某一步把 T0 从摄氏转成了开尔文再套式子Ta 被抬高了 273 度——不对是抬高了 273 倍的 0.926反正是彻底跑飞了。提示Ta 这个参数的敏感性虽然不高但它的错误来源往往是单位而不是数值精度。用同一套换算函数处理所有输入不要在中途手工转换。5. 单窗算法公式落地C、D 先算主式后算5.1 公式结构拆解覃志豪单窗算法的主式Ts [a(1-C-D) (b(1-C-D) C D) × T6 - D × Ta] / C其中 C ε × τD (1-τ) × [1 (1-ε)τ]。a 和 b 是拟合系数跟温度区间有关0 到 30℃ 区间取 a -60.3263、b 0.434360 到 70℃ 区间取 a -67.355351、b 0.458606。选哪一组看研究区的温度范围一般区域研究用第一组地表温度可能很高的干旱区用第二组。这个选择不是随便的两个系数是从普朗克函数在不同温度区间做线性近似拟合出来的用错了区间在边缘温度上会偏。在 ENVI 里我的做法是先算 C 和 D 两个中间量各存一个文件再把主式写成一步。为什么不分步算更多中间量因为 1-C-D 这个量在分子里出现了两次算一次存下来能省一次运算而且万一分母的 C 出了负数或者零你从中间文件能立刻看出问题——比如 ε 被误赋成了 0C 就是 0主式会直接除零ENVI 给出的结果是一堆极大值图上是一片刺眼的白。C 的表达式b1 * b2b1 是比辐射率b2 是大气透射率。D 的表达式(1 - b2) * (1 (1 - b1) * b2)5.2 单位一致性一个能白白吃掉 2K 的坑主式里的 T6 和 Ta必须和 a、b 的单位区间一致也就是都用摄氏度。但亮温算出来天然是开尔文所以要先减 273.15。我用基准参数验证过这件事T6 300.15K、Ta 295.15K如果直接代进公式而不是先转成 27℃ 和 22℃得到的结果是 302.14K也就是 28.99℃而正确的答案是 26.96℃。差了整整 2.03K。这个误差不会让图看起来有任何异常只会让所有数值系统性偏高两度非常危险。反过来如果你用的是某个只给开尔文版本的公式有些文献会把公式和系数一起改成开尔文形式那就全程用开尔文别再减 273.15。关键是公式和系数、还有输入单位三者必须是一套的混搭就是找死。5.3 越界值、无效值和掩膜处理主式跑完先做两件事。第一件是看统计量地表温度应该在什么范围中纬度夏季白天植被覆盖区大约 20 到 35℃裸露地表和建筑屋顶到 45℃ 甚至更高都有可能水体在 20 到 28℃。如果你看到负值、或者看到 100 以上的数值别犹豫回头查参数。第二件是处理无效值。主式计算过程里会有一些像元的输入本身就是填充值或者被云污染的极端值这些像元算出来的结果会落在物理上不可能的范围。处理方式是先建掩膜再做统计而不是直接把异常像元删掉——保留原始结果用掩膜控制哪些像元参与统计和出图这样你随时能回溯。掩膜的来源可以是 QA_PIXEL 波段解析出来的云、云影、水体标志位也可以自己用亮温阈值加 NDVI 阈值组合一个。我的习惯是两条路都做取交集宁可保守一点也不要把云边缘的半透明像元放进结果里。5.4 结果转摄氏度、配色和归档单窗算法的输出本身就是摄氏度前提是你按 5.2 说的做了如果你用的是开尔文版本的公式最后减 273.15 得到摄氏度。ENVI 里转一次很方便b1 - 273.15出图配色不要用彩虹色带。温度是连续变量用彩虹色带会在数值梯度平缓的地方造出假的边界看图的读者容易误判。用单色系的渐变色带或者分位数断点都更稳妥。同时记得在图上标出像元值的取值范围不要只给一个色带。归档这件事我特别想强调。ENVI 的 Band Math 表达式在项目切换、软件升级之后经常找不回来我踩过这个坑。现在我的做法是所有表达式写在一个 txt 文件里按步骤编号所有中间结果存成 ENVI 格式文件名里带上关键参数比如LST_swa_tau0.85_Ta22_eps_ndvi.dat同一个项目建一个 README写清楚数据版本、MTL 系数、气象数据来源。半年后要复现十分钟就能跑通。6. 结果对不对三条校验路径和一张异常速查表6.1 和官方温度产品对拍最快的一道体检Landsat 的 Level-2 产品里直接带了地表温度波段 ST_B10量化方式是 DN × 0.00341802 149.0结果单位是开尔文另外还有一个 ST_QA 波段给出每个像元的精度估计。这个产品是用辐射传输加再分析资料大气廓线算出来的跟单窗算法的原理不完全一样但都是同一个传感器、同一时刻的数据拿来对拍非常方便。我的做法是把官方 ST_B10 转到摄氏度和我的单窗结果做逐像元差值看差值的均值和标准差。正常情况下两者应该差不多差值均值在 1K 以内标准差 1.5K 左右散点分布在植被、裸土、城区比辐射率估算的差异会体现出来。如果差值均值超过 3K那基本可以确定是我的某个参数出了问题而不是两种方法本来就有差异。要注意 ST_QA 这一步不能省。官方产品里质量差的像元本身误差就大拿这些像元去对拍等于用坏尺子校准好尺子。6.2 和 MODIS LST、地面观测对比的正确姿势和 MODIS 地表温度产品对比有两个必须处理的问题。第一是时间匹配Terra 星过境时间和 Landsat 接近都是上午十点半左右用 Terra 的 MOD11A1 日产品比用 Aqua 的 MYD11A1 更合理。第二是空间尺度MODIS 是 1 公里Landsat 是 30 米直接逐像元比毫无意义。正确做法是先把 30 米的结果按 1 公里网格求平均升尺度再和 MODIS 像元比。同时用 MODIS 的 QC 层筛掉质量等级差的像元。和地面观测对比是最容易出错的一环。常规气象站测的是 1.5 米高度的气温和地表温度是两个物理量白天晴朗天气下两者差 10K 以上很正常。拿气温验证地表温度得到的误差里绝大部分是物理量本身的不同不是反演误差。真正可用的地面数据是红外测温仪测的表面辐射温度、或者长波辐射表反演的地表温度而且要注意观测的视场角、表面类型和像元尺度是否匹配。如果手上只有常规气象站数据那就老老实实只做相对分析不要报绝对精度。6.3 异常现象速查表下面这张表是我自己攒的按现象反查原因比从头查参数快得多现象最可能的原因排查方向整体偏高 2K 左右单位混用开尔文和摄氏同时进入主式检查 T6、Ta 是否都已减 273.15整体偏高 5K 以上比辐射率被赋成常数且偏低检查 ε 图层的取值范围整体偏低定标系数用了旧版本 Collection 的参数核对 MTL 里的 MULT/ADD水体异常偏热水体像元 ε 赋成了陆地值检查 NDVI 负值区的处理逻辑城区异常偏冷混合像元 dε 修正项漏掉检查线性式里的 0.0668 系数边缘一圈高值填充像元 DN0 参与计算加无效值掩膜结果图上出现条带输入波段本身有条带或掩膜未对齐检查原始数据质量和各图层网格一致性研究区一半是负值NDVI 计算时波段映射错位核对红光和近红外的波段编号倒数第二条值得展开说一句。我在做多景时间序列的时候遇到过一次某一年的结果上有规则条带查了半天最后发现是当年那景影像本身在传感器异常期采集的数据质量就有问题。所以在做趋势分析之前一定要按景看数据质量别把传感器问题当成地表变化。6.4 什么情况下不该硬用单窗算法单窗算法有它明确的假设前提地表近似朗伯体、大气在垂直方向均匀、天气晴朗无云。这几个假设在某些场景下会明显失效。地形起伏剧烈的山区坡度和坡向会改变局部入射角和大气路径长度单窗算法里没有对应的修正项结果会有系统性偏差通常表现为阳坡偏冷、阴坡偏暖。如果研究区是山区要么换用带地形修正的方法要么至少用 DEM 做一个分区分析别把整片山当成平地处理。水汽含量很高的地区比如夏季的南方τ 会比较低而 τ 越低公式对 τ 误差的敏感性越高前面那张灵敏度表里 0.05 的 τ 误差对应 0.44K如果 τ 本身只有 0.6这个放大效应会更明显。这种情况下要么把水汽估计做得更细比如用逐像元的水汽产品而不是单一值要么考虑换用对水汽不那么敏感的方法。还有一种情况是数据本身的限制。如果你手上只有单热红外波段那单窗算法几乎是唯一选择前提是把参数做实。如果你手上有两个热红外波段并且质量可靠劈窗算法在理论上更有优势因为它用两个通道的差异直接消掉了一部分大气影响减少了对辅助大气参数的依赖。但要注意前提是质量可靠——Landsat 8 的第 11 波段就不满足这个前提硬做劈窗反而更差。我自己在实际操作中的体会是单窗算法的门槛不在数学在数据。把 MTL 系数核对清楚、把比辐射率按像元算出来、把水汽含量找对来源、把单位统一好这四件事做到位结果基本就能用。剩下那部分误差更多来自数据本身的时空代表性不是靠调参数能解决的。每次跑完做一遍敏感性测试知道自己的结果对哪个参数最敏感比反复调参数追求看起来合理的图要靠谱得多。
