IES87507Y1 - Nmf image unmixing method in wavelet domain based on adaptive local neighborhood constraints - Google Patents

Abstract

The invention provides an NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints, including: Biorthogonal wavelet packet transform is performed on the hyperspectral data matrix and the endmember data matrix respectively; the coefficient matrix at each node of the wavelet packet tree corresponding to the endmember data matrix is obtained, and the node with the minimum condition number is found; the correlation matrix data corresponding to the node with the minimum condition number is used as the input of the NMF model based on the adaptive local neighborhood weighting constraint, the abundance matrix, and endmember matrix are iteratively updated until the model converges, and the iterative updated abundance matrix and endmember matrix are output when the model converges. The wavelet packet tree nodes with the minimum condition number are used to reconstruct the endmembers.

Description

NMF IMAGE UNMIXING METHOD IN WAVELET DOMAIN BASED ON ADAPTIVE LOCAL NEIGHBORHOOD CONSTRAINTS TECHNICAL FIELD id="p-1" id="p-1" id="p-1"

[01] The invention relates to the technical field of hyperspeetral unmixing, in particular to an NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints.

BACKGROUND ART id="p-2" id="p-2" id="p-2"

[02] In recent years, hyperspeetral remote sensing imaging technology has gradually become a hot topic for scholars. The processing and analysis of hyperspeetral remote sensing images are still important bottlenecks restricting the application ofhyperspectral remote sensing, however, the research of hyperspeetral remote sensing image unmixing is a new data processing problem accompanied by the development of hyperspeetral imaging technology, and it is one of the important tasks of hyperspeetral data processing. id="p-3" id="p-3" id="p-3"

[03] Due to the limitations ofspatial resolution of hyperspeetral remote sensing image information and the influences ofcomplex ground objects, a large number ofmixed pixels appear in the image, that is, one pixel contains the spectral information ofseveral ground objects at the same time. Ifa pixel contains only one ground object, the pixel is called an endmember. Hyperspectral unmixing mainly includes endmember extraction and abundance estimation, that is, determining the types of ground objects that make up the mixed pixels and determining the proportion of each type. id="p-4" id="p-4" id="p-4"

[04] The existing methods of abundance estimation for hyperspeetral remote sensing images mainly use the fully constrained least squares method. However, when this method is used to process data with noise, the errors of abundance estimation are large and the accuracy of unmixing will be reduced. Or the non-negative matrix factorization (NMF) method is used. In order to avoid the NMF algorithm falling into the local optimal solution, scholars have proposed various constraints for NMF. Among them, the local neighborhood constraint using spatial information has achieved a good unmixing effect shortly. However, the algorithm is carried out in the time domain, but due to the ill- conditioned matrix problem, the accuracy of its unmixing is often affected.

SUMMARY id="p-5" id="p-5" id="p-5"

[05] The invention provides an NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints, which is used to solve the problem of low unmixing accuracy caused by the ill-conditioned matrix introduced by NMF constraints in the existing time domain. id="p-6" id="p-6" id="p-6"

[06] The invention provides an NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints, comprising: id="p-7" id="p-7" id="p-7"

[07] Step 1, obtaining hyperspeetral remote sensing images; id="p-8" id="p-8" id="p-8"

[08] Step 2, obtaining a hyperspeetral data matrix X, an endmember data matrix A, and an abundance data matrix S; id="p-9" id="p-9" id="p-9"

[09] Step 3, performing a d—layer biorthogonal wavelet packet transformation on the hyperspeetral data matrix and the endmember data matrix respectively to obtain a wavelet packet tree corresponding to the hyperspectral data matrix and a wavelet packet tree corresponding to the endmember data matrix respectively; id="p-10" id="p-10" id="p-10"

[10] Step 4, obtaining a coefficient matrix at each node of the wavelet packet tree corresponding to the endmember data matrix, calculating the condition number of each coefficient matrix, and finding the node with the minimum condition number which is denoted by min Cond; id="p-11" id="p-11" id="p-11"

