CN104463223A

CN104463223A - Hyperspectral image group sparse demixing method based on empty spectral information abundance restraint

Info

Publication number: CN104463223A
Application number: CN201410810715.3A
Authority: CN
Inventors: 张向荣; 焦李成; 吴健康; 马文萍; 马晶晶; 侯彪; 白静; 刘红英
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2014-12-22
Filing date: 2014-12-22
Publication date: 2015-03-25
Anticipated expiration: 2034-12-22
Also published as: CN104463223B

Abstract

The invention belongs to the technical field of image processing, and specifically discloses a hyperspectral image group sparse unmixing method based on spatial spectrum information abundance constraints. The implementation steps are: input hyperspectral image data sets and standard spectral libraries; Adaptive grouping of spectral image data; perform group sparse unmixing for each group of hyperspectral image data; use the abundance matrix corresponding to each group of hyperspectral image data as a basis to trim the spectral library; output hyperspectral image data sparse unmixing results. The present invention considers the structural characteristics of hyperspectral data and spectral library data, uses group sparse unmixing and spectral library pruning strategies for hyperspectral image unmixing, and improves the accuracy of hyperspectral data sparse unmixing in the prior art. The method of the invention is based on a semi-supervised hyperspectral image unmixing method, has the advantages of sparse unmixing and high unmixing precision, and is an effective hyperspectral image unmixing method.

Description

Sparse Unmixing Method for Hyperspectral Image Groups Based on Spatial Spectral Information Abundance Constraints

技术领域technical field

本发明属于图像处理技术领域，更进一步涉及遥感图像技术领域，具体是一种基于空谱信息丰度约束的高光谱图像组稀疏解混方法。The invention belongs to the technical field of image processing, and further relates to the technical field of remote sensing images, in particular to a hyperspectral image group sparse unmixing method based on the abundance constraints of spatial spectrum information.

背景技术Background technique

高光谱遥感图像技术在近几年来有很快的发展，其研究主要致力于寻找使计算机智能地学习和识别高光谱图像真实地物类的技术方法。高光谱图像在城市规划、环境检测、植被分类、军事目标探测以及矿物地质识别等诸多方面都有着巨大的应用前景。由于高光谱成像仪分辨率及地物地形复杂度的影响，使得图像中单个像素存在着多种物质的混合，从而形成混合像素，影响了对高光谱图像的进一步解译和应用，因此，高光谱图像解混技术成为当前遥感领域最有研究意义的一门课题。最为普遍的高光谱图像解混模型是线性混合模型，该模型描述高光谱图像中的混合像元是一组少量地物即端元的线性组合，并被噪声干扰，由于该模型简单实用、物理意义明确，常被作为基本模型来解决高光谱解混问题。高光谱解混的过程是：首先，提取高光谱图像中存在的端元光谱信息；其次，获得混合像元中不同端元的比例信息即丰度信息。另外，由于高光谱图像具有数据量大、冗余信息多、含有噪声等不利因素，因此要求在高光谱解混技术方法具有一定的抗干扰和高精度的能力。Hyperspectral remote sensing image technology has developed rapidly in recent years, and its research is mainly devoted to finding technical methods to enable computers to intelligently learn and identify real objects in hyperspectral images. Hyperspectral images have great application prospects in many aspects such as urban planning, environmental detection, vegetation classification, military target detection, and mineral geological identification. Due to the influence of the resolution of the hyperspectral imager and the complexity of the ground and terrain, there is a mixture of various substances in a single pixel in the image, thus forming a mixed pixel, which affects the further interpretation and application of the hyperspectral image. Therefore, high Spectral image unmixing technology has become the most research topic in the field of remote sensing. The most common hyperspectral image unmixing model is the linear mixture model, which describes that the mixed pixels in the hyperspectral image are a linear combination of a small number of ground objects, namely end members, and are disturbed by noise. It has a clear meaning and is often used as a basic model to solve the problem of hyperspectral unmixing. The process of hyperspectral unmixing is: firstly, extract the endmember spectral information existing in the hyperspectral image; secondly, obtain the proportion information of different endmembers in the mixed pixel, that is, the abundance information. In addition, because hyperspectral images have disadvantages such as large data volume, redundant information, and noise, it is required to have certain anti-interference and high-precision capabilities in hyperspectral unmixing technology.

目前，高光谱解混方法主要包括基于几何、统计和稀疏回归的三大类方法。基于稀疏回归的高光谱解混方法作为一种半监督的方法越来越得到很多人的认可和研究，它避免了端元数目的估计和端元谱信息的提取，利用光谱库可以同时实现获取端元信息和丰度估计。由于高光谱图像中每个像元包含的端元数目远远小于光谱库中的维度，得到的稀疏表示的系数是稀疏的，符合数据稀疏表示的特性，因此，利用基于稀疏回归的高光谱解混方法可以得到有效的结果。At present, hyperspectral unmixing methods mainly include three categories of methods based on geometric, statistical and sparse regression. As a semi-supervised method, the hyperspectral unmixing method based on sparse regression has been more and more recognized and researched by many people. It avoids the estimation of the number of endmembers and the extraction of endmember spectral information, and the spectral library can be used to simultaneously obtain Endmember information and abundance estimation. Since the number of endmembers contained in each pixel in the hyperspectral image is much smaller than the dimension in the spectral library, the coefficients of the obtained sparse representation are sparse, which conforms to the characteristics of sparse representation of data. Therefore, using hyperspectral solution based on sparse regression Mixed methods can get effective results.

