简介这份资源面向城市研究、区域经济与地理信息分析领域的学习者和科研人员提供1992至2020年广东省经过校正的长时间序列夜间灯光数据及配套矢量文件可用于城市建成区提取、GDP空间化、人口与社会经济指标的空间化分析等场景。压缩包共181个文件约4.18MB主要包含29个tif栅格影像及对应tfw坐标文件、30个dbf属性表、30个cpg编码说明、58个xml元数据以及shp、shx、prj、sbn、sbx等矢量与索引文件兼顾栅格与矢量两种数据形态。其中DMSP/OLS影像已完成相互校正与连续性校正NPP/VIIRS影像完成年度合成、去噪与连续性校正并将两套数据拟合为统一的长时序灯光影像省去自行预处理的繁琐步骤。目前已有209人学习下载适合需要直接开展建成区识别、经济空间化建模或跨年份对比分析的用户参考使用。1. 广东省 1992-2020 夜间灯光校正数据从 DMSP 到 NPP 的一条时间线手里拿到一份标着「1992-2020 年广东省完成校正夜间灯光数据和矢量数据」的压缩包第一反应通常不是兴奋而是警惕——因为跨 29 年的夜间灯光序列天然横跨了两代完全不同的传感器DMSP-OLS1992-2013和 NPP-VIIRS2012 至今。这两代数据在辐射定标、空间分辨率、光谱响应上根本不是一回事直接拼在一起做时间序列分析出来的趋势线大概率是假的。所谓「完成校正」核心工作就是把这两段数据拉到同一个可比口径上再裁出广东省边界。这份数据能解决的是长时序城镇化监测、GDP 空间化、建成区提取、电力消费估算这类需要连续夜间灯光序列的问题适合做区域经济、城市扩张、遥感指数分析的从业者。但前提是你得知道它怎么校正的、校正到什么程度、哪些年份仍然不能信。2. 两代传感器怎么接DMSP-OLS 与 NPP-VIIRS 的校正逻辑2.1 为什么原始数据不能直接用DMSP-OLS 的夜间灯光产品稳定灯光产品有几个先天缺陷像元 DN 值被压缩到 0-63 的整数区间高亮度城市核心区普遍饱和顶到 63 就不再增长不同年份的卫星F10 到 F18 多颗星交替之间没有星上定标同年不同星、同星不同年的 DN 值不可比像元被重采样到约 1km但实际瞬时视场是 2.7km 左右存在明显的地理定位误差。NPP-VIIRS 则是辐射定标过的单位是 nW·cm⁻²·sr⁻¹动态范围大得多空间分辨率约 500m但早期数据2012-2013含有大量杂散光影响且没有像 DMSP 那样的「稳定灯光」产品需要自己过滤背景噪声。把这两段直接拼接会出现两个典型问题一是 2013 到 2014 年之间出现一个巨大的数值跳变因为量纲和饱和特性不同二是城市核心区的相对亮度关系在两段之间不一致。所以「校正」这件事本质上要解决三件事DMSP 内部的年际连续性、DMSP 到 VIIRS 的跨传感器一致性、以及广东省边界内的空间一致性。2.2 常见的校正路线与参数选择业界做 DMSP-OLS 年际校正最常用的是 Elvidge 等人提出的一套互校正流程核心是选取参考年份和参考卫星用不变目标区域通常是稳定的小城镇或郊区做回归把其他年份的 DN 值拉到参考年份的口径上。具体操作上我一般会这样做第一步确定参考影像。通常选 F16 卫星 2007 年的数据作为 DMSP 内部的参考因为这一年覆盖完整、饱和程度相对可控。第二步选取不变目标区域。在广东省内选若干个人口和建成区在 1992-2013 年间变化很小的乡镇作为「伪不变区域」提取这些区域的 DN 值做逐年回归。第三步用二次多项式回归做年际校正把各年 DN 值映射到参考年份。跨传感器校正DMSP 到 VIIRS则更麻烦。常见做法有两种一种是用重叠年份2012-2013建立 DMSP DN 值与 VIIRS 辐射值之间的幂函数或多项式关系把 DMSP 拟合到 VIIRS 口径另一种是把 VIIRS 聚合到 DMSP 的 1km 网格和 0-63 的 DN 区间反向拟合。两种方向各有优劣前者保留了 VIIRS 的辐射信息但会放大 DMSP 的饱和误差后者则牺牲了 VIIRS 的细节。import numpy as np import pandas as pd from scipy.optimize import curve_fit # 假设已有重叠年份 2012、2013 的配对样本 # dmsp_vals: DMSP 校正后 DN 值0-63 # viirs_vals: VIIRS 聚合到 1km 后的辐射值nW·cm^-2·sr^-1 dmsp_vals np.array([...]) # 从不变区域提取 viirs_vals np.array([...]) # 幂函数拟合VIIRS a * DMSP^b def power_func(x, a, b): return a * np.power(x, b) popt, pcov curve_fit(power_func, dmsp_vals, viirs_vals, p0[1.0, 1.0], maxfev10000) a_fit, b_fit popt # 把 1992-2013 的 DMSP 全部映射到 VIIRS 口径 dmsp_corrected power_func(dmsp_vals, a_fit, b_fit)这段代码的关键在于curve_fit的初值p0。如果初值给得太随意幂函数拟合容易发散尤其是当样本里包含大量低 DN 值比如 0-3时。我一般会先把 DN 值为 0 的像元剔除只保留 DN≥5 的样本做拟合因为低值区的噪声会严重干扰幂函数的形状。参数a控制整体缩放b控制曲线的弯曲程度b通常在 1.2 到 2.0 之间如果拟合出来b小于 1说明样本选取有问题大概率是混入了 VIIRS 的杂散光噪声。提示跨传感器校正的精度高度依赖重叠年份的质量。2012 年的 VIIRS 数据本身含有较多杂散光建议用 2013 年或 2014 年的数据作为重叠基准或者对 VIIRS 先做杂散光过滤再参与拟合。2.3 广东省矢量边界的裁剪与投影统一拿到校正后的栅格数据下一步是用广东省矢量边界做裁剪。这里有一个容易被忽略的坑DMSP 和 VIIRS 的原始投影不同。DMSP 常用的是等面积投影或地理经纬度投影VIIRS 的月度产品通常是地理经纬度投影WGS84。如果直接裁剪再拼接会出现像元错位。正确做法是先把两套数据统一到同一个投影坐标系下再做裁剪。我一般用 Albers 等面积投影适合中国区域中央经线 105°E双标准纬线 25°N 和 47°N因为夜间灯光分析经常要做面积统计等面积投影能避免高纬度地区的面积畸变。裁剪时用rasterio的mask功能配合geopandas读取矢量边界。import geopandas as gpd import rasterio from rasterio.mask import mask from rasterio.warp import calculate_default_transform, reproject, Resampling # 读取广东省矢量边界 gdf gpd.read_file(guangdong_boundary.shp) # 统一到 Albers 等面积投影 gdf_albers gdf.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84) # 对每一年的栅格做重投影 裁剪 with rasterio.open(nightlight_2010.tif) as src: # 重投影到 Albers transform, width, height calculate_default_transform( src.crs, gdf_albers.crs, src.width, src.height, *src.bounds) kwargs src.meta.copy() kwargs.update({crs: gdf_albers.crs, transform: transform, width: width, height: height}) with rasterio.open(nightlight_2010_albers.tif, w, **kwargs) as dst: reproject(sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsgdf_albers.crs, resamplingResampling.bilinear) # 裁剪 with rasterio.open(nightlight_2010_albers.tif) as src_albers: out_image, out_transform mask(src_albers, gdf_albers.geometry, cropTrue) out_meta src_albers.meta.copy() out_meta.update({driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform}) with rasterio.open(nightlight_2010_gd.tif, w, **out_meta) as dest: dest.write(out_image)重采样方法选bilinear还是nearest取决于你的分析目的。如果做建成区提取二值化用nearest保留原始 DN 值如果做亮度趋势分析bilinear更平滑。但要注意DMSP 的 DN 值是整数双线性重采样后会变成浮点数后续如果要做 DN 值频率统计需要先取整或做分层处理。3. 拿到数据后怎么验证三步判断校正是否可信3.1 用统计量做快速体检打开数据第一件事不是画图而是算几个基本统计量。对广东省范围内每一年的栅格计算最大值、均值、DN 值等于 63 的像元占比DMSP 年份、以及非零像元数。如果某一年 DMSP 的饱和像元占比突然从 5% 跳到 20%说明这一年的校正可能有问题或者该年原始数据本身质量差。我一般会把这些统计量做成一张时间序列表肉眼扫一遍。正常的校正后序列均值应该呈现缓慢上升趋势城镇化推进饱和占比应该逐年下降或持平因为校正会缓解饱和非零像元数应该稳步增加。如果中间某一年出现断崖式下跌大概率是那一年卫星切换或校正参数出了问题。3.2 用不变区域做交叉验证选取广东省内几个已知变化很小的区域比如偏远乡镇、自然保护区提取这些区域在 1992-2020 年的灯光均值序列。如果校正做得好这些区域的灯光值应该基本平稳波动幅度在 ±10% 以内。如果某个区域在 2013-2014 年之间出现明显跳变说明跨传感器校正没有把这一段接好。这个验证方法的好处是直观坏处是需要你对广东省的地理有基本了解。如果不太熟悉可以用夜间灯光数据的「稳定灯光」掩膜——DMSP 有官方稳定灯光产品VIIRS 也有类似的产品用这些掩膜提取的区域理论上应该年际稳定。3.3 与统计年鉴做外部对照最硬核的验证是拿广东省各地级市的 GDP 或电力消费数据做相关性分析。夜间灯光总值SUM与 GDP 的相关性在省级尺度上通常能达到 0.9 以上如果算出来只有 0.7 甚至更低要么是校正有问题要么是矢量边界和统计口径不匹配。注意这里要用「灯光总值」而不是「灯光均值」因为均值会受行政区面积影响。做相关性时建议分两段算1992-2013 一段2014-2020 一段。如果两段的相关系数差异很大比如前段 0.95后段 0.75说明跨传感器校正可能引入了系统性偏差。这时候需要回头检查 VIIRS 的杂散光过滤和 DMSP 的饱和校正参数。4. 避坑与排查夜间灯光校正里最容易翻车的五件事4.1 现象2013 和 2014 年之间灯光总值突然翻倍原因DMSP 和 VIIRS 的量纲没有统一。DMSP 的 DN 值是 0-63 的相对值VIIRS 是绝对辐射值如果不做跨传感器校正直接拼接2014 年的数值会天然比 2013 年大一个量级。解决检查校正流程里是否有跨传感器拟合这一步。如果没有用重叠年份2012-2013 或 2013-2014建立 DMSP 到 VIIRS 的映射关系。拟合时注意剔除低 DN 值样本否则幂函数会被噪声带偏。4.2 现象城市核心区的灯光值在 2000 年后不再增长原因DMSP 的饱和效应。城市核心区 DN 值顶到 63 后无论实际亮度怎么增加数值都不变。如果校正时没有做饱和校正2000 年后的城市核心区灯光序列就是一条平线。解决饱和校正的常见做法是用非饱和像元DN63建立 DN 与实际辐射的回归关系然后外推到饱和像元。或者用 VIIRS 的高动态范围数据反向拟合 DMSP 的饱和区。如果数据里没有做这一步城市核心区的长时序分析基本不可用。4.3 现象裁剪后的栅格在边界处出现锯齿或空白原因矢量边界和栅格的投影不匹配或者裁剪时没有设置cropTrue和正确的nodata值。解决先统一投影再做裁剪裁剪时用rasterio.mask.mask的cropTrue参数并显式设置nodata值通常用 -9999。如果边界处仍有锯齿检查矢量边界的精度是否足够必要时用gdf.buffer(0)修复几何。4.4 现象VIIRS 早期数据2012-2013出现大片异常高值原因VIIRS 的杂散光效应。2012 年和 2013 年的部分月份VIIRS 数据在低纬度地区会出现非真实的亮值尤其是在夏季。解决用 VIIRS 的杂散光掩膜产品过滤或者直接用 2014 年之后的数据作为跨传感器校正的基准。如果必须用 2012-2013 年数据建议做中值滤波或分位数过滤把异常高值剔除。4.5 现象校正后的灯光序列与 GDP 相关性很低原因可能是矢量边界与统计口径不匹配比如统计年鉴用的是行政边界而矢量数据是地理边界也可能是校正参数只适用于全球或全国尺度在广东省这种局部区域不适用。解决先检查边界是否一致再检查校正参数是否用了区域自适应的方法。如果用的是全球统一的校正系数建议用广东省内的不变区域重新拟合参数。另外GDP 数据本身的口径变化比如 2004 年经济普查调整也会影响相关性做对照时要留意统计口径的断点。5. 进阶用法用校正后的灯光数据做建成区提取与趋势分析校正后的数据最大的价值在于可以做长时序的建成区提取。常见做法是阈值法对每一年的灯光栅格选取一个阈值比如 DN10 或辐射值 某个分位数把超过阈值的像元标记为建成区。但固定阈值在长时序里会出问题——1992 年的 DN10 和 2020 年的 DN10 代表的实际亮度完全不同。所以更稳妥的做法是动态阈值用每一年的灯光值分布的分位数比如 90% 分位数作为阈值或者用 Otsu 方法自动选取。我一般会先用动态阈值提取建成区再计算每一年的建成区面积最后和统计年鉴的建成区面积做对照。如果趋势一致但绝对值有偏差说明阈值需要微调如果趋势都不一致说明校正有问题。import numpy as np import rasterio from skimage.filters import threshold_otsu # 对每一年的栅格做 Otsu 阈值提取建成区 years range(1992, 2021) builtup_areas [] for year in years: with rasterio.open(fnightlight_{year}_gd.tif) as src: data src.read(1) # 剔除 nodata valid data[data ! -9999] # Otsu 阈值 thresh threshold_otsu(valid) builtup (data thresh) (data ! -9999) # 计算面积像元数 × 像元面积 pixel_area src.res[0] * src.res[1] # 平方米 area_km2 builtup.sum() * pixel_area / 1e6 builtup_areas.append(area_km2) # 输出趋势 for year, area in zip(years, builtup_areas): print(f{year}: {area:.2f} km²)Otsu 方法的好处是不需要手动设阈值但它假设灯光值的分布是双峰的建成区和非建成区。如果某一年数据质量差分布变成单峰Otsu 会失效。这时候需要人工检查那一年的灯光图像或者改用分位数阈值。另一个进阶用法是做灯光重心的迁移分析。计算每一年的灯光重心亮度加权的地理中心然后看重心在 1992-2020 年间怎么移动。广东省的灯光重心大概率是向珠三角核心区聚集但如果发现某几年重心突然跳变可能是校正或裁剪出了问题。这个分析方法对数据质量很敏感适合用来做最终的质量把关。注意做重心分析时要用投影后的平面坐标不要用经纬度坐标否则重心位置会受纬度影响。最后说一个我自己的习惯拿到任何长时序夜间灯光数据先画三张图——全省灯光总值的时间序列、珠三角和非珠三角的对比序列、以及灯光重心的迁移轨迹。这三张图如果都合理数据基本可用如果有一张不对劲就回头查校正流程。这个习惯帮我省了很多后悔药因为很多问题在数字上不明显一画图就露馅了。希望帮到你。本文还有配套的精品资源点击获取