An Elastic-Window-Based Method for the Underdetermined Problem in Linear Spectral Unmixing to Enhance the Spatial Resolution of the Normalized Difference Vegetation Index Time Series

Abstract

Inverting land cover reflectance or derived indices from low-spatial-resolution images to refine the spatial resolution of this data is cost-effective for land surface monitoring applications that face technical or budget limitations. Based on the linear spectral mixing model, many approaches have successfully unmixed coarse mixed pixels using high-spatial-resolution land cover maps in the past decades. However, in some cases, the solutions of linear systems composed of several mixed pixels may not be acquired due to the underdetermined problem. This study presents the causes of this problem and proposes an iterative approximation strategy to address it. An elastic-window-based algorithm was developed, where the initial size of the window was calculated based on the land cover of the mixed pixel. Mixed pixels of neighborhoods with similar land covers were then selected to form the unmixing linear system, which was examined through a simulation test to ensure it was not underdetermined. Otherwise, the window would expand to include more adjacent pixels. This process was repeated until a successful solution was obtained. A statistical analysis of sixty land cover maps from around the globe shows that the underdetermined problem exists at a low level but becomes more serious with an increase in mixed scale. The results demonstrate that the proposed algorithm effectively prevents the underdetermined problem for mixed pixels of different scales and can be integrated into the coarse NDVI downscaling procedure to refine spatial resolution. This study provides a reference for estimating underdetermined mixed pixels and benefits applications that require dealing with the inversion of land cover values directly.

1. Introduction

Natural hazards, such as floods, disasters, forest fires, and landslides, occur frequently around the world and cause serious damage, and monitoring the environment consistently can reduce these adverse effects [1,2,3]. Vegetation mapping, especially crop mapping, is an important issue for land cover mapping, which is related to carbon storage and production forecast [4,5,6]. With the help of high-temporal-resolution remote sensing images, an effective way can be provided for these land surface dynamic applications extensively [7,8,9,10]. Unfortunately, the temporal resolution of a satellite imaging system is roughly inversely proportional to the spatial resolution, leading to the dilemma that high-temporal resolution remote sensing images have degrading spatial resolution characteristics [11,12,13]. That means it is difficult to provide high-spatial and -temporal remote sensing images in large areas simultaneously. The AVHRR (Advanced Very-High-Resolution Radiometer) has a broad swath width, a short imaging period of less than a day, and about 1000 m of spatial resolution. The MODIS (Moderate Resolution Imaging Spectroradiometer) sensor in Terra and Aqua repeats daily global observation, but its highest image spatial resolution is only 250 m. The images from these high-temporal remote sensors are not suitable for monitoring land surface processes concerning to high-spatial variants [14,15]. Biotic communities and residential buildings usually have a small size in large proportion, resulting in the difficulty of monitoring the damage and assessing the effects caused by natural hazards [16]. The cultivated land field size is often smaller than 250 m in some agricultural regions [17]. Therefore, these high-temporal images seem to lack the ability to map detailed land cover [18]. The Landsat series sensors have a 30 m spatial resolution and 16-day temporal resolution for multispectral bands, which have proven useful for land surface applications [19,20]. However, about 35% of these Landsat images are cloud-cover-contaminated and cloud-free images can be as low as 10% in a given year [21]. According to a statistic, all Landsat images that can be used in the land cover mapping of the year 2010 on a global scale account for 41%. This shows even one available Landsat image can hardly be obtained for about 59% of the land surface in 2010. The current effective coverage and relatively low temporal resolution of Landsat images are limitations for the demands of high-frequency detailed land surface applications in large areas [14,22].

The generation of high-temporal-resolution image series with high-spatial-resolution is one popular solution to bridge the gap [23,24]. It utilizes the spatial detail information in high-spatial resolution images and the phonological information in high-temporal resolution images to create corresponding high-spatial-resolution image series, known as spatial and temporal data fusion [12,25,26,27]. As one of the most widely adopted spatial and temporal data fusion technologies, linear spectral unmixing-based approaches are downscaling the coarse pixel value (e.g., radiance, reflectance, or derived indices) of a low-spatial-resolution image to the inner parts (e.g., endmembers, components, or land cover types) value of the pixel to improve spatial resolution [28,29,30,31,32]. The coarse pixel is regarded as the mixed pixel. Based on the linear spectral mixing model, one mixed pixel is translated into a linear equation, in which the pixel value is denoted as a known parameter and the value of several inner parts is denoted as an unknown parameter. For obvious reasons, several unknown values cannot be solved within only one equation. More pixels must be involved to make up a linear system with equation numbers that are no less than unknown value numbers. Based on the hypothesis of Tobler’s First Law of Geography [33], the pixels around the mixed pixel are regarded as having the same inner parts as the mixed pixel, the same inner part spectrum, and a similar spectral mixing mode. These pixels in adjacent land surface regions contribute to the solving of inner part values. However, when near things are less related to each other on the land surface, the inner parts are usually different around the adjacent pixel regions. This phenomenon often leads to linear systems having no unique solution, called an “underdetermined problem” [34]. The downscaling approach is incapable in this situation because the pixel value cannot be decomposed into an inner part value [22]. Therefore, understanding the relationship between the mixed pixel and its region on the land surface and handling these cases may be the key to avoiding the underdetermined problem of linear spectral unmixing.