Marian-Daniel Iordache等人在论文“Sparse Unmixing of HyperspectralData”(GRS，2011)中提出了基本的稀疏解混，它是对高光谱图像中的每个混合像元分别独立的进行解混，同时获得端元信息和丰度信息，解决稀疏解混模型利用了变量分裂和增广拉格朗日的方法(SUnSAL),这个方法是基于交替方向乘子法的理论。Marian-Daniel Iordache et al. proposed the basic sparse unmixing in the paper "Sparse Unmixing of HyperspectralData" (GRS, 2011), which unmixes each mixed pixel in the hyperspectral image independently, and obtains Endmember information and abundance information, to solve the sparse unmixing model using the method of variable splitting and augmented Lagrangian (SUnSAL), this method is based on the theory of alternating direction multiplier method.

随后在此基础上，提出了一种改良的实现对高光谱图像所有像元同时施加稀疏约束的解混方法，由于高光谱图像中包含的所有端元数目远小于光谱库维数，得到的丰度矩阵是行稀疏的，也就是说光谱库中只有少量的谱特征用来拟合高光谱图像中的所有像元，通过解决协同稀疏回归的问题，提升了高光谱图像的解混效果，采用的是协同的变量分裂和增广拉格朗日的算法(CLSUnSAL)，它也是交替方向乘子法理论的具体实现。但是上述两种稀疏解混方法仍然存在的不足之处是：并未考虑高光谱图像数据的区域结构信息和进一步减弱光谱库中地物光谱间高相关性对稀疏解混的影响。Then, on this basis, an improved unmixing method is proposed to implement sparse constraints on all pixels of the hyperspectral image at the same time. Since the number of all endmembers contained in the hyperspectral image is much smaller than the dimension of the spectral library, the obtained abundant The degree matrix is row-sparse, which means that only a small number of spectral features in the spectral library are used to fit all pixels in the hyperspectral image. By solving the problem of collaborative sparse regression, the unmixing effect of the hyperspectral image is improved. is the cooperative variable splitting and augmented Lagrangian algorithm (CLSUnSAL), which is also a concrete realization of the theory of alternating direction multipliers. However, the shortcomings of the above two sparse unmixing methods are that they do not consider the regional structure information of the hyperspectral image data and further weaken the influence of the high correlation between object spectra in the spectral library on the sparse unmixing.

发明内容Contents of the invention

本发明的目的是克服上述现有技术的不足，提出一种基于空谱信息丰度约束的高光谱图像组稀疏解混方法，利用高光谱图像的区域结构信息，减弱光谱库中地物光谱高相关干性的影响，使高光谱图像稀疏解混效果进一步得到提高。The purpose of the present invention is to overcome the deficiencies of the above-mentioned prior art, and propose a hyperspectral image group sparse unmixing method based on the constraint of spatial spectral information abundance, which uses the regional structure information of hyperspectral images to weaken the spectral height of ground objects in the spectral library. The effect of correlation dryness further improves the hyperspectral image sparse unmixing effect.

本发明的技术方案是：Technical scheme of the present invention is:

(1)输入高光谱图像数据和标准光谱库其中为高光谱图像数据中的第i个像元的光谱特征，为标准光谱库中的第j个地物光谱，L为波段数，n为像元个数，m为光谱库中包含的地物光谱个数，表示实数域；(1) Input hyperspectral image data and standard spectral library in is the spectral feature of the i-th pixel in the hyperspectral image data, is the jth ground object spectrum in the standard spectral library, L is the number of bands, n is the number of pixels, m is the number of ground object spectra contained in the spectral library, represents the field of real numbers;

(2)利用均值漂移算法对高光谱图像数据Y进行自适应分割，分割为k个区域，每个区域即为一组，则Y＝{G₁,...G_i,...,G_k}，其中i＝1,…,k为高光谱图像数据分割得到的第i组像元集，|G_i|为第i组像元集内包含的像元个数；(2) Use the mean shift algorithm to adaptively segment the hyperspectral image data Y into k regions, each region is a group, then Y={G ₁ ,...G _i ,...,G _k }, where i=1,...,k is the i-th group of pixel sets obtained by segmentation of hyperspectral image data, |G _i | is the number of pixels contained in the i-th group of pixel sets;

(3)通过组稀疏解混的方法对步骤(2)中得到的每组像元集G_i，i＝1,…,k分别进行稀疏解混：(3) Perform sparse unmixing on each group of pixel sets G _i obtained in step (2), i=1,...,k by group sparse unmixing method:

3a)输入由高光谱图像数据分割得到的一组像元集G_i和光谱库A，设置初始迭代次数iter＝0，光谱库A中地物光谱个数t＝m；3a) Input a group of pixel sets G _i and spectral library A obtained by segmenting hyperspectral image data, set the initial iteration number iter=0, and the number of feature spectra in spectral library A t=m;

3b)利用组稀疏回归的模型求解像元集G_i对应的丰度矩阵X_i，数学模型为s.t X_i≥0，其中第一项是误差项，表示任意矩阵H的F范数，trace{·}为矩阵的迹，模型的第二项是矩阵的l_2,1范数，其约束了矩阵行稀疏，x^j表示丰度矩阵X_i的第j行，m为丰度矩阵的行数，λ为正则项参数；3b) Use the group sparse regression model to solve the abundance matrix X _i corresponding to the pixel set G _i , the mathematical model is st X _i ≥ 0, where the first term is the error term, Indicates the F norm of any matrix H, trace{ } is the trace of the matrix, the second item of the model is the l _2,1 norm of the matrix, which constrains the sparseness of the matrix rows, x ^j represents the jth row of the abundance matrix X _i , m is the number of rows of the abundance matrix, and λ is the regular term parameter;

3c)令迭代次数iter＝iter+1，利用步骤3b)获得像元集G_i对应的丰度矩阵X_i对光谱库A进行修剪，更新光谱库A，光谱库修剪的具体步骤为：3c) Make the number of iterations iter=iter+1, use step 3b) to obtain the abundance matrix X _i corresponding to the pixel set G _i to prune the spectral library A, update the spectral library A, the specific steps of spectral library pruning are:

第1步，输入丰度矩阵X_i＝[x¹,...x^j,...,x^m]^T和光谱库A＝[a₁,...a_j,...,a_m]，x^j为丰度矩阵的第j行向量，a_j为光谱库中的第j列光谱特征；Step 1, input abundance matrix X _i =[x ¹ ,...x ^j ,...,x ^m ] ^T and spectral library A=[a ₁ ,...a _j ,...,a _m ], x ^j is the jth row vector of the abundance matrix, and a _j is the jth column spectral feature in the spectral library;

第2步，判断行向量x^j中的所有元素是否均大于阈值τ，初始阈值τ大小一般取2×e^-3，若满足条件，则保留光谱库中对应的a_j，否则，从光谱库A中剔除a_j；The second step is to judge whether all the elements in the row vector x ^j are greater than the threshold τ. The initial threshold τ is generally 2×e ^-3 . If the condition is met, the corresponding a _j in the spectral library is retained; otherwise, the spectral library Eliminate a _j from A;

第3步，得到新的光谱库A，更新阈值τ＝iter*τ和参数t＝|A|，|A|为新的光谱库A中包含的地物光谱个数；Step 3, get a new spectral library A, update the threshold τ=iter*τ and parameter t=|A|, where |A| is the number of object spectra contained in the new spectral library A;

3d)重复步骤3b)–3c)，直至满足终止条件，终止条件为最大迭代次数iter＝20或者光谱库中剩余地物光谱个数t≤T，阈值T的取值范围为p＜T＜2p，p为高光谱图像数据中包含的端元数目，可以由基于最小误差的高光谱信号识别算法估计得到；3d) Repeat steps 3b)–3c) until the termination condition is met. The termination condition is the maximum number of iterations iter=20 or the number of remaining object spectra in the spectral library t≤T, and the value range of the threshold T is p<T<2p , p is the number of endmembers contained in the hyperspectral image data, which can be estimated by the hyperspectral signal recognition algorithm based on the minimum error;

(4)获得k组像元集{G₁,...G_i,...,G_k}对应的丰度矩阵{X₁,...X_i,...,X_k}，其中X_i为第i组像元集G_i对应的丰度矩阵，输出高光谱图像数据解混结果。(4) Obtain the abundance matrix {X ₁ ,...X _i ,...,X _k } corresponding to k groups of pixel sets {G ₁ ,...G _i ,...,G _k }, where X _i is the abundance matrix corresponding to the i-th pixel set G _i , and outputs the hyperspectral image data unmixing result.

上述步骤3b)中数学模型的具体求解步骤为：The specific solution steps of the mathematical model in the above step 3b) are:

第1步，对于数学模型s.t X_i≥0中已知的变量G_i，A，λ做预处理，令G_i＝G_i/const，A＝A/const，λ＝λ/const，const表示一个中间变量，|G_i|表示像元集G_i内包含的像元个数；Step 1, for the mathematical model The known variables G _i , A, and λ in st X _i ≥ 0 are preprocessed, so that G _i =G _i /const, A=A/const, λ=λ/const, const represents an intermediate variable, |G _i | represents the number of pixels contained in the pixel set G _i ;

第2步，引入辅助变量U，使得U＝X_i,则数学模型等价于如下形式：In the second step, the auxiliary variable U is introduced so that U=X _i , then the mathematical model is equivalent to the following form:

$\underset{U u,, {V V}_{11},, {V V}_{22},, {V V}_{33}}{min min} \frac{11}{22} {| | | | {V V}_{11} - - G G | | | |}_{F f}^{22} + + λ λ {| | | | {V V}_{22} | | | |}_{2,1 2,1} + + {ι ι}_{R R + +} (({V V}_{33}))$

s.t V₁＝AU (1)st V ₁ =AU (1)

V₂＝UV ₂ =U

V₃＝UV ₃ =U

其中为指示函数，|V₃|为矩阵V₃的列数，V_3i为矩阵V₃的第i列，ι_R+(V_3i)的数学表达式如下所示：in For indicator function, | V ₃ | is the column number of matrix V ₃ , V _3i is the i-th column of matrix V ₃ , and the mathematical expression of ι _R+ (V _3i ) is as follows:

${ι ι}_{R R + +} (({V V}_{33 i i})) = = \{\begin{matrix} 00 & {V V}_{33 i i} &GreaterEqual; &Greater Equal; 00 \\ \infty \infty & {V V}_{33 i i} < < 00 \end{matrix} - - - - - - ((22))$

第3步，根据交替方向乘子法，引入增广拉格朗日乘子D₁/μ，D₂/μ，D₃/μ，μ为常数，公式(1)等价转换为如下形式：In the third step, according to the method of alternating direction multipliers, the augmented Lagrangian multipliers D ₁ /μ, D ₂ /μ, D ₃ /μ are introduced, and μ is a constant, and the formula (1) is equivalently transformed into the following form:

$\begin{matrix} Γ Γ ((U u,, {V V}_{11},, {V V}_{22},, {V V}_{33},, {D D.}_{11},, {D D.}_{22},, {D D.}_{33})) \\ = = \frac{11}{22} {| | | | {V V}_{11} - - G G | | | |}_{F f}^{22} + + λ λ {| | | | {V V}_{22} | | | |}_{2,1 2,1} + + {ι ι}_{R R + +} (({V V}_{33})) + + \\ \frac{μ μ}{22} {| | | | AU AU - - {V V}_{11} - - {D D.}_{11} | | | |}_{F f}^{22} + + \frac{μ μ}{22} {| | | | U u - - {V V}_{22} - - {D D.}_{22} | | | |}_{F f}^{22} + + \\ \frac{μ μ}{22} {| | | | U u - - {V V}_{33} - - {D D.}_{33} | | | |}_{F f}^{22} \end{matrix} - - - - - - ((33))$

第4步，设迭代次数η＝0，常数μ≥0，初始化U⁽⁰⁾，令初始值为零，固定其他变量，分别求变量U，V₁，V₂，V₃，D₁，D₂，D₃的值；Step 4, set the number of iterations η=0, constant μ≥0, initialize U ⁽⁰⁾ , let The initial value is zero, other variables are fixed, and the values of variables U, V ₁ , V ₂ , V ₃ , D ₁ , D ₂ , and D ₃ are calculated respectively;

1)固定变量V₁，V₂，V₃，D₁，D₂，D₃，求变量U的值1) Fixed variables V ₁ , V ₂ , V ₃ , D ₁ , D ₂ , D ₃ , find the value of variable U

