一种边缘先验引导的单帧遥感图像超分辨率处理方法

摘要

一种边缘先验引导的单帧遥感图像超分辨率处理方法，将待处理图像的边缘细节作为先验，对经过插值生成的中间图像进行高频信息恢复，以避免细节丢失或模糊。该方法基于正则优化理论，提供了一种单帧图像超分辨率处理机制，同时能够有效地提升图像细节描述能力、改善影像质量，对数据量大的遥感图像并行能力好，执行效率高。

说明书

技术领域

本发明属于图像处理领域，涉及一种单帧遥感图像的超分辨率处理方法。

背景技术

分辨率是衡量光学卫星影像成像质量的一个重要评价指标，遥感图像中包 含的细节信息与分辨率密切相关，其直接影响着对目标的分析、识别和定位能 力。由于像元尺寸的限制，再加上高分辨率相机体积庞大、造价昂贵，往往只 能获取到低分辨率图像。因此在已有成像环境(无法改变成像系统硬件条件) 和数据源条件下，如何根据信息理论，采用数字图像处理技术来提高遥感图像 的分辨率，是卫星地面分析与处理的重要环节之一。

常用的超分辨率方法归纳起来可以分为两类：一类是多帧超分辨率方法， 即获取相对运动关系已知的图像，利用冗余的采样信息构建重叠区域的高分辨 率细节；另一类为单帧超分辨率方法，即只有目标区域的一幅影像。大多数情 况下，卫星成像系统不支持微动模式，无法得到多帧图像，或者多帧间的相对 运动关系难以求解，因此单帧超分辨率技术的突破更为迫切。

目前主要的单帧超分辨率方法包括：插值、重构和统计学习的方法。

典型的插值方法包括：最近邻插值、线性插值、双三次插值(Keys，Hou) 等。此类方法通常模糊了高频信息，尤其是边缘细节，且会产生严重的块状效 应，在遥感应用中难以忍受。但其计算速度快、实时性好，一般作为复杂超分 方法的基础。为了克服传统方法的不足，Li，Allebach等提出了基于最小二乘 法的边缘内插算法，对插值的边缘有所改善，但在大幅遥感图像应用上无法克 服分块间的条纹效应。

基于重构的方法通常需要操作者给出某种先验，对超分辨率的结果与原始 图像进行先验一致投影，实用性较高。这类算法主要有迭代反投影法、凸集投 影法。迭代反投影法由Irani等人提出并引入到超分辨率问题中，该方法将反投 影误差均匀地累加到重建图像上，导致重建图像在边缘存在锯齿效应且边缘信 息模糊。Stark等人基于集合理论提出了凸集投影法，试图将解投影到约束集上。 Morse等人利用水平集实现了重构。这些算法数学描述直观，但是约束集往往 难以定义，而且对于整幅卫星图像的处理，计算速度需改善。

基于学习的方法需要构造低分辨率和高分辨率图像样本库，通过学习样本 库得到两者内在的联系，完成重建。早期，Chang等人提出了邻域嵌入算法， 使基于学习的算法得到了很大的关注。Yang等人利用稀疏表示中过完备字典的 思想，取得了很好的效果。但需要注意的是，任何一种学习方法都会严重依赖 学习的样本，且计算量庞大。

发明内容

本发明解决的技术问题是：克服现有单帧超分辨率方法应用于遥感图像时 复杂度高、块状效应及地物边缘模糊等不足，针对大幅遥感图像，提供了一种 基于边缘先验约束的单帧超分辨率并行处理方法，融合梯度场转移和迭代后投 影处理方法，使得最终重建的图像边缘细节得以保护，同时信噪比等质量指标 未大幅下降且处理速度快。

本发明的技术解决方案是：一种边缘先验引导的单帧遥感图像超分辨率处 理方法，包括如下步骤：

(1)获取原始图像I_L，并对I_L进行分块，对I_L的分块图像进行索引并记为 SubI_L；

(2)根据超分辨率倍数m，对每一个SubI_L进行插值，得到分块图像 SubI_Htmp，SubI_Htmp的集合构成过程图像I_Htemp；

(3)针对任意SubI_L及其对应插值后的SubI_Htmp，均执行以下的操作：

(3.1)按超分辨率倍数m，形成模糊滤波器h，利用h对SubI_L和SubI_Htmp分别进行卷积，并在卷积的基础上分别计算SubI_L的行方向边缘图SubL_ex，列方 向边缘图SubL_ey，以及SubI_Htmp的行方向边缘图SubHtmp_ex，列方向边缘图 SubHtmp_ey；

(3.2)对SubL_ex按像素求绝对值，并按最大最小值进行线性拉伸得到 线性边缘图|SubL_ex|_lin，

${|\mathrm{SubL}}_{\mathrm{ex}}{|}_{\mathrm{lin}}=\frac{\left|{\mathrm{SubL}}_{\mathrm{ex}}\right|-{\mathrm{Min}}_{\mathrm{ex}}}{{\mathrm{Max}}_{\mathrm{ex}}-{\mathrm{Min}}_{\mathrm{ex}}}$

其中，像素绝对值最大值Max_ex＝Max(|SubL_ex|)，像素绝对值最小值 Min_ex＝Min(|SubL_ex|)；然后对SubL_ey、SubHtmp_ex、SubHtmp_ey执行相同的 操作，得到对应的线性边缘图|SubL_ey|_lin、|SubHtmp_ex|_lin和|SubHtmp_ey|_lin；

(3.3)分别计算线性边缘图|SubL_ex|_lin，|SubL_ey|_lin、|SubHtmp_ex|_lin和， |SubHtmp_ey|_lin的方差L_vx、L_vy、Htmp_vx、Htmp_vy，并求取梯度场转移参数L_var与H_var，其中：

L_vx＝std(|SubL_ex|_lin)

L_vy＝std(|SubL_ey|_lin)

Htmp_vx＝std(|SubHtmp_ex|_lin)

Htmp_vy＝std(|SubHtmp_ey|_lin)

L_var＝0.5(L_vx+L_vy)

H_var＝0.5(Htmp_vx+Htmp_vy)

其中std表示方差计算；

(3.4)采用一阶微分计算SubI_Htmp的行方向梯度场Htmp_Gx和列方向梯 度场Htmp_Gy，完成梯度场转移，其中：

${\mathrm{Htmp}}_{\mathrm{Gx}}={\mathrm{SubI}}_{H{\mathrm{tmp}}_x}-{\mathrm{SubI}}_{H{\mathrm{tmp}}_{x-1}}$

${\mathrm{Htmp}}_{\mathrm{Gy}}={\mathrm{SubI}}_{H{\mathrm{tmp}}_y}-{\mathrm{SubI}}_{H{\mathrm{tmp}}_{y-1}}$

梯度场转移后

${\mathrm{Htmp}}_{\mathrm{Gx}}={\mathrm{Htmp}}_{\mathrm{Gx}}*\mathrm{varf}*e^{\frac{-{\mathrm{dist}}_x^2}{L_{\mathrm{var}}}}/e^{\frac{-{\mathrm{dist}}_x^2}{H_{\mathrm{var}}}}$

$\mathrm{varf}=\frac{H_{\mathrm{var}}}{L_{\mathrm{var}}},{\mathrm{dist}}_x=0.5\left(\right|{\mathrm{SubHtmp}}_{{\mathrm{ex}}_x}{|}_{\mathrm{lin}}+\left|{\mathrm{SubHtmp}}_{{\mathrm{ex}}_{x-1}}{|}_{\mathrm{lin}}\right)$

${\mathrm{Htmp}}_{\mathrm{Gy}}={\mathrm{Htmp}}_{\mathrm{Gy}}*\mathrm{varf}*e^{\frac{-{\mathrm{dist}}_y^2}{L_{\mathrm{var}}}}/e^{\frac{-{\mathrm{dist}}_y^2}{H_{\mathrm{var}}}}$

${\mathrm{dist}}_y=0.5\left(\right|{\mathrm{SubHtmp}}_{{\mathrm{ex}}_y}{|}_{\mathrm{lin}}+\left|{\mathrm{SubHtmp}}_{{\mathrm{ey}}_{y-1}}{|}_{\mathrm{lin}}\right)$

式中角标x和x-1表示行方向的相邻两行，角标y和y-1表示列方向的相邻两 列；

