CN103679210A

CN103679210A - Ground object recognition method based on hyperspectral image unmixing

Info

Publication number: CN103679210A
Application number: CN201310647509.0A
Authority: CN
Inventors: 杨淑媛; 焦李成; 黄春海; 马晶晶; 马文萍; 侯彪; 刘芳; 程时倩; 马永刚
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2013-12-03
Filing date: 2013-12-03
Publication date: 2014-03-26
Anticipated expiration: 2033-12-03
Also published as: CN103679210B

Abstract

The invention discloses a ground feature recognition method based on hyperspectral image unmixing, which mainly solves the problem of inaccurate judgment of the feature category to which mixed pixel points belong in the existing method. The implementation steps are as follows: input a hyperspectral image, arrange the mixed pixels in the hyperspectral image into a matrix by columns to form a data matrix; use the manifold constraints of the data matrix, the sparse constraints of the abundance matrix and the endmembers Constraints formed by matrix smoothing constraints are added to the objective function of the NMF algorithm to form a new objective function; the new objective function is optimized and unmixed to obtain the endmember matrix and abundance matrix after unmixing the hyperspectral image ; According to the endmember matrix and abundance matrix after unmixing, judge the object category of all mixed pixels in the hyperspectral image. The invention can improve the accuracy of the endmember value and the abundance value obtained by unmixing, thereby improving the accuracy of the hyperspectral image ground object recognition, and can be used for target tracking.

Description

Object Recognition Method Based on Hyperspectral Image Unmixing

技术领域 technical field

本发明属于遥感图像处理技术领域，是一种基于高光谱图像解混的地物识别方法，该方法可用于高光谱图像的分析，将一个混合像素点分解为端元和对应的丰度值。 The invention belongs to the technical field of remote sensing image processing, and relates to a ground feature recognition method based on hyperspectral image unmixing. The method can be used in the analysis of hyperspectral images, and decomposes a mixed pixel point into endmembers and corresponding abundance values. the

背景技术 Background technique

高光谱图像是利用成像光谱仪对同一地表区域的几十乃至几百个波段同时成像而获得的三维图像，由二维空间信息和一维光谱信息组成。利用这些丰富的光谱信息对地物进行细分和鉴别，在多领域得到了广泛应用。因为高光谱图像的光谱传感器的光谱分辨率高，所以形成的谱段多，但是每个谱段接受的能量小，所以只能提高接受光谱的地面面积，即降低空间分辨率。所以由于高光谱图像的空间分辨率不高还有自然界地物的复杂性从而形成混合像素点。混合像素点的普遍存在不仅影响地物的识别和分类精度，而且是遥感技术向定量化发展的重要障碍。因此如何有效地进行混合像素点分解是高光谱遥感图像应用的关键问题之一。 A hyperspectral image is a three-dimensional image obtained by simultaneously imaging dozens or even hundreds of bands of the same surface area with an imaging spectrometer, and consists of two-dimensional spatial information and one-dimensional spectral information. Using these rich spectral information to subdivide and identify ground features has been widely used in many fields. Because the spectral resolution of the spectral sensor of the hyperspectral image is high, there are many spectral segments formed, but the energy received by each spectral segment is small, so the ground area for receiving the spectrum can only be increased, that is, the spatial resolution is reduced. Therefore, due to the low spatial resolution of hyperspectral images and the complexity of natural objects, mixed pixels are formed. The ubiquity of mixed pixels not only affects the recognition and classification accuracy of ground objects, but also is an important obstacle to the quantitative development of remote sensing technology. Therefore, how to effectively decompose mixed pixels is one of the key issues in the application of hyperspectral remote sensing images. the

高光谱图像中混合像素点的模型一般采用线型混合模型，它的优点是算法简单，物理意义明确。该模型的数学过程简要描述如下：一个具有L个谱段的像素点表示为X_ij∈R^L×1，具有P个端元的端元矩阵表示为M∈R^L×P，M对应的丰度矩阵表示为S_ij∈R^P×1。则有：X_ij＝MS_ij+n，其中n为噪声。在实际环境中该模型受两个条件的限制：①M_up≥0,(1≤u≤L,1≤p≤P)②

Figure 2013106475090100002DEST_PATH_IMAGE001

上面的两个式子分别表示光谱的能量不存在负值和混合的能量的大小是一定的，不可能无限大。上述的模型和限制条件都满足非负矩阵分解（Nonnegative Matrix Factorization，NMF）的数学模型，所以可以用NMF算法来进行解混。 The model of mixed pixels in hyperspectral images generally adopts the linear mixed model, which has the advantages of simple algorithm and clear physical meaning. The mathematical process of the model is briefly described as follows: a pixel point with L spectral segments is expressed as X _ij ∈ R ^L×1 , an endmember matrix with P endmembers is expressed as M∈R ^L×P , and M corresponds to a rich The degree matrix is denoted as S _ij ∈ R ^P×1 . Then: X _ij =MS _ij +n, where n is noise. In the actual environment, the model is limited by two conditions: ①M _up ≥ 0, (1≤u≤L, 1≤p≤P)②

The above two formulas respectively indicate that there is no negative value in the energy of the spectrum and the magnitude of the mixed energy is certain and cannot be infinite. The above models and constraints all satisfy the mathematical model of Nonnegative Matrix Factorization (NMF), so the NMF algorithm can be used for unmixing.