${U u}^{((η η + + 11))} &LeftArrow; &LeftArrow; {arg arg min min}_{U u} Γ Γ ((U u,, {V V}_{11}^{((η η))},, {V V}_{22}^{((η η))},, {V V}_{33}^{((η η))},, {D D.}_{11}^{((η η))},, {D D.}_{22}^{((η η))},, {D D.}_{33}^{((η η))}))$

得：U^(η+1)←(A^TA+2I)^-1(A^Tξ₁+ξ₂+ξ₃)Get: U ^(η+1) ←(A ^T A+2I) ^-1 (A ^T ξ ₁ +ξ ₂ +ξ ₃ )

其中： $ξ_{1} = V_{1}^{(η)} + D_{1}^{(η)}, ξ_{2} = V_{2}^{(η)} + D_{2}^{(η)}, ξ_{3} = V_{3}^{(η)} + D_{3}^{(η)}$ in: $ξ_{1} = V_{1}^{(η)} + {D.}_{1}^{(η)}, ξ_{2} = V_{2}^{(η)} + {D.}_{2}^{(η)}, ξ_{3} = V_{3}^{(η)} + {D.}_{3}^{(η)}$

2)固定变量U，V₂，V₃，D₁，D₂，D₃，求变量V₁的值2) Fixed variables U, V ₂ , V ₃ , D ₁ , D ₂ , D ₃ , and find the value of variable V ₁

${V V}_{11}^{((η η + + 11))} &LeftArrow; &LeftArrow; {arg arg min min}_{{V V}_{11}} Γ Γ (({U u}^{((η η))},, {V V}_{11},, {V V}_{22}^{((η η))},, {V V}_{33}^{((η η))}))$

得： $V_{1}^{(η + 1)} &LeftArrow; \frac{1}{1 + μ} [G + μ ({AU}^{(η)} - D_{1}^{(η)})]$ have to: $V_{1}^{(η + 1)} &LeftArrow; \frac{1}{1 + μ} [G + μ ({AU}^{(η)} - {D.}_{1}^{(η)})]$

3)固定变量U，V₁，V₃，D₁，D₂，D₃，求变量V₂的值3) Fixed variables U, V ₁ , V ₃ , D ₁ , D ₂ , D ₃ , and find the value of variable V ₂

${V V}_{22}^{((η η + + 11))} &LeftArrow; &LeftArrow; {arg arg min min}_{{V V}_{22}} Γ Γ (({U u}^{((η η))},, {V V}_{11}^{((η η))},, {V V}_{22},, {V V}_{33}^{((η η))}))$

此处求解需分别计算矩阵的每行其中，为矩阵的第r行，函数vect_soft(·,·)是行向量软阈值函数，计算公式为：The solution here needs to calculate the matrix separately for each line of in, is the rth row of the matrix, the function vect_soft(·,·) is a row vector soft threshold function, and the calculation formula is:

4)固定变量U，V₁，V₂，D₁，D₂，D₃，求变量V₃的值4) Fixed variables U, V ₁ , V ₂ , D ₁ , D ₂ , D ₃ , and find the value of variable V ₃

${V V}_{33}^{((η η + + 11))} &LeftArrow; &LeftArrow; {arg arg min min}_{{V V}_{33}} Γ Γ (({U u}^{((η η))},, {V V}_{11}^{((η η))},, {V V}_{22}^{((η η))},, {V V}_{33}))$

得： $V_{3}^{(η + 1)} &LeftArrow; \max (U^{(η)} - D_{3}^{(η)}, 0)$ have to: $V_{3}^{(η + 1)} &LeftArrow; \max (u^{(η)} - {D.}_{3}^{(η)}, 0)$

5)更新变量D₁，D₂，D₃ 5) Update variables D ₁ , D ₂ , D ₃

${D D.}_{11}^{((η η + + 11))} &LeftArrow; &LeftArrow; {D D.}_{11}^{((η η))} - - {AU AU}^{((η η + + 11))} + + {V V}_{11}^{((η η + + 11))}$

${D D.}_{22}^{((η η + + 11))} &LeftArrow; &LeftArrow; {D D.}_{22}^{((η η))} - - {U u}^{((η η + + 11))} + + {V V}_{22}^{((η η + + 11))}$

${D D.}_{33}^{((η η + + 11))} &LeftArrow; &LeftArrow; {D D.}_{33}^{((η η))} - - {U u}^{((η η + + 11))} + + {V V}_{33}^{((η η + + 11))}$

第5步，令迭代次数η＝η+1，重复上述对变量U，V₁，V₂，V₃，D₁，D₂,D₃的求解过程，直至满足终止条件，终止条件为最大迭代次数η＝200或者变量U前后两次的误差阈值ε一般取值大小为1×e^-5，输出变量U的值，即丰度矩阵X_i。Step 5, set the number of iterations η=η+1, repeat the above process of solving the variables U, V ₁ , V ₂ , V ₃ , D ₁ , D ₂ , D ₃ until the termination condition is satisfied, which is the maximum iteration The number of times η=200 or the error of the variable U twice before and after The threshold ε generally takes a value of 1×e ^-5 , and outputs the value of the variable U, that is, the abundance matrix X _i .

