A Bayesian unmixing method for multi-temporal hyperspectral images in wavelet domain

Abstract

The present invention discloses a Bayesian unmixing method for multi-time hyperspectral images in the wavelet domain. The method utilizes the advantages of wavelet transform in signal detection, the Bayesian model, and the characteristics of wavelet coefficients of endmembers and endmember variability, and establishes prior models for high- and low-frequency wavelet coefficients of hyperspectral images respectively, so as to obtain endmembers and abundances that are closer to the true values. The core of the method is to make full use of the properties of wavelet coefficients and perform prior modeling. Based on the traditional Bayesian unmixing model, the present invention establishes prior models for its high- and low-frequency wavelet coefficients according to the characteristics of wavelet coefficients of endmembers and endmember variability. Finally, the obtained Bayesian unmixing model is solved by the MCMC method, and the obtained result is subjected to an inverse wavelet transform. Experimental results show that the proposed method has better approximation in terms of endmembers and endmember variability.

Description

Technical Field

[0001] The invention belongs to the technical field of hyperspectral image unmixing, and in particular to a multi-time hyperspectral image Bayesian unmixing method in wavelet domain. Background Art

[0002] Hyperspectral imagers can acquire data on hundreds of continuous bands and obtain hyperspectral images with high spectral resolution. They are widely used in many fields such as earth observation, military reconnaissance, mineral exploration, precision agriculture, etc. Due to the spatial resolution limitation of spectral imagers and the complexity and diversity of surface materials, the image pixels obtained by hyperspectral imaging technology often contain multiple substances. These pixels are called mixed pixels. The existence of mixed pixels will have a certain impact on hyperspectral data analysis. Therefore, it is necessary to use hyperspectral unmixing technology to separate mixed pixels and decompose them into the component (end member) spectra of the mixed pixels and their proportional coefficients (abundance).

[0003] Most of the classical unmixing methods assume that the endmember spectra are constant, ignoring the fact that the reflectance of the material may vary with factors such as light, atmosphere, and season. Therefore, it is necessary to further study the spectral unmixing problem that takes into account the endmember variability. According to the possible causes of spectral changes, spectral variation needs to consider several aspects such as spectral feature selection, spectral weighting, spectral transformation, and spectral modeling.

[0004] Roberts first proposed an unmixing method under endmember variation. By constructing spectral beams and defining standard features for each pure endmember, a limited number of spectral features can capture most of the spectral variability of each material. However, the disadvantage of spectral beams is that they rely on the diversity of the extracted clusters. A simple solution is to randomly select spectra from the spectral dataset. In this category of methods, dictionary learning methods and sparse representation methods are often used to accurately estimate endmembers and abundances.

[0005] From a physical perspective, endmember variability can also be modeled directly. An extended linear mixed model (ELMM) has been proposed, which simulates the effect of illumination on reflectance by multiplying endmembers by scaling factors. However, the ELMM model uses a fixed scale ratio for different wavelengths, which lacks the necessary flexibility when endmembers change in more complex environments. For example, the experimental measurement results of vegetation spectra are significantly different from seasonal changes. Based on the ELMM model, the generalized linear mixed model (GLMM) uses a band-based scaling factor to make the new model adaptable to arbitrary changes in endmember spectra. The augmented linear mixed model (ALMM) applies data-driven learning strategies to hyperspectral inverse problems. By introducing a spectral variability dictionary, it models other spectral variability caused by environmental conditions, instrument configuration, and material nonlinear mixing effects, while learning the spectral variability dictionary and estimating abundance. By assuming that the spectral changes are caused by additive perturbations, a perturbed linear mixing model (PLMM) is established. This model considers endmembers and their variability separately, and different constraints or regularizations can be easily imposed to simulate more complex environmental changes.