In general, the mixed pixel and its adjacent land surface region are considered homogeneous, as Tobler’s First Law of Geography describes [33]. It implies that the inner parts of mixed pixels are included in the pixels of adjacent regions to a certain extent, and the adjacent region affords enough pixels to the linear system of unmixing pixels in inner part value solving. It can be deduced that the acquisition of pixels from a certain-sized adjacent region is an orthodox way to accomplish the unmixing process. On this account, the image scene with a certain extent size was adopted to access the pixels around the adjacent region. For instance, Kerdiles et al. (1995) introduced two sets of 3-by-3-pixel windows to extract the land cover given via TM classification and the AVHRR pixel response [29]. Asner et al. (1997) adopted a sub-scene with 3-by-3 pixels in the AVHRR image to construct the system of equations for each pixel [31]. Zhukov et al. (1999) performed the unmixing of coarse pixels in a moving window mode with a selected size, in which the contextual information of surrounding coarse pixels was used [35]. Zurita-Milla et al. (2008), Amorós-López et al. (2013), and Zhang et al. (2013) applied a sliding window to the classified image to select the coarse neighborhood pixels [32,36,37]. These approaches provide the neighborhood pixels for the unmixing process directly and effectively and have a wide range of applications. With the development of research, the underdetermined problem caused by the unbalance of inner part numbers and neighborhood pixel numbers in the unmixing process draws attention, and two major means are taken to address this problem. On one hand, it is achieved by reducing the number of inner parts in the window. Gevaert et al. (2015) considered the relation between the numbers of the inner part k and the sliding window size n [22]. The inner part with a low proportion in the window is discarded when k is larger than n² to allow a more accurate decomposition. On the other hand, it is achieved by raising the number of neighborhood pixels by enlarging the window size. Wu et al. (2015 and 2016) deemed the window size for each land cover class to be different and developed a correlation-coefficient-based method to obtain the best window size [38,39]. These manners are important for the unmixing process and show great improvement. However, they are not the fundamental solutions.

During the complex distribution of objects on the land surface, various land cover types contribute to each coarse pixel, and the relationship between the inner parts of the pixel and adjacent region pixels exceeds the qualitative definition of homogeneity or heterogeneity. The unmixing process is based on the pixel unit, aiming to solve the inner part values within the local land region. Therefore, setting the optimal window size for each land cover type is still not enough for the whole image pixels to unmix, which limits the performance of recent research facing the underdetermined problem [40]. For each mixed pixel under the unmixing process, the composition of the inner parts is known, but the composition of adjacent pixels is unknown. The distribution and extent of adjacent pixels cannot be estimated based on the inner parts and abundance of central pixels. In other words, the solution to the underdetermined problem in the linear spectral unmixing process might lie in the approximation in the adjacent land region instead of the estimation of the local region. The proposed elastic-window-based method in this paper adjusts its window size dynamically based on the land covers of mixed pixels, ensuring the unmixing linear system is not underdetermined. The rest of the paper is organized as follows. In Section 2, the theoretical basis of inner part inversion using a linear spectral mixing model is reviewed, and the imperfection in approximating the adjacent land region is revealed. In Section 3, an iterative procedure to approximate the local land surface is discussed, and an elastic-window-based algorithm is introduced. In Section 4, a case study is conducted to validate the proposed method based on the GlobeLand30 dataset. Finally, the conclusion and future research are presented.

2. Problem Analysis of Linear Spectral Unmixing

2.1. Inversion of Inner Part Values

Various land surface objects within the remote sensing pixel size are recorded in these low-spatial-resolution pixels. At the pixel scale, these various land surface objects are considered the internal components of the pixel, possessing a finer spatial resolution compared to the pixel itself. In different studies, these objects within the pixel have been referred to as endmembers, components, land cover types, class types, or sub-pixels. In this study, considering the spatial relation between the pixel and its internal components as well as the nature of these objects, we refer to them as the inner parts of the pixel. The main objective of unmixing approaches is to estimate the values of these inner parts based on the pixel value. The mathematical foundation of linear spectral unmixing approaches is the linear spectral mixing model. According to this model, the value of a coarse pixel can be expressed as a linear combination of several inner parts and their respective abundances, as shown in Formula (1), where R represents the pixel value, r denotes the inner part value, f represents the abundance of the inner part, n denotes the total number of inner parts in the pixel, and ε refers to the residual error of the model. Typically, the abundance of the inner parts is derived from a land cover map or land use database by calculating the area proportion of each class type relative to the pixel.

$$R={\displaystyle \sum _{i=1}^n{f}_i\cdot }{r}_i+\epsilon ,$$

(1)

The pixel value and the abundance of the inner parts are considered known parameters, while the values of several inner parts are treated as unknown parameters. Each coarse pixel is represented using a linear equation. If more pixels are available, a linear system can be constructed by combining these linear equations, as shown in Formula (2). Taking into account that the residuals remain constant within each equation, the coefficient matrix of this linear system can be expressed as Formula (3). According to the computational principles of linear systems, the number of equations m should be equal to or greater than the number of unknown parameters n. The optimal solution of the linear system can be obtained using the least squares method under this condition [41]. Otherwise, the linear system would have no unique solution, which is referred to as an “underdetermined problem”. The solvability of a linear system depends largely on its coefficient matrix. Therefore, the crucial aspect of inverting the values of inner parts for mixed pixels in linear spectral unmixing lies in creating linear systems with appropriate coefficient matrices.

