Deep Hyperspectral Unmixing using Transformer Network

Preetam Ghosh,
Swalpa Kumar Roy, Student Member, IEEE,
Bikram Koirala, Member, IEEE,
Behnood Rasti, Senior Member, IEEE,
and Paul Scheunders, Senior Member, IEEE

This work is supported by the Research Foundation-Flanders - Project G031921N.

P. Ghosh and S. K. Roy are with the Department of Computer Science and Engineering, Jalpaiguri Engineering College, West Bengal 735102, India (e-mail<EMAIL_ADDRESS>swalpa@cse.jgec.ac.in).

B. Koirala and P. Scheunders are with Imec-Visionlab, University of Antwerp (CDE) Universiteitsplein 1, B-2610 Antwerp, Belgium (e-mail<EMAIL_ADDRESS>paul.scheunders@uantwerpen.be).

B. Rasti is with Helmholtz-Zentrum Dresden-Rossendorf, Helmholtz Institute Freiberg for Resource Technology, Machine Learning Group, Chemnitzer Straße 40, 09599 Freiberg, Germany; (e-mail: b.rasti@hzdr.de).

Submission to arXiv
Ghosh et al.: Bare Demo of IEEEtran.cls for IEEE Journals

Transformers have intrigued the vision research community with their state-of-the-art performance in natural language processing. With their superior performance, transformers have found their way in the field of hyperspectral image classification and achieved promising results. In this article, we harness the power of transformers to conquer the task of hyperspectral unmixing and propose a novel deep unmixing model with transformers. We aim to utilize the ability of transformers to better capture the global feature dependencies in order to enhance the quality of the endmember spectra and the abundance maps. The proposed model is a combination of a convolutional autoencoder and a transformer. The hyperspectral data is encoded by the convolutional encoder. The transformer captures long-range dependencies between the representations derived from the encoder. The data are reconstructed using a convolutional decoder. We applied the proposed unmixing model to three widely used unmixing datasets, i.e., Samson, Apex, and Washington DC mall  and compared it with the state-of-the-art in terms of root mean squared error and spectral angle distance. The source code for the proposed model will be made publicly available at <https://github.com/preetam22n/DeepTrans-HSU>.

Hyperspectral image, unmixing, convolutional neural network, deep learning, transformer network, abundance map, endmember extraction, blind unmixing

§ INTRODUCTION

Advances in remote sensing technology improved environmental monitoring, e.g., for tracking rapid environmental changes and take precautionary actions. In particular, hyperspectral imaging (HSI) has attracted much attention in recent years. Its tasks include but are not limited to land used and land cover classification [1, 2, 3], forest applications [4, 5] and target detection [6] etc. In hyperspectral remote sensing, each spectral pixel might cover several pure materials on the ground due to its limited spatial resolution. The acquired spectral reflectance is then a mixture of the pure spectra (endmembers) of the materials within the pixel [7, 8]. Spectral unmixing techniques estimate the relative proportions (fractional abundances) of the endmembers within spectral pixels. The primary goal of spectral unmixing methods is to extract/estimate endmembers and their fractional abundances in each pixel by only utilizing the observed hyperspectral image. However, this often relies on the presence of a spectral library or the estimation/extraction of endmembers, i.e., pure spectral pixels that span the abundance subspace.

In remote sensing applications, it is generally  assumed that the spectra of the pure materials are mixed linearly and several linear unmixing techniques have been developed [9]. When the endmembers of the hyperspectral image are available, the fractional abundances can be estimated by minimizing the least squared errors between the actual reflectance spectra and the ones, reconstructed by the linear model. To have a physical interpretation of the estimated fractional abundances, one must assume that no endmember can have a negative abundance. This constraint is often described as the abundance non-negativity constraint (ANC). The second constraint is the abundance sum-to-one constraint (ASC), i.e., the observed reflectance spectrum is completely composed of endmember contributions. The fully constrained least squares unmixing algorithm (FCLSU) [10] obeys both ANC and ASC. The hyperspectral pixels that follow the fully constrained linear mixing model lie on a linear simplex whose corners (vertices) are given by the endmembers. As a result, many endmember extraction algorithms have been proposed to maximize the volume enclosing simplex in the hyperspectral dataset [11, 12, 13, 14]. When endmembers are not available in the hyperspectral image (no pure pixel-scenario), virtual endmembers can be estimated by seeking the minimum volume linear simplex, which encloses the data points [15, 16].

Spectral unmixing techniques that can simultaneously estimate the endmembers and the abundances are referred to as blind unmixing techniques [17, 18, 19, 20, 21]. These methods formulate the unmixing problem as a nonconvex optimization problem with respect to both endmembers and abundances. A common practice is to induce a geometrical penalty term in the fully constrained least squares method. In [22], the Euclidean distances between the estimated endmembers and the center of the hyperspectral pixels were selected to form a geometrical penalty term. In [21], the Euclidean distances between the estimated endmembers and endmembers extracted by Vertex Component Analysis (VCA) were selected for the penalty term. The total variation (TV) of all estimated endmembers was considered in [23] as a geometrical penalty. The optimization equation of these methods contains a regularization parameter, which denotes the trade-off between the geometrical penalty term and the fidelity term. This parameter is data-dependent, and selecting a proper parameter for each hyperspectral image is a highly complex problem. To tackle this challenge, in [24] an automatic parameter selection technique was proposed.

Blind unmixing methods can accurately estimate endmembers, if sufficient hyperspectral pixels are available on the facets of the data simplex. When the spectral pixels are highly mixed, the estimated endmembers are not satisfactory, which leads to poor abundance maps. To deal with highly mixed scenarios, sparse unmixing techniques have been proposed [25, 26, 27]. Sparse unmixing utilizes a rich and well-designed library of pure spectra and applies sparse regression for the abundance estimation. A major challenge is to correct mismatches between the real reflectance spectra and the library spectra, caused by differences in the acquisition conditions of the two data types.

Due to the success of deep learning-based networks in machine learning and computer vision applications [28, 29], recently, a variety of deep neural networks has been proposed for hyperspectral unmixing. These networks are mainly based on variations of deep encoder-decoder networks. The inputs of these networks are the reflectance spectra, while the outputs are the reconstructed spectra. The encoder transforms the input spectra to the fractional abundances while the decoder transforms the abundances to the reconstructed spectra using linear layers, with the endmembers as the weights. In [30], autoencoders that have been used for hyperspectral unmixing are grouped into five different categories: (a) Sparse nonnegative autoencoders (a stack of nonnegative sparse autoencoders (SNSA)) [31] (b) Variational autoencoders (Deep AutoEncoder Network (DAEN) [32], Deep Generative Unmixing algorithm (DeepGUn) [33] (c) Adversarial autoencoders (Adversarial autoencoder network (AAENet)) [34, 35, 36] (d) Denoising autoencoders (an untied Denoising Autoencoder with Sparsity (uDAS)) [37], and (e) Convolutional autoencoders [38, 39, 40]. Although the advantage of incorporating the spatial information for hyperspectral unmixing has been demonstrated in the literature (especially for  homogeneous regions), in SNSA, DAEN, DeepGUn, and uDAS, the spatial information is ignored. Several convolutional autoencoder-based unmixing techniques have been proposed to effectively incorporate the spatial correlation between adjacent pixels. In [41], a supervised hyperspectral unmixing method (i.e., the endmembers are assumed to be known) was proposed using a 3D convolutional autoencoder. The method referred to as unmixing using deep image prior (UnDIP) [42] utilizes endmembers extracted by a simplex volume maximization (SiVM) technique. Although several deep learning-based unmixing techniques have been specifically designed for blind unmixing, most of the methods fail when pure pixels are not available in the hyperspectral image. This is because they do not exploit the geometrical properties of the linear simplex. Recently, a minimum simplex convolutional network (MiSiCNet) [43] was proposed to incorporate both the spatial correlation between adjacent pixels and the geometrical properties of the linear simplex.

§.§ Contributions and Novelties

HSI, being complex in nature, pose a big challenge for Convolutional Neural Networks (CNN). As a convolution operation is limited to local features determined by the dimension of the kernel size, a significant amount of contextual information present in the original HS image is lost. Most autoencoders (AEs) are purely based on CNN networks and therefore fail to preserve a substantial portion of the original information due to the limited dimensionality of the latent space. That poses an even more significant problem in the case of HSI unmixing because the final number of endmembers is considerably lower than the initial number of spectra, causing a lot of contextual information to be lost. To address this issue, a transformer [44, 45] will be utilized that can recover some of the lost information, owing to its ability to capture global contextual feature dependencies [46].
For this, the AE output is rearranged in terms of patches. Inspired by [47], we propose a new attention mechanism, called Multihead Self-Patch Attention to calculate the long-range dependencies between these patches.

This leads to better quality abundance maps and an overall better unmixing result, which in turn helps the decoder to better reconstruct the HSI. Since the weights of the decoder are used to obtain the endmember spectra, a better quality of extracted endmembers is obtained. The contribution of the proposed methodology to this end is summarized below:

* We propose a new unmixing method based on a combination of a convolutional autoencoder and a transformer. The transformer is applied to the latent space of the autoencoder to enhance the feature extraction and to ensure a better estimation of abundances and endmembers. For this, the AE output is rearranged into patches.

* Inside the transformer encoder, we propose a new attention mechanism which is referred to as Multihead Self-Patch Attention. The attention modules of the multi-head self-patch attention find the global contextual feature dependencies by determining the long-range relationship between the image patches.

* To estimate the endmembers, we apply a single convolution layer whose weights are initialized by VCA. The weights are learned and improved during the training of the model to obtain  endmember spectra of superior quality.

The remaining of the paper is organized as follows: Section <ref> introduces the components of the proposed method including the novel Multihead Self-Patch Attention for transformer based deep HS image unmixing. In Section <ref>, extensive experiments are conducted with three benchmark datasets, and a hyperparameter sensitivity analysis and discussions are provided. Finally, comprehensive conclusions are drawn in Section <ref>.

§ PROPOSED METHODOLOGY

Graphical representation of the proposed deep hyperspectral unmixing model.

Let the HSI of spatial dimensions $H \times W$ with $B$ spectral bands be denoted by $\mathbf{I}\in \mathbb{R}^{B\times{H}\times{W}}$. The HSI can be reshaped to produce the matrix $\mathbf{Y} \in  \mathbb{R}^{B\times n}$, where $n = {H}\cdot{W}$ is the number of hyperspectral pixels. The endmember matrix will be denoted by $\mathbf{E} \in \mathbb{R}^{B \times R}$ where $R$ represents the number of endmembers present in the HSI. The corresponding abundance cube (i.e., the stack of $R$ abundance maps) is represented by $\mathbf{M}\in \mathbb{R}^{R\times{H}\times{W}}$. The abundance cube can be reshaped to produce the matrix $\mathbf{A} \in \mathbb{R}^{R\times{n}}$.

§.§ Problem formulation

In the Linear Mixing Model (LMM), the observed spectral reflectance is formulated as:

\begin{equation}
\mathbf{Y = EA + N}
\label{equ:lmm}
\end{equation}

where $\mathbf{N} \in \mathbb{R}^{B\times n}$ is the additive noise present in $\mathbf{Y}$.
Generally, three physical constraints should be satisfied: 1) the endmember matrix should be non-negative $\mathbf{E} \geq 0$; 2) ANC (Eq. (<ref>));  and 3) ASC  (Eq. (<ref>)):

\begin{equation}
\mathbf{A} \geq 0
\label{equ:ANC}
\end{equation}

\begin{equation}
{\bf 1}_{R}^{T}{\bf A}={\bf 1}_{n}^{T}
\label{equ:ASC}
\end{equation}

where ${\bf 1}_n$ indicates an $n$-component column vector of ones.

Since spectal unmixing is a reconstruction problem, in which abundance maps are reconstructed from the given HSI, AEs can be applied. AEs are quite capable at reconstructing and extracting information from the given inputs. In this work, the performance of an AE is complemented by the use of a transformer, to significantly improve the quality of the generated abundance maps and consequently the extracted spectral signatures of the endmembers. Fig. <ref> illustrates the proposed model for deep HSI unmixing. The components of the model
are discussed in detail in subsections <ref> through <ref>.

