基于传播滤波的单幅图像去雾方法与流程

本发明涉及计算机图像处理领域，具体涉及一种单幅图像去雾技术，尤其是一种结合最小颜色通道图与传播滤波的单幅图像去雾的方法。

背景技术：

近几年空气质量退化严重，雾霾等恶劣天气出现频繁，pm2.5值越来越引起人们的广泛关注。在有雾天气下拍摄的图像模糊不清，清晰度不够，细节不明显，色彩保真度下降，出现严重的颜色偏移和失真，达不到满意的视觉效果。由此限制和影响了室外目标识别和跟踪、智能导航、公路视觉监视、卫星遥感监测、军事航空侦查等系统效用的发挥，给生产与生活等各方面都造成了极大的影响。

以公路监控为例，由于大雾弥漫，道路的能见度大大降低，司机通过视觉获得的路况信息往往不准确，进一步影响对环境的判读，很容易发生交通事故，此时高速封闭或公路限行，给人们的出行带来了极大的不便。故现在迫切需要对雾天图像进行有效的去雾处理。

图像去雾算法目前已取得较大的进展，主要可分为两大类：第一类是基于图像增强的图像去雾方法，如直方图均衡化算法、基于小波变换的图像增强算法和基于retinex理论的图像增强算法等，这类算法并不考虑图像退化模型及降质原因，仅从图像处理的角度提高图像对比度，因此容易造成图像的信息丢失，并不能从根本上去雾；第二类是基于物理模型的图像去雾方法，该类方法主要基于大气散射模型，研究雾天图像降质的物理过程，再反推恢复景物无雾图像，例如，fattal等假设大气透射率和场景目标局部区域的反射率两者不相关，从而估算大气透射率及场景的辐照度，取得了一定的复原效果，但该方法不适用于浓雾天气下的图像去雾。he等提出了基于暗原色先验的图像去雾方法，该方法首先根据暗原色先验信息估算大气透射率，再采用软抠图算法对大气透射率进行细化，从而复原图像，但该算法时间复杂度大。因此gibson等提出用中值滤波估计大气散射函数，从而实现图像恢复，但该方法容易导致边缘信息丢失，造成去雾后图像的黑斑效应。he等后续提出采用引导滤波代替软抠图实现透射率的优化，该图像复原的速度大大加快，但因采用原图作为引导图像，容易导致图像去雾残留较大。此外，zhang等则提出采用双区域滤波代替中值滤波来估算大气透射率，该方法能在一定程度上能改善黑斑效应，但通过实验可观察到在图像景深突变密集的边缘区域，该算法的透射率估值容易出现偏差，边缘区域去雾不彻底，如图9(c)所示。

综上，现有的图像去雾方法无法兼顾图像复原速度快和图像复原质量高的优点，具有一定的局限性。

技术实现要素：

本发明的目的是克服上述现有技术的缺点，提供了一种复原速度快、质量高的图像去雾方法，使用该方法的图像去雾彻底、无黑斑、无残留、清晰度高。

基于传播滤波的单幅图像去雾方法，首先通过双区域滤波法估计初始大气透射率，再以最小颜色通道图作为参考图像，结合传播滤波法得到优化后的大气透射率，最后采用自适应方法优化的大气光强度，基于大气散射模型实现无雾图像的恢复。

优选地，上述方法的步骤如下：

步骤a)对有雾图像进行双区域滤波处理，结合雾、霾环境下的大气散射模型，得到大气初始透射率t′(x)；

步骤b)通过最小颜色通道图运用传播滤波器将步骤a所得大气初始透射率t′(x)优化后得到大气透射率t″(x)；

步骤c)采用大气光强度的自适应恢复算法，实现对大气光a的估计，结合步骤b所得t″(x)得到无雾图。

更优选地，步骤a的具体操作是：根据暗原色先验理论得到暗原色先验图，对暗原色先验图像进行双区域滤波操作，得到暗原色双区域值i′dark(x)，根据雾、霾环境下的大气散射模型，结合i′dark(x)，得到大气初始透射率t′(x)。

进一步，建立雾、霾环境下的大气散射模型公式为：

i(x)＝j(x)t(x)+(1-t(x))a(1)

