CN107770526B - 一种基于量化ica的超光谱大气红外遥感图像无损压缩方法 - Google Patents

Abstract

本发明提供一种基于ICA的超光谱大气红外遥感图像无损压缩方法，本发明为了解决现有的遥感图像压缩方法不适用于超光谱大气红外遥感图像处理过程的问题，而提出一种基于ICA的超光谱大气红外遥感图像无损压缩方法，包括：将三维超光谱数据转换成二维矩阵，并进行处理。求出ICA变换的分离矩阵，再求得ICA变换系数矩阵；然后对变换系数矩阵和独立成分矩阵进行量化得到矩阵AQ和YQ，再对矩阵AQ和YQ进行反量化和逆ICA得到的结果与原图做差，得到残差矩阵D，对AQ和YQ进行预测得到残差矩阵AQP和YQP。最后对三个残差矩阵D、AQP和YQP进行区间编码得到压缩码流。本方法能够对超光谱大气红外遥感图像进行有效的压缩，取得较高的压缩比。

Description

一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法

技术领域

本发明涉及遥感图像无损压缩领域，具体涉及一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法。

背景技术

随着超光谱大气红外遥感探测技术的发展，对大气的探测越来越精细，探测周期越来越短，从而探测信息的数据量也随之越来越大，无论是在星上还是在星下，对于探测信息的存储和传输是在数据应用过程当中必然面对的问题。因此，为达到快速传输超光谱大气红外遥感图像数据，对其进行压缩以减少数据存储空间是必要的。近年来，对超光谱大气红外遥感图像数据的有效压缩研究大多数依旧沿用已经成熟的对地高光谱图像数据的压缩方法。但是，超光谱大气红外图像有其相应的特点，其光谱分辨率要高于对地高光谱图像一个数量级，但其空间分辨率却比对地高光谱要低，并且其具有明显的分段特性等，而以往的技术未能考虑到超光谱图像自身空谱特点和物理特性。

发明内容

本发明的目的是为了解决现有的遥感图像压缩方法不适用于超光谱大气红外遥感图像处理过程的问题，而提出一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法。

一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法，包括：

步骤一、三维超光谱大气红外遥感图像数据大小为M×N×L，波段通道总数为L个，将其转化为L行M×N列的矩阵X，其中矩阵X的每一行代表一个光谱维的图像数据，是将M×N按行展开的结果；再求数据矩阵X各行向量中所有元素的平均值，并对X进行零均值处理得到

步骤二、对

进行白化处理，得到白化后的结果矩阵Z，矩阵Z的大小为L×(M×N)；

步骤三、采用基于极大似然估计的不动点算法求解分离矩阵U，通过ICA变换去除超光谱图像数据的谱间冗余，得到独立成分矩阵Y以及ICA的变换系数矩阵A；

步骤四、对独立成分矩阵Y以及ICA的变换系数矩阵A进行量化，分别得到YQ矩阵和AQ矩阵，设定s表示量化级别的数目，量化公式为：

其中，q表示量化索引，d为量化输入的实际数据值，d_max和d_min分别是实际输入的数据中的最大值和最小值，(d_max-d_min)/(2^s-1)表示量化的步长大小；round表示取整；

由矩阵A量化得到矩阵AQ的过程为：将矩阵A中的各个元素作为d带入至量化公式中得到量化索引q，再将每一个元素计算得到的量化索引q构成矩阵得到YQ矩阵；由矩阵Y得到矩阵YQ的过程同理。

在对量化后的独立成分矩阵Y和ICA的变换系数矩阵A进行反量化得到矩阵

和

再通过逆ICA变换，即计算

得到恢复后的数据矩阵

再将其与白化后的矩阵Z做差，得到残差矩阵D；反量化为量化公式的逆运算；

步骤五、将量化后得到的YQ矩阵和AQ矩阵按照步骤一中相反的方式还原成三维数据形式；还原后的数据大小为n_x×n_y×n_l，用i,j,k分别代表数据所在的行、列以及光谱位置，则得到还原后的数据用c_i′_,j,k表示，通过以下的公式对其进行预测得到预测结果