目前提出的基于非负矩阵分解（Nonnegative Matrix Factorization，NMF）的高光谱图像的解混方法都是在NMF算法的目标函数上加上正则项。因为该方法没有充分考虑高光谱图像的特性，使得解混效果差，从而导致地物识别的精度不高。 The currently proposed methods for unmixing hyperspectral images based on Nonnegative Matrix Factorization (NMF) all add regular terms to the objective function of the NMF algorithm. Because this method does not fully consider the characteristics of hyperspectral images, the unmixing effect is poor, resulting in low accuracy of object recognition. the

发明内容 Contents of the invention

本发明的目的在于针对现有方法的不足，提出一种基于高光谱图像解混的地物识别方法，本方法将高光谱图像的端元矩阵的平滑性，丰度矩阵的稀疏性和数据矩阵的流形假设结构结合在一起，使得高光谱图像解混效果更好，从而使得地物识别的精度更高。 The purpose of the present invention is to address the deficiencies of the existing methods, and propose a ground object recognition method based on hyperspectral image unmixing. This method combines the smoothness of the endmember matrix of the hyperspectral image, the sparsity of the abundance matrix and the data matrix The combination of the manifold hypothetical structure makes the hyperspectral image unmixing effect better, so that the accuracy of object recognition is higher. the

实现本方法的技术方案是：输入一幅高光谱图像，将该高光谱图像中的混合像素点按列排成一个矩阵，组成数据矩阵，用数据矩阵的流形约束，丰度矩阵的稀疏约束和端元矩阵的平滑约束构成的约束项，加入到NMF算法的目标函数里，构成新的目标函数，然后对这个新的目标函数进行优化解混，得到该高光谱图像的端元矩阵和丰度矩阵，然后根据该解混出的端元矩阵和丰度矩阵判断该高光谱图像中所有混合像素点的地物类别。具体步骤包括如下： The technical solution to realize this method is: input a hyperspectral image, arrange the mixed pixels in the hyperspectral image into a matrix by columns, form a data matrix, use the manifold constraints of the data matrix, and the sparse constraints of the abundance matrix Constraints composed of the smooth constraints of the endmember matrix and the endmember matrix are added to the objective function of the NMF algorithm to form a new objective function, and then the new objective function is optimized and unmixed to obtain the endmember matrix and the richness of the hyperspectral image. Degree matrix, and then according to the unmixed endmember matrix and abundance matrix, the object categories of all mixed pixels in the hyperspectral image are judged. The specific steps include the following:

（1）输入一幅高光谱图像X∈R^M×N×L，并将该高光谱图像中混合像素点X_ij∈R^1×L按列排列，构成数据矩阵Z∈R^L×B，其中M和N为二维图像的行和列，i和为j为二维图像的横坐标和纵坐标，L为谱段数，B为高光谱图像中混合像素点总数，B＝M×N，R表示实数集合； (1) Input a hyperspectral image X ∈ R ^M×N×L , and arrange the mixed pixels X _ij ∈ R ^1×L in the hyperspectral image in columns to form a data matrix Z ∈ R ^L×B , where M and N are the rows and columns of the two-dimensional image, i and j are the abscissa and ordinate of the two-dimensional image, L is the number of spectral segments, B is the total number of mixed pixels in the hyperspectral image, B=M×N, R represents a set of real numbers;

（2）根据流形假设理论，构造数据矩阵Z的流形约束项： (2) According to the manifold hypothesis theory, construct the manifold constraint item of the data matrix Z:

$\frac{11}{22} {Σ Σ}_{i i,, j j = = 11}^{B B} {| | | | {s the s}_{i i} - - {s the s}_{j j} | | | |}^{22} {W W}_{ij ij} = = {Σ Σ}_{i i = = 11}^{B B} {s the s}_{i i}^{T T} {s the s}_{i i} {D D.}_{ii i} - - {Σ Σ}_{i i,, j j = = 11}^{B B} {s the s}_{i i}^{T T} {s the s}_{j j} {W W}_{ij ij}$

$= = Tr Tr (({SYS SYS}^{T T}))$

其中z_i是Z的第i列，z_j是Z的第j列，且z_i是z_j的k个近邻中的一个，S是丰度矩阵，s_i是S的第i列，W是Z的权值矩阵，W_ij是W的一个元素，

Figure 2013106475090100002DEST_PATH_IMAGE004

为z_i和z_j的权值，Tr(·)表示矩阵的迹，T表示矩阵的转置，D是Z的对角线权值矩阵，D_ii是D对角线上的一个元素，D_ii＝Σ_jW_ij，Y是流形系数矩阵，Y＝D-W； where z _i is the i-th column of Z, z _j is the j-th column of Z, and z _i is one of the k nearest neighbors of z _j , S is the abundance matrix, s _i is the i-th column of S, W is The weight matrix of Z, W _ij is an element of W,

is the weight of z _i and z _j , Tr( ) represents the trace of the matrix, T represents the transpose of the matrix, D is the diagonal weight matrix of Z, D _ii is an element on the diagonal of D, D _ii =Σ _j W _ij , Y is the manifold coefficient matrix, Y=DW;

