企业官网建设流程全解析

简介：本资源为2012—2020年NPP-VIIRS夜间灯光遥感数据（VNL2）的预处理成品，面向从事城市扩张、经济活动评估、碳排放与人口空间化等研究的高校师生及科研人员，可省去自行下载、裁剪与清洗的繁琐环节。数据分辨率为500m，已用WGS84中国矢量边界裁剪出全国范围，并转换为Asia_Lambert_Conformal_Conic平面投影；针对负值及gas flares导致的极端异常值，参考Kang Wu等（2019）以472.86为全国像元辐射阈值进行剔除，可直接用于时序分析与建模。压缩包共36个文件，以tif栅格为主体，配套tfw坐标文件、xml元数据与ovr金字塔文件，整体约176.73MB，按年份命名便于逐年调用。目前已有3505人学习下载，适合需要现成、规范夜光数据支撑论文或课题的研究者。

1. 预处理后的2012-2020年NPP-VIIRS夜光遥感数据：从原始栅格到能直接跑模型的干净面板

手里拿到一份“预处理后的2012-2020年NPP-VIIRS夜光遥感数据”，很多人第一反应是直接丢进模型跑回归，结果发现城市边缘的亮度值忽高忽低，同一座城市2013年和2018年的灯光面积对不上，甚至有些月份整片区域是空的。这不是数据坏了，而是夜光遥感从传感器原始DN值到可比较的时序面板，中间隔着一整套预处理链路：去背景噪声、掩膜火光、剔除月光污染、跨传感器校正、年际合成、重采样对齐。这套链路没走通，后面所有分析都是沙上建塔。这篇笔记面向需要把夜光数据用起来的人——做城市扩张监测的、算区域经济代理变量的、跑碳排放与灯光相关性模型的——把2012到2020年这九年NPP-VIIRS数据的预处理逻辑、可复现步骤和踩坑点讲清楚，让你拿到一份能直接进模型的面板数据，而不是一堆需要反复返工的半成品。

2. 夜光遥感预处理到底在预处理什么：从DN值到辐射定标的四层过滤

2.1 原始NPP-VIIRS月度数据里混进了哪些非灯光信号

NPP-VIIRS搭载的VIIRS传感器有日间/夜间波段（DNB），其月度合成产品（常见为VCM或VCMSLCFG系列）在生成时已经做了一部分筛选，但直接拿来用仍然会碰到四类污染。第一类是火光，包括天然气火炬、森林火灾、工业燃烧，这些在夜间亮度极高，如果不剔除，会把一个油田区域误判成超大城市。第二类是月光，满月前后的夜间云层反射月光，导致无灯光区域出现大面积低亮度值，尤其在热带和干旱区明显。第三类是极光和高能粒子噪声，高纬度地区常见。第四类是背景噪声，传感器暗电流和大气散射带来的本底值，通常在0.5 nW/cm²/sr以下，但如果不设阈值，会让人误以为荒漠里也有经济活动。

常见做法是先用月度产品自带的QF（Quality Flag）波段做第一轮筛选，把云覆盖、月光、火光标记的像元剔掉，再对剩余像元做阈值截断。我一般会把阈值设在0.5到1.0 nW/cm²/sr之间，具体看研究区背景噪声水平——城市密集区可以取高一点，荒漠和山区取低一点。这一步做完，数据量会减少大约15%到30%，但剩下的像元才真正代表稳定的人造灯光。

2.2 跨传感器校正：为什么2012和2013不能直接比

NPP-VIIRS在2012年刚上线时，传感器增益设置和后续年份不同，导致2012年的DN值整体偏高。更麻烦的是，2014年之后部分月份的数据经过了重新处理，辐射定标系数有调整。如果直接把2012年和2018年的月度数据放在一起算增长率，会得到一个虚高的结果。业内常用的校正方法是参考已发表的跨传感器校正系数，把2012-2013年的数据调整到与2014年之后一致的辐射基准上。另一种做法是使用年度合成产品中已经做过交叉定标的版本，但要注意不同版本之间的差异。

实际操作中，我会先检查数据提供方是否附带了定标系数文件。如果没有，就用同区域稳定灯光像元（比如城市核心区）做相对辐射归一化：选一个在九年里没有明显变化的参考区域，计算每年该区域的平均DN值，然后以某一年为基准做线性缩放。这个方法虽然粗糙，但在没有官方系数时能救急。注意，归一化之后要重新检查阈值，因为缩放会改变背景噪声的绝对水平。

2.3 年际合成与重采样：把月度栅格变成可比面板