[11] Step 5, inputting the coefficient matrix of the hyperspectral data matrix corresponding to the node min Cond of the minimum condition number, the coefficient matrix of the endmember data matrix corresponding to the node, and the abundance matrix, the three matrices are used as the inputs of NMF model based on the adaptive local neighborhood weighting constraints, iteratively updating the abundance matrix and the endmember matrix until the model converges, and outputting an iterative updated abundance matrix and endmember matrix when convergence occurs; id="p-12" id="p-12" id="p-12"

[12] Step 6, using the coefficient matrices of other nodes in the same layer as the minCond node and the iterative updated endmember matrix as the input of a reconstruction function, and outputting a reconstructed endmember matrix. id="p-13" id="p-13" id="p-13"

[13] In Step 5, the NMF model based on adaptive local neighborhood weighting constraints is: .5) = mm H x — Asufi + AZ 2 wt,- IISz - sillz [l4] i=1jEN(i) id="p-15" id="p-15" id="p-15"

[15] s.t. A > 0,52 0,1,3: 1; id="p-16" id="p-16" id="p-16"

[16] where A is an adjustable parameter to balance the approximation error and the constraint term, Wij is a contribution of the local neighborhood pixel x]- of a given pixel xi to the weight of the pixel xi, xi is the i th column of the hyperspectral data matrix, x]- is the j th column of the hyperspectral data matrix, N(i) is a column set of the local neighborhood pixels of the given pixel xi, jEN(i); s1 denotes the i th column of the abundance matrix S, Sj denotes the j th column of the abundance matrix S, P is a number of endmembers, and N is the total number ofpixels; id="p-17" id="p-17" id="p-17"

[17] in Step 5, updating the abundance matrix S and the endmember matrix A iteratively. id="p-18" id="p-18" id="p-18"

[18] A<—A - *(XST ) ./(ASST) id="p-19" id="p-19" id="p-19"

[19] s <— s -* (AJIXf + asw)./(A;Afs + AW} 5) id="p-20" id="p-20" id="p-20"