本发明的有益效果：本发明可以通过对高光谱图像的自适应分组，合理的利用高光谱图像的区域结构信息，在此基础上，利用组稀疏解混的优势以及光谱库修剪的优点，实现对高光谱图像中实际存在的混合像元分离，达到对图像的进一步解译和应用。与现有技术相比，本发明具有以下优点：Beneficial effects of the present invention: the present invention can rationally utilize the regional structure information of hyperspectral images through self-adaptive grouping of hyperspectral images, and on this basis, utilize the advantages of group sparse unmixing and spectral library pruning to realize Separation of the mixed pixels that actually exist in the hyperspectral image can achieve further interpretation and application of the image. Compared with the prior art, the present invention has the following advantages:

第一，本发明把均值漂移分割算法应用于高光谱图像解混领域，它主要是作为高光谱解混的预处理手段，现有的高光谱图像稀疏解混和空谱信息的结合将有助于提高解混效果，利用高光谱图像的空谱先验信息，认为邻域内的高光谱像素具有相同的端元和相近的丰度值，而且实际高光谱图像是由不同的区域组成，将高光谱图像进行自适应分割有助于将要进行的组稀疏解混，使得本发明更加贴近实际应用。First, the present invention applies the mean shift segmentation algorithm to the field of hyperspectral image unmixing. It is mainly used as a preprocessing means for hyperspectral image unmixing. The combination of existing hyperspectral image sparse unmixing and spatial spectral information will help To improve the unmixing effect, use the spatial spectrum prior information of the hyperspectral image to consider that the hyperspectral pixels in the neighborhood have the same end members and similar abundance values, and the actual hyperspectral image is composed of different regions. The self-adaptive segmentation of the image is helpful for the sparse unmixing of the group to be performed, which makes the present invention closer to practical application.

第二，本发明利用组稀疏解混方法进行高光谱图像解混，将对高光谱图像自适应分割的不同区域作为不同像素组，利用组内高光谱像素拥有相同的组成成分和相近的比例的信息，进一步发挥了组稀疏解混的优势，提高了高光谱图像解混的精度，使得本发明拥有更好的解混效果。Second, the present invention uses the group sparse unmixing method to unmix the hyperspectral image, and uses the hyperspectral pixels in the group to have the same composition and similar proportions as different pixel groups. information, further exerts the advantages of group sparse unmixing, improves the accuracy of hyperspectral image unmixing, and makes the present invention have a better unmixing effect.

第三，本发明利用了丰度约束的光谱库修剪方法，剔除光谱库中部分光谱特征，降低了光谱库谱特征高相关性对稀疏解混的影响，使得本发明进一步提高了高光谱图像的解混精度。Third, the present invention utilizes an abundance-constrained spectral library pruning method to remove some spectral features in the spectral library, reducing the impact of high correlation of spectral features of the spectral library on sparse unmixing, making the present invention further improve the accuracy of hyperspectral images. Unmixing accuracy.

以下将结合附图对本发明做进一步详细说明。The present invention will be described in further detail below in conjunction with the accompanying drawings.

附图说明Description of drawings

图1为本发明方法的流程图；Fig. 1 is the flowchart of the inventive method;

图2为本发明仿真实验中高光谱图像；Fig. 2 is a hyperspectral image in the simulation experiment of the present invention;

图3为本发明仿真实验中高光谱图像中的一块区域；Fig. 3 is an area in the hyperspectral image in the simulation experiment of the present invention;

图4为本发明方法真实数据解混结果与其他方法的对比；Fig. 4 is the comparison of the real data unmixing result of the method of the present invention and other methods;

图4(a)是算法NCLS丰度估计图；Figure 4(a) is the algorithm NCLS abundance estimation diagram;

图4(b)是算法SUnSAL丰度估计图；Figure 4(b) is the abundance estimation diagram of the algorithm SUnSAL;

图4(c)是算法CLSUnSAL丰度估计图；Figure 4(c) is the abundance estimation diagram of the algorithm CLSUnSAL;

图4(d)是本发明方法丰度估计图。Fig. 4(d) is an abundance estimation diagram of the method of the present invention.

具体实施措施Specific implementation measures

实施例1：Example 1:

结合附图1对本发明的具体步骤描述如下：Concrete steps of the present invention are described as follows in conjunction with accompanying drawing 1:

模拟数据集大小为224×30×30，由九块三行三列同样大小的区域组成，每个区域大小为10×10，每个区域内包含的端元种类和数目不同，模拟数据集中包含的端元个数为9，端元从光谱库中随机选取，丰度服从狄利克雷分布。The size of the simulated data set is 224×30×30, and it consists of nine regions of the same size with three rows and three columns, each region is 10×10 in size, and the types and numbers of endmembers contained in each region are different. The simulated data set contains The number of endmembers is 9, the endmembers are randomly selected from the spectral library, and the abundance obeys the Dirichlet distribution.

真实数据集：由机载可见光/红外成像光谱仪AVIRIS获得的美国内达华地区矿物数据集(如附图2所示)，其中在实际仿真过程中我们利用其中一块区域，大小为250×191，每个像素拥有188个谱波段(如附图3所示)。Real data set: The mineral data set in the Nevada region of the United States obtained by the airborne visible/infrared imaging spectrometer AVIRIS (as shown in Figure 2), in which we use one of the areas in the actual simulation process, the size is 250×191, Each pixel has 188 spectral bands (as shown in Figure 3).

输入标准光谱库数据A，源自于美国地质调查局USGS在2007年发布的splib06光谱库，光谱库中所有地物光谱的光谱波段数为224。The input standard spectral library data A is derived from the splib06 spectral library released by the United States Geological Survey USGS in 2007. The number of spectral bands of all object spectra in the spectral library is 224.