其中median函数表示取中值；1＜i≤n_x,1＜j≤n_y,1≤k≤n_l；

然后得到预测误差e_i,j,k：

按照以下三个公式进行计算和更新：

e_i,1,k←e_i,1,k-e_i-1,1,k

e_1,j,k←e_1,j,k-e_1,j-1,k

e_i,j,k←e_i,j,k-e_i,j,k-1

其中，1＜i≤n_x，1＜j≤n_y，1＜k≤n_l；

最后将所得的三维预测误差按照步骤一中的方式重新构造成二维形式，得到YQ矩阵预测后的结果矩阵YQP矩阵，AQ矩阵预测后的结果矩阵AQP矩阵；

步骤六、将待编码数据，包括残差矩阵D、YQP矩阵和AQP矩阵进行区间编码，产生的码流中加入4个辅助信息，分别是ICA的独立成分数num、量化级别的数目s和最小值d_min、以及各个光谱段的均值数据形成最终的压缩码流，完成压缩。

优选地，步骤一中，零均值处理的具体过程为：

令变量

为数据矩阵X各行向量中所有元素的平均值，即：

其中1≤i≤L；均值向量

由M×N个

元素组成，表达式为：

然后对X进行以下零均值处理得到

优选地，步骤二具体为：

步骤二一、令C为零均值数据矩阵

的协方差矩阵，将C进行以下分解

其中，C为L阶矩阵，Λ为对角阵，即Λ＝Diag[λ₁,λ₂,...,λ_L]，其中λ_i为协方差矩阵C的特征根，λ₁≥λ₂≥...≥λ_L≥0，V＝[v₁,v₂,...v_L]，其中v_i为相应的特征向量；

步骤二二、把W＝Λ^-1/2V^T作为变换算子，左乘到零均值数据矩阵

上，得到：

则矩阵Z为对数据阵

白化后的结果，W＝Λ^-1/2V^T称为白化矩阵，并且矩阵Z中各行相互正交，ZZ^T＝I，I为L阶单位矩阵。

优选地，步骤三具体为：步骤三一、确定独立成分数num，并初始化分离矩阵U，U＝[u₁,u₂,...u_i,...,u_num]^T，其中，u_i＝[u₁,u₂,...,u_L](1≤i≤num)，元素u₁,u₂,...,u_L能够使得|u_i|＝1且互相正交；

步骤三二、计算Y＝UZ，(Y＝[y₁,y₂,...,y_num]^T，大小为num×(M×N))，设定N维对角矩阵K＝diag[k₁,k₂,...,k_num]，其中：

其中1≤i≤num

E表示求括号内数据的数学期望；

sech(x)是双曲正割函数，公式为：

e为自然常数；

sign(x)是符号函数，公式为：

步骤三三、用公式

更新分离矩阵；U←Φ(U)表示将计算得到的Φ(U)赋值给U；

步骤三四、根据U←(UU^T)^-1/2U将矩阵U对称正交化；

步骤三五、判步骤三四中分离矩阵U中的u_i和步骤三二中分离矩阵U中的u_i的内积是否等于1，如果等于1，则保留步骤三四中的分离矩阵U，其大小为num×L，然后计算Y＝UZ得到白化后数据矩阵Z的独立成分矩阵Y，大小为num×(M×N)，Y中各个行向量是ICA变换之后的得到的独立成分，独立成分数目为num。通过实验设定不同的独立成分数num，以得到相对较高的压缩率，然后计算ICA的变换系数矩阵A＝U^-1；否则，如果不相等，返回步骤三二，直到满足内积等于1。

本发明的有益效果为：1、具有较高的压缩比；2、可有效地在高阶统计量上取出相关性，更适合超光谱大气红外遥感图像的处理过程；3、平均压缩率可达3.37，高于其他压缩算法。

附图说明

图1是本发明所述的基于量化ICA的超光谱大气红外遥感图像无损压缩方法原理示意图；

图2是本发明所采用的ICA算法的变换过程示意图；

图3是本发明在实验验证时所用的AIRS超光谱大气红外遥感图像数据；其中Granule表示“景”，用于研究的标准超光谱探测仪AIRS的数据集包含10景，5个白天，5个夜晚，同时代表了不同的地理区域；一天的数据共有240景，景后面的数字，例如Granule60，其中的60为序号。

