Transformer for Multitemporal Hyperspectral Image Unmixing

For multitemporal hyperspectral images, it is represented as $\mathrm{\mathbf{Y}\in\mathbb{R}^{T\times L\times N}}$ , Where $\mathrm{T}$ represents the number of time phases, $\mathrm{L}$ represents the number of bands of the hyperspectral image, and $\mathrm{N}$ represents the total number of the hyperspectral image pixles, respectively. For linear unmixing of multitemporal hyperspectral images, the formula for time $t$ is as follows.

Where $\mathrm{\mathbf{M_{t}}\in{\mathbb{R}^{L\times P}}}$ represents the matrix of endmembers at time $t$, while $\mathrm{\mathbf{A_{t}}\in{\mathbb{R}^{P\times N}}}$ denotes the abundance matrix corresponding to the same timeframe, and $\mathrm{\mathbf{\varepsilon_{t}}\in{\mathbb{R}^{L\times N}}}$ represents the additional random noise added to simulate the real scenario, where $\mathrm{P}$ denotes the number of endmembers. (1) can also be refined to the pixel-level form, which satisfies the linear unmixing model for every pixel in every time phase. For the hyperspectral image at each time phase, the abundance nonnegativity constraint (ANC) and the abundance sum-to-one constraint (ASC) are satisfied.

In the multitemporal hyperspectral image unmixing task, the hyperspectral images in different time phases have spatial-temporal correlation, and the contribution of hyperspectral bands in each time phase to unmixing is not the same, so it is necessary to jointly consider the temporal-spatial-spectral information to design the unmixing model.

In this section, we introduce the overall structure of the proposed MUFormer. The framework of the MUFormer model is shown in Fig. 1.

First, we input a sequence of hyperspectral images $\mathbf{Y_{in}}$ to the encoder, where $\mathrm{\mathbf{Y_{in}}\in{\mathbb{R}^{T\times L\times H\times W}}}$, $\mathrm{H}$ and $\mathrm{W}$ represent the height and width of the images, the sequence represents the hyperspectral images acquired at the same location at different time instants. The encoder is mainly composed of convolutional layer, BatchNorm layer, and LeakyReLU layer. In order to prevent the model from overfitting, we add Dropout layer, which makes the neural network more robust and generalization ability [45]. Compared with transformer, CNN structure has more advantages in extracting local image features, and multitemporal feature maps containing richer semantic information can be obtained through multiple convolution operations [53, 54, 55]. At this point, the model outputs the multitemporal feature map $\mathbf{Y_{conv}}$, where $\mathrm{\mathbf{Y_{conv}}\in{\mathbb{R}^{T\times C\times H\times W}}}$, $\mathrm{C}$ represents the number of channels retained after convolution processing, and it is a hyperparameter that can be set. Compared with the hundreds of channels of the original hyperspectral image, the channels retained here contain more important band information.

Then, we send the obtained multitemporal feature map into the GAM. With the advantages of transformer in sequence information modeling, we divide the feature map into patch blocks of the same size to calculate attention of each dimension and fuse it. Then, we use the proposed CEM module to refine the changes of adjacent phases, and send the processed feature maps into softmax to obtain the abundance maps of each phase.

In the decoder part, we establish a linear decoder for each phase, incorporating a linear layer. The weight of this linear layer is initialized using the VCA algorithm. Upon completion of the training iteration, the weight of the decoder represents the endmember matrix of the phase. Leveraging the linear decoder, we achieve the realization of the multitemporal linear unmixing model.

In order to calculate temporal-spatial-spectral attention from a global perspective and to deeply integrate information from multiple dimensions, we design the GAM. This module mainly contains temporal attention, spatial attention, spectral attention and forward parts. Feature maps processed by convolutional layers remove redundant bands and contain more critical information. Following the design of ViT [42], the feature map sequence is divided into non-overlapping patches $\mathbf{Y_{patch}}$, where $\mathrm{\mathbf{Y_{patch}}\in{\mathbb{R}^{p\times p\times C}}}$, $\mathrm{p}$ represents the side length of the patch and $\mathrm{C}$ represents the number of channels processed by the encoder. We concatenate the $\mathbf{Y_{cls}}$ with the original patch matrix, $\mathbf{Y_{cls}}$ aggregates global features to avoid bias toward specific tokens in the sequence, and then add the position embedding vector $\mathbf{Y_{pos}}$, $\mathrm{\mathbf{Y_{pos}}\in{\mathbb{R}^{D}}}$ represents the temporal-spatial position information of the patch that can be learned. The new matrix $\mathbf{Y^{{}^{\prime}}}$ representation is given in 3.

