Landsat9地表温度反演实战：ENVI Band Math与USGS数据的高效结合

1. 从“-100度”的坑说起：为什么你的Landsat9温度图会出错？

如果你最近在用USGS的Landsat Collection-2地表温度产品，配合ENVI的Band Math做反演，很可能遇到过一种让人摸不着头脑的情况：辛辛苦苦下载了数据，套用了官方公式，结果图一出来，大片区域显示温度是零下100多度。这显然不符合常理，北极圈内都没这么冷，更别说你研究的可能是长三角或者华北平原了。

我刚开始用Landsat9的Collection-2数据时，也踩过这个坑。当时我还以为是公式输错了，或者数据下载有问题，反复核对了半天。后来仔细研究产品说明和像元值才发现，问题出在数据本身。USGS提供的这个地表温度产品（ST_B10波段），里面是包含缺失值的。这些缺失值在数据文件中通常用一个特定的数字（比如0）来标记，表示该位置由于云层覆盖、传感器问题或缺乏必要的发射率数据等原因，无法计算出有效的地表温度。

然而，当我们直接使用那个经典的转换公式 LST = DN * 0.00341802 + 149 - 273.15 时，ENVI的Band Math会忠实地对每一个像元（包括那些标记为缺失值的像元）进行计算。你把DN值（像元值）为0代入公式：0 * 0.00341802 + 149 - 273.15 = -124.15°C。看，这就是那片诡异的超低温区域的来源。它并不是真实温度，而是公式对“无效数据”进行数学运算后产生的“垃圾值”。

所以，我们面临的核心挑战就变成了：如何在使用便捷的Band Math和官方预处理数据的同时，智能地识别并剔除这些缺失值，让最终的温度图只反映真实、有效的地表信息。这不仅仅是让图看起来更“好看”，更是保证后续分析，比如城市热岛效应评估、农业旱情监测等，结果准确性的关键一步。下面，我就把自己摸索出来的几套实战解决方案分享给你，从快速修复到深度处理，总有一款适合你。

2. 方案一：Band Math 增强版——一招屏蔽缺失值

最直接、最快上手的方法，就是在我们原有的Band Math公式里，加入一个逻辑判断条件。ENVI的Band Math语法支持类似编程的逻辑运算符，我们可以利用这个特性来“过滤”数据。

核心思路：在计算温度之前，先判断像元值是否在有效范围内。USGS的文档指出，ST_B10波段的有效DN值范围通常大于0。因此，我们可以设定一个条件：如果DN值小于或等于一个很小的阈值（比如0），就输出一个特定的“无效值”（例如-9999）；否则，才进行正常的温度转换计算。

具体操作步骤：

1. 打开Band Math工具：在ENVI中，依次点击 Basic Tools -> Band Math。

2. 输入增强版公式：在公式输入框中，键入以下代码：

   (b1 le 0) * -9999 + (b1 gt 0) * (0.00341802 * b1 + 149 - 273.15)
   

   我来拆解一下这个公式：

   • b1 le 0：判断波段b1（即ST_B10）的值是否小于等于0。如果成立，结果为1（真），否则为0（假）。
   • (b1 le 0) * -9999：如果上述判断为真（1），则乘以-9999，结果就是-9999；如果为假（0），结果就是0。这相当于把无效像元赋值为-9999。
   • b1 gt 0：判断b1是否大于0。
   • (b1 gt 0) * (0.00341802 * b1 + 149 - 273.15)：如果b1大于0，则进行温度计算，否则结果为0。
   • 最后将两部分相加。因为对于任何一个像元，(b1 le 0)和(b1 gt 0)必然只有一个为1，所以最终每个像元要么是-9999，要么是计算出的正确温度值。

3. 指定输入波段：点击“OK”后，会弹出变量绑定窗口。将b1指定为你加载的ST_B10波段数据。

4. 设置输出：指定输出文件的路径和文件名，点击“OK”执行计算。

效果与优化：执行完毕后，你得到的新图像中，原来的“-100多度”区域会变成-9999。你可以在ENVI中，通过右键图层选择“Edit Color Mapping”或“Edit Raster Color Slices”，将-9999这个值的显示颜色设置为透明或者你喜欢的背景色，这样在出图时它们就“消失”了。

进阶技巧：有时候，缺失值可能不是简单的0。你可以通过查看ST_B10的直方图（Quick Stats），或者查阅该景数据的元数据文件（MTL.txt），找到标记缺失值的具体数值（可能是0，也可能是其他值，比如65535）。然后只需将公式中的0替换成那个具体的缺失值标记即可。这个方法最大的优点是快，几乎不增加操作步骤，就能得到干净的可视化结果。

