简介这份ENVI遥感影像处理实验报告面向地信、遥感及相关专业学生与实习人员围绕Landsat8多光谱影像的完整处理流程展开可用于课程实验、实习报告撰写与软件操作复盘。压缩包内仅含1个docx文档大小约7.92MB以图文结合方式组织实验目的、数据概况、图像裁剪、辐射定标与反射率求算、大气校正及比较、图像融合等模块并给出Brovey、HSV、主成分变换的目视解译与定量对比结果。报告还涵盖水体提取、植被覆盖度、温度反演大气校正法、单窗算法、监督与非监督分类及土地利用变化分析的步骤图和结果图文末附实习心得体会便于读者对照理解ENVI操作逻辑与成果表达。已有1246人学习适合需要一份较完整遥感实习范例、并希望按目录查漏补缺的初学者和进阶学习者参考。1. 一份ENVI实习报告能省多少摸索Landsat 8从裁剪到变化监测的全链路前阵子有人拿着一份ENVI的Landsat 8影像处理实习报告来问目录看着齐全照着点却在大气校正那步卡了两天出来的反射率曲线还是平的。这份资源的价值不在结果图多漂亮而在于它把一条完整的遥感影像处理链路摊开了——波段组合与矢量裁剪、辐射定标、FLAASH大气校正、四种融合方法横评、水体与植被覆盖度提取、温度反演、监督与非监督分类最后落到两期影像的土地利用变化监测步骤、参数、结果图、定量表格都在。地信和遥感方向的学生可以拿它当操作底稿刚接触影像预处理的从业者可以拿它对照自己漏了哪一步。下面按我复现时的顺序拆重点讲参数怎么定、哪一步容易翻车。2. 波段堆栈与矢量裁剪把8个波段、1个ROI对齐到同一套坐标2.1 波段顺序决定后面所有公式的生死Landsat 8 OLI 的 Level-1 产品里B1 到 B7 是多光谱B8 是全色B9 是卷云B10、B11 是热红外。ENVI 的 Layer Stacking 按你导入的顺序给波段编号Band Math 里的 b4、b5 指的就是堆栈后的第4、第5个波段。一旦顺序乱了NDVI 会算出全负值FLAASH 的光谱匹配也会错位后面所有指数、分类、反演全部作废。这份报告里把 B1 到 B8 正序堆栈把全色波段也放进去理由是融合那一步要用 B8 当高分辨率输入。这个做法没问题但要注意两件事一是堆栈时别把 B9 卷云波段混进去它会打乱编号二是如果确定不做全色融合堆栈时只放 B1 到 B7后面 Edit Header 写波长能少写一组。波段中心波长μm原生分辨率典型用途B1 海岸/气溶胶0.44330 m浅水、气溶胶B2 蓝0.48230 m水体穿透B3 绿0.56130 m植被健康、水体B4 红0.65530 m叶绿素吸收B5 近红外0.86530 m植被、水体边界B6 SWIR11.60930 m湿度、地质B7 SWIR22.20130 m岩性、云检测B8 全色0.59015 m融合、提升空间细节这张表是我自己复现时贴在显示器边上的因为 FLAASH 的 multispectral settings 里要手选波段范围写错一个数字水汽反演就会跑到别的窗口去。2.2 矢量范围先筛属性再导出别整包往里塞报告里的做法是在 ArcMap 里加载行政区 shp打开属性表只保留目标县域那一条记录导出成新的 shapefile。这一步看着琐碎但不做的话ENVI 里导 ROI 会把整张图的所有多边形一起转进去裁剪出来的边界是一堆碎块。属性表筛选有两个细节值得说。第一导出时坐标系选“数据框的坐标系”还是“原数据坐标系”取决于你后面要不要投影。Landsat 8 的 L1TP 产品是 UTM/WGS84 投影而下载到的行政区 shp 常常是地理坐标经纬度。两者不在一个坐标系下ENVI 裁剪时会弹提示选自动投影也能跑但基准面转换一旦选错边界会整体偏移几百米。我一般会在 ArcGIS 里先做一次 Project把矢量投到影像所在的 UTM 带再导入 ENVI这样后面所有步骤都不用再操心对齐问题。第二个细节是属性表里字段的编码。中文县名在某些 shapefile 里会变成乱码筛选时按 FID 或者按面积最大的那条记录选比按名字选稳。2.3 Layer Stacking 与 ROI 裁剪的完整操作序列下面是我按报告步骤整理出来的执行顺序每一步都对应一个明确的输出物中途出问题能立刻定位在哪一环。打开 ENVIFile Open选中 B1 到 B8 共 8 个 TIFF。Basic Tools Layer StackingImport File 依次加入 B1 到 B8确认列表顺序与波段号一致输出文件名带_stack后缀。File Open选 shapefile导入时坐标系选与影像一致的那一套。在 Layer Manager 里右键矢量图层打开属性表确认只剩目标区域一条记录。显示堆栈影像的 RGB常用 5-4-3让影像出现在主视图。File Export Active Layer to ROIs选择导出的记录生成 ROI。Basic Tools Subset Data via ROIs输入影像选堆栈结果Select Input ROIs 选上一步的 ROIMask pixels outside of ROI设为 Yes。输出文件命名带_subset检查统计值范围是否与原图一致。第 7 步的Mask pixels outside of ROI是关键开关。设成 No 时ENVI 只按外接矩形裁边界外的像元照样保留设成 Yes 才是真正按多边形掩膜。掩膜后的背景值默认是 0如果后面要跑 Band Math0 会参与运算并把指数拉偏所以要么在 Band Math 里加判断要么把背景设成 NaN 或 -999。这是我在处理边界水体时踩过的坑后面第 5 章还会展开。2.4 头文件里的波长和单位必须手工补原始 TIFF 被 ENVI 读进来时头文件里没有波长信息。不补的话Spectral Profile 的横轴是波段序号 1、2、3FLAASH 和光谱分析都拿不到物理意义。报告里的做法是右键影像 Edit Header Wavelength把 8 个中心波长按顺序填进去单位从默认的 Nanometers 改成 Micrometers。这里有两个容易漏的点。一是 B8 全色的波长0.590也要填虽然它不参与反射率计算但融合时 ENVI 会读这个值做重采样判断。二是填完之后要打开一次 Spectral Profile 验证横轴显示的是 0.4 到 0.9 之间的波长值才算成功如果还是 1 到 8说明头文件没写进去多半是文件被其他程序占用或者输出时选了只读目录。3. 辐射定标与FLAASH把DN值推到反射率的7个参数3.1 定标公式与MTL里该抄哪几个字段辐射定标做的是把无量纲的 DN 值转成有物理单位的辐亮度再进一步转成大气顶层反射率。Landsat 8 Collection 1 产品在 MTL 头文件里直接给了增益和偏置公式是辐亮度Lλ RADIANCE_MULT_BAND_x × DN RADIANCE_ADD_BAND_x表观反射率ρ (REFLECTANCE_MULT_BAND_x × DN REFLECTANCE_ADD_BAND_x) / sin(SUN_ELEVATION)太阳高度角从 MTL 的SUN_ELEVATION字段取单位是度要先转弧度再取正弦。这一步如果偷懒不除高纬度或者冬季影像的反射率会整体偏高和大气校正后的结果对不上。报告里是在 ENVI 的扩展工具里用 Apply Gain and Offset 手动输入一个一个抄 7 组乘数和加数。手动抄容易串行我一般会先用脚本把 MTL 里的定标参数抽出来对一遍import re from pathlib import Path # MTL 文件名按实际产品替换路径行号用占位符表示 mtl Path(LC08_L1TP_PPPRRR_20200220_20200225_01_T1_MTL.txt) kv {} for line in mtl.read_text(encodingutf-8).splitlines(): m re.match(r\s*([A-Z0-9_])\s*\s*(.?)\s*$, line) if m: kv[m.group(1)] m.group(2).strip() # 逐波段打印定标系数方便和 ENVI 面板里的输入值核对 for b in range(1, 8): print( fB{b}:, Radiance, kv.get(fRADIANCE_MULT_BAND_{b}), kv.get(fRADIANCE_ADD_BAND_{b}), | Reflectance, kv.get(fREFLECTANCE_MULT_BAND_{b}), kv.get(fREFLECTANCE_ADD_BAND_{b}), ) print(太阳高度角:, kv.get(SUN_ELEVATION)) print(成像时间:, kv.get(DATE_ACQUIRED), kv.get(SCENE_CENTER_TIME))脚本只做一件事把每个波段的乘数、加数、太阳高度角和成像时间打出来。参数说明上RADIANCE_MULT_BAND_x是辐亮度乘数量级在 1e-2 左右REFLECTANCE_MULT_BAND_x是反射率乘数量级在 2e-5 左右两者差三个数量级抄错一位结果就完全不对。SUN_ELEVATION用来做太阳高度角校正SCENE_CENTER_TIME是后面 FLAASH 要填的成像时间注意不是FILE_CREATION_TIME。3.2 为什么定标之后必须转成BILFLAASH 是按行扫描读取数据的BSQ 格式波段顺序存放在它眼里会把整幅影像按波段切开读进来的“行”其实是同一个波段的一整块导致校正结果出现条带或者直接报错。报告里专门写了“改变辐射定标文件类型从 BSQ 到 BIL”这一步不是可选项。转换有两种做法一种是在 Basic Tools Convert Data 里选 BSQ to BIL另一种是 Edit Header 里直接改 Interleave 字段。前者会生成新文件后者只改元数据、不动数据本体FLAASH 读的时候还是按原格式解析所以不推荐。转完之后可以用File View Metadata确认 Interleave 一栏显示 BIL再往下走。3.3 FLAASH参数逐项定值错一个整幅图作废FLAASH 的面板参数多但真正影响结果的只有七八个。报告里给出的取值和我复现时用的基本一致整理成表更清楚参数取值依据与说明Input Radiance Image定标后的 BIL 文件必须是辐亮度不能是反射率Sensor TypeLandsat-8 OLI选错会导致光谱响应函数不匹配Sensor Altitude705 km传感器轨道高度默认值可用Pixel Size30 m多光谱原生分辨率热红外另算Ground Elevation0.97 km研究区平均高程从 DEM 统计Flight Date影像获取日期取 DATE_ACQUIREDFlight Time GMT成像时刻取 SCENE_CENTER_TIMEUTCAtmospheric ModelMid-Latitude Summer中纬度夏季按纬度和季节选Aerosol ModelUrban研究区有城镇和建设用地Initial Visibility40 km能见度初值晴朗天气可用Scale Factor10抵消产品自带的 0.1 缩放Scale Factor这个参数最容易填错。Landsat 8 Collection 1 的辐亮度在 MTL 里已经按 0.1 的系数预缩放ENVI 读进来后数值范围是 0 到 1 左右FLAASH 需要的是真实物理量级所以要把输入值乘 10 还原。填小了校正后的反射率整体偏低填大了会溢出。Atmospheric Model的选择依据是纬度和季节不是行政区划。研究区在亚热带夏季影像选 Mid-Latitude Summer 合理如果换成冬季影像还选这个水汽含量会被高估近红外波段反射率偏低。Aerosol Model选 Urban 是因为研究区内有建成区如果整幅图以农田林地为主选 Rural 更合适。3.4 校正结果怎么验别只看图好不好看报告里放了三组对比图原图、定标后、校正后并排还附了 Z-Profile 曲线。这套验证方法是对的但判读标准要明确植被像元红光波段B4反射率应该明显低于近红外B5典型值红光 0.03 到 0.08近红外 0.25 到 0.45曲线呈“红谷近红外峰”。水体像元近红外和短波红外反射率趋近于 0可见光略高整条曲线贴地。裸土像元从可见光到短波红外单调上升没有明显吸收谷。如果校正后植被的红光反而高于近红外先查波段顺序再查定标系数有没有串行。如果整条曲线在蓝波段异常抬高0.2 以上大概率是气溶胶模型选错或者能见度初值给得太小。如果曲线形状对了但整体量级偏低回头看 Scale Factor。4. 四种融合方法横评Brovey、HSV、PC、Gram-Schmidt的取舍与定量打分4.1 融合之前先统一量级全色波段要拉伸多光谱波段经过大气校正后数值范围是 0 到 1 的反射率全色波段 B8 如果只做过辐射定标数值可能是 0 到 1 或者 0 到 10000取决于有没有乘 Scale Factor。两者量级不一致直接融合全色波段会几乎全黑或者全白融合结果的空间细节根本进不来。报告里的做法是对 B8 做反射率拉伸把它拉到和多光谱一致的范围。常见做法有两种用 Band Math 把 B8 乘 10000和多光谱的 0 到 10000 对齐或者用 Stretch Data 把 B8 线性拉伸到 0 到 1。我一般选前者因为整数运算不引入额外的重采样误差。拉伸完用统计面板看一眼均值B8 的均值应该和 B5 在一个数量级差太多说明拉伸系数不对。4.2 Brovey与HSV参数少但光谱代价不一样Brovey 变换的思路是把 RGB 三个波段按比例分配给全色波段的高频信息公式上每个彩色波段乘以“全色 / 三波段之和”。它只支持三个波段输入所以报告里用 B5、B4、B3 合成 RGB全色选 B8重采样选最近邻。最近邻不产生新像元值融合后的光谱畸变最小这也是报告里反复强调“最邻近法效果最清晰”的原因。HSV 变换先把 RGB 转到色度、饱和度、亮度三个通道用全色波段替换亮度通道再转回 RGB。色度和饱和度通道在重采样时也会损失精度所以 HSV 的光谱保真度天然不如 Brovey但纹理提升更明显。报告里的目视结论是 HSV 融合后居住区边界更清晰这和 HSV 对亮度通道的强替换有关。两种方法的共同限制是只能吃三个波段。如果要做多波段融合比如把 7 个多光谱波段全部提升到 15 mBrovey 和 HSV 都用不了只能上 PC 或者 Gram-Schmidt。4.3 主成分变换与Gram-Schmidt多波段融合的正解主成分变换的做法是先对多光谱做 PCA把信息集中到第一主成分然后用全色波段替换第一主成分再做逆变换。替换之前必须把全色波段的均值和方差匹配到第一主成分上否则逆变换出来的影像色调会整体偏。报告里提到 PC 融合“视觉上反而变模糊、颜色更鲜艳”这就是匹配没做好或者第一主成分信息过于集中的典型表现。Gram-Schmidt 是 ENVI 里光谱保真度最好的一档。它的思路是用多光谱波段模拟出一幅低分辨率全色再和真实全色做正交化把差异注入回多光谱。ENVI 5.x 的面板里要选传感器类型Landsat 8 OLI选 Low Res 还是 High Res 取决于你手里的全色波段是不是已经重采样过。报告里没展开这一档但如果要做 7 个波段全部融合GS 是首选。4.4 定量评价四个指标怎么读怎么算报告里裁了一块区域做定量对比算了均值、标准差、信息熵和平均梯度。这四个指标各有指向均值反映辐射亮度整体偏移标准差反映对比度信息熵反映信息量平均梯度反映边缘清晰度。原文的数据整理如下取 B5 红光通道为例影像均值标准差信息熵平均梯度融合前 B510.944.654.1180.158Brovey R0.1990.1873.8880.094HSV R63.1213.855.3390.207PCA R11.433.653.4470.049读这张表要注意两点。第一Brovey 的均值掉到 0.2是因为它的输出被归一化到 0 到 1和融合前的量级不可比只能横向比同量级的方法。第二HSV 的信息熵和平均梯度都最高但均值从 10.94 跳到 63.12说明辐射亮度整体被抬升光谱失真明显。PCA 的信息熵反而低于融合前这和目视看到的“变模糊”一致说明第一主成分替换时信息丢失了。信息熵和平均梯度如果不想手算可以用脚本读 ENVI 导出的 ASCII 或者直接读栅格import numpy as np from osgeo import gdal def entropy(arr, bins256): # 按灰度直方图统计概率分布再算香农熵 hist, _ np.histogram(arr, binsbins, range(arr.min(), arr.max())) p hist / hist.sum() p p[p 0] return float(-(p * np.log2(p)).sum()) def mean_gradient(arr): # 水平和垂直方向一阶差分的平均绝对值 gx np.abs(np.diff(arr, axis1)).mean() gy np.abs(np.diff(arr, axis0)).mean() return float((gx gy) / 2) ds gdal.Open(fusion_hsv_b5.img) band ds.GetRasterBand(1).ReadAsArray().astype(float32) print(信息熵, round(entropy(band), 4), 平均梯度, round(mean_gradient(band), 4))bins256是灰度级数8 位影像用 256 合适如果数据是浮点反射率要先把值域归一化或者加大 bins 数量否则直方图过于稀疏熵值会偏低。mean_gradient用的是简单差分没有做高斯平滑对噪声敏感比较同一景影像的不同融合结果时没问题跨影像比较要谨慎。5. 避坑与排查五个让结果白做的细节5.1 元数据与量级类坑一FLAASH 出来灰蒙蒙植被光谱还是平的。现象是校正后的反射率曲线在蓝波段异常抬高红光和近红外几乎贴在一起目视影像像蒙了一层雾。原因多半是 Flight Time 填成了 MTL 里的FILE_CREATION_TIME而不是SCENE_CENTER_TIME。这两个时间差几个小时到几天太阳天顶角和大气查找表全部对不上水汽和气溶胶反演直接跑偏。解决办法是从 MTL 里取DATE_ACQUIRED加SCENE_CENTER_TIME格式是YYYY-MM-DD和HH:MM:SS注意这个时间已经是 UTC不要再手动加时区偏移。填完重跑一次对比植被光谱的红谷深度正常应该能看出明显吸收。坑二融合时全色波段几乎全黑空间细节进不来。现象是融合结果和原多光谱看不出分辨率差异B8 的统计均值接近 0。原因是量级没统一。大气校正后的多光谱是 0 到 1 的反射率而 B8 如果只做了辐射定标、又乘过 Scale Factor数值范围可能在 0 到 10000。两者相差四个数量级Brovey 的比值计算直接把全色压成 0。解决办法是融合前先看两个输入的统计值把 B8 乘 10000 对齐到多光谱量级或者把多光谱除以 10000。我习惯在文件名里标上量级比如_refl_x10000避免隔一天回来自己都忘了。5.2 空间参照与掩膜类坑三裁剪结果是全黑或者只有角上一小块有数据。现象是 Subset Data via ROIs 跑完之后打开的影像大部分像元值为 0。原因是矢量和影像坐标系不一致。Landsat 8 L1TP 是 UTM 投影下载的行政区 shp 常常是经纬度地理坐标ENVI 自动投影时如果基准面或者 UTM 带号选错多边形会落到影像范围之外。解决办法是先在 ArcGIS 或者 QGIS 里把矢量 Project 到影像的 UTM 带确认两者 Extent 有交集再导入 ENVI。判别方法很直接把矢量和影像同时加载看多边形是不是压在影像上压不上就是坐标系问题不用怀疑裁剪参数。坑四水体提取时边缘河流消失掩膜后和背景混在一起。现象是按报告里的阈值 0.2256 做密度分割主体水体提出来了但研究区边缘那条河整段变黑和掩膜背景分不开。原因有两个。一是阈值照抄。不同时相、不同大气校正结果的 NDWI 分布不一样0.2256 在原始影像上合适换一期影像可能偏高或偏低。二是掩膜背景值设成了 0而暗水体的 NDWI 也接近 0两者在显示上无法区分。解决办法是用 Density Slice 看 NDWI 直方图取双峰之间的谷底作为阈值而不是记一个固定数字。掩膜输出时把背景值设成 -999 或者 NaN后续统计和制图都能把背景排除掉。5.3 评价与精度类坑五分类精度 0.95 以上换一期影像就崩。现象是监督分类的混淆矩阵 Kappa 接近 0.95各类生产者精度都在 0.9 以上但把同一套样本应用到另一期影像分类结果大面积错分。原因是训练样本和验证样本来自同一批 ROI空间上高度重叠空间自相关把精度抬高了。另一种情况是样本选在了两期影像都发生变化的区域模型学到的是变化本身而不是地物特征。解决办法是分层随机抽样训练样本和验证样本按 7:3 分开并且空间上不重叠。看精度不能只看总体精度要拆开看单类的生产者精度和用户精度哪一类低就回去补哪一类的样本。报告里的 ISODATA 非监督分类可以作为对照如果监督分类的类别数和 ISODATA 差太多说明样本的代表性有问题。6. 亮温与单窗算法收尾K值、比辐射率与阈值的三处校验温度反演这一段报告里给了大气校正法和单窗算法两条路。大气校正法本质是用辐射传输方程把大气影响减掉再解算地表温度步骤多、依赖大气参数单窗算法只需要亮温、地表比辐射率和大气透射率三个量落地更快但每个量都得校验。亮温计算用 Planck 反函数ENVI 的 Band Math 里写成K2/alog(K1/L 1) - 273.15。K1、K2 必须和波段对应B10 是 774.8853 和 1321.0789B11 是 480.8883 和 1201.1442。我见过最常见的错是拿 B10 的亮温去配 B11 的 K 值结果温度整体偏高十几度图上看不出来一和历史气象数据对比就露馅。单窗算法的三个输入里地表比辐射率按 NDVI 分段估算植被覆盖高的像元取 0.98 到 0.99裸土取 0.96 到 0.97水体取 0.99大气透射率从 NASA 的水汽产品或者经验公式取研究区夏季一般取 0.6 到 0.75。这两个值一旦拍脑袋定反演结果就是玄学所以我在报告基础上加了一步校验挑几个已知地物深水体、密集植被、裸地取像元值看反演温度是否落在这类地物的合理区间水体夏季表面温度应该低于周边建成区 3 到 5 度。水体提取的阈值也别照抄。NDWI 的密度分割要在直方图上找谷底报告里给的 0.2256 只对那一期影像成立。我现在的习惯是每换一景影像先画直方图、再定阈值、最后用掩膜把背景设成 -999三步缺一不可。从那以后我每次做温度反演都会先把 K 值、比辐射率、透射率三个数写在图旁边跑完对着地物类型抽查一遍才敢出图。希望帮到你。本文还有配套的精品资源点击获取