CN116757925B - 一种生成高时空谱分辨率卫星遥感影像的方法和装置 - Google Patents

一种生成高时空谱分辨率卫星遥感影像的方法和装置 Download PDF

Info

Publication number
CN116757925B
CN116757925B CN202310546379.5A CN202310546379A CN116757925B CN 116757925 B CN116757925 B CN 116757925B CN 202310546379 A CN202310546379 A CN 202310546379A CN 116757925 B CN116757925 B CN 116757925B
Authority
CN
China
Prior art keywords
image
matrix
time
multispectral
resolution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310546379.5A
Other languages
English (en)
Other versions
CN116757925A (zh
Inventor
孙伟伟
周骏
孟祥超
杨刚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ningbo University
Original Assignee
Ningbo University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ningbo University filed Critical Ningbo University
Priority to CN202310546379.5A priority Critical patent/CN116757925B/zh
Publication of CN116757925A publication Critical patent/CN116757925A/zh
Priority to PCT/CN2024/094956 priority patent/WO2024208373A1/zh
Application granted granted Critical
Publication of CN116757925B publication Critical patent/CN116757925B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10036—Multispectral image; Hyperspectral image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
    • Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种生成高时空谱分辨率卫星遥感影像的方法和装置，包括：通过融合同一场景下具有高时间、空间分辨率但低光谱分辨率多光谱影像和具有高光谱分辨率但低时间、空间分辨率高光谱影像生成具有高时间、空间和光谱分辨率的融合影像。本发明的有益效果是：本发明相比当前一体化时空谱融合方法，提出的方法更符合当前主流星载遥感数据的分辨率特点，并能实现高保真的融合表现。

Description

一种生成高时空谱分辨率卫星遥感影像的方法和装置
技术领域
本发明涉及光学遥感图像领域，更确切地说，它涉及一种生成高时空谱分辨率卫星遥感影像的方法和装置。
背景技术
遥感具有实时、准确、连续地获取大范围的地表信息等特点，具有巨大的应用价值，在各个领域得到了广泛的应用，例如在自然资源监测，环境监测、矿物识别、精细农业等领域中都有着重要的应用。然而受限于成像传感器技术和成本的限制，卫星传感器在时间，空间和光谱分辨率间相互权衡，同时获取高时间、空间与光谱分辨率数据非常困难，限制了当前遥感数据在诸多领域的进一步应用，使得遥感数据呈现出“又多又少”的特点，可用的遥感数据越来越来越多，但真正能够被使用的数据相对总量比例仍然较低。
当前的卫星遥感成像中，高光谱影像具有较高的光谱分辨率，然而高光谱影像的空间和时间分辨率相对多光谱影像较低。而多光谱影像的光谱分辨率较低，波段数一般低于十个。我国的高光谱卫星ZY-1 02D于2020年发射，其高光谱传感器能够捕获空间分辨率为30米的166个光谱波段数据，重访周期为55天。低的空间分辨率(通常为30m)以及长的重访周期(长于15天)使得现有的高光谱卫星数据对于快速变化场景精细监测困难。相比于高光谱传感器，多光谱传感器通常具有更高的空间分辨率和时间分辨率，例如2015年和2017年发射的Sentinel-2a/b双星可以以5天的重访周期对地表进行访问，提供10米空间分辨率的4波段数据。多光谱卫星快速的重访周期和精细空间分辨率为长时序的快速监测提供了可能，但少量的波段数量难以进行精确的地物识别。遥感影像融合作为一种信息处理技术能够“以软补硬”，通过系统的综合同一区域具有互补信息的两幅或多幅遥感影像，生成一幅比可获得的单个影像提供更多视觉感知和计算机处理信息的影像，已成为当前提升单源卫星的时间/空间/光谱分辨率的主要手段之一。然而当前的遥感影像融合方法主要为针对遥感传感器三个分辨率属性中两者的融合(空间-光谱融合，时间-空间融合)，仅有的几种时空谱融合方法多是针对高时间分辨率高光谱影像的融合，未考虑到当前星载遥感影像的分辨率特点，即高光谱影像空间和时间分辨率相对多光谱影像更低的特点。
因此，基于当前遥感数据分辨率特点，开展多源遥感时-空-谱协同融合方法研究，集成多源遥感影像的时间、空间、光谱分辨率间的互补优势，融合生成高时-空-谱分辨率影像数据，对于遥感数据的进一步应用具有重要的意义。
发明内容
本发明的目的是针对当前遥感影像融合方法对当前遥感数据的分辨率特点考虑不足，不能充分发挥当前多源遥感数据的协同优势生成高时空谱分辨率数据，提出了一种生成高时空谱分辨率卫星遥感影像的方法和装置，实现高时间、空间和光谱分辨率数据的获取。
第一方面，提供了生成高时空谱分辨率卫星遥感影像的方法，包括：
步骤1、获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
步骤2、基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到预测时刻的上采样后的高光谱影像
步骤3、获取初始化的目标融合影像时间变化影像，表示为：
步骤4、基于光谱线性混合模型，将T2时刻的目标影像X2分解为端元矩阵E2和丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像；
步骤5、建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型；
步骤6、引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵；
步骤7、获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T；
步骤8、目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2。
作为优选，步骤1包括：
步骤1.1、获取T1时刻的高光谱影像和多光谱影像以及T2时刻的多光谱影像其中，l、L为波段数量，w、W为像素数量，l<L，w<W；
步骤1.2、对高光谱影像H1和多光谱影M1、M2进行归一化，并对归一化后的高光谱影像H1上采样至多光谱影像空间大小，得到
步骤1.3、将多光谱影像M1、M2和上采样后的高光谱影像沿光谱维折叠为三维形式，并分割为大小的全波段块，其中全波段块的宽和高为光谱维为l。
作为优选，步骤5中，多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型表示为：
Mi＝RXi＝REiAi
Hi＝XiFS＝EiAiFS
其中，i∈{1,2}，分别为在Ti时刻的高时空谱分辨率影像、高光谱影像以及多光谱影像；分别表示在Ti时刻的高时空谱分辨率影像的端元矩阵和丰度矩阵，Q为端元的个数，Q<<L；表示光谱降采样矩阵，和分别表现空间模糊矩阵和空间下采样矩阵。
作为优选，步骤7包括：
步骤7.1、分别引入关于目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T的正则化项，获得融合问题的能量函数；
步骤7.2、分割能量函数为对应于丰度矩阵A2、时间变化矩阵P和时间变化影像T变量的三个子函数，并为每个子函数引入若干辅助变量分裂子函数为子问题，为每个子问题构建增广拉格朗日函数；
步骤7.3、利用交替方向乘子法，对增广拉格朗日函数进行优化求解，获得目标影像的丰度矩阵A2，时间变化矩阵P和时间变化影像T。
作为优选，步骤7.1中，能量函数表示为：
s.t.D＝M2-M1
其中D为多光谱影像M2和M1的差分影像。作为优选，步骤7.1包括：
步骤7.1.1、引入关于目标影像的丰度矩阵A2的正则化项，该正则化项为空间平滑约束，表示为：
其中，Hh和Hv是线性算子，分别用于计算二维信号中相邻像素的分量之间的水平和垂直梯度，||·||2,1为混合L2,1范数，λA是丰度矩阵A2的正则化参数；
步骤7.1.2、引入时间变化矩阵P的正则化项，该正则化项为光谱平滑约束，表示为：
其中，H1是光谱维方向的差分算子，为L行Q列元素全为1的矩阵，λ1和λ2是时间变化P的正则化参数；
步骤7.1.3、引入时间变化影像T的正则化项，该正则化项为低秩约束，表示为：
其中，为T的三阶张量形式，||·||TNN为张量核范数，有其中表示第j个奇异值，通过对T进行傅里叶变换得到，λT为时间变化影像T的正则化参数。
作为优选，步骤7.2中，目标影像的丰度矩阵A2的子问题为：
其中η为时间变化影像项参数，用来控制时间变化影像对模型的贡献。
引入辅助变量B1＝A2，B2＝Hh(A1)，B3＝Hv(A1)，令u＝vec(A2)，符号vec(·)为矢量化操作，给出丰度A2的子问题的增广拉格朗日函数为：
其中，矩阵Vi,i＝1,2,3为拉格朗日函数的对偶变量，ρ为步长，有ρ>0，repM(V)为矩阵V的块对角矩阵，表示将矩阵V沿对角线复制M次。
作为优选，7.2中，目标影像的时间变化矩阵P的子问题为：
引入辅助变量B＝PeΕ1，时间变化矩阵P子问题的增广拉格朗日函数为：
目标影像的时间变化影像T的子问题为：
通过引入辅助变量给出时间变换影像T子问题的增广拉格朗日函数为：
作为优选，步骤2中，使用多光谱影像M1相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示：其中，ωi,j,k表示Ωi,j,k中第k块影像块的权重，ωi,j,k通过基于局部线性约束求取：
min||M2(i,j)-M1(Ωi,j)ωi,j||2+λ||di,j e ωi,j||2
s.t.1Tωi,j＝1
其中，ωi,j表示对应于M2(i,j)影像块的权重，||di,je ωi,j||2为局部线性约束正则化项，di,j为相似像素距离M1(i,j)中心的欧几里得距离，λ表示正则项参数；对于此局部线性约束的优化问题的解为：
ωi,j＝((Ci,j+λdiag(di,j))\1)(1T(Ci,j+λdiag(di,j))\1)-1
其中Ci,j＝(M1(Ωi,j)-M2(i,j)T)(M1(Ωi,j)-1M2(i,j)T)T表示影像块的协方差矩阵，diag(di,j)表示以di,j为对角线元素的对角矩阵；之后从多光谱影像对求取的权重ωi,j将被共享给上采样后高光谱影像以生成预测时刻的上采样的高光谱影像
第二方面，提供了一种生成高时空谱分辨率卫星遥感影像的装置，用于执行第一方面任一所述的生成高时空谱分辨率卫星遥感影像的方法，包括：
第一获取模块，用于获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
搜寻模块，用于基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到预测时刻的上采样后的高光谱影像
第二获取模块，用于获取初始化的目标融合影像时间变化影像，表示为：
分解模块，用于基于光谱线性混合模型，将T2时刻的目标影像X2分解为端元矩阵E2和丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像；
建立模块，用于建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型；
表示模块，用于引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵；
第三获取模块，用于获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T；
得到模块，用于令目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2。
本发明的有益效果是：
1.本发明主要针对于低时间、空间分辨率高光谱遥感影像与高时间、空间分辨率多光谱遥感影像的融合，相比现有遥感融合方法，本发明能够一次进行时间、空间和光谱三个分辨属性的融合，生成高时空谱分辨率融合影像。此外，本发明针对的融合任务更符合当前星载卫星遥感数据的分辨率特点。
2.本发明建立的时空谱融合方法提供了一种广义的融合框架，即包括了时间变化模型与时间变化影像分辨率增强重建模型，其中时间变化模型用于对影像间时间变化的刻画和建模，分辨率增强重建模型通过重建分辨率增强的时间变化影像模型以实现对时间变化的精确估计求取，通过时间变化模型与时间变化影像分辨率增强重建模型的耦合，从而实现了高保真的时空谱融合表现。
3.本发明利用了局部线性约束使用T1观测时间观测到的多光谱影像M1对T2预测时间的多光谱影像M2进行表示，将影像分块后，利用前一时刻M1(i，j)的全波段块周围搜索相似块M1(Ωi,j,k)来表示目标块M2(i,j)，得到权重系数，并将此学习到的权重系数共享于T1(观测时间)观测到的上采样后的高光谱影像上，以生成预测时刻的上采样高光谱影像充分的利用了遥感影像中具有的局部自相似性对时间变化特征进行学习。
4.本发明建立了时间变化影像模型，由和生成的初始化时间变化影像中，包含了两个时刻间目标影像的时间变化特征，将为之后的目标影像的估计提供精确的时间变化信息，以及高光谱分辨率的信息。
5.本发明基于广义线性混合模型来对时间变化进行刻画，之后与时间变化影像分辨率增强重建模型进行耦合提出了时空谱融合模型，将目标影像的估计，转为了估计目标影像的丰度矩阵A2，时间变化P，以及时间变化影像T。广义线性混合模型中的时间变化P能够对目标影像的时间变化进行捕获，时间变化影像T的估计能够作为对时间变化求取的约束。
6.本发明通过为目标影像丰度引入了空间平滑先验，为时间变化引入了光谱平滑先验，为时间变化影像模型引入了低秩先验，先验信息的引入使得模型可以得到精确的求解，从而得到高保真的融合影像X。
附图说明
图1为本发明的局部线性约束生成初始化时间变换影像示意图；
图2为本发明提出的时空谱一体化融合模型示意图。
具体实施方式
下面结合实施例对本发明做进一步描述。下述实施例的说明只是用于帮助理解本发明。应当指出，对于本技术领域的普通人员来说，在不脱离本发明原理的前提下，还可以对本发明进行若干修饰，这些改进和修饰也落入本发明权利要求的保护范围内。
实施例1：
本申请针对的时-空-谱融合问题，为从T1观测时间观测到的多光谱影像(有l个波段和W个像素)，和高光谱影像(有L个波段和w个像素)，与T2预测时间观测到的多光谱影像预测T2具有高空间、时间和光谱分辨率的目标影像其中l<L，w<W。本申请提供的方法包括以下步骤：
步骤1、获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
步骤1包括：
步骤1.1、获取T1时刻的高光谱影像和多光谱影像以及T2时刻的多光谱影像其中，l、L为波段数量，w、W为像素数量，l<L，w<W；
步骤1.2、对高光谱影像H1和多光谱影M1、M2进行归一化，并对归一化后的高光谱影像H1上采样至多光谱影像空间大小，得到高光谱影像H1上采样方式可以为双线性内插；
步骤1.3、将多光谱影像M1、M2和上采样后的高光谱影像沿光谱维折叠为三维形式，并分割为的全波段块。
步骤2、基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到预测时刻的上采样后的高光谱影像
步骤2中，使用多光谱影像M1相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示：其中，ωi,j,k表示Ωi,j,k中第k块影像块的权重，ωi,j,k通过基于局部线性约束求取：
min||M2(i,j)-M1(Ωi,j)ωi,j||2+λ||di,jeωi,j||2
s.t.1Tωi,j＝1
其中，ωi,j表示对应于M2(i,j)影像块的权重，||di,j e ωi,j||2为局部线性约束正则化项，di,j为相似像素距离M1(i,j)中心的欧几里得距离，λ表示正则项参数；对于此局部线性约束的优化问题的解为：
ωi,j＝((Ci,j+λdiag(di,j))\1)(1T(Ci,j+λdiag(di,j))\1)-1
其中Ci,j＝(M1(Ωi,j)-M2(i,j)T)(M1(Ωi,j)-1M2(i,j)T)T表示影像块的协方差矩阵，diag(di,j)表示以di,j为对角线元素的对角矩阵；之后从多光谱影像对求取的权重ωi,j将被共享给上采样后高光谱影像以生成预测时刻的上采样的高光谱影像
步骤3、获取初始化的目标融合影像时间变化影像，表示为：
步骤4、基于光谱线性混合模型，将T2时刻的目标影像X2分解为保留有光谱信息的端元矩阵E2和保留有目标影像空间信息的丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像。
步骤5、建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型。
步骤5中，多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型表示为：
Mi＝RXi＝REiAi
Hi＝XiFS＝EiAiFS
其中，i∈{1,2}，分别为在Ti时刻的高时空谱分辨率影像、高光谱影像以及多光谱影像；分别表示在Ti时刻的高时空谱分辨率影像的端元矩阵和丰度矩阵，Q为端元的个数，Q<<L；表示光谱降采样矩阵，和分别表现空间模糊矩阵和空间下采样矩阵。
步骤6、引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵。
步骤6中，引入广义线性混合模型后，目标影像X2的估计问题转为了估计目标影像的丰度矩阵A2以及目标影像的时间变化矩阵P。
步骤7、获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T。
步骤7中，T的初始化由得到；E1的初始化为利用端元提取算法从高光谱影像H1中提取得到；时间变化P的初始化为目标影像丰度A2的初始化由对多光谱影像M2分类或通过光谱解混得到。
步骤7包括：
步骤7.1、分别引入关于目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T的正则化项，获得融合问题的能量函数。
步骤7.1中，能量函数表示为：
s.t.D＝M2-M1
步骤7.1包括：
步骤7.1.1、引入关于目标影像的丰度矩阵A2的正则化项，该正则化项为水平和垂直方向的空间平滑约束，表示为：
其中，Hh和Hv是线性算子，分别用于计算二维信号中相邻像素的分量之间的水平和垂直梯度，||·||2,1为混合L2,1范数，λA是丰度矩阵A2的正则化参数；
步骤7.1.2、引入时间变化矩阵P的正则化项，该正则化项为光谱平滑约束，表示为：
其中，H1是光谱维方向的差分算子，λ1和λ2是时间变化P的正则化参数；
步骤7.1.3、引入时间变化影像T的正则化项，该正则化项为低秩约束，表示为：
其中，为T的三阶张量形式，||·||TNN为张量核范数，有其中表示第j个奇异值，通过对T进行傅里叶变换得到，λT为时间变化影像T的正则化参数。
在步骤7.1中，本申请将广义光谱线性混合模型和时间变化影像重建模型进行耦合，并引入关于丰度矩阵A2，目标影像的时间变化矩阵P，以及时间变化重建影像T的先验项构建了时空谱融合模型。
步骤7.2、分割能量函数为对应于丰度矩阵A2、时间变化矩阵P和时间变化影像T变量的三个子函数，并为每个子函数引入若干辅助变量分裂子函数为子问题，为每个子问题构建增广拉格朗日函数。
在步骤7.2中，本申请使用的分割策略为在固定其它两个变量情况下，单独的求解其中一个变量，并使用交替最小二乘方法来最小化对应于每一个变量的损失函数，以获得一个局部的驻点。
步骤7.2中，目标影像的丰度矩阵A2的子问题为：
引入辅助变量B1＝A2，B2＝Hh(A1)，B3＝Hv(A1)，令u＝vec(A2)，符号vec(·)为矢量化操作，给出丰度A2的子问题的增广拉格朗日函数为：
其中，矩阵Vi,i＝1,2,3为拉格朗日函数的对偶变量，ρ为步长，有ρ>0，repM(V)为矩阵V的块对角矩阵，表示将矩阵V沿对角线复制M次。
本申请将丰度A2拉格朗日函数的求解，分裂为关于变量u,B1,B2,B3,V1,V2,V3的子问题
关于丰度A2的变量u子问题为：
对上式的求解有：其中Ω1和Ω2由以下给出：
Ω1＝2(R(PeΕ1))TR(PeΕ1)+η(PeΕ1)T(PeΕ1)+ρI
Ω2＝(R(PeΕ1))T(M2+D+H1)+(PeΕ1)T(T+E1A1)+ρ(B1+V1)
关于丰度A2的变量B1子问题为：
对上式的求解有：
上式中符号vec(·)表示将矩阵转为向量的算子。
关于丰度A2的变量B2子问题为：
对上式的求解等同于L2范数的近端算子。该解可以通过使用soft阈值算子得到，表示为：
上式中，soft阈值运算符被定义为：
关于丰度A2的变量B3子问题为：
与B2子问题的优化策略相同,求解由soft阈值得到：
关于丰度A2变量的对偶变量的更新规则如下：
V1←V1+vec(B1)-u
V2←V2+vec(B2)-Hhvec(B1)
V3←V3+vec(B3)-Hvvec(B1)
目标影像的时间变化矩阵P的子问题为：
引入辅助变量B＝PeΕ1，时间变化矩阵P子问题的增广拉格朗日函数为：
本申请将时间变换P的增广拉格朗日函数的求解，分裂为了关于变量B,P,V的子问题。
关于时间变换P的变量B子问题为：
对上式的求解为：
(2RTR+1)BA2A2 T+ρB＝RT(M2+D+RE1A1)A2 T+(T+E1A1)A2 T+ρ(PeΕ1-V)
上式构成了西尔韦斯特方程，在matlab中使用命令X＝dlyap(A,B,C)来求解方程。
关于时间变换P的变量P子问题为：
上式改写为向量的形式进行求解:
关于时间变换P的对偶变量的更新规则如下：
V←V+B-PeΕ1
目标影像的时间变化影像T的子问题为：
通过引入辅助变量给出时间变换影像T子问题的增广拉格朗日函数为：
本申请还将时间变换影像T拉格朗日函数的求解，分裂为关于变量的子问题,其中分别为矩阵T,B,V的三阶张量形式。
关于时间变换影像T的变量子问题为：
将展开为矩阵，并令等于0，从而得到关于T求解的方程：
关于时间变换影像T的变量子问题为：
对上式的求解有：其中是奇异值收缩算子，UΣVT为奇异值分解过程。
关于时间变换T的对偶变量的更新规则如下：
步骤7.3、利用交替方向乘子法，对增广拉格朗日函数进行优化求解，获得目标影像的丰度矩阵A2，时间变化矩阵P和时间变化影像T。
步骤8、目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2。
实施例2：
一种生成高时空谱分辨率卫星遥感影像的装置，包括：
第一获取模块，用于获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
搜寻模块，用于基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到预测时刻的上采样后的高光谱影像
第二获取模块，用于获取初始化的目标融合影像时间变化影像，表示为：
分解模块，用于基于光谱线性混合模型，将T2时刻的目标影像X2分解为端元矩阵E2和丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像；
建立模块，用于建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型；
表示模块，用于引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵；
第三获取模块，用于获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T；
得到模块，用于令目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2。
综上所述，本申请提出了对影像对时间变化特征进行刻画和建模的时间变化模型，以及为实现时间变化保真优化求取提出了时间变化影像分辨率增强重建模型，并通过耦合时间变化模型和时间变化影像分辨率增强重建模型，建立了一种广义的时空谱融合框架。时间变化模型对影像间的时间变化特征进行建模和刻画，时间变化影像分辨率增强重建模型对时间变化进行保真约束求取。该方法能够实现低时间、空间分辨率高光谱影像与高时间、空间分辨率多光谱影像的时空谱融合，生成空间和光谱高保真的高时空谱分辨率融合影像。

Claims (2)

1.一种生成高时空谱分辨率卫星遥感影像的方法，其特征在于，包括：
步骤1、获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
步骤2、基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到T2时刻的上采样后的高光谱影像
步骤3、获取初始化的目标融合影像时间变化影像，表示为：
步骤4、基于光谱线性混合模型，将T2时刻的目标影像X2分解为端元矩阵E2和丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像；
步骤5、建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型；
步骤6、引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵；
步骤7、获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T；
步骤8、目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2；
步骤1包括：
步骤1.1、获取T1时刻的高光谱影像和多光谱影像以及T2时刻的多光谱影像其中，l、L为波段数量，w、W为像素数量，l<L，w<W；
步骤1.2、对高光谱影像H1和多光谱影M1、M2进行归一化，并对归一化后的高光谱影像H1上采样至多光谱影像空间大小，得到
步骤1.3、将多光谱影像M1、M2和上采样后的高光谱影像沿光谱维折叠为三维形式，并分割为宽和高皆为光谱维为l的全波段块；
步骤5中，多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型表示为：
Mi＝RXi＝REiAi
Hi＝XiFS＝EiAiFS
其中，分别为在Ti时刻的高时空谱分辨率影像、高光谱影像以及多光谱影像；分别表示在Ti时刻的高时空谱分辨率影像的端元矩阵和丰度矩阵，Q为端元的个数；表示光谱降采样矩阵，和分别表现空间模糊矩阵和空间下采样矩阵；步骤7包括：
步骤7.1、分别引入关于目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T的正则化项，获得融合问题的能量函数；
步骤7.2、分割能量函数为对应于丰度矩阵A2、时间变化矩阵P和时间变化影像T变量的三个子函数，并为每个子函数引入若干辅助变量分裂子函数为子问题，为每个子问题构建增广拉格朗日函数；
步骤7.3、利用交替方向乘子法，对增广拉格朗日函数进行优化求解，获得目标影像的丰度矩阵A2，时间变化矩阵P和时间变化影像T；
步骤7.1中，能量函数表示为：
s.t.D＝M2-M1
其中D为多光谱影像M2和M1的差分影像；
步骤7.1包括：
步骤7.1.1、引入关于目标影像的丰度矩阵A2的正则化项，该正则化项为空间平滑约束，表示为：
其中，Hh和Hv是线性算子，分别用于计算二维信号中相邻像素的分量之间的水平和垂直梯度，‖·‖2,1为混合L2,1范数，λA是丰度矩阵A2的正则化参数；
步骤7.1.2、引入时间变化矩阵P的正则化项，该正则化项为光谱平滑约束，表示为：
其中，H1是光谱维方向的差分算子,为L行Q列元素全为1的矩阵，λ1和λ2是时间变化P的正则化参数；
步骤7.1.3、引入时间变化影像T的正则化项，该正则化项为低秩约束，表示为：
其中，为T的三阶张量形式，‖·‖TNN为张量核范数，有其中表示第j个奇异值，通过对T进行傅里叶变换得到，λT为时间变化影像T的正则化参数；
步骤7.2中，目标影像的丰度矩阵A2的子问题为：
其中η为时间变化影像项参数，用来控制时间变化影像对模型的贡献；
引入辅助变量B1＝A2，B2＝Hh(A1)，B3＝Hv(A1)，令u＝vec(A2)，符号vec(·)为矢量化操作，给出丰度A2的子问题的增广拉格朗日函数为：
其中，矩阵Vi,i＝1,2,3为拉格朗日函数的对偶变量，ρ为步长，有ρ>0，repM(V)为矩阵V的块对角矩阵，表示将矩阵V沿对角线复制M次；
步骤7.2中，目标影像的时间变化矩阵P的子问题为：
引入辅助变量B＝Pe Ε1，时间变化矩阵P子问题的增广拉格朗日函数为：
目标影像的时间变化影像T的子问题为：
通过引入辅助变量给出时间变换影像T子问题的增广拉格朗日函数为：
步骤2中，使用多光谱影像M1相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示：其中，ωi,j,k表示Ωi,j,k中第k块影像块的权重，ωi,j,k通过基于局部线性约束求取：
min||M2(i,j)-M1(Ωi,j)ωi,j||2+λ||di,jeωi,j||2
s.t.1Tωi,j＝1
其中，ωi,j表示对应于M2(i,j)影像块的权重，||di,je ωi,j||2为局部线性约束正则化项，di,j为相似像素距离M1(i,j)中心的欧几里得距离，λ表示正则项参数；对于此局部线性约束的优化问题的解为：
ωi,j＝((Ci,j+λdiag(di,j))\1)(1T(Ci,j+λdiag(di,j))\1)-1
其中Ci,j＝(M1(Ωi,j)-M2(i,j)T)(M1(Ωi,j)-1M2(i,j)T)T表示影像块的协方差矩阵，diag(di,j)表示以di,j为对角线元素的对角矩阵；之后从多光谱影像对求取的权重ωi,j将被共享给上采样后高光谱影像以生成T2时刻的上采样的高光谱影像
2.一种生成高时空谱分辨率卫星遥感影像的装置，其特征在于，用于执行权利要求1所述的生成高时空谱分辨率卫星遥感影像的方法，包括：
第一获取模块，用于获取T1时刻的高光谱影像H1和多光谱影像M1，以及T2时刻的多光谱影像M2，并进行预处理，获取上采样后高光谱影像
搜寻模块，用于基于局部线性约束，搜寻多光谱影像M1中的相邻相似全波段块M1(Ωi,j,k)对多光谱影像M2中心位于(i,j)处的影像块进行表示，并求取相似全波段块的表示权重函数；将表示权重函数共享给上采样后高光谱影像得到T2时刻的上采样后的高光谱影像
第二获取模块，用于获取初始化的目标融合影像时间变化影像，表示为：
分解模块，用于基于光谱线性混合模型，将T2时刻的目标影像X2分解为端元矩阵E2和丰度矩阵A2，表示为：X2＝E2A2；目标影像X2为高时空谱分辨率影像；
建立模块，用于建立多光谱影像、高光谱影像与高时空谱分辨率影像的观测模型；
表示模块，用于引入广义线性混合模型，将端元矩阵E2表示为E2＝Pe E1，符号e表示逐元素乘的哈达玛积，P表示目标影像的时间变化矩阵，E1为T1时刻的端元矩阵；
第三获取模块，用于获取目标影像的丰度矩阵A2、时间变化矩阵P和时间变化影像T；
得到模块，用于令目标影像最终通过T1时刻的端元矩阵E1逐元素乘以T2时刻的时间变化矩阵P后再乘以高分辨率丰度矩阵A2得到：X＝(Pe E1)A2。
CN202310546379.5A 2023-05-16 2023-05-16 一种生成高时空谱分辨率卫星遥感影像的方法和装置 Active CN116757925B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202310546379.5A CN116757925B (zh) 2023-05-16 2023-05-16 一种生成高时空谱分辨率卫星遥感影像的方法和装置
PCT/CN2024/094956 WO2024208373A1 (zh) 2023-05-16 2024-05-23 一种生成高时空谱分辨率卫星遥感影像的方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310546379.5A CN116757925B (zh) 2023-05-16 2023-05-16 一种生成高时空谱分辨率卫星遥感影像的方法和装置

Publications (2)

Publication Number Publication Date
CN116757925A CN116757925A (zh) 2023-09-15
CN116757925B true CN116757925B (zh) 2024-07-02

Family

ID=87948580

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310546379.5A Active CN116757925B (zh) 2023-05-16 2023-05-16 一种生成高时空谱分辨率卫星遥感影像的方法和装置

Country Status (2)

Country Link
CN (1) CN116757925B (zh)
WO (1) WO2024208373A1 (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110148103A (zh) * 2019-04-29 2019-08-20 中国科学院西安光学精密机械研究所 基于联合优化的高光谱和多光谱图像融合方法、计算机可读存储介质、电子设备
CN112581365A (zh) * 2020-11-10 2021-03-30 北京拙河科技有限公司 一种跨尺度自适应信息映射成像方法及装置、介质

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7491944B1 (en) * 2005-07-14 2009-02-17 Sandia Corporation Method to analyze remotely sensed spectral data
US10685230B2 (en) * 2018-09-06 2020-06-16 National Central University Method of top-of-atmosphere reflectance-based spatiotemporal image fusion using aerosol optical depth
GB202001065D0 (en) * 2020-01-24 2020-03-11 Secr Defence Computer implemented method of spectral un-mixing
CN112396029B (zh) * 2020-12-03 2022-02-18 宁波大学 一种聚类分割与耦合端元提取协同的高光谱滨海湿地亚像元变化检测方法
CN113205453B (zh) * 2021-04-06 2022-03-15 武汉大学 一种基于空间-光谱全变分正则化的高光谱融合方法
CN113392790B (zh) * 2021-06-24 2022-03-18 哈尔滨工业大学 全色/多光谱遥感图像与高光谱遥感图像的融合方法
CN115272144A (zh) * 2022-06-22 2022-11-01 昆明理工大学 面向高光谱图像和多光谱图像的时空谱融合方法
CN115439344A (zh) * 2022-08-02 2022-12-06 武汉大学 联合双低秩和空谱全变差的混合噪声高光谱影像复原方法
CN115984155A (zh) * 2022-12-12 2023-04-18 武汉大学 一种基于光谱解混的高光谱、多光谱和全色图像融合方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110148103A (zh) * 2019-04-29 2019-08-20 中国科学院西安光学精密机械研究所 基于联合优化的高光谱和多光谱图像融合方法、计算机可读存储介质、电子设备
CN112581365A (zh) * 2020-11-10 2021-03-30 北京拙河科技有限公司 一种跨尺度自适应信息映射成像方法及装置、介质

Also Published As

Publication number Publication date
CN116757925A (zh) 2023-09-15
WO2024208373A1 (zh) 2024-10-10

Similar Documents

Publication Publication Date Title
CN111145131B (zh) 一种基于多尺度生成式对抗网络的红外和可见光图像融合方法
CN110119780B (zh) 基于生成对抗网络的高光谱图像超分辨重建方法
Dong et al. Generative dual-adversarial network with spectral fidelity and spatial enhancement for hyperspectral pansharpening
He et al. A new pansharpening method based on spatial and spectral sparsity priors
He et al. Deep convolutional neural network framework for subpixel mapping
CN111161199A (zh) 一种空谱融合的高光谱影像混合像元低秩稀疏分解方法
Ge et al. Improved semisupervised unet deep learning model for forest height mapping with satellite sar and optical data
Wei et al. Deep residual learning for remote sensed imagery pansharpening
CN114529830B (zh) 基于混合卷积网络遥感图像时空融合方法
Long et al. Dual self-attention Swin transformer for hyperspectral image super-resolution
CN112686830B (zh) 基于图像分解的单一深度图的超分辨率方法
Wang et al. Local–global feature-aware transformer based residual network for hyperspectral image denoising
CN118212127A (zh) 一种基于未配准的物理指导生成式对抗高光谱超分辨率方法
CN116309136A (zh) 一种基于sar先验知识引导的遥感影像云区重建方法
Wang et al. Multiresolution analysis pansharpening based on variation factor for multispectral and panchromatic images from different times
Karthick et al. Deep RegNet-150 Architecture for Single Image Super Resolution of Real-time Unpaired Image Data
Pereira et al. An end-to-end framework for low-resolution remote sensing semantic segmentation
Wang et al. Poissonian blurred hyperspectral imagery denoising based on variable splitting and penalty technique
Li et al. EMLM-Net: An Extended Multilinear Mixing Model-Inspired Dual-Stream Network for Unsupervised Nonlinear Hyperspectral Unmixing
CN116757925B (zh) 一种生成高时空谱分辨率卫星遥感影像的方法和装置
Picone et al. Spectro-spatial hyperspectral image reconstruction from interferometric acquisitions
CN116563103A (zh) 一种基于自适应神经网络的遥感图像时空融合方法
Cui et al. Meta-TR: meta-attention spatial compressive imaging network with swin transformer
Rani et al. An efficient block based feature level image fusion technique using wavelet transform and neural network
Cao et al. CFMB-T: A cross-frequency multi-branch transformer for low-quality infrared remote sensing image super-resolution

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant