企业官网建设流程全解析

简介：面向从事城市扩张、经济活动空间分布与区域发展研究的学者及遥感数据处理初学者，这份2012—2020年NPP-VIIRS夜间灯光数据集提供了可直接开展分析的基础底图。数据以500m分辨率覆盖全国，基于WGS84中国矢量边界裁剪，并统一转换为Asia_Lambert_Conformal_Conic投影，便于面积量算与跨年度比较；针对原始产品中的负值和油气燃烧（gas flares）引起的极端异常，参照Kang Wu等（2019）以472.86为像元辐射阈值进行剔除，使各年度栅格在时间序列上更具可比性。压缩包共36个文件，以9个tif栅格为主体，配套tfw坐标信息、xml元数据与ovr金字塔文件，整体约176.73MB，解压后可按年份快速定位。已有3511人学习下载，适合用于夜间灯光总量提取、时空变化分析、城市化监测及论文复现等场景，省去裁剪、投影与异常值处理的重复劳动。

1. 2012-2020年NPP-VIIRS夜光遥感数据不是“越亮越真”：先看它被预处理成了什么

2012-2020年NPP-VIIRS夜光遥感数据的预处理版本，最容易被当成“开箱即用”的面板变量：有人直接对栅格求均值，往经济指标上一回归，发现2013到2014年一条硬跳变，才回头怀疑数据。预处理后的2012-2020年NPP-VIIRS夜光遥感数据，核心价值不是“亮”，而是把辐射值、背景噪声、异常像元和年际口径尽量拉到可比；它适合做区域经济活跃度、城市扩张、基础设施覆盖的长期观察，不适合拿单年极值讲结论。先弄清楚它被处理到什么程度，再谈模型。

2. 预处理后的NPP-VIIRS夜光数据改了什么：辐射定标、杂散光与年际可比性

2.1 NPP-VIIRS DNB的原始量纲为什么不能直接跨年比较

夜光遥感不是普通可见光影像。NPP-VIIRS DNB 在低照度条件下工作，原始数字值要经过辐射定标才变成辐亮度，常见单位是 nW·cm⁻²·sr⁻¹。同一像元在两年里数值不同，可能来自地表灯光明暗，也可能来自定标系数、杂散光校正、月光与云掩膜、合成窗口。预处理后的2012-2020年NPP-VIIRS夜光遥感数据如果只做了负值截断和异常值剔除，没有做跨年辐射基准拉齐，拿它做时间序列就像拿两把刻度不同的尺子量同一张桌子。

杂散光是一个绕不开的词。早期低照度通道容易受到非目标光源污染，海岸线、极区、云边缘经常出现成片伪亮。常见预处理会做杂散光校正、背景噪声扣除、火点剔除和气辉去除，但不同处理链的力度不一样。2012-2020年跨了九个年份，传感器在轨状态、合成规则、掩膜策略都可能调整。判断一份数据能不能跨年用，不能只看它叫“预处理后”，要看它有没有把辐射基准和合成口径写清楚。

我一般先问三个问题：单位是辐亮度还是无量纲数字值；年度值由月度合成还是年均直接合成；异常值是被置为缺失、被截断，还是被替换成邻域插值。前两个决定可比性，后一个决定你会不会把“处理痕迹”当成“经济波动”。如果数据说明里只写了“已去噪”，没有写年际一致性处理，就要自己补验证。

2.2 拿到一份“预处理后”数据，先查四个字段

不要急着写循环读栅格。先做一张核对表，把数据说明、文件头、栅格属性放在一起看。下面四项是最低门槛，缺一项都可能在后面翻车。

这张表看起来笨，但能省掉后面很多返工。尤其注意 NoData 和 0 的区别：夜光数据里 0 可能代表真实无灯，也可能代表填充。把填充值当 0 参与均值，区域越靠边界，均值被拉得越低。预处理后的数据如果已经把 NoData 写成 NaN，读栅格时要用masked=True或显式设置nodata，不要让它悄悄变成 0。

