基于混合范数追踪算法的高光谱图像重构方法

摘要

本发明公开了一种基于混合范数追踪算法的高光谱图像重构方法：结合高光谱图像性质分析现有压缩感知重构算法，选择相应的0范数算法和1范数算法，证明存在能将二者结合的混合策略并寻找最佳的混合策略，得到混合范数追踪算法，进而对稀疏采样的高光谱图像进行重构，最后输出经稀疏逆变换获得的重构高光谱图像。该方法具有人机交互接口模块、混合策略分析模块、高光谱图像重构模块、重构图像输出模块这四个功能模块。混合范数追踪算法相比于传统的重构方法，在速度和精度上都有明显的改善，减轻了高光谱大数据传输、存储带来的硬件压力，提高了数据传输的抗干扰能力。

说明书

技术领域

本发明涉及一种基于混合范数追踪算法的高光谱图像重构方法，适用于高光谱数据处理系统中，属于高光谱数据处理领域。 

背景技术

高光谱遥感是一项近几十年来迅速发展的空天对地观测技术，它图谱合一的性质有重大的研究意义和广泛的应用前景，已经在商业、军事和民间领域得到广泛应用。相比较于常规遥感图像，高光谱图像能够更精细的刻画地物特征，探测到常规遥感中无法探测到的物质，为后期地物分类以及探测提供前提条件。然而，它良好的性质建立在庞大的数据量上，一张高光谱图像的大小往往是常规图像的几百倍，大量的波段数带来了巨大的信息冗余。这为高光谱图像的采集、传输和存储带来不必要的麻烦。 

压缩感知是近年来提出的一个寻找欠定系统稀疏解的技术，该技术通过开发信号的稀疏特性在远小于奈奎斯特采样率的条件下，利用随机采样获取信号的离散样本，然后通过非线性重构算法完美重构信号。压缩感知理论一经提出便受到信息论、图像处理、地球科学、光学、微波成像、模式识别、无线通讯等领域的高度关注，并被广泛的应用于电子工程尤其是信号处理中，进行低采样率采样和高精度重构稀疏或可压缩的信号。目前广泛采用的压缩感知重构算法主要有两类：基于0范数的贪婪算法和基于1范数的凸优化算法。基于0范数的贪婪算法包括匹配追踪算法、正交匹配追踪算法、硬阈值迭代算法等，基于1范数的凸优化算法包括基追踪算法、梯度投影算法、内点法、同伦法等。这两类算法各有所长但也均有其不足之处：0范数对应的算法速度快但是精度不高；1范数对应的算法精度高但速度不快。 

高光谱图像数据庞大且获取不易，为了减轻硬件传输、存储压力需尽可能低的进行采样和高精度重构；另一方面，低速度的重构算法的缺点会因大数据量而被急剧的放大，造成的时间消耗往往让人无法忍受。因此，需要寻找一种重构精度足够高并且重构速度足够快的算法，满足高光谱图像实时处理的需求。 

发明内容

本发明旨在提供一种基于混合范数追踪算法的高光谱图像重构方法，具体是一种结合现有0范数算法和1范数算法优点的新型压缩感知重构方法。本方法具有较强的鲁棒性，能够 在保证重构精度的情况下快速的重构出原始的高光谱图像，对不同地物环境的高光谱图像均表现出良好的实验效果。 

本发明所涉及的方法流程具体包括以下四个步骤：1、获得初始信息及相关初始化操作；2、混合策略分析；3、高光谱图像重构；4、重构结果输出。下面对该方法流程各步骤进行详细说明： 

步骤一获得初始信息及相关初始化操作 

利用人工交互接口模块输入高光谱图像稀疏采样结果和两类算法范例，并设置如公式(1)到(3)所示的算法迭代终止条件，对是否终止算法进行判断。 

|norm(x_k-1，0)-norm(x_k，0)|＜ε₀  (1) 

$\frac{\left|\right|x_k-x_{k-1}\left|\right|}{\left|\right|x_k\left|\right|}\le {ϵ}_1---\left(2\right)$

