Hyperspectral remote sensing image abundance estimation method and device based on condition number weighting

Abstract

The present invention provides a method and device for estimating abundance of hyperspectral remote sensing images based on condition number weighting, the method comprising: obtaining a hyperspectral data matrix and an end metadata matrix; using biorthogonal wavelets to perform wavelet packet decomposition on the hyperspectral data matrix and the end metadata matrix respectively to obtain a corresponding wavelet packet tree; obtaining a coefficient matrix at each node of the wavelet packet tree corresponding to the end metadata matrix; calculating the condition number of each coefficient matrix; calculating the weight of each node according to each condition number; using the hyperspectral data matrix and the end metadata matrix corresponding to each node as inputs of the least squares method, multiplying the result output by each least squares method by the weight of the corresponding node, and finally summing each node to obtain an abundance estimation result. The present invention uses the condition number weighting of the coefficient matrix after wavelet packet decomposition to maximize the use of effective information in the data, thereby improving the least squares estimation and improving the unmixing accuracy of the hyperspectral remote sensing image.

Description

Technical Field

[0001] The present invention relates to the technical field of hyperspectral unmixing, and in particular to a hyperspectral remote sensing image abundance estimation method and device based on condition number weighting. Background Art

[0002] Due to the limitation of spatial resolution of hyperspectral remote sensing image information and the influence of complex objects, a large number of mixed pixels appear in the image, that is, one pixel contains the spectral information of several objects at the same time. If a pixel contains only one object, the pixel is called an endmember. Hyperspectral unmixing mainly includes endmember extraction and abundance estimation, that is, determining the types of objects that make up the mixed pixels and determining the proportion of each type.

[0003] The existing method for abundance estimation of hyperspectral remote sensing images mainly adopts the fully constrained least squares method. However, when this method processes noisy data, the error of abundance estimation is large and the accuracy of unmixing will be reduced. Summary of the invention

[0004] In order to solve the technical problems of large abundance estimation error and low unmixing accuracy of the fully constrained least squares method, the technical solution adopted by the present invention is to provide a hyperspectral remote sensing image abundance estimation method and device based on condition number weighting.

[0005] According to one aspect of the present invention, a method for estimating abundance of hyperspectral remote sensing images based on condition number weighting is provided, comprising the following steps:

[0006] Obtain hyperspectral data matrix and end metadata matrix;

[0007] Using biorthogonal wavelets to perform wavelet packet decomposition on the hyperspectral data matrix to obtain a wavelet packet tree corresponding to the hyperspectral data matrix;

[0008] Using biorthogonal wavelets to perform wavelet packet decomposition on the end metadata matrix to obtain a wavelet packet tree corresponding to the end metadata matrix;

[0009] Obtain the coefficient matrix at each node of the wavelet packet tree corresponding to the terminal metadata matrix;

[0010] calculating a condition number for each of the coefficient matrices;

[0011] Calculating the weight of each node according to each of the condition numbers;

[0012] The hyperspectral data matrix and end metadata matrix corresponding to each node are used as the input of the least squares method. The result of each least squares output is multiplied by the weight of the corresponding node to obtain the result corresponding to the node. Finally, the results corresponding to each node are summed to obtain the abundance estimation result.

[0013] Preferably, the step of performing wavelet packet decomposition on the hyperspectral data matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the hyperspectral data matrix comprises:

[0014] The size of the hyperspectral data matrix is ​​x×y×z, where x and y are image data and z is wavelength;

[0015] When the hyperspectral data matrix is ​​subjected to wavelet packet decomposition using biorthogonal wavelets, 1×1×z data are decomposed each time, and after the decomposition is completed, x×y wavelet packet trees are obtained;

[0016] The size of the hyperspectral data matrix corresponding to each node is k×1, where k is determined by the number of decomposition layers and the selected wavelet packet type.

[0017] Preferably, the step of performing wavelet packet decomposition on the end metadata matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the end metadata matrix comprises:

[0018] The end member data matrix is ​​z×e, where z is the wavelength, e is the number of end members, and each z×1 vector represents the spectral data of a ground feature;

