Multiple Endmember Hyperspectral Sparse Unmixing Based on Improved OMP Algorithm

，，

（College of Information and Communication Engineering，Harbin Engineering University，Harbin 150001，China）

Multiple Endmember Hyperspectral Sparse Unmixing Based on Improved OMP Algorithm

Chunhui Zhao∗，Haifeng Zhu，Shiling Cui and Bin Qi

（College of Information and Communication Engineering，Harbin Engineering University，Harbin 150001，China）

In conventional linear spectral mixture analysis model，a class is represented by a single endmember. However，the intra-class spectral variability is usually very large，which makes it difficult to represent a class，and in this case，it leads to incorrect unmixing results.Some proposed algorithms play a positive role in overcoming the endmember variability，but there are shortcomings on computation intensive，unsatisfactory unmixing results and so on.Recently，sparse regression has been applied to unmixing，assuming each mixed pixel can be expressed as a linear combination of only a few spectra in a spectral library.It is essentially the same as multiple endmember spectral unmixing.OMP（orthogonal matching pursuit），a sparse reconstruction algorithm，has advantages of simple structure and high efficiency.However，it does not take into account the constraints of abundance non-negativity and abundance sum-to-one（ANC and ASC），leading to undesirable unmixing results.In order to solve these issues，this paper presents an improved OMP algorithm（fully constraint OMP，FOMP）for multiple endmember hyperspectral sparse unmixing.The proposed algorithm overcomes the shortcomings of OMP，and on the other hand，it solves the problem of endmember variability. The ANC and ASC constraints are firstly added into the OMP algorithm，and then the endmember set is refined by the relative increase in root-mean-square-error（RMSE）to avoid over-fitting，finally pixels are unmixed by their optimal endmember set.The simulated and real hyperspectral data experiments show that FOPM unmixing results are ideally comparable and abundance RMSE reduces much lower than OMP and simple spectral mixture analysis（sSMA），and has a strong anti-noise performance.It proves that multiple endmember spectral mixture analysis is more reasonable.

hyperspectral image；sparse representation；multiple endmember spectral unmixing；OMP；ANC and ASC

1 Introduction

In recent years，hyperspectral remote sensing technology develops rapidly.It has been used successfully in food safety，quality control，and forensic application［1-3］.Mixed pixels are widespread in hyperspectral images due to the diversity and complexity of ground and the limitations of the sensor’s spatial resolution.Mixed pixels have a serious impact on the accuracy of subsequent hyperspectral image processing.Therefore，how to solve mixed pixels effectively is one of the most important problems in hyperspectral image analysis［4］.Simple linear spectral mixture analysis（sSMA）［5］models a mixed spectrum as a linear combination of pure spectral signatures，weighted by their sub-pixel fractional cover.Once the endmembers have been extracted，sub-pixel cover distribution maps can be generated using approaches like singular value decomposition，Gramm-Schmidt orthogonalization，maximum likelihood or least-squares regression analysis.

However，the lack of ability to account for sufficient temporal and spatial variability between and among endmember spectra has been acknowledged as a major shortcoming of sSMA approach with fixed endmembers.Over the past decades numerous efforts have been made to circumvent this issue［6］.

One of the first attempts to address endmember variability is Multiple Endmember Spectral Mixture Analysis（MESMA）［7-8］model，which is an iterative mixture analysis algorithm.Monte Carlo spectral unmixing model（AutoMCU）［9］and Bayesian Spectral Mixture Analysis model（BSMA）［10］is in parallel withMESMA.These iterative mixture analysis techniques explore all possible endmember combinations in search for the best solutions and the high number of iterations strongly slows down the analysis.Different spectral feature selections and spectral transformations algorithms have also been proposed to reduce the effect of endmember variability［11-12］，for example AutoSWIR and Stable Zone Unmixing（SZU）.The SZU has been proved to be good，but the parameters variability affects the unmixing accuracy greatly.For spectral transformations，there are Normalized Spectral Mixture Analysis（NSMA）and Derivative Spectral Unmixing（DSU）in which second derivative endmember spectra are used as inputs into SMA to emphasize the inter-class variability，while reducing the intra-class variability.The physical meaning of these two algorithms is not clear and the unmixing results are unsatisfactory.The other one to reduce the effect of endmember variability is Weighted Spectral Mixture Analysis（wSMA），prioritizing spectral bands less sensitive to endmember variability［13］.