$$\begin{array}{c}{f}_{1,1}\cdot r_1+f_{1,2}\cdot r_2+\cdots +f_{1,\mathsf{n}}\cdot r_n+{\epsilon }_1=R_1\\ f_{2,1}\cdot r_1+f_{2,2}\cdot r_2+\cdots +f_{2,\mathsf{n}}\cdot r_n+{\epsilon }_2=R_2\\ ⋮\\ f_{\mathsf{m},1}\cdot r_1+f_{\mathsf{m},2}\cdot r_2+\cdots +f_{\mathsf{m},\mathsf{n}}\cdot r_n+{\epsilon }_m=R_m\end{array},$$

(2)

$$F=\left[\begin{array}{cccc}{f}_{1,1}& f_{1,2}& \cdots & f_{1,\mathsf{n}}\\ f_{2,1}& f_{2,2}& \cdots & f_{2,\mathsf{n}}\\ ⋮& ⋮& ⋮& ⋮\\ f_{\mathsf{m},1}& f_{\mathsf{m},2}& \cdots & f_{\mathsf{m},\mathsf{n}}\end{array}\right],$$

(3)

2.2. Imperfection in Approximating Adjacent Land Region

The values of coarse pixels in different bands are distinct from one another. Therefore, the values of inner parts in different bands cannot be assumed to be the same unknown parameters in a linear system. The inversion of inner part values for a pixel needs to be performed within the current band. Considering the variations caused by natural or artificial influences, the same land cover type may exhibit inconsistent characteristics in different regions. For instance, crop growth can be affected by factors such as agrotype, surface water, and cultivation, resulting in spatially variable behavior. This phenomenon, known as spectral-spatial variability [42], can cause the inner part values of pixels in different locations to differ significantly from each other across the entire image band. Based on the assumption that “Everything is related to everything else, but near things are more related than distant things”, known as Tobler’s First Law of Geography [33], coarse pixels in adjacent regions are logically selected to form the linear system for the central pixel. Typically, the adjacent pixels are chosen using a regular window of a specific size. The window is defined using Formula (4), where S represents the odd length or width of the rectangular area and N represents the number of inner parts on the land surface.

$$S^2\ge N,$$

(4)

In most cases, this rectangular area is square-shaped, with a coarse pixel size as the unit—for example, a window region of 3 × 3 or 5 × 5 pixels centered on a coarse pixel. By using this formula, the window size S can be calculated based on the number of inner parts N in the mixed pixel. Then, the adjacent pixels within this window are selected to form the linear system, with pixel values and inner part abundances as known parameters and the inner part values as unknown parameters.

An example is provided to illustrate the problem by combining a land cover map and a low-spatial-resolution image band. In Figure 1, the grid represents the pixels in the low-spatial-resolution image band, and specific pixels are labeled as A, B, C, and D. Pixel A encompasses six land cover types: crop, built-up area, bare land, grass, forest, and shrub. Pixel B includes five land cover types: crop, forest, grass, water, and wetland. Pixel C consists of seven land cover types: crop, built-up area, bare land, grass, forest, shrub, and wetland. Finally, pixel D contains two land cover types: crop and grass.

According to Formula (4), the window size is 3 pixels by 3 pixels. A total of nine pixels inside the window are selected to form the linear system for the inversion of inner part values for each of these four pixels separately. However, the linear systems of pixels A, B, and C encounter the underdetermined problem.

While bare land and shrubs are present in pixel A, they do not exist in the pixels within the window. The same situation occurs for pixel C, where built-up areas and shrubs cannot be separated in the unmixing process. Both cases, A and C, involve situations where undividable inner parts are present within the window. The pixels around pixel B can handle this unmixing scenario, but due to a limited number of adjacent pixels (e.g., pixel c9) included in the linear system of pixel B, the underdetermined problem arises (although a weighted method can help improve this case). The situation with pixel B involves undividable inner parts in the adjacent region within the specific window. When it comes to selecting enough pixels to decompose a given pixel, the appropriate window sizes for pixels A, B, and C are 5, 3, and 7, respectively. The inner part numbers of these pixels satisfy Formula (4), but the window sizes adopted are ill-suited. Thus, it appears that opting for a larger window size might be a suitable choice. However, the case of pixel D indicates that this assumption is incorrect. Pixel D encompasses grass and crop, which can be easily decomposed using a 3 × 3 window. Unfortunately, employing a window size of 5 includes pixel C, which is challenging to decompose due to its adjacent region with pixel D. The abundances corresponding to these 3 × 3 windows for pixels A, B, C, and D are listed in

Table 1. The inner part types and abundance of the coarse pixels (1-1, 1-2, 1-3, and 1-4 show A, B, C, and D in Figure 1, respectively, unit: %).

