CN108827470A - 基于全约束最小二乘法的自适应光谱解混合的光学实现方法及系统 - Google Patents

Abstract

为了解决传统高光谱解混合算法处理数据量庞大的技术问题，本发明提供了一种基于全约束最小二乘法的自适应光谱解混合的光学实现方法及系统，利用计算光谱成像技术中的压缩编码原理，将光谱解混算法—全约束最小二乘法，与光谱分光&合光系统相结合，以空间目标场景为分析处理对象，能够直接输出解混合结果，不需要进行光谱数据立方体的重构过程，从而避免了大数据量问题，大幅减轻了系统的数据存储、传输、处理的负担，且为光谱解混合处理的实时应用提供有效的技术支撑。

Description

基于全约束最小二乘法的自适应光谱解混合的光学实现方法 及系统

技术领域

本发明属于光学技术领域，涉及一种高光谱解混合的光学实现方法及系统，尤其涉及一种基于光谱成像仪直接进行光谱解混合的光学实现方法及系统。

背景技术

高光谱成像技术是将成像技术和光谱技术相结合的多维信息获取技术，它的出现续写了光学影像从黑白全色影像经由多光谱到高光谱的影像信息链。高光谱成像技术能够同时探测目标的二维几何空间与一维光谱信息，具有“图谱合一”的特点；且实现了光谱分辨率的突破性提高，使得本来在宽波段多光谱遥感中不可探测的地物在高光谱遥感中能够被探测。

然而，在高光谱遥感中，由于受遥感器的空间分辨率限制以及自然界地物的复杂性、多样性，普遍存在混合像元。混合像元是指每个像元对应区域内往往包含多种不同的地物，它们有着不同的光谱响应特征，导致每个混合像元记录的信息是对应区域内全部特征地物光谱信息的综合叠加。因此，高光谱图像分析往往需要在亚像元级别进行，这就需要首先进行光谱解混合，即根据遥感图像提供的信息判断出每个混合像元是由哪些纯像元、以多少比例混合的。

而高光谱成像获取目标丰富的空间、辐射和光谱信息，提高精细信息表达能力的同时，也带来了新的挑战。高光谱成像仪输出并重构后的数据立方体具有两维空间坐标和一维光谱坐标，包含体量巨大的数据信息。以我国OMIS-1成像光谱仪为例，行像元数为512，总波段数为128，扫描率为20线/s，则1min内将产生7.86×10⁷个信息单元。如此庞大的数据量不仅给数据处理、传输、存储系统带来了极大的压力，同时使得高光谱成像仪无法实现实时或接近实时的应用。

近些年，国内外学者们针对高光谱解混合提出了很多经典算法，这些算法所采用的光谱混合模型分为两种：线性光谱混合模型和非线性光谱混合模型。而无论是哪种模型，目前实现方式均是以高光谱成像仪输出并重构的三维数据立方体为研究对象，先后经过端元提取、丰度反演两步骤，提取每个空间像元包含的端元和对应的端元丰度系数，这就必然要面对高光谱成像带来的大数据量问题。

发明内容

为了解决传统高光谱解混合算法处理数据量庞大的技术问题，本发明提供了一种基于全约束最小二乘法的自适应光谱解混合的光学实现方法及系统，利用计算光谱成像技术中的压缩编码原理，将光谱解混算法—全约束最小二乘法，与光谱分光&合光系统相结合，以空间目标场景为分析处理对象，能够直接输出解混合结果。

本发明的技术解决方案如下：

基于全约束最小二乘法的自适应光谱解混合的光学实现方法，其特征在于，每一行空间像元的光谱解混合包括以下步骤：

步骤1：将一行空间像元的辐射/反射光，分光色散成二维光谱；

步骤2：利用已知参考光谱库中的端元构建端元矩阵E，并根据将E扩展为E′，初始化R＝{1,2,…,m}，A_i＝0，s＝0，r＝1；P和R分别是设定的两个集合；A_i和s是1×b大小的一维矩阵；r为用于计数的中间变量；b为光谱波段数；m为端元数；

