1 引言

图像复原是图像研究领域中的一个基本问题，在生物医学成像、天文、遥感等多种领域均有应用。图像复原通常主要包括图像去模糊和图像去噪，然而在遥感领域，线阵推扫式遥感成像方式下，遥感器在成像过程中的轻微振动即会导致每次线阵成像的采样间距各不相同，从而使得生成的遥感图像各像素点的采样间距呈不规则性，我们将此种现象称之为不规则采样。低分辨率遥感成像时期，由于遥感图像的质量相对较差，不规则采样导致的图像退化效果并不明显，所以没有引起科研人员的足够重视。随着遥感领域技术的发展，遥感器地面分辨率的逐渐提高，不规则采样所造成的遥感图像退化成为了一个不可忽视的因素，尤其在遥感测绘领域，微小的采样偏差即会带来较大的配准误差，从而导致巨大的结果误差,影响遥感图像质量与测绘精度。因此需要研究行之有效的方法消除此种退化因素，将不规则采样图像恢复为采样间距相同的规则采样图像。不规则采样复原方法在SAR成像^[1]复原以及MR成像^[2]复原等领域也有着重要作用。

针对图像的不规则采样问题，学者们已经提出了许多解决方法。传统方法主要基于插值和迭代算法^{[3, 4, 5, 6]}，此种算法主要由一维信号的不规则采样复原方法衍生而来，但由于需要迭代，算法速度较慢。文献[7, 8]提出了一种快速有效去除图像不规则采样的方法，该方法使用了自适应权重(adaptive weights)、共轭梯度法(conjugate gradient)以及 Toeplitz分块矩阵性质，故称之为ACT算法。此方法无需迭代，且利用了Toeplitz分块矩阵的特点进行快速计算。文献[9, 10]随后将ACT算法应用于带有多种退化因素的图像复原框架中，通过最小化总变差泛函以达到图像同时消除不规则采样、去模糊、去噪声等的目的。之后也有一些针对特殊噪声情况的改进方法^[11]。虽然该方法的复原效果良好，但依然无法摆脱总变差方法的固有缺陷，具体主要为复原图像带有比较严重的阶梯效应且原图像的纹理细节信息的保持较差。

针对遥感图像的去噪与复原问题，国内外学者已经提出了众多方法^{[12, 13, 14]}，能够有效地提高遥感图像质量。在众多的图像复原方法中，总变差(total variation,TV)方法最早由Rudin,Osher和Fatemi提出，由于其解空间为有界变差函数空间，而该空间中的函数允许存在边缘等不连续的信息，因此在图像复原中拥有较好的边缘保持能力。然而总变差方法会产生“阶梯”效应^[15]，且不能较好地保持图像中的纹理细节信息。非局部均值(non local means,NLM)去噪方法最早由Buades等人提出^[16]，它由Yaroslavsky邻接滤波器演化而来，对周围具有相似结构的像素灰度值进行平均以达到去噪的目的。试验表明，NLM方法在有效去噪的同时能够很好地保持纹理和结构，这个性质恰是总变差方法需要解决的问题。

本文的目的是将非局部均值(NLM)方法融入总变差图像复原框架下，使得复原后的遥感图像具有更好的纹理和结构保持能力，从而达到更佳的遥感图像消除不规则采样、去模糊、去噪复原结果。本文第2部分回顾相关研究成果并提出基于NLTV(non local total variation)的遥感图像去不规则采样、去卷积、去噪复原方法；第3部分给出所提出方法的试验结果，并对结果进行分析；第4部分对全文进行总结。

2 基于NLTV的遥感图像消除不规则采样图像复原方法 2.1 成像模型与退化因素分析

通常的图像获取系统的成像模型可以表示为

式中，Λ={λ_k}^N2_k=1⊆R²为采样点集合，它可以是规则采样也可以是不规则采样；u:R²→R为理想的无退化图像；h:R²→R为模糊核，其傅里叶频谱的能量主要集中在[－1/2,1/2]×[－1/2,1/2]区间内；n通常为均值为零，方差为σ的高斯白噪声；g:Λ→R为得到的采样后退化图像。

在不规则采样情况下，令Ω_N表示区间[0,N]×[0,N]，则不规则采样点集合可视为在规则采样点集合基础上加入随机偏移，可表示为

式中，ε:R²→R²为一个平滑且微小的震动函数。在遥感成像领域，卫星的微小震动以及遥感器的位置偏移都将导致遥感图像的不规则采样，通常将其建模为