具体实施方式

具体实施方式一、结合图1说明本发明的具体实施方式，一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法：

步骤一、三维超光谱大气红外遥感图像数据大小为M×N×L，波段通道总数为L个，将其转化成二维，得到L行M×N列的矩阵X＝[x₁,x₂,...,x_L]^T，其中，x_i(1≤i≤L)为M×N维的行向量。

再对其进行零均值处理，具体方式如下：

令变量

为数据矩阵X各行向量中所有元素的平均值，即：

均值向量

由M×N个

元素组成，可写成如下形式：

然后对X进行零均值处理后得到

步骤二、对

进行白化处理，得到白化后的结果矩阵Z，具体方式如下：

令C为零均值数据矩阵

的协方差矩阵，可将C进行如下分解

其中，C为L阶矩阵，Λ为对角阵，即Λ＝Diag[λ₁,λ₂,...,λ_L]，其中λ_i(1≤i≤L)为协方差矩阵C的特征根，λ₁≥λ₂≥...≥λ_L≥0，V＝[v₁,v₂,…v_L]，其中v_i(1≤i≤L)为相应的特征向量。

把W＝Λ^-1/2V^T作为变换算子，左乘到零均值数据矩阵

上，可得到：

则矩阵Z为对数据阵

白化后的结果，W＝Λ^-1/2V^T称为白化矩阵，并且矩阵Z中各行相互正交，ZZ^T＝I，I为L阶单位矩阵。

步骤三.采用基于极大似然估计的不动点算法求解ICA的分离矩阵U，通过ICA变换去除超光谱图像数据的谱间冗余，得到独立成分矩阵Y以及ICA的变换系数矩阵A；

步骤三一.首先，确定独立成分数num，然后，初始化分离矩阵U，U＝[u₁,u₂,…u_i,...,u_num]^T，其中，u_i＝[u₁,u₂,...,u_L](1≤i≤num)，|u_i|＝1且互相正交。

步骤三二.计算Y＝UZ，(Y＝[y₁,y₂,...,y_num]^T，大小为num×(M×N))，设定num维对角矩阵K＝diag[k₁,k₂,...,k_num]，其中(E指的是求括号内数据的数学期望)：

sech(x)是双曲正割函数，公式为：

e为自然常数；

sign(x)是符号函数，公式为：

步骤三三.用公式

更新分离矩阵。

步骤三四.根据U←(UU^T)^-1/2U将矩阵U对称正交化。U←Φ(U)表示将计算得到的Φ(U)赋值给U。以下出现的箭头也是相同的含义，表示将箭头起始端的变量计算完后赋给箭头末端的变量，以此来完成对箭头末端变量的更新。

步骤三五.判步骤三四中分离矩阵U中的u_i和步骤三二中分离矩阵U中的u_i的内积是否等于1，如果等于1，则可以得到分离矩阵U，大小为num×L，然后计算Y＝UZ得到白化后数据矩阵Z的独立成分矩阵Y，大小为num×(M×N)，Y中各个行向量是ICA变换之后的得到的独立成分，其独立成分数目为num。为了得到相对较高的压缩率，通过实验分析，我们采用的独立成分的数目num＝70，然后计算ICA的变换系数矩阵A＝U^-1。否则，如果不相等，返回步骤三二，直到满足内积等于1为止。

步骤四、对独立成分矩阵Y以及ICA的变换系数矩阵A进行量化，分别得到YQ矩阵和AQ矩阵，设定s表示量化级别的数目，均匀量化方式如下式：

其中，q表示量化索引，d为量化输入的实际数据值，d_max和d_min分别是实际输入的数据中的最大值和最小值，(d_max-d_min)/(2^s-1)表示量化的步长大小。通过实验分析，对ICA的变换系数矩阵A进行量化处理时的量化级别为11，对独立成分矩阵Y进行量化处理时选取的量化级别为15时，结果最优。round表示取整；

然后，对量化后的独立成分矩阵Y和ICA的变换系数矩阵A进行反量化得到矩阵