[20] where ° * denotes a multiplication of the corresponding elements of the matrix, ./ denotes a division of the corresponding elements of the matrix, the matrix X1 X2 XN W = le Xzz {N X2 XN PXN , Xi denotes a sum of the local neighborhood weights of xi, " 31' = ZieN(i)Wij and the i th column of SW denotes Si Xf and Af are the augmented matrices of X and A, respectively. id="p-21" id="p-21" id="p-21"

[21] Compared with the existing technology, the beneficial effects of the invention are as follows: id="p-22" id="p-22" id="p-22"

[22] A NMF model with adaptive local neighborhood constraints is proposed, the local neighborhood of a given pixel can be adaptively determined according to the characteristics ofthe abundance matrix, the calculated weight makes full use ofthe spatial information of the given pixel, and the pixels in the neighborhood, which makes the unmixing accuracy higher. id="p-23" id="p-23" id="p-23"

[23] Furthermore, different from the traditional time—domain NMF image unmixing method, the NlVIF model with adaptive local neighborhood constraints is combined with the biorthogonal wavelet packet transform to obtain the NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints. The biorthogonal wavelet is used as the wavelet basis to represent the data. The linear phase and short support set of the wavelet basis can compress the data and reduce the influence of noise on the unmixing results in a better way, then the ill-conditioned matrix problem in the unmixing process is overcome, and further, the unmixing accuracy is improved. id="p-24" id="p-24" id="p-24"

[24] At the same time, the calculation selects the minimum condition number node for subsequent calculation and reconstruction, which reduces the amount of data processing effectively, and the image information has been more completely retained.

BRIEF DESCRIPTION OF THE DRAWINGS id="p-25" id="p-25" id="p-25"

[25] In order to explain the technical solution of the invention or the existing technology more clearly, the drawings that need to be used in the embodiment or the existing technology descriptions are introduced briefly in the following. Obviously, the drawings in the following descriptions are some embodiments of the invention. For ordinary technicians in this field, other drawings can be obtained according to these drawings without paying creative labor. id="p-26" id="p-26" id="p-26"

[26] Fig. l is a flow diagram of the NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints in the invention; id="p-27" id="p-27" id="p-27"

[27] Fig. 2 is a local region division diagram provided in the invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS id="p-28" id="p-28" id="p-28"

[28] In order to make the purpose, technical solution, and advantages of the invention easier to understand, the technical solution of the invention is described clearly and completely in the following in combination with the attached drawings of the invention.

Obviously, the drawings in the following descriptions are some embodiments of the invention. For ordinary technicians in this field, other drawings can be obtained according to these drawings without paying creative labor belonging to the protection scope of the invention. id="p-29" id="p-29" id="p-29"

[29] The NMF unmixing method based on adaptive local neighborhood constraints of hyperspectral remote sensing images in wavelet domain is described in combination with Fig. 1. Fig. l is a flow diagram of the NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints in the invention, as shown in Fig.1. The method includes: id="p-30" id="p-30" id="p-30"

[30] Step 1, obtaining hyperspectral remote sensing images; id="p-31" id="p-31" id="p-31"

[31] specifically, obtaining a hyperspectral remote sensing image X with a size ofN X M. id="p-32" id="p-32" id="p-32"

[32] Step 2, obtaining a hyperspectral data matrix, an endmember data matrix, and an abundance data matrix; id="p-33" id="p-33" id="p-33"

[33] specifically, the VCA algorithm and the FCLS algorithm are used to initialize the hyperspectral remote sensing image to obtain the endmember matrix A_init and the abundance matrix S_init. The purpose of initialization is to limit the endmember matrix A and the abundance matrix S to a certain range before using the NMF algorithm to unmix, therefore, the convergence of the algorithm can be accelerated, the running time can be reduced, and the accuracy of the algorithm can be improved during the irritation. id="p-34" id="p-34" id="p-34"

[34] NMF is a non—negative decomposition ofa large matrix under the condition that all elements of the matrix are non—negative, and two non—negative low—rank matrices are obtained, the mathematical model is: id="p-35" id="p-35" id="p-35"

[35] V:WH (l) id="p-36" id="p-36" id="p-36"

[36] where V is the original data; w is the base matrix; h is the coefficient matrix. And , where m, n, r are constants representing the size of the matrix, V, W, and H are non—negative matrices. id="p-37" id="p-37" id="p-37"

[37] If the influence ofnoise is not considered, the linear mixed model is the same as the NMF decomposition model. The remote sensing image X denotes the matrix V in the NMF, and the endmember spectral matrix A denotes W in the NMF decomposition matrix, then the abundance matrix S can be represented by II in the NMF, and the ‘the sum is l’ constraint is also added to 11. Therefore, the endmember and abundance in the mixed pixels can be extracted simultaneously when using NMF for unmixing remote sensing images. id="p-38" id="p-38" id="p-38"

[38] X:AS+E (2) id="p-39" id="p-39" id="p-39"

[39] Step 3, performing a d-layer biorthogonal wavelet packet transform on the hyperspectral data matrix and the endmember data matrix respectively to obtain a wavelet packet tree corresponding to the hyperspectral data matrix and a wavelet packet tree corresponding to the endmember data matrix respectively. id="p-40" id="p-40" id="p-40"

[40] Hyperspectral image data is compressible in a suitable transform domain, such as Fourier, wavelet, etc. Fourier transform has good frequency localization characteristics but behaves poorly in spatial localization. In order to make a compromise between spa id="p-41" id="p-41" id="p-41"

[41] Wavelet packet analysis can provide a more detailed analysis method for signals, it divides the frequency band into multiple levels and then decomposes the high- frequency part that is not subdivided in the multi-resolution analysis, and the corresponding frequency bands can be selected adaptively according to the characteristics of the analyzed signal, so the bands match the signal spectrum, thereby improving the time-frequency resolution. id="p-42" id="p-42" id="p-42"

[42] Wavelet packet transform is developed based on wavelet transform, which overcomes the shortcomings of low frequency resolution of wavelet transform in high frequency part. A time domain signal x(t) is given, it can be written as afterj layer wavelet packet decomposition: id="p-43" id="p-43" id="p-43"

[43] x(t) = Wild") (3) id="p-44" id="p-44" id="p-44"

[44] where: id="p-45" id="p-45" id="p-45"

[45] x}(t) = 2k Cf} rpm-(t) (4) id="p-46" id="p-46" id="p-46"

