简介本资源面向遥感与GIS方向的初学者及环境监测从业者提供利用高分一号影像识别重庆部分区域黑臭水体与水华的完整技术方案覆盖从数据获取到成图的全流程涉及ArcGIS与ENVI两大工具。压缩包内共1个docx文档约13.95MB以图文步骤形式记录操作要点便于对照复现。内容从陆地观测卫星数据服务网站检索下载高分一号WFV数据入手依次讲解ENVI中的辐射定标、FLAASH大气校正、RPC正射校正与黑背景值去除再通过Build Mask与Apply Mask完成水域裁剪最后用Band Math计算NDWI、MNDWI等指数反演水体并在ArcGIS中符号化区分正常水体、黑臭水体与水华区导出专题图。已有557人学习适合需要掌握遥感水体监测完整链路、积累项目实操经验的读者参考。1. 高分一号影像做黑臭水体与水华反演从数据到出图的完整链路手里只有一景高分一号GF-1多光谱影像要同时回答两个问题哪些河段属于黑臭水体哪些湖库区域出现了水华。这不是两套独立流程而是同一套预处理管线分叉出两条反演支路。黑臭水体靠的是水质参数与反射率的经验关系水华靠的是藻类色素在近红外波段的强反射特征。两者共用辐射定标、大气校正、水体掩膜这三步之后才各走各的指数计算。适合做这件事的人是手头有GF-1 L1级数据、会一点Python或ENVI、但没跑通过完整反演链路的遥感从业者。下面按我实际跑通的顺序拆开讲参数和阈值都给具体值拿到数据就能照着复现。2. 数据准备与预处理GF-1 L1到地表反射率的四步走2.1 为什么GF-1的波段配置适合做这两类反演GF-1搭载的PMS传感器有四个多光谱波段蓝0.45-0.52μm、绿0.52-0.59μm、红0.63-0.69μm、近红外0.77-0.89μm空间分辨率8m多光谱。黑臭水体的核心判别依据是水体中悬浮物、有色可溶性有机物CDOM和营养盐的变化这些会改变水体在绿-红波段的反射率斜率。水华则是藻类叶绿素a在近红外波段的强反射GF-1的近红外波段刚好覆盖这个特征区间。16m幅宽单景多光谱对城市河段和中小型湖库够用但做大型湖泊需要多景拼接。常见做法是先把L1级数据做辐射定标转成表观辐亮度再用FLAASH或6S做大气校正。我一般用FLAASH因为GF-1的波段响应函数在ENVI里有现成配置。2.2 辐射定标与大气校正的具体参数辐射定标用ENVI的Radiometric Calibration工具或者用Python的rasterio手动算。定标公式是L Gain × DN OffsetGain和Offset从GF-1的XML元数据里读。下面这段Python代码直接读元数据完成定标import rasterio import xml.etree.ElementTree as ET import numpy as np def radiometric_calibration(img_path, xml_path, band_index): 对GF-1 L1级单波段做辐射定标 img_path: 单波段tif路径 xml_path: 对应的元数据XML band_index: 波段序号1蓝2绿3红4近红外 tree ET.parse(xml_path) root tree.getroot() # 从XML里读定标系数GF-1的标签名通常是Gain和Offset gains [float(e.text) for e in root.iter(Gain)] offsets [float(e.text) for e in root.iter(Offset)] with rasterio.open(img_path) as src: dn src.read(1).astype(np.float32) # 定标表观辐亮度 L dn * gains[band_index-1] offsets[band_index-1] return L逻辑说明GF-1的元数据XML里Gain和Offset是按波段顺序排列的band_index从1开始对应蓝波段。定标后的单位是W/(m²·sr·μm)。参数说明如果XML里读不到Gain检查标签名是否被改成了“CalibrationGain”之类不同地面站出的元数据标签有差异。定标完的辐亮度再送FLAASH设置传感器高度、地面高程、气溶胶模型。气溶胶模型选Rural能见度按当天天气填一般取20-40km。大气校正后得到地表反射率值域在0-1之间水体反射率通常低于0.1。2.3 水体掩膜的三种做法与选择依据水体掩膜是后面所有反演的前提。GF-1影像里水体、植被、建筑混杂不掩膜会把建筑阴影误判成黑臭水体。常用三种方法NDWI阈值法、单波段阈值法、面向对象分割。我一般用NDWI因为GF-1的绿波段和近红外波段对水陆差异最敏感。NDWI (Green - NIR) / (Green NIR)阈值取0.1到0.2之间。下面代码做NDWI并生成掩膜def water_mask(green, nir, threshold0.15): 基于NDWI生成水体掩膜 green: 绿波段地表反射率 nir: 近红外波段地表反射率 threshold: 经验阈值GF-1常用0.15 ndwi (green - nir) / (green nir 1e-10) mask ndwi threshold return mask, ndwi逻辑说明分母加1e-10防止除零。阈值0.15是GF-1在长江中下游地区跑出来的经验值北方干旱区水体浑浊阈值要降到0.05-0.1。参数说明如果掩膜后水体破碎用形态学闭运算填小孔如果掩膜把水华区域漏掉了因为水华的NDWI会偏低需要把阈值降到0.05再补一次。这一步的坑在于水华覆盖区NDWI可能接近0单纯靠NDWI会漏所以后面水华反演要单独处理。3. 黑臭水体反演从指数计算到等级划分3.1 黑臭水体的光谱响应机理与指数选择黑臭水体的“黑”来自悬浮物和CDOM对可见光的强吸收“臭”来自厌氧分解产生的硫化物和氨氮。遥感能捕捉的是“黑”对应的反射率变化绿波段反射率下降红波段相对抬升近红外被水体强烈吸收。所以黑臭水体的反射率曲线在绿-红区间斜率变陡。常用的指数有三个黑臭水体指数BOI、归一化黑臭指数NDBI的变体、以及基于绿红波段比值的指数。我一般用BOI (Red - Green) / (Red Green)因为GF-1的红绿波段对黑臭水体的响应最直接。BOI值越高黑臭程度越重。但BOI受悬浮物浓度干扰需要配合透明度或悬浮物指数做校正。3.2 用Python计算BOI并划分黑臭等级下面代码在掩膜后的水体区域计算BOI并按阈值分三级def calculate_boi(red, green, water_mask): 计算黑臭水体指数BOI red: 红波段反射率 green: 绿波段反射率 water_mask: 水体掩膜 boi (red - green) / (red green 1e-10) boi[~water_mask] np.nan # 非水体区域置NaN return boi def classify_black_odor(boi): 按BOI值划分黑臭等级 返回0非黑臭1轻度2中度3重度 level np.zeros_like(boi, dtypenp.uint8) level[(boi 0.05) (boi 0.15)] 1 level[(boi 0.15) (boi 0.25)] 2 level[boi 0.25] 3 return level逻辑说明BOI的阈值0.05、0.15、0.25是我在几个城市河段用实测数据回归出来的。轻度黑臭对应BOI 0.05-0.15这个区间的水体感官上已经发暗但还没到明显臭味。参数说明如果研究区水体本底浑浊阈值整体上移0.05如果水体清澈但有机污染重阈值下移0.03。分级后做众数滤波去掉椒盐噪声窗口取3×3。这一步的坑在于桥梁阴影的BOI会异常高因为阴影在红波段反射率也低但绿波段更低导致BOI虚高。解决办法是用NDWI掩膜后再叠加一个亮度阈值把反射率整体低于0.02的像元剔除。3.3 用实测数据校正BOI阈值的回归方法如果手头有实测的黑臭等级数据比如溶解氧、氨氮、透明度可以用线性回归把BOI和实测指标挂钩。下面代码做一元线性回归并输出校正后的阈值from scipy import stats import numpy as np def calibrate_boi_threshold(boi_samples, do_samples): boi_samples: 采样点BOI值数组 do_samples: 对应溶解氧实测值 返回回归方程和校正后的分级阈值 slope, intercept, r_value, p_value, std_err stats.linregress(boi_samples, do_samples) # 按溶解氧阈值反推BOI阈值DO2为重度黑臭2-4为中度4-6为轻度 boi_thresholds [(do - intercept) / slope for do in [6, 4, 2]] return slope, intercept, r_value, boi_thresholds逻辑说明溶解氧低于2mg/L是重度黑臭的国标阈值2-4是中度4-6是轻度。用回归方程反推BOI阈值比拍脑袋定阈值靠谱。参数说明r_value低于0.6说明BOI和溶解氧的相关性不够需要换指数或加其他波段。采样点至少15个才做回归少了容易过拟合。这一步的坑在于采样时间和卫星过境时间差超过3天水质已经变了回归出来的阈值不能用。4. 水华反演藻类指数与覆盖度估算4.1 水华在GF-1波段上的特征与指数选择水华蓝藻、绿藻爆发的遥感特征是近红外波段反射率急剧升高因为藻类细胞内的叶绿素和藻胆素在近红外有强反射。同时蓝波段因为藻蓝蛋白的吸收会下降。所以水华区域的NDVI会异常高甚至超过陆地植被。常用指数有两个NDVI (NIR - Red) / (NIR Red)以及FAI浮游藻类指数 NIR - [Red (SWIR - Red) × (λNIR - λRed)/(λSWIR - λRed)]。GF-1没有短波红外波段所以FAI用不了只能用NDVI或单波段近红外阈值。我一般用NDVI因为GF-1的NDVI对水华的响应足够敏感。水华的NDVI通常在0.1以上重度水华能到0.4。4.2 计算NDVI并提取水华覆盖度下面代码计算NDVI并用混合像元分解估算水华覆盖度def calculate_ndvi(nir, red, water_mask): 计算NDVI并提取水华区域 nir: 近红外反射率 red: 红波段反射率 water_mask: 水体掩膜 ndvi (nir - red) / (nir red 1e-10) ndvi[~water_mask] np.nan return ndvi def algae_coverage(ndvi, ndvi_water0.05, ndvi_algae0.35): 用线性混合模型估算水华覆盖度 ndvi_water: 纯水体的NDVI取0.05 ndvi_algae: 纯藻类的NDVI取0.35 coverage (ndvi - ndvi_water) / (ndvi_algae - ndvi_water) coverage np.clip(coverage, 0, 1) return coverage逻辑说明纯水体的NDVI在0附近取0.05是GF-1在清洁水体的经验值。纯藻类的NDVI取0.35是重度水华区的实测值。覆盖度0-0.2为轻度0.2-0.5为中度0.5以上为重度。参数说明如果研究区水华以绿藻为主ndvi_algae可以降到0.3以蓝藻为主升到0.4。这一步的坑在于水华边缘的NDVI受水体混合影响覆盖度估算偏低需要做3×3的均值滤波平滑。另外水生植被的NDVI和水华重叠需要用形状指数区分——水华呈片状连续分布水生植被呈条带状。4.3 水华与黑臭水体的重叠区处理同一景影像里水华区域和黑臭水体区域可能重叠因为水华死亡分解会消耗溶解氧导致水体发黑发臭。处理原则是先判水华再判黑臭。如果NDVI大于0.1优先归为水华NDVI小于0.1且BOI大于0.05归为黑臭。下面代码做优先级合并def merge_results(boi_level, algae_cov, ndvi): 合并黑臭和水华结果 优先级水华 黑臭 final np.zeros_like(boi_level, dtypenp.uint8) # 水华覆盖度大于0.2的区域标为水华用10表示 final[algae_cov 0.2] 10 # 非水华区域按黑臭等级赋值 mask_no_algae algae_cov 0.2 final[mask_no_algae] boi_level[mask_no_algae] return final逻辑说明用10表示水华0-3表示黑臭等级这样出图时用不同色带区分。参数说明如果水华覆盖度阈值取0.2导致水华面积偏大调到0.3偏小调到0.15。这一步的坑在于水华区域的BOI也会偏高因为藻类在红波段的反射率不低如果不做优先级合并水华会被误判成重度黑臭。5. 避坑与排查跑GF-1反演最容易翻车的五个地方5.1 大气校正后水体反射率出现负值现象FLAASH校正后绿波段反射率出现-0.01到-0.05的负值。原因气溶胶模型选错或者能见度填得过高导致大气程辐射扣除过度。解决把气溶胶模型从Rural改成Urban能见度降到15-25km。如果还有负值用ENVI的QUAC快速校正做对比QUAC不依赖气溶胶模型但精度略低。负值像元在后续指数计算前统一置为0.001避免除零和指数异常。5.2 NDWI掩膜把水华区域漏掉现象水华覆盖区的NDWI低于0.15被掩膜剔除导致水华反演结果空白。原因水华在近红外反射率高NDWI (Green - NIR)/(Green NIR)被拉低。解决做两次掩膜第一次用NDWI0.15提取清洁水体第二次用NDVI0.1提取水华两者取并集作为最终水体掩膜。代码里先算NDVI再算NDWI合并掩膜。5.3 BOI阈值在不同季节漂移现象夏季跑出来的黑臭等级比冬季高一个级别但实测水质没变。原因夏季水温高藻类活动强水体反射率整体抬升BOI基线漂移。解决按季节分别定阈值夏季阈值上移0.03-0.05冬季下移0.02。或者用温度归一化但GF-1没有热红外波段只能靠季节经验值。我一般建一个阈值查找表按月份索引。5.4 水华覆盖度在影像边缘异常高现象影像边缘的水华覆盖度突然跳到0.8以上但目视检查没有水华。原因边缘像元受邻近地物光谱混合影响NDVI虚高。解决做边缘缓冲把影像外围10个像元剔除。或者用众数滤波窗口取5×5把边缘异常值平滑掉。如果边缘是山区地形阴影也会导致NDVI异常需要叠加坡度掩膜。5.5 分级结果出现大量椒盐噪声现象黑臭等级图里单个像元跳变重度黑臭旁边紧挨着非黑臭。原因GF-1的8m分辨率下水体像元混合了岸边建筑或植被的光谱。解决先做3×3众数滤波再做连通域分析剔除面积小于10个像元的斑块。如果噪声还在检查大气校正是否引入了条带噪声用ENVI的Destriping工具处理。6. 出图与验证从分级栅格到可交付成果出图不是把栅格往ArcGIS里一扔就完事。我一般分三层底层是GF-1真彩色影像做底图中层是黑臭等级的半透明色带顶层是水华覆盖度的等值线。色带用红-橙-黄表示黑臭轻度到重度用绿-蓝表示水华轻度到重度。下面代码用matplotlib出图import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap def plot_results(rgb, boi_level, algae_cov, output_path): rgb: 真彩色影像数组形状(3, H, W) boi_level: 黑臭等级0-3 algae_cov: 水华覆盖度0-1 fig, axes plt.subplots(1, 3, figsize(18, 6)) # 真彩色 axes[0].imshow(rgb.transpose(1, 2, 0)) axes[0].set_title(GF-1 True Color) # 黑臭等级 cmap_boi ListedColormap([none, #FFFF00, #FFA500, #FF0000]) axes[1].imshow(boi_level, cmapcmap_boi, vmin0, vmax3) axes[1].set_title(Black-Odor Level) # 水华覆盖度 axes[2].imshow(algae_cov, cmapYlGnBu, vmin0, vmax1) axes[2].set_title(Algae Coverage) for ax in axes: ax.axis(off) plt.tight_layout() plt.savefig(output_path, dpi300, bbox_inchestight)逻辑说明真彩色用rgb数组的transpose把波段维度换到最后一维。黑臭等级用ListedColormap0值设为透明避免覆盖底图。水华覆盖度用YlGnBu色带0值偏黄1值偏蓝。参数说明dpi取300保证出图精度bbox_inchestight去掉白边。如果要做成ArcGIS的.lyr文件把分级栅格存成GeoTIFF在ArcGIS里加载后手动配色带色带值参考上面的十六进制码。验证方法分两种一是用同期实测数据做混淆矩阵计算Kappa系数Kappa大于0.7说明分级可信二是用高分辨率影像如Google Earth截图做目视比对随机选20个点检查。我一般两种都做实测数据不够时以目视比对为主。最后说一个血泪经验GF-1的L1级数据里不同日期的辐射定标系数可能不一样跑批量处理时一定要逐景读XML不要用固定系数。这个坑我踩过导致一景影像的黑臭等级整体偏高一个级别排查了两天才发现是定标系数用错了。希望帮到你。本文还有配套的精品资源点击获取
