CN103164850A - 一种基于压缩感知的多聚焦图像融合方法及装置 - Google Patents

Abstract

本发明公开了一种基于压缩感知的多聚焦图像融合方法，属于图像信号处理技术领域。本发明首先对待融合的两幅图像分别进行压缩采样，得到两幅图像的观测向量；对两幅图像的观测向量进行融合，得到融合观测向量；利用融合观测向量重构出待融合的两幅图像的融合图像。本发明采用自适应局部能量测度(ALEM)融合准则对观测向量进行融合，并利用快速持续线性增广Lagrangian法(FCLALM)进行图像重构。本发明还公开了一种基于压缩感知的多聚焦图像融合装置。相比现有技术，本发明可获得更高的图像融合质量，且具有更高的信噪比和更快速的收敛速度。

Description

一种基于压缩感知的多聚焦图像融合方法及装置

技术领域

本发明公开了一种基于压缩感知的多聚焦图像融合方法及装置，属于图像信号处理技术领域。

背景技术

对某一场景成像时，由于光学系统的聚焦范围有限，一般的光学成像系统难以对同一场景中不同焦距上的物体都形成清晰的像。当成像系统的焦点聚焦在某个物体上时，它可以在像平面上形成一个清晰的像。这时，位于其它位置上的物体，在像平面上所成的像将呈现出不同程度的模糊。因此，仅依靠成像系统难以得到同一场景中所有物体都清晰的图像。为了更加全面、真实地反映某一场景的信息，希望得到一幅该场景中所有物体都清晰的图像。解决这一问题的一种方法就是对场景中的物体分别进行聚焦，得到多幅该场景的多聚焦图像，然后对这些多聚焦的图像进行融合。图像融合是将传感器在同一时间或不同时间获取的关于某个场景的两个或两个以上的图像或图像序列加以综合，以获得信息更加丰富、更适合于视觉感知或更适合于进一步处理的新图像的信息处理过程，其应用已遍及医学图像、遥感和计算机视觉等。多聚焦图像融合是图像融合的一个重要分支。

在过去的十年间，提出了很多图像融合的方法，主要分为基于单个像素的融合准则和基于区域特征的融合准则。基于单个像素的融合准则虽然比基于区域特征的融合准则简单，但前者有对边缘敏感等缺点。此外，所有的样本都是基于Nyquist采样定理进行采样的，这势必会增加数据存储、计算和处理的负担。因此，近年提出了基于压缩感知(Compressed Sensing,CS)的图像融合方法。CS可以减少计算和传输的损耗，成为图像融合的有效方法。

近年，提出了许多基于CS的图像融合方法。Wan等提出了用2维Fourier矩阵进行观测并采用最大值选取的方法融合观测向量。Li等提出了用2维Fourier矩阵进行观测并采用线性加权平均的方法融合观测向量。然而，Fourier矩阵仅与时域稀疏的信号不相干，应用范围较窄。Luo等提出用相似性分类的思想融合观测向量，但这种方法的计算量很大。Zebhi等提出了基于DCT的稀疏采样方法的图像融合，但当采样矩阵不是正交矩阵时，这种方法的计算量非常大。Luo等提出了一种有效的线性加权平均的融合准则，但该方法没有考虑多聚焦图像的自身特性。

基于CS的图像融合方法有两大关键技术：(1)设计合适的融合准则融合观测向量；(2)设计一种重构算法解决l₁-最小化问题。一种好的融合准则应保留源图像的所有显著特征，并且尽可能消除块效应。例如最大值选取(Maximum Selection,MS)准则、标准差加权平均(Standard Deviation Weighted Average,SDWV)准则、熵测度加权平均(EntropyMetrics Weighted Average,EMWV)准则、简单加权平均(Simple Weighted Average,SWV)准则等一些基于CS的多聚焦图像融合方法被提出，但这些融合准则均存在缺点。此外，增广Lagrangian法(Augmented Lagrangian Method,ALM)虽可以有效解决l1-最小化问题，但其收敛速度太慢。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于压缩感知的多聚焦图像融合方法及装置，采用自适应局部能量测度(ALEM)融合准则对观测向量进行融合，能够提供更高的融合质量。

本发明的基于压缩感知的多聚焦图像融合方法，首先对待融合的两幅图像分别进行压缩采样，得到两幅图像的观测向量；对两幅图像的观测向量进行融合，得到融合观测向量；利用融合观测向量重构出待融合的两幅图像的融合图像；所述对两幅图像的观测向量进行融合，具体按照以下方法：

首先将待融合的两幅图像的观测向量y₁=(u₁,u₂,…,u_M)、y₂=(v₁,v₂,…,v_M)分段表示为y₁=(u₁,u₂,…,u_P)^T和y₂=(v₁,v₂,…,v_P)^T；其中，u_j=(u_(j-1)·q+1,u_(j-1)·q+2,…,u_j·q)^T，v_j=(v_(j-1)·q+1,v_(j-1)·q+2,…,v_j·q)^T，j=1,2,…P,P×q=M；

然后按照下式分别计算分段向量u_j和v_j的能量匹配度D_j：

$D_j=\frac{2E_{12j}}{E_{1j}+E_{2j}},$

式中，E_1j、E_2j分别为u_j、v_j的局部能量测度，E_12j为u_j和v_j的混合能量测度，根据下式得到： $E_{12j}=\underset{(j-1)\·q+1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}\left|u_j{v}_j\right|/q;$