1-1
Num.	Pixel Code	Crop	Forest	Grass	Shrub	Wetland	Water	Built-Up	Bare Land
1	A	42	17	3	9	0	0	7	22
2	b2	78	2	20	0	0	0	0	0
3	b3	82	0	0	0	0	0	18	0
4	b4	88	4	0	0	0	0	8	0
5	c2	83	0	0	0	0	0	0	0
6	c4	1	61	5	0	28	3	1	0
7	c2	90	0	0	0	0	0	10	0
8	c3	85	15	0	0	0	0	0	0
9	c4	38	13	1	0	8	40	0	0
1-2
Num.	Pixel Code	Crop	Forest	Grass	Shrub	Wetland	Water	Built-Up	Bare Land
1	B	5	12	69	2	4	8	0	0
2	b7	35	10	26	0	9	20	0	0
3	b8	32	25	31	0	0	12	0	0
4	b9	69	12	12	0	0	8	0	0
5	c7	3	14	71	0	5	7	0	0
6	c9	50	19	15	0	0	0	4	13
7	d7	18	4	78	0	0	0	0	0
8	d8	75	0	25	0	0	0	0	0
9	d9	87	5	4	0	2	2	0	0
1-3
Num.	Pixel Code	Crop	Forest	Grass	Shrub	Wetland	Water	Built-Up	Bare Land
1	C	26	2	0	13	40	1	10	8
2	d5	0	0	0	0	36	64	0	0
3	d6	5	0	10	0	70	15	0	0
4	d7	19	4	77	0	0	0	0	0
5	e5	17	0	0	0	44	39	0	0
6	e7	43	2	56	0	0	0	0	0
7	f5	97	0	2	0	1	0	0	0
8	f6	95	0	5	0	0	0	0	0
9	f7	64	1	35	0	0	0	0	0
1-4
Num.	Pixel Code	Crop	Forest	Grass	Shrub	Wetland	Water	Built-Up	Bare Land
1	D	81	0	19	0	0	0	0	0
2	e7	43	2	56	0	0	0	0	0
3	e8	94	1	5	0	0	0	0	0
4	e9	80	0	2	0	5	1	12	0
5	f7	64	1	35	0	0	0	0	0
6	f9	85	0	8	0	0	0	7	0
7	g7	93	0	7	0	0	0	0	0
8	g8	77	0	22	0	0	0	1	0
9	g9	93	0	3	0	0	0	5	0

.

Based on the land cover map, the value of abundance in each coarse pixel can be derived by calculating the proportion of area occupied by every land cover type to the total coarse pixel area. The abundance value of adjacent pixels located in windows of pixels A, B, C, and D has been shown in Table 1, each row corresponding to the distribution of land convert types and their abundance within a pixel. These arrangements of land cover types and their abundance can be seen as the coefficient matrix of the linear systems for pixels A, B, C, and D. By combining the principle of the existence of solutions to linear systems, the solutions of pixels A (Table 1(1-1)), B (Table 1(1-2)), and C (Table 1(1-3)) are unavailable.

3. Solution

3.1. Iterative Approximation Procedure

The diversity of land surface types and their various distributions result in significant complexity on the land surface. This complexity extends from large to local areas and is reflected in remote sensing images with a discrete pattern. As a result, many pixels cover multiple land cover types, and it is difficult to ensure that the same land cover types are present in other pixels within the adjacent region. Expanding the extent around the adjacent region is important and effective for selecting pixels to form the linear system with the same land cover types, but it is not essential. The presence of the required pixels to such an extent is uncertain. Additionally, the size of the extent calculated based on the number of land cover types in the unmixing pixel is only an estimated value and is not always reliable. A larger extent can introduce more pixels compared to a smaller extent, which can help reduce the possibility of the underdetermined problem to some extent. However, as shown in the case of pixel D in Figure 1, it is not a fundamental solution. Moreover, since each linear system determines a least squares solution by involving the same land cover type but different values to a large extent, the outputs can reduce the spectral variability. Therefore, more attention should be paid to the complexity of the adjacent region.

Theoretically, each pixel has other pixels that can assist in the successful inversion of inner part values. These pixels are scattered in the adjacent region around the mixed pixel at various distances, implying a spatial distribution pattern of scattered pixels. Unmixing a pixel is a process that approximates this spatial distribution pattern using certain methods. However, due to the unknowns of these spatial distribution patterns, it is almost impossible to select all the pixels at once within a certain extent or specific shapes. Nevertheless, continuously approximating the spatial distribution pattern in the adjacent region is a preferable solution until knowledge of these patterns becomes available. Because of this idea, a strategy, as shown in Figure 2, is introduced to guide the approximation procedure. It involves utilizing the mixed pixel to estimate the initial spatial approximating mode and selecting pixels from the adjacent pixels within the estimation region. Based on the adjacent pixels, an optimization option can be chosen to adjust the initial approximating mode. The adjustment may involve assigning higher weights to similar or nearby pixels or directly excluding opposite pixels. Through iterative selection and optimization, the spatial approximating mode approaches the spatial distribution pattern of theoretical pixels. At the same time, it is ensured that the number of selected pixels is not less than the inner part numbers N of the mixed pixel, and the rank of the coefficient matrix F, which is composed of the abundances of these pixels, should be equal to N (as shown in Formula (5)). The strategy estimates the spatial distribution pattern of pixels based on the mixed pixel according to Tobler’s First Law of Geography. Furthermore, the composition of land cover types and their abundance among adjacent pixels is utilized to optimize the previous spatial approximation mode to prevent failure when the assumption is not valid. The iterative procedure modifies the spatial approximation mode continuously and ensures the approximation of the spatial distribution pattern properly in the adjacent region. In extreme cases, an excessively large range may be traversed for the approximation process.

$$rankF=N$$

(5)

3.2. The Elastic-Window-Based Algorithm

As an implementation of the iterative approximating strategy, an elastic-window-based algorithm was developed to approach the adjacent pixels around the coarse pixel and address the underdetermined problem of linear spectral unmixing. The entire workflow of the algorithm is shown in Figure 3.