[0019] Use biorthogonal wavelets to perform a-layer wavelet packet decomposition on e z×1 spectral data, and get e wavelet packet trees. Each wavelet packet tree has nodes;

[0020] The data of the corresponding nodes of the e wavelet packet trees are reassembled into a k×e coefficient matrix, where k is determined by the number of decomposition layers and the selected wavelet packet type.

[0021] Preferably, the step of calculating the weight of each node according to the condition number comprises:

[0022] Define the initial weight calculation function: y i =(x i -m) n ;

[0023] Where m is the maximum condition number, n is an adjustable parameter. The larger n is, the lower the weight of the matrix with the larger condition number will be. i represents the coefficient matrix of the i-th node, y i represents the initial weight of the i-th node;

[0024] According to the initial weight calculation function, the initial weight y corresponding to each node is calculated. i ;

[0025] For all y i Normalization is performed to obtain the final weight, that is, Where W iis the weight corresponding to each node.

[0026] According to another aspect of the present invention, there is also provided a hyperspectral remote sensing image abundance estimation device based on condition number weighting, comprising the following modules:

[0027] An acquisition module is used to acquire a hyperspectral data matrix and an end metadata matrix;

[0028] A wavelet packet decomposition module, used for performing wavelet packet decomposition on the hyperspectral data matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the hyperspectral data matrix;

[0029] The wavelet packet decomposition module is further used to perform wavelet packet decomposition on the end metadata matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the end metadata matrix;

[0030] The acquisition module is further used to acquire the coefficient matrix at each node of the wavelet packet tree corresponding to the terminal metadata matrix;

[0031] A calculation module, used for calculating the condition number of each coefficient matrix;

[0032] The calculation module is further used to calculate the weight of each node according to each condition number;

[0033] The least squares module is used to take the hyperspectral data matrix and end metadata matrix corresponding to each node as the input of the least squares method, multiply the result of each least squares output by the weight of the corresponding node to obtain the result corresponding to the node, and finally sum the results corresponding to each node to obtain the abundance estimation result.

[0034] The technical solution provided by the present invention has the following beneficial effects:

[0035] The present invention adopts biorthogonal wavelet transform to effectively extract features of hyperspectral data and end metadata, and ensures the linear relationship among hyperspectral data, end metadata and abundance estimation; using the condition number weighting of the coefficient matrix after wavelet packet decomposition, it can maximize the utilization of effective information in the data, thereby improving the least squares estimation and improving the unmixing accuracy of hyperspectral remote sensing images. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] The present invention will be further described below with reference to the accompanying drawings and embodiments, in which:

[0037] Figure 1 It is a flow chart of a method for estimating abundance of hyperspectral remote sensing images based on condition number weighting in an embodiment of the present invention;

[0038] Figure 2 is a schematic diagram of a wavelet packet tree having three decomposition layers in an embodiment of the present invention;

[0039] Figure 3 Schematic diagram of wavelet packet decomposition of a hyperspectral data matrix in an embodiment of the present invention;

[0040] Figure 4 is a wavelet packet tree corresponding to the terminal metadata matrix when the number of decomposition layers is four in the embodiment of the present invention;

[0041] Figure 5 is a comparison chart of the abundance estimation accuracy of the method of the present invention and the traditional method for five types of land objects in an embodiment of the present invention;

[0042] in, Figure 5 (a) is the abundance estimation result of hyperspectral remote sensing image using the traditional least squares method. Figure 5 (b) is the abundance estimation result of hyperspectral remote sensing image using the least squares method with only the minimum condition number matrix. Figure 5 (c) is the least squares hyperspectral remote sensing image abundance estimation result based on condition number weighting of the present invention, Figure 5 (d) is the real abundance map;

[0043] Figure 6 It is a structural diagram of a hyperspectral remote sensing image abundance estimation device based on condition number weighting in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to have a clearer understanding of the technical features, purposes and effects of the present invention, specific embodiments of the present invention are now described in detail with reference to the accompanying drawings.

[0045] like Figure 1 As shown, an embodiment of the present invention provides a method for estimating abundance of hyperspectral remote sensing images based on condition number weighting, comprising the following steps:

