CN103177431A - 一种多源遥感数据时空融合方法 - Google Patents

Abstract

本发明公开了一种多源遥感数据时空融合方法，特别是利用中空间分辨率影像与高时间分辨率影像融合出高时间分辨率和高空间分辨率影像的融合方法。该方法是中、高空间分辨率影像成像机理研究基础上，基于线性光谱混合模型，从时序低空间分辨率影像中提取地物反射率的时间变化特征；利用地物反射率的时间变化特征，结合前期高空间分辨率影像的纹理特征，采用反射率融合模型，融合出空间分辨率与中分辨率影像相同，时间分辨率与高时间分辨率影像相同的高时空分辨率融合影像；利用真实中分辨率影像，采用相关系数、方差、平均绝对误差、RMSE和偏差等参数对融合影像，对融合影像进行精度评价。所述方法包括：丰度矩阵提取、地物反射率的时间变化特征提取、高时空分辨率融合影像生成和融合影像评价体系等四个部分。

Description

一种多源遥感数据时空融合方法

技术领域

一种多源遥感数据时空融合方法，属于数字图像处理技术领域，特别涉及遥感图像及遥感图像分类技术。

背景技术

多源遥感数据融合是利用不同空间分辨率、不同时间分辨率和不同光谱分辨率的光学遥感数据，甚至是光学遥感数据与雷达数据所各自具有的不同数据特点，融合出能够综合多源数据特征的融合影像的方法。它对提高遥感影像利用效率和提高遥感应用效果都具有重要的意义。

多源遥感数据时空融合是多源遥感数据融合领域的前沿研究内容之一。目前多源遥感数据融合主要是全色高空间分辨率影像与多光谱影像间的融合，融合的目的只要为提高多光谱影像的空间分辨率，融合的方法主要有代数法和基于分量替换的方法等。

代数法是通过对待融合影像进行代数运算来实现融合，如利用乘法运算进行SAR影像和光学影像的融合，利用相关系数加权方法进行全色影像和多光谱影像的融合等。

基于分量替换的融合方法是将影像先进行某种变换，然后利用高空间分辨率影像替换变换空间中的一个分量，在进行逆变换来实现融合。分量替换的融合方法主要包括HIS变换融合和小波变换融合等。

HIS变换融合(参见文献：赵英时等.遥感应用分析原理与方法.北京：科学出版社，2003)是将图像彩色空间RGB空间变换到HIS空间，再用全色波段替换HIS空间中的强度分量I分量，再进行HIS逆变换实现影像的融合。其目的主要为提高多光谱影像的空间分辨率。

小波变换融合(参见文献：赵英时等.遥感应用分析原理与方法.北京：科学出版社，2003)是利用小波变换具有变焦性、信息保持性和小波基选择的灵活性等优点，将图像分解为一些具有不同空间分辨率、频率特性和方向持性的子信号(图像)，其中它的分频特征，相当于高、低双频滤波器，能够将一信号分解为低频信息(图像)和高频细节/纹理信息(图像)，同时又不失原信号所包含的信息，将高空间分辨率影像的小波变换后的高频部分代替多光谱影像的小波变换高频部分，再进行小波逆变换来实现融合。小波变换融合可以用于以非线性的对数映射方式融合不同类型的图像数据，使融合后的图像既保留原高分辨率迢感影像的结构信息，又融合多光谱影像丰富的光谱信息，提高影像的解译能力、分类精度。

基于分量替换的融合方法还包括基于G-S变换的影像融合、基于主成分分析的影像融合和基于高通滤波影像融合等。

发明内容

本发明提供一种多源遥感时间的时空融合方法，用以实现利用时序低空间分辨率影像和高空间分辨率影像，融合生成高时空分辨率影像，解决目前由于遥感卫星数据获取能力不足和云等天气条件的影响造成的遥感数据缺乏问题。

本发明技术方案如下：

一种多源遥感时间的时空融合方法，其特征在于包括一下具体步骤：

