我要提问
ARTICLE DETAIL

资讯详情

前沿编程新知与开发实战干货的深度解读。

高分一号WFV水体反演全流程:从预处理到专题图输出

高分一号WFV水体反演全流程:从预处理到专题图输出 简介这是一份以高分一号影像为数据源、借助ArcGIS与ENVI完成重庆局部区域水体黑臭与水华反演的实操型技术文档适合遥感入门者、环境监测相关专业学生以及需要系统掌握ENVI预处理和ArcGIS制图流程的工程人员。文档完整覆盖了从陆地观测卫星平台检索并下载影像数据开始到辐射定标、FLAASH大气校正、RPC正射校正、掩膜裁剪、水体指数计算NDWI/MNDWI再到ArcGIS中符号化出图与面积统计的各个阶段。每个环节不仅写明操作工具和参数设置还特别提醒Data Ignore Value、光谱响应函数导入、Build Mask与Apply Mask等细节可有效降低新手复现时的踩坑概率。资源包内共1个docx说明文档格式便于随时查阅压缩包整体约13.95MB。目前已有557人学习浏览适合想完整走通黑臭水体与水华反演链路、并直接对照操作说明出图的读者使用。1. 黑臭水体与水华反演方案从一景高分一号影像到一张可交付的水体分布图做水体反演项目最容易踩的坑不是你不会波段运算而是你下载完影像、打开ENVI之后面对一整条预处理链路不知道从哪一步开始。很多人拿着高分一号WFV数据直接去算MNDWI结果算出来全是噪点——因为GF-1 WFV传感器根本没有短波红外波段教科书里那套MNDWI公式在这里压根跑不通。这套方案的核心思路是先用辐射定标、大气校正、正射校正把影像从DN值还原成地表反射率再用掩膜把研究区的水域裁剪出来最后用NDWI加黑臭指数和水华指数做分级出图。整个流程从数据下载到ArcGIS出图全部覆盖适合正在做水体遥感方向毕设或课设的学生也适合需要快速出水体监测专题图的从业者参考。整条链路跑通的关键在于每一处参数背后为什么这么设。2. 数据获取与预处理链路GF-1 WFV影像的下载、辐射定标与FLAASH大气校正2.1 数据下载通用卫星载荷检索的参数怎么设才能不白等高分一号影像的下载是在陆地观测卫星数据服务网站上完成的。项目实际操作里走的是“通用卫星载荷检索”入口这个入口比普通检索更灵活可以同时约束时间、位置、云量和传感器。时间范围一般选在研究区水体较平静且晴空少云的季节比如秋冬枯水期避免夏季多云雾干扰。云覆盖量直接关系到后续大气校正的效果我一般会在检索条件里把云量阈值卡在10%以下超过这个值FLAASH校正出来的反射率经常对不上实测值后面反演阶段还得返工。传感器这一栏选择高分一号的WFV载荷也就是四波段多光谱相机。这里有一个容易被忽略的细节GF-1搭载的WFV传感器有多个相机下发的数据文件命名里会带上WFV1到WFV4的标识不同相机的光谱响应曲线略有差异后期在大气校正里选择光谱曲线文件时必须和实际传感器匹配。检索出来后不要急着下载全部文件先看缩略图确认研究区覆盖完整、无明显条带或坏线再按景下载。数据下载完成后解压会得到一个XML元数据文件加多个TIFF波段文件ENVI直接打开XML即可自动识别波段组合。2.2 辐射定标把传感器DN值还原成物理辐射量辐射定标的目的是把传感器记录的原始DN值转换为大气顶层的辐亮度或反射率消除传感器响应差异和太阳光照条件带来的影响。这一步是后面FLAASH大气校正的输入前提跳过去直接做大气校正FLAASH会因为输入数据格式不对而报错或者校正结果失真。在ENVI工具栏搜索rad调出Radiometric Calibration工具选择刚打开的GF-1影像。参数设置上要注意几点定标类型选择辐亮度输出数据类型建议选Float避免Byte型截断导致暗部细节丢失。ENVI新版本已经内置了GF-1各传感器的定标系数工具会自动读取XML元数据里的增益和偏置值不需要手动填表前提是打开影像时正确选择了XML文件。定标完成后输出结果此时影像的像元值从0-255的整型变成了带物理单位的浮点辐亮度值你可以用Cursor Value查看数值范围正常情况下辐亮度值应该落在0到20之间的合理区间如果出现数量级异常说明定标类型选错了。2.3 大气校正FLAASH参数与gf1-wfv3光谱响应曲线的匹配大气校正的作用是消除大气分子和气溶胶对地表反射信号的散射、吸收影响把辐射定标得到的辐亮度进一步还原为地表真实反射率。水体目标的信号本身偏暗反射率低大气路径辐射占比高不做大气校正水体提取的阈值会整体偏移黑臭水体的微弱光谱差异很容易被淹没。FLAASH工具的启动方式是在ENVI工具栏搜索fla选择FLAASH Atmospheric Correction。输入辐射定标结果设置输出路径。参数设置是整个链路中最容易翻车的环节参考配置如下参数项设置值说明Sensor TypeUnknown未识别时选由输入文件属性决定Flight Date / Time从XML元数据读取必须精确到小时影响太阳角度计算Atmospheric ModelMid-Latitude Winter / Summer按影像拍摄季节和纬度选择Aerosol ModelRural / Urban城市研究区选UrbanAerosol Retrieval2-band (K-T)水体区域常用Initial Visibility40-60 km晴空影像设40以上Multispectral SettingsFilter Function File必须手动加载gf1-wfv3光谱曲线重点在Multispectral Settings这一步。点击Filter Function File后通过Open打开Spectral Library选择gf1-wfv3光谱响应曲线文件。这个文件对应的是GF-1 WFV3传感器的光谱响应函数。如果下载的影像来自WFV1或WFV2相机选择错了曲线大气校正后的反射率会出现偏绿或偏红的现象。某些ENVI版本里加载光谱曲线后还需要确认每个波段对应的中心波长是否从曲线文件中正确读取未识别时手动填写的中心波长和FWHM值一并写入面板。高级设置不需要动保持默认即可点击Apply开始计算一景影像大约运行十几分钟到半小时视计算机性能而定。大气校正完成后建议立即对比校正前后的植被和水体光谱曲线。水体在近红外波段的反射率应显著下降植被在近红外出现高反射这两个特征没出现说明校正参数有问题回到FLAASH面板检查传感器型号和大气模型。3. 正射校正与水域裁剪消除地形畸变并用掩膜提取研究区水体3.1 正射校正RPC模型、参考影像与高程数据如何组合大气校正做完后影像还存在几何畸变。GF-1 WFV的侧摆角加上地形起伏会导致影像上的地物位置与实际经纬度有几十到上百米的偏移这种偏移在后续叠加研究区水域shp做裁剪时会直接让水体边界错位。正射校正需要用到影像自带的RPC有理多项式系数文件ENVI的RPC Orthorectification Using Reference Image工具让这个过程变成半自动的。在ENVI工具栏搜索RPC调用正射校正工具。输入文件选择大气校正结果参考影像选择同一区域、空间分辨率更高且几何精度可靠的影像数据如果没有现成的参考影像可以用ENVI自带的全球底图或者公开的高分辨率影像服务。高程数据一般选择同区域SRTM或ASTER GDEM。参数上重采样方法选Cubic Convolution输出分辨率参考GF-1 WFV的16米原始分辨率不要设置过小计算量会增加很多但精度提升有限。正射校正的输出结果周围会出现一圈黑色背景区域这是因为校正过程中影像旋转和重投影产生的无数据区域。这圈黑边不处理后续掩膜裁剪会把背景值一起算进去波段计算时0值会严重干扰水体指数的统计范围所以必须继续做去黑边。3.2 去除黑背景值Data Ignore Value设置的正确时机去黑边的标准做法是给影像设置数据忽略值Data Ignore Value让ENVI在后续所有计算中把指定数值当作无效值处理而不是参与运算。注意这个步骤必须在正射校正输出之后、波段运算之前完成顺序反了的话噪声会一路带进反演结果。具体操作是双击正射校正结果打开图层在图层属性中找到Edit Metadata打开元数据编辑窗口后点击Add Metadata items在下拉列表里选择Data Ignore Value并添加。这个操作的关键在于设置值必须是0因为正射校正产生的背景区域像元值就是0。点击OK保存后影像显示时黑色背景会变成透明或统一的无数据颜色。这里有一个很隐蔽的坑如果你在设置Data Ignore Value之前先执行了掩膜裁剪背景0值会进入裁剪结果导致后期统计水体面积时把整幅图的无效像素也统计进去返工成本很高。3.3 水域裁剪Build Mask与Apply Mask的标准流程正射校正和去黑边完成后就需要把研究区的水体范围从整景影像中切出来。这一步的核心是使用前期准备好的研究区水域矢量边界shp文件。如果手里没有现成的shp通常的做法是在ArcGIS里基于同区域的水体分布数据或土地利用数据按行政区边界或流域边界手动裁剪出水域范围再另存为多边形shp。操作分两步走。第一步在ENVI工具栏搜索mask选择Build Mask工具。输入文件选择正射校正去黑边后的影像用矢量数据包含的方式加载研究区水域shp设置好输出路径和文件名生成掩膜文件mask。Build Mask会以shp多边形为界生成一个内部值为1、外部值为0的单波段掩膜。第二步再次搜索mask调用Apply Mask工具。输入影像选择正射校正结果掩膜文件选择刚才生成的mask输出得到水体裁剪结果。Apply Mask会把掩膜中值为0的区域全部置为背景值保留多边形内部的水体像元。输出后的影像里水域部分保留地表反射率值非水域区域全部变成0。建议执行完裁剪后在Display里叠加shp检查一次水体边界与shp边界应该重合且无明显错位如果边界偏移明显说明第3步的正射校正没做准需要返回去检查参考影像和RPC设置。此外裁剪得到的water.dat就是后续波段计算的标准输入文件。4. 反演建模与制图出图NDWI波段计算、分级阈值与ArcGIS符号化4.1 Band Math波段计算NDWI公式的编写与波段映射反演阶段的核心任务是把预处理好的水体反射率影像转换成能区分“正常水体—黑臭水体—水华”的指数图。针对GF-1 WFV的波段配置蓝、绿、红、近红外四个波段标准做法是使用NDWI归一化差异水体指数。GF-1 WFV没有短波红外波段导致常用的MNDWI无法计算这里只能用NDWI提取水域再结合其他指数识别污染类型。NDWI的公式是(Green - NIR) / (Green NIR)在ENVI中通过Band Math实现。在工具栏搜索band打开波段计算器输入表达式(float(b2) - float(b4)) / (float(b2) float(b4))其中b2对应绿光波段中心波长约555nmb4对应近红外波段中心波长约830nm。用float()强制转换是为了防止整型除法导致结果全为0。表达式输入后Band Math会弹出一个映射窗口点击Map Variable to Input File把b2和b4映射到水域裁剪结果water.dat对应的波段上。映射时必须确认波段顺序GF-1 WFV的波段顺序是标准的蓝、绿、红、近红外如果你的数据是Landsat或其他传感器波段顺序完全不同这也是跨项目复用公式时最常见的一个坑。除了NDWI黑臭水体的识别还会用到BOI黑臭水体指数公式为(Blue Green - Red) / (Blue Green Red)水华识别则常用FAI浮游藻类指数或者用NDWI提取水域后叠加近红外高反射特征来圈定水华密集区。我在实际项目中一般会一次性在Band Math里把NDWI和BOI都算出来输出为两个单独的单波段文件便于后面在ArcGIS里做符号化叠加。这里补一个用Python快速验证NDWI公式的思路适合需要批量处理或做算法对比的同学参考# 基于GDAL读取预处理后的GF-1影像计算NDWI import numpy as np from osgeo import gdal ds gdal.Open(water.dat) b2 ds.GetRasterBand(2).ReadAsArray().astype(np.float32) # 绿光 b4 ds.GetRasterBand(4).ReadAsArray().astype(np.float32) # 近红外 ndwi (b2 - b4) / (b2 b4 1e-10) # 加极小值避免除零 # 输出GeoTIFF driver gdal.GetDriverByName(GTiff) out driver.Create(ndwi.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) out.GetRasterBand(1).WriteArray(ndwi) out.SetGeoTransform(ds.GetGeoTransform()) out.SetProjection(ds.GetProjection()) out.FlushCache()ENVI的Band Math适合单次交互计算脚本适合批量处理多景影像原理完全一致。波段运算完成后检查结果影像的直方图正常水体的NDWI值应为正值集中在0到1之间非水体和背景值应为负值或0。如果整幅影像的NDWI值全部为负几乎可以肯定是b2和b4映射反了。4.2 黑臭水体与水华的指数分级与阈值设定NDWI提取出水体像元后接下来要区分正常水体、黑臭水体和水华区域。基于GF-1的波段范围常用的分级策略是组合NDWI与BOI、FAI两个指数设定阈值后进行密度分割。黑臭水体的光谱特征表现为绿光和红光波段的反射率较高蓝光波段受悬浮物和有机质影响反射率相对低因此BOI指数Blue Green - Red/Blue Green Red在黑臭区域会显著偏高。水华的识别则利用藻类叶绿素在近红外波段的强反射特征FAI指数通过红光、近红外和短波红外波段的差值计算但GF-1没有短波红外波段实战里用NDWI高值加近红外反射率突出来近似替代。项目中的分级阈值建议如下水体类型判别条件备注正常水体NDWI 0 且 BOI 0.05NDWI正值无明显黑臭特征疑似黑臭水体NDWI 0 且 BOI ≥ 0.05结合绿地河岸阴影做目视修正疑似水华NDWI 0 且 NIR反射率 0.08NIR突增对应藻类暴发在实际项目中阈值一般是通过统计研究区样本点的指数值分布来确定先抽取几十个已知水体状态的点位分别计算NDWI和BOI值再做散点图找分界线。如果手头没有实测数据先按上表经验值试算一版再叠加真彩色合成影像人工目视调整。比如某条河流岸边有阴影遮挡的区域BOI容易虚高这类误检测需要结合空间位置过滤掉。4.3 ArcGIS出图符号化分级、图例与专题地图输出反演结果最终要通过ArcGIS制成一幅可读的水体状态专题图。入图数据包括NDWI分级结果、黑臭水体分布图层和水华分布图层底图建议选研究区的行政区划或河流水系矢量。在ArcGIS中加载反演得到的tif或dat文件后打开图层属性调整符号化。水体分级图用分类色带比连续渐变更直观正常水体设为蓝色疑似黑臭设为棕黑色疑似水华设为黄绿色。具体操作是Layer Properties → Symbology → Classified将栅格按阈值范围分为三类手动修改每类的Label和颜色。注意在Symbology里设置NoData为透明显示否则黑色背景会盖住底图。制图要素不能缺比例为尺、指北针、图例、标题、坐标网格和影像日期。标题写明“研究区黑臭水体与水华遥感反演分布图”日期精确到影像获取日这在答辩时是数据时效性的证明。出图前用Export Map导出分辨率设150dpi以上格式选TIFF或PDF避免直接截图导致文字模糊。还可以在ArcGIS里做一步面积统计在Spatial Analyst工具里使用Zonal Geometry或Tabulate Area统计三类水体各自的像元数和面积这个数字写进论文的结论部分非常加分。5. 常见问题与避坑指南预处理链路上最容易翻车的5个环节5.1 大气校正后影像色彩异常水体发红或植被发黄现象是FLAASH校正完输出影像后查看水体反射率出现异常升高整个画面偏红或偏黄正常的地物光谱特征全部消失。这个问题的根源几乎都是光谱响应曲线和传感器型号不匹配。gf1-wfv3光谱曲线对应的是WFV3相机如果你的影像来自WFV1或WFV2校正时各波段的权重系数就加载错了大气校正等于做了个反向运算。解决方法是回到Multispectral Settings里核对影像的传感器标识。在ENVI打开XML元数据文件找到Sensor ID字段与光谱响应曲线文件的名称对照。网上有些教程习惯性地让人统一选gf1-wfv3这是不对的WFV1到WFV4四台相机的光谱响应曲线都有单独文件。不确定的情况下可以分别用两个相机曲线各跑一次对比水体在近红外波段的反射率正常结果该值应接近0。5.2 正射校正后影像全黑或边缘残留黑边正射校正完成后输出影像发现整个图面全黑只有边缘的黑色背景正常显示或者黑边无法去掉。全黑的情况多半是输出数据类型设置成了Byte型正射校正重采样过程中插值产生了负值被Byte型截断后全部变成0。边缘黑边去不掉则是没有正确设置Data Ignore Value后续裁切时的背景0值被当作有效像元保留。全黑问题好解决把输出数据类型改成Float再跑一次正射校正。黑边问题的关键是要同步设置Data Ignore Value为0然后在使用Raster Management里的Build Mask之前确保工具能识别到该忽略值。有个笨办法可以验证在波段计算器里用b1 gt 0做一次筛选输出结果中背景区域应该全部是0值如果还有非0像素说明忽略值没生效。5.3 掩膜裁剪后输出影像全白或全黑Build Mask和Apply Mask两个工具跑完后打开结果影像发现整幅图要么全黑要么全白。全黑表示掩膜的外部区域被全部赋0内部水体也一起被置0了全白则正好相反掩膜和水体位置错开裁切的是一块空白区域。原因是掩膜生成时矢量与栅格没有对齐shp文件与研究区影像的投影坐标系不一致或者Build Mask里选择的区域和输入影像的空间范围不匹配。解决顺序是先检查shp文件的坐标系与影像的不一致时在ArcGIS里用Project工具统一再回到ENVI重新生成掩膜。另一个常见原因是Apply Mask时选择掩膜文件和输入文件时顺序颠倒工具要求输入文件选待裁影像掩膜文件选mask反过来选必然全黑。操作前记得用cursor查看掩膜文件的像元值分布正常应该是内部1、外部0。5.4 Band Math计算结果全是0或数值溢出NDWI表达式计算完结果影像打开看着像噪声图像元值全部为0或者出现大量NaN。除零是最常见的原因当b2加b4恰好等于0时会产生NaN或无穷值。GF-1影像定标后的辐亮度很少出现这个情况但如果你输入的是整型数据且裁切后背景值不是忽略值而是真实的0像素那么背景区域的b2加b4就等于0除法直接异常。解决办法是在公式中加一个极小偏移量(float(b2)-float(b4))/(float(b2)float(b4)0.001)或者先通过Band Math计算一个水体掩膜把背景置为NaN后再做归一化除法。另外注意Band Math表达式里变量名不要用系统保留字符直接b1、b2这种最省事。5.5 成果图里反演结果与shp边界错位入图叠加时发现水体反演结果的边界和shp里的河网边界相差好几十米甚至上百米看起来像是对不齐。这通常是两类原因叠加在一起一是正射校正的参考影像精度不够造成几何偏差二是ArcGIS中底图和反演图层的投影坐标系统一工作没做到位。处理上先确认所有图层在ArcGIS中设为同一投影坐标系例如CGCS2000或WGS84再检查正射校正使用的参考影像分辨率如果参考影像是30米分辨率的公开底图16米分辨率GF-1的几何校正就能做到半个像元左右这已经不错了。真正要做高精度只能在正射校正阶段使用更高分辨率的DOM作为参考影像并添加地面控制点。项目精度要求不高的话用RPC自带模型配合SRTM就已经可以满足水体反演出图的要求。6. 成果验收与答辩用交叉验证和阈值迁移说服评审验收一套水体反演结果光靠“看着像”是不够的。我的建议是准备一套混交矩阵式的抽样验证方案在ArcGIS里创建50到100个随机验证点叠加到反演结果和同期真彩色影像上逐个点目视判定水体状态再与反演分类结果对比计算出总体精度和Kappa系数。验证点实测类型反演分类一致数反演分类不一致数用户精度正常水体38490.5%黑臭水体15383.3%水华区域9281.8%如果拿不到实测水质数据就改用目视解译替代实测把RGB合成影像放到最大比例尺人工圈出明显的黑臭河段和水华暴发区域再与反演结果做空间重合度分析重合率超过85%就可以认为反演结果可信。这个验证过程建议在答辩PPT里放一张图左边是原始影像右边是反演分类结果中间用箭头标注几个典型点位。答辩时有几个高频问题要提前准备。第一个问题是为什么用NDWI而不用MNDWI回答要点是GF-1 WFV没有短波红外波段MNDWI依赖的SWIR波段缺失NDWI是当前数据条件下的合理选择。第二个问题是阈值是怎么定的回答思路是先用经验初值试算再结合研究区典型水体的光谱曲线和目视解译做局部调整并说明阈值适用范围限于本项目影像时相。第三个问题是精度能达到多少用抽样验证的总体精度和Kappa系数说话。第四个问题是这种方法能不能批量推广可以回答同一流程适用于GF-1全系列WFV数据但需要重新匹配光谱曲线和调整阈值。答辩前我还会把ENVI的处理流程和ArcGIS出图的每一个关键步骤截图按顺序整理成一页附录方便老师侧目一页就看到完整链路。一个容易被忽略的习惯是把辐射定标、大气校正、正射校正每个环节输出的中间文件单独存档命名贴进论文作为数据预处理说明的支撑材料这部分内容往往是答辩中加分最明显的细节。从那以后我每次做水体反演项目都强制自己走一遍“先验证传感器型号—再跑预处理链路—最后做交叉验证”的流程中间任何一步发现异常先停下来定位问题而不是继续往下跑。这套方案如果能帮你少走一次弯路希望帮到你。本文还有配套的精品资源点击获取
返回列表