CN105550993B - 基于多重变换域的超分辨率重建方法 - Google Patents

Abstract

本发明涉及一种基于多重变换域的超分辨率重建方法，属于图像超分辨率技术领域，包括下述内容：对K幅输入的低分辨率图像分别进行非相关的变换域分析；根据目标超分辨率图像和输入的低分辨率图像尺寸确定下采样矩阵D；通过对代价函数minλ₁‖W₁(z₁)‖₁+λ₂‖W₂(z₂)‖₂+……+λ_j‖W_j(z_j)‖₁在z＝z₁+z₂+…+z_j和‖DH_iM_iz‑g_i‖₂≤n_i i＝1,2,…K约束下采用多重变换域下L₁范数最小化的优化算法重建超分辨率图像。对比现有技术，本发明方法解决了在超分辨率重建中的先验概率模型的选取问题，即基于多重变换域的遥感图像稀疏性表达；同时克服了低分辨率遥感图像间由于载荷平台抖动和大气扰动带来的非规则的几何畸变和错位的问题，重建出具有良好效果的超分辨率图像。

Description

基于多重变换域的超分辨率重建方法

技术领域

本发明涉及一种超分辨率重建方法，特别涉及一种基于多重变换域的超分辨率重建方法，即利用多时相、同一场景的低分辨率图像采用后处理的方式重建出高分辨率图像的方法，属于图像超分辨率技术领域。

背景技术

高分辨率成像对遥感、军事、天文、生物医疗等诸多方面具有重要意义。人类对高分辨率图像的追求从未停止过，由于成像设备在理论上和客观条件的限制，成像的空间分辨率往往会遭遇瓶颈。例如在保持其它条件不变的情况下，最直接提高分辨率的方法就是提高CCD或CMOS的密度，即减小每个感光单元尺寸，然而当CCD或CMOS的感光单元小到一定程度时，图像的质量将开始下降，这是因为随着感光单元的减小,光通量也随之减少,在曝光过程中，每个感光单元所收集到的光子将被热噪声所掩盖。提高分辨率也可以通过相应加大光学系统设计和CCD/CMOS阵列来实现，然而这就必然带来成像设备体积和重量的增加，这在某些对成像设备体积、功耗、重量要求极为苛刻的情况下是不可接受的，例如卫星遥感成像领域。相比常规的成像设备，除了光学系统、CCD或CMOS阵列密度等因素外，还有一些特殊的因素决定卫星遥感成像分辨率，例如：观测幅宽、观测波段范围、卫星重返周期等。一般来说，星载成像仪的高时间分辨率和高空间分辨率是一对矛盾，在同样视场角的前提下，重访周期短的卫星，往往成像的空间分辨率较低。卫星一旦发射，星载相机的成像分辨率是固定的，但卫星以一定的周期绕地球旋转，会重复拍摄同一个区域，这就为超分辨率技术的应用提供了契机。超分辨率重建技术作为一种投资少、收益高的后处理技术，可以综合地利用现有的航天/航空遥感载荷和观测数据资源，从而满足资源探测、灾害监测、国防建设等行业对高空间分辨率遥感图像的要求。

从超分辨率重建技术被第一次提出已有近三十年的历史,在涉及公共安全的视频监控领域、医学成像领域和图像通信等诸多领域获得巨大的成功，例如：基于超分辨率技术可以把一般的NTSC格式或PAL格式电视信号转换为高清电视信号(简称HDTV)而不失真地在HDTV电视上播放。但在对地观测领域还远没有达到实用的程度,主要原因如下:常规超分辨率方法无法克服低分辨率遥感图像间由于载荷平台抖动和大气扰动带来的非规则的几何畸变和错位；缺乏针对遥感图像的高效重建方法；并且超分辨率重建技术在航天/航空遥感领域应用中会因客观条件而受限。由于一直以来获取多时相低分辨率遥感图像的时间跨度往往较大(例如美国陆地卫星LANDSAT-7的重访周期是16天),因此很难保证地物在这个时间段内没有发生变化。但是随着我国近年来航天/航空事业的迅猛发展，对地观测载荷和相应的观测数据呈爆炸的增长态势，极大地缩短了获取多时相低分辨率遥感图像的时间，例如，我国资源三号卫星的回归周期是59天，但是在特殊需要的情况下，采用侧摆功能，则它的重访周期可达5天；再比如我国的环境减灾卫星HJ-1A和HJ-1B各自的重访周期是4天，结合HJ双星，重访周期可达2天。这些在短时期内获取的在同一星下观测点的遥感影像，为此时开展超分辨率技术研究提供了前所未有的发展契机。