[46] The wavelet packet coefficient 'Cl'k is: id="p-47" id="p-47" id="p-47"

[47] Ch = Li" "Widow (5) id="p-48" id="p-48" id="p-48"

[48] where wiykyi (t) is a wavelet packet with a scale index j, a position index k and a frequency index i. Since lbj M (t) is a set of standard orthogonal bases, when min, there id="p-50" id="p-50" id="p-50"

[50] The total energy of the signal x(t) is: id="p-51" id="p-51" id="p-51"

[51] E: = ffgxz (0dr = 2m 2311 EM" (t)x}1(r)dt (7) id="p-52" id="p-52" id="p-52"

[52] From the orthogonality of wavelet packet, there is: .j ~ id="p-53" id="p-53" id="p-53"

[53] EX 2 25:1 5} ‘8) id="p-54" id="p-54" id="p-54"

[54] where: E? = f "’[x-"(tflz dt (9) id="p-55" id="p-55" id="p-55"

[55] J —r 1 id="p-56" id="p-56" id="p-56"

[56] In the formula: Ej is the wavelet packet node energy in the i band, and the total energy of the wavelet packet node of the signal is equal to the sum of the signal node energy of each frequency band. id="p-57" id="p-57" id="p-57"

[57] The wavelet packet basis used in the biorthogonal wavelet packet transform is orthogonal, which can preserve the linear structure of hyperspectral data and move hyperspectral data to the wavelet domain, the wavelet packet transform of hyperspectral data X is an inner product, which is defined as: id="p-58" id="p-58" id="p-58"

[58] x]: = id="p-59" id="p-59" id="p-59"

[59] In the formula, X is hyperspectral data, Xw is hyperspectral data after WPT . 1b . . . transform, d-b 1s a wavelet packet base, and 0 is an Inner product operator. The transformation of the iterative rules in the wavelet domain can be obtained by formula (10): x", =A,,.s, +12". s.t. oss, $1.25, =1 (1 1) id="p-60" id="p-60" id="p-60"

[60] '4 id="p-61" id="p-61" id="p-61"

[61] In the formula, XW denotes a hyperspectral data matrix with a size of KXN, Ew denotes a Gaussian noise matrix with a size of KXN, for the e endmembers, the endmembcr matrix Aw in the wavelet domain is KXe, and Si denotes an abundance matrix with a size of cXN. Where K is the data dimension after wavelet packet decomposition, and the specific size of the value is determined by the selection of decomposition nodes and wavelet packet bases. id="p-62" id="p-62" id="p-62"

[62] The LMM model in the wavelet domain is applied to the NMF algorithm to realize the combination of wavelet transform and NMF algorithm. Generally, the square expression ofEuclidean distance is more intuitive and more widely used, the optimization problem of NMF in wavelet domain on the above basis can be expressed as: minF(AW,S) = "Xw _ AWSHE‘ 12 id="p-63" id="p-63" id="p-63"

[63] AWS ( ) id="p-64" id="p-64" id="p-64"

[64] In the formula, S is the abundance matrix. id="p-65" id="p-65" id="p-65"