（3）根据高光谱图像成像理论，在丰度矩阵S中加入L_1/2范数，得到稀疏约束表达式‖S‖_1/2，以作为丰度矩阵S的稀疏约束项； (3) According to the hyperspectral image imaging theory, the L _1/2 norm is added to the abundance matrix S, and the sparse constraint expression ‖S‖ _1/2 is obtained as the sparse constraint item of the abundance matrix S;

（4）根据高光谱图像成像理论，在端元矩阵M中加入Frobenius范数，得到平滑约束表达式

Figure 2013106475090100002DEST_PATH_IMAGE005

以作为端元矩阵M的平滑约束项； (4) According to the hyperspectral image imaging theory, the Frobenius norm is added to the endmember matrix M to obtain the smooth constraint expression

Take as the smooth constraint term of the endmember matrix M;

（5）将步骤(2)-(4)得到的三个约束项添加到NMF算法的目标函数

Figure 2013106475090100002DEST_PATH_IMAGE006

里，以构成新的目标函数： (5) Add the three constraints obtained in steps (2)-(4) to the objective function of the NMF algorithm

, to form a new objective function:

$f f'' ((M m,, S S)) = = \frac{11}{22} {| | | | Z Z - - MS MS | | | |}_{F f}^{22} + + α α {| | | | S S | | | |}_{11 / / 22} + + β β {| | | | M m | | | |}_{F f}^{22} + + γTr γTr (({SYS SYS}^{T T}))$

其中，α为丰度矩阵S的稀疏约束正则参数，β为端元矩阵M的平滑约束正则参数，γ为数据矩阵Z的流形约束正则参数； Among them, α is the sparse constraint regularization parameter of the abundance matrix S, β is the smooth constraint regularization parameter of the endmember matrix M, and γ is the manifold constraint regularization parameter of the data matrix Z;

（6）对步骤（5）得到的目标函数f’（M，S）用迭代乘法进行优化求解，得到高光谱图像X∈R^M×N×L的端元矩阵M和丰度矩阵S； (6) Optimize and solve the objective function f'(M, S) obtained in step (5) by iterative multiplication, and obtain the endmember matrix M and abundance matrix S of the hyperspectral image X∈R ^M×N×L ;

（7）将上述高光谱图像X∈R^M×N×L中混合像素点X_ij用步骤（6）求解得到的端元矩阵M和丰度向量s_i表示，即混合像素点X_ij＝Ms_i； (7) Express the mixed pixel point X _ij in the above hyperspectral image X∈R ^M×N×L by the endmember matrix M and abundance vector s _i obtained from the solution in step (6), that is, the mixed pixel point X _ij = Ms _i ;

（8）根据高光谱图像统计分布理论，由步骤（7）中的丰度向量s_i对混合像素点X_ij进行地物类别判断，即当max(s_i)=s_ai时，则判混合像素点X_ij属于第a类，得到该混合像素点的类别标签为v_ij=a，其中max(·)表示取向量中的最大值，a＝1,2,...,P表示该高光谱图像中相应的地物类别编号，P表示该高光谱图像中地物类别总数，s_ai是s_i的第a个元素； (8) According to the statistical distribution theory of hyperspectral images, the abundance vector s _i in step (7) is used to judge the category of mixed pixels X _ij , that is, when max(s _i )=s _ai , the mixed pixel point X ij is judged to be mixed The pixel point X _ij belongs to the class a, and the category label of the mixed pixel point is v _ij =a, where max( ) means taking the maximum value in the vector, and a=1,2,...,P means the height The number of the corresponding feature category in the spectral image, P represents the total number of feature categories in the hyperspectral image, s _ai is the ath element of s _i ;

（9）对上述高光谱图像X∈R^M×N×L中所有混合像素点用步骤（8）的操作进行地物类别判断，得到该高光谱图像X∈R^M×N×L的地物类别矩阵V∈R^M×N。 (9) Use step (8) to judge the object category for all the mixed pixels in the above hyperspectral image X∈RM ^×N×L , and obtain the object category of the hyperspectral image X∈RM ^×N×L Class matrix V ∈ R ^{M × N} .

本发明与现有技术相比具有以下优点： Compared with the prior art, the present invention has the following advantages:

1.本发明相比现有的方法更充分的考虑了高光谱图像的结构特性，将数据矩阵的流形约束加入到解混，与未加入流形约束的NMF算法，CNMF算法和分段光滑NMF算法相比，得到的丰度矩阵的精度更高，使得判断混合像素点所属地物类别更准确，从而达到较高的地物识别精度。 1. Compared with the existing method, the present invention fully considers the structural characteristics of the hyperspectral image, adds the manifold constraint of the data matrix to the unmixing, and does not add the NMF algorithm of the manifold constraint, the CNMF algorithm and the piecewise smooth Compared with the NMF algorithm, the accuracy of the obtained abundance matrix is higher, which makes it more accurate to judge the object category to which the mixed pixel points belong, so as to achieve a higher accuracy of object recognition. the

2.本发明同时将端元矩阵的平滑约束和丰度矩阵的稀疏约束加入到解混，与现有的GLNMF算法相比，更充分考虑了高光谱图像特性，大大提高了端元矩阵的精度，使得计算出的地物类别更接近真实地物类别，更加符合实际应用的要求。 2. The present invention adds the smooth constraint of the endmember matrix and the sparse constraint of the abundance matrix to the unmixing at the same time. Compared with the existing GLNMF algorithm, it fully considers the characteristics of the hyperspectral image and greatly improves the accuracy of the endmember matrix , so that the calculated object category is closer to the real object category, and more in line with the requirements of practical applications. the

