简介这份文档面向GIS从业者与遥感数据分析学习者系统讲解DMSP/OLS夜间灯光数据在ArcGIS Desktop中的校正操作流程帮助解决传感器辐射性能差异、年际数据不连续、F18突变及DN值0-63天花板效应等实际问题。资源包内含1个docx文件约592KB以图文步骤形式呈现便于对照软件界面逐步操作。内容涵盖中国区域亮值像元影像提取、兰伯特方位角等面积投影转换、NEAREST重采样以及基于伪不变区域与最小二乘回归的传感器依次校正方法并说明无重合年份时如何借助相邻年份数据建立校正方程。目前已有5973人学习下载适合需要提升夜间灯光数据连续性与可比性、开展城市扩展监测或区域经济差异评估的研究人员参考可据此掌握从数据预处理到传感器相互校正的完整技术路线。1. 夜间灯光数据校正在ArcGIS里的真实门槛为什么你导出的DN值总对不上夜间灯光遥感数据比如常见的NTL产品拿到手很多人第一反应是直接拖进ArcGIS出图结果发现两个区域亮度差了三倍或者同一景影像在不同月份拼接后出现明显色阶断层。这不是软件的问题而是原始DN值本身带着传感器增益、大气散射、月相周期和饱和像元四层干扰。不做校正就直接做统计分析出来的GDP相关性、城市扩张指数基本是玄学。这篇笔记面向的是已经会用ArcGIS做基础裁剪和投影转换、但一碰到夜间灯光校正就卡住的从业者。我会把整个流程拆成可复现的步骤从DN值转辐射亮度、去饱和、年内多期合成到最终在ArcGIS里用栅格计算器落地。中间会给出具体参数、Python脚本和踩坑记录。如果你手头有DMSP-OLS或NPP-VIIRS类数据这套流程可以直接套用。ArcGIS 10.8和ArcGIS Pro 3.x在栅格计算器语法上略有差异我会分别标注。2. 校正前的数据准备投影、重采样与无效值处理2.1 为什么第一步不是打开栅格计算器夜间灯光原始数据通常以地理坐标系WGS84分发像元大小是30弧秒或15弧秒。直接做栅格计算ArcGIS会按地理坐标的度数计算导致高纬度区域像元面积严重变形。常见做法是先投影到适合研究区的等面积投影比如Albers或Lambert。我一般会先确认三件事数据是否已经过几何校正、无效值NoData填充的是什么、以及是否需要重采样到统一分辨率。以某跨平台系统的夜间灯光处理Demo为例原始数据是GeoTIFFNoData值为-9999。如果直接参与计算-9999会被当成真实DN值结果全错。所以第一步是用“设为空函数”或栅格计算器把无效值剔除。# ArcGIS Pro Python窗口将NoData值替换为真正的NoData import arcpy from arcpy.sa import * arcpy.env.workspace rC:\NTL_Project\raw arcpy.env.overwriteOutput True # 输入原始栅格 in_raster NTL_2020.tif # 用Con函数把-9999设为NoData其余保留原值 out_raster Con(Raster(in_raster) -9999, , Raster(in_raster)) out_raster.save(rC:\NTL_Project\processed\NTL_2020_clean.tif)这段代码的逻辑是Con函数逐像元判断如果值等于-9999就输出空否则输出原值。参数上注意表示NoData不要写成0否则后续统计会把0当成有效暗背景。ArcGIS 10.8里对应的是Spatial Analyst工具箱下的“条件函数”操作路径是Spatial Analyst 条件分析 条件函数表达式写NTL_2020.tif -9999真值为空假值为原栅格。2.2 投影转换与重采样的参数怎么设投影转换用“投影栅格”工具输出坐标系选Albers中央经线按研究区定。重采样方法选“双线性”还是“最近邻”夜间灯光是连续型栅格双线性更平滑但会改变DN值分布最近邻保留原始值但可能产生锯齿。我的经验是如果后续要做辐射定标用最近邻如果只是做可视化或趋势面分析双线性可以接受。重采样分辨率建议统一到500米或1公里。DMSP-OLS原始分辨率约1公里NPP-VIIRS约500米。如果混用两种数据必须重采样到同一网格。这里有个细节重采样时输出像元大小要写投影后的单位米不要写度数。# 投影并重采样到1公里Albers arcpy.ProjectRaster_management( in_rasterrC:\NTL_Project\processed\NTL_2020_clean.tif, out_rasterrC:\NTL_Project\processed\NTL_2020_albers.tif, out_coor_systemarcpy.SpatialReference(102025), # Albers for China resampling_typeNEAREST, cell_size1000 1000 )102025是某区域Albers投影的WKID实际用时换成你研究区对应的。cell_size写成1000 1000表示X和Y方向都是1000米。如果只写一个1000ArcGIS会自动应用为正方形像元。2.3 裁剪与掩膜别让背景值污染统计裁剪用“按掩膜提取”掩膜可以是研究区矢量边界。注意裁剪后边缘像元可能被重采样建议先裁剪再投影或者投影后裁剪时勾选“保持裁剪范围”。如果研究区跨多景影像先做镶嵌再裁剪镶嵌时重叠区域选“最大值”还是“平均值”夜间灯光重叠区通常取最大值因为灯光不会因为多景平均而变暗。# 按研究区边界裁剪 out_extract ExtractByMask( in_rasterrC:\NTL_Project\processed\NTL_2020_albers.tif, in_maskrC:\NTL_Project\boundary.shp ) out_extract.save(rC:\NTL_Project\processed\NTL_2020_clip.tif)到这里数据已经干净、投影统一、无效值处理完毕。接下来才是真正的校正环节。3. DN值转辐射亮度公式、参数与ArcGIS栅格计算器写法3.1 为什么不能直接用DN值做跨年比较DMSP-OLS的DN值是0-63的整数NPP-VIIRS是0-255左右的浮点。不同传感器、不同年份的增益设置不同直接比较DN值等于拿不同尺子的刻度对比。校正的第一步是转成物理量辐射亮度radiance或亮度温度。DMSP-OLS常用公式是L DN^(3/2) * 10^-6某版本定标公式NPP-VIIRS则用L DN * 10^-9具体系数看元数据。这里要强调不要背公式去查你下载数据时附带的元数据文件。元数据里会有radiance_calibration或gain字段。我见过有人把DMSP的公式套到VIIRS上结果整幅图亮度差了三个数量级。3.2 栅格计算器里的幂运算与浮点精度ArcGIS栅格计算器的幂运算用**不是^。^在Python里是异或在栅格计算器里可能被解释为其他含义。写公式时注意浮点精度DN是整数先转浮点再运算否则DN^(3/2)在整数运算下会截断。# ArcGIS Pro栅格计算器DMSP-OLS DN转辐射亮度 # 假设DN范围0-63公式 L DN^(3/2) * 10^-6 out_radiance (Float(NTL_2020_clip.tif) ** 1.5) * 0.000001 out_radiance.save(rC:\NTL_Project\processed\NTL_2020_radiance.tif)Float()函数把栅格转为浮点型避免整数幂运算截断。** 1.5就是DN的3/2次方。* 0.000001是乘以10的负6次方。如果你在ArcGIS 10.8里操作栅格计算器界面直接输入Float(NTL_2020_clip.tif) ** 1.5 * 0.000001注意文件名要带引号。对于NPP-VIIRS公式通常是L DN * 10^-9但有些产品已经提供了辐射亮度波段不需要再转。判断方法看元数据里units字段如果是nW/cm2/sr说明已经是辐射亮度如果是DN才需要转。3.3 去饱和DMSP-OLS的63阈值怎么破DMSP-OLS最头疼的是饱和城市中心DN值卡在63导致亮度被低估。常见做法是用NPP-VIIRS数据做参考对DMSP饱和像元进行替换或拟合。ArcGIS里可以用“栅格计算器”配合“Con”函数实现如果DMSP DN等于63就用VIIRS对应像元的辐射亮度按比例替换。# 去饱和用VIIRS辐射亮度替换DMSP饱和像元 # 先重采样VIIRS到与DMSP同一网格 viirs_resampled VIIRS_2020_radiance_resampled.tif dmsp_radiance NTL_2020_radiance.tif # 计算替换值VIIRS辐射亮度乘以一个经验系数需根据研究区拟合 # 这里假设系数为1.2实际应用时用回归分析确定 out_desaturated Con(Raster(dmsp_radiance) 63 * 0.000001, Raster(viirs_resampled) * 1.2, Raster(dmsp_radiance)) out_desaturated.save(rC:\NTL_Project\processed\NTL_2020_desaturated.tif)Con函数的第一个参数是条件DMSP辐射亮度是否达到饱和阈值63对应的辐射亮度值。第二个参数是真值用VIIRS辐射亮度乘以系数。第三个参数是假值保留原DMSP辐射亮度。系数1.2不是固定的需要用你研究区内未饱和像元做回归得到DMSP和VIIRS的线性关系再取斜率。注意去饱和只对城市中心有效如果研究区没有VIIRS数据可以用“饱和像元邻域均值”替代但效果差很多。4. 年内多期合成与跨年校正把12个月压成一张可比较的图4.1 月度数据合成的三种策略夜间灯光月度数据受月相、云层、气溶胶影响单月影像噪声大。常见合成策略有三种最大值合成MVC、平均值合成、中值合成。MVC保留最亮像元适合城市范围提取平均值合成平滑噪声适合趋势分析中值合成抗异常值适合长时间序列。我一般用MVC做城市扩张用中值合成做GDP相关性。ArcGIS里用“像元统计”工具统计类型选MAXIMUM或MEDIAN输入12个月栅格。# 月度数据最大值合成 arcpy.gp.CellStatistics_sa( in_rasters[NTL_2020_01.tif, NTL_2020_02.tif, ..., NTL_2020_12.tif], out_rasterrC:\NTL_Project\processed\NTL_2020_MVC.tif, statistics_typeMAXIMUM, ignore_nodataDATA )ignore_nodataDATA表示如果某月是NoData其他月参与统计。如果选NODATA则任一月为NoData则输出NoData会丢失大量像元。4.2 跨年校正用不变目标区域做相对辐射归一化不同年份的传感器增益不同即使都转了辐射亮度跨年比较仍有系统偏差。常用方法是选取“不变目标区域”如稳定城市中心或沙漠暗背景建立年份间的线性回归模型然后对整幅影像做归一化。步骤先在ArcGIS里用“创建随机点”在不变区域生成样本点再用“提取多值至点”获取各年份辐射亮度导出到Excel做回归得到斜率和截距最后用栅格计算器应用。# 假设回归得到2020年相对于2015年的校正L_2015 a * L_2020 b # a0.85, b0.02示例值实际用回归结果 out_corrected Raster(NTL_2020_radiance.tif) * 0.85 0.02 out_corrected.save(rC:\NTL_Project\processed\NTL_2020_corrected.tif)这里a和b必须来自你的回归分析不要用示例值。回归时注意剔除饱和像元和NoData否则斜率会被拉偏。4.3 用ArcGIS动态表格模块做校正质量检查校正后怎么验证我习惯用“动态表格模块”或“波段集统计”对比校正前后均值、标准差和直方图。如果校正后均值偏移超过10%说明回归模型有问题。另一个技巧在不变区域上计算校正前后的差值理想情况下差值应接近0。# 计算校正前后在不变区域的差值 diff Raster(NTL_2020_corrected.tif) - Raster(NTL_2020_radiance.tif) diff.save(rC:\NTL_Project\processed\diff_check.tif) # 然后用分区统计获取不变区域的均值 arcpy.gp.ZonalStatisticsAsTable_sa( in_zone_datainvariant_region.shp, zone_fieldID, in_value_rasterrC:\NTL_Project\processed\diff_check.tif, out_tablerC:\NTL_Project\processed\diff_stats.dbf, statistics_typeMEAN )如果MEAN绝对值大于0.05辐射亮度单位说明校正系数需要重新拟合。5. 避坑与排查夜间灯光校正里最容易翻车的5个地方5.1 现象栅格计算器报错“无法打开栅格”原因文件路径含中文或空格或者栅格被其他程序占用。ArcGIS对中文路径支持不稳定尤其是10.8版本。解决把所有数据放在纯英文路径下关闭其他打开该栅格的窗口。如果还是报错用“复制栅格”工具先转成Esri Grid格式再计算。5.2 现象校正后影像出现大面积0值原因NoData被当成0参与运算或者Con函数的假值写成了0。解决检查原始NoData值用SetNull或Con显式处理。另外栅格计算器里Float()转换时如果原栅格有NoData转换后仍是NoData不会变0。5.3 现象跨年校正后城市中心反而变暗原因回归斜率a小于1且截距b太小导致高值被压缩。解决检查回归样本是否包含饱和像元。如果包含剔除后重新拟合。另外如果研究区城市扩张明显不变区域选取要避开新城区。5.4 现象MVC合成后边缘出现条带原因不同月份影像的覆盖范围不一致边缘像元只有部分月份有值。解决在合成前用“镶嵌至新栅格”统一范围或者合成时选ignore_nodataDATA但边缘仍可能因月份数不足而偏低。建议至少保留6个月以上的有效像元。5.5 现象ArcGIS Pro 3.x里栅格计算器语法不兼容原因Pro 3.x默认使用Python 3Float()函数名没变但Con函数的参数顺序有调整。解决在Pro里用arcpy.sa.Con参数顺序是Con(in_conditional_raster, in_true_raster_or_constant, in_false_raster_or_constant)。如果从10.8迁移把旧表达式里的Con(条件, 真值, 假值)直接搬过来通常没问题但注意Raster()对象要显式声明。6. 进阶技巧用分区统计和动态表格模块做校正后验证校正做完不是终点验证才是。我一般会做两件事一是用“分区统计”计算校正前后各行政区的灯光总量看排名是否合理二是用“动态表格模块”生成时间序列曲线检查年际变化是否平滑。分区统计的代码前面已经给过这里补充一个技巧统计时用ALL会输出均值、标准差、最大值、最小值等但夜间灯光更关注SUM和MEAN。如果某区域SUM校正后反而下降说明该区域有大量饱和像元被错误替换。动态表格模块在ArcGIS Pro里叫“图表”功能可以绑定栅格的时间序列。操作路径右键栅格图层 创建图表 时间序列。如果数据没有时间字段先用“添加时间字段”工具把文件名里的年份提取出来。# 为多期栅格添加时间字段以文件名年份为例 arcpy.AddField_management(NTL_2020_corrected.tif, Year, LONG) arcpy.CalculateField_management(NTL_2020_corrected.tif, Year, 2020, PYTHON3)最后说一个我自己的习惯每次校正完先不做任何分析把校正前后影像并排打开用“卷帘”工具来回拉。如果肉眼能看到明显的亮度突变或色阶断层说明校正参数有问题。这个土办法比任何统计指标都直观。夜间灯光校正没有一劳永逸的公式每个研究区的大气条件、传感器状态都不同回归系数必须自己拟合。希望帮到你。本文还有配套的精品资源点击获取