[65] Specifically, in this embodiment, the d—layer biorthogonal wavelet packet transform is performed on X and A_init respectively, and the wavelet packet tree of X and the wavelet packet tree of A_init are obtained. The number of layers of packet decomposition d determines the degree of spectral band compression. When the value of d is larger, the degree of spectral band compression is higher. However, in practical applications, the spectral bands will be compressed to about 1/2 of the original (the specific value is affected by the selection ofthe wavelet basis) when there is an additional layer of the decomposition layers. Therefore, the size of the original data band should be considered when taking the value of d. lfthe value ofd is too small, the compression is not sufficient, which is easy to cause data redundancy and the ill-conditioned matrix problem cannot be solved well. If the value ofd is too large, the spectral data compression is too serious, and the original high—dimensional data cannot be effectively expressed by low—dimensional data, resulting in a decrease in the unmixing accuracy. id="p-66" id="p-66" id="p-66"

[66] Step 4, obtaining a coefficient matrix at each node of the wavelet packet tree corresponding to the endmembcr data matrix, calculating the condition number of each coefficient matrix, and finding the node with the minimum condition number which is denoted by min Cond. id="p-67" id="p-67" id="p-67"

[67] for the linear model X:AS, the condition number of matrix A is used to measure the sensitivity of matrix multiplication or inverse output to the input error, when the condition number is larger, the sensitivity becomes worse. For a general matrix A, there exists an inverse matrix A'1 ofA (i.e., A is a nonsingular matrix), then the condition number is: id="p-68" id="p-68" id="p-68"

[68] cond (A) =IIAH - llA"|| (13) id="p-69" id="p-69" id="p-69"

[69] In the formula, cond (A) denotes the condition number of matrix A, and H . denotes the norm of matrix. According to the definition, any matrix norm can be used to define the condition number. Then the condition number of the corresponding 2-norm is: id="p-70" id="p-70" id="p-70"

[70] cond, (A) =nAn2 - "Ad", (14) NAME:0 ,IIA'IHfl/O Illil X id="p-71" id="p-71" id="p-71"

[71] It has also been proved that "'i" , so the condition number of2-norm can be calculated by singular value decomposition of the matrix: id="p-72" id="p-72" id="p-72"

[72] condzm) = "‘— (15) id="p-73" id="p-73" id="p-73"

