简介：针对GIS领域夜间灯光提取与强度计算的完整操作指南，依托ArcGIS平台，面向需要处理遥感栅格与行政矢量数据的城市规划、环境监测、经济发展评估等方向学习者。内容覆盖VIIRS/DMSP-OLS夜间灯光数据获取、坐标统一、按属性选择行政区域、按掩膜提取、分区统计以及符号化制图全流程，并以太原地区为例演示具体参数设置与连接汇总方法。资源包为单个docx文档，整体6.52MB，文档内以步骤化方式组织，含菜单路径、工具位置和字段选择细节，便于跟随实操和排查属性表连接时的常见问题。已有3651人学习下载，适合希望系统掌握ArcGIS夜间灯光数据处理流程、提升空间分析能力的研究生或相关从业者。通过分区统计得到的均值、总和等灯光强度指标，可直接用于城镇化对比、经济活跃度分析等后续研究，能帮助读者少走弯路并快速复用整套操作模板。
1. 夜间灯光提取与强度计算：城市研究者手里的"透视镜"
夜间灯光提取与强度计算，说白了就是从遥感影像里把灯点亮的地方准确抠出来，再把这些像元的亮度值换算成可比较的强度指标。很多人第一次拿到夜间灯光影像时，觉得它不过是一张"世界夜景图"——亮点越亮就是城市越发达。可一旦你真去处理，就会发现提取这一步就够呛：月光云层反射、传感器噪声、饱和像元、不同年份亮度不可比，每一条都能让结果翻车。这篇笔记要讲的，正是从原始夜间灯光影像到可信强度指标的完整链路：数据源怎么选、阈值怎么定、强度怎么算、坑在哪里。适合遥感数据处理者、城市规划与区域经济研究者，以及所有想用夜间灯光数据做实证分析的从业者。
2. 先说数据源：不同夜间灯光影像的适用边界
夜间灯光数据不是只有一种。不同代际的传感器在空间分辨率、辐射动态范围和时间跨度上有本质差异，选错数据源后面再精细的提取都是白搭。我一般把数据源分成三类，按需求匹配。
2.1 三代夜间灯光数据：分辨率、动态范围和适用场景
第一代夜间灯光数据从 1990 年代开始积累，空间分辨率约 1 公里，特点是时间序列长，能回溯二十多年建成区扩张。但它的辐射动态范围很窄，城市中心会迅速饱和——所有像元都顶到最大值，亮度差异完全消失。做长时序趋势分析时它几乎不可替代，但拿它做城市内部强度差异研究，结果会严重失真。
第二代产品大约从 2011 年前后开始，分辨率提升到约 500 米，最大的改进是不再饱和。城市核心区和外围的亮度差异能被记录，动态范围提升一到两个数量级。代价是时间序列不够长，且数据中带有一层微弱的非灯光背景信号，需要额外去噪。适合做经济空间估算、逐级亮度分析和区域对比。
第三类是近几年出现的高分辨率夜间灯光数据，分辨率达到几十米级，能看到道路灯、建筑轮廓甚至单栋楼的光。细节丰富，但覆盖范围小、时间不连续，获取门槛也高。适合街区级光环境分析和局地细查，不适合做宏观区域比较。
| 数据代际 | 空间分辨率 | 辐射动态范围 | 时间跨度 | 短板 | 适合场景 |
|---|---|---|---|---|---|
| 第一代 | 约 1km | 窄，易饱和 | 约 20 年以上 | 城市中心饱和 | 长时序建成区扩张 |
| 第二代 | 约 500m | 宽，无封顶 | 约 10 年以上 | 带背景噪声 | 经济空间估算、强度分级 |
| 高分辨率 | 几十米 | 较宽 | 短且间断 | 覆盖小、获取贵 | 街区级照明分析 |
选型时先问自己一个问题：你要回答的是"这个区域二十年来怎么变的"，还是"这个城市内部哪里更强、哪里更弱"？前者选第一代，后者选第二代。硬拿第一代做强度分层，城市中心一片全白，分析就失去了意义。反过来，拿第二代数据做长时序，序列长度不够，结论的置信度也撑不住。
2.2 数据预处理：投影统一、范围裁剪、重采样
拿到影像后，第一步不是急着提取灯光，而是把数据整理到统一的空间基准上。夜间灯光影像的原始投影通常是地理坐标，也就是经纬度网格。经纬度下像元面积随纬度变化，直接计算面积和密度必然失真。我一般先转成等积投影，国内区域常用 Albers 等积圆锥投影，再按研究区边界裁剪，最后统一重采样像元大小。
# 1. 统一投影到 Albers 等积投影（国内常用参数） gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84" \ input_nightlight.tif nightlight_aea.tif -r bilinear # 2. 按研究区矢量边界裁剪 gdalwarp -cutline study_area.shp -crop_to_cutline \ nightlight_aea.tif nightlight_city.tif # 3. 统一像元大小（例如重采样到 500m） gdalwarp -tr 500 500 -r bilinear nightlight_city.tif nightlight_city_500m.tif三个命令对应三个预处理动作。第一条里的-r bilinear是重采样方法，双线性内插会保留一定的亮度渐变信息，适合连续型辐亮度数据。如果你的数据是分类结果或掩膜，应该改用-r near最近邻，避免插值造出虚假类别。第二条用-cutline加-crop_to_cutline裁剪，比手动算行列号再切片更稳妥，矢量边界在哪，影像就切到哪。第三条把像元统一成 500m 时，要注意-tr 500 500的单位和投影一致。投影到以米为单位的等积投影后，这里就是 500 米。
预处理还有个容易被忽略的点：重采样会改变 DN 值的分布。双线性内插会对相邻像元做加权平均，相当于做了一次平滑，直方图会比原始数据更集中在中间。如果你的后续提取使用固定阈值，重采样方法的选择会直接影响提取结果。我通常先把预处理做完整，然后停下来看一眼直方图，确认重采样没有把双峰结构抹平，再继续往下做。
3. 灯光提取的核心操作：阈值、去噪与像元筛选
数据准备好后，真正的"提取"才开始。灯光提取本质上是一个图像分割问题：把影像像元分成"灯光"和"非灯光"两类。看起来简单，但夜间灯光影像没有清晰的边界，灯光强度从城市中心向外围连续衰减，再加上背景噪声，分割阈值怎么定就成了最玄学的一步。
3.1 为什么阈值是关键：灯光的直方图形态与一刀切问题
打开一景夜间灯光影像的直方图，通常会看到这样的形态：大量像元集中在接近零的低值区，这是非灯光背景；一小部分像元拖出长长的尾巴，这是灯光像元。理想情况下直方图呈现双峰，两个峰之间的谷底就是最佳分割点。但实际数据几乎不会这么干净——背景峰很宽，灯光尾巴又和它交叠，谷底不明显。
固定阈值是最常用的方法，设定一个 DN 值，高于它的算灯光。问题在于，城市边缘的灯光强度可能只有城市中心的几个百分点，一个固定阈值要么把外围弱灯光全部扔掉，要么把背景噪声全放进来。这就是为什么有人说阈值设定是"玄学"——同一个阈值，换一景影像、换一个城市、换一代传感器，结果都可能完全不同。
自适应阈值是更可靠的做法。常见的有 Otsu 大津法：在直方图上遍历候选阈值，计算每个阈值下背景和灯光两类的类间方差，取类间方差最大的那个值作为分割点。它的逻辑是"灯光和非灯光的区别最大化的分界线"，不依赖人眼判断。Otsu 不是万能的，如果背景噪声很强、双峰结构不明显，它会把阈值抬得偏高，损失边缘灯光。所以我的习惯是：用 Otsu 求出初值，再结合影像空间信息做人工复核，必要时按区域做局部阈值。总之，阈值不是随手填的数字，它决定了你后续所有强度统计的底层边界。
3.2 从影像到掩膜的完整流程：阈值、去噪、连通域清理
下面这个函数是我做单景灯光提取的固定流程，依次处理无效值、背景去噪、阈值分割、形态学清理和连通域筛选。可以直接复制运行，替换成本地文件路径即可。
import rasterio import numpy as np from scipy import ndimage def otsu_threshold_from_img(img, bins=200): """对影像中大于 0 的像元做 Otsu，返回阈值""" data = img[img > 0].astype(np.float32) if data.size == 0: return 1.0 hist, edges = np.histogram(data, bins=bins) total = hist.sum() if total == 0: return float(data.mean()) sum_all = float(np.sum(hist * edges[:-1])) sum_b = 0.0 w_b = 0.0 best_t = edges[0] best_var = -1.0 for i in range(len(hist)): w_b += hist[i] if w_b == 0: continue w_f = total - w_b if w_f == 0: break sum_b += float(hist[i] * edges[i]) m_b = sum_b / w_b m_f = (sum_all - sum_b) / w_f var = w_b * w_f * (m_b - m_f) ** 2 if var > best_var: best_var = var best_t = edges[i] return float(best_t) def extract_light_strength(src_path, threshold=None, min_area=30, kernel=3): with rasterio.open(src_path) as src: img = src.read(1).astype(np.float32) profile = src.profile nodata = src.nodata # 1. 无效值清理：nodata 和负值统一置 0 if nodata is not None: img[img == nodata] = 0 img = np.nan_to_num(img, nan=0.0, posinf=0.0, neginf=0.0) img[img < 0] = 0 # 2. 背景噪声抑制：3x3 中值滤波 img_f = ndimage.median_filter(img, size=kernel) # 3. 阈值分割：未指定阈值时用 Otsu 自动求 if threshold is None: threshold = otsu_threshold_from_img(img_f) mask = img_f >= threshold # 4. 形态学开运算，去掉零散亮点 mask = ndimage.binary_opening(mask, structure=np.ones((3, 3))) # 5. 删除过小连通域（孤立亮点） labeled, num = ndimage.label(mask) sizes = ndimage.sum(mask, labeled, range(num + 1)) mask = np.isin(labeled, np.where(sizes >= min_area)[0]) # 6. 强度统计 light_values = img * mask n_pixels = int(mask.sum()) total_strength = float(light_values.sum()) mean_strength = total_strength / n_pixels if n_pixels > 0 else 0.0 light_area_ratio = n_pixels / img.size # 7. 输出掩膜结果 out_profile = profile.copy() out_profile.update(dtype=rasterio.uint8, count=1) with rasterio.open("light_mask.tif", "w", **out_profile) as dst: dst.write(mask.astype(np.uint8) * 255, 1) stats = { "threshold": float(threshold), "n_pixels": n_pixels, "total_strength": total_strength, "mean_strength": mean_strength, "light_area_ratio": light_area_ratio, } return mask, stats # 运行示例 mask, stats = extract_light_strength("nightlight_city_500m.tif") print(stats)逻辑说明：第 1 步把无效数据清零，避免它们干扰统计；第 2 步中值滤波是专门对付椒盐式背景噪声的，它能去掉孤立的亮暗点，同时不过度模糊灯光区域的边界；第 3 步是关键分类决策，Otsu 自动求得的阈值通常落在背景和灯光两个峰之间的谷底附近；第 4 步的开运算是先腐蚀再膨胀，用来切断灯光区域之间细小的虚假连接；第 5 步的连通域面积过滤删除那些面积小于阈值的孤立亮点，通常是噪声或偶发光源。
参数说明：kernel=3代表中值滤波窗口为 3×3 像元，窗口越大去噪越强但边缘越模糊，一般不要超过 5。min_area=30是最小连通域面积（像元个数），在 500 米分辨率下约等于 7.5 平方公里，在 1 公里分辨率下约等于 30 平方公里，需要根据你的数据分辨率调整。bins=200控制 Otsu 直方图分箱数，分箱太少阈值精度差，太多会被噪声带偏。还有一点需要注意：如果某景影像里城市之外还有大型工业火光或油气田燃烧，这些亮区会被算法当成灯光保留下来。做提取前最好人工查看一次影像，把这类非城市灯光先用掩膜排除掉。
4. 强度计算：从像元 DN 到区域灯光指数
提取出灯光掩膜之后，下一步就是"强度计算"。很多刚上手的人以为强度就是"把亮度值加起来"，但严格来说，DN 值本身只是一个无量纲的传感器响应等级，不同数据、不同年份之间不能直接比。要做可靠的强度分析，得先把概念理清楚。
4.1 DN、辐亮度与灯光强度指标：三者的关系
DN 值是影像原始记录值，数字越大代表传感器接收到的辐射越强。它受传感器增益、定标状态和大气条件影响，同一盏灯用不同传感器拍出来 DN 值可以差一倍。辐亮度则是把 DN 值经过辐射定标后的物理量，单位是 W/(m²·sr·μm)，表示地表物体实际向外辐射的能量。这才是可以做物理意义上比较的值。
数据处理中，我见过最多的误用是"拿 DN 直接比大小"。同一景影像内部比，没问题；跨年份、跨传感器比，必须做定标或归一化。如果你的数据源提供定标系数，常见形式是辐亮度 = a × DN 的幂函数，把系数代入即可。如果不提供系数，我的建议是：只做同一数据内部的相对比较，不要声称亮度绝对值。
基于提取的灯光掩膜和像元强度，常用三个区域聚合指标：
| 指标 | 计算方式 | 含义 |
|---|---|---|
| 总灯光强度 | 区域内灯光像元 DN 值求和 | 反映区域灯光总量 |
| 平均灯光强度 | 总强度 / 灯光像元数 | 反映灯光集中程度 |
| 灯光面积比 | 灯光像元面积 / 区域总面积 | 反映灯光覆盖范围 |
总灯光强度最常用，但它对面积敏感，面积大的区域天然数值更大。平均灯光强度消除了面积影响，反映的是"亮得有多集中"。灯光面积比则回答"灯光覆盖了多大范围"。三种指标组合使用，才能说清楚一个区域的灯光特征。
4.2 区域聚合统计：从像元强度到面状指标
有了灯光强度栅格和一个行政区边界，最常见的需求是：给每个行政区算出上面的三个指标。下面这段代码用矢量栅格化的方式做分区统计，比逐区域裁剪循环快得多。
import numpy as np import geopandas as gpd import rasterio from rasterio.features import rasterize def zonal_light_stats(region_shapefile, light_tif, id_col="id"): """按行政区统计灯光强度指标，id_col 必须是整数列""" with rasterio.open(light_tif) as src: img = src.read(1).astype(np.float32) transform = src.transform crs = src.crs # 1. 读取矢量并转为栅格坐标系 gdf = gpd.read_file(region_shapefile).to_crs(crs) shapes = [(row.geometry, int(row[id_col])) for _, row in gdf.iterrows()] # 2. 矢量栅格化：每个像元值 = 其所在区域的 id zone_raster = rasterize( shapes, out_shape=img.shape, transform=transform, fill=0 ) # 3. 按区域统计 results = [] for z in np.unique(zone_raster): if z == 0: continue zone_mask = zone_raster == z vals = img[zone_mask] light_mask = vals > 0 total_light = float(vals[light_mask].sum()) light_count = int(light_mask.sum()) pixel_area_m2 = abs(transform.a * transform.e) light_area_km2 = light_count * pixel_area_m2 / 1e6 results.append({ id_col: int(z), "total_light": total_light, "mean_light": total_light / light_count if light_count > 0 else 0.0, "light_area_km2": light_area_km2, }) return results stats = zonal_light_stats("city_boundary.shp", "nightlight_city_500m.tif") print(stats)逻辑说明：第 1 步把行政区矢量对齐到灯光影像的坐标系，第 2 步用rasterize把矢量转成 ID 栅格，每个像元的值为所在区域的编号，第 3 步利用布尔索引一次性取出每个区域的灯光像元值做统计。这里不再循环裁剪，速度快很多。
注意几个参数：id_col必须是整数列，如果原始数据里是字符串编号，先gdf[id_col] = gdf[id_col].astype(int)再传入。pixel_area_m2 = abs(transform.a * transform.e)的前提是影像已经是等积投影且单位为米，假如这一步算错了，后面的灯光面积和密度都会跟着错。vals > 0的前提是强度栅格中非灯光像元已经是 0，如果你在提取环节没有把背景清零，这里就需要先做一次阈值判断。
4.3 多时相可比性处理：简单交叉定标
多时相分析是夜光数据的核心应用场景，但不同年份影像之间往往存在系统性亮度偏移。直接比较会出现"所有像元同步变大或变小"，这会让你误判为城市在扩张或收缩。交叉定标的标准做法是用伪不变像元建立年份之间的回归关系：找一批在一段时间内灯光强度稳定不变的像元，比如大型机场、长期稳定运行的工业区，在基准年份影像和待校正影像上分别取值，然后做线性回归。
import numpy as np def cross_calibrate(img_target, img_ref, stable_mask): """用稳定灯光像元做线性校正，返回校正后影像和系数""" x = img_ref[stable_mask].flatten() y = img_target[stable_mask].flatten() valid = np.isfinite(x) & np.isfinite(y) & (x > 0) & (y > 0) x, y = x[valid], y[valid] # 拟合 y = k*x + b A = np.vstack([x, np.ones_like(x)]).T k, b = np.linalg.lstsq(A, y, rcond=None)[0] corrected = k * img_target + b return corrected, k, b # 示意调用 # corrected_img, k, b = cross_calibrate(img_2020, img_2015, stable_mask)逻辑说明：stable_mask是关键输入，它标记哪些像元是"稳定"的。一个简单可靠的构建方法是取三期影像中始终超过某个高阈值的像元，这些像元大概率是持续高强度光源。拟合出的k和b反映了待校正年份相对基准年份的系统性偏移。校正后再去做强度比较，结果才是可信的。
没有稳定的像元样本怎么办？一个后备方案是取区域总体均值做比例归一化：把所有年份影像按某一基准年份的均值缩放到相同尺度。这种方法的精度低一些，但至少能消除整体增益差异。
5. 避坑指南：夜间灯光处理中最容易翻车的 5 个细节
这部分是我踩过最多坑的地方。每一条都是实际项目中反复出现的问题，现象、原因、解决一次说清。
5.1 提取结果全是孤点噪声：背景没有被压住
现象是提取出的灯光掩膜上布满密密麻麻的孤立亮点，完全没有城市轮廓，更像一张撒满芝麻的饼。原因是新一代夜间灯光影像存在微弱背景信号，包括月光云层反射、大气辉光和传感器热噪声，它们本身就有一定 DN 值，直接做阈值分割时被当成灯光放进来。解决方法是先做背景抑制：用 3×3 或 5×5 中值滤波把孤立的噪声点压掉，再进入阈值分割。如果中值滤波还不够，就提高阈值，观察直方图中背景峰的右侧边界，把阈值设到贴着背景峰边缘的位置。
5.2 城市灯光连不成片：阈值定太高
现象是城市中心的掩膜是完整的，但到了城市外围，原本连续的灯光带断成一截一截的碎块，看起来像切碎的香肠。原因是全局阈值对城市边缘不友好——外围灯光强度低，离城市中心越远越低，被同一个高阈值切掉后就产生了断裂。解决方法是改用局部自适应阈值，比如 Otsu 分区域计算，或者用滑动窗口局部阈值。如果只想快速修正，可以先把阈值调低 10% 到 20%，再做一次形态学闭运算，把小缺口填上。但要注意，闭运算也会把原本不连接的灯光区粘在一起，填缺口要适度。
5.3 强度算出来异常偏低：过饱和与无效值混入
现象是某区域总灯光强度数值明显不合理，城市核心甚至比外围还低。原因有两个方向：第一代数据在城市中心过饱和，所有像元顶在最大值，看不出差异，但平均值并非异常低；真正让数值低到离谱的通常是无效值没清理干净，比如 nodata 或负值被当成 0 参与平均，拉低了统计结果。解决方法是做强度统计前先检查数据：把 nodata、负值和异常大的值单独列出来看分布，确认无效值已经掩膜掉，再计算指标。对于第一代数据的饱和问题，只能在结论里说明这种低估，或者改用宽动态范围的数据源。
5.4 跨年份影像直接比大小：可比性被忽略
现象是某年灯光总量突然比相邻年份暴涨 30%，但该区域并没有大规模建设。原因是不同年份的影像可能来自不同传感器、不同增益设置或不同处理版本，整体亮度水平本来就不一致。解决方法是先做交叉定标，用我第 4 章的cross_calibrate函数把各年份数据校正到同一基准，再比较。如果你不想做完整定标，至少要算一下各年份影像在稳定像元上的均值，确认没有系统性偏移后再下结论。
5.5 跨区域统计面积失真：投影与像元面积没处理好
现象是不同纬度区域的灯光面积对比结果不合理，纬度越高显得面积越大。原因是原始地理坐标下像元面积不恒定，纬度越高一个像元对应的实际地面面积越小，如果按像元个数当面积用，当然会失真。解决方法是全程使用等积投影。另外，矢量栅格化时边界像元会被栅格化网格截断，一个行政区边缘的灯光像元可能被划给相邻区，导致面积和强度统计产生偏差。对边界比较敏感的话，可以在统计时加一层缓冲区，或者用更高分辨率的栅格做统计后再聚合。
6. 验证与进阶：让灯光强度真正服务于你的分析
算完强度指标，事情还不算完。一个没有人验证过的提取结果，不能直接进入论文或决策，至少要做一次相关性检验：把你的区域灯光强度数据和该区域可获取的独立指标对接，比如用电量、生产总值、常住人口，算一下皮尔逊相关系数。相关系数在 0.8 以上说明提取结果基本可信；低于 0.6 就要回头检查提取流程——大概率是阈值或去噪出了问题，而不是数据本身不行。
进阶方向上有三个做法值得尝试。第一个是做局部自适应阈值替代全局阈值：from skimage.filters import threshold_local; local_thresh = threshold_local(img_f, block_size=15, method="gaussian")，用滑动窗口在每个局部区域单独求阈值，城市边缘弱灯光会被保住。第二个是构建年度合成影像：把一年内多期夜间灯光数据取最大值或均值，消除云遮挡和临时光源干扰，得到更稳定的灯光基准。第三个是把灯光掩膜叠加路网或建成区边界做提纯，去掉农田里的温室灯光、海上渔排灯光等非城市光源。
我最早做灯光提取时，死活调不出一份稳定的城市掩膜，后来发现自己的问题是拿一个固定阈值打天下，换了传感器数据就翻车。现在我的习惯是先看直方图，再让 Otsu 给初值，最后用城市边界和路网做空间复核，三步缺一不可。希望帮到你。
本文还有配套的精品资源，点击获取