和

再通过逆ICA(IICA)变换，即计算

得到恢复后的数据矩阵

再将其与白化后的矩阵Z做差，得到残差矩阵D。

步骤五、将量化后得到的YQ矩阵和AQ矩阵按照步骤一中相反的方式还原成三维数据形式。还原后的数据大小为n_x×n_y×n_l，用i,j,k分别代表数据所在的行、列以及光谱位置，则得到还原后的数据用c′_i,j,k表示，通过下面的公式对其进行预测得到预测结果

(median函数表示取中值)：

然后可以得到预测误差e_i,j,k：

按照如下三个公式进行计算和更新：

e_i,1,k←e_i,1,k-e_i-1,1,k(1＜i≤n_x)

e_1,j,k←e_1,j,k-e_1,j-1,k(1＜j≤n_y)

e_i,j,k←e_i,j,k-e_i,j,k-1(1＜k≤n_l)

最后将所得的三维预测误差按照步骤一中的方式重新构造成二维形式，得到YQ矩阵预测后的结果矩阵YQP矩阵，AQ矩阵预测后的结果矩阵AQP矩阵。

步骤六、将待编码数据，包括残差矩阵D、YQP矩阵和AQP矩阵进行区间编码，产生的码流中加入4个辅助信息，分别是ICA的独立成分数num、量化级别的数目s和最小值d_min、以及各个光谱段的均值数据

形成最终的压缩码流，完成压缩。辅助信息可以具体加在帧的头部中。

验证实验采用的数据是美国国家航天太空总署的标准超光谱红外大气数据，共分为十个景，其中五个为白天时间，五个为夜晚时间，代表了地球气象的典型地理区域。每一个数据块被保存为一个二进制文件，共2107个波段通道，135个辐射扫描面，每个波段通道的扫描线上共90个像素，每个数据块的存储大小为41249250bytes。图3是在波数800.01cm^-1下的十个AIRS图像数据块。

表1最后一行是本发明方法(PICA)，即量化ICA的超光谱大气红外遥感图像无损压缩方法，对图3中十景AIRS数据块进行无损压缩的结果，并与其他典型的压缩方法的对比结果。

表1各种算法的压缩率对比结果

从实验结果可以看出，本发明方法压缩这十景图像均比传统方法具有更高的压缩率，平均压缩率可达到3.37。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (4)

1.一种基于量化ICA的超光谱大气红外遥感图像无损压缩方法，其特征在于，

步骤一、三维超光谱大气红外遥感图像数据大小为M×N×L，波段通道总数为L个，将其转化为L行M×N列的矩阵X，其中矩阵X的每一行代表一个光谱维的图像数据；再求数据矩阵X各行向量中所有元素的平均值，并对X进行零均值处理得到

步骤二、对

进行白化处理，得到白化后的结果矩阵Z，矩阵Z的大小为L×(M×N)；

步骤三、求解分离矩阵U，并进一步得到独立成分矩阵Y以及ICA的变换系数矩阵A；

步骤四、对独立成分矩阵Y以及ICA的变换系数矩阵A进行量化，分别得到YQ矩阵和AQ矩阵，设定s表示量化级别的数目，量化公式为：

其中，q表示量化索引，d为量化输入的实际数据值，d_max和d_min分别是实际输入的数据中的最大值和最小值，(d_max-d_min)/(2^s-1)表示量化的步长大小；round表示取整；

由矩阵A量化得到AQ矩阵的过程为：将矩阵A中的各个元素作为d带入至量化公式中得到量化索引q，再将每一个元素计算得到的量化索引q构成矩阵得到YQ矩阵；

在对量化后的独立成分矩阵Y和ICA的变换系数矩阵A进行反量化得到矩阵

和

再通过逆ICA变换，即计算

得到恢复后的数据矩阵

再将其与白化后的矩阵Z做差，得到残差矩阵D；反量化为量化公式的逆运算；

步骤五、将量化后得到的YQ矩阵和AQ矩阵按照步骤一中相反的方式还原成三维数据形式；还原后的数据大小为n_x×n_y×n_l，用i,j,k分别代表数据所在的行、列以及光谱位置，则得到还原后的数据用c′_i,j,k表示，通过以下的公式对其进行预测得到预测结果

