Bayesian Subpixel Mapping of Hyperspectral Imagery via Discrete Endmember Variability Mixture Model and Markov Random Field

Although Bayesian methods have been very effective for spatial–spectral analysis of hyperspectral imagery (HSI), they had not been fully explored for enhanced subpixel mapping (SPM) by simultaneously considering several key issues, i.e., endmember variability, the discrete nature of subpixel class labels, and the spatial information in HSI. Therefore, we propose a new Bayesian SPM method based on the discrete endmember variability mixture model (DEMM) and Markov random field (MRF), which has three main characteristics. First, DEMM allows the advanced SPM by completely accounting for the endmember–abundance patterns of each pixel to accommodate the endmember variability, the discrete hidden class label field of subpixels, while taking into account the noise heterogeneity effect. Second, the discrete class label fields modeled by MRF together with the DEMM, which can be integrated into a novel Bayesian model to better exploit the spatial contextual and spectral information. Third, the resulting Bayesian model can be efficiently solved by a designed expectation–maximization iteration, where E-step estimates the subpixel class label field using a simulated annealing algorithm and M-step estimates the endmembers for each pixel in HSI using the alternating non-negative least squares approach. The experimental results on three HSI datasets demonstrate that the proposed approach outperforms previously available SPM techniques.


I. INTRODUCTION
H YPERSPECTRAL imagery (HSI) has been utilized in various Earth observation applications, thanks to its rich spectral information. However, owing to the influence of the instantaneous field of view of the remote sensor as well as the complexity and diversity of the land cover types, the measured spectrum of each pixel in HSI is always a mixture of the reflectance of multiple distinct materials [1]. To address the challenge of mixed pixels in HSI, numerous techniques have been developed, such as soft classification [2] and spectral unmixing [3]. Nevertheless, these approaches fail to locate the subpixels, and therefore cannot estimate the class label of subpixels for further applications [4], [5]. To deal with this issue, the subpixel mapping (SPM) [6] method, which can classify HSI at the subpixel level and obtain the optimal subpixel spatial distribution according to spatial information and abundance fractions, is typically considered to be an alternative solution.
There are many SPM algorithms that have been developed and extensively used in a number of remote sensing applications [7], [8], [9], [10], [11], including the pixel/subpixel spatial attraction model (SAM) [12], the pixel-swapping approach (PSA) [5], the Hopfield neural network [13], the Kriging and indicator co-Kriging [14], [15], [16], the spatial-spectral interpolation [17], [18], the class boundaries-based subpixel model [19], and the super-resolution convolutional neural network [20]. These algorithms usually contain the following two parts [18]. First, the soft attribute values (between 0 and 1) of each subpixel is estimated to generate initial class label of the subpixels. Second, the initial subpixel class labels are updated according to the land cover class fractions of mixed pixels and the soft attribute values of subpixels. In general, the SPM models can be solved by some optimization program, such as the genetic algorithms [21], [22], particle swarm optimization [23], and simulating annealing [24], [25]. For these algorithms, the estimate of abundance (i.e., land cover class fractions of mixed pixels) and endmember is predefined and not allowed to be modified in second step, and therefore the accuracy of SPM depends strongly on the accuracy of predefined abundances and endmembers.
Several SPM algorithms are designed to allow better estimates of the abundance and endmember, such as the Markov random field (MRF)-based SPM algorithm [26], the spectral-spatial integration SPM model [27], the adaptive MAP-based class This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ determination SPM [28], and the genetic algorithm-based approach [29]. The objective function of these methods mainly involves the likelihood term and the class label field prior (i.e., spatial regularization terms) [30]. The class label prior term makes it more likely that spatially close subpixels belong to the same land cover class membership, and the likelihood term assures that the number of subpixels within a coarse pixel coincides with the predefined abundances of land cover class.
However, the endmember variability issue is not investigated in the aforementioned methods, which can therefore reduce the accuracy of endmembers and their corresponding abundances, thus leading to unsatisfactory SPM performances. Endmember variability refers to a variation of the measured spectral signature for a single material in HSI due to the variable illumination, atmospheric, and temporal conditions [31], [32], [33]. As a consequence, it is necessary to explore the impact of endmember variability when performing the SPM method. For example, the joint sparse SPM model takes into account the endmember variability in SPM [34]. Nevertheless, the endmember variability information in this algorithm requires to being manual-defined as prior knowledge.
In this article, a novel Bayesian SPM for HSI is designed, which integrates discrete endmember variability mixture model (DEMM) and MRF into the Bayesian framework for improved SPM result. The main contributions and innovations of this study are introduced as follows.
1) The novel DEMM is proposed to enhance the performance of SPM by fully taking into account the patterns of endmember-abundance in HSI for accommodating the endmember variability and the discrete hidden subpixel class label information in each coarse pixel while considering the heterogeneous noise effect. Instead of estimating a group of endmembers for all given pixels, the proposed method estimates a group of endmembers for each pixel. In such a case, by considering a group of endmembers within a pixel using the linear mixed model (LMM), DEMM is capable of accommodating endmember variability when performing SPM. 2) A discrete hidden class label field modeled by MRF together with the DEMM, which can be seamlessly integrated into the newly developed Bayesian model to effectively utilize the spatial contextual and spectral reflectance information contained in HSI. By accommodating the discrete class label information of the subpixels, DEMM attempts to link SPM with LMM and allows the use of MRF to exploit the spatial information in the discrete class label field.
3) The resulting Bayesian model can be efficiently solved by a designed expectation-maximization (EM) iteration, which treats the subpixel class label field as missing observations and iteratively alternates the estimation of subpixel class labels based on the endmembers and the update of the endmembers given the subpixel class labels. In this EM optimization, the E-step intends to estimate the subpixel class label field using a simulated annealing (SA) algorithm, and M-step is able to estimate the endmembers for each pixel in HSI using the alternating non-negative least squares (ANLS) approach. The rest of this article is organized as follows. The proposed Bayesian framework for SPM taking into account the endmember variability is described in Section II. Section III provides the model optimization of the proposed method in detail. Section IV conducts the experimental analysis on both synthetic and real HSI imageries. Finally, Section V concludes this article.

