“1992-2024年全国、分省、分市的夜间灯光数据还要经过矫正”——我每次看到这个描述第一反应都是夜光遥感不是美国NASA和NOAA公开下载吗折腾矫正意义在哪里真自己跑一遍数据就知道DMSP到VIIRS之间那道坎能直接把长时序趋势干崩塌。这套数据的价值不在“图片好看”而在“可比性”。今天我就把长时序灯光数据的加工逻辑、矫正流程和落地用法完整捋一遍手把手把踩过的坑都标出来做城市研究、区域经济、能源负荷的朋友可以直接抄作业。1. 为什么1992-2024年夜间灯光数据“非矫正不可”1.1 DMSP和VIIRS的底层差异有多大现在能公开获取的夜间灯光影像主要来自两代传感器先说清楚家底DMSP-OLS覆盖1992-2013年空间分辨率约2.7公里量化位数只有6比特DN值范围0-63。它的优点是历史长缺点是灵敏度低、城市中心极易饱和。NPP/NOAA-20 VIIRS DNB覆盖2012年至今分辨率约500米-750米量化位数高动态范围大能捕捉微弱灯光但也因此混入了大量火灾、渔船、闪电这类瞬时事件噪声。两代传感器在Quantization、过采样、杂散光处理上都不同。拿同一座城市2013年前后的影像直接对比大概率会看到一夜之间“变暗”或“变亮”这不是城市真实变化是仪器切换的假象。如果你直接把这两段拼接成一个长序列任何回归分析都会算出莫名其妙的断点。1.2 所谓“矫正”到底解决什么问题我理解的“经过矫正”至少要处理三类问题第一是连续性矫正。DMSP和VIIRS量纲不同需要建立跨传感器转换模型把两段数据统一到同一亮度参考系。第二是饱和矫正。DMSP城市中心像元DN值大量触及63后亮度不再增长中心城区看起来像个“平地”。如果不修正后续的灯光增长率会被系统性低估。第三是噪声清理。VIIRS原始数据里约0.3 nW·cm⁻²·sr⁻¹以下的像元经常混有背景噪声城市边缘、乡镇区域的伪灯光非常影响统计数据。只有把这三个问题解决掉数据才能用于年份之间的比较才能做出“1992-2024”这种长时序面板。1.3 全国、分省、分市三级粒度是怎么定的这套数据在空间上做了三级输出全国栅格、分省汇总、分市汇总。这么做不是拍脑袋是因为不同研究尺度的需求完全不一样全国总量用于宏观趋势比如灯光增长与GDP增长是否匹配分省数据适合区域差异分析、空间收敛性研究分市数据是最常用、最棘手的部分可用于城市扩张、夜间经济活动密度、市政能源规划等。从栅格到统计表需把不同年份的灯光影像统一重投影到等积投影再叠加对应的行政区边界。这里有一个关键决策如果行政区划有调整是采用“最新边界回溯统计”还是“逐年边界统计”我的建议是主产品用最新边界回溯确保行政区单元不变时灯光变化不受边界变动污染同时对有区划调整的年份额外给出说明列。这样才能保证面板数据的可比性。2. 长时序夜间灯光矫正全链路拆解2.1 数据源选择与口径确认做矫正前先把原料选对。DMSP数据一般下载NOAA NGDC的Version 4 Stable Lights影像它已做过去除云层和火光处理。VIIRS数据推荐使用VNP46A2月合成或NCEI的Nighttime Lights月度产品。注意月度产品里仍然包含部分极光、船舟灯光和气体燃烧信号需要后续掩膜。拿到原始数据后第一步是统一“口径”将所有年份影像裁剪至同一研究边界重投影至等积投影比如Albers等积投影或Mollweide统一分辨率。常见做法是把DMSP重采样至1公里VIIRS也聚合到1公里这样计算量适中噪声也降低。我见过不少人直接拿500米VIIRS和2.7公里DMSP对比数值差一大截那不是数据问题是分辨率没对齐。重采样时注意聚合方法建议用“近邻”或“平均值”不要用双线性插值否则会在城市边缘产生人为渐变区域。2.2 跨传感器桥接这是最核心的一步跨传感器桥接的常规思路是找两代传感器重叠年份2012、2013或2014视产品而定同一地区的DMSP DN值和VIIRS辐亮度值一一对应建立回归模型。具体操作提取重叠年份影像对剔除DMSP饱和像元DN63剔除VIIRS零值像元和异常高值像元在非城市的均匀背景区域再取一部分零点样本避免模型只看亮区拟合线性或二次多项式常见近似关系为DN_v a b * DN_d或指数型方程。举例某版产品在2013年样本中拟合得到VIIRS亮度 -0.13 0.74 × DMSP DN线性近似这里只是示例系数拟合R²在0.8左右。用这个系数就可以把1992-2013期间的DMSP值统一换算到VIIRS亮度口径。反向操作也可以但为了保留高分辨率、大动态范围的信息我建议统一到VIIRS口径。做完回归之后要做“反推验证”把2013年DMSP原始值代入模型生成的VIIRS估计值与真实VIIRS影像做差值图。差值应集中在少数噪声区域如果出现系统性偏移说明回归模型需要加入城市饱和修正项。2.3 饱和校正和背景噪声滤除的方法DMSP饱和校正常用的有三类方法局部最大值补偿对饱和簇中心邻域用周围未饱和像元的亮度梯度外推内部亮度不变目标法选取亮度常年稳定的区域如沙漠、高海拔荒原里的零星居民点作为辐射定标场用其多年DN值变化来修正社会经济辅助法结合人口密度、不透水面比例、GDP等数据在饱和区域重新分配灯光总量。这个方法精度高但需要额外数据集适合区域级研究。VIIRS噪声滤除不能直接一刀切。夜间灯光的有效信号在0.1-10 nW·cm⁻²·sr⁻¹范围但不同区域背景噪声不同沿海地区渔船灯光偏多、西部无人区偶有火点信号。稳妥的办法是用多期影像的“中值合成”来压制瞬态光源——我的经验是至少取三年同期月度数据做三分位数然后对每个像元计算时间稳定性指数瞬态光源在时间维度上往往忽高忽低很容易识别。2.4 分省市统计里的面积权重与NoData陷阱生成分省、分市统计值最常见的错误是直接把栅格数值在边界范围内求平均殊不知这个平均没有考虑NoData和单位面积载荷。正确的统计链路是把行政区边界作为mask裁剪灯光栅格排除该范围内的NoData像元否则无数据区域会被当成0拉低均值如果是统计“平均灯光亮度”直接计算非空像元均值即可如果是统计“灯光密度或总通量”需要按像元面积加权。因为投影变形会导致不同纬度像元实际面积不同所以不能简单总亮度除以像元数。实际操作中我更倾向于先输出三个字段灯光总量、灯光均值、点火像元数亮度阈值的像元个数后面分析时按需组合。这三个字段分别对应经济总量、平均夜间活动强度、空间覆盖范围互不替代。3. 自己动手矫正并生成省市指标的实操示例Python3.1 为什么放弃纯ArcGIS改走Python批处理ArcGIS的Zonal Statistics很好用但做1992-2024三十多年分市统计手动操作会崩溃。QGIS的批量也有限。Python的rasteriogeopandasregionmask组合能完全自动化而且代码可复现适合日后数据更新。我会把流程写成三个脚本预处桥接校正、噪声滤除、分区统计。这里分享分区统计的核心片段完整版根据自己数据微调。3.2 从矫正后tif到分市统计表的代码流程假设你已有矫正后的年度栅格corrected_lights_YYYY.tif和行政区划city_boundary.gpkgimport geopandas as gpd import rasterio import numpy as np from rasterio.mask import mask import pandas as pd # 读取行政区划一般为市界 cities gpd.read_file(city_boundary.gpkg) # 打开灯光栅格 raster_path corrected_lights_2024.tif src rasterio.open(raster_path) # 计算分区统计 records [] for idx, city in cities.iterrows(): out_image, out_transform mask(src, [city.geometry], cropTrue, nodata0) arr out_image.astype(float32) # 排除NoData vals arr[arr 0] if len(vals) 0: continue total float(vals.sum()) mean float(vals.mean()) lit_area float((vals 0.3).sum() * (abs(out_transform[0]) * abs(out_transform[4]))) records.append({ 省_code: city[省_code], 市_code: city[市_code], year: 2024, light_total: total, light_mean: mean, lit_area_km2: lit_area / 1e6 }) result pd.DataFrame(records) result.to_csv(city_lights_2024.csv, indexFalse)代码很简单但三个细节要注意mask函数的nodata0要与矫正后数据的有效值约定一致通常我们把无效区域改成0或NaN如果用0就要在求均值时arr 0过滤lit_area计算时要乘以栅格分辨率并且要检查分辨率单位。如果栅格已投影为米transform[0]就是像元宽单位是米边界文件和栅格坐标系不一致时必须cities cities.to_crs(src.crs)否则裁出的全是空值。3.3 快速排查校正异常的小技巧跑完全部分区统计后不要直接进模型先画一遍诊断图。我更常用的做法对每个城市计算“相邻年份灯光均值比值”画出时间序列图。正常城市应该是平滑波动突然升降很可能是校正参数没覆盖到或行政区边界里有废弃地块。如果发现某市2012到2013年比值超出常规范围优先检查该市是否包含DMSP饱和区未修正、VIIRS噪声滤除是否不稳定而不是马上怀疑经济变化。还可以做“空间残差图”用2013年桥接模型预测的VIIRS值与真实VIIRS值相减把残差叠加到城市边界上看。残差较大且成片出现的位置就是后续统计需要特别谨慎的区域。4. 矫正后的夜间灯光数据实际应用到底怎么“接得住”4.1 城市扩张研究从连续灯光面积看“真扩张”城市扩张研究里最常用的是“灯光面积阈值法”。矫正数据统一量纲后我们可以设定一个统一阈值把超过阈值的像元视为夜间活动覆盖区然后计算每座城市逐年的“点亮面积”。比如取阈值0.5 nW·cm⁻²·sr⁻¹某市2000年点亮面积120平方公里2020年达到340平方公里这比直接用DMSP的亮度变化更稳健。但必须强调阈值的高低直接改变扩张速率建议在做城市间比较时固定同一个阈值同时报告不同阈值下的面积作为稳健性检验。4.2 经济空间分析均值灯光不等于经济活力很多论文用“夜间灯光总量”或“均值”作为GDP的代理变量。经过矫正的数据能提升拟合优度但你不要忽视两个问题均值会被城市内部高亮核心“带飞”如果只研究县域尺度最好使用灯光总量或加权稀疏指数灯光增长与经济增长并非严格同步不同产业类型的夜间能耗差别很大在服务业为主的城市里夜间灯光可能被低估。我个人的处理习惯是把灯光数据作为模型的一个核心解释变量而不是替代实测GDP。比如研究“产业结构转型对经济增长的影响”时用灯光增长率作为经济活跃度的稳健性替换能有效缓解统计口径调整导致的GDP连续性问题。4.3 面板数据使用中的动态区划匹配分市面板数据最容易翻车的地方是行政区划代码。一个地级市历史上可能经历县级市合并、撤县设区、新设地级市如果直接用2010年的行政区划代码链接到2020年的边界统计值会“漂移”。解决办法有三种统一用最新边界回溯统计全部年份构造“1992年不变行政区划”剔除后来源变区域的样本对发生边界变化的城市做“重叠加权重”换算把历史年度灯光总量按空间比例拆分归并。第三种最精确但需准备逐时代的历史边界。很多公开数据并不提供历史边界版本所以你说数据“经过矫正”一定要在元数据里标注使用的是哪一版行政区划否则用户做历史比较就是盲人摸象。5. 常见问题与排查技巧实录5.1 六个高频“翻车”场景速查表问题现象可能原因解决办法2012/2013年灯光值整体突变DMSP到VIIRS没有桥接或桥接系数不合适用重叠年份重新拟合回归检查饱和样本剔除城市中心亮度多年不变、形成“平台”DMSP饱和没修正使用局部最大值补偿或社会经济数据去饱和某市灯光总量异常偏低行政边界内大量NoData被当0统计统计时过滤NoData使用有效像元均值或面积权重偏远地区出现随机高亮像元火烧迹地、加油站瞬态光源未滤除多期影像中值合成或按时间稳定性指数掩膜前后两年灯光面积突然翻倍行政区划调整或栅格投影不匹配用统一最新边界回溯检查边界坐标系城市边缘灯光过渡太陡重采样用了双线性插值且裁剪顺序不对改用平均聚合或最近邻法先重投影再裁剪5.2 我怎么验证矫正后的数据质量验证矫正质量不能只看相关系数还要做三件事与外部数据交叉验证把灯光总量与同年度城市非农业人口、用电量、高速公路里程做相关性分析。正常情况下相关系数应在0.8以上如果某个年份骤降到0.5就需要检查该年份影像是否有质量问题或校正出偏。与同期夜间影像目视对照随机抽取30个城市逐年看灯光分布与其他地表覆盖产品Sentinel-1夜间雷达或 GlobeLand不透水面是否吻合。计算残差的Moran指数若残差在空间上高度集聚说明校正模型缺少了空间异质性需要考虑分区域桥接或加入局部修正项。5.3 补充一个小经验绝对亮度别迷信增量比较更稳做长时序数据分析这么多年我最想强调的一点是经过矫正的夜间灯光数据适合做增量和比例比较不适合做绝对亮度的“精确测量”。原因很简单任何矫正模型都会引入误差城市中心的绝对亮度误差可能还有15%-20%。但相邻年份的亮度变化、同一城市在不同时间的灯光面积扩张率在矫正后通常能够保留稳健的排序关系。所以我在写结论时会优先报告“灯光增长率”“扩张面积变化”“亮度排名变化”这类相对指标而不是说“某城市2024年灯光是1992年的XX倍”这种绝对结论。这不是数据不够好而是遥感产品本身的物理性质决定的。理解这个边界才能把长时序夜间灯光数据用得既体系化又不下次翻车。