其中median函数表示取中值；1＜i≤n_x,1＜j≤n_y,1≤k≤n_l；

然后得到预测误差e_i,j,k：

按照以下三个公式进行计算和更新：

e_i,1,k←e_i,1,k-e_i-1,1,k

e_1,j,k←e_1,j,k-e_1,j-1,k

e_i,j,k←e_i,j,k-e_i,j,k-1

其中，1＜i≤n_x，1＜j≤n_y，1＜k≤n_l；箭头表示将箭头右侧式子的运算结果赋值给箭头左侧的值；

最后将所得的三维预测误差e_i,j,k所构成的矩阵按照步骤一中的方式重新构造成二维形式，得到YQ矩阵预测后的结果矩阵YQP矩阵，AQ矩阵预测后的结果矩阵AQP矩阵；

步骤六、将待编码数据进行区间编码，所述待编码数据包括残差矩阵D、YQP矩阵和AQP矩阵，产生的码流中加入4个辅助信息，分别是ICA的独立成分数num、量化步长q和最小值d_min、以及各个光谱段的均值数据形成最终的压缩码流，完成压缩。

2.根据权利要求1所述的基于量化ICA的超光谱大气红外遥感图像无损压缩方法，其特征在于，步骤一中，零均值处理的具体过程为：

令变量

为数据矩阵X各行向量中所有元素的平均值，即：

其中1≤i≤L；x_ij为矩阵第i行第j列的元素；均值向量

由M×N个

元素组成；

然后对X进行以下零均值处理得到

3.根据权利要求1所述的基于量化ICA的超光谱大气红外遥感图像无损压缩方法，其特征在于，步骤二具体为：

步骤二一、令C为零均值数据矩阵

的协方差矩阵，将C进行以下分解

其中，C为L阶矩阵，Λ为对角阵，即Λ＝Diag[λ₁,λ₂,...,λ_L]，其中λ_i为协方差矩阵C的特征根，λ₁≥λ₂≥...≥λ_L≥0，V＝[v₁,v₂,...v_L]，其中v_i为相应的特征向量；

步骤二二、把W＝Λ^-1/2V^T作为变换算子，左乘到零均值数据矩阵

上，得到：

则矩阵Z为对数据阵

白化后的结果，W＝Λ^-1/2V^T称为白化矩阵，并且矩阵Z中各行相互正交，ZZ^T＝I，I为L阶单位矩阵。

4.根据权利要求1所述的基于量化ICA的超光谱大气红外遥感图像无损压缩方法，其特征在于，步骤三具体为：

步骤三一、确定独立成分数num，并初始化分离矩阵U，U＝[u₁,u₂,...u_i,...,u_num]^T，其中，u_i＝[u₁,u₂,...,u_L]，1≤i≤num，元素u₁,u₂,...,u_L能够使得|u_i|＝1且互相正交；

步骤三二、计算Y＝UZ，Y＝[y₁,y₂,...,y_num]^T，大小为num×(M×N)，设定对角矩阵K＝diag[k₁,k₂,...,k_num]，其中：

sech(x)是双曲正割函数，公式为：

e为自然常数；

sign(x)是符号函数，公式为：

E指的是求括号内数据的数学期望；

步骤三三、用公式

更新分离矩阵；U←Φ(U)表示将计算得到的Φ(U)赋值给U；

步骤三四、根据U←(UUT)^-1/2U将矩阵U对称正交化；

步骤三五、判步骤三四中分离矩阵U中的u_i和步骤三一中分离矩阵U中的u_i的内积是否等于1，如果等于1，则保留步骤三四中的分离矩阵U，其大小为num×L，然后计算Y＝UZ得到白化后数据矩阵Z的独立成分矩阵Y，大小为num×(M×N)，Y中各个行向量是ICA变换之后的得到的独立成分，独立成分数目为num；通过实验设定不同的独立成分数num，以得到相对较高的压缩率，然后计算ICA的变换系数矩阵A＝U^-1；否则，如果不相等，返回步骤三二，直到满足内积等于1。

Images