企业官网建设流程全解析

简介：这是一份面向GIS与遥感分析人员的ArcGIS夜间灯光处理指南，系统讲解夜间灯光数据的提取与强度计算方法。文档以DMSP/OLS、VIIRS等栅格数据为输入，围绕行政区域矢量数据准备、属性筛选、按掩膜提取、分区统计、字段连接与分级符号化等环节展开，并突出坐标系匹配、区域选择逻辑和统计字段设置等易错细节，读者可按步骤在ArcGIS中完成从原始灯光影像到区域灯光强度图的全流程操作。整包仅1个docx文档，共6.52MB，步骤配有界面说明与参数提示，适合独立学习或课程实验参考。目前已有3651人学习下载，对从事城市规划、环境监测或区域经济研究的GIS使用者尤其实用。

1. 夜间灯光的提取及强度计算：两个动作把黑夜变成城市活动数据

夜间灯光数据的提取，是把夜间遥感影像里“有人在用电”的亮像元从暗背景切出来；切完后把亮像元的 DN 值或辐射亮度汇总成指标，就是强度计算。两个动作连起来，就得到一张比统计年鉴更新更快的活动强度图。

城市边界画到哪、经济重心往哪偏，都能从亮斑变化里读出来。这篇文章写给三类人：做长时序分析的从业者、把夜间灯光当城市研究代理变量的规划研究者，以及刚入坑想跑通全流程的初学者。

下面按数据选型、预处理、提取、强度计算、避坑的顺序展开，每步都给可复现的代码和参数说明。里面的阈值和经验值都是项目里反复调过的，直接抄基本能跑，但你自己的数据难免要微调。

2. 夜间灯光数据选型与预处理：三种星源怎么选，裁剪重投影怎么做

2.1 三种主流夜间灯光数据源的特性与选型逻辑

做夜间灯光提取，第一步不是写代码，而是选对数据源。目前从业者用得最多的是长时序的 DMSP/OLS、覆盖当前的 NPP/VIIRS，以及近几年出现的国产高分辨率夜间灯光数据。这三类产品在时间跨度、空间分辨率、数值含义上差异很大，选错了源，后面的阈值和强度统计整体作废。

DMSP/OLS 覆盖 1992-2013 年，空间分辨率约 1 公里，DN 值范围 0-63。它采用 6 位量化，动态范围小，大城市中心很容易饱和，整片像元直接顶到 63。NPP/VIIRS 从 2012 年延续至今，分辨率约 500 米，保存的是辐射亮度值，单位是 nW/cm²/sr，没有饱和问题，但背景噪声明显比 DMSP 多。国产高分辨率夜间灯光数据时间跨度短、覆盖不连续，不过胜在分辨率高，适合城市内部精细研究。

选型逻辑其实就三条。要做 1992-2013 的完整历史趋势，只能选 DMSP/OLS；要做近十年的定量统计或连续年度对比，优先 NPP/VIIRS；要细化到单个城市内部街区，再考虑国产高分辨率数据。没有特殊必要，不要把三种源混进同一条时间序列，因为数值口径完全不同。下载时 VIIRS 优先选 vcmcfg 版本，它在合成阶段去掉了火光和背景值影响，比 vcm-orm 更适合城市灯光分析。

注意：网上不少教程把 DMSP 的 DN 值和 VIIRS 的辐射亮度直接拿来拟合回归，这是最常见的数据源错误。数值含义不同，必须先统一量纲再做趋势分析。

2.2 预处理两步走：矢量裁剪、等积投影与异常值清理

拿到原始影像后，我一般先做两步预处理：用研究区矢量边界裁剪，再把影像投影到等积坐标系（如 UTM 或 Albers）。裁剪的目的是把大量无意义的背景 0 值剔除，否则后续 Otsu 阈值会把直方图淹没在 0 值里；投影的目的则是保证每个像元对应固定的地面面积，后面的面积和强度统计才有物理意义。