式中，q≥1；a_k(x)为光滑调制函数；ξ_k为振动频率，其值在采样率的奈奎斯特频率以下。依据此模型得到的采样点位置扰动|ε(x)|不超过几个像素，并且扰动变化率|ε′(x)|不超过1/10个像素。在试验中，将式(3)进行简化，将震动函数ε=(ε₁,ε₂)建模为离散彩色噪声，且依然能够反映出震动函数的特点，对于ξ∈(Z²∩Ω_N) 的取值使得ε_i(x)的标准差为A，即，因此震动函数由振幅A和最小震动周期T_ε决定，其具体取值将在试验部分给出。

2.2 消除不规则采样的ACT算法

ACT算法是现今最有效的消除不规则采样方法之一。它的基本思想是将图像{u(i,j)}^N-1_i,j=0表示为离散的三角多项式，故不规则采样的插值图像可表示为

用矩阵形式可表示为 式中，a为规则采样图像的频谱。如果采样点在某个区域过于密集或过于稀疏，则以上系统将有较差的平衡性。为了改进条件数，将式(5)中的等式乘以权重函数 并且，ACT算法并不直接求解式(6)，而是将问题转化为求解如下等式 由于N²×N²矩阵T=S^*WS(W=diag({ω_k}_k=1,…,N²))为Toeplitz分块矩阵，因此Ta可以利用傅里叶方法进行快速计算，其时间复杂度为N²log₂(N²)。式(8)的解可由范数最小化求得^[8]，求得a后对其实行逆傅里叶变换即可得到不规则采样恢复为规则采样后的图像信号。

2.3 TV+ACT消除不规则采样图像复原模型

消除不规则采样图像复原的目标是根据式(1)的成像模型，在已知模糊核h、不规则采样点位置集合Λ、高斯白噪声方差σ以及带有不规则采样、模糊和噪声的退化图像g的情况下，尽可能地恢复原图像u。由式(1)的成像模型及上述ACT算法，基于总变差的消除不规则采样图像复原^[9]的无约束离散模型可表示为

式中，为成像系统的MTF；F为傅里叶变换矩阵；S为不规则采样矩阵；W为ACT算法中的加权矩阵；λ为总变差正则化参数；J_d(u)为总变差的离散形式。此模型解的存在性可以简单地证明，但由于无法确定不规则采样矩阵S是否为单射，因此模型解的唯一性无法保证。对式(9)的欧拉-拉格朗日方程进行整理，并使用坎贝尔投影算法迭代求解即可求得复原结果图像^[9]。

TV+ACT的图像复原模型能够有效地去除退化图像中的不规则采样、模糊、噪声等退化因素，但是其拥有TV复原模型的固有缺点，复原图像中存在较为明显的阶梯效应，并且对原图像纹理细节信息恢复较差。

2.4 基于NLTV的消除不规则采样图像复原方法

由于上节提到的TV+ACT方法的固有缺陷，希望可以在TV+ACT方法的主框架下，寻找有效的方法使得复原图像能够拥有较好的纹理细节保持能力，且减小“阶梯效应”等人为退化因素。

非局部均值(NLM)方法能够在有效去除噪声的同时保持图像的纹理细节信息，因此近年来引起众多学者的关注。文献[17]在图论和NLM方法的启发下提出了NL算子，NL算子不但拥有NLM方法的特点，并且能够融入正则化图像处理框架中，以结合两者各自的优势得到更好的处理结果。

首先，给出基于NLTV的消除不规则采样的图像复原模型

式中，基于NL算子的总变差正则化方程为 式中，NL算子的主要思想是将梯度和散度两个传统的局部定义，通过图论的相关思想扩展到非局部，使用NLM方法计算各像素间的相似度以得到各像素相互间的权重，从而构造NL梯度算子与NL散度算子。令，x,y∈Ω，NL算子定义为 权重系数w(x,y)的计算使用NLM中的形式 式中，v为已知图像；G_a为标准差为a的高斯核函数；h为滤波器系数，通常根据噪声的标准差确定；v(x+·)、v(y+·)分别表示图像v中以x点、y点为中心的图像块灰度值矩阵。由式(13)可知，两个像素的周围结构越相似，则权重系数越大。在试验实现中，基于存储空间与计算速度的考虑，通常采用半局部方法计算权重系数，即对于每个像素，只计算其周围一定窗口大小内像素的权值，认为其余像素对该像素的权值为0，这个窗口即为查找区域。因此，NLM方法中还需要使用相似性度量块大小与查找区域块大小两个参数，相似性度量块尺寸为对比相似度的各图像片的尺寸；查找区域块尺寸则确定了计算权值的区域，当查找区域块尺寸为原图像尺寸时，即为原始的全局部NLM方法。

