ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

高分一号水体反演全流程:从辐射定标到黑臭水/水华分级制图

高分一号水体反演全流程:从辐射定标到黑臭水/水华分级制图

简介:面向遥感与GIS学习者的高分一号影像水体反演完整流程资源,以重庆部分区域为案例,重点演示黑臭水体与水华识别方法,全程基于ENVI和ArcGIS两大软件操作。包体内含1个docx文档,压缩包约13.95MB,文档按数据下载、预处理、反演及出图三大阶段组织。数据下载部分介绍陆地观测卫星数据服务网站的参数检索要点;预处理涵盖辐射定标、FLAASH大气校正、RPC正射校正及水体掩膜裁剪;反演阶段给出NDWI与MNDWI波段计算公式,并借助ArcGIS完成符号化制图与空间统计分析。适合需要系统掌握遥感水质监测流程或准备相关课程设计、毕业设计的学生与从业者。该资源已有557人学习下载,文档内附各工具的具体参数设置与操作截图,可帮助读者避开常见坑点,快速复现从原始影像到最终水体分布图的完整链路。

1. 高分一号水体反演全流程:从原始影像到黑臭水与水华分级图

重庆部分区域的黑臭水体和水华反演,不是一个“跑通就行”的活。很多人在 ENVI 里一步一步操作完,到了 ArcGIS 出图阶段才发现:水体指数算出来的结果要么噪点爆炸,要么大片是黑色无效值,要么水域边界把江心洲都吞了。问题十有八九不是反演公式写错,而是前面的数据下载、辐射定标、大气校正、正射校正和掩膜裁剪某一个环节埋了雷。这篇笔记把从数据获取到最终成图的完整流程拆开讲透,覆盖陆地观测卫星数据服务检索、高分一号 WFV 数据下载、ENVI 预处理参数设置、FLAASH 大气校正的细节、RPC 正射校正与去黑边操作,最后落到 ENVI 波段计算器跑 NDWI/MNDWI 以及 ArcGIS 出图分级。适合正在做遥感水质课程设计、本科毕业论文或刚开始接触环境遥感反演的从业者照着复现。

2. 数据下载与产品选型:检索参数、云覆盖阈值与 GF-1 WFV 传感器选择

2.1 为什么选高分一号 WFV 传感器

高分一号搭载了两类载荷,一类是 2 米分辨率全色加 8 米多光谱的 PMS,另一类是 16 米分辨率四波段多光谱的 WFV 宽幅相机。水体反演用的是 WFV,不是 PMS。原因很简单:水体面积大但光谱差异细微,需要的是覆盖范围广、重访周期短、能捕捉水色变化的宽幅数据;16 米分辨率对城市河道和大面积湖库来说已经够用,而且单景覆盖范围可以达到 800 公里幅宽,做区域水体监测一次成像能把大片水域收进来。PMS 的 8 米多光谱虽然在空间细节上更好,但覆盖范围小、单景面积有限,用来做整个区域的连续水体监测需要镶嵌很多景,工作量翻倍而且拼接缝容易干扰水质反演。

GF-1 的 WFV 传感器有 WFV1、WFV2、WFV3、WFV4 四个相机,各相机光谱范围基本一致,都是四个波段:蓝光约 450-520nm、绿光约 520-590nm、红光约 630-690nm、近红外约 770-890nm。这次重庆部分区域的案例用的是 WFV3 相机,光谱响应曲线文件在 ENVI 里可以直接调取,后面的 FLAASH 大气校正需要指定 gf1-wfv3 这个滤波函数文件。验证一下数据头文件里 Sensor Type 是否写对,WFV1 和 WFV3 的光谱响应略有差异,选错光谱函数会让大气校正后的反射率出现系统性偏移,表面看不出问题,NDWI 阈值一变就全乱了。

2.2 陆地观测卫星数据服务检索参数