import rasterio import numpy as np import geopandas as gpd from rasterio.mask import mask # 1. 读取 VIIRS 年度合成影像，先看元信息 src_path = "VNL_v2_npp_2022_global_vcmcfg.tif" with rasterio.open(src_path) as src: profile = src.profile.copy() print(f"原始 CRS: {src.crs}，范围: {src.bounds}") # 2. 用研究区矢量裁剪 aoi = gpd.read_file("study_area.shp") # 矢量坐标系必须和影像一致，否则 mask 会报错或裁空 aoi = aoi.to_crs(profile["crs"]) with rasterio.open(src_path) as src: out_image, out_transform = mask( src, aoi.geometry, crop=True, nodata=0 ) out_profile = profile.copy() out_profile.update({ "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform, "nodata": 0 }) with rasterio.open("study_ntl_2022.tif", "w", **out_profile) as dst: dst.write(out_image)

这段代码里最关键的是aoi.to_crs(profile["crs"])，很多初学者跳过这步，导致矢量边界和影像坐标系不一致，裁剪出来是空的或者错位。nodata=0的意思是把区域外像元统一置为 0，这样后面所有“有效像元”的判断口径就统一了。VIIRS 原始数据是 float32 浮点类型，读取时不要用 int 承接，否则辐射亮度的小数部分全部被截掉。

如果不想写 Python，用 GDAL 命令行也能完成同样的事，适合在服务器上一键批处理：

gdalwarp -cutline study_area.shp -crop_to_cutline -dstnodata 0 \ VNL_v2_npp_2022_global_vcmcfg.tif study_ntl_2022.tif

裁剪完先别急着提取，打开影像看一眼数值分布。VIIRS 年度合成影像的辐射亮度正常是大于 0 的浮点数，如果出现负值或极端大值，多半是原始数据质量问题，先用分位数截断再进提取环节。我一般会顺手打印 min、max 和直方图分位数，确认 95 分位和最大值之间的差距没有超过一个数量级，再继续往下走。

重投影放在裁剪之后做。UTM 投影下像元尺寸是固定米数，后面的面积统计和灯光重心计算都依赖这个前提。连续辐射值用 bilinear 或 cubic 重采样，不要用 nearest，否则提取出来的边界会有明显锯齿。

from rasterio.warp import calculate_default_transform, reproject, Resampling with rasterio.open("study_ntl_2022.tif") as src: transform, width, height = calculate_default_transform( src.crs, "EPSG:32650", src.width, src.height, *src.bounds ) kwargs = src.profile.copy() kwargs.update({ "crs": "EPSG:32650", "transform": transform, "width": width, "height": height }) with rasterio.open("study_ntl_utm50.tif", "w", **kwargs) as dst: for i in range(1, src.count + 1): reproject( source=rasterio.band(src, i), destination=rasterio.band(dst, i), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs="EPSG:32650", resampling=Resampling.bilinear )

3. 灯光区域提取：Otsu 自动阈值、经验阈值与形态学后处理

3.1 阈值提取的前提假设：为什么直方图形状决定选型

灯光提取的核心假设是：夜间灯光影像的像元值分布可以分成“亮”和“暗”两类，亮的代表有电力活动，暗的代表没有。理想情况下直方图呈双峰分布，一个峰是背景噪声和黑暗地表，另一个峰是城市灯光，两个峰之间的谷底就是最佳分割点。

但真实影像远没有这么理想。DMSP 大量像元饱和在 DN=63，VIIRS 在乡村和沿海有大量低值噪声，直方图往往是一个巨大的背景峰加上一条长长的亮值拖尾，而不是漂亮的双峰。这时手选阈值就成了一个高风险动作：同一个地区，2000 年影像和 2010 年影像，手动选的阈值不可能一致，时间序列就被人为打断了。所以我定了一个规矩：单张影像做一次性可视化展示时可以用经验阈值；凡是做多年度对比或定量统计，一律用自动阈值，并把阈值选择方法写进处理记录。

常用经验阈值可以作为参考：DMSP 灯光区域一般取 DN >= 10 或 DN >= 12；VIIRS 年度合成影像取辐射亮度大于 1 或 2 nW/cm²/sr。这两个数字不是万能药，不同地区、不同产品版本差异很大，一定要结合你研究区的直方图去验证。

3.2 Otsu 自动阈值代码与两个必调参数

Otsu 的本质是穷举所有可能的分割阈值，让分割后两类像元值的类间方差最大。它不需要人工干预，对单张影像和批量处理都很友好，缺点是默认数据分布大致是双峰的，遇到严重偏态分布时，Otsu 会把阈值切在主要峰的一侧，导致灯光区域要么过大要么过小。