步骤3：提取E′^T的第r行系数矩阵，标记为R，并用R的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的；

步骤4：根据步骤3中给出的二维控制信号，利用空间光调制器对步骤1分光后的二维光谱进行光谱调制；

步骤5：将调制后光谱经反向合光系统聚合后投射到探测器上，探测器获取聚合后光谱；

步骤6：将探测器获取的数据与步骤3提取的一维系数矩阵R的第b+1个系数值相加，得到1×m矩阵E′^Tx_i′中第r个数值；

步骤7：如果r＝m，则已获取1×m矩阵E′^Tx_i′，转入步骤8；否则，更新r＝r+1，转入步骤3；

步骤8：根据公式w＝E′^T(x_i′-E′A_i)，计算w；

步骤9：如果且max(w_k)＞0(k∈R)，转入步骤10；否则，转入步骤20；

步骤10：求出t＝argmax(w_k)(k∈R)，更新R＝R-{t}和P＝P∪{t}；

步骤11：根据更新的集合P，更新E′^P，并利用公式E^P#＝[(E′^P)^TE′^P]^-1(E′^P)^T自适应地计算E^P#,pr＝1；pr为用于计数的中间变量；

步骤12：提取E^P#的第pr行系数矩阵，标记为PR，并用PR的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的；

步骤13：根据步骤12中给出的二维控制信号，利用空间光调制器对步骤1分光后的二维光谱进行光谱调制；

步骤14：将调制后光谱经反向合光系统聚合后投射到探测器上，探测器获取聚合后光谱；

步骤15：将探测器获取的数据与步骤12提取的一维系数矩阵PR的第b+1个系数值相加，得到1×m矩阵E^P#x_i′中第pr个数值；

步骤16：如果pr＝m，则已获取1×m矩阵s^P＝E^P#x_i′，并更新s^R＝0，转步骤17；否则，更新pr＝pr+1，转步骤12；

步骤17：如果min(s^P)≤0，转入步骤18；否则，转入步骤19；

步骤18：根据公式A_i＝A_i+β(s-A_i)，计算β和A_i，并把所有满足A_ik＝0的指标k从P移动到R，转入步骤11；

步骤19：更新A_i＝s，w＝E′^T(x_i′-E′A_i)，转入步骤9；

步骤20：当前行的光谱解混合完成，此时A_i为当前行的丰度系数矩阵，进入下一行的光谱解混合过程。

本发明同时还提供了一种基于上述光学实现方法的自适应光谱解混合系统，其特殊之处在于：包括CPU，以及依次设置的分光元件、二维空间光调制器、反向合光系统和探测器；

所述分光元件用于将一行空间像元的辐射/反射光色散成二维光谱；二维空间光调制器用于根据CPU发出的二维控制信号对经分光元件分光后的光谱进行光谱调制；所述反向合光系统用于对调制后光谱进行聚合；所述探测器用于获取聚合后光谱数据，并发送给所述CPU，该CPU与计算机可读存储介质相连，其上存储有计算机程序；所述计算机程序被执行时实现上述方法中的步骤2至步骤3、步骤6至步骤12、步骤15至步骤20。

与现有技术相比，本发明具有如下有益效果：

本发明直接对进入系统的目标场景分光后光谱进行处理，并直接输出光谱解混合结果，不需要进行光谱数据立方体的重构过程，从而避免了大数据量问题，大幅减轻了系统的数据存储、传输、处理的负担，且为光谱解混合处理的实时应用提供有效的技术支撑。

附图说明

图1是现有的计算光谱成像技术流程图；

图2是本发明自适应光谱解混合系统的原理框图。

具体实施方式

光谱解混合过程分为端元提取和丰度反演两个步骤。其中，端元是指仅包含一种地物类型的纯像元，可用参考光谱库中已知样本光谱作为端元光谱。因此，本发明主要研究光谱解混合中丰度反演过程的光学实现方法，而高光谱解混合狭义的理解就是指丰度反演。

具体发明构思如下：

假设一幅具有b个波段的高光谱图像有n个像元，每个像元是由m个端元{e_j}(j＝1,2,…,m)，按照一定比例混合而成，且第j个端元在第i个像元x_i中的比例为a_i,j，则可写为如下公式：