The above spectral unmixing methods have played a positive role in reducing the effect of endmember variability，but there are always limitations.In recent years，with the development of sparse representation theory，scholars have a more profound understanding of sparse unmixing.Beginning in 2009，some scholars try to add sparse constraint to the sSMA，and have achieved certain results，but to our best knowledge，there are almost no studies of adding sparse constraint to multiple endmember spectral mixture analysis. Moreover，most sparse unmixing methods are based on convex relaxation methods，while little attention has been paid to the use of greedy algorithms（GAs）. OMP［14-15］is one of the representatives of GAs，with a simple structure and high efficiency，which finds the best matches for each pixel by iteration from the existing library.When applying OMP in hyperspectral image unmixing，every pixel can get its optimal endmember set and their corresponding abundance. Obviously，it is essentially the same as multiple endmember spectral unmixing which iteratively produces the optimization for each pixel.However，there are big drawbacks about OMP that the abundance does not satisfy the ANC and ASC，besides the optimal endmember set is always larger than real situation，leading to accurate unmixing results.

Given the shortcomings of the above overcoming intra-class spectral variability algorithms and the consistency of the iterative multiple endmember spectral unmixing and the OMP unmixing，this paper presents an improved OMP algorithm（called fully constraint OMP，FOMP）in multiple endmember spectral unmixing，which adds two constraints to the original OMP algorithm，and then finds the optimal number of endmembers for each pixel by the relative increase in RMSE to prevent over-fitting and increase the unmixing accuracy.

2 OMP Sparse Unmixing of Hyperspectral Image

Spectral unmixing aims at estimating the fractional abundances of pure spectral signatures in each mixed pixel.The linear mixing model assumes that the observed spectrum of a pixel can be expressed as a linear combination of the spectra of the endmembers presented in the respective pixel［16-17］.It can be expressed mathematically as follows：

whereyis an observation mixed pixel；Ais the endmember collection and the total amount isM；xis the corresponding fractional abundance，andnrepresents the error term.ANC and ASC are usually applied to the abundance.

The linear sparse unmixing model［18-19］assumes that the observed spectrum of a mixed pixel can be expressed as a linear combination of only a few spectral signatures selected from a library known in advance.As the number of endmembers involved in a mixed pixel is usually very small compared with the dimensionality of the spectral library，the vector of fractional abundancesxis sparse.

Without considering the two constraints in Eqs.（4）and（5），the optimization problem of sparse unmixing is then：

Three kinds of methods can solve this problem，namely，greedy algorithms（GAs），convex relaxation methods，and sparse Bayesian methods.Convex relaxation methods and sparse Bayesian methods are far more complicated than the GAs.OMP，a typical GA，has been utilized to solve theP0problem.OMP is simple，with high efficiency，making it more practical in real applications.OMP iteratively selects at each step the member in the spectral library best correlated with the residual part of the mixed pixel.Then，it produces a new approximation by projecting the spectral vector of the mixed pixel onto the potential endmembers which have already been selected.

3 Fully Constraint OMP Algorithm

In this section，we firstly propose the concrete mathematical expression of linear multiple endmember spectral mixture analysis model，and then analyze the unmixing process in detail.The endmember matrix in the sSMA is fixed，however，in fact，the number and types of endmember are allowed to vary for different pixels.If we unmix the pixels by fixed endmember set，the unmixing abundance results are inaccurate. Apparently multiple endmember spectral mixture analysis is more reasonable.

The linear multiple endmember spectral mixture analysis is as follows：

wherenis the total number of hyperspectral image pixels andqdenotes a constant，whereq＝1 denotes the endmember belongs to the pixel，otherwiseq＝0 denotes the endmember does not belong to the pixel；is thek-th pixel andMis the total number of endmembers；Nidenotes the number of spectra in thei-th intra-class；ei，jdenotes thej-th spectrum in thei-th intra-class；denotes the spectral abundance of thej-th endmember in thei-th intra-class，andεkdenotes an additive perturbation.

OMP iteratively selects at each step the member in the spectral library best correlated with the residual part of the mixed pixel.It finds the optimal set of endmember for each pixel from the spectral library，and then unmixes each pixel by the endmember set. We can see that it is the same as the multiple endmember spectral mixture analysis.So we combine the two algorithms.However，OMP does not take into account the two constraints，which leads to negative abundance.

The unmixing residuals can be written as Eq.（8）：

Then，Eq.（8）can be restated as Eq.（10）：

