A Spatial–Temporal Bayesian Deep Image Prior Model for Moderate Resolution Imaging Spectroradiometer Temporal Mixture Analysis

Abstract

Time-series remote sensing images are important in agricultural monitoring and investigation. However, most time-series data with high temporal resolution have the problem of insufficient spatial resolution which cannot meet the requirement of precision agriculture. The unmixing technique can obtain the object abundances with richer spatial information from the coarse-resolution images. Although the unmixing technique is widely used in hyperspectral data, it is insufficiently researched for time-series data. Temporal unmixing extends spectral unmixing to the time domain from the spectral domain, and describes the temporal characteristics rather than the spectral characteristics of different ground objects. Deep learning (DL) techniques have achieved promising performance for the unmixing problem in recent years, but there are still few studies on temporal mixture analysis (TMA), especially in the application of crop phenological monitoring. This paper presents a novel spatial–temporal deep image prior method based on a Bayesian framework (ST-Bdip), which innovatively combines the knowledge-driven TMA model and the DL-driven model. The normalized difference vegetation index (NDVI) time series of moderate resolution imaging spectroradiometer (MODIS) data is used as the object for TMA, while the extracted seasonal crop signatures and the fractional coverages are perceived as the temporal endmembers (tEMs) and corresponding abundances. The proposed ST-Bdip method mainly includes the following contributions. First, a deep image prior model based on U-Net architecture is designed to efficiently learn the spatial context information, which enhances the representation of abundance modeling compared to the traditional non-negative least squares algorithm. Second, The TMA model is incorporated into the U-Net training process to exploit the knowledge in the forward temporal model effectively. Third, the temporal noise heterogeneity in time-series images is considered in the model optimization process. Specifically, the anisotropic covariance matrix of observations from different time dimensions is modeled as a multivariate Gaussian distribution and incorporated into the calculation of the loss function. Fourth, the "purified means" approach is used to further optimize crop tEMs and the corresponding abundances. Finally, the expectation–maximization (EM) algorithm is designed to solve the maximum a posterior (MAP) problem of the model in the Bayesian framework. Experimental results on three synthetic datasets with different noise levels and two real MODIS datasets demonstrate the superiority of the proposed approach in comparison with seven traditional and advanced unmixing algorithms.

1. Introduction

With the development of remote sensing technology, different Earth observation satellite sensors are widely used. Remote sensing data acquired in the same area show multi-platform, multi-sensor, multi-resolution, multi-spectral, and multi-temporal characteristics [1]. Time-series remote sensing images are significant for agricultural investigation [2], change detection applications [3], and environmental monitoring [4,5]. In general, high temporal remote sensing images usually cover a larger spatial range with a large image footprint but exhibit lower spatial resolution [1]. Therefore, mixed pixels in medium and low-resolution remote sensing data are an inevitable problem [6,7]. Each hybrid pixel is usually a combination of multiple land cover classes [8], while both the spectral reflectances and temporal growth curves of different classes are inconsistent. The spectral signal and temporal signal of a mixed pixel are composed of different “heterogeneous” components, while the spectral signal and temporal signal of a pure pixel represent only one class [9]. Compared with high-resolution remote sensing data, the coarse spatial resolution remote sensing data after the unmixing operation can provide more abundant temporal information and a wider coverage area. Therefore, searching for more effective unmixing methods to address the “mixed” problem in medium and low-resolution images remains a key concern.

Recently, traditional unmixing algorithms have usually been based on the linear spectral mixture model (LSMM) [10,11,12], which models the spectral reflectance of a single pixel as a linear combination of endmember reflectances, weighted by the abundance of each endmember within the pixel [6]. Liu et al. used LSMM for the unmixing of cotton, corn, tomato, and soil, and the experiments showed that LSMM achieves high efficiency, high precision, and strong adaptability [11]. The important process of unmixing can be divided into two tasks: endmember extraction and abundance estimation [13]. At present, most geometrical-based methods are commonly used for endmember extraction, including the vertex component analysis (VCA) [14], N-FINDR [15], and pixel purity index (PPI) [16]. Some researchers used cluster centers of Kmeans or KPmeans as endmembers, which employed the “Purified means” method, and also showed good performance by effective constraints on the derivation of endmembers [17,18]. However, it is difficult for these geometrical-based methods to extract accurate endmembers from highly mixed data because pure spectral features are not always available in reality [19]. In addition, fully constrained least squares (FCLS) [20] and sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL) [21] algorithms are commonly used for abundance inversion [22], which usually need to assume that the pure endmember signatures exist in hyperspectral images. These algorithms are poorly inclusive of advanced modeling to solve the large-scale spatial correlation effects in remote sensing data. In recent years, deep learning (DL) methods have shown great capability and potential for spectral unmixing (SU) applications. Most of these DL-based methods are based on autoencoder (AE) networks [23,24,25] or convolutional neural networks (CNNs) [26,27,28]. Furthermore, the U-Net network, which uses fully convolutional layers, can better capture the large-scale spatial correlation effect of remote sensing data [29,30]. However, if there is ambiguous noise in images or the number of endmembers is incorrectly estimated, the unmixing performance usually drops drastically. Many DL-based methods consider the denoising algorithm as a module of the iterative optimization process in order to address the problem. For example, an untied denoising autoencoder with sparsity (uDAS) is proposed to formulate the unmixing problem, which improves the performance of abundance estimation [24]. And Zhao et al. propose a robust plug-and-play (PnP) priors framework for hyperspectral unmixing, which plugs a variety of denoisers to replace prior models and avoids designing regularizers [31]. However, the limited spectral information of the images and the restriction of sensor spatial resolution may directly lead to the low accuracy of ground object recognition and extraction [32]. And most of them are designed for the remote sensing data of a single day, mainly applied in SU application currently [27,33,34,35]. Few studies have focused on the temporal unmixing of continuous temporal remote sensing data.