月度数据做完清洗后，下一步是合成年度数据。常见做法有两种：年均值合成和年最大值合成。年均值适合分析灯光强度的长期趋势，年最大值适合提取城市建成区边界。我一般会同时生成两套，因为后续分析可能既要强度又要范围。合成时要注意，如果某个月份整月缺失（比如传感器故障或数据未发布），不能简单跳过，否则年均值会被剩余月份拉偏。处理方法是先统计每个像元在一年内有效观测的月份数，少于6个月的像元标记为缺失，不参与后续分析。

重采样是对齐不同年份栅格的关键。NPP-VIIRS原始分辨率约500米，但不同年份的产品可能经过不同处理，栅格对齐有微小偏移。我会用最近邻重采样统一到同一网格，因为夜光数据是连续值，双线性插值会引入虚假的平滑过渡。重采样后，用研究区的矢量边界做掩膜提取，得到每个空间单元的年度灯光总值和平均值。这一步的输出就是一份面板数据：行是空间单元，列是年份，值是灯光指标。

import rasterio import numpy as np import pandas as pd from rasterio.mask import mask import geopandas as gpd # 读取年度合成栅格和矢量边界 with rasterio.open('annual_composite_2018.tif') as src: # 用研究区边界做掩膜提取 gdf = gpd.read_file('study_area.shp') geoms = [geom for geom in gdf.geometry] out_image, out_transform = mask(src, geoms, crop=True) # 将背景值（通常为0或负值）设为NaN out_image = out_image.astype(float) out_image[out_image <= 0] = np.nan # 计算区域内灯光总值和有效像元均值 total_light = np.nansum(out_image) mean_light = np.nanmean(out_image) valid_pixels = np.sum(~np.isnan(out_image)) print(f"2018年灯光总值: {total_light:.2f}, 均值: {mean_light:.4f}, 有效像元: {valid_pixels}")

这段代码做了三件事：用矢量边界裁剪栅格、把非灯光背景设为NaN、统计灯光总值和均值。参数上要注意，mask函数的crop=True会缩小输出范围，如果后续要做多年度对齐，建议保持crop=False以维持原始网格。out_image <= 0这个阈值可以根据实际数据调整，有些产品背景值是负数，有些是0。统计时用nansum和nanmean避免背景值干扰。

2.4 灯光指标怎么选：总值、均值还是中位数

面板数据建好后，选哪个指标进模型取决于研究问题。灯光总值适合衡量区域经济总量，因为它同时反映了灯光面积和强度。灯光均值适合分析灯光强度变化，但容易被少数极亮像元拉高。中位数对极端值不敏感，适合分析灯光分布的集中程度。我一般会同时计算总值、均值和有效像元数三个指标，建模时先做相关性分析，看哪个指标与目标变量关系最稳定。

还有一个容易被忽略的指标是灯光面积，即超过某个阈值的像元数量。这个指标对城市扩张监测特别有用，因为它直接对应建成区范围。阈值的选择会影响面积计算结果，常见做法是用背景噪声的3倍标准差作为阈值，或者用固定值如5 nW/cm²/sr。我倾向于用动态阈值，因为不同区域的背景噪声水平差异很大。

3. 用Python把2012-2020年夜光数据跑成面板：从栅格到CSV的完整链路

3.1 批量读取九年栅格并对齐网格

拿到九年的年度合成栅格后，第一步是检查它们的网格是否一致。用rasterio打开每个文件，比较transform、width、height和crs。如果不一致，需要先统一到同一网格。我一般以2015年为基准年，因为这一年数据质量相对稳定，然后把其他年份重采样到2015年的网格上。重采样用最近邻，避免改变灯光值的分布。

import rasterio from rasterio.warp import reproject, Resampling import numpy as np # 以2015年为基准网格 with rasterio.open('annual_composite_2015.tif') as ref: ref_transform = ref.transform ref_width = ref.width ref_height = ref.height ref_crs = ref.crs def align_raster(input_path, output_path): with rasterio.open(input_path) as src: # 创建输出文件，网格与基准年一致 kwargs = src.meta.copy() kwargs.update({ 'transform': ref_transform, 'width': ref_width, 'height': ref_height, 'crs': ref_crs }) with rasterio.open(output_path, '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=ref_transform, dst_crs=ref_crs, resampling=Resampling.nearest ) # 对2012-2020年逐个对齐 for year in range(2012, 2021): align_raster(f'annual_composite_{year}.tif', f'aligned_{year}.tif')

这段代码的核心是reproject函数，它把源栅格重投影到目标网格。Resampling.nearest保证不引入新的灯光值，但会导致像元位置有半个像元的偏移，对大多数区域分析影响可忽略。如果研究区跨越多个UTM带，建议先统一到地理坐标系再做分析，避免投影变形。

3.2 逐像元提取灯光值并处理缺失月份