A. Discrete Endmember Variability Mixture Model
This article presents a novel generative model, named DEMM, which describes how the observed hyperspectral pixels are associated with the discrete hidden class labels of subpixels, the endmember variability of each coarse pixel, and the heterogeneous noise effect across spectral channels in the HSI.
Let I (for I = 1, 2, . . . , N) denote the positions of coarse pixels within the hyperspectral image. Each coarse pixel x i is encoded in a P × 1-dimensional vector, where P is the number of spectral channels. In the following, a collection of positions of the subpixel in coarse pixel is expressed as J. Let the scale factor of the SPM be d, then including d 2 subpixels within each pixel. The l i,j represents the class label of the subpixel at jth spatial position in the ith coarse pixel. For the entire HSI, the image is described as X = {x i |i ∈ I}, and the hidden subpixel class labels field can be described as L = {l i,j |i ∈ I, j ∈ J}. We assume that the HSI shares K common endmembers, then s k i (for k = 1, 2, . . . , K) denotes the abundance coefficient of kth land cover class at the ith coarse pixel.
Traditionally, the LMM does not account for the endmember variability issue and assumes that different pixels in the image have the same set of endmembers [35]. According to the LMM, a given spectrum of the ith pixel, i.e., x i , can be obtained as a linear combination of the group of endmembers {a k } for the whole image weighted by the abundance coefficients s k i , plus the Gaussian noise n.
where a k is the same for all pixels in X, and therefore (1) is not expected to effectively tackle the problem of endmember variability. However, due to the environmental, atmospheric, and temporal factors, each coarse mixed pixel x i (for i = 1, 2, . . . , N) has its own specific group of endmembers a k i ; they fluctuate within a range in HSI [33]. By considering that a k i varies for different x i (endmember variability), the endmember can be better captured and more accurately estimated. Therefore, estimating a k i for each pixel x i is supposed to facilitate the SPM implementation.
Here, our newly proposed DEMM model takes into consideration both the discrete nature of the subpixel class labels and noise variance heterogeneity. The observed pixel spectrum x i in DEMM can be formulated as a linear combination of K endmembers with spectral variability {a k i } weighted by the discrete version of the fractional abundance s k i , plus the band-dependent noise term n resulting from the varying physical properties.
where an additive noise n is assumed to be Gaussian-distributed, and the probability density can thus be formulated as where Λ is described as a noise covariance matrix. The proposed DEMM model defined by (2) and (3) has the following characteristics. 1) Endmember Variability Mixture Model: Comparing with a k in (1), the a k i in (2) is specific to the ith pixel, and therefore it allows each pixel to have its own group of endmembers. Therefore, in (2), instead of estimating a group of endmembers A = {a k |k ∈ K} for the whole image X, it estimates a group of endmembers A i = {a k i |k ∈ K, i ∈ I} for a single pixel x i . Although ignoring the endmember variability in HSI may compromise the SPM performance, most existing SPM methods cannot well cope with this problem. By estimating the group of endmembers A i for each pixel, DEMM is expected to achieve favorable SPM result by addressing endmember variability. 2) Discrete Hidden Class Label Field: According to the generative model defined by (2) and (3), the number of subpixels belonging to the kth land cover class within the ith coarse spectral pixel can be given by s k i as follows: where δ(·) denotes the Kronecker delta function, with δ(a, b) is equal to 1 when a = b, otherwise, δ(a, b) can be set as 0. From (4), s k i is generally treated as a discrete version of the fractional contribution of the kth land cover class a k i , and it also reflects the discrete nature of the subpixel class labels. In this sense, the discrete class label field allows the use of some spatial models, such as MRF, to capture the spatial information in HSI.
3) Spatial Correlation Effect: Through an MRF prior distribution, the hidden class label field of subpixels can be seamlessly integrated into the DEMM model. It enables DEMM to use spatial contextual information of subpixel class labels, i.e., adjacent subpixels in HSI are supposed to have the same class label. 4) Band-Weighted Noise Model: Due to the model and measurement errors, the different physical properties of the varying spectral channels as well as the presence of junk bands, each spectral band tends to have different noise levels in the HSI [36]. Therefore, it is important to consider the noise heterogeneity effect for better suppress the impact of noise when performing SPM. It is worth noting that the noise distribution is characterized by Λ, which can be expressed as the following diagonal matrix [37]: where σ 2 j (for j = 1, 2, . . . , P ) in Λ represents the noise variance of the jth band. Using such a band-dependent noise model has great capability to handle the heterogeneous noise problem in the spectral domain for enhanced SPM results. Consequently, the novel developed DEMM model can improve SPM results by better characterize endmember variability, discrete nature of the subpixel class labels, and spatial correlation effect.