根据能量匹配度D_j与预设阈值T∈(0.5,1)的比较结果，按照以下方法计算融合观测向量y=(w₁,w₂,…,w_P)^T的分段向量w_j，j=1,2,…P：

若D_j<T，则：

$w_j=\left\{\begin{array}{cc}{u}_j& \mathrm{if}{E}_{1j}\&GreaterEqual;E_{2j}\\ v_j& \mathrm{if}{E}_{1j}<E_{2j}\end{array}\right.(j=\mathrm{1,2},...,P),$

若D_j≥T，则：

w_j=ωu_j+(1-ω)v_j(j=1,2,…,P)

其中ω为自适应加权因子，按照下式计算：

$\mathrm{\&omega;}=\frac{E_{1j}}{E_{1j}+E_{2j}};$

最终得到融合观测向量y=(w₁,w₂,…,w_P)^T。

作为本发明技术方案的进一步改进方案，在利用融合观测向量重构出待融合的两幅图像的融合图像时，所述融合图像的稀疏系数Θ通过求解以下的l₁-最小化问题得到：

$\begin{array}{cccc}\underset{\mathrm{\&Theta;}}{\min }& {\left|\right|\mathrm{\&Theta;}\left|\right|}_1& s.t.& y=\mathrm{A\&Theta;}\end{array},$

式中，y为融合观测向量，A为对所述压缩采样所使用的观测矩阵；

求解采用快速持续线性增广Lagrangian法，具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数μ₀>0，最大惩罚参数μ_max>μ₀，参数α>0，初始变量t₀=1，放大因子η>1，初始稀疏系数Θ₀，初始辅助变量

初始乘子d₀，初始迭代次数k=0；

步骤2、更新辅助变量 $\stackrel{\&OverBar;}{\mathrm{\&Theta;}}:{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k});$

步骤3、更新变量 $t:t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤4、更新稀疏系数 $\mathrm{\&Theta;}:{\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k);$

步骤5、更新惩罚参数μ：μ_k+1=min{ημ_k,μ_max}；

步骤6、更新乘子 $d:d_{k+1}=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k);$

步骤7、若满足终止条件则算法终止，否则，令k=k+1转步骤2。

优选地，所述压缩采样利用离散小波变换进行稀疏表示；在利用融合观测向量重构出待融合的两幅图像的融合图像时，通过对融合图像的稀疏系数进行反离散小波变换得到。

优选地，所述压缩采样采用随机高斯矩阵进行观测。

根据同一发明构思可得到本发明的基于压缩感知的多聚焦图像融合装置，包括压缩采样模块、融合模块以及图像重构模块；所述压缩采样模块用于对待融合图像进行压缩采样，得到待融合图像的观测向量；所述融合模块用于对两幅待融合图像的观测向量进行融合，得到融合观测向量；所述图像重构模块用于对融合观测向量进行图像重构，得到融合图像；所述融合模块对两幅待融合图像的观测向量进行融合，具体按照以下方法：

首先将待融合的两幅图像的观测向量y₁=(u₁,u₂,…,u_M)、y₂=(v₁,v₂,…,v_M)分段表示为y₁=(u₁,u₂,…,u_P)^T和y₂=(v₁,v₂,…,v_P)^T；其中，u_j=(u_(j-1)·q+1,u_(j-1)·q+2,…,u_j·q)^T，v_j=(v_(j-1)·q+1,v_(j-1)·q+2,…,v_j·q)^T，j=1,2,…P,P×q=M；

然后按照下式分别计算分段向量u_j和v_j的能量匹配度D_j：

$D_j=\frac{2E_{12j}}{E_{1j}+E_{2j}},$

式中，E_1j、E_2j分别为u_j、v_j的局部能量测度，E_12j为u_j和v_j的混合能量测度，根据下式得到： $E_{12j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}\left|u_j{v}_j\right|/q;$

根据能量匹配度D_j与预设阈值T∈(0.5,1)的比较结果，按照以下方法计算融合观测向量y=(w₁,w₂,…,w_P)^T的分段向量w_j，j=1,2,…P：

若D_j<T，则：

$w_j=\left\{\begin{array}{cc}{u}_j& \mathrm{if}{E}_{1j}\&GreaterEqual;E_{2j}\\ v_j& \mathrm{if}{E}_{1j}<E_{2j}\end{array}\right.(j=\mathrm{1,2},...,P),$

若D_j≥T，则：

w_j=ωu_j+(1-ω)v_j(j=1,2,…,P)

其中ω为自适应加权因子，按照下式计算：

$\mathrm{\&omega;}=\frac{E_{1j}}{E_{1j}+E_{2j}};$

最终得到融合观测向量y=(w₁,w₂,…,w_P)^T。

进一步地，所述图像重构模块通过快速持续线性增广Lagrangian法(以下简称FCLALM)求解以下的l₁-最小化问题得到融合图像的稀疏系数Θ：

$\begin{array}{cccc}\underset{\mathrm{\&Theta;}}{\min }& {\left|\right|\mathrm{\&Theta;}\left|\right|}_1& s.t.& y=\mathrm{A\&Theta;}\end{array},$

式中，y为融合观测向量，A为对所述压缩采样所使用的观测矩阵；

具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数μ₀>0，最大惩罚参数μ_max>μ₀，参数α>0，初始变量t₀=1，放大因子η>1，初始稀疏系数Θ₀，初始辅助变量

初始乘子d₀，初始迭代次数k=0；

步骤2、更新辅助变量 $\stackrel{\&OverBar;}{\mathrm{\&Theta;}}:{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k});$