3. 方案二：利用QA波段——做更精细的数据“保洁”

方案一虽然快，但它有点“粗放”，它只是基于DN值是否大于0来做判断。实际上，Collection-2产品提供了一个强大的“质检”工具——QA_PIXEL波段。这个波段里存储了每个像元丰富的质量信息，比如是否是云、云阴影、雪、冰、水体，以及置信度等级。利用QA波段，我们能进行更专业、更精细的数据清洗。

理解QA_PIXEL：QA_PIXEL是一个16位的整型波段，每一位（bit）或每几位组合，都代表一种特定的质量属性。例如，根据USGS的文档，第3位（从0开始数）表示“云”，第4位表示“云阴影”。如果某像元的“云”位被置为1，就说明这个像元被算法判定为云。

操作流程：先掩膜，后计算

1. 创建云掩膜文件：

   • 打开你的Landsat9 L2SP数据，加载QA_PIXEL波段。
   • 再次打开Band Math工具。这次，我们要创建一个云和云阴影的掩膜。公式可以这样写：

     (b1 and 8) eq 0 and (b1 and 16) eq 0
     

     • b1 and 8：8是二进制1000（第3位为1），这个操作是检查“云”位是否为1。
     • (b1 and 8) eq 0：如果“云”位是0，表示不是云，结果为真（1）。
     • 同理，(b1 and 16) eq 0检查“云阴影”位（16是二进制10000，第4位）。
     • 整个公式的意思是：当像元既不是云，也不是云阴影时，输出1（好像元），否则输出0（坏像元）。
   • 将b1绑定到QA_PIXEL波段，生成一个二值掩膜图像（值只有0和1）。

2. 应用掩膜计算温度：

   • 现在，我们结合掩膜和ST_B10来计算温度。输入一个新的Band Math公式：

     (b2 eq 1) * (0.00341802 * b1 + 149 - 273.15) + (b2 eq 0) * -9999
     

     • b1：绑定为ST_B10波段。
     • b2：绑定为上一步生成的云掩膜图像。
     • 公式逻辑：如果掩膜b2等于1（好像元），就计算温度；如果等于0（云或云阴影），就赋值为-9999。

这个方法的优势：它不仅仅处理了原始缺失值，还主动剔除了被云层覆盖的区域。云层顶部的温度非常低，会严重干扰真实的地表温度分析。通过QA波段掩膜，你得到的是晴空条件下的地表温度，这对于绝大多数生态、环境研究来说，数据质量更高，结论也更可靠。你可以根据需要，修改掩膜公式，把雪、冰（位5）或者低置信度像元也剔除掉，实现定制化的数据清洗。

4. 方案三：回归辐射传输方程——应对无ST产品的场景

前面两个方案都基于一个前提：你下载的是包含ST_B10波段的L2SP产品。但有时候，你下载的数据可能只是L2SR产品（只有地表反射率，没有地表温度），或者你希望从最原始的辐射值开始，完全自己控制反演流程，这时就需要回归到经典的辐射传输方程法。好消息是，即使过去常用的大气参数网站（如atmcorr.gsfc.nasa.gov）不可用，我们依然有办法。

新思路：从L2SP数据中“借用”大气参数 Collection-2的L2SP产品包里，除了ST_B10，还贴心地提供了反演过程中用到的几个关键中间波段：

• ST_ATRAN：大气热红外波段透过率（τ）
• ST_URAD：大气上行辐射亮度（L↑）
• ST_DRAD：大气下行辐射亮度（L↓）

完整操作步骤：

1. 数据准备：下载一景Landsat9的L2SP数据并解压。你需要用到以下文件：

   • *_SR_B*.TIF：地表反射率波段（用于计算植被指数和比辐射率）。
   • *_ST_ATRAN.TIF
   • *_ST_URAD.TIF
   • *_ST_DRAD.TIF
   • *_MTL.txt：元数据文件。

2. 计算地表比辐射率（ε）：

   • 使用反射率波段（如B4红波段，B5近红外波段）计算NDVI。
   • 通过NDVI估算植被覆盖度（Fv）。一个常用的经验公式是：Fv = ((NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil))^2，其中NDVI_soil可取0.2，NDVI_veg可取0.86。
   • 根据地物类型（水体、自然表面、城镇）计算比辐射率ε。一个广泛使用的模型是：对于自然植被覆盖区，ε = 0.004 * Fv + 0.986。

