基于非局部全变差和低秩稀疏的高光谱压缩感知重建方法
【技术领域】
[0001 ]本发明属于图像处理技术领域，主要涉及一种高光谱数据压缩感知重建方法，可 用于高光谱成像。
【背景技术】
[0002] 高光谱数据由上百个非常窄的波段组成，其空间和波段方向的高分辨率特性使得 高光谱数据具有很大的数据维度，庞大的信息量对高光谱数据的存储、传输以及后续处理 带来了困难。传统的压缩采样方法使用奈奎斯特采样速率对信号进行均匀或非均匀采样 后，通过预测、变换、矢量量化等算法进行压缩。这种高冗余采样再压缩的过程造成了极大 的资源浪费，给低功耗和资源有限的机载或星载应用带来巨大压力。
[0003] 压缩感知作为一种新颖的信号获取理论，融合了传统的采样与压缩过程，以远低 于奈奎斯特采样率的方式直接获取测量数据，降低了采样成本，减少了存储资源。压缩感知 方法有效地解决了高光谱数据计算和存储量大的问题。在现有的高光谱数据压缩感知重建 方法中，Wagadarikar A等人提出将全变差模型独立应用于每一个波段的数据，用来约束空 间平滑性，但该方法忽略了谱间的相关性。Aravind N V等人以分布式压缩感知理论为基础 提出了Simultaneous 0ΜΡ方法，该方法利用了高光谱数据的联合稀疏特性，虽然在重建时 间上具有一定优势，但重建图像的质量较低，难以达到后续应用的要求。
[0004] James E·Fowler提出的基于压缩-投影主成分分析CPPCA算法将主成分分析PCA算 法用于高光谱数据的压缩重建。该算法首先将原高光谱数据投影到随机选择的低维子空间 获得采样数据，然后对采样数据进行PCA变换。通过采样数据的PCA变换基，利用凸集交替投 影优化算法近似恢复出原数据的PCA变换系数和变换基，进而得到重建的高光谱数据。该方 法虽然能将PCA降维步骤由编码端转移到解码端，并且重建时间较短，但是重建图像的质量 较低，不足以为后续的图像处理提供可靠的数据来源。
【发明内容】
[0005] 本发明的目的在于针对上述已有技术的不足，提出一种基于非局部全变差和低秩 稀疏的高光谱压缩感知重建方法，以提高高光谱数据压缩采样后的重建质量，满足后续应 用的要求。
[0006] 为实现上述目的，本发明的技术方案包括如下：
[0007] (1)输入高光谱数据4' = {為制，^1",.开，,/ = 1，."，户,1 = 1，".3}￡]^必'，其中 ?表示实数空间，Η和P分别表示高光谱数据在空间域的水平和垂直方向像元点数，N表示波 段总数，xn= {Xi,j,n，i = l，. . .，H, j = l，. . .，Ρ}表示高光谱数据Xori的第η个波段，Xi,j = {Xi,j,l, . . .，Xi,j,N}表示高光谱数据的一个像元;将高光谱数据Xcxri进行向量化处理，得到向 量化的高光谱数据I e肢?~;
[0008] (2)对数据使用块对角采样矩阵φ进行采样，得到采样数据FeIR7?''1， 其中T表示单个高光谱像元的采样数据维度；
[0009] (3)对采样数据Υ桉下式讲行初始重津：
[0010]
[0011]其中，Xl表不重建数据的低秩分量，Xs表不重建数据的稀疏分量，Χι/表不初始重建 数据X'的低秩分量，Xs'表示初始重建数据X'的稀疏分量，mat( ·)表示取向量的矩阵形式， I I · I I*表示取矩阵核范数的操作，I I · I |ι表示取向量1范数的操作，μι表示约束低秩项的 参数，μ2表示约束稀疏项的参数，arg min表示取使目标函数达到最小值时的自变量值；
[0012] 按照上式计算得到初始重建数据f，Xr为向量XJ +Χ^的三维张量形式；
[0013] (4)利用K均值聚类法，对初始重建的高光谱数据X'进行聚类操作，得到Μ类数据 Xclu(l)，Xclu(2)，· · ·，Xclu(m)，· · ·，Xclu(M)，中Xclu(m)表弟Π 1 类数据，Π 1 - 1，· · ·，Μ;
[0014] (5)根据像元Xw的类别，对采样数据Υ进行分类，得到Μ类采样数据Υ⑴，Υ⑵，…， Y(m), · · ·，Υ(Μ)，其中Y(m)表示第m类像元的采样数据，m=l, · · ·，Μ;
[0015] (6)根据向量化的高光谱数据X、块对角采样矩阵φ以及采样数据Υ，构造二次重建 模型如下：
[0016]
[0017]
[0018] 其中，Xl(m)表不X(m)的低秩分量，X2(m)表不X(m)的稀疏分量，X(m)表不第m类的重建数 据，m=l, ...，M，M表示聚类总数，X表示估计高光谱数据，Y表示采样数据，Φ表示块对角采 样矩阵，▽?:尤表不X中第η个波段Xn的非局部全变差，X#表不二次重建的最优高光谱数据， μι表示约束低秩项的参数，μ2表示约束稀疏项的参数，μ表示约束非局部全变差项的参数。
[0019] (7)根据步骤(5)得到的各类采样数据Υω，对步骤(6)中的二次重建模型进行迭代 求解，并将二次重建的最优高光谱数据浐作为最终的重建数据。
[0020] 本发明与现有技术相比，具有以下优点：
[0021] 第一，本发明考虑到高光谱数据的非局部结构特性，将非局部全变差引入低秩稀 疏重建模型，不仅利用了高光谱较强的谱间相关性，而且结合数据中的空间非局部信息，克 服了传统的TV重建易造成变化剧烈的边缘处重建模糊的缺点，使得图像中细节和边缘的重 构效果得到提升。
[0022] 第二，本发明考虑到高光谱数据中同类像元间的高相似度，引入聚类算法，将低秩 稀疏重建模型应用于同类的像元中，大大提高了重建结果的信噪比，使重建效果得到明显 提尚。
【附图说明】
[0023]图1是本发明的实现流程图；
[0024]图2是本发明仿真所用高光谱数据IndianPines第10波段图像；
[0025]图3是本发明与现有算法对高光谱数据IndianPines在采样率分别为0.2、0.3、0.4 时进行重建得到的信噪比曲线图；
[0026] 图4是本发明与现有算法对高光谱数据IndianPines在采样率为0.2时进行重建得 到的第10波段图像。
【具体实施方式】
[0027] 下面结合附图，对本发明作进一步的详细描述。
[0028]参照附图1，本发明的实现步骤如下：
[0029] 步骤1，输入高光谱数据XOTi，得到向量化的高光谱数据X。
[0030] 输入高光谱数据= = 1,…丑，i = 1S,..，/U 二 1，…e RiixPxA',其中 R 表 示实数空间，Η和P分别表示高光谱数据在空间域的水平和垂直方向像元点数，N表示波段总 数，Xn={Xi,j,n，i = l，· · ·，H，j = l, · · ·，Ρ}表不尚光谱数据Xori 的弟II个波段，Xi , j = {Xi,j,l, . . .，Xi,j,N}表示高光谱数据的一个像元；
[0031 ]将高光谱数据XOT1的各像元逐个摞列，得到向量化的高光谱数据X 。
[0032]步骤2，对向量化的高光谱数据X采样。
[0033]按照下式，对向量化的高光谱数据X使用块对角采样矩阵Φ进行采样，得到采样数 据Y:
[0034] γ = ΦΧ ,
[0035] 其中，块对角采样矩阵φ是由HP个像元的采样矩阵Φ2, . . .，〇i，. . .，Φηρ构 成，其表示为：
[0036]
:，Φ i表示第i个像元所用的采样矩阵，其大小为Τ X Ν 的随机高斯矩阵，i的取值范围为U，. . .，HP}，采样数据Y是大小为THPX1的向量，T表示单 个高光谱像元的采样数据维度。
[0037]步骤3，对采样数据Y进行初始重建。
[0038]本步骤的实现是通过对如下所示的低秩稀疏重建模型求解进行：
[0039]
[0040]其中，Xl表不重建数据的低秩分量，Xs表不重建数据的稀疏分量，Χι/表不初始重建 数据X'的低秩分量，Xs'表示初始重建数据X'的稀疏分量，mat( ·)表示取向量的矩阵形式， I I · I I*表示取矩阵核范数的操作，I I · I |ι表示取向量1范数的操作，μι表示约束低秩项的 参数，μι的值设置为1500，μ2表示约束稀疏项的参数，μ2的值设置为500，arg min表示取使目 标函数达到最小值时的自变量值；
[0041 ]具体的初始重建步骤如下：
[0042] 3a)设置迭代次数初始值k=l，最大迭代次数kmax = 5,分别初始化低秩分量X〗、稀 疏分量，以及中间变量Yk;
[0043] 3b)利用下式对中间变量Yk进行更新：
[0044]
[0045] 其中，Yk+1表示第k+Ι次迭代得到的中间变量，k表示用于迭代求解XL和Xs所用的迭 代次数，xf和I】分别表示第k次迭代得到的低秩分量和稀疏分量；
[0046] 3c)利用下式对低秩分量和稀疏分量进行更新：
[0047]
[0048] 3 d )判断是否满足迭代停止条件：若k + 1 = k m a x，停止迭代，此时 ，X/ 初始重建数据X'为向量Xl'+Xs'的三维张量形式;否则，k = k+l，返回 3b) 〇
[0049] 步骤4,利用K均值聚类法，对初始重建的高光谱数据X'进行聚类操作，得到Μ类数 据。
[0050] 4a)从初始重建数据X'中随机选取Μ个像元点作为初始的聚类中心；
[0051] 4b)从f中任意选取一个未聚类的像元点，分别计算所选取的像元点与Μ个聚类中 心的欧氏距离，找出欧式距离最小值对应的聚类中心，将所选取出的像元点与欧氏距离最 小值对应的聚类中心作为同一类数据；
[0052] 4c)判断初始重建数据f中的全部像元点是否都完成聚类，若是，执行4d)，否则， 返回4b);
[0053] 4d)计算聚类后每一类像元点的均值，并将该均值作为更新后的聚类中心；
[0054] 4e)按照下式，分别计算出Μ个聚类中心更新前后的残差：
[0055] Cresi= | | fi_hi | 12
[0056] 其中，CreSl表示第i类聚类中心更新前后的残差，i的取值范围为{1，2，. ..，M}，M 表示聚类总数，Μ · I |2表示取向量2范数操作，h表示第i类更新后的聚类中心，lu表示第i 类更新前的聚类中心；
[0057] 4f)判断计算出的Μ个聚类中心更新前后的残差中最大值是否小于预设阈值 0.001，若是，则聚类过程结束，得到Μ类数据Xca u⑴，Xcau(2)，. . .，Xcau(m)，. . .，XcauW，且 = jr，七表示取集合并操作，m=l，. . .，M;否则，将更新后的Μ个聚类中心作为新的 m-1 聚类中心返回4b)。
[0058] 步骤5,根据像元的类别，对采样数据Y进行分类。
[0059] 根据Μ类数据Xcau(l)，Xcdu(2)，· · ·，XcdU(m), · · ·，XcdU(M)中像元Xi,j的类别，确定采样数 据Y中对应空间位置上的像元Xu的采样数据类别；
[0060] 将采样数据Y中属于同一类的采样数据进行组合，得到Μ类采样数据Y⑴，Y⑵，...， Y(m), · · ·，Υ(Μ)，其中Y(m)表示第m类像元的采样数据，m=l, · · ·，Μ。
[0061]