CN103576132A - 一种遥感影像的处理方法及*** - Google Patents

Abstract

本发明公开了一种遥感影像的处理方法，根据空间分辨率不同的两组影像求得某时刻的所求影像A0。首先根据空间分辨率不同的两组影像计算目标区域的反射率变化量在这两组影像之间的转换系数，其次根据所述转换系数将低空间分辨率影像中的变化量转换成高空间分辨率影像中的变化量，最后由此根据已知的高空间分辨率影像得到所求影像A0的目标像元的反射率。还相应的提供了遥感影像的处理***。本发明能生成将细节保持的更好的影像。

Description

一种遥感影像的处理方法及***

技术领域

本发明涉及遥感影像处理的技术领域，特别是指通过使用多于一幅遥感影像生成新的遥感影像的处理方法及***。

背景技术

综合应用多传感器遥感数据在空间、时间、光谱分辨率上的不同优势的技术，即多源遥感数据融合技术从20世纪80年代开始发展，目的是通过对多种传感器数据的处理，提高遥感数据的空间和光谱分辨率。目前较为成熟的融合方法有：IHS变换、主成分变换、小波变换以及Brovey变换、高通滤波等，但是这些方法主要用于融合同一或相近时刻的高空间分辨率全色波段影像和较低空间分辨率的多光谱波段影像，其结果图像为同一时刻的高空间分辨率多光谱影像。由于高空间分辨率的影像并不能和所有低空间分辨率影像同时获取，这些传统的多源数据融合方法并不能同时提高遥感数据的空间和时间分辨率。即，无法生成高时空分辨率的反射率数据。为此，Gao等人在2006年提出了STARFM方法，该方法能够得到较为准确的高时空分辨率反射率影像，但其生成的影像的细节保持的不好。

发明内容

有鉴于此，本发明的主要目的在于提供一种遥感影像的处理方法，以实现能生成将细节保持的更好的影像。

为解决上述技术问题，本发明提供了一种遥感影像的处理方法，设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域，获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像，按以下步骤分别计算各目标像元的反射率以生成A0影像：a，根据所选的影像计算所述目标区域在两组影像中的反射率变化量之间的转换系数山，从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像：根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量：c，根据所述转换系数将所述各反射率的变化量转换为在第一组影像中的各反射率变化量：和d，从第一组影像中选择与所述其它时刻对应的影像，根据转换得到的各反射率变化量和目标区域在所选影像中的反射率计算所述目标区域的反射率。

由至少两对时刻相对应的影像就可以得到时刻对应的两组影像中地表区域反射率的变化信息之间的关系，即转换系数。根据所挑的影像经由步骤a得到目标区域的转换系数。因为变化信息之间的关系在较短时间范围内可以认为是稳定的。所以只要所选的各影像的时刻和影像A0时刻相差不大，则可以根据转换系数经由步骤c将步骤b所得的目标区域在第二组影像中的反射率变化量(简称第二变化量)转换为目标区域在第一组影像中的反射率变化量(简称第一变化量)。由此，由步骤d所得的反射率是基于更精细的第一组影像中的反射率和转换后的变化量生成的。这样生成的A0影像就很好的保持了第一组影像中的细节，避免所生成的影像在视觉上有被平滑的效果。

优选的是，所述步骤a包括：a1，将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元：a2，从所述中心像元各自所在影像中分别筛选出至少一个与所述中心像元属于同类地物的相似像元：和a3，确定各相似像元所对应的地表区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述转换系数。

混合像元反射率的变化量为各端元反射率的变化量的综合反映。所以端元可以提供比混合像元更为准确的变化信息。又因为作为端元的相似像元与中心像元属于同类地物，所以，相似像元能够提供最为准确的变化信息，由此可以得到更为准确的转换系数。其中，中心像元也算是相似像元之一，在未能筛选出其它相似像元时，中心像元是唯一的相似像元。尤其是对于目标区域中存在破碎斑块和小地物的地表区域，其影像在第二组影像中为混合像元。其提供的变化信息(第二变化量)被由各相似像元得到的转换系数转换修正成与目标像元属于同类地物的端元的变化信息(第一变化量)。由此，目标像元的反射率准确度更高了，尤其是对破碎斑块和小地物具有更好的融合效果。

优选的是，所述步骤a3之前还包括：a21，根据所选的影像计算各相似像元的权重：所述步骤b包括：b1，根据B0影像和所述第二组的其他时刻的影像计算各所述相似区域在第二组影像中的反射率变化量山2，根据得到的所述权重将步骤b1得到的各反射率变化量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。

目标区域的第二变化量不是直接由目标区域在B0影像和所述第二组的其他时刻的影像的反射率求差得到，而是由各相似像元按照各自权重共同来提供。由此修正了B0影像和第二组的其他时刻的影像中目标区域对应的混合像元所提供的变化信息。从而得到更准确的第二变化量。

优选的是，所述步骤a3中所述转换系数的计算是将各所述相似区域在所选影像中的反射率以在第二组影像中的反射率的变化量为自变量，在第一组影像中的反射率的变化量为因变量进行回归分析。

利用回归分析可以消除因噪声等原因引起的误差，从而得到更稳健的变换系数，减少可能的误差。

优选的是，所述步骤a3之前还包括：a22，根据所选的影像计算各相似像元的权重：所述回归分析是根据各相似像元的权重进行的加权回归分析。

各相似像元距离其中心像元的空间距离，以及其的光谱特征和中心像元的光谱特征的一致性等原因都导致了各相似像元提供变化信息时，存在贡献大小的差异。在回归分析时，利用各相似像元的权重进行加权回归分析，可以体现出这些差异。进而得到更准确的转换系数，以进一步提高目标像元的反射率准确度。