对齐网格后，下一步是逐像元提取灯光值。如果之前已经做了月度清洗和年度合成，这一步直接读年度栅格即可。但如果要控制质量，建议从月度数据重新合成，因为年度合成产品可能已经做了一些不可逆的处理。我一般会保留月度数据，自己合成年度指标，这样能控制缺失月份的处理方式。

import numpy as np import rasterio def extract_annual_metrics(monthly_files, output_csv): # monthly_files是12个月的栅格路径列表 # 读取所有月份数据 arrays = [] for f in monthly_files: with rasterio.open(f) as src: arr = src.read(1).astype(float) # 将负值和背景值设为NaN arr[arr < 0.5] = np.nan arrays.append(arr) stack = np.stack(arrays, axis=0) # 形状: (12, height, width) # 统计每个像元的有效月份数 valid_count = np.sum(~np.isnan(stack), axis=0) # 少于6个月有效的像元标记为缺失 stack[:, valid_count < 6] = np.nan # 计算年均值和年最大值 annual_mean = np.nanmean(stack, axis=0) annual_max = np.nanmax(stack, axis=0) # 保存为CSV，每行是一个像元 height, width = annual_mean.shape rows = [] for i in range(height): for j in range(width): if not np.isnan(annual_mean[i, j]): rows.append({ 'row': i, 'col': j, 'mean_light': annual_mean[i, j], 'max_light': annual_max[i, j], 'valid_months': valid_count[i, j] }) import pandas as pd pd.DataFrame(rows).to_csv(output_csv, index=False)

这段代码的关键参数是arr < 0.5这个阈值和valid_count < 6这个缺失判定。阈值0.5是经验值，如果研究区背景噪声高，可以提高到1.0。有效月份数少于6个月就标记缺失，是为了避免用少数月份代表全年。如果研究区云覆盖严重，可以放宽到4个月，但要在论文里说明。

3.3 把像元值聚合到行政区或网格单元

逐像元数据量很大，直接进模型不现实。下一步是按行政区或规则网格聚合。如果研究区有矢量边界，用geopandas和rasterio的mask函数做分区统计。如果没有边界，可以生成规则网格，比如1公里×1公里的渔网，然后统计每个网格内的灯光总值和均值。

import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np def zonal_stats(raster_path, vector_path, output_csv): gdf = gpd.read_file(vector_path) results = [] with rasterio.open(raster_path) as src: for idx, row in gdf.iterrows(): geom = [row.geometry] try: out_image, _ = mask(src, geom, crop=True) out_image = out_image.astype(float) out_image[out_image < 0.5] = np.nan total = np.nansum(out_image) mean = np.nanmean(out_image) valid = np.sum(~np.isnan(out_image)) results.append({ 'region_id': row.get('id', idx), 'total_light': total, 'mean_light': mean, 'valid_pixels': valid }) except Exception as e: results.append({ 'region_id': row.get('id', idx), 'total_light': np.nan, 'mean_light': np.nan, 'valid_pixels': 0 }) import pandas as pd pd.DataFrame(results).to_csv(output_csv, index=False)

分区统计时要注意，如果行政区边界很小，裁剪后的栅格可能只有几个像元，统计值波动很大。我一般会设置一个最小像元数阈值，比如少于10个有效像元的区域不参与后续分析。另外，mask函数在边界复杂时会比较慢，如果区域数量多，可以考虑用rasterstats库，它底层做了优化。

3.4 输出面板CSV并做基础校验

所有年份的行政区统计做完后，把九年结果合并成一张面板表。行是区域，列是年份，值是灯光指标。合并后要做三项校验：第一，检查每个区域九年里是否有缺失年份，如果某年整年缺失，需要标记并在建模时处理；第二，检查灯光总值的年际变化是否有异常跳变，比如某年突然下降50%以上，可能是数据问题而非真实变化；第三，检查区域之间的灯光值分布，看是否有极端离群值，这些可能是火光未剔除干净。

import pandas as pd # 读取九年统计结果 years = range(2012, 2021) dfs = [] for year in years: df = pd.read_csv(f'zonal_stats_{year}.csv') df['year'] = year dfs.append(df) panel = pd.concat(dfs, ignore_index=True) # 透视成面板：行是区域，列是年份 panel_total = panel.pivot(index='region_id', columns='year', values='total_light') panel_mean = panel.pivot(index='region_id', columns='year', values='mean_light') # 校验：检查缺失年份 missing = panel_total.isnull().sum(axis=1) print("有缺失年份的区域数量:", (missing > 0).sum()) # 校验：检查年际跳变 for region in panel_total.index: series = panel_total.loc[region].dropna() if len(series) > 1: pct_change = series.pct_change().abs() if (pct_change > 0.5).any(): print(f"区域 {region} 存在超过50%的年际跳变")