After obtaining $\mathbf{Y^{{}^{\prime}}}$, inspired by the work of [56], we treat the hyperspectral image sequence as frames with different time intervals, and then send them to the attention module for processing. Firstly, it is sent to the temporal attention module. In the temporal attention module, $\mathbf{Y}$ will be divided into $\mathbf{Q,K,V}$ three matrices after the linear layer, which represent query, key and value respectively, where $\mathrm{\mathbf{Q}\in{\mathbb{R}^{h_{dim}\times h}}}$, $\mathrm{h_{dim}}$ represents the dimension of each head, $h$ represents the number of heads,the dimensions of $\mathbf{K}$ and $\mathbf{V}$ are the same as $\mathbf{Q}$. After obtaining the above three matrices, the self-attention weight size is calculated through the softmax function, and the calculation formula is as shown in 4.

Where $\mathrm{l}$ represents the $\mathrm{lth}$ patch block and $\mathrm{t}$ represents the hyperspectral image at time $\mathrm{t}$. By computing attention in the time dimension, we can assign reasonable weights to patches at the same location but at different moments, this enables the model to identify regions that have undergone changes at different times. By allocating higher attention to patches that exhibit significant changes throughout the time series, the model can concentrate more effectively on these key areas of change. Similarly, for patches at different positions at the same time, we also use the same way to calculate the attention weights of spatial dimensions, as shown in 5.

In order to solve the problem of gradient disappearance and gradient explosion in deep networks, we adopt residual connection [57]. The feature maps processed by temporal self-attention and spatial self-attention are fed into the forward module, which contains an MLP. This module can be set to multiple continuous modules, in order to prevent overfitting of the model here the model depth is set to 2. Different channels in hyperspectral images may correspond to different ground object components or features, and some channels may be affected by noise, atmospheric interference and other factors, resulting in low information quality. At the same time, the contribution of different channels to the unmixing task is not consistent [58, 59]. Combined with the above reasons, we send it to the spectral self-attention calculation module. In the spectral self-attention module, we respectively use Average pooling and Max pooling to process the channels, and then generate the channel attention vector through the full connection operation, which can adaptively adjust the weight of each channel to highlight important channel information. The ability to focus on channels with higher unmixing contributions is enhanced to improve the unmixing accuracy.The spectral self-attention mechanism is calculated as shown in 6. $\mathbf{Y_{map}}$ represents the feature map outputted by FF, and $\sigma$ represents the sigmoid function. The schematic diagram of GAM is shown in Fig. 2.

Through the GAM, we dynamically assign weights to the significance of the time, space, and spectral dimensions, enabling the more precise capture of complex features within image data. This holistic attention mechanism aids models in more effectively comprehending the temporal evolution of objects or phenomena, alongside their unique spatial distributions and spectral characteristics. Concurrently, the robustness of the model is enhanced as it emphasizes crucial features while diminishing the impact of irrelevant or disruptive information. This enhancement makes our model highly skilled in dealing with the complex and fluctuating environmental conditions prevalent in hyperspectral images.

Through the above GAM, we can obtain the temporal-spatial-spectral feature information from a global perspective. However, the acquisition time interval of hyperspectral images in different time phases is not certain. Therefore, how to adaptively deal with the feature changes of hyperspectral images between adjacent time phases is also an important issue. In order to solve this problem, we propose the CEM from the feature-level perspective, and its structure diagram is shown in the bottom left corner of Fig. 1. The CEM input is the feature map of two adjacent time phases after GAM processing, which size is $\mathrm{{\mathbb{R}}^{C\times H\times W}}$. At this time, the feature map has gone through the weight allocation of global attention. In order to accurately capture the difference between adjacent temporal hyperspectral images from the feature level better, we propose to perform multi-scale convolution on the feature map at time $\mathrm{T_{k+1}}$, where the convolution kernel size is $3\times 3$ and $7\times 7$, respectively. Then the feature maps after convolution are subtracted from the feature graphs at $\mathrm{T_{k}}$ time. The subtracted feature maps are fed into the convolution layer and processed by the $Sigmoid$ function, where the kernel size is set to $1\times 1$. The weight vector $\mathbf{A}$ is obtained by multiplying the obtained weights $\mathbf{A1}$ and $\mathbf{A2}$, and the final weight vector is obtained by multiplying $\mathbf{A}$ by $\alpha$. We add the obtained weight vector to $\mathbf{Y_{cls}}$, $\mathbf{Y_{cls}}$ tags can provide a high-level representation of the entire sequence to facilitate learning and inference by the model. It allows the model to extract important global information from the input sequence rather than just focusing on local segments [60]. We modify $\mathbf{Y_{cls}}$ by combining the output of CEM, so that the model not only allocates global information weights through the attention mechanism, but also uses convolutional neural network to reasonably deal with the change regions of adjacent time phases in the local information extraction. By this way of “coarse tuning” (GAM) plus “fine tuning” (CEM), we can reasonably simulate the changes of multitemporal features, so that our model can adaptively unmix multitemporal hyperspectral images.

