基于均衡化与同态滤波的光学图像去云模型

费子豪 刘勇 李虎 吕鑫

摘   要：近年来对太阳耀斑的研究受到大多数研究人员的青睐，但获得的影像资料容易受到云层影响。在图像处理中，云雾的存在不仅使太阳耀斑信息被削弱，影响准确性，而且对太阳耀斑的分析、处理、识别和预测造成严重影响。因此，太阳耀斑图像中云雾的去除显得十分必要，对后期的研判、分析及预测都有着重要的意义。首先对原始图像进行像素值提取，并进行高斯降噪、椒盐降噪处理和小波分析，分解高低频信号。在此基础上，建立了基于高频分量直方图均衡化，低频分量同态滤波处理的小波变换模型来进行去云处理。同时，将新模型与未经分解的图像进行同态滤波、小波去云进行对比去云效果，并对其进行分析，组合去云处理的去云效果较明显。针对该算法，通过改变反射参数、边缘锐度等4个指标进行鲁棒性检验。

关键词：OpeCV处理  同态滤波  小波重构  鲁棒性分析

Abstract： In recent years， research on solar flares has been favored by most researchers， but the images obtained are susceptible to cloud imagery. In image processing， the presence of clouds not only weakens the solar flare information， but also has image accuracy， and it has a serious impact on the analysis， processing， recognition and prediction of solar flares. Therefore， the removal of clouds in the solar flare image is very necessary， and it has important significance for the later research， analysis and prediction. Firstly， the original image is extracted with pixel values， and Gaussian noise reduction， salt and pepper noise reduction processing and wavelet analysis are performed to decompose the high and low frequency signals. On this basis， a wavelet transform model based on high-frequency component histogram equalization and low-frequency component homomorphic filtering is established to perform de-cloud processing. At the same time， the new model and the undecomposed image are homomorphic filtered， and the wavelet is compared with the cloud to compare the cloud effect， and the cloud effect is combined. For the algorithm， the robustness test is carried out by changing the four parameters of reflection parameters and edge sharpness.

太阳耀斑（Solar flare）是指发生在太阳表面局部区域中突然和大规模的能量释放过程[1]。它是空间环境的主要扰动源， 对地球空间环境造成很大影响。平均而言，一个正常发展的黑子群（Sunspot groups）几乎几小时就会产生一个耀斑。

近年来对太阳耀斑的研究也受到大多数研究人员的青睐，但是在观测时受到云的遮挡，拍到的耀斑图片需要经过一定的处理才能被用作研究。

去云处理（Go cloud processing），一般是通过对同一地区不同时间的多幅光学图像进行融合[2-3]，消除云覆盖，而多幅图像的取得并不容易。而且，成像过程经常受到云的干扰，获得的光学图像对比度下降，清晰度和分辨率无法达到使用要求，图像数据可利用率大大降低。因此，本研究对太阳耀斑图像去云处理对于提高数据的有效利用率具有重要意义。

理想的去云算法是仅对有云的区域进行处理，而对无云区域保持灰度值不变，在达到去云效果的同时，保持无云区域的低频成分[4]。遥感影像中薄云的分布是渐变的，并无清晰边界，去云算法只能面向整幅图像，而无法针对特定的区域，在滤去图像低频成分的同时，也会损失一部分图像的高频部分，从而引起了无云区的改变。

现阶段已有不少学者对去云方法进行研究，通過资料查询，本文将典型的各类方法做了优化组合，得出最适合太阳耀斑的去云方法。

通过对实测例子的分析，将实测一年的太阳耀斑图片进行去云处理，但是观察所给出的数据，发现其中含有较多噪声，因此考虑进行降噪预处理。

一幅图像可被划分为高频部分与低频部分，云区信息位于含云遥感图像中的低频部分。用小波变换的方法去除云雾便是建立在此依据之上：对含有云区的太阳耀斑图像进行小波分解，本文将高通滤波与直方图均衡化图像增强的方法引入小波变换去除云区方法之中，构建了一种新的算法。再基于高频分量直方图均衡化处理与低频分量同态滤波处理的小波变换去云处理。

初步分析步骤如下。

（1）选取db4小波，采取Mallat算法对原始图像进行小波分解，得到一组高层细节系数、低层细节系数和近似系数。

（2）将处理后高、低频分量及近似系数小波重构，得到去云图像。

（3）若分界层n不当，增大或减小n值，重复（2）。

（4）确定分界层n，对低频细节系数重构得到图像高频分量，对高频细节系数重构得到图像低频分量。

（5）高频分量直方图均衡化处理，低频分量同态滤波处理。

（6）将处理后高、低频分量及近似系数小波重构，得到去云图像。