[0046] S1: Obtain hyperspectral data matrix and end-point metadata matrix;

[0047] Among them, the hyperspectral data matrix is ​​the hyperspectral image to be unmixed;

[0048] The endmember data matrix can be obtained by some classical algorithms, such as vertex component analysis (VCA) algorithm, spatial spectral endmember extraction (SSEE) algorithm, etc., or directly extracted from the endmember spectral library.

[0049] S2: using biorthogonal wavelets, performing wavelet packet decomposition on the hyperspectral data matrix and the end data matrix respectively, and obtaining corresponding wavelet packet trees respectively;

[0050] Specifically, the wavelet packet tree is Figure 2 As shown. (0,0) represents the original data, and each node represents the result of each layer of decomposition. Figure 2The figure shows a wavelet packet tree with three decomposition levels.

[0051] In step S2, the biorthogonal wavelet is first used to perform wavelet packet decomposition on the hyperspectral data matrix:

[0052] The size of the hyperspectral data matrix is ​​x×y×z, which is equivalent to a three-dimensional data. Among them, x and y are image data, and z is the wavelength, such as Figure 3 shown.

[0053] When performing wavelet packet decomposition on the hyperspectral data matrix, 1×1×z data is decomposed each time, that is, Figure 3 After the decomposition is completed, x×y wavelet packet trees can be obtained, and the size of the hyperspectral data matrix corresponding to each node is k×1, where k is determined by the number of decomposition layers and the type of selected biorthogonal wavelet. In this embodiment, biorthogonal wavelets such as bior3.3 and haar can be selected.

[0054] Then, the same wavelet packet decomposition is performed on the peer metadata matrix:

[0055] The end member data matrix is ​​z×e, where z is the wavelength, e is the number of end members, and each z×1 vector represents the spectral data of a feature (so there are e data in total). After performing a-layer transformation on the e z×1 data, a wavelet packet tree is obtained. a can be any positive integer according to the actual situation. In this embodiment, a=4 layers are taken, and 30 nodes can be obtained in each wavelet packet tree (excluding the root node, because the root node is the original data, and for the 4-layer decomposition, the number of nodes in each layer is 2 1 , 2 2 , 2 3 , 2 4 ).like Figure 4 shown.

[0056] S3: Obtain the coefficient matrix at each node of the wavelet packet tree corresponding to the end data matrix;

[0057] Specifically, the data of the corresponding nodes of the wavelet packet tree corresponding to the e-end metadata matrices are reassembled into a k×e matrix, wherein, because the data will be reduced in dimension when the wavelet packet is decomposed, the reassembled matrix is ​​no longer z×e. In this way, 30 coefficient matrices can be obtained.

[0058] S4: Calculate the condition number of each coefficient matrix;

[0059] The calculation formula of the condition number is the norm of the coefficient matrix multiplied by the norm of the inverse matrix of the coefficient matrix, that is: cond(A) = ||A||·||A -1 ||, where A represents the coefficient matrix and cond(A) represents the condition number of A.

[0060] S5: Calculate the weight of each node according to each of the condition numbers;

[0061] S5 specifically includes:

[0062] Define the initial weight calculation function: y i =(x i -m) n ;

[0063] Where m is the maximum condition number, n is an adjustable parameter. The larger n is, the lower the weight of the matrix with the larger condition number will be. You can choose according to the actual situation. i represents the coefficient matrix of the i-th node, y i represents the initial weight of the i-th node;

[0064] According to the initial weight calculation function, the initial weight y corresponding to each node is calculated. i ;

[0065] For all y i Normalization is performed to obtain the final weight, that is, Where W i is the weight corresponding to each node. In this embodiment, 4 layers of wavelet packet decomposition are performed, and the corresponding i=1,2,....30.

[0066] S6: The hyperspectral data matrix and end metadata matrix corresponding to each node are used as the input of the least squares method. The result of each least squares method output is multiplied by the weight of the corresponding node to obtain the result corresponding to the node. Finally, the results corresponding to each node are summed to obtain the abundance estimation result.