其中，ε_i为像元x_i的误差项，为了简化问题，可忽略；E＝[e₁,e₂,…，e_m]是端元矩阵，A_i＝[a_i,1,a_i,2,…,a_i,m]^T是丰度系数矩阵。

丰度系数有两个约束条件，即“合为1”约束和“非负”约束，分别如公式(2)和(3)所示：

根据满足约束条件的程度，分为四种不同的最小二乘法：

无约束最小二乘法(UCLS)、“合为1”约束最小二乘法(SCLS)、“非负”约束最小二乘法(NCLS)、全约束最小二乘法(FCLS)。

一、对于UCLS法，可直接求解公式(1)的方程组，得到如公式(4)的解。

A_i＝(E^TE)^-1E^Tx_i＝E^#x_i (4)

二、对于SCLS法，可通过公式(5)和(6)对E和x_i进行扩展，用扩展后得到的E′和x_i′带入替代公式(4)中的E和x_i即可。

公式(5)中I_M为全1(m维)行向量。

三、对于NCLS法，算法流程如下：

Step0、初始化集合R＝{1,2,…,m}，初始化1×b矩阵A_i＝0，s＝0，计算w＝E^T(x_i-EA_i)；

Step1、如果且max(w_k)＞0(k∈R)，转step2；否则，转step7；

Step2、t＝argmax(w_k)(k∈R)，更新R＝R-{t}，P＝P∪{t}；等号右边的P、R是更新之前的集合，等号左边的P、R是更新之后的集合；

Step3、s^P＝[(E^P)^TE^P]^-1(E^P)^Tx_i，s^R＝0；其中s^P(s^R)表示向量s中与P(R)中元素对应的分量，E^P表示E中与P中对应的列；

Step4、如果min(s^P)≤0，转step5；否则，转step6；

Step5、A_i＝A_i+β(s-A_i)，把所有满足A_ik＝0的指标k从P移动到R，转step3；等号右边的A_i是更新之前的矩阵，等号左边的A_i是更新之后的矩阵；

Step6、更新A_i＝s，重新计算w＝E^T(x_i-EA_i)，转step1；

Step7、此时A_i就是求得的丰度系数矩阵，算法结束。

上述算法中，P和R分别是设定的两个集合；A_i和s是1×b大小的一维矩阵；A_ik是一维矩阵A_i的第k个系数；w是根据E、x_i、A_i计算得出的中间值，后续步骤中使用；T表示转置，是一种运算符号；m是端元数目；以上参数在step0分别进行初始化。

四、对于FCLS算法，可通过将公式(5)和(6)带入NCLS法中的E和x_i得到。

计算光谱成像技术是在传统的色散型光谱成像技术基础上，通过引入具有特殊数学形式的编码模板，对目标的空间信息和光谱信息进行调制，从而获取若干个压缩编码图像，并通过合理的重构算法，最终形成目标的三维光谱数据立方体，如图1所示。

本发明将计算光谱成像技术中光谱分光、光谱编码、反向合光部分与全约束最小二乘法(FCLS)中对像元光谱的数学变换过程相结合，从而达到光谱解混合算法光学实现的目的。

根据FCLS法，整个解混合过程中，待分解像元光谱有三处参与计算：

1)从x_i扩展为在FCLS法中，对像元光谱x_i的数学变换，均需用对x_i′的数学变换来代替。但是这并不影响其光学实现方式，如变换系数为Q＝[q₁,q₂,…,q_b,q_b+1]，则 即可用变换系数的前b个系数矩阵进行光谱编码，然后将编码后获取的数据与第b+1个系数值相加即可。

2)计算s^P＝[(E′^P)^TE′^P]^-1(E′^P)^Tx_i′。这里将变换矩阵简化表示为E^P#＝[(E′^P)^TE′^P]^-1(E′^P)^T。E^P#的每一行是1×(b+1)大小的系数矩阵。因此，根据每次求得的P，更新E′^P，自适应地计算E^P#；逐行提取E^P#的系数矩阵，每次提取的是1×(b+1)大小的系数矩阵，根据1)中描述，用其中前b个系数矩阵，对分光后光谱进行编码，然后将编码后获取的数据与第b+1个系数值相加，可完成s^P＝[(E′^P)^TE′^P]^-1(E′^P)^Tx_i＇计算的光学实现。

