基于弥散张量引导的高分辨率遥感影像融合评价方法与流程

本发明属于遥感影像融合技术领域，尤其涉及一种基于弥散张量引导的高分辨遥感影像上采样融合方法。

背景技术：

多光谱遥感影像在农业、森林、矿藏、环境等领域发挥着重要作用，然而由于辐射传输过程和传感器工艺的限制，其空间分辨率较低，无法较好地获取地物细节。另一方面，高分辨率遥感卫星能够获取亚米级的对地观测，其影像却只能提供全色波段的信息。因此，卫星平台往往同时搭载多光谱和全色高分传感器，通过影像融合的手段对地面的多光谱信息和细节信息进行有效综合，为后续应用提供更有价值的数据。

影像融合一般流程包括多光谱影像上采样、影像变换、灰度匹配、成分替换和影像反变换。现有的研究表明，目前绝大多数影像融合算法在尽量保持融合影像与原始多光谱影像光谱特征一致的同时，其空间细节的表现力会有一定程度下降。原因在于，大多数方法需要在预处理步骤中将多光谱影像上采样到与全色影像一致的分辨率，然而常用的上采样方法没有顾及多光谱影像和全色影像在几何上的一致性，导致重叠的影像并没有严格对应，在地物边缘和细小地物附近经常出现虚假结构。

现有遥感影像获取平台向着多星、多传感器、高空间分辨率、高光谱分辨率和短回访周期发展，导致获取的遥感数据呈几何增长，传统的CPU并行编程已经无法满足应用系统对数据处理的效率和精度的要求。

因此本领域亟待能够实用的相关技术方案出现。

技术实现要素：

针对现有技术存在的问题，本发明提出一种基于散张量引导的高分辨率遥感影像融合评价方法。

本发明技术方案提供一种基于散张量引导的高分辨率遥感影像融合评价方法，包括以下步骤：

步骤1，将经过影像配准的多光谱影像与全色影像分块；

步骤2，上采样，包括以下子步骤，

步骤2.1，从步骤1所得结果中选择一对待融合的影像块，用于训练参数；根据预设的参数α、β和γ初始值，采取各向异性模型进行弥散张量引导的上采样，自适应检测待融合影像边缘并采用最小二乘法拟合；

采取各向异性模型进行弥散张量引导的上采样时，利用以下方程描述多光谱遥感影像的点扩散和上采样过程，

M(k)＝X(k)·D·A(k)

k＝1,…,K

其中，变量k为多光谱影像的波段序号，K为参与融合的波段数，M为原始的多光谱影像，X为多光谱影像的高分辨率估计值，A为二维点扩散函数对应的变换，D为二维上采样过程对应的变换；

采用高斯滤波器平滑原始的全色影像，计算得到每个像素的梯度方向n＝d/||d||和梯度强度||d||，其中d为计算的梯度向量；对于梯度为零的情况，规定n＝[1,0]^T，避免梯度方向的未定性；

定义每个像素的弥散张量t如下，

$t={\mathrm{wdd}}^T+d_{\⊥}{d}_{\⊥}^T$

其中，d_⊥定义为与梯度向量d垂直，即灰度变化最小方向的向量，d^T、分别为向量d、d_⊥的转秩，边缘保护函数w为

w＝exp(-β||d||^γ)

将张量模型用于能量函数优化，得到

$X^{*}\left(k\right)=\underset{X\left(k\right)}{\arg }min\left(\right|\left|S\right(k)\·X(k)-M(k\left)\right|{|}^2+\mathrm{\α}\left|\right|T^{1/2}\·G\·X\left(k\right)\left|{|}^2\right)$

k＝1,…,K

其中，S＝D·A，G为梯度算子对应的矩阵变换，T为各像素的弥散张量t对角排列组成的灰度张量矩阵，α为各向异性正则约束的权值，X为待解求影像，M为原始的多光谱影像；X^*为所得融合影像，K为参与融合的波段数；

基于线性最小二乘方式，计算如下，

X^*＝(S^TS+αG^TTG)^-1S^TM

通过逐波段求解二次优化问题，得到上采样的多光谱影像；

步骤2.2，进行客观指标评价，如果是精度最优则进入步骤3，否则进入步骤2.3；

步骤2.3，调节参数α、β和γ的当前取值，返回步骤2.1；

步骤3，选择遥感影像融合算法；