B. Bayesian SPM Formulation
Under SPM, we intend to achieve the subpixel class labels l i,j of the jth subpixel within the ith coarse pixel given the observations {x i }. From the point of view of the Bayesian theorem, the subpixel class labels can be estimated by maximizing the posteriori probability of {l i,j } in the presence of {x i } as follows: where p({x i }|{l i,j }) denotes the data likelihood term that describes the probability density of observed spectrum {x i } conditioned on {l i,j }, which indicates the modeling of the error in the observation model defined in (2), and p({l i,j }) represents the a priori probability of subpixel class labels, which models the spatial contextual information for the estimation of subpixel class labels.

C. Data Likelihood
According to the proposed DEMM, the data likelihood term p({x i }|{l i,j }) in (6) can be reformulated as where z represents a scaling parameter and · 2 2 is the Euclidean measure.
Note that, the best assignment of the {l i,j }, that is, supposed to produce the shortest Euclidean distance measure between x i and the endmember-abundance reconstruction, will result in the highest probabilistic value of p({x i }|{l i,j }).

D. MRF-Based Class Label Prior
Assuming that the discrete class label field has MRF property, the land cover class occupying adjacent subpixels has a greater chance of belonging to the same land cover class membership [26]. In this sense, in order to encourage the fact that the subpixel l i,j and its neighborhood tend to have the same class label, the MRF-based multilevel logistic (MLL) prior can be adopted to model the spatial context information, which can be expressed as where v denotes a constant for normalization and N i,j represents the neighborhood centered at the subpixel l i,j . p({l i,j }) encourages the neighboring subpixels to belong to the same class labels.
In this way, resort to MLL in (8) that promotes the subpixel class label homogeneity of the local neighborhood, since p( demonstrates the fact that adjacent subpixels tend to have the same class labels.

A. MAP Model Estimation
The MAP framework has been proven to be an effective way to efficiently address the ill-posed SPM problem. It can achieve an optimal discrete class label field of the subpixels {l i,j } by maximizing the posteriori probability of L on the basis of the observed HSI imagery X as follows: where Θ = {Λ, {a k i }, L} is the model parameter, which includes the group of endmembers for each pixel a k i , the noise covariance matrix Λ, and the subpixel class label field L. The equation p(L|X, Θ) can be split by the multiplication of the data likelihood function in (7) and the MRF-based MLL class label prior in (8), and is simply denoted as Maximizing the MAP objective of p(L|X, Θ) is equivalent to taking the negative symbol on its minimum logarithm likelihood function, i.e., Then, the objective function can also be simplified by eliminating the normalizing coefficient of class label prior and adding an extra coefficient η as follows: where η and ω are the two important weighing coefficients, which can be used for determining the relative weights between the class label prior and the likelihood term.
As we can see, unknown parameters in Θ include the set of endelements a k i of each pixel, the noise distribution matrix Λ, and the class label field L of subpixels. Since the number of unknown quantities in Θ of the equation is bigger than the number of observations, and as such, offering a plausible estimation of Θ can be seen as an ill-posed issue, which is typically addressed by the EM optimization strategy.

B. EM Optimization
The EM optimization is an extensively utilized technique for dealing with the incomplete data issue by treating some potential parameters as missing observations and iteratively alternates between two-stage optimization process (i.e., the estimation of the model parameters based on missing observations and the updates of missing observations given the model parameters) until convergence to yield optimal estimates [38]. Accordingly, the proposed algorithm treats the class label field {l i,j } as the missing observations and Λ and {a k i } as the model quantities to derive iterative estimation of all unknown parameters. Then, we finally provide the main steps of estimating {l i,j }, Λ, and {a k i } using the designed EM approach as follows.
where var(·) denotes the bandwise variance function, which calculates the variances of the reconstruction error for measured pixels.

C. {a k i } Estimation in M-step
To simplify the estimation, we decompose {a k i } into some reference endmembers a k 0 , they are pixel-independent, and some pixel-dependent scale factors ψ k i can be adopted for endmember variability [33] as follows: Algorithm 1: Finding a Local Minimum of (15) Using the ANLS Approach. Input: spectral stack X, endmembers A 0 , and iteration numbers Kiter According to (14), a criterion (minimization of energy) can be designed to perform spectral unmixing using (2) and to achieve the a k i as where Ψ i ∈ R K×K represents a diagonal matrix with elements ψ k i ≥ 0, A 0 denotes the reference endmember matrix that can be predefined, whose columns are the a k 0 , A i ∈ R P ×K is the collection of pixel-dependent endmember matrices with values a k i , and S i ∈ R K×1 are the abundance vectors for pixel x i . The parameter λ s determines the relative weight of two terms in (15). In this article, we use the ANLS method to estimate the (15), which is summarized in Algorithm 1. Note that, here λ s is defined as 0.5, the iteration numbers Kiter is set to 100, the initial abundance maps used are those of scaled version of non-negative constrained least squares unmixing (SCLSU) [39] and the initial endmember variability scale factors will be set to 1.

D. {l i,j } Estimation in E-step
The SA technique [40] has been developed to estimate the class label {l i,j }, which is called the state of the system. Moreover, the new state, i.e., L new , in the designed SA algorithm can be achieved by randomly choosing and updating the subpixel class label fields, and whether the new state of L new is accepted depends on the acceptance probability prob(e new , e, T ) defined in (16). It should be noted that the optimal state L best with the lowest energy over all iterative steps will be obtained as the final output when a predefined maximum number of iteration (i.e., M max) is reached.
where e new and e represent the energies of L new and L, respectively, and the parameter T denotes the current system temperature. Note that here, if a value randomly produced by a uniform distribution within the range of [0,1] is not larger than the prob(e new , e, T ) in the SA approach, then the new state of L new is recorded. Otherwise, the current system still favors L. Moreover, the possibility of acceptance is expected to increase with the cooling of system temperature T according to i }) end while the following: where m and m − 1 are the current and the last iteration, respectively. More specific, T 0 is assigned to 3. The SA algorithm has been detailed in our previous article [41]; it is worth noting that we use a distinct endmember set for each pixel in the designed SA algorithm to account for the endmember variability in this article, rather than a single endmember set for whole HSI imagery discussed in [41].

