一种细节增强的遥感图像数据时空谱一体化融合方法与流程

本发明属于遥感图像融合技术领域，具体涉及一种遥感图像数据时空谱一体化融合方法。

背景技术：

现有的遥感图像融合方法采用三种不同的层级，分别为数据级别、特征级别和决策级别，根据目标的不同，遥感图像融合可分为多视空间融合，空谱融合和时空融合。

大多数现有的遥感图像数据融合方法只能应用于整合空间、时间和光谱分辨率其中的两个指标，因此不能获得最高的时间、空间、光谱分辨率融合图像。此外，一些遥感图像数据融合方法旨在合并来自一个或两个传感器的补充信息，不能充分利用来自更多传感器的互补观测信息。在一体化融合理论方法发展中，国内学者做出了较大贡献。从2010年开始，武汉大学在国家自然科学基金的资助下开始相关研究，首次提出了时空谱一体化融合的概念与基于最大后验概率理论的融合框架，并取得了一定的研究成果。香港中文大学对遥感图像数据的时、空、谱关联模型进行了进一步探索，同样基于最大后验估计理论，提出了一种新的时空谱一体融合方法，但仍然是针对两个传感器数据的融合方面。meng等人、wu等人在3个传感器融合方面取得了重要进展，但没有全面考虑时、空、谱特征。之后meng等人进一步提出了一个统一的融合框架，既能够有效融合时、空、谱互补信息，同时消除了传感器数量的限制，然而该方法嵌入的先验模型为gauss-markov模型，虽然能够有效的平滑图像噪声，但该模型也会平滑掉图像细节信息，进而影响了融合图像的细节信息，造成图像拍摄的视觉效果不佳。

技术实现要素：

本发明的目的是为解决现有的遥感图像数据融合方法在平滑图像噪声的同时，也会平滑掉图像细节信息的问题。

本发明为解决上述技术问题采取的技术方案是：

一种细节增强的遥感图像数据时空谱一体化融合方法，该方法包括以下步骤：

步骤一、采用拉普拉斯算子对多视图像y和辅助多源观察图像z进行边缘增强，得到增强后的多视图像y′和辅助多源观察图像z′；

步骤二、建立融合图像x与y′之间的空间退化模型和融合图像x与z′之间的时空谱关系模型，计算出y′对融合图像x的一致性约束p(y′|x)、z′对融合图像x的一致性约束p(z′|x)和图像空间关系的描述p(x)；

步骤三、根据步骤二计算出的p(y′|x)、p(z′|x)和p(x)，得到目标函数f(x)的表达式；

步骤四、通过共轭梯度算法求解目标函数f(x)，并利用目标函数f(x)的求解结果计算融合图像x。

本发明的有益效果是：本发明提供了一种细节增强的遥感图像数据时空谱一体化融合方法，本发明首先采用拉普拉斯算子对多视图像和辅助多源观察图像进行边缘增强，因此能够在平滑图像噪声的同时保留图像细节信息；然后利用融合图像与增强后的多视图像之间的空间退化模型和融合图像与辅助多源观察图像之间的时空谱关系模型，计算增强后的多视图像对融合图像的一致性约束、辅助多源观察图像对融合图像的一致性约束和图像空间关系的描述，并进一步得到目标函数；最后通过共轭梯度算法求解目标函数，利用目标函数的求解结果计算出融合图像。与现有方法相比，本发明的方法在平滑掉图像噪声的同时，可以保留95％以上的图像细节信息。

附图说明

图1是本发明的一种细节增强的遥感图像数据时空谱一体化融合方法流程图；

具体实施方式

具体实施方式一：如图1所示，本实施方式所述的一种细节增强的遥感图像数据时空谱一体化融合方法，该方法包括以下步骤：

步骤一、采用拉普拉斯算子对多视图像y和辅助多源观察图像z进行边缘增强，得到增强后的多视图像y′和辅助多源观察图像z′；

步骤二、建立融合图像x与y′之间的空间退化模型和融合图像x与z′之间的时空谱关系模型，计算出y′对融合图像x的一致性约束p(y′|x)、z′对融合图像x的一致性约束p(z′|x)和图像空间关系的描述p(x)；

步骤三、根据步骤二计算出的p(y′|x)、p(z′|x)和p(x)，得到目标函数f(x)的表达式；