In the algorithm, the inner part types and numbers within the mixed pixel were extracted from the fine spatial resolution land cover map. The initial window size was estimated using Formula (4). Each mixed pixel created its own initial window centered around it to select the adjacent pixels. Then, pixels containing different land cover types compared to the central mixed pixel were discarded. The abundances of all inner parts among these selected pixels were extracted, and a corresponding linear system was constructed. A simulated examination was also conducted to assess the presence of the underdetermined problem.

In the simulated examination, all inner part values were assumed to be known by assigning constant values. The values of the pixels could be calculated using the linear spectral mixing model with abundances and zero residual error. The linear system, treating the inner part values as unknown parameters, was then inverted using the ordinary least squares method. If the inverted values matched the constant values, the selected pixels could be used to estimate the inner part values of the mixed pixel. If not, the window size would be expanded by two-pixel units, and the additional pixels within the expanded window would be included in the new linear system. The simulated examination would be repeated until the inverted values matched the constant values. This simulated examination serves as a practical alternative to Formula (5).

The matrix form of the algorithm’s linear system is given in Formula (6), where F represents the coefficient matrix, r represents the matrix of inner part values, ε represents the residual error matrix of the linear system, and R_E represents the matrix of coarse pixel values. The limiting conditions for R_E are shown in Formula (7), where R_i represents the value of a coarse pixel located within the spatial range of the elastic window S_W, and k represents a coarse pixel that contains different inner part types compared to the central pixel.

$$F\cdot r+\epsilon =R_E,$$

(6)

$$R_E=\left\{R_{\mathrm{i}}|R_{\mathrm{i}}\in S_W,i\ne k\right\},$$

(7)

The least-squares solution of Formula (6) is given in Formula (8). The superscript F^T is the transformation of the matrix, and (F^T·F)⁻¹ is the inverse of the matrix.

$$r={(F^T\cdot F)}^{-1}\cdot F^T{R}_E$$

(8)

It should be noted that achieving an exact equivalency between the inverted values and constant values during the simulated examination may be unlikely due to the inherent computational errors of computers. Therefore, considering a minor difference between the values is a more practical approach. To prevent the situation where more than two inner parts of the mixed pixel never appear in the entire image, leading to the process getting stuck in a loop, it is necessary to set a maximum number of iterations for the loop. This ensures that the algorithm does not continue indefinitely and allows for a predetermined number of iterations before reaching a stopping condition.

4. Experiment

4.1. Study Area and Data

To investigate the theoretical existence of the underdetermined problem in linear spectral unmixing, a simulation experiment was conducted within a local neighborhood of coarse pixels. The local region had a size of 3 × 3 coarse pixels, and it was assumed that various land cover types appeared randomly around the center coarse pixel. The land cover layout of the entire region was ensured to be non-repetitive throughout the experiment. The underdetermined results were recorded for local regions with land cover types ranging from 2 to 6. Several issues of land cover types were taken into account, including the huge calculation workload among the global simulation, the relationship between the spatial resolution of coarse pixels and the spatial scale of land cover types within the corresponding coarse pixel extent, and the practical application of linear spectral unmixing.

Following the simulation experiment, a statistical experiment was conducted on a global scale using the land cover data from GlobeLand30 (www.ngcc.cn (accessed on 18 June 2023)). The land cover types considered in this dataset included water, wetland, artificial surface, cultivated land, forest, shrubland, grassland, bare land, tundra, permanent snow, and ice. The world terrestrial ecoregion map (www.worldwildlife.org (accessed on 18 June 2023)) was used as a reference to manually select map sheets from GlobeLand30. The sampling process aimed to achieve qualitative balance in terms of continent and ecological region factors, ensuring representative map sheets with both heterogeneous and homogeneous landscapes from around the globe. Most map sheets had a size of approximately 6 degrees in longitude and 5 degrees in latitude, while a few map sheets in high latitudes spanned 12 degrees in longitude. A total of sixty map sheets were used in the experiment, as shown in Figure 4 and

Table 2. The UTM zone and row of the sample map sheets of GlobeLand30 in 2010.

Num.	Name	Num.	Name	Num.	Name	Num.	Name	Num.	Name	Num.	Name
1	N07_65	11	S23_15	21	N33_40	31	N38_05	41	N45_30	51	N47_60
2	N11_60	12	S20_00	22	N32_50	32	N32_10	42	N43_40	52	N51_70
3	N15_55	13	N18_05	23	N35_50	33	N29_15	43	N47_40	53	N50_25
4	N11_45	14	N12_30	24	N36_40	34	N35_25	44	N47_50	54	N50_40
5	N14_40	15	N17_45	25	S34_30	35	N30_30	45	N45_55	55	N49_50
6	N13_20	16	N38_20	26	S35_20	36	S51_00	46	S50_30	56	N51_55
7	N16_30	17	N42_25	27	S38_15	37	N49_00	47	S55_25	57	N57_55
8	S18_45	18	N40_35	28	S34_10	38	N48_05	48	S52_20	58	S59_40
9	S21_30	19	N41_50	29	S33_00	39	N48_15	49	S54_15	59	N43_10
10	S19_15	20	N37_65	30	N35_05	40	N47_25	50	S54_00	60	N54_35

. For each map sheet, ten different mixed pixel scales were chosen to simulate varying mixed situations within the same landscape region. The inner part values were then inverted based on these mixed pixels and their abundances, following the methodology described in the simulation test of Section 2.2. The spatial resolution ratio between the mixed pixel and the land cover pixel ranged from 5 to 50, with an increment of 5 (referred to as the mixed scale).