附图说明 Description of drawings

图1是本发明的实现流程图； Fig. 1 is the realization flowchart of the present invention;

图2是本发明使用的高光谱图像示意图； Fig. 2 is the hyperspectral image schematic diagram that the present invention uses;

图3是本发明与现有几种算法解混出来的端元的光谱曲线和实际光谱曲线的对比示意图； Fig. 3 is the comparative schematic diagram of the spectral curve and the actual spectral curve of the end member that the present invention and several existing algorithms are unmixed;

图4是本发明与现有几种算法解混出的端元的光谱角距离SAD和解混出的丰度的均方根误差RMSE的平均30次误差棒对比图。 Fig. 4 is a comparison chart of the average 30 times error bars of the spectral angular distance SAD of the endmembers unmixed by the present invention and several existing algorithms and the root mean square error RMSE of the unmixed abundances. the

具体实施方式 Detailed ways

参照图1，本发明的具体实施步骤如下： With reference to Fig. 1, the specific implementation steps of the present invention are as follows:

步骤1：输入高光谱图像，构建数据矩阵，获得该高光谱图像的真实的地物类别矩阵，真实的端元矩阵和真实的丰度矩阵。 Step 1: Input a hyperspectral image, construct a data matrix, and obtain the real object category matrix, real endmember matrix and real abundance matrix of the hyperspectral image. the

1.1)输入如图2所示的高光谱图像，该图像大小为145×145，共有16类地物，该图像中的每一个混合像素点可看做是由200个谱段的光谱信息构成的光谱向量； 1.1) Input the hyperspectral image shown in Figure 2. The size of the image is 145×145, and there are 16 types of ground objects. Each mixed pixel in the image can be regarded as composed of spectral information of 200 spectral segments spectral vector;

1.2)将该高光谱图像X∈R^M×N×L中混合像素点X_ij∈R^1×L按列排列，构成数据矩阵Z∈R^L×B，其中，M和N为二维图像的行和列，i和j分别为二维图像的横坐标和纵坐标，L为谱段数，P为地物类别数，B为高光谱图像中混合像素点总数，B＝M×N，R表示实数集合，在本实例中L取值为200，B取值为21025，P取值为16； 1.2) The mixed pixels in the hyperspectral image X∈R ^M×N×L are arranged in columns X _ij ∈ R ^1×L to form a data matrix Z∈R ^L×B , where M and N are the two-dimensional image Row and column, i and j are the abscissa and ordinate of the two-dimensional image respectively, L is the number of spectral segments, P is the number of object categories, B is the total number of mixed pixels in the hyperspectral image, B=M×N, R represents A set of real numbers. In this example, the value of L is 200, the value of B is 21025, and the value of P is 16;

1.3)获得该高光谱图像的真实类别矩阵为V∈R^M×N，真实端元为M∈R^L×P，真实丰度为S∈R^P×B。 1.3) Obtain the true category matrix of the hyperspectral image as V∈R ^M×N , the true endmember as M∈R ^L×P , and the true abundance as S∈R ^P×B .

步骤2：构造约束项： Step 2: Construct constraints:

1.1）根据流形假设理论，构造数据矩阵Z的流形约束项： 1.1) According to the manifold hypothesis theory, construct the manifold constraint item of the data matrix Z:

1.1a）用热核函数计算数据矩阵Z的第i列z_i和数据矩阵Z的第j列z_j的权值W_ij： 1.1a) Use the thermokernel function to calculate the weight W ij _of the i-th column z _i of the data matrix Z and the j-th column z _j of the data matrix Z:

${W W}_{ij ij} = = {e e}^{- - \frac{{| | | | {z z}_{i i} - - {z z}_{j j} | | | |}^{22}}{σ σ}},,$

其中，σ为热核参数，取值为1； Among them, σ is the thermonuclear parameter with a value of 1;

1.1b）用L₂范数计算丰度矩阵S的第i列s_i和丰度矩阵S的第j列s_j的距离： 1.1b) Calculate the distance between the i-th column s _i of the abundance matrix S and the j-th column s _j of the abundance matrix S using the _L2 norm:

‖s_i-s_j‖²； ‖s _i -s _j ‖ ² ;

1.1c）根据流形假设理论，将步骤1.1a）的式子和步骤1.1b）的式子结合起来构成流形约束项： 1.1c) According to the manifold hypothesis theory, combine the formula of step 1.1a) and the formula of step 1.1b) to form the manifold constraint item:

$= = Tr Tr (({SYS SYS}^{T T}))$

其中z_i是z_j的k个近邻中的一个，W是Z的权值矩阵，Tr(·)表示矩阵的迹，T表示矩阵的转置，D是Z的对角线权值矩阵，D_ii是D对角线上的一个元素，D_ii＝Σ_jW_ij，Y是流形系数矩阵，Y＝D-W； where z _i is one of the k nearest neighbors of z _j , W is the weight matrix of Z, Tr( ) represents the trace of the matrix, T represents the transpose of the matrix, D is the diagonal weight matrix of Z, D _ii is an element on the diagonal of D, D _ii =Σ _j W _ij , Y is the manifold coefficient matrix, Y=DW;