Eq.（9）is unconstrained least squares abundance inversion，which leads to negative abundance. According to convex geometry theory，the mixed pixels locate inside the simplex composed of the related endmembers.If a mixed pixel is inside the simplex， unconstrained least squares and fully constraint least squares get the same optimal solution.The opposite condition，if a mixed pixel lives outside the simplex，unconstrained least squares will get unreasonable consequence compared to fully constraint least squares. So we take full constraint least squares instead.By doing so，residuals will not orthogonal to the selected set of endmembers，to solve which，we project the residuals to the orthogonal subspace of the already selected set of endmembers.

Eq.（8）can be rewritten as follows：

However，OMP often leads to more endmembers than actually present within a pixel.If more endmembers are used than required，the RMSE values may indeed be smaller，yet the fractional error will be higher.To deal with this，we refine the endmember set by calculating the relative increase in RMSE，reducing the number of endmembers by one.For each pixel，K-，（K-1）-，…，2-and 1-endmember reconstruction RMSE are compared［20-21］.The relative increase in RMSE by going fromK-to（K-1）-endmember and from 2-to 1-endmember are calculated as follows Eqs.（14）and（15）：

whereRMSE（k-1_EM）is the reconstruction RMSE of（k-1）-endmember andRMSE（k_EM）is the reconstruction RMSE ofk-endmember；Lis the number of bands；yiis a pixel andis its estimation.If theΔRMSEis less than a predefined threshold，it means that the extra endmember makes little contribution to the pixel，so the less endmembers are selected to unmix the corresponding pixels，otherwise，the more endmembers are selected.

FOMP algorithm procedure is as follows：

Input：

Dictionary：A＝［a1，a2，…，aM］

Signals to represent：y

Sparsity：K

Initialization：

Residual：r0＝y

Iteration counter：k＝1；

Procedure：

Iteration，in thekth cycle，run the following 1）-6）.

If the maximum occurs for multiple indices，break the tie deterministically.

2）Augment the index set：

3）Update residuals：

Of which hyperFcls means fully constraint least squares abundance inversion.

The orthogonal subspace of the selected endmember is as follows：

4）ComputeΔRMSEfor each pixel.

5）Find the optimal number of endmembers byΔRMSEkfor each pixel.

6）Incrementk＝k+1，and return to Step 1）ifk＜K．

Output：

Index set：Λ

Sparse coefficient：α＝hyperFcls（E，y）

4 Experiments

4.1 Experiments of Simulated Hyperspectral Data

In order to verify the effectiveness of the algorithm，two experiments are made in this section. We select five mineral spectra from the United States Geological Survey（USGS）library to form simulated hyperspectral data.These spectra are Almandine WS477（A1），Almandine WS478（A2），Nontronite GDS41（N3），Buddingtonite GDS85 D-206（B4）and Buddingtonite NHB2301（B5）and the dictionary is the USGS library.

4．1．1 The first experiment

The first set of experimental data is constituted randomly of the fiver mineral spectra above and the numbers is 10 000.Besides，35 dB Gaussian white noise is added.

Table 1 gives a quantitative evaluation about the abundance RMSE of the five endmembers.We can find that the RMSE of FOMP is the lowest among all，besides，FOMP and OMP perform beyond the sSMA，which means they can address the endmember variability issue to some degree.

Fig.1 gives the abundance error histograms of the three methods at 35 dBSNR sceneries.It gives a direct comparison between these three methods.FOMP algorithm abundance error is mostly clustered around the value 0，while the other algorithms error distributions are scattered.

Table 1 Unmixing abundance RMSE contrast

Fig.1 Histograms of abundance error

4．1．2 The second experiment

In thesecond experiment，we form 3 000 mixed pixels in which the 1 th-1 000 th pixels are made up ofA1andN3by random proportion satisfying the ANC and ASC.The 1 001th-2 000th pixels are constituted byA1andB5and the remaining pixels are composed byA2，N3andB4in which the random proportion is made in the same manner as the former.In order to verify theanti-noise performance，30-55 dB Gaussian white noise are added.

Table 2 gives a quantitative evaluation of the RMSE for the five endmembers between the real endmember abundance and the estimated endmember abundance at different SNR（signal to noise ratio）. Longitudinal view，we can find that the RMSE of FOMP is the lowest among all，where abundance RMSE is reduced about 90%and 40%compared to sSMA and OMP，besides，FOMP and OMP perform beyond the sSMA，which means they can address the endmember variability issue to some degree.Lateral view，the abundance error ofN3is much smaller than the others，we know that synthetic data only involve one spectrum of this mineral class，so we can get that intra-class spectral variability indeed affect the accuracy of unmixing results.

Table 2 Unmixing abundance RMSE contrast at different SNR

