本发明涉及一种基于中频信号调制和补偿的遥感图像融合技术,属于遥感影像处理数据融合技术领域。
背景技术：
卫星遥感技术可以对地球表面进行大规模全覆盖的观测，在测绘、气象、海洋、农业、自然资源调查、灾害监测、国防安全等诸多方面均发挥巨大作用，世界各国也都十分重视。由于传感器信噪比的限制，光学遥感卫星所获取图像的光谱分辨率和空间分辨率之间存在矛盾，因此目前光学遥感卫星一般配置一个高分辨率全色波段和多个低分辨率多光谱波段。
遥感图像融合技术是指将高空间分辨率低光谱分辨率的全色图像与低空间分辨率高光谱分辨率的多光谱图像集成，生成兼具高空间分辨率和高光谱分辨率的遥感图像的图像处理技术，又称为像素级遥感数据融合技术。
以ihs、pca、bt变换为代表的第一代遥感图像融合主要侧重于对多光谱图像进行视觉效果增强，所采用的理论假设较强，光谱扭曲比较严重；以小波变换、hpf等为代表的第二代遥感图像融合面向计算机图像处理和定量遥感应用，采用从全色图像提取空间细节，经调制后注入多光谱图像的策略，光谱保持性能远胜于第一代方法，但容易引起频谱混叠，表现为图像发虚以及空间分辨率增强不足等，全色与多光谱图像之间的配准偏差对融合结果亦会造成较大影响。目前学术界致力于构建第三代遥感图像融合框架，利用稀疏表达、压缩传感等新理论，以及机器学习等新技术，将遥感图像融合问题视为优化问题予以解决。目前第三代融合方法距离达到实用水平还有相当大的差距。
第一代融合方法多属于分量替换方法(componentsubstitution,cs)，即将多光谱数据投影变换到某一特征空间，利用全色数据或其某种变换形式替换其中某一相关性较高的特征，并进行逆变换获得融合图像。第二代融合方法则多采用多分辨率分析技术(multi-resolutionanalysis,mra)，对全色和多光谱数据分别进行空间频谱的分解，然后依据一定的波段间空间结构模型，将全色数据中包含的空间细节成分注入到低分辨率多光谱数据中。
cs方法可以较好的保持全色数据的空间细节信息，但会严重扭曲多光谱数据的光谱信息，因此融合图像的空间性能好，光谱性能差；mra方法则相反，光谱性能好，空间性能差。有相当数量的研究试图将cs和mra方法结合起来，构造混合模型以吸取二者的优势，但通过理论分析可知混合模型通常都与某种单一方法等价。因此也有学者认为，光谱性能和空间性能是图像融合中的一对矛盾。虽然这一猜想源于空间分辨率和光谱分辨率的矛盾，理论上确实存在对二者共同提升的上限，但现有方法呈现出这种矛盾平衡距离理论上限还有相当大的提升空间。提升mra方法的空间性能，或改善cs方法的光谱保持性能，是目前比较可行的解决方案。目前对mra方法的改进，主要着眼于基于mtf的滤波对mra方法几何性能的提高，但实际的mtf值受到诸多因素影响，并且基于mtf滤波的方法并不能解决配准误差的问题。如何构建稳健的、对配准偏差较不敏感的mra方法是一个需要解决的问题。
技术实现要素：
本发明为解决现有基于mra的遥感图像融合方法空间性能较差的问题，在采用适配mtf的滤波器基础上，通过添加空间/光谱补偿项，降低mtf的不确定性对融合结果的影响，减弱图像配准偏差造成的虚轮廓效应；在保持一定的光谱保持能力的基础上，将融合图像的空间分辨率提高到与cs方法相近的程度。
本发明的技术方案为一种基于遥感图像中频信息进行空间细节调制和空间补偿的遥感影像融合方法，包含以下步骤：
步骤1：对原始多光谱图像ms通过立方卷积法上采样至与全色图像p相同的尺寸，得到插值图像msl；
步骤2：根据多光谱传感器相对全色传感器的相对mtf构造高斯低通滤波器g；
步骤3：图像成分提取，包括以下步骤，
步骤3.1：利用g对p进行高通滤波，获取高频成分ph；
步骤3.2：对ph进行低通滤波，获取中频成分pm；
步骤3.3：对msl进行高通滤波，获取高频成分
步骤4：根据
步骤5；根据w将ph注入msl，得到未经空间补偿的融合图像
步骤6：根据
步骤7：根据
步骤1中可采用其他图像插值核，但必须保证采用的是图像插值方法，而不是信号插值方法，否则会造成半个像元的偏移。
步骤2中使用相对mtf值构造高斯滤波器。令全色与多光谱数据的mtf值分别为mp和mm,则相对mtf值m＝mp/mm.由于大部分数据厂商提供的全色图像已经过mtf补偿，所以m可直接采用mm.如果数据厂商未提供数据的mtf值，可通过其他方法对图像mtf值进行测量，或采用典型值，一般为0.15到0.3.
对于给定的m值，滤波器g定义为
其中m为多光谱传感器相对全色传感器的相对mtf，r为多光谱和全色图像像元尺寸的比值，r＞1.给定滤波器抽头数p,则归一化系数为
步骤3.1采用步骤2构造的滤波器进行高通滤波，通过下式实现
步骤3.2和3.3分别对进行高通滤波和相应的低通滤波，给定低通滤波核h,步骤3.2的滤波过程如下：
步骤3.3的滤波过程如下：
h推荐采用g,也可采用其他自定义滤波核。
步骤4根据
其中σm为
步骤5根据遥感图像像素级融合统一理论框架进行细节注入，得到未经空间补偿的融合图像，细节注入方法如下式所示：
步骤6根据
步骤7根据
其中fused是融合图像，k为自定义的补偿强度系数，推荐值为0.5，取值范围为(-1,1).配准偏差较大时可采用较大的k值，光谱性能要求很高时可采用较小的k值。
本发明与现有技术相比，改进主要体现在：首先，本发明充分保留了mra方法的优势，具有良好的光谱保持性能；其次，采用中频空间信息构建波段间关系，提高了空间细节调制的准确度；再次，利用中频空间信息进行信息补偿，采用缺省参数可以有效降低mtf的不确定性对融合结果的影响，提高了融合方法的稳健度；采用小于0的补偿系数可进行光谱信息增强；最后，本发明总体效果可以减弱图像配准偏差造成的虚轮廓效应，在保持良好的光谱保持能力的基础上，将融合图像的空间分辨率提高到与cs方法相近的程度，从而使得本发明的方法在空间性能和光谱性能上达到了较高水平的平衡。
具体实施方式
下面实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。采用真实国产遥感卫星高分二号(gf-2)星载遥感多光谱和全色光图像，多光谱图像包含四个波段(mss1～mss4)，全色图像(pan)为单波段。多光谱图像的空间分辨率为4.0m，大小为2500行×2500列。全色图像空间分辨率为1.0m，大小为10000行×10000列。实施本发明包括以下步骤：
步骤一：上采样多光谱图像：对原始多光谱图像4个波段(mss1～mss4)通过立方卷积法4倍上采样至全色图像pan相同的尺寸，得到插值图像(lms1～lms4)；
步骤二：构造高斯低通滤波器g.多光谱传感器mtf采用典型值0.15，由于全色数据已经过mtf补偿，则相对mtf取值为0.15.根据下式计算滤波器尺度参数σ
其中m＝0.15为多光谱传感器相对全色传感器的相对mtf，多光谱和全色图像像元尺寸的比值r＝4.计算可得σ＝2.4801.由于高斯滤波器是方向可分离的，所以可以先计算一维高斯滤波器，然后相乘得到二维滤波器。如下matlab代码构造了一个19抽头的二维高斯滤波器：
m＝0.15；p＝19；r＝4；
s＝r*sqrt(-2*log(m))/pi；
x＝1:p；
x＝x-(p+1)/2；
h＝1/(sqrt(2*pi)*s)*exp(-x.^2/(2*s^2))；
g＝h'*h；
g＝g/sum(g(:))；
步骤三：图像成分提取，包括以下步骤：
首先，利用步骤二构造的滤波器g对pan滤波，并用pan减去滤波结果，得到pan的高频成分图像panh；同样的方式，用lms1～lms4分别减去各自低通滤波结果，得到多光谱图像的高频成分图像lmsh1～lmsh4.；最后，利用g对panh进行滤波，获得pan的中频成分图像panm.由于高斯滤波器的方向分离特性，这一步骤中的二维滤波步骤也可用步骤二中构造的一维滤波器h在行和列两个方向上分别进行。
步骤四：计算各波段细节注入的调制参数w1～w4.分别统计panm的标准差sp和lmsh1～lmsh4的标准差s1～s4，用sp分别除s1～s4得到w1～w4.
步骤五：空间细节注入获取未经空间补偿的融合图像msu.对lms1～lms4，分别加上panh与对应波段的调制参数w1～w4的乘积，得到初步融合结果fmsu1～fmsu4.
步骤六：构造空间补偿项msc.各多光谱波段的空间补偿项为该波段的高频成分减去该波段对应的调制参数w与全色中频成分panm的乘积，以第一波段为例,msc1＝lmsh1-w1*panm
步骤七：对初步融合结果进行空间补偿，空间补偿方法如下：从未经补偿的初步融合结果fmsu中减去k倍的空间补偿项msc，k为用户自定义的性能参数，缺省值为0.5，较大的k值可以提高空间增强的效果，k<0时可以增强光谱信息削弱空间结构信息。
步骤八：对融合结果进行异常处理并输出。