目前实现图像稀疏性表达的方法主要分成两大类：一种是通过训练或分析的方法建立稀疏字典，这类方法能够根据图像本身的特点自适应地选取稀疏字典中的元素来实现图像的稀疏性表达，但是基于分析方法建立的图像数据表达模型，自适应能力和可移植性较差；基于训练方法建立的稀疏字典，缺乏正则性、存在冗余、计算效率低。另一类方法则是基于多尺度分析理论又称为单一基方法，建立适合表达图像边缘和纹理几何结构的基来实现图像的稀疏性表达，代表方法是：小波变换、曲波变换、脊波变换等，这类方法的特点是简单且更具普适性，因而本专利立足于后者，挖掘遥感图像在多重变换域的稀疏性表达，进而采用MAP来解决超分辨率这样一个病态问题的求解。常规的超分辨率方法在重建遥感图像时效果往往不是很理想，主要原因是它们的先验概率模型并不适合遥感图像。以经典方法Huber-MAP为例，把Huber马尔科夫模型作为最大后验概率方法MAP中的先验概率模型，该模型的本质是对图像的二阶导数进行Huber边缘惩罚，它比通用的高斯先验概率模型能较好地保持图像边缘信息，但是在重建细节信息丰富的遥感图像时，很难平衡去除噪声和保持边缘的矛盾。本发明借鉴压缩感知理论和最大后验概率方法MAP，充分挖掘遥感图像在多重变换域的稀疏性表达，并把这种稀疏性表达作为MAP中的先验概率模型，形成以L₁范数最小化为核心的超分辨率重建代价函数。由于地表特征的复杂性，很难找到单一的通用变换域实现遥感图像的稀疏性表达。这里我们将一幅遥感图像理解成由包含不同纹理特征的子图构成，不同纹理特征的子图分别由不同的变换域来实现其稀疏性表达，进而通过组合多种不同的变换域即多重变换域，来实现整体遥感图像的稀疏性表达。结合L₁范数最小化的优化算法实现超分辨率的重建，采用L₁范数最小化主要是因为，它是一种高效保护目标信号稀疏性并且鲁棒特性较好的凸集优化方法。随着近年来压缩感知的迅猛发展，很多基于L₁范数最小化的优化方法被相继提出，这就为我们构建高效的超分辨重建方法提供了很好的技术保障。

发明内容

本发明的目的是针对现有技术存在的不足，提供一种基于多重变换域的超分辨率重建方法，该方法利用压缩感知的重建方法来解决超分辨率的重建问题，以克服现有遥感图像模糊、不精确和畸变等问题。

本发明的目的是通过以下技术方案实现的：

一种基于多重变换域的超分辨率重建方法，包括以下内容：

对K幅输入的低分辨率图像分别进行非相关的变换域分析；

根据目标超分辨率图像和输入的低分辨率图像尺寸确定下采样矩阵D；

通过对如下代价函数采用多重变换域下L₁范数最小化的优化算法重建超分辨率图像：

minλ₁||W₁(z₁)||₁+λ₂||W₂(z₂)||₁+...λ_j||W_j(z_j)||₁