Moderate Resolution Imaging Spectroradiometer (MODIS) data are widely used in large-scale monitoring of land cover types [36,37,38]. Many studies analyze phenological characteristics of remote sensing vegetation indices (VIs) of different crops in different seasons to realize the recognition of crop types [39,40,41,42]. However, the application of MODIS data is limited due to the low spatial resolution from multiple interferences, intimate mixtures, and resolution trade-offs. In recent years, the temporal mixture analysis (TMA) method [43] and TMA-based methods [44,45,46] are proposed to provide a new unmixing pattern for hypertemporal images (HTIs), which treat the temporal profiles of the mixed pixel as temporal endmembers (tEMs). TMA is an extension of spectral mixture analysis (SMA). And the algebraic expressions of TMA and SMA are identical in theory, except that TMA applies to the temporal spectra rather than the reflectance spectra. Different from traditional spectral curves as unmixing elements, the curves of VIs as the unmixing elements of the mixed pixels have developed and provided higher accuracies of endmember extraction and abundance estimation [5]. Many TMA methods based on the time-series data have achieved satisfactory results [6,46], which assume fractions of the main ground land cover types in the study area would not change significantly in a period of time [47]. Chi et al. apply a machine learning technique for TMA to study the temporal characteristics of Antarctic sea ice [48]. And an improved phenology-based temporal mixture analysis (PTMA) method is proposed to identify early, middle and late rice planting patterns and has been proven to be effective for rice localization in hilly areas of small-size paddy rice fields [49]. The TMA of the MODIS data can fully utilize the advantages of spatial, spectral, and temporal features to obtain high spatial and temporal resolution data. Therefore, it is significant to develop TMA methods for time-series data to monitor the dynamic changes in large areas while the global surface is changing every day [5].

Although analyzing HTIs with DL-based methods can effectively alleviate the limitations of single views with unsatisfactory accuracy for ground object recognition, it also poses new challenges. The variation in light conditions and the difference in spectral absorption properties from different times when remotely sensed images are generated cause the noise levels of the images at different temporal dimensions to vary considerably. Therefore, time-series images have different noise variances from different time periods [50], which affects the correct characterization of time series, and limits the temporal application of remote sensing images. However, most studies have not considered the problem of heterogeneous noises in the temporal dimension.

To address the low spatial resolution of HTIs and the problems summarized above, this paper develops a novel s-temporal Bayesian deep image prior (ST-Bdip) method based on a Bayesian framework that combines the physically-driven TMA model and the DL-driven model for TU of MODIS crop normalized difference vegetation index (NDVI) data. This method has the following main characteristics:

• The U-Net convolutional architecture is specifically designed to estimate crop coverage of the MODIS time series NDVI data. To be specific, this model uses multiple downsampling layers and upsampling layers to obtain the spatial and temporal context information of MODIS NDVI data, so the abundance of tEMs could be estimated efficiently. And the deep image prior (DIP) is utilized to account for spatial correlation within the abundance field.
• The TMA model is incorporated into the U-Net training process, which can effectively utilize the prior knowledge of the physical imaging process in the forward model. The linear relationship between tEMs and corresponding abundances is modeled, which improves the U-Net model to better understand the information on different land cover types in the HTIs.
• To solve the temporal noise variance in HTIs, the heterogeneous noise in NDVI data distributed in different periods is modeled by the multivariate Gaussian distribution.
• The ST-Bdip approach estimates the tEMs using the “Purified means” method, which can be interpreted as a conditional distribution of the tEMs by the obtained abundances. And the above components are integrated into the Bayesian framework. The expectation–maximization (EM) algorithm is used to solve the maximum a posteriori (MAP) problem.
• We compare different traditional unmixing methods such as N-FINDR, PPI, SUnSAL, FCLS, and KPmeans, and some DL-based unmixing methods such as uDAs, PnP with the proposed ST-Bdip method in our experiments. And the performances of models are evaluated from the extracted tEMs and the estimated abundances.

The rest of this paper is organized as follows. In Section 2, dataset and preprocessing are introduced. In Section 3, we discuss details of our model for time-series image unmixing. Results and analysis of experiments are detailed in Section 4, followed by the discussion in Section 5 and the conclusion in Section 6.

2. Dataset and Preprocessing

2.1. Simulated Dataset

We simulate the synthetic time-series image using four endmember reference curves, with the size of $102\times 102$ pixels. First, referring to the corresponding cropland data layer (CDL) product, we use 250 m MOD13Q1 NDVI data to select NDVI curves of some relatively purer pixels of four crops and calculate the average values as the reference tEMs using the ENVI tool. The tEMs here represent the NDVI profile curves, rather than the traditional spectral curves. Each year, the MOD13Q1 NDVI data comprise 23 images (days 1–363, at 16-day intervals). Figure 1a–d visually presents the selected endmember references for four distinct crops, specifically the NDVI growth curves corresponding to Grass, Corn, Winter Wheat, and Alfalfa. Second, with the rapid development of remote sensing technology, HTIs will be more and more widely used. Therefore, considering the continuity of the time-series images and the abrupt errors of the data, it is reliable to simulate the data with higher temporal resolution and satisfies the realistic scene. Most studies demonstrated that the spline interpolation algorithm has stronger adaptability and parameter simplicity than other algorithms for reconstructing NDVI time series [51,52,53]. Thus, four tEMs with 115 temporal bands are generated symbolically, as shown in Figure 1e–h, in which “Grass(1–115)” represents that the temporal dimension of the original one-year grass NDVI curves is increased to 115 bands after interpolation processing, and the other crops are also adopted using this method. Third, abundance maps are generated by first being divided into $14\times 14$ sized homogeneous blocks of one of the four tEMs, then degrading the blocks by applying a spatial low-pass filter of $11\times 11$. And Figure 1i–l shows four simulated abundance maps with the size of $102\times 102$ pixels. Finally, each pixel in the simulated sequence data is a high blend of four tEMs based on the linear TMA model using the product of tEMs and abundance created previously.