还要看年度合成规则。月度合成里，云掩膜会改变有效观测数；有的产品用平均值，有的用中位数，有的只在无云像元足够多时才输出。2012-2020年跨度长，某一年有效观测少，年度值就可能整体偏暗。这个偏差不是地表变暗，而是观测窗口变窄。看到某年整体偏低，先查有效观测数和合成方法，不要先加“年份虚拟变量”把它盖过去。

2.3 年际可比性的最小验证：用稳定灯光样本做桥接

在没有完整定标文档的情况下，我通常做一个最小验证：在研究区A里选一块连续建成区，要求2012-2020年没有大拆大建、没有大型临时工程、没有频繁火点，取每年同一批像元的中位数，做年度比值。这个样本不需要很大，但必须稳定。稳定样本的年度中位数如果出现超过10%-15%的跳变，就要回头查处理规则，而不是直接进入回归。

具体步骤可以这样落地。第一步，用2012年和2020年两期数据做交集，先去掉任一年为负值、零值或 NoData 的像元。第二步，在交集里选建成区连片区域，按研究区A的边界裁一块，避免把农田、港口、机场跑道混进来。第三步，计算每年稳定样本的中位数M_t，再算比值R_t = M_t / M_2012。第四步，算同一像元年际变异系数CV = std / mean，把 CV 过大的像元剔掉，剩下的再重复一次。第五步，把R_t序列和年度有效像元占比放在一张表里看，跳变年份如果同时伴随有效像元下降，优先怀疑掩膜和合成，而不是灯光本身。

注意：稳定样本不是越亮越好。过亮的核心商业区容易饱和，过暗的郊区又容易被背景噪声吞掉，取中间亮度、连片、变化小的建成区更可靠。

这个验证不能替代正式定标，但能快速判断“这版数据能不能直接跨年”。如果稳定样本的年度比值平滑，区域统计的年度趋势才有解释力；如果稳定样本自己都跳，后面所有回归系数都要打问号。我的习惯是先把这条曲线画出来，再决定要不要做年际校正。校正方法可以很简单：以2012年为基准，用稳定样本的年比值做线性缩放；也可以更保守，只在模型里加入年份固定效应并报告敏感性。关键是别把处理痕迹当成真实变化。

3. 用xarray+rioxarray把九年夜光栅格读成时间立方体：坐标系、掩膜和内存三个关口

3.1 目录结构与年份命名：先把文件清单变成可核对的表

第一步不是读数据，是把文件清单管住。2012-2020年九个年份，文件名里可能带月份、合成方式、版本后缀，直接按字符串排序很容易把2012和2012_01混在一起。我一般用正则提年份，生成年份 -> 文件路径的字典，后面所有读取都从这张表走。

from pathlib import Path import re root = Path("data/npp_viirs_2012_2020") items = [] for f in root.glob("*.tif"): # 只认文件名里的四位年份，避免月份和后缀干扰 m = re.search(r"(20\d{2})", f.stem) if m: items.append((int(m.group(1)), f)) year_files = dict(sorted(items)) print(year_files.keys()) print(len(year_files))

这段代码的逻辑很直接：遍历目录下的 TIFF，用正则抓四位年份，按年份排序后转成字典。参数root指向你放九年栅格的目录；glob("*.tif")只匹配 TIFF，如果数据是 NetCDF，把后缀改掉即可。打印年份键是为了核对有没有缺年、重年。若同一年出现两个文件，先不要继续，回到数据说明确认哪一个是年度合成、哪一个是月度合成。这个检查看起来琐碎，但能挡住“把12个月当12个样本”的低级错误。

3.2 读取单年栅格并统一到同一网格：rio.open_rasterio 与 reproject_match

九年数据最常见的问题是网格不一致：有的年份像元大小一样但原点偏了半个像元，有的年份换了投影，有的年份多了一个波段维度。直接拼接会得到错位的结果，区域统计就会把边界像元算到隔壁单元。我的做法是先选2012年作为参考网格，其他年份全部reproject_match到它上面。

import rioxarray as rxr ref = rxr.open_rasterio(year_files[2012], masked=True).squeeze(drop=True) print(ref.rio.crs) print(ref.rio.resolution()) print(ref.rio.nodata) print(ref.shape) def load_aligned(path): da = rxr.open_rasterio(path, masked=True).squeeze(drop=True) if da.rio.crs is None: raise ValueError("栅格缺少CRS，先补投影信息") if da.rio.crs != ref.rio.crs: da = da.rio.reproject(ref.rio.crs) # 统一到参考网格：范围、分辨率、原点一次对齐 da = da.rio.reproject_match(ref) return da

