做城市遥感或者长时序社科研究的人大概率都绕不开一个数据源夜间灯光。我指的不是手机夜景拍照而是卫星从高轨道拍下来的地表夜间照明也就是大家常说的 NPP/VIIRS 夜间灯光数据。从城市扩张、GDP空间化、人口分布模拟到电力消耗估算、贫困指数推测这套“逐年全球夜间灯光”几乎是宏观研究里的万能底图。很多像我一样的GIS或遥感从业者从2012年第一次拿到VIIRS数据开始一路更新到2024年。在这个过程中我下载过各个机构发布的版本、踩过不同产品之间的兼容坑也在年度合成和逐日数据之间反复横跳。这篇就当成我处理这套数据的一份操作笔记从产品分类、下载路径、预处理到常见坑点排查一次说清。1. 这份数据到底记录了什么NPP/VIIRS夜间灯光的来龙去脉1.1 为什么2012年成了一个绕不开的起点在做长时序夜间灯光研究时2012年是一个非常特殊的时间节点。2011年10月Suomi NPP卫星发射升空上面搭载了可见光红外成像辐射计VIIRS。它的DNB波段专门设计为对极微弱可见光敏感能够捕捉夜间地表的人造光源。从2012年4月开始VIIRS DNB数据逐步稳定发布这正好接替了当时已经服役多年的DMSP-OLS传感器。DMSP-OLS的问题在于传感器动态范围只有6位量化城市中心很容易饱和而且同一时期有多颗卫星在飞卫星之间还存在辐射定标不一致的麻烦。相比之下VIIRS DNB采用了14位量化动态范围极大提升城市核心区的灯光不再是一团纯白夜间月光、杂散光、云层等干扰也可以通过后处理进行识别和剔除。所以在研究城市内部结构、夜间经济活动和精细尺度的地表覆盖变化时VIIRS数据几乎是必须选择的基础影像。那为什么很多人会强调“2012-2024”这个区间严格来说VIIRS可用的夜间灯光记录从2012年4月才开始。如果你要做1992到2020的长序列分析2012年和2013年正好是DMSP-OLS与VIIRS有重叠观测的两年这个重叠期是做传感器交叉定标的关键窗口。而2024年则是目前能拿到的比较新的完整年序列很多研究机构已经发布了相应年份的年度合成产品。1.2 “逐年”和“全球范围”并不像字面上那么简单刚接触这套数据的人最容易被“逐年全球”四个字迷惑。我在第一次下载时就以为每年一个全球TIFF文件拿回来直接叠栅格窗口分析就行。但实际上这个看似简单的“逐年全球”背后藏着几个容易踩雷的细节。第一同一年的数据存在不同的“自然年口径”。官方提供的年度合成产品通常按1月到12月统计但也有些早期产品会把2012年定义为2012年4月到2013年3月甚至直接用2012年4月至12月的数据做一版“首年合成”。如果你拿DMSP-OLS的2012年数据来衔接就很容易出现三个月到一年的错位。第二全球范围文件并不总是单个文件。很多网站把全球切成非洲、亚洲、美洲、欧洲等子区域或者按纬度分成若干块你需要下载后再拼接。NASA的Black Marble系列虽然也提供全球栅格但不同版本的获取方式不一样。若不注意很容易下载得到一堆没有空间参考信息的局部Tile。第三“逐年”不等于“逐个月份都很完整”。VIIRS工作以来传感器偶尔会进入校准模式或遭遇重大云覆盖有些逐日数据会有大面积缺测。年度产品是在此基础上做集合平均的所以不同年份之间的有效观测次数并不完全一致。这种由观测条件导致的年份间噪音在做短时段变化分析时是需要重点关注的。1.3 同一个“夜间灯光”为什么不同平台数值不一样我常被问到一个问题为什么在同一区域、同一年份从NASA官网下载的VIIRS夜间灯光和从EOG网站下载的VIIRS夜间灯光数值差别很大原因很简单它们不是同一种处理产线。VIIRS DNB的原始观测值会经过一系列校正包括去云、去月球光照影响、去杂散光、大气校正等。不同处理流程对噪声像元、边界像元和极低值像元的处理策略不同最终输出的合成平均值自然就有差异。另外许多官方合成产品为了压缩存储空间会先乘以一个缩放因子后再保存为整型。比如数据真实辐射值是0.37 nW/cm2/sr文件里可能存的是3或4。如果你没读元数据中的Scale Factor直接拿文件里的整型值去做统计结果就会差一个数量级。这是初学者常犯的错误后面预处理章节我会详细展开。2. 产品家族横向对比先搞清你手里的数据属于哪一类2.1 NASA Black Marble VNP46系列目前做逐年合成最方便的产品之一是NASA的Black Marble系列。它基于VIIRS DNB数据通过HDR高动态范围自适应处理剔除了月光影响、云污染和大气散射输出逐日、逐月、逐年等不同时间粒度产品。Black Marble系列的核心产品编号是VNP46。我在实际工作中常用到的是VNP46A1逐日原始大气顶部辐射值VNP46A2逐日经过云掩膜和质量筛选的夜间灯光产品VNP46A3逐月合成VNP46A4逐年合成。其中VNP46A4是全球范围的年度积温型产品像元值可以理解成全年有效夜光观测的集合统计结果。它最大的优势是不需要用户自己去逐日合成省掉了云掩膜、月光筛选、杂散光校正等一系列麻烦事。在做宏观长时序变化时直接调取A4或A3会明显提高生产效率。但要注意的是Black Marble产品虽然干净它的像元分辨率大约为500米卫星原生观测经过重采样在使用时需要将原始DNB值换算成辐射亮度。换算方法和传统VNL产品不完全一样后面我会专门说。另一个问题是A4年度数据通常要到次年年中甚至更晚才会发布如果你急着分析前一年的情况可能只能自己拿A2逐日合成或者先用月度产品临时拼接。2.2 NOAA/NCEI与EOG发布的VNL合成产品除了NASA Black Marble另一条经典产线是NOAA国家环境信息中心最早推出的VNL年度产品后续由科罗拉多矿业学院的EOG团队延续发布。这类产品在论文和行业报告中经常直接被称为“NPP/VIIRS夜间灯光年度数据”。如果你进过EOG数据下载页会看到产品名里有很强的命名规则常出现这些关键词vcm只保留无云像元未额外处理杂散光vcmsl在无云的基础上进一步剔除了杂散光影响vcm-orm额外做了离群值去除比纯vcm更干净vcmsl-orm同时做了杂散光剔除和离群值去除。这些缩写看着绕但选型逻辑很清晰。如果你的研究区在中低纬度尤其是赤道附近那么卫星成像时容易受到季节性的杂散光污染建议用vcmsl版本。如果研究区偏中高纬度、杂散光影响不显著用vcm版本问题也不大。不过为了稳妥我建议优先下载带orm后缀的版本因为数据更平滑单年抖动更小。早期NOAA/NCEI发布的版本更新到2019年前后就基本停止后续年份由EOG团队延续到目前已经发布了VNL V2、V3等新版本。不同版本之间的年际数值存在系统性偏移如果做长序列研究要么从头到尾只用同一个版本要么做好版本间的校正不建议在不同版本间混用。2.3 逐日产品与合成产品怎么选逐日产品和合成产品之间的选择本质上是在“时间分辨率”和“数据处理成本”之间做权衡。逐日数据如VNP46A1/A2能精确反映某一天的灯光状态适合研究节假日灯光变化、突发事件影响等场景。但它的问题也很明显每天全球数据量非常大且云层覆盖、月光条件、传感器扫描角变化都会带来巨大噪声。如果不做严格的质量筛选直接拿逐日数据做年平均结果很可能被几朵云或几个月夜污染。合成产品则是在时间维度上做统计平均。月合成能反映月度整体水平年合成适合趋势分析。我的建议是如果你的研究不需要逐日尺度就直接用官方月合成或年合成。不要自己从逐日数据重新做平均除非你愿意花大量时间在质量波段筛选上并且对产品的物理含义有清楚把握。从城市扩张与灯光趋势分析这个最常见的需求来看VNP46A4或EOG的年度VNL是性价比最高的选择。它们已经把云掩膜、杂散光、噪声等常见问题处理过了你拿到的数据在时间序列上相对稳定能够比较真实地反映人造光源的逐年变化。2.4 一个简单但常被忽略的产品选择流程我在实际处理时一般按下面这个流程来选产品大家可以参考先确认研究时间范围。如果研究跨1992-2020优先考虑DMSP-OLS与VIIRS的衔接问题尽量选择能提供“月度”或“逐年”且命名的VNL产品因为这类产品二次校正的参考资料最多。如果只研究2012年以后直接优先用NASA Black Marble系列或EOG的V2版本产品。我个人更倾向于Black Marble的月度/年度合成质量字段更完整。判断研究区是否位于低纬度决定要不要选vcmsl版本。如果后续要统计面积、做聚类最好一开始就用较高空间分辨率的A4产品如果只是全国尺度看总体趋势EOG年度产品也够用。下面是几个常见产品系列的对比表方便快速查阅产品系列时间粒度空间分辨率优点需要注意的点VNP46A2逐日约500米可自由合成、质量字段丰富噪声大、需自行处理云和月光VNP46A3月约500米月度对比方便质量较好单月有效观测不足时会有空值VNP46A4年约500米直接给出年合成省时发布有滞后需自行缩放EOG VNL v2cm年约500米经典传统产线易和DMSP衔接低纬度杂散光残留EOG VNL v2cmsl年约500米低纬地区更干净数据质量不太稳定个别年份有跳变3. 实际下载与快速上手GEE和官方下载两条路3.1 在GEE里快速查看2012-2024年灯光趋势如果你只是先看个趋势、做个快速出图最推荐的是Google Earth Engine。GEE里已经有专门的VIIRS Nighttime Light集合不需要在本地下载大量TIFF直接在云端做筛选和计算。我在这里给出一段可以直接运行的示例代码用来按年生成年度平均灯光影像。代码里选择了VNP46A2逐日产品并且通过日期窗口和DNB波段做简单合成。实际正式出结果时建议引入质量波段做进一步筛选。// 定义研究区以某东部城市群为例 var roi ee.Geometry.Rectangle([118.0, 30.0, 122.0, 33.0]); // 逐年合成函数 function yearlyNTL(year) { var start ee.Date.fromYMD(year, 1, 1); var end ee.Date.fromYMD(year, 12, 31); var collection ee.ImageCollection(NASA/VIIRS/BlackMarble/VNP46A2) .filterDate(start, end) .select(DNB_At_Sensor_Radiance); // 这里做简单均值再把负值归零减少噪声影响 var annual collection.mean().where(collection.mean().lt(0), 0); return annual.clip(roi).set(year, year).rename(NTL); } // 以2016年为例 var nl2016 yearlyNTL(2016); Map.centerObject(roi, 8); Map.addLayer(nl2016, {min: 0, max: 50, palette: [black, purple, yellow]}, NTL 2016); // 导出到个人Google Drive Export.image.toDrive({ image: nl2016, description: NTL_2016, scale: 500, crs: EPSG:4326, maxPixels: 1e10 });运行前注意几个点2012年只有4月以后才有VIIRS数据所以2012年的“年度合成”并不完整参考时需要留意。这段代码只做了非常简单的均值合成没有严格筛掉月光和低质量观测。在正式科研分析中最好使用官方发布的月合成或年合成产品而不是自己从逐日产品简单平均。Export导出时scale推荐设为500米和产品原生网格尽量一致避免重采样引入误差。如果追求更精确的像元面积统计建议投影到Albers或UTM后再导出。3.2 官方下载源与批量下载思路如果需要本地存档或者研究区比较固定、要放到传统桌面软件里做深度处理我通常选择官方数据仓库下载。NASA的Earthdata Search和LAADS DAAC是下载VNP46系列的主渠道。你需要先注册一个Earthdata账号然后根据产品编号、时间范围和空间范围进行查询。下载时可以直接在网页上选择子区域但如果您需要全中国或全球逐年数据更建议用批量下载工具。官方推荐的~/.netrc配置方式可以配合wget或curl实现自动下载但具体脚本需要掌握一定基础。EOG的下载页面则更直观一些。它按照年份列出了全球不同分幅文件常见文件命名类似SVDNB_npp_20150101-20151231_75N180W_vcmslcfg.tif。下载时需要注意每个年份对应文件名里的日期范围不要只改年份数字而不改日期。如果你想批量下载可以复制某一年的文件链接把年份和日期字段批量替换成其他年份再利用wget工具连续执行。举例如下# 示例思路实际链接以页面提供为准 wget -c https://eogdata.mines.edu/xxx/viirs_ntl/annual/v2/2016/SVDNB_npp_20160101-20161231.75N180W.vcmslcfg.tif在批量下载时我还会加两个习惯文件下载到一半断掉用wget -c继续续传。下载完成后马上比对文件大小凡是明显小于平均值的优先删除重新下。EOG页面里的分幅往往是全球切成小块因此如果你只需要某个区域先拼后裁可能更稳妥。部分平台也会提供全图tif可以直接在页面上筛选具备“global”字样的文件。3.3 拿到文件后先核对哪几件事下载完影像别急着拿去算先用GIS软件打开确认下面几个信息。第一是像元值范围。由于辐射亮度通常很低许多产品会把数值扩大十倍或一百倍保存在整型栅格里。如果文件里的最大值是几十或几百说明可能是带缩放因子的存档。此时你需要去产品说明里查找Scale Factor等后续处理时把数值还原回物理辐射单位。第二是空间参考。有些发布源提供的GeoTIFF自带WGS84经纬度坐标有的则没有完整投影信息。如果打开后发现位置偏移到海里多半是元数据里缺少投影定义需要手工指定。第三是无效值。不同产品对背景值和填充值的定义不同有的用0表示无灯光有的用-999或65535等填充值。如果不把无效值先提取成NoData后面的统计计算会把大片海洋当成“零值”区域直接拉低城市灯光指数。4. 预处理和指数计算从栅格到能写进论文的指标4.1 负值、零值和单位换算最常见的第一道坎VIIRS DNB原始信号中包含一定的背景噪声在低照度区域经过校正后的辐射值可能变成负值。年度合成产品也可能残留负值直接进行求均值、求和会影响结果。我习惯的第一步都是统一做“非负处理”。在ArcGIS栅格计算器或QGIS里可以这样写ntl 0 → 0如果是在Python的rasterio或numpy里可以直接用np.where(array 0, 0, array)。做完负值归零后还要确认数据单位。以EOG年度VNL产品为例很多早期产品提供的DN值乘以0.1才是真实的辐射值单位是nW/cm2/sr。但这个系数不是所有版本都相同。如果你想比较不同年份或不同产品的绝对亮度必须先统一转换到物理辐射量纲。一个常见误区是很多人直接对原始DN值取对数。如果DN里存在0或负值取对数后会出现无穷大或空值。推荐做法是先加一个偏移量比如log(DN 1)再去做变换。想更平滑地压缩动态范围也可以尝试arcsinh变换这种方法在夜间灯光研究里越来越常用。4.2 投影、裁剪、年度对比的标准化全球范围的VIIRS夜间灯光产品默认多为经纬度网格每个像元约0.0041667度差不多是500米。做小范围城市分析时直接按经纬度计算面积会产生误差所以需要投影到等积投影上再统计。国内做区域分析时比较常见的做法是投影到Albers等积投影或UTM分带投影。对于500米像元的夜间灯光我通常使用UTM投影保证面积精度。如果你处理的是全国大范围数据用Albers等积投影会更好。投影后还要注意裁剪边界。如果裁剪矢量边界比较精确裁出来的图像可能包含边界外明显过亮的区域例如邻省城市灯光被切割进来一半也有可能把城市边缘低值像元切掉。我习惯在裁剪时做一次1公里左右的缓冲区外扩在统计时再利用矢量范围做严格裁剪避免边界效应影响指标。跨年对比时标准化处理同样关键。每一年的影像不能单独分别做归一化因为这样会让所有年份都变成0到1失去真实变化信息。正确的做法是先确定整个时间范围内的最大值或99.9%分位数作为统一归一化基准然后所有年份都用同一个基准去缩放这样年份间的差异才能保留下来。4.3 与DMSP-OLS数据衔接时怎么校正如果你想从1992年一直分析到2024年就避不开DMSP-OLS和NPP/VIIRS之间的衔接问题。DMSP-OLS的DN值范围是0到63而VIIRS夜间灯光是连续辐射值两者不仅量纲不同饱和特性也不同。常见的校正思路是利用2012年至2013年的重叠期做统计回归。在重叠年份找到同一区域内的DMSP像元和VIIRS像元进行配对后拟合关系总体来看DMSP与VIIRS之间存在较强的非线性关系幂函数或分段线性都是常用的拟合形式。不同研究区的拟合参数差别较大所以不能直接套用论文里的公式最好自己在研究区内采样拟合一次。我个人的建议是不要试图把VIIRS低值区完全映射到DMSP的低值区因为VIIRS对微弱灯光的灵敏度远高于DMSP。在低照度区域用简单线性回归很容易造成大面积虚高。更稳定的做法是先设定一个低亮度阈值只对城市建成区等较高亮度区做转换或者直接用“灯光面积占比”这类相对指标做长时序分析而不是比绝对亮度。4.4 城市尺度灯光指数TNL与CNLI的计算示例在实际研究中最常用的几个指标包括单年灯光总量、灯光面积、平均灯光强度和综合灯光指数。下面是它们的常见计算方式TNL单一年份研究区内所有像元灯光强度之和可以理解为区域总灯光亮度。计算公式为TNL ΣDN_i。灯光面积超过某一亮度阈值的像元数量乘以单个像元面积。阈值常用经验值确定也可根据DMSP/VIIRS的夜间灯光背景噪声确定。CNLI综合灯光指数通常定义为平均相对灯光强度和灯光面积占比的乘积。以QGIS栅格计算器为例计算某年研究区内TNL时可以把裁剪后的灯光栅格乘上单个像元面积后再求和。如果数据已经是辐射亮度值那TNL的单位就是“辐射亮度×面积”代表该区域夜间总辐射通量。如果计算CNLI可以采用如下思路I (ΣDN_i) / (N_light × DN_max) S N_light / N_total CNLI I × S这里N_light是亮灯像元数N_total是研究区总像元数DN_max是归一化基准。更严谨的做法是DN_max取多年最大值的99.9%分位数而不是直接取文件中可能出现的异常最大值避免个别噪声像元把指数压得过低。用这些指数可以很清晰地观察到城市发展过程。例如某个老牌城市群在2012年到2019年间TNL持续上升但CNLI上升速度逐渐放缓这说明扩张主要发生在外围低密度地区2020年以后核心区灯光亮度趋于稳定亮灯面积继续扩展则反映城市空间结构从中心集聚走向多中心扩散。这些判断如果只用单一指标很难完全显现。5. 我踩过的坑和排查思路5.1 常见问题速查表表格是我整理了几年使用中频率最高的问题按现象、原因和处理方式列出来遇到问题可以直接查问题现象最常见原因处理方式整片影像看起来全黑1. 没按实际拉伸范围显示2. 没乘缩放因子3. 背景用0填充先看像元统计确认最大值和最小值再做2%线性拉伸城市中心亮度骤降或出现黑洞没有正确筛选杂散光或产品版本不匹配低纬地区低纬用vcmsl版本处理逐日数据时加入质量波段的杂散光掩膜跨年份数据突跳严重不同年份使用了不同版本产品或2012年首年日期口径不同统一产品版本核对各年份文件名内的日期范围负值特别多DNB极低照度区的背景噪声没被剔除负值归零或者在做对数变换时加偏移量海岸线周围有零星亮斑船、海上气井、油气燃烧都会形成夜光先通过掩膜文件剔除水体或通过辅助数据识别统计面积和实际城市面积差异巨大投影选择不当直接用经纬度角度计算面积投影到Albers或UTM后再做面积计算分辨率不是500米而是约463米全球经纬度网格在低纬度1个角秒约463米按需重采样并注意重采样方法对总量统计影响5.2 杂散光为什么某些低纬城市在“消失”杂散光问题是VIIRS DNB数据处理中比较隐蔽的坑。VIIRS传感器在昼夜明暗界线附近时太阳光可能通过光学系统内部散射进入探测器导致DNB信号被污染。这种影响不是均匀分布在全球的在特定季节和特定卫星轨道条件下低纬度地区受影响尤为突出。我第一次处理东南亚区域月度数据时就发现某个月度合成影像里城市灯光变得断断续续部分城市核心甚至出现亮度大幅下降。排查了很久才发现问题出在杂散光处理上。如果使用普通vcm版本个别低纬地区的月度效果确实可能很差这时需要转用vcmsl版本。VCMSL中的“sl”就是杂散光处理标记它专门针对低纬度的受污染像元做了校正。对大多数中纬度用户来说这个坑影响有限。但如果你研究赤道附近的国家比如东南亚、南美北部或非洲中部建议直接选择vcmsl版本进行常规处理不要等到发现问题再返工。5.3 版本断档和滞后更新怎么应对做逐年长时间序列的人最怕的不是数据本身噪声而是不同年份的数据来自不同处理版本。官方会在积累更多校准数据后发布新版本并且通常只保留最新版本或提供“Version History”。一旦你前期使用了旧版本后期为了获得最新年份被迫切换到新版本就会面临一个巨大的断层。解决这个问题我只能提供几条“血泪经验”第一项目一开始就锁定一个主产品主版本。不要因为看到某一年出了新版本就去更新单独一年的数据。版本更新带来的系统性差异通常比真实灯光变化还要大。第二如果必须切换版本要留出重叠期数据来做校正。比如旧版本只到2020年新版本从2015年开始那么你可以用2015-2020年这段重叠期做线性调整把新版本回归到旧版本的统计口径下。第三产品发布的滞后问题也需要提前规划。很多年度合成产品在次年年中才会更新如果你想在2024年初就用2023年的灯光数据很可能只能得到不完整版本。此时建议使用截至某个日期的部分月度合成做临时估算等年度产品正式发布后再做最终修正。5.4 几个能让你少走弯路的操作习惯这里分享几个我在长期处理过程中沉淀下来的小习惯不一定属于高深技术但确实能提升效率。我习惯对每一份下载的原始数据都保留一份“零处理”存档。不管后续怎么投影、裁剪、去负值、重采样原始文件都不会再改动。数据处理链条长了以后如果忘了哪一步导致结果异常还可以随时回到起点重新来一遍不用再重新下载。在处理灯光数据时维护一个“数据说明笔记”也很有用。笔记里记录产品名称、网站版本号、文件名日期范围、有效值范围、坐标系统、缩放系数等甚至可以把下载时间一并记上。版本追踪主要靠这些细节不然两年后回来看自己的数据目录很可能会搞不清每一版到底对应什么产品。最后一个建议做变化检测时不要依赖单年的异常值。比如某年城市核心区的灯光密度因为云覆盖或数据缺失出现明显下跌直接纳入趋势分析会影响整体判断。我通常会对逐年序列做一次简单的异常值检查比如计算每年的区域总亮度和平均值标记出偏离历史趋势超过3个标准差的年份再去检查原始影像原因。在我自己的项目里这套“下载-检查-统一-处理-验证”的流程已经用了很多年2012年到2024年的夜间灯光数据我前前后后处理过好几轮。每次换新版本或新研究区时还是会碰到一些小问题但只要产品选型正确、预处理流程规范、年份口径核对清楚这套数据依然是大规模宏观研究里非常值得信赖的空间数据源。如果你刚开始接触NPP/VIIRS夜间灯光数据我的建议很直接不要贪多求全先选定一条产品线把一个年份、一个小区域完整跑通再逐步拓展到全时间序列。这样既不会在数据下载上浪费时间也能在一次次实际运行中更快找到适合自己的处理套路。
NPP/VIIRS夜间灯光数据全解析：从产品选型到预处理实战
做城市遥感或者长时序社科研究的人大概率都绕不开一个数据源夜间灯光。我指的不是手机夜景拍照而是卫星从高轨道拍下来的地表夜间照明也就是大家常说的 NPP/VIIRS 夜间灯光数据。从城市扩张、GDP空间化、人口分布模拟到电力消耗估算、贫困指数推测这套“逐年全球夜间灯光”几乎是宏观研究里的万能底图。很多像我一样的GIS或遥感从业者从2012年第一次拿到VIIRS数据开始一路更新到2024年。在这个过程中我下载过各个机构发布的版本、踩过不同产品之间的兼容坑也在年度合成和逐日数据之间反复横跳。这篇就当成我处理这套数据的一份操作笔记从产品分类、下载路径、预处理到常见坑点排查一次说清。1. 这份数据到底记录了什么NPP/VIIRS夜间灯光的来龙去脉1.1 为什么2012年成了一个绕不开的起点在做长时序夜间灯光研究时2012年是一个非常特殊的时间节点。2011年10月Suomi NPP卫星发射升空上面搭载了可见光红外成像辐射计VIIRS。它的DNB波段专门设计为对极微弱可见光敏感能够捕捉夜间地表的人造光源。从2012年4月开始VIIRS DNB数据逐步稳定发布这正好接替了当时已经服役多年的DMSP-OLS传感器。DMSP-OLS的问题在于传感器动态范围只有6位量化城市中心很容易饱和而且同一时期有多颗卫星在飞卫星之间还存在辐射定标不一致的麻烦。相比之下VIIRS DNB采用了14位量化动态范围极大提升城市核心区的灯光不再是一团纯白夜间月光、杂散光、云层等干扰也可以通过后处理进行识别和剔除。所以在研究城市内部结构、夜间经济活动和精细尺度的地表覆盖变化时VIIRS数据几乎是必须选择的基础影像。那为什么很多人会强调“2012-2024”这个区间严格来说VIIRS可用的夜间灯光记录从2012年4月才开始。如果你要做1992到2020的长序列分析2012年和2013年正好是DMSP-OLS与VIIRS有重叠观测的两年这个重叠期是做传感器交叉定标的关键窗口。而2024年则是目前能拿到的比较新的完整年序列很多研究机构已经发布了相应年份的年度合成产品。1.2 “逐年”和“全球范围”并不像字面上那么简单刚接触这套数据的人最容易被“逐年全球”四个字迷惑。我在第一次下载时就以为每年一个全球TIFF文件拿回来直接叠栅格窗口分析就行。但实际上这个看似简单的“逐年全球”背后藏着几个容易踩雷的细节。第一同一年的数据存在不同的“自然年口径”。官方提供的年度合成产品通常按1月到12月统计但也有些早期产品会把2012年定义为2012年4月到2013年3月甚至直接用2012年4月至12月的数据做一版“首年合成”。如果你拿DMSP-OLS的2012年数据来衔接就很容易出现三个月到一年的错位。第二全球范围文件并不总是单个文件。很多网站把全球切成非洲、亚洲、美洲、欧洲等子区域或者按纬度分成若干块你需要下载后再拼接。NASA的Black Marble系列虽然也提供全球栅格但不同版本的获取方式不一样。若不注意很容易下载得到一堆没有空间参考信息的局部Tile。第三“逐年”不等于“逐个月份都很完整”。VIIRS工作以来传感器偶尔会进入校准模式或遭遇重大云覆盖有些逐日数据会有大面积缺测。年度产品是在此基础上做集合平均的所以不同年份之间的有效观测次数并不完全一致。这种由观测条件导致的年份间噪音在做短时段变化分析时是需要重点关注的。1.3 同一个“夜间灯光”为什么不同平台数值不一样我常被问到一个问题为什么在同一区域、同一年份从NASA官网下载的VIIRS夜间灯光和从EOG网站下载的VIIRS夜间灯光数值差别很大原因很简单它们不是同一种处理产线。VIIRS DNB的原始观测值会经过一系列校正包括去云、去月球光照影响、去杂散光、大气校正等。不同处理流程对噪声像元、边界像元和极低值像元的处理策略不同最终输出的合成平均值自然就有差异。另外许多官方合成产品为了压缩存储空间会先乘以一个缩放因子后再保存为整型。比如数据真实辐射值是0.37 nW/cm2/sr文件里可能存的是3或4。如果你没读元数据中的Scale Factor直接拿文件里的整型值去做统计结果就会差一个数量级。这是初学者常犯的错误后面预处理章节我会详细展开。2. 产品家族横向对比先搞清你手里的数据属于哪一类2.1 NASA Black Marble VNP46系列目前做逐年合成最方便的产品之一是NASA的Black Marble系列。它基于VIIRS DNB数据通过HDR高动态范围自适应处理剔除了月光影响、云污染和大气散射输出逐日、逐月、逐年等不同时间粒度产品。Black Marble系列的核心产品编号是VNP46。我在实际工作中常用到的是VNP46A1逐日原始大气顶部辐射值VNP46A2逐日经过云掩膜和质量筛选的夜间灯光产品VNP46A3逐月合成VNP46A4逐年合成。其中VNP46A4是全球范围的年度积温型产品像元值可以理解成全年有效夜光观测的集合统计结果。它最大的优势是不需要用户自己去逐日合成省掉了云掩膜、月光筛选、杂散光校正等一系列麻烦事。在做宏观长时序变化时直接调取A4或A3会明显提高生产效率。但要注意的是Black Marble产品虽然干净它的像元分辨率大约为500米卫星原生观测经过重采样在使用时需要将原始DNB值换算成辐射亮度。换算方法和传统VNL产品不完全一样后面我会专门说。另一个问题是A4年度数据通常要到次年年中甚至更晚才会发布如果你急着分析前一年的情况可能只能自己拿A2逐日合成或者先用月度产品临时拼接。2.2 NOAA/NCEI与EOG发布的VNL合成产品除了NASA Black Marble另一条经典产线是NOAA国家环境信息中心最早推出的VNL年度产品后续由科罗拉多矿业学院的EOG团队延续发布。这类产品在论文和行业报告中经常直接被称为“NPP/VIIRS夜间灯光年度数据”。如果你进过EOG数据下载页会看到产品名里有很强的命名规则常出现这些关键词vcm只保留无云像元未额外处理杂散光vcmsl在无云的基础上进一步剔除了杂散光影响vcm-orm额外做了离群值去除比纯vcm更干净vcmsl-orm同时做了杂散光剔除和离群值去除。这些缩写看着绕但选型逻辑很清晰。如果你的研究区在中低纬度尤其是赤道附近那么卫星成像时容易受到季节性的杂散光污染建议用vcmsl版本。如果研究区偏中高纬度、杂散光影响不显著用vcm版本问题也不大。不过为了稳妥我建议优先下载带orm后缀的版本因为数据更平滑单年抖动更小。早期NOAA/NCEI发布的版本更新到2019年前后就基本停止后续年份由EOG团队延续到目前已经发布了VNL V2、V3等新版本。不同版本之间的年际数值存在系统性偏移如果做长序列研究要么从头到尾只用同一个版本要么做好版本间的校正不建议在不同版本间混用。2.3 逐日产品与合成产品怎么选逐日产品和合成产品之间的选择本质上是在“时间分辨率”和“数据处理成本”之间做权衡。逐日数据如VNP46A1/A2能精确反映某一天的灯光状态适合研究节假日灯光变化、突发事件影响等场景。但它的问题也很明显每天全球数据量非常大且云层覆盖、月光条件、传感器扫描角变化都会带来巨大噪声。如果不做严格的质量筛选直接拿逐日数据做年平均结果很可能被几朵云或几个月夜污染。合成产品则是在时间维度上做统计平均。月合成能反映月度整体水平年合成适合趋势分析。我的建议是如果你的研究不需要逐日尺度就直接用官方月合成或年合成。不要自己从逐日数据重新做平均除非你愿意花大量时间在质量波段筛选上并且对产品的物理含义有清楚把握。从城市扩张与灯光趋势分析这个最常见的需求来看VNP46A4或EOG的年度VNL是性价比最高的选择。它们已经把云掩膜、杂散光、噪声等常见问题处理过了你拿到的数据在时间序列上相对稳定能够比较真实地反映人造光源的逐年变化。2.4 一个简单但常被忽略的产品选择流程我在实际处理时一般按下面这个流程来选产品大家可以参考先确认研究时间范围。如果研究跨1992-2020优先考虑DMSP-OLS与VIIRS的衔接问题尽量选择能提供“月度”或“逐年”且命名的VNL产品因为这类产品二次校正的参考资料最多。如果只研究2012年以后直接优先用NASA Black Marble系列或EOG的V2版本产品。我个人更倾向于Black Marble的月度/年度合成质量字段更完整。判断研究区是否位于低纬度决定要不要选vcmsl版本。如果后续要统计面积、做聚类最好一开始就用较高空间分辨率的A4产品如果只是全国尺度看总体趋势EOG年度产品也够用。下面是几个常见产品系列的对比表方便快速查阅产品系列时间粒度空间分辨率优点需要注意的点VNP46A2逐日约500米可自由合成、质量字段丰富噪声大、需自行处理云和月光VNP46A3月约500米月度对比方便质量较好单月有效观测不足时会有空值VNP46A4年约500米直接给出年合成省时发布有滞后需自行缩放EOG VNL v2cm年约500米经典传统产线易和DMSP衔接低纬度杂散光残留EOG VNL v2cmsl年约500米低纬地区更干净数据质量不太稳定个别年份有跳变3. 实际下载与快速上手GEE和官方下载两条路3.1 在GEE里快速查看2012-2024年灯光趋势如果你只是先看个趋势、做个快速出图最推荐的是Google Earth Engine。GEE里已经有专门的VIIRS Nighttime Light集合不需要在本地下载大量TIFF直接在云端做筛选和计算。我在这里给出一段可以直接运行的示例代码用来按年生成年度平均灯光影像。代码里选择了VNP46A2逐日产品并且通过日期窗口和DNB波段做简单合成。实际正式出结果时建议引入质量波段做进一步筛选。// 定义研究区以某东部城市群为例 var roi ee.Geometry.Rectangle([118.0, 30.0, 122.0, 33.0]); // 逐年合成函数 function yearlyNTL(year) { var start ee.Date.fromYMD(year, 1, 1); var end ee.Date.fromYMD(year, 12, 31); var collection ee.ImageCollection(NASA/VIIRS/BlackMarble/VNP46A2) .filterDate(start, end) .select(DNB_At_Sensor_Radiance); // 这里做简单均值再把负值归零减少噪声影响 var annual collection.mean().where(collection.mean().lt(0), 0); return annual.clip(roi).set(year, year).rename(NTL); } // 以2016年为例 var nl2016 yearlyNTL(2016); Map.centerObject(roi, 8); Map.addLayer(nl2016, {min: 0, max: 50, palette: [black, purple, yellow]}, NTL 2016); // 导出到个人Google Drive Export.image.toDrive({ image: nl2016, description: NTL_2016, scale: 500, crs: EPSG:4326, maxPixels: 1e10 });运行前注意几个点2012年只有4月以后才有VIIRS数据所以2012年的“年度合成”并不完整参考时需要留意。这段代码只做了非常简单的均值合成没有严格筛掉月光和低质量观测。在正式科研分析中最好使用官方发布的月合成或年合成产品而不是自己从逐日产品简单平均。Export导出时scale推荐设为500米和产品原生网格尽量一致避免重采样引入误差。如果追求更精确的像元面积统计建议投影到Albers或UTM后再导出。3.2 官方下载源与批量下载思路如果需要本地存档或者研究区比较固定、要放到传统桌面软件里做深度处理我通常选择官方数据仓库下载。NASA的Earthdata Search和LAADS DAAC是下载VNP46系列的主渠道。你需要先注册一个Earthdata账号然后根据产品编号、时间范围和空间范围进行查询。下载时可以直接在网页上选择子区域但如果您需要全中国或全球逐年数据更建议用批量下载工具。官方推荐的~/.netrc配置方式可以配合wget或curl实现自动下载但具体脚本需要掌握一定基础。EOG的下载页面则更直观一些。它按照年份列出了全球不同分幅文件常见文件命名类似SVDNB_npp_20150101-20151231_75N180W_vcmslcfg.tif。下载时需要注意每个年份对应文件名里的日期范围不要只改年份数字而不改日期。如果你想批量下载可以复制某一年的文件链接把年份和日期字段批量替换成其他年份再利用wget工具连续执行。举例如下# 示例思路实际链接以页面提供为准 wget -c https://eogdata.mines.edu/xxx/viirs_ntl/annual/v2/2016/SVDNB_npp_20160101-20161231.75N180W.vcmslcfg.tif在批量下载时我还会加两个习惯文件下载到一半断掉用wget -c继续续传。下载完成后马上比对文件大小凡是明显小于平均值的优先删除重新下。EOG页面里的分幅往往是全球切成小块因此如果你只需要某个区域先拼后裁可能更稳妥。部分平台也会提供全图tif可以直接在页面上筛选具备“global”字样的文件。3.3 拿到文件后先核对哪几件事下载完影像别急着拿去算先用GIS软件打开确认下面几个信息。第一是像元值范围。由于辐射亮度通常很低许多产品会把数值扩大十倍或一百倍保存在整型栅格里。如果文件里的最大值是几十或几百说明可能是带缩放因子的存档。此时你需要去产品说明里查找Scale Factor等后续处理时把数值还原回物理辐射单位。第二是空间参考。有些发布源提供的GeoTIFF自带WGS84经纬度坐标有的则没有完整投影信息。如果打开后发现位置偏移到海里多半是元数据里缺少投影定义需要手工指定。第三是无效值。不同产品对背景值和填充值的定义不同有的用0表示无灯光有的用-999或65535等填充值。如果不把无效值先提取成NoData后面的统计计算会把大片海洋当成“零值”区域直接拉低城市灯光指数。4. 预处理和指数计算从栅格到能写进论文的指标4.1 负值、零值和单位换算最常见的第一道坎VIIRS DNB原始信号中包含一定的背景噪声在低照度区域经过校正后的辐射值可能变成负值。年度合成产品也可能残留负值直接进行求均值、求和会影响结果。我习惯的第一步都是统一做“非负处理”。在ArcGIS栅格计算器或QGIS里可以这样写ntl 0 → 0如果是在Python的rasterio或numpy里可以直接用np.where(array 0, 0, array)。做完负值归零后还要确认数据单位。以EOG年度VNL产品为例很多早期产品提供的DN值乘以0.1才是真实的辐射值单位是nW/cm2/sr。但这个系数不是所有版本都相同。如果你想比较不同年份或不同产品的绝对亮度必须先统一转换到物理辐射量纲。一个常见误区是很多人直接对原始DN值取对数。如果DN里存在0或负值取对数后会出现无穷大或空值。推荐做法是先加一个偏移量比如log(DN 1)再去做变换。想更平滑地压缩动态范围也可以尝试arcsinh变换这种方法在夜间灯光研究里越来越常用。4.2 投影、裁剪、年度对比的标准化全球范围的VIIRS夜间灯光产品默认多为经纬度网格每个像元约0.0041667度差不多是500米。做小范围城市分析时直接按经纬度计算面积会产生误差所以需要投影到等积投影上再统计。国内做区域分析时比较常见的做法是投影到Albers等积投影或UTM分带投影。对于500米像元的夜间灯光我通常使用UTM投影保证面积精度。如果你处理的是全国大范围数据用Albers等积投影会更好。投影后还要注意裁剪边界。如果裁剪矢量边界比较精确裁出来的图像可能包含边界外明显过亮的区域例如邻省城市灯光被切割进来一半也有可能把城市边缘低值像元切掉。我习惯在裁剪时做一次1公里左右的缓冲区外扩在统计时再利用矢量范围做严格裁剪避免边界效应影响指标。跨年对比时标准化处理同样关键。每一年的影像不能单独分别做归一化因为这样会让所有年份都变成0到1失去真实变化信息。正确的做法是先确定整个时间范围内的最大值或99.9%分位数作为统一归一化基准然后所有年份都用同一个基准去缩放这样年份间的差异才能保留下来。4.3 与DMSP-OLS数据衔接时怎么校正如果你想从1992年一直分析到2024年就避不开DMSP-OLS和NPP/VIIRS之间的衔接问题。DMSP-OLS的DN值范围是0到63而VIIRS夜间灯光是连续辐射值两者不仅量纲不同饱和特性也不同。常见的校正思路是利用2012年至2013年的重叠期做统计回归。在重叠年份找到同一区域内的DMSP像元和VIIRS像元进行配对后拟合关系总体来看DMSP与VIIRS之间存在较强的非线性关系幂函数或分段线性都是常用的拟合形式。不同研究区的拟合参数差别较大所以不能直接套用论文里的公式最好自己在研究区内采样拟合一次。我个人的建议是不要试图把VIIRS低值区完全映射到DMSP的低值区因为VIIRS对微弱灯光的灵敏度远高于DMSP。在低照度区域用简单线性回归很容易造成大面积虚高。更稳定的做法是先设定一个低亮度阈值只对城市建成区等较高亮度区做转换或者直接用“灯光面积占比”这类相对指标做长时序分析而不是比绝对亮度。4.4 城市尺度灯光指数TNL与CNLI的计算示例在实际研究中最常用的几个指标包括单年灯光总量、灯光面积、平均灯光强度和综合灯光指数。下面是它们的常见计算方式TNL单一年份研究区内所有像元灯光强度之和可以理解为区域总灯光亮度。计算公式为TNL ΣDN_i。灯光面积超过某一亮度阈值的像元数量乘以单个像元面积。阈值常用经验值确定也可根据DMSP/VIIRS的夜间灯光背景噪声确定。CNLI综合灯光指数通常定义为平均相对灯光强度和灯光面积占比的乘积。以QGIS栅格计算器为例计算某年研究区内TNL时可以把裁剪后的灯光栅格乘上单个像元面积后再求和。如果数据已经是辐射亮度值那TNL的单位就是“辐射亮度×面积”代表该区域夜间总辐射通量。如果计算CNLI可以采用如下思路I (ΣDN_i) / (N_light × DN_max) S N_light / N_total CNLI I × S这里N_light是亮灯像元数N_total是研究区总像元数DN_max是归一化基准。更严谨的做法是DN_max取多年最大值的99.9%分位数而不是直接取文件中可能出现的异常最大值避免个别噪声像元把指数压得过低。用这些指数可以很清晰地观察到城市发展过程。例如某个老牌城市群在2012年到2019年间TNL持续上升但CNLI上升速度逐渐放缓这说明扩张主要发生在外围低密度地区2020年以后核心区灯光亮度趋于稳定亮灯面积继续扩展则反映城市空间结构从中心集聚走向多中心扩散。这些判断如果只用单一指标很难完全显现。5. 我踩过的坑和排查思路5.1 常见问题速查表表格是我整理了几年使用中频率最高的问题按现象、原因和处理方式列出来遇到问题可以直接查问题现象最常见原因处理方式整片影像看起来全黑1. 没按实际拉伸范围显示2. 没乘缩放因子3. 背景用0填充先看像元统计确认最大值和最小值再做2%线性拉伸城市中心亮度骤降或出现黑洞没有正确筛选杂散光或产品版本不匹配低纬地区低纬用vcmsl版本处理逐日数据时加入质量波段的杂散光掩膜跨年份数据突跳严重不同年份使用了不同版本产品或2012年首年日期口径不同统一产品版本核对各年份文件名内的日期范围负值特别多DNB极低照度区的背景噪声没被剔除负值归零或者在做对数变换时加偏移量海岸线周围有零星亮斑船、海上气井、油气燃烧都会形成夜光先通过掩膜文件剔除水体或通过辅助数据识别统计面积和实际城市面积差异巨大投影选择不当直接用经纬度角度计算面积投影到Albers或UTM后再做面积计算分辨率不是500米而是约463米全球经纬度网格在低纬度1个角秒约463米按需重采样并注意重采样方法对总量统计影响5.2 杂散光为什么某些低纬城市在“消失”杂散光问题是VIIRS DNB数据处理中比较隐蔽的坑。VIIRS传感器在昼夜明暗界线附近时太阳光可能通过光学系统内部散射进入探测器导致DNB信号被污染。这种影响不是均匀分布在全球的在特定季节和特定卫星轨道条件下低纬度地区受影响尤为突出。我第一次处理东南亚区域月度数据时就发现某个月度合成影像里城市灯光变得断断续续部分城市核心甚至出现亮度大幅下降。排查了很久才发现问题出在杂散光处理上。如果使用普通vcm版本个别低纬地区的月度效果确实可能很差这时需要转用vcmsl版本。VCMSL中的“sl”就是杂散光处理标记它专门针对低纬度的受污染像元做了校正。对大多数中纬度用户来说这个坑影响有限。但如果你研究赤道附近的国家比如东南亚、南美北部或非洲中部建议直接选择vcmsl版本进行常规处理不要等到发现问题再返工。5.3 版本断档和滞后更新怎么应对做逐年长时间序列的人最怕的不是数据本身噪声而是不同年份的数据来自不同处理版本。官方会在积累更多校准数据后发布新版本并且通常只保留最新版本或提供“Version History”。一旦你前期使用了旧版本后期为了获得最新年份被迫切换到新版本就会面临一个巨大的断层。解决这个问题我只能提供几条“血泪经验”第一项目一开始就锁定一个主产品主版本。不要因为看到某一年出了新版本就去更新单独一年的数据。版本更新带来的系统性差异通常比真实灯光变化还要大。第二如果必须切换版本要留出重叠期数据来做校正。比如旧版本只到2020年新版本从2015年开始那么你可以用2015-2020年这段重叠期做线性调整把新版本回归到旧版本的统计口径下。第三产品发布的滞后问题也需要提前规划。很多年度合成产品在次年年中才会更新如果你想在2024年初就用2023年的灯光数据很可能只能得到不完整版本。此时建议使用截至某个日期的部分月度合成做临时估算等年度产品正式发布后再做最终修正。5.4 几个能让你少走弯路的操作习惯这里分享几个我在长期处理过程中沉淀下来的小习惯不一定属于高深技术但确实能提升效率。我习惯对每一份下载的原始数据都保留一份“零处理”存档。不管后续怎么投影、裁剪、去负值、重采样原始文件都不会再改动。数据处理链条长了以后如果忘了哪一步导致结果异常还可以随时回到起点重新来一遍不用再重新下载。在处理灯光数据时维护一个“数据说明笔记”也很有用。笔记里记录产品名称、网站版本号、文件名日期范围、有效值范围、坐标系统、缩放系数等甚至可以把下载时间一并记上。版本追踪主要靠这些细节不然两年后回来看自己的数据目录很可能会搞不清每一版到底对应什么产品。最后一个建议做变化检测时不要依赖单年的异常值。比如某年城市核心区的灯光密度因为云覆盖或数据缺失出现明显下跌直接纳入趋势分析会影响整体判断。我通常会对逐年序列做一次简单的异常值检查比如计算每年的区域总亮度和平均值标记出偏离历史趋势超过3个标准差的年份再去检查原始影像原因。在我自己的项目里这套“下载-检查-统一-处理-验证”的流程已经用了很多年2012年到2024年的夜间灯光数据我前前后后处理过好几轮。每次换新版本或新研究区时还是会碰到一些小问题但只要产品选型正确、预处理流程规范、年份口径核对清楚这套数据依然是大规模宏观研究里非常值得信赖的空间数据源。如果你刚开始接触NPP/VIIRS夜间灯光数据我的建议很直接不要贪多求全先选定一条产品线把一个年份、一个小区域完整跑通再逐步拓展到全时间序列。这样既不会在数据下载上浪费时间也能在一次次实际运行中更快找到适合自己的处理套路。