Furthermore, due to the influence of the atmosphere and sensor noise in the process of data acquisition [15], zero-mean Gaussian noise with different noise variances in each temporal band is added to the synthetic data using different signal-to-noise ratio (SNR) values from 10 to 30 dB. These highly mixed simulated time-series images are obtained to verify the TU performance of the proposed algorithm under different noise levels, as shown in Figure 1m–o (the 100th band selected as an example). Figure 1p shows the hard label of the simulated data.

2.2. MODIS Dataset

2.2.1. Study Areas

Two real datasets are used in our experiments: respectively located in Scott County and Sherman County, the MODIS and NDVI datasets are publicly available, as shown in Figure 2. Scott County and Sherman County are both located in the western of Kansas, between ${38}^{\circ }{14}^{\prime }\sim {38}^{\circ }{42}^{\prime }$ north latitude and ${100}^{\circ }{39}^{\prime }\sim {101}^{\circ }{7}^{\prime }$ west longitude, and ${38}^{\circ }{14}^{\prime }\sim {38}^{\circ }{42}^{\prime }$ north latitude and ${100}^{\circ }{39}^{\prime }\sim {101}^{\circ }{7}^{\prime }$ west longitude, respectively. The main crops in these two counties are corn, soybean, winter wheat, etc. And the NDVI curve of each crop reflects unique phenological characteristics, with specific planting dates and seasonal growth cycles. Five typical land cover types are selected in our experiments, including corn, winter wheat, grass, fallow cropland, and sorghum. And all other crops are uniformly classified into the “Others” category.

2.2.2. Data Preprocessing

We downloaded the MODIS NDVI time-series data of Scott County and Sherman County in the 2018 year from Level-1 and the Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC) website, which is a sinusoidal projection coordinate system, with a spatial resolution of 250 m and a temporal resolution of 16 days, i.e., 23 available images each year [46,54,55].

Table 1. The data and pertinent information utilized in this study.

Data	Spatial Resolution	Time Period	Projection	Download Source
MOD13Q1	250 m	2018/01/01	Sinusoidal projection	https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 1 June 2023)
CDL	30 m	2018–2018/12/31	Albers Equal Area Conic projection	https://nassgeodata.gmu.edu/CropScape/ (accessed on 1 June 2023)

provides a comprehensive summary of the datasets utilized in this study, including MODIS NDVI time-series data and CDL reference data. The MODIS Reprojection Tool (MRT) was used for data preprocessing, including projection transformation, region of interest extraction, and superposition analysis. All MODIS tiles here were reprojected from Sinusoidal projection to Albers Conical Equal Area projection, to achieve alignment with the CDL dataset [56]. Subsequently, the transformed MODIS data were superimposed, clipped, and processed using the shapefiles corresponding to the two counties. Finally, the trimmed size of Scott County is $188\times 143$ pixels, and the size of Sherman County is $185\times 216$ pixels.

2.2.3. Reference Selection

The creation of reference labels of the datasets mainly includes endmember selection and abundance calculation. Firstly, the prior knowledge of the image, CDL, and high-resolution images of Google Earth were used to effectively determine the number of tEMs in the study areas. Secondly, tEMs were mainly obtained from the NDVI profile curves of pixels composed of a single coverage type [6]. We collected uniform and homogeneous pixels representing a similar object in the MODIS NDVI image of the study area, while the pixel reflects the same crop in the neighborhood $8\times 8$ patch from the corresponding CDL map. The average values of the homogeneous pixels of the same land cover acquired from different regions were calculated as the final reference tEMs, shown in Figure 3. Thirdly, the reference abundance was mainly calculated based on the CDL map, as shown in Figure 4. For example, a MODIS NDVI pixel with 250 m spatial resolution roughly corresponds to 64 CDL pixels with 30 m spatial resolution, and the probabilities of respective tEMs were calculated, so that the reference abundances corresponding to different land cover types in the study area can be obtained.

3. Methodology

In this section, we propose an temporal unmixing method ST-Bdip based on the Bayesian framework, which innovatively combines the knowledge-driven TMA model and the DL-driven model. This approach considers the influence of spatial correlation and temporal noise variance simultaneously. The specific architecture is depicted in Figure 5. The ST-Bdip method uses the TMA model as the forward model for temporal unmixing, which can mainly be divided into two steps: (1) identify the tEMs of crops; (2) decompose each pixel into crop tEMs and corresponding crop coverage. In addition, the ST-Bdip method uses the inverse model based on U-Net architecture [57] as the principal part for a prior estimate of abundances. The deep image prior model of U-Net architecture can effectively learn the spatial context information, which enhances the representation ability of abundance modeling compared with the traditional non-negative least squares algorithm. Subsequently, the “Purified means” concept and the expectation–maximization (EM) algorithm [58] are used for iterative optimization of tEMs and abundances.

3.1. Temporal Mixture Analysis

With a fully constrained linear TMA, the hybrid pixel can be regarded as the mixture of several specific crops, that is, the multi-temporal NDVI profiles, instead of the traditional reflectance spectra, can be identified for the linear combination of these crop tEMs and the corresponding crop coverages. By selecting the crop types and the number of crop tEMs, TMA can be mathematically modeled as Equation (1).

$$NDV{I}_{mix}=\sum _{i=1}^N{f}_i\times NDV{I}_i+\epsilon$$

(1)

$$\mathrm{Subject}\mathrm{to}\sum _{i=1}^N{f}_i=1\mathrm{and}{f}_i\ge 0$$

(2)

where $NDV{I}_{mix}$ is the NDVI value of a mixed pixel in the time-series image, N is the number of crop tEMs, $NDV{I}_i$ represents the NDVI value of the ith crop tEM, and $f_i$ is the fraction of crop tEM i, and $\epsilon$ is the residual. As defined in Equation (2), the abundances of the crop tEMs should satisfy two constraints: (1) non-negative and (2) sum-to-one.