步骤3、更新变量 $t:t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤4、更新稀疏系数 $\mathrm{\&Theta;}:{\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k);$

步骤5、更新惩罚参数μ：μ_k+1=min{ημ_k,μ_max}；

步骤6、更新乘子 $d:d_{k+1}=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k);$

步骤7、若满足终止条件则算法终止，否则，令k=k+1转步骤2。相比现有技术，本发明具有以下有益效果：

（1）采用自适应局部能量测度(ALEM)融合准则对观测向量进行融合，能够提供更高的融合质量。

（2）采用快速持续线性增广Lagrangian法来重构出融合图像的稀疏系数，具有更高的信噪比和更快速的收敛速度。

附图说明

图1为本发明的基于压缩感知的多聚焦图像融合方法流程框图；

图2a、图2b、图2c分别为参考图像Lab、Pepsi、Clock在基于DCT基和DWT基的FCLALM算法的PSNRs；

图3a、图3b分别为采用不同算法重构Pepsi参考图像的PSNR和迭代次数；

图4a～图4h依次为Lab参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像；

图5a～图5h依次为Pepsi参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像；

图6a～图6h依次为Clock参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的基于压缩感知的多聚焦图像融合装置，包括压缩采样模块、融合模块以及图像重构模块；所述压缩采样模块用于对待融合图像进行压缩采样，得到待融合图像的观测向量；所述融合模块用于对两幅待融合图像的观测向量进行融合，得到融合观测向量；所述图像重构模块用于对融合观测向量进行图像重构，得到融合图像。

正如背景技术部分所述，基于CS的图像融合方法有两大关键技术：(1)设计合适的融合准则融合观测向量；(2)设计一种重构算法解决l₁-最小化问题。一种好的融合准则应保留源图像的所有显著特征，并且尽可能消除块效应。例如最大值选取(MS)准则、标准差加权平均(SDWV)准则、熵测度加权平均(EMWV)准则、简单加权平均(SWV)准则等一些基于CS的多聚焦图像融合方法被提出，但这些融合准则均存在缺点。局部能量通常比平均梯度更体现融合图像的高分辨率等特性，具有大的局部能量的图像通常包含更多的信息。此外，基于区域特征的融合准则可以改善图像的清晰度。基于上述考虑，本发明为了获得更高的融合质量，提出了一种自适应局部能量测度(ALEM)融合准则对观测向量进行融合。

本发明的基于压缩感知的多聚焦图像融合方法，包括以下步骤：

步骤A、对待融合的两幅图像分别进行压缩采样，得到两幅图像的观测向量。

对待融合图像进行压缩采样，可采用现有各种压缩采样方法，例如首先对源图像在DCT域或DFT域进行稀疏表示，然后进行观测。本发明优选在小波域进行稀疏表示，并优选随机高斯矩阵进行观测。具体如下：在采样阶段，将待融合的源图像x₁∈R^n×n和x₂∈R^n×n分别列向量化，并用离散小波变换（DWT）进行稀疏表示，得到各自的长度为N=n²的稀疏系数Θ₁和Θ₂；然后，用随机高斯矩阵A分别对Θ₁和Θ₂进行观测，得到观测向量y₁和y₂。

步骤B、对两幅图像的观测向量进行融合，得到融合观测向量。

局部能量通常比平均梯度更体现融合图像的高分辨率等特性，具有大的局部能量的图像通常包含更多的信息。此外，基于区域特征的融合准则可以改善图像的清晰度。基于上述考虑，本发明提出了一种新的称为自适应局部能量测度（ALEM）的融合准则。具体如下：

首先将待融合的两幅图像的观测向量y₁=(u₁,u₂,…,u_M)、y₂=(v₁,v₂,…,v_M)分段表示为y₁=(u₁,u₂,…,u_P)^T和y₂=(v₁,v₂,…,v_P)^T；其中，u_j=(u_(j-1)·q+1,u_(j-1)·q+2,…,u_j·q)^T，v_j=(v_(j-1)·q+1,v_(j-1)·q+2,…,v_j·q)^T，j=1,2,…P,P×q=M，M为观测向量的维数。

分段向量u_j和v_j的传统局部能量测度分别定义为：

$E_{1j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}{u}_j^2/q---\left(1\right)$

$E_{2j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}{v}_j^2/q---\left(2\right)$

分段向量u_j和v_j的传统混合局部能量测度定义为：

$E_j=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}{u_{j}v}_j/q---\left(3\right)$

局部能量测度为现有技术，具体可参见文献[Local energy based multi-focus image fusion methodon curvelet transforms]。

为了避免u_jv_j的值为负数，本发明用|u_jv_j|代替u_jv_j，则得到分段向量u_j和v_j的混合局部能量测度：

$E_{12j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}\left|u_j{v}_j\right|/q---\left(4\right)$

定义分段向量u_j和v_j的能量匹配度D_j：

$D_j=\frac{2E_{12j}}{E_{1j}+E_{2j}}---\left(5\right)$