A land cover map tile from the N11_60 region, consisting of 4000 × 4000 pixels, was selected to analyze the performance of the proposed elastic-window-based algorithm for mixed scales ranging from 5 to 50. As a point of comparison, regular windows with fixed sizes of 3, 7, 11, 15, 21, 27, 33, and 45 were also selected. Additionally, an experiment was conducted to downscale MODIS-NDVI (Normalized Difference Vegetation Index) data in the Shandong Province of China using the proposed algorithm. Fourteen Landsat images from the year 2010 (shown in

Table 3. The Landsat-8 OLI multi-spectral images covering the Shandong Province.

Num.	WRS Path and Row	Day of Year	Num.	WRS Path and Row	Day of Year
1	119034	2010256	8	121036	2010014
2	119035	2010256	9	122034	2010117
3	120034	2010312	10	122035	2010069
4	120035	2010123	11	122036	2010117
5	120036	2010231	12	123034	2010172
6	121034	2010254	13	123035	2010172
7	121035	2010014	14	123036	2010172

), covering the local region, were selected to generate a fine-resolution land cover map with a spatial resolution of 30 m. The ISO-DATA classification method in ENVI was employed for this purpose. The 16-day composite of MODIS-NDVI data (MOD13Q1 h27v05) for the Shandong Province in 2010 was collected to obtain the coarse pixel values. A total of 23 NDVI products were resampled to a resolution of 240 m to align with the resolution of the land cover map.

4.2. Results and Discussion

4.2.1. Local Simulation Results of the Coarse Pixel Neighborhood

For each coarse pixel, the number of class types for its inner parts could be 1, 2, 3, and so on. When the number of class types is greater than 2, the possible combinations of these class types increase significantly.

Table 4. Simulated results within an adjacent region of 3 × 3 pixels.

Class Types	Combinations of Class Types	Layouts of Combinations	Underdetermined Layouts	Underdetermined Proportion (‰)
2	3	45	0	0
3	7	3003	3	9.99
4	15	319,770	342	10.70
5	31	48,903,492	44,528	9.11
6	63	9,440,350,920	8,098,340	8.58

illustrates the layouts of total combinations for each coarse pixel within a 3 × 3 pixel region, and it can be observed that the number of layout combinations increases substantially with the growth of class types. A coarse pixel covered by two class types would never encounter the underdetermined problem within a 3 × 3 pixel adjacent region. However, when three or more class types are present in a coarse pixel, the likelihood of encountering underdetermined layouts becomes more diverse. Consequently, the number of underdetermined layouts increases rapidly. Despite this increase, the proportion of underdetermined problems remains stable, typically at a level of one thousandth.

The results presented in Table 4 demonstrate a theoretical random state, and it is worth noting that the size of the 3 × 3 region may not be large enough. However, these results still clearly indicate the presence of the underdetermined problem and its trend in linear spectral unmixing. It provides statistical averages regarding the likelihood of encountering an underdetermined problem within a certain size region. Even though the adjacent region is a small area, it still encompasses a wide variety of objects present on the land surface and their complex distribution. The local region includes diverse class types and distributions of these objects. Consequently, the linear system composed of adjacent pixels within the region cannot guarantee unique solutions. This implies that the underdetermined problem is inherently associated with the adjacent pixels to some extent (Figure 5).

4.2.2. Statistics for Underdetermined Pixels in Globe Scale

An unmixing test was conducted on all sixty land cover maps at ten different mixed scales, as described in Section 4.1. The test aimed to count the number of mixed pixels that encountered the underdetermined problem and calculate the proportion of these pixels among all mixed pixels. Figure 6 illustrates the proportions of underdetermined pixels for the ten mixed scales. The results indicate that the underdetermined problem is present in all land surface areas of the study when performing mixed-pixel unmixing based on the linear spectral mixing model. Additionally, Figure 6 shows that the underdetermined problem exists to varying degrees across the ten mixed scales.

For the mixed scale of 5, half of the land cover tiles have more than 21 underdetermined pixels in every 10,000 mixed pixels. As for the mixed scales ranging from 15 to 50, these numbers are 51, 92, 99, 119, 139, 170, 176, 193, and 194 underdetermined pixels in every ten thousand mixed pixels. Considering the spatial resolution of the land cover map is 30 m, this indicates that approximately 4.67 acres per 100-by-100 mixed pixels area cannot be inverted based on the linear spectral unmixing theory when the image pixel spatial resolution is 150 m. As the pixel spatial scale increases, more and more pixels become mixed pixels, and larger portions of the land surface become unsolvable. These unsolvable areas amount to 11.34, 20.45, 22.01, 26.46, 30.91, 37.80, 39.14, 42.92, and 43.14 acres when the image spatial resolutions are 300, 450, 600, 750, 900, 1050, 1200, 1350, and 1500 m, respectively, for half of the land cover maps over the study area. These results can serve as a reference for estimating unsolvable areas when unmixing land surfaces with mixed pixels using remote sensors with low spatial resolutions.