Temporal unmixing aims to extract sequential tEMs $NDV{I}_i$ and estimate the corresponding abundance $f_i$ for a given time-series NDVI value $NDV{I}_{mix}$ of a mixed pixel. The proposed ST-Bdip method is a deep convolution model based on the Bayesian theorem, which exploits s-temporal deep image prior information. According to the Bayesian theorem, a posteriori distribution of unknown quantities can be established using the original prior distribution of noiseless data, which can be formulated as

$$p(NDV{I}_i, f_i\mid NDV{I}_{mix})\propto p(NDV{I}_{mix}\mid NDV{I}_i, f_i)p(NDV{I}_i\mid f_i)$$

(3)

Then, the posterior distribution can be maximized, i.e., the maximum a posterior (MAP) problem, so that the $NDV{I}_i$ and $f_i$ can be estimated. We write $\left\{NDV{I}_i, f_i\right\}$ as $\mathit{\theta }$ in the following formula:

$$\mathit{\theta }=\underset{\mathit{\theta }}{argmax}\left\{p\left(NDV{I}_i, f_i\mid NDV{I}_{mix}\right)\right\}$$

(4)

Theoretically, the maximization of $p\left(NDV{I}_i,f_i\mid NDV{I}_{mix}\right)$ is also the minimization of the negative log-likelihood representation, i.e.,

$$\begin{array}{c}\hfill \mathit{\theta }=\underset{\mathit{\theta }}{argmin}\left\{-{log}_{}p\left(NDV{I}_i, f_i\mid NDV{I}_{mix}\right)\right\}\propto \\ \hfill \underset{\mathit{\theta }}{argmin}\left\{-{log}_{}p\left(NDV{I}_{mix}\mid NDV{I}_i, f_i\right)-{log}_{}p\left(NDV{I}_i\mid f_i\right)-{log}_{}p\left(f_i\right)\right\}\end{array}$$

(5)

3.2. Noise Heterogeneity

Generally, most researchers assume that the noise existing in the images is homogeneous noise [59,60]. However, remote sensing time-series images are affected by different levels of noise in both temporal and spectral dimensions due to the illumination variation and difference acquisition time of the dataset. Therefore, we design the model to adapt to the impact of noise heterogeneity in different temporal periods. Based on the least squares optimization model, the term of the data likelihood function is modeled and solved by anisotropic multivariate Gaussian distribution, and the formula $p\left(NDV{I}_{mix}\mid NDV{I}_i, f_i\right)$ can be expressed as follows:

$$p\left(NDV{I}_{mix}\mid NDV{I}_i, f_i\right)={\prod }_1^{N}exp\left({\left(NDV{I}_{mix}-NDV{I}_i{f}_i\right)}^T{\Lambda }^{-1}\left(NDV{I}_{mix}-NDV{I}_i{f}_i\right)\right)$$

(6)

where ${\left(NDV{I}_{mix}-NDV{I}_i{f}_i\right)}^T{\Lambda }^{-1}\left(NDV{I}_{mix}-NDV{I}_i{f}_i\right)$ represents the reconstruction error of image, and $\Lambda$ represents the diagonal matrix of noise covariance $\epsilon$, the formula of which is as follows:

$$\Lambda =\left(\begin{array}{cccc}{\epsilon }_1^2& & & \\ & {\epsilon }_2^2& & \\ & & \dots & \\ & & & {\epsilon }_T^2\end{array}\right)$$

(7)

where ${\epsilon }_t^2,t=1,2,\dots ,T$ represents the noise variance of the tth temporal band. Compared with most existing unmixing methods which assume the same noise level in different temporal or hyperspectral bands, our proposed method considers the heterogeneity of noise in different temporal bands, so as to obtain better noise estimation and unmixing results.

3.3. ST-Bdip Framework

As shown in Figure 5, the U-Net framework combining skip connections in the model is used to efficiently extract the prior spatial information of abundance in the large-scale range area. The model consists of convolution, upsampling, downsampling, and nonlinear activation layers. In addition, the skip connection is used to restrain overfitting and promote the convergence of the model.

The last activation layer of the inverse model is set to the softmax layer because the obtained abundance needs to satisfy the “sum-to-one” and “non-negative” constraints, as shown in Equation (2). In conclusion, the model for estimating the tEMs and abundances contains the following major steps:

• The estimated abundances are acquired from an inverse model based on the U-Net framework combining skip connections using the MODIS NDVI time-series data.
• The initial tEMs are obtained based on the VCA algorithm.
• According to the TMA, the time-series image is reconstructed using the abundances from step 1 and tEMs from step 2.
• The multivariate Gaussian distribution with an anisotropic covariance matrix to represent the conditional temporal distribution is used. And the Mahalanobis distance-based loss function is calculated for the model optimization to solve the noise heterogeneity effect in remote sensing images. Thus, the model parameters are optimized, and the new abundance is obtained.
• The “Purified means” method is used to build and further optimize tEMs.
• The tEMs and abundances are constantly updated through the above step 4 and step 5 until the optimal solution is obtained from EM iteration.

4. Experiments and Results

4.1. Comparison Methods and Parameters