数据下载步骤在陆地观测卫星数据服务网站完成,选“通用卫星载荷检索”入口。检索时有几个参数直接决定拿到的数据能不能用:

检索参数推荐设置说明
时间范围反演目标时段前后 15-30 天时间太长水体状态变化大,时间太短可能一景都没有
空间范围重庆目标区域行政区或经纬度框用经纬度矩形框更精确,避免下载无关区域浪费硬盘
云覆盖量≤ 10%,最多放宽到 20%重庆多雾,20% 以上的云量基本没法做大气校正
传感器GF-1 WFV(注明 WFV1/WFV3 等相机号)不同相机光谱响应不同,FLAASH 参数设置时要用对应文件
产品级别L1A 级预处理产品这是几何粗校正产品,适合接下来自己做辐射定标和大气校正

检索出来的结果里,同一个时间段往往有多景候选影像,不要只看文件名,点开查看每景的云覆盖百分比和缩略图。云覆盖量是整个文件的统计值,容易出现“整体云量 5%,目标河道正好被一片云盖住”的情况,这时候老老实实换一景重复操作,不然大气校正时被云覆盖的水体反射率会被强行插值,反演结果会出现一条条“云影假水体”。

2.3 下载后的第一步:先做数据体检

下载下来的高分一号数据一般是压缩包,解压后里面包含一个 XML 元数据文件、一个 TIFF 或 GeoTIFF 影像文件以及 RPC 系数文件。我一般会先用 ENVI 打开影像检查三件事:影像波段数量是否为 4 个、空间分辨率是否为 16 米、影像最上方是否有明显条带或坏线。检查元数据里的成像时间也很重要,FLAASH 大气校正需要精确到小时的成像时刻,如果元数据时间和下载列表里的时间对不上,输出结果在校正时用的太阳高度角会不准,水体反射率会整体偏高或偏低。这个体检过程只有 5 分钟,但能避开后面两小时的白费功夫。

3. 辐射定标与 FLAASH 大气校正:参数表、光谱响应函数与常见错误

3.1 辐射定标:把 DN 值变成亮度值

高分一号下载的 L1A 级产品存的是原始 DN 值,也就是传感器记录的数字量化值。DN 值本身没有物理意义,不同时间、不同光照条件下同一个地物的 DN 值可能差很多,所以必须先做辐射定标,把 DN 值转换成辐亮度或表观反射率。在 ENVI 的工具栏搜索 rad,打开 Radiometric Calibration 工具,选择打开的高分影像。这里有几个关键参数:

定标参数推荐设置原因
Calibration TypeRadiance输出辐亮度,供 FLAASH 使用
Output Data TypeFloat浮点型,避免整型截断导致波段差值失真
Scale Factor1.0保持原始量级,FLAASH 能正常读取
Output StoredInterleaveBIL 或 BSQ 均可以,FLAASH 推荐 BSQ

点击 OK 后输出辐亮度文件,文件名我用 rad.dat 这种明显标识,方便后面找文件。定标结果可以快速拉一个波段查看亮度和对比度,重庆地区植被覆盖率高,正常影像应显示绿色植被区域明显偏亮、水体区域较暗,如果整体亮度异常那大概率是数据类型设置不对。

3.2 FLAASH 大气校正参数设置

辐射定标完成后,在 ENVI 工具栏搜索 fla,点击 FLAASH Atmospheric Correction。这里的大气校正是整个流程里参数最多、最容易出问题的一步,也是决定反演精度上限的一步。

首先是 Radiance Image 选择刚才的辐射定标结果,Output Reflectance File 设置输出反射率文件名,这里我习惯命名为 ref.dat。接着需要填写中心点经纬度、传感器类型、飞行时间和大气模型。中心点的经纬度直接读 FLAASH 界面提供的 RFL 信息,或者在影像几何信息里查看,注意坐标要填对,填到别的城市去了,大气模型和太阳几何角度全部错位。