若D_j<T(其中T∈(0.5,1)为预设阈值)，则融合观测向量y=(w₁,w₂,…,w_P)^T可以由下式得到

$w_j=\left\{\begin{array}{cc}{u}_j& \mathrm{if}{E}_{1j}\&GreaterEqual;E_{2j}\\ v_j& \mathrm{if}{E}_{1j}<E_{2j}\end{array}\right.(j=\mathrm{1,2},...,P)---\left(6\right)$

若D_j≥T，定义自适应加权因子ω为

$\mathrm{\&omega;}=\frac{E_{1j}}{E_{1j}+E_{2j}}---\left(7\right)$

则融合观测向量y=(w₁,w₂,…,w_P)^T可以由下式得到

w_j=ωu_j+(1-ω)v_j(j=1,2,…,P)    (8)

当D_j<T时，ALEM采用选取最大绝对值的融合准则；当D_j≥T时，ALEM采用基于局部能量测度的自适应选取准则。

因为x=ΨΘ且Ψ为正交基，故有

${\left|\right|x\left|\right|}_2^2={\left|\right|\mathrm{\&Psi;\&Theta;}\left|\right|}_2^2={\left(\mathrm{\&Psi;\&Theta;}\right)}^T\mathrm{\&Psi;\&Theta;}={\mathrm{\&Theta;}}^T{\mathrm{\&Psi;}}^T\mathrm{\&Psi;\&Theta;}={\mathrm{\&Theta;}}^T\mathrm{\&Theta;}={\left|\right|\mathrm{\&Theta;}\left|\right|}_2^2--\left(9\right)$

也就是说稀疏系数Θ含有x的全部能量。此外，观测矩阵A的约束等距条件(RestrictedIsometry Property,RIP)为

$(1-{\mathrm{\&delta;}}_K){\left|\right|\mathrm{\&Theta;}\left|\right|}_2^2\&le;{\left|\right|\mathrm{A\&Theta;}\left|\right|}_2^2\&le;(1+{\mathrm{\&delta;}}_K){\left|\right|\mathrm{\&Theta;}\left|\right|}_2^2---\left(10\right)$

即有

$(1-{\mathrm{\&delta;}}_K){\left|\right|\mathrm{\&Theta;}\left|\right|}_2^2\&le;{\left|\right|y\left|\right|}_2^2\&le;(1+{\mathrm{\&delta;}}_K){\left|\right|\mathrm{\&Theta;}\left|\right|}_2^2---\left(11\right)$

其中δ_K∈(0,1)为常数。由RIP条件可以得到观测向量y含有Θ的几乎全部能量。从而，观测向量y含有x的几乎全部能量。因具有大的局部能量的图像通常包含更多的信息。所以，相比其他基于CS的融合准则，ALEM是一种更好的融合准则。

步骤C、利用融合观测向量重构出待融合的两幅图像的融合图像。

利用得到的融合观测向量y=(w₁,w₂,…,w_P)^T，即可利用的现有的各种重构算法重构出融合图像的稀疏系数Θ，进而重构出融合图像。图像重构的应用已经遍及到医学图像、天文图像、核磁共振成像、图像和视频编码等领域。CS技术是近几年出现的一种新兴的信号采样和压缩技术，基于该理论所获得的原始信号采样值，不仅数量大大低于基于传统Nyquist准则的采样值，而且CS技术还具有对未知信号边采样边压缩的特性。基于CS的图像重构问题是一个病态的欠定问题，可以利用一些先验信息对其进行求解。l₁-最小化问题便是一种标准的小波域图像重构问题。基于CS的图像重构问题假设图像x∈R^N可以用一个正交基

表示为

其中，为投影系数。当Θ满足||Θ||₀=K(K<<N)时，称图像x是K-稀疏的，其中||Θ||₀表示向量Θ中非零元素的个数。

若图像x在Ψ域上是K-稀疏的，我们直接用矩阵A∈R^M×N对其稀疏系数Θ进行线性测量，得到观测向量y∈R^M为

y=AΘ    (13)

其中A称为观测矩阵。

稀疏系数Θ可以通过求解如下的l₁-最小化问题

$\begin{array}{cccc}\underset{\mathrm{\&Theta;}}{\min }& {\left|\right|\mathrm{\&Theta;}\left|\right|}_1& s.t.& y=\mathrm{A\&Theta;}\end{array}---\left(14\right)$

进而通过式(12)重构出图像x。

l₁-最小化问题可以转化为线性问题通过经典的线性规划(Linear Programming,LP)的方法对其进行求解。但是，当观测矩阵为大规模、稠密的随机高斯矩阵时，这种方法的收敛速度将非常慢。因此，本发明提出了采用快速持续线性增广Lagrangian法（FCLALM）来解决l₁-最小化问题，该算法重构的图像与基于LP算法重构的图像具有同样好的精度，但该算法具有更快速的收敛速度。

式(14)的l₁-最小化问题的乘子函数为

$\mathrm{\&Gamma;}(\mathrm{\&Theta;},\mathrm{\&lambda;},\mathrm{\&mu;})={\left|\right|\mathrm{\&Theta;}\left|\right|}_1-{\mathrm{\&lambda;}}^T(\mathrm{A\&Theta;}-y)+\frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{A\&Theta;}-y\left|\right|}_2^2={\left|\right|\mathrm{\&Theta;}\left|\right|}_1+\frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{A\&Theta;}-y-\frac{\mathrm{\&lambda;}}{\mathrm{\&mu;}}\left|\right|}_2^2---\left(5\right)$

