企业官网建设流程全解析

这个系列更新到第五篇了。前几篇我们把遥感影像的基础处理讲了一遍，包括波段合成、辐射定标和大气校正，今天这篇正好进入大家最关心的“算指数”环节——用一景Landsat 5 TM影像，一次性把NDVI、MNDWI、NDBI三个指数和地表温度反演全部跑通。如果你正在做土地覆盖分类、城市热岛、水环境污染或者植被长势监测，这四个参数基本就是你日常工作中绕不开的基础功，学会之后可以平移到Landsat 7、Landsat 8，甚至Sentinel-2上。

这篇我不只列公式，还会把每一步为什么要这么做、一般会踩什么坑都写清楚，适合有一定遥感基础、但还没完整跑通过这个流程的同学参考。文章会从原理、数据准备、逐个指数的实操计算，一路讲到地表温度反演的完整链路，最后把我在实际项目里遇到过的典型问题和排查思路一并整理出来。

1. 项目需求拆解：三个指数和LST到底在算什么

1.1 为什么这个流程偏偏选Landsat 5 TM

先说数据选择的问题。Landsat 5从1984年发射，到2011年退役，2012年才正式停止数据采集，服役时间接近三十年，是全球对地观测领域里运行时间最长的卫星之一。它搭载的TM传感器虽然属于上世纪八十年代的技术，但波段设计和现在的Landsat 8/9有很好的延续性，这就在长时序研究中形成了一个难得的优势：我可以用Landsat 5做2000年前后的历史影像分析，再用Landsat 8接着做现在的部分，两代数据之间的指数结果可以互相比较，不需要担心波段错位。

TM传感器一共有7个波段，波段1到波段5和波段7的空间分辨率是30米，波段6是热红外波段，原始分辨率为120米，但在标准产品中通常会被重采样到30米。具体波段参数如下：

这三个指数加一个温度反演，其实都建立在一个朴素的光谱规律上：不同地物在不同波段上的反射/辐射特性不一样。植被在近红外高反射、红光低反射；水体在绿光高反射、近红外和短波红外几乎全吸收；建筑物则在短波红外反射较强。Landsat 5 TM恰好把这三个关键波段都覆盖到了，所以可以一次算全。

1.2 四个参数各自的计算逻辑和用途

NDVI、MNDWI、NDBI这三个都是归一化差值类指数，公式结构看起来很像，但背后的地物光谱逻辑完全不同。

**NDVI（归一化植被指数）**用于反映植被覆盖和生长活力。它的逻辑是：健康植被在红光波段被叶绿素大量吸收，在近红外波段被叶片细胞结构强烈反射，所以用这两个波段的差除以和，可以把“植被强”的像元顶到高值。范围理论上在-1到1之间，一般水体为负值，裸土接近0，稀疏植被在0.2-0.4，茂密植被可以超过0.6。

**MNDWI（改进归一化差异水体指数）**是在NDWI基础上的改进版本。传统NDWI用的是绿光和近红外，而MNDWI把近红外换成了短波红外，这样做的关键好处是：水体在短波红外的吸收比在近红外更强，而建筑物和土壤在短波红外的反射比较高，所以MNDWI对建筑物的抑制作用比NDWI明显得多。城市水体提取的时候，这个改进几乎是决定性的。

**NDBI（归一化差异建筑指数）**用于识别城市建成区和不透水面。它的思路正好和NDVI相反：建筑区在短波红外反射强、近红外反射弱，所以用（SWIR1-NIR）/（SWIR1+NIR）可以把建筑区域凸显为正的高值，植被和清洁水体则表现为负值。但要注意，NDBI单独使用的误提率其实不低，裸地和部分干燥土壤在NDBI里也容易是正值，所以实际项目中最好把NDBI、MNDWI、NDVI放在一起做规则判断，而不是只用一个指数就下结论。

**地表温度反演（LST）**是所有遥感反演参数里对数据质量最敏感的之一。Landsat 5 TM的波段6接收的是地表和大气共同作用后的热红外辐射，要从这里面把“真实地表温度”解算出来，就必须想办法去除大气的影响，同时准确估计地表比辐射率。这个过程比前面三个指数复杂得多，所以我会在第四章专门展开。

2. 数据准备与预处理：这一步做不好，后面全白搭

2.1 数据获取、波段组成与头文件读取

做这个项目的起点是拿一景合适的Landsat 5 TM数据。我这次用的是2007年夏季某地区的一景L1T产品。L1T产品已经做过几何校正和地形校正，对大多数区域分析场景来说，几何精度足够了，不需要再做配准。