对于带有不规则采样与模糊的退化图像，直接计算退化图像的NL算子往往无法得到较好的复原效果，因为带有不规则采样与模糊的退化图像已经改变了原图像的结构信息，无法反应原图像各像素间的相似程度，因此，在本文中使用算子分割的方法，将复原分为两个步骤：第1步对不规则采样与模糊进行复原，得到带噪复原图像；第2步使用带噪图像计算NL算子并进行去噪。具体方法如下，首先将求解式(10)改为求解下式

显然，当α足够大时，式(15)与式(10)等价。式(15)又可使用如下算子分裂方法计算 其中，式(16)是一个可微问题，可以使用共轭梯度法、Gauss-Seidel方法、傅里叶方法等求解。式(17)则是一个NLTV去噪问题，可以使用基于NL算子的扩展的坎贝尔投影算法求解。算子分割由于将一个复杂问题分解为两个较易解决的子问题，并且减少了计算量，因此可以同时提高算法的运行效率。在本文方法中，对式(16)、(17)进行迭代，每次迭代首先计算得到消除不规则采样、去模糊的带噪图像f_m，然后使用此去模糊带噪图像f_m进行NLTV方法去噪以得到迭代结果u_m。虽然f_m会放大原图像中的噪声，但NLM方法对于噪声具有较强的抗差性，因此能够获得更有效的权重系数，以得到更好的复原结果。对于初值u₀的选取，当仅为图像复原(无不规则采样)问题且H为列满秩矩阵时，式(15)为严格凸问题，能够证明算子分割方法的解序列{u_i}收敛到式(15)的不动点，也即收敛到式(15)的全局最小值点，因此无论u₀如何选取都不会影响最终的迭代结果^[18]。然而，在不规则采样问题中，由于不规则采样矩阵S的不确定性，无法得到上述结论。在试验中，当不规则采样偏移较小时，初值u₀的选取对迭代结果影响较小；当不规则采样偏移较大时，将使用Tikhonov正则化或TV正则化等其他方法预计算的复原图像作为初值能够得到更好的迭代结果，也能够减少收敛所需的迭代次数。

综上所述，基于NLTV的去不规则采样图像复原方法可归纳为：

步骤1 使用NUFFT计算T和。

步骤2 取u₀=0或更好的初始值，m=0，进行如下迭代

当 ，且m＜max iteru时

(1) 使用共轭梯度法求解f_m=F^－1(T′+αI)^－1F(F^－1+αu_m)；

(2) 对f_m计算NLM权重系数w(x,y)；

(3) 使用坎贝尔投影算法求解u_m+1，具体步骤为

① 取ξ⁰=0，p=0,进行如下循环

当，且p＜maxiterp时，计算

② 使用共轭梯度法求解u_m+1=f_m+α^－1λdiv_NL,其中，为①的迭代结果；

③ p=p+1；

(4)m=m+1。

步骤3 迭代结果即为所求的复原图像。

可以证明，u_m,m→∞是式(15)的解，具体的证明方法类似于已有方法^{[9, 19]}，在此不再赘述。但同样由于上文提到的不规则采样矩阵S的不确定性，无法证明解的唯一性，因此无法从理论上证明此解必为全局最优解。在实际试验中，当不规则采样偏移在2个像素以内时，通常能够得到良好的复原结果。

3 试验结果与分析

为了验证所提方法的性能并与其他方法进行比较，需要生成具有各种退化因素的试验图像。首先对原始图像加入模糊，此处使用SPOT-5的系统调制传递函数(MTF)模型

式中，α=0.58；β=0.14。使用式(14)对原始图像进行滤波操作，得到具有模糊的退化图像。然后根据式(4)，由A和T_ε的取值得到不规则采样点的位置，并使用NUFFT计算图像在不规则采样点处的像素灰度值，得到具有模糊和不规则采样的退化图像。本文试验中，为了测试算法的有效性，设置A=1，T_ε=10，因此退化图像将表现出明显的不规则采样效果。最后，在此基础上加入σ=2的高斯白噪声。