式中，x是像素点索引，i(x)表示有雾图，j(x)表示无雾图，t(x)为大气透射率，a为大气光强度。

通过对大量无雾图像统计发现，对于清晰的户外图像，在绝大多数非天空的局部区域内，某些像素总存在至少一个颜色通道亮度值很低并趋近于0，即暗原色先验图，计算公式如下：

式中，上标c表示r、g、b某个通道，ω(x)表示以像素x为中心的邻域，则j^c(y)为对应ω(x)中任一像素的r、g、b某个颜色通道值，则jdark(x)是两次最小值滤波结果。

假设ω(x)邻域内t(x)保持不变，则对式(1)两边进行两次最小值运算：

将式(2)代入式(3)，得大气透射率t(x)初始值：

为了使图像恢复结果更加逼真，引入因子θ，θ＝0.95，因此，式(4)中的大气初始透射率为：

随后采用双区域滤波优化透射率：以对无雾图j(x)进行透射率优化为例，对每个像素x的邻域内三原色通道最小值排序，取其中值，若中值大于像素点x三原色通道的最小值，暗原色双区域值则取三原色通道的最小值，否则取像素点邻域内三原色通道最小值的中值，具体定义如下：

式中，j^c(x)为j(x)中像素x的r、g、b某个颜色通道值，j′dark(x)为像素x的暗原色双区域值。

类似地，对有雾图i(x)进行双区域滤波处理，得到暗原色双区域滤波图i′dark(x)；将所得i′dark(x)代入式(5)，得到重新定义的大气透射率t′(x)：

更优选地，步骤b的具体操作是：由有雾图i(x)得到最小颜色通道图w(x)，运用传播滤波器对步骤a中所得大气透射率t′(x)进行优化，得到优化后的大气透射率t″(x)。

进一步，由有雾图i(x)得到最小颜色通道图w(x)的公式为：

式中，i^c(x)为i(x)中像素x的r、g、b某个颜色通道值，则w(x)为像素x中三个颜色通道的最小值。

计算传播滤波器输出v^p(x)：

式中，v(y)代表像素值，为归一化系数，且

wx，y＝g(d^a(v(x)，v(y))；σa)g(d^r(v(x)，v(y))；σr)(10)

式中，为高斯函数，σ为方差，s、s+1表示区域邻接路径φ上相邻的两个像素点，v(s)、v(s+1)、v(x)分别表示像素s、s+1、x对应的像素值，d^a(v(x)，v(y))定义了邻接路径φ上相邻像素间的关联程度，d^r(v(x)，v(y))定义了邻接路径上φ每个像素与中心像素x的关联程度，φ则表示两个像素x与y之间相连的区域邻接路径上所有像素点集合，而且由相邻的像素迭代生长而成，其生长规则如下：

(a)若两像素在水平或垂直方向，如ys-2与x，则邻接路径为两像素直接相连的线段；

(b)若两像素不在水平或垂直方向，如y与x，则通过两个像素间的曼哈顿距离决定路径φ，具体定义如下：

式中，(a1，b1)和(a2，b2)代表像素y与x的坐标。

则优化后的大气透射率t″(x)为：

更优选地，步骤c的具体操作是：采用大气光强度的自适应恢复算法，实现对大气光a的估计，根据步骤b所得t″(x)，结合所得a对有雾图i(x)进行大气散射模型反解，得到无雾图j(x)。

进一步，由于无雾图j(x)可分解为图像的明暗系数和表面反射系数两个变量的乘积，即

j(x)＝l(x)r(x)(15)

式中，l(x)为图像明暗系数，r(x)为图像表面的反射系数，其矢量值与rgb颜色空间有关，且||r(x)||＝1；

将有雾图像i(x)分块，基于式(1)和式(15)，第i个图块的去雾模型ii(x)可以重新修改如下：

ii(x)＝l(x)ri(x)ti(x)+(1-ti(x))a(16)

式中，ri(x)，ti(x)分别代表第i个图块中像素的表面反射系数和透射率，ii(x)由矢量l(x)ri(x)ti(x)和矢量(1-ti(x))a合成。根据图像具有局部区域平滑的特性，则第i个图块内的表面反射系数ri(x)和大气透射率ti(x)基本保持不变，因此ii(x)值仅随l(x)的变化而变化。