大气模型的选择以影像获取季节和纬度带为基准,重庆地处中低纬度,气候湿润,夏季选 Tropical,冬季选 Mid-Latitude Summer 或 Mid-Latitude Winter。重庆的最大问题是水汽含量高,大气水汽反演经常报错,所以在 Multispectral Settings 里打开 Filter Function File 时,一定要确认选择的 gf1-wfv3 光谱响应函数和图像传感器一致,否则水汽反演这一列数据会出现大量 0 值。

气溶胶模型,城区及工业区选 Urban,郊区植被密布选 Rural。重庆主城区水体周边建筑密集,气溶胶成分复杂,选 Urban 更稳。初始能见度 Initial Visibility 通常设置为 20-40 km,重庆多雾天气,能见度很可能只有 10-20 km,如果设置太高,大气校正输出的反射率在蓝绿波段容易出现负值。

3.3 光谱响应函数加载方法

在 FLAASH 界面点击 Multispectral Settings,打开后按 Open,然后进入 Spectral Library 找 gf1-wfv3 光谱响应曲线文件。这个文件是 ENVI 内置的高分一号 WFV3 传感器光谱响应函数,误差不要自己拍脑袋去修改。选择完成后 FLAASH 会读取每个波段的响应函数来计算大气透过率和散射贡献。这里要清楚一点:如果找不到 gf1-wfv3,直接选了 gf1-wfv1,大气校正虽然能跑完,但输出的反射率在绿光波段和近红外波段会有一个系统性偏差,NDWI 的值整体被抬高或压低,后面阈值分割时用固定值就不准了。

设置好后点击 Apply,FLAASH 开始大气校正过程。处理时间取决于影像大小和机器性能,16 米分辨率的 WFV 数据一般需要 15-30 分钟。如果中途报错或者输出的反射率影像出现大面积负值,优先检查两个地方:一是输入影像是否真的是辐亮度而非反射率,二是能见度和气溶胶模型是否过于激进。负值偶尔出现在个别像元是正常现象,出现在水体区域则说明初始能见度设置偏小或气溶胶模型偏滞后。

3.4 大气校正结果的验证方法

大气校正做没做准,不用等到最后反演阶段。快速验证方法是对比水体在近红外波段的反射率,深水体的近红外反射率通常很低,一般在 3% 到 5% 以下。如果近红外反射率在 8% 以上,要么是水体本身含沙量很高,要么是大气校正没做干净。另一个方法是看水体在蓝光波段的反射率是否比绿光低,正常清水体符合这个规律,如果蓝光高于绿光且数值异常偏大,说明大气校正的瑞利散射部分没有消除干净。

4. 正射校正、去黑边与水域裁剪:RPC 模型参数与掩膜处理细节

4.1 RPC 正射校正参数选择

大气校正后的反射率影像还存在地形起伏引起的几何形变。重庆是典型山地城市,地形高差大,河道两侧有大量梯坎和山体,如果不做正射校正,山体阴影和水体会互相错位,在 ArcGIS 里叠加水域边界时会出现边界错开的地方。ENVI 的正射校正工具是 RPC Orthorectification Using Reference Image,在工具栏搜索 RPC 找到它。

参数设置上,输入文件是大气校正后的反射率影像,参考影像可以用已有的经过精确校正的高分辨率影像或具有准确地理定位的底图,高程数据用覆盖重庆区域的 DEM,我用的是 ASTER GDEM 30 米分辨率高程数据,对 16 米分辨率的 WFV 数据足够。重采样方法选三次卷积,输出像元大小设置为 16 米,保持原始分辨率。

需要注意 RPC 正射校正时需要读取影像自带 RPC 系数文件,下载的数据包里的 RPC 文件是关键,没有它正射校正是跑不动的。如果 ENVI 提示找不到 RPC,检查影像文件是否和 .rpc 同名且放同一目录,某些下载数据解压后会丢失 RPC 文件导致无法处理,这时候只能重新下载。