1  建模准备与定性、定量分析

1.1 OpeCV图形处理

均值、标准差和平均梯度是验证图像质量的常用指标。我们利用OpenCV提供的函数辅助计算图像指标。其中：（1）均值反映了图像的亮度，均值越大说明图像亮度越大，反之越小。（2）标准差反映了图像像素值与均值的离散程度，标准差越大说明图像的质量越好。（3）平均梯度反映了图像的清晰度和纹理变化，平均梯度越大说明图像越清晰。（4）信息嫡是反映了图像的信息量丰富程度，图熵值越大其信息量也就越大。

（1）计算图像的均值和标准差。

式中，M×N为图像的大小，P（i，j）为第i行、第j列的像素值，μ表示均值。使用meanStdDev函数计算均值和标准差。

（2）计算图像的平均梯度。

平均梯度的计算公式如下：

式中，M×N为图像的大小，为水平方向的梯度;为垂直方向的梯度。

（3）计算图像的信息熵。

信息嫡是用来衡量融合影像信息量丰富程度的重要指标，在影像的表示上为偏离影像直方图高峰值灰度区的大小，通过对图像信息熵的比较可以对比出图像的细节表现能力，融合影像的熵值越大其信息量也就越大。其中，对于灰度范围的图像直方图，其熵的定义式如公式：

式中，pi为第i个灰度的出现概率。计算过程由Matlab编程实现

1.2 测试图像

选取附件中的5张图片为例进行函数计算测试，测试图像为如图1所示。

测试结果如表1所示。

通过计算图像的指标值来定量分析去云处理后图像与原始图像的变化，其结果反映了图像的信息含量在太阳耀斑图像去云前后的变换情况。因此，本文主要采用影像的均值、标准差以及平均梯度法指标对处理结果进行定量评。

1.3 原始图像像素值处理

图2中，X轴表示为0～2048的像素值，Y轴表示为0～2048的像素值，而Z轴则表示为在Z轴上的灰度值。我们可以看出;在数据主要集中在中间部分;而X轴的像素值在向2048的正方向上扩展;Y轴的像素值在向2048的正方向上扩展;Z轴则在0正方向上扩展，扩展幅度较X、Y轴要大。

2  模型的建立

问题二要求将附件中一年的太阳耀斑图片进行去云处理，但是观察所给出的数据，发现其中含有较多噪声，因此考虑进行降噪预处理。再基于高频分量直方图均衡化处理与低频分量同态滤波处理的小波变换去云处理。

2.1 模型预处理

2.1.1 降噪算法

对于问题二，在处理图像的过程中，我们发现不论何种重建算法，不可避免会带入伪影点，同时由于数据含有噪声，成像质量并不令人满意。本文使用NLM（非局部均值）降噪方式对图像进行处理[5]。NLM算法是图像降噪领域非常有效的算法之一，效率较高，实现简单。其思路是对像素的某邻域窗口内的像素灰度值做加权平均，且像素越相似，权重越大。J.Huang等人将其运用到CT图像迭代重建过程中，显示其具有良好的性能[6]。

具体算法公式如下：

2.1.2 去噪预处理

我们取附件任意一个太阳耀斑的原始图像进行去噪预处理，并分别使用高斯算法，椒盐算法进行重建。结果如图3所示。

从图3中可以看出利用算法对于重建图像的去噪有着较好的效果，并且保留了较多的边缘信息。

2.2 问题二模型的建立

2.2.1 高頻分量直方图处理

结合建模准备和问题一提出的均匀化处理，针对问题二我们还是选取直方图匹配来描述图形的变化。

设Pi（i）和Pj（j）分别表示源太阳耀斑图像与去云处理后的图像的灰度分布概率密度函数，直方图匹配的目的就是通过调整太阳耀斑图像的直方图，使之具有Pj（j）所表示的形状。下面我们建立Pi（i）和Pj（j）之间的关系，首先对源太阳耀斑图像直方图进行均衡化处理，即求变换函数：

上式的逆变换为j=G-1（n）在对源太阳耀斑图像和希望云处理后的图像都做了均衡化处理后，两者便具有相同的概率密度函数，所以可以用源图像均衡化处理后的灰度级m来代替希望图像中的n，即有：

这就是说希望图像的灰度级j可由源图像均衡化后图像的灰度值来求算。还是选取问题一中的用到的20110106_075847.fts，即2011年1月获取的太阳耀斑影像数据给出的实例，对其进行均衡化后的直方图处理如图4所示。

2.2.2 低频分量同态滤波法处理

同态滤波是一种传统的去云方法 ，它是把频率过滤和灰度变换结合起来的之中处理方法[7]。图像模型变化为：

