图像去雾方法可以分为两大类：基于图像处理的雾天图像增强方法[1]和基于物理模型的雾天图像复原方法[2]。图像增强方法不考虑图像降质的原因，能提高雾天图像对比度，突出图像细节，然而对图像信息会造成一定损失。图像复原方法针对雾天图像降质过程，建立退化物理模型，反演退化过程，进而获得无雾图像[3]。基于物理模型的图像去雾方法得到的去雾效果更自然，信息损失较少，该类方法成为当前图像去雾领域研究的热点。
基于物理模型的图像去雾方法首先估计大气光和透射率，然后计算得到无雾图像。当前在估计介质透射率上已有大量的研究成果，而这些方法在大气光向量通常选取图像上的“最亮点”或者利用简单的规则来计算[4]。然而，不准确的大气光强会导致去雾后图像色彩的失真[5]，因此如何准确估计雾天图像中的大气光越来越受到国内外研究人员的重视。其中，Narasimhan和Nayar[6]提出了一种用户交互的去雾方法，该方法需要用户选取一组颜色相同但景深不同的区域。Fattal[7]提出了一种用户辅助的方法，用户需要选取若干个不同反照率的图像块，将包含所选图像块内像素的RGB平面相交估计得到大气光强。此后，Fattal等[8]又提出了一种自动估计大气光强的方法，该方法继承了文献[7]所述的方法，通过设置选取条件自动选取不同反照率的图像块，并通过最小化透射率和表面反照率之间的相关性来估计大气光的幅度。Tan[9]假设图像中最亮的像素为饱和状态并依此估计大气光。在He等[10]提出的暗通道方法中，大气光从最亮区域的像素同时也是暗通道中较亮区域的像素估计得到，这种不准确的大气光估计方法直接导致该方法得到的介质透射率是不准确的。Tarel和Hautiere[11]根据灰色实际假设首先对图像进行白平衡，然后利用纯白大气光向量A=[1, 1, 1]进行去雾。
对于无人机拍摄的航空图像而言，由于飞行高度大多在几百米或上千米[12]，地物高度差几乎可以忽略不计，因此获取的图像景深可以看为常数，当前方法的诸多假设条件已不再成立。另外，无人机对地观测图像通常不包含天空，因此不能通过简单选取像素的方法估计大气光。本文根据无人机对地观测时的成像特点，提出了自动估计单幅有雾无人机图像的大气光估计方法，进而计算得到无雾图像。与现有的户外图像去雾方法相比，本文方法对于无人机图像去雾具有更好的效果。
1 大气光方向的估计 1.1 大气光方向估计原理 基于物理模型的图像去雾方法常采用以下模型表示雾天图像[3]：
 | (1) |
式中：I为雾天图像；J为无雾图像；A为大气光向量；t为介质透射率；x为像素点的位置。
对于无人机成像而言，大多数情况是垂直向下或近似垂直向下成像的[13]。在这种情况下，可以认为无人机图像中各个像素点的介质透射率相同，定义这一常数介质透射率为图像的全局透射率，记为t0。因此，无人机图像雾天模型可表示为
 | (2) |
可以通过设置规则[8]自动选取图像中若干个不同反照率[5]的图像块，满足如下方程：
 | (3) |
对比式(2) 和式(3) 可知，无雾图像J(x)被分解为J(x)=l′(x)Ri，l′(x)为反射光强度的标量，Ri为第i个图像块的反照率向量。图像块的反照率可以认为是该图像块主体颜色的RGB向量，同一物体通常具有相同的反照率，其不同光照部分(如阴影部分和光照部分)具有不同的反射光强度。
式(3) 将l′(x)和t0合并为l(x)，则J(x)=l(x)Ri。可知，通过适当选取可以得到如下无雾图像块：该图像块上的大部分像素点处在同一条直线上，这条直线的方向为该图像块反照率向量Ri的方向。进一步分析可知，满足式(3) 的图像块I(x)上的像素点处在反照率向量Ri和大气光向量A张成的RGB平面上。这一原理可以估计大气光的方向：从一幅图像中选取2个或以上满足式(3) 的图像块，大气光向量A的方向
以图 1的户外图像为例说明该方法。图 1中：图像各个方框部分为选取不同反照率的图像块。图 2为各图像块(选取4个)中的归一化像素值及其拟合成的直线，其坐标为归一化像素值。图 3为各图像块拟合成的直线，选取各条直线拟合交点(垂直视图)的方向作为大气光方向的估计值。
 |