步骤4，采用GPU并行计算加速，以块的形式生成和保持融合结果，包括读取与进程数目相同的数据块至内存，每个进程分别处理一对待融合的图形块；

步骤5，遥感影像融合评价指标相关性分析，包括因子聚类分析客观评价指标相关性，选取独立的评价指标从多方面综合评价融合结果精度。

而且，步骤3中，

设将配准后的多光谱影像M上采样至分辨率与全色影像相同的处理表达如下，

M^L＝usp(M)

其中，usp()表示上采样算法，M^L为上采样的结果，用和分别表示融合前后的多光谱影像，用pan_(i,j)表示参与融合的全色影像，用k表示参与融合的波段，(i,j)表示影像像素位置；

提出遥感影像像素级融合算法的通用模型如下，

$X_{(k,i,j)}^H=M_{(k,i,j)}^L+F_{(k,i,j)}\·S_{(i,j)}$

其中，融合策略用F_(k,i,j)表示，全色影像pan_(i,j)各像素在(i,j)位置的高频空间细节信息用S_(i,j)表示。

和现有技术相比，本发明具有如下优点和有益效果：

(1)上采样过程中可以通过调节β、γ、α三个独立参数的值，调节弥散张量模型对影像像素的作用强度，使得边缘保护和去除噪声根据实际需求达到平衡。上采样不受原始影像亮对比度等成像质量的影响，可省去遥感影像预处理过程。所得的上采样结果顾及多光谱影像和全色影像在几何上的一致性，使重叠的影像严格对应，避免了在边缘附近像素的虚假结构现象。

(2)通用模型概括的融合策略可清晰的反映各融合算法对光谱保持和空间细节增强的不同偏重。通过相关性分析选取的评价指标相互独立，可从不同维度综合评价所得融合影像质量。

(3)充分发挥现有设备性能，不同于传统方法，本实施例将多光谱影像的上采样和融合过程分配给GPU，利用GPU强度大的计算优势多线程处理，可大幅减少运行时间。且影像分块可避免大幅面影像无法读入内存问题。

附图说明

图1是本发明实施例的流程图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

参见图1，本发明实施例提供一种散张量引导的高分辨率遥感影像融合评价方法，包括以下步骤：

步骤1，将经过影像配准的多光谱影像与全色影像分块，多光谱影像的影像块尺寸为512×512像素，全色影像的影像块尺寸为2048×2048像素。多光谱影像的影像块和全色影像的影像块一一对应。实施例中的输入数据已经严格配准，可以创建存储待融合文件，将多光谱与全色影像各波段数据按分块结果存储于文件中，便于后续融合策略模块化复用。

步骤2，上采样，包括以下子步骤：

步骤2.1，从步骤1所得结果中选择一对待融合的影像块，以便训练参数；根据预设的参数α、β和γ初始值，采用弥散张量引导的上采样方法(disffusion tensor driven，DTD)，采用各向异性假设，自适应检测待融合影像边缘并采用最小二乘法拟合，在保持多光谱影像和全色影像的边缘一致性同时，提高了后续融合影像的质量和细节表现力。

具体实施时，可从步骤1分割所得的多光谱影像的影像块和全色影像的影像块中，由用户指定或随机选取一对位置相应的影像块。本实施例推荐的参数初始值α＝0.5，β＝50，γ＝1.5。

步骤2.2，进行客观指标评价，如果是精度最优则进入步骤3，否则进入步骤2.3；

具体实施时可根据实际需求采用SAM，RMSE，UIQI，PSNR，SD等客观评价指标。

步骤2.3，调节参数α、β和γ的当前取值，返回步骤2.1。具体实施时，本领域技术人员可自行设定调节方式，例如α、β和γ的调节步长分别为0.1、1、0.1，可根据实际情况增大或减小。

实施例提供的DTD上采样方法，实现说明如下：

采取各向异性模型进行弥散张量引导的上采样时，利用以下方程描述多光谱遥感影像的点扩散和上采样过程。

M(k)＝X(k)·D·A(k)

k＝1,…,K

其中，变量k为多光谱影像的波段序号，K为参与融合的波段数，M为配准后原始的多光谱影像，X为多光谱影像的高分辨率估计值，A为二维点扩散函数对应的变换，D为二维上采样过程对应的变换。两个变换都是线性的，该线性方程系统用来描述高分辨率场景在传感器中发生点扩散并以较低分辨率采样后得到多光谱影像的过程。