1.2）根据高光谱图像成像理论，在丰度矩阵S中加入L_1/2范数，得到稀疏约束表达式‖S‖_1/2，以作为丰度矩阵S的稀疏约束项； 1.2) According to the hyperspectral image imaging theory, the L _1/2 norm is added to the abundance matrix S, and the sparse constraint expression ‖S‖ _1/2 is obtained as the sparse constraint item of the abundance matrix S;

1.3）根据高光谱图像成像理论，在端元矩阵M中加入Frobenius范数，得到平滑约束表达式

Figure 2013106475090100002DEST_PATH_IMAGE011

以作为端元矩阵M的平滑约束项。 1.3) According to the hyperspectral image imaging theory, the Frobenius norm is added to the endmember matrix M to obtain a smooth constraint expression

Take as the smoothness constraint term of the endmember matrix M.

步骤3：构造新的目标函数 Step 3: Construct a new objective function

1.1）将步骤2得到的三个约束项添加入到NMF算法的目标函数

Figure 2013106475090100002DEST_PATH_IMAGE012

里，以构成新的目标函数： 1.1) Add the three constraints obtained in step 2 to the objective function of the NMF algorithm

, to form a new objective function:

$f f'' ((M m,, S S)) = = \frac{11}{22} {| | | | Z Z - - MS MS | | | |}_{F f}^{22} + + α α {| | | | S S | | | |}_{11 / / 22} + + β β {| | | | M m | | | |}_{F f}^{22} + + γTr γTr (({SYS SYS}^{T T})),,$

其中，α为丰度矩阵S的稀疏约束正则参数，取值为2，β为端元矩阵M的平滑约束正则参数，取值为1，γ为数据矩阵Z的流形约束正则参数，取值为0.6。 Among them, α is the sparse constraint regularization parameter of the abundance matrix S, and the value is 2; β is the smooth constraint regularization parameter of the endmember matrix M, and the value is 1; γ is the manifold constraint regularization parameter of the data matrix Z, and the value is 0.6. the

步骤4：对步骤3得到的新的目标函数f’（M，S）用迭代乘法进行优化求解，得到端元矩阵M和丰度矩阵S。 Step 4: The new objective function f'(M, S) obtained in step 3 is optimized and solved by iterative multiplication, and the endmember matrix M and abundance matrix S are obtained. the

1.1）用[0,1]之间随机数初始化端元矩阵M和丰度矩阵S，将S的每列进行归一化处理； 1.1) Initialize the endmember matrix M and the abundance matrix S with random numbers between [0,1], and normalize each column of S;

1.2）输入端元矩阵M和丰度矩阵S； 1.2) Input endmember matrix M and abundance matrix S;

1.3）通过对新的目标函数f’（M，S）中的端元矩阵M和丰度矩阵S求导，分别得到端元矩阵M和丰度矩阵S的迭代乘法的计算公式： 1.3) By deriving the endmember matrix M and the abundance matrix S in the new objective function f'(M, S), the calculation formulas of the iterative multiplication of the endmember matrix M and the abundance matrix S are obtained respectively:

M’=M.*(ZS^T-βM)./MSS^T， M'=M.*(ZS ^T -βM)./MSS ^T ,

$S S . . * * (({M m}^{T T} Z Z + + 22 γSW γ SW)) . . / / (({M m}^{T T} MS MS + + \frac{11}{22} α α {S S}^{- - \frac{11}{22}} + + 22 γSD γSD))$

其中，（.*）表示矩阵的元素乘法，（./）表示矩阵的元素除法； Among them, (.*) represents the element multiplication of the matrix, and (./) represents the element division of the matrix;

1.4）用步骤1.3）的式子计算出的新的丰度矩阵S’和新的端元矩阵M’，将新的丰度矩阵S’和新的端元矩阵M’作为步骤1.2）的输入； 1.4) Use the new abundance matrix S' and new endmember matrix M' calculated by the formula in step 1.3), and use the new abundance matrix S' and new endmember matrix M' as the input of step 1.2) ;

1.5）重复执行步骤1.2）～1.4）共n次，输出最终需要求出的端元矩阵M和丰度矩阵S，计算结束，其中，n为执行次数，取值为3000。 1.5) Repeat steps 1.2) to 1.4) for a total of n times, output the endmember matrix M and the abundance matrix S that need to be obtained finally, and the calculation is completed, where n is the number of executions, and the value is 3000. the

步骤5：判断步骤1输入的高光谱图像X∈R^M×N×L中所有混合像素点的地物类别。 Step 5: Determine the object category of all mixed pixels in the hyperspectral image X ∈ R ^M×N×L input in step 1.

1.1）对步骤1输入的高光谱图像X∈R^M×N×L中的混合像素点X_ij用步骤3解混得到的端元矩阵M和丰度向量s_i表示，即X_ij=Ms_i； 1.1) For the mixed pixel point X _ij in the hyperspectral image X∈R ^M×N×L input in step 1, use the endmember matrix M and abundance vector s _i obtained by unmixing in step 3 to represent, that is, X _ij = Ms _i ;

