去除DMSP/OLS影像模糊的RTSVD算法

于志文，任 超,2，邓开元，毕旋旋

(1. 桂林理工大学测绘地理信息学院，广西 桂林 541004; 2. 广西空间信息与测绘重点实验室，广西 桂林 541004)

DMSP/OLS(the defense meteorological satellite program operational linescan system)影像是一种应用广泛的夜光遥感数据，在社会经济数据挖掘方面具有重要价值。该影像容易获取且可重复监测，但不当的数据采集与管理机制导致影像模糊，一般通过选取阈值或结合其他数据来抑制模糊。文献[1]提出89%累计发光阈值滤波，有效缓解了灯光溢出问题,但仅针对城市面积提取，且丢失了很多细节信息。文献[2—3]通过对DMSP/OLS影像设定3种不同阈值来获取城市边界，得到了比较精确的城市边界，但是处理过程过于烦琐。文献[4]利用DMSP/OLS辐射定标夜间灯光强度数据校正,有效扩展了强度值的动态范围，但并未改善模糊现象。文献[5]将NDVI和地表温度结合起来，在削弱城市中心DN值饱和的同时缓解了模糊现象，但辅助数据需经过复杂的计算。文献[6]提出了Pct影像滤波方法。Pct影像容易获取且结合Pct影像滤波的去模糊处理方法便于实现，但仅依靠Pct影像滤波不足以消除灯光边缘部分的模糊及噪声，还需结合其他方法。

综上，本文提出一种结合Pct影像滤波处理的正则化截断奇异值分解(RTSVD)去模糊算法。首先，根据光源像素发光频率必定高于非光源像素的规律，利用Pct数据比较像素发光频率，关闭图像中所有非光源像素，并设定15%发光频率阈值去除偶然影响。然后，利用高斯分布拟合点扩散函数(PSF)，并用自反边界条件[7]构造模糊矩阵。最后，绘制L曲线计算截断参数，采用RTSVD法对图像进行去模糊处理。

1 DMSP/OLS数据介绍

DMSP/OLS数据始于20世纪70年代美国国防部开展的美国军事气象卫星计划。DMSP卫星系统由6颗卫星组成，采用双星运行机制，两颗卫星同时运行。卫星运行在高度为833 km的太阳同步轨道上，运行周期为101 min，扫描宽幅为3000 km，每天收集14个轨道信息，可提供4次全球覆盖。OLS传感器是DMSP卫星的主要传感器，拥有可见光近红外(VNIR:0.4～1.1 μm)和热红外(TIR:10～13 μm)两个通道[8]，近地点分辨率可达500 m，扫描带边缘分辨率大于1 km。

本文以平均灯光影像和Pct影像两种DMSP/OLS数据为研究对象，数据来源于美国国家地球物理数据中心。DMSP/OLS非辐射定标夜间灯光强度数据包括无云观测频数影像、平均灯光影像、稳定灯光影像及Pct影像。其中，平均灯光影像是未经滤波处理的可见光强度平均值影像，是模糊问题的主要研究对象；稳定灯光影像是去除偶然噪声且经过降噪滤波的影像，图像质量优于平均灯光影像，可以用于影像去模糊效果评估；Pct影像是无云观测强度与灯光探测百分比的乘积，可用于监测灯光的发光频率。

2 平均灯光影像模糊原因

2.1 数据采集部分

DMSP卫星正下方的视场是正圆[8]，扫描运动使视场拉伸为椭圆。位于视场质心处的精细像素大小为0.56 km×0.56 km，远小于视场面积。在视场移动过程中，同一固定光源会落在多个重叠的视场范围内，形成以光源为中心的视场范围大小的模糊椭圆，并且由于卫星误算，模糊椭圆的质心存在随机偏移，导致DMSP影像中光源的地理位置偏移。文献[8]通过在全球不同地区移动便携式光源，测算出每夜DMSP影像中光源位置与GPS观测位置的距离与方位，所测的28个点全部向北偏移，偏移量均值为2.9 km,大约为一个平滑像素的宽度。文献[5]进一步推测出偏移的标准差为1～1.12 km。

2.2 数据管理部分

DMSP卫星的技术有限，不能存储浮点型小数，因此每个精细像素都是取整后的结果。同时，DMSP卫星也不能直接存储精细像素的高分辨率信息，因此存储时对精细像素进行平滑，把5×5像素的精细像素作为一个平滑像素。但DMSP卫星对精细像素平滑时每夜的起始计算位置并不统一，导致每夜的平滑像素图像不会完全重叠，这也是年平均影像模糊的原因之一。

3 基于Pct影像滤波的RTSVD算法

3.1 Pct影像滤波

Pct影像是由无云探测的平均DN值乘以探测百分比频率得到的，可以代表灯光的百分比频率。依据灯光影像成像原理，光源像素的发光频率必定高于非光源像素，可利用Pct影像关闭发光频率小于周围像素的像元，以去除非光源像素误差。