§.§ Hyperspectral feature extraction using AE

AEs encode the input into a latent space with a lower dimensionality, learning only the salient features within the input image while avoiding unnecessary details. Owing to CNNs ability to extract high-level abstract features, using them in the encoder part of an AE provides a twofold benefit. Firstly, it heavily reduces the large number of spectral bands of a HSI and secondly, it extracts discriminative high-level features that form the base for the transformer in the next step.

The CNN applied in the encoder block of the proposed model contains three layers. Each layer progressively reduces the number of spectral bands of the HSI until $C$ spectral bands remain. The value of $C$ is a hyperparameter to be set.

As the convolutional layer is primarily used to reduce the number of channels of the input HSI, a kernel size of $1\times{1}$ is used to keep the number of parameters low and to facilitate a faster training of the model. All three layers use a 2D convolution operation followed by a batch normalization (BN). To mitigate the vanishing gradient problem of the network, the first layer uses a dropout function. To introduce non-linearity, Leaky ReLU is used in the output of the first two layers of the AE. Table <ref> summarizes the structure of the encoder.

Layerwise summary of the Encoder block where $B$ represents the number of spectral bands and $C$ is the number of output bands.

Layers     Composition     Kernel     Bands in      Bands out
Layer 1     c]@c@Conv 2D
Leaky ReLU     3*$(1\times{1})$     B     128
1-2 4-5
Layer 2     c]@c@Conv 2D
Leaky ReLU          128     64
1-2 4-5
Layers 3     c]@c@Conv 2D
BN         64     C