同态滤波法的依据是：图像分频后，云雾信息主要集中于低频部分这一特性。

主要原理是：将图像通过Fourier变换变换到频率域，选用合适的滤波器提取出薄云相对集中低频部分，然后设置滤波器的截止频率将其滤除，对滤波后的结果进行必要的背景增强，再通过Fourier逆变换把图像还原到空间域，最终得到去除薄云信息的图像。

其方法步骤流程图如图5所示。

（1）取对数。

对于两个函数的乘积来说，其Fourier变换是不可分开计算的，所以照度和反射的频率部分是无法直接分开进行操作的[8]。我们可对原始图像f（x，y）取对数，以分开照度分量和反射分量，得到：

（2）Fourier变换。

通过Fourier变换将它们转换到频域。

如果用频率分量表示，上式即可轉换成低频成分与高频成分的加性运算。上式可简单记为：

（3）选择合适的同态滤波器。

选择合适的同态滤波器，压缩低频段，拉伸高频段，从而使占据低频成分的云的信息从影像信息中剔除出去或抑制。薄云的去除效果取决于高通滤波时滤波函数的选择，为了减少乃至消除振铃效应，滤波器频率响应应该具有光滑的、缓慢变化的特征，应用中我们选用巴特沃思高通滤波器[9]。

2.2.3 小波变换去云处理

一幅图像含有低频与高频两个部分，并可被划分，对图像做小波分解，经一次分解后得到LL1、LH1、HL1、HH1四个部分。L代表低通分量，H代表高通分量，每部分的前后分量分别对应着竖直方向与水平方向。其中LL1为近似系数，LH1、HL1、HH1与各自小波系数相对应。经再一次小波分解后，LL1又得到了LL2、LH2、HL2、HH2      4个部分，依次类推，每一次小波分解后均可得到3组小波系数，过程如图8所示。

经过3次分解后得到的最高层近似系数LL3频率最低，其中集中了图像90%以上的能量[15]。图9为原始信号经小波变换以后，分解成不同层次系数的波形图。

由图9分析可知：图像经小波变换以后，可以分解成不同层次的高层细节系数、低层细节系数和近似系数，其中高层细节系数频率较低，低层细节系数频率较高。云雾图像中的云信息主要位于图像低频部分，对应小波分解后的高层细节系数;信息集中于图像高频部分，对应小波分解后的高层细节系数。可见，小波分解可有效分离云雾信息与地物信息。减小高层细节系数，增大低层细节系数后再做小波重构便可实现去云。

2.2.4 模型的确定

综上，可得基于高频分量直方图均衡化处理与低频分量同态滤波处理的小波变换去云处理的模型可总结为以下几点。

（1）选取db4小波，采取Mallat算法对原始图像进行小波分解，得到一组高层细节系数、低层细节系数和近似系数。

（2）将处理后高、低频分量及近似系数小波重构，得到去云图像。

（3）若分界层n不当，增大或减小n值，重复（2）。

（4）确定分界层n，对低频细节系数重构得到图像高频分量，对高频细节系数重构得到图像低频分量。

（5）高频分量直方图均衡化处理，低频分量同态滤波处理。

（6）将处理后高、低频分量及近似系数小波重构，得到去云图像。

3  结论

采用数据为20110106_075847.fts，即2011年1月获取的太阳耀斑影像数据，研究区域太阳耀斑图像灰度均值为749.666，标准差为49.288，平均梯度为9.457。图像中部出现明显的太阳耀斑，云层遮挡也十分明显。

根据问题二模型，进行了高频分量直方图均衡化处理与低频分量同态滤波处理的小波变换去云处理后的太阳耀斑对比图像结果如图10所示。

根据前面OpeCV图形所采用的实验数据为20110106_075847.fts，即2011年1月获取的太阳耀斑影像数据处理的定性分析参数计算结果如表2所示。

通过表2得出：

（1）去云处理后的标准差较大，说明太阳耀斑影像细节表现损失较小，处理后的图像清晰度较高。

（2）3类方法处理后的平均梯度较小，说明使图像变模糊了，而基于组合模型处理后的结果稍大，说明图像的清晰度改进，图像纹理变得更加清晰。

（3）前面3类方法处理后的图像信息熵较原始影像都有所减少，说明图像信息量都有所减少;而基于组合算法处理得到的信息熵增大，说明信息量增加。

（4）将新算法与未经分解的图像进行同态滤波、小波去云进行对比去云效果，并对其进行分析，后者的去云效果较明显。

从目视的角度可以很直观、很明晰地看到本文所使用的组合模型算法对太阳耀斑去除云层处理是较为有效的。