其中，λ为Lagrange乘子，μ为惩罚参数。

对于给定的λ_k有

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\underset{\mathrm{\&Theta;}}{\min }\mathrm{\&Gamma;}(\mathrm{\&Theta;},{\mathrm{\&lambda;}}_k,\mathrm{\&mu;})=\underset{\mathrm{\&Theta;}}{\min }{\left|\right|\mathrm{\&Theta;}\left|\right|}_1+\frac{\mathrm{\&mu;}}{2}{\left|\right|\mathrm{A\&Theta;}-(y+\frac{{\mathrm{\&lambda;}}_k}{2})\left|\right|}_2^2---\left(16\right)$

优化问题(16)很难求解，借助于线性化的思想，提出用线性化增广Lagrangian法(LALM)来求解此优化问题。通过线性化的思想，二次项

可以近似表示为

$\frac{1}{2}{\left|\right|\mathrm{A\&Theta;}-(y+\frac{{\mathrm{\&lambda;}}_k}{\mathrm{\&mu;}})\left|\right|}_2^2\&ap;\frac{1}{2}{\left|\right|A{\mathrm{\&Theta;}}_k-(y+\frac{{\mathrm{\&lambda;}}_k}{\mathrm{\&mu;}})\left|\right|}_2^2+<g_k,\mathrm{\&Theta;}-{\mathrm{\&Theta;}}_k>+\frac{1}{2\mathrm{\&alpha;}}{\left|\right|\mathrm{\&Theta;}-{\mathrm{\&Theta;}}_k\left|\right|}_2^2---\left(17\right)$

其中α>0为近似参数，<g_k,Θ-Θ_k>表示向量g_k和向量Θ-Θ_k的内积，且

$g_k=A^T(A{\mathrm{\&Theta;}}_k-(y+\frac{{\mathrm{\&lambda;}}_k}{\mathrm{\&mu;}}))---\left(18\right)$

为

在Θ_k的梯度。将式(18)代入式(16)得到

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\underset{\mathrm{\&Theta;}}{\min }{\left|\right|\mathrm{\&Theta;}\left|\right|}_1+\mathrm{\&mu;}(<g_k,\mathrm{\&Theta;}-{\mathrm{\&Theta;}}_k>+\frac{1}{2\mathrm{\&alpha;}}{\left|\right|\mathrm{\&Theta;}-{\mathrm{\&Theta;}}_k\left|\right|}_2^2)$

$=\underset{\mathrm{\&Theta;}}{\min }\frac{\mathrm{\&alpha;}}{\mathrm{\&mu;}}{\left|\right|\mathrm{\&Theta;}\left|\right|}_1+\frac{1}{2}{\left|\right|\mathrm{\&Theta;}-({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{g}_k)\left|\right|}_2^2$

$=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{g}_k,\frac{\mathrm{\&alpha;}}{\mathrm{\&mu;}})---\left(19\right)$

$=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-(y+\frac{{\mathrm{\&lambda;}}_k}{\mathrm{\&mu;}})),\frac{\mathrm{\&alpha;}}{\mathrm{\&mu;}})$

其中soft(·,Th)为阈值是Th的软阈值函数。

在经典的ALM算法中，惩罚参数μ是固定的。大量研究表明基于固定参数的ALM收敛速度非常慢，LALM也具有此特点。因此，选取动态变化的惩罚参数是很有必要的。采用持续策略来选取μ以此来加速LALM的收敛速度。采用持续策略式(19)可以变为

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-(y+\frac{{\mathrm{\&lambda;}}_k}{{\mathrm{\&mu;}}_k})),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k})---\left(20\right)$

持续惩罚参数μ可以定义为

μ_k+1=min{ημ_k,μ_max}    (21)

其中η>1为放大因子。

设 $d_k=y+\frac{{\mathrm{\&lambda;}}_k}{{\mathrm{\&mu;}}_k},$ 则有

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k})---\left(22\right)$

为了加速上述迭代的收敛速度，采用快速思想对变量Θ进行如下的二次更新

${\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k)---\left(23\right)$

其中 $t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2},$ t₀=1， $\frac{t_k-1}{t_{k+1}}\&Element;\left[\mathrm{0,1}\right]$ 为步长。

此外，使用下式更新Lagrange乘子

λ_k+1=λ_k-μ_k(AΘ_k+1-y)     (24)

又因为 $d_k=y+\frac{{\mathrm{\&lambda;}}_k}{{\mathrm{\&mu;}}_k},$ 故有

$d_k=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k)---\left(25\right)$

综上所述，本发明的FCLALM的算法步骤为：

步骤1、初始化：设置初始惩罚参数μ₀>0，最大惩罚参数μ_max>μ₀，参数α>0，初始变量t₀=1，放大因子η>1，初始稀疏系数Θ₀，初始辅助变量

初始乘子d₀，初始迭代次数k=0；

步骤2、更新辅助变量 $\stackrel{\&OverBar;}{\mathrm{\&Theta;}}:{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k});$

步骤3、更新变量 $t:t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤4、更新稀疏系数 $\mathrm{\&Theta;}:{\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k);$