3b)利用组稀疏回归的模型求解像元集G_i对应的丰度矩阵X_i，数学模型为s.t X_i≥0，其中第一项是误差项，表示任意矩阵H的F范数，trace{·}为矩阵的迹，模型的第二项是矩阵的l_2,1范数，其约束了矩阵行稀疏，使得非零行尽量少，x^j表示丰度矩阵X_i的第j行，m为丰度矩阵的行数，λ为正则项参数，一般通过手动调节的方式使最终的解混结果达到所需要求；3b) Use the group sparse regression model to solve the abundance matrix X _i corresponding to the pixel set G _i , the mathematical model is st X _i ≥ 0, where the first term is the error term, Indicates the F norm of any matrix H, trace{ } is the trace of the matrix, the second item of the model is the l _2,1 norm of the matrix, which constrains the sparseness of the matrix rows, so that there are as few non-zero rows as possible, x ^j represents the jth row of the abundance matrix X _i , m is the number of rows of the abundance matrix, and λ is the regular term Parameters, generally through manual adjustment to make the final unmixing result meet the required requirements;

根据上述基于空谱信息丰度约束的高光谱图像组稀疏解混方法，步骤3b)中数学模型的具体求解步骤为：According to the above hyperspectral image set sparse unmixing method based on the constraint of spatial spectral information abundance, the specific solution steps of the mathematical model in step 3b) are:

s.t V₁＝AU (1)st V ₁ =AU (1)

V₂＝UV ₂ =U

V₃＝UV ₃ =U

此处求解需分别计算矩阵的每行其中，为矩阵的第r行，函数vectsoft(·,·)是行向量软阈值函数，计算公式为：The solution here needs to calculate the matrix separately for each line of in, is the rth row of the matrix, and the function vectsoft(·,·) is a row vector soft threshold function, and the calculation formula is:

5)更新变量D₁，D₂，D₃ 5) Update variables D ₁ , D ₂ , D ₃

实施例2：Example 2:

下面结合附图3和附图4对本发明的效果做进一步描述。The effect of the present invention will be further described in conjunction with accompanying drawings 3 and 4 below.

本发明的仿真实验是在Intel Core(TM)2Duo CPU、主频2.00GHz，内存2G，Windows 7平台上的MATLAB R2011b上实现的。Simulation experiment of the present invention is realized on the MATLAB R2011b on the Windows 7 platform on Intel Core (TM) 2Duo CPU, main frequency 2.00GHz, internal memory 2G.

本发明的仿真是在模拟数据集和真实数据集上做的实验仿真，模拟数据是由10×10的九个小块组成，每块包含的端元数目不同，端元从九个谱特征中随机选取，所有小块的丰度服从狄利克雷分布，模拟数据集，大小为224×900，数据被不同级别高斯白噪声所干扰，信噪比SNR(dB)＝Ε||Ax||²/Ε||n||²分别为：20dB，30dB和40dB。The simulation of the present invention is an experimental simulation done on a simulated data set and a real data set. The simulated data is composed of nine small blocks of 10×10, and each block contains different numbers of endmembers. Randomly selected, the abundance of all small blocks obeys the Dirichlet distribution, the simulated data set, the size is 224×900, the data is interfered by different levels of Gaussian white noise, the signal-to-noise ratio SNR (dB) = Ε||Ax|| ² /Ε||n|| ² are: 20dB, 30dB and 40dB respectively.

真实数据是美国内达华地区的矿区数据集(如附图2所示)，其中在实际仿真过程中我们利用其中一块区域，大小为250×191，去除波段中一些被水蒸气和噪声严重干扰的波段，剩余拥有188个谱波段(如附图3所示)The real data is the mining data set in the Nevada region of the United States (as shown in Figure 2). In the actual simulation process, we use one of the areas with a size of 250×191 to remove some of the bands that are severely interfered by water vapor and noise. The remaining bands have 188 spectral bands (as shown in Figure 3)

模拟数据集解混精度评价指标是信号与重建误差信号比(SRE)，数学表达式如下所示：The evaluation index of the unmixing accuracy of the simulated data set is the signal-to-reconstruction error signal ratio (SRE), and the mathematical expression is as follows:

$SRE SREs ((dB dB)) = = 2020 log log ((E E. [[{| | | | X x | | | |}_{F f}^{22}]] / / E E. [[{| | | | X x - - \overset{&OverBar; &OverBar;}{X x} | | | |}_{F f}^{22}]]))$

其中：X为真实丰度矩阵，为估计丰度矩阵。SRE值越大，表示解混效果越好。Where: X is the true abundance matrix, is the estimated abundance matrix. The larger the SRE value, the better the unmixing effect.

本发明在模拟数据集上做的高光谱解混效果与非负约束最小二乘(NCLS)、SUnSAL、CLSUnSAL得到的结果做了对比，得到的高光谱解混效果对比如表1所示。The hyperspectral unmixing effect of the present invention on the simulated data set is compared with the results obtained by non-negative constrained least squares (NCLS), SUnSAL, and CLSUnSAL, and the obtained hyperspectral unmixing effect is compared as shown in Table 1.

表1不同方法在模拟数据集上的解混精度Table 1 Unmixing accuracy of different methods on simulated datasets

SNRSNR NCLSNCLS SUnSALSUnSAL CLSUnSALCLSUnSAL 本发明方法The method of the invention 20dB20dB 4.3704dB4.3704dB 5.5499dB5.5499dB 6.2458dB6.2458dB 7.0340dB7.0340dB 30dB30dB 9.6123dB9.6123dB 11.3869dB11.3869dB 12.5738dB12.5738dB 14.9755dB14.9755dB 40dB40dB 16.0141dB16.0141dB 19.3443dB19.3443dB 21.8499dB21.8499dB 27.3691dB27.3691dB

从表1可以看出，在模拟数据集上的仿真结果，在不同噪声下本发明的解混精度比已有的稀疏解混方法高。It can be seen from Table 1 that the simulation results on the simulated data set show that the unmixing accuracy of the present invention is higher than that of existing sparse unmixing methods under different noises.