[0006] In order to solve the problem of unmixing of spectral variation, the Bayesian method incorporates meaningful prior information into the modeling process through statistical means, and can model the variability and uncertainty in spectral data, abundance and endmembers, thereby forming a more robust estimate. For endmembers and endmember variability, a commonly used model is the Normal Compositional Model (NCM), in which endmembers and endmember variability are modeled as random vectors with Gaussian distributions of different variances. However, in practice, the distribution of reflectance values ​​of materials has a certain offset phenomenon, and the Gaussian distribution cannot well express this offset of endmembers. The Beta distribution is an offset distribution. By modeling the endmembers with the Beta Composition Mode (BCN), the estimation of endmember variability is improved. However, the above endmember distributions are all assumed for each pixel, and using the same distribution for all endmembers in the hyperspectral image obviously has certain limitations. Therefore, the more complex Gaussian Mixture Model (GMM) of endmembers has achieved a good unmixing effect, but the disadvantage is that the amount of calculation is increased. In addition, using the Bayesian method, on the basis of the PLMM model, a priori models for abundance, endmembers and endmember variability were established and applied to multi-temporal hyperspectral unmixing. Although this method makes a good estimate of the endmember variability caused by temporal changes, in the modeling of abundance and endmember variation, only the temporal changes are taken into account, ignoring the spectral and spatial characteristics of the parameters.

[0007] Wavelet transform is effective in extracting endmember curve features. The low-frequency data of the endmembers obtained by wavelet transform reflects the main trend of the endmember curve, and the hyperspectral data restored from it roughly coincides with the trend of the original endmember curve. The high-frequency data reflects the changes in the endmember curve and the degree of change. The extreme value part of the high-frequency data can well reflect the information that the endmember data changes sharply in certain bands (absorption bands). Different chemical compositions and physical structures of different surface covers produce spectral curves with different spectral characteristics, while the same land object has similar spectral characteristics, and this similarity is mainly reflected in the spectral absorption characteristics. Therefore, the extraction of the spectral absorption characteristics of the endmember curve has an important influence on the unmixing results. Summary of the invention

[0008] The present invention provides a Bayesian unmixing method for multi-temporal hyperspectral images in wavelet domain.

[0009] The technical solution to achieve the purpose of the present invention is: a multi-time hyperspectral image Bayesian unmixing method in wavelet domain, which performs wavelet transform on hyperspectral data to obtain corresponding high-frequency coefficients and low-frequency coefficients, and then combines the obtained high-frequency and low-frequency coefficients into vector form, and then based on the hierarchical Bayesian method, the unmixing problem is converted into a maximum a posteriori probability solution, and then a Bayesian unmixing model is established according to the modeling method given in this paper, and the posterior distribution of the parameters is sampled and calculated by the GIBBS sampler in the MCMC method. Finally, wt(M) (q) ,wt(dM) (q) Perform the inverse transform to obtain M,dM.

[0010] Compared with the prior art, the present invention has the following significant advantages: by performing wavelet transform and modeling rules given in this article to model, the obtained end member matrix and abundance coefficient matrix can be made closer to the true value while ensuring that the peak value and extreme value of the curve are close to the true value. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 This is the flow chart of the Bayesian unmixing method for multi-temporal hyperspectral images in wavelet domain in this paper.

[0012] Figure 2 It is an image sequence of real data taken at different times in the same area.

[0013] Figure 3 It is the end members and end member variability curves obtained by different methods of real data experiments.

[0014] Figure 4 It is a graph of the abundance of the first end member (vegetation) extracted from the real data over time.

[0015] Figure 5 It is a graph of the abundance of the second end member (water) extracted from the real data over time.

[0016] Figure 6 It is a graph showing the abundance of the third end member (soil) extracted from real data over time. DETAILED DESCRIPTION

[0017] The present invention deeply analyzes the priors of abundance and endmember variation in time, space and spectral dimensions and the effectiveness of wavelet transform in endmember curve feature extraction, and proposes a hierarchical Bayesian unmixing model and algorithm in the wavelet domain for multi-time hyperspectral images. The method obtains unmixing results that are closer to the real data through wavelet transform and Bayesian unmixing method. First, the hyperspectral data is wavelet transformed, and the high and low frequency data obtained by wavelet transform are combined into vector form, and then based on the hierarchical Bayesian method, the unmixing problem is converted into a maximum a posteriori probability solution, and then a Bayesian unmixing model is established for the wavelet domain data according to the modeling method given in this article, and the MCMC method is used to solve it. Finally, the unmixing result is obtained by inverse wavelet transform. Based on the traditional hyperspectral unmixing method, this paper makes full use of the advantages of wavelet transform in signal detection, the Bayesian model, and the wavelet coefficient characteristics of endmembers and endmember variability, and establishes prior models for high and low frequency wavelet coefficients respectively, and obtains unmixing results that are closer to the real value.

[0018] The present invention is further described in detail below in conjunction with the accompanying drawings.

