Nonlinear extended blind end-member and abundance extraction for hyperspectral images

Highlights

• Hyperspectral images are essential tools in different application frameworks.

• Nonlinear blind unmixing is studied by a multilinear mixture model.

• Three subproblems: end-members, abundances, and nonlinear interaction levels.

• A cyclic coordinate descent algorithm and quadratic constrained optimizations are used.

• Synthetic and remote sensing datasets, and a biomedical imaging application were evaluated.

Abstract

Hyperspectral images had become an essential tool in different application frameworks, such as mineral exploration, food inspection, and medical assessment, among others. However, the interpretability of these images involves an initial processing stage to model the optical interaction and analyze the spectral information. In this work, we study nonlinear unmixing of hyperspectral images by a multilinear mixture model (MMM). In this sense, we propose a nonlinear version of the extended blind end-member and abundance extraction (NEBEAE) method for blind unmixing, i.e. estimation of the end-members, their abundances, and the nonlinear interaction levels. In the problem formulation, we include a normalization step in the hyperspectral measurements for the end-members and abundances to improve robustness. The blind unmixing process can be separated into three estimation subproblems for each component in the model, which are solved by a cyclic coordinate descent algorithm and quadratic constrained optimizations. Each problem is mathematically formulated and derived to construct the general nonlinear iterative unmixing technique. We evaluated our proposal with synthetic and experimental datasets from the remote sensing literature (Cuprite, Urban and Pavia University scene datasets) and a biomedical hyperspectral imaging application. In the validation stage, we compared NEBEAE with three state-of-the-art methods to show its advantages in terms of precision and computational time.

Introduction

Hyperspectral images (HSIs) have become a valuable tool to identify qualitative and quantitative mappings of the constitutive elements of a scene, which are indistinguishable in plain sight [1]. The applications of these imaging modalities are quite diverse: mineral exploration, earth observation, seafloor exploration and monitoring, food inspection, biological tissue characterization, artwork investigation, etc. [2], [3], [4], [5], [6], [7], [8], [9], [10]. However, for a proper evaluation of hyperspectral information, a processing stage is required to estimate the spectral signatures of the constitutive elements or end-members, and their fractional contributions or abundances in each pixel of the image [11], [12]. The physical relation between end-members and their abundances is represented by a mathematical model, which describes the optical paths and interactions of the reflected light by the objects in the scene [13].

The basic mathematical model assumes a linear mixing model (LMM), so each measurement is represented by a linear combination of end-members plus an uncertainty/noise component [14], [15]. However, under some conditions, the light captured by the camera shows multipath scattering and nonuniform reflection patterns, so a nonlinear mixture model (NMM) should be pursued. Since late 2000’s and early 2010’s, diverse research groups have proposed NMM’s to highlight different nonlinear physical phenomena: Fan et al. model [16], generalized bilinear model (GBM) [17], polynomial post-nonlinear model (PPNM) [18], and recently multilinear mixing model (MMM) [19]. In all previous NMMs, the model preserves a linear component and includes a nonlinear term depending on the assumed physical interaction. In particular, the MMM shows an evolution of the GBM and PPNMM frameworks to describe all possible optical interactions in the scene. One advantage of the MMM is the addition of a single parameter to quantify the nonlinear phenomena per pixel.

Quantitative characterization of HSIs involves the estimation of end-members and their abundances, and the specific parameters of the NMM [13]. In some approaches, the end-members are assumed known in advance by early studies of the materials in the scene and their spectral signatures. This scenario is called a supervised strategy. On the other hand, if the end-members are unknown, and they have to be jointly estimated, this condition is denoted as an unsupervised approach, or a blind unmixing methodology (BUM) [20].

Recently, the joint application of artificial intelligence (AI) and non-convex modeling in hyperspectral imaging was reviewed in Hong et al. [21], where the main challenges and open problems in image restoration, dimensionality reduction, spectral unmixing, data-fusion and enhancement, and cross-modality learning were discussed. In this way, a deep learning perspective was considered for spectral unmixing in Hong et al. [22], where an end-member guided unmixing network (EGU-NET) was proposed that relies on a two-stream Siamese deep network for blind unmixing. In this work, the deviations of a LMM are denoted as spectral variability and by the structure of the EGU-NET, a nonlinear unmixing of the studied HSI is achieved. A different perspective is followed in Yao et al. [23], where a tensor-based approach was used for spectral unmixing to have a more adequate modeling in both the spatial and spectral domains. The unmixing methodology in Yao et al. [23] is called sparsity-enhanced convolutional decomposition, where the HSI is modeled as a third-order tensor, but no high-order interactions are taken into account in the end-members. Meanwhile, in Yao et al. [24], sparseness of the abundance maps and non-negative matrix factorization (NNMF) are taken into account for spectral unmixing, but once more without high-order components. Recently, in Hong et al. [25], a general multimodal deep learning (MDL) framework was proposed for classification and spatial information modeling in remote sensing datasets. Finally, a new training scheme was suggested in Hong et al. [26] for graph convolutional networks in HSI classification tasks.