优选的是，所述各相似像元的权重的计算包括以下步骤：a221，根据所选影像计算各相似像元所有波段光谱向量和第二组影像中各所述相似区域的像元所有波段光谱向量的一致程度：a222，根据所述一致程度计算各相似像元的权重。

所述第一组影像的像元为高空间分辨率像元，所述第二组影像的像元为低空间分辨率像元，高空间分辨率像元对应地表区域所属地物在对应的低空间分辨率像元对应地表区域的各种地物中所占比例即为高空间分辨率像元的纯度。所以相似像元的纯度越大，也就意味着其所提供的变化信息越准确。从而通过步骤a221和步骤a222可以提高能提供更准确变化信息的纯度大的相似像元的权重。而纯度大小的计算则是由所有波段光谱向量的一致程度计算得到。

优选的是，所述步骤a2还包括：将各筛选出的相似像元所对应的地表区域取交集，以取交集后的地表区域所对应的像元为相似像元。

有的地物随时间会发生变化，其光谱特征也相应发生变化，只利用一组时刻对应的影像筛选出的相似像元可能会因此筛选出已和中心像元不属于同一地物的相似像元。所以取交集后可以筛选出的各时刻的均和中心像元都是同一地物的相似像元。从而得到更可靠的目标像元。

本发明还相应的提供了一种遥感影像的处理***，根据所述两组影像求得高空间分辨率的A0影像，所述A0影像不属于第一组影像，并由至少一个待计算的像元组成，所述处理***包括：影像选择模块101，用于获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像：目标像元设定模块102，用于设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域：目标像元反射率计算模块104，用于计算目标像元的反射率：和影像生成模块108，用于根据目标像元反射率计算模块104计算得到的各目标像元的反射率生成A0影像：所述目标像元反射率计算模块104包括以下子模块：转换系数模块103，根据所选的影像计算所述目标区域在两组影像中的反射率变化量之间的转换系数：第二变化量计算模块105，用于从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量：变化量转换模块106，用于根据所述转换系数将所述各第二变化量计算模块105得到的各变化量转换为在第一组影像中的各反射率变化量：和反射率计算模块107，用于从第一组影像中选择与所述其它时刻对应的影像，根据转换得到的各反射率变化量和目标区域在所选影像中的反射率计算所述目标区域的反射率。

优选的是，所述转换系数模块103包括：中心像元设定模块1031，用于将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元：相似像元筛选模块1032，用于从所述中心像元各自所在影像中分别筛选出至少一个与所述中心像元属于同类地物的相似像元：和转换系数计算模块1033，用于确定各相似像元所对应的地表区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述转换系数。

优选的是，所述目标像元反射率计算模块104还包括：权重模块，用于根据所选的影像计算各相似像元的权重：所述第二变化量计算模块105包括：相似像元第二变化量计算模块1051，用于从第二组影像中选择与所述某时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算各所述相似区域在第二组影像中的反射率变化量：第二变化量加权平均模块1052，根据得到的所述权重将相似像元第二变化量计算模块1051得到的各反射率变化量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。

附图说明

图1为遥感影像的处理方法的实施例流程图：

图2为筛选相似像元的流程图：

图3为取交集筛选相似像元示意图：

图4为本实施例中计算相似像元的权重的流程图：

图5为得到转换系数的流程图：

图6为转换系数的计算示意图：

图7为计算预测时刻的目标像元的高空间反射率的流程图：

图8为处理***的结构图：

图9为转换系数模块103结构图：

图10为第二变化量计算模块105结构图：

图11为一个实施例中权重模块109结构图：

图12为另一个实施例中权重模块109结构图：

图13为相似像元筛选模块1032结构图：

图14为调整模块1034结构图

图15为模拟物候变化的情况：

图16为模拟线状地物的情况：

图17为MDCM和STARFM对图16(b)中线状地物的反射率的预测结果：

图18为模拟小面积地物的情况：

图19为MDCM和STARFM对图18(b)中各圆形小地物的预测误差：

具体实施方式

在对实施方式进行说明之前，先对本发明的理论基础进行说明。

对于同一区域来自不同传感器的遥感数据经过辐射定标、几何配准和大气纠正后，这些数据之间具有一定的可比性和相关性。但是，由于不同传感器在波段宽度、光谱响应函数、获取时刻的大气状况等方面存在差异，使得这些多源数据间存在一定的***偏差。在对***误差进行纠正的同时，如何利用它们之间的相关性实现多源数据的融合构成了本发明的核心。以下分两种情况介绍本发明理论基础。

对于均一地表情况：假设低空间分辨率、高时间分辨率(略写为低空间分辨率)和高空间分辨率、低时间分辨率(略写为高空间分辨率)传感器有类似的光谱波段设置。当低空间分辨率影像重采样为与高空间分辨率影像相同的空间分辨率(即相同的像元大小)和坐标***时，对于同一且均一的地物而言，低空间分辨率像元内仅为一种地物类型纯像元，此时在波段B上闻空间分辨率反射率与低空间分辨率反射率之间的关系即为闻、低空间分辨率影像上反射率之间的相对辐射定标关系，该关系仅由两种传感器特性差异(波段宽度及光谱响应函数)和成像时刻大气状况差异所决定，一般可视为以下线性关系：

F(Xi，Yj，tk，B)＝aXC(Xi，jptk，B)+b1

其中，F、C分别代表高、低空间分辨率的像元反射率，Xi，，yj)是像元位置，B为光谱波段，tk代表了影像获取时刻(以日为单位)，a、b分别代表高空间分辨率和低空间分辨率影像相对辐射定标的增益和偏离值。对于同一且均一的地物而言，1式所示关系在任何影像获取时刻都成立。

如果已经获取了、时刻波段B的高、低空间分辨率的影像和tp时刻同一波段的低空间分辨率影像，那么tp时刻该波段高空间分辨率的像元反射率可由式2来表示：F(Xi，yj，tp，B)＝aXC(Xi，yj，tp，B)+b2