[0019] The present invention is a multi-time hyperspectral image Bayesian unmixing method in wavelet domain. First, the hyperspectral image data in the spatial domain is converted to the wavelet domain, and the hyperspectral data is subjected to wavelet transform, and the high and low frequency coefficients obtained by the wavelet transform are combined into a vector form. Then, according to the method given in this article, a Bayesian unmixing model is established for the wavelet domain data to obtain a linear perturbation model (Perturbed Linear Mixing Model, PLMM) in the wavelet domain. Based on the hierarchical Bayesian method, the unmixing problem is converted into a maximum posterior probability solution, and then the end members, end member changes and noise are modeled a priori according to the method given in this article. With the prior distribution of end members and end member variability, we can obtain the posterior distribution of parameters in the wavelet domain through the Bayesian formula, and then use the GIBBS sampler in the MCMC method to sample and calculate the posterior distribution of the parameters, and finally perform an inverse wavelet transform on the obtained results to obtain M, dM.

[0020] The specific steps to achieve the above are:

[0021] Step 1: Convert the hyperspectral image data in the spatial domain to the wavelet domain, perform a one-dimensional wavelet transform on the spectral dimension data of the hyperspectral image, and obtain a linear perturbation model in the wavelet domain:

[0022] wt(y n,t )=wt((M+dM t )*a n,t +b n,t )=(wt(M)+wt(dM t ))*a n,t +wt(b n,t )

[0023] Where wt(*) represents the wavelet transform operator, y n,t ∈R L×1 represents the nth pixel at time t, M = [m1,...,m R ] represents the end member matrix of size L×R, dM t =[dm 1,t ,...,dm R,t ] represents the endmember variability matrix of size L×R at time t, a n,t ∈R L×1 represents the abundance coefficient of the nth pixel at time t, b n,t ∈R L×1 Represents the error or additive noise generated in the data acquisition and modeling process.

[0024] Step 2: Convert the unmixing problem of the hierarchical Bayesian method into a maximum a posteriori probability problem, that is, model the original problem through the Bayesian formula:

[0025]

[0026] Among them, wt(Y t )=[wt(y 1,t ),...,wt(y N,t )] represents the hyperspectral data at time t, wt(M)=[wt(m1),...,wt(m R )] represents the end member matrix, wt(dM t )=[wt(dm 1,t ),...,wt(dm R,t )] represents the endmember variability matrix at time t, A t =[a 1,t ,...,a N,t ] represents the abundance coefficient of the entire hyperspectral data at time t, represents the noise variance at time t, represents the set of parameters and hyperparameters, Ψ 2 The matrix form representing the variance components of the endmember variation terms.

[0027] Step 3: A priori modeling of endmembers, endmember variations, and noise:

[0028] Step 3.1: Likelihood function p(wt(Y t )|Θ) prior modeling:

[0029]

[0030] Among them, λ represents the scale of wavelet transform, ||*|| Frepresents the Frobenius norm.

[0031] Step 3.2: Prior modeling of the endmember prior distribution p(wt(M)) in the wavelet domain:

[0032]

[0033] in, represents the high-frequency data of the rth end member in the lth band, represents the low-frequency data of the rth end member in the lth band, ξ low,high is a sufficiently large number to ensure an uninformative prior.

[0034] Step 3.3: Prior distribution of endmember variability p(wt(dM) in the wavelet domain t )|wt(M),Ψ 2 ) Modeling:

[0035]

[0036]

[0037] in Represents the high-frequency and low-frequency data vectors of the endmember variability in the wavelet domain at time t=1.

[0038]

[0039]

[0040] in Represents the high and low frequency data of end member variability in the wavelet domain, express The variance of .

[0041] It should be noted that, in the calculation process, we extract the first k modulus maximum coefficients from the posterior distribution of the high-frequency wavelet coefficients of the endmembers and the high-frequency wavelet coefficients of the endmember variations for calculation.

[0042] Step 4: Calculate the posterior distribution of the required parameters in the wavelet domain:

[0043] Step 4.1: With the prior distribution of endmembers and endmember variability, we can obtain the posterior distribution of parameters in the wavelet domain through the Bayesian formula:

[0044]

[0045]

[0046]

[0047] Step 4.2: Posterior distribution of abundance:

[0048]

[0049]

[0050]

[0051] Among them, δ() represents the characteristic function, σ 2 (a r,n,t ) is a spatial data adaptive variance defined using local differences:

[0052]

[0053] and represent the vertical and horizontal difference operators respectively, and α is used to adjust the confidence.

