ArcGIS夜间灯光数据校正实战:DMSP-OLS去饱和与NPP-VIIRS年际归一化
简介这份文档面向GIS从业者与遥感数据分析学习者系统讲解DMSP/OLS夜间灯光数据在ArcGIS Desktop中的校正流程解决传感器辐射性能差异、年际数据不连续、F18突变及DN值0-63天花板效应等导致数据不可比的问题。资源包共1个docx文件约592KB内容涵盖校正缘由剖析、中国区域亮值像元提取、兰伯特方位角等面积投影与NEAREST重采样设置以及基于伪不变区域与最小二乘回归的传感器依次校正方法并说明相邻传感器无重合年份时的处理思路。已有5973人学习下载适合需要掌握夜间灯光数据预处理、开展城市化与GDP匹配研究的中高级GIS用户参考可帮助读者理解校正模型构建逻辑与关键参数选择提升长时序灯光数据的连续性与可比性。1. 夜间灯光数据校正到底在修什么从一张“亮得离谱”的城区图说起如果你手上有 DMSP-OLS 或 NPP-VIIRS 这类夜间灯光影像直接拿原始 DN 值去做建成区提取或 GDP 空间化十有八九会翻车。我最早做某城市扩张分析时把两期灯光影像一叠加发现同一个城区边缘的亮度差了两倍多当时以为是城市真的变亮了后来才发现是传感器饱和和年际定标差异在作怪。夜间灯光数据校正本质上就是解决三件事饱和像元去饱和、年际数据可比化、以及跨传感器一致性。在 ArcGIS 里做这件事不需要写复杂脚本但每一步的栅格计算和掩膜逻辑必须搞清楚否则出来的结果只是“看起来像校正过”。这篇笔记面向的是已经会用 ArcGIS 基本栅格工具、但被灯光数据校正卡住的从业者从操作步骤到参数设置再到我踩过的坑全部拆开讲。热词里常出现的 arcgis 统计分析、arcgis 裁剪影像、arcgis 做坡度图这些操作在灯光校正里都会以变体形式出现但逻辑完全不同别混用。2. 校正前的数据准备与 ArcGIS 环境确认别让坐标系和像元对齐毁掉一切2.1 灯光数据校正需要哪几类输入以及为什么不能直接拿原始 DN 开算夜间灯光校正不是单一操作它依赖三类输入待校正的灯光栅格、参考掩膜或辅助数据、以及目标年份的定标参数。以 DMSP-OLS 稳定灯光产品为例常见做法是先用一份高分辨率建成区掩膜比如从 Landsat 提取的不透水面来界定“哪些像元属于真实城市灯光”再对掩膜内的像元做去饱和。NPP-VIIRS 则更依赖年度合成产品和杂散光校正标志位。在 ArcGIS 里这些数据必须满足两个硬条件所有栅格必须统一到同一投影坐标系且像元大小完全一致。我见过太多人直接把 WGS84 地理坐标的灯光图和投影坐标的掩膜丢进栅格计算器结果 ArcGIS 不报错但输出栅格偏移了几百个像元校正完全失效。常见做法是先用“投影栅格”工具把灯光数据转到与掩膜一致的投影如 Albers 等积投影再用“重采样”把像元对齐到 1km 或 500m。注意重采样方法选“双线性”还是“最近邻”取决于你的灯光数据是连续型还是离散型——DMSP-OLS 的 DN 值是整数建议用最近邻避免产生非整数 DN 导致后续阈值判断混乱。2.2 在 ArcGIS 里统一坐标系与像元对齐的具体命令和参数假设你手头有一份地理坐标的 NPP-VIIRS 月度合成栅格和一份投影坐标的行政区矢量。第一步不是裁剪而是投影。打开 ArcToolbox → Data Management Tools → Projections and Transformations → Raster → Project Raster。输入栅格选灯光图输出坐标系选与行政区一致的投影重采样技术选 NEAREST输出像元大小填 500单位与投影一致。这一步完成后用“栅格转点”或“识别”工具抽查几个已知城市中心的坐标确认没有整体偏移。第二步是像元对齐如果掩膜像元是 1000m而灯光是 500m不要直接重采样掩膜而是用“重采样”工具把灯光聚合到 1000m聚合方法选“平均值”或“最大值”——做去饱和时通常用最大值保留亮区峰值。第三步用“按掩膜提取”把灯光裁到研究区但注意这个工具默认输出会保留掩膜范围如果掩膜有孔洞灯光也会被挖掉所以掩膜最好先做“栅格转面→消除→面转栅格”清理一遍。这些步骤听起来繁琐但少一步后面栅格计算器里就会出现 NoData 蔓延你以为是校正公式错了其实是数据没对齐。2.3 用“栅格计算器”做初步统计先看清 DN 分布再动手在正式校正前我习惯先用栅格计算器算几个统计量最大值、最小值、平均值、以及大于某个阈值比如 DN50的像元数。ArcGIS 里没有直接的“栅格统计”按钮但可以用“分区统计”或“Zonal Statistics as Table”以研究区为分区统计灯光栅格的 MAX、MEAN、STD。更直接的办法是打开栅格属性 → 源 → 统计值但那只对全图有效。如果你要按行政区统计用 Zonal Statistics as Table输入分区数据选行政区输入值栅格选灯光统计类型勾选 ALL。输出的表里MAX 列能告诉你饱和像元大概在什么量级MEAN 列能看出整体亮度水平。这一步的意义在于校正不是盲目套公式而是根据你的数据实际分布决定阈值。比如 DMSP-OLS 的饱和阈值通常在 DN63但不同年份、不同传感器版本会有差异你得先看统计表再定。另外热词里常有人搜“arcgis统计分析”在灯光校正里最实用的就是分区统计和栅格直方图别去折腾复杂的空间自相关先把基础分布摸清。3. DMSP-OLS 去饱和校正从阈值掩膜到回归调整的完整 ArcGIS 操作链3.1 为什么 DMSP-OLS 必须做去饱和以及 ArcGIS 里怎么构建饱和掩膜DMSP-OLS 的 DN 值上限是 63城市核心区往往一大片全是 63这就是饱和。饱和像元不携带内部差异信息直接用来做城市内部结构分析会得到“铁板一块”的假象。去饱和的核心思路是用辅助数据如植被指数、不透水面比例建立饱和像元 DN 与真实亮度的回归关系再把 63 替换成预测值。在 ArcGIS 里第一步是生成饱和掩膜栅格计算器输入Con(light.tif 63, 1, 0)输出一个二值栅格1 代表饱和。注意有些版本的数据最大值是 63但经过重采样后可能出现 63.0 浮点所以条件写成light.tif 63更稳妥。生成掩膜后用“栅格转面”把饱和区转成矢量面再与不透水面数据做相交得到“饱和且不透水面比例高”的区域作为回归样本区。这一步的坑在于如果掩膜范围太大回归样本会混入非城市灯光如油气田火炬导致校正后农村也变亮。我一般会把饱和掩膜再与人口密度栅格做一次叠加只保留人口密度大于一定阈值的像元。3.2 用栅格计算器实现经典去饱和公式参数怎么设、NoData 怎么避常见的去饱和公式是线性回归DN_corrected a * NDVI b其中 a、b 由饱和区样本回归得到。在 ArcGIS 里你可以先用“采样”工具Spatial Analyst → Extraction → Sample在饱和区随机采样导出 NDVI 和 DN 值到表格然后在 Excel 或 Python 里做回归得到 a 和 b。回到栅格计算器输入# 去饱和校正仅对饱和像元应用回归非饱和像元保留原值 Con(light.tif 63, 63 (0.85 * Float(ndvi.tif) 2.3), light.tif)这里0.85和2.3是假设的回归系数实际要用你的采样结果替换。逻辑说明Con函数先判断饱和条件满足则用回归预测值替换不满足则保留原始 DN。注意Float()转换因为 NDVI 是浮点不转换会导致整数截断。参数说明回归系数 a 通常为负值NDVI 越高灯光越暗不在城市内部NDVI 低而灯光高所以 a 可能为负具体符号取决于你的样本。如果回归 R² 低于 0.5建议换辅助数据比如用不透水面比例代替 NDVI。另一个坑是 NoData如果 NDVI 在饱和区有 NoDataCon会输出 NoData导致校正后出现空洞。解决办法是在栅格计算器里加一层IsNull判断或者提前用“焦点统计”填充 NDVI 的小空洞。3.3 校正后验证用分区统计对比校正前后 DN 均值变化校正完不能直接出图得验证。我通常用 Zonal Statistics as Table 分别统计校正前后各行政区的灯光均值然后算变化率。如果某个区的均值变化超过 30%要么是回归系数不合理要么是该区饱和像元占比过高导致外推过度。另一个验证方法是看直方图校正后 DN 最大值应该超过 63但不应出现极端异常值比如 200否则说明回归斜率过大。在 ArcGIS 里右键栅格图层 → 属性 → 符号系统 → 拉伸观察直方图形态。如果校正后直方图在 63 处仍有尖峰说明饱和掩膜没覆盖全检查条件是否用了而不是。这一步的血泪经验是去饱和不是越亮越好而是让饱和区的相对差异恢复出来如果校正后城市核心区反而比边缘暗那肯定是回归系数符号反了。4. NPP-VIIRS 年际校正与跨传感器一致性ArcGIS 里的相对辐射定标操作4.1 NPP-VIIRS 为什么需要年际校正从杂散光到传感器衰减NPP-VIIRS 的夜间灯光产品比 DMSP-OLS 动态范围大得多但它有另一个问题年际之间的绝对辐射值不可直接比较。原因包括传感器衰减、杂散光校正版本更新、以及月度合成中云污染残留。常见做法是选取一个参考年份比如 2015 年把其他年份的灯光影像通过线性回归调整到参考年份的辐射尺度上。在 ArcGIS 里这步叫“相对辐射归一化”操作上就是栅格计算器加回归系数。但前提是你要有稳定不变的目标区作为回归样本比如远离城市的沙漠或深海区域——这些区域灯光应该接近 0如果某年份在这些区域出现高值说明杂散光没除干净。我一般会选研究区内的几个大型公园或水库作为“暗目标”统计其 DN 均值然后计算年份间的偏移量。4.2 用“栅格计算器”做线性拉伸增益和偏移量的确定方法假设你已经通过暗目标统计得到 2016 年相对于 2015 年的增益 gain1.12偏移 offset-0.35。在 ArcGIS 栅格计算器里输入# 年际相对辐射归一化将 2016 年拉伸到 2015 年尺度 Float(viirs_2016.tif) * 1.12 - 0.35逻辑说明先转浮点避免整数溢出再乘增益加偏移。参数说明gain 通常接近 1如果偏离超过 20%说明两年数据版本差异太大建议换年份或改用官方年度合成产品。offset 一般为负值因为暗目标在后期年份可能因杂散光校正残留而略高。注意这个公式是全局应用但城市核心区可能因饱和VIIRS 也有饱和只是阈值更高导致拉伸后过亮。解决办法是先用“按掩膜提取”把城市核心区单独处理或者用分段拉伸对 DN 小于某阈值的区域用一套系数大于阈值的用另一套。ArcGIS 里可以用Con嵌套实现但代码会很长建议用 Python 脚本批量跑。4.3 跨传感器一致性把 DMSP-OLS 和 NPP-VIIRS 放到同一尺度如果你要做长时间序列比如 1992-2020必然遇到 DMSP-OLS 和 NPP-VIIRS 衔接问题。常见做法是找重叠年份2012-2013建立两种数据在同区域的回归关系然后把 DMSP 调整到 VIIRS 尺度或反之。在 ArcGIS 里步骤是先分别提取重叠年份的灯光栅格用“采样”工具在建成区随机取点导出 DN 对做回归得到斜率 k 和截距 b。然后对 DMSP 全系列应用Float(dmsp.tif) * k b。这里最大的坑是DMSP 的饱和像元在回归中会拉低斜率所以回归前必须先去饱和或者只用非饱和像元做样本。另一个坑是空间分辨率差异DMSP 是 1kmVIIRS 是 500m直接回归会受尺度效应影响。我一般先把 VIIRS 聚合到 1km用“重采样”选“平均值”再采样回归。这样得到的系数更稳健。校正后用同一套行政区边界做分区统计看两种数据在重叠年份的均值是否接近如果差异仍大于 15%说明回归样本有偏需要增加暗目标样本。5. 避坑与排查夜间灯光校正中最容易翻车的 5 个操作5.1 现象栅格计算器输出全为 NoData原因像元对齐或掩膜范围不匹配解决先做投影和重采样这是最高频的翻车现场。你写了一个看起来完美的公式点运行结果输出一片空白。原因通常有两个一是输入栅格的像元大小或对齐方式不一致ArcGIS 在栅格计算器里会取交集只要有一个栅格在某像元是 NoData输出就是 NoData二是掩膜范围比灯光范围小裁剪后外围全是 NoData。解决办法在栅格计算器之前用“投影栅格”和“重采样”把所有输入统一到同一坐标系和像元大小并用“栅格转面”检查掩膜是否有孔洞。如果只是小范围 NoData可以用“焦点统计”的“均值”填充但注意这会平滑灯光只适合非饱和区。5.2 现象校正后城市核心区反而变暗原因回归系数符号错误或样本混入非灯光像元解决检查回归样本的 NDVI 与 DN 散点图去饱和回归中如果辅助数据是 NDVI城市核心区 NDVI 低DN 高回归斜率应为负。但如果你采样时混入了水体或云影NDVI 异常低而 DN 也低会导致斜率变正校正后核心区被压低。排查方法把采样点导出到 Excel画 NDVI 和 DN 的散点图看趋势是否合理。如果散点图一团乱说明样本区不纯需要重新定义饱和掩膜比如加入不透水面比例阈值。另一个可能是回归截距过大导致非饱和区也被误改但你的公式只对饱和像元生效所以问题一定在样本。5.3 现象年际校正后暗目标区域出现负值原因偏移量过大或数据版本不一致解决限制输出最小值并检查暗目标统计NPP-VIIRS 年际拉伸时如果 offset 设得太大比如 -2.0暗目标区域原本 DN 接近 0减去后变成负数。ArcGIS 栅格计算器不会自动截断负值输出栅格会出现负 DN后续做对数或比值运算直接报错。解决办法在公式外层加Con(result 0, 0, result)或者用“栅格计算器”的Con函数限制最小值。但更根本的是检查暗目标统计如果暗目标区域在参考年份的均值是 0.5而目标年份是 2.0offset 应该是 -1.5而不是 -2.0。另外如果两年数据版本不同比如 2015 是 v1 杂散光校正2016 是 v2暗目标差异会很大建议统一使用同一版本。5.4 现象跨传感器校正后 DMSP 和 VIIRS 在重叠年份均值差 30% 以上原因回归样本未去饱和或尺度未统一解决先聚合 VIIRS 到 1km再只用非饱和 DMSP 像元回归跨传感器回归最容易忽略的是 DMSP 饱和。如果你直接用全部 DMSP 像元包括 63和 VIIRS 做回归斜率会被饱和像元拉低导致校正后 DMSP 整体偏暗。正确做法先用Con(dmsp.tif 63, dmsp.tif, 0)把饱和像元设为 0 或 NoData只保留非饱和像元做采样。同时VIIRS 要先聚合到 1km否则 500m 的 VIIRS 和 1km 的 DMSP 在空间上不对应采样点会错位。聚合时用“重采样”选“平均值”不要用“最大值”因为最大值会放大城市核心导致回归斜率偏高。做完这两步重叠年份的均值差异通常能降到 10% 以内。5.5 现象校正后栅格文件巨大ArcGIS 卡死原因输出格式未压缩或浮点精度过高解决输出为 TIFF 并设置压缩或转成整数夜间灯光校正涉及多次栅格计算如果每一步都输出浮点栅格文件会迅速膨胀。我见过一个研究区校正后单幅栅格超过 2GBArcGIS 打开就崩。解决办法在“环境设置”里把输出栅格格式设为 TIFF并勾选“压缩”为 LZW。如果最终结果不需要浮点精度可以在最后一步用Int()转成整数但注意转整数前先乘以 100 保留两位小数否则去饱和的细微差异会被抹掉。另外中间步骤可以用“栅格转整型”临时降低精度但只建议在非饱和区使用。血泪经验是每做完一步就检查文件大小超过 500MB 就考虑压缩或分块处理。6. 进阶技巧用 ArcGIS 模型构建器把校正流程自动化以及一个验证校正效果的小方法如果你要处理多期灯光数据手动点栅格计算器会疯掉。我后来用 ArcGIS 的模型构建器ModelBuilder把整个流程串起来输入灯光栅格和掩膜自动投影、重采样、去饱和、年际拉伸、输出校正后栅格。关键节点是“栅格计算器”工具在模型里可以插入变量把回归系数作为参数传入这样换年份只需改参数不用重连。模型构建器里还有一个“前提条件”设置确保投影完成后再执行计算避免顺序错乱。具体操作在模型里拖入“投影栅格”和“重采样”用连接线串到“栅格计算器”右键计算器 → 创建变量 → 从参数 → 数据类型选“栅格”然后双击计算器在表达式里用%变量名%引用。这样你就能批量跑 10 年数据晚上挂着早上收结果。验证校正效果除了分区统计我还会用一个简单方法计算校正前后灯光栅格与不透水面比例的相关性。校正前由于饱和城市核心区相关性会被压低校正后如果去饱和合理相关性应该提升。在 ArcGIS 里用“波段集统计”或“多元分析”里的“波段集统计”工具输入校正前后栅格和不透水面栅格输出相关系数矩阵。如果校正后相关系数反而下降说明回归系数有问题或者辅助数据选错了。另一个技巧是看边缘梯度校正后城市边缘的 DN 过渡应该更平滑而不是突然从 63 跳到 0。你可以沿一条穿过城市的剖面线用“堆栈剖面”工具画校正前后的 DN 曲线如果校正后曲线在核心区有起伏而不是平顶说明去饱和起作用了。最后说个我自己的习惯每次校正完我都会把关键参数回归系数、增益、偏移、阈值记在一个文本文件里和输出栅格放一起。过半年回头看没有这些记录你根本不知道当时为什么设了 0.85 而不是 0.9。夜间灯光校正不是一劳永逸的事数据版本更新、研究区变化都会让旧参数失效但有了记录至少能快速定位问题。希望帮到你。本文还有配套的精品资源点击获取