Fig.2 intuitively shows the average unmixing abundance RMSE as a function of SNR.As shown，all the algorithms get worse performance with the SNR getting lower，but our proposed algorithm can obtain the best estimation of the average abundance at different noise scenarios among all the methods. However，when introducing the random noise，although the accuracy of FOMP decreases，its curve is the most flat one，meaning that it possesses the best anti-noise performance.Meanwhile，OMP and sSMA are sensitive to noise.

Fig.3 gives the abundance error histograms of the three methods at 30 dB，40 dB，50 dB SNR sceneries. It gives a direct comparison between these three methods.FOMP algorithm abundance error is mostly clustered around the value 0，while the error distributions of the other algorithms are scattered.With the decline of SNR，all these methods’number of pixels deviation gradually from 0 value increases，but FOMP algorithm is always better than the other algorithms.

Fig.2 Average unmixing abundance RMSE as a function of SNR

Fig.3 Histograms of abundance error

4.2 Experiments of Real Hyperspectral Data

Our proposed method has been applied to a hyperspectral image taken over northwest Indiana’s Indian Pine Test Site shown in Fig.4.The data consist of 145×145 pixels with 220 bands which are initially reduced to 200 by removing bands，covering water absorption as well as noisy bands and the total number of endmember class is 16.Three typical types of crops：hong-windrowed，woods and corn，respectively 489，1 294，834 sample points are chosen because they can best verify the effectiveness of the improved algorithm. Their spatial distribution of the label figures are shown in Fig.5（a）.We randomly select ten spectra from the total sixteen classes to form the dictionary.So the dictionary contains 160 spectra.

Table 3 gives a quantitative evaluation about the abundance RMSE of the ten endmembers for the three crops.It can be seen that FOMP performs the lowest RMSE compared with others.In the unmixing distribution maps，the brighter gray represents the abundance of the corresponding categories in the mixed pixels the greater，and vice versa.

Fig.4 Pseudocolor composed map

Table 3 Unmixing abundance RMSE contrast

Fig.5 shows a more intuitive representation.FOMP algorithm represents the three crops better，and the unmixing results are much closer to the true spatialdistribution.Meanwhile，the other algorithms unmixing results are vague and cannot distinguish different crops effectively.However，FOMP and OMP are better than the sSMA in both the average RMSE and the unmixing distribution maps.Above all，except for sSMA，they all have played a positive role in overcoming the intra-class spectral variability and implying multiple endmember spectral mixture analysis is more reasonable.

Fig.6 shows the abundance error histograms of the three methods.It gives a direct comparison between these three methods.FOMP algorithm abundance errors are mostly clustered around the value 0，while the error distribution of other algorithms is scattered.

4.3 Analysis of Unmixing Results

By comparing the experimental results above，we can get that FOMP algorithm is better than the other algorithms in unmixing abundance error，which takes ANC and ASC into account and adoptsΔRMSEto refine the optimal endmember set to avoid over-fitting where the number of endmembers is greater than the actual number in each pixel.sSMA gets the worst unmixing results without considering the intra-class spectral variability.OMP algorithm does not consider ANC and ASC，causing negative abundance，besides，the iterative optimization process of OMP results in over-fitting.Although the reconstruction error is cut down，the unmixing abundance error increases.OMP may not accurately unmix some signals with low SNR or measured insufficiently.Besides，it is strict to build a dictionary because of its use of orthogonalization.

Fig.5 Unmixing distribution maps of hong-windrowed，woods and corn

5 Conclusions

This paper has proposed a multiple endmember hyperspectral sparse unmixingmethod based on the improved OMP algorithm.It firstly adds ANC and ASC to the OMP algorithm to receive reasonable abundance，and then finds the optimal number of endmembers for each pixel by the relative increase in RMSE to avoid over-fitting.The simulated and real hyperspectral image data experimental results have showed that FOMP is superior to the other two and is effective in overcoming intra-class spectral variability.We can see that multiple endmember spectral mixture analysis is more reasonable.Since the introduction of fully least squares abundance inversion，the computation is large，how to improve its efficiency is the focus of our future work.

［1］Zhao Chunhui，Cheng Baozhi，Yang Weichao.Algorithm for hyperspectral unmixing using constrained nonnegative matrix factorization.Journal of Harbin EngineeringUniversity，2012，33（3）：377-382.（in Chinese）

［2］Cheng Baozhi，Zhao Chunhui.A particle swarm optimization clustering-based approach for hyperspectral image anomaly targets detection.Journal of Optoelectronics Laser，2013，24（10）：2047-2054.（in Chinese）

