遥感影像的处理方法及系统的制作方法

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

发明内容
有鉴于此，本发明的主要目的在于提供一种遥感影像的处理方 法，以实现能生成将细节保持的更好的影像。
为解决上述技术问题，本发明提供了一种遥感影像的处理方法， 设定某时刻的所求影像AO的一像元为目标像元，确定目标像元对应 的地表区域为目标区域，获取同一地表区域的具有高空间分辨率的 第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中 对应的选择至少两对不同时刻的影像，按以下步骤分别计算各目标 像元的反射率以生成AO影像a，根据所选的影像计算所述目标区 域在两组影像中的反射率变化量之间的转换系数；b,从第二组影像
6中选择与所述影像A 0时刻对应的影像B 0和至少 一 个其他时刻的影 像；根据所选择的影像计算所述目标区域在各所述其他时刻的影像 中的反射率相对于其在影像BO中的反射率的变化量；c，根据所述 转换系数将所述各反射率的变化量转换为在第 一组影像中的各反射 率变化量；和d,从第一组影像中选择与所述其它时刻对应的影像， 根据转换得到的各反射率变化量和目标区域在所选影像中的反射率 计算所述目标区域的反射率。
由至少两对时刻相对应的影像就可以得到时刻对应的两组影像 中地表区域反射率的变化信息之间的关系，即转换系数。根据所挑 的影像经由步骤a得到目标区域的转换系数。因为变化信息之间的 关系在较短时间范围内可以认为是稳定的。所以只要所选的各影像 的时刻和影像AO时刻相差不大，则可以根据转换系数经由步骤c将 步骤b所得的目标区域在第二组影像中的反射率变化量(简称第二 变化量)转换为目标区域在第一组影像中的反射率变化量(简称第 一变化量)。由此，由步骤d所得的反射率是基于更精细的第一组 影像中的反射率和转换后的变化量生成的。这样生成的AO影像就很 好的保持了第一组影像中的细节，避免所生成的影像在视觉上有被 平滑的效果。
优选的是，所述步骤a包括al，将所选的第一组影像中的与所 述目标区域对应的像元i殳定为中心像元；a2， 乂人所述中心像元各自 所在影像中分别筛选出至少 一 个与所述中心像元属于同类地物的相 似^f象元；和a3,确定各相似4象元所对应的地表区域为相似区域，根 据所述相似区域在所选影像中的反射率计算所述转换系数。
混合像元反射率的变化量为各端元反射率的变化量的综合反映。 所以端元可以提供比混合像元更为准确的变化信息。又因为作为端 元的相似像元与中心像元属于同类地物，所以，相似像元能够提供 最为准确的变化信息，由此可以得到更为准确的转换系数。其中， 中心像元也算是相似像元之一，在未能筛选出其它相似像元时，中 心像元是唯一的相似4象元。尤其是对于目标区域中存在石皮碎斑块和小地物的地表区域，其影像在第二组影像中为混合像元。其提供的 变化信息(第二变化量)被由各相似像元得到的转换系数转换修正 成与目标像元属于同类地物的端元的变化信息(第一变化量)。由 此，目标像元的反射率准确度更高了 ，尤其是对破碎斑块和小地物 具有更好的融合效果。
优选的是，所述步骤a3之前还包括a21,根据所选的影像计算 各相似像元的权重；所述步骤b包括bl,根据BO影像和所述第二 组的其他时刻的影像计算各所述相似区域在第二组影像中的反射率 变化量；b2，根据得到的所述权重将步骤bl得到的各反射率变化量 加权平均计算所述目标区域在各所述其他时刻的影像中的反射率相 对于其在影像B0中的反射率的变化量。
目标区域的第二变化量不是直接由目标区域在B0影像和所述第 二组的其他时刻的影像的反射率求差得到，而是由各相似像元按照 各自权重共同来提供。由此修正了 BO影像和第二组的其他时刻的影 像中目标区域对应的混合像元所提供的变化信息。从而得到更准确 的第二变化量。
优选的是，所述步骤a3中所述转换系数的计算是将各所述相似 区域在所选影像中的反射率以在第二组影像中的反射率的变化量为 自变量，在第一组影像中的反射率的变化量为因变量进行回归分析。
利用回归分析可以消除因噪声等原因引起的误差，从而得到更稳 健的变换系数，减少可能的误差。
优选的是，所述步骤a3之前还包括a22，根据所选的影像计算 各相似像元的权重；所述回归分析是根据各相似像元的权重进行的 力口权回归分析。
各相似像元距离其中心像元的空间距离，以及其的光谱特征和中 心像元的光谱特征的一致性等原因都导致了各相似像元提供变化信 息时，存在贡献大小的差异。在回归分析时，利用各相似像元的权 重进行加权回归分析，可以体现出这些差异。进而得到更准确的转 换系数，以进一步提高目标像元的反射率准确度。
8优选的是，所述各相似像元的权重的计算包括以下步骤a221， 根据所选影像计算各相似像元所有波段光谱向量和第二组影像中各 所述相似区域的像元所有波段光谱向量的一致程度；a222,根据所 述一致程度计算各相似像元的权重。
所述第一组影像的像元为高空间分辨率像元，所述第二组影像的 像元为低空间分辨率像元，高空间分辨率像元对应地表区域所属地 物在对应的低空间分辨率像元对应地表区域的各种地物中所占比例 即为高空间分辨率像元的纯度。所以相似像元的纯度越大，也就意 味着其所提供的变化信息越准确。从而通过步骤a221和步骤a222 可以提高能提供更准确变化信息的纯度大的相似像元的权重。而纯 度大小的计算则是由所有波段光谱向量的 一 致程度计算得到。
优选的是，所述步骤a2还包括将各筛选出的相似像元所对应 的地表区域取交集，以取交集后的地表区域所对应的像元为相似像 元。
有的地物随时间会发生变化，其光谱特征也相应发生变化，只利 用 一组时刻对应的影像筛选出的相似像元可能会因此筛选出已和中 心像元不属于同一地物的相似像元。所以取交集后可以筛选出的各 时刻的均和中心像元都是同 一地物的相似像元。从而得到更可靠的 目标像元。
本发明还相应的提供了 一种遥感影像的处理系统，根据所述两组 影像求得高空间分辨率的A 0影像，所述A 0影像不属于第 一 组影像， 并由至少一个待计算的像元组成，所述处理系统包括影像选择模 块101 ,用于获取同一地表区域的具有高空间分辨率的第一组影像 和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择 至少两对不同时刻的影像；目标像元设定模块102 ,用于设定某时 刻的所求影像AO的一像元为目标像元，确定目标像元对应的地表区 域为目标区域；目标像元反射率计算模块104 ,用于计算目标像元 的反射率；和影像生成模块108 ,用于根据目标像元反射率计算模 块104计算得到的各目标像元的反射率生成AO影像；所述目标像元反射率计算模块104包括以下子模块转换系数模块103 ，根 据所选的影像计算所述目标区域在两组影像中的反射率变化量之间 的转换系数；第二变化量计算模块105 ,用于从第二组影像中选择 与所述影像AO时刻对应的影像BO和至少一个其他时刻的影像，根 据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反 射率相对于其在影像BO中的反射率的变化量；变化量转换模块 106 ,用于根据所述转换系数将所述各第二变化量计算模块105得 到的各变化量转换为在第一组影像中的各反射率变化量；和反射率 计算模块107 ，用于从第一组影像中选择与所述其它时刻对应的影
像，根据转换得到的各反射率变化量和目标区域在所选影像中的反 射率计算所述目标区域的反射率。
优选的是，所述转换系数模块103包括中心像元设定模块 1031 ，用于将所选的第一组影像中的与所述目标区域对应的像元设 定为中心像元；相似像元筛选模块1032 ,用于从所述中心像元各 自所在影像中分别筛选出至少 一个与所述中心像元属于同类地物的 相似像元；和转换系数计算模块1033 ，用于确定各相似像元所对 应的地表区域为相似区域，根据所述相似区域在所选影像中的反射 率计算所述转换系数。
优选的是，所述目标像元反射率计算模块104还包括权重模 块，用于根据所选的影像计算各相似像元的权重；所述第二变化量 计算模块105包括相似像元第二变化量计算模块1051 ,用于从 第二组影像中选择与所述某时刻对应的影像BO和至少一个其他时 刻的影像，根据所选择的影像计算各所述相似区域在第二组影像中 的反射率变化量；第二变化量加权平均模块1052 ,根据得到的所 述权重将相似像元第二变化量计算模块1051得到的各反射率变化 量加权平均计算所述目标区域在各所述其他时刻的影像中的反射率 相对于其在影像BO中的反射率的变化量。


