☰
国产宽幅多光谱影像黑臭水体与水华反演:ENVI/ArcGIS全流程避坑
2026/10/31 7:12:06 网站建设 项目流程

简介:面向遥感与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 2>&1 && 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 Value=0。

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 = r"D:\water_project\05_mask\城区A_water.shp" raster = r"D:\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)/(b2+b4)提取水体范围大于 0.1城区阴影可能误入
NDVI(b4-b3)/(b4+b3)水华候选大于 0.3岸边植被会抬高
绿红差异指数(b2-b3)/(b2+b3)黑臭候选初筛大于 0.05需按河段标定
蓝红差异指数(b1-b3)/(b1+b3)辅助判断暗色水体结合绿红指数受气溶胶残余影响

分级时不要一上来做全城区统一阈值。某西南山地城区的干流、支流和湖库背景差异很大,上游清水段和下游城区段的绿红差异指数不是一个量级。我的做法是先把水体掩膜内的指数分布拉出来,按河段或行政区做分位数统计,再定“正常水体、疑似黑臭、水华”三类。ArcGIS 里用重分类做分级,ENVI 里也可以用 Band Math 直接输出分类栅格,但分类前一定保留原始指数栅格,后面调阈值还要用。

import arcpy from arcpy.sa import * ndvi = Raster(r"D:\water_project\06_index\NDVI.tif") ndwi = Raster(r"D:\water_project\06_index\NDWI.tif") # 水华候选:NDVI 高且在水体范围内;2 为候选,0 为非候选 algae = Con((ndvi > 0.3) & (ndwi > 0.1), 2, 0) algae.save(r"D:\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 Value=0,再重新做掩膜和指数计算;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 Value=0
水体掩膜河段边界贴合,无大面积锯齿查矢量坐标系和掩膜值
NDWI 统计水体为正,背景不混入查波段映射和分母除零
水华候选斑块集中在湖库和缓流区调整 NDVI 阈值并复核岸边植被
面积统计分类面积与栅格像元数对得上查像元大小和 NoData 设置

6.2 两段脚本:面积统计与成果命名

面积统计我一般用 ArcPy 先把分类栅格转成面积栅格,再按类别汇总,避免在 ArcGIS 属性表里手算。下面这段脚本按像元面积把水华候选转成平方米,实际像元大小从栅格元数据读取,不要写死。

import arcpy from arcpy.sa import * arcpy.env.workspace = r"D:\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(r"D:\water_project\06_index\algae_area.tif") print("pixel_area:", pixel_area)

逻辑说明:meanCellHeight和meanCellWidth从栅格元数据读取像元大小,Con只保留候选类并乘像元面积。参数上,工作空间换成项目06_index目录,cls == 2对应水华候选编码;如果影像坐标系是地理坐标系,面积单位是平方度,出报告前要投影到投影坐标系再算平方米。成果命名统一带日期和阈值版本,比如城区A_水华候选_NDVI030_日期.tif,后面改阈值还能追溯。

从那以后我每次交付前都强制走一遍验收表,先看黑边和 NoData,再看阈值分级,最后才打开 ArcGIS 排版;参数记录不写清楚,过两周自己都说不清这版阈值是怎么来的。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询