E. Complete Algorithm
The newly developed SPM framework based on DEMM and MRF is summarized in Algorithm 2.

IV. EXPERIMENT AND ANALYSIS
The proposed algorithm is evaluated on three HSIs, in comparison with several popular SPM methods, including the subpixel/pixel SAM [12], PSA [5], the spectral and spatial integration SPM model (SPMLM) [27],the MRF-based SPM model (MRFSPM) [26], as well as the genetic algorithm-based approach (GAAI) [29]. For all experimental data, we use the SCLSU to obtain the initial abundance maps for all SPM methods in this article. To determine the best relative weight of η value in this approach, we perform the hyperparameter searching with interval of 0.1 within the range of 0.1-0.9. In addition, we present a detailed quantitative evaluation using the kappa coefficient and the overall accuracy (OA).

A. Experiment 1: Simulated HSI
The simulated HSI is generated by using the (2) with some Gaussian noise. Six reference endmembers, which has 188 spectral bands ranging from visible to near-infrared, are randomly picked from the U.S. geological survey digital spectral library, as shown in Fig. 1(a). The simulated class label for the subpixel with size 144 × 144 is designed, as shown in   (14)] are generated. Each map is associated with one endmember, because the pixels in the map have different values, and therefore different pixels tend to have different endmembers. Fig. 1(b). To simulated endmember variability, we generated endmember variability scale factors maps for each endmember, using mixtures of Gaussian with the range from 0.75 to 1.25, as shown in Fig. 2.
In addition, the pixel-dependent endmember instances can be acquired through multiplying the selected reference endmembers by their corresponding endmember variability scale factors, then the pixel-dependent endmember combined with the class label field map are used to reconstruct a fine-resolution image. Moreover, by averaging the original image with higher resolution, the coarse resolution imagery with downsampling scale coefficients of 2, 3, and 4 can be achieved, respectively. In this sense, the reconstruction scale factor is set to 2, 3, and 4, respectively. Consequently, the simulated HSI data are degraded using the zero-mean Gaussian. To simulate the noise heterogeneity effect, we assign different SNR values to different bands. The band-dependent SNR values used for simulation are estimated from the benchmark Indian Pines image. Suppose that the estimated SNR vector q has been centralized and normalized; then the simulated SNR r can be obtained according to the following rule: where α is the amplitude that determines the magnitude of fluctuation of band-dependent SNR and c is the center value that defines the overall SNR of all bands. In this article, the SNR r is calculated with α = 9 and c = 25.
In this experiment, the reference endmembers in each coarse pixel directly use the selected endmembers for data simulation. Fig. 3 shows the SPM outcomes obtained using distinct approaches over various scale factors. As can be seen, the proposed SPM method outperforms the other referenced algorithms across all downsampling scales and produces more consistent SPM results that closely match the ground truth map shown in Fig. 1(b). We can see that there exist numerous noise artifacts in the SAM and PSA results. Although the GAAI, MRFSPM, and SPMLM methods achieve smoother result, some details are also lost. The proposed method not only can resist the noise influence, but also well preserve the details (see the details within the red rectangles in Fig. 3), since it takes full account of endmember variability and band-dependent noise. Table I presents the quantitative results of the proposed algorithm and the baseline methods, which indicate the consistent results with the abovementioned visual comparison. Overall, it is obvious that the proposed algorithm achieves the best performance among all the SPM methods, with the highest numerical evaluation indexes of OA and kappa coefficients over different scale factors. Specifically, compared with the GAAI, the accuracy of the proposed algorithm increased dramatically, with the OA rising by 2.95%, 8.26%, and 7.93% for scale factors of 2, 3, and 4, respectively. The OA values of the proposed algorithm are 1.19%, 3.24%, and 4.61% higher than those of SPMLM for scale factors of 2, 3, and 4, respectively. The OA values of the proposed algorithm are 1.37%, 2.35%, and 2.64% higher than those of MRFSPM for scale factors of 2, 3, and 4, respectively.

