基于变分正则化的大气扰动图像恢复方法

摘要

基于变分正则化的大气扰动图像恢复方法，涉及图像处理。获取对固定目标的受到大气干扰的视频帧或者用模拟软件对一张图片附加模拟扰动以生成一组视频帧；对视频帧进行低秩分解，得到初始参考图像；使用基于非局部全变差正则项与可控核回归正则项的优化模型对参考图像进行优化，并用分离的Bregman算法实现加快该优化过程；用优化后的参考图像对视频帧进行B样条插值配准，得到经过配准后的视频帧；使用空间加权核范数最小化，将配准后的视频帧融合成一张近衍射极限图片；对近衍射极限图片进行去卷积处理，得到最后的去除了模糊和噪声的图片。提升扰动去除结果，得到了视觉清晰以及细节丰富的恢复图像，可用于高空对地观测、远距离监控等。

说明书

技术领域

本发明涉及图像处理，尤其是涉及可用于高空对地观测、远距离监控的一种基于变分正 则化的大气扰动图像恢复方法。

背景技术

随着社会的发展和科技的进步，越来越多的对地监控设备得到使用。例如无人机、高空 对地观测和远距离监控等技术得到广泛关注。与普通的拍摄行为不同的是，因为在远距离拍 摄中，拍摄目标与拍摄设备之间的距离较远，拍摄效果极容易受到大气扰动的影响，严重影 响拍摄效果。而这种扰动无法通过对硬件设备的提升来去除，为此如何从软件技术层面，运 用一些算法将由于大气扰动引起的视频模糊去除就十分重要了。

大气扰动造成的图像畸变主要原因是空气扰动造成时空变化的模糊及噪声的影响。因此 如何去掉时空变化的模糊及去噪，成为解决受大气扰动造成的图像畸变恢复的关键。现有的 解决图像畸变的流行工作可分为两大类，一类是基于参考图像配准的方法，先找到参考图像， 然后将视频中的帧与参考图像配准，这类方法并不对模糊核进行估计，而是使用固定的高斯 核来逼近模糊核，因此，该方法很难得到较高的图像恢复质量。另外一类方法是基于幸运图 的图像恢复方法，该类方法认为，小尺寸的图像块可以看作是3D场景在某一等平面角度形 成的，所以，小尺寸的图像块可以看作是等平面的。因此，这类方法首先把图像划分成图像 块序列，先在视频流中找到比较锐化的图像块，以此作为衍射极限图像块。然而，实验发现， 这类方法虽然压缩了大气扰动的影响，然而却增加了噪声的影响。

