简介本资源面向遥感与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个波段蓝B1450-520nm、绿B2520-590nm、红B3630-690nm、近红外B4770-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(dtyperasterio.float32, count1) 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(dtyperasterio.float32, count1) 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, threshold0.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(dtyperasterio.uint8, count1) 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(dtyperasterio.float32, count1) 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的区域说明水华加重。变化检测的结果比单时相更有说服力也更容易被采纳。我一般会选枯水期和丰水期各一期枯水期黑臭问题更突出丰水期水华更容易爆发。两期影像的时相尽量接近间隔不超过一个月否则地表覆盖变化会干扰结果。做这个方向最深的体会是指数计算本身不难难的是把预处理做扎实。定标系数、大气校正、水体掩膜这三步里任何一步出问题后面的反演结果都不可信。我现在的习惯是每做完一景先看清洁水体的反射率曲线正常了再往下走。这个习惯帮我省了很多返工的时间。希望帮到你。本文还有配套的精品资源点击获取