import numpy as np from skimage.filters import threshold_otsu from scipy import ndimage def extract_lit_area(data: np.ndarray, method: str = "otsu") -> tuple: """提取灯光像元掩膜。 参数: data: 单波段夜间灯光影像，float32 method: otsu 或 empirical 返回: (掩膜数组, 阈值) """ # 只取正值为有效像元，避免背景 0 淹没直方图 valid = data[data > 0].ravel() if valid.size < 100: raise ValueError("有效像元过少，检查裁剪范围或数据质量") if method == "otsu": thres = threshold_otsu(valid) else: # 经验法：均值加 1 倍标准差，适用于噪声较多的 VIIRS thres = valid.mean() + 1.0 * valid.std() lit_mask = (data >= thres).astype(np.uint8) return lit_mask, float(thres) with rasterio.open("study_ntl_utm50.tif") as src: data = src.read(1).astype(np.float32) lit_mask, thres = extract_lit_area(data) print(f"Otsu 自动阈值: {thres:.4f} nW/cm²/sr")

代码里两个要点需要解释。第一，data[data > 0]只取正值像元参与阈值计算，因为裁剪时把区域外统一设成了 0，这些背景像元数量巨大，会把直方图完全压扁。第二，Otsu 返回的是像元值的分界点，VIIRS 影像单位是辐射亮度，所以阈值也带 nW/cm²/sr 单位，打印出来看一眼，如果小于 0.1 或大于 10，大概率是数据分布有问题。

提示：如果 Otsu 阈值明显偏低（比如只有 0.3），通常意味着数据里混入大量低值噪声。这时不要急着调参数，先回预处理环节检查，看是不是裁剪范围包含了太多未供电区域。

如果直方图偏态特别严重，可以先把有效像元取对数再跑 Otsu，得到的阈值在对数空间里，需要再指数还原。这个 trick 在处理大范围农村区域时很有效，因为大片低值像元会把阈值往低处拉。

3.3 形态学后处理：孤立亮斑清理与空洞填充的尺度控制

阈值得到的二值掩膜通常很碎：有孤立的单像元亮斑（传感器噪声），也有城市内部的黑洞（公园、水体、农田）。我用两步形态学处理来清理，但会刻意控制尺度，避免把真实的城市结构削掉。

# 1. 清除面积小于阈值的孤立亮斑 min_pixels = 9 # 500m 分辨率下约 2.25 平方公里 label_img, num_labels = ndimage.label(lit_mask) sizes = ndimage.sum(lit_mask, label_img, range(num_labels + 1)) remove_mask = sizes < min_pixels remove = remove_mask[label_img] lit_clean = lit_mask.copy() lit_clean[remove] = 0 # 2. 填充城市内部的小空洞 structure = np.ones((7, 7), dtype=bool) lit_filled = ndimage.binary_fill_holes( lit_clean, structure=structure ).astype(np.uint8)

min_pixels = 9的经验含义是：在 500 米分辨率 VIIRS 数据上，小于 9 个像元的孤立亮斑几乎不可能是城市，多数是单次火点或传感器噪声。binary_fill_holes里我用 7x7 结构元素，它能填掉城市内部尺度小于 3 个像元的空洞，但保留大尺度的水体隔断。结构元素开得太大，会把两片不相连的城区强行连通，小城市研究里特别容易翻车。

清理完之后做一次人工目视检查：把掩膜叠加到同期高分辨率影像上，看亮斑边界是否落在可见的城市建成区内。这一步花不了几分钟，能拦住 80% 的提取逻辑错误。掩膜存储建议输出为单波段 GeoTIFF，方便后续直接和强度统计代码对接。

4. 强度计算：灯光总量、平均强度与灯光重心怎么算

4.1 三个核心指标的物理含义与选型场景

提取完成后，下一步是把亮像元的数值汇总成有意义的强度指标。最常用的三个指标是灯光总量、平均强度和灯光面积占比，它们在实际应用中是互补关系，不是互相替代。