[0054] Step 4.3: Posterior distribution of endmember variability in wavelet domain:

[0055]

[0056]

[0057]

[0058] in It means wt(dM after removing the rth vector t ), is the lth row of wt(M), represents wt(Y t ), Yes A t The transpose of the vector formed by the rth row of A \r,t Represents the matrix A t Eliminate row r.

[0059] Step 4.4 Wavelet domain endmember posterior distribution:

[0060]

[0061]

[0062]

[0063] in It means wt(M after removing the rth vector t ), is wt(dM t )'s first row.

[0064] Step 4.5: Prior distribution of noise variance and variability variance:

[0065]

[0066]

[0067] Where IG(*) represents the inverse gamma distribution, parameter a σ , b σ , a ψ , b ψ =10 -3 is used to ensure a weakly informative prior.

[0068] Step 4.6: Calculate the posterior distribution using the Bayesian formula:

[0069]

[0070]

[0071] Where IG(*) represents the inverse gamma distribution, ||.|| F represents the Frobenius norm.

[0072] Step 5: Use the GIBBS sampler in the MCMC method to sample and calculate the posterior distribution of the parameters. The algorithm is as follows:

[0073]

[0074]

[0075] Step 6: wt(M) (q) ,wt(dM) (q) Perform the inverse transform to obtain M,dM.

[0076] Step 7: The real data used is a series of AVIRIS hyperspectral images taken in the Lake Tahoe area (California, USA) from 2014 to 2015. The scene size is 50×50 and consists of a lake and a nearby site. The scene mainly contains three types of substances: vegetation, lake, and soil. After removing these heavily polluted bands and water absorption bands, 169 of the 224 bands were used in the experiment. The method proposed in this paper (denoted as WB_SU) is compared with VCA, SISAL, OU, PLMM, HB, ELMM and MIX_HB. The relevant parameter settings used in the experiment are shown in Table 1.

[0077] Table 1 Relevant parameter settings required in the experiment

[0078]

[0079] Figure 3 The spectral curves of the three endmembers and their changes at different times are given. In general, the spectral curves obtained by these unmixing methods reflect the endmember changes to a certain extent. Among them, VCA and SISAL do not take into account the endmember variation, so the curves obtained are chaotic, while the spectra obtained by other methods have certain rules. On the second endmember (water), the endmembers and endmember variability of the PLMM and OU methods have negative values. At t = 3, the curve values ​​of HB and MIX_HB are larger than those of the curves obtained by other methods, which means that a large variability has occurred, and this result is consistent with Figure 2 The real image shown is consistent. In the image at t=3, it can be seen that the color of the water is very different from that of the images at other times. For the WB_SU method, the results obtained by this method on the second end member are somewhat different from those of other methods. The results obtained on the first and third end members are relatively stable and have better effects.

[0080] In real data (such as Figure 2 ) at t = 1, 2, and 5, the second end member (water) accounts for a large proportion of the image. At t = 3, 4, and 6, the first end member (vegetation) and the third end member (soil) are more mixed, resulting in a larger spectral change. The abundance diagram of the three end members is shown in Figure 4-6 shown. Figure 5 In 4 and 5, the water abundance estimated by SISAL and VCA methods is smaller than the actual abundance value and has a large difference with the real image. Figure 6 In the case of t = 3 and 5, the estimation error of the VCA method is large and even the end members cannot be identified. The ELMM cannot correctly identify the distribution of water, especially at t = 3 and 5 ( Figure 5 ). At the same time, it cannot effectively identify soil and vegetation, which is manifested in the high abundance of soil and the low abundance of vegetation. The OU method also has the same shortcomings in distinguishing vegetation from water, such as Figure 4-5 As shown, especially at t = 1, 2, 3. Figure 6 In the PLMM, HB and MIX_HB methods, very similar soil abundance maps were obtained. However, PLMM Figure 6 The abundance coefficient in the upper right corner cannot be accurately estimated when t=6. Compared with the HB method, the abundance result obtained by MIX_HB is smoother and more accurate. The same conclusion can be drawn from Table 2. The reconstruction errors of PLMM, HB, MIX_HB and the method proposed in this paper are relatively low. For the WB_SU method, although the secondary reconstruction error is large (as shown in Table 2), the obtained abundance image is roughly equivalent to the distribution of substances in the real image, and the effect is better.

[0081] Table 2 Numerical results obtained from real data experiments

[0082]