图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)中各圓形小地物的预测 误差；
图20为MODIS上排和Landsat下排假彩色合成影像；
图21为分别由STARFM和MDCM方法生成的7月11日Landsat
分辨率的假彩色合成影像；
图22为真实影像与重建影像的各波段散点图。
具体实施例方式
在对实施方式进行说明之前，先对本发明的理论基础进行说明。 对于同一区域来自不同传感器的遥感数据经过辐射定标、几何配 准和大气纠正后，这些数据之间具有一定的可比性和相关性。但是， 由于不同传感器在波段宽度、光谱响应函数、获取时刻的大气状况等方面存在差异，使得这些多源数据间存在一定的系统偏差。在对 系统误差进行纠正的同时，如何利用它们之间的相关性实现多源数 据的融合构成了本发明的核心。以下分两种情况介绍本发明理论基 础。
对于均一地表情况假设低空间分辨率、高时间分辨率(略写为 低空间分辨率)和高空间分辨率、低时间分辨率(略写为高空间分 辨率)传感器有类似的光谱波段设置。当低空间分辨率影像重采样 为与高空间分辨率影像相同的空间分辨率(即相同的像元大小)和坐 标系统时，对于同一且均一的地物而言，低空间分辨率像元内仅为 一种地物类型纯像元，此时在波段S上高空间分辨率反射率与低空 间分辨率反射率之间的关系即为高、低空间分辨率影像上反射率之 间的相对辐射定标关系，该关系仅由两种传感器特性差异(波段宽度 及光谱响应函数)和成像时刻大气状况差异所决定， 一般可视为以下 线性关系
,力,用="x C(x,.,力,^， 5) + 6 i
其中，F、 C分别代表高、低空间分辨率的像元反射率，a,，力) 是像元位置，S为光谱波段，？k代表了影像获取时刻(以日为单位)， "、6分别代表高空间分辨率和低空间分辨率影像相对辐射定标的增 益和偏离值。对于同一且均一的地物而言，l式所示关系在任何影像 获耳又时刻都成立。
如果已经获取了 Z。时刻波段B的高、低空间分辨率的影像和/p 时刻同一波段的低空间分辨率影像，那么/p时刻该波段高空间分辨 率的像元反射率可由式2来表示~'，乃^，5) = """，^ )" 2
而/。时刻高、低空间分辨率影像的像元反射率关系为
F(、 ， ^, ^ ， B) = a x C(x,,力，~, B) + 6 3
由式2和式3相减并移项整理可得到式4: ^",力,^，5) = 7^,，3;乂0,5) + "><((:"，力,^,￡)一(:",3;-。,5)) 4
上式可以理解为/p时刻高空间分辨率的像元反射率由？o时刻高 空间分辨率的反射率加上 相对？。时刻的反射率变化量。式中fl代表高空间分辨率和低空间分辨率影像相对辐射定标的增益，也可理 解为高空间分辨率和低空间分辨率像元反射率之间的转换系数，该 系数由两种传感器特性差异和成像时刻大气状况差异所决定。如果 两种传感器成像时刻相近，可视为大气状况相同，则该系数仅由两 种传感器特性差异所决定，且不随时间发生变化。
如果获取任意两个时刻Zm、 /n的两对高、低空间分辨率影像，那 么通过对这两个时刻的高、低空间分辨率图象上的各波段纯像元反
射率进行线性回归即可求出各波段的转换系数a,通过该系数和式4
即可由 一 系列时间连续的低空间分辨率影像生成一 系列时间连续的 高空间分辨率的反射率影像。
对于非均一地表情况由于传感器空间分辨率的限制以及地物的
复杂多样性，混合像元普遍存在于低空间分辨率影像上，对地表地
物类型多样、分布杂乱的区域尤其如此。此时，低空间分辨率像元 与其对应的高空间分辨率像元的反射率并不是对同 一 种地物的反
应，二者之间的差异不仅仅是辐射定标的差异，此时基于均一地表
纯像元的式4并不适合非均一地表混合像元的情况。
混合像元的光谱特征是其内部各类地物光谱端元光谱的综合反 映，所以两个时刻/m、 ^之间其反射率的变化为各类端元光谱变化的 综合反映。目前混合像元光谱混合模型最常用的是线性光谱混合模 型，即利用 一 个线性关系来表达混合像元和各端元光谱之间的关系。
假设低空间分辨率影像上fm、 ？n时刻的某波段混合像元反射率分别为
4、 K,根据线性光谱混合模型，IV K可由下式表达
其中，M为该混合像元包含的端元数目，《M和；^分别表示^、
^时刻第i种端元的光谱，S为误差项，且考虑S在Zm、 Zn时刻不变。
由式5可得到变化量其中下式中的力是第/个端元的比例；
13,=i 6
假设/m、 ^的时间间隔在一定范围内时，各端元反射率的变化可
以近似为线性变化，则端元光谱《 可由式7表示其中下式中的n' 是第/个端元反射率的线性变化率；X, =。xA"X,.w 7
其中，Vm,将式7代入式6可得到下式^="￡〃 8 已知第A:端元/m、/n时刻的反射率，则由式7可求出△〖& = ^"_义" 将式9代入式8，可得%^ = + 10
由于假设Zm、 ^的时间间隔在一定范围内，各端元比例/和各端 元反射率的线性变化率r在,m、 Zn时刻之间可视为稳定，则式10的 右边为一稳定变量，用V表示，其意义为A类地物端元反射率变化 量与混合像元反射率变化量的比值变化比例系数。
对于非均一地表，低空间分辨率像元为混合像元，其对应的高空
间分辨率影像上的各像元视为端元。根据式10， Zm、 Z。时刻的高空间 分辨率像元反射率尸(JC,.，》，^， B)和F(X,.， 乂， ^， 5)与对应的低空
间分辨率像元反射率C(x,.， 乂.， l 5)和C(x,.，为，"，B)满足以下关系..
通过式11可知，/m、 ^时刻的各波段高空间分辨率反射率与对应 的低空间分辨率反射率为线性关系，进行线性回归即可求出各波段 的变化比例系数v。
如果在fm到rn内，已经获取了 /。时刻波段B的高、低空间分辨 率的影像和fp时刻同一波段的低空间分辨率影像。根据以上的假设， ^相比"时刻该波段高空间分辨率的像元与低空间分辨率像元的反
射率变化量之间的比例系数仍为v则c^，乂，/p，旬-c"，乂，/。，5) 12 通过式12即可求解出^时刻的高空间分辨率反射率F"，力，,p,5卜F(x";^。,万)+ vx(C(x"力，^5) —C"，力J。,5)) 13
虽然公式13与公式4的形式相同，但其意义和适用情况不同。 式4表示的是纯像元不同分辨率反射率之间相对定标的关系，在任 何时刻都成立，其计算的结果最为准确。而式13表示的是混合像元 中不同分辨率反射率变化量之间的比例关系，根据假设只在^和& 之间或其附近时刻才成立，而且是基于短时间范围内各地物反射率 的变化为线'f生的假设，所以其计算的结果是近似解。
基于上述的理论基础，当低空间分辨率像元为纯像元时，可以利 用式4求解出^时刻的高空间分辨率反射率，而当低空间分辨率像 元为混合像元时利用式13求解。但纯像元计算的结果较混合像元准 确可靠。为尽量利用纯像元的信息，以及减小影像受到云、大气等 污染给计算结果带来的不确定性，本实施例中采取窗口内求解的办 法。即，以像元XW2， ， ^W2)为中心，取宽度为W的窗口，先筛选出
与中心像元属同类地物的像元相似像元，然后利用式14计算中心像 元^时刻的高空间分辨率反射率
F(Uw2，^5",( 2,尺/"^) + iXxrX(C(H^S)-C(《，乂，^邻
m 14 其中，iV为窗口内与中心4象元属于同类地物的相似^象元数目，『 为各相似像元的权重，F为高、低空间分辨率的转换系数。式14的 意义是/p时刻的反射率由已知时刻的反射率加上^时刻相对已知/0 时刻的变化量决定，而该变化量由窗口内的相似像元共同预测。
在本实施例之前，通常需要先对获取的高、低空间分辨率影像进 行预处理，使得多源数据之间具有相似的波段设置，相同的像元大 小和图幅尺寸以及相同的坐标系统。采用常用的软件和方法即可实 现。下文将结合附图对遥感影像的处理方法的实施方式进行详细说 明。
本实施方式中遥感影像为对同一地表摄像而得到的两组影像，其 中，第一组影像具有高空间分辨率，第二组影像具有低空间分辨率， 根据所述两组影像求得tp时刻的高空间分辨率的A影像，所述A影 像由至少一个待计算的像元组成。如图1所示，遥感影像的处理方法包括以下步骤步骤100，从所述遥感影像中挑选至少两对时刻
对应的影像。本实施例中，选用的是已知的tm、 tn时刻的两对影像。
从第二组影像中挑选其时刻与tp时刻所对应影像；步骤101，设定 一个待计算的像元为目标像元，确定目标像元的应的地表区域为目 标区域；步骤102，以目标区域对应的像元为中心像元，筛选出已知 时刻的高空间分辨率影像中所述中心像元的窗口内的同类地物像元 得到相似像元，其中，中心像元也算是相似像元之一，在未能筛选 出其它相似像元时，中心像元是唯一的相似像元；步骤104，计算所 述相似像元的权重。本实施例中，根据/m、 /w时刻的两幅低空间分辨 率影像，计算出步骤102中筛选出的相似像元的权重；步骤106,根 据所述权重采用加权算法对已知时刻的高、低空间分辨率反射率进 行线性回归计算得到转换系数；步骤108,计算预测时刻的目标像元 的高空间分辨率反射率，本实施例中，根据转换系数和预测时刻tp 的低空间分辨率影像计算已知时刻高空间分辨率影像的时间权重， 并根据已知时刻高空间分辨率影像按照所述时间权重加权计算得到 预测时刻的目标像元的高空间分辨率反射率；步骤IIO,逐次计算所 有目标像元并最终生成A影像。下文将结合附图对上述各个步骤进
行详细说明o
图2为本实施例中筛选相似像元的流程图，如图所示包括 步骤1022,以目标区域对应的像元为中心4象元，筛选出各已知 时刻的高空间分辨率影像中所述中心像元的窗口内的与目标像元属 于同类地物的像元得到相似像元。在目标像元( )为中心
的窗口内，与要计算的中心像元属于同类地物的高空间分辨率像元 才能提供比较正确可靠的反射率变化信息。可以使用以下两种方法
筛选出相似像元(a)对高空间分辨率影像进行非监督分类，和中 心像元类型相同的像元为相似像元；(b)4吏用阈值判断，以窗口内 其它像元反射率与中心像元反射率之差来判断是否属于同类地物。 这两种筛选相似像元的方法具有共同的意义，即是以光谱相似的像 元为同一类地物。但是又有不同之处阈值法是在窗口内寻找与中心像元的光谱相似的像元，中心像元作为了搜寻的中心，而且其筛
选条件是局部适用的，随着中心像元的位置变化而变；而分类是在 全影像进行，如果分类有误差，那么判断出的相似像元可能具有很 大的错误。所以在本实施例中，使用阈值法来筛选相似像元。在本 实施例中，判断窗口内某像元所有波段都满足式15的条件，则该像 元被确 定为中 心像元的相似像元 。
其中，o(S)为&时刻高空间分辨率影像所有像元S波段的标准差， m为预估的地物类别数，比如影像覆盖的区域主要有植被、棵地、 水体、岩石4类地物，则可将m值设为4。当然，也可以用所有波 段中的部分波段进行筛选或者采用方法(a)进行筛选。
步骤1024，将筛选出的两组已知时刻的相似像元取交集。有的 地物随时间会发生变化，其光谱特征也相应发生变化，只利用一个 时相的影像筛选出的相似像元可能会出现错误。比如窗口内有棵地 和农作物两种地物类型，如果中心像元是农作物，^时刻农作物并 未生长，其光谱特征与棵地相同，此时筛选出来的相似像元是棵地， 如果"时刻农作物已生长，则在"时刻的影像上能筛选出正确的农 作物像元。鉴于此，本实施例中分别利用fM、 /"时刻两幅高空间分辨 率影像筛选相似像元，并将筛选结果取交集，以此提高筛选结果的 准确性。如图3所示为取交集筛选相似像元示意图，"*"表示中心 像元，"O"表示中心像元的相似像元。
图4为本实施例中计算相似像元权重的流程，如图所示包括
步骤10402,计算相似像元的纯度大小。权重『决定了窗口内筛 选出的各相似像元对计算目标像元反射率的贡献。在本实施例中权 重『由各相似像元的纯度大小和离中心像元的空间距离远近共同决 定。纯度越高、空间距离越近的相似像元权重越大，这是因为均一 地表的纯像元能提供最准确的变化信息，而且认为空间上离中心像 元越近的相似像元的变化情况和中心像元越 一 致。
高空间分辨率影像上的某像元代表了某种地物类型，而其对应位
17置低空间分辨率影像的像元包含了该类地物，纯度即是指低空间分 辨率 <象元中该种地物所占比重。如果在均 一 地表情况下低空间分辨
率像元全由该种地物覆盖，则认为是纯像元。比如Landsat影像上某
像元所覆盖的30mx 30m的区域为小麦，其对应位置的MODIS影像
上像元所覆盖的500mx 500m的区域也全为小麦，则认为该像元为
小麦纯像元。由于不同的地物具有不同的光谱特征，所以可以通过
比较高空间分辨率像元和其对应位置的低空间分辨率像元的光谱曲
线的一致程度来判断纯度大小。如果高空间分辨率像元和其对应的
低空间分辨率像元的光谱一致程度越高，则认为纯度越大。而高、
低空间分辨率像元的光谱曲线之间的 一致性可由相关系数来刻画，
所以筛选出的各相似像元的纯度i 可由式16计算 及=譜(F,，C,)
'一V^y,(C,) 16
F, ={F(x,,_y,，,m,JB1),...,F(x,.,;;,,,m,JBJ,F(x,j,,, ,JB1),...,JF(x,，j;,， ,5 )} 其中向量Fi， Ci分别表示第/个相似像元高、低空间分辨率的光
镨向量，cov为求协方差，D为求方差。i 值的范围为-l到l, iM直
越大表示纯度越高。之所以将G、 4时刻的光谱曲线放在一起计算相 关系数，是考虑到地物随物候的变化，光谱特征也会发生变化，只 使用一个时相的光谱可能会出现错误。比如小麦在返青之前其光谱 特征和棵地一致，如果某小麦高空间分辨率像元对应的低空间分辨 率像元除了小麦外还包含棵地，则用小麦返青前的获取的影像计算 会得到较高的纯度，显然是不正确的。当加入小麦返青后获取的影 像进行判断时，就会得到更加正确的纯度。由上可以看出，只要获 得至少两组已知时刻的高空间分辨率像元的所有波段光谱向量和在 其对应位置同时的低空间分辨率像元的所有波段光谱向量，就可以 根据其一致程度来计算纯度大小。
利用纯度大小作为计算权重的一个要素，较之目前采用同时刻 高、低空间分辨率反射率之差表示的光谱特性作为计算权重的要素将更加有效的利用相似像元并且得到的结果也更加准确。
步骤10404,判断相似像元中是否有纯像元，有则进入步骤 10406，没有则进入步骤10408。
步骤10406，设定纯像元的权重为1。在本实施例中，权重W的 取值范围是0到1，所有相似像元的权重之和为1。 R=l的像元为 纯像元。为了尽量利用纯像元，将纯像元的权重设为1。即，中心像 元的反射率变化信息完全由这些纯像元提供。在另一个实施例中， 当相似像元中有P个纯像元存在即R=l的像元时，规定这些纯像元 的权重W为1/P，其他相似像元的权重为0。当然，纯像元的判断也 可以是在计算权重过程中的其它步骤中实现。例如在步骤10402之 前或者在步骤10408之后。甚至不判断是否有纯像元也可以。
步骤10408,计算相似像元距离中心像元的空间距离。相似像元 ，")与窗口中心像元;w2, ， JV2)的空间距离^由式17计算
《=1 + V"W2 - A )2 + - X )2 2) "
上式包含了对空间距离的归一化，在窗口 w内各相似像元的空间距 离《的取值范围为1到l+2G's,值越大，该相似像元离中心像元越远。 步骤10410，计算相似像元的权重。接下来将纯度和空间距离这 两个不相关的量联系在一起，来确定各相似像元的权重『。上文已 经提到纯度越高、距离越近的像元能提供更准确的信息，即A越大， d越小，其权重『越大。本实施例中，通过公式18，将纯度和空间 距离结合为总距离D: A=(l-^>《 18
Z)越大的相似像元为目标像元反射率的计算贡献越小，再将D的倒
『'=(i/A)/i](i/A)
数归一化得到各相似像元的权重『式19: '=1 19
当然，建立纯度、空间和权重的关系可以采用多种形式，只要能
实现纯度越大，距中心像元越近的相似像元的权重越大即可。 图5为得到转换系数的流程，如图所示包括 步骤1062,计算转换系数。由公式4和公式13可知，转换系数
K包含了相对定标的增益a和高低空间分辨率变化量的比例系数v,
19二者都可由已经获取的fw、"时刻高、低空间分辨率各波段反射率进 行线性回归得到。为了减小噪声污染等对转换系数计算带来的不确 定性，将所有相似像元^、 ^时刻高、低空间分辨率反射率进行线性
回归，回归直线的斜率即是转换系数F。为了利用更可靠的相似像元 的信息计算转换系数，采用加权最小二乘法进行线性回归，权重为 各相似像元的『。如图6所示为转换系数的计算示意图。某窗口内 有12个像元和中心像元属于同类地物包括中心像元本身，虚线框分 别标注出"、4时刻的反射率，可见从^到,"时刻，高空间分辨率 反射率和其对应的低空间分辨率反射率都有所增加，但是高空间分
辨率反射率增加幅度较大，是低空间分辨率的6.448倍^=0.915, 尸<0.001。这说明该类地物该波段反射率在/w、"时刻之间发生了较 大的变化，比如植被生长导致的近红外反射率急剧增加。根据理论 基础的论述，各相似像元都反映了低、高空间分辨率反射率之间的 转换系数，但由于噪声等影响，通常不能只利用某一个相似像元来 计算V，这样的结果带有很大的不确定性，而线性回归能利用所有 相似像元的信息得到 一个比较稳健的转换系数。又因为各相似像元 存在贡献大小的差异，为了体现这些差异，在线性回归之中进一步 还考虑各样本的权重，本实施例采用的是加权最小二乘法。当然， 也可以釆用其它的回归方法。
步骤1064，根据不确定性调整转换系数。由于几何配准、大气 修正等数据处理过程都给反射率带来一 定的误差，如果低空间分辨 率像元^、"时刻之间的反射率变化太小，在误差范围内，则进行加 权回归得到的斜率J/具有很大的不确定性。这可能是以下两种情况 低空间分辨率像元内反射率发生变化的地物所占面积比例很小，在 整个低空间分辨率像元尺度上体现不出这种变化；也可能是低空间 分辨率像元内有的地物反射率增加，有的地物反射率减少，而增加 和减少的量相等，使得在低空间分辨率像元尺度上体现不出变化来。 为解决这个问题，在本实施例中，在所有相似像元的低空间分辨率 反射率变化均值小于至少两个已知时刻反射率的变化量的不确定性时将转换系数设为1。假设反射率各波段的不确定性U为该波段最大 值的1%,而且这种不确定性在各影像各波段间是互相独立的，那么
两个已知时刻反射率的变化量的不确定性为u = AX 20
其中，Utm是Um分别是^、"时刻反射率的不确定性。如果所有 相似像元的低空间分辨率反射率变化均值小于式20计算出的不确定 性，则只能将第二变化量平均分配到其内部的每一个高空间分辨率 像元上，即变化量的转换系数r为l。当然，已知时刻反射率的变 化量的不确定性也可以通过其它统计方法得到。对式20的不同表达 方式，以及其它的变换的计算方式都应包含在本发明的保护范围之 内。这样就避免了两个时刻间反射率的变化太小时，因为反射率自 身带有一定的不确定性，也就是误差，而导致计算出不真实的转换 系数V。
步骤1066，修正奇异的转换系数。由于影像自身噪声的存在， 导致绝对值较大的转换系数出现。为修正个别较奇异的转换系数， 对全幅影像所有像元的转换系数进行统计，在本实施例中将均值上 下2倍标准差以外的转换系数r设定为1。当然，也可以根据需要将 其它预设倍数的标准差以外的转换系数设定为1。
图7为计算预测时刻的目标像元的高空间反射率的流程，如图所 示包括步骤10802，根据预测时刻的低空间分辨率影像计算各相似 像元相对已知时刻的第二变化量；步骤10804,根据所述转换系数将 所述第二变化量转换为第一变化量；步骤10806，将各相似像元的第 一变化量按所述权重加权平均后得到目标像元的第一变化量；和， 步骤10808 ，将所述目标像元的第 一变化量与已知时刻的高空间分辨 率反射率相加得到基于不同已知时刻的预测时刻目标像元的高空间 分辨率反射率。通过上述4个步骤完成了式14的计算。以两组分别 为"、^的已知时刻所获取的时刻之间或附近的任意^时刻的低空间 分辨率影像，计算出^时刻的高空间分辨率反射率为例。根据式14, 由已知的Zw、"时刻高、低空间分辨率影像即可计算出 时刻的高空间分辨率反射率，此时/w、 ^时刻的高、低空间分辨率反射率都已知， 所以可以分别由/w、 时刻的高空间分辨率反射率计算出Zp时刻的高
空间分辨率反射率Fjcw/2, — 5,将这两个计算结果分别记
作Fwxw/2, ~， B、 i^W2， JV2， ~, 5。当然，只要已知转
换系数、预测时刻的低空间分辨率影像和各相似单元的权重就可以 在已知时刻的高空间分辨率影像的基础上计算得到预测时刻的目标
像元的高空间分辨率反射率。对式14的不同表达方式，以及其它的
变换的计算方式都应包含在本发明的保护范围之内。
步骤10810,加权计算得到预测时刻的高空间分辨率反射率。为 充分利用G、 /"时刻的高空间分辨率反射率的信息，将两个计算结果 加权求和得到最后的预测结果。假设/p时刻的低空间分辨率反射率 相对已知时刻变化越小，则其计算的结果越准确可靠。这是因为如 果低空间分辨率影像上没有体现出反射率的变化，则认为高空间分 辨率影像也不会有反射率的变化。基于以上假设，两个计算结果的 时间权重T由式21计算
<formula>formula see original document page 22</formula>得到时间权重T后，^时刻的高空间分辨率反射率F:cw/2，少w， ^， S可由下式计算
F(^2,3^2，,p,5^7;xFm(^2,;^2,/p，5) + 7;xF"(^2,:^2^,5) 22
当然，也可以采用其它方法来计算基于不同已知时刻计算得到的
高空间分辨率反射率的权重。例如，如式23，用预测时刻/。和已知
时刻4低空间分辨率反射率之差来计算权重大小。
<formula>formula see original document page 22</formula>
至此得到了最终的预测时刻的目标像元的高空间分辨率反射率。
仅是单独应用转换系数配合步骤108就能改善对小地物和线状地物 的计算精度。而相似像元的筛选以及权重计算的改进则有助于整体
精度的提高。为了更好的对本实施方式更进一步说明和验证，下面对本实施方 式分别进行模拟数据和真实数据的检验。
模拟数据检验。为了更好的检验方法的精确性和可靠性下面将本
实施例运用于简单的模拟数据。并且为了便于和GAO的STARFM 方法进行比较，利用GAO文中的模拟数据进行检验。在下文中本实 施方式简称为MDCM。
物候变化的情况。物候是引起反射率变化的主要原因，尤其在植 被覆盖率较大的区域。如图15所示为模拟物候变化的情况。其中 (a)-(c)中的圆形区域的反射率固定为0.05,周围区域的反射率从O.l 变为0.2再到0.4;(d)-(f)分别由(a)-(c)聚合而成的低空间分辨率影像； (g)和(h)分别是利用MDCM和STARFM对(b)的预测结果。筒单假设 只存在两类地物，分别为水体和植被，而且形状固定不变。假设圆 形的水体具有固定的反射率0.05，而水体四周植被的反射率从0.1 图15(a)增加到0.2图15(b)再到0.4图15(c),代表随着物候变化，植 被的生长导致的近红外反射率增加。图15(d)- (f)是将图15(a)- (c)按 17*17聚合成低空间分辨率的影像(a)- (c)类比于30米分辨率的 landsat, (d)- (f)类比于500米分辨率的MODIS。利用图15(a)和图 15(c)-(f)预测图15(b)。图15(g)和(h)分别是MDCM方法和STARFM 的结果，图15(g)和(h)与图15(b)完全一致，说明两种方法生成的高 空间分辨率反射率都能准确地反映由于物候引起的反射率变化。物 候变化所引起的反射率变化是真实影像中最普遍存在的现象，说明 MDCM方法和STARFM都能运用于真实影像，生成高时空分辨率反 射率影像。
线状地物的情况。在土地覆盖或土地利用类型中，线状地物相当 普遍，比如道路、桥梁、河流等，而人眼对线状地物的敏感程度较 高，生成的高分辨率影像中线状地物的正确与否直接影响影像的视 觉效果。如图16所示为模拟线状地物的情况。图16。 (a)-(c)中的圆 形区域的反射率固定为0.05,直线的反射率固定为0.5,周围区域的 反射率从O.l变为0.2再到0.4; (d)-(f)分别由(a)-(c)聚合而成的低空
23间分辨率影^象；(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)中的圆形区域的反射率从O.l变为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个小面积地物， 圆形地物的反射率从O.l图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方法如果在窗口内找 不到纯像元时，通过这些小面积地物反射率变化量和其所在的低空 间分辨率像元反射率变化量之间的比例系数，将非纯像元提供的变 化信息加以了修正，所以能准确地生成这些小面积地物的反射率。
真实数据检验。将MDCM方法运用于真实的Landsat-7和MODIS 遥感影像。同模拟数据一样，为了与Gao等的STARFM进行比较， 使用 Gao 等人文中的遥感数据进行检验数据可在 http:〃ledaps.nascom.nasa.gov/ledaps/Tools/StarFM.htm下载。Landsat-7 ETM+数据的空间分辨率为30米，MODIS数据的空间分辨率为500 米，波段为绿光、红光、近红外3个波段MODIS的4、 1、 2波段， 分别对应ETM+的2、 3、 4波段。
数据所覆盖区域的地理位置是54° N， 104° W,大小为12kmx 12km,主要地表覆盖类型为植被、水体、棵地、道路等，植被的生 长季较短，物候变化较快。影像共有3组，获取时刻分别为2001年 5月24日、2001年7月11日、2001年8月12日。如图20为MODIS 上排和Landsat下排假彩色合成影像近红-红-绿波段对应红-绿-蓝彩 色合成。从左到右获取时刻依次是2001年5月24日、2001年7月 11日、2001年8月12日。上面一排为MODIS影像，下面一排为 Landsat影像。利用5月24日和8月12日的两组影像和7月11曰 的MODIS影像生成7月11日的Landsat分辨率的影像，并将结果与 真实的7月11日Landsat影像进行比较，来评价MDCM方法的精度。 根据该区域特点，在本实施例中将窗口 w设置为3 x 3个MODIS像元即50 x 50个Landsat像元，地物类型数量w设置为4。图21显示 了分别由STARFM和MDCM方法生成的7月11日Landsat分辨率 的假彩色合成影像。如图21所示为真实的7月11日Landsat影像左 和由STARFM中和MDCM方法右生成的Landsat分辨率影像。
由图21可见，MDCM方法生成的影像和真实影像视觉上很接近， 细节特征保持较好，而STARFM的结果有些细节特征没有突出，在 视觉上有影像被平滑的效果。这主要是因为微小地物在对应的低空 间分辨率上没有对应的纯像元存在，对于这些小地物^象元，STARFM 不能生成正确高分辨率反射率，而MDCM方法能克服这个问题。
为定量比较两种方法结果的精度，现计算出预测反射率和真实反 射率之间的绝对误差，即预测反射率和真实反射率之差的绝对值。 如图22所示为真实影像与重建影像的各波段散点图。上一排为 STARFM计算的结果，下一排为MDCM计算的结果。R是相关系数， AD为平均绝对误差。
图22显示了两种方法的结果各波段与真实值之间的散点图。对 于绿光波段green band, MDCM方法的结果相比真实值有所偏高， 使得整幅影像的绝对误差均值AD为0.0039,略高于STARFM的 0.0033;对于红光波段red band, MDCM方法的结果比STARFM略 好一点，预测值与真实值的相关系数R为0.90,绝对误差为0.0034, 而STARFM的R为0.89,绝对误差为0.0037,从散点图来看，MDCM 方法的结果偏高，而STARFM的结果偏低；对于近红外波段NIR band, MDCM方法的结果明显优于STARFM,预测值与真实值的相 关系数R为0.95，绝对误差为0.0103,而STARFM的R为0.85， 绝对误差为0.0155,从散点图来看，MDCM方法的结果和真实值基 本一致，而STARFM有很多像元的反射率值比真实值低，而这些低 估的像元正是影像中的小面积地物板块。
本发明还相应的提供了 一种遥感影像的处理系统。遥感影像为对 同一地表摄像而得到的两组影像，其中，第一组影像具有高空间分 辨率，第二组影像具有低空间分辨率，根据所述两组影像求得tp时刻的高空间分辨率的A影像，所述A影像由至少一个待计算的像元 组成。下面结合图8~图14对遥感影像的处理系统的实施方式进行介 绍。
如图8所示，遥感影像的处理系统包括以下模块影像选择模块 101 ，用于获取同 一地表区域的具有高空间分辨率的第 一组影像和具 有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少 两对不同时刻的影像；目标像元设定模块102,用于设定某时刻的所 求影像AO的一像元为目标像元，确定目标像元对应的地表区域为目 标区域；目标像元反射率计算模块104,用于计算目标像元的反射率； 和，影像生成模块108,用于根据目标像元反射率计算模块104计算 得到的各目标像元的反射率生成AO影像。
所述目标像元反射率计算模块104包括以下子模块
转换系数模块103，根据所选的影像计算所述目标区域的反射率 变化量在两组影像中之间的转换系数。如图9所示，在本实施例中， 所述转换系数模块103包括中心像元设定模块1031，用于将所选 的第一组影像中的与所述目标区域对应的像元设定为中心像元；相 似像元筛选模块1032，用于从所述中心像元各自所在影像中分别筛 选出至少 一 个与所述中心像元属于同类地物的相似像元；转换系数 计算才莫块1033,用于确定各相似像元所对应的地表区域为相似区域， 根据所述相似区域在所选影像中的反射率计算所述转换系数。调整 模块1034，用于调整误差导致的可能不正确的转换系数和修正奇异 的转换系数。调整模块1034将在下文结合图14进行详细介绍。在 本实施例中，所述转换系数计算模块1033中所述转换系数的计算是 将所述所获得的反射率以第二变化量为自变量，第一变化量为因变 量根据各相似像元的权重进行加权回归分析。相似像元的篩选的具 体实现可以参见上文中步骤1022的说明，在此不再赘述。
权重模块109，用于根据所选的影像计算各相似像元的权重。下 文将结合图11和图12进行详细介绍。
第二变化量计算模块105，用于从第二组影像中选择与所述影像AO时刻对应的影像B0和至少一个其他时刻的影像，根据所选择的 影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于 其在影像BO中的反射率的变化量。如图10所示，所述第二变化量 计算模块105包括相似像元第二变化量计算模块1051,用于从第 二组影像中选择与所述某时刻对应的影像BO和至少一个其他时刻 的影像，根据所选择的影像计算各所述相似区域在第二组影像中的 反射率变化量；第二变化量加权平均模块1052,根据得到的所述权 重将相似像元第二变化量计算模块1051得到的各反射率变化量加权
平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于 其在影像BO中的反射率的变化量。
变化量转换模块106，用于根据所述转换系数将所述各第二变化 量计算模块105得到的各变化量转换为在第一组影像中的各反射率 变化量。
反射率计算模块107,用于根据变化量转换模块106得到的各反 射率变化量和目标区域在第 一组相应时刻的各影像中的反射率计算 所述目标区域的反射率。在本实施例中经由第二变化量加权平均模 块1052、变化量转换模块106和反射率计算模块107完成了式14 的计算，得到了目标像元的反射率。影像生成模块108,用于逐次计 算所有待计算的像元的反射率生成AO影像。
在另一个实施例中，目标像元反射率计算模块104还包括时间 权重计算模块111，用于根据第二变化量计算模块105所得到的变化 量计算所得到的所述目标区域的各反射率的时间权重。反射率计算
模块107还包括反射率加权计算模块，用于根据所述时间权重， 将各反射率加权平均得到所求目标像元的反射率。时间权重的计算 和最终目标像元反射率的计算的具体实现可以参见上文中步骤 10810中的说明，在此不再赘述。
如图ll所示，所述权重模块109包括纯度计算模块1091,用 于根据所选影像计算各相似像元所有波段光谱向量和第二组影像中 各所述相似区域的像元所有波段光谱向量的 一致程度；纯像元判断
28模块1093,用于判断相似像元中是否有纯像元，并当有纯像元时将 该纯像元权重设定为最大值，即，所有变化信息均取自该纯像元， 当没有纯像元时不调整各相似像元的权重；空间距离计算模块1094, 用于计算各相似像元距离各自中心像元的空间距离；权重计算模块 1095，用于根据各相似像元的纯度和距各自中心像元的空间距离计 算各相似像元的权重。权重模块109中各模块的具体实现可以参见上 文中步骤104的说明，在此不再赘述。
如图12所示，在另一个实施例中，所述权重模块109包括纯度 计算模块1091，用于根据所述至少两对时刻对应的影像计算各相似像 元所有波段光谱向量和与其时刻对应的第二组影像中地表区域对应的 像元所有波段光谱向量的一致程度；相似像元权重计算模块1092，用 于根据所述一致程度计算各相似像元的权重。具体纯度计算和权重计 算可以参考权重模块109的上一个实施例，在此不再赘述。
如图13所示，在另一个实施例中，所述相似4象元筛选才莫块1032 还包括交集筛选模块10321,用于将各筛选出的相似像元所对应的地 表区域取交集，以取交集后的地表区域所对应的像元为相似像元。交 集筛选模块10321的具体实现可以参见上文中步骤1024的说明，在此 不再赘述。
如图14所示，所述调整模块1034包括均值判断模块10341,用 于判断各相似像元对应的第二变化量的均值是否小于反射率自身误 差；误差调整模块10342，用于当均值判断模块(10341 )判断为小于 时，转换系数设定为使第一变化量等于第二变化量，当均值判断模块 (10341 )判断为不小于时，转换系数保持不变；奇异转换系数判断模 块10343，用于判断转换系数是否奇异；奇异转换系数修正模块10344, 用于将奇异的转换系数修正成。调整模块1034的具体实现可以参见上 文中步骤1064和步骤1066的说明，在此不再赘述。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明， 凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等， 均应包含在本发明的保护范围之内。
29
权利要求
1. 一种遥感影像的处理方法，其特征在于，设定某时刻的所求影像A0的一像元为目标像元，确定目标像元对应的地表区域为目标区域，获取同一地表区域的具有高空间分辨率的第一组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像，按以下步骤分别计算各目标像元的反射率以生成A0影像a，根据所选的影像计算所述目标区域的反射率变化量在两组影像中之间的转换系数；b，从第二组影像中选择与所述影像A0时刻对应的影像B0和至少一个其他时刻的影像；根据所选择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像B0中的反射率的变化量；c，根据所述转换系数将所述各反射率的变化量转换为在第一组影像中的各反射率变化量；和d，从第一组影像中选择与所述其它时刻对应的影像，根据转换得到的各反射率变化量和目标区域在所选影像中的反射率计算所述目标区域的反射率。
2. 根据权利要求l所述方法，其特征在于，所述步骤a包括 al,将所选的第一组影像中的与所述目标区域对应的像元设定为中心像元；a2,从所述中心像元各自所在影像中分别筛选出至少 一个与所述 中心像元属于同类地物的相似像元；和a3，确定各相似像元所对应的地表区域为相似区域，根据所述相 似区域在所选影像中的反射率计算所述转换系数。
3. 根据权利要求2所述方法，其特征在于，所述步骤a3之前还 包括a21,根据所选的影像计算各相似像元的权重；所述步骤b包括b 1 ，根据B0影像和所述第二组的其他时刻的影像计算各所述相 似区域在第二组影像中的反射率变化量；b2,根据得到的所述权重将步骤bl得到的各反射率变化量加权 平均计算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影像BO中的反射率的变化量。
4. 根据权利要求2所述方法，其特征在于，所述步骤a3中所述 转换系数的计算是将各所述相似区域在所选影像中的反射率以在第 二组影像中的反射率的变化量为自变量，在第一组影像中的反射率 的变化量为因变量进行回归分析。
5. 根据权利要求4所述方法，其特征在于，所述步骤a3之前还 包括a22,根据所选的影像计算各相似像元的权重；所述回归分析是才艮据各相似像元的权重进行的加权回归分析。
6. 根据权利要求3或5所述方法，其特征在于，所述各相似像 元的权重的计算包括以下步骤a221 ，根据所选影像计算各相似像元所有波段光谱向量和第二组 影像中各所述相似区域的像元所有波段光谱向量的 一 致程度； a222,根据所述一致程度计算各相似像元的权重。
7. 根据权利要求2所述方法，其特征在于，所述步骤a2还包括 将各筛选出的相似像元所对应的地表区域取交集，以取交集后的地表区域所对应的^f象元为相似4象元。
8. 遥感影像的处理系统，其特征在于，根据所述两组影像求得 高空间分辨率的AO影像，所述AO影像不属于第一组影像，并由至 少一个待计算的像元组成，所述处理系统包括影像选择模块(101),用于获取同一地表区域的具有高空间分 辨率的第 一 组影像和具有低空间分辨率的第二组影像，从所述两组影像中对应的选择至少两对不同时刻的影像；目标像元设定模块(102),用于设定某时刻的所求影像A0的一像元为目标-像元，确定目标像元对应的地表区域为目标区域；目标像元反射率计算模块(104),用于计算目标像元的反射率；和影像生成模块(108 ),用于根据目标像元反射率计算模块(104 ) 计算得到的各目标像元的反射率生成AO影像；所述目标像元反射率计算模块(104)包括以下子模块转换系数模块(103)，根据所选的影像计算所述目标区域的反 射率变化量在两组影像中之间的转换系数；第二变化量计算模块(105)，用于从第二组影像中选择与所述 影像A 0时刻对应的影像B 0和至少 一 个其他时刻的影像，根据所选 择的影像计算所述目标区域在各所述其他时刻的影像中的反射率相 对于其在影像BO中的反射率的变化量；变化量转换模块(106),用于根据所述转换系数将所述各第二 变化量计算模块(105)得到的各变化量转换为在第一组影像中的各 反射率变化量；和反射率计算模块(107)，用于从第一组影像中选择与所述其它 时刻对应的影像，根据变化量转换模块(106)得到的各反射率变化 量和目标区域在所选影像中的反射率计算所述目标区域的反射率。
9. 根据权利要求8所述系统，其特征在于，所述转换系数模块 (103 )包括中心像元设定模块(1031),用于将所选的第一组影像中的与所 述目标区域对应的像元设定为中心像元；相似像元筛选模块(1032),用于从所述中心像元各自所在影像 中分别筛选出至少 一个与所述中心像元属于同类地物的相似像元；和转换系数计算模块(1033 ),用于确定各相似像元所对应的地表 区域为相似区域，根据所述相似区域在所选影像中的反射率计算所述 转换系数。
10. 根据权利要求9所述系统，其特征在于，所述目标像元反射 率计算模块(104)还包括权重模块，用于根据所选的影像计算各相似像元的权重； 所述第二变化量计算模块(105)包括相似像元第二变化量计算模块(1051 ),用于从第二组影像中选 择与所述某时刻对应的影像B0和至少一个其他时刻的影像，根据所 选择的影像计算各所述相似区域在第二组影像中的反射率变化量；第二变化量加权平均模块(1052)，根据得到的所述权重将相似 像元第二变化量计算模块(1051 )得到的各反射率变化量加权平均计 算所述目标区域在各所述其他时刻的影像中的反射率相对于其在影 像BO中的反射率的变化量。
全文摘要
本发明公开了一种遥感影像的处理方法，根据空间分辨率不同的两组影像求得某时刻的所求影像A0。首先根据空间分辨率不同的两组影像计算目标区域的反射率变化量在这两组影像之间的转换系数，其次根据所述转换系数将低空间分辨率影像中的变化量转换成高空间分辨率影像中的变化量，最后由此根据已知的高空间分辨率影像得到所求影像A0的目标像元的反射率。还相应的提供了遥感影像的处理系统。本发明能生成将细节保持的更好的影像。
文档编号G06K9/62GK101482929SQ20091007937
公开日2009年7月15日 申请日期2009年3月9日 优先权日2009年3月9日
发明者朱孝林, 晋 陈, 陈仲新, 陈学弘 申请人:中国农业科学院农业资源与农业区划研究所;北京师范大学