技术领域
[0001] 本
发明属于
图像处理领域，具体涉及一种高分辨率遥感图像辐射重建方法。
背景技术
[0002] 当前的高空间分辨率遥感影像辐射重建方法一般是利用基于反射率的
时空融合技术来处理得到，时空融合技术的本质是通过将基于多个采集时相的低空间-高
时间分辨率遥感影像所建立的辐射关联模式应用到对应时相的高空间-低时间分辨率遥感影像，进而利用已知时相的高空间-低时间分辨率影像以及多时相的低空间-高时间分辨率数据来预测得到目标时相下的高空间-低时间分辨率遥感影像。在上述数据预测过程中，地物的纹理特征与
光谱特征是获得高
质量融合结果的两个关键要素，当前已发展出一些能够有效合成光谱特征的融合方法(如STARFM(The Spatial and Temporal Adaptive Reflectance Fusion Model)
算法)，但鉴于高空间分辨率影像和低空间分辨率影像之间的空间尺度差异问题，针对纹理特征的融合过程描述尚不成熟，一般需依赖多时序(两个及以上)的高空间分辨率数据(如ESTARFM(Enhanced STARFM)算法)或处理效率较低(如基于稀疏学习的时空融合算法)。
发明内容
[0003] 本发明克服
现有技术存在的不足，所要解决的技术问题为：提供一种更准确的高分辨率遥感图像辐射重建方法。
[0004] 为了解决上述技术问题，本发明采用的技术方案为：一种高分辨率遥感图像辐射重建方法，其特征在于，包括以下步骤：
[0005] S1：获得待重建区域t1时刻下的高分辨率遥感图像L1，以及t1、t2时刻下的低分辨率遥感图像M1、M2，其中，L1、M1、M2三幅图像的波段数量相同，均为b，图像尺寸也相同，均为n×n，待重建的图像为该区域t2时刻下的高分辨率遥感图像L2；
[0006] S2：通过降
采样方法或者像元聚合技术将L1处理成为尺寸为(n/2)×(n/2)的图像L′1；并对L′1的全波段图像进行非监督分类，得到包含m个类别的分类结果图像K；
[0007] S3：分别对原始的高分辨率遥感图像L1和低分辨率遥感图像M1、M2进行逐波段的离散
小波变换处理，对高分辨率遥感图像L1处理后得到高频分量L1(bi,HH)、低频分量L1(bi,LL)以及高低频的混合分量L1(bi,HL)和L1(bi,LH)，对低分辨率遥感图像M1、M2处理后分别得到高频分量M1(bi,HH)、低频分量M1(bi,LL)、高低频的混合分量M1(bi,HL)和M1(bi,LH)以及高频分量M2(bi,HH)、低频分量M2(bi,LL)、高低频的混合分量M2(bi,HL)和M2(bi,LH)；其中bi代表波段；
[0008] S4：根据高频分量L1(bi,HH)、M1(bi,HH)、M2(bi,HH)计算L2的高频分量L2(bi,HH)；根据低频分量M1(bi,LL)和M2(bi,LL)以及高低频的混合分量M1(bi,HL)、M1(bi,LH)、M2(bi,HL)和M2(bi,LH)计算L2的低频分量L2(bi,LL)、高低频的混合分量L2(bi,HL)和L2(bi,LH)；
[0009] S5：利用步骤S4所求解得到的L2的高频分量L2(bi,HH)、低频分量L2(bi,LL)、高低频的混合分量L2(bi,HL)和L2(bi,LH)进行上述小波变换的逆变换，重建出t2时刻尺寸为n×n的高分辨率遥感图像L2。
[0010] 所述步骤S4中，计算得到L2的高频图像L2(bi,HH)的具体步骤为：
[0011] S401：在高频分量图像L1(bi,HH)中设置一个以待处理像元p(x,y)为中心的z×z个像元大小的滑动窗口，其中，z＝INT(r/s)，r为低分辨率遥感图像M1的分辨率，s表示高分辨率遥感图像L1的分辨率，bi表示波段；
[0012] S402：计算待处理像元p(x,y)与滑动窗口内的坐标为(e,f)的相似像元j的相对空间距离d′，计算公式为：
[0013]
[0014] 其中，d表示像元p(x,y)与其相似像元j的欧几里得距离，其计算公式为：
[0015] S403：计算待处理像元p(x,y)的相似像元j的权重函数Wj(p)，计算公式为：
[0016]
[0017] S404：分波段计算得到L2的高频分量图像L2(bi,HH，p)，计算公式为：
[0018]
[0019] 其中，j表示待处理像元p(x,y)在滑动窗口内的相似像元的编号，J为滑动窗口内的光谱相似像元数量，其通过借助分类图像K在此滑动窗口内搜寻与待处理像元p(x,y)属于同一类别的邻近像元得到。
[0020] 所述步骤S4中，计算得到L2的低频分量图像L2(bi,LL)、高低频的混合图像和L2(bi,HL)，L2(bi,LH)的具体步骤为：根据关系式M2(bi,LL,k)＝u×M1(bi,LL,k)+v，进行最小二乘回归方法计算得到分类结果图像K中的类别k的线性模型的斜率u和截距参数v，然后根据线性模型的斜率u和截距参数v计算得到类别k的低频分量图像L2(bi,LL,k)、高低频的混合图像L2(bi,HL,k)和L2(bi,LH,k)，计算公式为：
[0021] L2(bi,LL,k)＝u×L1(bi,LL,k)+v；
[0022] L2(bi,HL,k)＝u×L1(bi,HL,k)+υ；
[0023] L2(bi,LH,k)＝u×L1(bi,LH,k)+v；
[0024] 针对分类结果图像K中的每一个类别均采用上述求解步骤则可得到所有类别的低频分量图像L2(bi,LL)、高低频混合图像L2(bi,HL)和L2(bi,LH)。
[0025] 所述步骤2中，采用的ISODATA或K-Means对L′1的全波段图像进行非监督分类。
[0026] 所述步骤3中，分别对原始的高分辨率遥感图像L1和低分辨率遥感图像M1、M2进行逐波段的
离散小波变换处理时，小波变换处理过程为一层小波变换或二层小波变换或二层以上的小波变换。
[0027] 本发明与现有技术相比具有以下有益效果：
[0028] 本发明首先利用小波变换将所有时相下的对应波段数据转换为1个高频分量(HH)与3个包含低频的分量(HL、LH、LL)来提升融合算法对于纹理特征的描述
精度，再采用邻域光谱贡献法来处理高频分量，以及利用逐地类的线性回归建模来处理3个包含低频的分量，以此来最大程度的保留高、低空间分辨率影像间的辐射光谱融合精度；在利用上述两种处理模式分别处理得到目标时相下高分辨率影像的1个高频分量与3个包含低频的分量数据后，再利用小波逆变换来合成目标时相下的高空间分辨率反射率数据，本方法只需单个时相的高空间分辨率遥感影像以及两个时相的低空间分辨率遥感影像即可处理得到高空间分辨率的遥感影像，本方法的核心是将基于单数据对的时间-空间反射率融合数据分解成高频与包含低频的两类分量并进行针对性的时空融合处理，从而可以更准确的描述、表达和处理多源遥感影像的时间、空间(尤其是纹理特征)以及光谱特征信息，进而可以更加准确的预测地物在空间维的变化信息。
具体实施方式
[0029] 为使本发明
实施例的目的、技术方案和优点更加清楚，下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例；基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0030] 本发明实施例提供了一种高分辨率遥感图像辐射重建方法，包括以下步骤：
[0031] S1：获得某地表区域在t1时刻下的高分辨率遥感图像L1，以及t1、t2时刻下的低分辨率遥感图像M1、M2，其中，L1、M1、M2三幅图像的波段数量相同，均为b，图像尺寸也相同，均为n×n，假设待重建的图像为该区域t2时刻下的高分辨率遥感图像L2。
[0032] S2：通过降采样方法或者像元聚合技术将L1处理成为尺寸为(n/2)×(n/2)的图像L′1；采用ISODATA(Iterative Selforganizing Data Analysis，
迭代自组织数据分析算法)或K-Means(k均值)对L′1的全波段图像进行非监督分类，得到包含m个类别的分类结果图像K；
[0033] S3：分别对原始的高分辨率遥感图像L1和低分辨率遥感图像M1、M2进行逐波段的一层离散小波变换处理，例如，对高分辨率遥感图像L1的波段bi进行小波变换处理后，得到高频分量L1(bi,HH)、低频分量L1(bi,LL)以及高低频的混合分量L1(bi,HL)和L1(bi,LH)，对低分辨率遥感图像M1、M2的波段bi进行小波变换处理后，处理后分别得到高频分量M1(bi,HH)、低频分量M1(bi,LL)、高低频的混合分量M1(bi,HL)和M1(bi,LH)以及高频分量M2(bi,HH)、低频分量M2(bi,LL)、高低频的混合分量M2(bi,HL)和M2(bi,LH)；其中bi代表波段；
[0034] S4：分别计算L2的高频分量L2(bi,HH)、低频分量L2(bi,HH)、高低频的混合分量L2(bi,HL)和L2(bi,LH)。
[0035] 其中，图像的高频部分主要是纹理结构特征，该特征一般可采用差值图像来表示，则高分辨率图像的高频分量在t1、t2时刻间的变化可通过两个时刻间低分辨率图像高频分量的差值图像来有效表达，即
[0036] L2(bi,HH)＝L1(bi,HH)+(M2(bi,HH)-M1(bi,HH))； (1)
[0037] 为了得到更为鲁棒的小波系数值，这里加入一个基于空间相对
位置关系的权重函数Wj(p)来描述与待处理像元p(x,y)周围邻近的某光谱相似像元j对像元(x,y)的光谱贡献，像元之间的光谱相似条件通过分类结果图像K来确定。加入权重函数后，L2(bi,HH)图像中任一像元p(x,y)的值可通过下式计算得到：
[0038]
[0039] 其中，光谱相似像元是通过在高频分量图像中设置一个以待处理像元p(x,y)为中心的z×z个像元大小的滑动窗口，然后借助分类图像K在此滑动窗口内搜寻与待处理像元(x,y)属于同一类别的邻近像元得到的，J为滑动窗口内的光谱相似像元数量，j＝1,2,3，……J，表示相似像元的编号；若原始高分辨率遥感图像L1的空间分辨率为s、原始低分辨率遥感图像M1和M2的空间分辨率为r，则z一般可设为INT(r/s)个像元(其中INT为取整操作)。若在此滑动窗口内的相似像元j的位置坐标为(e,f)，则定义一个描述待处理像元(x,y)与该相似像元(e,f)的欧几里得距离为：
[0040]
[0041] 其相对空间距离d′定义为：
[0042]
[0043] 最终的权重函数Wj(p)可定义为：
[0044]
[0045] 因此，如果在高频分量图像L1(bi,HH)中设置一个以待处理像元p(x,y)为中心的z×z个像元大小的滑动窗口，则通过式(3)和式(4)可以计算待处理像元p(x,y)与滑动窗口内的坐标为(e,f)的相似像元j的相对空间距离d′，代入式(5)，可以计算得到待处理像元p(x,y)的相似像元j的权重函数Wj(p)；将权重函数Wj(p)代入式(2)，可以得到高频分量L2(bi,HH)中任一像元p(x,y)的值L2(bi,HH,p)，计算出所有像元的值L2(bi,HH,p)即可以得到L2的高频分量L2(bi,HH)。
[0046] 一般认为，某地物类别k在t1、t2时刻的高分辨率图像的低频分量(光谱特征)间以及高低频混合分量间均呈线性相关，并且该地物类别k在t1、t2时刻的低分辨率图像的低频分量间也呈线性相关，并且它们的线性参数基本相同，即：
[0047] L2(bi,LL,k)＝u×L1(bi,LL,k)+v； (6)
[0048] L2(bi,HL,k)＝u×L1(bi,HL,k)+v； (7)
[0049] L2(bi,LH,k)＝u×L1(bi,LH,k)+v； (8)
[0050] M2(bi,LL,k)＝u×M1(bi,LL,k)+v； (9)
[0051] 其中，u、v为线性模型的斜率和截距参数，对于分类结果图K中的任意类别k，由式(9)的关系，可以根据已知的M1(bi,LL,k)、M2(bi,LL,k)利用最小二乘回归方法得到其对应的u、v，然后代入式(6)～式(8)，即可以计算得到类别k的低频分量L2(bi,LL,k)、高低频的混合分量L2(bi,HL,k)和L2(bi,LH,k)，针对分类结果图像K中的每一个类别均采用上述求解步骤则可得到所有类别的低频分量L2(bi,LL)、高低频混合分量L2(bi,HL)和L2(bi,LH)。
[0052] S5：利用步骤S4所求解得到的L2的高频分量L2(bi,HH)、低频分量L2(bi,LL)、高低频的混合分量L2(bi,HL)和L2(bi,LH)进行上述小波变换的逆变换，重建出t2时刻尺寸为n×n的高分辨率遥感图像L2。
[0053] 其中，本实施例中，分别对原始的高分辨率遥感图像L1和低分辨率遥感图像M1、M2进行逐波段的离散小波变换处理时，进行的是一层小波变换处理，得到高频分量L2(bi,HH)、低频分量L2(bi,LL)、高低频的混合分量L2(bi,HL)和L2(bi,LH)，进行一层小波逆变换，即可以重建出t2时刻尺寸为n×n的高分辨率遥感图像L2。本领域的技术人员在上述
说明书的教导下，也可以根据对原始的高分辨率遥感图像L1和低分辨率遥感图像M1、M2进行二层或二层以上的小波变换，当采用二层或二层以上小波变换时，需要对所有经小波分解得到的高频分量HH和低频分量LL均重复进行上述步骤的求解，但由于尺度效应与图像信息的损失，层数过高的小波变换会降低图像融合精度；小波变换方法可以采用多种小波模型，但均在本
专利的诉求范畴。
[0054] 此外，上述步骤S2中采用的ISODATA或K-Means等非监督分类方法，既可以换为最大似然等监督分类方法、
阈值分析法，也可以通过额外提供高精度的分类图像得到。
[0055] 进一步地，上述步骤S3中，对于高频分量L2(bi,HH)和低频分量L2(bi,LL)的求取，既可以采用上述方法，也可以采用现有的时空融合方法来各自独立求解，如时空自适应的反射率融合算法(Spatial-Temporal Adaptive Reflectance Fusion Model：STARFM等)。在计算过程中，只需要将时空自适应的反射率融合算法中的反射率参数换成本实施例中的高频分量HH，低频分量LL以及高低频混合分量HL和LH即可以分别求解。在求取过程中，L2(bi,HH)的计算需利用L1(bi,HH)、M1(bi,HH)、M2(bi,HH)三个分量图像，L2(bi,LL)的计算则需利用L1(bi,LL)、M1(bi,LL)、M2(bi,LL)三个分量图像，高低频分量L2(bi,LH)和L2(bi,HL)的计算则分别需利用L1(bi,LH)、M1(bi,LL)、M2(bi,LL)和L1(bi,LL)、M1(bi,LL)、M2(bi,LH)。
[0056] 最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行
修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