采用一个高斯滤波器平滑原始全色影像，以减弱噪声对梯度计算的影响。计算得到每个像素的梯度方向n＝d/||d||和梯度强度||d||，其中d为计算的梯度向量。对于梯度为零(完全均质)的情况，规定n＝[1,0]^T，避免梯度方向的未定性。

对于边缘锐利的区域，自适应平滑能够相应地降低约束强度。各向异性扩散张量作为偏微分方程遥感影像处理的重要概念，对于这一问题有很强的适用性。定义每个像素的弥散张量t：

$t={\mathrm{wdd}}^T+d_{\⊥}{d}_{\⊥}^T$

其中，d_⊥定义为与梯度向量d垂直，即灰度变化最小方向的向量，d^T分别为向量d d_⊥的转秩，边缘保护函数w为

w＝exp(-β||d||^γ)

当边缘两侧灰度差越大，w的取值越小，即平滑约束越弱。通过设置参数β和γ的取值可以使约束项对超过规定强度的边缘才发生有效响应。

将张量模型用于能量函数优化，得到：

$X^{*}\left(k\right)=\underset{X\left(k\right)}{\arg }min\left(\right|\left|S\right(k)\·X(k)-M(k\left)\right|{|}^2+\mathrm{\α}\left|\right|T^{1/2}\·G\·X\left(k\right)\left|{|}^2\right)$

k＝1,…,K

其中S＝D·A，D、A分别为上文中的点扩散和上采样函数，G为梯度算子对应的矩阵变换，T为各像素的弥散张量t对角排列组成的灰度张量矩阵，α为各向异性正则约束的权值，X为待解求影像，M为原始的多光谱影像。X^*为所得融合影像，K为参与融合的波段数。

上式为经典的二次规划问题，有全局最优解，使求得的残差最小，即线性最小二乘问题，计算过程由下式给出：

X^*＝(S^TS+αG^TTG)^-1S^TM

定义的灰度张量矩阵T满秩，因此逆矩阵存在。通过逐波段求解二次优化问题，得到上采样的多光谱影像。该融合影像结果既顾及了点扩散过程的模糊效应，也保持了与高分辨率全色波段的边缘一致性。可较大程度减少边缘处的虚假结构现象，避免后续融合操作中由于边缘未对齐导致的高频全色影像信息注入错误现象。

实施例用各项异性的弥散过程来描述多光谱遥感影像像素与高分辨率全色影像像素间的一对多关系。通过计算多光谱遥感影像像素灰度梯度来检测边缘，构建弥散张量，用于高分辨率遥感影像边缘细节的检测，为后续边缘处采取不同的上采样策略提供判定依据。将张量模型用于能量函数优化，获得全波段残差最小的解，提高影像边缘拟合的精度。

步骤3，遥感影像融合算法的选取：分析常用融合算法的数学模型，以及各项客观评价指标之间的相关性，组成了一套面向对象的遥感影像融合系统，可按实际需求对光谱保持和空间增强的不同偏重灵活选择融合方法。

采用通用模型概括融合策略，根据实际情况选择光谱信息保持和空间细节增强的不同偏重。

设经步骤2将配准后原始的多光谱影像M上采样至分辨率与全色影像相同的处理，可表达为：

M^L＝usp(M)

usp()表示上采样算法，M^L为上采样的结果，用和分别表示融合前后的多光谱影像，用pan_(i,j)表示参与融合的全色影像，用k表示参与融合的波段，(i,j)表示影像像素位置。可提出遥感影像像素级融合算法的通用模型：

$X_{(k,i,j)}^H=M_{(k,i,j)}^L+F_{(k,i,j)}\·S_{(i,j)}$

上式中，融合策略用F_(k,i,j)表示，全色影像pan_(i,j)各像素在(i,j)位置的高频空间细节信息用S_(i,j)表示。

通用模型的建立可以清晰的比较各种融合策略的优势与局限，在选取融合方案时思想清晰。最终融合结果由采用的融合策略系数F_(k,i,j)与提取的高频空间细节信息参数S_(i,j)两者共同决定。

分析融合参数F_(k,i,j)和S_(i,j)的过程中可能需要使用的信息有：参与融合的多光谱遥感影像各波段信息；全色影像(i,j)像素的位置信息；以及可能用到的(i,j)邻域像素信息和基于全局的统计信息。