B. Experiment 2: Jasper Ridge Imagery
Jasper Ridge imagery with a size of 512 × 614 pixels was acquired by the AVIRIS sensor in central California, USA. Each pixel contains 224 spectral channels. Due to the original Jasper Ridge datasets are too complicated to collect the ground truth map, we achieve a subimage composed of 100 × 100 pixels, sampled on 198 wavelengths by removing the bands severely affected by the dense water vapor and ever-changing atmospheric effects. The 2 × 2 and 4 × 4 mean filters are applied to the fine imagery to obtain coarse resolution imagery with the size of 25 × 25 and 50 × 50. Remove one row and one column to get an imagery with the size of 99 × 99 pixels. The 3 × 3 mean average is applied to this imagery to obtain coarse resolution imagery with the size of 33 × 33. Therefore, the reconstruction scales are assigned to 2, 3, and 4. The fine and coarse resolution imagery with the scale factor of 2, 3, and 4, are depicted in Fig. 4(a)-(d), respectively. The study region contains four land cover classes of interest: 1) road; 2) soil; 3) water; and 4) tree. Note that, the reference endmembers of the four land covers used in the experiment came from the previous publication in [42], as shown in Fig. 4(f). The ground truth is also produced by [42], as shown in Fig. 4(e).     5 shows the SPM results obtained by different approaches on the Jasper Ridge imagery with the scale factor of 2, 3, and 4; it also indicates that the proposed approach can better estimate the subpixel class label field with rich detail information with less noise influence (see the details within red rectangles in Fig. 5). It can be seen that SAM and PSA contain the noise artifacts in the results. Although GAAI, SPMLM, and MRFSPM can better resist the noise impact, they cannot recover detail very well, due to the lack of consideration of endmember variability and heterogeneous noise effect. Table II gives that the proposed method for scale factor 2, 3, and 4 achieve OA of 88.29%, 86.54% and 84.90%, respectively. The OA values achieved by the proposed algorithm are 2.77%, 4.71%, and 7.21% higher than that of GAAI for scale factors of 2, 3, and 4, respectively. The OA values of the proposed algorithm are 1.00%, 1.20%, and 1.46% higher than those of SPMLM for scale factors of 2, 3, and 4, respectively. Compared with the proposed algorithm, the value of OA decreases by 0.88%, 0.84%, and 1.06% when it comes to MRFSPM for scale factors of 2, 3, and 4, respectively.  Fig. 6(a) shows the Urban dataset (300 × 300 pixels), which was recorded by the hyperspectral digital imagery collection experiment in 1995 over an urban area at Copperas Cove, TX, U.S [42]. The Urban hyperspectral data contain six land covers of: 1) asphalt; 2) grass; 3) tree; 4) roof; 5) dirt; and 6) metal, and its ground truth map are generated by [42], as shown in Fig. 6(e). Each original pixel is observed at 210 bands. There are 162 spectral reflectance bands are collected in this study after excluding the channels that are severely influenced by the water vapor and the atmospheric conditions. The coarse resolution imagery of 150 × 150, 100 × 100, and 75 × 75 can be obtained by applying the 2 × 2, 3 × 3, and 4 × 4 mean filters operation, respectively, as shown in Fig. 6(b)-(d), then the reconstruction scale is set to 2, 3, and 4, respectively. It is worth noting that the reference endmembers of the six land cover classes used in the experiment came from the previous publication in [42], as can be seen in Fig. 6(f). Fig. 7 shows the SPM results produced by different approaches on the Urban dataset imagery with different scale factors d = 2, 3, 4. Overall, it shows the consistent results with the experiments 1 and 2, which reveal that the proposed method can yield higher accuracy outcomes than that of the other five popular SPM algorithms. Specifically, the GAAI, SPMLM, and MRFSPM can eliminate certain unmixing errors and generate smoother subpixel spatial distribution than the PSA and SAM, especially on the homogeneous regions. Moreover, we also noticed that it is difficult for GAAI, SPMLM, and MRFSPM to reconstruct some subpixel distribution details information compared with the proposed method. For example, the proposed method is able to recover more details of the dirt at different scales within the red rectangles in Fig. 7, which is more consistent with the ground truth label map displayed in Fig. 6(e). Table III lists the OA and kappa indexes of the different SPM algorithms for quantitative evaluation; it indicates the consistent results with Fig. 7. As can be seen from Table III, the proposed method provides the greatest performance among the methods, the accuracies of the proposed method for scale factors of 2, 3, and 4 are equal to 82.38%, 79.05%, and 76.19%, respectively, which is an improvement in the SPM accuracies of about 2.80%, 3.73%, and 3.93% when compared with the SPMLM, 2.30%, 7.76%, and 6.54% when compared with the GAAI, and 2.31%, 1.53%, and 0.83% when compared with the MRFSPM. With the increase of scale factors, the accuracy of all methods becomes lower, because the spatial distribution of mixed pixels in large scale is more complicated. These experimental results demonstrate that the proposed algorithm outperforms the other popular methods for the Urban dataset.