首先筛选出像素xi,j周围发光频率最高的像素xmax，比较xi,j与xmax为

(1)

把经过式(1)处理后的Pct影像与模糊影像相乘，得到光源像素影像。但Pct影像滤波不能完全关闭灯光边缘部分的模糊像素，而且无法去除影像中的噪声。因此，本文在Pct影像滤波的基础上提出RTSVD算法进一步处理平均灯光影像。

3.2 RTSVD算法

RTSVD算法是融合了正则化思想的截断奇异值分解法，目的在于去除影像中的噪声并恢复图像质量。对于m行n列的模糊图像，一般有线性模糊模型

Ax=b

(2)

式中，x、b分别为真实图像与模糊图像排列成的列向量，长度为N=mn；A为模糊矩阵，A∈RN×N。由此可以得出基础去模糊模型

x=A-1b

(3)

模糊矩阵A由点扩散函数与边界条件两个成分决定。其中，点扩散函数可以利用高斯分布来拟合[7，9]。根据模糊原因，模糊矩阵A的边界条件设为自反边界条件[10]。然后，利用奇异值分解进一步完善去模糊模型。对于m×n阶模糊矩阵A，存在一个分解使得

A=USVT

(4)

式中，U和V为正交矩阵，分别为左奇异向量阵和右奇异向量阵，满足UTU=Im，VTV=In；S为奇异值矩阵，为对角矩阵，σ1、σ2、…、σm为对角线上的奇异值，满足

σ1≥σ2≥…≥σm≥0

(5)

式(3)可写成

(6)

式中，ui为左奇异向量；vi为右奇异向量。模糊矩阵A为病态矩阵，条件数cond(A)=σ1/σn非常大，奇异值衰减的速度非常快，如图1所示。

与小奇异值对应的误差分量被放大，干扰图像真实信息，公式描述如下[11]

(7)

式中，Δ为误差。为了削弱误差Δ的干扰，运用RTSVD算法，用截取的奇异值矩阵Sk来替代原奇异值矩阵Sm，公式如下

(8)

式中，k为截断参数。RTSVD算法的误差来源于正则化与噪声两个方面，在选择截断参数时，应考虑在抑制噪声的同时尽可能缩小正则化误差。

3.3 L曲线

截断参数的选择对于RTSVD算法至关重要，一般采用GCV或L曲线来确定截断参数。但针对图像模糊问题时，GCV难以计算出截断误差的最小值。因此，本文采用L曲线确定截断参数，利用三次样条曲线拟合离散点。L曲线是以log-log为尺度，在以正则化解范数为纵轴，对应的残差范数为横轴的直角坐标系中，对一系列截断参数绘制的曲线[12]。L曲线上曲率最大的点对应的截断参数为最优解。令ρ、η分别为残差与解的对数，L曲线的曲率K公式表达为

(9)

在截断奇异值分解中，L曲线是由若干离散点组成的，要计算出最优解，必须在保证曲线整体形状不变的条件下定义一条与离散点相关联的可微、光滑的曲线。因此，采用三次样条曲线逼近L曲线离散点,计算连续曲线上曲率最大的点从而确定截断参数。

4 算例分析

为了验证RTSVD算法对DMSP图像去模糊的有效性,裁剪大小为417×223像素的台湾地区F142000期平均灯光影像进行去模糊试验。为了评价本文方法效果，选取稳定灯光影像和89%累计发光阈值滤波影像与试验结果进行比较,如图2所示。

除了对试验结果进行目视观察评价外，引入峰值信噪比(PSNR)、图像信息熵、灰度平均梯度(GMG)与边缘强度作为客观评价标准。由表1可以看出，原始的平均灯光影像(图2(a))存在严重的模糊现象。依据L曲线选取截断参数k=67，处理得到图像(图2(d))与原始影像相比，峰值信噪比增大2.505 3，信息熵增加0.197 4，灰度平均梯度增加1.023 9，边缘强度增大9.855 3，与稳定灯光影像相比，各项指标均有所增长。对于直接频率滤波处理的影像(图2(c))，峰值信噪比提高3.443 5，信息熵降低0.482 8，表明该方法虽然有效降低了噪声与模糊，但严重丢失了影像的细节信息。这也说明结合Pct影像滤波的RTSVD算法在消除模糊的同时保留了影像的细节信息。

表1 去模糊效果评价

5 结 论

本文采用结合Pct影像滤波的RTSVD算法来解决DMSP/OLS数据的模糊问题，通过F142000期台湾地区的平均灯光数据试验，得出以下结论：

(1) 通过比较像素的发光频率可以去除非光源像素。

(2) 利用样条曲线拟合L曲线离散点可以计算出RTSVD的截断参数。

(3) RTSVD算法在去除模糊的同时能够保留影像的细节信息。