1.2）根据高光谱图像统计分布理论，由步骤1.1）中得到的丰度向量s_i对混合像素点X_ij进行地物类别判断，即当max(s_i)=s_ai时，则判混合像素点X_ij属于第a类，得到该混合像素点的类别标签为v_ij=a，其中max(·)表示取向量中的最大值，a＝1,2,...,P表示该高光谱图像中相应的地物类别编号，P表示该高光谱图像中地物类别总数，s_ai是s_i的第a个元素； 1.2) According to the statistical distribution theory of hyperspectral images, the abundance vector s _i obtained in step 1.1) is used to judge the mixed pixel point X _ij of the ground feature category, that is, when max(s _i )=s _ai , the mixed pixel is judged Point X _ij belongs to category a, and the category label of the mixed pixel point is v _ij =a, where max( ) means taking the maximum value in the vector, and a=1,2,...,P means the hyperspectral The number of the corresponding ground object category in the image, P represents the total number of ground object categories in the hyperspectral image, s _ai is the ath element of s _i ;

1.3）对该高光谱图像X∈R^M×N×L中所有混合像素点用步骤1.2）的操作进行地物类别判断，得到该高光谱图像X∈R^M×N×L的解混之后的地物类别矩阵V∈R^M×N。 1.3) Use the operation of step 1.2) to judge the object category of all mixed pixels in the hyperspectral image X∈RM ^×N×L , and obtain the hyperspectral image X∈RM ^×N×L after unmixing Object category matrix V∈RM ^×N .

本发明的效果通过以下仿真实验进一步说明： Effect of the present invention is further illustrated by following simulation experiments:

（1）实验仿真条件： (1) Experimental simulation conditions:

本实验使用的高光谱图像是典型的AVIRIS高光谱图像：取自1992年6月拍摄的美国印第安纳州西北部印第安遥感试验区，地物类别共计16类，图像的大小为145×145。原始数据共有220个谱段，除去被噪声污染和水域污染的20个谱段，仅保留剩下的200个谱段。本实验在CPU为Intel(R)Core(TM)i5-2450、主频2.5GHz，内存为4G的WINDOWS7系统上采用软件MATLAB2009a进行仿真。 The hyperspectral image used in this experiment is a typical AVIRIS hyperspectral image: it was taken from the Indian Remote Sensing Experimental Area in Northwest Indiana, USA in June 1992. There are 16 types of ground objects in total, and the size of the image is 145×145. There are 220 spectral segments in the original data, 20 spectral segments polluted by noise and water are removed, and only the remaining 200 spectral segments are retained. This experiment uses the software MATLAB2009a to simulate on the WINDOWS7 system with Intel(R) Core(TM) i5-2450 CPU, main frequency 2.5GHz, and 4G memory. the

（2）评价指标： (2) Evaluation indicators:

1.1）对于端元，用光谱角距离SAD来比较真实端元M_t和估计端元M_t的相似性。光谱角距离SAD的值越小则两个端元的光谱曲线越接近。对于丰度，用均方根误差RMSE来比较真实丰度S_t和估计丰度S_t之间的差异。均方根误差RMSE的值越小则两个丰度值越接近。上述两个评价标准的公式分别为： 1.1) For endmembers, the spectral angular distance SAD is used to compare the similarity between the real endmember _Mt and the estimated endmember _Mt. The smaller the value of the spectral angle distance SAD is, the closer the spectral curves of the two end members are. For abundance, the root mean square error RMSE was used to compare the difference between the true abundance S _t and the estimated abundance S _t . The smaller the value of the root mean square error RMSE, the closer the two abundance values are. The formulas of the above two evaluation criteria are:

${SAD SAD}_{t t} = = {cos cos}^{- - 11} ((\frac{{M m}_{t t}^{T T} {M m}_{t t}}{| | | | {M m}_{t t}^{T T} | | | | | | | | {M m}_{t t} | | | |})),,$

${RMSE RMSE}_{t t} = = {((\frac{11}{B B} {| | {S S}_{t t} - - {S S}_{t t} | |}^{22}))}^{\frac{11}{22}},,$

其中，M_t为真实端元矩阵M的第t列，M_t为估计端元矩阵M的第t列，S_t为真实丰度矩阵S的第t列，S_t为估计丰度矩阵S的第t列，t表示地物的类别， t=1,2,...,P，P为地物类别总数； Among them, M _t is the tth column of the real endmember matrix M, M _t is the tth column of the estimated endmember matrix M, S _t is the tth column of the real abundance matrix S, and S _t is the tth column of the estimated abundance matrix S In the tth column, t represents the category of the feature, t=1,2,...,P, P is the total number of categories of the feature;

（3）实验仿真内容： (3) Experimental simulation content:

实验一 experiment one

利用本发明对步骤（1）中所述高光谱图像进行优化解混，得到端元矩阵M和丰度矩阵S，再分别用NMF算法，CNMF算法，GLNMF算法和分段光滑NMF算法对上述高光谱图像进行优化解混，与本发明解混得到的结果进行对比，实验结果如图3所示，其中: Utilize the present invention to optimize and unmix the hyperspectral image described in step (1), obtain the endmember matrix M and the abundance matrix S, and then respectively use the NMF algorithm, the CNMF algorithm, the GLNMF algorithm and the piecewise smooth NMF algorithm to process the above hyperspectral images The spectral image is optimized and unmixed, and compared with the results obtained by the unmixing of the present invention, the experimental results are as shown in Figure 3, wherein:

