基于正交匹配追踪的超光谱图像预测压缩方法与流程

1.本发明涉及图像处理技术领域，具体为基于正交匹配追踪的超光谱图像预测压缩方法。


背景技术：

2.超光谱大气探测为气象领域提供了重要的大气资料。探测技术的提高使得光谱分辨率极大提升，使得探测仪能够获得包含更加丰富和精细的大气状态参数信息，大量探测器每天产生庞大数据给地面数据的传输和存储带来了巨大困难，目前应用较为广泛的iasi数据，在3.62-15.5um的波段范围内具有8461个探测通道，每个通道空间分辨率为60
×
1530
×
16bits，这样一幅iasi超光谱图像的数据量大约为1.45gbytes，如此巨大的数据量给超光谱图像的存储传输和应用处理带来巨大困难，另外，超光谱图像主要用于反演大气温湿度和微量气体等，这些地球物理参数的物理反演涉及到求解辐射传输积分方程，这类方程通常存在非线性问题，微小的数据扰动会带来巨大影响，因此，十分必要研究一种针对于超光谱大气红外遥感图像的无损压缩方法。
3.超光谱图像不仅具有类似于普通图像的空间相关性，而且具有极强的光谱相关性，其压缩方法主要是去除空间上和光谱上的冗余性达到减小熵的目的，结合熵编码减小数据量。就国内外公开的文献来看，现有的压缩方法跟三维地物高光谱图像及其他传统图像的研究类似，主要可以分为三类，基于预测技术、基于变换技术以及基于矢量量化技术的压缩方法。基于预测和矢量量化的技术都是利用超光谱图像的局部相关性例如c-dpcm-rlso方法，基于变换的方法更加注重超光谱图像整体上的数据去相关，例如目前被认为去除光谱间线性相关性最优的pca(主成分分析)方法。而基于小波的方法，例如3d-spiht方法，都是先对空间二维和光谱一维上分别进行小波变换，在对小波系数和误差量化编码。而目前效果较好的多级聚类rklt m-calic方法虽然采用了可逆kl变换和谱间线性预测，但也没有很好的利用已预测的光谱信息和空间数据分布对光谱维度的影响。
4.目前已有超光谱数据压缩方法没有综合考虑超光谱图像在光谱维度上的局部相关性和整体数据分布以及选取多个预测数据时的互相影响。


技术实现要素：