从官方数据源下载下来的压缩包解压之后，会看到两类核心文件：一类是各波段的TIFF影像，文件名一般形如“LT05_L1TP_xxx_20070701_20161201_01_T1_B1.TIF”，B1到B7分别对应波段1到波段7；另一类是MTL开头、后缀为txt的头文件。这个MTL文件非常关键，辐射定标需要的增益、偏置参数、太阳高度角、成像时间都记录在里面。我习惯先把MTL文件打开看一遍，确认影像的云量、采集日期和定标参数是否齐全。只有Level-1以上级别的产品才带完整的定标参数，如果拿到的是原始Level-0数据，那就得自己去查USGS官方发布的定标系数表，麻烦很多。

拿到数据后我通常还会先做一步“目视检查”：把影像用假彩色合成（RGB分别放波段4、3、2）看一眼，确认研究区有没有大片云层覆盖、有没有条带或坏行。如果云量大，与其后面做一堆无效计算，不如直接换一景数据。云不仅会让NDVI偏低，还会让地表温度反演结果产生严重的伪高温区。

2.2 辐射定标和大气校正到底要不要做

这是整个流程里争议最大的前置步骤。我的结论是：做指数计算，必须做辐射定标，但大气校正视情况而定。

辐射定标是必须的。因为TM影像记录的DN值只是传感器响应的数字量化值，不同时间、不同太阳条件下，同一个地物的DN值可能相差很大。定标就是把DN值转换成具有物理意义的表观反射率（可见光和近红外波段）或辐射亮度（热红外波段）。如果跳过这一步直接算NDVI，虽然比值在一定程度上能抵消大气和太阳高度的影响，但算出来的数值在不同影像之间没有可比性，后期要做时序分析或阈值跨影像使用时会出大问题。

大气校正则看用途。NDVI、MNDWI、NDBI这种比值型指数，因为分子分母同时受大气影响，一定程度上能够自我抵消，所以很多研究中直接在表观反射率上算，也被主流期刊接受。但如果你要把指数结果用来做定量分类、阈值跨区域迁移，或者要和实测光谱数据对比，那还是建议做一遍完整的大气校正。ENVI里的FLAASH模块、QUAC模块都可以用，前提是先完成辐射定标，将数据转换为反射率，然后才能运行大气校正模型。

需要特别提醒的是：热红外波段一般不要用FLAASH那套可见光大气校正流程去处理。热红外波段的大气效应主要体现为大气上行辐射、下行辐射和大气透过率对地表辐射信号的削弱与叠加，处理方式是不同的。这个我在第四章会专门说明。

2.3 影像裁剪、坏值处理与云掩膜

预处理里最容易忽略的是坏值和云掩膜。

下载的Landsat影像边缘通常有值为0的区域，这些0值像元在参与指数计算时会被当成真实的极低反射率，从而在归一化指数里产生大量-1附近的异常值。如果不先进行掩膜，最终成图会出现一片诡异的黑色边框，还会把统计直方图拉歪。我通常会在计算之前把所有0值像元统一设置为NoData，在ENVI里可以通过Build Mask工具，在ArcGIS里则用SetNull函数处理。

云掩膜我建议根据影像的QA波段来做。Landsat 5的L1T产品中有一个质量评估波段，记录了云、云阴影、雪等像元标记。实际操作中也可以做一个简单的“蓝波段阈值法”：因为云在蓝光波段的反射率非常高，把Band 1大于某个经验阈值（比如0.2）的像元标记为云，效果通常也不错。这个方法虽然粗糙，但在没有QA数据的情况下足够实用。

3. NDVI、MNDWI、NDBI的实操计算

3.1 三指数公式的统一记忆法

很多初学者一次性记三个公式容易混，我提供一个理解之道。这三个指数的分子分母结构都是“（特征波段A - 特征波段B）/（特征波段A + 特征波段B）”，关键是记住每个地物在哪两个波段表现得“反差最大”。

NDVI表示植被，特征是红光弱、近红外强，所以公式是（B4-B3）/（B4+B3）。

MNDWI表示水体，特征是绿光反射中等、短波红外极弱，所以公式是（B2-B5）/（B2+B5）。

NDBI表示建筑，特征和植被正好相反，短波红外比近红外强，所以公式是（B5-B4）/（B5+B4）。

你可以看到，MNDWI和NDBI用的都是B5波段，只是一个和B2做差，一个和B4做差。理解了光谱逻辑，就不需要死背公式，遇到Landsat 8或Sentinel-2的时候，也能自己推导出对应波段的指数表达式。

3.2 ENVI Band Math与ArcGIS栅格计算器实现