步骤一、丰度矩阵提取。选取2期的中空间分辨率影像，计算其NDVI，选定遥感影像分类方法，对时序NDVI影像进行逐窗口分类，根据低空间分辨率影像空间分辨率计算低空间分辨率影像的丰度矩阵；其中窗口大小采用融合效果最佳的窗口，本发明中为40×40MODIS像元大小。

步骤二、地物反射率的时间变化特征曲线提取。利用时序低空间分辨率影像和丰度矩阵，采用最小二乘法，逐窗口解算各个类别的时序平均反射率，作为地物的时间变化特征曲线。

步骤三、融合影像生成。利用地物反射率的时间变化特征，结合前期高空间分辨率影像的纹理特征，采用反射率融合模型，融合出空间分辨率与中分辨率影像相同，时间分辨率与高时间分辨率影像相同的高时空分辨率融合影像；其中，反射率融合模型是基于类内像元反射率的时间变化与类均值反射率时间变化相似假设的，其重点在于确定同类像元对目标像元反射率的贡献，即权值。

步骤四、融合效果精度评价。利用真实中分辨率影像，基于相似性原理，采用相关系数、方差、平均绝对误差、RMSE和偏差等参数对融合影像，对融合影像进行精度评价。

本发明与现有技术相比所具有的优点：本发明是一种基于像元反射率时间变化特征的融合方法；与传统方法相比，本发明不仅可以提高融合影像的空间分辨率，还可以提高影像的时间分辨率，从而可以解决目前由于卫星传感器数据获取能力不足、以及光学卫星容易受到云的影响而导致的遥感数据缺失问题。

附图说明

图1：本发明的一种多源遥感时间的时空融合方法流程图；

图2：转换系数v计算图；

图3：融合结果图；

图4：融合结果精度评价图；

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图及具体实施方式详细介绍本发明。

本发明一种多源遥感时间的时空融合方法，该方法主要包括以下几个步骤：

1、丰度矩阵提取；

2、地物反射率的时间变化特征曲线提取；

3、融合影像生成；

4、融合效果精度评价；

本发明的具体实现流程如图1所示，各部分具体实施细节如下：

1、丰度矩阵提取

根据线性混合模型，每一个低空间分辨率像元的反射率R(i，t)可以看做各组分平均反射率

的线性组合：

$R(i,t)=\underset{c=0}{\overset{k}{\mathrm{\Σ}}}{f}_c(i,c)\×\stackrel{\‾}{r}(c,t)+\mathrm{\ξ}(i,t)---\left(1\right)$

$\underset{c=0}{\overset{k}{\mathrm{\&Sigma;}}}{f}_c(i,c)<1;$ f_c(i，c)≥0

其中：f_c(i，c)是组分C在t时刻的丰度矩阵；ξ(i，t)是残差；k是组分总数量。在可以从中高分辨率分类图中获得丰度矩阵f_c(i，c)的情况下，通过等式(1)，可以利用最小二乘法解求t时刻组分平均反射率

要进行低空间分辨率影像与高空间分辨率影像融合，丰度矩阵f_c(i，c)是一个必须已知的量。传统方法中，丰度矩阵f_c(i，c)都是从中高分辨率影像分类图中获取的，并假设其不随时间发生变化。然而，考虑到历史时期的一景中空间分辨率影像是较容易获取的，若能够获取中空间分辨率影像的时间变化规律，则可从前期的中空间分辨率影像推算出新时期的中空间分辨率影像。

为获取中空间分辨率影像的时间变化规律，本文用中空间分辨率影像的变化分类图，因此在本文中，类别表示从期初(T_l)到期末(T_n)具有相同变化的像元集合，并假设类别不随时间变化。

首先，我们计算两期影像的NDVI，选定分类方法，如ISODATA等，设定类别初始值，对时序NDVI影像在窗口内逐窗口进行分类，获得分类图，再设置与低空间分辨率影像空间分辨率相同的格网，统计提取出丰度矩阵f_c(i，c)。