而、时刻高、低空间分辨率影像的像元反射率关系为：

F(Xi，yj，t0，B)＝aXC(Xi，yj，t0，B)+b3

由式2和式3相减并移项整理可得到式4：

F(Xi，yj，tp，B)＝F(Xi，jPt0，B)+aX(C(Xi，y』，tp，B)-C(Xi，jPt0，B))4

上式可以理解为tp时刻高空间分辨率的像元反射率由、时刻高空间分辨率的反射率加上tp相对、时刻的反射率变化量。式中a代表高空间分辨率和低空间分辨率影像相对辐射定标的增益，也可理解为高空间分辨率和低空间分辨率像元反射率之间的转换系数，该系数由两种传感器特性差异和成像时刻大气状况差异所决定。如果两种传感器成像时刻相近，可视为大气状况相同，则该系数仅由两种传感器特性差异所决定，且不随时间发生变化。

如果获取任意两个时刻I、tn的两对高、低空间分辨率影像，那么通过对这两个时刻的高、低空间分辨率图象上的各波段纯像元反射率进行线性回归即可求出各波段的转换系数a，通过该系数和式4即可由一系列时间连续的低空间分辨率影像生成一系列时间连续的高空间分辨率的反射率影像。

对于非均一地表情况：由于传感器空间分辨率的限制以及地物的复杂多样性，混合像元普遍存在于低空间分辨率影像上，对地表地物类型多样、分布杂乱的区域尤其如此。此时，低空间分辨率像元与其对应的高空间分辨率像元的反射率并不是对同一种地物的反应，二者之间的差异不仅仅是辐射定标的差异，此时基于均一地表纯像元的式4并不适合非均一地表混合像元的情况。

混合像元的光谱特征是其内部各类地物光谱端元光谱的综合反映，所以两个时刻tm、tn之间其反射率的变化为各类端元光谱变化的综合反映。目前混合像元光谱混合模型最常用的是线性光谱混合模型，即利用一个线性关系来表达混合像元和各端元光谱之间的关系。假设低空间分辨率影像上tm、tn时刻的某波段混合像元反射率分别为Ym、Yn，根据线性光谱混合模型，Ym、Yn可由下式表达：

M

Ym＝YJfiXim+S

/＝1

M

Yn＝Yuf＜Xin+s

′＝1  5

其中，M为该混合像元包含的端元数目，Xim和Xin分别表示tm、tn时刻第i种端元的光谱，e为误差项，且考虑e在tm、tn时刻不变。由式5可得到变化量Yn-Ym：其中下式中的fi是第i个端元的比例：

Yn-Ym＝2]/：(Xjn-Xim)

，.＝i6

假设tm、tn的时间间隔在一定范围内时，各端元反射率的变化可以近似为线性变化，则端元光谱Xin可由式7表示：其中下式中的ri是第i个端元反射率的线性变化率：Xin＝AXAt+Xim7

其中，At＝tn-tm，将式7代入式6可得到下式-r^A^：众8已知第

/＝1

k端元tm、tn时刻的反射率，则由式7可求出At＝A^Xkn～Xkm9将式9代入式8，可得：

rk

XknXkm_1

/＝1

由于假设tm、tn的时间间隔在一定范围内，各端元比例f·和各端元反射率的线性变化率r在tm、tn时刻之间可视为稳定，则式10的右边为一稳定变量，用v表示，其意义为k类地物端元反射率变化量与混合像元反射率变化量的比值变化比例系数。

对于非均一地表，低空间分辨率像元为混合像元，其对应的高空间分辨率影像上的各像元视为端元。根据式10，tm、tn时刻的高空间分辨率像元反射率F(Xi，Yi，tffi，B)和F(Xi，Yi，tn，B)与对应的低空间分辨率像元反射率C(Xi，Yi，tm，B)和C(Xi，Yi，tn，B)满足以下关系：

F(xi，yi，tn，B)-F(xi，yl，tm，B)

Cix^t^-C^y^t^B)′11

通过式11可知，tm、tn时刻的各波段闻空间分辨率反射率与对应的低空间分辨率反射率为线性关系，进行线性回归即可求出各波段的变化比例系数V。

如果在tm到tn内，已经获取了tQ时刻波段B的高、低空间分辨率的影像和％时刻同一波段的低空间分辨率影像。根据以上的假设，tp相比、时刻该波段高空间分辨率的像元与低空间分辨率像元的反射率变化量之间的比例系数仍为v则：F{xi，y，，tB)-F(xt，t0，B)

C(xi，yi，tp，B)-C(xi，yi，t0，B)^

通过式12即可求解出tp时刻的高空间分辨率反射率：

F(Xi，yj，tp，B)＝F(Xi，y』，t0，B)+vX(C(Xi，y』，tp，B)-C(Xi，y』，t0，B))13

虽然公式13与公式4的形式相同，但其意义和适用情况不同。式4表示的是纯像元不同分辨率反射率之间相对定标的关系，在任何时刻都成立，其计算的结果最为准确。而式13表示的是混合像元中不同分辨率反射率变化量之间的比例关系，根据假设只在之间或其附近时刻才成立，而且是基于短时间范围内各地物反射率的变化为线性的假设，所以其计算的结果是近似解。

基于上述的理论基础，当低空间分辨率像元为纯像元时，可以利用式4求解出tp时刻的高空间分辨率反射率，而当低空间分辨率像元为混合像元时利用式13求解。但纯像元计算的结果较混合像元准确可靠。为尽量利用纯像元的信息，以及减小影像受到云、大气等污染给计算结果带来的不确定性，本实施例中采取窗口内求解的办法。即，以像元xw/2，，yw/2)为中心，取宽度为w的窗口，先筛选出与中心像元属同类地物的像元相似像元，然后利用式14计算中心像元tp时刻的高空间分辨率反射率：