$\frac{\left|\right|f\left(x_k\right)-f\left(x_{k-1}\right)\left|\right|}{\left|\right|f\left(x_k\right)\left|\right|}\le {ϵ}_2---\left(3\right)$

其中，norm(x，0)代表x中非零元素个数，ε₀，ε₁，ε₂是选定的正常数，x是第k步迭代结果，f是目标函数。 

步骤二混合策略分析 

假设原始信号是k稀疏的n维信号，它的稀疏变换结果为x，求解的近似结果为x′。实际上，无论是基于0范数的算法还是基于1范数的算法，它们的迭代过程都是n维空间上逼近目标的一条折线。假定n维空间的基为

基于0范数的算法满足下面的约束优化问题： 

$\underset{x}{\min }{\left|\right|x\left|\right|}_{l_0},s.t.\mathrm{\Phi \Psi x}=y---\left(4\right)$

其中，y代表观测值，x是我们需要求解的结果，Φ是观测矩阵Ψ是稀疏矩阵。基于0范数的贪婪算法的每一步都会获得某个结果维度上的近似值，按照下面的步骤进行迭代： 

第1步，从原点出发沿坐标轴获取第一个分量的近似值

第2步，从x₁′出发沿坐标轴获取第二个分量的近似值

第n步，从x_n-1′出发沿坐标轴获取第n个分量的近似值

基于1范数的算法满足下面的凸优化问题： 

$\underset{x}{\min }{\left|\right|x\left|\right|}_{l_1},s.t.\mathrm{\Phi \Psi x}=y---\left(5\right)$

它可以利用拉格朗日方法写成如下形式： 

$\underset{x}{\min }({\left|\right|x\left|\right|}_{l_1}+\tau ||y-\mathrm{\Phi \Psi x}|\left|\right)---\left(6\right)$

其中，τ代表一个用来平衡信号稀疏度和重构误差的非负常量。凸优化算法的每一步迭代都会使迭代结果离目标更近一些。为了更直接的说明，利用梯度引导的算法作为一个例子，这类算法以下面的步骤进行着迭代过程： 

第1步，从原点出发，沿此时的梯度方向前进α⁰步长

第2步，从x₁′出发，沿此时的梯度方向前进α¹步长

第m步，从x_m-1′出发，沿此时的梯度方向前进α^m-1步长

其中，参数αⁱ代表第i步的步长，代表目标函数。 

基于步骤二的分析，将混合范数追踪算法确定为由0范数获得近似的结果，再由1范数算法进行调整得到精确结果的两步迭代过程。混合范数追踪算法有效的结合了0范数的速度和1范数的精度，能够高速准确的重构信号。现在存在的问题有两个：1、是否存在最优的混合策略；2、如何寻找最优混合策略。 

首先分析最优的混合策略是否存在，也就是说交换点t是否存在。定义0范数算法的速度时间函数为v₀(t)，1范数算法的速度时间函数为v₁(t)。为了分析方便，假设两类速度时间函数都足够光滑。 

如图2所示，将一次迭代的初始点到目标点的距离记做d₁，终止点到目标点的距离记做d₂，将这次迭代所走的路程记做d₁-d₂。那么这两个时间点之间的算法的速度时间函数满足公式(7)。 

${\int }_{t_1}^{t_2}v\left(t\right)\mathrm{dt}=d_1-d_2---\left(7\right)$

选择时间节点t作为两类算法的交替节点。在节点t之前，由0范数算法获得结果的近似值，之后由1范数算法对近似值进行调整，具体情况如图3所示。由于两类算法的速度不一样，相同时间里0范数算法前进的距离大于1范数算法，所以1范数对应的迭代初始时间应该大于t，记做f(t)。由前进距离相同可知，那么两者满足公式(8)所示的关系。 

${\int }_{t_0}^t{v}_0\left(t\right)\mathrm{dt}={\int }_{t_0}^{f\left(t\right)}{v}_1\left(t\right)\mathrm{dt}---\left(8\right)$

定义新算法的平均速度为那么它满足公式(9)所示的等式关系。 

