ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

DMSP-OLS夜间灯光数据在ArcGIS中的三步校正实操指南

DMSP-OLS夜间灯光数据在ArcGIS中的三步校正实操指南 简介夜间灯光数据作为衡量城市活动空间格局的重要遥感指标被广泛应用于城市化监测与经济社会估算。然而DMSP-OLS传感器搭载于多颗卫星不同年份、不同传感器的DN值存在增益差异且城市中心易出现饱和值直接使用原始影难以支撑长时序比较。为此学界普遍采用相互校正、饱和校正与连续校正三步流程以参考影像统一尺度、修复极亮像元、消除年际波动。ArcGIS作为常用桌面GIS平台提供了栅格计算器、焦点统计、掩膜提取等完整工具链能够实现上述校正操作并保持投影与像元对齐。工程实践中需要注意环境捕捉设置、背景值处理和表达式规范才能真正产出可比的夜间灯光序列。本文以详细步骤方式梳理了DMSP-OLS夜间灯光数据在ArcGIS中的校正流程与避坑要点适用于城镇化研究、GDP空间化等应用场景。1. 夜间灯光数据校正为什么ArcGIS里这步没人替你省做过DMSP-OLS夜间灯光数据的人都有个共同感受数据刚从NGDC官网下下来的时候影像上城市中心亮得刺眼边缘乡镇却暗到几乎看不见不同年份之间同一座城市的DN值还忽高忽低——这不是数据坏了而是原始产品没做过任何校正。DMSP-OLS本身是面向气象观测的传感器不同卫星、不同年份的传感器灵敏度存在差异直接拿来做城市扩张或GDP拟合结果基本不能用。本文要讲的这份《夜间灯光数据校正在arcgis中的具体操作步骤非常详细.docx》就是把相互校正、饱和校正和连续校正这三步落到ArcGIS桌面端的完整操作文档从新建工作空间到栅格计算器写回归公式每步都有界面截图对应的文字描述。适合正在写论文、做城镇化研究或者第一次拿到F182013数据不知道从哪下手的从业者照着文档点一遍比自己瞎试一星期省事得多。2. 先搞懂为什么要校正DN值不是拿来就能用的2.1 DMSP-OLS数据的三宗罪不连续、饱和、过拟合DMSP-OLS夜间灯光影像的原始DN值范围是0到63其中63代表饱和。这意味着上海陆家嘴和北京国贸的DN值都是63但实际上这两个区域的夜间辐射强度显然不在一个量级。更麻烦的是DMSP卫星有F10、F12、F14、F15、F16、F18等多颗卫星在轨运行同一时期至少有两颗卫星同时成像而不同卫星的传感器增益设置完全不同。比如1999年F12卫星和2000年F14卫星拍摄的同一座城市灯光DN值可能差出30%以上这个差异纯粹是传感器差异不是真实的灯光变化。这就带来三个连锁问题一是多期影像无法直接比较做不了时间序列分析二是城市中心饱和区DN值恒定在63掩盖了真实的城市内部结构差异三是用原始数据做回归分析时饱和区会成为异常高杠杆点把模型系数带偏。校正的本质就是找一个基准把不同年份、不同传感器的影像统一到一个可比尺度上。2.2 三步校正的底层逻辑先定基准再削峰最后归一业内通用的校正流程分三步顺序不能乱。第一步是相互校正Mutual Correction做法是选一颗卫星、某一年的影像作为参考影像通常是F121999然后用其余影像的DN值对该参考影像做二次多项式回归把传感器差异拉齐。第二步是饱和校正Saturation Correction针对DN值等于63的像元做处理常见做法是用NDVI或城市建成区边界做掩膜对饱和像元赋值一个更高、符合实际辐射水平的DN值。第三步是连续校正Continuous Correction目的是消除不同年份影像之间的异常波动做法是统计各年份影像的累计DN值设定一个基准年份将其他年份影像乘以一个比例系数让整体亮度水平连续变化。从ArcGIS操作角度来看三步分别对应栅格计算器、按掩膜提取和栅格重分类工具都不复杂难点在于参数设置和中间数据的组织。文档里给的参数表很关键如果只是顺手点一下对话框很容易漏掉环境设置里的捕捉栅格和像元大小这一步漏了后面所有栅格计算的坐标系和分辨率都会乱掉。3. 准备工作数据下载与ArcGIS环境设置3.1 数据下载与文件整理不要把所有影像堆在一个文件夹DMSP-OLS数据从NGDC官网下载文件命名规则是FHHHHHHVVVV比如F152001.v4b_cf_cvg.tgzF代表卫星后四位是年份。注意有v4b和v4c两个版本v4c是去除了云覆盖影响的版本做长时间序列分析建议用v4cv4b保留了云覆盖信息但噪声较多。下载后是tgz压缩包解压后每个年份包含三个文件稳定灯光影像stable_lights、平均可见光影像avg_vis、云覆盖影像cf_cvg。做校正只用稳定灯光影像其他两个可以归档不管。拿到数据后在本地磁盘建一个工作目录我一般按这个结构组织E:\nightlight\ ├── raw\ # 原始下载的tgz包和解压后的tif ├── clip\ # 裁剪后的中国区影像或研究区影像 ├── reproject\ # 重投影后的影像 ├── correction\ # 校正中间结果 ├── mask\ # 夜间灯光以外的掩膜数据 └── output\ # 最终校正结果注意raw目录下只放原始文件任何中间处理结果都不要写进去。ArcGIS的栅格计算器在计算时如果输出路径包含中文字符或空格偶尔会报“Unable to save raster”英文路径能避开绝大多数这类问题。数据准备好后在ArcMap里打开ArcToolbox确认Spatial Analyst扩展模块已勾选栅格计算器、按掩膜提取、重分类这几个工具都属于这个模块没启用的话工具全灰。3.2 统一的投影与范围研究区一致是校正的地基所有校正操作开始前先把所有影像统一到一个投影坐标系和一个研究区范围。中国区域的项目我一般用Krasovsky_1940_Albers等积圆锥投影这是国内地理学界用得较多的投影方式如果是全球范围的研究用WGS_1984_Web_Mercator_Auxiliary_Sphere也行但要注意面积变形问题。具体操作分两步。第一步用投影栅格工具先指定输出坐标系为Albers投影同时设置地理变换从WGS84到Krasovsky_1940需要选相应的地理变换参数国内通常用“Krasovsky_1940_To_WGS_1984”这个变换重采样方法选双线性内插。第二步用按掩膜提取或栅格裁剪工具把研究区边界矢量作为掩膜统一影像范围。掩膜矢量必须和影像位于同一投影坐标系否则裁剪结果会出现边缘偏移或空白条带。这里有一个很容易忽略的细节投影栅格时ArcGIS对话框下方有一个“环境”按钮点击后需要设置捕捉栅格为某一期基准影像。这样所有年份的影像在投影后像元边界完全对齐后续栅格计算不会出现因为边界错位导致的像元错位和重采样噪声。文档里专门有一页讲这个环境设置我看完才发现自己之前一直漏了这一步。4. 核心操作把三步校正落实到ArcGIS对话框4.1 相互校正栅格计算器里写回归公式相互校正是整个流程里最核心的一步。先在ArcToolbox里打开Spatial Analyst工具→地图代数→栅格计算器把参考影像通常选F121999和待校正影像都加载进来用“栅格计算器”的界面写回归表达式。参考影像的DN值作为自变量待校正影像的DN值作为因变量用ArcGIS的波段集统计工具获取两个栅格之间的回归系数。如果嫌波段集统计给出的线性参数不够精确可以用ArcToolbox里的“波段集统计”配合Excel做二次多项式拟合得到系数a、b、c后写回栅格计算器。假设参考影像为ref待校正影像为img提取的回归系数为a0.0755、b0.0924、c2.983校正表达式为# 栅格计算器表达式输出校正后的影像 Con(img 0, a * img * img b * img c, 0)这里用Con条件函数是为了把背景0值排除在回归校正之外。0值像元本身没有灯光信号如果参与二次多项式计算会得到一个非零的背景值导致暗区出现大量伪灯光。参数a、b、c来自回归拟合不同年份的影像系数不同需要逐对影像单独计算不能复用同一组系数。校正完成后再做一次波段集统计确认校正后的影像均值与参考影像的均值差异在5%以内这一步能快速判断回归是否到位。4.2 饱和校正用建成区掩膜修复DN63的像元饱和校正是最容易出问题的一步。直接用栅格计算器把DN值63改成某个常量比如80是最粗暴的做法会让影像看起来像贴了一块补丁。文档推荐的做法是用夜间灯光数据自带的城市边界信息或者用NDVI数据辅助。一个可行的方案是先对校正后的影像做重分类将DN值大于60的像元提取出来作为饱和区掩膜然后基于周边非饱和像元的均值进行插值再替换饱和像元的值。在ArcGIS里可以这样操作先用栅格计算器提取饱和区# 提取饱和区DN值等于63 Con(img_corrected 63, 1, 0)然后使用焦点统计工具对非饱和区像元计算3×3窗口的均值再用栅格计算器把饱和区像元的值替换为周边均值# 用邻域均值替换饱和像元 Con(img_corrected 63, FocalStatistics(img_corrected, NbrRectangle(3, 3, CELL), MEAN), img_corrected)注意焦点统计窗口的选择3×3窗口适合城市尺度较小的区域如果研究区是长三角这种超大范围用5×5窗口会让边界更平滑但计算量会变大。另外焦点统计会把背景0值也算进均值里在饱和区边缘会产生暗色阴影。解决的办法是先做一次零值掩膜让背景区域不参与统计或者在Con条件里先判断周边像元是否大于0。这个细节文档里有专门提示实际操作中我踩过两次坑一次是城市边缘出现暗色环带一次是饱和区替换后亮度比周边低了20%都是窗口和背景值没处理好。4.3 连续校正让不同年份的DN值真正可比连续校正是把各年份影像拉到同一个亮度基准上。具体原理是统计每一期影像的累计DN值即所有像元DN值之和以某一期通常是中间年份或基准年份的累计DN值为基准计算其他年份的校正系数然后将各年份影像乘以对应系数。操作上分成两步先用分区统计或栅格计算器的合计功能统计各年份影像的DN值总和然后在Excel里计算系数最后再回到栅格计算器做乘法。各年份累计DN统计结果如下年份累计DN值校正系数基准2000199928631.451.1523200032989.211.0000200135721.880.9238200230119.471.0953200338275.620.8618在栅格计算器里执行乘法# 连续校正将2002年影像乘以系数1.0953 img_2002_corr * 1.0953这一步看起来简单但有一个隐藏问题乘以系数后DN值不再是整数会变成带小数的浮点数。如果后续要进行栅格转矢量或叠加分析建议先对结果做取整用栅格计算器的Int函数把浮点转为整型。另外如果研究区内有DN值为0的背景区域乘法不会改变0值但会对极低DN值比如1的小像元产生放大效应导致噪声被放大。文档建议在连续校正前先对影像做一次低通滤波把噪声压制掉再做乘法实际操作中用焦点统计里的多数值滤波即可窗口大小为3×3。5. 避坑指南夜间灯光校正的五个典型翻车现场5.1 影像裁剪后背景值不是0导致回归系数失真现象裁剪出的影像里背景区域不是0而是255或NaN栅格计算器直接报错或者回归拟合的R²极低。原因原始影像的背景值是0但经过投影、重采样或裁剪后背景值可能被重采样为255或异常值。尤其是双线性内插会把边缘的0值插值成非零小数而按掩膜提取时默认的NoData值在不同工具间并不统一。解决裁剪后的影像先用栅格计算器做一次条件赋值把大于100的像元强制设为0确保背景统一为0。常用表达式是Con(img 100, 0, img)处理后再进行后续校正。5.2 栅格计算器表达式超过长度限制报错现象写的Con嵌套表达式很长点击OK后ArcGIS弹出“Invalid expression”或“Expression is too long”。原因栅格计算器的文本框并不是无限长的而且中文字段名、带空格的字段名比如“F152000”和“F15 2000”都会增加字符解析负担。更隐蔽的是如果栅格图层名里带有文件路径中的空格计算器会自动加引号表达式长度暴增。解决把影像重命名为简短无空格的文件名比如f152000.tif存放在英文路径下表达式拆成两步先算中间结果再算最终结果如果表达式仍然过长改用Python窗口运行Raster Calculator的ArcPy版本比如arcpy.gp.RasterCalculator_sa(expression, output)这个接口对表达式长度更宽容。5.3 多期影像像元边界没有对齐相减后出现条纹噪声现象两期校正后影像做栅格相减结果图上一块亮一块暗像格纹状。原因各期影像投影时没有设置捕捉栅格输出像元网格的起始点不同导致不同影像的同一位置实际覆盖的地面范围有微小偏移。栅格相减时边缘像元错位产生了虚假差异。解决重新执行投影栅格在环境设置里打开“捕捉栅格”选择基准年份的影像作为捕捉参考同时在“像元大小”中勾选“与捕捉栅格相同”所有影像写到一个空的文件地理数据库里避免各期影像另行指定输出范围。5.4 饱和区校正后出现“黑洞”现象对饱和区用焦点统计替换后城市中心不再是亮点反而比周边暗像城市中心被挖掉一块。原因焦点统计计算均值时把饱和区旁边刚校正过的非饱和暗像元也包含在内导致均值偏低。另外饱和区本身尺度大于3×3窗口时窗口内全是饱和像元用均值替换后饱和像元被低估。解决加大窗口尺寸到5×5或7×7但优先采用迭代替换每次只替换饱和区外扩一圈的像元下一轮再替换更靠内的一圈。在ArcGIS里用“条件函数焦点统计”循环三次能明显改善暗洞现象。循环可以用模型构建器完成或者手动执行三次栅格计算器每轮用上一轮的输出作为输入。5.5 连续校正后整体DN值下降反而不如原始数据现象做完连续校正统计累计DN值发现比基准年份小了30%空间分布也失真了。原因校正系数计算错误。常见的错误是把待校正影像自身的累计DN值作为基准来算系数或者把参考影像的范围没有做掩膜统一导致累计DN值统计口径不一致。另外如果某些年份影像被裁剪成不同范围累计DN值天然不可比。解决所有参与连续校正的影像必须确保范围和像元大小完全一致——统一用同一份研究区矢量掩膜裁剪计算校正系数前把所有影像的背景0值区域排除在统计范围之外用分区统计排除NoData的选项校正系数建议用中间年份的影像做基准不要用首尾年份首尾年份往往是传感器状态最不稳定的时期。6. 验证校正效果三个必须做的小实验校正做完不能直接拿去分析花十五分钟做三个验证实验能避免你拿着一套错误数据跑完整个研究。第一个实验是抽几个城市的DN值量级检查。在ArcMap里用识别工具或值提取至点工具把校正前后同一点位的DN值列出来对比。拿上海、北京、成都各选五个点校正前市中心DN值63郊区10到30校正后市中心应该在63以上饱和区被修复后通常会调到80左右郊区DN值应略有浮动但幅度不应超过10%。如果郊区的DN值在校正后掉到1以下说明回归系数有问题。第二个实验是热点对比。选一个典型城市如西安把校正后的影像做成密度图看城市中心的热点范围是否与已知的建成区边界吻合。用栅格计算器提取DN值大于30的区域和Landsat影像或土地利用数据做目视对比重合度应高于90%。这一步能直观暴露饱和校正的质量如果热点范围明显小于建成区说明饱和区替换不足如果明显大于建成区说明背景噪声被放大了。第三个实验是趋势合理性检查。把校正后各年份的累计DN值画成折线图看是否符合研究区的社会经济规律。比如某城市要做2000到2013年的城镇化分析期内的累计DN值应该逐年稳步上升不应出现“V”形或“倒V”形的剧烈波动。如果累计DN值在2003年突然比2002年还低15%而这段时期并没有大规模停电或经济倒退那就说明连续校正系数算错了。做完这三个检查后再进回归分析或者空间统计就比较放心了。我自己的习惯是校正结果保存一个备份原始数据、中间文件、最终结果分别存放万一哪一步做错还能重来——夜间灯光校正在ArcGIS里每一步都有后悔药但前提是你把每一步的中间结果都留好。这份文档教的流程并非最短路径却是最不容易翻车的一种。从F121999的选取到栅格计算器里的Con嵌套再到连续校正的系数表整条线走完大概需要一下午。希望这篇文章能帮你绕开我当年踩过的那些坑把时间花在分析本身而不是数据清理上。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进