In the encoder, the HSI $\mathbf{I} \in \mathbb{R}^{B \times H \times W}$ is transformed by the three consecutive layers of the encoder block into $\mathbf{I'} \in \mathbb{R}^{H \times W \times C}$:

\begin{equation}
\begin{aligned}
\mathbf{I_{1}} &= f_1(\mathbf{W_1I} + \mathbf{U_1}) \\
\mathbf{I_{2}} &= f_2(\mathbf{W_2I_1} + \mathbf{U_2}) \\
\mathbf{I_{3}} &= f_3(\mathbf{W_3I_2} + \mathbf{U_3}) \\
\mathbf{I'} &= \mathbf{I_{3}^T}
\end{aligned}
\end{equation}

where $f_1(\cdot), ~f_2(\cdot)$ and $f_3(\cdot)$ denote the three encoder layers and $\mathbf{W_1,W_2,W_3}$ and $\mathbf{U_1,U_2,U_3}$ are the weights and biases, respectively of each layer. The superscript T denotes the matrix transpose operation.

§.§ Patch and Position Embeddings

To efficiently capture the long range feature dependencies, the AE output is rearranged in terms of patches.
The output of the AE encoder is the cube $\mathbf{I'}$ of dimension (${H}\times{W}\times{C}$) where $H,W$ are the spatial dimensions and $C$ represents the reduced number of bands of the output. These features are grouped in patches $({({m}\cdot{p})}\times{({n}\cdot{p})}\times{C})$ where $p$ is the patch size and $m \cdot n$ is the total number of patches. Then the cube is reshaped to a matrix $\mathbf{X_{patch}}$ of size $({({m}\cdot{n})}\times{({p}\cdot{p}\cdot{C})})$ = ($N' \times D$) where $N'$ is the total number of patches and $D$ is the dimension of each patch embedding. As an example, for the Samson dataset (Section <ref>),  with $p=5$ and $C=24$, the rearrangement is given as:

\begin{equation}
\begin{aligned}
\mathbf{ I'}&= ({95}\times{95}\times{24})\\
&= ({({19}\cdot{5})}\times{({19}\cdot{5})}\times{24})\\
& \to \\
\mathbf{X_{patch}} &=  ({({19}\cdot{19})}\times{({5}\cdot{5}\cdot{24})})\\
&= (361 \times 600)
\end{aligned}
\nonumber
\end{equation}

In a next step, learnable class tokens $\mathbf{X_{cls}}$ of dimensions ($1 \times D$) are defined, in which the transformer encoder will capture the long range semantic information of the patch tokens.
Moreover, positional tokens $\mathbf{X_{pos}}$ of shape $(N \times D)$, with $N=N'+1$ are generated to retain patch positional information. Rather than providing pixel and patch positional information, the positional tokens will be learned by the transformer encoder as well. Both are randomly initialized.

$\mathbf{X_{cls}}$ is appended as an extra row to the matrix $\mathbf{X_{patch}}$ and  $\mathbf{X_{pos}}$  is added to the feature embedding:

\begin{equation}
\mathbf{X' = (X_{cls} \parallel X_{patch}) + X_{pos} \\
= (X'_{cls} \parallel X'_{patch})}
\label{equ:patch_embed}
\end{equation}

with $\parallel$ the concatenation operation.

Transformer Encoder with Multihead Self-Patch Attention.

§.§ Transformer Encoder with Multihead Self-Patch Attention

$\mathbf{X'}$ is the input of the next phase, which is composed of one or several transformer encoders. Each transformer encoder contains a Multihead Self-Patch Attention network [44]. The goal of this network is the exchange of information within the patch tokens to capture their long range contextual information and to feed this into the class token. To preserve the overall patch structure, the patch tokens are appended again to the learned class token.

Fig. <ref> depicts the proposed Multihead Self-Patch Attention network. A detailed description of each step is given below.

Step 1: In a fist step, the overall patch matrix $\mathbf{X}'$ enters the self attention block of the transformer after going through a layer normalisation step. Attention is calculated by three linear layers. One layer works on the class token only (weight $\mathbf{W_q}$ and output $\mathbf{q}$ of size $(1\times{D})$). The other 2 layers work on the entire patch matrix (weights $\mathbf{W_k}$ and  $\mathbf{W_v}$ and outputs $\mathbf{k}$ and $\mathbf{v}$, both of size $(N\times{D})$):

\begin{equation}
\mathbf{q=W_qX'_{cls},\hskip 10pt k=W_kX',\hskip 10pt v=W_vX'}
\nonumber
\end{equation}

Step 2: In the next step, the attention weight ($\mathbf{A}$) is calculated by computing the pairwise similarity between $\mathbf{q}$ and $\mathbf{k}$ and applying a softmax function:

\begin{equation}
\mathbf{A} = \texttt{softmax}\mathbf{(qk^T/\sqrt{D})}
\nonumber
\end{equation}

The scaling term ($1/\sqrt{D}$) counteracts the small gradients of the softmax function. The self-patch attention () is then computed as:

\begin{equation}
\texttt{PA}(\mathbf{X'}) = \mathbf{Av}
\end{equation}

To further enhance the relationships among the different patches, self-patch attention with multiple heads is applied.
For this, $\mathbf{q}$, $\mathbf{k}$ , and $\mathbf{v}$ have to reshape into matrices $\mathbf{q'}$, $\mathbf{k'}$ , and $\mathbf{v'}$ of size $(h_n\times{D/h_n})$, $((N\cdot h_n)\times{D/h_n})$, $((N\cdot h_n)\times{D/h_n})$ respectively, where $h_n$ denotes the number of heads (attention modules).  Then, the attention weight becomes:

\begin{equation}
\mathbf{A'} = \texttt{softmax}\mathbf{(q'k'^T/\sqrt{D/h_n})}
\nonumber
\end{equation}

The self-patch attention with multiple heads () is then computed as:

\begin{equation}
\texttt{MPA}(\mathbf{X'}) = \mathbf{A'v'}
\label{equ:MPA}
\end{equation}

Step 3: The output of  is a matrix of size $(h_n\times{D/h_n})$, and is then reshaped back to a matrix of size $(1\times{D})$. This matrix is further passed through a linear layer (weights $\mathbf{W_l} \in \mathbb{R}^{D \times D}$) and added up with the original class token $\mathbf{X'_{cls}}$ to obtain the class token $\mathbf{y_{cls}}$:

\begin{equation}
\mathbf{y_{cls}} = \texttt{MPA}(\mathbf{X'})\mathbf{W_l}+\mathbf{X'_{cls}}
\label{equ:token_out}
\end{equation}

Step 4: Finally, $\mathbf{y_{cls}}$ is concatenated with the layer normalised patch tokens to obtain the output of the attention network $\mathbf{X}''$:

\begin{equation}
\mathbf{X''} = \mathbf{y_{cls}} \parallel \texttt{LN}(\mathbf{X'_{patch}})
\label{equ:attn_out}
\end{equation}

As the output of the Multihead Self-Patch Attention network, the feature embedding $\mathbf{X}''$ is passed through a normalization layer and then fed into an Multi Layered Perceptron () block along with a residual connection to obtain the final output of the transformer encoder block (see bottom right of Fig. <ref>):

\begin{equation}
\mathbf{X'''} = \mathbf{X''} + \texttt{MLP}(\texttt{LN}(\mathbf{X''}))
\label{equ:trans_out}
\end{equation}

Any number of such transformer encoders can be applied sequentially. In this work, two encoders have been applied. The output of the final block is used for further processing down the line.

The pseudo code of the Transformer Encoder with Multihead Self-Patch Attention, is shown in Algorithm <ref>.

$\mathbf{X'}$, $\mathbf{X'_{cls}}$, $\mathbf{X'_{patch}}$, $D$, $h_n$

Multihead Self-Patch Attention (Begin)

10pt Step 1. $\mathbf{q=W_qX'_{cls},\hskip 10pt k=W_kX',\hskip 10pt v=W_vX'}$,

42pt $ \mathbf{q' = \textup{reshape}(q), k' = \textup{reshape}(k)}$,

42pt $ \mathbf{v' =\textup{reshape}(v)}$

10pt Step 2. $\mathbf{A'} = \texttt{softmax}\mathbf{(q'k'^T/\sqrt{D/h_n})}$,

45pt      $\texttt{MPA}(\mathbf{X'}) = \mathbf{A'v'}$ (<ref>)

Multihead Self-Patch Attention (End)

Step 3. $\mathbf{y_{cls}} = \textup{reshape}(\texttt{MPA}(\mathbf{X'}))\mathbf{W_l}+\mathbf{X'_{cls}}$ (<ref>)

Step 4. $\mathbf{X''} = \mathbf{y_{cls}} \parallel \texttt{LN}(\mathbf{X'_{patch}})$ (<ref>),

32pt      $\mathbf{X'''} = \mathbf{X''} + \texttt{MLP}(\texttt{LN}(\mathbf{X''}))$ (<ref>),

32pt      $\mathbf{X'''_{cls}} = \mathbf{X'''}(1,:)$

Transformer Encoder with Multihead Self-Patch Attention

§.§ Unmixing with decoder

The transformer produces the results $\mathbf{X'''}\in \mathbb{R}^{N\times D}$, where $N$ is the total number of tokens and $D$ is the dimension of each token. However, for the purpose of unmixing, only the class token $\mathbf{X'''_{cls}}$ (i.e., the first row of $\mathbf{X'''}$) of size $(1 \times D)$ is considered and forwarded to the upsampling block.  To do so, we reshape $\mathbf{X'''_{cls}}$ to a matrix of size $R \times (D/R)$, and then upscale it to size $R\times (H \cdot W)$. Upscaling from a relatively small dimension of $D/R$ to the dimensions $H\cdot W$ introduces noise in the final output. To solve this issue, a convolution operation with parameters $kernel\_size=(3 \times 3),~stride=1,~padding=1$ is used. Finally, a reshaping operation is carried out to convert the output to the shape of the adundance cube $\mathbf{M}$ i.e., $(R \times H \times W)$. To ensure that the ANC and ASC constraints (Eqs. (<ref>) and (<ref>)) are satisfied, a softmax layer is used along the $R$ dimension.

To calculate the endmembers, the abundance matrix $\mathbf{M}$ is passed through the decoder block of the AE which consists of a single convolutional layer. This convolution operation increases the number of bands in $\mathbf{M}$ from $R$ to $B$, to obtain the reconstructed HSI $\widehat{\mathbf{I}}$. The weights of the convolution layer, which are initialized with the endmembers obtained from VCA, are optimized to estimate the final endmembers $\mathbf{\hat{E} \in \mathbb{R}^{B \times R}}$.

§.§ Losses and Optimization functions

In order to train the proposed model, a combination of two different losses: Reconstruction Error (RE) loss and Spectral Angle Distance (SAD) loss were applied:

\begin{equation}
L_{RE}({\bf I},\hat{\bf I}) = \frac{1}{H\cdot W}\sum\limits_{i=1}^H\sum\limits_{j=1}^W({\mathbf{\hat{I}_{ij}} - \mathbf{I_{ij}}})^2
\label{equ:mse_loss}
\end{equation}

\begin{equation}
L_{SAD}({\bf I},\hat{\bf I})=\frac{1}{R}\sum_{i=1}^R\arccos{\left(\frac{ \left\langle{\bf I}_{i}, \hat{\bf I}_{i}\right\rangle}{\|{\bf I}_{i}\|_2\|{\bf \hat I}_{i}\|_2}\right)}
\label{equ:sad_loss}
\end{equation}

The RE loss is calculated by the Mean Squared Error (MSE) objective function and helps the encoder part to learn only the essential features of the input HSI while discarding non-essential details. The SAD loss is a scale invariant objective function. MSE discriminates between  endmembers, based on their absolute magnitude which is not desirable in case of HSI unmixing. Including SAD loss helps to counter this drawback of the MSE objective function and makes the overall model converge much faster. The total loss is calculated as the weighted sum of these two losses:

\begin{equation}
L = \beta L_{RE} + \gamma L_{SAD}
\label{equ:total_loss}
\end{equation}

with regularization parameters $\beta$ and $\gamma$.

§ EXPERIMENTAL RESULTS

§.§ Hyperspectral Data Description

We performed experiments on three datasets. The description of the datasets are given below.

§.§.§ Samson

(a)     (b)

Samson image: (a) True-color image (Red: 571.01 nm, Green: 539.53 nm, and Blue: 432.48 nm) (b) Endmembers.

The Samson hyperspectral dataset ([48]) (Fig. <ref>(a)) utilized in this work contains 95$\times$95 hyperspectral pixels. Each hyperspectral pixel contains reflection values from 156 bands covering the wavelength range [401–889] nm. In this hyperspectral image, there are three endmembers (i.e., Soil, Tree, and Water). The ground truth endmember spectra (see Fig. <ref>(b)) were manually selected from the image, and ground truth abundance maps were produced by applying FCLSU.

§.§.§ Apex

(a)     (b)

Apex image: (a) True-color image (Red: 572.2 nm, Green: 532.3 nm, Blue: 426.5 nm); (b) Endmembers.

Fig. <ref>(a) shows a cropped image of the Apex dataset ([49]), as used in this work. This image contains 110$\times$110 hyperspectral pixels. Each hyperspectral pixel contains reflection values from 285 bands covering the wavelength range [413–2420] nm. There are four endmembers (i.e., Water, Tree, Road, and Roof) in this hyperspectral image. The ground truth endmember spectra (see Fig. <ref>(b)) were manually selected from the image, and ground truth abundance maps were produced by applying FCLSU.

§.§.§ Washington DC Mall

(a)     (b)

Washington DC Mall image: (a) True-color image (Red: 572.7 nm, Green: 530.1 nm, Blue: 425.0 nm);  (b) Endmembers.

This hyperspectral image is acquired over the Washington DC Mall using the HYDICE sensor [https://engineering.purdue.edu/ biehl/MultiSpec/hyperspectral.html]. Fig. <ref>(a) shows the cropped data used in this paper that contains 290 $\times$ 290 pixels. Each hyperspectral pixel contains reflection values from 191 bands covering the wavelength range [400–2400] nm. There are six endmembers (i.e., Grass, Tree, Roof, Road, Water, and Trail) in this hyperspectral image. The ground truth endmember spectra (Fig. <ref>(b)) were manually selected from the image, and ground truth abundance maps were produced by applying FCLSU.

§.§ Experimental Setup

The performance of the proposed model is evaluated and compared to six different unmixing techniques from different categories: Geometrical unmixing method  [10] using VCA [12] for endmember extraction, Geometrical and blind unmixing method  [24], Sparse unmixing method Collaborative LASSO () [50] and Deep unmixing methods  [37],  [42], and  [40].

§.§ Hyperparameters

In deep unmixing models, the produced results are typically largely dependent on the hyperparameter settings. Choosing proper values for the hyperparameters can significantly improve results. Table <ref> shows the hyperparameters used for training the proposed model, which are further discussed below:

Samson dataset: The patch size $p$ was selected to be $(5 \times 5)$

and the transformer input dimensionality $C$ was chosen to be $24$.
The values 4e3 and 5e-3 were used for the regularization parameters $\beta$ and $\gamma$ respectively. The model was trained during $200$ epochs with an initial learning rate of 6e-3, which was reduced by $20\%$ after every $15$ epochs. A weight decay rate of 4e-5 was incorporated in the optimization function to keep the losses in check.

Apex dataset: The patch size $p$ was selected to be $(5 \times 5)$

and the transformer input dimensionality $C$ was chosen to be $32$. Values of  4e3 and 5e-2 were used for the regularization parameters $\beta$ and $\gamma$ respectively. The model was trained during $200$ epochs with an initial learning rate of 9e-3, which was reduced by $20\%$ after every $15$ epochs. A weight decay rate of 4e-5 was incorporated in the optimization function to keep the losses in check.

Washington DC Mall dataset: The patch size $p$ was selected to be $(10 \times 10)$

and the transformer input dimensionality $C$ was chosen to be $24$. Values 5e3 and 2e-3 were used for the regularization parameters $\beta$ and $\gamma$ respectively. The model was trained during $150$ epochs with an initial learning rate of 6e-3, which was reduced by $20\%$ after every $15$ epochs. A weight decay rate of 3e-5 was incorporated in the optimization function to keep the losses in check.

Hyperparameters used for training the proposed model.

Hyperparameters     Samson     Apex     WDC Mall
$p$     $(5\times{5})$     $(5\times{5})$     $(10\times{10})$

$C$     24     32     24
$\beta$     5e3     5e3     5e3

$\gamma$     3e-2     5e-2     1e-4
Epoch     200     200     150

Learning rate     6e-3     9e-3     6e-3

Weight decay     4e-5     4e-5     3e-5

§.§ Quantitative Performance Measures

Quantitative results are provided by the root mean squared error (RMSE) between the estimated and ground truth abundance fractions:

\begin{equation}
\textnormal{RMSE}({\bf M},\hat{\bf M}) =\sqrt{\frac{1}{RHW}\sum_{k=1}^R\sum_{i=1}^H\sum_{j=1}^W\left({\hat{\bf M}}_{kij}-{\bf M}_{kij}\right)^2}
\end{equation}

and the spectral angle distance (SAD) in degree between the estimated and ground truth endmembers:

\begin{equation}
\mbox{SAD}({\bf S},\hat{\bf S})=\frac{1}{R}\sum_{i=1}^R\arccos{\left(\frac{ \left\langle{\bf s}_{(i)}, \hat{\bf s}_{(i)}\right\rangle}{\left\|{\bf s}_{(i)}\right\|_2\left\|\hat{\bf s}_{(i)}\right\|_2}\right)},
\end{equation}

where $\langle . \rangle$ denotes the inner product  and ${\bf s}_{(i)}$ indicates the $i$th column of the ground truth endmembers matrix ${\bf S}$.

§.§ Unmixing Experiments: Quantitative Results

Samson Dataset: Quantitative results on the Samson dataset can be found in Tables <ref> and <ref>. The results confirm that the proposed model outperforms the other techniques in terms of both abundance and endmember estimation with a mean RMSE of $0.0783$ showing a $48.02\%$ improvement to the next best method and a mean SAD value of $0.0608$ which amounts to a $35.93\%$ improvement.

-2 pt
RMSE (Samson Dataset). The best performances are shown in bold.

CyCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Soil     0.2417     0.1506     0.1766     0.2011     0.1778     0.1799     0.0712

Tree     0.1386     0.0607     0.0653     0.1466     0.1330     0.1383     0.0683

Water     0.2654     0.1181     0.1492     0.2063     0.2096     0.2303     0.0930

Overall     0.2222     0.1159     0.1387     0.1866     0.1763     0.1867     0.0783

-2 pt
SAD (Samson Dataset). The best performances are shown in bold.

CyCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Soil     0.1144     0.0155     0.0391     0.0259     0.0259     0.0358     0.0128

Tree     0.1517     0.0832     0.1239     0.0748     0.0961     0.0960     0.0674

Water     0.2081     0.1402     1.5201     0.1554     0.1554     0.1527     0.0729

Overall     0.1581     0.0796     0.5610     0.0854     0.0925     0.0948     0.0510

GT     CyCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Samson dataset - Visual comparison of the abundance maps obtained by the different unmixing techniques.

CyCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Samson dataset - Visual comparison of the endmembers obtained by the different unmixing techniques. Blue: ground truth endmembers;  Orange: estimated endmembers.

Apex Dataset: Quantitative results on the Apex dataset can be found in Tables <ref> and <ref>. The endmember “Road" of the Apex dataset is found to be quite a challenge for the other methods, while the proposed method estimates this endmember satisfactorily. The proposed method outperforms the other unmixing techniques with a mean RMSE value of $0.1264$ and a mean SAD value of $0.0867$.  Additionally, it provides the best endmember estimation for Road and Water in terms of SAD.

-2 pt
RMSE (Apex Dataset). The best performances are shown in bold.

CyCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Road     0.2921     0.3078     0.2331     0.1806     0.1737     0.1973     0.1776

Tree     0.2020     0.1907     0.0944     0.2468     0.2154     0.1419     0.0993

Roof     0.1630     0.1483     0.1201     0.2359     0.2554     0.2303     0.1200

Water     0.1213     0.0797     0.1327     0.3751     0.4170     0.2887     0.0902

Overall     0.2046     0.1997     0.1543     0.2692     0.2809     0.2210     0.1264

-2 pt
SAD (Apex Dataset). The best performances are shown in bold.

CyCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Road     0.4543     0.6772     0.4003     0.0907     0.6915     0.4551     0.0836

Tree     0.0850     0.2063     0.2710     0.1339     0.2644     0.1405     0.1295

Roof     0.1298     0.1002     0.1753     0.0689     0.1471     0.0860     0.0903

Water     0.6223     0.5137     1.8417     0.5040     0.5176     0.2251     0.0434

Overall     0.3228     0.3744     0.6721     0.1994     0.4052     0.2267     0.0867

GT     CyCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Apex dataset - Visual comparison of the abundance maps obtained by the different unmixing techniques.

CYCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Apex dataset - Visual comparison of the endmembers obtained by the different unmixing techniques. Blue: ground truth endmembers;  Orange: estimated endmembers.

Washington DC Mall Dataset: Quantitative results on the Washington DC Mall dataset can be found in Tables <ref> and <ref>. Among all the considered datasets, the similarity between the spectral signatures of its six endmembers provides the greatest challenge. The “Tree" and “Grass" endmembers have almost identical spectral signatures, and most methods struggle to find the difference. The proposed model successfully separated these two endmembers due to its ability to find long-range dependencies among the image patches, thus leading to a RMSE value of $0.1661$ and $0.0963$ for the “Grass" and “Tree" endmembers, respectively. In terms of overall RMSE and SAD, the proposed model outperforms the closest competitor method by $43.71\%$ and $52.11\%$ respectively.

-2 pt
RMSE (Washington DC Mall Dataset). The best performances are shown in bold.

CYCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Grass     0.4104     0.2901     0.3090     0.3624     0.2978     0.3780     0.1661

Tree     0.2824     0.4167     0.4025     0.2761     0.3514     0.3351     0.0963

Road     0.2545     0.2263     0.1757     0.2351     0.2436     0.2497     0.1353

Roof     0.4157     0.0437     0.0380     0.0862     0.0493     0.0463     0.0863

Water     0.3957     0.3102     0.2921     0.2076     0.3812     0.5156     0.1326

Trail     0.2072     0.1875     0.1230     0.1011     0.2360     0.1769     0.1492

Overall     0.3379     0.2715     0.2550     0.2322     0.2814     0.3206     0.1307

-2 pt
SAD (Washington DC Mall Dataset). The best performances are shown in bold.

CYCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Grass     0.0895     0.3171     0.1952     0.1851     0.3170     0.1897     0.2379

Tree     0.2704     0.3335     0.4507     0.7258     0.2883     0.4251     0.1225

Road     0.4642     0.3439     0.2243     0.8608     0.2316     0.6585     0.0781

Roof     0.9500     0.0331     0.2078     0.2826     0.0343     0.1992     0.3352

Water     0.4205     0.0305     0.6736     0.9495     0.7766     0.2328     0.0533

Trail     0.7906     0.3446     0.0615     0.1754     0.6472     0.0940     0.0951

Overall     0.4975     0.2338     0.3022     0.5299     0.3825     0.2999     0.1537

GT     CYCU     Collab     FCLSU     NMF     UnDIP     uDAS     Proposed

Washington DC Mall dataset - Visual comparison of the abundance maps obtained by the different unmixing techniques.

Overall observations:
From Tables <ref>, <ref> and <ref>, one can conclude that the overall performance of the proposed method beats the other competing methods by a significant margin in terms of RMSE.  , , and  also shows decent performance on the  Samson, Apex and WDC Mall datasets, but their performance is not consistent across the datasets.  and  were unable to beat any of the methods for any given class; however, their performance was consistent throughout the different datasets used in the experiments.  produced mixed results within a given dataset, with good performance on particular endmembers and significantly worse on other endmembers.

Tables <ref>, <ref> and <ref> make it clear that obtaining good spectral signatures for the endmembers is more difficult than producing a good abundance map. The proposed model considerably outperforms all the other competing methods. On the  Apex and WDC Mall datasets, the proposed model obtains SAD values of $0.0867$ and $0.1537$, respectively, about half of the next best method.

It is worth mentioning that a good SAD value does not necessarily guarantee good abundance maps, because SAD removes the norm of the endmember spectra. In other words, it ignores endmember scaling factors, caused by multiple reflections of the light and continuously variable illumination conditions in practical situations. However, such scaling factors can considerably affect the abundance estimation. As the proposed method provides the best results in both SAD and RMSE, one can conclude that it overcomes the mentioned problem at least to some degree.

CYCU     Collab     NMF     SiVM     VCA     uDAS     Proposed

Washington DC Mall dataset - Visual comparison of the endmembers obtained by the different unmixing techniques. Blue: ground truth endmembers;  Orange: estimated endmembers.

§.§ Visual Analysis of Abundance Maps and Endmembers

The abundance maps and spectral signatures of the endmembers provide a way to visually compare the generated results of the different unmixing methods. Figs. <ref>, <ref>, and <ref> show the abundance maps obtained from the various competing methods. It can be inferred that the abundance maps obtained from the proposed method are visually most similar to the ground truth abundance maps. The methods  and  fail to properly represent the endmember “Water" across all the experimental datasets. Decent results are obtained by the methods , , and , but their performance suffers from inconsistencies in RMSE values from one endmember to another. This causes the models to lose in terms of overall performance, even if they manage to obtain  good results on a particular endmember. For example, none of the competing methods was able to correctly produce the “Road" endmember in the Apex dataset. The success of the proposed model on this endmember can be attributed to the ability of the transformer encoder block with self-patch attention to find the long distance feature dependencies, which are otherwise lacking in the abundance maps obtained from the output of the convolutional network.

Figs. <ref>, <ref> and <ref> depict the extracted endmembers.  It was observed that the methods using  as initialization could not further improve the  extracted endmembers by much, leading to higher values of SAD later on. The proposed method however is also initialized by  analysis, but modifies the spectral signatures in a way that they more closely resemble the ground truth endmembers, with much lower SAD errors.

(a)     (b)     (c)

Effect of the hyperparameter $\gamma$ on RMSE and SAD error for (a) Samson dataset (b) Apex dataset and (c) Washington DC Mall dataset

§.§ Sensitivity Analysis to Hyperparameters

The hyperparameters $\beta$ and $\gamma$ play essential roles in determining the model's overall performance. In order to keep the training process simple, the value of $\beta$ was kept constant at 5e3 for all the datasets. Fig. <ref> depicts the sensitivity of the proposed unmixing model to the hyperparameter $\gamma$. Both SAD and RMSE values are correlated, and changing  $\gamma$ affects both of them similarly in most cases. The figure suggests that $\gamma$ can be set in the range 1e-4 to 1e-2, with a higher number of endmembers favouring a lower $\gamma$ value.

(a)     (b)     (c)

Effect of learning rate and weight decay on RMSE and SAD values for (a) Samson dataset (b) Apex dataset and (c) Washington DC Mall dataset

Apart from the hyperparameters mentioned above, the learning rate and the weight decay were also found to have a significant impact on the obtained results, as can be seen in Fig. <ref>. Learning rates were tested in the range from $0.001$ to $0.009$, and the best results were obtained in the  range from $0.006$ to $0.009$,  with images having lower spatial dimensions preferring a slightly lower learning rate. The weight decay was tested in the range from $\num{1e-5}$ to $\num{9e-5}$. Fig. <ref> suggests an optimal value  around $\num{3e-5}$. It was observed that the quality of the abundance maps quickly deteriorates with increasing weight decay.

The optimal parameters were selected using a grid search-based approach on the sample space [51], and the combination of parameter values which resulted in the minimal value of the loss function in Eq. (<ref>) was finally applied to obtain the reported results.

§ CONCLUSION

In this article we proposed  a novel HSI unmixing model that uses a convolutional autoencoder combined with a transformer. We demonstrated the viability of the novel Multihead Self-Patch Attention mechanism used in the encoder block of the transformer. The experiments were carried out on three real datasets, each with its unique set of challenges, and were successfully handled by the proposed model with consistent performance across the range of endmembers. The accuracy and consistency of the proposed model can be credited to the use of the transformer block which captures the long range feature dependencies that are otherwise not reachable by a CNN based architecture. This enables our model to achieve superior unmixing results, which are significantly better than the competing methods.

§ ACKNOWLEDGMENT

The authors thank Ganesan Narayanasamy who is leading IBM OpenPOWER/POWER enablement and ecosystem worldwide for his support to get the IBM AC922 system's access. B. Koirala is funded by the Research Foundation-Flanders - Project G031921N.

[1]
S. K. Roy, G. Krishna, S. R. Dubey, and B. B. Chaudhuri, “HybridSN:
Exploring 3-D-2-D CNN feature hierarchy for hyperspectral image
classification,” IEEE Geoscience and Remote Sensing Letters, vol. 17,
no. 2, pp. 277–281, 2020.

[2]
B. Rasti, D. Hong, R. Hang, P. Ghamisi, X. Kang, J. Chanussot, and J. A.
Benediktsson, “Feature extraction for hyperspectral imagery: The evolution
from shallow to deep: Overview and toolbox,” IEEE Geoscience and
Remote Sensing Magazine, vol. 8, no. 4, pp. 60–88, 2020.

[3]
S. K. Roy, R. Mondal, M. E. Paoletti, J. M. Haut, and A. Plaza, “Morphological
convolutional neural networks for hyperspectral image classification,”
IEEE Journal of Selected Topics in Applied Earth Observations and
Remote Sensing, vol. 14, pp. 8689–8702, 2021.

[4]
B. Koetz, F. Morsdorf, S. Van der Linden, T. Curt, and B. Allgöwer,
“Multi-source land cover classification for forest fire management based on
imaging spectrometry and lidar data,” Forest Ecology and Management,
vol. 256, no. 3, pp. 263–271, 2008.

[5]
M. Ahmad, S. Shabbir, S. K. Roy, D. Hong, X. Wu, J. Yao, A. M. Khan,
M. Mazzara, S. Distefano, and J. Chanussot, “Hyperspectral image
classification—traditional to deep models: A survey for future prospects,”
IEEE Journal of Selected Topics in Applied Earth Observations and
Remote Sensing, vol. 15, pp. 968–999, 2022.

[6]
W. Li, Q. Du, and B. Zhang, “Combined sparse and collaborative representation
for hyperspectral target detection,” Pattern Recognition, vol. 48,
no. 12, pp. 3904–3916, 2015.

[7]
J. M. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders,
N. Nasrabadi, and J. Chanussot, “Hyperspectral remote sensing data
analysis and future challenges,” IEEE Geoscience and Remote Sensing
Magazine, vol. 1, no. 2, pp. 6–36, 2013.

[8]
P. Ghamisi, N. Yokoya, J. Li, W. Liao, S. Liu, J. Plaza,
B. Rasti, and A. Plaza, “Advances in hyperspectral image and signal
processing: A comprehensive overview of the state of the art,” IEEE
Geoscience and Remote Sensing Magazine, vol. 5, no. 4, pp. 37–78, 2017.

[9]
J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du,
P. Gader, and J. Chanussot, “Hyperspectral unmixing overview:
Geometrical, statistical, and sparse regression-based approaches,”
IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 5, no. 2,
pp. 354–379, April 2012.

[10]
D. C. Heinz and Chein-I-Chang, “Fully constrained least squares linear
spectral mixture analysis method for material quantification in hyperspectral
imagery,” IEEE Transactions on Geoscience and Remote Sensing,
vol. 39, no. 3, pp. 529–545, 2001.

[11]
R. Heylen, D. Burazerovic, and P. Scheunders, “Fully constrained least
squares spectral unmixing by simplex projection,” IEEE Transactions on
Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4112–4122, Nov 2011.

[12]
J. Nascimento and J. Bioucas-Dias, “Vertex component analysis:
A~fast algorithm to extract endmembers spectra from
hyperspectral data,” in Pattern Recognition and Image Analysis, F. J.
Perales, A. J. C. Campilho, N. P. de la Blanca, and A. Sanfeliu, Eds.1em plus 0.5em minus 0.4emBerlin, Heidelberg: Springer Berlin
Heidelberg, 2003, pp. 626–635.

[13]
T.-H. Chan, C.-Y. Chi, Y.-M. Huang, and W.-K. Ma, “A convex analysis-based
minimum-volume enclosing simplex algorithm for hyperspectral unmixing,”
IEEE Transactions on Signal Processing, vol. 57, no. 11, pp.
4418–4432, 2009.

[14]
M. E. Winter, “N-FINDR: an algorithm for fast autonomous spectral end-member
determination in hyperspectral data,” in Imaging Spectrometry V,
M. R. Descour and S. S. Shen, Eds., vol. 3753, International Society for
Optics and Photonics.1em plus 0.5em minus 0.4emSPIE, 1999, pp.
266 – 275. [Online]. Available: <https://doi.org/10.1117/12.366289>

[15]
W. Full, R. Ehrlich, and J. Klovan, “Extended qmodel—objective definition of
external end members in the analysis of mixtures,” Journal of the
International Association for Mathematical Geology, vol. 13, pp. 331–344,
08 1981.

[16]
M. Craig, “Minimum-volume transforms for remotely sensed data,” IEEE
Transactions on Geoscience and Remote Sensing, vol. 32, no. 3, pp. 542–552,

[17]
J. Li and J. M. Bioucas-Dias, “Minimum volume simplex analysis: A fast
algorithm to unmix hyperspectral data,” in IGARSS 2008 - 2008 IEEE
International Geoscience and Remote Sensing Symposium, vol. 3, 2008, pp. III
– 250–III – 253.

[18]
J. M. Bioucas-Dias, “A variable splitting augmented lagrangian approach to
linear spectral unmixing,” in 2009 First Workshop on Hyperspectral
Image and Signal Processing: Evolution in Remote Sensing, 2009, pp. 1–4.

[19]
N. Dobigeon, S. Moussaoui, M. Coulon, J. Tourneret, and A. O. Hero,
“Joint bayesian endmember extraction and linear unmixing for hyperspectral
imagery,” IEEE Transactions on Signal Processing, vol. 57, no. 11,
pp. 4355–4368, 2009.

[20]
L. Miao and H. Qi, “Endmember extraction from highly mixed data using
minimum volume constrained nonnegative matrix factorization,” IEEE
Transactions on Geoscience and Remote Sensing, vol. 45, no. 3, pp. 765–777,

[21]
J. Li, J. M. Bioucas-Dias, A. Plaza, and L. Liu, “Robust collaborative
nonnegative matrix factorization for hyperspectral unmixing,” IEEE
Transactions on Geoscience and Remote Sensing, vol. 54, no. 10, pp.
6076–6090, 2016.

[22]
J. Li, J. M. Bioucas-Dias, and A. Plaza, “Collaborative nonnegative
matrix factorization for remotely sensed hyperspectral unmixing,” in
IEEE International Geoscience and Remote Sensing Symposium (IGARSS),
July 2012, pp. 3078–3081.

[23]
M. Berman, H. Kiiveri, R. Lagerstrom, A. Ernst, R. Dunne, and J. Huntington,
“Ice: a statistical approach to identifying endmembers in hyperspectral
images,” IEEE Transactions on Geoscience and Remote Sensing, vol. 42,
no. 10, pp. 2085–2095, 2004.

[24]
L. Zhuang, C. Lin, M. A. T. Figueiredo, and J. M. Bioucas-Dias,
“Regularization parameter selection in minimum volume hyperspectral
unmixing,” IEEE Transactions on Geoscience and Remote Sensing,
vol. 57, no. 12, pp. 9858–9877, 2019.

[25]
M. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Sparse unmixing of
hyperspectral data,” IEEE Transactions on Geoscience and Remote
Sensing, vol. 49, no. 6, pp. 2014–2039, 2011.

[26]
——, “Total variation spatial regularization for sparse hyperspectral
unmixing,” IEEE Transactions on Geoscience and Remote Sensing,
vol. 50, no. 11, pp. 4484–4502, 2012.

[27]
B. Rasti and B. Koirala, “Suncnn: Sparse unmixing using unsupervised
convolutional neural network,” IEEE Geoscience and Remote Sensing
Letters, pp. 1–5, 2021.

[28]
He, Kaiming and Zhang, Xiangyu and Ren, Shaoqing and Sun, Jian, “Identity
mappings in deep residual networks,” in European conference on
computer vision.1em plus 0.5em minus 0.4emSpringer, 2016, pp.

[29]
S. K. Roy, P. Kar, D. Hong, X. Wu, A. Plaza, and J. Chanussot, “Revisiting
deep hyperspectral feature extraction networks via gradient centralized
convolution,” IEEE Transactions on Geoscience and Remote Sensing,

[30]
B. Palsson, J. R. Sveinsson, and M. O. Ulfarsson, “Blind hyperspectral
unmixing using autoencoders: A critical comparison,” IEEE Journal of
Selected Topics in Applied Earth Observations and Remote Sensing, pp. 1–1,

[31]
Y. Su, A. Marinoni, J. Li, J. Plaza, and P. Gamba, “Stacked
nonnegative sparse autoencoders for robust hyperspectral unmixing,”
IEEE Geoscience and Remote Sensing Letters, vol. 15, no. 9, pp.
1427–1431, 2018.

[32]
Y. Su, J. Li, A. Plaza, A. Marinoni, P. Gamba, and S. Chakravortty,
“Daen: Deep autoencoder networks for hyperspectral unmixing,” IEEE
Transactions on Geoscience and Remote Sensing, vol. 57, no. 7, pp.
4309–4321, 2019.

[33]
R. A. Borsoi, T. Imbiriba, and J. C. M. Bermudez, “Deep generative
endmember modeling: An application to unsupervised spectral unmixing,”
IEEE Transactions on Computational Imaging, vol. 6, pp. 374–384,

[34]
Q. Jin, Y. Ma, F. Fan, J. Huang, X. Mei, and J. Ma, “Adversarial autoencoder
network for hyperspectral unmixing,” IEEE Transactions on Neural
Networks and Learning Systems, pp. 1–15, 2021.

[35]
M. Tang, Y. Qu, and H. Qi, “Hyperspectral nonlinear unmixing via generative
adversarial network,” in IGARSS 2020-2020 IEEE International
Geoscience and Remote Sensing Symposium.1em plus 0.5em minus
0.4emIEEE, 2020, pp. 2404–2407.

[36]
S. K. Roy, J. M. Haut, M. E. Paoletti, S. R. Dubey, and A. Plaza, “Generative
adversarial minority oversampling for spectral–spatial hyperspectral image
classification,” IEEE Transactions on Geoscience and Remote Sensing,
vol. 60, pp. 1–15, 2021.

[37]
Y. Qu and H. Qi, “udas: An untied denoising autoencoder with sparsity for
spectral unmixing,” IEEE Transactions on Geoscience and Remote
Sensing, vol. 57, no. 3, pp. 1698–1712, 2019.

[38]
X. Zhang, Y. Sun, J. Zhang, P. Wu, and L. Jiao, “Hyperspectral
unmixing via deep convolutional neural networks,” IEEE Geoscience and
Remote Sensing Letters, vol. 15, no. 11, pp. 1755–1759, 2018.

[39]
B. Palsson, M. O. Ulfarsson, and J. R. Sveinsson, “Convolutional
autoencoder for spectral-spatial hyperspectral unmixing,” IEEE
Transactions on Geoscience and Remote Sensing, pp. 1–15, 2020.

[40]
L. Gao, Z. Han, D. Hong, B. Zhang, and J. Chanussot, “Cycu-net:
Cycle-consistency unmixing network by learning cascaded autoencoders,”
IEEE Transactions on Geoscience and Remote Sensing, pp. 1–14, 2021.

[41]
F. Khajehrayeni and H. Ghassemian, “Hyperspectral unmixing using deep
convolutional autoencoders in a supervised scenario,” IEEE Journal of
Selected Topics in Applied Earth Observations and Remote Sensing, vol. 13,
pp. 567–576, 2020.

[42]
B. Rasti, B. Koirala, P. Scheunders, and P. Ghamisi, “UnDIP: Hyperspectral
unmixing using deep image prior,” IEEE Transactions on Geoscience and
Remote Sensing, pp. 1–15, 2021.

[43]
B. Rasti, B. Koirala, P. Scheunders, and J. Chanussot, “Misicnet: Minimum
simplex convolutional network for deep hyperspectral unmixing,” IEEE
Transactions on Geoscience and Remote Sensing, pp. 1–1, 2022.

[44]
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez,
Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in
Advances in neural information processing systems, 2017, pp.

[45]
A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai,
T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly et al.,
“An image is worth 16x16 words: Transformers for image recognition at
scale,” arXiv preprint arXiv:2010.11929, 2020.

[46]
J. Guo, K. Han, H. Wu, C. Xu, Y. Tang, C. Xu, and Y. Wang, “Cmt: Convolutional
neural networks meet vision transformers,” arXiv preprint
arXiv:2107.06263, 2021.

[47]
C.-F. R. Chen, Q. Fan, and R. Panda, “Crossvit: Cross-attention multi-scale
vision transformer for image classification,” in Proceedings of the
IEEE/CVF International Conference on Computer Vision, 2021, pp. 357–366.

[48]
F. Zhu, Y. Wang, B. Fan, S. Xiang, G. Meng, and C. Pan, “Spectral unmixing via
data-guided sparsity,” IEEE Transactions on Image Processing,
vol. 23, no. 12, pp. 5412–5427, 2014.

[49]
M. E. Schaepman, M. Jehle, A. Hueni, P. D'Odorico, A. Damm, J. Weyermann, F. D.
Schneider, V. Laurent, C. Popp, F. C. Seidel, K. Lenhard, P. Gege,
C. Küchler, J. Brazile, P. Kohler, L. De Vos, K. Meuleman, R. Meynart,
D. Schläpfer, M. Kneubühler, and K. I. Itten, “Advanced radiometry
measurements and earth science applications with the airborne prism
experiment (apex),” Remote Sensing of Environment, vol. 158, pp.
207–219, 2015.

[50]
L. Drumetz, T. R. Meyer, J. Chanussot, A. L. Bertozzi, and C. Jutten,
“Hyperspectral image unmixing with endmember bundles and group sparsity
inducing mixed norms,” IEEE Transactions on Image Processing,
vol. 28, no. 7, pp. 3435–3450, 2019.

[51]
J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl, “Algorithms for
hyper-parameter optimization,” Advances in neural information
processing systems, vol. 24, 2011.


# Lubrication-mediated rebounds off fluid baths

K. A. Phillips1,2  P. A. Milewski1,3  1 Department of Mathematical Sciences,
University of Bath, Bath BA2 7AY, United Kingdom 3 Department of Mathematics,
Pennsylvania State University, State College, PA 16802, USA

###### Abstract

We present herein the derivation of a lubrication-mediated (LM) quasi-
potential model for droplet rebounds off deep liquid baths, assuming the
presence of a persistent dynamic air layer which acts as a lubricating
pressure transfer. We then present numerical simulations of the LM model for
axisymmetric rebounds of solid spheres and compare quantitatively to current
results in the literature, including experimental data in the low speed impact
regime. In this regime the LM model has the advantage of being far more
computationally tractable than DNS and is also able to provide detailed
behaviour within the micro-metric thin lubrication region. The LM system has
an interesting mathematical structure, with the lubrication layer providing a
free boundary elliptic problem mediating the drop and bath free-boundary
evolutionary equations.

###### keywords:

##  1 Introduction

The role of the entrained air layer present throughout droplet rebounds has
been well studied. Since the work of Worthington (1881), where it was first
questioned, it has been well established that rebounds are only able to occur
in the presence of a lubricating layer of air preventing the two liquids to
have contact. In the case of immiscible fluids, or one of the impactors being
solid, the same assumption can be made in the pre-impact time; that the layer
of air helps aid the rebound to occur (Moore, 2021). Recent studies of drops
that bounce exist in several configurations, for example off solid hydrophobic
surfaces Kolinski et al. (2014) and off liquid films Tang et al. (2019).

There are fewer studies of drops rebounding off deep baths, particularly in
the regime where the underlying fluid is inviscid (or nearly so) and capillary
waves are long lived after the droplet impact. These models gained interest in
part because of the work of Couder on indefinitely bouncing drops on vibrating
baths Couder et al. (2005) and its fascinating quantum mechanical analogies
(Bush (2015) and references therein). The studies of (Galeano-Rios et al.,
2017, 2019) and Alventosa et al. (2023) imposed a so-called kinematic match
between the solid sphere or drop and the bath to avoid modelling the air layer
and obtained accurate results when compared to Faraday-pilot wave studies of
bouncing and walking droplets. In (Galeano-Rios et al., 2021) the kinematic
match model was validated against experiments and DNS simulations of the
Navier-Stokes equations of single rebounds of solid hydrophobic spheres for
low Weber number impacts.

Current reduced models for solid or droplet rebounds on deep baths do not
capture the dynamics of the air layer which is present throughout the drop-
bath dynamics. This paper builds off the work on the impact of two-dimensional
drops in Phillips et al. (2024), where the lubrication air layer is considered
as an additional fluid region coupling the drop deformation and the wave
equations for the bath. We extend that model to three-dimensions in generality
(Section 2), and compare its results to the kinematic match, direct numerical
simulations and experiments for solid hydrophobic sphere impacts on water
(Section 3).

##  2 Model

Consider an almost spherical droplet falling vertically downwards through the
air, towards a deep liquid bath initially at rest. Let the subscripts
$\alpha\in\\{a,b,d\\}$ denote the air, bath, and droplet respectively, and
introduce notation $\Omega_{\alpha}$ for each fluid domain,
$\mathbf{u}_{\alpha}$ for the fluid velocity, with the density, dynamic
viscosity, kinematic viscosity, and pressure of the fluids are given as
$\rho_{\alpha},\mu_{\alpha},\nu_{\alpha},p_{\alpha}$ respectively. Finally we
use $\beta\in\\{b,d\\}$ to denote the bath-air and droplet-air boundaries
$\partial\Omega_{\beta}$, defined by a level set of $F_{\beta}$ and the
surface tension coefficient $\sigma_{\beta}$ at the air-liquid interfaces. In
each fluid region we begin from the Navier-Stokes equations, and at the
interfaces the usual kinematic condition, stress continuity and velocity
continuity boundary conditions $\mathbf{u}_{a}=\mathbf{u}_{\beta}$ are taken.
Hence,

$\displaystyle\nabla\cdot\mathbf{u}_{\alpha}$ $\displaystyle=0,\qquad$
$\displaystyle\mathbf{x}$ $\displaystyle\in\Omega_{\alpha},$ (1)
$\displaystyle\rho_{\alpha}(\partial_{t}\mathbf{u}_{\alpha}+(\mathbf{u}_{\alpha}\cdot\nabla)\mathbf{u}_{\alpha})$
$\displaystyle=-\nabla
p_{\alpha}+\mu_{\alpha}\nabla^{2}\mathbf{u}_{\alpha}-\rho_{\alpha}g\mathbf{e}_{z},\qquad$
$\displaystyle\mathbf{x}$ $\displaystyle\in\Omega_{\alpha},$ (2)
$\displaystyle D_{t}F_{\beta}$ $\displaystyle=0,\qquad$
$\displaystyle\mathbf{x}$ $\displaystyle\in\partial\Omega_{\beta},$ (3)
$\displaystyle
p_{\beta}\mathbf{n}_{\beta}+\mu_{\beta}\tau_{\beta}\mathbf{n}_{\beta}$
$\displaystyle=p_{a}\mathbf{n}+\mu_{a}\tau_{a}\mathbf{n}_{\beta}-\sigma_{\beta}\kappa_{\beta},\qquad$
$\displaystyle\mathbf{x}$ $\displaystyle\in\partial\Omega_{\beta}.$ (4)

Where
$\tau_{\alpha}=(\nabla\mathbf{u}_{\alpha}+\nabla\mathbf{u}_{\alpha}^{T})$ is
the stress tensor for each fluid, $\kappa_{\beta}$ denotes curvature of the
free surface with outward unit normal vector $\mathbf{n}_{\beta}$, and $D_{t}$
is the material derivative. These equations are supplemented by imposing decay
of velocities in the far field or appropriate boundary conditions in the case
of finite domains. We proceed with several modelling approximations which will
reduce the system to a tractable state.

###  2.1 Linear quasi-potential model for the liquid bath

We assume viscous effects are small except in the lubrication layer, and
follow methodology from Dias et al. (2008), and take a Helmholtz decomposition
of the liquid as small perturbation from potential flow
$\mathbf{u}_{l}=\nabla\phi+\nabla\times\psi$. The $\psi$ terms are then
written in terms of $\phi$ through a boundary layer argument. The free surface
of the bath is given by $z=\eta_{b}(\mathbf{x},t)$ where $z=0$ is its
undisturbed free surface. The governing linear system then becomes

$\displaystyle 0$ $\displaystyle=\Delta_{H}\phi+\partial_{z}^{2}\phi,\qquad$
$\displaystyle z$ $\displaystyle\leq 0,$ (5) $\displaystyle\partial_{t}\phi$
$\displaystyle=-\frac{1}{\rho_{b}}p_{a}+g\eta_{b}+\frac{\sigma}{\rho_{b}}\Delta_{H}\eta_{b}+2\nu_{l}\Delta_{H}\phi,\qquad$
$\displaystyle z$ $\displaystyle=0,$ (6) $\displaystyle\partial_{t}\eta_{b}$
$\displaystyle=\partial_{z}\phi+2\nu_{l}\Delta_{H}\eta,\qquad$ $\displaystyle
z$ $\displaystyle=0,$ (7)

where $\Delta_{H}$ is used to denote the horizontal Laplacian, and velocity
continuity $\mathbf{u}_{b}=\mathbf{u}_{a}$ is imposed at the interface
$\eta_{b}$ between the bath and air layer.

###  2.2 Air layer lubrication approximation

The lower part of the droplet is described by
$z=\eta_{d}(\mathbf{x}_{H},t)=Z(t)-S(\mathbf{x}_{H},t)$, composed of the shape
of the droplet $S$, and the vertical height of its centre of mass $Z$. The
domain of $S$ in $\mathbf{x}_{H}$ is chosen such that $\nabla_{H}$ is well
defined. The dynamics of the air layer are described through a lubrication
approximation as in Phillips et al. (2024). Consider a lubrication region
$\Omega_{a}=\Omega_{L}\times[\eta_{b}(r,\theta,t),\eta_{d}(r,\theta,t)]$ where
its footprint is described in cylindrical coordinates as

$\Omega_{L}=\\{(r,\theta):r\in[0,r^{*}(\theta,t)],\theta\in[0,2\pi)\\},$ (8)

where $r^{*}$ is the edge of the lubrication region. Denoting
$h=\eta_{d}-\eta_{b}$, the criterion for the edge of the lubrication layer is
given by a criterion $h(r^{*},\theta,t)=\varepsilon$. Assuming that the layer
between the droplet and the bath is thin, the leading order balance of the
Navier-Stokes equations result in the lubrication balance where
$\partial_{z}p_{a}=0$ with

$\displaystyle\nabla_{H}\cdot\mathbf{u}_{a}^{H}+\partial_{z}w_{a}$
$\displaystyle=0,\qquad$ $\displaystyle\mathbf{x}$
$\displaystyle\in\Omega_{a},$ (9)
$\displaystyle\mu\partial_{z}^{2}\mathbf{u}_{a}^{H}$
$\displaystyle=\nabla_{H}p_{a},\qquad$ $\displaystyle\mathbf{x}$
$\displaystyle\in\Omega_{a},$ (10)

and velocity continuity holds at both interfaces:
$\mathbf{u}_{a}=\mathbf{u}_{d}$ at $z=\eta_{d}$ and
$\mathbf{u}_{a}=\nabla_{H}\phi$ at $z=\eta_{b}$. Here $\nabla_{H}$ is the
horizontal gradient operator, and $\mathbf{u}_{a}^{H}=(u_{r},u_{\theta})$ is
the horizontal velocity of the air layer, such that
$\mathbf{u}_{a}=(\mathbf{u}_{a}^{H},w_{a})$. As expected, the pressure
$p_{a}=P(r,\theta,t)$ is independent of $z$ throughout the air layer. Further,
we can solve the equation for $\mathbf{u}_{a}^{H}$ in terms of $P$:

$\mathbf{u}_{a}^{H}=\frac{\nabla_{H}P}{2\mu_{a}}(z-\eta_{b})(z-\eta_{d})+\frac{\nabla_{H}\phi_{b}|_{\eta_{b}}}{\eta_{b}-\eta_{d}}(z-\eta_{d})+\frac{\mathbf{u}^{H}_{d}|_{\eta_{d}}}{\eta_{d}-\eta_{b}}(z-\eta_{b}).$
(11)

To obtain the thin film equation we integrate vertically the conservation of
mass (incompressibility) equation resulting in

$\nabla_{H}\cdot\int_{\eta_{b}}^{\eta_{d}}\mathbf{u}_{a}^{H}\;\text{d}z+\left(\mathbf{u}_{a}^{H}\Big{|}_{\eta_{b}}\cdot\nabla_{H}\eta_{b}-\mathbf{u}_{a}^{H}\Big{|}_{\eta_{d}}\cdot\nabla_{H}\eta_{d}\right)+\left(w_{a}\Big{|}_{\eta_{d}}-w_{a}\Big{|}_{\eta_{b}}\right)=0.$
(12)

All terms after the integral can expressed exactly as
$\partial_{t}\eta_{d}-\partial_{t}\eta_{b}=\partial_{t}h$ from the kinematic
conditions at both interfaces. Finally, evaluating the integral in the first
term using (11) to obtain the form of the flux

$\mathbf{Q}=\int_{\eta_{b}}^{\eta_{d}}\mathbf{u}_{a}^{H}\;\text{d}z=-\frac{1}{12\mu_{a}}h^{3}\nabla_{H}P+\frac{1}{2}h\left(\nabla_{H}\phi_{b}+\mathbf{u}_{d}^{H}\right),$
(13)

results in the thin film equation

$\partial_{t}h+\nabla_{H}\cdot\mathbf{Q}=0,\qquad\mathbf{x_{H}}\in\Omega_{L}.$
(14)

As we shall see below this should be considered a nonlinear free boundary
elliptic equation for $P$ where $-h_{t}$ is considered as a forcing term, with
the boundary condition $P=0$ at the evolving boundary curve
$\partial\Omega_{L}$.

###  2.3 Droplet Model

####  2.3.1 Droplet Deformation

Linear damped droplet oscillations are well understood (Lamb (1924)) and the
results have been used in a variety of contexts e.g. (Aalilija et al., 2020).
In this section we write equivalent equations for these oscillations using the
same quasi-potential theory used above for the bath. Consider a capillary
scale droplet of liquid which when unperturbed is a sphere of radius $R_{0}$.
Describing the problem using spherical polar coordinates with azimuthal angle
$\theta$ and polar angle $\varphi$, we can take a Helmholtz decomposition for
the liquid velocity $\mathbf{u}_{d}=\nabla\phi_{d}+\nabla^{\perp}\psi$, and
using a boundary layer argument (Phillips et al. (2024)) obtain similar
governing equations to the bath equations (5)-(7)

$\displaystyle\Delta\phi_{d}$ $\displaystyle=0,$ $\displaystyle r<R_{0},$ (15)
$\displaystyle\partial_{t}\phi_{d}$
$\displaystyle=-\frac{1}{\rho_{d}}p_{a}+\frac{\sigma}{\rho_{d}}\kappa_{d}+2\nu_{d}\partial_{r}^{2}\phi_{d},$
$\displaystyle r=R_{0},$ (16) $\displaystyle\partial_{t}\zeta_{d}$
$\displaystyle=\partial_{r}\phi_{d}+2\nu_{d}\left(\frac{1}{r^{2}}\Delta_{s}\zeta_{d}-\frac{1}{r^{3}}\int\Delta_{s}\phi_{d}\;\text{d}t\right),$
$\displaystyle r=R_{0},$ (17)

where surface of the droplet is given by $r=R_{0}+\zeta_{d}(\theta,\varphi)$,
and $\kappa_{d}$ is the mean curvature

$\kappa_{d}=-\left(\frac{2}{R_{0}^{2}}\zeta_{d}+\frac{1}{R_{0}^{2}}\Delta_{s}\zeta_{d}\right),\qquad\Delta=\frac{1}{r^{2}}\partial_{r}(r^{2}\partial_{r})+\frac{1}{r^{2}}\Delta_{s},$

with $\Delta_{s}$ the Laplace-Beltrami operator (i.e. the surface Laplacian).
The equations are valid in the small oscillation limit and have been truncated
to contain terms that are linear in the Ohnesorge number
$Oh=\frac{\mu}{\sqrt{\rho R_{0}\sigma}}$. One may proceed with an
eigenfunction decomposition of $\zeta_{d}$ and $\phi_{d}$ as

$\phi_{d}=\sum_{l=2}^{\infty}\sum_{m=-l}^{l}a_{lm}(t)r^{l}Y_{l}^{m},\qquad\zeta_{d}=\sum_{l=2}^{\infty}\sum_{m=-l}^{l}c_{lm}(t)Y_{l}^{m},$

where the $Y_{l}^{m}(\theta,\varphi)=e^{im\varphi}P_{l}^{m}(cos\theta)$ are
spherical harmonics satisfying

$-\Delta_{s}Y_{l}^{m}=(l)(l+1)Y_{l}^{m}.$

The $P_{l}^{m}$ are the associated Legendre polynomials of degree $l$ and
order $m$. The terms $l=0,1$ have been omitted as $l=0$ corresponds to a
dilation of the sphere and $l=1$ corresponds to translational and rotational
symmetries which need to be considered separately as they will be affected by
the external forces due to the air layer. Substitution into the boundary
conditions (16)-(17) leads to the differential equation (where terms of order
$Oh^{2}$ have been omitted)

$\ddot{c}_{l,m}+2\lambda_{l}\dot{c}_{l,m}+\omega_{l}^{2}c_{l,m}=-\frac{l\hat{p}_{l,m}}{\rho
R_{0}}\;,$ (18)

where $\hat{p}_{l,m}$ is the $Y_{l,m}$ coefficient of of the lubrication air
pressure $p_{a}$. The damping coefficients and modal frequencies are given by

$\lambda_{l}=(2l+1)(l-1)\frac{\mu}{\rho
R_{0}^{2}},\qquad\omega_{l}^{2}=l(l-1)(l+2)\frac{\sigma}{\rho R_{0}^{3}},$
(19)

as seen in the literature (Lamb, 1924; Aalilija et al., 2020; Alventosa et
al., 2023).

####  2.3.2 Equations for the centre of mass

The location of the droplet is determined through updating the position of the
centre of mass, given by $\mathbf{X}(t)$. We shall disregard drag effects from
the air, instead only considering contributions from the lubrication layer,
then the equation of motion is

$m\ddot{\mathbf{X}}=F=-\int_{\Omega_{L}}p_{a}\;\mathbf{n_{d}}\;\text{d}a-mg\;\hat{\mathbf{z}},$
(20)

where $\mathbf{n_{d}}$ is the outward normal of the droplet $m=4\rho\pi
R_{0}^{3}/3$ is its mass, $g$ is the gravitational constant and
$\hat{\mathbf{z}}$ is the unit vector in the vertical direction. The boundary
of the droplet is now given by
$\mathbf{X}+(R_{0}+\eta_{d}(\theta,\varphi))\mathbf{\hat{r}}$, where
$\mathbf{\hat{r}}$ is a radial unit vector relative to the position
$\mathbf{X}$.

We note that we have disregarded global rotational motion of the droplet
(forced by both asymmetric pressure and shear stresses in the lubrication
layer) which would significantly add to the modelling complexity, and which
will not be relevant for the vertical (axisymmetric) impacts we shall consider
next. We also expect these effects to be small for quasi-normal impacts (e.g.
Galeano-Rios et al. (2019)).

###  2.4 The lubrication-mediated droplet-bath system

We can now propose a closed system of equations for the evolution of the thin
film, liquid bath and droplet. Within the bath and the droplet one may reduce
the system further using the framework of Dirichelet to Neumann (DtN) maps.
Consider the “trace” of the potentials on the boundary of the domains in which
Laplace’s equation are solved

$\Phi_{b}(\mathbf{x_{H}},t)=\phi_{b}|_{z=0},\qquad\Phi_{d}(\theta,\varphi,t)=\phi_{d}|_{r=R_{0}}.$

The DtN map, denoted $D$, is the linear map that provides the normal
derivative of the potential given its value on the boundary. In particular we
denote

$(\partial_{z}\phi_{b})|_{z=0}=D_{b}\Phi_{b},\qquad(\partial_{r}\phi_{d})|_{r=R_{0}}=D_{d}\Phi_{d}.$

We shall see that these maps are easily expressed in terms of eigenfunction
expansions. Making use of DtN maps we now have the evolution system

$\displaystyle\partial_{t}\Phi_{b}$
$\displaystyle=-\frac{1}{\rho_{b}}p_{a}+\frac{\sigma}{\rho_{b}}\Delta_{H}\eta_{b}+2\nu_{l}\Delta_{H}\Phi_{b},$
(21) $\displaystyle\partial_{t}\eta_{b}$
$\displaystyle=D_{b}\Phi_{b}+2\nu_{l}\Delta_{H}\eta_{b}\doteq F_{b},$ (22)
$\displaystyle\partial_{t}\Phi_{d}$
$\displaystyle=-\frac{1}{\rho_{d}}p_{a}+\frac{\sigma}{\rho_{d}}\kappa_{d}+2\nu_{d}\partial_{r}^{2}\Phi_{d},$
(23) $\displaystyle\partial_{t}\eta_{d}$
$\displaystyle=D_{d}\Phi_{d}+2\nu_{d}\left(\frac{1}{r^{2}}\Delta_{s}\eta_{d}-\frac{1}{r^{3}}\Delta_{s}D_{d}^{-1}\Phi_{d}\right)\doteq
F_{d},$ (24) $\displaystyle m\ddot{\mathbf{X}}$
$\displaystyle=F=-\int_{\Omega_{L}}p_{a}\;\mathbf{n_{d}}\;\text{d}a-mg\;\hat{\mathbf{z}},$
(25) $\displaystyle\nabla_{H}\cdot\mathbf{Q}$
$\displaystyle=W+F_{d}-F_{b},\qquad\mathbf{x_{H}}\in\Omega_{L}.$ (26)

In the surface potential equations, the pressure $p_{a}$ is nonzero only in
$\Omega_{L}$ and $W=\frac{d\mathbf{X}}{dt}$ is the vertical velocity of the
droplet’s centre of mass. Equation (26), where $Q$ is given by (13) is an
inhomogeneous elliptic equation for the pressure with the boundary condition
that $p_{a}=0$ on the free boundary $\partial\Omega_{L}$.

##  3 Axisymmetric solid impacts and numerical results

We now consider the simplest case in which this framework can be used and
compared to prior simulations and experiments: the rebound of solid
(hydrophobic) spheres Galeano-Rios et al. (2021). This simplifies considerably
the system, reducing the dimension and eliminating droplet oscillations.
Considering the axisymmetric domain for the bath to of finite radius $r=L$,
and imposing Neumann boundary conditions there for both $\eta_{b}$ and
$\Phi_{b}$ we expand these in terms of Bessel functions of the first kind,

$\eta_{b}(r,t)=\sum_{j=1}^{\infty}a_{j}(t)J_{0}(k_{j}r),\qquad\Phi_{b}(r,t)=\sum_{j=1}^{\infty}b_{j}(t)J_{0}(k_{j}r),$
(27)

where the $k_{j}$ satisfy $J_{0}^{\prime}(k_{j}L)=0$. We note that in this
basis $\phi_{b}=\sum_{j=1}^{\infty}b_{j}J_{0}(k_{j}r)e^{k_{j}z}$ to satisfy
Laplace’s equation and hence the DtN map is expressed simply as

$D_{b}\Phi=\partial_{z}\phi_{b}|_{z=0}=\sum_{j=1}^{\infty}k_{j}b_{j}J_{0}(k_{j}r).$

Hence the DtN operator in this basis corresponds to multiplication by $k_{j}$.
Similarly, the horizontal Laplacian $\Delta_{H}$ results in a Bessel
multiplier $-k_{j}^{2}$. The system can now be written for the Bessel
coefficients as

$\displaystyle\dot{a}_{j}$ $\displaystyle=k_{j}b_{j}-2\nu_{l}k_{j}^{2}a_{j},$
(28) $\displaystyle\dot{b}_{j}$
$\displaystyle=-\frac{1}{\rho_{l}}\hat{p}_{j}-\frac{\sigma}{\rho_{l}}k^{2}a_{j}-2\nu_{l}k^{2}b_{j},$
(29) $\displaystyle m\ddot{Z}$
$\displaystyle={2\pi}\int_{0}^{r^{*}}p_{a}\;r\text{d}r-mg,$ (30)

with the elliptic equation for pressure

$-\frac{1}{r}\partial_{r}\left(\frac{r}{12\mu_{a}}h^{3}\partial_{r}p_{a}+\frac{r}{2}h\partial_{r}\Phi_{b}\right)=\partial_{t}h=W-(D_{b}\Phi_{b}+2\nu_{l}\Delta_{H}\eta_{b}),$
(31)

where $W=\dot{Z}$, $h=Z-\sqrt{R_{0}^{2}-r^{2}}-\eta_{b}$ for $r\leq R_{0}$,
$r^{*}$ is defined by $h(r^{*},t)=\epsilon$ and where $\hat{p}_{j}$ are the
Bessel coefficients of the Pressure provided by

$\hat{p}_{j}=\frac{2}{(LJ_{0}(k_{j}L))^{2}}\int_{0}^{L}p_{a}(r,t)J_{0}(k_{j}r)r\text{d}r\;.$
(32)

Similar expressions apply for obtaining $a_{j}$ from $\eta_{b}$ and $b_{j}$
from $\Phi_{b}$. The pressure can be integrated from (31) using
$p_{a}(r^{*},t)=0$ and the symmetry condition $\partial_{r}p_{a}(r^{*},t)=0$:

$p_{a}=-\int_{r^{*}}^{r}\left[\frac{12\mu_{a}}{r^{\prime}h(r^{\prime},t)^{3}}\left(\int_{0}^{r^{\prime}}r^{\prime\prime}\partial_{t}h(r^{\prime\prime},t)\text{d}r^{\prime\prime}\right)+\frac{6\mu_{a}}{h(r^{\prime},t)^{2}}\partial_{r}\Phi_{b}(r^{\prime},t)\right]\text{d}r^{\prime}\;.$

###  3.1 Numerical implementation

We shall compute a solutions approximating the system (28)-(3) by truncating
the Bessel expansions at a large $N$ and integrating the resulting
differential equations (28)-(30) with a fourth order Runge-Kutta scheme. The
equation for the pressure (3) is calculated using a trapezium rule. Since
Bessel expansions are ill-conditioned on regular grids, we use an oversampled
nonuniform grid in $r$ with $M$ points distributed on Chebyshev collocation
points on $[0,L]$ using $\theta_{j}=(j-1)\pi/(M-1)$ with $j=1\ldots M$ and
$r_{j}=\frac{L(1+\cos(\theta_{j}))}{2}$. Hence we compute two matrices: an
$M\times N$ matrix which evaluates an N-term Bessel series at $r_{j}$ and an
$N\times M$ matrix corresponding to calculating the projection formula (32)
for the Bessel coefficients of a function given at $r_{j}$.

###  3.2 Model Results

Parameter | Drop | Value | Air | Value | Bath | Value | Units
---|---|---|---|---|---|---|---
Density | $\rho_{d}$ | 250-1200 | $\rho_{a}$ | 1.225 | $\rho_{b}$ | 1000 | kgm-3
Viscosity | $\mu_{d}$ | 9.78$\times 10^{-4}$ | $\mu_{a}$ | 1.825$\times 10^{-4}$ | $\mu_{b}$ | 9.78$\times 10^{-4}$ | kgm-1s-1
Radius | $R_{0}$ | 2.5-8.3$\times 10^{-4}$ | - | - | - | - | m
Initial Velocity | $W_{0}$ | 0.1-0.45 | - | - | - | - | ms-1
Table 1: Parameters used in simulations of a water bath and various solid
spheres.

We first present results of the rebound of two representative spheres of
different radii and densities rebounding off a deep water bath against a sweep
of initial downward velocities. The choice of parameters is displayed in table
1, and were chosen to correspond to the data for low velocity impacts in
Galeano-Rios et al. (2021). This will permit the results of the lubrication-
mediated model (LM) derived within this paper to be compared to the kinematic
match (KM) model Galeano-Rios et al. (2017) and direct numerical simulation of
the Navier-Stokes equations (DNS). In that study experiments were also
performed using hydrophobically coated solid spheres.

For such studies, the coefficient of restitution is usually taken as the
negative ratio of the speeds between the time of impact $t_{\text{imp}}$,
taken to be when the south pole of the sphere first crosses $z=0$
(correspondingly the centre of mass crosses $z=R_{0}$), and the time of
liftoff $t_{\text{lift}}$, which is when the sphere leaves passes the same
height on an upwards trajectory. The time spent below $z=0$ is the contact
time $t_{c}=t_{\text{lift}}-t_{\text{imp}}$. In Galeano-Rios et al. (2017) for
low impact velocities, the square coefficient of restitution $\alpha^{2}$,
which is the ratio of mechanical energies before and after impact (using the
$z=0$ crossing as reference height) is used, as it can capture very small
rebounds where $\alpha^{2}<0$ where the sphere detaches but it’s centre of
mass does not reach the $z=R_{0}$ after rebound. Thus

$\alpha^{2}=\frac{E_{\text{out}}}{E_{\text{in}}}=\frac{W_{\text{imp}}^{2}}{W_{\text{detach}}^{2}-2g(Z-R_{0})},$
(33)

where $W_{\text{imp}}$ and $W_{\text{detach}}$ are the vertical velocities at
the moment of impact $t_{\text{imp}}$ and separation $t_{\text{detach}}$, when
$\text{min}(h)>\epsilon$ is first achieved after impact. In Figure 1,
$\alpha^{2}$, the maximum penetration depth $\delta=-(Z_{\text{min}}-R_{0})$,
and the pressing time $t_{p}=t_{\text{detach}}-t_{\text{imp}}$, are displayed
for the LM and KM models and DNS. In Galeano-Rios et al. (2017) $t_{p}$ is
favoured over $t_{c}$ to capture the rebounds which don’t return to $z=0$.

Figure 1: Comparison of $\alpha^{2}$, $\delta$, and $t_{p}$ for two parameter
regimes, with LM ($\circ$), compared to DNS ($\Diamond$) and KM ($*$) from
Galeano-Rios et al. (2021).

Figure 1 demonstrates that the agreement is good between the lubrication
mediated model against both the results from the KM and DNS, with general
trends being upheld. We note here that the latter two methods also involve
certain approximations. The KM model assumes no lubrication layer (the method
matches the position and velocities of the free surface to those of the
sphere) and uses a tangential contact between the sphere and the bath
throughout the impact. In addition, for the KM model we are comparing with,
the evolution equations for the bath are linear and with the same
quasipotential approximation that we use here. The DNS results in Galeano-Rios
et al. (2021) were performed on the open source Gerris where the solid drop is
modelled by using artificially high viscosity and surface tension
coefficients. Gerris uses a volume of fluid formulation and interfaces are
approximated using continuous density changes (see Popinet (2009)).

System | $\alpha$ | $\delta$ | $t_{c}$
---|---|---|---
Lubrication-mediated | 0.115-0.133 | 8.2-8.3$\times 10^{-4}$ m. | 0.0164-0.0168
Kinematic match | 0.15 | 9.5$\times 10^{-4}$m. | 0.0175
Direct numerical simulation | 0.24 | 9$\times 10^{-4}$m. | 0.0155
Experimental | 0.12-0.18 | 9-10$\times 10^{-4}$m. | 0.0155-0.0175
Table 2: Results for a falling sphere of radius $R_{0}=0.83\times 10^{-4}$ m,
density $\rho=1200$kgm-3, and initial velocity $W_{0}=-0.35$ms-1, onto a deep
water bath. The range in the lubrication-mediated data arise from varying bath
radius $L$.

We conclude by providing numerical results for a configuration which can be
compared to the KM model, DNS, and also experimental values, as shown in table
2. The experimental configuration used a rectangular container while our
simulation is cylindrical. In table 2 we varied the radius of the bath over
comparable values to their length and width dimensions resulting in a range of
values for our simulations.

The results for the lubrication mediated model restitution coefficient
$\alpha$ is at the lower end of the experimental range, the depth of impact is
somewhat below the experimental range and the impact time is well within the
experimental range. In addition to modelling assumptions, one potential source
of error is that the physics of the hydrophobic coating used in experiments is
not considered in any of the models.

(a)

(b)

Figure 2: (a) Thickness of the air layer, $h(r^{*},t)$. (b) Free surface of
the bath, $\eta^{*}_{b}(r^{*},t)$ shows the sinking depth of the sphere
relative to the resultant waves on the free surface.

Figures 2(a) and 3(b), show the rapid inertially-driven expansion of the
lubrication layer, lasting $\mathcal{O}(1\;\text{ms})$ followed by a gradual
retraction of the layer, lasting $\mathcal{O}(10\;\text{ms})$. Finally a short
“suction” event lasting $\mathcal{O}(1\;\text{ms})$ precedes separation. A
highlight of the lubrication mediated model is the ability to capture and
quantify the behaviour within the air layer during the rebound. In particular,
we can study the evolution of the pressure profile and corresponding height
layer within the system, quantities which are challenging to extract even in
the DNS setting. In what follows, we have nondimensionalised all lengths with
the radius of the sphere (e.g. $r^{*}=r/R_{0}$), the pressure with the surface
tension pressure jump across a spherical interface $P^{*}=P/(2\sigma
R_{0}^{-1})$, and the Force with $F^{*}=F/(2\sigma R_{0})$. Throughout the
impact, the height of the air layer retains thickness $\mathcal{O}(1)\;\mu$m,
as expected from the literature (Tang et al., 2019). However there is some
variation as shown in Figure 2(a), with a thickness between 0.5-2 $\mu$m in
the bulk of the layer and a narrowing on its retraction boundary to
0.15-0.5$\mu$m. In effect, for the bulk of the impact the air layer forms a
static cushion constricted on its circular boundary. This narrowing effect is
enhanced when the lubrication region boundary “turns around” and immediately
prior to detachment, where the lubrication layer has shrunk to the centre of
the sphere. Turning to the profile of the pressure within the air layer, shown
in figure 3(b), we find a range of 0.9-1.5 $2\sigma R_{0}^{-1}$ over the vast
majority of the region. The maximum pressure is at impact with 7.6 $2\sigma
R_{0}^{-1}$ and the minimal pressure being $-$3.8 $2\sigma R_{0}^{-1}$ at
liftoff. A pressure spike occurs where the layer thins at its edge, aligning
with results found in Hendrix et al. (2016). Figure 3(a) captures the pressure
profiles during expansion and retraction. Small oscillations observed in the
pressure profile at two intermediate times are numerical artifacts, as they
are lessened by increased mesh resolution. They have minimal effect on the
overall dynamics of the system.

Figure 4 shows a temporal slice of pressure, under the south pole of the
sphere $P(0,t)$, and the total pressure force acting on the sphere. The
pressure has clearly fine scale behaviour at the initial impact and detachment
but these features do not substantially affect the overall force on the
system. Most of the work on the drop is due to the quasi-uniform pressure at
intermediate times. Figure 2(b) shows the full wavefield during impact.

All computations presented were performed in matlab with ranges
$N=2^{11}-2^{12}$, $M=2^{8}-2^{9}$. For lower impact velocities this was an
over-resolution, whereas the finer grids were needed for stronger impacts. The
domain length varied $L=12R_{0}-16R_{0}$ the latter ensuring the impact is
unaffected by waves reflecting from the boundary, even for longer impacts. The
longest simulations took $\mathcal{O}(10\;\text{hours})$ on a standard
personal computer, although the use of an adaptive timestep would probably
reduce this considerably. The only free parameter in the system is $\epsilon$,
the cutoff for lubrication and this was set to 40$\mu$m.

(a)

(b)

Figure 3: Lubrication layer pressure (a) Profiles $P^{*}(r^{*},t)$ at selected
times. Dark and light curves correspond to expanding and contracting layer
respectively. (b) Colormap showing near uniform pressure in the interior and
its rapid expansion and slow contraction. Figure 4: Comparison between the
pressure force felt by the droplet, the pressure along $r=0$ during impact,
and the pressed area. Note initial (inertial) and final (suction) spikes.
During most of the impact the force is proportional to the lubrication layer
area.

##  4 Conclusions

The lubrication-mediated model herein efficiently resolves certain features of
droplet rebounds heretofore unattainable computationally in reduced models.
Future work include computational studies including drop deformation and non-
axisymmetric impact.

This work follows several (e.g. Galeano-Rios et al. (2019, 2017)) which were
motivated by Faraday bouncing droplet dynamics (Bush, 2015). We expect as
future development to use this model with Faraday pilot-wave systems, as an
accurate reproduction of the vertical motion of the droplet is thought to be
crucial to the observed statistical features resembling quantum mechanics
(Durey et al., 2020). The present formulation of a lubrication layer mediating
the impact of two bodies of fluid also has potential interest in the
mathematics community which has recently addressed rigorously questions of
whether fluid surfaces can intersect (Fefferman et al., 2016).

## References

  * Aalilija et al. (2020) Aalilija, A, Gandin, Ch-A & Hachem, Elie 2020 On the analytical and numerical simulation of an oscillating drop in zero-gravity. Computers & Fluids 197, 104362.
  * Alventosa et al. (2023) Alventosa, Luke FL, Cimpeanu, Radu & Harris, Daniel M 2023 Inertio-capillary rebound of a droplet impacting a fluid bath. Journal of Fluid Mechanics 958, A24.
  * Bush (2015) Bush, John WM 2015 Pilot-wave hydrodynamics. Ann. Rev. of Fluid Mechanics 47, 269–292.
  * Couder et al. (2005) Couder, Yves, Protiere, Suzie, Fort, Emmanuel & Boudaoud, Arezki 2005 Walking and orbiting droplets. Nature 437 (7056), 208–208.
  * Dias et al. (2008) Dias, Frederic, Dyachenko, Alexander I & Zakharov, Vladimir E 2008 Theory of weakly damped free-surface flows: a new formulation based on potential flow solutions. Physics Letters A 372 (8), 1297–1302.
  * Durey et al. (2020) Durey, Matthew, Milewski, Paul A & Wang, Zhan 2020 Faraday pilot-wave dynamics in a circular corral. Journal of Fluid Mechanics 891, A3.
  * Fefferman et al. (2016) Fefferman, Charles, Ionescu, Alexandru D & Lie, Victor 2016 On the absence of splash singularities in the case of two-fluid interfaces .
  * Galeano-Rios et al. (2021) Galeano-Rios, Carlos A, Cimpeanu, Radu, Bauman, Isabelle A, MacEwen, Annika, Milewski, Paul A & Harris, Daniel M 2021 Capillary-scale solid rebounds: experiments, modelling and simulations. Journal of Fluid Mechanics 912, A17.
  * Galeano-Rios et al. (2017) Galeano-Rios, Carlos A, Milewski, Paul A & Vanden-Broeck, J-M 2017 Non-wetting impact of a sphere onto a bath and its application to bouncing droplets. Journal of Fluid Mechanics 826, 97–127.
  * Galeano-Rios et al. (2019) Galeano-Rios, Carlos A, Milewski, Paul A & Vanden-Broeck, J-M 2019 Quasi-normal free-surface impacts, capillary rebounds and application to faraday walkers. Journal of Fluid Mechanics 873, 856–888.
  * Hendrix et al. (2016) Hendrix, Maurice HW, Bouwhuis, Wilco, van der Meer, Devaraj, Lohse, Detlef & Snoeijer, Jacco H 2016 Universal mechanism for air entrainment during liquid impact. Journal of fluid mechanics 789, 708–725.
  * Kolinski et al. (2014) Kolinski, John M, Mahadevan, L & Rubinstein, Shmuel M 2014 Drops can bounce from perfectly hydrophilic surfaces. Europhysics Letters 108 (2), 24001.
  * Lamb (1924) Lamb, Horace 1924 Hydrodynamics. University Press.
  * Moore (2021) Moore, Matthew R 2021 Introducing pre-impact air-cushioning effects into the wagner model of impact theory. Journal of Engineering Mathematics 129 (1), 6.
  * Phillips et al. (2024) Phillips, Katie A., Cimpeanu, Radu & Milewski, Paul A. 2024 Modelling droplet rebounds off fluid baths. ArXiv .
  * Popinet (2009) Popinet, Stéphane 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228 (16), 5838–5866.
  * Tang et al. (2019) Tang, Xiaoyu, Saha, Abhishek, Law, Chung K & Sun, Chao 2019 Bouncing drop on liquid film: Dynamics of interfacial gas layer. Physics of Fluids 31 (1).
  * Worthington (1881) Worthington, AM 1881 On impact with a liquid surface. Proceedings of the Royal Society of London 33, 347–349.