minλ₁||W₁(z₁)||₁+λ₂||W₂(z₂)||₁+...λ_j||W_j(z_j)||₁

s.t.z＝z₁+z₂+...z_j

||DH_iM_iz-g_i||₂≤n_i i＝1,2,...K

其中，z表示原始图像，z₁、z₂……z_j分别表示超分辨率图像中包含不同纹理特征的j个子图，W₁、W₂……W_j分别表示针对不同特征子图的j个非相关多重变换域，λ₁、λ₂……λ_j分别表示对应子图的权重，||·||₁表示L₁范数，g_i表示输入的第i个低分辨率图像，M_i表示第i幅图像的位置错位和几何变形的矩阵，H_i表示第i幅图像的模糊降质过程的一个块循环矩阵，n_i表示第i幅图像的加性噪声，||·||₂表示对·取绝对值，s.t.表示约束条件，min f表示f的最小值。

作为优选，所述多重变换域下L₁范数最小化的优化算法在优化时采用收缩循环迭 代法，其收缩运算符为

有益效果

对比现有技术，本发明方法针对常规超分辨率方法无法克服低分辨率遥感图像间由于载荷平台抖动和大气扰动带来的非规则的几何畸变和错位的问题，基于遥感图像的特点，本发明借鉴压缩感知理论和最大后验概率方法MAP，充分挖掘遥感图像在多重变换域的稀疏性表达，并把这种稀疏性表达作为MAP中的先验概率模型，形成以L₁范数最小化为核心的超分辨率重建代价函数；并通过采用收缩循环迭代法对其求解从而获得效果良好的超分辨率图像。

本发明方法的原理如下：

超分辨率重建是一个典型的病态问题，针对这类病态问题，一般情况下，我们需要一个先验模型来指导重建过程，最终重建出原始的高分辨率图像。由于卫星遥感图像细节丰富、地表覆盖复杂，因而很难构建一种适用于所有遥感图像的单一稀疏表达方式，本专利提出一种结合多种变换域，来挖掘卫星遥感图像的稀疏性表达。一幅遥感图像可以理解为由几幅不同纹理特征的子图构成。不同纹理特征的子图在不同的变换域体现出不同的稀疏性，通过结合不同的变换域，来挖掘遥感图像的稀疏性表达；再把这种多重变换域稀疏性表达作为先验模型引入到高分辨率图像的重建过程。

这里的变换域分析主要是用来挖掘信号的稀疏性表达。正如小波变换在图像压缩标准JPEG2000中获得极大的成功，主要是因为把常规图像变换到小波域中，大多数的小波系数由零或接近零的系数构成，即常规图像在小波域具有稀疏性表达。每种变换域具有不同的特性，适合不同信号的稀疏性表达。例如小波变换比较适合平滑图像中边缘的稀疏性表达，总变差适合分段线性图像的稀疏性表达，曲波变换比较适合卡通图片和文字的稀疏性表达，每种变换都有其针对的特定图像。由于遥感图像细节丰富、地表覆盖复杂，现实中很难找到一种通用的变换域适合所有遥感图像的稀疏性表达。在组合变换域的研究中，需要考虑到不同变换域的非相关性的问题，这里的“非相关”的概念主要是参照压缩感知理论中“Incoherence”，例如傅里叶变换域和空域本身是非相关的，小波变换域和噪波变换(Noiselet)域是非相关的。这里强调非相关的概念主要是避免遥感图像在多重变换域中，稀疏性被重复投影的问题。

由于超分辨率重建本质上是基于退化模型反演出原始图像，因而分析原始遥感图像的退化模型是一个不可回避的问题。成像过程中的降质因素一般包括：平台漂移和大气扰动带来的局部挤压、运动模糊和光学系统带来的模糊、下采样和噪声污染等不良影响。拟通过下列的数学表达来描述原始遥感图像的退化模型：

g_i＝DH_iM_iz+n_i,i＝1,2,...,K