图3（a）是用NMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； Fig. 3(a) is a schematic diagram of the comparison between the spectral curve of a surface object chlorotopaz obtained by unmixing with the NMF algorithm and the actual surface object chlorotopaz spectral curve;

图3（b）是用CNMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； Figure 3(b) is a schematic diagram of the comparison between the spectral curve of a surface object chlorotopaz obtained by unmixing with the CNMF algorithm and the actual surface object chlorotopaz spectral curve;

图3（c）是用GLNMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； Figure 3(c) is a schematic diagram of the comparison between the spectral curve of a surface object chlorotopaz obtained by unmixing with the GLNMF algorithm and the actual surface object chlorotopaz spectral curve;

图3（d）是用分段光滑NMF算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； Figure 3(d) is a schematic diagram of the comparison between the spectral curve of a surface object chlorotopaz and the actual surface object chlorotopaz spectral curve unmixed by the piecewise smoothing NMF algorithm;

图3（e）是用本发明的算法解混出的一种地物氯黄晶光谱曲线与实际地物氯黄晶光谱曲线的对比示意图； Fig. 3 (e) is the comparison schematic diagram of a kind of surface object chlorotopaz spectral curve and actual surface object chlorotopaz spectral curve that the algorithm of the present invention unmixes;

由图3可见，本发明相比于其他现有方法，得到的地物的光谱曲线更接近真实地物地物的光谱曲线。 It can be seen from FIG. 3 that compared with other existing methods, the spectral curve of the ground object obtained by the present invention is closer to the spectral curve of the real ground object. the

实验二 Experiment 2

利用步骤（2）对实验一得到的端元矩阵M和丰度矩阵S分别计算端元矩阵M的光谱角距离SAD值和丰度矩阵S的均方根误差RMSE值，实验结果如图4所示，其中: Use step (2) to calculate the spectral angular distance SAD value of the endmember matrix M and the root mean square error RMSE value of the abundance matrix S for the endmember matrix M and abundance matrix S obtained in Experiment 1. The experimental results are shown in Figure 4 show, of which:

图4（a）是本发明与上述4种现有算法解混出的端元的光谱角距离SAD值的平均30次结果误差棒对比图； Figure 4(a) is a comparison chart of the average 30-time result error bars of the spectral angular distance SAD values of the endmembers unmixed by the present invention and the above four existing algorithms;

图4（b）是本发明与上述4种算法解混出的丰度的均方根误差RMSE值的平均30次结果误差棒对比图； Figure 4(b) is a comparison chart of the error bars of the average 30 times of the root mean square error RMSE value of the abundance unmixed by the present invention and the above four algorithms;

由图4可见，本发明相比于其他现有方法，解混得到的端元值和丰度值的精度更高。 It can be seen from FIG. 4 that compared with other existing methods, the accuracy of the endmember value and abundance value obtained by unmixing is higher in the present invention. the

实验三 Experiment three

利用实验一得到的端元矩阵M和丰度矩阵S对该高光谱图像中所有混合像素点进行类别判断，得到解混之后的类别矩阵V，再利用支撑矢量机SVM对解混之后得到的类别矩阵V和真实类别矩阵V进行地物识别精度的测量，将本发明测量所得的地物识别精度与上述4种现有方法测量所得的地物识别精度进行对比，如表1所示。 Use the endmember matrix M and abundance matrix S obtained in Experiment 1 to judge the categories of all mixed pixels in the hyperspectral image, and obtain the category matrix V after unmixing, and then use the support vector machine SVM to classify the categories obtained after unmixing The matrix V and the real category matrix V are used to measure the accuracy of feature recognition, and the accuracy of feature recognition measured by the present invention is compared with the accuracy of feature recognition measured by the above four existing methods, as shown in Table 1. the

表1 地物识别精度数值指标对比 Table 1 Comparison of Numerical Indexes of Object Recognition Accuracy

从表1可见，利用本发明进行地物识别得到的地物识别精度已经明显的优于上述4种现有算法的地物识别精度。 It can be seen from Table 1 that the accuracy of feature recognition obtained by using the present invention is significantly better than that of the above four existing algorithms. the

综上所述，本发明可以大幅度提高高光谱图像中端元值和丰度值的精度，从而更好地提高了地物识别的精度，所以本发明作为一种基于高光谱图像解混的地物识别方法在高光谱图像识别领域具有广阔的应用前景。 In summary, the present invention can greatly improve the accuracy of endmember values and abundance values in hyperspectral images, thereby better improving the accuracy of ground object recognition, so the present invention is used as a hyperspectral image unmixing based Object recognition methods have broad application prospects in the field of hyperspectral image recognition. the

Claims

1. A ground feature recognition method based on hyperspectral image unmixing, comprising the steps of:

(1) Input a hyperspectral image X ∈ R ^M×N×L , and arrange the mixed pixels X _ij ∈ R ^1×L in the hyperspectral image in columns to form a data matrix Z ∈ R ^L×B , where M and N are the rows and columns of the two-dimensional image, i and j are the abscissa and ordinate of the two-dimensional image, L is the number of spectral segments, B is the total number of mixed pixels in the hyperspectral image, B=M×N, R represents a set of real numbers;

(2) According to the manifold hypothesis theory, construct the manifold constraint item of the data matrix Z:

