本发明属于图像处理技术领域，更进一步涉及遥感图像处理技术领域中的一种基于相位一致性的遥感图像质量评价方法。本发明可应用于对遥感图像的质量进行客观评价，用于保证获取遥感图像数据的有效性和准确性。
背景技术：
遥感图像在成像、传输和存储的过程中，受到光学系统、大气环境以及其它因素的影响，最终会使得地面获得的图像存在条带噪声、一般噪声、模糊等各种图像质量下降的问题，遥感图像质量是比较各种遥感图像处理算法性能优劣以及优化遥感光学系统参数的重要指标，因此在遥感图像处理领域建立有效的图像质量评价机制具有重要的意义。
目前，图像质量评价一般包括主观评价和客观评价。主观评价是以人眼视觉效果作为质量评价准则，需要多名专业观测者参加评分，没有统一标准，需要多人参与，繁琐、费时费力、不能实时处理，不利于海量数据的自动判读。图像质量的客观评价方法根据需要参考图像信息量的多少可以分为全参考评价、部分参考评价和无参考评价三类。常用客观评价方法主要有峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)、结构相似度(Structure Similarity,SSIM)、信息熵、平均梯度等方法。
张飞艳、谢伟等人在论文“No-reference remote sensing image quality assessment based on natural scene statistical in wavelet domain”(Journal of Electronics and Information Technology,2011,33(11):2742-2747)中用遥感图像小波域统计特性，通过噪声、模糊对影像线性关系的影响方式及破坏程度的定量分析，来判定影像的噪声和模糊强度，最后，利用二者的加权和作为遥感影像的最终质量评价指标。该方法是一种无参考的图像质量评价方法，该方法仍然存在的不足是，综合评价模型比较简单，没有考虑人眼视觉的特点，评价结果不能更好的符合主观评价的结果。
北京空间机电研究所申请的专利“一种基于遥感图像相位相关性的图像质量判别方法”(公开号：CN102013100A，公开日：2011年4月13日，申请日：2010年11月26日)中公开了一种基于遥感图像相位相关性的图像质量判别方法。该方法先对遥感图像进行预处理得到遥感图像的周期分量频谱，然后向遥感图像的周期分量频谱中加入随机相位形成的N幅加入随机相位形成多幅遥感图像,再进行离散傅里叶逆变换得到N幅时域遥感图像，再计算出N幅时域遥感图像peru1,peru2,…,peruN的全变差值和相位相关值，最后利用图像相位相关值对多幅遥感图像进行质量评判。该方法是一种无参考的数字图像质量评价方法，优点是不需要任何参考图像数据，更适合于评价遥感图像。但是，该方法仍然存在的不足是，但是，该方法仍然存在的不足是，只考虑了图像的锐度，仅以单个参数表征图像质量优劣，没有综合考虑遥感图像中各个降质参数之间的相互关系，评价结果不能更好的符合主观评价的结果。
技术实现要素：
本发明的目的是针对上述现有技术的不足，提出一种基于相位一致性的遥感图像质量评价方法。本发明综合考虑遥感图像中条带噪声、一般噪声和模糊等因素的相互关系，利用相位一致性能够检测图像中边缘、阶跃和马赫带等特征点的特性，通过相位一致性的计算方式检测出遥感图像中的条带噪声、一般噪声和模糊，利用相位一致性表征对应的强度指标并得到评价参数，最后利用Choquet非线性函数将三个评价参数结合，得到一个能够对遥感图像进行全面评价的综合评价值，评价结果能更好地符合人的主观评价。
本发明实现的具体步骤包括如下：
(1)输入一幅待评价的遥感图像；
(2)获得相位一致性值：
采用相位一致性公式，计算遥感图像在条带噪声方向的相位一致性，得到相
位一致性值；
(3)二值化处理：
采用二值化方法，对相位一致性值进行二值化处理，得到0-1模板；
(4)获得条带噪声强度矩阵：
用遥感图像矩阵乘以0-1模板，得到条带噪声，将条带噪声的像素值进行归一化处理，得到条带噪声强度矩阵；
(5)按照下式，计算遥感图像的条带噪声值：
其中，N1表示遥感图像的条带噪声值；∑表示求和操作；j表示遥感图像像素的列坐标；i表示遥感图像像素的行坐标；S表示条带噪声强度矩阵；P表示遥感图像的相位一致性值；M表示遥感图像的行数；N表示遥感图像的列数；
(6)获得一般噪声信息：
(6a)采用边缘提取方法，提取遥感图像的边缘，得到边缘图；
(6b)采用形态学中的膨胀运算，对边缘图进行膨胀，得到膨胀后的边缘图；
(6c)采用局部相位算法，计算遥感图像的局部相位值，得到局部相位图；
(6d)从膨胀后的边缘图中找出像素值为1的坐标，将找出的坐标所对应的局部相位图中像素值赋值为0，得到去掉边缘信息的局部相位图；
(7)按照下式，计算遥感图像的一般噪声值：
其中，N2表示遥感图像一般噪声的值；∑表示求和操作；p表示去掉边缘信息的局部相位图的行标号；M表示遥感图像的行数；l表示去掉边缘信息的局部相位图的列标号；N表示遥感图像的列数；P(p,l)表示去掉边缘信息的局部相位图像素点(p,l)的像素值；(p,l)表示去掉边缘信息的局部相位图像素坐标；
(8)获得遥感图像的模糊度值：
(8a)采用局部相位算法，计算遥感图像的局部相位值，得到局部相位图；
(8b)利用matlab软件中的reshape函数，将局部相位图矩阵拉伸成为一维数组，得到局部相位图的一维数组；
(8c)将局部相位图的一维数组按照从小到大进行排序；
(8d)采用权重函数，计算权重值；
(8e)按照下式，计算遥感图像的模糊度值:
其中，B表示遥感图像的模糊度值；∑表示求和操作；n表示局部相位图的一维数组从1到M×N共M×N个数的标号；M表示遥感图像的行数；N表示遥感图像的列数；Wn表示局部相位图的一维数组中第n个元素的权重值；Pn表示局部相位图的一维数组第n个元素的像素值；n表示数组从1到M×N共M×N个数的标号；
(9)计算遥感图像质量的综合评价值：
(9a)采用极差法，分别对遥感图像的条带噪声值、一般噪声值和模糊度值进行归一化处理；
(9b)将归一化处理后的遥感图像的条带噪声值、一般噪声值和模糊度值按从大到小进行排序；
(9c)根据遥感图像中条带噪声、一般噪声和模糊度在遥感图像质量中的重要程度，确定综合评价因子中各个值得的权重系数；
(9d)按照模糊测度公式，计算按降序排列后条带噪声值、一般噪声值和模糊度值的模糊测度；
(9e)按照下式，计算遥感图像质量的综合评价值：
Q＝I1g1+I2g2+I3g3
其中，Q表示遥感图像质量的综合评价值，Q越大，表示遥感图像的质量越差；I1表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最大值；g1表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最大值的模糊测度；I2表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的中间值；g2表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的中间值的模糊测度；I3表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最小值；g3表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最小值的模糊测度；
(10)输出遥感图像质量的综合评价值。
与现有技术相比，本发明具有以下优点：
第一，由于本发明在获得遥感图像的条带噪声值时，应用遥感图像的相位信息，从含有条带噪声的遥感图像中准确定位条带噪声的位置，克服了现有遥感图像评价技术中没有考虑条带噪声的问题，使得本发明具有在评价遥感图像条带噪声时结果更加精确的优点。
第二，由于本发明在遥感图像综合指标中，将导致遥感图像质量下降的条带噪声、一般噪声和模糊度综合起来，克服了现有技术单纯由人眼判决产生的主观性和仅以单个参数表征遥感图像质量优劣的片面性的问题，使得本发明在评价遥感图像质量时的结果更加全面优点。
第三，由于本发明在遥感图像综合评价指标的建立中采用Choquet模糊积分函数，充分考虑了人的视觉感知，克服了现有技术中遥感图像质量评价综合模型简单，权重系数没有考虑人眼视觉特点的问题，使得本发明在综合评价遥感图像质量时符合人的主观评价的优点。
附图说明
图1是本发明评价遥感图像质量的流程图。
具体实施方式
参照图1，对本发明具体实施方式做进一步的描述。.
步骤1，输入一幅待评价的遥感图像。
在计算机中应用matlab软件分别读取存储在计算机硬盘空间中待评价遥感图像。
步骤2，获得相位一致性值。
采用下述的相位一致性公式，计算遥感图像在条带噪声方向的相位一致性，得到相位一致性值。
其中，P(x,y)表示在遥感图像像素点(x,y)的相位一致性；P(x,y)越接近于1，表示遥感图像在像素点(x,y)的相位一致性越好；∑表示求和操作；o表示Log-Gabor小波滤波器的方向；Eo(x,y)表示对遥感图像进行在Log-Gabor小波滤波方向o上的滤波，得到的偶对称频率响应和奇对称频率响应的平方和，再对其平方和开方；Ado(x,y)表示在Log-Gabor小波滤波器的尺度d和Log-Gabor小波滤波器的方向o上遥感图像像素点(x,y)与Log-Gabor小波滤波器卷积的幅值；ε表示0.01的常数。
步骤3，二值化处理。
采用下述的二值化方法，对相位一致性值进行二值化处理，得到0-1模板。
第一步，用相位一致性值的最大值，减去相位一致性值的最小值，用其差值乘以0.146，将乘积结果作为二值化操作的阈值。
第二步，将相位一致性值中大于二值化阈值的值赋值为1，将相位一致性值中小于二值化阈值的值赋值为0，完成相位一致性值的二值化处理。
步骤4，获得条带噪声强度矩阵。
用遥感图像矩阵乘以0-1模板，得到条带噪声，将条带噪声中的每个像素值除以条带噪声的最大像素值进行归一化处理，得到条带噪声强度矩阵。
步骤5，采用下述的条带噪声公式，计算遥感图像的条带噪声值。
其中，N1表示遥感图像的条带噪声值，N1越大，遥感图像受条带噪声污染越严重，∑表示求和操作，j表示遥感图像像素的列坐标，i表示遥感图像像素的行坐标，S表示条带噪声强度矩阵，P表示遥感图像的相位一致性值，M表示遥感图像的行数，N表示遥感图像的列数。
步骤6，获得一般噪声信息。
采用下述的边缘提取方法，提取遥感图像的边缘，得到边缘图。
第一步，按照下式，计算噪声能量：
To＝2μ1+2σ1
其中，To表示在Log-Gabor小波滤波器o方向上遥感图像的噪声能量，o表示Log-Gabor小波滤波器方向，μ1表示通过Log-Gabor小波滤波器滤波后遥感图像像素的平均值，σ1表示通过Log-Gabor小波滤波器后遥感图像像素的标准差。
第二步，按照下式，计算遥感图像的边缘值：
其中，E表示遥感图像的边缘值，∑表示求和操作，o表示Log-Gabor小波滤波器方向，Eo表示在Log-Gabor小波滤波器o方向上遥感图像和Log-Gabor小波滤波器卷积的幅值，To表示在Log-Gabor小波滤波器o方向上的噪声能量，表示向下取整操作，cos(·)表示余弦操作，Φo表示在Log-Gabor小波滤波器o方向上遥感图像和Log-Gabor小波滤波器卷积后的相位值。
采用形态学中的膨胀运算，对边缘图进行膨胀，得到膨胀后的边缘图。
采用下述的局部相位算法，计算遥感图像的局部相位值，得到局部相位图。
局部相位算法的步骤如下:
第一步，按照下式，计算噪声能量：
Tk＝μ2+2σ2
其中，Tk表示在Log-Gabor小波滤波器k方向上遥感图像的噪声能量，k表示Log-Gabor小波滤波器的方向，μ2表示通过Log-Gabor小波滤波器后遥感图像像素点的平均值，σ2表示通过Log-Gabor小波滤波器后遥感图像像素点的标准差。
第二步，按照下式，计算遥感图像的局部相位值：
其中，M2表示遥感图像的局部相位值，k表示Log-Gabor小波滤波器的方向，∑表示求和操作，Ek表示在Log-Gabor小波滤波器k方向上遥感图像和Log-Gabor小波滤波器卷积后的幅值，Tk表示在Log-Gabor小波滤波器k方向上，遥感图像的噪声能量，表示向下取整操作，cos(·)表示求余弦值操作，Φk表示在Log-Gabor小波滤波器k方向上遥感图像和Log-Gabor小波滤波器卷积后的相位值。
从膨胀后的边缘图中找出像素值为1的坐标，将找出的坐标所对应的局部相位图中像素值赋值为0，得到去掉边缘信息的局部相位图。
步骤7，采用下述的一般噪声公式，计算遥感图像的一般噪声值。
其中，N2表示遥感图像一般噪声的值，N2越大，表示遥感图像受一般噪声污染越严重，∑表示求和操作，p表示去掉边缘信息的局部相位图的行标号，M表示遥感图像的行数，l表示去掉边缘信息的局部相位图的列标号，N表示遥感图像的列数，P(p,l)表示去掉边缘信息的局部相位图像素点(p,l)的像素值；(i,j)表示去掉边缘信息的局部相位图像素坐标。
步骤8，获得遥感图像的模糊度值。
采用下述的局部相位算法，计算遥感图像的局部相位值，得到局部相位图。
局部相位算法具体步骤如下：
第一步，按照下式，计算噪声能量：
Tk＝μ2+2σ2
其中，Tk表示在Log-Gabor小波滤波器k方向上遥感图像的噪声能量，k表示Log-Gabor小波滤波器的方向，μ2表示通过Log-Gabor小波滤波器后遥感图像像素点的平均值，σ2表示通过Log-Gabor小波滤波器后遥感图像像素点的标准差。
利用matlab软件中的reshape函数将局部相位图矩阵拉伸成为一维数组，得到局部相位图的一维数组。
将局部相位图的一维数组按照从小到大进行排序。
采用下述的权重函数，计算权重值。
其中，Wq表示按照升序排列后一维数组中第q个元素的权重值，exp(·)表示以自然常数e为底的指数操作，q表示数组从1到M×N共M×N个数的标号，M表示遥感图像的行数，N表示遥感图像的列数。
按照下述的模糊度公式，计算遥感图像的模糊度值：
其中，B表示遥感图像的模糊度值，∑表示求和操作，n表示局部相位图的一维数组从1到M×N共M×N个数的标号，M表示遥感图像的行数，N表示遥感图像的列数，Wn表示局部相位图的一维数组中第n个元素的权重值，Pn表示局部相位图的一维数组第n个元素的像素值；n表示数组从1到M×N共M×N个数的标号；
步骤9，计算遥感图像质量的综合评价值。
采用极差法，分别对遥感图像的条带噪声值、一般噪声值和模糊度值进行归一化处理。
将归一化处理后的遥感图像的条带噪声值、一般噪声值和模糊度值按从大到小进行排序。
根据遥感图像中条带噪声、一般噪声和模糊度在遥感图像质量中的重要程度，确定综合评价因子中各个值得的权重系数。
根据choquet模糊积分中，权重系数的确定方法，专家按照遥感图像中条带噪声、一般噪声和模糊度在遥感图像质量中的重要程度，得出语气算子，最后获得权重系数。
按照模糊测度公式，计算按降序排列后条带噪声值、一般噪声值和模糊度值的模糊测度。
模糊测度公式如下：
其中，m表示遥感图像中的条带噪声值、一般噪声值和模糊度的标号；g1(Im)表示遥感图像中的条带噪声值、一般噪声值和模糊度值按降序排列后第m位的模糊测度值，∑表示求和操作，g(Im)表示降序排列后第m位的权重系数。
利用choquet模糊积分模型，计算遥感图像质量的综合评价值。
Q＝I1g1+I2g2+I3g3
其中，Q表示遥感图像质量的综合评价值，Q越大，表示遥感图像的质量越差；I1表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最大值，g1表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最大值的模糊测度，I2表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的中间值，g2表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的中间值的模糊测度，I3表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最小值，g3表示归一化处理后的条带噪声值、一般噪声值和模糊度值中的最小值的模糊测度，
步骤10，输出遥感图像质量的综合评价值。