masked=True会让 NoData 变成掩膜，后面做统计时不会把填充值当有效值。squeeze(drop=True)去掉单波段维度，避免出现band=1的冗余维。reproject只处理投影差异，reproject_match进一步把范围、分辨率和原点对齐到参考栅格。这里把2012年当参考，是因为它是序列起点，后面所有年份都跟它对齐，口径统一。若你的研究区在2012年边界不完整，也可以选覆盖最全的一年做参考，但要在文档里写清楚。

reproject_match默认会做重采样，连续型辐射值建议用双线性或平均，不要用最近邻。最近邻适合分类栅格，用在夜光辐亮度上会产生块状伪影。这个参数在rioxarray里可以通过resampling控制，常见做法是Resampling.bilinear。如果你不确定，先拿一年做对比：最近邻和双线性的区域均值差多少，差异大就说明网格错位明显，需要回头查原始网格。

3.3 九年数据堆栈：用 xarray.concat 而不是 for 循环反复读盘

单年读取完成后，把九年堆成一个year, y, x的时间立方体。这样后面做年度统计、区域掩膜、时间序列都能一次完成，不用每年重新读盘。数据量大时，直接 concat 会吃内存，所以我一般加chunk，让 dask 按年份分块。

import numpy as np import xarray as xr years = sorted(year_files) das = [load_aligned(year_files[y]) for y in years] stack = xr.concat(das, dim="year").assign_coords(year=years) # 负值在夜光产品里多为噪声或处理残留，先置为缺失 stack = stack.where(stack >= 0) # 已安装 dask 时按年份分块，避免一次性把九年数据全塞进内存 stack = stack.chunk({"year": 1, "x": 2048, "y": 2048}) print(stack)

xr.concat按year维拼接，assign_coords把年份坐标补上，后面stack.sel(year=2020)才不会拿到错位数据。stack >= 0是最保守的负值处理：负值置为 NaN，不参与均值。若数据说明明确保留负值作为噪声估计，就不要一刀切，改为在统计时用min_count控制有效像元数。chunk里的2048不是固定值，按你的内存和像元总数调；原则是单块能放进内存，且年份维不要合并，方便按年计算。

注意：stack.where(stack >= 0)只处理负值，不处理极端亮斑。火点、气辉、异常反光需要额外的 QA 或分位数掩膜，不要指望负值过滤能解决所有异常。

到这一步，你手里应该有一个带 CRS、带年份坐标、负值已掩膜的时间立方体。接下来才轮到区域统计。如果这一步的网格没对齐，后面所有区域均值都会带系统偏差，而且这种偏差很难在回归里被诊断出来。

4. 区域夜光指标怎么算：掩膜、重采样、分区统计的四个必调参数

4.1 先定分析单元：行政边界、渔网网格还是自定义缓冲区

夜光数据本身是栅格，但大多数分析要落到区域上。分析单元可以是行政边界、规则渔网、缓冲区，也可以是自定义的交通走廊。选哪种取决于问题：看区域经济活力，行政边界更贴近统计口径；看城市扩张，渔网能保留空间梯度；看基础设施覆盖，缓冲区更合适。不要一上来就用最细的行政单元，夜光像元有空间溢出，过细的单元会把邻区灯光算进来。

我一般先做两套：一套是分析单元，用来出主结果；一套是粗网格，用来做空间稳健性。比如主结果用区县边界，稳健性用5 km渔网。两套结果的趋势方向一致，才敢说结论不是边界画出来的。边界数据要统一到栅格 CRS，不要用经纬度直接算面积。下面代码以zones.gpkg为例，字段zone_id是区域唯一标识。

4.2 重采样不是越小越好：500 m、1 km 与聚合网格的取舍

NPP-VIIRS DNB 的原始像元接近750 m，常见预处理产品会重采样到500 m或1 km。很多人以为网格越细越好，其实夜光的光溢出和配准误差不会因为重采样变细而消失。500 m网格适合看城市内部形态，但文件大、噪声像元多；1 km网格更适合区域面板，统计更稳；再聚合到5 km，适合跨区域比较和减少边界效应。

