一种基于L0范数平滑滤波的遥感影像底图匀色算法

赵薇薇， 王 艳， 丁一帆， 薛晓伟， 孙 赜

(北京市遥感信息研究所， 北京 100192)

0 引言

遥感影像镶嵌产品被广泛应用于城市规划、建设、管理、应急响应等多个场景中[1]。镶嵌产品生产时，需要将区域内不同来源、不同时间获取的影像进行平差校正和色彩处理，如图1所示为区域镶嵌产品生产过程中色彩处理前后的对比。色彩一致化在遥感影像镶嵌产品中扮演重要的角色。

图1 色彩处理前后对比(a)色彩处理前的安徽卫星遥感数据，由快鸟[2]卫星拍摄；(b)色彩处理后的结果

然而受到卫星、载荷、时相、季节等因素差异的影响，拍摄的影像中通常会存在色彩差异[2]。如果色彩差异过大，会对影像的质量产生负面影响，如图2所示的细节损失和伪影现象等。为了解决这些问题，研究者在这方面做了大量的研究。在这些研究中，一种代表性的算法为基于底图的匀色算法。在该算法中主要存在两个注意点：①从底图及目标片中提取所需要的信息。具体的，首先从底图中提取低频信息，即全局的色彩信息，然后从目标片提取高频信息，即局部的细节信息;②底图低频信息与目标片高频信息的融合。提取底图低频信息与目标片高频信息后，需要通过恰当的方法，将二者进行融合。

图2 负面效果举例

上述两点直接影响匀色算法的效果，如果方法不当，会导致细节丢失、伪影或者曝光等负面结果。

为此，本研究在底图匀色算法的基础上，做了一些改进。

第一，引入了L0范数平滑滤波[3]，用于提取参考底图和目标片的低频信息，并进一步计算目标片的高频信息。

第二，在融合算法上，使用IHS融合算法[4]，并进行改进，用以消除曝光、细节丢失等问题。

总之，本文的贡献有如下几个方面：

(1)引入L0范数平滑滤波作为低通滤波。与其他滤波相比，L0范数平滑滤波可以获得更好的保边效果。该效果可以帮助消除伪影以及帮助计算出更好的目标片高频信息。

(2)改进了IHS融合算法。传统的IHS算法是将一张单通道的全色影像与一张多通道的多光谱影像进行融合，得到最终结果。在本研究中，对该方法进行改进，实现多通道对多通道的融合，并更好保留了细节信息。

为了更好阐述本研究所涉及的内容，本文将按下述方式进行组织：在第一节中，对已有的匀色方法进行简要描述；在第二节中，介绍了本研究中所用的方法及其改进；在第三节中，从主观以及量化评价两个角度，与其他方法进行实验对比与分析；在第四节中，对整个研究进行总结，并提出相应建议；在文章末尾，列举了涉及的参考文献。

1 相关研究

为了做好影像匀色处理，研究者在此领域投入了大量的研究。总的来说，这些方法可以被分为两类：第一类是直方图匹配法；第二类是数学模型法[5]。

第一类方法典型算法包括Wallis滤波[6]以及Reinhard色彩迁移等。该类方法有两点优势：一是计算高效；二是提供参数化的质量控制。Wallis方法提出之际主要用于图像增强。它可以增强高分辨率影像的对比度，并且可以压制影像中的噪点。在文献[7]中，Wallis滤波被用于遥感影像的色彩平衡处理。由于其高效性与有效性，文献[8]中采用了Wallis滤波算法，并在陆地和水域的处理中获得了较好的效果。在文献[9]中，研究者发现在LAB色彩空间中，各分量存在更小的相关性，这对色彩迁移来说很有利。在文献[10]中，对比了三种不同的色彩处理算法，分别为直方图规定化、Wallis滤波以及Reinhard色彩迁移算法。实验结果显示，当影像间色彩差异较大时，Reinhard色彩迁移算法可以获得更好的效果。在文献[11]中，Reinhard色彩迁移算法在一个局部窗口中进行计算，并获得了更好的效果。该类方法常选定测区内的某一影像作为参考，适用于景物内容相似的影像处理，当覆盖范围广、场景内容差异大时，容易造成整体或者局部偏色。

第二类方法典型算法有基于重叠区的直方图规定化和伽马校正等。在数字影像处理中，直方图规定化是最为广泛应用的方法之一。通过调整影像灰度分布，可以实现不同影像间的色彩一致化处理。在文献[12]中，累计直方图方法被用于影像间的色彩映射，由此消除不同影像帧之间的光照以及饱和度差异。然而，在文献[13]中，作者指出直方图规定化会破坏影像自身内部的色彩对比，造成色彩失真或者偏色问题。为此，在文献[14]中，作者提出了使用伽马校正实现影像间的色彩处理，但伽马系数会直接影响色彩处理质量。文献[15]提出了一种自动伽马系数计算方法，通过统计影像的色彩信息，结合内置质量评估算法，可以自动计算出最优系数，用于最后的色彩处理中。该类方法需要考虑影像之间色彩传递顺序，计算相对复杂；传递过程极有可能带来误差累积，导致原理基准位置的影像出现色差，不适用于覆盖范围比较大的影像色彩一致性处理。