[73] In the formula, condg (A) denotes the 2-norm condition number of matrix A, (EmaX denotes the maximum singular value of the matrix, and 0mm denotes the minimum singular value of the matrix. id="p-74" id="p-74" id="p-74"

[74] There is a certain similarity between the pixels in the homogeneous region, but when considering the similarity between the pixels, it is inconvenient to judge the similarity between the central pixel and all other pixels considering the computational complexity of the algorithm. Therefore, it is necessary to solve this problem using the idea of the local neighborhood, only the pixels in the local neighborhood N(i) of xi are considered when considering the similarity ofpixels xi, and Xi denotes the i th column of the hyperspectral data matrix Y E RLXX. id="p-75" id="p-75" id="p-75"

[75] There are many ways to select the local neighborhood N(i), the common local neighborhood structures are square and circular. These structures are easy to operate, but the spatial information cannot be fully used when dealing with transition regions and boundary regions. Combined with the characteristics of abundance data, the algorithm proposes a local neighborhood determination algorithm that can make full use of spatial information. id="p-76" id="p-76" id="p-76"

[76] The above method is based on the principle that the whole hyperspectral image is divided into several different similar regions, and each similar region is close to a ground object, that is, an endmember. Fig. 2 is a partition map of a local region in hyperspectral images. Assuming that the surface features in region 1 are close to the endmember K, the proportion of endmember K in region 1 is relatively large, while the abundance values of other endmembers in region 1 are relatively small. Therefore, by comparing the abundance values of different abundance maps at the same position, the endmembers similar to the pixels at that position can be determined. Subsequently, the local neighborhood is determined by using the points with similar abundance values of the corresponding endmembers in the similar region. id="p-77" id="p-77" id="p-77"

[77] Firstly, the spatial coordinates of each pixel are determined. It is assumed that the hyperspectral image is arranged on an r> corresponds to an L—dimensional pixel vector, that is, a pixel. The N pixels of the hyperspectral data matrix Y E R are arranged into the grid in the order from left to right and from top to bottom. In this way, the relationship between the pixel yr and its spatial position (p, q) can be obtained. id="p-78" id="p-78" id="p-78"

[78] 1: (p-l)*1‘+q (l6) id="p-79" id="p-79" id="p-79"

[79] id="p-80" id="p-80" id="p-80"

[80] Secondly, the representation of abundance data is changed to obtain the similarity between pixels. Let each row of abundance matrix S be sorted into the corresponding matrix according to the spatial position of its pixels. Taking the abundance of any endmember as an example, let it be organized as , where r and e represent the number of rows and columns of the matrix, respectively. The abundances of p rows corresponding to p endmembers are organized to is expressed as , the value of the element is equal to the abundance value. id="p-81" id="p-81" id="p-81"

[81] Then, it is necessary to determine which endmember the pixel is close to before determining the local neighborhood of a given pixel with a spatial coordinate of(i, j). By comparing the abundance values corresponding to each abundance matrix with spatial coordinates of(i,j), the maximum is found, it can bejudged that the pixel is close to the endmember K. id="p-82" id="p-82" id="p-82"

[82] Finally, according to the characteristics that the abundance values corresponding to the similar regions are relatively close, in the abundance matrix corresponding to the endmember K, the positions ofall abundance values in the range of is found in the 3X3 window with as the center as the local neighborhood N (i), ’E is a controllable parameter. The local neighborhood region of each pixel can be obtained by traversing each pixel. id="p-83" id="p-83" id="p-83"

[83] In order to determine the contribution of local neighborhood pixels to a given pixel x], two factors should be considered: spatial distance and similarity of local neighborhood pixels. Because the corresponding abundance ofpixels in the local neighborhood is quite close, the similarity of pixels can be judged by the similarity of abundance. id="p-84" id="p-84" id="p-84"

[84] Based on the above analysis, the contribution of the local neighborhood pixel x]- to the weight of xi can be defined as follows: id="p-85" id="p-85" id="p-85"

[85] id="p-86" id="p-86" id="p-86"