步骤5、更新惩罚参数μ：μ_k+1=min{ημ_k,μ_max}；

步骤6、更新乘子 $d:d_{k+1}=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k);$

步骤7、若满足终止条件则算法终止，否则，令k=k+1转步骤2。

通过上述算法即可重构出融合图像的稀疏系数Θ，对稀疏系数Θ进行反离散小波变换（IDWT）即可得到最终的融合图像。

本发明方法的基本流程如图1所示。

FCLALM算法的复杂度主要在步骤2和步骤5，其计算量为O(MN)，而步骤3、4、6的计算量均仅为O(1)。FCLALM算法的优点在于它可以直接对问题(3)进行求解，而其他算法仅解决的是该问题的近似问题。此外，线性化、持续、快速的思想可以充分加速该算法的收敛速度。因此，FCLALM算法比其他重构图像的算法具有更好的性能和更快速的收敛速度。

接下来，我们将通过LALM、FPC(Fixed-Point Continuation)和FISTA(Fast IterativeShrinkage/Thresholding Algorithm)的收敛性给出FCLALM算法的收敛性。

引理1设

其中ρ_max(A^TA)表示矩阵A^TA的最大特征值。对于任意固定的μ>0，给定初始迭代点Θ₀，LALM算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。

引理2对于任意的μ_max>μ₀>0，给定初始迭代点Θ₀，FPC算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。

引理3给定初始迭代点Θ₀，FISTA算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。

定理1设

对于任意的μ_max>μ₀>0，给定初始迭代点Θ₀，FCLALM算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。

通过引理1，2和3很容易证明定理1。因优化问题(14)是凸优化问题，所以该问题的解若存在，则唯一。由引理1和式(19)可以得到，LALM算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。通过引理2和3可以得到，基于持续和快速思想的LALM算法也收敛于Θ^*。因此，FCLALM算法生成的序列{Θ_k}收敛于优化问题(14)的解Θ^*。

为了验证本发明所提出的FCLALM算法的性能，进行了以下实验。实验对象为：Lab(512×512)、Pepsi(512×512)和Clock(256×256)。首先将n×n的图像转化为N×1的向量，然后将N×1的向量分成n段，每段的长度为n。压缩比定义为

其中n为每段向量的维数，m为其相应的观测向量的维数。

将本发明的FCLALM重构算法与LALM、FPC、FISTA和OMP(Orthogonal MatchingPursuit)4种经典的重构算法进行了比较，并用仿真实验验证了FCLALM算法的有效性和普适性。

FCLALM算法的终止准则为

$\frac{{\left|\right|{\mathrm{\&Theta;}}_{k+1}-{\mathrm{\&Theta;}}_k\left|\right|}_2}{{\left|\right|{\mathrm{\&Theta;}}_k\left|\right|}_2}<\mathrm{\&epsiv;}---\left(26\right)$

其中ε>0是一个预先给定的较小的常量。

FCLALM算法的参数设置如表1所示：

表1FCLALM算法的参数设置

其中0表示零向量，即其所有元素均为零。

理论上，惩罚参数μ越大，重构图像的性能越好。但在实际中，当

时，重构图像的PSNR将是一个常数不再变化。因此，设

${\mathrm{\&mu;}}_{\max }={10}^8\&CenterDot;\frac{1}{{\mathrm{\&rho;}}_{\max }\left(A^{T}A\right)}.$ 假设 $0<\mathrm{\&alpha;}\&le;\frac{1}{{\mathrm{\&rho;}}_{\max }\left(A^{T}A\right)}$ 是为了保证FCLALM算法的收敛，此外，α越大，FCLALM算法的收敛速度越快，因此我们取

取η=2，取

Θ₀=0,λ₀=0。因为令

故取d₀=y。较小的阈值(例如ε=10^-4)并不一定能改进算法的重构精度，但却会到来迭代次数的增加，所以我们取ε=10^-3。

基于DCT基和DWT基的FCLALM算法重构3组参考图像在不同正交基下的PSNRs如图2所示。很显然，在相同压缩比下，基于DWT基的FCLALM算法重构图像的PSNRs值大于基于DCT基的FCLALM算法重构图像的PSNRs值。因此，本发明优选DWT基作为图像的稀疏正交基。

为了进一步验证FCLALM算法的性能，采用PSNR和迭代次数作为衡量指标。

各种重构算法在不同压缩比下重构Lab图像的PSNR如表2所示

表2重构Lab图像的PSNR

由表2可以看出，各种重构算法均可以较好的重构Lab图像。特别地，当压缩比r≥0.5时，所有的PSNR值均超过40dB。此外，随着压缩比数值的增大，所有算法重构图像的PSNR值均增大。也就是说，通过增大观测数目，可以提高重构图像的质量。最重要的是，在相同压缩比条件下，FCLALM算法重构图像的质量高于LALM、FPC、FISTA和OMP重构图像的质量。例如，当压缩比r=0.5，FCLALM重构图像的PSNR值为43.527dB，而LALM、FPC、FISTA和OMP重构图像的PSNR值分别为42.631dB、42.407dB、42.012dB、40.578dB。此时，FCLALM重构图像的PSNR值比OMP重构图像的PSNR值高大约3dB。特别地，当压缩比r=1时，提出的FCLALM重构图像的PSNR值比OMP重构图像的PSNR值高10.892dB。由表2可以看出，在相同压缩比条件下，FCLALM算法的性能优于LALM、FPC、FISTA和OMP算法的性能。

