Landsat 8遥感影像ENVI全链路处理：从DN值到专题图复现

简介这份ENVI遥感影像处理实验报告面向地信、遥感专业学生及实习实训者围绕Landsat8多光谱影像的完整处理流程展开适合课程实验、软件操作复盘与实习报告撰写参考。内容按实验目的、数据与实验区概况、图像裁剪、辐射定标与反射率求算、大气校正及比较、图像融合、水体提取与植被覆盖度、温度反演、分类与土地利用变化等模块组织融合部分给出Brovey、HSV、主成分变换的目视解译与定量对比温度反演涵盖大气校正法和单窗算法并附结果图与实习心得体会。资源包为1个docx文档约7.92MB图文结合便于按目录查阅操作步骤与成图。目前已有1245人学习可作为遥感影像处理实验的全流程范例帮助读者理解从预处理到专题信息提取、分类及变化分析的组织方式。1. 一份 ENVI 实验报告能复现到什么程度从 DN 值到专题图的全链路带过几届遥感综合实习最常见的场面是实验报告翻了三遍步骤全认识真打开 ENVI 却卡在第一步——波段叠完发现顺序错了定标完发现单位没换算FLAASH 跑完反射率全是 0.0x最后交上来的专题图连自己都不敢信。这份以 Landsat 8 为数据源的 ENVI 影像处理实验报告价值不在于它记录了多少张结果图而在于它把「预处理—大气校正—融合—专题提取—变化监测」串成了一条真能跑通的链路每一步都留了参数、中间产物和对比图。它适合地信、遥感方向需要交一份完整实习报告的人也适合刚接手卫星影像、想找一套可对照流程模板的从业者。下面按我复现时的顺序拆参数怎么填、哪里容易翻车一并写清。2. 预处理三件套波段组合、ROI 裁剪与辐射定标拿到一景 Landsat 8 的 Level-1 产品解压出来是一堆独立 GeoTIFF外加一个 MTL.txt 元数据文件。很多人上来就点 Layer Stacking叠完直接往下做结果后面每一步都在给这一步还债。预处理的顺序和口径一旦定错后面的定标、大气校正、融合全是连锁反应返工成本比一开始多花二十分钟高得多。2.1 先定波段顺序再谈 Layer StackingLandsat 8 OLI 的 Level-1 产品给的是 11 个独立波段文件ENVI 里用 Basic Tools Layer Stacking 把它们叠成一个多波段文件。操作本身不复杂选输入文件、调顺序、设输出路径。翻车点全在「调顺序」上——ENVI 默认按文件名字母序排B1、B10、B11、B2、B3 这么排下来你点完 OK 得到的是一个波段号和物理含义完全错位的文件中心波长、增益偏置、FLAASH 的波段索引全部对不上。先把每个波段的物理含义摆清楚再决定叠哪些波段波长范围μm中心波长μm分辨率主要用途B10.43–0.450.44330 m海岸/气溶胶B20.45–0.510.48230 m蓝水体穿透B30.53–0.590.56130 m绿植被反射峰B40.64–0.670.65530 m红叶绿素吸收B50.85–0.880.86530 m近红外生物量B61.57–1.651.60930 m短波红外湿度B72.11–2.292.20130 m短波红外地质B80.50–0.680.59015 m全色融合用B1010.60–11.1910.895100 m热红外地表温度B1111.50–12.5112.005100 m热红外温度反演我的习惯是分三次叠第一次只叠 B1–B7给大气校正和专题提取用第二次把 B8 单独留着做融合第三次把 B10、B11 单独叠出来做温度反演。原因很直接——FLAASH 只认多光谱反射波段你把全色和热红外一股脑塞进去它要么报错要么按错误的波长解方程出来一片玄学数值。B9 卷云波段一般不进主流程留着做云检测参考就行。手工排序之外我更常用脚本做一遍省得每次开 GUI 数格子import glob import re from osgeo import gdal def band_key(path): # 从文件名里抠出波段号按整数排序避免 B10 排在 B2 前面 m re.search(r_B(\d)\.TIF$, path, re.IGNORECASE) return int(m.group(1)) if m else 999 # 只取 1-7 波段全色和热红外单独处理 files [f for f in glob.glob(LC08_*_B*.TIF) if 1 band_key(f) 7] files.sort(keyband_key) # separateTrue 表示每个输入文件成为一个独立波段而不是拼成一幅大图 vrt gdal.BuildVRT(stack_1_7.vrt, files, separateTrue) gdal.Translate(stack_1_7.tif, vrt)这段代码的关键就两处band_key用正则取波段号转成整数解决字符串排序的坑separateTrue保证每个 GeoTIFF 独立成层而不是空间拼接。输出成 GeoTIFF 之后 ENVI 可以直接打开头文件里的波段数自动识别。如果你习惯全在 ENVI 里做那就手动把 B1 到 B7 依次排好排完在 Layer Stacking 的对话框里再核对一遍顺序别嫌麻烦。2.2 矢量裁剪shp → ROI → Subset 的三步与两个开关研究区一般给的是行政单元矢量第一步不在 ENVI 里做而在 ArcMap 里把目标单元筛出来。打开 shp 的属性表找到行政区名称字段用 Select By Attributes 选中研究区右键图层 → Data → Export Data导出成新的 shp。这一步别偷懒直接拿全省的 shp 往 ENVI 里塞后面 ROI 会大到没法用。回到 ENVI流程是三步File Open Vector File 打开 shp在图层管理器里右键该矢量层 → Export Active Layer to ROIs把它转成 ROI再走 Basic Tools Subset Data via ROIs选影像、选 ROI、设输出。Subset Data via ROIs 关键参数 Select Input File : stack_1_7.tif Select ROI/Vector Layer : 研究区 ROI Mask pixels outside of ROI : Yes Output Background Value : -999或 0看后续处理习惯Mask pixels outside of ROI这个开关是整步的核心。选 No裁出来是一个矩形研究区外面全是背景值选 Yes研究区外的像元被掩掉后续统计和分类不会被背景拉偏。我一般选 Yes同时把背景值设成 -999 而不是 0因为 0 在反射率影像里是有物理含义的比如深水体设成 0 会让后面的阈值分割分不清「背景」和「真水体」这是后面水体提取那一节会反复踩的坑。裁剪本身不改变 DN 值只做空间子集所以像元值口径不变。但有一件事必须在这一步确认矢量和影像的坐标系要一致。不一致就先做 Vector Reproject Vector别指望 ENVI 自动帮你对齐它不会。2.3 辐射定标增益偏置、中心波长与 WGS-84辐射定标是把 DN 值转成有物理单位的辐射亮度公式不复杂Lλ Gain × DN Offset增益和偏置不在影像里在 MTL.txt 里。打开文件搜RADIANCE_MULT_BAND_和RADIANCE_ADD_BAND_每个波段一组。ENVI 扩展工具里的 Apply Gain and Offset 就是干这个的输入影像、按波段填 mult 和 add、输出。也可以用 Radiometric Calibration 面板选 Radiance 输出效果一样。先把系数抠出来核对一遍省得在 GUI 里一个个手输import re mtl open(LC08_L1TP_120038_20200220_20200225_01_T1_MTL.txt).read() for b in range(1, 8): mult float(re.search(rfRADIANCE_MULT_BAND_{b} ([\d.E-]), mtl).group(1)) add float(re.search(rfRADIANCE_ADD_BAND_{b} ([\d.E-]), mtl).group(1)) k1 re.search(rfK1_CONSTANT_BAND_{b} ([\d.E-]), mtl) print(fB{b}: gain{mult:.6f}, offset{add:.6f})RADIANCE_MULT_BAND_x是乘性增益RADIANCE_ADD_BAND_x是加性偏置两个一起决定 DN 到辐射亮度的线性映射。Landsat 8 的偏置通常是个很小的负数看到负号别慌那是正常的。定标输出的单位是 W/(m²·sr·μm)这个单位记牢下一章 FLAASH 会因为它差出十倍。定标完还有两件事。一是改中心波长右键影像 → Edit Header在 Wavelength 里按 B1 到 B7 的顺序填 0.443、0.482、0.561、0.655、0.865、1.609、2.201单位选 Micrometers。这一步错了FLAASH 的波谱响应匹配和后面的温度反演都会偏。二是确认坐标系统一改成 WGS-84跟研究区矢量保持一致。3. FLAASH 大气校正六个参数决定成败大气校正是整条链路上最像「黑匣子」的一步。同样一景数据参数填错一个出来的反射率能差出几倍而且表面上不报错图看起来还挺正常。想让它可控得先知道它在解什么再谈参数怎么填。3.1 大气校正解的是什么方程遥感器收到的辐射亮度不等于地表反射的辐射亮度。中间隔着大气分子的散射、气溶胶的散射、水汽和臭氧的吸收还有地物之间的多次反射。简化后的辐射传输方程长这样L Aρ / (1 - ρeS) Bρe / (1 - ρeS) Laρ 是地表反射率ρe 是像元周围的平均反射率S 是大气球面反照率La 是大气程辐射A、B 是大气透过率和光照几何相关的系数。FLAASH 干的事就是用 MODTRAN 辐射传输模型根据你给的传感器、大气模式、气溶胶、高程、能见度这些参数把这几个系数算出来再从观测辐亮度里反解 ρ。所谓「参数决定成败」本质是这些参数决定了 A、B、La 的取值参数偏一点反解出来的 ρ 就整体平移或缩放。这也解释了一个常见误解大气校正不是调色。它不会让影像「更好看」只会让像元值更接近地表真实反射率。校正完看着比原图灰、对比度低是正常的因为原图那股「鲜亮」里有相当一部分是大气程辐射贡献的假信号。3.2 参数表逐项填传感器、大气模型、气溶胶、高程与波谱响应打开 Spectral FLAASH Atmospheric Correction第一件事是把定标后的文件转成 BIL 格式。FLAASH 对输入的交错方式有要求BSQ 进去可能直接报错Basic Tools Convert Data 里转一下输出 BIL一分钟的事。进主面板之后参数按下面这张表填参数填法依据Input Radiance Image定标后的 BIL 文件单位 W/(m²·sr·μm)Output Reflectance File自定义路径输出浮点型Scene Center Location自动填充从影像头文件读经纬度Sensor TypeLandsat-8 OLI决定波谱响应函数Ground Elevation0.97 km研究区平均海拔约 970 米Flight Date影像获取日期MTL 里 DATE_ACQUIRED不是文件生成日期Atmospheric ModelMid-Latitude Summer按纬度和成像季节选Aerosol ModelUrban研究区周边有城镇和人为气溶胶Water RetrievalYes启用近红外水汽反演Scale Factor10单位换算见下最容易出错的是 Flight Date 和 Scale Factor。Flight Date 要的是影像成像时间不是 MTL 文件的生成时间两者能差好几天大气模式按错日期选水汽含量估计就偏了。Scale Factor 填 10 的原因在单位定标输出是 W/(m²·sr·μm)FLAASH 输入期望 μW/(cm²·sr·nm)两者数值上正好差 10 倍不填这个系数出来的反射率会整体小一个量级看着就像「所有地物都变暗了」。再往下是 Multispectral Settings 对话框要指定水汽反演和气溶胶反演用的波段。报告里的填法是第一个框选第 7 波段第二个框选第 1 波段。这里不同 ENVI 版本的默认值不完全一样别照抄看清楚 GUI 上每个框旁边的提示文字再填填完点 OK 之前把参数页从头到尾过一遍FLAASH 一旦跑起来就是几十分钟参数错了重跑一次很肉疼。3.3 结果验收Z-Profile、光谱曲线与三种典型翻车FLAASH 跑完之后别急着往下做融合和分类先做验收。验收分两层数值分布和光谱形态。数值分布看反射率的取值范围。写几行代码把分位数打出来比肉眼扫图靠谱import numpy as np from osgeo import gdal ds gdal.Open(flaash_reflectance.dat) arr ds.ReadAsArray() # 形状 (波段, 行, 列) # 只统计有效像元排除背景填充值 valid arr[arr -0.1] print(分位数 1/5/50/95/99:, np.percentile(valid, [1, 5, 50, 95, 99])) print(最大值:, valid.max(), 最小值:, valid.min())合理的地表反射率在 0 到 1 之间绝大部分地物落在 0.02 到 0.6。如果 95% 分位数只有 0.05 左右八成是 Scale Factor 没填如果大量像元大于 1先排查云和定标系数再看气溶胶模型是不是选得太离谱。arr -0.1这个阈值是为了排除裁剪时设的背景值背景值设成 -999 的话会被这一步滤掉。光谱形态看典型地物的曲线。在植被区拉一条 Z-Profile正常应该能看到 0.67 μm 附近的红光吸收谷和 0.76 μm 之后的近红外陡坡。校正前这条曲线是平的或者斜率很弱校正后陡坡明显说明大气影响被压下去了。如果曲线整体贴地、陡坡消失回到上一段查 Scale Factor如果曲线在短波红外段剧烈抖动多半是波段顺序或中心波长填错了。三种典型翻车我按出现频率排第一反射率整体偏小十倍Scale Factor 漏填第二影像一半亮一半暗定标时增益偏置按错了波段第三云区反射率溢出成负值或大于 2那是云顶反射太强FLAASH 的反演假设失效这类像元在后续分类里直接掩掉就行不用纠结。4. 融合与专题提取一条链上跑出五类成果预处理和大气校正做完手上是一景 30 米分辨率的反射率影像。接下来的活分两条线一条是把 15 米全色波段融进来提升空间细节另一条是从反射率里抠出水体、植被、温度、分类和变化这几类专题信息。两条线共用同一套像元值口径所以放在一章里说顺序也是我实际操作的顺序。4.1 全色融合四法Brovey、HSV、PC、Gram-Schmidt 的数学路径融合的目的是把全色波段的高空间分辨率和高光谱波段的光谱信息合到一张图上。ENVI 里 Transform Image Sharpening 下面挂了一排方法常用的有四种数学路径完全不同方法核心操作波段限制光谱保真空间增强Brovey各波段乘以 PAN 与波段和的比值通常 3 波段较好中等HSVRGB 转 HSV用 PAN 替换亮度分量再转回必须 3 波段一般较强PC主成分变换用 PAN 替换第一主成分再反变换无硬限制中等中等Gram-Schmidt正交化后用 PAN 替换第一分量无硬限制最好较强Brovey 的公式是逐波段的R R × PAN / (R G B)。它假设全色波段覆盖的波谱范围约等于参与融合的几个多光谱波段之和所以参与融合的波段最好在波谱上连续、范围接近全色。这也是它受波段限制的原因你拿 B5、B4、B3 去做近红外和可见光混在一起光谱畸变就压不住。HSV 是把 RGB 转到色度、饱和度、亮度三个分量用全色替换亮度再转回 RGB。空间细节提升最直观但色调容易偏因为全色的亮度分布和多光谱合成的亮度分布本来就不是一回事。PC 变换先对多光谱做统计正交化第一主成分集中了大部分方差也集中了大部分空间信息。把全色匹配到第一主成分的均值和方差之后替换进去再反变换。它没有波段数限制但第一主成分被换掉之后整体色调变化大报告里那句「PC 融合后反而更模糊、颜色更鲜明本文还有配套的精品资源点击获取