尧图网页设计 
简介基于Python实现的Landsat 8影像地表温度反演算法资源包面向遥感、地理信息系统及环境监测方向的初学者和研究者帮助理解热红外波段反演地表温度的完整流程。压缩包体积仅2KB共包含2个文件主程序脚本负责从影像读取、辐射校正、大气校正、地理校正到亮度温度计算以及单窗法和分裂窗口算法转换地表温度的核心步骤文本辅助文件记录元数据、算法参数或运行日志。已有356人学习/下载。这一实现可实际用于城市热岛效应、植被健康状况或水体温度监测等场景尤其适合希望借助Python生态中的rasterio、numpy等工具快速上手遥感数据处理的学习者。通过拆解这份代码能够掌握热红外遥感反演中的预处理、算法选择、结果验证与比辐射率估算等关键细节是一份轻量但完整的入门参考。1. 拿到Landsat8 zip包真正要解的问题是什么做城市热岛、农业旱情监测、地表能量平衡研究的人几乎都会遇到同一个卡点从USGS或地理空间数据云下载下来的Landsat8 L1级数据是一个zip压缩包解压后是一堆TIF文件——一个波段一个文件外加一个MTL.txt元数据文件。你盯着_B10.TIF这个文件手里拿着python却不知道下一步该按哪个键。这个标题拆开看本质就三件事第一Landsat8的热红外波段Band10/Band11存的是DN值而不是温度需要辐射定标、亮温转换、地表发射率校正才能把“数字”变成“摄氏度”第二反演算法不是只有一种单窗算法、劈窗算法、辐射传输方程法各有适用条件和参数要求选错了结果偏差几度是常事第三这些都是能用pythonnumpyrasterio在本地跑通的计算流程不需要ENVI或者ArcGIS的专用模块。这篇笔记按照“数据准备 → 算法选型 → 代码实现 → 踩坑排查 → 结果验证”的顺序把一个能复现的地表温度反演流程讲透。适合已经会用python读栅格、但没系统做过热红外反演的从业者。2. 从zip到可计算的数组预处理与数据组织2.1 解压、目录规划与MTL元数据的读取Landsat8 L1级数据包解压后典型结构里包含11个波段的TIF文件、一个*_MTL.txt、一个*_ANG.txt几何角度文件以及质量评估文件。地表温度反演真正用到的只有三个波段Band4红、Band5近红外、Band10热红外。Band11在理论上可以做劈窗算法但USGS官方对Band11的定标稳定性有保留意见这里先放一放后面避坑章节会具体说。拿到zip后的第一步不是急着读TIF而是把目录整理成固定结构。我一般会用glob把三景数据放到一个工作目录下然后写一个解析MTL的小函数。MTL文件是文本格式里面以GROUP PRODUCT_METADATA之类的分组形式存放定标系数。用python读取时正则或简单的字符串分割都可以关键是拿到以下字段import re from pathlib import Path def read_mtl(mtl_path): 读取Landsat8 MTL文件提取反演需要的定标参数 params {} with open(mtl_path, r, encodingutf-8) as f: for line in f: # MTL每行格式形如: KEY VALUE 或 KEY 数字 match re.match(r\s*([A-Z_0-9])\s*\s*(.), line.strip()) if match: key, value match.groups() params[key] value.strip().strip() return params mtl_path sorted(Path(data).glob(*_MTL.txt))[0] mtl read_mtl(mtl_path) # 热红外波段的辐射定标系数和热红外常数 ml10 float(mtl[RADIANCE_MULT_BAND_10]) al10 float(mtl[RADIANCE_ADD_BAND_10]) k1_10 float(mtl[K1_CONSTANT_BAND_10]) k2_10 float(mtl[K2_CONSTANT_BAND_10])这段代码并不复杂关键是字段名不带数字后缀容易写成RADIANCE_MULT_BAND_10。我在实际项目中见过有人用python处理时错用了Band11的系数_BAND_11去算Band10的亮温导致最后结果差好几度且毫无规律。MTL里这些系数对同一景影像是固定的但对不同采集时间、不同产品版本会有区别所以每次处理新数据都必须重新读MTL不要硬编码。2.2 地理坐标保持与波段对齐Landsat8的L1级TIF是自带地理参考的GeoTIFFEPSG一般是326xx系列的UTM投影区。用rasterio读取时需要保留profile信息尤其是transform和crs这样最终写出的温度栅格才能和原影像完美叠加。这里有一个新手容易翻车的点用rasterio.open读出来的数组是二维numpy数组但多个波段的坐标系和像元尺寸必须严格对齐否则NDVI和发射率计算时数组错位表面上能算出数值实际空间位置已经对不上了。import rasterio import numpy as np with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B10.TIF) as src: b10 src.read(1).astype(np.float32) profile src.profile.copy() transform src.transform crs src.crs with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B4.TIF) as src: b4 src.read(1).astype(np.float32) with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B5.TIF) as src: b5 src.read(1).astype(np.float32) print(fB10 shape: {b10.shape}, B4 shape: {b4.shape}, B5 shape: {b5.shape}) print(fTransform (B10): {transform})我一般在读取完三个波段后先做一次shape断言不一致就停下来检查不急着往下算。数据预处理到这个程度就够了不需要做大气校正——地表温度反演里的“大气校正”多数情况下是通过大气透过率和大气平均作用温度参数来实现的不是像反射率反演那样先做6S或FLAASH校正。3. 反演算法选型单窗、劈窗还是辐射传输方程3.1 三种主流算法的边界和适用条件地表温度反演算法不是越复杂越好而是要看输入参数能拿到什么。辐射传输方程法大气校正法理论上最严谨公式为LST BT / (1 (λ·BT/ρ)·ln(ε))其中λ是热红外波段的有效波长Band10取10.895μmρ h·c/σ约1.438×10⁻² m·K。但真实应用时它需要大气上行辐射、下行辐射、大气透过率这组参数通常要借助NASA的在线大气参数计算器按影像中心经纬度和成像时间查询。这导致两个问题批量处理多景影像时每一景都要手动查参数自动化程度低而且查询结果针对的是整景影像的平均大气状态局部水汽不均时也会引入误差。劈窗算法利用Band10和Band11两个热红外波段的亮温差异来消除大气影响公式形式有很多种常见的有Price、Ulivieri、覃志豪等版本。它不需要实时大气剖面数据可以说是自动化批处理的理想选择。但问题是Landsat8的Band11定标存在不确定性USGS在Collection 1时代就发布过相关说明提示用户慎重使用Band11做定量反演。这个“理论最优、实践翻车”的属性让劈窗算法在Landsat8上的实际落地率并不高。单窗算法覃志豪2001在实践中最常用。它只需要三个输入热红外波段亮温、地表发射率、大气透过率和大气平均作用温度。其中大气透过率可以根据影像所在纬度带和成像月份按中纬度夏季/冬季的经验值表查取大气平均作用温度可以用近地表气温近似估算。整套流程不依赖外部在线服务python代码完全可以跑通。3.2 单窗算法的参数表和取值逻辑单窗算法的核心公式是LST [a(1 - C - D) (b(1 - C - D) C D)·T - D·Ta] / C其中C ε·τD (1 - τ)·[1 (1 - ε)·τ]。T是亮温KTa是大气平均作用温度Kε是地表发射率τ是大气透过率。a和b是回归系数在0~70℃地表温度区间内通常取a -67.355351b 0.458606。大气平均作用温度Ta的估计有几种经验公式我常用的是中纬度夏季的近似式Ta 17.9769 0.91715·T0其中T0是成像时刻的近地表气温K。这个气温可以从气象站观测数据取也可以从ERA5再分析资料里提取。如果都没有就用“16.0110 0.92621·T0”这个美标大气版本两者在夏季差别很小。大气透过率τ是单窗算法里最影响结果的参数。按覃志豪的查找表中纬度夏季、水汽含量1.0~3.0 g/cm²时Landsat8 Band10对应的τ区间大约在0.7~0.9。我的操作惯例是有探空或再分析水汽数据就按水汽含量插值没有就按夏季取0.75~0.83、冬季取0.85~0.91。宁可在这个参数上多花时间不要用默认值硬算。参数符号取值/来源说明回归系数aa-67.3553510~70℃区间拟合回归系数bb0.4586060~70℃区间拟合大气透过率夏季τ0.75~0.83中纬度水汽1~3g/cm²大气透过率冬季τ0.85~0.91中纬度干燥大气大气平均作用温度Ta17.9769 0.91715·T0T0为近地表气温热红外有效波长λ10.895μmBand10这里需要注意的是a和b的取值文献版本很多有的文章用角度参数或地表温度分段拟合但我复现过多个版本在0~50℃区间差异都在0.3℃以内。真正让结果拉开差距的是τ和Ta尤其是τ。4. 核心实现辐射定标、亮温、NDVI、发射率到单窗反演的完整代码4.1 从DN值到亮温Landsat8的DN值是16位整型辐射定标公式是L ML·DN ALL的单位是W/(m²·sr·μm)。亮温转换公式是T K2 / ln(K1/L 1)T的单位是开尔文。这里有个细节必须说清楚亮温是“黑体等效温度”它没有考虑地表发射率的影响所以直接生成的温度在白天通常比真实地表温度低几度必须走完发射率校正才算完。def dn_to_bt(dn, ml, al, k1, k2): DN值转亮度温度: 先辐射定标再亮温转换 # dn可能包含0值和饱和值先掩膜掉无效像元 valid (dn 0) (dn 65535) radiance ml * dn al bt np.where(valid, k2 / np.log(k1 / radiance 1), np.nan) return bt, valid bt10, valid_mask dn_to_bt(b10, ml10, al10, k1_10, k2_10) print(f亮温范围: {np.nanmin(bt10):.2f} K ~ {np.nanmax(bt10):.2f} K)掩膜这一步不是可有可无。Landsat8影像的边界外、云区域、饱和像元DN值往往为0或65535如果不处理np.log计算时会出inf或nan后面整景结果都会报错。我在项目里用valid_mask来标记有效像元所有中间结果都乘这个掩膜最后写文件时把无效位置设置为NoData。4.2 NDVI与地表发射率估计地表发射率是地表温度反演里最不能拍脑袋的部分。常见做法是Sobrino提出的NDVI阈值法它把地表分成三类水体NDVI 0自然地表NDVI介于0.2~0.5之间全植被NDVI 0.5。对Landsat8的Band10发射率经验公式为ε 0.004·PV 0.986其中PV是植被覆盖度。def calc_emissivity(b4, b5, ndvi_thres(0.2, 0.5)): 基于NDVI阈值法估算地表发射率返回发射率和植被覆盖度 # 计算归一化植被指数 ndvi (b5 - b4) / (b5 b4 1e-10) # 植被覆盖度NDVI线性拉伸到0-1 ndvi_soil, ndvi_veg ndvi_thres pv np.clip((ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil), 0, 1) # 地表发射率取0.986为裸土基值随植被覆盖度线性增加 emissivity 0.004 * pv 0.986 # 水体像元单独处理NDVI0 的像元发射率取 0.995 water ndvi 0 emissivity[water] 0.995 return emissivity, ndvi, pv emissivity, ndvi, pv calc_emissivity(b4, b5) print(f发射率范围: {np.nanmin(emissivity):.4f} ~ {np.nanmax(emissivity):.4f})这段代码里有一个参数值得注意ndvi_thres(0.2, 0.5)。这两个阈值在文献里通常写作NDVI_soil 0.2、NDVI_veg 0.5但对稀疏植被区比如半干旱地区实际影像的NDVI往往整体偏低阈值需要下调到0.15/0.4。我在处理西北地区的Landsat8影像时如果按0.2/0.5来算PV会被大量截断为0导致城市周边地表温度整体偏高。做之前先看NDVI直方图再决定阈值这个步骤省不得。水体像元的发射率取0.995是一种近似。更严格的做法是用MODIS的水体发射率谱库但单窗算法本身对发射率的敏感度大约是“发射率误差0.01 → 温度误差0.5℃左右”所以0.99和0.995的差异可以忽略不需要在这一项上过度纠结。4.3 单窗算法主函数把前面所有中间量汇到一起就是单窗算法的实现。这里我把τ和Ta设为函数参数方便后面批量处理时按景设置不同的值。def single_window_lst(bt_k, emissivity, tau, ta_k): 单窗算法反演地表温度 bt_k: 亮温数组 (K) emissivity: 地表发射率数组 tau: 大气透过率 (标量或数组) ta_k: 大气平均作用温度 (K) a -67.355351 b 0.458606 # C和D是大气影响参数 c emissivity * tau d (1 - tau) * (1 (1 - emissivity) * tau) # 单窗算法主公式 lst (a * (1 - c - d) (b * (1 - c - d) c d) * bt_k - d * ta_k) / c # 输出单位转换为摄氏温度 lst_c lst - 273.15 return lst_c T0 302.15 # 近地表气温 29℃ tau 0.78 # 中纬度夏季水汽约2.0 g/cm² ta_k 17.9769 0.91715 * T0 lst_c single_window_lst(bt10, emissivity, tau, ta_k) lst_c[~valid_mask] np.nan print(f地表温度范围: {np.nanmin(lst_c):.2f} ℃ ~ {np.nanmax(lst_c):.2f} ℃)参数选择上T0我用的是成像时刻的气象站气温。如果手头没有可以用同一天的MOD11 LST产品在影像内平均气温代替误差在可接受范围内。tau 0.78这个值是查表中纬度夏季水汽含量2.0 g/cm²对应的Band10大气透过率如果影像区域当天有降雨或明显锋面过境水汽含量变化会很大此时去NASA查一下再插值才是负责任的。4.4 输出GeoTIFF算完温度只是完成了“反演”把结果写成带地理参考的GeoTIFF才算真正交付。输出时我一般同时写两个文件一个保持开尔文单位的float32原始值一个转成摄氏度后缩放到0~100℃区间的int16乘以100方便在ArcGIS或QGIS里直接做渲染。def write_geotiff(output_path, data, profile, transform, crs, nodata-9999): 写GeoTIFF保留地理参考 profile.update( driverGTiff, heightdata.shape[0], widthdata.shape[1], count1, dtyperasterio.float32, nodatanodata, transformtransform, crscrs ) with rasterio.open(output_path, w, **profile) as dst: dst.write(data.astype(np.float32), 1) # 先把无效值替换为nodata再写文件 lst_out lst_c.copy() lst_out[~valid_mask] -9999 write_geotiff(output/LST_20230601.tif, lst_out, profile, transform, crs)这里有一个容易忽略的操作profile是从Band10的读句柄里拿来的里面可能包含count1、dtypeuint16这些原影像的属性所以我用profile.update把关键字段覆盖掉避免写出一个“类型还是uint16”但数据是float的温度栅格那会直接导致精度丢失。5. 避坑指南5个真实场景里的数据偏差与处理5.1 亮温整体偏高3~5K问题出在没做辐射定标现象用DN值直接代入亮温公式T K2 / ln(K1/DN 1)反演出的地表温度在白天普遍比实测地温高4~8℃。原因Landsat8 L1级产品的DN值和辐射亮度是线性关系但DN值本身量级在几千到三万之间直接代入亮温公式相当于把辐射亮度放大了上百倍ln里面根本不在合理区间。所有反演教程里的K1_CONSTANT_BAND_10应用前提都是“先完成定标到辐射亮度”。解决严格按L ML·DN AL先把DN转辐射亮度再算亮温。判断自己是否踩了这个坑的办法很简单检查亮温数组的值域。Landsat8夏季白天星上亮温一般在290~310K之间如果你的结果出现350K以上的像元且占比不小大概率定标系数没读到或波段用错了。5.2 单窗算法结果比辐射传输方程法低6℃大气透过率选错了现象同一景影像用辐射传输方程法算出地表温度35℃的区域单窗算法只算出29℃。原因在单窗算法中τ取值偏大比如用了冬季的0.9而实际夏季水汽含量高会直接改变C和D的值进而压低LST。τ每偏大0.1LST大约会变化2~3℃这是单窗算法对参数最敏感的环节。解决先用影像中心的经纬度和成像时间在NASA的ATMOSPHERIC PARAMETER CALCULATOR查询反演时刻的大气透过率再把这个τ代回单窗算法两者结果一般能控制在1.5℃以内。如果确实没有网络查询条件就按“夏季保守取0.75冬季取0.88”的惯例并在结果文档中注明参数的不确定性。5.3 zip包解压后Band10和Band11文件混用现象代码报数组维度不匹配或算出结果在空间分布上出现横向条纹。原因批量下载Landsat8数据时一个zip包里同时包含*_B10.TIF和*_B11.TIF而热红外两个波段的像元尺寸虽然都是100m但文件名容易拼错。还有新版Collection 2的产品命名里热红外波段有时写为*_B10.TIF有时又出现在*_SR_B10.TIF表面反射率产品里后者是大气校正后的反射率产品不能直接用于亮温计算。解决每次处理新数据集时先把解压后的文件全部列出来用glob(*B10.TIF)拿到完整匹配再打印文件名确认不是SR_B10。不要手动输入文件名复制粘贴是血泪教训。import glob b10_files glob.glob(data/*B10.TIF) print(匹配到B10文件:, b10_files) # 确认输出中不包含SR这样的表面反射率产品标记5.4 写入GeoTIFF后ArcGIS显示全黑现象用rasterio写出的LST结果文件在ArcGIS或QGIS打开整景黑屏但属性表能查到值。原因地表温度float32数组的值域在290~310K开尔文附近而GIS的默认拉伸范围是0~255所有像元的值都落在255之上因此显示为纯黑。这不是数据出错是渲染拉伸问题。解决写文件时同时输出一个LST_C.tif并将数组乘以100后转为int16值域变为约-3000~4000配合GIS的百分比拉伸就能正常显示。或者直接在python里用matplotlib把结果渲染成PNG标注好色标再交付给别人看。5.5 多景拼接后边缘温度异常现象相邻两景影像分别反演后拼接接缝处出现明显温度台阶有时一侧温度整体高3℃。原因两景影像成像时间不同大气状态、太阳高度角、近地表气温都不同各自代入的τ和Ta也不同。如果用了同一个常数τ去算两景表面看是统一了参数实际等于把东景的大气条件强加到西景上。解决正确的思路是“每景独立估算τ和Ta再拼”。拼接时用gdal.Warp或rasterio.merge并设置nodata-9999让重叠区的边缘像元做加权融合而不是简单覆盖。6. 进阶用MOD11A1交叉验证你的反演结果温度算完不代表工作结束。最怕的是流程全程跑通、数值看着“正常”但放到真实地理场景里对不上。我建议第一次做某个区域的反演时至少用MOD11A1MODIS地表温度产品1km分辨率做一次横向对比。方法是在研究区内随机选取200~500个样本点提取MOD11A1的温度值和你的LST结果做线性回归看相关系数的R²和平均偏差。对比时有三个点要注意。第一MOD11A1是1km分辨率你的Landsat8结果是100m采样时不要单点对比用3x3窗口均值去对应MOD11A1像元。第二MOD11A1白天的过境时间和Landsat8差几小时下午地表温度变化很快允许存在2~4℃的偏差但如果偏离超过6℃就要回头检查发射率和τ。第三云掩膜必须一致MOD11A1自带QC层Landsat8这边你需要用质量评估波段或目视检查剔除云影响。最后批量处理多景时我会在循环里把每隔几景的tau、T0、成像时间打印出来人工扫一遍看看有没有离群值宁可多花一分钟确认参数也不要在没有校验的前提下直接出一批结果。这个习惯让我在好多次数据质量参差不齐的任务里免受返工之苦。地表温度反演这条路本质上就是跟大气和地表发射率两件事较劲。参数调好了精度能贴近实测参数糊弄流程再漂亮也是白搭。希望这篇笔记能帮你少踩几个我已经踩过的坑。本文还有配套的精品资源点击获取 Landsat8地表温度反演：基于Python单窗算法的完整实现
简介基于Python实现的Landsat 8影像地表温度反演算法资源包面向遥感、地理信息系统及环境监测方向的初学者和研究者帮助理解热红外波段反演地表温度的完整流程。压缩包体积仅2KB共包含2个文件主程序脚本负责从影像读取、辐射校正、大气校正、地理校正到亮度温度计算以及单窗法和分裂窗口算法转换地表温度的核心步骤文本辅助文件记录元数据、算法参数或运行日志。已有356人学习/下载。这一实现可实际用于城市热岛效应、植被健康状况或水体温度监测等场景尤其适合希望借助Python生态中的rasterio、numpy等工具快速上手遥感数据处理的学习者。通过拆解这份代码能够掌握热红外遥感反演中的预处理、算法选择、结果验证与比辐射率估算等关键细节是一份轻量但完整的入门参考。1. 拿到Landsat8 zip包真正要解的问题是什么做城市热岛、农业旱情监测、地表能量平衡研究的人几乎都会遇到同一个卡点从USGS或地理空间数据云下载下来的Landsat8 L1级数据是一个zip压缩包解压后是一堆TIF文件——一个波段一个文件外加一个MTL.txt元数据文件。你盯着_B10.TIF这个文件手里拿着python却不知道下一步该按哪个键。这个标题拆开看本质就三件事第一Landsat8的热红外波段Band10/Band11存的是DN值而不是温度需要辐射定标、亮温转换、地表发射率校正才能把“数字”变成“摄氏度”第二反演算法不是只有一种单窗算法、劈窗算法、辐射传输方程法各有适用条件和参数要求选错了结果偏差几度是常事第三这些都是能用pythonnumpyrasterio在本地跑通的计算流程不需要ENVI或者ArcGIS的专用模块。这篇笔记按照“数据准备 → 算法选型 → 代码实现 → 踩坑排查 → 结果验证”的顺序把一个能复现的地表温度反演流程讲透。适合已经会用python读栅格、但没系统做过热红外反演的从业者。2. 从zip到可计算的数组预处理与数据组织2.1 解压、目录规划与MTL元数据的读取Landsat8 L1级数据包解压后典型结构里包含11个波段的TIF文件、一个*_MTL.txt、一个*_ANG.txt几何角度文件以及质量评估文件。地表温度反演真正用到的只有三个波段Band4红、Band5近红外、Band10热红外。Band11在理论上可以做劈窗算法但USGS官方对Band11的定标稳定性有保留意见这里先放一放后面避坑章节会具体说。拿到zip后的第一步不是急着读TIF而是把目录整理成固定结构。我一般会用glob把三景数据放到一个工作目录下然后写一个解析MTL的小函数。MTL文件是文本格式里面以GROUP PRODUCT_METADATA之类的分组形式存放定标系数。用python读取时正则或简单的字符串分割都可以关键是拿到以下字段import re from pathlib import Path def read_mtl(mtl_path): 读取Landsat8 MTL文件提取反演需要的定标参数 params {} with open(mtl_path, r, encodingutf-8) as f: for line in f: # MTL每行格式形如: KEY VALUE 或 KEY 数字 match re.match(r\s*([A-Z_0-9])\s*\s*(.), line.strip()) if match: key, value match.groups() params[key] value.strip().strip() return params mtl_path sorted(Path(data).glob(*_MTL.txt))[0] mtl read_mtl(mtl_path) # 热红外波段的辐射定标系数和热红外常数 ml10 float(mtl[RADIANCE_MULT_BAND_10]) al10 float(mtl[RADIANCE_ADD_BAND_10]) k1_10 float(mtl[K1_CONSTANT_BAND_10]) k2_10 float(mtl[K2_CONSTANT_BAND_10])这段代码并不复杂关键是字段名不带数字后缀容易写成RADIANCE_MULT_BAND_10。我在实际项目中见过有人用python处理时错用了Band11的系数_BAND_11去算Band10的亮温导致最后结果差好几度且毫无规律。MTL里这些系数对同一景影像是固定的但对不同采集时间、不同产品版本会有区别所以每次处理新数据都必须重新读MTL不要硬编码。2.2 地理坐标保持与波段对齐Landsat8的L1级TIF是自带地理参考的GeoTIFFEPSG一般是326xx系列的UTM投影区。用rasterio读取时需要保留profile信息尤其是transform和crs这样最终写出的温度栅格才能和原影像完美叠加。这里有一个新手容易翻车的点用rasterio.open读出来的数组是二维numpy数组但多个波段的坐标系和像元尺寸必须严格对齐否则NDVI和发射率计算时数组错位表面上能算出数值实际空间位置已经对不上了。import rasterio import numpy as np with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B10.TIF) as src: b10 src.read(1).astype(np.float32) profile src.profile.copy() transform src.transform crs src.crs with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B4.TIF) as src: b4 src.read(1).astype(np.float32) with rasterio.open(data/LC08_L1TP_118038_20230601_20230609_02_T1_B5.TIF) as src: b5 src.read(1).astype(np.float32) print(fB10 shape: {b10.shape}, B4 shape: {b4.shape}, B5 shape: {b5.shape}) print(fTransform (B10): {transform})我一般在读取完三个波段后先做一次shape断言不一致就停下来检查不急着往下算。数据预处理到这个程度就够了不需要做大气校正——地表温度反演里的“大气校正”多数情况下是通过大气透过率和大气平均作用温度参数来实现的不是像反射率反演那样先做6S或FLAASH校正。3. 反演算法选型单窗、劈窗还是辐射传输方程3.1 三种主流算法的边界和适用条件地表温度反演算法不是越复杂越好而是要看输入参数能拿到什么。辐射传输方程法大气校正法理论上最严谨公式为LST BT / (1 (λ·BT/ρ)·ln(ε))其中λ是热红外波段的有效波长Band10取10.895μmρ h·c/σ约1.438×10⁻² m·K。但真实应用时它需要大气上行辐射、下行辐射、大气透过率这组参数通常要借助NASA的在线大气参数计算器按影像中心经纬度和成像时间查询。这导致两个问题批量处理多景影像时每一景都要手动查参数自动化程度低而且查询结果针对的是整景影像的平均大气状态局部水汽不均时也会引入误差。劈窗算法利用Band10和Band11两个热红外波段的亮温差异来消除大气影响公式形式有很多种常见的有Price、Ulivieri、覃志豪等版本。它不需要实时大气剖面数据可以说是自动化批处理的理想选择。但问题是Landsat8的Band11定标存在不确定性USGS在Collection 1时代就发布过相关说明提示用户慎重使用Band11做定量反演。这个“理论最优、实践翻车”的属性让劈窗算法在Landsat8上的实际落地率并不高。单窗算法覃志豪2001在实践中最常用。它只需要三个输入热红外波段亮温、地表发射率、大气透过率和大气平均作用温度。其中大气透过率可以根据影像所在纬度带和成像月份按中纬度夏季/冬季的经验值表查取大气平均作用温度可以用近地表气温近似估算。整套流程不依赖外部在线服务python代码完全可以跑通。3.2 单窗算法的参数表和取值逻辑单窗算法的核心公式是LST [a(1 - C - D) (b(1 - C - D) C D)·T - D·Ta] / C其中C ε·τD (1 - τ)·[1 (1 - ε)·τ]。T是亮温KTa是大气平均作用温度Kε是地表发射率τ是大气透过率。a和b是回归系数在0~70℃地表温度区间内通常取a -67.355351b 0.458606。大气平均作用温度Ta的估计有几种经验公式我常用的是中纬度夏季的近似式Ta 17.9769 0.91715·T0其中T0是成像时刻的近地表气温K。这个气温可以从气象站观测数据取也可以从ERA5再分析资料里提取。如果都没有就用“16.0110 0.92621·T0”这个美标大气版本两者在夏季差别很小。大气透过率τ是单窗算法里最影响结果的参数。按覃志豪的查找表中纬度夏季、水汽含量1.0~3.0 g/cm²时Landsat8 Band10对应的τ区间大约在0.7~0.9。我的操作惯例是有探空或再分析水汽数据就按水汽含量插值没有就按夏季取0.75~0.83、冬季取0.85~0.91。宁可在这个参数上多花时间不要用默认值硬算。参数符号取值/来源说明回归系数aa-67.3553510~70℃区间拟合回归系数bb0.4586060~70℃区间拟合大气透过率夏季τ0.75~0.83中纬度水汽1~3g/cm²大气透过率冬季τ0.85~0.91中纬度干燥大气大气平均作用温度Ta17.9769 0.91715·T0T0为近地表气温热红外有效波长λ10.895μmBand10这里需要注意的是a和b的取值文献版本很多有的文章用角度参数或地表温度分段拟合但我复现过多个版本在0~50℃区间差异都在0.3℃以内。真正让结果拉开差距的是τ和Ta尤其是τ。4. 核心实现辐射定标、亮温、NDVI、发射率到单窗反演的完整代码4.1 从DN值到亮温Landsat8的DN值是16位整型辐射定标公式是L ML·DN ALL的单位是W/(m²·sr·μm)。亮温转换公式是T K2 / ln(K1/L 1)T的单位是开尔文。这里有个细节必须说清楚亮温是“黑体等效温度”它没有考虑地表发射率的影响所以直接生成的温度在白天通常比真实地表温度低几度必须走完发射率校正才算完。def dn_to_bt(dn, ml, al, k1, k2): DN值转亮度温度: 先辐射定标再亮温转换 # dn可能包含0值和饱和值先掩膜掉无效像元 valid (dn 0) (dn 65535) radiance ml * dn al bt np.where(valid, k2 / np.log(k1 / radiance 1), np.nan) return bt, valid bt10, valid_mask dn_to_bt(b10, ml10, al10, k1_10, k2_10) print(f亮温范围: {np.nanmin(bt10):.2f} K ~ {np.nanmax(bt10):.2f} K)掩膜这一步不是可有可无。Landsat8影像的边界外、云区域、饱和像元DN值往往为0或65535如果不处理np.log计算时会出inf或nan后面整景结果都会报错。我在项目里用valid_mask来标记有效像元所有中间结果都乘这个掩膜最后写文件时把无效位置设置为NoData。4.2 NDVI与地表发射率估计地表发射率是地表温度反演里最不能拍脑袋的部分。常见做法是Sobrino提出的NDVI阈值法它把地表分成三类水体NDVI 0自然地表NDVI介于0.2~0.5之间全植被NDVI 0.5。对Landsat8的Band10发射率经验公式为ε 0.004·PV 0.986其中PV是植被覆盖度。def calc_emissivity(b4, b5, ndvi_thres(0.2, 0.5)): 基于NDVI阈值法估算地表发射率返回发射率和植被覆盖度 # 计算归一化植被指数 ndvi (b5 - b4) / (b5 b4 1e-10) # 植被覆盖度NDVI线性拉伸到0-1 ndvi_soil, ndvi_veg ndvi_thres pv np.clip((ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil), 0, 1) # 地表发射率取0.986为裸土基值随植被覆盖度线性增加 emissivity 0.004 * pv 0.986 # 水体像元单独处理NDVI0 的像元发射率取 0.995 water ndvi 0 emissivity[water] 0.995 return emissivity, ndvi, pv emissivity, ndvi, pv calc_emissivity(b4, b5) print(f发射率范围: {np.nanmin(emissivity):.4f} ~ {np.nanmax(emissivity):.4f})这段代码里有一个参数值得注意ndvi_thres(0.2, 0.5)。这两个阈值在文献里通常写作NDVI_soil 0.2、NDVI_veg 0.5但对稀疏植被区比如半干旱地区实际影像的NDVI往往整体偏低阈值需要下调到0.15/0.4。我在处理西北地区的Landsat8影像时如果按0.2/0.5来算PV会被大量截断为0导致城市周边地表温度整体偏高。做之前先看NDVI直方图再决定阈值这个步骤省不得。水体像元的发射率取0.995是一种近似。更严格的做法是用MODIS的水体发射率谱库但单窗算法本身对发射率的敏感度大约是“发射率误差0.01 → 温度误差0.5℃左右”所以0.99和0.995的差异可以忽略不需要在这一项上过度纠结。4.3 单窗算法主函数把前面所有中间量汇到一起就是单窗算法的实现。这里我把τ和Ta设为函数参数方便后面批量处理时按景设置不同的值。def single_window_lst(bt_k, emissivity, tau, ta_k): 单窗算法反演地表温度 bt_k: 亮温数组 (K) emissivity: 地表发射率数组 tau: 大气透过率 (标量或数组) ta_k: 大气平均作用温度 (K) a -67.355351 b 0.458606 # C和D是大气影响参数 c emissivity * tau d (1 - tau) * (1 (1 - emissivity) * tau) # 单窗算法主公式 lst (a * (1 - c - d) (b * (1 - c - d) c d) * bt_k - d * ta_k) / c # 输出单位转换为摄氏温度 lst_c lst - 273.15 return lst_c T0 302.15 # 近地表气温 29℃ tau 0.78 # 中纬度夏季水汽约2.0 g/cm² ta_k 17.9769 0.91715 * T0 lst_c single_window_lst(bt10, emissivity, tau, ta_k) lst_c[~valid_mask] np.nan print(f地表温度范围: {np.nanmin(lst_c):.2f} ℃ ~ {np.nanmax(lst_c):.2f} ℃)参数选择上T0我用的是成像时刻的气象站气温。如果手头没有可以用同一天的MOD11 LST产品在影像内平均气温代替误差在可接受范围内。tau 0.78这个值是查表中纬度夏季水汽含量2.0 g/cm²对应的Band10大气透过率如果影像区域当天有降雨或明显锋面过境水汽含量变化会很大此时去NASA查一下再插值才是负责任的。4.4 输出GeoTIFF算完温度只是完成了“反演”把结果写成带地理参考的GeoTIFF才算真正交付。输出时我一般同时写两个文件一个保持开尔文单位的float32原始值一个转成摄氏度后缩放到0~100℃区间的int16乘以100方便在ArcGIS或QGIS里直接做渲染。def write_geotiff(output_path, data, profile, transform, crs, nodata-9999): 写GeoTIFF保留地理参考 profile.update( driverGTiff, heightdata.shape[0], widthdata.shape[1], count1, dtyperasterio.float32, nodatanodata, transformtransform, crscrs ) with rasterio.open(output_path, w, **profile) as dst: dst.write(data.astype(np.float32), 1) # 先把无效值替换为nodata再写文件 lst_out lst_c.copy() lst_out[~valid_mask] -9999 write_geotiff(output/LST_20230601.tif, lst_out, profile, transform, crs)这里有一个容易忽略的操作profile是从Band10的读句柄里拿来的里面可能包含count1、dtypeuint16这些原影像的属性所以我用profile.update把关键字段覆盖掉避免写出一个“类型还是uint16”但数据是float的温度栅格那会直接导致精度丢失。5. 避坑指南5个真实场景里的数据偏差与处理5.1 亮温整体偏高3~5K问题出在没做辐射定标现象用DN值直接代入亮温公式T K2 / ln(K1/DN 1)反演出的地表温度在白天普遍比实测地温高4~8℃。原因Landsat8 L1级产品的DN值和辐射亮度是线性关系但DN值本身量级在几千到三万之间直接代入亮温公式相当于把辐射亮度放大了上百倍ln里面根本不在合理区间。所有反演教程里的K1_CONSTANT_BAND_10应用前提都是“先完成定标到辐射亮度”。解决严格按L ML·DN AL先把DN转辐射亮度再算亮温。判断自己是否踩了这个坑的办法很简单检查亮温数组的值域。Landsat8夏季白天星上亮温一般在290~310K之间如果你的结果出现350K以上的像元且占比不小大概率定标系数没读到或波段用错了。5.2 单窗算法结果比辐射传输方程法低6℃大气透过率选错了现象同一景影像用辐射传输方程法算出地表温度35℃的区域单窗算法只算出29℃。原因在单窗算法中τ取值偏大比如用了冬季的0.9而实际夏季水汽含量高会直接改变C和D的值进而压低LST。τ每偏大0.1LST大约会变化2~3℃这是单窗算法对参数最敏感的环节。解决先用影像中心的经纬度和成像时间在NASA的ATMOSPHERIC PARAMETER CALCULATOR查询反演时刻的大气透过率再把这个τ代回单窗算法两者结果一般能控制在1.5℃以内。如果确实没有网络查询条件就按“夏季保守取0.75冬季取0.88”的惯例并在结果文档中注明参数的不确定性。5.3 zip包解压后Band10和Band11文件混用现象代码报数组维度不匹配或算出结果在空间分布上出现横向条纹。原因批量下载Landsat8数据时一个zip包里同时包含*_B10.TIF和*_B11.TIF而热红外两个波段的像元尺寸虽然都是100m但文件名容易拼错。还有新版Collection 2的产品命名里热红外波段有时写为*_B10.TIF有时又出现在*_SR_B10.TIF表面反射率产品里后者是大气校正后的反射率产品不能直接用于亮温计算。解决每次处理新数据集时先把解压后的文件全部列出来用glob(*B10.TIF)拿到完整匹配再打印文件名确认不是SR_B10。不要手动输入文件名复制粘贴是血泪教训。import glob b10_files glob.glob(data/*B10.TIF) print(匹配到B10文件:, b10_files) # 确认输出中不包含SR这样的表面反射率产品标记5.4 写入GeoTIFF后ArcGIS显示全黑现象用rasterio写出的LST结果文件在ArcGIS或QGIS打开整景黑屏但属性表能查到值。原因地表温度float32数组的值域在290~310K开尔文附近而GIS的默认拉伸范围是0~255所有像元的值都落在255之上因此显示为纯黑。这不是数据出错是渲染拉伸问题。解决写文件时同时输出一个LST_C.tif并将数组乘以100后转为int16值域变为约-3000~4000配合GIS的百分比拉伸就能正常显示。或者直接在python里用matplotlib把结果渲染成PNG标注好色标再交付给别人看。5.5 多景拼接后边缘温度异常现象相邻两景影像分别反演后拼接接缝处出现明显温度台阶有时一侧温度整体高3℃。原因两景影像成像时间不同大气状态、太阳高度角、近地表气温都不同各自代入的τ和Ta也不同。如果用了同一个常数τ去算两景表面看是统一了参数实际等于把东景的大气条件强加到西景上。解决正确的思路是“每景独立估算τ和Ta再拼”。拼接时用gdal.Warp或rasterio.merge并设置nodata-9999让重叠区的边缘像元做加权融合而不是简单覆盖。6. 进阶用MOD11A1交叉验证你的反演结果温度算完不代表工作结束。最怕的是流程全程跑通、数值看着“正常”但放到真实地理场景里对不上。我建议第一次做某个区域的反演时至少用MOD11A1MODIS地表温度产品1km分辨率做一次横向对比。方法是在研究区内随机选取200~500个样本点提取MOD11A1的温度值和你的LST结果做线性回归看相关系数的R²和平均偏差。对比时有三个点要注意。第一MOD11A1是1km分辨率你的Landsat8结果是100m采样时不要单点对比用3x3窗口均值去对应MOD11A1像元。第二MOD11A1白天的过境时间和Landsat8差几小时下午地表温度变化很快允许存在2~4℃的偏差但如果偏离超过6℃就要回头检查发射率和τ。第三云掩膜必须一致MOD11A1自带QC层Landsat8这边你需要用质量评估波段或目视检查剔除云影响。最后批量处理多景时我会在循环里把每隔几景的tau、T0、成像时间打印出来人工扫一遍看看有没有离群值宁可多花一分钟确认参数也不要在没有校验的前提下直接出一批结果。这个习惯让我在好多次数据质量参差不齐的任务里免受返工之苦。地表温度反演这条路本质上就是跟大气和地表发射率两件事较劲。参数调好了精度能贴近实测参数糊弄流程再漂亮也是白搭。希望这篇笔记能帮你少踩几个我已经踩过的坑。本文还有配套的精品资源点击获取