高分一号水体反演全流程:从辐射定标到黑臭水华分类
简介这份资源面向遥感与GIS方向的初学者及环境监测从业者围绕高分一号影像在重庆部分区域的水体监测展开完整覆盖从数据获取、预处理到反演出图的全流程重点解决黑臭水体与水华现象的识别问题适合作为ArcGIS与ENVI联合实操的入门范例。压缩包内共1个docx文档约13.95MB以图文步骤形式记录操作流程便于对照软件界面逐步复现。内容涵盖陆地观测卫星数据服务网站的数据检索与下载、ENVI中的辐射定标、基于FLAASH的大气校正、RPC正射校正与黑背景值去除、水域掩膜裁剪以及利用Band Math计算NDWI、MNDWI等水体指数完成反演最后在ArcGIS中符号化分类并输出专题图。已有559人学习读者可借此掌握遥感水体信息提取的完整技术链路与常见参数设置思路。1. 高分一号水体反演从一景原始影像到黑臭与水华分布图拿到一景高分一号 WFV 影像很多人第一反应是直接丢进 ENVI 算 NDWI结果水体边界糊成一片黑臭段和水华段根本分不开。问题不在指数本身而在前面的辐射定标、大气校正、正射校正三步没做扎实——DN 值没转成反射率大气散射把蓝波段抬高一截正射残差又让水体边界和 shp 对不上。这套流程把数据下载、ENVI 预处理、Band Math 反演、ArcGIS 出图串成一条完整链路适合做水环境遥感监测的从业者也适合刚接触 GF-1 影像、想跑通一次端到端流程的新手。下面按我实际拆解的步骤走参数和坑都标出来。2. 数据获取与预处理链路GF-1 WFV 的四步定标逻辑高分一号 WFV 传感器幅宽大、重访周期短做区域水体监测性价比高但它不是定量遥感载荷原始产品只给 DN 值直接算指数等于拿没校准的尺子量长度。整条预处理链路的目的是把 DN 变成有物理意义的地表反射率并落到统一的地理坐标系里后面反演才有可比性。2.1 数据检索与下载的参数取舍在陆地观测卫星数据服务平台检索时几个参数直接决定后面能不能用。时间窗口按水华高发期选一般夏秋两季云覆盖量我一般卡在 10% 以下超过 20% 的水体区域基本被云吃掉传感器选 WFV它有蓝、绿、红、近红外四个波段正好够算 NDWI 和 MNDWI。下载下来通常是 tar 包解压后一个文件夹里带多光谱 tiff、全色 tiff、RPC 文件和元数据 xml别只拷 tiffRPC 文件后面正射校正要用。提示下载前先看检索结果里的快视图确认目标水域没被云或者条带噪声盖住省得下完几十景再筛。2.2 辐射定标把 DN 转成辐亮度ENVI 里搜 rad 打开 Radiometric Calibration选多光谱数据。参数上Calibration Type 选 RadianceOutput Interleave 保持 BILScale Factor 不勾保持默认 1.0输出数据类型选 Float。定标系数从元数据 xml 里自动读如果软件没读到手动填 Gain 和 Offset。# ENVI 命令行方式做辐射定标IDL 环境 envi_open_file, GF1_WFV_multispectral.tiff, r_fidfid envi_doit, RadiometricCalibration, $ FIDfid, $ CAL_TYPE0, $ ; 0Radiance, 1Reflectance OUTPUT_INTERLEAVE1, $ ; BIL SCALE_FACTOR1.0, $ OUT_NAMEradiance.dat逻辑说明CAL_TYPE0 输出辐亮度这是 FLAASH 大气校正的输入要求如果直接选 ReflectanceFLAASH 会报错说输入不是辐亮度。SCALE_FACTOR 保持 1.0 是为了不引入额外缩放后面 FLAASH 自己会处理量纲。定标完检查一下直方图辐亮度值应该在合理范围内如果全是 0 或者异常大多半是定标系数没读对。2.3 FLAASH 大气校正光谱曲线文件别选错搜 fla 打开 FLAASH Atmospheric Correction输入辐射定标结果。关键参数Sensor Type 选 UNKNOWN-MSIGround Elevation 填研究区平均高程重庆部分区域按 0.3 km 左右填Flight Date 和 Time 从元数据抄Atmospheric Model 选 Mid-Latitude SummerAerosol Model 选 RuralInitial Visibility 填 40 km 左右。多光谱设置里点 Multispectral SettingsFilter Function File 点 Open → Spectral Library选 gf1-wfv3 光谱曲线文件。这一步是血泪经验选错传感器光谱响应文件反演出来的反射率整体偏移水体指数阈值全乱。高级设置不动直接 Apply。# FLAASH 调用示例IDL envi_doit, FLAASH, $ FIDfid_rad, $ OUT_NAMEflaash_out.dat, $ SENSOR_TYPEUNKNOWN-MSI, $ GROUND_ELEV0.3, $ FLIGHT_DATE2023-08-15, $ FLIGHT_TIME03:20:00, $ ATMOSPHERIC_MODELMid-Latitude Summer, $ AEROSOL_MODELRural, $ INITIAL_VISIBILITY40.0, $ FILTER_FUNCTION_FILEgf1-wfv3.sli参数说明GROUND_ELEV 单位是 km填错会导致大气路径长度估算偏差FLIGHT_TIME 用 UTC 时间不是北京时间差 8 小时会让太阳天顶角算错反射率整体偏亮或偏暗。校正完对比一下校正前后的光谱曲线水体在蓝波段应该明显下降如果反而升高检查光谱曲线文件是不是选成了别的传感器。2.4 正射校正与黑边去除搜 RPC 打开 RPC Orthorectification Using Reference Image输入大气校正结果、参考影像用同区域已正射的底图或 DEM、高程数据。参数DEM 选 SRTM 或 ASTER 都行Resampling 选 Bilinear输出像元大小按 WFV 原始 16 m 填。输出后双击图层 → Edit Metadata → Add Metadata Items找 Data Ignore Value设为 0黑边就透明了。这一步不做后面裁剪水域时黑边会被当成有效值参与计算Band Math 结果在边界处出现异常高值。2.5 水域掩膜裁剪搜 mask 打开 Build Mask载入正射结果和重庆水域 shp输出 mask 文件再搜 mask 打开 Apply Mask输入正射结果掩膜选刚生成的 mask输出 water.dat。参数上Mask Value 设 0 表示掩掉非水域Outside Value 也设 0。裁剪完检查一下水域范围是否完整如果 shp 坐标系和影像不一致掩膜会整体偏移这是最常见的翻车点。3. 黑臭与水华反演Band Math 公式与阈值怎么定预处理做完影像已经是地表反射率接下来用波段运算提取水体指数再在黑臭和水华两个维度上做区分。这里的关键不是公式多复杂而是阈值怎么定、波段顺序别搞反。3.1 NDWI 与 MNDWI 的公式选择ENVI 里搜 band 打开 Band Math输入公式。NDWI 用绿和近红外# NDWI绿波段(b2) 和 近红外(b4) (float(b2) - float(b4)) / (float(b2) float(b4))MNDWI 用绿和短波红外但 GF-1 WFV 没有短波红外波段所以这个项目里只能用 NDWI或者用绿-红组合做替代水体指数。这是载荷限制不是流程缺陷。公式里 float() 必须加否则整型相减会截断结果全是 0 或 1。映射波段时点 Map Variable to Input File选 water.datb2 对应绿波段b4 对应近红外。GF-1 WFV 波段顺序是蓝、绿、红、近红外别按 Landsat 的习惯去映射。3.2 黑臭水体与水华的区分逻辑NDWI 只能分出“是不是水”分不出“水好不好”。黑臭水体在可见光波段反射率整体偏低、波形平缓水华则在近红外抬升明显、绿波段出现反射峰。实操里我用两个辅助指数# 黑臭参考指数红/绿比值黑臭段偏高 float(b3) / float(b2) # 水华参考指数近红外与红波段差值归一化 (float(b4) - float(b3)) / (float(b4) float(b3))逻辑说明黑臭水体悬浮物和有机物高红光相对绿光反射增强比值大于 1 的像元重点排查水华区域近红外反射被藻类抬升归一化差值明显高于普通水体。两个指数叠加 NDWI 做决策树分类比单阈值靠谱。3.3 阈值标定与分类后处理阈值不能拍脑袋。我一般在水域范围内选几个已知采样点正常水、黑臭段、水华段各若干统计各指数均值取类间中点做阈值。ENVI 里用 ROI 工具画点再 Statistics 看均值。分类结果用 Classification → Post Classification → Majority/Minority Analysis 做 3×3 众数滤波去掉椒盐噪声。滤波核别开太大5×5 以上会把小面积水华斑块抹掉。4. ArcGIS 出图与统计符号化、面积量算与版面反演结果在 ENVI 里是栅格分类图直接出图不好看也不方便量算面积。导到 ArcGIS 里做符号化和统计是这套流程的收尾环节。4.1 分类栅格导入与符号化ENVI 分类结果存成 tiff 或 ENVI 标准格式ArcGIS 直接拖进来。图层属性 → Symbology → Unique Values按类别值赋色正常水体用蓝色黑臭水体用深灰或褐色水华用绿色。色带别用默认彩虹色审图的人看着累。如果分类值是 1/2/3 这种编码在 Symbology 里手动改 Label别让图例显示成数字。4.2 面积统计与变化对比打开属性表加字段 Area_km2用 Calculate Geometry 按投影坐标系算面积别用地理坐标系算出来是平方度。统计各类面积占比导出表格。# ArcPy 批量统计分类面积 import arcpy arcpy.env.workspace rD:\water_result.gdb fc water_class with arcpy.da.SearchCursor(fc, [Class, SHAPEAREA]) as cursor: stats {} for cls, area in cursor: stats[cls] stats.get(cls, 0) area for k, v in stats.items(): print(f类别 {k}: {v/1e6:.2f} km²)逻辑说明SHAPEAREA 返回投影坐标系下的面积平方米除以 1e6 转平方公里。如果要素类是地理坐标系这个值不可信先 Project 再统计。4.3 版面输出与出图检查Layout 视图里加图例、比例尺、指北针图例文字改成中文类别名。导出分辨率至少 300 dpi格式选 tiff 或 pdf。出图前检查三件事水域边界和底图是否套合、图例类别是否和栅格值一一对应、面积统计数字和属性表是否一致。5. 避坑与排查五个让反演结果翻车的细节5.1 大气校正后反射率整体偏亮或偏暗现象FLAASH 输出反射率曲线形状不对水体在蓝波段不降反升。原因Flight Time 填了北京时间而不是 UTC太阳天顶角算错或者光谱曲线文件选成了其他传感器。解决从元数据 xml 里核对成像时间确认是 UTC重新选 gf1-wfv3 光谱曲线文件校正后对比已知水体的光谱形态。5.2 正射校正后影像和水域 shp 对不上现象掩膜裁剪出来的水域范围整体偏移几百米。原因shp 坐标系和影像坐标系不一致或者 RPC 校正时参考影像选错。解决统一转到同一投影坐标系如 CGCS2000 高斯克吕格检查 RPC 文件的角点坐标是否和影像匹配。5.3 Band Math 结果全是 0 或异常值现象NDWI 计算结果大片为 0或者出现极大值。原因波段映射顺序错了把红波段当绿波段或者没加 float() 导致整型截断。解决核对 GF-1 WFV 波段顺序蓝绿红近红外公式里所有波段变量套 float()映射后先看直方图确认数值范围在 -1 到 1 之间。5.4 黑边参与运算导致边界异常现象水域边缘出现一圈高值或低值异常带。原因正射校正后的 Data Ignore Value 没设成 0黑边被当成有效反射率参与 Band Math。解决正射结果元数据里把 Data Ignore Value 设为 0或者在 Band Math 前用掩膜把黑边去掉。5.5 分类结果椒盐噪声严重现象分类图上一堆孤立小斑块水华和黑臭混在一起。原因单像元阈值分类没有空间约束传感器噪声被放大。解决先做 3×3 众数滤波再结合 ROI 采样点重新标定阈值如果噪声集中在条带区域检查原始影像是否有条带缺陷。6. 进阶技巧用 ROI 采样把阈值从玄学变成可复现阈值定得好不好直接决定反演结果能不能用。我早期做这类项目阈值全靠试换个区域就得重调后来改成 ROI 采样标定才算把这件事从玄学拉回工程。具体做法在 ENVI 里对预处理后的影像开 ROI 工具按正常水体、黑臭段、水华段三类各画 15 到 20 个样本点尽量分散在不同位置避开岸边混合像元和云影。然后对每个 ROI 统计 NDWI、红绿比值、近红外归一化差值三个指数的均值与标准差。类别NDWI 均值红绿比值均值近红外归一化均值正常水体0.35 左右0.85 左右负值黑臭水体0.20 左右1.10 以上接近 0水华区域0.30 左右0.90 左右明显正值阈值取两类均值的中点比如 NDWI 水体/非水体阈值取正常水体和背景的中点黑臭判定用红绿比值大于两类中点。这样换一景影像只要重画 ROI 重算均值阈值就跟着数据走不用凭感觉调。验证环节我一般留 20% 的样本点不参与标定分类完拿这些点算混淆矩阵总体精度低于 85% 就回去检查预处理或者重新选样本。这套习惯是从一次返工里逼出来的——当时没留验证点图交上去被问“你怎么证明分类是对的”答不上来。从那以后我每次做反演都强制走一遍 ROI 标定加留出验证哪怕时间紧也不省这一步。希望帮到你。本文还有配套的精品资源点击获取