我最常用的还是ENVI的Band Math工具，因为它在处理多波段影像时更顺手。假设我已经在ENVI里打开了一景完成定标和裁剪的TM影像，波段顺序和原始文件一致：b1对应Band 1，b2对应Band 2，以此类推。

计算NDVI时，输入表达式：

(float(b4)-b3)/(float(b4)+b3)

有人会问，为什么要写float(b4)？这是因为Landsat 5 TM的原始波段是16位整型。整型除以整型，在某些软件环境里结果还是整型，导致大量介于0和1之间的小数被直接截断成0，最后还是变成0和1两个极端值，图像一片花白。这个坑我早年踩过，浪费了一整个晚上。后来养成了习惯，涉及比值计算一律先转float，并且顺手把b3也一起转float，避免运算中不同类型变量参与的隐式转换问题。

同样的道理，MNDWI的表达式：

(float(b2)-b5)/(float(b2)+b5)

NDBI的表达式：

(float(b5)-b4)/(float(b5)+b4)

在Band Math里把对应波段加载为变量后，点OK就能得到结果。输出的时候建议选择浮点型。完成后用Link工具把三个指数结果和原始假彩色影像联动浏览，一眼就能看出哪些地物被正确识别出来了。

如果你更习惯用ArcGIS，在栅格计算器里写类似的公式，记得用Float()函数包裹波段名，例如NDVI：

Float("Band_4") - Float("Band_3") / Float("Band_4") + Float("Band_3")

等等，这里必须加括号。正确写法是：

(Float("Band_4") - Float("Band_3")) / (Float("Band_4") + Float("Band_3"))

ArcGIS栅格计算器对运算优先级是敏感的，少一个括号结果完全变样。我见过有同事把NDVI算成“B4 - B3 / B4 + B3”，出来的结果和真实NDVI差了十万八千里。公式本身不复杂，但写进计算器的时候一定多检查两遍括号。

3.3 GEE和R的轻量级替代方案

如果你的研究区范围较大，或者需要一次性处理多期影像，桌面软件一个个算效率太低，这时候Google Earth Engine和R是更合适的选择。

GEE里如果直接用Landsat 5表面反射率产品，三行代码就能得出三个指数：

var l5 = ee.Image("LANDSAT/LT05/C01/T1_SR") .filterBounds(roi) .filterDate('2007-05-01', '2007-09-30') .first(); var ndvi = l5.normalizedDifference(['B4', 'B3']).rename('NDVI'); var mndwi = l5.normalizedDifference(['B2', 'B5']).rename('MNDWI'); var ndbi = l5.normalizedDifference(['B5', 'B4']).rename('NDBI');

GEE会替你把辐射定标和大气校正都处理好，直接使用表面反射率产品。但要注意，GEE里也有一个问题：Collection 1和Collection 2的表面反射率产品，波段缩放系数不一样。Collection 2的SR产品默认是原始整数，使用前需要乘上0.0001，否则指数会被放大一万倍。我建议养成习惯，在提交数据之前先print看一下波段值和属性。

R语言侧重数据分析和可视化，配合raster或terra包也能完成类似的工作。思路是先读取各波段为栅格对象，再用overlay函数逐像元计算指数。R的好处在于后续统计分析一条流水线走完，比如算完NDVI马上做植被覆盖度分级、统计面积，都不用换工具。

4. 地表温度反演全流程详解

4.1 方法选型：辐射传输方程法还是单窗算法

地表温度反演的方法很多，但对Landsat 5 TM来说，实际应用中主要就两条路线：辐射传输方程法（也叫大气校正法）和覃志豪的单窗算法。

辐射传输方程法的思路非常直观：传感器接收到的热红外辐射亮度，由三部分组成——地表自身辐射经过大气衰减后的部分、大气上行辐射、大气下行辐射经地表反射后的部分。把后两个“大气污染”减掉，再除以大气透过率，就得到地表真实辐射亮度。这个方法需要三个大气参数：大气透过率τ、大气上行辐射L↑、大气下行辐射L↓。

单窗算法是国内遥感圈非常熟悉的方法，只需要大气透过率和大气平均作用温度两个参数，计算更加简便。它的形式看起来稍复杂，但原理类似，都是在大气参数支持下把大气影响从星上辐射中剥离出来。

我的习惯是优先用辐射传输方程法，因为大气参数可以从遥感大气校正参数计算网站直接获取，整个过程更透明、更好复现。如果研究对象是大区域多期影像、每一景都去网站查参数不现实，那就退一步用单窗算法，配合MODIS水汽产品估算大气透过率。下面我以辐射传输方程法为主线，把每一步的操作细节讲透。