The experiments are conducted with one synthetic dataset and two real MODIS NDVI images. Several traditional and state-of-the-art unmixing algorithms are used for comparison with the proposed ST-Bdip method, e.g., N-FINDR, PPI, VCA-S, VCA-F, uDAs, PnP, and KPmeans. N-FINDR and PPI algorithms are two unsupervised endmember extraction methods widely used in the field of hyperspectral unmixing [61,62,63,64]. The N-FINDR algorithm iteratively searches for endmembers by seeking the direction to maximize the abundance of the selected pixel in the data hyperspace, until all the required endmembers are identified, and then updates the abundances of the selected pixel according to the identified endmembers [15]. Differently, the PPI algorithm directly calculates the pixel purity index of each pixel in the image, selects the pixel with the highest purity index as the endmember, and estimates the corresponding abundances [65]. In addition, both VCA-S and VCA-F methods use the VCA algorithm to extract tEMs, but VCA-S uses the SUnSAL method to estimate the abundance, while VCA-F uses the FCLS method to estimate the abundance. Specifically, SUnSAL is a sparse unmixing algorithm, which initiates the abundance matrix either randomly or based on specific criteria, and subsequently updates the abundance matrix using sparse optimization algorithms [66]. The FCLS method solves the fully constrained least squares equation to minimize the error between the observed data and the linear combination of endmembers, yielding the abundance matrix [22]. Besides conventional unmixing methods, we also compare DL-based approaches that use deep learning techniques to estimate the results in the unmixing application, including uDA and PnP. The uDAs method leverages unsupervised learning algorithms like autoencoders or variational autoencoders to acquire the intrinsic representation of hyperspectral data, enabling the estimation of endmembers and abundances [24]. The PnP method integrates endmember information with prior constraints and employs iterative algorithms to solve optimization problems to obtain abundance maps [31]. KPmeans is a variant of the traditional K-means algorithm that uses kernel functions to process nonlinear data. The conversion data of the K-means clustering algorithm is used to identify the endmembers, and then the corresponding abundances are estimated according to the concept of “Purified means” [17]. The proposed method ST-Bdip refers to the “Purified means” idea of the KPmeans method for iterative updating.

For the proposed ST-Bdip method, we empirically set the number of EM iterations as 50 for simulated data and 1000 for real MODIS NDVI datasets. The experiments are implemented through the MATLAB R2018a platform using codes of the PnP and uDAs methods run with an Intel Xeon Silver 4110 CPU with 2.10 GHz. Other algorithms run on the Pytorch toolbox with an NVIDIA GeForce RTX 2080 Ti GPU.

4.2. Numerical Measures

We employ the following seven evaluation metrics in the experiments: the temporal angle distance (TAD) and the temporal information divergence (TID) to test the performance of tEMs extraction, while the expression formulas are similar to spectral angle distance and the spectral information divergence; the abundance angle distance (AAD), the abundance information divergence (AID), structural similarity (SSIM), and root mean square error between predicted abundances and reference abundance (RMSE-S) to test the performance of abundance estimation; root mean square error of reconstruction image and observation image (RMSE-X) to test the performance of image reconstruction. Their specific representations are as follows:

$$TAD/AAD={cos}^{-1}\left(\left(x^T\widehat{x}\right)/(\parallel \widehat{x}\parallel \parallel x\parallel )\right)$$

(8)

$$TID/AID=D(x/\widehat{x})+D(\widehat{x}/x)$$

(9)

$$RMSE=\sqrt{\frac{1}{N}\sum _{n=1}^N{∥x_n-\widehat{x_n}∥}^2}$$

(10)

where x is the reference tEM, true abundance, or observation HTI image, and $\widehat{x}$ represents reconstruction tEM, reconstruction abundance, or reconstruction HTI image. $D(a/b)$ calculates the relative entropy between a and b. SSIM [67] is an evaluation index related to visual effects, which evaluates image quality from contrast, brightness, and structure dimensions. We calculate the SSIM metric using the compare_ssim module from the scikit-image package of Python.

4.3. Test on Simulated Dataset

Table 2. Mean TAD, AAD, TID, AID, 1-SSIM, and the variances of them obtained from 10 independent runs by different methods using the simulated data over SNR from 10 to 30 dB. The best results are in bold.

SNR = 10, 4 tEMs, 102 × 102 size, 115 bands
Method	TAD	AAD	TID	AID	1-SSIM
N-FINDR	0.3096 ± 0.02	0.5818 ± 0.03	0.1320 ± 0.02	4.5547 ± 0.36	0.5473
PPI	0.3006 ± 0.01	0.7055 ± 0.03	0.1101 ± 0.01	5.5562 ± 0.33	0.7019
VCA-S	0.1108 ± 0.03	0.6669 ± 0.18	0.0170 ± 0.01	5.1622 ± 4.00	0.4683
VCA-F	0.0953 ± 0.04	0.5936 ± 0.13	0.0137 ± 0.01	6.3974 ± 3.47	0.3889
uDAs	0.1531 ± 0.04	0.5420 ± 0.09	0.0355 ± 0.02	4.1086 ± 1.22	0.3990
PnP	0.0936 ± 0.04	0.5186 ± 0.13	0.0123 ± 0.01	4.4523 ± 1.97	0.2872
KPmeans	0.1286 ± 0.04	0.5098 ± 0.09	0.0241 ± 0.01	3.8117 ± 1.30	0.3343
ST-Bdip	0.0816 ± 0.04	0.3339 ± 0.15	0.0110 ± 0.01	2.1410 ± 1.10	0.1387
SNR = 20, 4 tEMs, 102 × 102 size, 115 bands
Method	TAD	AAD	TID	AID	1-SSIM
N-FINDR	0.1849 ± 0.01	0.5896 ± 0.05	0.0404 ± 0.00	5.0966 ± 0.72	0.5711
PPI	0.1073 ± 0.01	0.2756 ± 0.02	0.0132 ± 0.00	1.7968 ± 0.25	0.1669
VCA-S	0.0672 ± 0.03	0.3803 ± 0.19	0.0056 ± 0.00	1.7978 ± 2.55	0.1965
VCA-F	0.0648 ± 0.03	0.3151 ± 0.10	0.0052 ± 0.01	2.0104 ± 1.23	0.1961
uDAs	0.0612 ± 0.04	0.2932 ± 0.12	0.0072 ± 0.01	1.9552 ± 1.68	0.1531
PnP	0.0655 ± 0.03	0.3035 ± 0.10	0.0055 ± 0.01	2.1320 ± 1.34	0.1763
KPmeans	0.0376 ± 0.02	0.2341 ± 0.13	0.0025 ± 0.00	1.5301 ± 1.88	0.1217
ST-Bdip	0.0354 ± 0.02	0.1640 ± 0.08	0.0032 ± 0.00	1.0483 ± 0.53	0.0658
SNR = 30, 4 tEMs, 102 × 102 size, 115 bands
Method	TAD	AAD	TID	AID	1-SSIM
N-FINDR	0.1758 ± 0.00	0.6864 ± 0.08	0.0391 ± 0.00	9.4483 ± 2.97	0.4656
PPI	0.0577 ± 0.00	0.1632 ± 0.03	0.0036 ± 0.00	0.5543 ± 0.24	0.0777
VCA-S	0.0546 ± 0.02	0.2612 ± 0.18	0.0036 ± 0.00	1.0269 ± 2.36	0.1553
VCA-F	0.0537 ± 0.02	0.2264 ± 0.07	0.0036 ± 0.00	1.2240 ± 0.87	0.1064
uDAs	0.0554 ± 0.03	0.2874 ± 0.12	0.0053 ± 0.01	1.8998 ± 1.79	0.1278
PnP	0.0530 ± 0.01	0.2255 ± 0.05	0.0031 ± 0.00	1.2114 ± 0.57	0.1435
KPmeans	0.0306 ± 0.01	0.1044 ± 0.04	0.0013 ± 0.00	0.3667 ± 0.28	0.0539
ST-Bdip	0.0107 ± 0.01	0.0536 ± 0.02	0.0003 ± 0.00	0.2663 ± 0.21	0.0120

