夜光遥感被公认为大范围同步监测社会经济活动的良好数据来源，已在人口估算、经济估计、城市化监测、环境与能源消费研究等领域得到广泛应用^[1-8]。与美国国防气象卫星(defense meteorological satellite program，DMSP)搭载的业务型线性扫描传感器(operational linescan system，OLS)获取的数据相比，美国国家极地轨道合作卫星(national polar-orbiting partnership，NPP)搭载的可见光红外成像辐射仪(visible infrared imaging radiometer suite，VIIRS)获取的日夜波段(day and night band，DNB)继承并优化了前者的微光探测能力，消除了像元饱和，减少了溢出现象^[9]。由于该传感器对夜间微光的优质探测灵敏度，影像易受噪声干扰，影响了应用的准确性。目前，针对NPP/VIIRS影像校正的研究尚不多见，大多停留在掩膜去噪^[10]、固定阈值分割^[11]等操作。利用DMSP/OLS掩膜剔除VIIRS噪声的方法仅适用于2012和2013年，且随着影像更新，掩膜去噪的弊端越来越明显；固定阈值分割法主观性强且通用性差，分割效果不稳定。此外，国内外针对NPP/VIIRS的校正和研究侧重于宏观尺度的年际数据，影像月度发布的优势很少体现。

本文在研究NPP/VIIRS月度影像辐射特征的基础上，提出多时相夜光遥感影像校正方法。该方法基于无光源区域样本剔除背景噪声，使用峰值特征化放大辐射差异，通过卷积阈值分离异常像元与稳定像元；对于卷积阈值边界附近的像元，按最大熵阈值和八邻域中值滤波重分类，实现对异常像元的动态剔除。

1 NPP/VIIRS夜光遥感影像数据分析

美国国家地球物理数据中心利用NPP/VIIRS的DNB波段获取的夜间灯光数据生产了月度影像集，提供像元辐射值表征灯光强度，灯光辐射值的大小反映灯光的强弱，单位为nW·cm^-2·sr^-1，像元辐射值0值及以下视为无灯光。使用2012年4月至2017年12月的69期NPP/VIIRS VCMCFG影像产品，以江苏省为研究区，行政区划矢量数据来自国家基础地理信息中心，GDP数据来自国家统计局。

由于影像显著提高了清晰度，火光、渔船和油气井等短暂光源及河流湖泊等反射弱光极易对数据造成局部扰动。表 1中，背景噪声表现为影像中存在的部分低辐射及负值像元，依据“无灯光区域数值为0”的原则^[12]，辐射值应为0；异常像元则表现为像元辐射值在相邻月份影像中的无序波动，2016年8—10月的最大灯光辐射值及辐射总值差距较大，根据近几年研究区持续稳定的发展状况和无重大灾害的自然情况，不应存在大量像元的异常波动。背景噪声和异常像元的存在影响了夜光数据应用的准确性，需对其进行校正。

2 NPP/VIIRS数据校正方法

选取研究区69期月度夜光遥感数据，进行预处理、背景像元校正及异常像元校正等处理，流程如图 1所示。

将月度数据作为待校正影像，编号M_ij(i为年份，j为月份)，使用行政区划数据批量裁剪研究区。为避免形变，将影像重投影为Albers等积坐标系，重采样为500 m×500 m的网格。为降低影像的月变化敏感度，利用均值法对重采样后的M_ij求得季度图像Q_ik(i为年份，k为季度)和年度图像Y_i(i为年份)。引用Google Earth图像，从各地级市的Y_i中选择19处样本区域，类型包括湖泊、大型水库和农田等无光源区。在VIIRS影像中寻找并记录该位置的辐射值，计算样本辐射均值作为最小阈值，将最小阈值以下的像元统一替换为0值，达到去除背景噪声的目的。

2.1 峰值特征化卷积运算

NPP/VIIRS影像在获取时已进行星上辐射定标，影像间可比，因此相互校正的峰值特征化法可行。以季度图像Q_ik为参考，与M_ij作差得月度峰值特征图像M′_ij，异常像元在M′_ij中呈现为零散分布的高辐射亮点目标，稳定像元成为低辐射背景。

为明确区分M′_ij的目标与背景，在峰值特征图像区域中设置3×3的卷积图像函数T, i和j代表像元所在行与列，设置与T区域相同的卷积权窗口Q，如下

此处权矩阵模板为高通滤波器，可在保持高频信息的同时，弱化图像中的低频成分。将T与Q作加权求和得卷积值T′。卷积权窗口由左至右，由上至下顺次滑动一个像元直至遍历完所有像元，按式(2)求取每个单元的卷积值，产生M′_ijm的卷积运算结果M″_ijm

经过卷积处理，M″_ij中的目标形成相对独立的像元群，在图像中表现为卷积正值区，背景像为负值，二者明显区分。将M″_ij中的卷积负值区统一赋值为0，设置大于0的区域为掩膜，提取峰值特征图像M′_ij的高频信息区域。

2.2 基于最大熵阈值的自适应滤波