N

Hk’/2，yw/2，5)＝f(xw/2＞yw/2A^B)+xVx(c(uJp’B、-c(u，h，B))

＜＝114

其中，N为窗口内与中心像元属于同类地物的相似像元数目，W为各相似像元的权重，V为高、低空间分辨率的转换系数。式14的意义是tp时刻的反射率由已知时刻的反射率加上tp时刻相对已知h时刻的变化量决定，而该变化量由窗口内的相似像元共同预测。

在本实施例之前，通常需要先对获取的高、低空间分辨率影像进行预处理，使得多源数据之间具有相似的波段设置，相同的像元大小和图幅尺寸以及相同的坐标***。采用常用的软件和方法即可实现。下文将结合附图对遥感影像的处理方法的实施方式进行详细说明。

本实施方式中遥感影像为对同一地表摄像而得到的两组影像，其中，第一组影像具有高空间分辨率，第二组影像具有低空间分辨率，根据所述两组影像求得tp时刻的高空间分辨率的A影像，所述A影像由至少一个待计算的像元组成。如图1所示，遥感影像的处理方法包括以下步骤：步骤100，从所述遥感影像中挑选至少两对时刻对应的影像。本实施例中，选用的是已知的tm、tn时刻的两对影像。从第二组影像中挑选其时刻与tp时刻所对应影像：步骤101，设定一个待计算的像元为目标像元，确定目标像元的应的地表区域为目标区域：步骤102，以目标区域对应的像元为中心像元，筛选出已知时刻的高空间分辨率影像中所述中心像元的窗口内的同类地物像元得到相似像元，其中，中心像元也算是相似像元之一，在未能筛选出其它相似像元时，中心像元是唯一的相似像元：步骤104，计算所述相似像元的权重。本实施例中，根据tm、、时刻的两幅低空间分辨率影像，计算出步骤102中筛选出的相似像元的权重：步骤106，根据所述权重采用加权算法对已知时刻的高、低空间分辨率反射率进行线性回归计算得到转换系数：步骤108，计算预测时刻的目标像元的高空间分辨率反射率，本实施例中，根据转换系数和预测时刻tp的低空间分辨率影像计算已知时刻高空间分辨率影像的时间权重，并根据已知时刻高空间分辨率影像按照所述时间权重加权计算得到预测时刻的目标像元的高空间分辨率反射率：步骤110，逐次计算所有目标像元并最终生成A影像。下文将结合附图对上述各个步骤进行详细说明。

图2为本实施例中筛选相似像元的流程图，如图所示包括：

步骤1022，以目标区域对应的像元为中心像元，筛选出各已知时刻的高空间分辨率影像中所述中心像元的窗口内的与目标像元属于同类地物的像元得到相似像元。在目标像元(xw/2，，yw/2)为中心的窗口内，与要计算的中心像元属于同类地物的高空间分辨率像元才能提供比较正确可靠的反射率变化信息。可以使用以下两种方法筛选出相似像元：(a)对高空间分辨率影像进行非监督分类，和中心像元类型相同的像元为相似像元：(b)使用阈值判断，以窗口内其它像元反射率与中心像元反射率之差来判断是否属于同类地物。这两种筛选相似像元的方法具有共同的意义，即是以光谱相似的像元为同一类地物。但是又有不同之处：阈值法是在窗口内寻找与中心像元的光谱相似的像元，中心像元作为了搜寻的中心，而且其筛选条件是局部适用的，随着中心像元的位置变化而变：而分类是在全影像进行，如果分类有误差，那么判断出的相似像元可能具有很大的错误。所以在本实施例中，使用阈值法来筛选相似像元。在本实施例中，判断窗口内某像元所有波段都满足式15的条件，则该像元被确定为中心像元的相似像元。

|F(Xi，tk，B)-F(xw/2，yw/2，tk，B)|(o(B)X2/m15

其中，o(B)为tk时刻高空间分辨率影像所有像元B波段的标准差，m为预估的地物类别数，比如影像覆盖的区域主要有植被、裸地、水体、岩石4类地物，则可将m值设为4。当然，也可以用所有波段中的部分波段进行筛选或者采用方法(a)进行筛选。

步骤1024，将筛选出的两组已知时刻的相似像元取交集。有的地物随时间会发生变化，其光谱特征也相应发生变化，只利用一个时相的影像筛选出的相似像元可能会出现错误。比如窗口内有裸地和农作物两种地物类型，如果中心像元是农作物，tm时刻农作物并未生长，其光谱特征与裸地相同，此时筛选出来的相似像元是裸地，如果tn时刻农作物已生长，则在、时刻的影像上能筛选出正确的农作物像元。鉴于此，本实施例中分别利用tm、tn时刻两幅高空间分辨率影像筛选相似像元，并将筛选结果取交集，以此提高筛选结果的准确性。如图3所示为取交集筛选相似像元示意图表示中心像元表示中心像元的相似像元。

图4为本实施例中计算相似像元权重的流程，如图所示包括：

步骤10402，计算相似像元的纯度大小。权重W决定了窗口内筛选出的各相似像元对计算目标像元反射率的贡献。在本实施例中权重W由各相似像元的纯度大小和离中心像元的空间距离远近共同决定。纯度越高、空间距离越近的相似像元权重越大，这是因为均一地表的纯像元能提供最准确的变化信息，而且认为空间上离中心像元越近的相似像元的变化情况和中心像元越一致。

高空间分辨率影像上的某像元代表了某种地物类型，而其对应位置低空间分辨率

影像的像元包含了该类地物，纯度即是指低空间分辨率像元中该种地物所占比重。如果在均一地表情况下低空间分辨率像元全由该种地物覆盖，则认为是纯像元。比如Landsat影像上某像元所覆盖的30mX30m的区域为小麦，其对应位置的M0DIS影像上像元所覆盖的500mX500m的区域也全为小麦，则认为该像元为小麦纯像元。由于不同的地物具有不同的光谱特征，所以可以通过比较高空间分辨率像元和其对应位置的低空间分辨率像元的光谱曲线的一致程度来判断纯度大小。如果高空间分辨率像元和其对应的低空间分辨率像元的光谱一致程度越高，则认为纯度越大。而高、低空间分辨率像元的光谱曲线之间的一致性可由相关系数来刻画，所以筛选出的各相似像元的纯度R可由式16计算：分cov(F，，C，)

^D(ct)16

Fi＝{F(Xi，Yi，tm，Bj)，…，F(Xi，tm，Bn)，F(x^tn，Bj)，…，F(Xi，tn，Bn)}

Ci＝{C(Xi，tm，Bj)，C(Xi，y”tm，Bn)，C(Xi，y”tn，BD，…，C(Xi，tn，Bn)}

其中向量FyCi分别表示第i个相似像元高、低空间分辨率的光谱向量，cov为求协方差，D为求方差。R值的范围为-1到1，R值越大表示纯度越高。之所以将tm、tn时刻的光谱曲线放在一起计算相关系数，是考虑到地物随物候的变化，光谱特征也会发生变化，只使用一个时相的光谱可能会出现错误。比如小麦在返青之前其光谱特征和裸地一致，如果某小麦高空间分辨率像元对应的低空间分辨率像元除了小麦外还包含裸地，则用小麦返青前的获取的影像计算会得到较高的纯度，显然是不正确的。当加入小麦返青后获取的影像进行判断时，就会得到更加正确的纯度。由上可以看出，只要获得至少两组已知时刻的高空间分辨率像元的所有波段光谱向量和在其对应位置同时的低空间分辨率像元的所有波段光谱向量，就可以根据其一致程度来计算纯度大小。

利用纯度大小作为计算权重的一个要素，较之目前采用同时刻高、低空间分辨率反射率之差表示的光谱特性作为计算权重的要素将更加有效的利用相似像元并且得到的结果也更加准确。

步骤10404，判断相似像元中是否有纯像元，有则进入步骤10406，没有则进入步骤10408。

步骤10406，设定纯像元的权重为1。在本实施例中，权重W的取值范围是0至IJ1，所有相似像元的权重之和为1。R＝1的像元为纯像元。为了尽量利用纯像元，将纯像元的权重设为1。即，中心像元的反射率变化信息完全由这些纯像元提供。在另一个实施例中，当相似像元中有P个纯像元存在即R＝1的像元时，规定这些纯像元的权重W为1/P，其他相似像元的权重为0。当然，纯像元的判断也可以是在计算权重过程中的其它步骤中实现。例如在步骤10402之前或者在步骤10408之后。甚至不判断是否有纯像元也可以。

步骤10408，计算相似像元距离中心像元的空间距离。相似像元Xi，，Yi)与窗口中心像元xw/2，，yw/2)的空间距离屯由式17计算：

dt＝1+yl(xw/2-X，.)2+(yw/2-兄)2/(w/2)17

上式包含了对空间距离的归一化，在窗口w内各相似像元的空间距离屯的取值范围为1到1+20°_5，值越大，该相似像元离中心像元越远。

步骤10410，计算相似像元的权重。接下来将纯度和空间距离这两个不相关的量联系在一起，来确定各相似像元的权重W。上文已经提到纯度越高、距离越近的像元能提供更准确的信息，即R越大，d越小，其权重W越大。本实施例中，通过公式18，将纯度和空间距尚结合为总距尚D：Dj＝(1-Rj)Xdj18

D越大的相似像元为目标像元反射率的计算贡献越小，再将D的倒数归一化得到各相似像元的权重W式19

，＝i19

当然，建立纯度、空间和权重的关系可以采用多种形式，只要能实现纯度越大，距中心像元越近的相似像元的权重越大即可。

图5为得到转换系数的流程，如图所示包括：

步骤1062，计算转换系数。由公式4和公式13可知，转换系数V包含了相对定标的增益a和高低空间分辨率变化量的比例系数V，二者都可由已经获取的tm、tn时刻高、低空间分辨率各波段反射率进行线性回归得到。为了减小噪声污染等对转换系数计算带来的不确定性，将所有相似像元tm、tn时刻高、低空间分辨率反射率进行线性回归，回归直线的斜率即是转换系数V。为了利用更可靠的相似像元的信息计算转换系数，采用加权最小二乘法进行线性回归，权重为各相似像元的W。如图6所示为转换系数的计算示意图。某窗口内有12个像元和中心像元属于同类地物包括中心像元本身，虚线框分别标注出、、tn时刻的反射率，可见从tm到tn时刻，高空间分辨率反射率和其对应的低空间分辨率反射率都有所增加，但是高空间分辨率反射率增加幅度较大，是低空间分辨率的6.448倍R＝0.915，P＜0.001。这说明该类地物该波段反射率在tm、tn时刻之间发生了较大的变化，比如植被生长导致的近红外反射率急剧增加。根据理论基础的论述，各相似像元都反映了低、高空间分辨率反射率之间的转换系数，但由于噪声等影响，通常不能只利用某一个相似像元来计算V，这样的结果带有很大的不确定性，而线性回归能利用所有相似像元的信息得到一个比较稳健的转换系数。又因为各相似像元存在贡献大小的差异，为了体现这些差异，在线性回归之中进一步还考虑各样本的权重，本实施例采用的是加权最小二乘法。当然，也可以采用其它的回归方法。

步骤1064，根据不确定性调整转换系数。由于几何配准、大气修正等数据处理过程都给反射率带来一定的误差，如果低空间分辨率像元tm、tn时刻之间的反射率变化太小，在误差范围内，则进行加权回归得到的斜率V具有很大的不确定性。这可能是以下两种情况：低空间分辨率像元内反射率发生变化的地物所占面积比例很小，在整个低空间分辨率像元尺度上体现不出这种变化：也可能是低空间分辨率像元内有的地物反射率增加，有的地物反射率减少，而增加和减少的量相等，使得在低空间分辨率像元尺度上体现不出变化来。为解决这个问题，在本实施例中，在所有相似像元的低空间分辨率反射率变化均值小于至少两个已知时刻反射率的变化量的不确定性时将转换系数设为1。假设反射率各波段的不确定性u为该波段最大值的1％，而且这种不确定性在各影像各波段间是互相独立的，那么两个已知时刻反射率的变化量的不确定性为：U＝如J+Utn2

其中，utD1是utn分别是、、、时刻反射率的不确定性。如果所有相似像元的低空间分辨率反射率变化均值小于式20计算出的不确定性，则只能将第二变化量平均分配到其内部的每一个高空间分辨率像元上，即变化量的转换系数V为1。当然，已知时刻反射率的变化量的不确定性也可以通过其它统计方法得到。对式20的不同表达方式，以及其它的变换的计算方式都应包含在本发明的保护范围之内。这样就避免了两个时刻间反射率的变化太小时，因为反射率自身带有一定的不确定性，也就是误差，而导致计算出不真实的转换系数V。

步骤1066，修正奇异的转换系数。由于影像自身噪声的存在，导致绝对值较大的转换系数出现。为修正个别较奇异的转换系数，对全幅影像所有像元的转换系数进行统计，在本实施例中将均值上下2倍标准差以外的转换系数V设定为1。当然，也可以根据需要将其它预设倍数的标准差以外的转换系数设定为1。

图7为计算预测时刻的目标像元的高空间反射率的流程，如图所示包括：步骤10802，根据预测时刻的低空间分辨率影像计算各相似像元相对已知时刻的第二变化量：步骤10804，根据所述转换系数将所述第二变化量转换为第一变化量：步骤10806，将各相似像元的第一变化量按所述权重加权平均后得到目标像元的第一变化量：和，步骤10808，将所述目标像元的第一变化量与已知时刻的高空间分辨率反射率相加得到基于不同已知时刻的预测时刻目标像元的高空间分辨率反射率。通过上述4个步骤完成了式14的计算。以两组分别为tm、tn的已知时刻所获取的时刻之间或附近的任意tp时刻的低空间分辨率影像，计算出tp时刻的高空间分辨率反射率为例。根据式14，由已知的I、tn时刻高、低空间分辨率影像即可计算出tp时刻的高空间分辨率反射率，此时乜、、时刻的高、低空间分辨率反射率都已知，所以可以分别由tm、tn时刻的高空间分辨率反射率计算出tp时刻的高空间分辨率反射率Fxw/2，yw/2，tp，B，将这两个计算结果分别记作Fmxw/2，yw/2，tp，B、Fnxw/2，yw/2，tp，B。当然，只要已知转换系数、预测时刻的低空间分辨率影像和各相似单元的权重就可以在已知时刻的高空间分辨率影像的基础上计算得到预测时刻的目标像元的高空间分辨率反射率。对式14的不同表达方式，以及其它的变换的计算方式都应包含在本发明的保护范围之内。

步骤10810，加权计算得到预测时刻的高空间分辨率反射率。为充分利用tm、tn时刻的高空间分辨率反射率的信息，将两个计算结果加权求和得到最后的预测结果。假设tp时刻的低空间分辨率反射率相对已知时刻变化越小，则其计算的结果越准确可靠。这是因为如果低空间分辨率影像上没有体现出反射率的变化，则认为高空间分辨率影像也不会有反射率的变化。基于以上假设，两个计算结果的时间权重T由式21计算。

当然，也可以采用其它方法来计算基于不同已知时刻计算得到的高空间分辨率反射率的权重。例如，如式23，用预测时刻、和已知时刻tk低空间分辨率反射率之差来计算权重大小。

TiJk＝|C(Xi，yj，t0)-C(Xi，yj，tk)|23

至此得到了最终的预测时刻的目标像元的高空间分辨率反射率。仅是单独应用转换系数配合步骤108就能改善对小地物和线状地物的计算精度。而相似像元的筛选以及权重计算的改进则有助于整体精度的提高。

为了更好的对本实施方式更进一步说明和验证，下面对本实施方式分别进行模拟数据和真实数据的检验。

模拟数据检验。为了更好的检验方法的精确性和可靠性下面将本实施例运用于简单的模拟数据。并且为了便于和GA0的STARFM方法进行比较，利用GA0文中的模拟数据进行检验。在下文中本实施方式简称为MDCM。

线状地物的情况。在土地覆盖或土地利用类型中，线状地物相当普遍，比如道路、桥梁、河流等，而人眼对线状地物的敏感程度较高，生成的高分辨率影像中线状地物的正确与否直接影响影像的视觉效果。如图16所示为模拟线状地物的情况。图16。(a)-(c)中的圆形区域的反射率固定为0.05，直线的反射率固定为0.5，周围区域的反射率从0.1变为0.2再到0.4；(4)-(f)分别由(a)-(c)聚合而成的低空间分辨率影像：(g)、(h)为MDCM和STARFM对(b)中线状地物反射率的预测结果。

图16(a)(c)在图15(a)(c)的基础上在对角线上增加一线状地物，反射率固定为0.5，模拟一条道路陆地之上和桥梁水体之上，图16(d)-(f)分别由图16(a)-(c)聚合而成的低空间分辨率影像。利用图16(b)之外的5幅影像重新生成图16(b)，图16(g)显示了MDCM和STARFM重新生成的图16(b)，可见MDCM和STARFM都能定性地预测出该线状地物。如图17所示为MDCM和STARFM对图16(b)中线状地物的反射率的预测结果。从图17所示两种方法对图16(b)的线状地物预测的定量结果看，STARFM并不能完全正确地得到线状地物的反射率，错误主要出现在线状地物周围区域发生了变化的一段圆形地物之外，而MDCM方法能得到完全正确的线状地物的反射率。

小面积地物的情况。小面积地物在实际影像中也是非常普遍存在的，特别是在地物斑块比较破碎的区域，比如小块农田或小面积水体等。如果这些小面积地物的尺寸小于了低空间分辨率像元的尺度，则在低空间分辨率影像中不能清晰观察到它们的形状以及得到它们的反射率信息，而重新生成的高空间分辨率影像希望能准确地反映这些小地物的形状和反射率信息。如图18所示为模拟小面积地物的情况。(a)-(c)中的圆形区域的反射率从0.1变为0.2再到0.4，各图中4个圆形地物的直径依次为5、10、15、20个高空间分辨率单位，周围区域的反射率固定为0.05：(d)-(f)分别由(a)-(c)聚合而成的低空间分辨率影像：(g)和(h)分别为STARFM和MDCM方法对(b)的预测结果。

假设只存在两类地物，图18(a)-(c)中的圆形地物及周围背景。图中4个圆形地物的直径依次为5、10、15、20的4个小面积地物，圆形地物的反射率从0.1图18(a)增加到0.2图18(b)再到0.4图18(c)，周围区域的反射率固定为0.05。图18(d)-(f)是将图18(a)-(c)按17*17聚合成低空间分辨率的影像。图18(g)和(h)分别是STARFM和MDCM方法对图18(b)的重建结果，而如图19所示为MDCM和STARFM对图18(b)中各圆形小地物的预测误差。图19显示了二者对不同直径小地物的预测平均误差小地物各像元反射率预测值与真实值之差的平均值，可见STARFM只有直径为20的圆形地物的反射率是正确的，其他小于低空间分辨率像元尺度17的地物的反射率都比真实值偏小，而MDCM方法对所有尺寸的地物都能生成正确的反射率。原因主要是STARFM需要纯像元提供变化信息，而小于低空间分辨率像元尺度的小面积地物在低空间分辨率影像上不存在纯像元，所以生成的反射率出现误差，而MDCM方法如果在窗口内找不到纯像元时，通过这些小面积地物反射率变化量和其所在的低空间分辨率像元反射率变化量之间的比例系数，将非纯像元提供的变化信息加以了修正，所以能准确地生成这些小面积地物的反射率。

如图8所示，遥感影像的处理***包括以下模块：影像选择模块101，用于获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像：目标像元设定模块102，用于设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域：目标像

元反射率计算模块104，用于计算目标像元的反射率：和，影像生成模块108，用于根据目标像元反射率计算模块104计算得到的各目标像元的反射率生成A0影像。

所述目标像元反射率计算模块104包括以下子模块：

转换系数模块103，根据所选的影像计算所述目标区域的反射率变化量在两组影像中之间的转换系数。如图9所示，在本实施例中，所述转换系数模块103包括：中心像元设定模块1031，用于将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元：相似像元筛选模块1032，用于从所述中心像元各自所在影像中分别筛选出至少一个与所述中心像元属于同类地物的相似像元：转换系数计算模块1033，用于确定各相似像元所对应的地表区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述转换系数。调整模块1034，用于调整误差导致的可能不正确的转换系数和修正奇异的转换系数。调整模块1034将在下文结合图14进行详细介绍。在本实施例中，所述转换系数计算模块1033中所述转换系数的计算是将所述所获得的反射率以第二变化量为自变量，第一变化量为因变量根据各相似像元的权重进行加权回归分析。相似像元的筛选的具体实现可以参见上文中步骤1022的说明，在此不再赘述。

权重模块109，用于根据所选的影像计算各相似像元的权重。下文将结合图11和图12进行详细介绍。

第二变化量计算模块105，用于从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。如图10所示，所述第二变化量计算模块105包括：相似像元第二变化量计算模块1051，用于从第二组影像中选择与所述某时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算各所述相似区域在第二组影像中的反射率变化量：第二变化量加权平均模块1052，根据得到的所述权重将相似像元第二变化量计算模块1051得到的各反射率变化量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。

变化量转换模块106，用于根据所述转换系数将所述各第二变化量计算模块105得到的各变化量转换为在第一组影像中的各反射率变化量。

反射率计算模块107，用于根据变化量转换模块106得到的各反射率变化量和目标区域在第一组相应时刻的各影像中的反射率计算所述目标区域的反射率。在本实施例中经由第二变化量加权平均模块1052、变化量转换模块106和反射率计算模块107完成了式14的计算，得到了目标像元的反射率。影像生成模块108，用于逐次计算所有待计算的像元的反射率生成A0影像。

在另一个实施例中，目标像元反射率计算模块104还包括：时间权重计算模块111，用于根据第二变化量计算模块105所得到的变化量计算所得到的所述目标区域的各反射率的时间权重。反射率计算模块107还包括：反射率加权计算模块，用于根据所述时间权重，将各反射率加权平均得到所求目标像元的反射率。时间权重的计算和最终目标像元反射率的计算的具体实现可以参见上文中步骤10810中的说明，在此不再赘述。

如图11所示，所述权重模块109包括：纯度计算模块1091，用于根据所选影像计算各相似像元所有波段光谱向量和第二组影像中各所述相似区域的像元所有波段光谱向量的一致程度：纯像元判断模块1093，用于判断相似像元中是否有纯像元，并当有纯像元时将该纯像元权重设定为最大值，即，所有变化信息均取自该纯像元，当没有纯像元时不调整各相似像元的权重：空间距离计算模块1094，用于计算各相似像元距离各自中心像元的空间距离：权重计算模块1095，用于根据各相似像元的纯度和距各自中心像元的空间距离计算各相似像元的权重。权重模块109中各模块的具体实现可以参见上文中步骤104的说明，在此不再赘述。

如图12所示，在另一个实施例中，所述权重模块109包括：纯度计算模块1091，用于根据所述至少两对时刻对应的影像计算各相似像元所有波段光谱向量和与其时刻对应的第二组影像中地表区域对应的像元所有波段光谱向量的一致程度：相似像元权重计算模块1092，用于根据所述一致程度计算各相似像元的权重。具体纯度计算和权重计算可以参考权重模块109的上一个实施例，在此不再赘述。

如图13所示，在另一个实施例中，所述相似像元筛选模块1032还包括：交集筛选模块10321，用于将各筛选出的相似像元所对应的地表区域取交集，以取交集后的地表区域所对应的像元为相似像元。交集筛选模块10321的具体实现可以参见上文中步骤1024的说明，在此不再赘述。

如图14所示，所述调整模块1034包括：均值判断模块10341，用于判断各相似像元对应的第二变化量的均值是否小于反射率自身误差：误差调整模块10342，用于当均值判断模块(10341)判断为小于时，转换系数设定为使第一变化量等于第二变化量，当均值判断模块(10341)判断为不小于时，转换系数保持不变：奇异转换系数判断模块10343，用于判断转换系数是否奇异：奇异转换系数修正模块10344，用于将奇异的转换系数修正成。调整模块1034的具体实现可以参见上文中步骤1064和步骤1066的说明，在此不再赘述。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种遥感影像的处理方法，其特征在于，设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域，获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像，按以下步骤分别计算各目标像元的反射率以生成A0影像：

a，根据所选的影像计算所述目标区域的反射率变化量在两组影像中之间的转换系数：

b，从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像：根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量：

c，根据所述转换系数将所述各反射率的变化量转换为在第一组影像中的各反射率变化量：和

d，从第一组影像中选择与所述其它时刻对应的影像，根据转换得到的各反射率变化量和目标区域在所选影像中的反射率计算所述目标区域的反射率。

2.根据权利要求1所述方法，其特征在于，所述步骤a包括：

a1，将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元：

a2，从所述中心像元各自所在影像中分别筛选出至少一个与所述中心像元属于同类地物的相似像元：和

a3，确定各相似像元所对应的地表区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述转换系数。

3.根据权利要求2所述方法，其特征在于，所述步骤a3之前还包括：a21，根据所选的影像计算各相似像元的权重：

所述步骤b包括：

b1，根据B0影像和所述第二组的其他时刻的影像计算各所述相似区域在第二组影像中的反射率变化量：

b2，根据得到的所述权重将步骤b1得到的各反射率变化量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。

4.根据权利要求2所述方法，其特征在于，所述步骤a3中所述转换系数的计算是将各所述相似区域在所选影像中的反射率以在第二组影像中的反射率的变化量为自变量，在第一组影像中的反射率的变化量为因变量进行回归分析。

5.根据权利要求4所述方法，其特征在于，所述步骤a3之前还包括：a22，根据所选的影像计算各相似像元的权重：所述回归分析是根据各相似像元的权重进行的加权回归分析。

6.根据权利要求3或5所述方法，其特征在于，所述各相似像元的权重的计算包括以下步骤：

a221，根据所选影像计算各相似像元所有波段光谱向量和第二组影像中各所述相似区域的像元所有波段光谱向量的一致程度：

a222，根据所述一致程度计算各相似像元的权重。

7.根据权利要求2所述方法，其特征在于，所述步骤a2还包括：将各筛选出的相似像元所对应的地表区域取交集，以取交集后的地表区域所对应的像元为相似像元。

8.遥感影像的处理***，其特征在于，根据同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像求得高空间分辨率的A0影像，所述A0影像不属于第一组影像，并由至少一个待计算的像元组成，所述处理***包括：

影像选择模块(101)，用于获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像：目标像元设定模块(102)，用于设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域：

目标像元反射率计算模块(104)，用于计算目标像元的反射率：和影像生成模块(108)，用于根据目标像元反射率计算模块(104)计算得到的各目标像元的反射率生成A0影像：

所述目标像元反射率计算模块(104)包括以下子模块：

转换系数模块(103)，根据所选的影像计算所述目标区域的反射率变化量在两组影像中之间的转换系数：

第二变化量计算模块(105)，用于从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量：

变化量转换模块(106)，用于根据所述转换系数将所述各第二变化量计算模块(105)得到的各变化量转换为在第一组影像中的各反射率变化量：和反射率计算模块(107)，用于从第一组影像中选择与所述其它时刻对应的影像，根据变化量转换模块(106)得到的各反射率变化量和目标区域在所选影像中的反射率计算所述目标区域的反射率。

9.根据权利要求8所述***，其特征在于，所述转换系数模块(103)包括：

中心像元设定模块(1031)，用于将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元：

相似像元筛选模块(1032)，用于从所述中心像元各自所在影像中分别筛选出至少一个与所述中心像元属于同类地物的相似像元：和转换系数计算模块(1033)，用于确定各相似像元所对应的地表区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述转换系数。

10.根据权利要求9所述***，其特征在于，所述目标像元反射率计算模块(104)还包括：权重模块，用于根据所选的影像计算各相似像元的权重：

所述第二变化量计算模块(105)包括：

相似像元第二变化量计算模块(1051)，用于从第二组影像中选择与所述某时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的影像计算各所述相似区域在第二组影像中的反射率变化量：

第二变化量加权平均模块(1052)，根据得到的所述权重将相似像元第二变化量计算模块(1051)得到的各反射率变化量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量。

Images