2、地物反射率的时间变化特征曲线提取

输入丰度矩阵f_c(i，c)和从期初(T_l)到期末(T_n)的时序低空间分辨率反射率影像，利用最小二乘法，解算等式(2)，获得时序的中空间分辨率像元类别平均反射率，将类别评价反射率作为地物类别的时间变化特征曲线。

3、融合影像生成

通过等式(1)只能够获得中空间分辨率像元的类别平均反射率，要获得新时期的中空间分辨率反射率影像，需要从类别平均反射率进一步获得每个中空间分辨率像元的反射率。为此，本发明提出了一种从类别平均反射率获取像元反射率的模型，其假设同类像元的反射率时间变化一致，即：

r(c，t_i，k)-r(c，t_j，k)＝r(c，t_i，l)-r(c，t_j，l)，(2)

其中：r(c，t_i，k)，r(c，t_j，k)，r(c，t_i，l)和r(c，t_j，l)是分别为类别C中像元k和像元l在时间t_i和t_j的反射率。

由于类别平均反射率是属于该类别像元的反射率平均值：

$\stackrel{\&OverBar;}{r}(c,t_i)=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}r(c,t_i,k),---\left(3\right)$

其中：n是属于类别C的像元个数。因此，从t_i时刻到t_j时刻，类别C的平均反射率变化

为：

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{r}(c,t_j-t_i)=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;r}(c,t_j-t_i,k),---\left(4\right)$

其中Δr(c，t_j-t_i，k)是类别C中像元k从t_i时刻到t_j时刻的反射率变化。而

和Δr(c，t_j-t_i，k)可由下式计算：

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{r}(c,t_j-t_i)=\stackrel{\&OverBar;}{r}(c,t_j)-\stackrel{\&OverBar;}{r}(c,t_i),---\left(5\right)$

Δr(c，t_j-t_i，k)＝r(c，t_j，k)-r(c，t_i，k)，(6)

从而

Δr(c，t_j-t_i，k)＝Δr(c，t_j-t_i，l).        (7)

将等式(5)(6)(7)代入等式(4)得：

$\stackrel{\&OverBar;}{r}(c,t_j)-\stackrel{\&OverBar;}{r}(c,t_i)=r(c,t_j)-r(c,t_i).---\left(8\right)$

由于

和

可以通过最小二乘法，解算等式(1)获得，r(c，t_i，k)可以从前期的中空间分辨率影像获得，因此可以利用等式(8)计算新时期的中空间分辨率影像反射率r(c，t_j)。

由于几何配准误差的存在，直接采用等式(8)进行中低空间分辨率影像的融合会导致较大的误差。大量研究表明，由于同类像元的反射率时间变化特征相似，选择临近窗口S范围内的同类像元，并根据其与目标像元的光谱相似性和距离等赋予权重W，可以有效降低几何配准误差的影响，用公式表示为：

$r(c,t_j)=r(c,t_i)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{W}_i\&times;{\mathrm{\&upsi;}}_{c,i}\&times;(R_j-R_i)---\left(9\right)$

其中：N是窗口S中属于组份c的像元个数；W_i是窗口S属于组份c的每个像元的权重。

权重W决定了窗口内同类别像元对目标像元反射率解算的贡献。它主要由类别像元与目标像元之间的光谱相似性和距离决定。光谱相似性越高、距离越近的同类别像元具有更高的权重。同时由于像元反射率的时间变化特征对融合具有重要的意义，在计算光谱相似性时，本研究计算的是时序反射率影像间的光谱相似性，因而获得光谱相似性能够同时反应反射率的时间变化特征与光谱相似性，光谱相似性的计算公式为：

$S_i=\frac{E\left[(r_i-E\left(r_i\right))(R_i-E\left(R_i\right))\right]}{\sqrt{D\left(r_i\right)}\&CenterDot;\sqrt{D\left(R_i\right)}}---\left(10\right)$