其中，z是原始高分辨率图像；g_i是观测到的第i个低分辨率图像；K是低分辨率图像的数量；M_i是表示第i幅图像的位置错位和几何变形的矩阵；H_i是一个块循环矩阵用来表示模糊降质过程，它可以包括失焦、运动模糊或光学系统传递函数带来的模糊等；D是一个下采样矩阵；n_i表示第i幅图像的加性噪声。原始遥感图像的退化过程和参数描述的示意图如图1所示。超分辨率的目标是使用已知系列观测到的低分辨率图像，重建出原始高分辨率图像z。退化模型的建立是后续反演出高分辨率图像的基础。

本发明利用压缩感知的重建方法即L₁范数最小化的方法来解决超分辨率的重建问题。假设已知遥感图像在多重变换域中具有稀疏性即经多重变换域变换后具有最少的非零系数，因而可以采用L0范数最小化来求出最优解，很不幸的是L0范数最小化求解过程是一个NP-Hard问题。根据压缩感知理论，当测量矩阵满足约束等距特性条件下，L0范数最小化可以采用L₁范数最小化来代替。利用L₁范数最小化来求解可以追溯到当年Beurling基于部分傅里叶变换系数推测原始目标信号的时代。从上个世纪70年代到80年代开始，随着现代计算机技术的发展，计算能力呈爆炸式的增长，因而基于L₁范数最小化被广泛地应用于解决实际的问题。特别是在近几年，随着压缩感知理论的蓬勃发展，再一次引起了专家学者的广泛关注。正是由于L₁范数最小化方法具有保持原始信号的稀疏性和优良的鲁棒特性，这里拟采用把L₁范数最小化和最大后验概率(MAP)相结合的方法，即把遥感图像在多重变换域的稀疏性表达当作先验概率模型，再通过正则化的方法来完成超分辨率的重建。针对求解二次型问题时，我们采用收缩循环迭代法，也称为软门限法(Soft-Threshloding)，最早由David Donoho提出，它是一种经典的基于小波变换域的图像去噪声方法。其收缩运算符定义为：

采用的代价函数为：

minλ₁||W₁(z₁)||₁+λ₂||W₂(z₂)||₁+...λ_j||W_j(z_j)||₁

s.t.z＝z₁+z₂+...z_j

||DH_iM_iz-g_i||₂≤n_i i＝1,2,...K

其中，z是原始图像，z₁、z₂和z_j分别是遥感图像中包含不同纹理特征的j个子图，W₁、W₂和W_j分别表示针对不同特征子图的j个非相关多重变换域，λi表示第i个子图的权重。

附图说明

图1是卫星遥感图像在成像过程中的退化模型示例。

图2是本发明实施例遥感图像超分辨率重建流程示意图。

图3是多重变换域算法流程示意图。

图4是重建稀疏成分流程示意图。

图5是重建平滑成分流程示意图。

图6是重建曲波成分流程示意图。

具体实施方式

现将详细参照本发明的方法结合附图对本发明的实施例作进一步详细说明。如图2所示，本发明的超分辨率重建方法的步骤如下：

第一步，获取需要超分辨重建的多张低分辨率图g₁,g₂,...,g_k，并选定一幅低分辨率图像g_ref作为参考图像。

第二步，通过弹性配准的方法得到配准映射矩阵M₁,M₂,...,M_k，该弹性配准方法克服了载荷平台抖动和大气扰动带来的非规则的几何畸变和错位。

所述弹性配准可以采用目前任何一种能够克服非规则的几何畸变和错位的方法，如S.Periaswamy和H.Farid于2003年提出的一种弹性配准方法，为了达到提高配准效率的目的，本实施例使用的弹性配准方案如下：针对图像间包含非规则几何畸变的情形，虽然不能用全局的数学表达式来描述这种映射关系，但是我们可以粗略地理解为图像间包含刚性的全局表达式再加上图像局部区域各自不同的映射参数。待配准图像和参考图像首先通过刚性配准方法进行粗配，而后再通过分块的策略来精细地微调，逐像素滑动块的窗口，计算出每个块中心像素位置的微调参数。

