这个系列更新到第五篇了。前几篇我们把遥感影像的基础处理讲了一遍包括波段合成、辐射定标和大气校正今天这篇正好进入大家最关心的“算指数”环节——用一景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米。具体波段参数如下波段波长范围微米分辨率米主要用途Band 1 蓝0.45-0.5230水体穿透、叶绿素敏感Band 2 绿0.52-0.6030植被反射峰、水体识别Band 3 红0.63-0.6930叶绿素吸收、植被分类Band 4 近红外0.76-0.9030植被高反射、水体强吸收Band 5 短波红外1.55-1.7530植被含水量、土壤/建筑识别Band 6 热红外10.40-12.50120产品重采样30地表温度反演Band 7 短波红外2.08-2.3530矿物识别、干/湿区分这三个指数加一个温度反演其实都建立在一个朴素的光谱规律上不同地物在不同波段上的反射/辐射特性不一样。植被在近红外高反射、红光低反射水体在绿光高反射、近红外和短波红外几乎全吸收建筑物则在短波红外反射较强。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/SWIR1NIR可以把建筑区域凸显为正的高值植被和清洁水体则表现为负值。但要注意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/B4B3。MNDWI表示水体特征是绿光反射中等、短波红外极弱所以公式是B2-B5/B2B5。NDBI表示建筑特征和植被正好相反短波红外比近红外强所以公式是B5-B4/B5B4。你可以看到MNDWI和NDBI用的都是B5波段只是一个和B2做差一个和B4做差。理解了光谱逻辑就不需要死背公式遇到Landsat 8或Sentinel-2的时候也能自己推导出对应波段的指数表达式。3.2 ENVI Band Math与ArcGIS栅格计算器实现我最常用的还是ENVI的Band Math工具因为它在处理多波段影像时更顺手。假设我已经在ENVI里打开了一景完成定标和裁剪的TM影像波段顺序和原始文件一致b1对应Band 1b2对应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()函数包裹波段名例如NDVIFloat(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 TMK1 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估算思路是把每个像元近似看成“裸土”和“植被”的混合根据植被覆盖度线性推算比辐射率。先计算植被覆盖度FcFc (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↓中纬度夏季0.711.792.76中纬度冬季0.821.121.79热带0.662.213.31注意这些只是典型情况下的参考不是所有研究区都适用。精度要求高时还是应该用对应成像时刻的大气剖面产品或在线计算工具。有了ε、τ、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的K1774.89、K21321.08Landsat 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作为因变量做一个简单的线性回归或者按缓冲区统计温度梯度这些扩展都能直接从今天的结果接上。