$\overline{v\left(t\right)}=\frac{{\int }_{t_0}^t{v}_0\left(t\right)\mathrm{dt}+{\int }_{f\left(t\right)}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{t-t_0+t_{\mathrm{end}}-f\left(t\right)}---\left(9\right)$

联合公式(8)和公式(9)，可以得到下面的结果： 

$\overline{v\left(t\right)}=\frac{{\int }_{t_0}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{t_{\mathrm{end}}-t_0-(f\left(t\right)-t)}---\left(10\right)$

由公式(10)可知，新算法运行的距离为这意味着它能达到和1范数算法相同的精度。随着测试样本的稀疏度增加，0范数算法的迭代时间明显增加，但是1范数算法的时间变化并不明显。例如，当测试样本的稀疏度增加到200的时候，正交匹配追踪算法的CPU运行时间明显增加，但是梯度投影算法的CPU运行时间没有多大改变。两种算法的混合使得新算法对测试样本的CPU运行时间更加稳定。另一方面，大多数1范数算法对迭代的初始点有条件要求，但是0范数算法则很少有要求。利用0范数算法得到的近似值作为1范数算法的迭代初始值往往比随机选择的初始值更符合1范数算法的要求。终上所述，混合范数追踪算法较于原来的0范数算法和1范数算法具有更强的鲁棒性。 

由于新算法运行已经确定，关系新算法效率的仅仅是算法的时间节点选取。因此问题转化为下面的优化问题： 

$\underset{t}{\max }(f\left(t\right)-t)$

$s.t.{\int }_{t_0}^t{v}_0\left(t\right)\mathrm{dt}={\int }_{t_0}^{f\left(t\right)}{v}_1\left(t\right)\mathrm{dt}---\left(11\right)$

现在分别来证明公式(11)中的函数f(t)连续和最优混合策略在开区间(t₀，t_end)上存在。 

首先定义函数F(t)：$F\left(t\right)={\int }_{t_0}^t{v}_0\left(t\right)\mathrm{dt}={\int }_{t_0}^{f\left(t\right)}{v}_1\left(t\right)\mathrm{dt}$

对于选择开区间(t₀，t_end)中的任意一个点t_c，那么 

$\left|\right|F\left(t\right)-F\left(t_c\right)\left|\right|=\left|\right|{\int }_{t_c}^t{v}_0\left(t\right)\mathrm{dt}\left|\right|\le \left|\right|t-t_c\left|\right|*\underset{t}{\max }\left|\right|v_0\left(t\right)\left|\right|---\left(12\right)$

选取当||t-t_c||＜δ时，有||F(t)-F(t_c)||＜ε。也就是说F(t)在t_c点连续，由于t_c点的任意性，F(t)在开区间(t₀，t_end)上连续。 

所以，对于存在一个正常数δ＞0，当||t-t_c||＜δ时 

$\left|\right|F\left(t\right)-F\left(t_c\right)\left|\right|=\left|\right|{\int }_{f\left(t_c\right)}^{f\left(t\right)}{v}_1\left(t\right)\mathrm{dt}\left|\right|=\left|\right|f\left(t\right)-f\left(t_c\right)\left|\right|*\left|\right|v_1\left(\eta \right)\left|\right|<ϵ---\left(13\right)$

$\left|\right|f\left(t\right)-f\left(t_c\right)\left|\right|<\frac{ϵ}{\left|\right|v_1\left(\eta \right)\left|\right|}---\left(14\right)$

公式(13)中，η是一个常数满足η∈[f(t)，f(t_c)]。因此，f(t)在点开区间(t₀，t_end)上任意点t_c处连续，所以f(t)在开区间(t₀，t_end)处连续。 

对于优化问题(11)，函数f(t)-t在区间(t₀，t_end)上连续。因此，对于任意函数f(t)-t在区间[t₀+ε，t_end-ε]上都存在最大值。现在问题转换成函数f(t)-t在区间[t₀，t_end]上的最大值取值位置不在t₀和t_end处。 

由之前的分析知，0范数的速度比1范数快，这意味着存在一点t₀′，在这点之前0范数的速度都会大于1范数。与之对应的，1范数运行时间比0范数长，这意味着在0范数的运行停止时1范数速度仍然在运行。即存在一点t_end′，在该点之后1范数的速度大于0范数速度。那么有： 

$\overline{v\left({t_0}^{\prime }\right)}=\frac{{\int }_{t_0}^{{t_0}^{\prime }}-v_0\left(t\right)\mathrm{dt}+{\int }_{f\left({t_0}^{\prime }\right)}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{(t_{\mathrm{end}}-t_0)-(f\left({t_0}^{\prime }\right)-{t_0}^{\prime })}=\frac{{\int }_{t_0}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{(t_{\mathrm{end}}-t_0)-(f\left({t_0}^{\prime }\right)-{t_0}^{\prime })}>\frac{{\int }_{t_0}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{t_{\mathrm{end}}-t_0}=\overline{v\left(t_0\right)}---\left(15\right)$

$\stackrel{‾}{v\left({t_{\mathrm{end}}}^{\prime }\right)}=\frac{{\int }_{t_0}^{{t_{\mathrm{end}}}^{\prime }}-v_0\left(t\right)\mathrm{dt}+{\int }_{f\left({t_{\mathrm{end}}}^{\prime }\right)}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{(t_{\mathrm{end}}-t_0)-(f\left({t_{\mathrm{end}}}^{\prime }\right)-{t_{\mathrm{end}}}^{\prime })}=\frac{{\int }_{t_0}^{t_{\mathrm{end}}}{v}_1\left(t\right)\mathrm{dt}}{(t_{\mathrm{end}}-t_0)-(f\left({t_{\mathrm{end}}}^{\prime }\right)-{t_{\mathrm{end}}}^{\prime })}>\frac{{\int }_{t_0}^{t_{\mathrm{end}}}{v}_0\left(t\right)\mathrm{dt}}{t_{\mathrm{end}}-t_0}=\overline{v\left(t_{\mathrm{end}}\right)}---\left(16\right)$

即函数f(t)-t在区间[t₀，t_end]上的最大值取值位置不在t₀和t_end处。 

根据前面的结论，对最优混合策略的选取进行分析。对公式(11)进行求导，得 

[f(t)-t]′＝f′(t)-1＝0  (17) 

对公式(8)进行分析，有 

${\left[{\int }_{t_0}^t{v}_0\left(t\right)\mathrm{dt}\right]}^{\prime }=v_0\left(t\right)---\left(18\right)$

${\left[{\int }_{t_0}^{f\left(t\right)}{v}_1\left(t\right)\mathrm{dt}\right]}^{\prime }=v_1\left[f\left(t\right)\right]f^{\prime }\left(t\right)=v_1\left[f\left(t\right)\right]---\left(19\right)$

因此，可以得到： 

v₀(t)＝v₁[f(t)]  (20) 

所以，效率最高的时刻的必要条件是运行相同路程时速度相同。也就是说，运行相同路程是速度相同的点有很大可能就是最优的混合策略点。 

步骤三重构高光谱图像 

通过步骤二的分析，可以整理得到混合范数追踪算法的算法流程。为了方便理解，用梯度引导的算法作为1范数算法的例子： 

步骤1，利用0范数算法获取待求解信号的第一个分量的近似值

步骤2，利用0范数算法获取待求解信号的第二个分量的近似值

步骤p，利用0范数算法获取待求解信号的第p个分量的近似值

步骤p+1，利用1范数算法沿当前位置的梯度方向获取求解信号的位置 ${x_{p+1}}^{\prime }={x_p}^{\prime }-{\alpha }^{p+1}▿f\left({x_p}^{\prime }\right);$

步骤q，利用1范数算法沿当前位置梯度方向获取求解信号的位置 $x^{\prime }={x_q}^{\prime }={x_{q-1}}^{\prime }-{\alpha }^q▿f\left({x_{q-1}}^{\prime }\right).$

在上面的算法流程中p∈{1，2，…k}，k是原始信号的稀疏度，αⁱ是第i步迭代步进；新算法将会很大概率的达到最优状态，如果利用0范数算法迭代p步花的时间t满足(8)和(20)。 

利用混合范数追踪算法对输入的低采样率高光谱图像数据进行重构，按照输入的终止条件终止迭代步骤，获得需要的稀疏近似解。根据高光谱图像的特点，选择相应的稀疏基，通过混合策略分析模块得到的混合范数重构算法对人机交互接口模块获取的稀疏采样结果进行重构，利用稀疏逆变换得到高光谱图像的重构结果。 

步骤四重构结果输出 

通过重构结果输出模块，输出重构完成的高光谱图像。 

本发明是基于混合范数追踪算法的高光谱图像重构方法，其优点在于：能够在较短的时间内对低采样率的高光谱图像样本进行准确的重构，对于不同背景的地物目标和不同的采样率均能达到良好的重构效果，较原有的算法具有更高的鲁棒性。 

附图说明

图1所示为基于混合范数追踪算法的高光谱图像重构方法流程图 

图2三类重构算法比较 

图3所示为速度的定义 

图4所示为三类算法的速度时间关系 

图5所示为仿真实验数据 

图6所示为仿真实验结果 

具体实施方式

下面结合附图与仿真实验对该发明的技术方法进行进一步的说明。 

基于本发明开发了仿真原型系统，该系统包括：人机交互接口模块、混合策略分析模块、高光谱图像重构模块、重构图像输出模块这四个功能模块。 

通过人机交互接口模块获得稀疏采样后的高光谱数据。本实例采用Indiana Pine高光谱数据和Washington D.C.Mall高光谱数据截取的一部分，具体情况如图4所示。其中，Indiana Pine高光谱数据大小为145×145，波长范围为400～2400nm，除去水汽吸收波段和低信噪比波段后，保留220个波段，Washington D.C.Mall高光谱数据截取的一部分，大小为100×100，波长范围为0.4～2.4μm，除去水汽吸收波段和低信噪比波段后，保留191个波段，采样矩阵为Gaussian随机矩阵。 

对有关迭代终止条件的参数进行初始化处理，设置用来终止算法的参数ε₀，ε₁，ε₂，当算法满足公式(1)(2)(3)时，停止运行。 

结合高光谱图像稀疏采样结果，选择合适的0范数算法和1范数算法，分析它们的具体算法流程，寻找它们的最优混合时间节点，获得混合范数追踪算法的算法流程。 

假设原始信号是k稀疏的n维信号，它的稀疏变换结果为x，求解的近似结果为x′。实际上，无论是基于0范数的算法还是基于1范数的算法，它们的迭代过程都是n维空间上逼近目标的一条折线。假定n维空间的基为为了方便解释，以下是由梯度引导的算法作为1范数算法的例子： 

步骤1：利用0范数算法获取待求解信号的第一个分量的近似值：

步骤2：利用0范数算法获取待求解信号的第二个分量的近似值：

步骤p：利用0范数算法获取待求解信号的第p个分量的近似值：

步骤p+1：利用1范数算法沿当前位置的梯度方向获取求解信号的位置： ${x_{p+1}}^{\prime }={x_p}^{\prime }-{\alpha }^{p+1}▿f\left({x_{p+1}}^{\prime }\right);$

步骤q：利用1范数算法沿当前位置梯度方向获取求解信号的位置： $x^{\prime }={x_q}^{\prime }={x_{q-1}}^{\prime }-{\alpha }^{q-1}▿f\left({x_{q-1}}^{\prime }\right);$

在上面的算法流程中p∈{1，2，…k}，k是原始信号的稀疏度，αⁱ是第i步迭代步进；新算法将会很大概率的达到最优状态，如果利用0范数算法迭代p步花的时间t满足(8)和(20)。 

根据高光谱图像的特点，选择相应的稀疏基，通过混合策略分析模块得到的混合范数重构算法对人机交互接口模块获取的稀疏采样结果进行重构，利用稀疏逆变换得到高光谱图像的重构结果； 

通过重构结果输出模块，输出重构完成的高光谱遥感图像。 

本发明方法经过仿真系统的具体实施，实验结果表明能够高速高精度重构出不同采样率下不同背景的地物高光谱图像，具有良好的鲁棒性。