［3］Wang Liguo，Zhang Jing.An improved spectral unmixing modeling based on linear spectral mixing modeling.Journal of Optoelectronics Laser，2010，21（8）：1222-1226.（in Chinese）

［4］Zhang Bing，Gao Lianru.Hyperspectral Image Classification and Target Detection.Beijing：Science Press，2011.201-213.（in Chinese）

［5］WangLiguo，Liu Danfeng，Wang Qunming.Geometric method of fully constrained least squares linear spectral mixture analysis.IEEE Transactions on Geoscience And Remote Sensing，2013，51（6）：3558-3566.

［6］Wang Liguo，Zhao Chunhui.Processing Techniques of Hyperspectral Imagery.Beijing：National Defense Industry Press，2013.13-17.（in Chinese）

［7］SomersB，Asner G P，Tits L，et al.Endmember variability in spectral mixture analysis：a review.Remote Sensing of Environment，2011，115（7）：1603-1616.

［8］Roberts D A，Gardner M.Mapping chaparral in the Santa Monica Mountains using multiple endmember spectral mixture models.Remote Sensing of Environment，1998，65（3）：267-279.

［9］Asner G P，Lobell D B.A biogeophysical approach for automated SWIR unmixing of soils and vegetation.Remote Sensing of Environment，2000，74（1）：99-112.

［10］Rao Yuhan，Chen Jin，Chen Xuehong.Quantitative assessment of the different methods addressing the endmember variability.IEEE International Geoscience and Remote Sensing Symposium（IGARSS）.Piscataway：IEEE，2013.3317-3320.

［11］Somers B，Delalieux S，Verstraeten W W.An automated waveband selection technique for optimized hyperspectral mixture analysis.International Journal of Remote Sensing，2010，20（31）：5549-5568.

［12］ShiChen，Wang Le.Incorporating spatial information in spectral unmixing：A review.Remote Sensing Environment，2014，149：70-87.

［13］Song Meiping，Zhang Yongrong，An Jubai，et al.Effective endmember based bilinear unmixing mode.Spectroscopy and Spectral Analysis，2014，34（1）：196-200.（in Chinese）

［14］Raksuntorn N，Du Qian，Younan N，et al.Orthogonal matching pursuit for nonlinear unmixing of hyperspectral imagery.2014 IEEE China Summit＆International Conference on Signal and Information Processing（ChinaSIP）.Piscataway：IEEE，2014，7：157-161.

［15］IordacheM D，Bioucas-Dias J M，Plaza A.Sparse unmixing of hyperspectral data.IEEE Transactions on Geoscience And Remote Sensing，2011，49（6）：2014-2039.

［16］Ding Haiyong，Shi Wenzhong.Fast N-FINDR algorithm for endmember extraction based on chi-square distribution. Journal of Remote Sensing，2013，17（1）：130-137.（in Chinese）

［17］ZareA，Ho K C.Endmember variability in hyperspectral analysis.IEEE Signal Processing magazine，2014，31（1）：95-104.

［18］Gharavi-Alkhansari M，Huang T S.A Fast orthogonal matching pursuit algorithm.IEEE International Conference on Speech and Signal Processing，1998，3，1389-1392.

［19］Shi Zhenwei，Tang Wei，Duren Zhana，et al.Subspace matching pursuit for sparse unmixing of hyperspectral data. IEEE Transactions on Geoscience And Remote Sensing，2014（52）：3256-3274.

［20］Wu Bo，Zhou Xiaocheng，Zhao Yindi.Study on the relationships between endmember variance and decomposition accuracy of mixture pixel.Remote Sensing Information，2007（3）：3-7.（in Chinese）.

［21］Demarchi L，Canters F，Chan J C-W，et al.Mapping sealed surfaces from CHRIS／Proba data：a multiple endmember unmixing approach.Hyperspectral Image and Signal Processing：Evolution in Remote Sensing（WHISPERS），2010.1-4.

TN911.73

：

：1005-9113（2015）05-0097-08

10.11916／j.issn.1005-9113.2015.05.015

2014-07-07.

Sponsored by the National Natural Science Foundation of China（Grant No.61405041，61571145），the Key Program of Heilongjiang Natural Science Foundation（Grant No.ZD201216），the Program Excellent Academic Leaders of Harbin（Grant No.RC2013XK009003），the China Postdoctoral Science Foundation（Grant No.2014M551221）and the Heilongjiang Postdoctoral Science Found（Grant No.LBH-Z13057）.

∗Corresponding author.E-mail：zhaochunhui＠hrbeu.edu.cn.