3. 处理大气参数并计算温度：

   • 将ST_URAD、ST_DRAD、ST_ATRAN三个文件导入ENVI。
   • 关键步骤：单位换算。根据产品指南，ST_URAD和ST_DRAD的栅格值需要乘以0.001（即除以1000），ST_ATRAN需要乘以0.0001（即除以10000）。你可以在导入时用Band Math处理，例如对ST_URAD：b1 * 0.001。
   • 现在，运用辐射传输方程的核心公式，在Band Math中输入：

     (b2 - b3 - b4 * (1 - b1) * b5) / (b4 * b1)
     

     • b1：绑定为你计算出的地表比辐射率ε图像。
     • b2：绑定为热红外波段辐射亮度。对于Landsat9 TIRS Band 10，你需要从MTL文件中找到辐射定标参数，将原始DN值转换为辐射亮度值（单位：W/(m²·sr·μm)）。公式通常是：Lλ = ML * DN + AL，其中ML和AL在MTL文件中。
     • b3：绑定为处理后的上行辐射亮度ST_URAD（已乘0.001）。
     • b4：绑定为处理后的大气透过率ST_ATRAN（已乘0.0001）。
     • b5：绑定为处理后的下行辐射亮度ST_DRAD（已乘0.001）。
   • 这个公式计算出来的是黑体在表面温度Ts下的热辐射亮度B(Ts)。

4. 亮度温度转地表温度：

   • 最后，利用普朗克公式的反函数，将B(Ts)转换为地表温度Ts。对于Landsat9 Band 10，公式为：

     (1321.08) / alog(774.89 / b1 + 1) - 273.15
     

     • 其中b1绑定为上一步得到的B(Ts)图像。
     • 1321.08和774.89是Landsat9 Band 10的K2和K1常数。
     • alog是ENVI中的自然对数函数。
     • 减去273.15是为了将开尔文温度转换为摄氏度。

这个过程虽然步骤多一些，但它给了你最大的灵活性和透明度。你可以控制每一个参数，并且这个方法不依赖于现成的ST产品，适用性更广。计算出的结果理论上不会有方案一中那种大片的缺失值区域，因为你是从辐射亮度开始一步步推导的。

5. 实战心得：数据下载、处理与结果解读的避坑指南

走通了技术流程，我想分享几个在实战中直接影响效率和结果的关键点，这些都是我亲身经历后总结的。

关于数据下载：在USGS EarthExplorer上找数据时，一定要认准“Collection 2 Level-2”这个标签。在数据集选择中，明确勾选“Landsat 9 OLI/TIRS C2 L2”。下载后，检查文件列表。如果包含*_ST_B10.TIF，那就是L2SP产品，可以采用方案一或二。如果只有*_SR_B*.TIF而没有ST相关文件，那就是L2SR产品，你需要采用方案三。我建议优先下载L2SP，省去很多计算大气参数的麻烦。

关于ENVI Band Math的使用技巧：

• 公式保存：在Band Math窗口输入复杂公式后，别急着点OK执行。先点“Add to List”，然后点“Save”，可以把公式保存成.exp文件。下次遇到同样的计算，直接“Restore”加载，只需重新绑定波段变量即可，效率倍增。
• 中间结果：像计算NDVI、植被覆盖度、比辐射率这些步骤，都会产生中间文件。建议给它们起个清晰的名字，比如Area_20230915_NDVI.img。虽然会占用一些磁盘空间，但万一后续步骤出错，你可以从中间步骤开始排查，不用全部重来。
• 结果检查：温度图像生成后，不要只看渲染图。一定要用“Cursor Location/Value”工具，在地图上移动光标，查看典型地物（如水体、森林、城市建筑、农田）的像元温度值是否在合理范围内。也可以用“Quick Stats”查看整景图像的最小值、最大值、均值，判断是否有异常值污染。

如何解读你的第一张地表温度图：当你得到一张色彩斑斓的温度分布图时，可以结合谷歌地球或其它高分辨率影像进行解读。通常你会发现：

• 水体（湖泊、河流）呈现明显的低温（蓝色调）。水的比热容大，升温慢，在白天通常比陆地温度低。
• 茂密植被（森林、公园）温度也相对较低。植物的蒸腾作用会消耗热量，起到降温效果。
• 城市建成区（商业区、工业区、密集住宅区）往往是大片的高温区（红色或黄色调），这就是“城市热岛效应”的直观体现。混凝土、沥青等建筑材料吸热快、储热能力强，加上人类活动排放的热量，共同导致了这一现象。
• 农田的温度介于自然植被和城市之间，具体取决于作物类型和灌溉情况。

记住，地表温度反演只是一个工具，更重要的是你如何利用这个工具去发现和解释空间现象背后的规律。无论是监测城市扩张对热环境的影响，还是评估农作物旱情，一张准确、干净的温度图都是你做出可靠分析的基石。希望这些从踩坑到填坑的经验，能让你在利用Landsat9和ENVI探索地表温度的世界时，走得更顺畅一些。