本发明在美国内达华地区的矿区数据集上和不同算法NCLS、SUnSAL、CLSUnSAL的仿真实验，通过实验结果(如附图4所示)可以看出本发明能够获得较高的解混效果，得到的矿物(从左至右矿物质依次为alunit-明矾石、buddingtonite-水铵长石、chalcedony-玉髓)丰度图更加清晰明了，每种矿物质分布更加集中。The present invention is on the mining area data set of U.S. Nevada area and the emulation experiment of different algorithms NCLS, SUnSAL, CLSUnSAL, can find out that the present invention can obtain higher unmixing effect by experimental result (as shown in accompanying drawing 4), The obtained mineral abundance map (from left to right is alunit-alumite, buddingtonite- feldspar, chalcedony-chalcedony) is more clear, and the distribution of each mineral is more concentrated.

综上所述，本发明可以通过对高光谱图像的自适应分组，合理的利用高光谱图像的区域结构信息，在此基础上，利用组稀疏解混的优势以及光谱库修剪的优点，实现对高光谱图像中实际存在的混合像元分离，达到对图像的进一步解译和应用。与现有技术相比，本发明具有以下优点：In summary, the present invention can rationally utilize the region structure information of hyperspectral images through self-adaptive grouping of hyperspectral images. The separation of the mixed pixels actually existing in the hyperspectral image can achieve further interpretation and application of the image. Compared with the prior art, the present invention has the following advantages:

因此，本发明是一种比较有效的高光谱图像稀疏解混方法。Therefore, the present invention is a relatively effective hyperspectral image sparse unmixing method.

本实施方式中没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。The parts that are not described in detail in this embodiment are commonly known and commonly used means in this industry, and will not be described here one by one. The above examples are only illustrations of the present invention, and do not constitute a limitation to the protection scope of the present invention. All designs that are the same as or similar to the present invention fall within the protection scope of the present invention.

Claims

1. The hyperspectral image group sparse unmixing method based on the spatial spectrum information abundance constraint is characterized by comprising the following steps of: the method comprises the following steps:

(1) inputting hyperspectral image dataAnd a standard spectrum libraryWhereinIs the spectral feature of the ith pixel in the hyperspectral image data,is the jth surface feature spectrum in the standard spectrum library, L is the number of wave bands, n is the number of pixels, m is the number of surface feature spectra contained in the spectrum library,representing a real number domain;

(2) self-adaptive segmentation is carried out on the hyperspectral image data Y by using a mean shift algorithm, the hyperspectral image data Y is segmented into k regions, each region is a group, and then Y is equal to { G ═ G₁,...G_i,...,G_kTherein ofi is 1, …, k is the i-th group of image elements, | G, obtained by splitting the hyperspectral image data_iI is the number of pixels contained in the ith group of pixel sets;

(3) performing unmixing on each group of pixel sets G obtained in the step (2) by using a group sparseness unmixing method_iI-1, …, k are each subjected to sparse unmixing:

3a) inputting a group of image element sets G obtained by segmenting hyperspectral image data_iSetting the initial iteration number iter to be 0 and the number t of the surface feature spectrums in the spectrum library A to be m;

3b) solving a set of image elements G using a model of group sparse regression_iCorresponding abundance matrix X_iThe mathematical model is

<math> <mrow> <munder> <mi>min</mi> <msub> <mi>X</mi> <mi>i</mi> </msub> </munder> <mn>1</mn> <mo>/</mo> <mn>2</mn> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>AX</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>G</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>λ</mi> <msub> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>X</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mn>2,1</mn> </msub> <mi>s</mi> <mo>.</mo> <mi>t</mi> <msub> <mi>X</mi> <mi>i</mi> </msub> <mo>&GreaterEqual;</mo> <mn>0</mn> <mo>,</mo> </mrow> </math>

The first item ofIs the term for the error as a function of,f norm representing arbitrary matrix H, trace {. cndot } is the trace of the matrix, the second term of the modelIs the matrix of_2,1Norm, which constrains the matrix row sparsity, x^jRepresenting the abundance matrix X_iM is the number of rows of the abundance matrix, and lambda is a regular term parameter;

3c) let iter be iter +1, obtain primitive set G using step 3b)_iCorresponding abundance matrix X_iPruning the spectrum library A, and updating the spectrum library A, wherein the concrete steps of pruning the spectrum library are as follows:

step 1, inputting abundance matrix X_i＝[x¹,...x^j,...,x^m]^TAnd spectral library A ═ a₁,...a_j,...,a_m]，x^jIs the j-th row vector of the abundance matrix, a_jIs the jth column spectral feature in the spectral library;

step 2, judging the row vector x^jWhether all elements in (1) are larger than a threshold value tau, the size of the initial threshold value tau is generally 2 × e^-3If the condition is satisfied, the corresponding a in the spectrum library is reserved_jOtherwise, a is removed from the spectrum library A_j；

Step 3, obtaining a new spectrum library A, updating a threshold value tau and a parameter t | A |, wherein | A | is the number of the surface feature spectrums contained in the new spectrum library A;

3d) repeating the steps 3b) to 3c) until a termination condition is met, wherein the termination condition is that the maximum iteration time iter is 20 or the number T of the residual ground object spectrums in the spectrum library is less than or equal to T, the value range of the threshold value T is that p is less than T and less than 2p, and p is the number of end members contained in the hyperspectral image data and can be obtained by estimation through a hyperspectral signal identification algorithm based on the minimum error;

(4) obtaining k sets of image elements G₁,...G_i,...,G_kCorresponding abundance matrix { X }₁,...X_i,...,X_kIn which X is_iAs the ith group of image elements G_iAnd outputting a hyperspectral image data unmixing result by using the corresponding abundance matrix.

2. The sparse unmixing method for hyperspectral image groups based on empty spectral information abundance constraints as claimed in claim 1, wherein the specific solving step of the mathematical model in the step 3b) is:

step 1, for the mathematical model