[0067] In S6, it should be noted that when the hyperspectral data matrix and the end metadata matrix are used as input, the nodes should be matched, that is, the hyperspectral data matrix is ​​the result of the node, and the end metadata matrix should also be the result of the corresponding node. In this embodiment, since there are 1×1×z data and 30 nodes, this process needs to be performed x×y×30 times. Each result obtained by the least squares estimation must be multiplied by the corresponding weight, and the final abundance estimation result can be obtained by summing up.

[0068] In order to verify the estimation accuracy of the hyperspectral remote sensing image abundance estimation method based on condition number weighting of the present invention, the synthetic data was simulated using the software MatlabR2019a. In order to construct the synthetic data, five reflectance spectra were selected from the USGS spectral library: Asphalt (gds317), Brick (gds350), Fiberglass (gds374), Sheetmetal (gds352), and Vinylplastic (gds372). An end metadata matrix of size 431×5 was formed. Using the abundance vector of size 5×1, a hyperspectral data vector of size 431×1 was generated. Considering the five different abundance maps, a hyperspectral data matrix of 431 bands, i.e., 128×128×431, was generated.

[0069] In this embodiment, for the abundance estimation of five types of land features, the least square method based on condition number weighting of the present invention is compared with the traditional least square method and the least square method using only the minimum condition number matrix. The comparison results are as follows: Figure 5 As shown, Figure 5 (a) is the abundance estimation result of hyperspectral remote sensing image using the traditional least squares method. Figure 5 (b) is the abundance estimation result of hyperspectral remote sensing image using the least squares method with only the minimum condition number matrix. Figure 5 (c) is the least squares hyperspectral remote sensing image abundance estimation result based on condition number weighting of the present invention, Figure 5 (d) is the real abundance map. By comparison, it can be seen that the method of the present invention has the highest similarity with the real abundance map, so the method of the present invention is significantly superior to the other two methods in terms of unmixing accuracy.

[0070] In addition, in order to perform quantitative analysis, in this embodiment, the root mean square error (RMSE) is used as the evaluation standard, and the specific calculation formula is: where y i is the true value, is an estimated value. The calculation results are shown in Table 1:

[0071] Table 1 Comparison of abundance estimation accuracy of the method of the present invention and other methods

[0072]

[0073] In some embodiments, a hyperspectral remote sensing image abundance estimation device based on condition number weighting is also provided. Figure 6 As shown, it includes the following modules:

[0074] Acquisition module 1, used to acquire hyperspectral data matrix and end metadata matrix;

[0075] Wavelet packet decomposition module 2, used for performing wavelet packet decomposition on the hyperspectral data matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the hyperspectral data matrix;

[0076] The wavelet packet decomposition module 2 is further used to perform wavelet packet decomposition on the end metadata matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the end metadata matrix;

[0077] The acquisition module 1 is further used to obtain the coefficient matrix at each node of the wavelet packet tree corresponding to the terminal metadata matrix;

[0078] A calculation module 3, used for calculating the condition number of each coefficient matrix;

[0079] The calculation module 3 is further used to calculate the weight of each node according to each condition number;

[0080] The least squares module 4 is used to use the hyperspectral data matrix and the end metadata matrix corresponding to each node as the input of the least squares method, multiply the result of each least squares output by the weight of the corresponding node to obtain the result corresponding to the node, and finally sum the results corresponding to each node to obtain the abundance estimation result.

[0081] The present invention adopts biorthogonal wavelet transform to effectively extract features from hyperspectral data and end metadata, and ensures the linear relationship among hyperspectral data, end metadata and abundance estimation; using the condition number weighting of the coefficient matrix after wavelet packet decomposition, it can maximize the utilization of effective information in the data, thereby improving the least squares estimation and improving the accuracy of hyperspectral data unmixing.

[0082] It should be noted that, in this article, the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article or system including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or system. In the absence of further restrictions, an element defined by the sentence "comprises a ..." does not exclude the existence of other identical elements in the process, method, article or system including the element.

[0083] The serial numbers of the embodiments of the present invention are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments. In a unit claim that lists several means, several of these means may be embodied by the same hardware item. The use of the words first, second, and third, etc. does not indicate any order and these words may be interpreted as identifiers.

