简介：本资源面向遥感与GIS方向的初学者及环境监测从业者，提供利用高分一号影像对重庆部分区域开展黑臭水体与水华反演的完整技术方案，涉及ArcGIS与ENVI两大工具，适合作为水体遥感提取的实操参考。压缩包内共1个docx文档，约13.95MB，以图文步骤形式记录从数据获取到成图的全流程。内容涵盖在陆地观测卫星数据服务网站检索并下载高分一号影像，在ENVI中依次完成辐射定标、FLAASH大气校正、基于RPC模型的正射校正与黑背景值去除，再通过Build Mask与Apply Mask裁剪水域，最后利用Band Math计算NDWI、MNDWI等指数提取水体，并在ArcGIS中符号化区分正常水体、黑臭水体与水华区，导出专题图。目前已有557人学习，读者可据此掌握一套可复用的遥感水体监测流程与排错思路。
1. 高分一号影像做黑臭水体与水华反演：从数据到出图，一条能跑通的路
黑臭水体和水华，一个是岸上排口惹的祸，一个是水里营养盐爆的表，看着是两件事，但在遥感影像上它们经常同时出现在一幅图里——城区河段发黑发臭，下游湖库绿得发亮。高分一号（GF-1）这类国产宽幅多光谱数据，16米分辨率、4天重访，做城市尺度到流域尺度的水环境反演，性价比很高。这篇笔记讲的就是一条完整链路：从拿到GF-1的L1A级原始影像开始，做辐射定标、大气校正、水陆分离、黑臭水体指数计算、水华指数计算，最后出专题图。适合做环保遥感、水务监测、湖泊富营养化评估的从业者，也适合想拿国产卫星数据练手的学生。整条链路我按实际项目里跑通的顺序写，参数给到能直接抄的程度。
2. 数据准备与预处理：GF-1的PMS数据怎么变成能算的反射率
2.1 为什么必须做辐射定标和大气校正
GF-1的PMS传感器有4个波段：蓝（B1，450-520nm）、绿（B2，520-590nm）、红（B3，630-690nm）、近红外（B4，770-890nm）。L1A级产品是DN值，直接拿来算指数，不同时相之间没有可比性，同一时相不同景之间也没法拼。黑臭水体反演对蓝绿波段的差异很敏感，水华反演依赖近红外和红波段的比值，这两个应用都要求输入是地表反射率。
常见做法是：辐射定标用定标系数把DN转成表观辐亮度，再做FLAASH或6S大气校正得到地表反射率。如果只是做同一景内的相对比较，跳过大气校正用表观反射率也能出图，但跨时相对比就必须做。我一般会做完整的大气校正，因为黑臭水体的信号本来就弱，大气程辐射的干扰不能忽略。
2.2 用Python批量做辐射定标
GF-1的定标系数在元数据XML里，不同年份发射的星（GF-1有01、02、03、04星）系数不同，必须从元数据读，不能写死。
import xml.etree.ElementTree as ET import numpy as np import rasterio def radiometric_calibration(dn_path, xml_path, out_path): """ GF-1 PMS 辐射定标：DN -> 表观辐亮度 dn_path: 单波段DN影像路径 xml_path: 元数据XML路径 out_path: 输出辐亮度影像路径 """ # 从XML读取定标系数和太阳高度角 tree = ET.parse(xml_path) root = tree.getroot() # 不同星和不同年份的XML结构略有差异，这里取通用字段 gain = float(root.find('.//Gain').text) # 增益 offset = float(root.find('.//Offset').text) # 偏移 sun_elev = float(root.find('.//SunElevation').text) # 太阳高度角，度 with rasterio.open(dn_path) as src: dn = src.read(1).astype(np.float32) profile = src.profile # 辐亮度 = 增益 * DN + 偏移 radiance = gain * dn + offset # 转表观反射率：ρ = π * L * d^2 / (ESUN * cos(θ)) # d为日地距离（天文单位），GF-1常用近似值1.0 # ESUN为波段太阳辐照度，需按波段查表 esun = {'B1': 1958.0, 'B2': 1852.0, 'B3': 1562.0, 'B4': 1093.0} band = 'B3' # 示例用红波段 d = 1.0 theta = np.deg2rad(90.0 - sun_elev) reflectance = np.pi * radiance * d**2 / (esun[band] * np.cos(theta)) profile.update(dtype=rasterio.float32, count=1) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(reflectance, 1) return out_path这段代码的关键参数有三个：Gain和Offset必须从对应影像的XML里读，不能跨星混用；SunElevation影响反射率量级，填错会让整景偏亮或偏暗；ESUN是波段太阳辐照度，不同波段值不同，用错波段会导致反射率量级错误。实际项目中，我建议把4个波段循环处理，输出4个反射率文件，后续指数计算直接调用。
2.3 大气校正的替代方案与验证
如果手头没有FLAASH或6S的授权，可以用暗像元法做简化大气校正：在影像里找一块清洁水体（近红外反射率接近0），把它的反射率当作大气程辐射的估计值，从所有波段减去。这个方法在城区影像上效果一般，因为暗像元不好找，但在湖库场景够用。
验证校正结果是否合理，看两个地方：清洁水体的近红外反射率应该在0.02以下，植被的近红外反射率应该在0.3以上。如果清洁水体近红外超过0.05，说明大气校正不够或者定标系数有问题。
3. 黑臭水体反演：用指数法把发黑发臭的河段圈出来
3.1 黑臭水体的光谱特征与指数选择
黑臭水体的核心问题是：水体里有机污染物和悬浮物浓度高，导致蓝绿波段反射率升高，红波段和近红外被强烈吸收。表现到光谱上，就是蓝绿波段反射率相对清洁水体明显抬升，红波段到近红外急剧下降，整个可见光区间的反射率曲线比清洁水体更“平”。
常用的指数有两个：一是黑臭水体指数（BOI），用绿波段和红波段的差值除以绿红之和，再减去近红外；二是归一化黑臭指数（NDBWI），用蓝波段和绿波段的组合。我一般用BOI，因为它对悬浮物和有机物的响应更直接。
3.2 BOI计算与阈值分割
import numpy as np import rasterio def calculate_boi(green_path, red_path, nir_path, out_path): """ 黑臭水体指数 BOI BOI = (Green - Red) / (Green + Red) - NIR 值越高，黑臭可能性越大 """ with rasterio.open(green_path) as src: green = src.read(1).astype(np.float32) profile = src.profile with rasterio.open(red_path) as src: red = src.read(1).astype(np.float32) with rasterio.open(nir_path) as src: nir = src.read(1).astype(np.float32) # 避免除零 denom = green + red denom[denom == 0] = 1e-6 boi = (green - red) / denom - nir profile.update(dtype=rasterio.float32, count=1) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(boi, 1) return out_pathBOI的阈值需要按区域调。在城区河段，我一般先用0.05作为初筛，把BOI大于0.05的像元当作疑似黑臭，然后叠加水体掩膜（用NDWI提取的水体范围）做交集。如果初筛结果太多，把阈值提到0.08；如果漏掉了明显黑臭的河段，降到0.03。这个阈值没有全国统一值，必须用地面实测数据或者目视解译做校验。
3.3 水陆分离：先用水体指数把非水体剔掉
黑臭水体反演只在水体范围内做，所以要先提取水体。用NDWI（归一化水体指数）：
def calculate_ndwi(green_path, nir_path, out_path, threshold=0.1): """ NDWI = (Green - NIR) / (Green + NIR) 大于阈值判为水体 """ with rasterio.open(green_path) as src: green = src.read(1).astype(np.float32) profile = src.profile with rasterio.open(nir_path) as src: nir = src.read(1).astype(np.float32) denom = green + nir denom[denom == 0] = 1e-6 ndwi = (green - nir) / denom water_mask = (ndwi > threshold).astype(np.uint8) profile.update(dtype=rasterio.uint8, count=1) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(water_mask, 1) return out_pathNDWI阈值取0.1是经验值，城区水体因为悬浮物多，NDWI可能偏低，可以降到0.05。提取完水体后，用栅格计算器把BOI和水体掩膜相乘，非水体区域置为0，这样后续统计和出图都干净。
4. 水华反演：用NDCI和FAI区分藻类浓度
4.1 水华的光谱特征与指数选择
水华（蓝藻爆发）的光谱特征和黑臭水体正好相反：藻类细胞里的叶绿素在近红外波段有强反射，在红波段有吸收，导致红波段反射率下降、近红外反射率升高。表现到影像上，就是水体在近红外波段变亮，和清洁水体在近红外波段的暗信号形成鲜明对比。
常用指数有两个：NDCI（归一化叶绿素指数），用红波段和红边波段；FAI（浮游藻类指数），用近红外、红、短波红外三个波段。GF-1没有红边波段，所以NDCI只能用红波段和近红外近似，精度会打折扣。我一般用FAI的简化版：用近红外减去红波段和短波红外的线性插值。但GF-1的PMS没有短波红外，所以只能用NDCI的变体：NDCI = (NIR - Red) / (NIR + Red)。
4.2 NDCI计算与藻类浓度分级
def calculate_ndci(red_path, nir_path, out_path): """ 归一化叶绿素指数 NDCI NDCI = (NIR - Red) / (NIR + Red) 值越高，藻类浓度越大 """ with rasterio.open(red_path) as src: red = src.read(1).astype(np.float32) profile = src.profile with rasterio.open(nir_path) as src: nir = src.read(1).astype(np.float32) denom = nir + red denom[denom == 0] = 1e-6 ndci = (nir - red) / denom profile.update(dtype=rasterio.float32, count=1) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(ndci, 1) return out_pathNDCI的分级阈值：小于0.1为清洁水体，0.1到0.2为轻度水华，0.2到0.3为中度水华，大于0.3为重度水华。这个分级在不同湖区会有偏移，比如太湖的阈值可能比滇池低0.05左右。实际项目中，我会用同期地面采样点的叶绿素a浓度做回归，把NDCI转成叶绿素a浓度，这样出图更有说服力。
4.3 黑臭与水华的叠加分析
同一景影像里，黑臭水体和水华可能同时存在。黑臭水体在BOI上表现为高值，水华在NDCI上表现为高值，两者在空间上一般不重叠——黑臭水体多在排污口附近，水华多在湖库开阔水域。但下游河段如果流速慢、营养盐高，可能出现黑臭和水华交替出现的复合污染。
叠加分析的步骤：先把BOI大于阈值的像元标记为黑臭，再把NDCI大于阈值的像元标记为水华，最后用分类图输出。如果某个像元同时满足两个条件，说明该区域既有黑臭又有藻类，需要单独标注。这种复合污染区域在实际项目里不多，但一旦出现，往往是治理的重点。
5. 避坑与排查：GF-1反演里最容易翻车的五个地方
5.1 定标系数用错星号，整景反射率偏亮或偏暗
现象：同一区域不同时相的影像，清洁水体的反射率差异超过0.05，或者整景看起来比实际亮很多。
原因：GF-1有01、02、03、04星，每颗星的定标系数不同，XML里的字段名也可能不一样。如果批量处理时用了同一套系数，或者读错了字段，反射率就会系统性偏移。
解决：处理前先确认影像的卫星编号和传感器编号，从对应XML里读系数。批量处理时，每景影像单独读自己的XML，不要用全局变量。验证方法是看清洁水体的近红外反射率，正常应该在0.02以下。
5.2 大气校正后水体反射率出现负值
现象：大气校正后，部分水体像元的反射率小于0，尤其是蓝波段。
原因：大气校正模型对气溶胶光学厚度的估计偏大，或者暗像元选取不当，导致过度校正。城区影像因为气溶胶复杂，这个问题更常见。
解决：对负值像元做截断，置为0。如果负值比例超过5%，说明大气校正参数有问题，需要重新调整气溶胶光学厚度或者换用暗像元法。我一般会在校正后加一步检查：统计各波段的负值比例，超过3%就回退。
5.3 NDWI阈值取0.1，城区水体漏提
现象：城区河段在NDWI图上信号很弱，提取的水体范围比实际窄，黑臭河段被漏掉。
原因：城区水体悬浮物浓度高，绿波段和近红外的差异被压缩，NDWI值偏低。用0.1的阈值会漏掉这些水体。
解决：城区场景把NDWI阈值降到0.05，或者改用MNDWI（用绿波段和中红外），但GF-1没有中红外，所以只能用NDWI加阈值调整。另一个办法是叠加建筑阴影掩膜，把阴影误提的水体剔掉。
5.4 BOI阈值一刀切，不同河段结果不一致
现象：上游清洁河段被误判为黑臭，下游重污染河段反而漏判。
原因：BOI是相对指数，受悬浮物本底值影响。上游水体悬浮物低，BOI本底也低，用同一个阈值会误判；下游悬浮物高，BOI本底高，阈值需要相应提高。
解决：分河段设定阈值，或者用地面实测数据做回归，把BOI转成黑臭等级。如果没有实测数据，至少按上游、中游、下游分三段，每段单独定阈值。我一般会先做目视解译，选几个典型河段，看BOI的分布范围，再定阈值。
5.5 水华和云影混淆，NDCI高值区出现在云边缘
现象：NDCI图上，云边缘出现高值区，被误判为水华。
原因：云的近红外反射率高，云边缘的混合像元在近红外和红波段的比例关系和水华类似，导致NDCI偏高。
解决：加云掩膜。GF-1没有热红外波段，不能用亮温做云检测，只能用蓝波段和近红外的组合做简单云检测：蓝波段反射率大于0.2且近红外大于0.3的像元判为云。云掩膜做完后，再算NDCI，云边缘的误判就没了。
6. 出图与验证：把反演结果变成能交差的专题图
6.1 用QGIS做分级设色和出图
反演结果最终要出图，我一般用QGIS做分级设色。黑臭水体按BOI值分三级：轻度（0.05-0.08）、中度（0.08-0.12）、重度（大于0.12），分别用黄、橙、红。水华按NDCI值分三级：轻度（0.1-0.2）、中度（0.2-0.3）、重度（大于0.3），分别用浅绿、中绿、深绿。
出图时注意两点：一是底图用真彩色影像，方便对照；二是加比例尺和指北针，专题图没有这两个元素，交差时会被打回。QGIS的打印布局里可以直接加，不用额外插件。
6.2 用地面实测数据做精度验证
如果有地面实测数据，做混淆矩阵。黑臭水体的验证指标是总体精度和Kappa系数，一般要求总体精度大于85%，Kappa大于0.7。水华的验证指标是叶绿素a浓度的回归R²，要求大于0.6。如果达不到，先检查大气校正和阈值，再检查实测数据和影像的时相是否匹配——时相差超过3天，验证结果参考价值就有限。
6.3 一个我踩过的坑：出图时忘了统一坐标系
GF-1的L1A产品是WGS84地理坐标系，但有些项目要求投影到UTM或者高斯克吕格。如果出图时忘了投影转换，专题图的面积统计会偏，尤其是跨带区域。我现在的习惯是：预处理第一步就做投影转换，统一到项目要求的坐标系，后续所有计算都在投影坐标系下做。这样出图时不用再转，面积统计也准。
6.4 进阶技巧：用多时相影像做变化检测
单时相的反演只能看一个时刻的状态，实际项目里更关注变化。用两期GF-1影像，分别算BOI和NDCI，然后做差值：BOI差值大于0.05的区域，说明黑臭加重；NDCI差值大于0.1的区域，说明水华加重。变化检测的结果比单时相更有说服力，也更容易被采纳。
我一般会选枯水期和丰水期各一期，枯水期黑臭问题更突出，丰水期水华更容易爆发。两期影像的时相尽量接近，间隔不超过一个月，否则地表覆盖变化会干扰结果。
做这个方向，最深的体会是：指数计算本身不难，难的是把预处理做扎实。定标系数、大气校正、水体掩膜，这三步里任何一步出问题，后面的反演结果都不可信。我现在的习惯是，每做完一景，先看清洁水体的反射率曲线，正常了再往下走。这个习惯帮我省了很多返工的时间。希望帮到你。
本文还有配套的精品资源，点击获取