步骤四、通过共轭梯度算法求解目标函数f(x)，并利用目标函数f(x)的求解结果计算融合图像x。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一的具体过程为：

多视图像y＝{y1,...,yk,...}，yk是第k个观测图像，k是多视图像的总数；为了获取k时刻的高分辨率融合图像，前提是在给定的时间k时刻有观测图像yk，一般来说，图像y具有更高的光谱和时间分辨率，但空间分辨率较低。因此，我们认为y是空间退化的图像。相对应的的增强融合过程得到的x是具有较高空间分辨率的图像。

辅助多源观察图像z＝{z1,z2...,zn,...zn}，其中：zn代表第n个图像，n是辅助多源观察图像的总数；通常是一个全色图像或多光谱(ms)，具有更高的空间分辨率，但光谱或时间分辨率较低。一体化融合框架能完成不同类型的融合任务，包括多视图空间融合，空间谱融合，时空融合，时空谱融合。

拉普拉斯算子(laplace)是一个二阶线性微分算子，与一阶线性微分算子相比，二阶微分的边缘定位能力更强，锐化效果更好。假设多视图像y为二维图像，y的函数表达式为f(x,y)，x,y分别是二维图像的横坐标和纵坐标，laplace算子定义为：

对于任意阶微分算子都是线性算子，所以二阶微分算子和后面的一阶微分算子都可以用生成模板然后卷积的方式得出结果。根据微分的定义有：

将公式(2)、(3)与拉普拉斯算子定义相结合得到：

将拉普拉斯算子表示成模板的形式以便编程需要，则拉普拉斯算子的扩展模板为

从模板形式可以看出，如果在图像中一个较暗的区域中出现了一个亮点，那么使用拉普拉斯算子运算可以使这个亮点变得更亮。因为图像中的边缘就是那些灰度发生跳变的区域，所以拉普拉斯锐化模板在边缘的细节增强过程中有很大作用。

最后将计算像素与其周边像素的加权差值然后再叠加到原图中，得到细节增强后的图像。则增强后的多视图像y′的函数表达式g(x,y)为：

其中：c为权重系数；c的取值和上文的模板定义有关系，当模板中心数值取正时，c＝-1，相反c＝1；

同理，得到增强后的辅助多源观察图像z′。

具体实施方式三：本实施方式与具体实施方式二不同的是：步骤二的具体过程为：

融合图像x与y'之间的空间退化模型表示如下：

y′k,b＝dsk,bmkxb+vk,b,1≤b≤bx,1≤k≤k(6)

y′k,b表示y′的第k个图像的第b个波段的退化观测图像，xb代表融合图像x的第b个波段，bx为波段的总数，mk表示运动矩阵，sk,b表示光学模糊矩阵，d是下采样矩阵，vk,b是由传感器和外界环境引起的零均值高斯噪声；

对公式(6)进行简化，则空间退化模型表示为：

y′k,b＝ay,k,bxb+vk,b(7)

其中：ay,k,b表示融合图像x和y′之间的空间退化矩阵图像，ay,k,b＝dsk,bmk；

融合图像x与z′之间的时空谱关系模型表示如下：

z′n,q＝ψn,qcn,qaz,n,qxq+τn,q+vn,q,1≤q≤bz,n,1≤n≤n(8)

z′n,q表示z′的第n幅图像的第q个波段的退化观测图像，xq代表融合图像x的第q个光谱带，bz,n为光谱带的总数，az,n,q表示融合图像x和z′之间的空间退化矩阵图像，cn,q表示光谱相关性矩阵，ψn,q为时间相关矩阵，τn,q为时偏移量，vn,q代表零均值传感器噪声；

则融合图像x的估计值表示如下：

p(x|y′,z′)代表融合图像x对y'和z′的一致性约束；代表当p(x|y′,z′)取最大值时，对应的x的取值；

根据贝叶斯公式得：

其中：p(y′|x)表示y'对融合图像x的一致性约束，p(z′|x)表示z′对融合图像x的一致性约束，p(x)代表图像空间关系的描述；

假设传感器和外界环境引起的零均值高斯噪声vk,b服从随机高斯分布，则p(y′|x)表示如下：