C. Experiment 3: Urban Dataset Imagery
D. Sensitivity Analysis of η η is a very important parameter that controls the contribution weights of the likelihood and the class label prior (the spatial regularization term). The large value of η, indicating that the spatial regularization term plays a dominant role in the objective function, will result in the oversmooth of the results of SPM. In addition, the small value of η, indicating that the likelihood occupies the dominant role in the objective function, will lead to the SPM results affected by the unmixing error and the existence of noise. To determine the optimal value of η in the proposed approach, we perform the hyperparameter searching ranging from 0 to 0.9 with a step of 0.1 in the abovementioned experiments. Meanwhile, to investigate the robustness of the algorithm under different η values, we repeated the experiment of each η value for ten times and calculated its mean and variance values. The statistical summary results of each experiment are shown in Fig. 8. Fig. 8 shows the sensitivity analysis of the η value using three different datasets; it indicates that with the increase of η value, the accuracy of the algorithm increases and then gradually decreases, which is in-line with our expectation that the η value somewhere between 0 and 0.9 finds the best balance point between the spatial regularization term and the likelihood. Moreover, we can also conclude that the accuracy fluctuation of the algorithm is not obvious in the experiment of selecting the value of η, indicating that our method exhibits a satisfactory robustness.

E. Further Discussion
In this article, we propose a Bayesian SPM of HSI via DEMM and MRF, which establishes a Bayesian formula based on the DEMM and MRF and adopts iterative optimization by the EM algorithm. The E-step estimates the discrete subpixel class label modeled by MRF, and M-step estimates the endmembers for each pixel by ANLS. We compare with other methods in the three datasets. The advancement and effectiveness of the proposed method are verified from the SPM results, the local amplification map, and the index table. We assume that the better performance   of the proposed approach over the other approaches is mainly due to the following. 1) First, the designed DEMM model enables the proposed SPM to consider endmember variability, discrete hidden label fields of subpixels, and noise heterogeneity. 2) Second, the DEMM and the discrete label fields modeled by the MRF are integrated into a new Bayesian model with better use of both spatial and spectral information.
3) Third, the efficiency of the new optimization approach for updating endmembers for each pixel and the abundance iteratively. In the comparison method of this article, the PSA and SAM methods should first unmixing the mixed pixel, and then make the SPM after obtaining the abundance of each pixel. These two methods strictly adhere to abundance constraints, where the numbers of subpixels of each land cover class within the mixed pixel corresponds to the input abundance map. However,  since the spectral unmixing is still an open problem, the obtained abundance map does not necessarily accurately predict the true proportion of land cover classes within the mixed pixel, so there are a lot of noise artifacts in the SPM results, which affect the performance of SPM techniques. The GAAI, SPMLM, MRFSPM, and the proposed methods are all based on the linear spectral unmixing models and perform SPM in the form of constructing objective functions. Their objective functions both contain the following two parts, one part is the spectral term, so that the SPM satisfies the constraints of the linear spectral unmixing model, and the other part is the spatial term, so that the discrete subpixel class label field produces the maximum spatial dependence and tends to smooth. These two parts are combined according to a certain weighing coefficient to achieve a best balance, relaxing the constraints on abundance, and eliminates the uncertainty caused by spectral unmixing. However, the spectral terms of the GAAI and SPMLM methods do not consider endmember variability and noise heterogeneity, while the MRFSPM method only considers noise heterogeneity, not endmember variability, while the spectral terms of the proposed method are considered both, such as (12). Therefore, the proposed method is more suitable for the original spectral information in spectral reconstruction, and then in the iterative optimization of abundance information estimation is more accurate, thus more accurate to estimate the details of the subpixel label field (such as the red rectangles part of the three experimental results), obtained the high-quality SPM results.

V. CONCLUSION
In this article, we have proposed a Bayesian SPM for HSI by integrating DEMM and MRF to address the SPM problems. First, the DEMM framework that allows the enhanced SPM performances by fully taking into account the patterns of endmember-abundance within each pixel for accommodating the endmember variability as well as the discrete subpixel hidden class label field while considering the noise heterogeneity effects. Consequently, the novel SPM model was less susceptible to the impact of endmember variability and band-dependent noise present in HSI. Second, we used MRF to model the discrete class label field, and then combined MRF with DEMM in a novel Bayesian model to efficiently investigate the spatial contextual and spectral information. By accommodating the discrete class label information of the subpixels, DEMM links SPM with LMM and allows the use of MRF to exploit the spatial information in the discrete class label field. Third, we designed a new EM approach that solve the resulting Bayesian model, where E-step estimates the abundance using the SA algorithm and M-step estimates the endmembers for each pixel in HSI using the ANLS approach. The visual and quantitative assessments on three different HSIs demonstrate that the proposed algorithm is an efficient SPM method. The analysis η of parameters indicates the robustness of the proposed algorithms. Future research directions include the fusion of a more accurate spatial regularization term into the Bayesian framework and the improvement of the efficiency and practicability of the algorithm and its application to large-scale images.