(3.5)根据转移后的梯度场，利用稀疏矩阵求解器解算标准Poisson 图像方程，得到新的高分辨率图像子块SubI_C；

(3.6)采用迭代后向投影对SubI_C与SubI_L进行迭代投影，得到超分结 果图像子块记为SubI_H；

(4)对步骤(3)获得的全部SubI_H进行图像拼接，得到I_L的超分辨率结 果图像I_H。

所述步骤(3.1)中边缘图的计算采用归一化Sobel算子，算子的像素长和 宽均取[2R(m+1)+1]，R为取整操作。

所述步骤(3.6)中进行迭代投影的方法为：

SubI_C(i+1)＝SubI_C(i)*(1+β*RMS)，i＝1，2，…，N_s

迭代执行N_s步，其中$\mathrm{RMS}={\left|\right|{\mathrm{SubI}}_L-{d\left({\mathrm{SubI}}_{C\left(i\right)\otimes }h\right)}_m)\left|\right|}^2,$表卷积操作，d(I)_m表示将图像I按m倍进行降采样，β为拉格朗日平衡参数。β＝0.2。

本发明与现有技术相比的优点在于：

(1)本发明方法综合了插值与重构的优势，通过引入边缘与梯度场先验， 确保了超分辨出的结果图像地物边缘等高频信息保持性好；

(2)本发明方法采用了迭代后投影的处理，避免了传统重构方法和学习方 法中的块状效应；

(3)本发明方法并行执行效率高，计算复杂度小，同时信噪比等质量指标 未大幅下降，特别适用于单帧大幅遥感图像。

附图说明

图1为本发明方法的流程框图；

图2为本发明大幅遥感图像分块并行处理示意图。

具体实施方式

从遥感图像单帧超分辨趋势来看，保持边缘和细节的重构方法在快速性和 有效性上有巨大的优势。

边缘是指图像中像素灰度有阶跃变化的像素的集合，它存在于目标与背景、 目标与目标、区域与区域之间，并与图像亮度或图像亮度的一阶导数(梯度场) 的不连续性有关，从而表现为阶跃边缘和线条边缘。由于遥感图像中的边缘信 息往往对应着真实的地物，如道路、河流等。因此本发明基于重构的单帧超分 辨率处理方法的先验值将采用待处理图像的梯度场和边缘信息。

在此需要说明的是，本发明中所述的图像均为[0，1]空间浮点数据。

记原始低分辨率图像为I_L，分辨率提升倍数为m，本发明的超分辨率结果 图像为以下正则方程(最大似然问题)的解I_H：

即在已知低分辨率图像I_L及高分辨率图像的梯度场或边缘信息▽I_H情况下， 获取后验概率上最大的高分辨率图像结果。

上式可进一步采用拉格朗日方程表示为：

$I_H={\arg \min }_{I_{H*}}{\left|\right|I_L-d{(I_H^{*}\otimes h)}_m\left|\right|}^2+\beta {\left|\right|▿I_H-▿I_H^{*}\left|\right|}^2---\left(2\right)$

其中，表示的是I_H对I_L的后向投影误差，表示卷积操 作，h为模糊滤波器，为自变量，||·||²为标准二范数，表示梯度场一致性误差，d(u)_m表示将图像u按m倍进行降采样，β为拉格朗 日平衡参数，β∈(0，1)。即在已知I_L的情况下，获取高分辨率图像的梯度场 先验，并按后向投影和梯度场转移误差寻找迭代结果。

按目前使用的超分辨率使用方法，取I_L的插值结果作为I_H的最初解,鉴于 双三次插值的性能较高，本发明推荐采用双三次插值，以降低后续迭代处理的 复杂度。按式(1)和(2)的描述需要找到边缘和梯度场信息完成引导。结合 目前处理系统的多线程特性，可以将原始图像进行分块处理。通过这一系列操 作，既完整重现了图像的边缘细节且图像块之间不出现块状条纹，又使得到的 结果子块图像间也过度自然。

以下结合附图详细说明本发明的技术方案，具体流程如图1所示。

1)式(1)中的模糊滤波器h与具体的成像环境相关，通常考虑到图像的 带限及平滑特性，普遍使用归一化二维高斯卷积模板。本发明不限定对h的选取， 可以是任意的归一化二维模板，但其他的二维模板也容易导致边缘过渡锐化。