4.2 热红外波段辐射定标与亮度温度计算

第一步，对TM的波段6做辐射定标，得到热红外辐射亮度L6，单位是W/(m²·sr·μm)。在ENVI里直接用Radiometric Calibration工具，选择波段6，设置输出类型为Radiance即可。如果是自己手动计算，公式是：

L6 = gain * DN + bias

这里的gain和bias在MTL头文件里都有，不同影像会略有差异，本质是传感器定标参数随时间和卫星轨道变化而更新的结果。需要留意的是，Landsat 5在2004年之后有一套新的定标系数，和早期不同。如果你用的是旧版本软件或者网上找到的旧教程里的固定系数，可能会算出偏差。最稳妥的办法永远是从MTL文件里读当前影像的定标参数。

第二步，计算星上亮温T6，也就是“黑体等效温度”。公式是普朗克公式的反函数：

T6 = K2 / ln(K1 / L6 + 1)

对Landsat 5 TM，K1 = 607.76 W/(m²·sr·μm)，K2 = 1260.56 K。这两个常数和Landsat 8是不同的。很多人做完Landsat 8的项目直接切到Landsat 5，忘记换常数，结果温度差了十几K，这是非常典型的低级错误。

用ENVI的Band Math可以这样写：

1260.56 / alog(607.76 / b6 + 1)

注意Temperature输出默认是开尔文，后续如果要转摄氏温度，还要再减273.15。

4.3 地表比辐射率的三种估算方式

地表比辐射率是反演精度影响最大的地表参数，但也是最难精确获得的一个。在没有实测波谱数据的情况下，最常用的方法是借助NDVI估算，思路是把每个像元近似看成“裸土”和“植被”的混合，根据植被覆盖度线性推算比辐射率。

先计算植被覆盖度Fc：

Fc = (NDVI - NDVI_min) / (NDVI_max - NDVI_min)

这里的NDVI_min和NDVI_max可以取经验值0.05和0.7，也可以统计整景影像NDVI直方图的5%和95%分位数来代替，后者更贴合具体情况。实测下来，用影像自己的分位数比固定经验值更靠谱，因为不同季节、不同区域的NDVI分布差异不小。

得到Fc之后，按下面规则估算比辐射率：

• 当NDVI小于0时，判为水体像元，比辐射率取0.995；
• 当NDVI在0到0.7之间时，比辐射率取0.004 * Fc + 0.986；
• 当NDVI大于0.7时，判为完全植被像元，比辐射率取0.986。

这个规则来自Sobrino等人的研究，在大多数中低纬度研究区表现良好。也可以简化为三个固定类型：水体0.995、城镇0.970、自然地表0.986。如果你研究的是城市区域，我建议把城镇像元的比辐射率单独设置成0.970附近，而不是统一用0.986，因为水泥、沥青、屋顶材料在热红外波段的发射特性确实比植被低不少，统一的取值会低估城市像元的地表温度。

4.4 大气参数求解与最终温度合成

在地表偏好的估算和热红外波段的辐射亮度都到手之后，下一步获取三个大气参数。目前最方便的做法是登录NASA的大气校正参数计算网站，输入成像时间、影像中心经纬度、高程、大气模式等信息，系统会返回对应的τ、L↑和L↓。一般选择中纬度夏季或冬季大气模式，水汽和温度廓线用标准大气即可。

如果因为项目需要无法在线查参数，可以用MODTRAN典型大气的近似值。我列出几组常用参考值（以Landsat 5 TM的band 6为例）：

注意这些只是典型情况下的参考，不是所有研究区都适用。精度要求高时，还是应该用对应成像时刻的大气剖面产品或在线计算工具。

有了ε、τ、L↑、L↓、L6之后，先计算同温度下黑体辐射亮度：

B(Ts) = (L6 - L↑ - τ * (1 - ε) * L↓) / (τ * ε)

再代入普朗克反函数：

Ts = 1260.56 / ln(607.76 / B(Ts) + 1)

在ENVI Band Math里可以直接把这个公式写成一整条表达式：

1260.56 / alog(607.76 / ((b6 - 1.79 - 0.71 * (1 - e) * 2.76) / (0.71 * e)) + 1) - 273.15

这里的e就是上一步算出来的比辐射率栅格，b6是热红外波段的辐射亮度。减掉273.15之后，结果直接就是摄氏度。

算完之后一定要做一次简单的合理性检查。正常情况下，中等纬度夏季晴天的地表温度应该在15到50摄氏度之间。如果发现大范围出现70度以上的值，通常是比辐射率取值偏低、大气透过率偏大或者辐射定标单位不对；如果全是0到5度的低温值，则要检查是不是K1/K2常数用成了Landsat 8的，或者是不是忘记做辐射定标直接拿DN值代入公式了。

