国产宽幅多光谱影像黑臭水体与水华反演：ENVI/ArcGIS全流程避坑

简介面向遥感与GIS初学者及环境监测从业者的实操型文档围绕高分一号影像在重庆部分区域的水体监测展开用ArcGIS与ENVI串起从数据获取、预处理到反演出图的完整链路重点解决黑臭水体与水华识别中流程不清、参数难设的问题。压缩包内为1个docx文件约13.95MB以图文步骤形式记录陆地观测卫星数据服务检索下载、辐射定标、FLAASH大气校正、RPC正射校正、水域掩膜裁剪以及NDWI/MNDWI波段计算与专题制图等环节参数设置与工具入口均有交代。已有560人学习。读者可据此复现一套可落地的高分一号水体信息提取方案理解辐射定标与大气校正的作用、正射校正去黑边的处理细节并掌握从反演结果到ArcGIS符号化出图的衔接方法适合课程作业、项目练手或监测业务入门参考。1. 从一张发暗的河汊影像说起国产宽幅多光谱影像做黑臭水体与水华反演难在哪拿到一景国产宽幅多光谱影像WFV 传感器第一眼看到城区河汊发暗、支流拐弯处泛绿很多人会直接套 NDWI 提水体再拿阈值分黑臭和水华。真正跑一遍才发现卡住交付的往往不是指数公式而是检索时云量没控住、辐射定标缩放系数没对齐、FLAASH 光谱响应文件波段顺序接错、正射黑边没归零最后 ArcGIS 里面积统计对不上。这套流程面向手里有 ENVI 和 ArcGIS、需要把一景宽幅影像从原始 DN 做到专题图的环境监测从业者也适合第一次做水体异常筛查、想把步骤和参数一次钉死的新手。下面按数据获取、预处理、掩膜反演、出图验收的顺序拆开。2. 影像获取与ENVI/ArcGIS工程目录检索条件怎么设才不白跑2.1 为什么选宽幅多光谱而不是高分辨率单景做城区黑臭水体和水华筛查第一步不是比谁的分辨率高而是先看覆盖能力和波段组合。某西南山地城区的河网往往是“干流窄、支流多、湖库散”高分辨率单景幅宽小拼一景城区要十几景甚至更多接边、云影、时相不一致都会把后面的反演搅乱。宽幅多光谱影像单景覆盖大蓝、绿、红、近红外四个波段足够做水体提取、植被和藻类响应分析重访周期也更适合做汛期前后、藻华高发期的对比。黑臭水体识别依赖水体颜色异常、低反射和周边排污口分布宽幅影像的空间细节不如高分辨率但胜在能把河段放在同一个辐射基准下比较。水华识别依赖近红外波段的高反射特征WFV 的四波段组合能做 NDVI 初筛但细小河道会混入岸边植被和建筑阴影混合像元是这套数据的边界。我的做法是先用宽幅影像做全城区筛查把疑似黑臭河段和水华斑块标出来再用高分辨率影像或现场快检做重点复核。ENVI 负责辐射、大气、正射和指数计算ArcGIS 负责成果分级、面积统计和版面出图两个软件的分工在项目一开始就定好后面少很多反复。注意宽幅影像不是万能的城区小于两个像元宽的支流不要硬判黑臭容易把岸边阴影算进去。2.2 检索条件怎么设时间、云量、侧摆与轨道重叠数据获取在某陆地观测数据服务平台完成选用通用卫星载荷检索。时间窗口不要只卡一天我一般按“藻华高发期一景 枯水期一景”来配前者看水华后者看黑臭水体的底色差异。云量控制在 10% 以内城区小河道建议压到 5% 以内云影比云本身更麻烦快视图上看着没云河面上可能一片灰。侧摆角尽量小侧摆大时几何变形和太阳耀斑都会影响水体反射率。传感器选 WFV地点用行政区矢量或流域缓冲面不要只画一个矩形框否则下载回来一半是山体。检索条件建议起点原因常见翻车点时间窗口藻华高发期与枯水期各一景做异常对比不拿单景下结论两景相隔太久河岸施工和水位变化混入云量全图小于 10%河段小于 5%云影会压低水体反射率只看快视图不看云影方向侧摆角尽量小于 15 度减小几何变形和耀斑侧摆大时正射后河宽对不上轨道重叠优先同轨道连续景接边辐射一致不同轨道拼接后色差明显传感器WFV 宽幅多光谱覆盖城区河网误选其他载荷导致波段顺序不同下载前把快视图放大到河段级别看一遍重点看支流汇入口、湖库岸边和城区暗渠出口。文件名里带成像日期、轨道号、传感器和云量别等下载完再靠记忆整理。检索结果多的时候先下元数据和快视图确认河段没被云盖住再下全量能省很多解压和转存时间。2.3 下载后的文件体检与目录规范影像包下载下来先别急着往 ENVI 里拖先做一轮文件体检。压缩包损坏、元数据缺失、解压后少波段这些问题在预处理阶段才暴露会很耗时间。下面这段 bash 放在下载目录里跑逐景检查文件数、大小、元数据和压缩包完整性输出结果直接贴到项目记录里。# 在影像下载目录执行逐景检查文件数、大小、元数据和压缩包完整性 for d in */; do echo $d find $d -maxdepth 1 -type f | wc -l du -sh $d # 检查元数据 xml 是否存在 ls $d/*.xml 2/dev/null | head # 检查 zip 是否完整不是 zip 包会跳过 unzip -t $d*.zip /dev/null 21 echo zip ok || echo zip check skipped done逻辑说明for d in */遍历下载目录下的每个影像包find统计一级文件数量文件数明显偏少说明解压不完整du -sh看包大小是否和下载记录一致ls *.xml确认元数据在不在unzip -t只做完整性测试不改动文件。参数上-maxdepth 1避免把子目录里的临时文件算进去head只显示第一个 xml防止刷屏。目录我习惯按处理阶段分开避免 ENVI 临时文件、ArcGIS 工程文件和原始影像混在一起01_raw放原始包和解压文件02_rad放辐射定标结果03_flaash放大气校正结果04_rpc放正射校正结果05_mask放水域矢量和掩膜文件06_index放 NDWI、NDVI、黑臭候选和水华候选07_map放 ArcGIS 工程和出图08_doc放参数记录和检查表。命名统一成“城区A_WFV_成像日期_传感器_处理级别_步骤.tif”不要用“最终版”“改改改”这种名字。提示路径里不要带中文、空格和特殊符号ENVI 和 ArcGIS 在长路径下报错时排查成本比改目录高得多。3. ENVI预处理四连辐射定标、FLAASH、RPC正射与黑边归零3.1 辐射定标DN值到辐亮度别把缩放系数当装饰原始影像的 DN 值只代表传感器响应不同时间、不同景之间不能直接比。辐射定标把 DN 转成具有物理意义的辐亮度后面 FLAASH 才有输入基础。ENVI 里在工具栏搜索rad打开 Radiometric Calibration选择刚刚打开的宽幅影像。参数不要一路默认尤其是输出类型和缩放系数定标结果值域异常后面全盘皆输。参数项建议设置说明Input File原始多光谱影像确认波段顺序为蓝、绿、红、近红外Calibration TypeRadiance后续大气校正需要辐亮度Output InterleaveBIL兼顾读写速度和兼容性Output Data TypeFloat避免整数截断Scale Factor按元数据给的值常见有 0.1 或 1不能凭感觉填Output Filename城区A_WFV_日期_rad.dat与目录规范一致定标完成后先看统计值不要只看图像亮不亮。辐亮度值全为 0、全为极大值或者空间上出现条带通常是缩放系数和元数据没对上或者输入选成了已经拉伸过的快视图。我的习惯是把定标前后的统计值各记一行到08_doc的参数记录里后面 FLAASH 输出异常时能往前倒查。辐射定标不是“点一下 OK”的步骤它是后面大气校正的黑匣子入口这里错一位后面反射率就敢错一片。3.2 FLAASH大气校正WFV光谱响应文件怎么接大气校正用 FLAASH目的是削弱大气散射和吸收对水体反射率的影响。ENVI 工具栏搜索fla打开 FLAASH Atmospheric Correction输入辐射定标后的影像设置输出文件名和路径。关键在 Multispectral Settings 里的 Filter Function File点开 Filter Function File选择 Open再点 Spectral Library选中 WFV 的光谱响应函数文件。这个文件告诉 FLAASH 每个波段的光谱响应范围选错或者波段顺序对不上输出的反射率会出现整体偏蓝、偏灰或者负值。参数项取值原则常见问题大气模型按成像月份和区域纬度选选错会让蓝波段反射率整体偏移气溶胶模型城区选城市型郊区选乡村型城区气溶胶复杂能见度不要照搬默认初始能见度参考当天能见度记录没有记录时先用默认再做敏感性检查地面高程用 DEM 均值或区域平均高程山地城区高程差大填 0 会失真像元大小按元数据和正射校正的输出像元保持一致光谱响应文件WFV 对应文件波段顺序错会导致反射率曲线变形高级设置一般不动点 Apply 开始输出。FLAASH 跑完后先看反射率统计正常水体在蓝绿波段反射率低近红外更低水华斑块近红外明显抬高。如果整景集中在 0 到 10000 的整数区间常见做法是按 10000 缩放到 0 到 1 再参与指数计算如果出现大片负值先查光谱响应文件波段顺序和辐射定标缩放系数不要急着调大气模型。大气校正后的结果命名成城区A_WFV_日期_flaash.dat和辐射定标结果分目录存放。注意FLAASH 的输出缩放方式不同版本可能不一样指数计算前一定看统计值不要默认它就是 0 到 1。3.3 RPC正射校正与Data Ignore Value去黑边大气校正解决的是光谱问题正射校正解决的是几何问题。ENVI 工具栏搜索RPC打开 RPC Orthorectification Using Reference Image输入大气校正结果、参考影像和 DEM。参考影像选同区域已经正射好的底图DEM 用项目区的高程数据。重采样方式我一般用 Cubic输出像元按参考影像或原始像元大小不要为了“看起来清楚”重采样到很小后面面积统计会虚高。参数项建议设置说明Input File大气校正结果不要用辐射定标结果代替Reference Image同区域正射底图用于几何参考DEM项目区高程数据山地城区必须给ResamplingCubic水体边缘比 Nearest 平滑Output Pixel Size按元数据或参考影像保持面积统计可追溯Output Filename城区A_WFV_日期_rpc.dat与目录规范一致正射输出后别直接掩膜先处理黑边。双击输出结果图层点 Edit Metadata再点 Add Metadata items找到Data Ignore Value把它设为 0。这样旋转和重采样带出来的黑色背景在后续统计里会被识别为无效值而不是被当成低反射水体。黑边不归零是血泪经验里最常见的一条掩膜后水体面积虚高ArcGIS 里一算多了几平方公里回头查半天才发现是背景值混进去了。正射后还要看一眼河宽和桥的位置和参考影像对得上再往下走。如果整景往一个方向偏移查 RPC 参数和参考影像坐标系如果局部扭曲查 DEM 高程异常。正射结果命名成城区A_WFV_日期_rpc.dat黑边归零后的版本可以另存一份参数记录里写清楚Data Ignore Value0。4. 水域掩膜与指数反演NDWI、MNDWI和黑臭/水华阈值怎么落4.1 水域掩膜Build Mask Apply Mask的边界处理正射和黑边处理完下一步把研究区水域裁出来。ENVI 工具栏搜索mask找到 Build Mask选择正射校正结果载入正射结果和城区水域 shp 数据设置好输出路径和文件名生成掩膜文件mask。然后再搜索mask找到 Apply Mask输入文件选正射校正结果掩膜文件选刚生成的mask输出结果命名成water.dat。这两步的顺序不能反Build Mask 是生成掩膜Apply Mask 才是把掩膜套到影像上。参数项建议设置说明Input File正射校正结果确认黑边已归零Mask Vector城区水域 shp坐标系必须和影像一致Mask Value按掩膜规则设 0 或 1输出前预览水体边界Output Filewater.dat只保留水域像元掩膜前先用 ArcGIS 或 ENVI 看一眼 shp 和影像的坐标系是否一致不一致先投影不要直接掩膜。城区水域矢量如果有缓冲区确认缓冲区宽度不会把岸边建筑带进来黑臭水体往往在支流汇入口矢量裁得太干净反而漏掉异常段。掩膜后检查图像四周和河岸边缘出现锯齿状黑边说明掩膜值设反了出现整景全黑说明输入文件和掩膜文件对不上。import arcpy shp rD:\water_project\05_mask\城区A_water.shp raster rD:\water_project\04_rpc\城区A_rpc.tif print(shp:, arcpy.Describe(shp).spatialReference.name) print(raster:, arcpy.Describe(raster).spatialReference.name) # 坐标系名称不一致时先做投影不要直接进 ENVI 掩膜逻辑说明用 ArcPy 读取矢量和栅格的坐标系名称先做一致性检查。参数上shp换成项目里的水域矢量raster换成正射结果如果输出两个坐标系名称不同先在 ArcGIS 里投影矢量或栅格再回 ENVI 做掩膜。这一步花不了一分钟能省掉后面掩膜错位的大麻烦。4.2 Band Math写NDWI/MNDWI公式、波段对应与无效值ENVI 工具栏搜索band打开 Band Math输入水体指数公式。WFV 的波段顺序按蓝、绿、红、近红外排列NDWI 用绿波段和近红外波段(float(b2)-float(b4))/(float(b2)float(b4))。这里的b2是绿b4是近红外具体映射以打开影像后的波段列表为准不要凭记忆填。MNDWI 常见公式需要短波红外波段WFV 四波段里没有短波红外硬套 MNDWI 公式会算错波段这一点在宽幅数据上必须说清楚没有短波红外时用 NDWI 提水体用 NDVI 看水华黑臭用蓝绿红组合做初筛不要为了公式好看硬造一个 MNDWI。# ENVI Band Math 表达式模板b1蓝 b2绿 b3红 b4近红外 ndwi (float(b2)-float(b4))/(float(b2)float(b4)) ndvi (float(b4)-float(b3))/(float(b4)float(b3)) # 分母为 0 时置 0避免出现无效大值 ndwi_safe (b2 eq 0 or b4 eq 0) ? 0 : (float(b2)-float(b4))/(float(b2)float(b4)) ndvi_safe (b4 eq 0 or b3 eq 0) ? 0 : (float(b4)-float(b3))/(float(b4)float(b3)) print(ndwi_safe) print(ndvi_safe)逻辑说明这段 Python 只是把 ENVI Band Math 表达式打印出来复制到 Band Math 输入框即可。参数上b2、b3、b4对应打开影像后的绿、红、近红外波段eq 0判断分母波段是否为 0避免除零。点 OK 后在 Map Variable to Input File 里选中water.dat确认波段映射和表达式一致再输出NDWI.dat和NDVI.dat。指数结果先看统计值NDWI 水体通常为正NDVI 水华斑块明显高于周边水体。提示Band Math 里的波段变量名是 ENVI 按输入文件动态生成的换一景影像就要重新核对映射不要复制上次的表达式直接跑。4.3 黑臭水体与水华反演的指标组合与阈值分级NDWI 只能回答“哪里是水”回答不了“水黑不黑、有没有水华”。黑臭水体反演需要多指标组合先用 NDWI 或掩膜把水体范围定住再用蓝绿红波段的差异做黑臭候选初筛最后用水华指数把藻类高反射区单独分出来。下面这张表是我在项目里常用的初筛阈值组合阈值不是通用常数只是让第一版结果能看最终要拿现场快检、历史黑臭点位或者高分辨率复核来标定。指标表达式或来源用途初筛阈值边界说明NDWI(b2-b4)/(b2b4)提取水体范围大于 0.1城区阴影可能误入NDVI(b4-b3)/(b4b3)水华候选大于 0.3岸边植被会抬高绿红差异指数(b2-b3)/(b2b3)黑臭候选初筛大于 0.05需按河段标定蓝红差异指数(b1-b3)/(b1b3)辅助判断暗色水体结合绿红指数受气溶胶残余影响分级时不要一上来做全城区统一阈值。某西南山地城区的干流、支流和湖库背景差异很大上游清水段和下游城区段的绿红差异指数不是一个量级。我的做法是先把水体掩膜内的指数分布拉出来按河段或行政区做分位数统计再定“正常水体、疑似黑臭、水华”三类。ArcGIS 里用重分类做分级ENVI 里也可以用 Band Math 直接输出分类栅格但分类前一定保留原始指数栅格后面调阈值还要用。import arcpy from arcpy.sa import * ndvi Raster(rD:\water_project\06_index\NDVI.tif) ndwi Raster(rD:\water_project\06_index\NDWI.tif) # 水华候选NDVI 高且在水体范围内2 为候选0 为非候选 algae Con((ndvi 0.3) (ndwi 0.1), 2, 0) algae.save(rD:\water_project\06_index\algae_candidate.tif)逻辑说明用 ArcPy 的 Con 做栅格条件组合先把水华候选从水体范围里分出来。参数上ndvi 0.3和ndwi 0.1是初筛阈值实际项目要按河段调整2是候选类编码0是非候选。输出algae_candidate.tif后进 ArcGIS 符号化不要直接拿它当最终成果后面还要用高分辨率影像或现场记录复核。5. 避坑与排查黑臭水体反演最容易翻车的5个点5.1 大气校正后反射率一片灰或负数现象FLAASH 输出打开后整景发灰水体反射率曲线没有明显的蓝绿红差异统计值里还有大片负值。原因通常有三个辐射定标缩放系数和元数据不一致光谱响应文件波段顺序接错或者气溶胶模型和能见度选得离成像条件太远。解决时先回到辐射定标结果看统计值确认辐亮度没有全 0 或极大值再检查 Multispectral Settings 里的 WFV 光谱响应文件是否按蓝、绿、红、近红外顺序对应最后用当天能见度记录替换默认值跑一小块河段做敏感性对比。不要一看到负值就调大气模型先把输入链路查清楚。5.2 正射黑边没归零掩膜后面积虚高现象正射校正后图像四周有黑色背景Apply Mask 之后水体面积比 ArcGIS 里量的大出一截。原因是黑边像元值被当成低反射水体参与统计掩膜只裁了矢量范围没有把背景值排除。解决是在正射结果元数据里加Data Ignore Value0再重新做掩膜和指数计算ArcGIS 统计时确认栅格 NoData 已设为 0面积统计只算有效像元。这个坑在出图阶段才暴露最气人图和面积表对不上回头重跑一遍预处理。5.3 水域矢量坐标系不一致掩膜错位现象ENVI 里 Build Mask 看起来正常输出水体却整体偏到河岸上或者只掩到半条河。原因是水域 shp 和正射影像坐标系不一致或者 shp 有未知单位。解决是先用 ArcPy 或 ArcGIS 查两者坐标系名称不一致就投影到影像坐标系再检查 shp 的几何是否有自相交、空几何和重复要素。城区水域矢量经常从不同来源拼接坐标系混用很常见掩膜前做一次检查比掩膜后返工便宜得多。5.4 把NDWI当黑臭指数阈值一套用全城现象NDWI 提水体很漂亮但拿同一个阈值分黑臭上游清水段也被判成疑似黑臭下游暗色河段反而漏掉。原因是 NDWI 只区分水和非水不区分水质黑臭需要蓝绿红差异、周边排污口和现场快检共同判断。解决是把水体掩膜和黑臭候选分开做先定水体范围再在水体掩膜内按河段统计绿红差异指数分河段定阈值。全城一套阈值只适合做第一版筛查图不适合直接下结论。5.5 中文路径、长路径、临时文件导致ENVI报错现象FLAASH 跑到一半报错退出或者 RPC 正射输出文件写不进去换一台机器又正常。原因是路径里有中文、空格、特殊符号或者临时目录空间不足、权限不够。解决是把项目根目录改成短英文路径比如D:\water_project所有输入输出都放在这个根目录下ENVI 临时目录留足空间跑 FLAASH 前关掉占用同一目录的同步盘和压缩软件。这个坑不涉及算法但能让前面所有参数白调项目开始就把目录规范钉死。6. 验收与批量复现把单景流程压成一张检查表和两段脚本6.1 出图前验收表单景流程跑完后不要直接打开 ArcGIS 排版。先按下面这张表过一遍任何一项不通过就回到对应步骤别带着问题出图。检查项通过标准不通过处理黑边与 NoData四周背景值为 0 且已设为 NoData回正射结果加Data Ignore Value0水体掩膜河段边界贴合无大面积锯齿查矢量坐标系和掩膜值NDWI 统计水体为正背景不混入查波段映射和分母除零水华候选斑块集中在湖库和缓流区调整 NDVI 阈值并复核岸边植被面积统计分类面积与栅格像元数对得上查像元大小和 NoData 设置6.2 两段脚本面积统计与成果命名面积统计我一般用 ArcPy 先把分类栅格转成面积栅格再按类别汇总避免在 ArcGIS 属性表里手算。下面这段脚本按像元面积把水华候选转成平方米实际像元大小从栅格元数据读取不要写死。import arcpy from arcpy.sa import * arcpy.env.workspace rD:\water_project\06_index cls Raster(algae_candidate.tif) pixel_area arcpy.Describe(algae_candidate.tif).meanCellHeight * \ arcpy.Describe(algae_candidate.tif).meanCellWidth # 只统计候选类 2 的面积单位与影像坐标系一致 area Con(cls 2, pixel_area, 0) area.save(rD:\water_project\06_index\algae_area.tif) print(pixel_area:, pixel_area)逻辑说明meanCellHeight和meanCellWidth从栅格元数据读取像元大小Con只保留候选类并乘像元面积。参数上工作空间换成项目06_index目录cls 2对应水华候选编码如果影像坐标系是地理坐标系面积单位是平方度出报告前要投影到投影坐标系再算平方米。成果命名统一带日期和阈值版本比如城区A_水华候选_NDVI030_日期.tif后面改阈值还能追溯。从那以后我每次交付前都强制走一遍验收表先看黑边和 NoData再看阈值分级最后才打开 ArcGIS 排版参数记录不写清楚过两周自己都说不清这版阈值是怎么来的。希望帮到你。本文还有配套的精品资源点击获取