在rgb颜色空间中，像素值ii(x)的分布在同一条线上，即第i个图块的rgb线，并相交于矢量(1-ti(x))a；通过分析若干图像块，发现较多图块内像素点主要集中在同一条rgb线附近，且不经过原点。找到满足式(16)的图块，则每个图块构成的rgb线与原点将构成一个平面，然而每个平面都包含着矢量(1-ti(x))a，因此，将两两平面相交的线设定为侯选大气光强度的方向然后，计算第i个图块对应的rgb线到每条候选线的欧氏距离dij，则大气光强度a的方向定义如下：

由式(17)可知，大气光强度的方向即当所求的dij取中值时对应的第i个图块与所有侯选相交线距离最小的

为了准确的估算出大气光强度的模值，引入惩罚因子a，则式(17)可转化为：

i(x)＝la(x)ra(x)ta(x)+(1-ta(x))aa(18)

式中，la(x)、ra(x)与ta(x)分别代表当大气光强度为aa时，像素x的明暗系数、表面反射系数与透射率，由式(18)与式(16)可知，(1-ta(x))aa＝(1-t(x))a，得

同理：la(x)ra(x)ta(x)＝l(x)r(x)t(x)，其中ra(x)＝r(x)，

令因此，当a＞1时，则0＜σa(t)＜1，大气光强度估值偏大造成去雾后的图像偏暗；当a＜1时，则t＞1-a，σa(t)＞1，大气光强度估值偏小则使得去雾后的图像偏亮，这种趋势在低透射率区域(即远处)尤为明显。

由于图像的不同区域的透射率存在着一定的差别，则可将透射率划分为若干等级s，根据对自然图像的统计可知：

l^*(s)＝percentile{l(x)|t(x)＝s}(21)

式中，当透射率级别为s时，l^*(s)值则为去雾后图像的最大亮度的1％像素取平均值，而引入惩罚因子后的大气光强度模值aa，同样具有：

la^*(s)＝percentile{la(x)|ta(x)＝s}(22)

式中，代表采用大气光强度为aa时图像去雾后，其透射率为s的最亮的前1％像素点的平均值。

由式(19)、式(20)与式(22)可知：

将雾图模型分成三个颜色通道，模型也可写成以下形式：

i^c(x)＝j^c(x)t(x)+(1-t(x))a^c(24)

其中，以i^c(x)为纵轴，j^c(x)为横轴，t(x)为斜率，建立平面直角坐标系，若图像的像素值j^c＝a^c，i^c(x)＝a^c，对应于不同透射率的直线都相交于点(a^c，a^c)。

在不同等级的透射率条件下，处于j^c局部区域内像素值基本保持不变，假设矢量是一个方向模值为a||a||，随着透射率的变化，去雾后图像的最大明暗系数基本保持不变，则有因此，则式(23)可重新定义如下：

利用l-bfgs优化算法实现对式(25)的求解，得计算得到正确的大气光的模值||a||，然后根据正确大气光强度的模值||a||和方向得到矢量

根据已求得的t″(x)和对大气光的估计a，并结合雾、霾环境下的大气散射模型式(1)，得到复原后的无雾清晰图像j(x)：

与现有技术相比，本发明所用图像去雾方法首先通过双区域滤波法估计初始大气透射率，再以最小颜色通道图作为参考图像，结合传播滤波法得到优化后的大气透射率，最后采用l-bfgs自适应方法优化的大气光强度，基于大气散射模型实现无雾图像的恢复。其有益效果是：能得到更为精确的大气透射率，可以改善现有的图像去雾方法中景深突变边缘的去雾残留问题，同时边缘细节得到有效保持，局部纹理区域具有较好的空间平滑性。

附图说明

图1为本发明提供的基于传播滤波的单幅图像去雾方法的算法流程图；

图2为本发明实施例1提供的有雾图；

图3为图2的暗原色先验图；

图4为图2的双区域滤波图；

图5为图2的最小颜色通道图；

图6为实施例1的路径模型图；

图7为单个原色通道模型与透射率关系图；

图8为图2某图块模型的rgb线；

图9(a)为图2采用gibson算法得到的去雾图；

图9(b)为图2采用he算法得到的去雾图；