| 图 1 户外图像示例 Fig. 1 Sample outdoor image |
| 图选项 |
 |
| 图 2 图像块中的像素直线 Fig. 2 Pixel lines in image patches |
| 图选项 |
 |
| 图 3 各条像素直线及其交点 Fig. 3 Pixel lines and their intersection |
| 图选项 |
1.2 自动选取图像块 将一幅尺寸为M×N的图像划分为M×N/m2个m×m的图像块，可通过设置约束条件来自动从这些图像块中选取满足式(3) 的图像块。假设λ1≥λ2≥λ3为图像块经过PCA分解[14]后3个最大的特征值，这些约束条件如下[8]：
1) 最大主成分向量应为非负，即图像块表面反照率向量应为非负。
2) 图像块中不应包含边缘，可利用边缘检测算子(本文采用Sobel算子)滤除包含边缘的图像块。
3) 单一的最大特征值，即要求最大特征值λ1大于某一阈值τ1(λ1>τ1)。
4) PCA矩阵为一阶，即要求最大特征值λ1与第二大特征值λ2的比值大于某一阈值τ2(λ1/λ2>τ2)。
5) 图像块所在直线不应经过原点，即要求直线距原点的最小欧氏距离distance大于某一阈值τ3(distance>τ3)。
通过上述5个约束条件，可以筛选出满足式(3) 的图像块：该图像块的颜色基本一致，且不包含边缘。
在得到若干个这样的图像块后，计算各个图像块主体颜色的RGB向量，进而得到各向量所在直线的交点，原点到交点的向量方向即为大气光方向
2 大气光幅度的估计 2.1 全局透射率的范围 对于一般户外图像而言，每个像素点x对应的透射率为[15]
 | (4) |
式中：β为大气散射系数；d(x)为像素点x的景深。
对于无人机图像来讲，大气散射系数β随雾的浓度变化[15]，对于单张含雾无人机图像而言，β很难估计准确。由于景深d(x)数值较大且图像全部像素近似相等，全局透射率t0通常为数值接近0的小数，也有文献[8]直接假定距离较远处物体的介质透射率为0.01。不失一般性，无人机图像的全局透射率可设置为t0=ε(0＜ε≤0.01)。
2.2 计算大气光幅度 将大气光向量A写为其方向
 | (5) |
由式(5) 可知，图像块I(x)所在直线在向量Ri和
利用上述结论可以进行大气光幅度的估计：根据1.2节的约束条件，选取P个图像块的像素直线，分别计算各条像素直线在
2.3 获得去雾图像 本文方法的基本流程如图 4所示。对于一幅含雾图像I(x)，首先采用1.2节介绍的方法自动选取图像块，然后基于1.1节的大气光方向估计原理得到大气光方向
 | (6) |
 |