3)计算w＝E′^T(x_i′-E′A_i)。该式可分解为w＝E′Tx_i′-E′TE′A_i，其中E′^TE′A_i可根据已知的E′和已求得的A_i事先计算得出，因此只需光学实现E′^Tx_i′的计算。具体方法如2)中所述，逐行提取E′^T的一维系数矩阵，对分光后光谱进行编码。

基于以上发明构思，结合图2，以下对本发明提供的自适应光谱解混合系统及方法进行详述。

本发明中拟采用空见光调制器实现光谱编码。空间光调制器具有对光的调制特性，是可实现谱带编码功能的典型器件，也是本发明中将高光谱解混合算法光学实现的关键之一。空间光调制器含有许多独立单元，在空间上排列成一维或二维阵列，每个单元都可以独立地接收光学信号或电学信号的控制，并按此信号改变自身的光学性质，从而实现可编程编码模板的功能。图2中虚线框内部分为存储在计算机可读存储介质上的计算机程序，由CPU执行。

本发明采用推扫扫描方式，对目标空间逐行进行光谱解混合处理，如果端元数目是m，每一行空间像元的光谱解混合步骤如下：

步骤1：一行空间像元的辐射/反射光，进入光谱解混合系统后，经过分光元件，色散成二维光谱，投射到二维空间光调制器上。

步骤2：根据已知参考光谱库中的端元，构成端元矩阵E，并利用前述公式(5)将E扩展为E′，并初始化R＝{1,2,…,m}，A_i＝0，s＝0，r＝1。P和R分别是设定的两个集合；A_i和s是1×b大小的一维矩阵；r为用于计数的中间变量；b为光谱波段数；m为端元数。

步骤3：提取E′^T的第r行系数矩阵，标记为R，并用R的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的。

步骤4：根据步骤3中给出的二维控制信号，空间光调制器对步骤1分光后光谱进行光谱调制。

步骤5：调制后光谱经反向合光系统聚合投射到探测器上，探测器获取聚合后光谱。

步骤6：将探测器获取的数据与步骤3提取的一维系数矩阵R的第b+1个系数值相加，得到1×m矩阵E′^Tx_i′中第r个数值。

步骤7：如果r＝m，则已获取1×m矩阵E′^Tx_i′，转步骤8；否则，更新r＝r+1，转步骤3。

步骤8：根据公式w＝E′^T(x_i′-E′A_i)，计算w。

步骤9：如果且max(w_k)＞0(k∈R)，转步骤10；否则转步骤20。

步骤10：求出t＝argmax(w_k)(k∈R)，更新R＝R-{t}和P＝P∪{t}。

步骤11：根据更新的集合P，更新E′^P，并利用公式E^P#＝[(E′^P)^TE′^P]^-1(E′^P)^T，自适应地计算E^P#，pr＝1，pr为用于计数的中间变量。

步骤12：提取E^P#的第pr行系数矩阵，标记为PR，并用PR的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的。

步骤13：根据步骤12中给出的二维控制信号，空间光调制器对步骤1分光后光谱进行光谱调制。

步骤14：调制后光谱经反向合光系统聚合后投射到探测器上。

步骤15：将探测器获取的数据与步骤12提取的一维系数矩阵PR的第b+1个系数值相加，得到1×m矩阵E^P#x_i′中第pr个数值。

步骤16：如果pr＝m，则已获取1×m矩阵s^P＝E^P#x_i′，并更新s^R＝0，转步骤17；否则，pr＝pr+1，转步骤12。

步骤17：如果min(s^P)≤0，转步骤18；否则，转步骤19。

步骤18：根据公式A_i＝A_i+β(s-A_i)，计算β和A_i，并把所有满足A_ik＝0的指标k从P移动到R，转步骤11。