r_i＝{r(x_i，y_i，t_m，B₁)，…，r(x_i，y_i，t_m，B_n)，r(x_i，y_i，t_n，B₁)，…，r(x_i，y_i，t_n，B_n)}

R_i＝{R(x_i，y_i，t_m，B₁)，…，R(x_i，y_i，t_m，B_n)，R(x_i，y_i，t_n，B₁)，…，R(x_i，y_i，t_n，B_n)}

其中：S_i是位置(x_i，y_i)处的像元与目标像元的光谱相似性；r_i是高分辨率像元的时序反射率矢量；R_i是低分辨率像元的时序反射率矢量。

距离权重的计算公式为：

$d_i=1+\sqrt{{(x_o-x_i)}^2+{(y_o-y_i)}^2}/(w/2)---\left(11\right)$

其中：d_i为位置(x_i，y_i)处的像元与目标像元(x₀，y₀)的距离；w是窗口大小。

将光谱相似性权重和距离权重进行组合获得组合权重D_i：

D_i＝(1+S_i)×d_i        (12)

对组合权重进行归一化，获得最终的权重Wi：

$W_i=(1+D_i)/\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(1+D_i)---\left(13\right)$

对时序低分辨率像元反射率和高分辨率像元反射率进行回归分析可以获得转换系数v。由于几何配准误差的存在，严格根据几何位置进行转换系数v的回归会存在较大的误差。同时，在临近区域内，由于环境相似，相临像元间的转换系数v应相似，因此本发明是逐低空间分辨率像元进行转换系数v计算的，并认为在整个低空间分辨率像元内的v值相同(图2)。

生成的融合影像见图3，从图中我们可以看出，融合影像和真实影像非常相似。

4、融合效果精度评价

基于相似性原理，选择与融合影像同一天的真实高空间分辨率影像，利用相关性分析方法进行融合效果精度评价(图4)。由于相关系数表征了两者之间的相关程度，选用相关系数(r)、方差(variance)、平均绝对误差(MAD)、RMSE和偏差(bias)等参数对融合效果进行定量精度评价。定量评价结果见下表：

表1融合效果精度评价

Claims (3)

1.一种多源遥感数据时空融合方法，其特征在于包括以下步骤：

(1)选取2期的中空间分辨率影像，计算其NDVI，选定遥感影像分类方法，对时序NDVI影像进行逐窗口分类，根据低空间分辨率影像空间分辨率计算低空间分辨率影像的丰度矩阵；

(2)利用时序低空间分辨率影像和丰度矩阵，采用最小二乘法，逐窗口解算各个类别的时序平均反射率，作为地物的时间变化特征曲线；

(3)利用地物反射率的时间变化特征，结合前期高空间分辨率影像的纹理特征，采用反射率融合模型，逐窗口融合出空间分辨率与中分辨率影像相同，时间分辨率与高时间分辨率影像相同的高时空分辨率融合影像；

(4)利用真实中分辨率影像，采用相关系数、方差、平均绝对误差、RMSE和偏差等参数对融合影像，对融合影像进行精度评价。

2.根据权利要求1所述步骤(3)高时空分辨率融合影像生成，其特征在于1)反射率融合模型，该模型是根据时序中低空间分辨率影像间的关系模型，利用前期中分辨率影像反射率和地物时序反射率特征，计算出融合高空间分辨率影像；其中关键在于计算过程中，同窗口、同类别像元的权值的确定，本发明在确定同类像元权重时不仅考虑像元之间的光谱相似性和距离，还考虑像元间时间变化特征曲线的相似性。2)逐窗口运算，通过试算，选取融合效果最佳的窗口作为融合窗口，大大降低配准误差对融合效果的影响。

3.根据权利要求1所述步骤(4)融合影像精度评价，其特征在于区别于传统融合影像精度评价采用熵、清晰度和逼真度等指标，本发明采用基于相似性原理进行融合影像和真实影像的融合效果评价，采用相关系数、方差、平均绝对误差、RMSE和偏差等参数对融合影像，对融合影像进行精度评价。

Images