Application of spatially related MRF model in NMF hyperspectral unmixing

Aiming at Non-negative Matrix Factorization (NMF)’s problem of initialization and "local minima" in hyperspectral unmixing, a NMF linear unmixing algorithm with spatial correlation constrains (SCNMF) based on Markov Random Field (MRF) was proposed. Firstly, Hyperspectral Signal identification by minimum error (HySime) method was adopted to estimate the number of endmembers, initialized endmember matrix and abundance matrix by Vertex Component Analysis (VCA) and Fully Constrained Least Squares (FCLS) respectively. then established energy function to depict the spatial distribution characteristics of ground objects by MRF model. Finally, spatial correlation constraint based on MRF model and NMF standard objective function were combined in the form of altemating iteration to estimate endmember spectrum and abundance of hyperspectral image. Theoretical analysis and experimental results indicated that, the endmember decomposition precision of SCNMF is 10.6% higher than that of Minimum Volume Constrained NMF (MVC-NMF), 12.3% higher than that of Piecewise Smoothness NMF with Sparseness Constraints(PSNMFSC), 14.1% higher than that of NMF with Alternating Projected Subgradients(APS-NMF); the abundance decomposition precision of SCNMF is 14.4% higher than that of MVC-NMF, 15.9% higher than that of PSNMFSC, 15.3% higher than that of APS-NMF.The proposed SCNMF can remedy NMF's deficiency in describing spatial correlation characteristics, and decrease spatial energy distribution error.


Introduction
"Mixed pixel" is one of the key questions that affect the precision of Hyperspectral Remote Sensing Application. With the continuous expansion of hyperspectral data's application scope and improving demand of application precision, more and more unmixing methods emerge. Among them, Non-negative matrix factorization (NMF) proposed by Lee and Seung [1] is suitable for processing high dimensional data, and excludes meaningless negative elements. Its formula is similar to the linear spectral mixture model, which makes NMF based unmixing algorithm the hot spot of linear hyperspectral unmixing study.
Because of problems such as initialization and"local minima"caused by non-convexity, the unmixing effect of standard NMF is unsatisfactory, and numbers of NMF extension algorithms for hyperspectral unmixing emerge. These algorithms can reduce the solution space, speed up the iteration or improve the accuracy of the mixing, by adding some typical features of hyperspectral images into the standard NMF objective function as constraint terms. For example, Pauca [2] presented a NMF-based unmixing algorithm,which introduces constraints like abundance sparseness and endmember smoothness into the standard NMF algorithm; Miao [3] proposed a minimum volume constrained NMF (MVC-NMF) algorithm; Zymnis [4] proposed a NMF algorithm with altemating projected subgradients, APS-NMF); Jia [5] also proposed a piecewise smoothness NMF with sparseness constraint (PSNMFSC), which needs to specify the sparsity beforehand.
The problems of algorithms mentioned above include: treating hyperspectral data only as a collective of spectral or statistical information, which ignores the spatial features, especially spatial correlation features among different kinds of ground objects (endmembers) ； and the internal function terms of extended objective function may interfere with each other. Such problems to some extent restrict the further improvement of hyperspectral unmixing algorithms' accuracy.
In recent years, a large number of academic achievements have been achieved in the study of hyperspectral unmixing based on spatial correlation characteristics. For example, Du [6] proposed a semi supervised NMF algorithm based on Graph regularization(GSNMF), which overcomes the defect of ignoring the local geometric structure of sample data; Liu [7] put forward a new hyperspectral unmixing algorithm, by using training data to construct structured dictionary, which introduced space correlation constraints and spatial information of training data; Chen [8] used a priori probability density function to express the spatial correlation degree of two adjacent regions, and proposed a region related NMF unmixing algorithm; Jiang [9] proposed a semi supervised NMF algorithm based on sparseness constraint and graph regularization, which could maintain the geometry structure when processing low dimensional nonnegative decomposition.
This dissertation proposes a spatial correlation constrained NMF unmixing algorithm based on Markov Random Field (MRF) model, recorded as SCNMF. SCNMF can be classified as a non-supervised algorithm, which has low degree of dependence on prior knowledge. It uses NMF to realize unmixing procedure and ensure basic unmixing accuracy, adopts MRF model to depict the spatial correlation feature of hyperspectral image, rectifies unmixing error of NMF's standard objective function, and further improves unmixing accuracy.