Figure 7 displays the average proportions of underdetermined pixels and mixed pixels for the ten mixed scales across all sixty land cover maps. The graph shows that the number of mixed pixels increases with the mixed scales, and a larger portion of these mixed pixels gradually encounters the underdetermined problem. Generally, the proportions of underdetermined pixels remain at a low level, accounting for 0.1% to 2% of mixed pixels. However, when considering the enormous number of mixed pixels from various remote sensing systems that image the global land surface day and night, this issue becomes significant.

4.2.3. Performance of Elastic-Window-Based Algorithm

Figure 8 shows the underdetermined pixel proportions of the test between the regular window and elastic window approaches for the land cover map with a size of 4000 pixels by 4000 pixels from tile N11_60. The underdetermined pixel proportions of the regular window approach increase with the mixed scale. On the other hand, the results of the elastic window approach remain at a low level for all ten mixed scales and are almost unaffected by the mixed scale size. For a mixed scale of 5, only 0.003% of mixed pixels cannot be inverted using the elastic window algorithm. This proportion decreases to 0.001% for mixed scale 10 and becomes zero for mixed scales ranging from 15 to 50. This indicates that the elastic-window-based algorithm is effective and robust for all ten mixed scales when inverting the inner part values using the linear spectral mixing model.

In contrast, when inverting the inner part values based on regular windows with fixed sizes, the underdetermined pixel proportions vary with the window size. The lowest proportion is 0.223% with a window size of 3, while the highest proportion is 14.799% with a window size of 45. On each mixed scale, the number of underdetermined pixels increases with the expansion of the window size.

Traditionally, it is believed that a larger window size provides more adjacent pixels to construct the unmixing linear system, leading to improved accuracy with the help of the least squares method. However, this result indicates that a larger regular window does not help reduce the number of mixed pixels encountering the underdetermined problem. A larger regular window involves more mixed pixels that consist of land covers different from the unmixing mixed pixel. These different land covers introduce additional unknowns to the unmixing linear system, making the inversion more challenging. With the increase in the mixed scale, more and more mixed pixels are likely to encounter the underdetermined problem. Therefore, compared to a regular window with a fixed size, the elastic-window-based algorithm, which approximates the land surface iteratively, is preferable for preventing the underdetermined problem.

Figure 9 shows a comparison of time consumption between the elastic window and the fixed window (size increasing from 3 to 45) on 10 mixed scales separately. It shows that the unmixing process based on a regular window with a size of 3 is the fastest in terms of time consumption for all mixed scales. With an increase in the mixed scale, the number of mixed pixels in the study area decreases, leading to a reduction in the total unmixing time. The reason why the decrease in time consumption appears with the increase in mixed scales is the reduction in the total number of mixed pixels in a fixed experimental region. Therefore, the time consumption decreases sharply in exponent form for either a fixed window or an elastic window in each mixed scale. It can be seen that the elastic window takes a relatively small amount of time compared with other fixed windows in statistics. Until the mixed scale reaches 45, the time consumption of the elastic window stands at a low level, despite not being the lowest value. At the mixed scale 50, the elastic window spends extra time on iterative approximation of the local region compared with the least time-consuming fixed window. The elastic-window-based algorithm demonstrates high performance, indicating that the iterative approximation of the land surface does not require significant time, making it a practical approach.

4.2.4. MODIS-NDVI Downscaling Based on the Elastic Window

The experiment involved downscaling 23 coarse spatial resolution MOD13Q1 NDVI images to a finer resolution within the study area of Province Shandong in the year 2010 using the elastic window method described in Section 3.2. The land cover class maps were obtained through the ISO-DATA classification of Landsat images, as shown in Table 3. The results of the downscaling experiment are presented in Figure 10, where the images are listed in row-major order based on 16-day intervals.