quantitatively presents the performance of different unmixing methods using the mean TAD, AAD, TID, AID, and SSIM results of four crops on synthetic data. It is immediately observed that our proposed ST-Bdip method achieves optimal performance under all levels of noise. Table 2 also shows that the larger the SNR of the image, the better the result of unmixing. The SSIM values indicate that the abundance obtained by our proposed ST-Bdip method is most similar to the reference abundance. Different metrics are visually displayed in Figure 6 and Figure 7 with different SNR noises. In Table 2 and Figure 6, we can see that the variance in the ST-Bdip method does not perform the best when the input image has a small SNR, but the overall accuracy is still optimal. In addition, it shows better robustness when the SNR is large, which is similar to other methods. In Figure 7, we can see more intuitively that the ST-Bdip method has the smallest reconstruction error of the abundance, exceeding all other traditional and advanced methods, and the reconstruction errors of the image show few differences in all methods. However, the smaller the error value of the reconstructed image does not signify the better the result due to the existence of the temporal noise.

Figure 8 shows the tEM curves of four crops extracted by all methods when SNR noise equals 20 dB. It can be seen that the tEMs extracted by N-FINDR and PPI methods have more noise. Especially, the N-FINDR method deviates more from the reference tEMs of grass, corn, and winter wheat crops. There is little difference between the other methods, and ST-Bdip is closer to the ground truth for all crops. The abundance maps of the four materials for SNR noises equal to 10 dB and 20 dB are shown in Figure 9. We can find that the abundances of the VCA-S method have brighter results because they have no constraint of sum-to-one. In addition, the PPI method fails to unmix this dataset with small SNR, as shown in Figure 9b. Although the result is better when SNR equals 20 dB in Figure 9k, there is still no ideal for “Grass” and “Alfalfa” materials. Overall, the proposed ST-Bdip method shows a minimum salt and pepper noise of abundance, which is closest to the reference abundances.

4.4. Test on MODIS NDVI Dataset

The goal of the experiments is to obtain the abundances at 250 m resolution according to the extracted tEMs of the six main cover types in Scott and Sherman counties, which is helpful for the determination of crop spatial mapping and yield estimation in a higher resolution subsequently. We evaluated the average TAD, AAD, TID, AID, and RMSEs of abundances (represented by RMSE-S) and reconstruction images (represented by RMSE-X) from 20 independent runs by different methods using the MODIS NDVI products in two study areas, as shown in

Table 3. Mean TAD, AAD, TID, AID, their variances, the RMSE-S and RMSE-X obtained from 20 independent runs by different methods using the MODIS NDVI products in Scott County. The best results are in bold.

Methods	TAD	AAD	TID	AID	RMSE-S	RMSE-X
N-FINDR	0.2757 ± 0.05	0.9939 ± 0.04	0.3208 ± 0.01	15.0797 ± 0.99	0.2911	0.0576
PPI	0.3007 ± 0.06	1.0194 ± 0.09	0.9141 ± 0.44	16.3344 ± 2.42	0.2969	0.0671
VCA-S	0.2463 ± 0.04	1.0092 ± 0.09	0.0796 ± 0.03	49.9739 ± 16.47	0.2865	0.0359
VCA-F	0.2541 ± 0.05	0.9612 ± 0.05	0.0832 ± 0.03	14.8227 ± 1.90	0.3019	0.0473
uDAs	0.2943 ± 0.09	0.9452 ± 0.07	0.3228 ± 0.55	15.1977 ± 2.02	0.3082	0.0477
PnP	0.2644 ± 0.08	0.9668 ± 0.09	0.0871 ± 0.05	16.2128 ± 3.31	0.2804	0.0559
KPmeans	0.2557 ± 0.06	0.8888 ± 0.05	0.1128 ± 0.13	13.6755 ± 0.78	0.2606	0.0403
ST-Bdip	0.2300 ± 0.05	0.8833 ± 0.06	0.0710 ± 0.05	13.1397 ± 1.21	0.2382	0.0451

and

Table 4. Mean TAD, AAD, TID, AID, their variances, the RMSE-S and RMSE-X obtained from 20 independent runs by different methods using the MODIS NDVI products in Sherman County. The best results are in bold.