其中：p(y′k,b|xb)表示y′的第k个图像的第b个波段的退化观测图像对融合图像x的第b个波段的一致性约束；ay,k,b是vk,b的方差，bx是光谱波段号，φ1φ2表示y′k,b的空间维度，||·||2表示2范数；

假设vn,q服从随机高斯分布，则p(z′|x)表示如下：

其中：p(z′n,q|xq)表示z′的第n幅图像的第q个波段对融合图像x的第q个波段的一致性约束；az,n,q表示噪声vn,q的方差，h1h2表示z′n,q的空间维度；

第三个概率密度函数p(x)是图像先验知识，用来描述图像的空间关系。我们引入huber-markov先验模型。相比于gauss-markov模型，该模型可以在融合过程中消除噪声带来的不稳定现象，并保持影像边缘，即细节增强。基于拉普拉斯算子的三维空谱的自适应加权，则p(x)表示为：

其中：ax,b是随机高斯分布的噪声vx,b的方差，l1l2表示空间维度，ρ(·)为huber函数，qxb表示自适应加权三维拉普拉斯算子矩阵。且qxb表示为：

是根据重采样得到的融合图像的初始估计值，β是参数，由(21)表示：

其中μ为阈值参数。假设融合图像只有几个光谱波段，其光谱曲线是不连续的，因此令β＝0是可行的。相反，如果所需的图像是一个高光谱图像，可以假定光谱曲线是连续的，因此，自适应加权的先验项可以有效地保持光谱的曲线，降低光谱失真。假设不同的光谱波段之间的区别很小，光谱约束就越强。表示第b个波段的光谱维度的梯度。

具体实施方式四：本实施方式与具体实施方式三不同的是：步骤三的具体过程为：

将公式(11)、(13)和(15)代入公式(10)，根据对数函数的单调性进行处理，则目标函数f(x)表示为正则最小化问题，即根据对数函数的单调性并进行简化处理，很多参数都可以删除，最后的目标函数可以表示为正则最小化问题，融合图像x的估计值的表达式为：

其中：表示当f(x)取最小值时，对应的x的取值为估计值

其中第一项表示x和y′之间的一致性约束，第二项表示x和z′之间的关系，第三项表示先验图像。λ1和λ2分别表示各部分的权重系数，λ1和λ2与噪声方差相关，ωn,q表示z′n,q对融合图像x的贡献；它是自适应计算ωn,q＝λ′n,qun得到的。un是由z′n和x的相关性计算得到。它是假定相关性越大，贡献越大。辅助参数λ′n,q为自适应地调整融合图像各个光谱带的空间细节平衡，表示为每个波段的融合图像，表示为

λ′n,q＝f(zn,q,x)/min[f(z'1,2,x),...f(z'1,q,x)...f(z'n,q,x)]

其中f(z'n,q,x)表示融合的z'n,q光谱带数量。λ1和λ2与噪声方差相关，参数用来调节三个部分的比重。

具体实施方式五：本实施方式与具体实施方式四不同的是：步骤四的具体过程为：

通过共轭梯度算法求解目标函数f(x)的过程如下：

对目标函数f(x)的xb求导得：

表示对目标函数f(xb)求导，为的转置，表示融合图像x的第b个波段对应的光谱相关性矩阵；

融合图像x的第b个波段的估计值通过连续迭代得到：

xb,d+1＝xb+θdeb,d(19)

其中，xb,d+1表示第d+1次迭代后融合图像x的第b个波段，eb,d表示第d次迭代的搜索方向，搜索方向的初始值为即对于第1次迭代，搜索方向为接下来的搜索方向与当前迭代以及之前的搜索方向相关，θd为迭代步长；

f(xb)d为第d次迭代对应的目标函数；则第d+1次迭代的搜索方向eb,d+1为：

中间变量

利用第d+1次迭代的搜索方向eb,d+1计算出第d+2次迭代后融合图像x的第b个波段xb,d+2，利用xb,d+2计算出f(xb)d+2；

以此类推，计算出每次迭代对应的目标函数；利用目标函数f(x)的求解结果计算融合图像x的具体过程为：将每次迭代对应的目标函数分别代入公式(16)，计算出每次迭代的目标函数对应的融合图像x的第b个波段的估计值，直至时，则停止迭代，将作为融合图像x的第b个波段的估计值，其中ζ为迭代终止阈值；

同理，计算出融合图像x的其余波段的估计值，得到融合图像x。