4.2 去除黑色背景值:Data Ignore Value 设置

正射校正输出的结果影像边缘通常有一圈黑色背景,这些像元值为 0,但和真实水体里的 0 反射率完全不同。如果不处理,后续波段计算和水域裁剪会把黑色背景当成“超低反射率纯黑水体”,NDWI 计算时这些区域的值会被拉低,水域提取结果会出现从影像边缘延伸进内部的假水道。

去除黑边要双击输出结果图层,选择 Edit Metadata,然后点击 Add Metadata items,找到 Data Ignore Value 并添加,将值设置为 0,点击 OK 确认。这一步是把黑色背景像元标记为忽略值,后续统计和波段运算时 ENVI 自动跳过这些像元。做完这步再打开图层,黑边部分不会参与后面的数值计算。

4.3 Build Mask + Apply Mask 裁剪水体

反演黑臭水体和水华,只关心水体区域,陆地和植被的信号会干扰阈值提取。水域裁剪分两步:第一步用 Build Mask 构建掩膜,第二步用 Apply Mask 应用掩膜。

Build Mask 里输入文件选正射校正后的影像,然后载入重庆水域 shp 矢量数据。这里有个容易踩的坑:Build Mask 默认认为加载进来的矢量区域是“感兴趣区”,但需要确认掩膜输出的逻辑值方向。在面板里把 Shp 区域的值设为 1,其他区域设为 0,保存成 mask 文件。如果你的 shp 是线状水系或者包含很多岛礁,注意检查 shp 拓扑是否闭合,未闭合的 shp 生成掩膜时会出现锯齿状空洞。

Apply Mask 操作里输入文件选择正射校正结果,掩膜文件选择刚才生成的 mask,输出文件命名为 water.dat。通过这个操作,陆地区域全部被掩膜掉,只留下水体像元。打开 water.dat 检查一下,重庆的长江、嘉陵江以及周边水库和池塘应该清晰可见,而山地和城区全部变成忽略值。如果水体区域出现了断裂,说明输入的 bool mask 值方向反了,把陆地保留了下来而水体被裁掉,重新回到 Build Mask 检查。

5. 避坑排查:定标过暗、黑边残留、掩膜反向等五个实测问题

5.1 辐射定标后影像整体偏暗,水体细节看不见

现象:做完辐射定标,打开影像发现整体亮度很低,水体区域几乎全黑,拉直方图也救不回来。 原因:定标时输出数据类型选成了 Byte 整型,辐亮度数值范围大,整型存储会把大量中间灰度截断,导致暗部被压缩成一团。 解决:重新做辐射定标,把 Output Data Type 设为 Float,输出浮点型辐亮度。操作不麻烦,但很多人第一步顺手默认就选了 Byte,后续 FLAASH 读取这样的“假辐亮度”数据也会出错。记住一条:ENVI 做定标和大气校正,输出类型永远是 Float,宁可后面存储空间多占一点,也不能为了省空间把物理量级截断。

5.2 FLAASH 大气校正后输出的反射率出现大量负值

现象:大气校正后打开影像,水体区域大量像元反射率为负,整幅图像是打了马赛克。 原因:两种可能:一是初始能见度设置过大,把大气浑浊度估计得过于干净,导致校正算法把大气散射分量过度扣除;二是气溶胶模型选错,把乡村模型用到城市污染严重的区域,或者反过来,把城市模型用到干净山区。 解决:重新跑 FLAASH,把能见度从 40 km 调到 20 km,气溶胶模型选 Urban,然后查看水体区域在近红外波段的取值。通常把能见度调低一档后,负值大幅减少。注意重庆这种多雾低能见度城市,不要照抄平原地区的能见度参数。

5.3 正射校正后黑边没去掉,掩膜裁剪把黑边也包进来了