Zhu等人(Zhu X,Milanfar P.Removing atmospheric turbulence via space-invariant  deconvolution[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2013,35(1): 157-170.)提出了一种基于整体而非局部图像块的处理大气扰动造成的影响的方法，作者把视 频序列的时间平均作为参考图像，然后将视频中的每帧与参考图像进行配准、矫正，通过时 间轴的融合得到衍射极限图像。实验发现该方法在大气扰动造成的图像畸变矫正方面取得了 颇有希望的结果。然而，该方法将参考图像取为视频序列的时间平均，造成参考图像的细节 损失，图像序列的配准精度受到限制，从而使得矫正精度不高。另外，衍射极限图像的生成 也受限于融合时采用的时间轴平均的影响，该步骤本质上采用了时间轴去噪的方法。众所周 知，基于时间轴的平均这一大类去噪方法，去噪性能不高。从公式化的角度来看，Zhu的方 法没有从整体优化的角度来建模大气扰动造成的畸变，理论方面可以继续完善。

发明内容

本发明的目的在于提供可去除视频图像中因为大气扰动而产生的图像畸变的基于变分正 则化的大气扰动图像恢复方法。

本发明包括以下步骤：

(1)有一组在空中拍摄固定目标的视频帧序列{f_i}，它们受到了大气扰动的影响而产生 了模糊和畸变；

(2)对视频帧序列{f_i}进行低秩分解，得到初始参考图像I₀；

(3)初始化外层迭代，令外层循环变量Out_loop＝1，将I₀作为当前要优化的参考图像u， 设置最大外层循环次数max_out；

(4)初始化Bregman距离(Bregman distance)迭代，令中层循环变量Middle_loop＝1，将 初始帧序列{f_i}作为当前帧设置最大中层循环次数max_middle；

(5)对于当前帧序列中的每一帧利用B-spline配准方法计算参考图像到的形 变矩阵Φ_i；

(6)初始化向前-向后算子分离(forward-backward operator splitting)迭代，令Inter_loop＝1， 设置最大内层循环次数max_inter；

(7)利用{Φ_i}、和u进行向前步骤，得到迭代过程变量v；

(8)用非局部全变差(Non-local total variation)和可控核回归全变差(Steering kernel  regressiontotal variation)为正则项得到最优化式子，然后使用分离的Bregman迭代算法求解最 优化式子得到新的u；

(9)将Inter_loop>max_inter作为循环终止条件，若Inter_loop不满足大于max_inter的 条件，则在将Inter_loop的值增加1后返回步骤(7)，否则直接执行步骤(10)；

(10)利用{Φ_i}、u和初始帧序列{f_i}进行当前帧序列的更新；

(11)将Middle_loop>max_middle作为循环终止条件，若Middle_loop不满足大于 max_middle的条件，则在将Middle_loop值增加1后返回步骤(5)，否则直接执行步骤(12)；

(12)将优化后的参考图像作为当前参考图像u，即

(13)利用当前参考图像u，对当前视频帧进行B-spline近似配准，去掉几何形变， 得到配准后的视频帧序列{H_i}；

(14)对配准后的视频帧序列{H_i}进行低秩分解，得到结果作为新的参考图像u；

(15)将Out_loop>max_out作为循环终止条件，若Out_loop不满足大于max_out的条件， 则在将Out_loop的值增加1后返回步骤(4)，否则直接执行步骤(16)；

(16)使用空间加权核范数最小化(Space weighted nuclear norm minimization)将配准以后 的视频帧序列{H_i}融合成一张近衍射极限图片I_fusion；

(17)将融合图片I_fusion进行去卷积处理，得到最后的去除大气扰动后的图片I_final。

在步骤(7)中，所述向前步骤的公式如下：

$v=u-{\mathrm{\delta \Sigma }}_i{\Phi }_i^T({\Phi }_{i}u-{\tilde{f}}_i)$

式中，δ是参数，一般情况下设置其值为1。

在步骤(8)中，所述最优化式子如下：

$u=\underset{u}{\arg \text{><mi>m</mi><mi>i</mi><mi>n</mi>}}\{\frac{1}{2}{||u-v||}_2^2+{\lambda }_N\left|d_N\right|+{\lambda }_S\left|d_S\right|\}$

式中，λ_N和λ_S为参数，|d_N|表示非局部全变差正则化项，|d_S|表示可控核回归全变差正则 化项，表达公式如下：

|d_N＝(I-Q^N)u|

|d_S＝(I-Q^S)u|

其中，

$Q^N(i,j)=\left(\begin{array}{cc}{c}_{ij}^N& j\in S_i\\ 0& j\notin S_i\end{array}\right)$

$Q^S(i,j)=\left(\begin{array}{cc}{c}_{ij}^S& j\in L_i\\ 0& j\notin L_i\end{array}\right)$

$c_i^S=e_1^T{\left({\Psi }^T{W}_S\Psi \right)}^{-1}{\Psi }^T{W}_S$

$\Psi =\left(\begin{array}{ccc}1& {(u_i-u_j)}^T& {\mathrm{vech}}^T(u_i-u_j){(u_i-u_j)}^T\\ .& .& .\\ .& .& .\\ .& .& .\end{array}\right)$

式中，是由非局部权重w_N(i,j)串联得到的，S_i表示与u_i相似的图像块的位置集合，L_i表示与u_i相邻的图像块的位置集合，e₁＝[1 0 … 0]^T，W_S表示可控核回归权重，W_N表示非 局部权重，其表达式如下：

$w_N(i,j)=\exp (-\frac{{||R_{u_i}u-R_{u_j}u||}^2}{h_N^2})$

$w_S(i,j)=\frac{\sqrt{\det \left(C_i\right)}}{2{\mathrm{\pi h}}_S^2}\exp (-\frac{{(u_i-u_j)}^T{C}_i(u_i-u_j)}{2h_S^2})$

式中，表示从u中取出以u_i为中心的图像块，C_i是u_i处水平以及竖直方向梯度的下 方差矩阵，h_N，h_S是参数。

在步骤(8)中，所述使用分离的Bregman迭代算法求解最优化式子的迭代过程如下：

$u^{k+1}={[\underset{i}{\Sigma }{\Phi }_i^T{\Phi }_i+2{\lambda }_N\mu {(I-Q^N)}^T(I-Q^N)+2{\lambda }_S\mu {(I-Q^S)}^T(I-Q^S)]}^{-1}$

$*[{\mathrm{\Sigma \Phi }}_i^{T}v+2{\lambda }_N\mu {(I-Q^N)}^T(d_N-b_N^k)+2{\lambda }_S\mu {(I-Q^S)}^T(d_S-b_S^k)]$

$d_N^{k+1}=shrinkage\left(\right(I-Q^N)u^{k+1}+b_N^k,\frac{1}{\mu })$

$d_S^{k+1}=shrinkage\left(\right(I-Q^S)u^{k+1}+b_S^k,\frac{1}{\mu })$

$b_N^{k+1}=b_N^k+(I-Q^N)u^{k+1}-d_N^{k+1}$

$b_S^{k+1}=b_S^k+(I-Q^S)u^{k+1}-d_S^{k+1}$

其中，

$shrinkage(\chi ,\gamma )=\frac{\chi }{\left|\chi \right|}*max\left(\right|\chi |-\gamma ,0).$

在步骤(10)中，所述更新的公式如下：

${\tilde{f}}_i^{k+1}={\tilde{f}}_i^k+f_i-{\Phi }_i{u}^k$

在步骤(16)中，所述使用空间加权核范数最小化的具体步骤如下：

(16.1)将每张帧图片分割成许多大小相等的图像块，用H_ij表示第i个帧中以第j个像素 为中心的图像块，初始化像素位置j；

(16.2)认为对于每个图像块H_ij，在其他所有帧H_k上，同样中心在第j个像素的图像块 H_kj都是H_ij的相似块；将图像块H_ij的所有行一次连接成一个向量Z_ij，然后再将Z_ij与它的所 有相似块组合成一个矩阵Z_j＝[H_1j,H_2j,...,H_nj]；

(16.3)对Z_j进行奇异值分解得到F_j＝UΣV^T，再使用加权核范数最优化进行去噪，得到 其中S_W(Σ)＝max(Σ_ii-w_i,0)，其中σ_i(Z_j)为Z_j的第i个奇 异值，另外c＝2.8，ε＝10^-16；

(16.4)对的每一列求平均值，得到一个以位置j为中心的图像块

(16.5)若对于所有中心位置j都得到了则进行步骤(16.6)，否则将j移到下一个中 心位置，返回步骤(16.2)；

(16.6)将图像块序列拼凑成一副完整的图像I_fusion。

本发明从对参考图像进行优化入手，增加参考图像的细节，以使匹配后的图像序列更好， 从而提升扰动去除结果，实现了一种基于变分正则化的大气扰动图像恢复方法，得到了视觉 清晰以及细节丰富的恢复图像。

本发明具有以下突出优点：

1.本发明利用可控核回归全变差以及非局部全变差正则化项，对参考图像进行迭代优化， 其中可控核回归全变差正则项保留了参考图像的局部结构，非局部全变差正则项保留了非局 部的重复结构以及纹理。此方法可以很好地提高参考图像，本发明还使用了分离的Bregman 算法，加快了计算速度。

2.本发明使用低秩逼近的方法，使用空间加权原子核范数最小化模型对视频帧序列进行 融合优化，得到了一个去噪后的近衍射极限图像，最后通过去卷积得到了最后的清晰图片。 通过实验证实，本发明所使用的方法对真实视频帧序列以及模拟视频帧序列都有很好地效果， 并且优于当前其他最好的方法。

附图说明

图1是本发明的总体流程图；

图2是使用本发明和现有三种方法对烟囱的真实视频帧进行大气扰动去除后得到的图像；

图3是使用本发明和现有三中方法对城市的模拟扰动视频进行大气扰动取出后得到的图 像。

具体实施方式

参见图1，本发明实施例包括参考图像的生成与优化和视频帧的配准与融合两个部分。

一、参考图像的生成与优化：

步骤1，获取受到扰动的视频帧序列。

从网上下载远距离拍摄的对固定目标的视频，或者对一张清晰图片使用模拟软件，生成 加上扰动效果后的畸变帧序列。这样我们就得到了初始帧序列{f_i}。

步骤2，得到初始参考图像。

(2a)对初始帧序列{f_i}，每个帧f_i都用一个向量g_i表示，g_i的长度为f_i的像素数量。 这样我们可以用矩阵G＝[g₁,g₂,...,g_n]表示{f_i}，其中n为视频帧的数量。对G进行低秩分解 得到G的低秩部分L，L也是一个大小为n*m的矩阵，对L的每一列求平均，在恢复成一个 帧的矩阵形式，得到了初始图像I₀。

(2b)将I₀作为当前要优化的参考图像u。

步骤3，优化参考图像。

(3a)初始化Bregman距离(Bregman distance)迭代，令中层循环变量Middle_loop＝1， 将初始帧序列{f_i}作为当前帧

(3b)对于当前帧序列中的每一帧利用B-spline配准方法计算参考图像到的形 变矩阵Φ_i。

(3c)初始化向前-向后算子分离(forward-backward operator splitting)迭代，令Inter_loop＝ 1。

(3d)利用{Φ_i}、和u进行向前步骤得到v：

$v=u-{\mathrm{\delta \Sigma }}_i{\Phi }_i^T({\Phi }_{i}u-{\tilde{f}}_i)$

(3e)计算u的非局部权重以及可控核权重如下：

$w_N(i,j)=\exp (-\frac{{||R_{u_i}u-R_{u_j}u||}^2}{h_N^2})$

$w_S(i,j)=\frac{\sqrt{\det \left(C_i\right)}}{2{\mathrm{\pi h}}_S^2}\exp (-\frac{{(u_i-u_j)}^T{C}_i(u_i-u_j)}{2h_S^2})$

式中，表示从u中取出以u_i为中心的图像块，C_i是u_i处水平以及竖直方向梯度的下 方差矩阵，h_N＝20，h_S＝2。

(3f)计算u的非局部全变差正则化项和可控核回归全变差正则化项如下：

|d_N＝(I-Q^N)u|

|d_S＝(I-Q^S)u|

其中，

$Q^N(i,j)=\left(\begin{array}{cc}{c}_{ij}^N& j\in S_i\\ 0& j\notin S_i\end{array}\right)$

$Q^S(i,j)=\left(\begin{array}{cc}{c}_{ij}^S& j\in L_i\\ 0& j\notin L_i\end{array}\right)$

$c_i^S=e_1^T{\left({\Psi }^T{W}_S\Psi \right)}^{-1}{\Psi }^T{W}_S$

$\Psi =\left(\begin{array}{ccc}1& {(u_i-u_j)}^T& {\mathrm{vech}}^T(u_i-u_j){(u_i-u_j)}^T\\ .& .& .\\ .& .& .\\ .& .& .\end{array}\right)$

式中，λ_N和λ_S为参数，是由非局部权重w_N(i,j)串联得到的，S_i表示与u_i相似的图像 块的位置集合，L_i表示与u_i相邻的图像块的位置集合，

(3g)根据以下迭代公式进行10次迭代计算：

$u^{k+1}={[\underset{i}{\Sigma }{\Phi }_i^T{\Phi }_i+2{\lambda }_N\mu {(I-Q^N)}^T(I-Q^N)+2{\lambda }_S\mu {(I-Q^S)}^T(I-Q^S)]}^{-1}$

$*[{\mathrm{\Sigma \Phi }}_i^{T}v+2{\lambda }_N\mu {(I-Q^N)}^T(d_N-b_N^k)+2{\lambda }_S\mu {(I-Q^S)}^T(d_S-b_S^k)]$

$d_N^{k+1}=shrinkage\left(\right(I-Q^N)u^{k+1}+b_N^k,\frac{1}{\mu })$

$d_S^{k+1}=shrinkage\left(\right(I-Q^S)u^{k+1}+b_S^k,\frac{1}{\mu })$

$b_N^{k+1}=b_N^k+(I-Q^N)u^{k+1}-d_N^{k+1}$

$b_S^{k+1}=b_S^k+(I-Q^S)u^{k+1}-d_S^{k+1}$

其中，

$shrinkage(\chi ,\gamma )=\frac{\chi }{\left|\chi \right|}*max\left(\right|\chi |-\gamma ,0)$

(3h)若Inter_loop<＝3，在将Inter_loop的值增加1后返回步骤(3d)，否则直接执行步骤 (3i)。

(3i)更新${\tilde{f}}_i^{k+1}={\tilde{f}}_i^k+f_i-{\Phi }_i{u}^k$

(3j)若Middle_loop<＝3，则将Middle_loop值增加1后返回步骤(3b)，否则直接执行步 骤(4a)。

步骤4，配准新的视频帧。

(4a)将优化后的参考图像作为当前参考图像：

(4b)利用当前参考图像u，对当前视频帧进行B-spline配准得到新的视频帧序列 {H_i}。

步骤5，融合新视频帧。

(5a)将每张帧图片分割成许多大小相等的图像块，用H_ij表示第i个帧中以第j个像素为 中心的图像块。下面初始化像素位置j。

(5b)我们认为对于每个图像块H_ij，在其他所有帧H_k上，同样中心在第j个像素的图像 块H_kj都是H_ij的相似块。我们将图像块H_ij的所有行一次连接成一个向量Z_ij，然后再将Z_ij与 它的所有相似块组合成一个矩阵Z_j＝[H_1j,H_2j,...,H_nj]。

(5c)对Z_j进行奇异值分解得到F_j＝UΣV^T，再使用加权核范数最优化进行去噪，得到 其中S_W(Σ)＝max(Σ_ii-w_i,0)，其中σ_i(Z_j)为Z_j的第i个奇 异值，另外c＝2.8，ε＝10^-16。

(5d)对的每一列求平均值，得到一个以位置j为中心的图像块

(5e)若对于所有中心位置j都得到了则进行步骤(5f)，否则将j移到下一个中心 位置，返回步骤(5b)。

(5f)将图像块序列拼凑成一副完整的图像I_fusion。

步骤6，将融合图片I_fusion进行去卷积处理，得到最后的去除大气扰动后的图片I_final。

本发明通过以下实验进行优点和有效性的证明。

1.实验条件：

实验室用台式机参数：CPU为Inter(R)Core(TM)i7-2600，主频为3.40GHz，内存为4G， 操作系统为Win764位系统，实验平台为Matlab2014b。

实验室用本发明提出的基于变分正则化的大气扰动图像恢复方法，其中视频帧配准的方 法见参考文献“Zhu X,Milanfar P.Removing atmospheric turbulence via space-invariant  deconvolution[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,2013,35(1): 157-170.”

2.实验结果以及结果分析：

实验一，用本发明和现有的方法进行真实视频帧的大气扰动图像恢复。

从网络上下载一组从远处拍摄的烟囱的视频，它受到热空气的影响而产生了模糊。我们 选取其中的100帧进行实验，每一帧的大小为256×256，采用本发明和现有的三种方法分别 进行扰动的去除。去除扰动的结果如图2所示，其中图2(a)是没受到扰动的真实图片；图 2(b)是其中一张受到扰动的帧；图2(c)是选取幸运图片的方法得到的结果；图2(d)是 由O.Oreifej等人提出的TwoStage的方法得到的结果；图2(e)是由X.zhu等人提出的NDL 的方法得到的结果；图2(f)是本发明得到的结果。

图2的实验结果说明：从直观视觉上看，图2(c)、图2(d)和图2(e)的结果都很模 糊，细节不清晰，丢失了真实图像的很多细节信息。图2(f)中则显得清晰了许多，细节丰 富，整体上更接近真实图像。

实验二，用本发明和现有的方法进行模拟视频帧的大气扰动图像恢复。

从网络上下载一张从城市上空俯拍的图片，利用扰动模拟软件生成60帧的扰动模拟视频， 每一帧的大小为256×256，采用本发明和现有的三种方法分别进行扰动的去除。去除扰动的 结果如图3所示，每张图的表示和实验一中一致。表1给出了本发明与现有三种发明的PSNR 和SSIM的对比结果。

图3和表1的实验结果说明：图3(b)中受到扰动的帧有着很大的扭曲和模糊；图3(c) 和图3(e)虽然矫正了扭曲，但是还是很模糊，边缘杂乱，细节信息丢失严重；图3(d)的 结果相比前面要清晰，但还是略显模糊，桥的两头一边大一边小，与真实图片差距略大；图 3(e)的结果相比前面的图都要好，图像清晰，扭曲较少，细节信息更丰富，与真实图片最 为接近。表1中，本发明的PSNR和SSIM值都高于现有的三种方法。

表1

  幸运图片 TwoStage NDL 本发明 PSNR 20.3783 21.9903 21.2538 22.9349 SSIM 0.45 0.6448 0.53 0.7486

本发明主要解决了由于大气扰动造成视频时空变化的模糊及噪声的问题。其主要的实现 步骤是：1)获取对固定目标的受到大气干扰的视频帧或者用模拟软件对一张图片附加模拟扰 动以生成一组视频帧。2)对视频帧进行低秩分解，得到初始参考图像。3)使用基于非局部 全变差正则项与可控核回归正则项的优化模型对参考图像进行优化，并用分离的Bregman算 法实现加快该优化过程。4)用优化后的参考图像对视频帧进行B样条插值配准，得到经过 配准后的视频帧。5)使用空间加权核范数最小化，将配准后的视频帧融合成一张近衍射极限 图片。6)对近衍射极限图片进行去卷积处理，得到最后的去除了模糊和噪声的图片。本发明 可用于高空对地观测、远距离监控等。