拓冰网站定制 
简介一份以高分一号影像为例、面向遥感与GIS学习者的完整水体反演流程教程聚焦重庆部分区域的黑臭水体与水华识别覆盖数据获取、预处理、指数反演、专题制图到成果分析的全部环节并贯穿ArcGIS与ENVI两套软件的应用。包内仅含1个docx操作文档压缩包大小约13.95MB文档以分步操作说明配合参数设置的形式清晰拆解辐射定标、FLAASH大气校正、基于RPC的正射校正、掩膜裁剪等关键预处理步骤并详解利用Band Math计算NDWI/MNDWI水体指数、区分正常水体/黑臭水体/水华区以及在ArcGIS中完成符号化与地图输出的方法。除标准流程外还提供了数据下载时的传感器与云量选择技巧、数据无效值忽略设置等实用排错经验可直接用于复现或改造为自有区域的监测方案。目前已有556人学习适合需要系统上手遥感水体监测的本科生、研究生及环境监测相关从业者。1. 高分一号影像反演黑臭水体与水华一条能直接落地的重庆方案重庆的长江、嘉陵江干流水质尚可但次级支流和城市内河的返黑返臭一直没消停夏季库湾水华也常常一夜间冒出来。用高分一号16米分辨率的WFV影像做黑臭水体和水华反演是区县环境监测里性价比很高的一条路。高分一号的幅宽大、重返周期短在重庆这种云多雾重的地区能抓住晴空窗口本身就值钱。下面把整条流程写透从卫星影像下载、辐射定标、大气校正到水体提取、黑臭指数构建、水华判别最后落到专题出图。这里不走语义分割、也不指望transformer反演只用波段比值指数先把业务跑起来——反演不是玄学但参数设置确实有一点讲究。2. 高分一号WFV数据获取与预处理从影像下载到反射率的四个关键参数2.1 高分一号影像从哪里下载三个渠道和存档检索高分一号WFV影像在国内可以免费获取这是它能用来做常态化监测的前提。常见渠道有三个中国资源卫星应用中心的陆地观测卫星数据服务网站、地理空间数据云、国家对地观测科学数据中心。这些平台都需要先用单位或学校身份注册然后按传感器、载荷、时间范围检索。检索时选传感器GF-1、载荷WFV级别选L1A或L1级正射产品。时间范围上重庆地区优先选云量小于10%的景。不要只看一景就下载建议把目标区域前后一周内的条带全部浏览一遍选出云量最少、覆盖完整的那一景。历史卫星影像的存档检索特别适合做“返黑返臭”时间对比比如去年前年同一季节的影像能看出黑臭河段的动态变化。注意WFV有四台相机相邻条带拼接处会有辐射差异如果目标区域跨了两景尽量选单景覆盖实在不行就把两景都下载后面拼接时按云量和重叠区均值处理。2.2 辐射定标与大气校正ENVI里的参数配置和校验拿到L1A影像后第一步做辐射定标。高分一号的定标公式是辐亮度等于DN值乘以绝对定标系数再加偏移量这些系数会随影像文件一起给到有的平台打包在XML文件里有的以独立文本形式提供。ENVI里用Radiometric Calibration工具校准类型选Radiance输出单位为W/(m²·µm·sr)勾选“Scale factor”为1.0避免数值被压缩成0100区间后影响后续计算。第二步是大气校正。FLAASH是目前用得最多的方法参数设置有四个直接影响结果传感器类型选GF-1 WFV重庆主城平均海拔按200400米设置长江沿线过水面低一点问题不大大气模式按成像季节选Mid-Latitude Summer或Mid-Latitude Winter春秋季在两者之间选择靠近成像月份的模式气溶胶模式选Urban因为重庆城区和工业区气溶胶以人为源为主初始能见度在秋季晴天设为3040公里云雾天或春季设为20公里。参数项推荐值说明传感器类型GF-1 WFV四波段多光谱16米分辨率大气模式Mid-Latitude Summer/Winter按成像月份和纬度选择气溶胶模式Urban重庆城区人为气溶胶特征明显初始能见度2040 km雾天取小值晴天取大值水汽反演不开启WFV没有水汽波段反演会引入噪声大气校正完成后立刻做两个检查。一是看水体像元的近红外反射率是否接近0.010.03如果出现大面积负值多半是能见度给低了或者气溶胶模式选错。二是看植被像元的NDVI是否大于0.3如果植被指数偏低说明大气扣除过量需要回退参数重新跑。没有实测气象数据时改用QUAC做快速大气校正也能接受虽然精度略低但不容易出现负反射率。2.3 几何精校正与区域裁剪让反演结果对齐重庆政区边界L1A数据自带RPC有理多项式系数在ENVI里用RPC Orthorectification做正射校正DEM选30米SRTM。校正好坏取决于控制点的地形起伏重庆沟谷纵横影像边缘容易出现几像元的偏移如果后续要叠加河长制断面点位误差会很扎眼。校正后把影像叠到天地图影像上检查河道边界是否贴合如果偏差大于一个像元就在配准工具里手动加控制点重点选在桥梁、江心洲、库湾拐角这些清晰目标上。裁剪用重庆区县行政边界或者重点流域边界的shp文件。ENVI里用Subset Data from ROIs按矢量范围裁剪ArcGIS里用Extract by Mask也可以。这里有个容易被忽略的细节裁剪时不要正好卡在边界上向外保留30个像元的缓冲这样后续做水体矢量化时边界处的像元不会因为切片而出现不连续。ArcGIS裁剪影像的操作本身比ENVI直观但ENVI在处理十六位反射率数据时能保留更多信息推荐在ENVI里完成裁剪再导出。3. 水体信息提取重庆复杂山地下抑制阴影干扰的实用方法3.1 为什么重庆水体提取不能只靠NDWI单阈值重庆是立体城市长江、嘉陵江深切河谷主城区高差几十米到上百米太阳高度角稍低时山坡阴影直接投在水面附近。NDWI的公式是绿波段近红外/绿波段近红外它利用水体在绿光反射略高、近红外几乎全吸收的特征来增强水体。问题是阴影也具有同样的光谱形态——阴影里的绿光贡献比近红外高所以NDWI会把大片山体阴影误判成水体。这是重庆水体提取最容易翻车的地方。另一个麻烦是水体破碎。16米像元下长江干流没有问题但次级河流和塘坝只有几个像元宽单阈值分割要么把这些细小水体漏掉要么把周围山地一起圈进来。所以我的做法不是只设一个NDWI阈值而是组合NDWI、坡度、面积三个条件做掩膜。其中NDWI负责光谱识别坡度负责剔除山地阴影像元面积阈值负责去掉零碎孤立点。3.2 用GDAL计算水体掩膜NDWI、坡度与形态学后处理下面是一段可以直接跑的Python脚本。它读取预处理后的四波段GF-1影像计算NDWI结合DEM坡度剔除山体阴影最后输出水体掩膜GeoTIFF。波段顺序按ENVI输出默认为蓝、绿、红、近红外如果你自己做过波段重排先检查data.shape的波段顺序。from osgeo import gdal import numpy as np def extract_water(input_tif, output_tif, slope_tifNone, ndwi_threshold0.10, slope_threshold5.0): 从GF1_WFV预处理影像中提取水体掩膜。 波段顺序: 0Blue, 1Green, 2Red, 3NIR ds gdal.Open(input_tif) data ds.ReadAsArray().astype(np.float32) green, nir data[1], data[3] # NDWI1e-6防止分母为0 ndwi (green - nir) / (green nir 1e-6) mask np.ones_like(ndwi, dtypenp.uint8) mask[ndwi ndwi_threshold] 0 # 用坡度剔除山地阴影: 坡度阈值的像元不可能是静水水体 if slope_tif: ds_slope gdal.Open(slope_tif) slope ds_slope.ReadAsArray().astype(np.float32) ds_slope None mask[slope slope_threshold] 0 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(output_tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(mask) out_ds.FlushCache() out_ds None ds None if __name__ __main__: extract_water( input_tifGF1_WFV_prep.tif, output_tifwater_mask.tif, slope_tifslope.tif, ndwi_threshold0.10, slope_threshold5.0 )这段代码的逻辑分三块先算NDWI并做阈值分割再用坡度栅格把高坡地区的水体候选像元删掉最后把结果写成单波段Byte型GeoTIFF。参数方面NDWI阈值0.10是重庆秋季晴空影像的常见起点如果影像里江面大面积漏提说明阈值偏高降到0.05重跑如果掩膜里出现大量山腰斑块说明阈值太低配合坡度一起调整。坡度阈值5度看起来很理想但注意长江河谷两侧的阶地坡度也接近这个值所以坡度掩膜只用来剔除典型山体阴影不要贪心加大到10度以上。3.3 水体矢量化与细小沟渠合并从栅格掩膜到矢量工作底图水体掩膜栅格不能直接用于制图和空间统计需要矢量化。GDAL自带的gdal_polygonize.py一行命令就能完成gdal_polygonize.py water_mask.tif -f ESRI Shapefile water_polygons.shp生成的矢量面会包含很多细碎斑块。在ArcGIS里按面积字段筛选保留面积大于3个像元约768平方米的斑块小于这个面积的斑块绝大多数是阴影残留或混合像元。接着用“消除”工具把靠近大水面且边界共享的小斑块合并进相邻水体避免河道在视觉上断断续续。重庆还有一个特殊问题山区沟谷里的水体被山脊线隔开矢量结果里一个连续河流可能被切成几段。这一步把多个河段按“同一条河流”字段做溶解或者直接用河湖管理范围线叠加修边界。我个人习惯保留一份未融合的水体掩膜用于黑臭指数统计再保留一份融合后的矢量用于出图底图两份数据各有用途。4. 黑臭水体反演模型指数构建、阈值标定与实地数据校准4.1 黑臭水体的光谱机理与两类反演路线黑臭水体的光谱特征和普通水体差异很明确整个可见光波段的反射率明显偏低绿光波峰的抬升幅度被压扁红光和近红外区域的反射率曲线趋于平缓整体像一条被压低了的直线。原因是黑臭水体里含有大量有机碎屑、还原性物质和黑色悬浮颗粒对可见光的吸收强烈。这个特征给了遥感反演一个物理基础不是靠单个波段而是靠波段之间的形态差异来判别。实际项目里常见两类反演路线。第一类是经验统计模型用实测水质指标透明度、溶解氧、COD、氨氮和同步影像反射率做回归适合有监测任务的单位。第二类是光谱指数判别构造黑臭水体指数BOI阈值分割出水体黑臭等级。指数法的好处是无需大量实测数据就能出初判结果坏处是阈值需要标定各地各季节的反射率绝对水平不同。下面用的就是指数法为主、实测点校准为辅的混合路线。4.2 黑臭水体指数BOI的计算和分级黑臭水体指数BOI在文献里有多种变体有的用绿红比值有的用红和近红外组合我这里用很多项目里表现稳定的一种形式把红波段加近红外作为分子绿波段加近红外作为分母黑臭水体的绿光贡献相对更弱所以指数值偏高普通水体绿光贡献强指数值偏低。计算公式为# 读取预处理后的4波段反射率 blue, green, red, nir data[0], data[1], data[2], data[3] # 黑臭水体指数BOI boi (red nir) / (green nir 1e-6) # 先用0.85作为初判阈值方向是“越高越黑臭” odorous_level np.zeros_like(boi, dtypenp.uint8) odorous_level[boi 0.85] 3 # 重度黑臭 odorous_level[(boi 0.75) (boi 0.85)] 2 # 中度黑臭 odorous_level[(boi 0.65) (boi 0.75)] 1 # 轻度黑臭这里0.65、0.75、0.85只是一组示例初始值。凡是比值型指数绝对数值会随大气校正质量、季节和影像内部的辐射一致性变化。所以算完BOI后不要急着出专题图先看整个水体的BOI直方图找到两个峰之间的低谷位置把阈值设在谷底才符合这一景影像的实际情况。4.3 实地点位校准阈值用Excel和影像值做最简单的验证现场校准是黑臭反演里最值得花时间的一步。带着手机和GPS去现场记录点位每个点拍一张水色照片标注等级。回来后把点位的经纬度坐标做成shp文件用下面的代码提取对应像元的BOI值import geopandas as gpd import rasterio # 加载现场点位shp points gpd.read_file(field_samples.shp) # 读取BOI栅格采样点位对应的像元值 with rasterio.open(boi.tif) as src: sampled [v[0] for v in src.sample( [(geom.x, geom.y) for geom in points.geometry] )] points[boi_value] sampled points.to_file(field_samples_with_boi.shp)然后把BOI值和现场等级整理成Excel画一个箱线图。如果黑色水体的BOI值分布和普通水体明显分开说明指数方向是对的如果两组数据重叠严重说明该区域黑臭水体的光谱形态不符合这个指数假设就要换用绿红波段的其他比值形式。这里踩过坑越多越能体会“公式抄来的、阈值一定是自己标定的”这句话的分量。5. 反演全流程的五个典型坑云、负反射率、阴影、阈值漂移和混合像元5.1 云和云影没处理反演指数大面积异常现象水体提取结果里出现大片连续高值斑块黑臭指数异常偏高对照原始影像才发现是一团云及其影子。 原因高分一号WFV不下发标准云掩膜产品很多人在预处理时直接忽略云污染结果云顶反射率把比值指数推向极端值。 解决在计算NDWI和BOI之前先做一次云掩膜。经验方法是云在蓝波段反射率极高、植被指数极低可以用蓝波段阈值加NDVI阈值组合识别云影则表现为近红外被明显压低需要在影像上人工圈出可疑区域。ENVI的FLAASH工具里也有云掩膜选项跑大气校正时顺带输出Clound Mask这一步不要省。5.2 大气校正后水体出现负反射率指数计算全乱现象FLAASH跑完后检查水体像元近红外或红光反射率出现负值BOX计算时产生异常尖峰。 原因大气校正参数里能见度设置过小把大气贡献扣得太多水面本身的暗信号被扣成负数。重庆春季雾天多如果直接把能见度设成10公里基本会翻车。 解决把能见度提高到30公里以上重新校正或者改用QUAC。QUAC不需要输入气溶胶参数校正后的光谱形态比FLAASH稳虽然绝对值精度略低但对于比值型指数够用。实在不行在计算指数前列一个最小值钳制np.clip(data, 0.001, 1.0)能救回负数导致的除零和锯齿。5.3 山体阴影被当成水体重庆地形的经典误判现象NDWI掩膜里出现大量沿山脉走向的斑块面积比真实水体还大。 原因阴影和水体在绿、近红外两波段的光谱形态高度相似单靠NDWI根本分不开。 解决叠加坡度掩膜只能解决一部分还有一部分是缓坡上的建筑物阴影和桥梁阴影。我的组合方案是坡度大于5度直接剔除NDWI阈值不低于0.05最后再按水体形状的圆度或长宽比筛一遍。山体阴影的形状往往细长且不规则和河流的连续条带差异明显。把这两层叠加后误判能减少八成以上。5.4 黑臭阈值照搬文献换一张影像就翻车现象上一景影像用BOI大于0.75识别黑臭效果很好换到下一景0.75以上只有零星几个像元全区域没有黑臭。 原因比值型指数的绝对值会随太阳高度角、气溶胶、水体浑浊度季节性变化不同月份的影像之间阈值会漂移。文献发表的阈值是针对特定影像和特定季节的直接拿来用等于在赌运气。 解决每次处理新影像必须重新标定阈值。标定材料不需要很多现场拍几张水色照片就够了。先把明显清洁的水库和明显黑臭的沟渠在BOI上取值取两类点的中位数中间位置作为分界阈值。把历史影像的阈值记录在Excel里同季节的影像可以套用相近阈值跨季节必须重新来。5.5 16米分辨率下的混合像元小沟渠黑臭被平均掉现象现场看着发黑的窄沟渠在反演结果里完全没信号水体掩膜甚至没把它提出来。 原因高分一号WFV是16米像元重庆的黑臭沟渠很多只有5到15米宽一个像元里混合了水面、岸壁和植被光谱被平均成普通地表。 解决项目一开始就明确遥感反演的目标水体是水面宽度大于50米至少3个像元的河流、湖库和塘坝。窄沟渠交给高分二号2米影像或现场巡查来查漏。把这个限制写进报告结论里比事后解释为什么漏掉几条沟渠要体面得多。6. 水华反演与专题出图把反演结果变成能上报的图件6.1 用NDCI快速圈定藻华聚集区水华暴发时藻类叶绿素a浓度升高红光吸收增强、近红外反射抬升归一化叶绿素指数NDCI正好捕捉这个反差。公式是近红外红光/近红外红光GF-1没有710纳米中心波段不能用红边波段做精细反演但用来圈定中高浓度藻华聚集区已经足够。# 基于预处理后的反射率数据计算NDCI ndci (nir - red) / (nir red 1e-6) # 阈值0.01为常见起始值水华核心区常超过0.05 bloom np.where((ndci 0.01) (water_mask 1), 1, 0)阈值标定方法跟BOI一样先看NDCI直方图找双峰低谷。重庆夏季高温晴热天水华一般在库湾和回水区先聚集NDCI高值区的位置和实际水华照片对照后基本能对上。验证时拿当天或前一天的现场照片做时间匹配能对上就是有效反演。6.2 QGIS加载天地图核对空间分布并按制图规范出图反演结果必须叠加到真实地理底图上才具备说服力。QGIS里加载天地图影像作为参考底图把黑臭分级栅格和水华分布栅格按半透明叠加检查异常斑块是否落在河道内部、桥梁阴影是否误判。出图时按专题图标准配置色带用绿到红的分级表示黑臭程度水华用红色半透明图层叠加图例里写明指数名称、阈值和影像获取日期指北针、比例尺、坐标参考系一个都不能少。比例尺用1:50000到1:100000比较合适太细会暴露16米像元的锯齿。我习惯在出图前把阈值参数、验证点表、影像日期一起存档进项目文件夹。这样两个月后回访同一个水库时能直接调出上一期的阈值和色带保证两期图件可比。反演这件事做到最后拼的不是算法多先进而是阈值标定和过程记录的细致程度。希望帮到你。本文还有配套的精品资源点击获取 高分一号遥感反演黑臭水体与水华：重庆落地全流程
简介一份以高分一号影像为例、面向遥感与GIS学习者的完整水体反演流程教程聚焦重庆部分区域的黑臭水体与水华识别覆盖数据获取、预处理、指数反演、专题制图到成果分析的全部环节并贯穿ArcGIS与ENVI两套软件的应用。包内仅含1个docx操作文档压缩包大小约13.95MB文档以分步操作说明配合参数设置的形式清晰拆解辐射定标、FLAASH大气校正、基于RPC的正射校正、掩膜裁剪等关键预处理步骤并详解利用Band Math计算NDWI/MNDWI水体指数、区分正常水体/黑臭水体/水华区以及在ArcGIS中完成符号化与地图输出的方法。除标准流程外还提供了数据下载时的传感器与云量选择技巧、数据无效值忽略设置等实用排错经验可直接用于复现或改造为自有区域的监测方案。目前已有556人学习适合需要系统上手遥感水体监测的本科生、研究生及环境监测相关从业者。1. 高分一号影像反演黑臭水体与水华一条能直接落地的重庆方案重庆的长江、嘉陵江干流水质尚可但次级支流和城市内河的返黑返臭一直没消停夏季库湾水华也常常一夜间冒出来。用高分一号16米分辨率的WFV影像做黑臭水体和水华反演是区县环境监测里性价比很高的一条路。高分一号的幅宽大、重返周期短在重庆这种云多雾重的地区能抓住晴空窗口本身就值钱。下面把整条流程写透从卫星影像下载、辐射定标、大气校正到水体提取、黑臭指数构建、水华判别最后落到专题出图。这里不走语义分割、也不指望transformer反演只用波段比值指数先把业务跑起来——反演不是玄学但参数设置确实有一点讲究。2. 高分一号WFV数据获取与预处理从影像下载到反射率的四个关键参数2.1 高分一号影像从哪里下载三个渠道和存档检索高分一号WFV影像在国内可以免费获取这是它能用来做常态化监测的前提。常见渠道有三个中国资源卫星应用中心的陆地观测卫星数据服务网站、地理空间数据云、国家对地观测科学数据中心。这些平台都需要先用单位或学校身份注册然后按传感器、载荷、时间范围检索。检索时选传感器GF-1、载荷WFV级别选L1A或L1级正射产品。时间范围上重庆地区优先选云量小于10%的景。不要只看一景就下载建议把目标区域前后一周内的条带全部浏览一遍选出云量最少、覆盖完整的那一景。历史卫星影像的存档检索特别适合做“返黑返臭”时间对比比如去年前年同一季节的影像能看出黑臭河段的动态变化。注意WFV有四台相机相邻条带拼接处会有辐射差异如果目标区域跨了两景尽量选单景覆盖实在不行就把两景都下载后面拼接时按云量和重叠区均值处理。2.2 辐射定标与大气校正ENVI里的参数配置和校验拿到L1A影像后第一步做辐射定标。高分一号的定标公式是辐亮度等于DN值乘以绝对定标系数再加偏移量这些系数会随影像文件一起给到有的平台打包在XML文件里有的以独立文本形式提供。ENVI里用Radiometric Calibration工具校准类型选Radiance输出单位为W/(m²·µm·sr)勾选“Scale factor”为1.0避免数值被压缩成0100区间后影响后续计算。第二步是大气校正。FLAASH是目前用得最多的方法参数设置有四个直接影响结果传感器类型选GF-1 WFV重庆主城平均海拔按200400米设置长江沿线过水面低一点问题不大大气模式按成像季节选Mid-Latitude Summer或Mid-Latitude Winter春秋季在两者之间选择靠近成像月份的模式气溶胶模式选Urban因为重庆城区和工业区气溶胶以人为源为主初始能见度在秋季晴天设为3040公里云雾天或春季设为20公里。参数项推荐值说明传感器类型GF-1 WFV四波段多光谱16米分辨率大气模式Mid-Latitude Summer/Winter按成像月份和纬度选择气溶胶模式Urban重庆城区人为气溶胶特征明显初始能见度2040 km雾天取小值晴天取大值水汽反演不开启WFV没有水汽波段反演会引入噪声大气校正完成后立刻做两个检查。一是看水体像元的近红外反射率是否接近0.010.03如果出现大面积负值多半是能见度给低了或者气溶胶模式选错。二是看植被像元的NDVI是否大于0.3如果植被指数偏低说明大气扣除过量需要回退参数重新跑。没有实测气象数据时改用QUAC做快速大气校正也能接受虽然精度略低但不容易出现负反射率。2.3 几何精校正与区域裁剪让反演结果对齐重庆政区边界L1A数据自带RPC有理多项式系数在ENVI里用RPC Orthorectification做正射校正DEM选30米SRTM。校正好坏取决于控制点的地形起伏重庆沟谷纵横影像边缘容易出现几像元的偏移如果后续要叠加河长制断面点位误差会很扎眼。校正后把影像叠到天地图影像上检查河道边界是否贴合如果偏差大于一个像元就在配准工具里手动加控制点重点选在桥梁、江心洲、库湾拐角这些清晰目标上。裁剪用重庆区县行政边界或者重点流域边界的shp文件。ENVI里用Subset Data from ROIs按矢量范围裁剪ArcGIS里用Extract by Mask也可以。这里有个容易被忽略的细节裁剪时不要正好卡在边界上向外保留30个像元的缓冲这样后续做水体矢量化时边界处的像元不会因为切片而出现不连续。ArcGIS裁剪影像的操作本身比ENVI直观但ENVI在处理十六位反射率数据时能保留更多信息推荐在ENVI里完成裁剪再导出。3. 水体信息提取重庆复杂山地下抑制阴影干扰的实用方法3.1 为什么重庆水体提取不能只靠NDWI单阈值重庆是立体城市长江、嘉陵江深切河谷主城区高差几十米到上百米太阳高度角稍低时山坡阴影直接投在水面附近。NDWI的公式是绿波段近红外/绿波段近红外它利用水体在绿光反射略高、近红外几乎全吸收的特征来增强水体。问题是阴影也具有同样的光谱形态——阴影里的绿光贡献比近红外高所以NDWI会把大片山体阴影误判成水体。这是重庆水体提取最容易翻车的地方。另一个麻烦是水体破碎。16米像元下长江干流没有问题但次级河流和塘坝只有几个像元宽单阈值分割要么把这些细小水体漏掉要么把周围山地一起圈进来。所以我的做法不是只设一个NDWI阈值而是组合NDWI、坡度、面积三个条件做掩膜。其中NDWI负责光谱识别坡度负责剔除山地阴影像元面积阈值负责去掉零碎孤立点。3.2 用GDAL计算水体掩膜NDWI、坡度与形态学后处理下面是一段可以直接跑的Python脚本。它读取预处理后的四波段GF-1影像计算NDWI结合DEM坡度剔除山体阴影最后输出水体掩膜GeoTIFF。波段顺序按ENVI输出默认为蓝、绿、红、近红外如果你自己做过波段重排先检查data.shape的波段顺序。from osgeo import gdal import numpy as np def extract_water(input_tif, output_tif, slope_tifNone, ndwi_threshold0.10, slope_threshold5.0): 从GF1_WFV预处理影像中提取水体掩膜。 波段顺序: 0Blue, 1Green, 2Red, 3NIR ds gdal.Open(input_tif) data ds.ReadAsArray().astype(np.float32) green, nir data[1], data[3] # NDWI1e-6防止分母为0 ndwi (green - nir) / (green nir 1e-6) mask np.ones_like(ndwi, dtypenp.uint8) mask[ndwi ndwi_threshold] 0 # 用坡度剔除山地阴影: 坡度阈值的像元不可能是静水水体 if slope_tif: ds_slope gdal.Open(slope_tif) slope ds_slope.ReadAsArray().astype(np.float32) ds_slope None mask[slope slope_threshold] 0 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(output_tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(mask) out_ds.FlushCache() out_ds None ds None if __name__ __main__: extract_water( input_tifGF1_WFV_prep.tif, output_tifwater_mask.tif, slope_tifslope.tif, ndwi_threshold0.10, slope_threshold5.0 )这段代码的逻辑分三块先算NDWI并做阈值分割再用坡度栅格把高坡地区的水体候选像元删掉最后把结果写成单波段Byte型GeoTIFF。参数方面NDWI阈值0.10是重庆秋季晴空影像的常见起点如果影像里江面大面积漏提说明阈值偏高降到0.05重跑如果掩膜里出现大量山腰斑块说明阈值太低配合坡度一起调整。坡度阈值5度看起来很理想但注意长江河谷两侧的阶地坡度也接近这个值所以坡度掩膜只用来剔除典型山体阴影不要贪心加大到10度以上。3.3 水体矢量化与细小沟渠合并从栅格掩膜到矢量工作底图水体掩膜栅格不能直接用于制图和空间统计需要矢量化。GDAL自带的gdal_polygonize.py一行命令就能完成gdal_polygonize.py water_mask.tif -f ESRI Shapefile water_polygons.shp生成的矢量面会包含很多细碎斑块。在ArcGIS里按面积字段筛选保留面积大于3个像元约768平方米的斑块小于这个面积的斑块绝大多数是阴影残留或混合像元。接着用“消除”工具把靠近大水面且边界共享的小斑块合并进相邻水体避免河道在视觉上断断续续。重庆还有一个特殊问题山区沟谷里的水体被山脊线隔开矢量结果里一个连续河流可能被切成几段。这一步把多个河段按“同一条河流”字段做溶解或者直接用河湖管理范围线叠加修边界。我个人习惯保留一份未融合的水体掩膜用于黑臭指数统计再保留一份融合后的矢量用于出图底图两份数据各有用途。4. 黑臭水体反演模型指数构建、阈值标定与实地数据校准4.1 黑臭水体的光谱机理与两类反演路线黑臭水体的光谱特征和普通水体差异很明确整个可见光波段的反射率明显偏低绿光波峰的抬升幅度被压扁红光和近红外区域的反射率曲线趋于平缓整体像一条被压低了的直线。原因是黑臭水体里含有大量有机碎屑、还原性物质和黑色悬浮颗粒对可见光的吸收强烈。这个特征给了遥感反演一个物理基础不是靠单个波段而是靠波段之间的形态差异来判别。实际项目里常见两类反演路线。第一类是经验统计模型用实测水质指标透明度、溶解氧、COD、氨氮和同步影像反射率做回归适合有监测任务的单位。第二类是光谱指数判别构造黑臭水体指数BOI阈值分割出水体黑臭等级。指数法的好处是无需大量实测数据就能出初判结果坏处是阈值需要标定各地各季节的反射率绝对水平不同。下面用的就是指数法为主、实测点校准为辅的混合路线。4.2 黑臭水体指数BOI的计算和分级黑臭水体指数BOI在文献里有多种变体有的用绿红比值有的用红和近红外组合我这里用很多项目里表现稳定的一种形式把红波段加近红外作为分子绿波段加近红外作为分母黑臭水体的绿光贡献相对更弱所以指数值偏高普通水体绿光贡献强指数值偏低。计算公式为# 读取预处理后的4波段反射率 blue, green, red, nir data[0], data[1], data[2], data[3] # 黑臭水体指数BOI boi (red nir) / (green nir 1e-6) # 先用0.85作为初判阈值方向是“越高越黑臭” odorous_level np.zeros_like(boi, dtypenp.uint8) odorous_level[boi 0.85] 3 # 重度黑臭 odorous_level[(boi 0.75) (boi 0.85)] 2 # 中度黑臭 odorous_level[(boi 0.65) (boi 0.75)] 1 # 轻度黑臭这里0.65、0.75、0.85只是一组示例初始值。凡是比值型指数绝对数值会随大气校正质量、季节和影像内部的辐射一致性变化。所以算完BOI后不要急着出专题图先看整个水体的BOI直方图找到两个峰之间的低谷位置把阈值设在谷底才符合这一景影像的实际情况。4.3 实地点位校准阈值用Excel和影像值做最简单的验证现场校准是黑臭反演里最值得花时间的一步。带着手机和GPS去现场记录点位每个点拍一张水色照片标注等级。回来后把点位的经纬度坐标做成shp文件用下面的代码提取对应像元的BOI值import geopandas as gpd import rasterio # 加载现场点位shp points gpd.read_file(field_samples.shp) # 读取BOI栅格采样点位对应的像元值 with rasterio.open(boi.tif) as src: sampled [v[0] for v in src.sample( [(geom.x, geom.y) for geom in points.geometry] )] points[boi_value] sampled points.to_file(field_samples_with_boi.shp)然后把BOI值和现场等级整理成Excel画一个箱线图。如果黑色水体的BOI值分布和普通水体明显分开说明指数方向是对的如果两组数据重叠严重说明该区域黑臭水体的光谱形态不符合这个指数假设就要换用绿红波段的其他比值形式。这里踩过坑越多越能体会“公式抄来的、阈值一定是自己标定的”这句话的分量。5. 反演全流程的五个典型坑云、负反射率、阴影、阈值漂移和混合像元5.1 云和云影没处理反演指数大面积异常现象水体提取结果里出现大片连续高值斑块黑臭指数异常偏高对照原始影像才发现是一团云及其影子。 原因高分一号WFV不下发标准云掩膜产品很多人在预处理时直接忽略云污染结果云顶反射率把比值指数推向极端值。 解决在计算NDWI和BOI之前先做一次云掩膜。经验方法是云在蓝波段反射率极高、植被指数极低可以用蓝波段阈值加NDVI阈值组合识别云影则表现为近红外被明显压低需要在影像上人工圈出可疑区域。ENVI的FLAASH工具里也有云掩膜选项跑大气校正时顺带输出Clound Mask这一步不要省。5.2 大气校正后水体出现负反射率指数计算全乱现象FLAASH跑完后检查水体像元近红外或红光反射率出现负值BOX计算时产生异常尖峰。 原因大气校正参数里能见度设置过小把大气贡献扣得太多水面本身的暗信号被扣成负数。重庆春季雾天多如果直接把能见度设成10公里基本会翻车。 解决把能见度提高到30公里以上重新校正或者改用QUAC。QUAC不需要输入气溶胶参数校正后的光谱形态比FLAASH稳虽然绝对值精度略低但对于比值型指数够用。实在不行在计算指数前列一个最小值钳制np.clip(data, 0.001, 1.0)能救回负数导致的除零和锯齿。5.3 山体阴影被当成水体重庆地形的经典误判现象NDWI掩膜里出现大量沿山脉走向的斑块面积比真实水体还大。 原因阴影和水体在绿、近红外两波段的光谱形态高度相似单靠NDWI根本分不开。 解决叠加坡度掩膜只能解决一部分还有一部分是缓坡上的建筑物阴影和桥梁阴影。我的组合方案是坡度大于5度直接剔除NDWI阈值不低于0.05最后再按水体形状的圆度或长宽比筛一遍。山体阴影的形状往往细长且不规则和河流的连续条带差异明显。把这两层叠加后误判能减少八成以上。5.4 黑臭阈值照搬文献换一张影像就翻车现象上一景影像用BOI大于0.75识别黑臭效果很好换到下一景0.75以上只有零星几个像元全区域没有黑臭。 原因比值型指数的绝对值会随太阳高度角、气溶胶、水体浑浊度季节性变化不同月份的影像之间阈值会漂移。文献发表的阈值是针对特定影像和特定季节的直接拿来用等于在赌运气。 解决每次处理新影像必须重新标定阈值。标定材料不需要很多现场拍几张水色照片就够了。先把明显清洁的水库和明显黑臭的沟渠在BOI上取值取两类点的中位数中间位置作为分界阈值。把历史影像的阈值记录在Excel里同季节的影像可以套用相近阈值跨季节必须重新来。5.5 16米分辨率下的混合像元小沟渠黑臭被平均掉现象现场看着发黑的窄沟渠在反演结果里完全没信号水体掩膜甚至没把它提出来。 原因高分一号WFV是16米像元重庆的黑臭沟渠很多只有5到15米宽一个像元里混合了水面、岸壁和植被光谱被平均成普通地表。 解决项目一开始就明确遥感反演的目标水体是水面宽度大于50米至少3个像元的河流、湖库和塘坝。窄沟渠交给高分二号2米影像或现场巡查来查漏。把这个限制写进报告结论里比事后解释为什么漏掉几条沟渠要体面得多。6. 水华反演与专题出图把反演结果变成能上报的图件6.1 用NDCI快速圈定藻华聚集区水华暴发时藻类叶绿素a浓度升高红光吸收增强、近红外反射抬升归一化叶绿素指数NDCI正好捕捉这个反差。公式是近红外红光/近红外红光GF-1没有710纳米中心波段不能用红边波段做精细反演但用来圈定中高浓度藻华聚集区已经足够。# 基于预处理后的反射率数据计算NDCI ndci (nir - red) / (nir red 1e-6) # 阈值0.01为常见起始值水华核心区常超过0.05 bloom np.where((ndci 0.01) (water_mask 1), 1, 0)阈值标定方法跟BOI一样先看NDCI直方图找双峰低谷。重庆夏季高温晴热天水华一般在库湾和回水区先聚集NDCI高值区的位置和实际水华照片对照后基本能对上。验证时拿当天或前一天的现场照片做时间匹配能对上就是有效反演。6.2 QGIS加载天地图核对空间分布并按制图规范出图反演结果必须叠加到真实地理底图上才具备说服力。QGIS里加载天地图影像作为参考底图把黑臭分级栅格和水华分布栅格按半透明叠加检查异常斑块是否落在河道内部、桥梁阴影是否误判。出图时按专题图标准配置色带用绿到红的分级表示黑臭程度水华用红色半透明图层叠加图例里写明指数名称、阈值和影像获取日期指北针、比例尺、坐标参考系一个都不能少。比例尺用1:50000到1:100000比较合适太细会暴露16米像元的锯齿。我习惯在出图前把阈值参数、验证点表、影像日期一起存档进项目文件夹。这样两个月后回访同一个水库时能直接调出上一期的阈值和色带保证两期图件可比。反演这件事做到最后拼的不是算法多先进而是阈值标定和过程记录的细致程度。希望帮到你。本文还有配套的精品资源点击获取