图9(c)为图2采用zhang算法得到的去雾图；

图9(d)为图2采用本发明实施例1得到的去雾图；

图10(a)为图9(a)的局部放大图；

图10(b)为图9(b)的局部放大图；

图10(c)为图9(c)的局部放大图；

图10(d)为图9(d)的局部放大图；

图11(a)为对比例1的透射率图；

图11(b)为对比例2的透射率图；

图11(c)为对比例3的透射率图；

图11(d)为实施例1的透射率图；

图12(a)为图11(a)的局部放大图；

图12(b)为图11(b)的局部放大图；

图12(c)为图11(c)的局部放大图；

图12(d)为图11(d)的局部放大图；

图13为表1中mountain图像去雾效果比较图；

图14为表1中canyon图像去雾效果比较图；

图15为表1中train图像去雾效果比较图；

图16为表1中road图像去雾效果比较图；

图17为表1中swans图像去雾效果比较图。

具体实施方式

为了能够更清楚地描述本发明的技术内容，下面结合具体实施例来进行进一步的描述。

图1为本发明提供的基于传播滤波的单幅图像去雾方法的算法流程图，以对图2进行去雾为例，采用本发明提供的图像去雾方法，步骤如下：

1.建立雾、霾环境下的大气散射模型i(x)：

i(x)＝j(x)t(x)+(1-t(x))a(1)

式中，x是像素点索引，i(x)表示有雾图，j(x)表示无雾图，t(x)为大气透射率，a为大气光强度。

2.暗原色先验j^c(y)：

通过对大量无雾图像统计发现，对于清晰的户外图像，在绝大多数非天空的局部区域内，某些像素总存在至少一个颜色通道亮度值很低并趋近于0，即暗原色先验图，公式如下：

上标c表示r、g、b某个通道，ω(x)表示以像素x为中心的邻域，则j^c(y)为对应ω(x)中任一像素的r、g、b某个颜色通道值，则ldark(x)是两次最小值滤波结果，图2对应的暗原色先验如图3所示。

3.大气初始透射率t′(x)：

假设ω(x)邻域内t(x)保持不变，则对式(1)两边进行两次最小值运算：

将式(2)代入式(3)，得大气透射率t(x)初始值：

为了使图像恢复结果更加逼真引入因子θ，因此，式(4)中的大气初始透射率t′(x)为：

随后采用双区域滤波优化透射率：以对无雾图j(x)进行透射率优化为例，对每个像素x的邻域内三原色通道最小值排序，取其中值，若中值大于像素点x三原色通道的最小值，暗原色双区域值则取三原色通道的最小值，否则取像素点邻域内三原色通道最小值的中值，具体定义如下：

式中，j^c(x)为j(x)对应ω(x)中像素x的r、g、b某个颜色通道值，j′dark(x)为像素x的暗原色双区域值；

本实施例中，对有雾图i(x)进行双区域滤波处理(类似式(6))，得到暗原色双区域滤波图i′dark(x)，如图4所示；将所得i′dark(x)代入式(4)，得到重新定义的大气透射率t′(x)：

式中，已知θ，a为大气光强度，i′dark(x)为本实施例所得暗原色双区域值。

4.由有雾图得到最小颜色通道图w(x)：

式中，i^c(x)为像素x对应的r，g，b值，则w(x)为像素x中三个颜色通道的最小值，所得如图5所示。

5.计算传播滤波器输出v^p(x)：

式中，v(y)代表像素值，为归一化系数，且

wx，y＝g(d^a(v(x)，v(y))；σa)g(d^r(v(x)，v(y))；σr)(10)

式中，为高斯函数，σ为方差，s、s+1表示区域邻接路径φ上相邻的两个像素点，v(s)、v(s+1)表示像素s、s+1对应的像素值，d^a(v(x)，v(y))定义了邻接路径φ上相邻像素间的关联程度，d^r(v(x)，v(y))定义了邻接路径上φ每个像素与中心像素x的关联程度，φ则表示两个像素x与y之间相连的区域邻接路径上所有像素点集合，而且由相邻的像素迭代生长而成，其生长规则如下：

(a)若两像素在水平或垂直方向，如ys-2与x，则邻接路径为两像素直接相连的线段，如图6中的ys-2→x；