现象:Apply Mask 后输出的水域结果不是按 shp 边界裁剪,而是连同黑色背景一起形成了一个矩形区域,水体边缘全是锯齿状黑块。 原因:跳过 Edit Metadata 设置 Data Ignore Value 的步骤,黑色背景像元值为 0,在 Build Mask 的时候软件会把这个区域判定为有效数据,导致掩膜和黑边一起作用。 解决:回到正射校正结果,重新设置 Data Ignore Value=0,保存后再做 Build Mask 和 Apply Mask。如果已经做了掩膜,直接重新建掩膜重来,不要懒,黑边一旦参与掩膜计算,后期每步都会受影响。

5.4 波段计算时提示“Division by zero”或输出大量异常大值

现象:用 Band Math 计算 NDWI 时,输出结果里出现很多极端大正数或大负数,水体提取结果像雪花屏。 原因:分母为 0 或接近 0。正射校正结果里黑色背景像元值虽然是 0,但如果没设置 Data Ignore Value,波段计算器会正常参与运算,(b2+b4) 等于 0 的区域在浮点运算里就变成无穷大。 解决:两种方式,一是在波段计算前确保设置了 Data Ignore Value=0,ENVI 在波段运算时会自动跳过忽略值;二是把表达式写成 float 类型保护,例如 (float(b2)-float(b4))/(float(b2)+float(b4))。另外注意计算前先检查影像的坐标系和范围是否一致,不一致的影像做波段运算,结果也是错位相加。

5.5 水域裁剪后长江主河道断裂成很多小段

现象:水域 shp 本身连续,但裁剪后的河道出现多处断裂,尤其是江心洲两侧和桥梁附近。 原因:两类可能。一类是 shp 数据本身经过简化处理,在江心洲和桥梁位置存在多边形的自相交或空洞;另一类是掩膜生成时投影变换导致矢量边界和栅格像元不对齐,像元重采样后细窄河道被吞掉。 解决:对 shp 做拓扑修复,在 ArcGIS 里用 Repair Geometry 工具处理后再放进 ENVI;同时确认 shp 和影像的坐标系保持一致,影像通常是 UTM 投影,shp 也需要是同一个投影。投影不一致时先在 ArcGIS 里做投影转换,不要指望 ENVI 自动处理好。

6. 波段计算反演与 ArcGIS 制图:从水体指数到黑臭水/水华分级出图

6.1 NDWI 与 MNDWI 的波段计算器实现

在 ENVI 工具栏搜索 band,打开 Band Math 波段计算器,输入 NDWI 表达式并以 Map Variable to Input File 映射到 water.dat 裁剪后的水体影像。高分一号 WFV3 的波段顺序是:b1 蓝光、b2 绿光、b3 红光、b4 近红外。NDWI 公式用绿光和近红外波段:

# NDWI - 归一化差异水体指数(GF-1 WFV3 对应波段2绿光、波段4近红外) (float(b2)-float(b4))/(float(b2)+float(b4))

表达式里的 float() 运算符做浮点转换,避免整型相除时接近 0 的分母导致极值。NDWI 的物理含义是绿光反射率减去近红外反射率再除以二者之和,水体在绿光反射率高于近红外,所以水体区域的 NDWI 通常为正值,而植被和土壤为负值。长江和嘉陵江主河道应该呈现明显正值区,而支流和水库因为水深不同会有梯度差异。

MNDWI 的表达式则把绿光波段与中红外波段之间的差异换成了绿光与近红外,需要再确认一下 GF-1 WFV 没有中红外波段,所以严格意义上 MNDWI 在 WFV 数据上是没法直接算的。但很多教材项目里说的“MNDWI”其实是用绿光和红光波段的差值替代,表达式可以写成:

# MNDWI 替代表达式(用绿光和红光波段的差异增强水体信号) (float(b2)-float(b3))/(float(b2)+float(b3))