Where $\mathrm{f_{3\times 3}}$ and $\mathrm{f_{7\times 7}}$ represent convolution kernel sizes of 3 and 7, respectively, followed by BN layer and LeakyReLU layer. $\mathrm{f_{1\times 1}}$ represents a convolutional layer with a kernel size of 1 and $\sigma$ represents the $Sigmoid$ function, where $\alpha$ is a hyperparameter that can be set empirically. We perform multi-scale convolution on the feature map, which can capture features at different scales. The smaller $3\times 3$ kernels help to capture local details and edge features, while the larger $7\times 7$ kernels are more suitable to capture larger contextual information. By using convolution kernels at both scales, a more comprehensive multi-scale feature representation can be obtained. In the process of feature fusion, the information loss that may exist in a single scale can be reduced and the expression ability of the model can be improved. Moreover, the large-scale convolution kernel may be more effective for removing some noise and unimportant detail information in the image, while the small-scale convolution kernel helps to retain important local features, which also suppresses noise to a certain extent. The ablation study on multi-scale CEM can be seen in Section IV.

In the design of loss function, we choose the weighted sum of multiple loss functions as the total loss function to better train our proposed model. Reconstruction loss and Spectral Angle Distance (SAD) loss are commonly used loss functions in unmixing tasks. At the same time, it was demonstrated in [10] that the data simplex loss plays an important role in endmember estimation, the endmembers can be further constrained by adding the data simplex loss. The total training loss function is shown in 13, including $\mathrm{L_{RE}}$, $\mathrm{L_{SAD}}$ and $\mathrm{L_{E}}$, where $\beta$, $\gamma$ and $\lambda$ are artificially set hyperparameters to balance each loss function. We adopt Adam algorithm for optimization and set scheduler to dynamically adjust the learning rate.

In order to accurately test the effect of our proposed model, we expand the test on three datasets, including one real dataset and two synthetic datasets, which will be described in the following expansion.

1. 1.

   Lake Tahoe: The dataset was acquired by the airborne Visible Infrared Imaging Spectrometer (AVIRIS) between 2014 and 2015 and contains hyperspectral images of six time phases. The size of the hyperspectral image in each time phase is $150\times 110$ pixels, and there are 173 bands after removing the water absorption band. The image contains three types of endmembers: water, soil, and vegetation, each of which produces distinct changes over time [36]. This sequence of images is shown in Fig. 3.

2. 2.

   Synthetic data1: The dataset contains six temporal hyperspectral synthetic images, each of size $50\times 50$ pixels, containing three endmembers, and three features of $\mathrm{L=224}$ bands randomly sampled from the USGS library as the reference endmember matrix. To model the endmember variability, the endmembers of each pixel are obtained by multiplying the reference signatures with a piecewise linear random scaling factor of the amplitude interval [0.85,1.15]. Local pixel mutation is added to $t\in\left\{2,3,4,5\right\}$ to better match the endmember changes of real scene, To simulate realistic scenarios, Gaussian noise with a SNR of 30dB was added [41].

3. 3.

   Synthetic data2: The second synthetic dataset contains 15 temporal hyperspectral images, each $50\times 50$ pixels in size. In order to introduce realistic spectral variability, the endmember features of each pixel and phase are randomly selected from the artificially extracted pure pixels of water, vegetation, soil, and road in Jasper Ridge HI, which contains 198 bands. To simulate realistic scenarios, Gaussian noise with a SNR of 30dB was added [41].