(b)若两像素不在水平或垂直方向，如y与x，则通过两个像素间的曼哈顿距离决定路径φ，具体定义如下：

式中，(a1，b1)和(a2，b2)代表像素y与x的坐标，如图6所示，像素y与x对应的曼哈顿距离dy，x＝5，则邻接路径向垂直方向迭代生长，即y→ys；再计算像素ys与x对应的则邻接路径向水平方向迭代生长，即ys→ys-1。再依此类推因此邻接路径继续向垂直方向迭代生长，则ys→ys-1→ys-2，而ys-2与x保持在水平方向，则根据生长规则(a)组成完整的路径ys→ys-1→ys-2→x。

6.根据步骤4所得w(x)并采用步骤5所得v^p(x)的传播滤波器对双区域滤波处理后重新定义的大气透射率t′(x)进行优化，得到优化后的大气透射率t″(x)：

所得如图7所示。

7.大气光强度a的方向

由于无雾图j(x)可分解为图像的明暗系数和表面反射系数两个变量的乘积，即

j(x)＝l(x)r(x)(15)

式中，l(x)为图像明暗系数，r(x)为图像表面的反射系数，其矢量值与rgb颜色空间有关，且||r(x)||＝1；

将有雾图像i(x)分块，基于式(1)和式(15)，第i个图块的去雾模型ii(x)可以重新修改如下：

ii(x)＝l(x)ri(x)ti(x)+(1-ti(x))a(16)

式中，ri(x)，ti(x)分别代表第i个图块中像素的表面反射系数和透射率，ii(x)由矢量l(x)ri(x)ti(x)和矢量(1-ti(x))a合成。根据图像具有局部区域平滑的特性，则第i个图块内的表面反射系数ri(x)和大气透射率ti(x)基本保持不变，因此ii(x)值仅随l(x)的变化而变化。

图8为图2选取的第i个图块的rgb线，在rgb颜色空间中，像素值ii(x)的分布在同一条线上，即第i个图块的rgb线，并相交于矢量(1-ti(x))a；通过分析若干图像块，发现较多图块内像素点主要集中在同一条rgb线附近，且不经过原点。找到满足式(16)的图块，则每个图块构成的rgb线与原点将构成一个平面，然而每个平面都包含着矢量(1-ti(x))a，因此，将两两平面相交的线设定为侯选大气光强度的方向然后，计算第i个图块对应的rgb线到每条候选线的欧氏距离dij，则大气光强度a的方向定义如下：

由式(17)可知，大气光强度的方向即当所求的dij取中值时对应的第i个图块与所有侯选相交线距离最小的

8.估计大气光强度a的模值||a||

为了准确的估算出大气光强度的模值，引入惩罚因子a，则式(17)可转化为：

i(x)＝la(x)ra(x)ta(x)+(1-ta(x))aa(18)

式中，la(x)、ra(x)与ta(x)分别代表当大气光强度为aa时，像素x的明暗系数、表面反射系数与透射率，由式(18)与式(16)可知，(1-ta(x))aa＝(1-t(x))a，得

同理：la(x)ra(x)ta(x)＝l(x)r(x)t(x)，其中ra(x)＝r(x)，

令因此，当a＞1时，则0＜σa(t)＜1，大气光强度估值偏大造成去雾后的图像偏暗；当a＜1时，则t＞1-a，σa(t)＞1，大气光强度估值偏小则使得去雾后的图像偏亮，这种趋势在低透射率区域(即远处)尤为明显。

由于图像的不同区域的透射率存在着一定的差别，则可将透射率划分为若干等级s，根据对自然图像的统计可知：

l^*(s)＝percentile{l(x)|t(x)＝s}(21)

式中，当透射率级别为s时，l^*(s)值则为去雾后图像的最大亮度的1％像素取平均值，而引入惩罚因子后的大气光强度模值aa，同样具有：

la^*(s)＝percentile{la(x)|ta(x)＝s}(22)

式中，代表采用大气光强度为aa时图像去雾后，其透射率为s的最亮的前1％像素点的平均值。

由此可知，可将式(19)、式(20)与式(22)可知：

将雾图模型分成三个颜色通道，模型也可写成以下形式：