5. 常见问题与排查技巧实录

5.1 指数异常值：八成出在数据类型和NoData上

做指数计算时最常遇到的现象就是结果图边缘有一圈黑色、或者整个影像内出现大量-1和1的极值。遇到这种情况，先不要怀疑算法，先查数据类型。

有人把DN值没定标就直接丢进Band Math，是整型，相除之后小数全被截断，结果图自然不是0就是1。解决办法就是我前面强调的，表达式里用float()包一下再算。还有一类情况是影像本身带着NoData值或者0值像元，归一化指数的分母一旦遇到0值还是0，会出现除零错误。在ENVI中，零值像元参与运算后通常会变成无穷大或-9999之类，后续统计时会严重影响直方图。我习惯的做法是在计算指数之前，先做一次无效值掩膜，或者在指数计算完成后用条件函数把异常值设为NoData。

另外提醒一句，如果你在ArcGIS里做，先把原始影像复制出一个临时文件，把NoData值显式设置为某个数值（比如-9999），计算完再设置回来。这一步虽然琐碎，但能让后面的分类和制图省心很多。

5.2 温度反演结果偏差的排查清单

地表温度反演出问题，通常不是“某个环节错了”，而是多个参数共同造成的系统性偏差。我整理了一个排查顺序：

先检查输入的单位。辐射定标后L6的单位必须是W/(m²·sr·μm)，如果定标工具输出了W/(m²·sr·μm·nm)或者别的单位，数值会差好几个数量级，后面的计算就全乱了。

再检查K1、K2常数值。这是最容易被忽视的一步。每位从Landsat 8教程转到Landsat 5的人，几乎都踩过这个坑。Landsat 8的K1=774.89、K2=1321.08，Landsat 5是607.76和1260.56，完全不是一回事。

然后对比辐射率取值做敏感性测试。比如把研究区所有像元的比辐射率从0.98改成0.97，温度结果会升高约1到1.5度。如果你发现自己反演的温度比实测气象站数据整体偏高5度，很可能就是因为比辐射率设低了。

最后检查大气参数。在线工具拿到的大气参数是针对单一大气的，如果研究区面积非常大、跨越了不同气候区，最好分区域分别获取参数，而不是全图用同一组值。

5.3 结果验证与阈值选择的经验

指数和温度反演做完之后，一定要做验证，否则论文或报告里底气不足。

NDVI、MNDWI、NDBI这类指数，最简单的验证方式是在高分辨率影像上随机选取典型的植被、水体、建筑样本点，统计这些点在指数影像上的分布范围，确认不同地物的指数区间是否存在明显的可分性。我一般会选每类30到50个样本点，计算均值和标准差。如果某两类地物的指数直方图完全重叠，说明这个指数在你的研究区不适合单独使用，需要引入其他波段信息。

地表温度验证的常见做法是和气象站地表温度观测数据做线性回归分析，统计相关系数和均方根误差。要注意对比口径的问题：气象站测的是站点周边2米或10米的气温，遥感反演的是混合像元尺度的地表辐射温度，两者之间存在系统性差异，白天晴热天气下地表温度通常比气温高10度以上。所以对比时用回归校正而不是直接要求数值相等，均方根误差在2到3开尔文以内，就算比较理想。

至于阈值选择，我的经验是不要直接套用期刊论文里的固定阈值，应该针对研究区影像单独统计。操作方法是先画出指数直方图，观察双峰或多峰分布，再以峰谷位置作为初始阈值，结合高分辨率样本点反复调整。不同时相、不同季节、不同气候区的阈值很难通用，与其纠结“最科学”的阈值，不如把你的样本点检验结果写清楚，用可重复的统计过程来支撑阈值设定。

这组流程跑下来，NDVI、MNDWI、NDBI和地表温度反演四个结果就能放在同一个工作空间里做后续分析。我自己做了几轮之后最大的感触是，技术流程本身并不复杂，真正拉开差距的地方在于对每个参数的来源和误差范围心里有数。比如植被覆盖度估算公式里的NDVI阈值，换一景影像可能就要改；大气参数没查准，温度结果虽然有趋势但绝对值不可信。建议大家养成把每一步的参数截图存档的习惯，尤其是MTL文件里的定标参数和在线获取的大气参数，这样后期改任何一步都不需要重新从头推导。接下来如果你要做城市热岛分析，可以把NDBI作为自变量、LST作为因变量做一个简单的线性回归，或者按缓冲区统计温度梯度，这些扩展都能直接从今天的结果接上。