随着遥感影像处理范围的不断扩大，需要使用更多方式获取对应测区的影像数据，这导致影像间的色彩差异因素越发多样化，处理的难度也越发增大[16-18]。目前解决大测区匀色问题多采用基于参考底图的色彩处理方法。该方法所选取的底图本身拥有较好的色彩一致性，色彩差异巨大的待处理影像在底图中找到相应参考区域，然后进行色彩迁移，最终实现待处理影像间的色彩一致化处理。

文献[17-18]均使用了不同方法将底图色彩迁移到目标片中，实现了色彩的一致化操作。文献[17]使用高斯滤波提取底图低频信息，计算目标片高频信息，然后使用线性方法将高频信息与低频信息进行融合得到结果。文献[18]使用Wallis滤波，统计底图均值与方差信息，然后修改目标片像素值，使得均值与方差信息与底图一致，进而实现色彩的统一。但文献[17]所述方法虽然可以获得较好的匀色效果，但存在影像对比度较差、细节损失、有压平感的问题。这点在文献[18]中也有所提及。因此，为更好消除影像间的色彩差异，保持影像细节信息不丢失，本研究提出了一种基于L0范数平滑滤波的遥感影像匀色算法。

2 提出的方法

本文提出的一种基于L0范数平滑滤波的遥感影像匀色算法，首先选择待处理影像对应测区的参考底图，使用L0范数平滑滤波提取底图的低频色彩信息以及目标片的高频细节信息，然后使用改进的IHS融合算法，实现多通道对多通道数据的融合。最后，得到的结果不仅实现影像间色彩的统一，同时也消除了细节损失、曝光等负面效果。

2.1 L0范数平滑滤波

本研究采用L0范数平滑滤波[3]来获得底图的低频信息。该方法在提取影像低频信息的同时，可兼顾边缘保持效果，具有计算复杂度低、实用性强的特点。具体原理如下：

影像高频信息存在于地物边缘处，属于边缘特征。在L0范数平滑滤波中，通过限制非零梯度的数量来增强影像中对比度最高的边缘，同时在全局进行影像平滑。既能实现边缘区域梯度信息不丢失，又能达到保边效果。

假设输入影像为I，平滑后的影像为S，计算S中每个像素在x和y方向上的差分：

∇Sp=(∂xSp,∂ySp)T

(1)

式中，p为像素的位置，梯度L0的范数可以表达为：

C(∇S)=#{p||∇Sp||1≠0,p∈Ω}

(2)

式中，#{.}为计数符号; Ω 为位置的取值范围;C(∇S) 控制图像边缘的保真度。

考虑影像的梯度L0范数。从全局考虑显著边缘的保护问题，并加入图像相似性约束，从而获得下述能量方程：

(3)

式中，K为全局的梯度数量约束，K越大，保留的细节越多，反之越小。由于在实际操作中很难确定K的最优值，因此(3)式可以改写为：

(4)

式中，λ为平滑系数，图像越平滑，保留的高频信息越少，即保留的图像边缘较少。

通过控制平滑系数，构建上述能量方程并求解，可以获得L0范数平滑滤波的平滑结果。

得到的平滑结果为影像的低频信息，然后将影像与低频信息相减取绝对值，即可得到对应的高频信息。

2.2 改进的IHS融合算法

传统的图像信息融合方法主要是全色影像与多光谱影像的融合[4, 16, 19]。全色影像为单通道影像，具有较高的空间分辨率，细节信息好；多光谱影像为多通道影像，空间分辨率低，色彩信息丰富。该类方法本身是单通道影像与多通道影像之间的融合。

本研究中，由于底图为多通道数据，目标片也为多通道数据，需要将底图多通道低频信息与目标片多通道高频信息进行融合。因此，本文基于传统的IHS融合算法进行了改进和优化。

2.2.1 传统的IHS融合算法原理

将RGB色彩影像转为IHS色彩影像[20]，具体计算公式如下：

I=(R+G+B)/3

(5)

(6)

(7)

如果B>G, 则有H=2π-H, 对应的转换方式如下。

如果H满足 0≤H<2π/3, 则有:

(8)

如果H满足 2π/3≤H<4π/3, 则有:

(9)

如果H满足 4π/3≤H<2π, 则有:

(10)

通过上述公式，可以将RGB色彩转为IHS色彩。然而参考文献[20,21]的方法，将I通道替换为单通道的高频信息。

一种效果较好的替换方式为文献[22]所述方式：

(11)

式中，μ和σ为影像的均值和方差。

2.2.2 IHS融合算法改进

本研究中目标片计算得到的高频细节信息为多通道数据，若将其视为新的I通道，使用IHS融合算法融合到底图低频色彩信息中时，需要在IHS融合算法上做一些改进。

