简介:一份围绕高分一号(GF-1)影像的重庆部分区域水体监测实操方案,面向遥感、GIS学习者及从事环境监测的从业者,重点解决黑臭水体与水华反演问题。资源内含1个docx文档,压缩包约13.95MB,文档详细记录了从数据下载到ArcGIS出图的完整工作流:在陆地观测卫星数据服务网站设置时间、地点、云覆盖率及GF-1传感器参数;ENVI中按步骤完成辐射定标、FLAASH大气校正、RPC正射校正及Mask水域裁剪;通过Band Math输入NDWI/MNDWI公式提取水体信息,最后在ArcGIS中符号化制图并做面积统计。已有557人学习下载。这套资源不仅有每一步的操作说明与参数指导,还附有gf1-wfv3光谱曲线文件选择、Data Ignore Value置零去除黑边等易错点说明,适合需要完整复现水体遥感反演流程的初学者按文档步骤逐步实践。
1. 高分一号黑臭水体反演:不是有影像就能算出水体
我第一次拿到覆盖西南某市的高分一号 WFV 影像时,想法很简单:打开 ENVI,拉伸一下对比度,找个水体指数算一算,黑臭和水华不就浮出来了?结果整幅图花得没法看,云和建筑全被识别成水体,正经河面反而裂成一截一截的。问题不在反演公式,而在原始影像根本不能用——DN 值只是传感器记录的灰度,既没物理意义,也没做过几何纠正,拿它直接算任何比值指数,本质都是在跟噪声做运算。
这个项目就是把「从数据下载到出图」的完整链路拆给你看:陆地观测卫星数据服务页面选片、ENVI 里做辐射定标、FLAASH 大气校正、RPC 正射校正、用目标区域水系 shp 裁出水域,再在 Band Math 里算 NDWI,最后交给 ArcGIS 符号化出图。适合刚接触遥感反演、想跑通一个完整案例的从业者,也适合做水环境监测但不想在预处理上反复折腾的工程师。读完你能说清楚每一步在解决什么问题、参数该往哪个方向调,而不是照着截图瞎点。
2. 数据下载与选片:把检索条件当过滤器而不是填空题
2.1 通用卫星载荷检索的字段怎么填
数据下载这一步,大部分教程一笔带过,实际恰恰是翻车重灾区。登录陆地观测卫星数据服务页面后,进入通用卫星载荷检索,看到的是一排下拉框,很多人习惯性填「全部」「任意」,结果要么列表几千条没法看,要么一张能用的都没有。
我的习惯是把检索条件当过滤器用,而不是填空题。卫星选高分一号,载荷选 WFV。高分一号上有两类相机,PMS 是 2 米全色加 8 米多光谱,适合看细节;WFV 是 16 米多光谱,幅宽大,一景差不多能盖住大半个市辖区,做水体监测完全够用,而且数据量小、预处理快。时间范围不要只给一天,至少留出两周缓冲,高分一号虽然是高重访星座,但你要的是「研究区上空 + 云量低 + 成像质量好」三条件同时满足,窗口期没那么宽裕。地理位置直接填研究区中心经纬度,或者用平台自带的画框工具圈定范围。
云量这一项要特别说。默认是 0 到 100 的区间,新手一看「云量越少越好」,直接填 0。这么做大概率检索结果直接为空,或者只有寥寥几景、成像时间还对不上。云量是成像时真实记录的云覆盖比例,不是按需定制的,西南地区常年多云,想等一景云量 0 的高分一号,运气成分太大。我一般先按 10 到 20 搜,下载前再配合平台提供的快视图确认目标水域附近没云、没有云影。选片阶段的目标是「研究区干净」,不是「整景干净」,云不在水面上就完全可用。
2.2 云量、时相与影像替补:选片翻车现场
时相选择上,黑臭水体反演最好挑夏季晴好天气的影像,气温高、水体流动性差、藻类繁殖快,黑臭和水华的信号最明显。春秋季也能做,但水体温差小、分层弱,反演结果容易边界模糊。另外注意避开刚下过大雨的几天,强降雨会把河道冲一遍,浊度异常升高,这时候检出的「黑臭」可能只是泥沙。
快视图的坑也不小。16 米分辨率的 WFV 在缩略图上看着挺干净,下载回来边缘偶尔会出现拉花或者条带,这种问题在阶段很难筛掉。所以我一般一次下载两到三景同一时间段的数据留替补,尤其是跨越不同轨道号的景。多花的这点下载时间,远比重跑一遍预处理值得。补一句:入库后的影像文件名和产品说明要认真看,记录成像时间、云量、中心经纬度,后面 FLAASH 大气校正要填这些信息,遗漏了还得回去翻元数据。
2.3 下载产物与数据目录规划
下载解压后,一个时相的高分一号 L1 级产品通常会包含影像数据主体、XML 或 TXT 格式的元数据文件,以及 RPC 参数文件。RPC 文件别删,正射校正全靠它。我会在项目根目录下先把目录结构建好,免得预处理到一半文件乱套:
mkdir -p 项目目录/00_参考数据 项目目录/01_原始L1 mkdir -p 项目目录/02_辐射定标 项目目录/03_大气校正 mkdir -p 项目目录/04_正射校正 项目目录/05_水域裁剪 mkdir -p 项目目录/06_反演结果 项目目录/07_出图工程 # 把下载的压缩包解压后统一改名,命名规则:卫星_传感器_成像日期 mv GF1_WFV1_20240815_L1A 项目目录/01_原始L1/GF1_WFV1_20240815这段脚本没什么技术含量,但能帮你省大量时间。命名规则建议固定成「传感器_成像日期」这种一眼可读的格式,FLAASH 和正射校正都要填根文件名,名字越规整越不容易填错。水域 shp、高程 DEM、gf1-wfv3 光谱响应文件这类反复使用的外部数据统一放 00_参考数据,后面所有环节都从这个目录取,不会出现「咦这个文件放哪了」的情况。
3. 预处理三件套:辐射定标、FLAASH 大气校正、正射校正
3.1 辐射定标:先把 DN 值翻译成辐射亮度
高分一号的 L1 级产品,像元存的是 DN 值,也就是传感器记录的数字量化值,它跟地物实际反射的电磁波能量是线性关系,但中间隔着传感器响应、增益、偏置一堆因素。不消除这些影响,绿波段和近红外波段的 DN 值根本没有可比性,NDWI 这种比值指数算出来自然全是噪音。
在 ENVI 里操作不复杂:工具栏搜 rad,打开 Radiometric Calibration,选刚打开的影像。关键参数如下:
| 参数 | 推荐设置 | 说明 |
|---|---|---|
| Calibration Type | Radiance | 输出辐射亮度,单位是 W/(m²·µm·sr) |
| Output Interleave | BIL 或 BSQ | 保持默认即可,后续工具会自动识别 |
| Output Data Type | Float | 定标后数值可能带小数,用整型会截断精度 |
| Scale Factor | 1.0 | 默认值,不要随便填 0.1 或 10,除非你清楚自己在干嘛 |
新版 ENVI 的 Radiometric Calibration 工具会把高分一号的定标参数从元数据里自动读出来,不需要手动输入。点 OK 后会生成一个定标结果,文件名里带 rad 后缀,这个文件就是下一步大气校正的输入。定标完可以顺手打开 Quick Stats 看一眼辐射亮度值的量级,多在 0 到几十之间,DEM 之类的高程数据进 FLAASH 前,最好先把所有波段都确认一遍像元深度和坐标信息,避免参数填错。
3.2 FLAASH 大气校正的逐项参数与光谱响应函数
大气校正是整个预处理里最像黑匣子的一步,参数填错不会立刻报错,而是等到反演阶段给你一份整体偏移的结果。工具栏搜 fla,打开 FLAASH Atmospheric Correction,输入辐射定标后的影像,设置输出文件名和路径,然后按顺序处理参数。
FLAASH 主面板的输入参数每一格都有讲究:
| 参数 | 推荐设置 | 说明 |
|---|---|---|
| Input Radiance Image | 定标结果文件 | 后缀通常是 rad,确认是辐射亮度而不是 DN |
| Output Reflectance File | 自定义命名 | 命名包含 atm 或 ref,区分处理阶段 |
| Sensor Type | Multispectral | 多光谱传感器类型,选错会加载错误模型 |
| Flight Date / Time | 元数据里的成像时间 | 注意 UTC 与北京时间的差异,填错影响大气模型 |
| Scene Center Location | 自动读取后人工复核 | 务必确认中心点落在研究区附近 |
| Sensor Altitude | 保持默认 | GF1 的轨道高度 ENVI 一般能自动带出 |
| Ground Elevation | 高程数据中位数 | 用 SRTM 平均值,别填 0,平原区至少也有几十米 |
| Atmospheric Model | Sub-Arctic Summer / Mid-Latitude Summer | 按季节和纬度选,夏季研究区选亚热带夏季模型 |
| Aerosol Model | Rural | 城区/郊区通用,没有特殊污染源就用它 |
| Water Retrieval | None | 没有 820nm 波段时关闭水汽反演 |
Multispectral Settings 这一项不能跳过。点开后,Filter Function File 一栏选择 gf1-wfv3 光谱曲线文件。这个文件是高分一号 WFV 传感器各波段的响应函数,FLAASH 靠它把影像波段和大气传输模型对应起来,选错了等于传感器模型全部错位。文件可以通过 Open 按钮进入 Spectral Library 选择,我通常直接把工具包附带的 gf1-wfv3 文件拷到 00_参考数据目录,避免每次去翻库。
高级设置全部不动。点 Apply 开始运行,几分钟到十几分钟不等,取决于影像大小和机器性能。跑完输出的大气校正结果,像元值表示地表反射率,范围一般在 0 到 1 之间(有的版本乘了 10000 缩放系数,不影响后续比值计算)。这一步做完,影像才真正有物理意义,之后的比值指数才有对比基础。
3.3 正射校正与 Data Ignore Value:把黑边从统计里踢出去
大气校正后的影像仍然存在地形起伏和传感器姿态造成的几何形变,河流这种细长地物尤其明显,山区的河道可能整体平移几十米。需要做正射校正。工具栏搜 RPC,打开 RPC Orthorectification Using Reference Image,输入大气校正结果,参考影像选同区域的 DEM 或者已校正的影像,高程数据选 SRTM 30 米 DEM,重采样方法用 Bilinear。
正射校正输出的影像四周是黑的,这是无效像元,值为 0。如果直接拿去做反演,黑边会在水体指数图里变成一圈极低值或极高值的异常区,用 ArcGIS 出图时直方图拉伸还会把中间的正常水体压暗。处理办法是给影像设置 Data Ignore Value:在图层列表右键输出结果,选 Edit Metadata,点 Add Metadata Items,找到 Data Ignore Value,把值设为 0。
这一步很多人忽略,但它直接决定后续所有统计结果的可信度。设置完 Data Ignore Value 后,再去做水域裁剪、波段运算,ENVI 会自动忽略值为 0 的像元。需要提醒的是,Data Ignore Value 设置的是「元数据标记」,它告诉软件这些像元无效,但不会物理删除这些像元。因此后续 Band Math 计算时,黑边像元仍然可能参与像素级的运算,唯一的可靠做法是配合第 4 章的掩膜步骤,把无效区彻底排除。
4. 水域裁剪:用掩膜把研究区从影像里抠出来
4.1 Build Mask 前先统一坐标系
拿到正射校正结果后,下一步是只保留水域范围。这里涉及两个常见误区:一是拿 shp 直接用,不管坐标系;二是在 Build Mask 里选错数据源,导致掩膜和影像对不上。
水系 shp 的坐标系经常是通用的地理坐标系或者省级地方坐标系,而正射校正后的影像是投影坐标系(通常是 UTM 或高斯投影)。坐标系不统一,掩膜虽然在视窗里看着叠上了,实际像素级对位会偏几百米。我一般在 ArcGIS 里先查看 shp 的属性表或 .prj 文件,确定原始坐标系,然后用 Project 工具转换到跟影像一致的坐标系。如果 shp 没有 .prj 文件,先用 Define Projection 指定可能的坐标系,再投影转换。转换完成后,在 ENVI 里把影像和 shp 叠加看一眼,河流边界和影像上的水体轮廓贴合才算过关。
确认坐标系一致后,工具栏搜 mask,打开 Build Mask。在 Mask Definition 选项里选择 Vector,导入水域 shp 数据。这里有个细节:Build Mask 输出的是一个掩膜文件,里面只有 0 和 1 两种像元值,作用是把研究区内外区分开,输出路径和文件名自己定。操作时要注意,导入 shp 后,要检查一下掩膜文件覆盖的范围是否跟影像范围一致,有时候 shp 范围大于影像,会导致掩膜文件边缘被截断。
4.2 Apply Mask 参数与结果检查
掩膜文件生成后,工具栏再搜 mask,选 Apply Mask。输入文件选择正射校正结果,掩膜文件选择刚生成的 mask,输出文件名建议命名为 water.dat,方便后面波段运算时识别。Apply Mask 对话框里,默认会列出掩膜值,如果你希望保留的水域区域在掩膜文件里是 1,就确保 Mask Values 只勾选 1,把 0 值区域剔除。
输出后先不要急着做反演,做两件事。第一,用 Quick Stats 查看 water.dat 的像元统计,确认水域内部有有效值,四周被 0 值覆盖。第二,把 water.dat 加载到视窗,叠加水系 shp,检查边缘是否贴合。如果发现水域边缘被削掉一圈,或者掩膜外有残留,说明掩膜值选反了,回 Apply Mask 里改勾选值重跑一次。
需要说明的是,Apply Mask 只是把非目标区域置为无效值,并没有删除像元。这也是为什么前面一定要设 Data Ignore Value,前后配合才能让无效区在统计、制图、波段运算中被一致地忽略。如果偷懒只做 Apply Mask 不设 ignore value,ArcGIS 里拉直方图时会看到大量的 0 值像元参与运算。
4.3 掩膜黑匣子:0 和 1 的意义别靠猜
Build Mask 和 Apply Mask 是 ENVI 里最容易让人栽跟头的两个工具,因为它们对「哪个值有效」的语义在不同版本里可能相反。有的版本掩膜文件里 1 表示有效区域,0 表示无效区域;有的版本在 Apply Mask 里默认 Mask Value 是 0,你反倒要勾选 0 才能保留目标区域。
不要看教程怎么写就照抄,而是做完掩膜后自己验证。方法很简单:加载 mask 文件,用 Cursor Value 工具分别点一下研究区内和区外的像元,看它的值到底是 1 还是 0,再回到 Apply Mask 里勾选对应的值。这个确认动作三十秒,能帮你避免掩膜反了之后整个反演结果报废。掩膜做完,数据准备部分就结束了,接下来进入真正的反演环节。
5. 避坑:高分一号反演最常见的五个翻车现场
5.1 黑边漏进反演结果,水体指数图外围一圈假异常
现象:NDWI 结果图上,影像边缘出现一圈极高值或极低值,正常水体被压缩成一小块。
原因:正射校正后的 0 值黑边没有设置 Data Ignore Value,或者设置了但在 Build Mask 时没有勾选忽略 0 值的选项,导致黑边像元参与了波段运算。0 值像元在 (b2-b4)/(b2+b4) 这个公式里会得到 0/0 的未定义结果,ENVI 通常把它当作 0 处理,于是边缘出现大片「非水」假象。
解决:正射校正后立即通过 Edit Metadata 设置 Data Ignore Value 为 0。Build Mask 时,在 Mask Definition 面板里勾选 Ignore Value 选项,确保掩膜文件同样继承这个无效值标记。最后用 Quick Stats 验证边缘像元不参与统计。
5.2 水体指数算出来全是负值,或者整幅图整体偏移
现象:NDWI 的像元值几乎全在 -0.5 以下,水体和陆地拉不开差距。
原因:最常见的是直接用 L1 级产品的 DN 值算指数,没做辐射定标和大气校正。DN 值里绿波段和近红外波段的量级受传感器增益和大气散射影响,两者不在同一基准上,比值指数自然失真。其次是大气校正时大气模型选错,比如夏季影像选了中纬度冬季模型,导致近红外波段的反射率被系统性压低了。
解决:严格按照第 3 章的流程做完辐射定标和 FLAASH 大气校正,确认大气模型与成像季节、纬度匹配。检查方式是在 ENVI 里对比定标前后同一像元的绿波段数值,正常水体区域的大气校正后反射率应该在 0.03 到 0.08 之间,远低于 DN 值的几百几千。
5.3 FLAASH 报错或输出影像全黑/全白
现象:FLAASH 运行到一半报错,或者跑完输出的反射率影像不是全黑就是全白,直方图集中在一个值上。
原因:多数是 Flight Date / Time 填错,或者 Scene Center Location 自动读取的值偏离研究区太远。时间填错会导致大气模型的水汽含量和气溶胶光学厚度计算完全偏离,输出反射率不是负数就是溢出。
解决:回看下载影像的元数据文件,把成像时间精确到小时分钟,注意 UTC 和北京时间差 8 小时,国产卫星元数据里一般写的是北京时间,FLAASH 里要确认时区设置。Scene Center Location 如果自动读取的经纬度明显偏离,手动改成影像中心点的实际坐标。
5.4 Band Math 公式报错或者结果全为 0
现象:在 Band Math 里输入 (b2-b4)/(b2+b4) 后提示变量未定义,或者运行成功但结果全是 0。
原因:两个细节。第一,不要在 ENVI 里把波段号写成 b1、b2 这样没错,但如果你的影像波段顺序跟假设不同也会错;第二,整型数据直接做除法会截断小数,0 到 255 的整型相减后结果是整型,除法结果也是整型,几乎都被截成 0。
解决:做波段运算前,先用文件头信息确认波段顺序,设置好 Data Ignore Value。公式里用 float() 显式转换数据类型:
# ENVI Band Math 表达式,注意 float 强制转换避免整型截断 ndwi = (float(b2) - float(b4)) / (float(b2) + float(b4))如果公式还报错,检查括号是不是全角字符,从 PDF 或网页复制公式到 ENVI 时经常带入全角括号导致解析失败。
5.5 掩膜后水域像元大面积缺失,只剩零星几个点
现象:Apply Mask 之后 water.dat 里只有少数零散像元,河道主体没了,或者河岸线明显偏移。
原因:水系 shp 和影像坐标系不一致,没有做投影转换就直接裁剪。另一种情况是 shp 范围比影像小很多,或者矢量边界精度不够,16 米分辨率的影像上,河流这种几十米宽的线状地物本来就只有几列像元宽。
解决:裁剪前在 ArcGIS 中检查 shp 的坐标系,统一转换成影像所在投影坐标系。如果河流太窄,只有 1 到 2 个像元宽,反演结果的意义不大,建议换用更高分辨率影像或者对 WFV 影像做相邻像元合并,用「水域占比」而非「单个像元值」来分析。
6. 反演与 ArcGIS 出图:把指数变成别人看得懂的图
6.1 Band Math 里的 NDWI 公式
水域裁剪完成后,反演的核心就是把 NDWI 算出来。NDWI 全称是归一化差异水体指数,公式是 (绿波段 - 近红外波段) / (绿波段 + 近红外波段)。水体在绿波段反射强、近红外波段吸收强,所以正常水体的 NDWI 是正值,陆地植被和建筑多为负值。
高分一号 WFV 的波段顺序通常是 B1 蓝、B2 绿、B3 红、B4 近红外,所以公式写作:
# ENVI Band Math 表达式 ndwi = (float(b2) - float(b4)) / (float(b2) + float(b4))在 Band Math 对话框里输入表达式后,点 Map Variable to Input File,选中水域裁剪结果 water.dat,运行即可。结果文件里,水体部分 NDWI 大于 0,越接近 1 表示水体信号越强。需要提醒的是,MNDWI 需要短波红外波段,而 GF1 WFV 只有 4 个波段,没有短波红外,直接套模板写 MNDWI 公式会算不出来,这是这个传感器最典型的边界。
6.2 ArcGIS 符号化与黑臭水体初判
反演结果要在 ArcGIS 里打开才像样。把 NDWI 结果和 water.dat 一起加载进 ArcMap 或 ArcGIS Pro,先对 NDWI 做符号化,选择拉伸类型为「分类」或「分位数」,把数值范围映射到蓝到绿的色带。正常水体是蓝色,黑臭水体因为浑浊度高、绿光反射弱,NDWI 通常偏低,会显示成深色或暗色调。水华区域因为蓝藻在近红外波段反射较强,在假彩色合成影像上呈亮红色或橙色斑块。
出图时我习惯做一张标准假彩色图叠加 NDWI 半透明图层:假彩色合成用近红外、红、绿波段分别赋红、绿、蓝色通道,水体呈暗色,植被呈红色,黑臭水体呈黑色或深褐色,水华呈亮红色。把这张图加上经纬网、比例尺、图例,导出成 300dpi 的 JPG 或 PDF,反演成果就能拿出去汇报了。
做完这个流程,我给自己立了个规矩:任何反演结果在出图之前,必须把原始影像、掩膜文件、NDWI 结果的 Quick Stats 对一遍,确认数值范围和黑边都干净,再谈图。这个习惯后来帮我挡掉了好几次误判,也希望帮到你。
本文还有配套的精品资源,点击获取