[86] where 0t denotes the spatial distance between pixels, [3 denotes the similarity between pixels, and the value range ofB is (0,l], both of which are of the same order of magnitude. Assuming that th id="p-88" id="p-88" id="p-88"

[88] Where [ ' [ denotes the absolute value, < ' > denotes the inner product, and (1 denotes the Manhattan distance between pixel xi and pixel xi» [3 denotes the similarity between abundances s1 and s1, and s1 denotes the i th column of abundance matrix S. When the distance is relatively close and the abundance similarity is relatively high, the influence of thej th pixel on the i th pixel is greater. Obviously, Wij can measure the contribution of thej th pixel to the i th pixel in terms of spatial distance and similarity. id="p-89" id="p-89" id="p-89"

[89] For a given pixel xi, the influence of all pixels in the local neighborhood on it can be expressed as Xi: which can be expressed by calculating the sum of the contribution of (jEN(i)) each neighborhood pixel x-‘ i. id="p-90" id="p-90" id="p-90"

[90] Xizzievnwwizl’Z’M’N (21) id="p-91" id="p-91" id="p-91"

[91] Specifically, in this embodiment, (2 1'2) nodes can be obtained in the wavelet packet tree ofA, and the 2-norm condition number of each node coefficient matrix can be calculated to find the node with the minimum condition number, which is denoted as min Cond. id="p-92" id="p-92" id="p-92"

[92] Step 5, inputting the coefficient matrix of the hyperspeetral data matrix corresponding to the node min Cond of the minimum condition number, the coefficient matrix of the endmember data matrix corresponding to the node, and the abundance matrix, the three matrices are used as the inputs of NMF model based on the adaptive local neighborhood weighting constraints, iteratively updating the abundance matrix and the endmember matrix until the model converges, and outputting an iterative updated abundance matrix and endmember matrix when convergence occurs; id="p-93" id="p-93" id="p-93"

[93] There is a local minimum problem when the NMF algorithm is directly used for hyperspeetral image unmixing because of the non-convexity of the objective function of the NMF model. In order to estimate the endmemberA and abundance S in a better way, some auxiliary constraints must be added to the NMF model. According to the similarity between the abundance of a given pixel and the abundance of pixels in the local neighborhood, the prior information of the local neighborhood abundance is used as a constraint term and weighted to make full use of the spatial information of the hyperspeetral and the similarity of the local neighborhood pixels, the accuracy of unmixing is improved. id="p-94" id="p-94" id="p-94"

[94] The NMF model based on adaptive local neighborhood weighting constraint is: . . 2 .5) = mm H X — Aslli' + XXL ZjEN(i)Wij "Si — Sj"2 (22) id="p-95" id="p-95" id="p-95"

[95] s.t. A20,S>0,1',’;S=1',C id="p-96" id="p-96" id="p-96"

[96] where A is an adjustable parameter to balance the approximation error and the constraint term, Wij is the contribution ofthe local neighborhood pixel x,- ofthe given pixel xi to the weight of the pixel xi, xi denotes the i th column of the hyperspeetral data matrix, x]- denotes thej th column of the hyperspeetral data matrix, and N(i) denotes the column set of the local neighborhood pixels of the given pixel xi, j ENG); si denotes the i th column of the abundance matrix S, Si denotes thej th column ofthe abundance matrix S, P is the number of endmembers, and N is the total number of pixels; id="p-97" id="p-97" id="p-97"

[97] After derivation, the update rules of endmembers A and S are obtained: Ai—A - * (xs‘) ./ (ASST) (23) id="p-98" id="p-98" id="p-98"

[98] s <— s -* (A}'X, + 15W)./(A}'Afs + AW.* s) (24) id="p-99" id="p-99" id="p-99"

[99] where .* denotes a multiplication of the corresponding elements of the matrix, ./ denotes a division of the corresponding elements of the matrix, the matrix X1 X2 XN W = ){1 X12 {N X2 XN PXN , Xi denotes a sum of the local neighborhood weights of xi, and the i th column ofSw denotes Si ’ Si 2 EJENU) WU Si . id="p-100" id="p-100" id="p-100"

[100] Xf and Af are the augmented matrices, as formula (25): X :[ ..]lA :[ ] (25) id="p-101" id="p-101" id="p-101"

[101] f 51;, f 51;, id="p-102" id="p-102" id="p-102"

[102] In the formula, 5 is a positive constant, which is used to adjust the abundance sum to 1, when the constraint ofB is larger, the sum of each column ofS is closer to l. id="p-103" id="p-103" id="p-103"

[103] Specifically, in this embodiment, the coefficient matrix at the min Cond node in the wavelet packet tree of X, the coefficient matrix at the min Cond node in the wavelet packet tree of A_init, and the initialized abundance matrix S_init are used as the inputs of the NMF unmixing algorithm based on adaptive local neighborhood constraints, and the abundance matrix S and the endmember matrix Aw in the wavelet domain are obtained by iterative solution. id="p-104" id="p-104" id="p-104"

[104] Step 6: using the coefficient matrices of other nodes in the same layer as the min Cond node and the iteratively updated endmember matrix as the inputs of a reconstruction function, and outputting a reconstructed endmember matrix, that is, the wavelet packet tree node with the minimum condition number is used to reconstruct the endmember. id="p-105" id="p-105" id="p-105"

[105] The key points of the embodiment of the invention are as follows: 1. A NMF model with adaptive local neighborhood constraints is proposed. According to the characteristics of the abundance matrix, the local neighborhood of a given pixel can be adaptively determined, the calculated weights make full use of the spatial information of the given pixel and the pixels in the neighborhood, and the accuracy of unmixing is higher; . Secondly, different from the traditional time-domain NMF image unmixing method, the NMF model with adaptive local neighborhood constraints is combined with the biorthogonal wavelet packet transform to obtain the NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints. The biorthogonal wavelet is used as the wavelet basis to represent the data. The linear phase and short support set of the wavelet basis can compress the data and reduce the influence of noise on the unmixing results in a better way, then the ill-conditioned matrix problem in the unmixing process is overcome, and the unmixing accuracy is improved. 3. Not all nodes are selected to input the NMF model with adaptive local neighborhood constraints and participate in the reconstruction of the endmember matrix, but the minimum condition number node is selected. The application and the selection of this node reduces the amount of data processing effectively, and the image information has been preserved more completely, that is, the integrity of the image information is taken into account, and the amount of calculation is small and efficient. id="p-106" id="p-106" id="p-106"

[106] Finally, it should be explained that the above embodiment is only used to illustrate the technical solution of the invention, rather than to restrict it; although the invention is described in detail concerning the aforementioned embodiment, the general technical personnel in this field should understand that they can still modify the technical solutions recorded in the aforementioned embodiments, or replace some of the technical features with equivalents; these modifications or replacements do not make the essence of the corresponding technical solution separate from the spirit and scope of the technical solution of each embodiment of the invention.

Claims (1)

1. An NMF image unmixing method in wavelet domain based on adaptive local neighborhood constraints, comprising: Step 1, obtaining hyperspectral remote sensing images; Step 2, obtaining a hyperspectral data matrix X, an endmember data matrix A, and an abundance data matrix S; Step 3, performing a d-layer biorthogonal wavelet packet transformation on the hyperspectral data matrix and the endmember data matrix respectively to obtain a wavelet packet tree corresponding to the hyperspectral data matrix and a wavelet packet tree corresponding to the endmember data matrix respectively; Step 4, obtaining a coefficient matrix at each node of the wavelet packet tree corresponding to the endmember data matrix, calculating the condition number of each coefficient matrix, and finding the node with the minimum condition number which is denoted by min Cond; Step 5, inputting the coefficient matrix of the hyperspectral data matrix corresponding to the node min Cond of the minimum condition number, the coefficient matrix of the endmember data matrix corresponding to the node, and the abundance matrix, the three matrices being used as the inputs of the NMF model based on the adaptive local neighborhood weighting constraints, iteratively updating the abundance matrix and the endmember matrix until the model converges, and outputting an iterative updated abundance matrix and endmember matrix when convergence occurs; Step 6, using the coefficient matrices of other nodes in the same layer as the minCond node and the iteratively updated endmember matrix as the input of a reconstruction function, and outputting a reconstructed endmember matrix; in Step 5, the NMF model based on adaptive local neighborhood weighting constraints is: 101.5) = mm H x — Asufi + AZ 2 wt,- IISz - sillz i=1 jeN(i) s.t. A;0,s;0,1;s= 1; where K is an adjustable parameter to balance the approximation error and the constraint term, Wij is a contribution of the local neighborhood pixel Xi of a given pixel xi to the weight of the pixel xi, xi is the i th column of the hyperspectral data matrix, Xi is the j th column of the hyperspectral data matrix, N(i) is a column set of the local neighborhood pixels of the given pixel xi, jEN(i); si denotes the i th column of the abundance matrix S, s denotes the j th column of the abundance matrix S, P is a number of endmembers, and N is the total number of pixels; in Step 5, updating the abundance matrix S and the endmember matrix A iteratively; A<—A - *(xsT ) ./(AssT) s <— s -* (AJIXf + asw)./(A;Afs + aw_* s)