遥感影像融合通用模型在应用方面的优势表现为以下两点：一是有利于概括不同融合算法的数学思想，便于在后续的融合中针对各自优势按需选择融合策略；二是在算法的执行层面有利于模块的集成与复用。

例如以下三种融合方法：

A.基于变量替换的融合方法

基于变量替换的融合方法中，S_(i,j)可概括为提取的高分全色波段的高频细节信息，通常由全色波段直方图变化增强或经过正交变换后多光谱影像某个分量在数据空间的差构成，因此在融合策略中S_(i,j)起到主要作用，对后续融合结果的影像较大；F_(k,i,j)通常由矩阵变换构成。

B.基于调制信息的融合方法

S_(i,j)通常由高分辨率全色波段影像与其滤波后影像求差值计算得出，另一种情况时可由全色波段与其相似的合成波段求差计算得出；F_(k,i,j)可概括为全色影像与一些合成波段的比值；此类方式的融合结果可根据S_(i,j)、F_(k,i,j)的不同偏重满足获取多光谱或高空间分辨率的不同需要。

C.基于多尺度分析的融合方法

用F_(k,i,j)表示高分全色影像和其分级降采样后的近似分量的差值，级数越大，F_(k,i,j)相对于S_(i,j)权值越小，表现为空间细节保持更好，但分解级数的增加会造成多光谱信息的丢失。

步骤4，GPU并行计算加速：采用GPU并行计算加速，以块的形式生成和保持融合结果。整个流程中仅读取与进程数目相同的数据块至内存，每个进程分别处理一对待融合的图形块。与单次导入导出全幅影像相比，大大减少内存消耗。

统一计算架构CUDA(Compute Unified Device Architecture)为Nvidia公司开发的GPU并行编程平台，在CUDA环境下可同时启用远多于CPU的处理线程，极大的提高计算效率。在CUDA环境下，利用GPU强度大的计算优势多线程处理，可大幅减少运行时间。

多核计算平台下，一个计算核心对应并行环境中的一个进程，每个进程通过进程号0从k-1标识。输出融合影像按照行优先顺序划分并编号，行数为m，列数为n，按0到m×n-1编号，则并行计算过程中第i块影像分配的进程号为i％(k-1)+1，每个进程处理的影像块数为(m×n)/(k-1)。

步骤5，遥感影像融合评价指标相关性分析：因子聚类分析客观评价指标相关性，选取独立的评价指标从多方面综合评价融合结果精度。

通过计算各指标间的相关系数建立相关系数矩阵，用因子分析法计算出彼此相互独立因子进行聚类分析，最后将所得结果分类可得在质量评价某一方面具有代表性的独立客观指标，用于后续融合后遥感影像质量评价。根据该系统获取彼此独立的典型客观指标，建立基于综合得分客观评价系统，用于从各方面全方位定性综合评价融合结果质量。

5.1计算皮尔森相关系数

采用多元统计分析中常用的皮尔森相关系数反应变量的线性相关程度

${\mathrm{\ρ}}_{x,y}=\frac{\mathrm{cov}(x,y)}{\mathrm{\σ}x\mathrm{\σ}y}$

其中x、y为变量，cov(x,y)表示变量x、y的协方差，σxσy表示变量x、y标准差的积，由柯西-施瓦茨不等式可知，当两个变量的线性关系增强时，相关系数趋近于1或-1。

5.2构建相关系数矩阵

因子分析的数学模型可概括为

X＝AF+ε

其中X为原始变量向量，A为公因子荷载系数矩阵，由因子荷载a_ij构成，a_ij为第i个变量与第j个公共因子的相关系数。F为因子向量，ε为残差向量，每个变量可用公共因子的线性函数与残差向量之和表示。

5.3因子聚类分析获取评价指标

通过因子分析得到独立性较强的因子，当因子包含的指标较多是对因子进行聚类分析，通过细化分类求出每个类别的典型指标。本实施例采用平均结连法，用欧式距离得到类内紧凑、类间独立的聚类结果。平均结连法为现有技术，本发明不予赘述。

以上步骤所得的相互独立评价指标用于后续遥感影像融合结果各维度的定量评价。具体实施时，本领域技术人员可采用计算机软件方式实现以上流程的自动运行。

本发明中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。