本发明属于遥感图像处理
技术领域：
，具体来说，涉及一种改进的遥感多光谱与全色图像融合方法。
背景技术：
：由于传感器物理特性和数据传输能力等的限制，高空间分辨率和高光谱空间分辨率影像难以同时获取。因此，大量目前在轨的卫星（例如Landsat7ETM+、QuickBird、SPOT-5和WorldView-2/3等），同时提供高空间分辨率的全色波段（PAN）和低空间分辨率的多光谱波段（MS）遥感图像。由于大量的应用需要使用高空间分辨率的多光谱图像，因此破切地需要将MS与PAN图像进行融合处理，以获得空间分辨率增强的MS图像，以应用于遥感图像解译、地表覆盖分类、目标检测等应用中。近年来国内外研究人员发展了大量遥感融合技术来融合MS和PAN图像以得到高空间分辨率的MS图像。现有的MS与PAN融合方法可以分为三类：基于成分替代的方法、基于PAN调制的方法和基于多尺度分析的方法。基于成分替代融合的典型代表算法有Intensity-Hue-Saturation（IHS）变换、PrincipalComponentAnalysis（PCA）,Gram-Schmidt等，而基于PAN调制技术的融合算法主要有Brovey变换、Pradines’、合成变化比率（SyntheticVariableRatio）、SmoothingFilter-basedIntensityModulation、PANSHARP（PS）、Haze-andRatio-based（HR）等。基于成分替代和调制技术的融合算法的特点是快速且易于实现，但其融合结果会导致一定程度的光谱失真。基于多尺度变换的融合方法，有效地保持了多光谱图像的光谱信息。然而这类方法的融合图像可能会出现空间畸变，典型的现象有振铃效应、虚景混淆、边缘以及纹理模糊。为综合成分替代和多分辨率分析等方法分别在空间和光谱信息保留方面的优势，一些研究人员提出了成分替代（PCA、IHS等）或调整技术与多分辨率分析相结合的融合算法，这类方法能获得优于成分替代和标准多分辨率分析融合的结果，但增加了计算复杂度。有一些研究利用正则化优化的方法对融合结果进行迭代。尽管国内外研究人员提出了大量新的融合方法来降低融合图像光谱失真，如何进一步降低光谱失真仍然是遥感融合研究面临的主要挑战。基于PAN调制的融合方法（例如PCI软件中的PANSHARP算法）由于具有计算简单、鲁棒性强等优点，被广泛地应用于卫星数据的融合。基于PAN调制的MS与PAN融合基于假设：融合后MS波段与原始MS波段的比值等于PAN图像与合成的PAN图像的比值。但当图像受雾气（即Haze）影响的情况下，该假设并不成立。在JingandCheng（2009）的文献中，HR融合图像从目视比较和定量评价指标来看都优于PANSHARP和Gram-Schmidt融合方法，这证明了PAN调制融合方法中考虑雾气的必要性。然而，由于各个波段的雾气值决定了融合像元的光谱矢量的方向，影响了融合图像的光谱失真程度。因此，雾气（Haze）值的取值对融合图像的质量非常重要。此外，在考虑雾气影响的情况下，图像中暗色调像元（特别是水体、阴影等像元）的融合光谱的失真问题是待解决的关键问题。针对相关技术中的问题，目前尚未提出有效的解决方案。技术实现要素：针对相关技术中的上述技术问题，本发明提出一种改进的遥感多光谱与全色图像融合方法，能够显著降低融合图像的光谱失真，特别是图像中对应水体、阴影等暗色调地物像元的融合光谱的失真。为实现上述技术目的，本发明的技术方案是这样实现的：一种改进的遥感多光谱与全色图像融合方法，包括以下步骤：S1将原始MS图像采用立方卷积方式上采样到原始PAN图像的分辨率得到上采样MS图像；将原始PAN图像采用平均法下采样到MS图像空间分辨率，再采用立方卷积方法上采样到原始PAN空间分辨率得到合成PAN图像；S2利用原始PAN图像计算PAN波段的雾气值，利用原始MS图像计算各个MS波段的雾气值；根据PAN波段雾气值和各个MS波段的雾气值去除上采样MS图像、原始PAN图像和合成PAN图像中的雾气，得到图像Ir、图像Pr和图像Psr；S3对图像Psr中的接近0值的像元进行偏移，并根据图像Pr和图像Psr得到比值图像；利用图像Ir、比值图像和各个MS波段的雾气值得到融合图像。进一步的，在步骤S2中，PAN波段的雾气值为Hp=min（P），MS图像第i个波段的雾气值为Hi=min（MSi），其中，P为原始PAN图像，MSi为原始MS图像的第i波段。进一步的，在步骤S2中，根据PAN波段雾气值和各个MS波段的雾气值去除上采样MS图像、原始PAN图像和合成PAN图像中的雾气，得到图像Ir、Pr和Psr的计算公式为：；；，式中，Ii为上采样MS的第i波段，Ps为合成PAN图像。进一步的，在步骤S3中，对图像Psr中的接近0值的像元进行偏移具体包括：根据图像Psr的直方图确定阈值，小于阈值的像元即为图像Psr中接近0值的像元；对图像Psr中接近0值的像元，增加偏移量。进一步的，根据图像Psr的直方图确定阈值具体包括：统计图像Psr的直方图，并计算直方图累积分布，取累积比例Pt所对应的灰度值设置为阈值，其中，累积比例Pt的取值范围为0.01≤Pt≤0.03。进一步的，在步骤S3中，根据图像Pr和图像Psr得到比值图像的计算公式为：，式中，T为阈值，R为比值图像，S为偏移量，S可定义为S=αT,0.5≤α≤2。进一步的，在步骤S4中，利用图像Ir、比值图像和各个MS波段的雾气值得到融合图像的计算公式为：，式中，Fi为融合图像的第i波段。本发明的有益效果：本发明方法的融合图像跟HR方法融合图像相比，进一步降低了光谱失真、且进一步地增强了空间细节。而跟GSA、GLP、ATWT和ATWP等方法相比，本发明方法的融合图像在多数情况下都显著优于这些方法的融合图像。本发明通过对HR方法的改进，进一步降低了融合图像光谱失真并增强了空间细节，特别是能显著降低水体、阴影等像元的融合后的光谱失真，特别适用于城市高分辨率遥感影像的融合。此外，本发明方法具有算法简单、稳健高效等特点，适用于大图像的融合。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1是根据本发明实施例所述的改进的遥感多光谱与全色图像融合方法的流程图。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。如图1所示，根据本发明实施例所述的一种改进的遥感多光谱与全色图像融合方法，包括以下步骤：S1将原始MS图像采用立方卷积方式上采样到原始PAN图像的分辨率得到上采样MS图像；将原始PAN图像采用平均法下采样到MS图像空间分辨率，再采用立方卷积方法上采样到原始PAN空间分辨率得到合成PAN图像；S2利用原始PAN图像计算PAN波段的雾气值，利用原始MS图像计算各个MS波段的雾气值；根据PAN波段雾气值和各个MS波段的雾气值去除上采样MS图像、原始PAN图像和合成PAN图像中的雾气，得到图像Ir、图像Pr和图像Psr；S3对图像Psr中的接近0值的像元进行偏移，并根据图像Pr和图像Psr得到比值图像；利用图像Ir、比值图像和各个MS波段的雾气值得到融合图像。其中，在步骤S2中，PAN波段的雾气值为Hp=min（P），MS图像第i个波段的雾气值为Hi=min（MSi），其中，P为原始PAN图像，MSi为原始MS图像的第i波段。其中，在步骤S2中，根据PAN波段雾气值和各个MS波段的雾气值去除上采样MS图像、原始PAN图像和合成PAN图像中的雾气，得到图像Ir、Pr和Psr的计算公式为：；；，式中，Ii为上采样MS的第i波段，Ps为合成PAN图像。其中，在步骤S3中，对图像Psr中的接近0值的像元进行偏移具体包括：根据图像Psr的直方图确定阈值，小于阈值的像元即为图像Psr中接近0值的像元；对图像Psr中接近0值的像元，增加偏移量。其中，根据图像Psr的直方图确定阈值具体包括：统计图像Psr的直方图，并计算直方图累积分布，取累积比例Pt所对应的灰度值设置为阈值，其中，累积比例Pt的取值范围为0.01≤Pt≤0.03。其中，在步骤S3中，根据图像Pr和图像Psr得到比值图像的计算公式为：，式中，T为阈值，R为比值图像，S为偏移量，S可定义为S=αT,0.5≤α≤2。其中，在步骤S4中，利用图像Ir、比值图像和各个MS波段的雾气值得到融合图像的计算公式为：，式中，Fi为融合图像的第i波段。本发明是在一种HR方法（PAN调制融合方法）的基础上，提出一种改进的遥感多光谱（MS）与全色（PAN）图像融合方法。HR方法是一种考虑了大气程辐射的PAN调制融合方法，该方法基于假设：融合后图像与原始多光谱图像的比值等于PAN图像与合成的PAN图像（空间分辨率同原始MS图像）的比值。根据该方法，融合图像第i波段Fi的计算公式为：，式中，Ii为上采样MS图像的第i波段，P为PAN波段，Ps为合成MS空间分辨率的PAN图像。Hp和Hi分别为PAN波段和MS第i波段的雾气值。在JingandCheng提出的HR方法中，上采样MS图像通过对MS图像进行立方卷积插值得到；Ps则通过对原始PAN波段（即P）采用平均法下采样到原始MS图像空间分辨率，然后对下采样图像采用立方卷积插值得到；雾气值Hp和Hi分别PAN波段和MS图像第i波段的最小值确定。本发明对HR主要进行了两个方面的改进。首先，明确了雾气值Hi和Hp应分别取值为原始低分辨率MS图像的第i波段（即Ii）和高分辨率PAN波段（即P）的最小值；其次，为了避免式中分母（即Ps-Hp）出现零值，针对图像中暗像元的融合处理进行了改进。下面分别介绍这两方面的改进：改进一：雾气值的确定由于算法执行过程中涉及两个MS图像，即原始低分辨率MS图像和上采样MS图像。尽管原始HR方法中指出雾气值Hi和Hp根据波段最小值确定，但并没有明确指出Hi是根据原始MS的第i波段还是上采样MS的第i波段的最小值来确定。根据采用多个传感器高分辨遥感影像的融合实验，本发明确定，雾气值Hi应应取值为原始MS的第i波段的最小值。改进二：暗像元融合处理为了避免式中分母（即Ps-Hp）出现零值，对于去雾气后的图像Psr（Psr=Ps-Hp）中的低值像元进行了偏移。具体步骤为：如果某个像元(m,n)的值Ps（m,n）小于阈值T，则像元（m,n）的融合光谱的计算公式为：，式中，S为低值像元的偏移，可定义为S=αT,0.5≤α≤2。在本发明中，阈值T根据原始PAN图像的直方图累积分布曲线来确定。具体为，将T设置为累积比例为Pt所对应的灰度。其中，Pt的取值范围为0.01≤Pt≤0.03。由于雾气值决定了融合像元光谱矢量的光谱角，因此影响了融合图像的光谱失真程度。在使用来自多个传感器的高分辨率影像的融合实验的结果表明，在采用HR进行融合时，多光谱波段的雾气值Hi取原始低分辨率MS图像第i波段（即MSi）的最小值时，融合图像在质量评价指标上优于Hi取上采样MS图像第i波段（即Ii）的最小值时的融合图像。前者优于后者的原因在于，后者的部分像元是采用插值方式估算的，而前者是完全真实的MS影像，因此后者的波段最小值更为准确。在去除雾气后，在PAN图像去除雾气后的图像PR和合成PAN图像去除雾气后的图像PSR中会出现较多0值附近的像元，这些像元主要是图像中的阴影、水体像元。令，对于像元(m,n)，如果Ps（m,n）-Hp（m,n）的计算值特别接近0，将出现R(m,n)的值特别大或特别小的现象，从而导致融合像元的光谱失真。而Ps（m,n）-Hp（m,n）的计算值进行偏移，即，将能避免Rt（m,n）的值过高或过低，从而避免引入光谱失真。因此本发明能显著改进水体、阴影等融合像元的光谱失真现象。为了评价本发明方法的性能，我们设计了对比实验。实验数据包括来自WorldView-2、Pleiades、IKONOS等3个传感器的5个高分辨率遥感图像；对比方法选择了5种目前公认比较优秀的融合算法，具体包括HR、PANSHARP、AdaptiveGram-Schmidt（GSA）、GeneralizedGaussianGeneralized（GLP）、“àtrous”wavelettransform（ATWT）和AdditiveWaveletLuminanceProportional（AWLP）；融合图像质量评价指标选用了相对全局维数综合指标（ERGAS）、光谱角（SAM）、综合质量指数Q4/Q8、和空间相关系数（SCC）。其中，EASE反映了融合图像跟参考图像的偏差，值越小融合效果越好；ERGAS反映了融合图像跟参考图像的全局光谱辐射变形误差，越小越好；SAM反映融合图像跟参考图像的光谱差异，值越小融合效果越好；Q4与Q8是同时考虑了融合图像跟参考图像的局部均值偏差、对比度变化以及相关性丢失情况的综合质量指标，值越大越好；SCC是考虑了融合图像跟PAN图像空间细节相关性的指标，值越大越好。5个实验图像的融合图像质量评价指标的统计见表1-表5。表1WorldVIew-2卫星图像1的融合图像的质量评价指标统计IndexRASE(%)ERGASSAM(º)Q8SCC本发明方法6.5401.7082.2610.9370.936HR7.4931.9572.2810.9250.929GSA8.5652.2362.7960.8960.902GLP7.2831.9002.5500.9240.922ATWT7.5481.9662.5670.9210.919ATWP7.6701.9782.6110.9200.917表2WorldVIew-2卫星图像2的融合图像的质量评价指标统计IndexRASE(%)ERGASSAM(º)Q8SCC本发明方法4.7861.2591.7910.9760.978HR8.1782.1401.7850.9620.954GSA7.3291.8492.3250.9550.958GLP5.4571.4171.9510.9700.970ATWT5.8391.5081.9980.9670.968ATWP6.1811.5812.2460.9660.966表3Pleiades卫星图像1的融合图像的质量评价指标统计IndexRASE(%)ERGASSAM(º)Q4SCC本发明方法7.2281.7811.6370.8880.890HR7.3901.8221.6530.8840.886GSA8.0232.0711.6410.8460.857GLP6.9781.7531.4840.8830.885ATWT6.8221.7121.4660.8850.888ATWP6.8911.7021.4270.8850.887表4Pleiades卫星图像2的融合图像的质量评价指标统计IndexRASE(%)ERGASSAM(º)Q4SCC本发明方法8.7922.1681.8750.8630.867HR9.1202.2471.9200.8560.861GSA10.1612.6221.7240.7970.818GLP7.9011.9861.5470.8690.875ATWT7.7161.9391.5290.8710.877ATWP7.8751.9401.5230.8690.876表5IKONOS卫星图像的融合图像的质量评价指标统计IndexRASE(%)ERGASSAM(º)Q4SCC本发明方法5.681.491.820.89010.892HR5.721.501.810.88980.892GSA8.632.252.230.8230.833GLP6.191.621.990.8760.877ATWT6.321.662.010.8700.871ATWP6.331.662.050.8710.870上述表中的统计指标可见，从光谱质量评价指标（RASE、ERGAS、SAM和Q4/Q8等）和空间质量评价指标（SCC）来看，本发明方法都优于HR方法。本发明方法的融合图像跟HR方法融合图像相比，进一步降低了光谱失真、且进一步地增强了空间细节。而跟GSA、GLP、ATWT和ATWP等方法相比，本发明方法的融合图像在多数情况下都显著优于这些方法的融合图像。由此可见，借助于本发明的上述技术方案，通过对HR方法的改进，进一步降低了融合图像光谱失真并增强了空间细节，特别是能显著降低水体、阴影等像元的融合后的光谱失真，特别适用于城市高分辨率遥感影像的融合。此外，本发明方法具有算法简单、稳健高效等特点，适用于大图像的融合。以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页1 2 3