各种重构算法在不同压缩比下重构Lab图像的迭代次数如表3所示。

表3重构Lab图像的迭代次数

由表3可以看出，随着压缩比数值的增大，所有重构算法的迭代次数均增大。而且，在相同压缩比条件下，FCLALM算法的迭代次数小于LALM、FPC、FISTA和OMP算法的迭代次数。也即，此时，FCLALM算法的收敛速度快于其他4种算法的收敛速度。

为了验证FCLALM算法具有普适性，将该算法用来重构Pepsi参考图像，并与LALM、FPC和FISTA算法进行了比较。各种重构算法重构图像的PSNR和迭代次数分别如图3a、图3b所示。

因为重构图像的PSNR和迭代次数都与压缩比有关，有必要选择一种合适的压缩比，本实验取压缩比r=0.5。

为了验证本发明所提出的ALEM融合准则的效果，将本发明的ALEM融合准则与经典的MS、SWV、SDWV和EMWV融合准则进行了比较。重构算法采用本发明的FCLALM算法，压缩比为r=0.5，实验结果分别如图4a～图4h、图5a～图5h和图6a～图6h所示，其中，图4a～图4h依次为Lab参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像；图5a～图5h依次为Pepsi参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像；图6a～图6h依次为Clock参考图像、两幅不同聚焦的参考图像，以及分别采用ALEM、MS、SWV、SDWV和EMWV融合准则所得到的融合图像。

由图可以看出，压缩比r=0.5时，基于各种融合准则的融合图像均较好。但是，ALEM融合图像保留了参考图像的更多信息，更加接近于原参考图像。此外，ALEM融合图像比MS、SWV、SDWV和EMWV融合图像更清晰。也就是说，从视觉效果上来看，ALEM融合图像的效果最好。

为了定量地评价提出的图像融合与重构装置的性能，采用PSNR、扭曲度(WarpingDegree,WD)、平均梯度(Average Gradient,AG)和空间频率(Spatial Frequency,SF)作为评价指标。实验结果如表4所示。

表4融合结果的定量评价

由表4可以看出，对于评价指标PSNR和WD，ALEM融合准则的性能好于MS、SWV、SDWV和EMWV融合准则的性能。对于评价指标AG和SF，ALEM融合准则的性能稍好于MS、SWV、SDWV和EMWV融合准则的性能。从而有，提出的ALEM融合准则比其他4种经典的基于CS的图像融合准则具有更好的性能。

Claims (6)

1.一种基于压缩感知的多聚焦图像融合方法，首先对待融合的两幅图像分别进行压缩采样，得到两幅图像的观测向量；对两幅图像的观测向量进行融合，得到融合观测向量；利用融合观测向量重构出待融合的两幅图像的融合图像；其特征在于，所述对两幅图像的观测向量进行融合，具体按照以下方法：

首先将待融合的两幅图像的观测向量y₁=(u₁,u₂,…,u_M)、y₂=(v₁,v₂,…,v_M)分段表示为y₁=(u₁,u₂,…,u_P)^T和y₂=(v₁,v₂,…,v_P)^T；其中，u_j=(u_(j-1)·q+1,u(_j-1)·q+2,…,u_j·q)^T，v_j=(v_(j-1)·q+1,v_(j-1)·q+2,…,v_j·q)^T，j=1,2,…P,P×q=M；

然后按照下式分别计算分段向量u_j和v_j的能量匹配度D_j：

$D_j=\frac{2E_{12j}}{E_{1j}+E_{2j}},$

式中，E_1j、E_2j分别为u_j、v_j的局部能量测度，E_12j为u_j和v_j的混合能量测度，根据下式得到： $E_{12j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}\left|u_j{v}_j\right|/q;$

根据能量匹配度D_j与预设阈值T∈(0.5,1)的比较结果，按照以下方法计算融合观测向量y=(w₁,w₂,…,w_P)^T的分段向量w_j，j=1,2,…P：

若D_j<T，则：

$w_j=\left\{\begin{array}{cc}{u}_j& \mathrm{if}{E}_{1j}\&GreaterEqual;E_{2j}\\ v_j& \mathrm{if}{E}_{1j}<E_{2j}\end{array}\right.(j=\mathrm{1,2},...,P),$

若D_j≥T，则：

w_j=ωu_j+(1-ω)v_j(j=1,2,…,P)

其中ω为自适应加权因子，按照下式计算：

$\mathrm{\&omega;}=\frac{E_{1j}}{E_{1j}+E_{2j}};$

最终得到融合观测向量y=(w₁,w₂,…,w_P)^T。

2.如权利要求1所述基于压缩感知的多聚焦图像融合方法，其特征在于，在利用融合观测向量重构出待融合的两幅图像的融合图像时，所述融合图像的稀疏系数Θ通过求解以下的l₁-最小化问题得到：

$\begin{array}{cccc}\underset{\mathrm{\&Theta;}}{\min }& {\left|\right|\mathrm{\&Theta;}\left|\right|}_1& s.t.& y=\mathrm{A\&Theta;}\end{array},$

式中，y为融合观测向量，A为对所述压缩采样所使用的观测矩阵；

求解采用快速持续线性增广Lagrangian法，具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数μ₀>0，最大惩罚参数μ_max>μ₀，参数α>0，初始变量t₀=1，放大因子η>1，初始稀疏系数Θ₀，初始辅助变量

初始乘子d₀，初始迭代次数k=0；

步骤2、更新辅助变量 $\stackrel{\&OverBar;}{\mathrm{\&Theta;}}:{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k});$