<math> <mrow> <munder> <mi>min</mi> <msub> <mi>X</mi> <mi>i</mi> </msub> </munder> <mn>1</mn> <mo>/</mo> <mn>2</mn> <msubsup> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>AX</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>G</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>F</mi> <mn>2</mn> </msubsup> <mo>+</mo> <mi>λ</mi> <msub> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>X</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mn>2,1</mn> </msub> <mi>s</mi> <mo>.</mo> <mi>t</mi> <msub> <mi>X</mi> <mi>i</mi> </msub> <mo>&GreaterEqual;</mo> <mn>0</mn> </mrow> </math>

Known variable G in_iA, lambda is pretreated to obtainG_i＝G_iA/const, a ═ a/const, λ ═ λ/const, const denotes an intermediate variable, | G_iI represents the image element set G_iThe number of pixels contained in the image;

step 2, introducing an auxiliary variable U, so that U is equal to X_iThen the mathematical model is equivalent to the following form:

s.t V₁＝AU (1)

V₂＝U

V₃＝U

whereinTo indicate the function, | V₃I is the matrix V₃Number of rows of (V)_3iIs a matrix V₃Column i, iota of_R+(V_3i) The mathematical expression of (a) is as follows:

<math> <mrow> <msub> <mi>ι</mi> <mrow> <mi>R</mi> <mo>+</mo> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>V</mi> <mrow> <mn>3</mn> <mi>i</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>V</mi> <mrow> <mn>3</mn> <mi>i</mi> </mrow> </msub> <mo>&GreaterEqual;</mo> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mo>∞</mo> </mtd> <mtd> <msub> <mi>V</mi> <mrow> <mn>3</mn> <mi>i</mi> </mrow> </msub> <mo><</mo> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> </math>

step 3, introducing an augmented Lagrange multiplier D according to an alternate direction multiplier method₁/μ，D₂/μ，D₃Mu, mu is a constant, equation (1) is equivalently transformed into the following form:

step 4, setting the iteration number eta to be 0 and the constant mu to be more than or equal to 0, and initializing U⁽⁰⁾Let us order Setting the initial value to zero, fixing other variables, and respectively solving the variables U and V₁，V₂，V₃，D₁，D₂，D₃A value of (d);

1) fixed variable V₁，V₂，V₃，D₁，D₂，D₃Evaluating the value of the variable U

<math> <mrow> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <mo>&LeftArrow;</mo> <mi>arg</mi> <msub> <mi>min</mi> <mi>U</mi> </msub> <mi>Γ</mi> <mrow> <mo>(</mo> <mi>U</mi> <mo>,</mo> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> </math>

Obtaining: u shape^(η+1)←(A^TA+2I)^-1(A^Tξ₁+ξ₂+ξ₃)

Wherein:

2) fixed variables U, V₂，V₃，D₁，D₂，D₃Calculating the variable V₁Value of (A)

<math> <mrow> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <mi>arg</mi> <msub> <mi>min</mi> <msub> <mi>V</mi> <mn>1</mn> </msub> </msub> <mi>Γ</mi> <mrow> <mo>(</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msup> <mo>,</mo> <msub> <mi>V</mi> <mn>1</mn> </msub> <mo>,</mo> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> </math>

Obtaining:

<math> <mrow> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <mfrac> <mn>1</mn> <mrow> <mn>1</mn> <mo>+</mo> <mi>μ</mi> </mrow> </mfrac> <mo>[</mo> <mi>G</mi> <mo>+</mo> <mi>μ</mi> <mrow> <mo>(</mo> <msup> <mi>AU</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msup> <mo>-</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>]</mo> </mrow> </math>

3) fixed variables U, V₁，V₃，D₁，D₂，D₃Calculating the variable V₂Value of (A)

<math> <mrow> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <mi>arg</mi> <msub> <mi>min</mi> <msub> <mi>V</mi> <mn>2</mn> </msub> </msub> <mi>Γ</mi> <mrow> <mo>(</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msub> <mi>V</mi> <mn>2</mn> </msub> <mo>,</mo> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> </math>

Here, the solution requires separate computation of the matricesEach row ofWherein, for the r-th row of the matrix, the function vectsoft (·, · s) is a row vector soft threshold function, and the calculation formula is:

4) fixed variables U, V₁，V₂，D₁，D₂，D₃Calculating the variable V₃Value of (A)

<math> <mrow> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <mi>arg</mi> <msub> <mi>min</mi> <msub> <mi>V</mi> <mn>3</mn> </msub> </msub> <mi>Γ</mi> <mrow> <mo>(</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msub> <mi>V</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </math>

Obtaining:

<math> <mrow> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <mi>max</mi> <mrow> <mo>(</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msup> <mo>-</mo> <msubsup> <mi>D</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </math>

5) updating variable D₁，D₂，D₃

<math> <mrow> <msubsup> <mi>D</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>-</mo> <msup> <mi>AU</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <mo>+</mo> <msubsup> <mi>V</mi> <mn>1</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> </mrow> </math>

<math> <mrow> <msubsup> <mi>D</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <msubsup> <mi>D</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>-</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <mo>+</mo> <msubsup> <mi>V</mi> <mn>2</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> </mrow> </math>

<math> <mrow> <msubsup> <mi>D</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>&LeftArrow;</mo> <msubsup> <mi>D</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>)</mo> </mrow> </msubsup> <mo>-</mo> <msup> <mi>U</mi> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <mo>+</mo> <msubsup> <mi>V</mi> <mn>3</mn> <mrow> <mo>(</mo> <mi>η</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> </mrow> </math>

Step 5, making the iteration number eta equal to eta +1, and repeating the pair of variables U and V₁，V₂，V₃，D₁，D₂,D₃Until a termination condition is satisfied, the termination condition is that the maximum iteration number eta is 200 or the error of the variable U is twice before and afterThe threshold value is generally 1 × e^-5The value of the output variable U, i.e. the abundance matrix X_i。