做热红外遥感的人多多少少都被人问过一句这张图能不能直接给我一个地表温度。早年我也觉得这事不难不就是 DN 值转辐射亮度、再转亮度温度、最后套个公式嘛。真到自己拿 ENVI5.3 上手才发现从原始影像到一张靠谱的 LST 地表温度反演数据中间埋的坑比想象中多得多辐射定标漏了一个系数、比辐射率用了固定值 0.95、水汽含量拍脑袋填了个 1.5最后出图看着挺漂亮一验证跟地面实测差了四五度。这篇就说单窗算法这条路——它是国内做 Landsat 热红外温度反演最常用的方案之一参数少、公式封闭、在 ENVI 里纯靠波段运算就能跑完全程特别适合刚接触热红外的朋友和需要批量出图的人。我把参数怎么定、表达式怎么写、哪里最容易翻车尽量一次讲透。需要说明的是下面涉及的参数取值和程序步骤一部分来自公开文献里被反复引用的通用做法一部分是我自己在实操中的选择具体项目请按自己的研究区条件调整。1. 单窗算法凭什么成为Landsat热红外反演的常客1.1 从DN值到地表温度中间隔着三层先把物理链条捋直。卫星传感器在热红外通道记录到的其实是大气顶端的辐射亮度这个辐射亮度由三部分叠加而成地表自身的热辐射、大气上行辐射、以及被大气反射回传感器的下行辐射。你真正想要的只是第一项可传感器把三项一起打包送回来了。所以反演的本质就是剥洋葱。第一层剥掉传感器自身的响应特性得到大气顶的辐射亮度这一步叫辐射定标第二层要从辐射亮度反推传感器的入瞳亮温靠的是普朗克函数的反函数也就是那个K2/ln(K1/L1)的式子第三层才是最难的——把大气上行、下行辐射和地表比辐射率的影响一并剥离得到真正的地表温度。单窗算法的聪明之处就在第三层。它不去逐个求解大气上行下行辐射而是用一个简化的辐射传输近似把大气影响压缩成三个可测可算的参数大气透射率 τ、大气平均作用温度 Ta、地表比辐射率 ε。只要这三个数搞得定一个封闭公式就能直接算出 Ts。这也是它比辐射传输方程直接解更接地气的原因——后者需要同步的大气廓线数据多数人根本拿不到。我第一次做的时候卡在最朴素的地方以为亮温就是地表温度。后来拿亮温图和气象站数据一对比夏天中午差了七八度才明白亮温是传感器看到的温度包含了比辐射率小于 1 和大气削弱两个因素跟真实地表温度差了整整一个量级。1.2 单窗、劈窗、温度发射率分离该选哪条路热红外反演不是只有单窗一条路。常见的还有劈窗算法、多波段温度发射率分离TES、以及基于辐射传输方程的查找表法。各自的门槛和精度不一样选错了工具会白费很多力气。方法需要的数据主要参数适用情况大致精度单窗算法单个热红外波段ε、τ、Ta只有 Band 10 可用或只需单波段±1.0~1.5 K劈窗算法两个相邻热红外波段ε 差异、水汽有 Band 1011 或 Band 3132±0.8~1.2 K温度发射率分离多热红外波段波段间发射率关系MODIS 类多通道传感器依赖先验约束辐射传输方程单/多波段同步大气廓线有探空或再分析廓线参数敏感对 Landsat 8/9 来说情况有点特殊。TIRS 虽然有两个热红外波段但 USGS 官方建议不要用 Band 11 做地表温度定量反演因为它受杂散光影响明显定标精度不稳定。这个建议一出来等于把大多数人推回了单窗算法——你手上真正能信的只有 Band 10。如果只是做区域热环境分析、城市热岛分级、或者若干期的温度变化趋势对比单窗算法完全够用。它的问题不在精度上限而在参数取值的稳定性。劈窗算法虽然理论精度高一点但它需要两个波段的比辐射率差值而这个差值恰恰是热红外里最难确定的量。在 Landsat 8 只有 Band 10 可靠的前提下纠结劈窗意义不大。1.3 开干之前先确认数据够不够这一步经常被跳过然后做到一半发现材料不全。清单如下热红外波段影像Landsat 8 的 Band 10必须是 L1TP 级别的原始产品附带的 MTL.txt 文件要完整里面存着定标系数。光学波段影像红波段与近红外波段用来算 NDVI进而推植被覆盖度 Pv 和比辐射率 ε。成像时间与中心经纬度估算大气水汽和大气平均作用温度都要用。近地面气温反演当天、成像时刻附近的气温用于计算大气平均作用温度。相对湿度或可降水量用于估算大气透射率。气温和湿度这两个数据气象站资料、再分析数据集、或者公开的气象数据服务都能拿到。我建议至少取到成像时刻前后一小时的观测值别直接拿当天日均气温凑数——单窗算法对 Ta 的敏感性不低日均值和时间瞬时值能差三五度。另外提醒一句Landsat 8 在 2017 年前后有过一次产品重处理定标系数有微调。拿到影像后一定要打开 MTL.txt 核一遍RADIANCE_MULT_BAND_10和RADIANCE_ADD_BAND_10到底是多少别默认套用记忆里的 0.0003342 和 0.1。2. 三个大气参数和两个地表参数到底怎么定2.1 亮度温度从DN到开尔文的完整链条亮度温度是整个反演的中间量算错了后面全废。链条分两步。第一步DN 值转大气顶辐射亮度L RADIANCE_MULT_BAND_10 * DN RADIANCE_ADD_BAND_10第二步辐射亮度转亮温。这一步用的是热红外定标常数T6 K2 / alog(K1 / L 1)ENVI 的波段运算里自然对数写作alog不是ln这是新手最容易踩的一个书写陷阱。Landsat 8 Band 10 常用的 K1 774.8853 W/(m²·sr·µm)K2 1321.0789 K算出来的 T6 单位是开尔文。ENVI 里的表达式长这样b1 是定标后的 Band 10 辐亮度1321.0789 / alog(774.8853 / b1 1)结果记得减去 273.15 才是摄氏度。但注意减 273.15 这一步建议放到最后因为单窗算法公式内部的 T6 和 Ta 都要求用开尔文中途换成摄氏度会把 a、b 两个经验系数彻底搞乱。实测下来这个环节最常见的错误是把 DN 和辐亮度搞混。ENVI 的辐射定标工具Radiometric Calibration默认输出 BIL 格式的浮点辐亮度影像文件名里带_rad或者带定标类型后缀。做波段运算前用统计功能看一眼数据范围——辐亮度的数量级一般在个位数到几十之间如果看到的是几千上万的整数那它还是 DN 值。2.2 地表比辐射率NDVI阈值法的三种地表类型比辐射率 ε 是单窗算法里最敏感的参数之一。文献里常能看到取 0.95 就行的说法我劝你别这么干。Landsat 8 Band 10 对应 10.6~11.2 微米不同地物的比辐射率在这里能差 0.02 以上折算成温度误差接近 1 K足够让研究结论反过来。我采用的是 NDVI 阈值法把地表分成三类分别赋值。先算 NDVILandsat 8 的红波段是 Band 4、近红外是 Band 5(float(b5) - b4) / (float(b5) b4)这里float()不能省。整数波段相除在 ENVI 里会截断得到的是 0 或 1画面会变成二值图。我第一次做的时候就吃了这个亏盯着满屏黑白图愣了半天。接着算植被覆盖度 Pv((b1 - NDVImin) / (NDVImax - NDVImin))^2NDVImin 和 NDVImax 建议从影像自身的统计值里取一般取累积概率 5% 和 95% 对应的值比较稳比硬套 0.05 和 0.7 更贴合实际。但 Pv 必须限制在 [0,1] 区间ENVI 波段运算没有现成的 clamp 函数可以用关系运算符拼出来(b1 lt 0.05) * 0 (b1 ge 0.05 and b1 le 0.7) * ((b1 - 0.05) / 0.65)^2 (b1 gt 0.7) * 1ENVI 里关系运算返回 1 或 0三个分支相乘再相加等价于分段函数。这个写法比反复裁剪要省事得多值得记下来。有了 Pv比辐射率按地表类型分三种情况地表类型判定条件比辐射率算式水体NDVI 00.995自然表面0 NDVI 0.1570.9625 0.0614·Pv − 0.0461·Pv²城镇建成区NDVI ≥ 0.1570.9589 0.086·Pv − 0.0671·Pv²把三支拼成一个表达式(b1 lt 0) * 0.995 (b1 ge 0 and b1 lt 0.157) * (0.9625 0.0614*b2 - 0.0461*b2*b2) (b1 ge 0.157) * (0.9589 0.086*b2 - 0.0671*b2*b2)其中 b1 是 NDVIb2 是 Pv。水体单独给 0.995 是因为水体在热红外波段的比辐射率接近黑体用自然表面的算式会低估。这套系数来自公开文献中对自然表面和城镇下垫面的拟合结果属于通用取值。如果你的研究区是沙漠、雪盖或者大面积水体务必查一遍该区域专门的比辐射率文献别照搬。2.3 大气透射率水汽含量是唯一的钥匙大气透射率 τ 反映的是热红外辐射穿过整层大气后剩下的比例主要由大气水汽含量决定。获取途径有三条按可靠程度排一下。最稳的是探空资料或者再分析数据直接取整层可降水量再通过查找表插值。这条路精度最高但不是每个项目都有条件。其次是公开的大气校正参数计算服务输入成像日期时间、中心经纬度、气压、气温、相对湿度会自动返回透射率、上行辐射、下行辐射。这个方法门槛低几秒钟出结果我大部分项目都用它。最后是经验公式。文献里针对热红外通道有一些拟合式比如在水汽含量 0.4~1.6 g/cm² 区间用τ 0.974290 − 0.08007w在 1.6~3.0 g/cm² 区间用τ 1.031412 − 0.11536ww 是整层水汽含量。**这些式子是针对特定传感器通道拟合出来的跨传感器套用时心里要有个问号。**如果研究区跨越干湿两季建议用同期探空数据校一遍。如果连水汽含量都拿不到只能退一步用地面气温和相对湿度先算近地面水汽压公式是e0 6.11 * exp(17.27 * T / (237.3 T)) * RH / 100再由水汽压通过经验系数换算整层水汽。这类换算式的系数按月、按区域标定差别不小用之前查一查本区域有没有现成的拟合结果。整层水汽含量 w 估算偏差 0.5 g/cm²对应的透射率偏差大概在 0.04 左右最后反映到地表温度上约 0.4~0.6 K不算致命但也不可忽略。2.4 大气平均作用温度一个气温回归式就够了大气平均作用温度 Ta 表示整层大气按辐射加权后的等效温度没法直接测得靠近地面气温回归。不同大气廓线模型对应的回归系数不一样常用的是这几组大气模型回归式T0 为近地面气温单位 K美国 1976 标准大气Ta 25.9396 0.88045·T0热带大气Ta 17.9769 0.91720·T0中纬度夏季Ta 16.0110 0.92621·T0中纬度冬季Ta 19.2704 0.91118·T0怎么选看研究区的气候带和成像季节。华北平原的七八月选中纬度夏季华南全年偏暖湿用热带大气更合适东北的一月那是妥妥的中纬度冬季。拿不准的时候可以几组都算一遍看差异通常各模型之间 Ta 差 1~3 K反映到最终温度上大约 0.2~0.5 K。这里有个细节值得说T0 必须是成像时刻的近地面气温而且要转成开尔文。不少人图省事取当天最高气温结果正午成像的影像还行上午十点成像的就偏高好几度。如果是多时相分析每个时相都要单独配气温别拿一个值从头用到尾。3. ENVI 5.3 里的完整操作链路3.1 辐射定标只对热红外波段做一次打开 ENVI 5.3走Toolbox → Radiometric Correction → Radiometric Calibration。输入原始影像在定标类型里选 Radiance辐射亮度输出格式选 BIL数据类型选浮点。有个容易忽略的点这个工具会把所有波段一起定标包括可见光近红外。热红外反演其实只需要 Band 10全波段定标会生成一个大文件白白占着硬盘还拖慢后续运算。我的做法是先做波段裁剪把 Band 10 单独提取出来再定标或者定标后立刻用Layer Stacking之外的方式裁掉多余波段。定标完之后用Statistics → Compute Statistics看一眼输出影像的最小值、最大值和均值。正常情况下辐亮度均值应该在 5~15 之间如果均值到了几千那就是在 DN 值上重复乘以了定标系数。关于热红外波段要不要做大气校正不需要。单窗算法已经把大气影响参数化进公式了如果你对 Band 10 再做一次 FLAASH 或 QUAC等于把大气效应剥离了两次结果会系统性偏低。光学波段做大气校正与否对 NDVI 有轻微影响但对 Pv 和 ε 的影响很小用表观反射率算 NDVI 也可以接受。3.2 光学波段把NDVI和Pv一步步算出来NDVI 的表达式前面给过这里补充一个实操细节ENVI 波段运算里输入波段时需要保证所有波段的空间范围和像元大小完全一致。如果红波段和近红外来自同一景影像这自然满足如果是拼接过的数据一定要先做Layer Stacking否则会报维度不匹配。算完 NDVI用统计功能取 5% 和 95% 分位数作为 NDVImin 和 NDVImax。ENVI 的统计面板里有直方图也可以把数据导出后用别的工具算我一般直接看Statistics里的 Percentile 结果。然后是 Pv再然后是按分段函数生成比辐射率 ε。第三步的表达式比较长建议先在文本编辑器里写好、检查括号配对再粘贴进 ENVI 的波段运算输入框。ENVI 的表达式输入框对中文标点极其不友好一个全角括号就能让你找半天错误。生成 ε 之后强烈建议做一次合理性检查用统计功能看最小值最大值应该落在 0.95~0.995 之间。如果出现了 0.98 以下或者大于 0.995 的值说明分段函数的判断条件写错了或者 NDVI 里出现了异常值比如云、云影、雪。3.3 把参数拼成单窗算法表达式单窗算法的核心公式是Ts [a·(1 − C − D) (b·(1 − C − D) C D)·T6 − D·Ta] / C其中C ε · τ D (1 − τ) · (1 (1 − ε) · τ)a −67.355351b 0.458606这两个系数是针对地表温度 0~70°C 区间拟合的国内大部分区域都适用。T6 和 Ta 都用开尔文算出来的 Ts 也是开尔文最后减 273.15 得摄氏度。**别把整个公式塞进一个表达式。**我第一次就这么干一行两百多个字符中间漏了个括号排查了四十分钟。正确的做法是分步生成中间波段第一步算 Cb1 * b2b1 是 εb2 是 τ。第二步算 D(1 - b2) * (1 (1 - b1) * b2)第三步算 Ts此时 b1 Cb2 Db3 T6Kb4 TaK(-67.355351 * (1 - b1 - b2) (0.458606 * (1 - b1 - b2) b1 b2) * b3 - b2 * b4) / b1 - 273.15分步的好处不只是好排查还方便做敏感性分析——想知道 τ 变 0.05 会怎样改一个中间波段重跑第三步就行不用从头再来。生成 Ts 之后立刻看统计值。中国中东部夏季地表温度合理范围大致在 15~55°C西北干旱区夏季白天可以到 60°C 以上。如果你看到 −30°C 或者 200°C 这种值基本可以确定是参数没对齐或者波段顺序搞错了。3.4 出图、裁剪与温度分级结果出来只是半成品。我最常用的几步后处理掩膜异常值Ts 小于 −20°C 或大于 80°C 的像元基本是云、水汽异常或者参数越界导致用Band Math配合关系运算做一个掩膜或者直接Masking → Build Mask。研究区裁剪用矢量边界做Subset Data via ROIs注意保持像元对齐别重采样。温度分级出图城市热岛分析常用等间距分级或者自然断点分级用Color Mapping → Density Slice或者Raster Color Slice分级阈值别硬套按你自己的数据分布定。统计导出需要算区域均温、热岛强度的时候用ROI Tool的统计功能直接读比导出到别的软件再算省事。出图配色上多说一句。热红外结果用彩虹色带是行业惯例但它对色觉障碍不友好而且高值区的细微差别会被压缩。我现在的习惯是用分级色带把关键阈值比如城市热岛强度的分界温度单独设一个醒目颜色读者一眼就能看出重点区在哪。4. 精度到底靠不靠谱怎么验证4.1 三种可行的交叉验证方案没人验证过的 LST 结果说服力是打折的。手头没有实测数据的时候我会用下面几种方式互相印证。第一种与亮温对比。地表温度应该高于或等于亮温因为比辐射率小于 1 会抬高反演温度大气削弱会降低亮温两者差值在植被茂密区通常 1~3 K在裸土区 3~6 K。如果反演结果比亮温还低一大截先把参数回头查一遍。第二种与同区域其他传感器产品对比。MODIS 的 LST 产品时间分辨率高同一天的过境时间差一两个小时可以把两个结果做散点回归看斜率和 R²。要注意空间分辨率差异巨大得先升尺度到同一网格且避开云污染像元。斜率接近 1、R² 超过 0.8就说明量级关系是对的。第三种与气象站地表温度对比。有些气象站有地表温度观测注意是地温不是气温可以取若干站点做验证。这一种最直接但站点的地温观测在空间代表性上跟卫星像元差得远一个站代表不了一平方公里所以别指望单点误差能小于 1 K看整体偏差趋势更有意义。4.2 参数敏感性哪些数错了最致命我在同一景影像上做过参数扰动测试把每个参数人为改一点看 Ts 变化多少结果大致是这样扰动参数扰动幅度Ts 变化敏感程度比辐射率 ε±0.01±0.4~0.6 K高大气透射率 τ±0.05±0.5~0.8 K高大气平均作用温度 Ta±2 K±0.3~0.5 K中近地面气温 T0±2 K±0.2~0.4 K中亮温 T6±0.5 K±0.5 K高结论很清楚ε 和 τ 是命门必须花力气搞准。ε 用固定值 0.95 的做法之所以危险是因为植被茂密区和裸土区的真实 ε 能差 0.02~0.03直接对应 1 K 以上的系统偏差而这个偏差在图上还看不出来——它不表现为斑点而是整片区域一致偏高或偏低非常隐蔽。τ 的敏感度也很高但它的麻烦在于难以获得高精度值。所以如果项目对绝对温度精度要求很高建议优先想办法拿到探空或再分析的水汽廓线别用经验公式硬凑。相比之下 Ta 和 T0 的敏感度温和一些取值的容错空间大用日均值还是瞬时值对最终结论的影响通常在 0.5 K 以内。这不意味着可以随便填而是说在精力有限时优先保证 ε 和 τ 的质量。5. 踩过的坑与速查表5.1 常见问题速查表现象可能原因排查方向结果出现负温度波段顺序错、ε 或 τ 为 0检查中间波段的统计值整片区域温度偏高ε 用了固定 0.95重算比辐射率看均值温度图呈黑白二值波段运算里漏了 float()检查 NDVI 表达式结果与亮温几乎相同τ 被误设为 1查透射率波段是否全为 1边缘出现直线异常影像拼接未对齐重新 Layer Stacking云区温度极低未做云掩膜用光学波段做云检测后掩膜波段运算报错表达式用了全角符号全文替换半角括号逗号多时相结果不可比各期用了不同参数来源统一参数获取方法5.2 几条不写在手册里的心得第一条把参数存成表格别只存在脑子里。每个时相的 ε 均值、τ、Ta、T0、定标系数都记在一个表格里。半年后回来复现或者审稿人问参数怎么取的你能立刻答上来。我吃过这个亏做了一年的多时相序列回头想统一参数口径发现有几期的取值记录早就丢了。第二条先在小区块上试跑。整景影像跑一遍波段运算可能要几分钟到十几分钟参数错了就得重来。用 ROI 裁一个 500×500 的小块先跑通全流程确认温度范围合理之后再上全图。这一步花五分钟能省半小时。第三条比辐射率的分段阈值要看研究区调。0.157 这个分界是文献里针对城镇和自然表面的经验值如果你的研究区是连片农田这个阈值可能把所有耕地都判成城镇比辐射率会系统性偏低。稳妥做法是先在影像上圈几个典型地物样本算一下它们的 NDVI 分布在什么区间再把阈值挪到合适的位置。第四条透射率的经验公式用之前先做个反向验证。如果你手上有几期探空资料用经验公式算一遍再和探空结果比偏差大于 0.05 就说明这套系数不适合你的区域得另找或者自己拟合。这一步很多人省了然后带着一个系统性偏差做完了整个项目。第五条单窗算法对水体的处理要留神。水体比辐射率给 0.995 没问题但水体在热红外波段的热惯性大、日夜温差小如果你的研究重点是水体温度单窗算法的简化假设会带来额外误差这时候考虑用专门的水体温度反演方案更合适。第六条Landsat 8 和 Landsat 9 的定标常数不完全一样。两者 K1、K2 数值接近但不相同做长时间序列的时候每一景都要从自己的 MTL.txt 里读系数别用一套常数通吃。这一点在跨传感器拼接分析里特别容易出错而且往往要等到结果出现系统性偏移才被发现。