Due to the large size of the study region, images with a spatial resolution of 30 m occupy a significant amount of computer memory. The Dell workstation (i5, 8G) with a Windows-based operating system was used for the experiment. The source code of the method was compiled with ENVI 5.4 desktop software (https://envi.geoscene.cn/ accessed on 18 June 2023). To overcome this limitation, the high-resolution NDVI data images were resampled and segmented into small domains, resulting in a striping phenomenon in the mosaic image. Nevertheless, the high-resolution NDVI data images exhibit clear details and provide spatial information. Isolated outlier pixels are almost non-existent in Figure 10, indicating that the elastic-window-based downscaling method can effectively handle various land surface conditions over a wide area and mitigate the underdetermined problem to a certain extent. A comparison of the 23 high-resolution NDVI images throughout the year reveals noticeable changes in the land surface during vegetation growth periods.

Figure 11 presents the high-resolution NDVI time series curves for several land cover types within the study area of Shandong Province. The dashed line in the figure represents the NDVI curve of MOD13Q1 products at the same geographic location. It can be observed that the range of NDVI values for typical land cover types is greater in high-resolution pixels compared to low-resolution pixels at the same scale. The fluctuation trends of the generated high-resolution NDVI curves throughout the year align with the original MOD13Q1 NDVI curves. Notably, for the built-up land cover type, significant differences in NDVI values can be observed at multiple time points between the original NDVI values and the downscaling values. However, the overall fluctuation trend remains consistent. For cropland, the generated NDVI curve exhibits a similar trend to the original curve, with higher NDVI values in the fine-resolution data. The NDVI values show a sudden decrease and subsequent increase in the 14th period, indicating that the generated NDVI curve captures the cultivation characteristics of two crops per year in the study area, similar to the original NDVI curve.

On forest land, the trend of the NDVI curve changes after construction is consistent with the original curve, and the values are generally higher. Similarly, on grassland, the constructed NDVI curve shows good consistency with the original NDVI curve in terms of trends and values. However, for wetlands, there is a noticeable difference between the generated NDVI curve and the original curve. The fluctuation amplitude of the generated curve is smaller and smoother, with values generally fluctuating around zero. In the case of water, the generated curve is generally lower than the original curve, consistently below zero, and exhibits relatively small variation compared to the original curve.

To evaluate the accuracy and reliability of the generated fine NDVI compared to the original MODIS-NDVI, an accuracy assessment was conducted. A stratified sampling strategy was employed, where 2 out of 14 Landsat scenes were selected and points were randomly scattered to obtain values of fine curves and original curves for different land cover types, as shown in

Table 5. The numbers of sampling points of two areas and the according land cover types.

Year	Land Cover	Sampling Point Number of 119034	Sampling Point Number of 122035
2010	Crop	455	622
2010	Forest	417	228
2010	Grass	193	252
2010	Water	348	240
2010	Built-up	565	418

. In total, there were 3738 sampling points in 2010.

The original values of the MODIS-NDVI curve were considered the true values. Three measures, namely variance, mean difference, and correlation coefficient, were used to assess the similarity between the constructed curve and the original curve. The variance calculation is shown in Formula (9), where n represents the total number of periods in the NDVI time series, and ΔNDVI represents the difference between the two curves during the corresponding period. The parameter Δμ represents the average value of the NDVI difference between the two curves, as shown in Formula (10). The correlation coefficient calculation is shown in Formula (11), where n represents the total number of periods in the temporal NDVI, a represents the generated NDVI curve, and b represents the low-spatial-resolution curve. The parameters μ_a and μ_b represent the NDVI means of the fine curve and the low-spatial-resolution curve, respectively.

$$\sigma =\sqrt{\frac{{\displaystyle \sum _{i=1}^n(\Delta \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_i-\Delta \mu )}}{n-1}}$$

(9)

$$\Delta \mu =\frac{1}{n}{\displaystyle \sum _{i=1}^n(\Delta \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_i)}$$

(10)

$$c=\frac{{\displaystyle \sum _{i=1}^n(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{(a,i)}-{\mu }_a)(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{(b,i)}-{\mu }_b)}}{\sqrt{{\displaystyle \sum _{i=1}^n{(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{(a,i)}-{\mu }_a)}^2}}\sqrt{{\displaystyle \sum _{i=1}^n{(\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{(b,i)}-{\mu }_b)}^2}}}$$

(11)

The statistical distribution of the constructed NDVI and the original NDVI for the main land cover categories at the sampling points in 2010 is presented in Figure 12. It can be observed that the mean difference of NDVI values (Figure 12A) in the sampled land cover types follows a normal distribution. Generally, the constructed NDVI curve has higher values than the original NDVI curve, with an average increase in approximately 0.16 points across all points. Except for forest land types, the standard deviation (Figure 12B) for all other types is around 0.1, and the correlation coefficients (Figure 12C) for each sampling point are generally above 0.9. Based on the local climate characteristics of the study area, it can be inferred that vegetation growth was stagnant during the specified period. Consequently, this image may not comprehensively reflect the accurate vegetation distribution in the area. When this image is used as the basis for generating classification maps, it may introduce errors and hurt the classification results of ground objects. Therefore, some errors were introduced into the generated high-spatial-resolution NDVI dataset in this specific scenic area.

5. Conclusions

In this study, we have described the underdetermined problem that arises when inverting the inner part values using a linear spectral mixing model. Through simulation studies and a case study with a class map, we have provided a detailed illustration of this problem. We have concluded that the main reason for encountering the most underdetermined cases is the imperfect approximation of land surface using current regular window-based methods. To address this issue, we have proposed an iterative approximation strategy that focuses on local mixed pixels and their adjacent regions. This strategy aims to provide a solution to the underdetermined problem by gradually approaching the solution instead of attempting to solve it in a single step when dealing with mixed pixels and their neighborhoods. As an example of this strategy, we have developed an elastic-window-based algorithm that can handle mixed pixels that may be involved in linear systems to prevent these systems from encountering the underdetermined problem.

Based on a statistical analysis of sixty land cover maps from around the world, the results have shown that the underdetermined problem exists at a relatively low level (ranging from dozens to hundreds in every ten thousand mixed pixels) in all of these areas. However, the severity of the problem tends to increase to some extent with the increase in mixed scales. A comparison between the single-window-based method and the proposed algorithm has demonstrated that the proposed algorithm can reduce the proportions of mixed pixels that encounter the underdetermined problem. However, it should be noted that the proposed algorithm requires more time compared to the single window-based method, mainly due to the iterative approximation of local mixed pixels. Ensuring that equations do not result in an underdetermined problem is crucial for successfully inverting linear systems.

The presented strategy and algorithm can have potential applications in various areas that require the inversion of inner part values directly—for example, downscaling high-temporal-resolution images or time series of derived indices (such as NDVI) to obtain high-spatial-resolution time sequences based on linear spectral unmixing. Additionally, it can provide an alternative approach for parameter inversions related to local land surface approximation, with mixed pixels as the computational unit. However, it is important to note that although the experiment has been tested with ten mixed scales based on 30 m spatial resolution land cover maps, further studies are needed considering various land cover maps with different spatial resolutions and classification systems. Future research should focus on improving the accuracy of unmixing in these scenarios.