第三步，根据低分辨率图g₁,g₂,...,g_k估算模糊矩阵H₁,H₂,...,H_k。

假定模糊图像模糊模型为高斯类型，可以采用根据图像中点目标或线状目标被模糊后的像素个数来估算出一个高斯卷积核H。

第四步，根据超分辨率的放大倍数算出下采样矩阵D。

对于一幅尺寸为M*N的图像I，对其进行s倍下采样，即得到(M/s)*(N/s)尺寸的分辨率图像，s是M和N的公约数。这个下采样像素点的值就是窗口内所有像素的均值。

第五步，对g₁,g₂,...,g_k分别估算出每幅低分辨率图像的噪声方差n₁,n₂,...,n_k。

可以利用第一级小波变换后的系数来估算出，详见David Donoho等人的去噪文章。

第六步，结合配准映射矩阵、模糊矩阵和噪声方差，通过多重变换域下的L₁范数的最小优化算法重建出以g_ref作为参考图像的超分辨率图像z。

通过多重变换域下的L₁范数的最小优化算法重建分辨率图像z可以采用多种方式实现，本实施例以3重变换域为例通过以下过程重建z：

首先，初始化参数变量，对每个低分辨率图像g_i重建其稀疏成分z_s；再对每个低分辨率图像g_i重建其平滑成分z_w；最后对每个低分辨率图像g_i，重建其边缘结构z_c。如图3所示即为本步骤所述的多重变换域算法流程示意图，具体流程如下：

1)初始化：

a)先选择软门限值τ(绝大多数情况下可以设置为1)和停止门限值δ(根据噪声强度设置)

b)计算总迭代数量l＝floor(τ/δ)

c)设置权重λ₁...λ₃

d)定义矩阵z_s、z_w和z_c为0矩阵z_s＝0；z_w＝0；z_c＝0；z＝z_s+z_w+z_c

2)循环计算l次：

a)对每个低分辨率图像g_i,重建其稀疏成分z_s，如图4所示即为本步骤重建稀疏成分的流程示意图，具体过程如下：

i.计算残差：

ii.计算梯度：

iii.更新z_s＝z_s-d

iv.当矩阵z_s中的元素的绝对值小于等于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z

b)对每个低分辨率图像g_i，重建其平滑成分z_w，如图5所示即为本步骤重建平滑成分的流程示意图，具体过程如下：

i.计算残差：

ii.计算梯度：

iii.z_w＝z_w-d

iv.计算z_w的DWT(离散小波)系数Wz_w

v.当矩阵Wz_w中的元素的绝对值小于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z

c)对每个低分辨率图像g_i，重建其边缘结构z_c，如图6所示即为本步骤重建曲波成分的流程示意图，具体过程如下：

i.计算残差：

ii.计算梯度：

iii.z_c＝z_c-d

iv.计算z_c的曲波变换系数Cz_c

v.当矩阵Cz_c z_c中的元素的绝对值小于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z

d)τ＝τ-δ

3)重建超分辨率图像z＝z_s+z_w+z_c

其中，·^*表示·的共轭转置矩阵。本领域技术人员通过上述重建过程可知，当多重变换域不只3重时可以同理为之，此处不再赘述。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于多重变换域的超分辨率重建方法，其特征在于，包括以下内容：

对K幅输入的低分辨率图像分别进行非相关的变换域分析；

根据目标超分辨率图像和输入的低分辨率图像尺寸确定下采样矩阵D；

通过对如下代价函数采用多重变换域下L₁范数最小化的优化算法重建超分辨率图像：

minλ₁||W₁(z₁)||₁+λ₂||W₂(z₂)||₁+...λ_j||W_j(z_j)||₁

s.t.z＝z₁+z₂+...z_j