步骤19：更新A_i＝s，并重新获取w＝E′^T(x_i′-E′A_i)，转步骤9。

步骤20：完成了当前行的光谱解混合，此时A_i就是当前行的丰度系数矩阵，然后进入下一行的光谱解混合过程。

Claims (2)

1.基于全约束最小二乘法的自适应光谱解混合的光学实现方法，其特征在于，每一行空间像元的光谱解混合包括以下步骤：

步骤1：将一行空间像元的辐射/反射光，分光色散成二维光谱；

步骤2：利用已知参考光谱库中的端元构建端元矩阵E，并根据将E扩展为E′，初始化R＝{1,2,…,m}，A_i＝0，s＝0，r＝1；P和R分别是设定的两个集合；A_i和s是1×b大小的一维矩阵；r为用于计数的中间变量；b为光谱波段数；m为端元数；

步骤3：提取E′^T的第r行系数矩阵，标记为R，并用R的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的；

步骤4：根据步骤3中给出的二维控制信号，利用空间光调制器对步骤1分光后的二维光谱进行光谱调制；

步骤5：将调制后光谱经反向合光系统聚合后投射到探测器上，探测器获取聚合后光谱；

步骤6：将探测器获取的数据与步骤3提取的一维系数矩阵R的第b+1个系数值相加，得到1×m矩阵E′^Tx_i′中第r个数值；

步骤7：如果r＝m，则已获取1×m矩阵E′^Tx_i′，转入步骤8；否则，更新r＝r+1，转入步骤3；

步骤8：根据公式w＝E′^T(x_i′-E′A_i)，计算w；

步骤9：如果且max(w_k)＞0(k∈R)，转入步骤10；否则，转入步骤20；

步骤10：求出t＝argmax(w_k)(k∈R)，更新R＝R-{t}和P＝P∪{t}；

步骤11：根据更新的集合P，更新E′^P，并利用公式E^P#＝[(E′^P)^TE′^P]^-1(E′^P)^T自适应地计算E^P#,pr＝1；pr为用于计数的中间变量；

步骤12：提取E^P#的第pr行系数矩阵，标记为PR，并用PR的前b个系数构成空间光调制器的二维控制信号，其中空间光调制器的各行控制信号是相同的；

步骤13：根据步骤12中给出的二维控制信号，利用空间光调制器对步骤1分光后的二维光谱进行光谱调制；

步骤14：将调制后光谱经反向合光系统聚合后投射到探测器上，探测器获取聚合后光谱；

步骤15：将探测器获取的数据与步骤12提取的一维系数矩阵PR的第b+1个系数值相加，得到1×m矩阵E^P#x_i′中第pr个数值；

步骤16：如果pr＝m，则已获取1×m矩阵s^P＝E^P#x_i′，并更新s^R＝0，转步骤17；否则，更新pr＝pr+1，转步骤12；

步骤17：如果min(s^P)≤0，转入步骤18；否则，转入步骤19；

步骤18：根据公式A_i＝A_i+β(s-A_i)，计算β和A_i，并把所有满足A_ik＝0的指标k从P移动到R，转入步骤11；

步骤19：更新A_i＝s，w＝E′^T(x_i′-E′A_i)，转入步骤9；

步骤20：当前行的光谱解混合完成，此时A_i为当前行的丰度系数矩阵，进入下一行的光谱解混合过程。

2.基于权利要求1所述光学实现方法的自适应光谱解混合系统，其特征在于：包括CPU，以及依次设置的分光元件、二维空间光调制器、反向合光系统和探测器；

所述分光元件用于将一行空间像元的辐射/反射光色散成二维光谱；二维空间光调制器用于根据CPU发出的二维控制信号对经分光元件分光后的光谱进行光谱调制；所述反向合光系统用于对调制后光谱进行聚合；所述探测器用于获取聚合后光谱数据，并发送给所述CPU，该CPU与计算机可读存储介质相连，其上存储有计算机程序；

所述计算机程序被CPU执行时实现权利要求1中的步骤2至步骤3、步骤6至步骤12、步骤15至步骤20。