校验通过后，把面板表保存为CSV，后续建模直接读取。如果要做空间分析，还可以把面板表与矢量边界合并，导出为GeoPackage或Shapefile。

4. 避坑指南：夜光数据预处理里最容易翻车的五个地方

4.1 火光未剔除导致城市灯光总值虚高

现象：某区域灯光总值在2015年突然翻倍，但该区域并没有大规模建设活动。原因：天然气火炬或工业燃烧被误判为城市灯光。解决：用月度产品的QF波段筛选火光标记像元，或者用热红外波段辅助判断。如果QF不可用，可以对灯光值做空间连续性检查，孤立的高亮像元很可能是火光。

4.2 月光污染让无灯区域出现虚假亮度

现象：荒漠或山区在满月月份出现大面积低亮度值，导致年均值偏高。原因：月光被云层反射，传感器接收到了非人造光源。解决：用月度产品的月光标记筛选，或者只使用新月前后的数据。如果数据已经合成，可以用时间序列滤波剔除异常月份。

4.3 跨年比较时忽略传感器定标差异

现象：2012年灯光总值明显高于2013年，但实际经济活动没有这么大波动。原因：2012年传感器增益设置不同，DN值整体偏高。解决：使用官方定标系数，或用稳定灯光区域做相对辐射归一化。归一化后要重新检查阈值，因为缩放改变了背景噪声水平。

4.4 重采样方法选错导致灯光值失真

现象：重采样后城市核心区灯光值下降，边缘区域灯光值上升。原因：用了双线性或三次卷积插值，把高值平滑到了低值区域。解决：夜光数据是连续值但空间突变明显，重采样必须用最近邻。如果必须改变分辨率，先用最近邻重采样，再做聚合统计。

4.5 缺失月份处理不当引入偏差

现象：某区域年均灯光值在数据缺失年份反而升高。原因：缺失月份被跳过，剩余月份恰好是灯光旺季。解决：统计每个像元的有效月份数，少于6个月的标记为缺失。如果缺失月份集中在某几个月，可以考虑用同区域其他年份的同月数据做插补，但要在论文里说明。

5. 进阶技巧：用夜光面板做城市扩张监测的阈值选择与验证

5.1 动态阈值法提取建成区边界

用夜光数据提取建成区边界，核心是选阈值。固定阈值如5 nW/cm²/sr在发达地区好用，但在发展中地区会漏掉大量灯光。我一般用动态阈值：先计算研究区背景噪声的均值和标准差，然后取均值加3倍标准差作为阈值。这个方法能适应不同区域的背景差异，但需要逐区域计算。

import numpy as np import rasterio def extract_urban_extent(raster_path, output_path): with rasterio.open(raster_path) as src: arr = src.read(1).astype(float) # 计算背景噪声：取灯光值最低的20%像元 valid = arr[arr > 0] if len(valid) == 0: return background = np.percentile(valid, 20) std = np.std(valid[valid <= background]) threshold = background + 3 * std # 提取建成区 urban = (arr > threshold).astype(np.uint8) # 保存结果 kwargs = src.meta.copy() kwargs.update({'dtype': 'uint8', 'count': 1}) with rasterio.open(output_path, 'w', **kwargs) as dst: dst.write(urban, 1) print(f"阈值: {threshold:.4f}, 建成区像元数: {np.sum(urban)}")

这段代码用最低20%像元估计背景噪声，然后取3倍标准差作为阈值。参数上，20%和3倍标准差是经验值，如果研究区灯光分布集中，可以调整。提取出的建成区可以进一步做形态学滤波，去掉孤立像元。

5.2 用多年度建成区做扩张速率验证

提取出九年建成区后，可以计算每年的扩张面积和扩张速率。验证方法是与已知的城市扩张数据对比，比如土地利用遥感产品中的不透水面数据。如果两者趋势一致，说明阈值选择合理。如果差异大，需要调整阈值或检查数据质量。

表格里的差异控制在5%以内，说明动态阈值法在这个研究区可用。如果差异超过10%，需要检查是否有火光未剔除或月光污染。

5.3 灯光面板与统计数据的相关性检验

最后一步是用灯光面板与GDP、人口等统计数据做相关性检验。如果相关性高，说明预处理后的数据能作为经济活动的代理变量。我一般会算皮尔逊相关系数和斯皮尔曼秩相关系数，前者对线性关系敏感，后者对单调关系敏感。如果两个系数都高于0.7，说明数据质量可靠。

做这个检验时，要注意统计数据的口径和年份匹配。有些区域的GDP数据是市辖区口径，而灯光数据是全市口径，直接比较会低估相关性。我踩过这个坑，后来统一了口径才得到合理结果。希望帮到你。

本文还有配套的精品资源，点击获取