||DH_iM_iz-g_i||₂≤n_i i＝1,2,...K

其中，z表示原始图像，z₁、z₂……z_j分别表示超分辨率图像中包含不同纹理特征的j个子图，W₁、W₂……W_j分别表示针对不同特征子图的j个非相关多重变换域，λ₁、λ₂……λ_j分别表示对应子图的权重，||·||₁表示L₁范数，g_i表示输入的第i个低分辨率图像，M_i表示第i幅图像的位置错位和几何变形的矩阵，H_i表示第i幅图像的模糊降质过程的一个块循环矩阵，n_i表示第i幅图像的加性噪声，||·||₂表示对·取绝对值，s.t.表示约束条件，min f表示f的最小值。

2.根据权利要求1所述的一种基于多重变换域的超分辨率重建方法，其特征在于，所述多重变换域下L₁范数最小化的优化算法在优化时采用收缩循环迭代法，其收缩运算符为

3.根据权利要求1所述的一种基于多重变换域的超分辨率重建方法，其特征在于，包括以下步骤：

第一步，获取需要超分辨重建的多张低分辨率图g₁,g₂,...,g_k，并选定一幅低分辨率图像g_ref作为参考图像；

第二步，通过弹性配准的方法得到配准映射矩阵M₁,M₂,...,M_k，该配准方法克服了载荷平台抖动和大气扰动带来的非规则的几何畸变和错位；

第三步，根据低分辨率图g₁,g₂,...,g_k估算模糊矩阵H₁,H₂,...,H_k；

第四步，根据超分辨率的放大倍数算出下采样矩阵D；

第五步，对g₁,g₂,...,g_k分别估算出每幅低分辨率图像的噪声方差n₁,n₂,...,n_k；

第六步，结合配准映射矩阵、模糊矩阵和噪声方差，通过多重变换域下的L₁范数的最小优化算法重建出以g_ref作为参考图像的超分辨率图像z。

4.根据权利要求3所述的一种基于多重变换域的超分辨率重建方法，其特征在于，所述弹性配准方法如下：待配准图像g_i和参考图像g_ref首先通过刚性配准方法进行粗配，而后再通过分块的策略来精细地微调，逐像素滑动块的窗口，计算出每个块中心像素位置的微调参数。

5.根据权利要求3-4任一所述的一种基于多重变换域的超分辨率重建方法，其特征在于，步骤六所述超分辨率图像z通过以下过程重建：

1)初始化如下参数：

a)设置门限值τ和停止门限值δ；

b)计算总迭代数量l＝floor(τ/δ)；

c)设置权重λ₁...λ₃；

d)定义矩阵z_s、z_w和z_c为0矩阵，z_s＝0；z_w＝0；z_c＝0；z＝z_s+z_w+z_c；

2)对下述过程循环计算l次：

a)对每个低分辨率图像g_i,通过以下步骤重建其稀疏成分z_s：

i.计算残差：

ii.计算梯度：其中·^*表示·的共轭转置矩阵；

iii.更新z_s＝z_s-d；

iv.当矩阵z_s中的元素的绝对值小于等于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z；

b)对每个低分辨率图像g_i，通过以下步骤重建其平滑成分z_w：

i.计算残差：

ii.计算梯度：

iii.z_w＝z_w-d；

iv.计算z_w的DWT系数Wz_w；

v.当矩阵Wz_w中的元素的绝对值小于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z；

c)对每个低分辨率图像g_i，通过以下步骤重建其边缘成分z_c：

i.计算残差：

ii.计算梯度：

iii.z_c＝z_c-d；

iv.计算z_c的曲波变换系数Cz_c；

v.当矩阵Cz_c中的元素的绝对值小于时，设置对应的元素值为0，根据公式z＝z_s+z_w+z_c更新z；

d)τ＝τ-δ；

3)重建超分辨率图像z＝z_s+z_w+z_c。