i^c(x)＝j^c(x)t(x)+(1-t(x))a^c(24)

其中，c为rgb任一颜色通道，以i^c(x)为纵轴，j^c(x)为横轴，t(x)为斜率，建立平面直角坐标系，若图像的像素值j^c＝a^c，i^c(x)＝a^c，对应于不同透射率的直线都相交于点(a^c，a^c)，如图8中q点。

由图8可知，q位置的j^c值不会随着透射率变化，且在不同等级的透射率条件下，处于j^c局部区域内像素值基本保持不变，若大气光强度的方向已知，假设矢量是一个方向为其模值为a||a||，随着透射率的变化，去雾后图像的最大明暗系数基本保持不变，则有因此，本文结合图像明暗系数的最大值与透射率的等级无关的这一约束条件，提出最小化图像明暗系数与透射率估计值的距离为目标函数，则式(23)可重新定义如下：

利用l-bfgs优化算法实现对式(25)的求解，得

再计算得到正确的大气光的模值||a||，然后根据正确大气光强度的模值||a||和方向即为矢量a。

综上，对图进行分块，根据每个图块的像素在rgb颜色空间分布在同一条线上，可求得大气光强度的方向然后对大气光强度模值||a||引入一个惩罚因子α实现图像去雾，对去雾后图像提出优化大气光强度模值的目标函数，根据图像明暗系数的最大值与透射率的等级无关的约束条件，利用l-bfgs优化从而得到正确的大气光强度模值||a||，最后根据大气光强度的方向和模值||a||得到矢量

9.根据步骤6已求得的t″(x)、步骤8对大气光的估计a，并结合雾、霾环境下的大气散射模型(式(1))，得到复原后的无雾清晰图像j(x)：

无雾清晰图像j(x)如图9(d)所示，局部放大图如图10(d)所示。

对比例1gibson算法去雾

使用gibson算法对图2进行去雾，得到去雾图，如图9(a)所示，局部放大图如图10(d)所示。

对比例2he算法去雾

使用he算法对图2进行去雾，得到去雾图，如图9(b)所示，局部放大图如图10(d)所示。

对比例3zhang算法去雾

使用zhang算法对图2进行去雾，得到去雾图，如图9(c)所示，局部放大图如图10(c)所示。

去雾效果评价

1.现对实施例和对比例各算法求解的透射图进行评价，图11为各算法透射率图，图12为图11的局部放大图，由图11和图12可知，实施例1有效改善了景深突变边缘区域的透射率被估高的情况(图中叶子边缘)，更好的保持了图像边缘细节信息，具有更好的清晰度，且透射图在局部均匀的区域更加平滑(墙体区域)。

2.现分别对mountain、canyon、train、road和swans图采用本发明提供的算法与gibson算法、he算法和zhang算法的去雾效果进行评价，主要评价指标为：iv、dsim、avegrad，其中：

iv表示细节强度，能有效表征图像边缘对比度，其公式定义如下：

式中，ls表示照度图像的总细节强度，n表示像素点邻域内的像素个数，利用canny算子检测图像边缘并求和；而去雾图像的halo强度ihalo表示图像亮通道的估计值。

dsim表示原图像和去雾后图像直方图的形状相似度，相似度越高则表征去雾图像的色彩还原能力较好，可用两幅图像相关系数来测量，其定义如下：

式中，h，h′为去雾图像前后的灰度值，为h，h′的均值，dsim(h，h′)越大，说明相似度越高，色调还原程度越好。

avegrad表示图像的平均梯度，要反映了图像细节的对比度，其定义如下：

其中，m，n为图像的尺寸，平均梯度值越高，则图像的清晰度越高。所得结果如表1所示。

由表1可知，采用本申请算法对图像去雾所得细节强度值iv更高，即图像的细节表达能力更强；采用本申请算法对图像去雾所得图像相关系数dsim更高，即图像相似度更高，色调还原程度更好；采用本申请算法对图像去雾所得平均梯度avegrad更高，即图像的清晰度更高。

在此说明书中，本发明已参照其特定的实施例作了描述。但是，很显然仍可以作出各种修改和变换而不背离本发明的精神和范围。因此，说明书和附图应被认为是说明性的而非限制性的。

表1