Methods	TAD	AAD	TID	AID	RMSE-S	RMSE-X
N-FINDR	0.3067 ± 0.02	1.0186 ± 0.08	0.1263 ± 0.01	15.6274 ± 2.61	0.3103	0.0795
PPI	0.3185 ± 0.02	1.0446 ± 0.09	0.1407 ± 0.02	17.5854 ± 8.49	0.3201	0.1177
VCA-S	0.3312 ± 0.10	1.0687 ± 0.13	0.1558 ± 0.11	47.1089 ± 26.49	0.3348	0.0415
VCA-F	0.3173 ± 0.07	0.9284 ± 0.06	0.1431 ± 0.06	14.5976 ± 1.16	0.2781	0.0543
uDAs	0.3296 ± 0.07	0.8975 ± 0.05	0.7204 ± 0.75	13.6823 ± 1.09	0.2799	0.0516
PnP	0.2915 ± 0.02	0.9298 ± 0.02	0.1201 ± 0.01	16.6188 ± 1.30	0.2914	0.0598
KPmeans	0.3878 ± 0.06	0.9071 ± 0.02	0.6087 ± 0.17	14.8144 ± 0.85	0.2986	0.0408
ST-Bdip	0.2612 ± 0.03	0.8639 ± 0.03	0.1016 ± 0.06	12.7759 ± 1.57	0.2420	0.0556

. The proposed ST-Bdip method shows optimal results in all metrics except RMSE-X. The TAD value of the ST-Bdip method is 0.2300 in Scott County and 0.2612 in Sherman County, and the AAD value is 0.8833 in Scott County and 0.8639 in Sherman County. Empirically, abundance estimation is more valuable and essential than tEM extraction. The results of abundance estimation in Table 3 and Table 4 show that the KPMeans method also has a good performance in Scott County with the AAD value of 0.8888. In Sherman County, the uDAs and KPMeans methods perform well with AAD values of 0.8975 and 0.9071, respectively.

For illustrative purposes, Figure 10, Figure 11, Figure 12 and Figure 13 show some of the classification maps corresponding to the experiments reported in Table 3 and Table 4. There is no exact tEM of the “Other” class while it is calculated by the average of all the other classes, so it does not represent a certain type of land class, and the results do not reflect the superiority of the methods. And five crops (“Corn”, “Sorghum”, “Winter Wheat”, “Fallow Cropland”, and “Grass”) are mainly presented. The comparison of the five tEMs extracted by the eight methods and the reference tEMs (red lines) for Scott and Sherman counties is shown in Figure 10 and Figure 12. Correspondingly, the abundance maps of the two counties are displayed in Figure 11 and Figure 13, respectively. As we can see, compared with simulated data, real MODIS data greatly increase the difficulty of temporal unmixing. The tEMs estimated by the ST-Bdip method are more consistent with the reference tEMs than other methods, which proves the effectiveness of the TMA for the real dataset. Although the NDVI curves of corn and sorghum crops are similar, most methods extract good tEMs of them except PPI, VCA-S, and PnP methods.

5. Discussion

This paper presented a deep s-temporal method ST-Bdip based on the Bayesian framework for temporal unmixing. At present, hyperspectral remote sensing data are mainly used in most crop unmixing applications [8,68]. However, hyperspectral images have some limitations such as data redundancy, small coverage area, and challenges in data acquisition [69]. Given the time-series data as a crucial indicator of crop growth trends, this paper focuses on expanding the endmembers in unmixing from spectral curves to time-series curves and extending the application of unmixing algorithm from hyperspectral images to time-series images, providing different exploration directions to innovatively calculate the individual crop fractions in coarse resolution pixels.

5.1. Advantages of the Proposed ST-Bdip Method

We designed the experiment of traditional and advanced algorithms of hyperspectral unmixing being innovatively applied to the temporal unmixing domain. Different methods show different generalization effects. As shown in Table 2, Table 3 and Table 4, our proposed ST-Bdip method achieves relatively better performance than other comparison methods, both in terms of the extraction of tEMs and the estimation of abundances. Surprisingly, although the “Corn” and “Sorghum” classes show similar growth curves, most methods obtain distinguishable results of them in Figure 10 and Figure 12a,b. It is found that the mistakes of tEMs extraction are mainly concentrated in the “Fallow Cropland” and “Grass” classes. In addition, most existing methods assumed that the noise in the image is homogeneous noise, which cannot reflect the noise distribution in the real remote sensing image. The corresponding solution was designed in our proposed method to solve the problem of heterogeneity noise in the different temporal dimensions. In this paper, three simulated data with different noise levels and two MODIS NDVI data are used to discuss the superiority of the ST-Bdip method in the application of temporal unmixing. The results of the experiments show that most of the traditional and advanced methods have a good performance on the simulated data, and the complexity of real data brought great difficulty to temporal unmixing.

Overall, the abundance maps obtained by the ST-Bdip method are closer to the ground truth maps in terms of both the intense brightness and the structural characteristics, both in Scott County and Sherman County. The proposed ST-Bdip model demonstrates promising results in all materials, particularly in the “Fallow Cropland” class of Scott County and the “Grass” class of Sherman County versus other methods. Essentially, all of the methods misclassify the “Winter Wheat" class of Scott County in the third row of Figure 11. Comparing the CDL map with the predicted abundance maps of this material, we find that most of the highlighted abundances identified by these methods appeared in CDL as the “Alfalfa” class, which is not selected in our experiments because it is a relatively small number in the study areas. It is possible that the extracted tEMs have been biased by these methods, which results in the misclassification of the abundances. Although the uDAs and PnP methods consider the effect of image noise, the results of MODIS data unmixing are still inferior to the proposed ST-Bdip method, which has simulated heterogeneous noise in the temporal dimension.

5.2. Data Complexity Analysis