NMF unmixing based on MRF model
Because of two-dimensional spatial image attributes, hyperspectral remote sensing images not only contain a lot of spectral and statistical information, but also contain abundant spatial information. For example, the natural objects has flaky and continuous feature of spatial distribution, which leads to similar values of adjacent image pixels, and the spatial distribution of endmembers has continuity and relevance.
The space energy of hyperspectral images depends on the frequency and amplitude of the ground objects' spatial variation. The lower the frequency and the smaller the amplitude, the smaller the space energy is. Because of the spatial correlation of endmembers, the variation frequency and amplitude of endmembers in local neighborhood are both small. In other words, in the ideal NMF decomposition results, the endmember's "space energy" should be as small as possible. Therefore, the smaller the energy function value (indicates space energy) is, the closer to the truth of the ground.
Dobrushin [10] pointed out that if the NMF unmixing algorithm does not consider the similarity among adjacent ground objects, and ignores the spatial correlation feature of endmember distribution, other types of endmember or noise will have more interference on endmember's spatial distribution in unmixing results, which may result in bigger space energy. A bigger space energy means bigger distribution proportion error of endmembers in mixed pixels.
A model which could describe the spatial correlation of endmember distribution need to be built, so that the endmember distribution's "space energy" will move more closer to the ideal state (minimum), and the coincidence degree between the endmembers' spatial distribution and the ground truth will be improved.

MRF model
MRF model is a powerful tool to establish joint probability distribution according to the local spatial correlation among adjacent pixels, its biggest advantage is that it cannot only make use of the characteristic information of the pixel itself, but also make use of the correlation information among adjacent pixels. Supposing 1 2 B x={x ,x ,...,x } as one realization of a random variable cluster ∈ is a neighborhood system defined on T.
If random field X meets condition (1)，then X is called a MRF with neighborhood system N: MRF model can provide statistical description of hyperspectral images, and depicts the conditional distribution (spatial correlation degree of endmember distribution) of each pixel on its adjacent pixels, and the conditional distribution is usually expressed as MRF's local characteristics. When the order of the neighborhood system is large enough, any hyperspectral image defined on T can be treated as a realization of MRF.

MRF model describing the spatial correlation characteristics
As there is no specific form of joint probability in the definition of MRF, it is inconvenient to use. According to Spitzer [11], Hammersley-Clifford theorem establishes the equivalence between MRF and Gibbs random fields, and the optimal solution problem of MRF model is transformed to the minimum value problem of Gibbs random field's energy function.
Different forms of energy function result in different MRF models. There are some common MRF models such as auto model and multi-level logic model (MLL). In this paper, the MRF model proposed by Deng [12] is used to describe the spatial correlation of endmembers, which divides the energy function into two parts, as shown in formula (2): In formula (2), F E is the energy function to represent the feature of the image itself, and only has relationship with the distribution of observational values, as shown in formula (3); R E is the energy function to represent the correlation between pixel and its neighborhood.
In formula (3), B represents band number； k i μ and k i σ represent the mean and standard deviation of F's kth feature. Supposing one image as a two-dimensional random field, the gray value of each pixel is affected by its neighborhood, and MLL is often used to characterize this type of correlation. In order to simplify the calculation, the model is usually assumed as homogeneous and isotropic,so the potential function is independent of the cluster's position and direction. R E 's calculation formula is shown as formula (4): In formula (4), P represents endmember number；if represents the order of neighborhood model. According to linear spectral mixed model (R represents hyperspectral data matrix ， M represents endmember spectral matrix ， S represents abundance matrix) ， suppose U as the separation matrix，Y as M's estimation，then： Substitutes formula (5) into formula (2), formula (3) and formula (4), and suppose W as the mean matrix of the eigenvectors corresponding to all pixels. The energy function of the MRF model with U as the independent variable could be inferred as ( ) E U , as shown in formula (6): As mentioned above, when the energy function takes the minimum value, the "space energy" described by the MRF model is the smallest, and the space distribution is optimal. Take the derivative of U in formula (6), and the simplified result is shown in formula (7).