重采样方法也要固定。连续辐射值用平均或双线性，不要用最近邻；聚合到粗网格时用均值，不要用总和后再除，否则像元数不同会引入尺度偏差。最关键的一条：所有年份必须用同一套重采样规则，不能2012年用最近邻、2020年用双线性，那样年度跳变可能来自算法而不是地表。

4.3 分区统计的四个参数：all_touched、nodata、统计量和最小有效比例

分区统计看起来只是“把栅格值按多边形求平均”，实际有四个参数直接决定结果。第一个是all_touched：设为False时只统计中心落在多边形内的像元，结果保守；设为True时边界像元也算，区域会被“放大”。第二个是nodata：必须显式设为 NaN 或数据说明里的填充值，不能让它变成0。第三个是统计量：均值反映强度，总和反映总量，二者不要混用。第四个是最小有效像元比例：有效像元太少时，均值不可信，应该置为缺失。

import numpy as np import geopandas as gpd import pandas as pd from rasterstats import zonal_stats zones = gpd.read_file("data/zones.gpkg").to_crs(stack.rio.crs) rows = [] for y in stack.year.values: da = stack.sel(year=y) stats = zonal_stats( zones, da.values, affine=da.rio.transform(), nodata=np.nan, stats=["mean", "sum", "count"], all_touched=False, ) for zid, s in zip(zones["zone_id"], stats): rows.append({ "zone_id": zid, "year": int(y), "light_mean": s["mean"], "light_sum": s["sum"], "valid_pixels": s["count"], }) panel = pd.DataFrame(rows) print(panel.head())

这段代码按年份循环，把每个区域的均值、总和和有效像元数收进长表。affine=da.rio.transform()把数组坐标转成地理坐标，nodata=np.nan防止填充值参与统计，all_touched=False让边界更保守。stats里同时要mean和sum，是因为很多分析会用均值代表强度、用总和代表规模，两者不能互相替代。valid_pixels不是装饰字段，后面要用它过滤掉有效像元过少的区域。

注意：count的含义在不同库版本里可能略有差异，拿到结果后先用一个已知区域手工核对一遍。不要直接拿count当总面积，像元面积要另算。

4.4 把年度统计变成面板表：区域ID、年份、均值、总和、有效像元占比

分区统计输出的是长表，进入模型前还要补两列：像元面积和有效像元占比。有效像元占比可以用valid_pixels除以该区域落在栅格范围内的总像元数。总像元数可以用同样的分区统计对一张全1栅格求count得到，也可以在GIS里预先算好。这个比例低于0.6的区域，我一般先标记，不直接删除，做敏感性分析时再决定是否剔除。

面板表至少保留zone_id、year、light_mean、light_sum、valid_pixels、valid_ratio。如果后面要做增长率，最好再加一列light_mean_lag，避免在代码里反复shift出错。区域ID不要用中文名，统一用编码，防止合并时编码不一致。年份用整数，不要用字符串，否则排序会出现"2020"排在"2019"前面的问题。

到这里，区域夜光指标才算落地。真正容易出错的不是公式，而是边界、NoData、重采样和有效像元比例这四个环节。任何一个环节口径不一致，年度趋势都可能被伪造出来。

5. 2012-2020年夜光数据排查清单：五类翻车现场与补救顺序

5.1 负值成片：先分清填充值、背景扣除和真实暗区

现象：区域均值出现大片负值，或者海岸线、沙漠边缘整片为负。原因通常有三种：NoData 被写成负值；背景扣除把暗区压到0以下；杂散光校正的残差没有清理。解决时先看栅格头里的nodata，再用masked=True读取；如果负值仍然成片，按数据说明决定是置为 NaN 还是截断为0。不要直接把负值当0参与均值，那会把暗区抬亮，区域趋势看起来更平滑，但这是人为的。

5.2 某一年整体偏暗或偏亮：先查有效观测数，不要先加年份虚拟变量