首先按波段顺序，从底图低频色彩信息中提取相应通道并生成新的多通道数据；然后在目标片高频信息中，按相同波段顺序，提取相应通道，与新生成的多通道数据进行IHS融合，再将融合结果转为RGB数据；最后从多个RGB数据中抽取不同波段数据，组合成最后的结果。

详细过程见下表。

算法一:色彩一致化处理流程输入:底图Iref和目标片IsrcStep-1:计算底图Iref和目标片Isrc的低频信息Lref和Lsrc。Lref=L0Filter(Iref)Lsrc=L0Filter(Isrc)Step-2:计算目标片Isrc的高频信息Hsrc 。Hsrc=abs(Isrc-Lsrc)Step-3:对于底图低频信息Lref以及目标片高频信息Hsrc中的每一个通道,进行以下操作。特别的,作为约束条件,底图低频信息与目标片高频信息的通道数应相同,均为N。定义变量n,n的取值从1到N,为整数。 Step-4.1:CnLref和CnHsrc分别为Lref和Hsrc第n个通道数据; Step-4.2:对于CnLref中的每一个元素,分别进行+2和-2的操作,得到CnLref+2和CnLref-2,并将计算结果限制在[0,255]之间。然后将CnLref+2、CnLref、CnLref-2分别作为R、G、B通道,合成新的影像Nref; Step-4.3:将Nref转换到IHS色彩空间中; Step-4.4:使用IHS算法,将Nref的I通道替换为CnHsrc; Step-4.5:将Nref转回至RGB色彩空间; Step-4.6:提取Nref的G通道数据,作为结果影像result第n个通道的数据。 重复Step4中的每个子步骤,直到Lref 和 Hsrc中的每个通道都被处理过。输出: 结果影像result

图3所示为各阶段输出结果的示例。

图3 本文方法输出展示

3 实验分析

3.1 实验环境

本次实验中，使用Landsat卫星数据作为参考底图。参考底图的分辨率为30 m。目标片为广东广州市的航拍数据。

测试电脑使用的CPU型号为：CPU 是Intel(R) Core(TM) i7-6800 3.40 GHz，内存32 GB。数据经过PixelGrid做过正射校正。通过实验对比，推荐在L0范数平滑滤波中，平滑系数为0.02。对比方法选用GeoDoging使用的Wallis方法[23]，该方法在遥感影像匀色处理中，最为被广泛使用。然后通过实验，从主观以及量化两个角度对两种色彩一致化算法进行比较。

3.2 实验结果评价与分析

3.2.1 主观评价

选用了4组典型场景采用不同方法进行处理前后对比。

本文方法与Wallis方法处理结果的细节对比如下。

如图4、图5所示，原始影像以及不同方法处理结果均叠加在底图上。两种方法均实现了影像间的色彩统一，但与本文方法相比，Wallis方法结果缺乏立体感，整体有压平感。

图4 4组案例不同方法处理结果

图5 不同方法处理细节对比结果

案例一拍摄于不同时间。因此光照存在差异。Wallis方法对于光照差异不能进行较好地处理，处理后的影像依旧存在差异。而本文方法可以较好地消除光照差异。

案例二受到Halo现象的污染。在影像边缘存在晕光。Wallis方法不能消除此种问题，而本文方法可以较好地消除。

案例三和四都存在一定的色彩和光照差异。本文方法仍然获得比Wallis方法好的结果，具体见图6中所示细节。

图6 本文方法细节丢失以及曝光现象得到压制(与Wallis方法相比)

3.2.2 量化评价

在实验中，使用均值、方差和均方差来衡量匀色效果。均值越接近，表示色彩差异越小。方差越大，表示影像越立体，压平感越小。均方差越小也表示色彩差异越小。

4组案例的评价指标计算结果见表1，通过对比，本文方法在均值、方差以及均方差中均获得了较好的结果。

表1 量化结果对照表

在案例1、案例3、案例4中，本文方法均值较为接近，均方差较小，表示色彩较为统一。同时，方差较大，图像没有压平感和目视效果一致。

在案例2中，本文优势不明显。

4 结论

本文提出了一种新的图像匀色算法。通过使用L0范数平滑滤波，得到影像的高低频信息，并使用改进的IHS算法，获得最终的融合后的结果。并通过实验，验证了本文方法的优越性。

然而在部分案例中，本文方法效果并不尽如人意：例如在案例二中，本研究所使用的方法效果优势并不明显。此外，在选取底图的过程中，没有考虑地物差异以及色温差异带来的影响。例如，当底图存在云层或者地物变化较大时，可能会导致匀色结果的异常，部分地物色彩在处理后与实际应有色彩不一致。或者目标片为夏季数据，底图为秋季或者冬季数据，导致匀色后结果影像质量在全局角度的下降。

因此，为了更好地实现影像匀色，后续应重点进行地物变化检测，防止变化地物对结果的干扰；其次是确立底图选取标准，实现色彩一致化效果的标准化处理，尽量减少个人审美因素对结果的干扰。