SCNMF Model
SCNMF calculates the separation matrix U by alternately running standard NMF and MRF based spatial correlation constraint model. The NMF procedure selects the standard objective function (8) based on Euclidean distance as the basis for ensuring the accuracy of unmixing. The iterative rule adopts the maximum likelihood method, and the multiplicative iterative formulas of the endmember spectral matrix and the abundance matrix are respectively shown in formula (9) and formula (10).
In the MRF-based spatial correlation constraint step, we use the energy function model and related formulas in section 2.2 to correct the endmember distribution of the NMF unmixing result, and make it gradually approach to the ideal distribution state as the minimal space energy.
The steps of the SCNMF algorithm can be summarized as follows: 1) The Hyperspectral Signal identification by minimum error (HySime) method was adopted to estimate endmember number P; 2) Endmember matrix M and abundances matrix S was initialized, based on the estimated endmember number P in step 1), combined with Vertex Component Analysis (VCA) and fully constrained least squares (FCLS); 3) The iterative results of the separation matrix U was calculated by iterative procedure of standard NMF; 4) Every column of U was normalized, and the mean matrix W of the corresponding eigenvectors was estimated, then U was calculated according to formula (7); Repeat step 3) and step 4) until each stop criterion is met at the same time, and obtains an estimated ground component, convert it into a matrix, then an endmember's distribution is obtained. Continue to iterate until the threshold conditions are satisfied, then the estimated results of all endmembers' distribution are obtained.
The endmember number estimation algorithm (HySime) used in step 1) is a method to estimate the hyperspectral signal subspace. Although its computation procedure is complex, it does not need any parameters, has high adaptability and high accuracy. The principle and detailed implementation of HySime were given by Bioucas [13] .
VCA in Step 2) is a kind of endmember identification algorithm. Its advantages include short running time, high efficiency and high accuracy, and its drawback is dependence of pure pixels. However, the experimental results of Yuan [14] show that, even if the pure pixel does not exist, the accuracy and comprehensive performance of VCA as an initialization algorithm for NMF are also satisfying.
SCNMF abandons the traditional way of adding new constraint items in NMF objective function, and builds an independent auxiliary function, which running alternately with NMF objective function, and avoids the interference among internal function items.
The convergence of SCNMF is discussed below. NMF objective function's convergence was testified by Lee [15] , then the convergence of separation matrix U's iterative results can be ensured, and the eigenvector mean matrix W calculated by normalizing each column of U in step 4) is convergent.
It is easy to prove that formula (7) is convergent on W, therefore the objective function in step 4) is convergent. As step 3) and step 4) form together a parallel alternating iterative algorithm with initial boundary value, this algorithm can be abstracted as a Dirichlet-Neuman (D-N) alternate iteration, and the analysis of D-N iterative algorithm and its convergence by Wu [16] can prove the convergence of SCNMF.

Experimental results and analysis
By calculating the global Moran 'I index, the spatial correlation degree of experimental data is quantitatively described. Two sets of true hyperspectral data with different spatial correlation degree were adopted in the experiment. The bigger the global Moran 'I index, the higher degree of the image' spatial correlation, and vice versa.

Experiment 1
Experimental 1 adopted the hyperspectral image acquired by Hyperspectral digital imagery collection experiment (HYDICE) from Washington, D.C., as shown in Figure 1  Through visual interpretation, the image mainly contained three kinds of ground components, including vegetation, bare soil and water, and other types of ground components, such as linear roads and small buildings. In order to simplify the description, the ground components with very small area were not considered. In other words, all pixels were treated as only containing three main objects: vegetation, bare soil and water.
As HYDICE was an airborne imaging spectrometer with high spatial resolution, each end-member could be considered as containing a large number of pure pixels. Therefore, the paper collected the spectrum of each endmember by selecting reference points in the original image as reference value, and adopted full constrained least squares(FCLS) method to calculate end-member's abundance reference value. The difference of spectral curves between water and other two kinds of objects was very significant.
Degraded the image's spatial resolution accuracy to 0.1 times of original value by mean-value resampling, and plenty of mixed pixels were produced. Then rectified SCNMF's effect and accuracy by implementing unmixing experiment of new produced image with a large number of mixed pixels. The abundance reference value after resampling was shown as figure 2.
In Figure 2, pure white represented that endmember's area proportion in the pixel is 1 (100%), pure black represented that endmember's area proportion in the pixel is 0 (0%) , and the rest of grayscale corresponded to different proportions between 0 and 1 respectively, the meaning of different grayscale's representation in other abundance figures was the same as Figure 2. To demonstrate the effectiveness of SCNMF algorithm in the unmixing of hyperspectral images, three representative NMF algorithms, MVC-NMF, PSNMFSC and APS-NMF, were selected as reference algorithms. Similar to SCNMF, the other three algorithms introduced hyperspectral image's spatial information in the process of NMF unmixing in different principles and forms, and literature [3][4] [5] showed that the three algorithms were effective and have fine accuracy. The result of SCNMF's abundance decomposition was shown in Figure 3, white parts in each subgraph represented the location of a certain ground object in the image. It can be seen that the three subgraphs were basically consistent with true values in Figure 2, which indicating that the SCNMF can effectively separate the three objects with different spectral features. Table 1 and table 2 illustrated the accuracy of endmember spectrum and abundance distribution's decomposition results respectively. Spectral Angle Distance (SAD) was adopted to calculate the endmember's spectral estimation accuracy, and the Root Mean Square Error(RMSE) was used to calculate endmember's abundance estimation accuracy.
Formula of SAD： Formula of whole image's RMSE: And î is the residual error of the i th pixel in the image. As can be seen from table 1~2, the accuracy of SCNMF was the highest in four algorithms , either endmembers' spectrum or abundance. Speaking specifically, the accuracy of SCNMF's endmember spectral decomposition was 10.6% higher than MVC-NMF, 12.3% higher than PSNMFSC, 14.1% higher than APS-NMF, and SCNMF's endmember abundance accuracy was 14.4% higher than MVC-NMF, 15.9% higher than PSNMFSC, and 15.3% higher than APS-NMF.