Figure 8 and Figure 9 reflect that the proposed ST-Bdip method and comparison methods can perform well in the simulation datasets without large deviation. However, the experimental results on the real MODIS data in Figure 10, Figure 11, Figure 12 and Figure 13 show that there are significant differences between different methods and different crops. We conducted a comprehensive analysis and concluded that the above results may be attributed to the complexity of real MODIS data, which can be discussed in the following aspects. First, the simulated datasets contain fewer crop categories and have lower data complexity than the real MODIS data. Second, real MODIS data are susceptible to encountering imbalanced distributions of different crops in the selected study areas, thereby introducing additional complexity to the task of crop unmixing. Third, real MODIS data are usually subject to more complex noise effects, such as Gaussian noise, impulse noise, speckle noise, etc., rather than just the Gaussian noise considered in the simulated data. Considering the above reasons, we implemented the additional data complexity experiment based on the number of crop classes, while the new simulated dataset involves six crop types for comparative assessment. The results are shown in

Table 5. Mean TAD, AAD, and their variances obtained from multiple independent runs by different methods using the simulated data of four crop categories and six crop categories, respectively, over SNR from 10 to 30 dB. The best results are in bold.

TAD	SNR = 10		SNR = 20		SNR = 30
	4 Classes	6 Classes	4 Classes	6 Classes	4 Classes	6 Classes
NFINDR	0.3096 ± 0.02	0.3339 ± 0.01	0.1849 ± 0.01	0.1365 ± 0.01	0.1758 ± 0.00	0.0821 ± 0.00
PPI	0.3006 ± 0.01	0.3287 ± 0.00	0.1073 ± 0.01	0.2432 ± 0.01	0.0577 ± 0.00	0.2354 ± 0.01
VCA-S	0.1108 ± 0.03	0.1406 ± 0.01	0.0672 ± 0.03	0.0849 ± 0.01	0.0546 ± 0.02	0.0784 ± 0.00
VCA-F	0.0953 ± 0.04	0.1590 ± 0.02	0.0648 ± 0.03	0.0811 ± 0.01	0.0537 ± 0.02	0.0710 ± 0.01
uDAs	0.1531 ± 0.04	0.1597 ± 0.01	0.0612 ± 0.04	0.0501 ± 0.01	0.0554 ± 0.03	0.0474 ± 0.00
PnP	0.0936 ± 0.04	0.1361 ± 0.01	0.0655 ± 0.03	0.0851 ± 0.01	0.0530 ± 0.01	0.0685 ± 0.01
KPmeans	0.1286 ± 0.04	0.1235 ± 0.01	0.0376 ± 0.02	0.0482 ± 0.01	0.0306 ± 0.01	0.0539 ± 0.01
ST-Bdip	0.0816 ± 0.04	0.1053 ± 0.01	0.0354 ± 0.02	0.0436 ± 0.01	0.0107 ± 0.01	0.0419 ± 0.01
AAD	SNR=10		SNR = 20		SNR = 30
	4 classes	6 classes	4 classes	6 classes	4 classes	6 classes
NFINDR	0.5818 ± 0.03	0.5728 ± 0.02	0.5896 ± 0.05	0.3780 ± 0.08	0.6864 ± 0.08	0.2179 ± 0.00
PPI	0.7055 ± 0.03	0.8999 ± 0.01	0.2756 ± 0.02	0.9445 ± 0.03	0.1632 ± 0.03	1.1779 ± 0.13
VCA-S	0.6669 ± 0.18	0.8490 ± 0.05	0.3803 ± 0.19	0.4406 ± 0.04	0.2612 ± 0.18	0.3074 ± 0.04
VCA-F	0.5936 ± 0.13	0.6685 ± 0.04	0.3151 ± 0.10	0.3573 ± 0.07	0.2264 ± 0.07	0.2666 ± 0.03
uDAs	0.5420 ± 0.09	0.5728 ± 0.02	0.2932 ± 0.12	0.2834 ± 0.02	0.2874 ± 0.12	0.2111 ± 0.01
PnP	0.5186 ± 0.13	0.5936 ± 0.05	0.3035 ± 0.10	0.3492 ± 0.05	0.2255 ± 0.05	0.2745 ± 0.05
KPmeans	0.5098 ± 0.09	0.5216 ± 0.02	0.2341 ± 0.13	0.2210 ± 0.02	0.1044 ± 0.04	0.2134 ± 0.13
ST-Bdip	0.3339 ± 0.15	0.3727 ± 0.06	0.1640 ± 0.08	0.2089 ± 0.01	0.0536 ± 0.02	0.1394 ± 0.02

, in which the results indeed demonstrate that temporal unmixing using time-series images with higher complexity leads to diminished accuracies. This is logical with the prevailing trends in existing applications of MODIS crop unmixing, where experiments are often based on the selection of highly distinguishable crops or a small number of crops [6]. For example, Mutlu Ozdogan used the ICA method for the temporal unmixing of winter crops, summer crops, and other classes, achieving satisfactory results [70].

6. Conclusions

In this paper, we proposed a novel DL-based ST-Bdip method for TMA of time-series images, and the EM algorithm was applied to solve the MAP problem in a Bayesian framework. This method consisted of the U-Net model and the EM algorithm, which were used for the iterative optimization of tEMs and abundances. The 3D convolution layers in the model effectively extracted the information both on the spatial and temporal dimensions. In addition, the noise heterogeneity effect in HTI was addressed by the conditional distribution of temporal observations on the M-distance loss. Three simulated data with different noise levels and two MODIS NDVI data were used for the detailed comparison experiments of seven unmixing algorithms with our proposed method, which proved that the proposed method performs better overall. The exploration of temporal unmixing methods is highly beneficial, which utilizes the important temporal information of crops for timely crop monitoring and environmental, economic, and policy decision-making. However, there are still significant challenges due to the complexity of MODIS data and the similarity of crop growth curves. In the future, we will explore how to simulate data with more accurate noise to be closer to the real data. And the advantages of time-series remote sensing data can be used more effectively when applied to a wider range of fields, particularly in the application of data fusion, where higher resolution remote sensing data (e.g., Landsat) can be synergistically utilized.