| 图 4 本文方法的基本流程 Fig. 4 Basic flow of proposed method |
| 图选项 |
3 提高方法鲁棒性的措施 本文方法的鲁棒性指对不同飞行高度、尺寸和地物特点的无人机图像均可以获得较准确的大气光估计。对上述方法流程加入如下措施以提高其鲁棒性。
3.1 自动调整图像块尺寸 从输入图像中选取满足式(3) 的图像块是估计大气光方向
3.2 自动调整条件阈值 为使本文方法鲁棒性增强，采取自适应方法调整1.2节中的阈值τ1、τ2和τ3，具体步骤如下：
1) 根据1.2节约束条件3)~约束条件5) 对上面获得的100个图像块分别单独排序，每一种条件选取50个图像块，形成一个图像块备选集，并记录初始的τ1、τ2和τ3。
2) 从图像块备选集中提取同时满足条件阈值τ1、τ2和τ3的图像块，排序得到前10个图像块(即P=10)。若同时满足的图像块数目少于10，则将阈值τ1、τ2和τ3分别降低5%再进行判断。
3) 根据获得的10个图像块计算其像素所在的直线，得10条直线的交点(两两相交后交点取均值)，大气光方向
4 实验结果与分析 4.1 数据源 为了验证本文方法的去雾性能，选用某型无人机在1 000、3 000和5 000 m飞行高度不同时间(天气条件不同)拍摄的3幅含雾数码照片作为实验图像源，其尺寸分别为512×500、1 392×1 040和1 392×1 040，如图 5(a)~图 5(c)所示。
 |
| 图 5 实验图像源 Fig. 5 Experimental images |
| 图选项 |
4.2 大气光估计 为了验证全局透射率设置对大气光估计的影响，实验中分别设置3组全局透射率t0为0.01、0.005和0.001。
从表 1可以看出，对于范围在(0，0.01) 之间的全局透射率，随着其取值的变化，估计得到的大气光的方向相同，而幅度变化微乎其微。这是由于在计算大气光幅度║A║时，利用到║A║(1－t0)的均值，当t0=ε(0＜ε≤0.01) 时，1－t0≈1，因此║A║(1－t0)可近似视为║A║。对于高空拍摄的无人机图像，通常可设置t0为数值较小的固定值，实验中取t0=0.01。
表 1 大气光向量估计结果 Table 1 Estimation results of atmospheric light vector
| t0 | 图 5(a) | 图 5(b) | 图 5(c) | |||||
| 方向 | 幅度 | 方向 | 幅度 | 方向 | 幅度 | |||
| 0.01 | 319.33 | 412.51 | 390.89 | |||||
| 0.005 | (0.081 5, 0.558, 0.523 5) | 319.67 | (0.618 9, 0.604 4, 0.502 4) | 412.69 | (0.653 8, 0.589 7, 0.474 4) | 391.15 | ||
| 0.001 | 319.91 | 413.16 | 391.24 | |||||
表选项
4.3 去雾效果评价 为评估本文方法的有效性，利用主观评价和客观评价2种方法分别对去雾效果进行评估。作为对比，选取3种经典的去雾方法(Fattal等[8]方法、Tan[9]方法和He等[10]方法)与本文方法进行比较, 实验结果如图 6~图 8所示。
 |
| 图 6 图 5(a)采用不同方法去雾效果对比结果 Fig. 6 Comparative results of Fig. 5(a) by dehazing with different methods |
| 图选项 |
 |
| 图 7 图 5(b)采用不同方法去雾效果对比结果 Fig. 7 Comparative results of Fig. 5(b) by dehazing with different methods |
| 图选项 |
 |
| 图 8 图 5(c)采用不同方法去雾效果对比结果 Fig. 8 Comparative results of Fig. 5(c) by dehazing with different methods |
| 图选项 |
在主观评价方面，从图 6~图 8的去雾结果可以看出，文献[9]方法的去雾结果出现明显的“光晕”效果，图像细节也变得模糊；文献[8]方法虽然增强了图像的对比度，但图像部分区域出现了颜色的过饱和现象；文献[10]方法虽然获得了较好的去雾效果，但由于简单地选取了“最亮像素”为大气光，使得去雾后图像亮度偏暗，图像细节变得模糊；相比而言，本文方法既有效地增强了图像的对比度，也较好地保护了图像细节，获得优良的去雾效果。
在客观评价方面，本文主要针对去雾图像的对比度和颜色保真度进行客观评估。其中，对比度评价指标采用Hautiere等[16]提出的新增可见边比e和图像去雾前后可见边梯度比r，其表达式如下：
 | (7) |