卷积运算放大了像元间的辐射差异，但在图像特征提取环节仍使用固定阈值，对卷积阈值边界附近的像元易造成漏警和虚警^[13]。熵是对随机变量不确定性的度量，若将影像的像元辐射值看作一组随机变量，则影像的熵为测量辐射级分布随机性的一种特征参数^[14]。图像分割时越靠近目标与背景的边界，熵越大，使目标与背景的熵值之和达到最大的分割阈值，记为最大熵阈值。相较于固定阈值，最大熵阈值综合考虑了像素点的辐射值和邻域分布信息，具有更好的分割效果和自适应性。设置阈值Thr将峰值特征图像的卷积结果M″_ij划分为目标A和背景B，其各像元概率分布如下

式中，Thr为阈值；i为像元辐射值，P_i为直方图中像元辐射值出现的概率；N(R)为图像中最小辐射值；M(R)为最大辐射值。其中A、B的概率P_A、P_B及其熵H_A、H_B的计算如下

图像的熵H_mm为目标A和背景B熵的和，如下

目标A和背景B每个概率分布对应的熵，记为分布信息熵H(A)与H(B)，如式(6)所示。使得目标与背景分布信息熵H(A)与H(B)之和最大的Thr为最大熵阈值。

利用搜索算法寻找每一幅M″_ijm的最大熵阈值并作二值化处理，生成0值与非0值掩膜，对M_ijm的0值掩膜区域保留，非0掩膜区域采用八邻域算法的中值滤波运算，去除图像中具有椒盐噪声特性的异常像元，并最大化保留原图像信息^[15]。两种掩膜的像元合并即完成了多时相NPP/VIIRS夜光遥感影像的校正工作，校正后部分月度影像示例如图 2所示。与校正前影像相比，去除了负值与低辐射的背景噪声，将月度数据的波动控制在合理范围。

2.3 影像的校正结果评定

为检验校正效果，构建表征整体辐射情况的灯光辐射总数(total number of radiation，TNR)及反映辐射分布的像元离散指数(pixel dispersion index，PDI)，在微观尺度上比较影像的变化趋势，TNR计算如下

式中，i为像元辐射值；N_i表示i辐射值对应的像元数。月度序列TNR变化情况如图 3所示。

由图 3可知，因数据缺失及噪声的存在，初始TNR随时间推移存在大幅度无序波动的现象；校正修补了因云层遮盖而缺失的影像，TNR的波动幅度缩小，相邻月度影像TNR差异平均下降了63.84%，数据呈现持续增长轨迹，符合研究区发展情况，实现了时间序列上影像变化更稳定的效果。

校正对异常像元的剔除也影响了辐射值与像元数的分布及变化。构建像元离散指数PDI表征研究区域像元分布的离散程度，计算公式为

式中，i为像元辐射值；i为影像辐射值的算术平均值；N_i表示i辐射值对应的像元数目。PDI变化情况见表 2，指数越小则像元越集聚。由于噪声影响，原始图像灯光辐射值的范围过大，数值离散，PDI在(6，10]的月份占89.85%；校正有效缩小了分布区间，PDI在(0，6]的月份占82.61%，像元分布呈现集聚状态，时间序列的像元离散指数平均降低26.98%，对地表照明分布的刻画更准确。

针对DMSP/OLS与社会经济参量进行宏观定量关联分析的研究表明，夜光数据和年度GDP数据存在显著的线性关系^[16-17]。由于掩膜去噪法不适用于2013年后的数据，本研究根据文献[11]的方法取0.30×10 nW·cm^-2·sr^-1作为固定阈值，得到江苏省13市2012—2016年的年度NPP/VIIRS校正数据。同时，基于本文的多时相夜光遥感影像校正方法(简称最大熵法)生产的月度影像求均得到年度数据。采用线性模型对两种方法模拟社会经济参量的潜力进行比较，由于数据值分布宽泛，对数据的显示和处理造成麻烦，统一采用对数变换处理，回归结果如图 4所示。

图 4(a)和图 4(b)分别表示固定阈值法和最大熵法校正的NPP/VIIRS年度数据与全市生产总值的回归结果，lnGDP表示季度GDP的对数，lnTNR表示校正灯光辐射总数的对数。回归结果表明，全市GDP与固定阈值法的R²值(0.909)低于其与最大熵法回归的R²值(0.936)。说明在宏观年度尺度上，本文方法生产的夜光数据在区域生产总值上的相关性高于固定阈值法，具备更大的模拟潜力。

3 结语

基于统计数据的定量分析和基于DMSP/OLS的定性分析已无法满足当前研究的需求，NPP/VIIRS数据以其影像清晰度、信息丰富度及发布频率克服了统计数据与DMSP/OLS的局限，可为多级区域尺度下的宏观动态监测提供更加强大的数据支撑。

多时相夜光遥感影像校正方法修补了缺失影像，实现对短暂光源和背景噪声的有效剔除，修正了灯光辐射及像元分布，微观上夜光数据的稳定程度和集聚状态均有提升，宏观上模拟社会经济参量更具优势，能满足月、季、年等多种时间分辨率的研究需求，为了解区域社会经济活动及城市化进程提供数据基础，为相关单位提供决策支持。