灯光总量等于所有亮像元数值之和，对应区域的总灯光产出，GDP 回归、用电量估算基本都用它。平均强度等于总量除以亮像元数量，反映光的集中程度，适合描述城市紧凑度。灯光面积占比等于亮像元数量除以区域总像元数量，近似城市建成区的相对规模，做城市扩张研究时优先看它。

4.2 区域尺度强度统计的代码实现与单位换算

实际分析里很少对整幅影像算一个数，通常是把研究区切成行政区、格网或缓冲区，逐个区域统计。下面按区域编号统计的实现中，区域编号可以从矢量面转栅格得到，也可以用固定大小的渔网格网。

import numpy as np def zonal_light_stats(lit_mask: np.ndarray, intensity: np.ndarray, zones: np.ndarray) -> list: """按区域统计灯光强度指标。 参数: lit_mask: 灯光提取结果的二值掩膜 intensity: 原始辐射亮度/DN 影像 zones: 区域编号栅格，0 表示无效区 返回: 每个区域的指标字典列表 """ results = [] for zid in np.unique(zones): if zid == 0: continue zone_mask = zones == zid lit_in_zone = lit_mask & zone_mask n_lit = int(np.count_nonzero(lit_in_zone)) n_zone = int(np.count_nonzero(zone_mask)) if n_lit == 0: results.append({ "zone_id": int(zid), "tnli": 0.0, "mean_intensity": 0.0, "lit_ratio": 0.0, "n_lit": 0 }) continue sum_light = float(np.sum(intensity[lit_in_zone])) results.append({ "zone_id": int(zid), "tnli": round(sum_light, 4), "mean_intensity": round(sum_light / n_lit, 4), "lit_ratio": round(n_lit / n_zone, 4), "n_lit": n_lit }) return results

lit_mask & zone_mask这一步在布尔数组上做按位与，性能很好，几千万像元的影像也能秒级完成。tnli的单位和原始数据一致，VIIRS 影像的 TNLI 是“辐射亮度 × 像元数”，不能直接对应物理上的总功率，但在回归模型里作为相对指标完全够用。如果要做面积统计，在 UTM 投影下把n_lit乘以单像元面积即可，比如重采样到 500 米后每个像元约 0.25 平方公里。

4.3 灯光重心计算：城市扩张方向的隐性指标

灯光重心是把每个亮像元的强度当作权重，计算区域灯光的加权中心。一次计算结果看不出什么，把 2000、2005、2010、2015 四期影像的灯光重心连成一条线，城市扩张的方向就一目了然。

def light_centroid(lit_mask: np.ndarray, intensity: np.ndarray, transform) -> tuple: """计算灯光重心，返回投影坐标系下的坐标。""" rows, cols = np.nonzero(lit_mask) if rows.size == 0: return None vals = intensity[rows, cols].astype(np.float64) total = vals.sum() if total <= 0: return None # 像元中心坐标：x = transform.c + col * transform.a xs = transform.c + cols * transform.a ys = transform.f + rows * transform.e cx = np.sum(vals * xs) / total cy = np.sum(vals * ys) / total return (cx, cy)

这个函数要求影像已经是投影坐标系，因为transform.a和transform.e在投影坐标下是米，算出来的重心坐标才有距离意义。如果直接拿经纬度影像算重心，得到的坐标差值是度数而不是实际距离，跨纬度地区完全不可比。

注意：灯光重心对极端亮值非常敏感，DMSP 中心饱和像元会把重心拉向中心。用 VIIRS 算重心比 DMSP 可靠得多，因为它没有饱和问题。

5. 夜间灯光处理避坑指南：饱和、噪声与投影翻车实录

5.1 DMSP 城市中心饱和：DN 封顶 63，强度被系统性低估

现象：大城市核心区 DN 值全部是 63，像被切了顶，灯光总量明显低于同等级城市。

原因：DMSP/OLS 探测器是 6 位量化，动态范围只有 0-63，城市中心超高亮度直接触顶。这不是数据下载错了，是传感器本身的物理限制。解决：历史趋势研究里，把 DMSP 强度指标定位成相对量而不是绝对量，只用来描述年份间的相对变化；近十年的定量分析一律改用 VIIRS。如果一定要用 DMSP 做强度回归，考虑在饱和区域用灯光面积替代强度，或把城市分成中心区和外围分别统计，中心区不做强度对比。