Due to the lack of real comparison results for the Lake Tahoe dataset, we only present its abundance map qualitatively, and the results are shown in Fig. 4. We alse show the results of multitemporal endmember estimation results in Fig. 7. It can be seen that the MUFormer proposed by us achieves better results in the unmixing of real multitemporal hyperspectral datasets. We can see from the results that FCLS has the worst performance, and the endmember is missing in multiple phases. There is also confusion between different endmembers. KalmanEM and OU have similar effects, but there are large errors in the abundance estimation of water and soil endmembers. Kalman filter has shortcomings in a large range of abundance changes and there are a lot of artifacts, which is because the Kalman filter sets the abundance to be continuous. When there is a large change in abundance, the effect will be affected. At the same time, in order to form a reasonable comparison, we add the DeepTrans model, which has superior performance on single-temporal hyperspectral image unmixing, as one of the comparison models. We can see that DeepTrans performs well in some individual temporal phases, but it does not consider the temporal information of multitemporal hyperspectral images, resulting in more erroneous unmixing regions. Finally, by comparing our proposed method with ReSUDNN, it can be seen that on the whole, both of them have achieved good results, and they are meticulous in the edge processing of endmembers. However, ReSUDNN has the phenomenon of misjudgment of endmembers in some phases, especially in the estimation of soil and vegetation.

Due to the lack of ground truth in Lake Tahoe, we can only qualitatively compare the abundance map results. In order to compare their unmixing accuracy more precisely, we also tested them on two synthetic datasets. We performed a quantitative comparison on synthetic datasets, the results of which are shown in Table I. The results of the abundance maps for each model are shown in Fig. 8. Since the synthetic dataset 2 contains 15 phases, we only compare the MUFormer estimated abundance plot with the real one, taking into account space constraints. The comparison results are shown in Fig. 9, and the quantitative results are shown in Table II.

1) Discussion: From Fig. 8 and 9, we can see that each model is able to capture the abundance change from time t=2 to time t=6, but all are affected by noise to varying degrees. FCLS has a large estimation error at time t=6, OU and KalmanEM are seriously polluted by noise, and DeepTrans, as a single-temporal unmixing model, performs better at the first two moments, but it does not consider the time information, so the effect is poor for the latter moments, and there is a case of wrong estimation. As a whole, ReSUDNN and our method are closer to the true abundance map. We can see that for long time sequence hyperspectral images, our model shows excellent results, which is because transformer has a good effect on long sequence information processing. At the same time, we also combine the GAM to process the abundance transformation information of adjacent time phases, which further enhances the effect of our model in processing multitemporal hyperspectral images. In Table I and II we compare the quantitative results of the individual algorithms, it can be seen that our proposed MUFormer achieves SOTA results in three indicators of $\mathbf{NRMSE_{A}}$, $\mathbf{NRMSE_{M}}$ and $\mathbf{SAM_{M}}$. Although OU has achieved good results on $\mathbf{NRMSE_{Y}}$, the reconstruction loss cannot explain the effect of abundance and endmember estimation, so we see that his results of $\mathbf{NRMSE_{A}}$, $\mathbf{NRMSE_{M}}$ and $\mathbf{SAM_{M}}$ are poor, and the abundance map obtained is also affected by a lot of noise. and at the same time, the running speed is several times higher than that of KalmanEM and ReSUDNN. It is worth mentioning that synthetic dataset 2 contains 15 temporal phases, which is more challenging than synthetic dataset 1. We can see that our method has a substantial improvement, which also highlights the potential of our model in processing long sequences of multitemporal hyperspectral images. In Fig. 5 and Fig. 6, we show the results of endmember estimation of synthetic dataset 1 and 2. It can be seen that our estimation results are very close to the real endmember, but we still need to further enhance the perception ability of the model for the change of phase at different times of the same endmember.

To fully verify the role played by our proposed CEM and GAM in the model, we launched an ablation study, the results of which are shown in Table III. From Table III, we can see that the model performance is improved to varying degrees when the CEM module and GAM are added separately, and the model effect is the best when the two modules are added at the same time. This is because through the GAM, we extract the spatial, global time and spectral information features of the multitemporal hyperspectral image. Combined with the CEM module to capture the changes between adjacent time phases, it is more conducive to our model to deal with multitemporal tasks. At the same time, as mentioned above, we processed images at different scales in CEM. In order to explore the influence of different scale processing on the unmixing effect, we carried out ablation experiments on the size of the convolution kernel, and the experimental results are shown in Table IV. From Table IV, we can see that the addition of convolution modules of different scales to $cls$ has different degrees of improvement, and the fusion of features at different scales can improve the robustness of the model. We compare the two fusion methods of addition and multiplication, and it can be seen that the multiplication effect of A1 and A2 is more obvious. We believe that multiplying them enables the model to pay more attention to the relationship between features at different scales, which helps to capture the information in the image more comprehensively, especially in the presence of multiple scale structures. Moreover, the correlation between two features can be emphasized by multiplication, and the weighted attention to different regions can also be realized, making the model more flexible to learn different degrees of attention to different locations.