试验结果将本文方法与TV+ACT图像复原模型以及扩展的ACT算法进行比较，TV+ACT图像复原模型在2.3节已作介绍，扩展的ACT算法则是将ACT算法扩展为能够去卷积、去噪，具体方法为计算 优化问题的最小范数解，其中，T与b的定义如2.2节。扩展的ACT算法原理等同于逆滤波方法，因此在去除不规则采样和模糊的同时会放大噪声，导致复原效果较差。

试验参数设置方面，TV方法与本文的NLTV方法的通用参数使用相同的设置以便于比较。具体为λ=0.005，α=200，max iterp=10，max iteru=20。NLTV中权重参数的方面，相似性度量块尺寸的设置需要对噪声具有抗差性，通常设置为5像素×5像素或7像素×7像素，在图像噪声水平不高的情况下也可以设置更小的块尺寸，在试验中笔者发现相似性度量块尺寸对于本文试验结果的影响较小，为减少算法运算时间，将其设置为3像素×3像素；查找区域块尺寸的设置与图像的噪声水平和图像内容有关，较小的查找区域无法显著消除噪声，较大的查找区域不但增加运算负担，也会影响图像中的细节信息，本文试验中将其设置为5像素×5像素；滤波器系数为h=8。试验选用两幅细节信息差异较大的遥感图像：第1幅(图 1中的原始图像，123像素×123像素)的细节信息较少，主要为阶梯的纹理信息； 第2幅(图 3中的原始图像，237像素×237像素)中既有较为细小的纹理信息，也有丰富的细节信息。评价指标使用PSNR与SSIM两种，其中SSIM(structural similarity)为结构相似度，数值越大表示复原图像相比与原图像结构特征保持较好。

第1组试验中，退化图像中仅有不规则采样与噪声两种退化因素，3种方法的试验结果见图 1。

从试验结果中可以看到，3种方法均能达到较好的消除不规则采样的作用，然而ACT方法复原图像中依然能够明显观察到存留的噪声，TV+ACT方法与本文方法的复原结果则更为理想。

第2组试验的退化图像中带有不规则采样、模糊、噪声3种退化因素，使用两幅图像进行试验，试验结果见图 2~图 4。

由试验结果可以看到，相较于TV+ACT方法，本文方法具有更好的纹理细节信息保持能力。图 2试验中，本文方法对于阶梯的纹理恢复效果已接近原图像；图 3细节信息较为丰富，复原结果难以完全保持所有信息，但NLTV方法对于阶梯纹理以及细节信息也达到了较好的保持作用；图 4的细节对比图也可明显看到，NLTV方法相对于TV方法在纹理细节信息的保持方面更具优势。

第3组试验，在第2组试验的退化因素基础上加入2倍下采样的混叠影响，对TV+ACT方法与本文方法进行比较，试验结果见图 5。

由上述复原结果可以看到，本文方法略好于TV+ACT方法的复原结果，然而两者在复原方法上均未考虑混叠退化因素的复原，所以在复原结果中，左侧阶梯依然带有较大的混叠效应，导致复原图像视觉效果较差。因此，今后的工作也将针对带有不规则采样、模糊、混叠、噪声等多种退化因素的退化图像进行研究，寻找有效同时去除上述退化因素的有效手段。

在评价指标方面，从表 1可以看到在单纯去不规则采样与去噪的试验中，TV+ACT方法与本文方法的评价指标相差不大，均显著高于ACT方法。但在加入模糊的试验中，本文方法的评价指标均优于TV+ACT方法，且在带有混叠的试验中本文方法也稍稍领先。

4 结论

本文提出了一种遥感图像消除不规则采样的图像复原方法，该方法结合了不规则采样复原的ACT算法和总变差图像复原模型以达到同时去除多种图像退化因素的目的。并且，针对总变差图像复原模型具有阶梯效应以及纹理细节信息保持较差等不足，引入非局部均值(NLM)的思想，使用NL算子结合总变差方法(NLTV)以提高复原图像的纹理细节信息保持能力，提出了基于NLTV的遥感图像消除不规则采样图像复原模型。在模型的求解方面，使用算子分裂将求解过程分解为去卷积(消除不规则采样和去模糊)与去噪两个步骤，使得计算出的NLM权重更为准确地反应原始图像各像素间的关系，且提高了模型求解速度。试验结果表明，此方法能够有效地去除不规则采样、模糊、噪声等退化因素，相较于ACT算法与总变差图像复原模型，具有更好的纹理细节保持能力。该方法可应用于遥感图像复原、MR图像重建、SAR图像重建、压缩感知、图像修补等多种领域。今后的工作将主要针对带有混叠的多种退化因素问题进行研究，寻求更好的解决方法。