步骤3、更新变量 $t:t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤4、更新稀疏系数 $\mathrm{\&Theta;}:{\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k);$

步骤5、更新惩罚参数μ：μ_k+1=min{ημ_k,μ_max}；

步骤6、更新乘子 $d:d_{k+1}=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k);$

步骤7、若满足终止条件则算法终止，否则，令k=k+1转步骤2。

3.如权利要求1或2所述基于压缩感知的多聚焦图像融合方法，其特征在于，所述压缩采样利用离散小波变换进行稀疏表示；在利用融合观测向量重构出待融合的两幅图像的融合图像时，通过对融合图像的稀疏系数进行反离散小波变换得到。

4.如权利要求1或2所述基于压缩感知的多聚焦图像融合方法，其特征在于，所述压缩采样采用随机高斯矩阵进行观测。

5.一种基于压缩感知的多聚焦图像融合装置，包括压缩采样模块、融合模块以及图像重构模块；所述压缩采样模块用于对待融合图像进行压缩采样，得到待融合图像的观测向量；所述融合模块用于对两幅待融合图像的观测向量进行融合，得到融合观测向量；所述图像重构模块用于对融合观测向量进行图像重构，得到融合图像；其特征在于，所述融合模块对两幅待融合图像的观测向量进行融合，具体按照以下方法：

首先将待融合的两幅图像的观测向量y₁=(u₁,u₂,…,u_M)、y₂=(v₁,v₂,…,v_M)分段表示为y₁=(u₁,u₂,…,u_P)^T和y₂=(v₁,v₂,…,v_P)^T；其中，u_j=(u_(j-1)·q+1,u_(j-1)·q+2,…,u_j·q)^T，v_j=(v_(j-1)·q+1,v_(j-1)·q+2,…,v_j·q)^T，j=1,2,…P,P×q=M；

然后按照下式分别计算分段向量u_j和v_j的能量匹配度D_j：

$D_j=\frac{2E_{12j}}{E_{1j}+E_{2j}},$

式中，E_1j、E_2j分别为u_j、v_j的局部能量测度，E_12j为u_j和v_j的混合能量测度，根据下式得到： $E_{12j}=\underset{(j-1)\&CenterDot;q+1}{\overset{\mathrm{jq}}{\mathrm{\&Sigma;}}}\left|u_j{v}_j\right|/q;$

根据能量匹配度D_j与预设阈值T∈(0.5,1)的比较结果，按照以下方法计算融合观测向量y=(w₁,w₂,…,w_P)^T的分段向量w_j，j=1,2,…P：

若D_j<T，则：

$w_j=\left\{\begin{array}{cc}{u}_j& \mathrm{if}{E}_{1j}\&GreaterEqual;E_{2j}\\ v_j& \mathrm{if}{E}_{1j}<E_{2j}\end{array}\right.(j=\mathrm{1,2},...,P),$

若D_j≥T，则：

w_j=ωu_j+(1-ω)v_j(j=1,2,…,P)

其中ω为自适应加权因子，按照下式计算：

$\mathrm{\&omega;}=\frac{E_{1j}}{E_{1j}+E_{2j}};$

最终得到融合观测向量y=(w₁,w₂,…,w_P)^T。

6.如权利要求5所述基于压缩感知的多聚焦图像融合装置，其特征在于，所述图像重构模块通过快速持续线性增广Lagrangian法求解以下的l₁-最小化问题得到融合图像的稀疏系数Θ：

$\begin{array}{cccc}\underset{\mathrm{\&Theta;}}{\min }& {\left|\right|\mathrm{\&Theta;}\left|\right|}_1& s.t.& y=\mathrm{A\&Theta;}\end{array},$

式中，y为融合观测向量，A为对所述压缩采样所使用的观测矩阵；

具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数μ₀>0，最大惩罚参数μ_max>μ₀，参数α>0，初始变量t₀=1，放大因子η>1，初始稀疏系数Θ₀，初始辅助变量

初始乘子d₀，初始迭代次数k=0；

步骤2、更新辅助变量 $\stackrel{\&OverBar;}{\mathrm{\&Theta;}}:{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}=\mathrm{soft}({\mathrm{\&Theta;}}_k-\mathrm{\&alpha;}{A}^T(A{\mathrm{\&Theta;}}_k-d_k),\frac{\mathrm{\&alpha;}}{{\mathrm{\&mu;}}_k});$

步骤3、更新变量 $t:t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤4、更新稀疏系数 $\mathrm{\&Theta;}:{\mathrm{\&Theta;}}_{k+1}={\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_{k+1}-{\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}_k);$

步骤5、更新惩罚参数μ：μ_k+1=min{ημ_k,μ_max}；

步骤6、更新乘子 $d:d_{k+1}=y-\frac{{\mathrm{\&mu;}}_k}{{\mathrm{\&mu;}}_{k+1}}(A{\mathrm{\&Theta;}}_{k+1}-d_k);$

步骤7、若满足终止条件则算法终止，否则，令k=k+1转步骤2。

Images