本发明实施例根据超分辨率倍数m，采用[2R(m+1)+1]×[2R(m+1)+1] 的归一化高斯模板h，其中R为取整操作。

2)根据超分辨率倍数m，计算超分辨率结果图像的大小，对I_L进行双三次 插值得到图像高分辨率过程图像I_Htemp。

设原始图像I_L中，以像素表示的图像宽和高分别为I_W和I_H，则超分辨率 结果图像I_H的像素宽和高分别为R(L_W*m)和R(L_H*m)，其中R为取整操作。

由于遥感图像幅宽很大，常常超过10000*10000像元，不能一次性整体的 进行内存分配，因此按附图2所示的分块方法并行处理。

取subw＝subh＝500(其中各图像分块在宽度方向和高度方向的长度分别 记为subw和subh)对图像进行分块，按如下多线程处理得到图像I_L双三次插值 后的高分辨率过程图像I_Htemp：

a)如果R(L_W*m)能被subw整除且R(L_H*m)能被subh整除，则将图像整 数分块，并进行双三次插值，结果直接写回；

b)如果R(L_W*m)不能被subw整除或R(L_H*m)不能被subh整除，则按能 整除数将图像整数分块，并进行双三次插值。不能整除的宽度部分向左扩充至 subw，高度部分向上扩充至subh处理，然后进行双三次插值，结果分步写回。

以下处理步骤针对的为每一个图像分块以及可能的扩展块，处理时可以独 立进行，最终线程合并时写出到最终结果文件。

3)对于原始图像I_L分块SubI_L和过程图像I_Htmp分块SubI_Htmp，与步骤1) 中的模板h卷积，以减小噪声的影响。分别计算卷积结果行方向和列方向边缘图。

记SubI_L的行方向边缘图记为SubL_ex，列方向边缘图记为SubL_ey。相应地， SubI_Htmp的为SubHtmp_ex、SubHtmp_ey。边缘图的计算采用归一化Sobel算子， 算子的像素长和宽对应地取为[2R(m+1)+1]。

4)对SubL_ex、SubL_ey和SubHtmp_ex、SubHtmp_ey按像素求绝对值，并按最 大最小值进行线性拉伸得到线性边缘图|SubL_ex|_lin，|SubL_ey|_lin和 |SubHtmp_ex|_lin，|SubHtmp_ey|_lin。

以|SubL_ex|_lin计算为例，像素绝对值最大值Max_ex＝Max(|SubL_ex|)，像素 绝对值最小值Min_ex＝Min(|SubL_ex|)，则：

${|\mathrm{SubL}}_{\mathrm{ex}}{|}_{\mathrm{lin}}=\frac{{|\mathrm{SubL}}_{\mathrm{ex}}|-{\mathrm{Min}}_{\mathrm{ex}}}{{\mathrm{Max}}_{\mathrm{ex}}-{\mathrm{Min}}_{\mathrm{ex}}}$

同理得到线性边缘图|SubL_ey|_lin和|SubHtmp_ex|_lin，|SubHtmp_ey|_lin。

5)计算线性边缘图|SubL_ex|_lin，|SubL_ey|_lin和|SubHtmp_ex|_lin，|SubHtmp_ey|_lin的方差L_vx、L_vy和Htmp_vx、Htmp_vy，求取梯度场转移参数L_var与H_var。

具体地：

L_vx＝std(SubL_ex|_lin)

L_vy＝std(|SubL_ey|_lin)

Htmp_vx＝std(|SubHtmp_ex|_lin)

Htmp_vy＝std(|SubHtmp_ey|_lin)

其中std表示方差计算。取方差L_vx、L_vy的均值L_var＝0.5(L_vx+L_vy)，Htmp_vx、 Htmp_vy的均值H_var＝0.5(Htmp_vx+Htmp_vy)作为步骤6)的输入。

6)采用一阶微分计算SubI_Htmp的行方向梯度场Htmp_Gx和列方向梯度场 Htmp_Gy，完成梯度场转移。