现象：2012-2020年序列里某一年整体下沉，所有区域同向变化，回归里年份虚拟变量显著。原因可能是年度合成时有效观测少、云掩膜过严、月度权重不同，也可能是辐射基准调整。解决顺序是：先统计每年的有效像元占比；再看稳定灯光样本的年度中位数；最后查数据说明里的合成规则。如果有效像元占比明显低，优先考虑重算年度值或做缺失标记，而不是用模型年份效应掩盖。

5.3 城市边缘出现异常亮斑或条带：火点、气辉和配准误差要分开处理

现象：城市边缘、港口、机场附近出现孤立极亮像元，或者沿条带方向出现规则亮线。原因可能是火点、气辉、船只灯光、云边缘残留，也可能是年份之间配准偏移导致亮区错位。解决时先用 QA 层或分位数阈值把极端值置为缺失，再检查两年影像的亮区边缘是否重合。配准误差不能靠阈值解决，必须回到reproject_match统一网格。极端值处理要做敏感性分析：去掉前后区域均值变化超过20%的区域，结论要谨慎。

5.4 区域统计对边界敏感：换个边界结论就变，先查缓冲区和 all_touched

现象：同一个城市，换一套边界后夜光均值变化很大，增长率方向甚至相反。原因通常是边界穿过了亮区边缘，all_touched=True把邻区灯光算进来，或者不同年份用了不同边界。解决时固定一套边界，所有年份共用；主结果用all_touched=False，稳健性用True对比；对边界附近区域做缓冲区敏感性分析。如果两种口径下趋势方向不一致，说明结论对边界太敏感，应该缩小研究范围或改用网格分析。

5.5 有效像元比例低但均值正常：最隐蔽的坑

现象：某区域夜光均值看起来合理，但有效像元只占区域面积的30%，而且每年比例波动很大。原因是云掩膜、负值掩膜或 NoData 把大部分像元去掉了，剩下的像元可能集中在亮区，均值被高估。解决时把valid_ratio作为硬门槛，低于阈值的区域标记为缺失；同时报告不同阈值下的结果。我的习惯是主结果用0.6，稳健性用0.5和0.7各跑一遍，如果结论只在0.5成立，就不敢把它当主结论。

6. 把2012-2020年夜光序列用进面板分析：先做同质区检验，再决定是否取对数

6.1 先做同质区检验，再决定是否取对数

夜光均值右偏，很多人第一反应是取对数。取对数没问题，但前提是同一区域内年度变化不是由掩膜和边界造成的。我一般先做同质区检验：把研究区A按稳定灯光样本分成高、中、低三组，分别看2012-2020年的中位数曲线。三组曲线平行，说明年际口径基本一致；如果只有高亮组跳变，低亮组平稳，问题多半在饱和或异常值处理，不在经济变化。检验通过后再取log1p，可以压缩极值，但不能修复口径问题。

import numpy as np panel["ln_light"] = np.log1p(panel["light_mean"].clip(lower=0)) panel["weight"] = panel["valid_ratio"].clip(0, 1)

clip(lower=0)防止负值进入对数，log1p在均值为0时也能返回0。weight用有效像元占比做权重，适合加权回归或加权平均。参数不在多，关键是每一步都能解释：为什么取对数、为什么加权、剔除阈值是多少。把这些写进分析日志，比事后补“稳健性检验”有用得多。

6.2 一个稳健技巧：用有效像元占比做权重，别让边界像元主导结果

如果非要用一个技巧收尾，我会选“有效像元占比加权”。夜光数据的误差不是均匀的：城市核心区有效像元多，边界和云掩膜边缘有效像元少，后者噪声更大。直接对区域均值做回归，等于默认每个区域同样可靠。用valid_ratio做权重，可以降低低质量区域的影响，又不至于把它们全删掉。权重不要用原始像元数，像元数大的区域会主导结果；用比例，范围固定在0到1之间，更好解释。

我现在拿到任何一版2012-2020年NPP-VIIRS夜光数据，都会先做三件事：核对年份和单位，画稳定样本年度曲线，检查有效像元占比。这套笨办法救过我两次，一次是发现某年合成规则变了，一次是发现边界像元被all_touched放大。夜光数据很亮，但坑也很亮。希望帮到你。

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