式中：n0和nr分别为有雾图像和去雾图像中可见边的数目；gr为去雾图像的平均梯度；g0为含雾图像的平均梯度。通常，e和r较大，表面去雾后的图像边缘强度较大，数量较多。
在去雾图像颜色保真度评价方面，本文采用色调保真度作为评价指标进行度量，这里色调保真度采用Jobson等[17]提出的基于图像统计特征的度量参数H，其代表色调变化，定义为
 | (8) |
式中：abs(·)表示对所有元素取绝对值；mean(·)表示对所有元素取均值；Hin(x)和Hout(x)分别为去雾前后彩色图像在HSV空间中的色调分量。通常，H值越小，去雾后的色调保真度越高。
表 2给出了4种方法去雾后的量化评价结果。可以看到，与其他3种方法相比，本文方法获得了最低的e和r值，以及最高的H值，说明本文方法获得了较强的纹理细节保护能力和颜色保真能力；文献[10]方法在H值上与本文方法相似，说明其颜色保真效果较好，但e和r值较本文方法有较大差距；文献[9]方法的量化指标最差，e和r均较低，H值较高；文献[8]方法虽然在量化指标上较文献[9]方法有所提升，但由于出现颜色的过饱和现象，与文献[10]方法和本文方法相比，各量化指标均不够理想。
表 2 几种方法的量化结果 Table 2 Quantitative results of several methods
| 图像源 | 文献[9]方法 | 文献[8]方法 | 文献[10]方法 | 本文方法 | |||||||||||
| e | r | H | e | r | H | e | r | H | e | r | H | ||||
| 图 5(a) | 0.02 | 1.28 | 0.123 1 | 0.04 | 1.33 | 0.0978 | 0.05 | 1.34 | 0.145 | 0.06 | 1.48 | 0.113 1 | |||
| 图 5(b) | 0.03 | 1.45 | 0.081 2 | 0.03 | 1.55 | 0.086 | 0.06 | 1.73 | 0.052 | 0.08 | 2.12 | 0.043 3 | |||
| 图 5(c) | 0.06 | 1.86 | 0.084 2 | 0.07 | 1.56 | 0.077 | 0.06 | 2.04 | 0.034 2 | 0.08 | 2.72 | 0.025 6 | |||
表选项
5 结论 1) 本文提出的无人机图像去雾方法获得了较好的去雾性能，与现有方法相比，无论从主观视觉上还是客观评价指标上均有较大提升。
2) 本文方法在计算大气光向量方向和幅度的过程中，充分利用了图像本身的信息，同时还考虑到不同图像的差异性，采取了自动调整图像块尺寸和条件阈值的方法来增强方法的鲁棒性。
3) 本文方法在计算大气光幅度时的约束条件是整幅无人机图像景深相似，但这一条件并不适用于无人机进行大倾斜角度成像[18]的情况，因此，笔者下一步的工作是将无人机姿态和飞行高度作为输入参数纳入方法设计中，以增强本文方法的适用范围。
参考文献
| [1] | 刘海波, 杨杰, 吴正平, 等. 基于暗通道先验和Retinex理论的快速单幅图像去雾方法[J].自动化学报, 2015, 41(7): 1264–1273. LIU H B, YANG J, WU Z P, et al. A fast single image dehazing method based on dark channel prior and Retinex theory[J].Acta Automatica Sinica, 2015, 41(7): 1264–1273.(in Chinese) |
| [2] | 张登银, 鞠铭烨, 王雪梅. 一种基于暗通道先验的快速图像去雾算法[J].电子学报, 2015, 43(7): 1437–1443. ZHANG D Y, JU M Y, WANG X M. A fast image haze removal algorithm using dark channel prior[J].Acta Electronica Sinica, 2015, 43(7): 1437–1443.(in Chinese) |
| [3] | NARASIMHAN S G, NAYAR S K. Vision and the atmosphere[J].International Journal of Computer Vision, 2002, 48(3): 233–254.DOI:10.1023/A:1016328200723 |
| [4] | 周雨薇, 陈强, 孙权森, 等. 结合暗通道原理和双边滤波的遥感图像增强[J].中国图象图形学报, 2014, 19(2): 313–321. ZHOU Y W, CHEN Q, SUN Q S, et al. Remote sensing image enhancement based on dark channel prior and bilateral filtering[J].Journal of Image and Graphics, 2014, 19(2): 313–321.(in Chinese) |
| [5] | LAI Y H, CHEN Y L, CHIOU C J, et al. Single-image dehazing via optimal transmission map under scene priors[J].IEEE Transactions on Circuits and Systems for Video Technology, 2015, 25(1): 1–14.DOI:10.1109/TCSVT.2014.2329381 |
| [6] | NARASIMHAN S G, NAYAR S K. Contrast restoration of weather degraded images[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 2003, 25(6): 713–724.DOI:10.1109/TPAMI.2003.1201821 |
| [7] | FATTAL R. Single image dehazing[J].ACM Transactions on Graphics, 2008, 27(3): 1–9. |
| [8] | SULAMI M, GLATZER I, FATTAL R, et al.Automatic recovery of the atmospheric light in hazy images[C]//2014 IEEE International Conference on Computational Photography (ICCP).Piscataway, NJ:IEEE Press, 2014:1-11. |
| [9] | TAN R T. Visibility in bad weather from a single image[C]//IEEE Conference Computer Vision and Pattern Recognition.Piscataway, NJ:IEEE Press, 2008:1-8. |
| [10] | HE K, SUN J, TANG X. Single image haze removal using dark channel prior[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(12): 2341–2353.DOI:10.1109/TPAMI.2010.168 |
| [11] | TAREL J P, HAUTIERE N.Fast visibility restoration from a single color or gray level image[C]//IEEE 12th International Conference on Computer Vision.Piscataway, NJ:IEEE Press, 2009:2201-2208. |
| [12] | 姚鹏, 王宏伦, 刘畅. 基于动态流体扰动原理的三维滚动航路规划[J].北京航空航天大学学报, 2015, 41(12): 2280–2287. YAO P, WANG H L, LIU C. Three-dimensional rolling path planning via dynamic fluid disturbance[J].Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(12): 2280–2287.(in Chinese) |
| [13] | LIU S, LI H, YUAN Y, et al. A method for uav real-time image simulation based on hierarchical degradation model[J].Advances in Intelligent Systems and Computing, 2014, 277: 221–232.DOI:10.1007/978-3-642-54924-3 |
| [14] | MATEOS G, GIANNAKIS G B. Robust PCA as bilinear decomposition with outlier-sparsity regularization[J].IEEE Transactions on Signal Processing, 2012, 60(10): 5176–5190.DOI:10.1109/TSP.2012.2204986 |
| [15] | MCCARTNEY E J. Optics of the atmosphere:scattering by molecules and particles[M].New York: Wiley, 1976. |
| [16] | HAUTIERE N, TAREL J P, AUBERT D, et al. Blind contrast enhancement assessment by gradient ratioing at visible edges[J].Image Analysis and Stereology, 2008, 27(2): 87–95. |
| [17] | JOBSON D J, RAHMAN Z, WOODELL G A.Statistics of visual representation[C]//Proceedings of the 2002 Visual Information Processing XI.Bellingham:SPIE, 2002:25-35. |
| [18] | HUANG Y, DING W, LI H. Haze removal for UAV reconnaissance images using layered scattering model[J].Chinese Journal of Aeronautics, 2016, 29(2): 502–511.DOI:10.1016/j.cja.2016.01.012 |