Figure 2013106475090100001DEST_PATH_IMAGE001

where _zi is the i-th column of Z, z _j is the j-th column of Z, and _zi is one of the k nearest neighbors of z _j , S is the abundance matrix, s _i is the i-th column of S, s _j Is the jth column of S, W is the weight matrix of Z, W _ij is an element of W,

Figure 2013106475090100001DEST_PATH_IMAGE003

is the weight of z _i and z _j , Tr( ) represents the trace of the matrix, T represents the transpose of the matrix, D is the diagonal weight matrix of Z, D _ii is an element on the diagonal of D, D _ii ＝Σ _j W _ij , Y is the manifold coefficient matrix, Y=DW;

(3) According to the hyperspectral image imaging theory, the L _1/2 norm is added to the abundance matrix S, and the sparse constraint expression ‖S‖ _1/2 is obtained as the sparse constraint item of the abundance matrix S;

(4) According to the hyperspectral image imaging theory, the Frobenius norm is added to the endmember matrix M to obtain the smooth constraint expression

Figure 2013106475090100001DEST_PATH_IMAGE004

Take as the smooth constraint term of the endmember matrix M;

(5) Add the three constraints obtained in steps (2)-(4) to the objective function of the NMF algorithm

Figure 2013106475090100001DEST_PATH_IMAGE005

, to form a new objective function:

Figure 2013106475090100001DEST_PATH_IMAGE006

Among them, α is the sparse constraint regularization parameter of the abundance matrix S, β is the smooth constraint regularization parameter of the endmember matrix M, and γ is the manifold constraint regularization parameter of the data matrix Z;

(6) Optimize and solve the objective function f'(M, S) obtained in step (5) by iterative multiplication, and obtain the endmember matrix M and abundance matrix S of the hyperspectral image X∈R ^M×N×L ;

(7) Express the mixed pixel point X _ij in the above hyperspectral image X∈R ^M×N×L by the endmember matrix M and abundance vector s _i obtained from the solution in step (6), that is, the mixed pixel point X _ij = Ms _i ;

(8) According to the statistical distribution theory of hyperspectral images, the abundance vector s _i in step (7) is used to judge the category of mixed pixels X _ij , that is, when max(s _i )=s _ai , the mixed pixel point X ij is judged to be mixed The pixel point X _ij belongs to the class a, and the category label of the mixed pixel point is v _ij =a, where max( ) means taking the maximum value in the vector, and a=1,2,...,P means the height The number of the corresponding feature category in the spectral image, P represents the total number of feature categories in the hyperspectral image, s _ai is the ath element of s _i ;

(9) Use step (8) to judge the object category for all the mixed pixels in the above hyperspectral image X∈RM ^×N×L , and obtain the object category of the hyperspectral image X∈RM ^×N×L Class matrix V ∈ R ^{M × N} .

2. The ground object recognition method based on hyperspectral image unmixing according to claim 1, wherein the manifold constraint item of constructing the data matrix Z according to the manifold hypothesis theory described in step (2) is carried out as follows :

2a) Calculate the weight W _ij of the i-th column z _i of the data matrix Z and the j-th column z _j of the data matrix Z using the thermokernel function:

Figure 2013106475090100001DEST_PATH_IMAGE007

Among them, σ is the thermonuclear parameter with a value of 1;

2b) Calculate the distance between the i-th column s _i of the abundance matrix S and the j-th column s _j of the abundance matrix S using the _L2 norm:

‖s _i -s _j ‖ ² ;

2c) According to the manifold hypothesis theory, combine the formula of step 2a) and the formula of step 2b) to form the manifold constraint term:

Figure 2013106475090100001DEST_PATH_IMAGE008

Among them, B is the total number of mixed pixels in the hyperspectral image.

3. The ground object recognition method based on hyperspectral image unmixing according to claim 1, wherein in the step (6), iterative multiplication is used to optimize and solve the objective function f'(M, S), as follows:

3a) Initialize the endmember matrix M and the abundance matrix S with random numbers between (0,1), and normalize each column of the abundance matrix S; input the endmember matrix M and the abundance matrix S;

3b) Deriving the abundance matrix S in the objective function f'(M, S)=0, and obtaining the iterative multiplication calculation formula of the abundance matrix S:

Figure 2013106475090100001DEST_PATH_IMAGE009

Among them (.*) represents the element multiplication of the matrix, (./) represents the element division of the matrix, W is the weight matrix of Z, D is the diagonal weight matrix of Z, Z is the data matrix, T represents the transformation of the matrix α is the sparse constraint regularization parameter of the abundance matrix S, and γ is the manifold constraint regularization parameter of the data matrix Z;

3c) Deriving the endmember matrix M in the objective function f'(M, S)=0, and obtaining the iterative multiplication calculation formula of the endmember matrix M:

M.*(ZS ^T -βM)./MSS ^T

Among them, β is the smooth constraint regularization parameter of the endmember matrix M;

3d) Use the new abundance matrix S' and new endmember matrix M' calculated by the formulas of step 3b) and step 3c), and use the abundance matrix S' and new endmember matrix M' as step 3a) input of;

3e) Steps 3a) to 3d) are repeated n times, and the endmember matrix M and abundance matrix S to be obtained are output, and the calculation ends, where n is the number of executions, and the value is 3000. the