4  模型鲁棒性分析

在问题二中，我们假设了一些参数，如采用的巴沃特斯滤波器，压缩低频段，拉伸高频段，从而使占据低频成分的云的信息从影像信息中剔除出去或抑制，为了减少乃至消除振铃效应，滤波器频率响应应该具有光滑的、缓慢变化的特征，我们选取滤波器频率为0.414;同态滤波模型中时，在考虑照射分量时，我们采用的是RL=2.0;在考虑反射分量时，我们采用的是RH=1.0;在考虑图片边缘锐度时，我们采用的是c=2;在考虑同态滤波频率时，我们采用的是d0=10。

对于这些主观选择的因素时，我们对其进行鲁棒性分析来观察当这些因素与我们假定因素不同时对模型结果的影响来观察模型是否鲁棒，计算出4个主观假定值与结果（图像的信息熵值）如表3所示。

为了观察效果，我们直接做出4个主观假定值与结果（图像的信息熵值）的函数图像，如图11、图12所示。

图11左图为照射分量假定参数选定对去云后图像的信息熵值的影响，我们可以发现从参数选定的前3个，图像的信息熵值均在上升，但上升的幅度在不断减小，达到最大值后呈线性大幅下降，即我们取照射分量参数为0.7合理。

图11右图为反射参数选择对去云后的图像信息熵值的影响，我们可以发现反射参数选定的前4个的信息熵值在波动，且波动范围不大，在反射分量为1.2时达到最大值，达到最大值后的信息熵值不断减小，即反射参数取最大值为1.2图像效果最好。

图12左图为锐度边缘假定参数选定对去云后的图像信息熵值的影响，可以直观地看出曲线在1.5处到达峰值，图像的信息熵值稳定保持在3.57左右，而后出现波动，整体呈下降趋势，即可判断锐度边缘假定参数 确定为1.5是最优。

图12右图为滤波频率假定参数选定对去云后的图像信息熵值的影响，可以直观地看出曲线波动较大，总体呈上升趋势，即滤波频率假定参数在一定范围内与信息熵值呈正相关。

5  结语

本文利用直方图对g（x，y）进行直方图规定划，使处理后的图像在视觉上达到更好的效果;选用巴特沃斯滤波器，发挥其图像模糊程度轻，且滤除噪声效果好的优势;模型兼顾了太阳耀斑图像的空域和频域特点，能综合考虑局部对比度增强性能和频域信息的高通处理。

本文提出了去降噪算法和噪预处理的简单分析，并将高通滤波器与直方图均衡化图像增强的方法引入小波变换去除云区方法之中，构建了一种新的算法。

同时，在新模型的建立过程中还存在亟需解决的问题。傅利叶变换及其反变换中的自然对数底e为无限不循环小数，计算机算法实现中的截断误差无法避免，从而导致图像失真;损失图像无云区域的低频部分，破坏了图像的边缘信息，对非云区的边界造成影响。

本文所提出的方法在去云过程基本是交互式实现，太阳耀斑图像上云区的选取、直方图匹配参考区域及参数选择等工作都需要人工参与，这势必造成速度慢、精度低，影响处理效果等问题，应对如何实现智能化进一步地研究探讨。

參考文献

[1] 林元章.太阳物理学导论[M].北京：科学出版社，2000.

[2] Zhao P， Zhou XJ， Liu G.Decadal-centennial-scale change in Asian-Pacific summer thermal contrast and solar activity[J]. Chin Sci Bull，2011（56）：3012-3018.

[3] Wu JY， Shao XH， Kong XG， et al. Imprint of solar activity on Nanjing stalagmite annual layer thickness sequence during the Last Glacial Maximum[J]. Chin Sci Bull，2006（51）：441-447.

[4] 赖格英，辜晓青.基于ERDAS IMAGINE的遥感图像去云方法[A].GIS协会第六界年会论文集[C].成都：GIS协会，2001.

[5] 谈建国，周红妹，陆贤.NOAA卫星云检测和云修复业务应用系统的研制和建立[J].遥感技术与应用，2000，15（4）：228-231.

[6] 张新明，沈兰荪.基于小波的同态滤波器用于图像对比度增强[J].电子学报，2001，29（4）：531-533.

[7] 冯益明，雷相东，陆元昌.应用空间统计学理论解译遥感影像信息“缺失”区[J].遥感学报，2004，8（4）：317-322.

[8] 冯春，马建文，戴芹，等.一种改进的遥感图像薄云快速去除方法[J].国土资源遥感，2004（62）：1-5.

[9] 宋小宁，赵英时.MODIS图像的云检测及分析[J].中国图像图形学报：A辑，2003（9）：1079-1083.