注意这是一种经验替代,不要把它当成严格意义上的 MNDWI 去发表论文。实际处理时我一般优先用 NDWI 做水体提取主判据,用绿光和红光差值来辅助区分高浊度水体,因为黑臭水体的典型特征是绿光和红光反射率关系发生改变。

6.2 黑臭水体和水华的阈值分割

NDWI 计算结果是一幅连续变化的灰度图,需要设定阈值把水体信息分类。水体区域的 NDWI 阈值定位在 0 附近,但重庆的长江水体含沙量较高,近红外反射率偏大,实际 NDWI 在 0 附近的像元很多。查看 NDWI 直方图,在 0 到 0.2 之间会看到两个波峰,第一个波峰对应含沙量较高的浑浊水体,第二个波峰对应较清澈的水体,阈值选在两个波峰之间的谷底处。

黑臭水体的判别在光谱上表现为蓝紫光反射率显著升高,绿光反射率相对降低,红光反射率升高,水体整体呈现黑暗甚至深褐色。实际操作中我习惯用两个波段组合来判别:一是用蓝光波段和绿光波段的差值 (b1-b2) 提取发黑发暗的水体,二是用红光波段和绿光波段的差值与 NDWI 值组合分级。在 ENVI Band Math 里输入分类表达式:

# 黑臭水体判别条件(波段值均为反射率,条件是经验阈值,需要按区域微调) ((float(b1)-float(b2)) > 0.02) AND ((float(b4)) < 0.05) AND (NDWI > -0.1)

水华区域则表现为叶绿素浓度高,近红外反射率异常升高,所以用 NDVI 更好:

# 水华区域判别(近红外高于红光,叶绿素富集) (float(b4)-float(b3))/(float(b4)+float(b3))

水华区域 NDVI 明显高于周围水体,用直方图分割找异常高值区,一般阈值在 0.2 以上。使用 ArcGIS 的栅格计算器配合条件判断,把正常水体、黑臭水体和疑似水华区域合并成一个三值栅格,方便出图和统计。建议保留原 NDWI 连续值图层用于后续定量分析,而分类结果只作为可视化展示。

6.3 ArcGIS 制图与面积统计技巧

反演结果最后要落到 ArcGIS 出图。先把 ENVI 输出的分类栅格导出成 GeoTIFF,然后用 ArcMap 加载,设置不同分类值的符号:正常水体用浅蓝色,黑臭水体用深褐色,疑似水华区域用黄绿色。图例、比例尺、指北针和文字注记逐一加上,重庆长江和嘉陵江交汇处的点位标注清楚,图名写成“重庆部分区域黑臭水体与水华分布反演图”,再导出成 JPEG 或 PDF 作为成果输出。

面积统计在 ArcGIS 里用栅格转面工具把分类栅格转成多边形要素类,然后新建字段计算 Area,按分类值汇总面积即可。注意利用栅格转面前先把分类栅格的 Data Ignore Value 或 NoData 值清理好,不然统计会把无效像元算进去。如果是论文或报告,建议同时在 ArcGIS 里生成采样点的 NDWI 值提取,用来和实地水质监测数据做相关性分析,这样反演成果的可靠性更高。

我自己的习惯是:无论项目要求多紧,从数据下载到反演出图每一步都保留过程文件。数据下载的 XML、定标后的 rad.dat、大气校正后的 ref.dat、正射校正后的 ortho.dat、裁剪后的 water.dat、NDWI 中间结果和最终分类栅格,全部按阶段存好,每一层的元数据里都备注当时用的参数。尤其是原始数据下载时那个云覆盖量和影像的规范性检查,从那以后我每次做重庆区域的水体反演都会强制走一遍“先体检再开工”的流程,绝不在没确认云量分布和传感器型号的情况下盲目开跑。做遥感反演,翻车往往不在最后一步,而在最开始那些你没耐心检查的参数上。希望这篇笔记能帮你少走点弯路,顺利把整套流程复现出来。

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

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进