Experiment 2
The natural scene of experimental data 1 was dominated by natural scenes, such as water, vegetation and bare soil, with significant spatial correlation. For comparing and verifying unmixing effectiveness of SCNMF of images with different degrees of spatial correlation, an Airborne Visible Infrared Imaging Spectrometer(AVIRIS) hyperspectral image acquired in July 1995, mainly mining area of Nevada Cuprite in American was selected in experiment 2, as shown in figure 4. The global Moran 'I index of experimental data 2 is 0.2247. Compared with experimental data 1, data 2's spatial correlation degree was obviously reduced, indicating that the distribution of minerals in the mining area was scattered and chaotic.
The words "A, B, C, K, M" in figure 4 represent the distribution of 5 widely distributed minerals, including Alunite, Buddingtointe, Calcite, Kaolinite and Muscovite respectively. The image had 400 columns, 350 rows, spatial resolution was 20m, wavelength range was 1.99~2.48 μ m , and spectral resolution was 10nm. It had 50 bands as the 172~221th bands of AVIRIS image's original spectrum. This area was located in the southern Nevada of the United States, its surface was mostly bare minerals without vegetation. Compared with the HYDICE hyperspectral data in experiment 1, its spatial correlation characteristic was not obvious, and the spatial distribution of main minerals in this image was random.
Because this AVIRIS image's spatial resolution (20m) was relatively high, it can be assumed that every kind of end-member had a certain amount of pure pixels.Therefore, the paper selected reference points manually in the original image to collect spectrum of each end-member as reference value, and use Full Constrained Least Squares(FCLS) method to calculate the reference value of abundance .
As for data 2, The experimental procedure of experimental 1 was repeated to compare and verify SCNMF's performance for different scenarios. The experiment was focused on the above 5 main mineral types, and corresponding treatment of other small objects was identical to data 1. Figure 5 was the result of abundance estimation corresponding to SCNMF method. Due to length limitations, the results of other reference methods were not listed one by one.
Compared SCNMF's abundance estimation of each end-member in Figure 5 with the abundance reference values calculated by FCLS, their similarity was over 85%, which indicating that SCNMF can effectively separate the 5 different land types and their general location.  Table 3 and table 4 listed the precision analysis results of 4 unmixing methods. It could be seen that all of the 4 algorithms can effectively decompose the 5 main types of minerals, and SCNMF had the best performance. The end-member spectrum accuracy of SCNMF was 7.82%, 12.4% and 10.1% higher than that of MVCNMF, PSNMFSC, and APS-NMF respectively. And the abundance estimation accuracy of SCNMF was increased by 8.34%, 12.6% and 10.1%, MVCNMF, PSNMFSC, and APS-NMF respectively.
The results listed above indicated that, comparing with MVCNMF, PSNMFSC, and APS-NMF, SCNMF still had a significant accuracy advantage when experimental data's space correlation degree decreased evidently, yet slightly lower than experimental data with significant spatial correlation.

Conclusion
Aiming at the spatial correlation characteristics of endmember distribution which were easy to be ignored in many NMF unmixing methods, the paper adopted MRF model to describe spatial correlation, and proposed a new NMF based hyperspectral unmixing method, by constructing energy function of endmember's spatial distribution to characterize spatial correlation characteristics, and the auxiliary function was independent with NMF's objective function, which abandoned traditional practice of constructing new internal constraints of NMF's objective function.
The experimental results showed that, for most true hyperspectral data, the decomposition accuracy of proposed method was better than other representative methods like MVCNMF, PSNMFSC, and APS-NMF.
Although spatial correlation features are common in hyperspectral images, not all hyperspectral images have significant spatial correlation. Experimental results also showed that the performance of proposed method decreased slightly when the spatial correlation was not high. It can be predicted that if the spatial correlation between pixels of experimental images is close to zero or negative correlation, the performance of SCNMF would have severe degradation. How to further stabilize SCNMF's performance will be the main content of future research work.