5.本发明的目的是：针对现有技术中已有超光谱数据压缩方法没有综合考虑超光谱图像在光谱维度上的局部相关性和整体数据分布以及选取多个预测数据时的互相影响进而导致预测准确性低的问题，提出基于正交匹配追踪的超光谱图像预测压缩方法。
6.本发明为了解决上述技术问题采取的技术方案是：
7.基于正交匹配追踪的超光谱图像预测压缩方法，包括以下步骤：
8.步骤一：设定参考通道数目上限s和空间聚类分组数k；
9.步骤二：获取待压缩的超光谱图像，将m
×n×
k维超光谱图像的前两个维度合并形成(m
×
n)
×
k维度的矩阵x，并记录原始数据维度m,n,k；
10.步骤三：使用聚类算法将矩阵x按行进行中心化后，得到均值mu，然后将矩阵x按行分为k组，得到分组索引idx，并根据分组索引idx得到子矩阵x
(1)
,x
(2)
...x
(k)
，最后记录均值mu和分组索引idx；
11.步骤四：对于子矩阵x
(j)
，1≤j≤k，将子矩阵x
(j)
第一列作为已知数据矩阵d的第一列且记录在残差矩阵e的第一列e(:,1)＝x
(j)
(:,1)，然后将系数矩阵c的第一行填充0，c(1,:)＝0，将子矩阵x
(j)
的第2列至第s列，执行步骤五，将子矩阵于x
(j)
第s 1列至第k列，执行步骤六至七；
12.步骤五：利用最小二乘法对已知数据矩阵d求解，得到d的伪逆，然后将x
(j)
的第i列记为y，求解方程解为c为y预测系数，然后根据c和已知数据矩阵d得到y预测残差e，表示为e＝y-[dc]，其中[]表示向下取整，最后将c用0补齐后并记录在系数矩阵c的第i行，将e记录在残差矩阵e的第i列，并根据y和e更新已知数据矩阵d，表示为
[0013]
步骤六：将x
(j)
的第i列记为y，并根据已知数据矩阵d利用稀疏最小二乘方法得到c；
[0014]
步骤七：选择没有残差数值溢出的迭代次数t，并将该次迭代中的c记录在c的第i行，计算取整残差e＝y-[dc]，将e记录在残差矩阵e的第i列，并根据y和e更新已知数据矩阵d，d表示为
[0015]
步骤八：压缩残差矩阵e得到残差压缩码流，获取系数矩阵c中非零元素位置，然后使用区间编码器对非零元素位置、均值mu、分组索引idx和m,n,k进行编码，得到区间编码器压缩码流，
[0016]
根据非零元素位置索引系数矩阵c得到非零元素的值，
[0017]
最后将非零元素的值、残差压缩码流和区间编码器压缩码流发送给解码端；
[0018]
步骤九：解码端将非零元素的值、残差压缩码流和区间编码器压缩码流解码得到残差矩阵e、系数矩阵c、分组索引idx、均值mu、原始数据维度标识m、n、k；
[0019]
步骤十：根据分组索引idx和原始数据维度标识m、n、k得到每一个子矩阵x
(j)
的维度，对于子矩阵x
(j)
的第一列，将残差矩阵e的第一列作为已知数据矩阵d的第一列，d表示为根据残差矩阵e的第一列得到x
(j)
第一列的原始数据，表示为x
(j)
(:,1)＝e(:,1)；
[0020]
步骤十一：对于子矩阵x
(j)
的第i列，使用x
(j)
(:,i)＝[dc] e(:,i)无损恢复，[]表示向下取整，2≤i≤k；
[0021]
步骤十二：使用更新d；
[0022]
步骤十三：利用x
(j)
、分组索引idx和原始数据维度m,n,k进行重构，得到矩阵x，然后使用均值mu将矩阵x逆中心化，得到原始超光谱图像。
[0023]
进一步的，所述步骤六的具体步骤为：
[0024]
步骤六一：将x
(j)
的第i列记为y，记迭代次数t＝1，初始化残差r
t
＝y，初始化参考矩阵ψ为空矩阵ψ＝null；
[0025]
步骤六二：计算r
t
与已知数据矩阵d的内积向量q，表示为q＝《r
t
,d》；
[0026]
步骤六三：记q中绝对值最大的元素索引为p，将已知数据矩阵d的第p列放入参考矩阵ψ的第t列，表示为ψ(:,t)＝d(:,p)；
[0027]
步骤六四：判断ψ
t
ψ是否可逆，若不可逆，则执行步骤六五，若可逆，则执行步骤六七；
[0028]
步骤六五：记q中绝对值第二大的元素索引为p，将已知数据矩阵d的第p列放入参考矩阵ψ的第t列，表示为ψ(:,t)＝d(:,p)，然后判断ψ
t
ψ是否可逆，若不可逆，则执行步骤六六，若可逆，则执行步骤六七；
[0029]
步骤六六：利用最小二乘法对ψ求解，得到ψ的伪逆表示为
[0030]
解为更新残差r
t
＝y-ψc，t＝t 1；
[0031]
步骤六七：利用最小二乘法求解方程解为c＝(ψ
t
ψ)-1
ψ
t
y，更新残差r
t
＝y-ψc，记迭代次数t＝t 1；
[0032]
步骤六八：判断迭代次数t是否小于s，若t小于s，则返回六二，若t不小于s，则执行步骤七。
[0033]
进一步的，所述步骤五中最小二乘法为svd方法或正规方程直接求逆法。
[0034]
进一步的，所述步骤六六中最小二乘法为svd方法，所述步骤六七中最小二乘法为正规方程直接求逆法。
[0035]
进一步的，所述参考通道数目上限s＝9，空间聚类分组数k＝8。
[0036]
进一步的，所述步骤三中聚类算法为k均值聚类算法。
[0037]
进一步的，所述残差压缩码流包括单通道残差压缩码流和多通道残差压缩码流。
[0038]
进一步的，所述单通道残差压缩码流通过压缩残差矩阵e的第一个通道得到，所述多通道残差压缩码流通过压缩残差矩阵e的剩余通道得到。
[0039]
进一步的，所述压缩残差矩阵e的第一个通道使用jpeg2000。
[0040]
进一步的，所述压缩残差矩阵e的剩余通道使用webp。
[0041]
本发明的有益效果是：
[0042]
1)通过空间聚类使得相似的光谱束一起预测，每组每个光谱通道共享一组线性预测系数。
[0043]
2)在已知通道数目多于最多预测系数数目时，使用正交匹配追踪的方法代替传统线性预测求解最小二乘的方法，传统最小二乘方法的参考通道一般选取当前通道的前s个
通道，其思想为临近光谱通道相关性高，它们之中包含的当前通道的信息多，但是，正因为临近光谱通道相关性高，当s大于1时，这些参考通道也包含了数据冗余，这会导致方程求解出现奇异值的同时也没有更好的利用上所有的已知数据信息。而正交匹配追踪方法利用迭代的思想逐一选择参考通道，每一次选择的通道都与当前预测残差相关性最强，产生的预测残差由于最小二乘的约束都与参考通道正交，这样可以保证每一次迭代中都是选择减小残差最有利的通道且考虑了之前选择的参考通道对本次选择的影响。
[0044]
3)考虑到应用svd方法求取矩形矩阵伪逆矩阵耗时较长，而计算机求满秩实对称矩阵的逆较为便捷，所以在运行次数最多的步骤九至十二，即根据参考矩阵求取当前通道的最小二乘时，应倾向于更多的直接使用正规矩阵求解而不是svd方法以保证算法运行效率，但是如果矩阵不可逆，只能使用svd方法保证算法能够顺利运行。
[0045]
4)本技术通过omp算法自动地在已知数据中根据每一次输入的待预测数据选择指定书目的最优预测参考通道，如果具有先前已知的数据(如过往已经传输过的超光谱图像)，可以利用这些数据初始化已知数据集合，omp算法会自动进行选择，也就是说增加更多的已知数据不仅不会带来更大的预测系数体积，并且提高了本技术的预测准确性。
附图说明
[0046]
图1为十五次实验中压缩码流各组分占用空间对比图；
[0047]
图2为编码端的流程框图；
[0048]
图3为解码端流程框图。
具体实施方式
[0049]
需要特别说明的是，在不冲突的情况下，本技术公开的各个实施方式之间可以相互组合。
[0050]
具体实施方式一：参照图1具体说明本实施方式，本实施方式所述的基于正交匹配追踪的超光谱图像预测压缩方法，包括以下步骤：
[0051]
步骤一：设定参考通道数目上限s和空间聚类分组数k；
[0052]
步骤二：获取待压缩的超光谱图像，将m
×n×
k维超光谱图像的前两个维度合并形成(m
×
n)
×
k维度的矩阵x，并记录原始数据维度m,n,k；
[0053]
步骤三：使用聚类算法将矩阵x按行进行中心化后，得到均值mu，然后将矩阵x按行分为k组，得到分组索引idx，并根据分组索引idx得到子矩阵x
(1)
,x
(2)
...x
(k)
，最后记录均值mu和分组索引idx；
[0054]
步骤四：对于子矩阵x
(j)
，1≤j≤k，将子矩阵x
(j)
第一列作为已知数据矩阵d的第一列且记录在残差矩阵e的第一列e(:,1)＝x
(j)
(:,1)，然后将系数矩阵c的第一行填充0，c(1,:)＝0，将子矩阵x
(j)
的第2列至第s列，执行步骤五，将子矩阵于x
(j)
第s 1列至第k列，执行步骤六至七；
[0055]
步骤五：利用最小二乘法对已知数据矩阵d求解，得到d的伪逆，然后将x
(j)
的第i列
记为y，求解方程解为c为y预测系数，然后根据c和已知数据矩阵d得到y预测残差e，表示为e＝y-[dc]，其中[]表示向下取整，最后将c用0补齐后并记录在系数矩阵c的第i行，将e记录在残差矩阵e的第i列，并根据y和e更新已知数据矩阵d，表示为
[0056]
步骤六：将x
(j)
的第i列记为y，并根据已知数据矩阵d利用稀疏最小二乘方法得到c；
[0057]
步骤七：选择没有残差数值溢出的迭代次数t，并将该次迭代中的c记录在c的第i行，计算取整残差e＝y-[dc]，将e记录在残差矩阵e的第i列，并根据y和e更新已知数据矩阵d，d表示为
[0058]
步骤八：压缩残差矩阵e得到残差压缩码流，获取系数矩阵c中非零元素位置，然后使用区间编码器对非零元素位置、均值mu、分组索引idx和m,n,k进行编码，得到区间编码器压缩码流，
[0059]
根据非零元素位置索引系数矩阵c得到非零元素的值，
[0060]
最后将非零元素的值、残差压缩码流和区间编码器压缩码流发送给解码端；
[0061]
步骤九：解码端将非零元素的值、残差压缩码流和区间编码器压缩码流解码得到残差矩阵e、系数矩阵c、分组索引idx、均值mu、原始数据维度标识m、n、k；
[0062]
步骤十：根据分组索引idx和原始数据维度标识m、n、k得到每一个子矩阵x
(j)
的维度，对于子矩阵x
(j)
的第一列，将残差矩阵e的第一列作为已知数据矩阵d的第一列，d表示为根据残差矩阵e的第一列得到x
(j)
第一列的原始数据，表示为x
(j)
(:,1)＝e(:,1)；
[0063]
步骤十一：对于子矩阵x
(j)
的第i列，使用x
(j)
(:,i)＝[dc] e(:,i)无损恢复，[]表示向下取整，2≤i≤k；
[0064]
步骤十二：使用更新d；
[0065]
步骤十三：利用x
(j)
、分组索引idx和原始数据维度m,n,k进行重构，得到矩阵x，然后使用均值mu将矩阵x逆中心化，得到原始超光谱图像。
[0066]
本技术实现步骤如下：
[0067]
步骤一：设定参考通道数目上限，空间聚类分组数，输入待压缩超光谱图像。
[0068]
步骤二：对超光谱图像重组为由光谱束按行构成的矩阵，不同行代表不同光谱束，其中光谱束指的是在同一个空间点的各个光谱通道下的像素构成的向量。
[0069]
步骤三：将上述矩阵按行分组划分为多个子矩阵，对每个子矩阵按列中心化，即使得每一个子矩阵的每一列均值为0。这里保证了以下2点：
[0070]
1)线性相关性与内积等价，通常用pearson线性相关系数表示变量间线性关系，其
计算公式如下：
[0071][0072]
其中fi,fj分别为第i,j通道图像拉伸成的向量，m为向量中元素个数，μi,μj为其均值，计算公式如下：
[0073][0074][0075]
当满足以下条件时
[0076]
μi＝μj＝0
[0077][0078]
考虑到fi,fj的内积为如果只考虑对于相同的i的r
ij
大小排序时，内积的计算与线性相关性的计算等价，即说明在步骤八中选择的参考通道就是线性相关性最高的通道。
[0079]
2)由于各列中心化，即在每一光谱束分组下每一个光谱通道的均值都是0，则该分组下不同光谱通道的线性组合的均值也是0，在计算线性预测系数时无需考虑常数项。
[0080]
步骤四：对每一个子矩阵，将其第一列作为已知数据并记录于残差，对于第2至s列，执行步骤五(使用全部已知数据进行线性预测)，对于第s 1至k列，执行步骤六至十四(优选s列已知数据进行线性预测)。
[0081]
步骤五：使用伪逆法利用已知数据对子矩阵的第i列线性表示，记录残差与系数并在预测完毕后将子矩阵的第i列和残差加入已知数据，使用svd分解给出moore-penrose伪逆矩阵。
[0082]
步骤六：将子矩阵的第i列作为待预测数据，初始化残差，迭代次数，参考矩阵。
[0083]
步骤七：计算残差与已知数据的内积。
[0084]
步骤八：记录内积最大的位置索引，将其对应的列加入参考矩阵，即认为与残差内积大的已知数据向量能更好的预测残差。
[0085]
步骤九：判断是否可以使用正规方程直接求解，如果不可以执行步骤十，可以则执行步骤十二。
[0086]
步骤十：将参考矩阵新加入的列替换为步骤八中内积第二大的列，判断是否可以使用正规方程直接求解，如果不可以执行步骤十一，可以则执行步骤十二。
[0087]
步骤十一：使用伪逆求解待预测数据在参考矩阵的线性表示系数，更新残差。
[0088]
步骤十二：直接使用正规方程求解待预测数据在参考矩阵的线性表示系数，更新残差。
[0089]
步骤十三：判断迭代次数是否大于阈值，若小于阈值则返回步骤七，否则执行步骤十四。
[0090]
步骤十四：寻找满足残差各元素都没有数值溢出的迭代，保存其系数，计算取整残差并记录，并在预测完毕后将子矩阵的第i列和残差加入已知数据。
[0091]
步骤十五：使用jpeg2000压缩预测残差的第一个通道，使用webp压缩残差剩余通道，使用区间编码器编码预测系数中非零元素位置，由于非零元素的值使用32位浮点数记录且数目较少，不进行压缩直接传输。
[0092]
实验数据来自于欧洲气象卫星开发组织metop系列卫星上搭载的iasi探测仪的l1c级数据(60行
×
1530列
×
8461通道
×
16bits)，实验一共对15个不同时间不同轨道的数据进行了无损压缩编码，数据来源于欧洲气象卫星应用组织的数据中心(http://catalogue.ceda.ac.uk/)。
[0093]
表1实验所用15个iasi完整轨道的超光谱探测数据具体信息及压缩比
[0094][0095]
其中压缩比计算方式为(原始数据大小)/(压缩后数据大小)，压缩比越大说明压缩效果越好。
[0096]
在参数参考通道数目上限s＝9，空间聚类分组数k＝8时，本发明平均压缩比为2.5386，与目前压缩效果最好的多级聚类rklt m-calic方法2.54的压缩比相当，且本方法支持光谱通道渐进传输，一方面可以使得接收端无需等待发送端将全部光谱通道压缩完毕
在传输数据，更加适应了目前超光谱图像光谱通道数目快速增长的发展趋势，另一方面当接收端也可以在数据传输一开始就进行解压操作。
[0097]
需要注意的是，具体实施方式仅仅是对本发明技术方案的解释和说明，不能以此限定权利保护范围。凡根据本发明权利要求书和说明书所做的仅仅是局部改变的，仍应落入本发明的保护范围内。