[0084] The above are only preferred embodiments of the present invention, and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A hyperspectral remote sensing image abundance estimation method based on condition number weighting, characterized in that: The following steps are involved: Obtain hyperspectral data matrix and end metadata matrix; Using biorthogonal wavelets to perform wavelet packet decomposition on the hyperspectral data matrix to obtain a wavelet packet tree corresponding to the hyperspectral data matrix; Using biorthogonal wavelets to perform wavelet packet decomposition on the end metadata matrix to obtain a wavelet packet tree corresponding to the end metadata matrix; Obtain the coefficient matrix at each node of the wavelet packet tree corresponding to the terminal metadata matrix; calculating the condition number of each of the coefficient matrices; Calculating the weight of each node according to each of the condition numbers; The hyperspectral data matrix and end metadata matrix corresponding to each node are used as the input of the least squares method. The result of each least squares output is multiplied by the weight of the corresponding node to obtain the result corresponding to the node. Finally, the results corresponding to each node are summed to obtain the abundance estimation result.

2. The method for estimating abundance of hyperspectral remote sensing images based on condition number weighting according to claim 1, characterized in that: The step of using biorthogonal wavelets to perform wavelet packet decomposition on the hyperspectral data matrix to obtain a wavelet packet tree corresponding to the hyperspectral data matrix includes: The size of the hyperspectral data matrix is ​​x×y×z, where x and y are image data and z is wavelength; When the hyperspectral data matrix is ​​subjected to wavelet packet decomposition using biorthogonal wavelets, 1×1×z data are decomposed each time, and after the decomposition is completed, x×y wavelet packet trees are obtained; The size of the hyperspectral data matrix corresponding to each node is k×1, where k is determined by the number of decomposition levels and the selected biorthogonal wavelet type.

3. The method for estimating abundance of hyperspectral remote sensing images based on condition number weighting according to claim 1, characterized in that: The step of using biorthogonal wavelets to perform wavelet packet decomposition on the end metadata matrix to obtain a wavelet packet tree corresponding to the end metadata matrix includes: The end member data matrix is ​​z×e, where z is the wavelength, e is the number of end members, and each z×1 vector represents the spectral data of a ground feature; Use biorthogonal wavelets to perform a-layer wavelet packet decomposition on e z×1 spectral data, and get e wavelet packet trees. Each wavelet packet tree has nodes; The data of the corresponding nodes of the e wavelet packet trees are reassembled into a k×e coefficient matrix, where k is determined by the number of decomposition levels and the selected biorthogonal wavelet type.

4. The method for estimating abundance of hyperspectral remote sensing images based on condition number weighting according to claim 1, characterized in that: The step of calculating the weight of each node according to the condition number comprises: Define the initial weight calculation function: y i =(x i -m) n ; Among them, m is the maximum condition number, n is an adjustable parameter, the larger the n, the lower the weight of the matrix with a larger condition number, x i represents the coefficient matrix of the i-th node, y i represents the initial weight of the i-th node; According to the initial weight calculation function, the initial weight y corresponding to each node is calculated. i ; For all y i The final weight is obtained by normalization, that is, Where W i is the weight corresponding to each node.

5. A hyperspectral remote sensing image abundance estimation device based on condition number weighting, characterized in that: Includes the following modules: An acquisition module is used to acquire a hyperspectral data matrix and an end metadata matrix; A wavelet packet decomposition module, used for performing wavelet packet decomposition on the hyperspectral data matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the hyperspectral data matrix; The wavelet packet decomposition module is further used to perform wavelet packet decomposition on the end metadata matrix using biorthogonal wavelets to obtain a wavelet packet tree corresponding to the end metadata matrix; The acquisition module is further used to acquire the coefficient matrix at each node of the wavelet packet tree corresponding to the terminal metadata matrix; A calculation module, used for calculating the condition number of each coefficient matrix; The calculation module is further used to calculate the weight of each node according to each condition number; The least squares module is used to take the hyperspectral data matrix and end metadata matrix corresponding to each node as the input of the least squares method, multiply the result of each least squares output by the weight of the corresponding node to obtain the result corresponding to the node, and finally sum the results corresponding to each node to obtain the abundance estimation result.

Images