Some recent nonlinear unmixing methodologies are reviewed in what follows [27]. As described in Hong et al. [21], nonlinearity in spectral unmixing is an open and pending problem in hyperspectral imaging. In this sense, a bilinear mixture model was considered in Yang et al. [28] to propose an unsupervised scheme that relies on geometric projection for nonlinear unmixing, and NNMF for end-members estimation. In this work, the concepts of virtual nonlinear vertex and nonlinear hyperplanes are adopted to derive the synthesis procedure. Also based on a bilinear mixture model, in Sigurdsson et al. [29], a blind sparse nonlinear hyperspectral unmixing (BSNHU) is suggested that relies on iterative cyclic descent algorithms and the ${\ell }_q$-regularizer to obtain sparse abundances. In [30], an augmented LMM is proposed to consider high-order interactions by incorporating a spectral variability dictionary based on prior knowledge. A different approach has also been suggested in Licciardi et al. [31], where first a nonlinear dimensionality reduction is applied to consider the subsequent estimation of end-members by the N-finder algorithm (N-FINDR) [32]. Finally, the abundances estimation was achieved by constrained quadratic optimization (CQO). In [33], a supervised unmixing was achieved by using continuum removal and natural logarithmic operations to eliminate the nonlinear effects in the HSIs, followed by a linear unmixing strategy. A nonlinear unmixing based on a version of PPNM was proposed in Dixit and Agarwal [34], where a Bayesian framework was used for abundances, nonlinearity factor, and noise variance estimation.

Meanwhile, the state-of-the-art in non-linear unmixing methodologies focused on the MMM is described next. In [19], a supervised scheme was pursued to estimate the abundances and nonlinear interaction terms in the MMM by numerical optimization. In [35], another supervised unmixing approach with a MMM was formulated by Laplacian graph regularizers and alternating direction method of multipliers. This proposal has two distinctive features: the sparsity of the abundances and the use of superpixels to handle large HSIs. To improve noise robustness, a robust unmixing scheme with MMM was proposed in Li et al. [36] by a supervised perspective, where ${\ell }_{2,1}$ norm-based loss function is adopted for the synthesis problem and also the sparseness of abundances is taken into account. To our knowledge, Wei et al. suggested the only unsupervised strategy for nonlinear unmixing based on the MMM [37]. In this work, a block coordinate descent methodology was used, where the abundances and end-members were updated by a gradient projection method, and the nonlinear interaction level by constrained least-squares.

In this context, the MMM provides a general framework to quantify the nonlinear optical interactions in the scene, but the research work on unsupervised unmixing strategies departing from this model is limited to the work by Wei et al. [37]. Moreover, new BUMs based on the MMM need to take into account precision, robustness, and computational time, as performance requirements. Our proposed BUM follows the estimation framework of the Extended Blind End-member and Abundance Extraction (EBEAE) methodology [20], but is now adapted to the MMM to address these performance issues. The new methodology is called nonlinear EBEAE (NEBEAE). One advantage of the original formulation of EBEAE was the measurements’ normalization to restrict the search space and to limit the variability of the estimation information. Now, we derive a similar condition for NEBEAE and establish a restriction on the abundances. The synthesis optimization objective function has three components: estimation error, entropy of abundances, and similarity of end-members. Due to the non-linear dependence on decision variables, a cyclic coordinate descent algorithm (CCDA) is used to derive an iterative scheme [38]. The implementation details of NEBEAE with respect to initialization, convergence, and measurements reconstruction are discussed in this work. Finally, we evaluate NEBEAE in three different scenarios: 1) synthetic datasets with different types and noise levels; 2) three experimental benchmarks in the remote sensing literature (Cuprite, Urban and Pavia University scene datasets [39]); 3) a biomedical hyperspectral imaging application. In the validation stage, we compared NEBEAE to Heylen and Scheunders [19] with the estimation of end-members by vertex component analysis (VCA), Sigurdsson et al. [29], Wei et al. [37] to show its advantages in terms of precision, robustness and computational time.

The mathematical notation used along this document is described next. The sets of real numbers, real $L$-dimensional vectors, and $L\times N$ real matrices are expressed as $\mathbb{R},{\mathbb{R}}^L$ and ${\mathbb{R}}^{L\times N}$, respectively. Scalars and vectors are described by italic and boldface lower-case letters, respectively. Matrices are denoted by upper-case boldface letters, and sets by an upper-case calligraphy font. ${\mathbf{1}}_L$ and ${\mathbf{1}}_{L\times N}$ refer to an $L$-dimensional vector and a $L\times N$ matrix both with unitary entries, respectively, and $\mathbf{I}$ denotes to the identity matrix. For vectors $\mathbf{x}$ and $\mathbf{y}$ with the same dimension, $\parallel \mathbf{x}\parallel$ describes its Euclidean norm, ${\mathbf{x}}^{\top }$ the transpose operation, $〈\mathbf{x},\mathbf{y}〉={\mathbf{x}}^{\top }\mathbf{y}$ the inner-product, and $\mathbf{x}\odot \mathbf{y}$ the Hadamard or element-wise product of both vectors. For a matrix $\mathbf{A}$, ${\parallel \mathbf{A}\parallel }_F=\text{Tr}\left({\mathbf{A}}^{\top }\mathbf{A}\right)$ refers to its Frobenius norm, where $\text{Tr}\left(\mathbf{A}\right)$ denotes the trace operation of $\mathbf{A}$, i.e., sum of diagonal elements in $\mathbf{A}$, and ${\lambda }_{\min }\left({\mathbf{A}}^{\top }\mathbf{A}\right)$ to the minimum eigenvalue of ${\mathbf{A}}^{\top }\mathbf{A}$. For a vector $\mathbf{x}$, $\text{diag}\left(\mathbf{x}\right)$ represents the squared matrix whose diagonal elements are in $\mathbf{x}$. For two sets $\mathcal{X}$ and $\mathcal{Y}$, the $K$-elements of $\mathcal{X}$ are described by $\mathcal{X}=\{x_1,\dots ,x_K\}$, $\text{card}\left(\mathcal{X}\right)=K$ denotes its cardinality, and $\mathcal{X}\setminus \mathcal{Y}$ the set difference, i.e., the elements in $\mathcal{X}$ that do not belong to $\mathcal{Y}$.