[0083] Therefore, the multi-temporal hyperspectral image Bayesian unmixing method in wavelet domain proposed in the present invention performs well in terms of endmembers, endmember variability and abundance estimation.

Claims

1. A Bayesian unmixing method for multi-temporal hyperspectral images in wavelet domain, characterized by: The hyperspectral image data in the spatial domain is converted to the wavelet domain to obtain the linear perturbation model in the wavelet domain. Then, based on the hierarchical Bayesian method, the unmixing problem is transformed into a maximum a posteriori probability solution. Then, the endmembers, endmember changes, and noise are modeled a priori. Finally, under the Bayesian framework, the posterior distribution of the parameters in the wavelet domain is obtained. The GIBBS sampler in the MCMC method is used to sample and calculate the posterior distribution of the parameters, and then the unmixing result is obtained through the inverse wavelet transform. The method for converting the spectral dimension data of the hyperspectral image in the spatial domain into the wavelet domain is as follows: a one-dimensional wavelet transform is performed on the spectral dimension data of the hyperspectral image to obtain a PLMM model in the wavelet domain. wt(y n,t )=wt((M+dM t )*a n,t +b n,t ) Where wt(*) represents the wavelet transform operator, y n,t ∈R L×1 represents the nth pixel at time t, M = [m1,...,m R ] represents the end member matrix of size L×R, dM t =[dm 1,t ,...,dm R,t ] represents the endmember variability matrix of size L×R at time t, a n,t ∈R L×1 represents the abundance coefficient of the nth pixel at time t, b n,t ∈R L×1 It represents the error or additive noise generated in the process of data acquisition and modeling; The unmixing problem based on the hierarchical Bayesian method is transformed into the maximum a posteriori probability problem, that is, the original problem is modeled using the Bayesian formula: Among them, wt(Y t )=[wt(y 1,t ),...,wt(y N,t )] represents the hyperspectral data at time t, wt(M)=[wt(m1),...,wt(m R )] represents the end member matrix, wt(dM t )=[wt(dm 1,t ),...,wt(dm R,t )] represents the endmember variability matrix at time t, A t =[a 1,t ,...,a N,t ] represents the abundance coefficient of the entire hyperspectral data at time t, represents the noise variance at time t, represents the set of parameters and hyperparameters, Ψ 2 The matrix form representing the variance components of the end member variation terms; The method for implementing a priori modeling of endmembers, endmember changes and noise is: making full use of the effectiveness of wavelet transform in extracting endmember curve features; The likelihood function p(wt(Y t )|Θ) prior modeling: Among them, λ represents the scale of wavelet transform, ||*|| F represents the Frobenius norm; The end member prior distribution p(wt(M)) in the wavelet domain is modeled as: in, represents the high-frequency data of the rth end member in the lth band, represents the low-frequency data of the rth end member in the lth band, ξ low,high is a sufficiently large number to guarantee an uninformative prior; The prior distribution of endmember variability p(wt(dM t )|wt(M),Ψ 2 ) Modeling: in Represents the high-frequency and low-frequency data vectors of the end-member variability in the wavelet domain at time t=1; in Represents the high and low frequency data of end member variability in the wavelet domain, express The variance of In the process of calculation, the first k modulus maximum coefficients are extracted from the posterior distribution of the high-frequency wavelet coefficients of the end member and the high-frequency wavelet coefficients of the end member variation for calculation; In the Bayesian framework, the posterior distribution of the parameters in the wavelet domain is obtained as follows: The posterior distribution of abundance: Among them, δ() represents the characteristic function, σ 2 (a r,n,t ) is a spatial data adaptive variance defined using local differences: and They represent the vertical and horizontal difference operators respectively, and α is used to adjust the confidence; The posterior distribution of the endmember variability in the wavelet domain is: in It means wt(dM after removing the rth vector t ), is the lth row of wt(M), represents wt(Y t ), Yes A t The transpose of the vector formed by the rth row of A \r,t Represents the matrix A t Eliminate row r; The posterior distribution of endmembers in the wavelet domain is: in It means wt(M after removing the rth vector t ), is wt(dM t )'s first row; Prior distributions of the noise variance and variability variance: Where IG(*) represents the inverse gamma distribution, parameter a σ , b σ , a ψ , b ψ =10 -3 is used to ensure a weakly informative prior; The posterior distribution is calculated by the Bayesian formula: Where IG(*) represents the inverse gamma distribution, ||.|| F represents the Frobenius norm.

Images