1. 这套夜间灯光数据，到底解决了什么问题
第一次拿到1992-2024年跨年度的夜间灯光数据集时，说实话我的第一反应不是兴奋，是头大。资料里厚厚一摞文件告诉我，这批栅格数据来自DMSP-OLS和NPP-VIIRS两代传感器，中间还隔了好几个卫星版本、多次辐射定标修正。如果直接拿原始影像拼成时间序列跑回归，基本等于在噪声里找信号。这篇文章把我整理全国、分省、分市尺度经过矫正的夜间灯光数据时踩过的坑、用顺的流程记录下来，希望能帮你少走一点弯路。
这套经过矫正的夜间灯光数据，核心价值在于“可比”。它把1992年到2024年跨三十多年的卫星观测统一到同一把尺子上，让不同年份、不同省份、不同城市之间的灯光强度可以直接相加、比较、回归。适合的人群很明确：做空间计量和区域经济研究的硕博生、城乡规划的研究人员、长期从事经济地理分析的工程师。如果你只是临时需要一个灯光总量变量，拿到矫正后的面板数据直接用就够；如果你想参与数据处理，那第二部分讲的矫正流程就是必读内容。
1.1 两代传感器，两条“不搭”的时间线
先交代背景。DMSP-OLS夜间灯光影像覆盖1992到2013年，主要产品是年度稳定灯光合成网格，数值范围通常整理成0到63的整数DN，空间分辨率约1公里。NPP-VIIRS从2012年4月开始发布月度合成产品，分辨率约500米，单位是辐射亮度值。两代传感器在光谱响应、定标方式、饱和特性上完全不同，数值一个无量纲一个带物理单位，直接拼接出来的序列会出现系统性断裂。
更麻烦的是DMSP这代本身也不省心。1992到2013年间有F10、F12、F14、F15、F16、F18等多颗卫星交替运行，某些年份同时有两颗卫星在采集，同一区域同一年的影像数值却不一致。这就是为什么行业内反复强调要先做“相互校正”，而不是简单把tif文件按年份排在一起。
1.2 “矫正”二字的完整含义
很多人把矫正理解成“去噪”，实际上这里的矫正动作要复杂得多。我拆解下来至少包含四层处理：第一是相互校正，消除不同卫星、不同年份之间的系统性亮度差异；第二是去饱和，把城市核心区被压到63的灯光信息重新展开；第三是连续性平滑，处理相邻年份之间因为卫星切换而产生的跳变；第四是跨传感器定标，把VIIRS辐射亮度统一到与DMSP可比的光谱尺度上。
只有这四步都跑完，才能说这份数据是“经过矫正的”。更严谨地讲，矫正的最终目标应当是：同一个像元在不同年份的数值变化，只反映地面灯光的真实变化，不反映传感器更替、卫星版本差异和大气路径差异。这句话是我做质量验收时反复问自己的，建议你也把这条标准贴在工位旁边。
1.3 谁需要它，能拿来做什么
从用途上看，这套数据主要解决四类问题：一是测度城市化进程，比如城市实体地域的扩张范围和速度；二是作为经济变量的代理指标，最常见的就是用夜间灯光总量拟合GDP、电力消费和碳排放；三是做城市蔓延和空间结构研究，比如识别开发区的灯光轮廓；四是为生态学提供人为活动的空间输入，评估光污染影响。
分省分市版本尤其受社科研究者欢迎，因为栅格数据对他们来说门槛太高，而按行政区汇总好的面板数据可以直接丢进Stata或者R跑固定效应模型。这也是为什么我在第三部分花了大量篇幅讲“从栅格到面板”的实操细节——这是整个流程里最容易被低估、也最容易出错的一个环节。
2. 矫正方法选型：我在实践中用的四步流水线
2.1 四步走的总体结构
我自己搭建的处理流水线长这样：原始影像 → 相互校正 → 去饱和 → 连续性修正 → 跨传感器桥接 → 掩膜与分级汇总。这个顺序不是随便排的，每一步之间都有逻辑依赖。
相互校正必须先做，因为多源影像之间还没有对齐时，任何后续拉伸和拟合都会把系统误差带到模型里。去饱和放在相互校正之后，原因是饱和拉伸会改变整体数值分布，如果先拉伸再校正，会把相互校正的基准打乱。跨传感器桥接必须放在最后，因为VIIRS的噪声和负值需要先处理干净，否则转换函数会把虚假偏误放大并传染给后续年份。
2.2 相互校正：参考区域和回归模型怎么选
相互校正的原理不复杂：在同一个时间断面下，不同卫星对同一地面亮度的观测结果存在系统性差异。要消除这个差异，就得找一块“不太变化”的区域作为基准，提取该区域在每一年、每一颗卫星影像上的DN值，然后建立回归方程把所有影像校正到同一参考基准。
实际操作中，我推荐选择某省内经济活动长期平稳的县域组合作为参考区域，而不是用全国或全省。原因很简单：如果用经济快速增长的省会作为参考区，拟合出来的校正曲线会把真实的灯光增长也一起“掰平”，最后得到一条被严重低估的时间序列。我在早期就犯过这个错误，教训很深刻。
回归模型用两种形式比较多：二次多项式 y = a·x² + b·x + c，适合DN值偏低、需要把低亮度区细节拉开的情形；幂函数 y = a·xᵇ，适合更严重的非线性压缩。拟合时要在参考区域内随机抽几万个像元，按DN十分位分层抽样，保证低亮度区也有足够的样本权重。
以二次多项式拟合为例，我实测下来二次项系数通常在0.002到0.006之间，一次项在0.75到0.95之间，常数项基本接近0。如果拟合结果出现二次项大于0.01或者一次项小于0.6，就需要立刻怀疑样本里是否混入了云层、火点等噪声像元，不能继续往下走。
2.3 饱和校正的三种落地思路
DMSP数据的0到63整数范围在城市核心区特别容易“顶格”，大量像元停在63，真实灯光强度差异被完全抹平。要解决这个问题，有几种实际可操作的路径。
第一种是基于辅助数据分解：利用NDBI建筑指数、人口密度或者建设用地数据，把饱和区的灯光潜力重新分配到各个像元。常见做法是先建立非饱和区灯光强度与建筑密度的回归关系，再把回归残差叠加到饱和像元上。这个办法精度最高，但太依赖同一年份的辅助数据，1990年代很多地方没有像样的土地利用资料，很难落地。
第二种是空间扩展法：假设饱和区外侧的灯光梯度符合某种衰减曲线，用城市边缘的非饱和像元去插值估计中心区强度。说人话就是“用城市边缘的灯光亮度推城市中心的亮度”。操作上可以用高斯核函数，但边界定义对结果影响很大，需要反复调试。
第三种是借用VIIRS反演：在2012到2013年重叠期内，用VIIRS同位置的辐射值来估计DMSP饱和区的数值分布，做局部回归或分位数映射。我个人经验是，饱和像元至少要超过50个才值得做局部回归，样本太少就直接用全局幂函数，别加戏。
去饱和做完之后，必须检查直方图。如果DN=63的频数仍然明显凸起，说明拉伸没有真正生效，要回头调整模型参数，而不是直接进入下一步。
2.4 连续性修正与跨传感器桥接
连续性修正常被忽略，但对面板数据来说至关重要。多源拼接之后，很容易出现“锯齿形”序列：比如某省2005年的灯光总量突然比2004年低20%，这不是经济收缩，大概率是卫星切换导致的。我用三年中心移动平均来处理明显断点，思路类似于S-G滤波。但移动平均窗口不能超过5年，否则真实的经济周期波动信息也会被一起磨平。
跨传感器桥接就是把DMSP和VIIRS两根线接成一根绳。我的做法是：取2012和2013年两个重叠年份，先把VIIRS月度合成聚合到1公里分辨率，做负值置零和火点掩膜，再到各区域分别拟合幂函数转换模型。
为什么选幂函数而不是线性？因为DMSP在高亮度区是压扁的，VIIRS则相对线性，两者关系天然是一条下凸的增长曲线。实际拟合时幂指数大多落在0.9到1.15之间，常数项很小。转换完成之后别急着做统计，先抽20个像元做人工对比，确认城市亮斑轮廓和DMSP原始高分影像一致，这一步能救回很多后期问题。
3. 从栅格到省市级面板：重投影、分区统计与阈值设定
3.1 数据源、分辨率和坐标系的基本配置
做分省分市数据，至少需要三类输入：已经矫正好的夜光栅格年度合成文件，最好是GeoTIFF格式；分年份的省市行政边界；辅助掩膜数据，包括水体、火点等。
坐标系选择是第一个容易踩的坑。DMSP和VIIRS的原始产品大多是WGS84地理坐标系，也就是经纬度网格。这种坐标系下的像元面积在高纬度地区会明显小于低纬度，直接用经纬度做面积统计会扭曲结果。我统一投影到Albers等积圆锥投影，中央经线105°E，标准纬线25°N和47°N，分辨率统一设为1公里，这样算出来的亮灯面积和灯光总强度才具备空间可比性。
3.2 栅格预处理：重投影、裁剪、重采样
预处理的核心动作是：先重投影，再按省界市界裁剪。顺序反过来的话，在地理坐标系下裁剪出的面元边界投影后会发生形变，边界像元的归属会错乱。
重采样方法推荐双线性或三次卷积，不要用最邻近。夜间灯光是一个连续变化的地表物理量，最邻近会带来明显的锯齿效应，尤其在行政边界上会出现“棋盘格”一样的伪影。分辨率统一到1公里，这个取值不是随便定的：DMSP原生约1公里，VIIRS原生约500米，统一到1公里既能抑制VIIRS高分辨率带来的椒盐噪声，又不会损失太多信息。
裁剪操作有个小技巧：先用行政边界向外扩一个像元，也就是buffer约1000米，然后再对栅格做裁剪。这样可以避免边界像元因为裁剪误差被重复计数或无意丢失。
3.3 分区统计：核心代码逻辑与误差控制
到这一步，核心问题是怎么把栅格值按行政区划汇总成表格。我平时最常用的是rasterio加geopandas这套组合。核心逻辑分两段：第一段先把行政区面要素栅格化为“像元ID数组”，第二段用numpy的bincount做加权汇总。
下面是一个简化但可运行的处理骨架：
import rasterio import numpy as np from rasterio.features import rasterize from rasterio.warp import calculate_default_transform, reproject, Resampling import geopandas as gpd # 统一投影到 Albers 等积投影 proj4 = "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" with rasterio.open("ntl_corrected.tif") as src: transform, width, height = calculate_default_transform( src.crs, proj4, src.width, src.height, resolution=1000) kwargs = src.meta.copy() kwargs.update({"crs": proj4, "transform": transform, "width": width, "height": height}) with rasterio.open("ntl_corrected_aea.tif", "w", **kwargs) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=proj4, resampling=Resampling.bilinear) # 分区统计：按面要素对应的像元ID汇总 shp = gpd.read_file("city_2020.shp").to_crs(proj4) with rasterio.open("ntl_corrected_aea.tif") as src: a = src.read(1).astype("float64") a[a < 0] = 0 id_raster = rasterize( [(g, i) for i, g in enumerate(shp.geometry)], out_shape=a.shape, transform=src.transform, fill=-1, dtype="int32" ) valid = id_raster >= 0 idx = id_raster[valid].ravel() vals = a[valid].ravel() city_sum = np.bincount(idx, weights=vals, minlength=len(shp)) city_area = np.bincount(idx, minlength=len(shp))这段代码最后能得到每个城市的灯光总强度和有效像元个数。误差控制我有三条经验：第一，不要直接拿全国栅格和边界做分区统计，必须先裁剪外扩后的区域；第二，对跨边界像元，要么用像元面积加权、要么按中心点归属，不能简单四舍五入，否则小城市的边界会“吃掉”大城市的边缘灯光；第三，如果城市数量超过几百个，一定给每个面要素加上FID，把统计中间表导出保留，方便日后回溯是哪一年、哪个边界出了问题。
3.4 指標定义与阈值选取
汇总之后要明确指标口径，否则不同文章之间的结果根本没有可比性。常用的指标有四类，直接列成表格：
| 指标 | 计算方法 | 适用场景 |
|---|---|---|
| 夜间灯光总强度（TNL） | 区域内所有大于零像元数值的总和 | 衡量城市经济总规模、空间扩张总量 |
| 平均灯光强度 | 总强度除以有效像元数 | 衡量灯光密度、平均发展水平 |
| 点亮面积 | 数值超过阈值的像元数乘以单位像元面积 | 识别城市实体地域范围、蔓延程度 |
| 灯光覆盖率 | 点亮面积除以行政区面积 | 城镇化率替代指标 |
阈值选取是最容易被糊弄的地方。文献里常见阈值有3、6、12、20，我不能说哪个绝对正确。我的处理方式是：把不同阈值下的点亮面积分别与同一年份的建成区面积、道路照明范围做相关性检验，选取相关性最高且边际变化趋于平稳的阈值。对于经过矫正的整合数据，我个人在省级尺度常用等效值≥3来识别亮灯；到了地级市尺度，考虑到边缘噪声，我会把阈值提高到5，宁可漏掉微弱灯光，也不要把噪声算进城市面积。
4. 夜间灯光数据处理常见问题与排查方法
4.1 一张速查表先解决问题
| 症状 | 可能原因 | 排查与处理 |
|---|---|---|
| 某省某年TNL突然断崖式下降 | 卫星切换、年度合成缺少月份 | 检查原始文件版本和覆盖月数，做三年中心平滑 |
| 城市边界出现一整圈高亮像元 | 投影不统一或缓冲区设置过大 | 改为Albers投影，重做裁剪，检查buffer |
| VIIRS数据出现大面积负值 | 传感器噪声、云层污染 | 用邻域中位数替换，不要直接删除像元 |
| DMSP数据城市核心全为63 | 饱和矫正未执行 | 回到去饱和流程，检查直方图 |
| 2012-2013年变化幅度超过150% | 跨传感器桥接缺失 | 检查转换模型是否拟合进异常样本 |
| 分区统计面积与官方面积明显不符 | 边界年份与栅格年份错配 | 核对行政区划变更记录，按年份就近匹配边界 |
| 小城市灯光为0 | 阈值过高或重采样丢边 | 降低阈值，扩展边界缓冲区 |
4.2 三个踩过的坑的详细复盘
第一个坑是参考区域选错。早期我直接用某经济增长很快的省会做相互校正基准区，结果整条时间序列的灯光增速被显著压低，后续所有面板模型跟着失真。后来改成“省内若干经济增长几乎停滞的县组合区域”，才把基准稳住。这条经验建议所有做矫正的人刻在脑子里：参考区一定要“稳”，不是越发达越好。
第二个坑是负值处理。VIIRS月度数据里的负值大多来自暗像元的随机噪声，但也确实包含微弱信号。我曾经统一置零，结果欠发达县域的灯光总和莫名其妙少了近两成。正确的做法是先做3×3邻域统计：如果邻域中位数大于等于0，就用中位数替换当前像元；如果邻域本身都是负值，再置零。这个逻辑既保留了微弱信号，又压掉了噪声。
第三个坑是行政边界错配。直辖市成立、撤地设市、地区改市这类区划调整，在三十多年里发生过很多次。拿同一份2020年边界去汇总1995年数据，会把老的“地区”和新的“地级市”混在一起，形成伪变化。我现在准备五套边界，分别是1995、2000、2010、2015、2020版，按年份就近匹配，同时生成一个“可比口径”变量：把所有历史区划归并到2020年行政区划后再重新统计一遍。这一份额外的面板在稳健性检验里非常有用。
5. 数据质量验证：我的质检套路
5.1 用“伪不变区域”检验连续性
自动流程跑完之后，我第一件事不是看统计结果，而是找“伪不变区域”做连续性验证。所谓伪不变区域，就是过去三十年几乎不会发生灯光建设的区域，比如大型自然保护区核心区、远离城镇的连片农田区。这些区域的灯光均值理论上应该是一条近乎水平的线。
检验方法很简单：把伪不变区域各年的灯光均值拉出来画时序折线，允许的波动我一般控制在±5%以内。如果哪一年突然抬升或者下降，说明矫正环节里有遗漏，要回到相互校正和连续性修正重新检查。
5.2 用影像轮廓目视复核
机器校验只能抓统计层面的异常，数据最终能不能用，还是要靠眼睛。我通常会随机抽20个不同规模的城市，把矫正后的影像做成逐年对比的图片序列，重点看城市亮斑的几何轮廓是否平滑扩张、有没有突然消失的孤立亮块、海岸线附近有没有不正常的亮线。
这种目视检查看起来土，但真的能发现自动流程发现不了的问题。比如某省会的亮斑某一年突然被“切掉”一半，多半是要素ID顺序错乱；某沿海城市出现一条细长亮线，很可能是漏掉了渔业灯光的掩膜。这些细节，写代码时很难提前想到。
5.3 与GDP、人口、用电量做交叉验证
最后是外部统计验证。以地级市为样本，把TNL与GDP、常住人口、全社会用电量分别做对数回归，拟合优度R²一般应该在0.7以上。如果某一年突然掉到0.5以下，那一年基本可以判定有问题。
需要说明的是，0.7只是我的经验值，不同地区的数据基础差异很大。这个方法适合当日常体检指标，不适合当严格的学术检定标准。一旦发现问题，我还会带上前几年的差值图一起定位，找到具体是哪个区域、哪个年份出了问题，而不是笼统地返工。
6. 一些经验和补充建议
所有中间版本文件都要保留。矫正流程每一步都会改变栅格数值，不保留中间结果，出了问题你根本不知道是引入在哪一步的。我的习惯是按“原始数据 / 相互校正后 / 去饱和后 / 桥接后 / 汇总表”五个目录存放，文件名里带上处理日期和版本号。
如果是非GIS背景的研究者，我的建议更直接：不必自己从头做整套矫正，优先使用公开的组织化数据产品。但拿到数据后，还是要把第二部分讲的四步流程当作检查清单过一遍，重点看数据描述文件里有没有说明跨传感器桥接是怎么做的、有没有用重叠年份做验证。没写清楚的，默认它没做过，自己再验一遍。
动手做面板时，建议把统计脚本封装成“输入栅格 + 输入边界 + 输出表格”的固定范式，并把每次统计的生成日期写入文件注释。这样半年后回看数据，你还记得当年是怎么处理的。
我自己现在跑长时序灯光数据，已经养成了一个死习惯：数据汇总完成后，先抽一年甚至连续两年做人工抽查。不是不信任流程，而是这种跨三十多年、覆盖两级行政区的数据，差一个像元都可能被模型放大成“显著结果”。多看两眼原始栅格，远比在回归结果里找p值有意义。