${\mathrm{Htmp}}_{\mathrm{Gx}}={\mathrm{SubI}}_{H{\mathrm{tmp}}_x}-{\mathrm{SubI}}_{{\mathrm{Htmp}}_{x-1}},$即相邻列差分；

${\mathrm{Htmp}}_{\mathrm{Gy}}={\mathrm{SubI}}_{H{\mathrm{tmp}}_y}-{\mathrm{SubI}}_{{\mathrm{Htmp}}_{y-1}},$即相邻行差分。

按以下方法完成梯度场转移，完成式(2)中的分项：

a)行方向，${\mathrm{Htmp}}_{\mathrm{Gx}}={\mathrm{Htmp}}_{\mathrm{Gx}}*\mathrm{varf}*e^{\frac{-{\mathrm{dist}}_x^2}{L_{\mathrm{var}}}}/e^{\frac{-{\mathrm{dist}}_x^2}{H_{\mathrm{var}}}},\mathrm{varf}=\frac{H_{\mathrm{var}}}{L_{\mathrm{var}}},$${\mathrm{dist}}_x=0.5\left(\right|{\mathrm{SubHtmp}}_{{\mathrm{ex}}_x}{|}_{\mathrm{lin}}+\left|{\mathrm{SubHtmp}}_{{\mathrm{ex}}_{x-1}}{|}_{\mathrm{lin}}\right),$即线性边缘图相邻列像素的 平均值；

b)列方向，${\mathrm{Htmp}}_{\mathrm{Gy}}={\mathrm{Htmp}}_{\mathrm{Gy}}*\mathrm{varf}*e^{\frac{-{\mathrm{dist}}_y^2}{L_{\mathrm{var}}}}/e^{\frac{-{\mathrm{dist}}_y^2}{H_{\mathrm{var}}}},\mathrm{varf}=\frac{H_{\mathrm{var}}}{L_{\mathrm{var}}},$${\mathrm{dist}}_y=0.5\left(\right|{\mathrm{SubHtmp}}_{{\mathrm{ey}}_y}{|}_{\mathrm{lin}}+\left|{\mathrm{SubHtmp}}_{{\mathrm{ey}}_{y-1}}{|}_{\mathrm{lin}}\right),$即线性边缘图相邻行像素的 平均值。

上述步骤就是根据低分辨图像与高分辨图像间的边缘及梯度场差异实现边 缘先验的转移，从而改善插值带来的边缘模糊及高频丢失。

7)根据转移后的梯度场，利用稀疏矩阵求解器解算标准Poisson图像方 程，得到新的高分辨率图像子块SubI_C。

Poisson方程的构建见Prez等文献，稀疏矩阵求解器可采用开源实现，如 Eigen。

此处理将在保证原始图像与目标图像梯度信息的同时，可以很好的融合背 景信息，避免分块处理时出现块状线。

8)为进一步避免大幅遥感图像分块处理后的块状效应，采用迭代后向投 影对SubI_C与SubI_L进行迭代投影。

迭代执行N_s步，

SubI_C(i+1)＝SubI_C(i)*(1+β*RMS)，i＝1，2，…，N_s

其中$\mathrm{RMS}={\left|\right|{\mathrm{SubI}}_L-{d\left({\mathrm{SubI}}_{C\left(i\right)\otimes }h\right)}_m)\left|\right|}^2,$表卷积操作，d(I)_m表示将图 像I按m倍进行降采样。β为公式2)中的平衡参数，根据视觉效果选取β＝0.2。

上述实现即为式(2)中的分项。

RMS与N_s的大小控制着迭代投影的结果，RMS越小，N_s越大则计算结果 的块状效应越不明显，但时间复杂度大大上升。对目前主要的10位量化遥感 数据，本发明实例残留误差RMS推荐阈值取为0.0001，其他量化位数的遥感 数据可按位数差进行转换。根据实际计算结果，迭代次数N_s一般设置为100。 即如果RMS＜0.0001或者i＞100，则停止迭代，否则，继续步骤8)。

最终得到的超分结果图像子块记为SubI_H。

9)收集各个子块对应的所有SubI_H，将结果写出到文件，为低分辨率图像 I_L的超分辨率结果I_H。

对发明实施例的分析发现，超分辨后的图像层次较为清晰，信噪比等各项 质量指标均有提高或未下降。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。