5.2 VIIRS 背景噪声：渔船、火点与临时光源混入亮斑

现象：沿海和乡村地区出现大量孤立小亮斑，面积不大但数量很多，灯光面积占比被明显抬高。

原因：VIIRS 灵敏度比 DMSP 高一个量级，海上渔船、天然气火把、夜间施工灯光都会进入影像。官方合成已经过滤了火光，但渔船灯光很难完全去除。解决：阈值提取后先做面积过滤，清除小于 3-9 个像元的孤立斑块；海边城市额外叠加陆地掩膜，把海洋像元先剔除。如果噪声仍顽固，把阈值从 Otsu 结果上调 20%-50%，代价是城市周边小村落可能被一起切掉，需要人工权衡。

5.3 坐标系不一致：经纬度影像直接算面积和距离会偏差

现象：同一个研究区，用 WGS84 和用 UTM 算出来的灯光面积差了 20% 以上。

原因：经纬度影像每个像元的实际地面面积在不同纬度不一样，直接统计像元数再乘以固定常数，数字自然不对。解决：所有统计分析统一在投影坐标系下进行。城市级研究选对应 UTM 带号，省域或全国范围选 Albers 等积投影。提取和统计全部基于投影后的影像，不要在经纬度影像上做面积和距离判断。

5.4 大范围处理内存溢出：分块读取与提前裁剪

现象：处理全省或全国范围影像时，直接src.read(1)把数组读进内存，程序直接崩掉。

原因：全球 VIIRS 年度合成是几万列乘几万行的浮点数组，一次性载入要占数 GB 内存，普通开发机扛不住。解决：用 rasterio 窗口读取，循环分块处理，每块算完立即释放。更省事的做法是先用gdalwarp命令行完成裁剪和重采样，把影像缩小到研究区范围再进 Python。研究区范围大时，输出也不要攒成一个大数组，每算完一个分块就把统计结果追加到文件。

5.5 跨传感器时间序列不可比：两个产品的数值口径

现象：把 2012-2014 年 DMSP 和 VIIRS 数据拼在一起做趋势，2013 到 2014 年灯光总量突然跳变。

原因：两种传感器的空间分辨率、量化位数、数值含义完全不同，直接拼接等于把两把不同刻度的尺子接在一起用。解决：要么全部用 VIIRS 做近十年统计，要么以 VIIRS 为基准对 DMSP 做回归校正后再拼接。常用做法是取 2012-2013 重叠年份的两套数据做幂函数拟合，用拟合系数把 DMSP 的 DN 值转换到 VIIRS 量纲。拟合样本要覆盖从暗到亮的所有亮度段，否则校正曲线在高端会失真。

6. 结果验证与跨期校正：拿统计数据回测灯光强度的三个动作

6.1 回归验证、目视检查与突变审查

强度算完不经过验证就写进报告，是分析翻车的根源。我验证一个夜间灯光结果通常做三件事：第一，把灯光总量和统计年鉴里的 GDP 或全社会用电量做线性回归，R² 高于 0.7 说明整体趋势靠谱，低于这个数就得回头查提取环节；第二，把灯光掩膜叠加到同期高分辨率遥感影像上，目视检查亮斑边界是不是落在可见城区范围内；第三，对时间序列数据，检查相邻年份的灯光总量突变，某一年突然翻倍，优先怀疑数据源更换或阈值漂移，而不是经济活动真的暴涨。

6.2 固定阈值的年度序列处理习惯与跨传感器校正

批量做多年份影像时，最忌讳每一年重新选一次阈值。我的习惯是先用三年影像做阈值稳定性测试，如果 Otsu 阈值在这三年里都稳定在同一个区间，就用这三年的均值作为整个时间序列的固定阈值。每跑完一年，记录影像版本、阈值、后处理参数和处理日期，存成 CSV 处理日志。这个习惯帮我拦住了无数次“第二年结果突然对不上”的排查。

这几年我越来越确认一个教训：夜间灯光分析里，强度计算永远比提取简单，真正决定结果质量的都是提取前后的细节——数据源选没选对、投影有没有统一、阈值有没有漂移、噪声有没有清干净。把这些流程固化下来，比换更复杂的算法有用得多。希望帮到你。

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