1.

Introduction

Moderate Resolution Imaging Spectroradiometer (MODIS), as well as later-developed hyperspectral sensors have made great breakthroughs in spectral channel settings compared with earlier remote sensors. There are 36 discrete channels, including 20 reflective spectral channels, in a MODIS image, and each pixel of the image acquires many bands of light intensity data from the spectrum, instead of just the three bands of the RGB color model, which makes it possible to accurately depict the spectrum characteristics of typical ground features using not only the wavelengths, ranges, and intensities of the peaks and valleys but also the integral area that is in the range enclosed by the spectral reflectance curves of the ground features and the $x$-axis (in Cartesian coordinates). The MODIS visits the globe once or twice per day with coarse resolution of 250 to 1000 m. However, the spatial resolution of MODIS images is not high enough to clearly distinguish different ground features. In many cases, a MODIS pixel is a mixed pixel that is covered by multiple land cover types, which has a significant influence on the information that can be derived.^1,2 Thus, the decomposition of mixed pixels in MODIS images is critically important for the application of MODIS data in many fields, such as mapping land cover distributions,³ evaluating vegetation/soil fractional cover,^4–6 monitoring and evaluating karst rocky desertification,⁷ flood mapping,^8,9 and retrieving fire temperature and area.¹⁰

The spectral characteristics of ground features are the basis not only for identifying them in remote sensing images but also for decomposing mixed pixels in images. The decomposition of mixed pixels is generally based on a linear spectral mixture model (LSMM) or a nonlinear spectral mixture model (NLSMM).¹¹ Although the NLSMM is more applicable when the multiple scattering among distinct endmembers is not negligible,¹² such as in intimate mineral mixtures and vegetation canopies,¹³ the LSMM is a mature and more widely used technique than the NLSMM.^14,15

To apply existing methods for decomposing mixed pixels, the endmembers must be obtained. Endmember extraction is the process of selecting a collection of pure signature spectra of ground features present in a remote sensing image.^16–18 The corresponding abundance of each endmember is usually estimated by using the fully constrained least squares (FCLS) method based on the LSMM.¹⁹ The endmember extraction is generally performed in two ways: (1) by deriving them directly from the remote sensing images, which is referred to as image endmember analysis;¹ or (2) from a spectral library that contains the spectra of known target features measured in the field or laboratory, which is referred to as library endmember analysis.²⁰ When considering the effect factors, such as the atmospheric interaction and remote sensor peculiarities and noise, image endmember analysis is now most widely used. Two major approaches are used to extract endmembers based on the LSMM. One approach uses geometrical methods, including the pixel purity index (PPI),²¹ the N-finder algorithm (N-FINDR),²² the sequential maximum angle convex cone (SMACC),²³ vertex component analysis (VCA),²⁴ etc., of which the PPI and SMACC methods are widely used for decomposing mixed pixels in remote sensing images due to their publicity and availability in the Environment for Visualizing Images software.²⁵ Another approach uses statistical methods, such as independent component analysis.²⁶

It is usually difficult to acquire pure pixels in a MODIS image because of its spatial resolution limit. Many researchers have suggested that there are no pure pixels in remote sensing images with low spatial resolution.^17,27,28 Some authors have tried to use nonnegative matrix factorization (NMF) for hyperspectral data unmixing.^29,30 Miao and Qi³¹ presented a constrained NMF (MVC-NMF) method without the pure-pixel assumption for unsupervised endmember extraction from highly mixed image data.

The accuracy of extracted endmembers has a great impact on the unmixing accuracy. To assure unmixing accuracy, an unmixing method for MODIS data that does not resort to extracting endmember spectra is taken into account.

Adjacent channels in multispectral/hyperspectral imagery have good correlation and often contain similar information, which produces redundancies in a multispectral/hyperspectral dataset.^32,33 Thus, many conventional unmixing methods, e.g., PPI,²¹ manual endmember selection tool,³² N-FINDR,²² spectral mixture analysis based on simulated annealing,³⁴ VCA,²⁴ simplex growing algorithm,³⁵ Gaussian elimination method,³⁶ etc., use statistical techniques such as principal component analysis (PCA) to reduce the dimensionality of the image data for both computational time saving and signal-to-noise improvement. Then, a set of uncorrelated variables (principal components) are generated, and those containing the most information from the original bands are selected to extract endmember spectra. Each endmember spectrum can be constructed as a linear combination of the principal components.³² As a statistical technique, the PCA transformation is highly dependent on the numerical characteristics of the image. Hence, the principal components vary with the images, and the difficulty of interpreting a priori the content of the principal components is an inherent problem of PCA.^33,37

A set of basis functions are independent of each other as well as principal components, and they are purely theoretical functions. In mathematics, a complex curve can be represented as a linear combination of a set of basis functions.^38,39 Similarly, the spectral curve made by mixing spectra with more than one ground cover type can also be represented as a linear combination of a set of basis functions. The basis functions can be employed to reduce the dimensionality of the image data and characterize the spectral curve of each pixel without redundant information. A comparison of basis functions with the principal components generated by using PCA shows that on one hand, the basis functions can be used to depict each endmember spectrum with a linear combination as well as the principal components do. On the other hand, there exists the difference that the basis functions are invariant and independent of image data. Thus, the coefficients of the basis functions for pixels in different images are comparable, and the coefficients can be employed to depict the spectral curves of mixed pixels with various combinations of ground feature abundance fractions. Thus, to ensure unmixing accuracy, an unmixing method for MODIS data based on a set of basis functions, which does not resort to extracting endmember spectra, is proposed and tested in our study.

This study exploits a set of third-order Bernstein basis functions to construct the pixel spectral unmixing indexes (PSUIs), i.e., the coefficients of the basis functions, for a MODIS image without resort to extracting endmember spectra, and then a higher spatial resolution image, such as a Landsat Enhanced Thematic Mapper Plus (ETM+)/Operational Land Imager (OLI) image from the same region and same day with the MODIS image, is utilized to calibrate these indexes, which then creates a calibration model. The calibration model indicates the relationship between the PSUIs and the component abundances and thus can be used for calculating the abundances of the mixed pixel’s components in MODIS images. This method was tested on MODIS and ETM+/OLI images in different scenes or at different times and was compared with other methods, such as the PPI, the N-FINDR, the SMACC, and VCA.

2.

Methodology

2.1.

Bezier Curve and Bernstein Basis Functions

Given a set of control points, $P_i$, $i=0,1,\dots ,n$, its $n$’th-order Bezier curve is defined as

Eq. (1)

$$P(t)=\sum _{i=0}^n{P}_i{B}_{i,n}(t),t\in [\mathrm{0,1}],$$

where $P_i$ is the control point, and $B_{i,n}(t)$ is known as the $n$’th-order Bernstein basis function.⁴⁰

For the $n$’th-order Bernstein basis function, the expansion terms of the binomial expression $1={[t+(1-t)]}^n$ are defined as

Eq. (2)

$$B_{i,n}(t)=\left(\begin{array}{c}n\\ i\end{array}\right)t^i{(1-t)}^{n-i},t\in [\mathrm{0,1}]\text{and}i=\mathrm{0,1},\dots ,n.$$

When $n=3$, it is known as a Bernstein basis function of order 3 (see Fig. 1), which may be defined as

Eq. (3)

$$\begin{gathered} B_{\mathrm{0,3}}(t)=\frac{3!}{0!(3-0)!}{t}^0{(1-t)}^{3-0}={(1-t)}^3, \\ {B}_{\mathrm{1,3}}(t)=\frac{3!}{1!(3-1)!}{t}^1{(1-t)}^{3-1}=3t{(1-\mathrm{t})}^2, \\ {B}_{\mathrm{2,3}}(t)=\frac{3!}{2!(3-2)!}{t}^2{(1-t)}^{3-2}=3t^2(1-t), \\ {B}_{\mathrm{3,3}}(t)=\frac{3!}{3!(3-3)!}{t}^3{(1-t)}^{3-3}=t^3. \end{gathered}$$

Fig. 1

Curves of third-order Bernstein basis functions (see the text).

In a plane or in a higher-dimensional space, the explicit form of this cubic Bezier curve with four control points can be written as

Eq. (4)

$$P(t)=P_0{B}_{\mathrm{0,3}}(t)+P_1{B}_{\mathrm{1,3}}(t)+P_2{B}_{\mathrm{2,3}}(t)+P_3{B}_{\mathrm{3,3}}(t),t\in [\mathrm{0,1}].$$

2.2.

Calculation of Pixel Spectral Unmixing Indexes

2.2.1.

Calculation principle

To calculate the PSUIs for a MODIS image, a set of samples needs to be taken from the image. Here, a MODIS image of the Pearl River Delta region of China was used [MOD021KM: level 1b calibrated, $1000\times 1000\mathrm{m}$ spatial resolution, date: 2001324, time: 03:10, derived from the Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC)]. Acquiring pure pixels containing only one ground object from a MODIS image is difficult because of its spatial resolution limit, but it is possible for each sampling pixel to be dominated by only one category of ground features. For the convenience of discussion, in this paper, such sampling points are called pseudo-MODIS pure pixels, and the ground features estimated from the pseudo-pure pixels are called quasiground features. There are four main kinds of spectral reflectance curves for quasiground features (e.g., water body, sediment-laden water, vegetation, and bare soil) derived from a MODIS image.

Figure 2 shows the spectral reflectance curves of the four types of quasiground features obtained from the sampling points for the above-mentioned MODIS image (a total of 280 samples, each category accounts for a quarter of the total sampling points) and the spectral reflectance curve of a random mixed pixel in the MODIS image, where the original reflectance curve obtained from the MODIS data has been normalized with respect to the total area enclosed by the curve and the $x$-axis. Normalization offers the advantage that it can reduce statistical fluctuations without losing any information. Each curve in Fig. 2 contains 13-channel reflectance data points distributed in a wavelength range from 405 to 2155 nm. Channels 13 to 18 and 26 are not used in Fig. 2, because channels 13 to 16 and 26 are invalid on land and the wavelength ranges of channels 17 and 18 overlap with that of channel 19. According to the spectral reflectance curves of the quasiground features derived from MODIS data shown in Fig. 2(a), different quasiground features reach high reflectance in different channels, e.g., water body in blue–green channels, sediment-laden water in red channels, vegetation in near-infrared (with shorter wavelengths) channels, and bare soil in near-infrared (with longer wavelengths) channels. The peak feature of spectral reflectance curves is important for identifying different ground features. Based on this point, the reflectance data on each spectral curve can be divided into four groups [see Fig. 2(a)], including $G_0$ at wavelengths from 405 to 565 nm, $G_1$ at wavelengths from 620 to 876 nm, $G_2$ at wavelengths from 915 to 1250 nm, and $G_3$ at wavelengths from 1628 to 2155 nm. This way of grouping reflectance data guarantees that high reflectance of the quasiground features appears in different groups. In addition, according to the property of Bernstein basis functions, $B_{i,n}(t)$ reaches a maximum when $t_{\mathrm{i}}=i/n$, which means that the peaks of different basis functions appear at different $t$-values. Both the Bernstein basis function curves and the spectral curves of the ground features have evident peak features. Thus, the third-order Bernstein basis functions with four curves (Fig. 1) are used to characterize the spectral signatures of mixed pixels in MODIS data by using their coefficients.

Fig. 2

(a) Spectral reflectance curves of four types of quasiground features that are derived from a MODIS image (280 samples, each category accounts for a quarter of the total sampling points). $G_0$, $G_1$, $G_2$, and $G_3$ are groupings of the spectral reflectance data. (b) Spectral curve of a mixed pixel in a MODIS image. Gray-shaded areas with the names $S_0$, $S_1$, $S_2$, and $S_3$ were used to show the spectral integral areas corresponding to these four groups of the mixed pixel. Red rectangles below the horizontal axis indicate the locations of MODIS channels.

A cubic Bezier curve from a linear combination of the third-order Bernstein basis functions consists of innumerable data points, whereas a spectral reflectance curve derived from MODIS data consists of 13 data points. Consequently, the spectral reflectance curve of each mixed pixel should be mapped to a cubic Bezier curve before employing the third-order Bernstein basis functions to characterize the spectral curve with their coefficients. A cubic Bezier curve mapped to the mixed spectral curve can be expressed as

Eq. (5)

$$f:F(\lambda )\to P(t),$$

where $F(\lambda )$ represents the spectral curve of a mixed pixel, and $P(t)$ represents the mapped cubic Bezier curve, which is determined by four control points.

Here, the spectral integral area ($S_0$, $S_1$, $S_2$, and $S_3$), which is in the range enclosed by the spectral curve for each group and the $x$-axis [see Fig. 2(b)], and the $t$-value ($t_{\mathrm{i}}=i/n$, $i=0$, 1, 2, 3 and $n=3$) are used together to generate four data points for the mapped cubic Bezier curve. The spectral integral area, which is a combination of sequentially related channels, is employed to replace the single reflectance value. This is because data points generated by 4 of the 15 valid channels of the MODIS sensor cannot fully reflect information about channel width and interrelation, whereas data points generated by the spectral integral areas can do so. Thereafter, four control points can be determined by these four data points. Thus, the cubic Bezier curve is determined.

The components in the LSMM are endmembers with physical meaning, and the abundances are nonnegative. The components in Eq. (4) are third-order Bernstein basis functions, namely $B_{\mathrm{0,3}}(t)$, $B_{\mathrm{1,3}}(t)$, $B_{\mathrm{2,3}}(t)$, and $B_{\mathrm{3,3}}(t)$, which have exact shapes. The coefficients in Eq. (4), namely $P_0$, $P_1$, $P_2$, and $P_3$, express the content of the four basic functions for the mixed spectrum and can be positive or negative. Because the geometric shapes of the four basic functions are invariant, $P_0$, $P_1$, $P_2$, and $P_3$ can objectively describe the complex spectral curves of mixed pixels in MODIS images. Here, these coefficients are called PSUIs.

2.2.2.

Calculation process

There are four steps used to generate PSUIs ($P_0$, $P_1$, $P_2$, and $P_3$) for each mixed pixel in a MODIS image.

• Step 1: Divide the spectral reflectance data of each pixel in the MODIS data into four groups ($G_0$, $G_1$, $G_2$, and $G_3$) according to the peak locations of the spectral curves of the four types of quasiground features [see Fig. 2(a)].

• Step 2: Calculate the spectral integral areas corresponding to these four groups for each pixel as $S_0$, $S_1$, $S_2$, and $S_3$, respectively.

  $$\begin{gathered} S_0=[(R_8+R_9)({\lambda }_9-{\lambda }_8)+(R_9+R_3)({\lambda }_3-{\lambda }_9)+(R_3+R_{10})({\lambda }_{10}-{\lambda }_3) \\ +(R_{10}+R_{11})({\lambda }_{11}-{\lambda }_{10})+(R_{11}+R_{12})({\lambda }_{12}-{\lambda }_{11})+(R_{12}+R_4)({\lambda }_4-{\lambda }_{12})]/2, \end{gathered}$$

  Eq. (6)

  $$\begin{gathered} S_1=(R_1+R_2)({\lambda }_2-{\lambda }_1)/2, \\ {S}_2=(R_{19}+R_5)({\lambda }_5-{\lambda }_{19})/2, \\ {S}_3=(R_6+R_7)({\lambda }_7-{\lambda }_6)/2, \end{gathered}$$
  where $R_i$ represents the reflectance (%) at the $i$’th channel of a pixel, and ${\lambda }_i$ represents the central wavelength (nm) of the $i$’th channel:

  Eq. (7)

  $$\begin{gathered} S=S_0+S_1+S_2+S_3, \\ {S}_0^{\prime }=S_0/S, \\ {S}_1^{\prime }=S_1/S, \\ {S}_2^{\prime }=S_2/S, \\ {S}_3^{\prime }=S_3/S. \end{gathered}$$

  To reduce statistical fluctuations without losing any information, $S_0$, $S_1$, $S_2$, and $S_3$ are normalized to be dimensionless [Eq. (7)]. Hereinafter, $S_0$, $S_1$, $S_2$, and $S_3$ represent the normalized values of the spectral integral areas, respectively.

• Step 3: According to the property of Bernstein basis functions that $B_{i,n}(t)$ reaches a maximum when $t_i=i/n$, and considering the importance of the peak feature of spectral curves, we set $t_0=0$, $t_1=1/3$, $t_2=2/3$, and $t_3=1$; then, four data points of a cubic Bezier curve are generated as $(t_0,S_0)$, $(t_1,S_1)$, $(t_2,S_2)$, and $(t_3,S_3)$.

• Step 4: Substituting $t=t_i$, $P(t)=S_i$, $i=\mathrm{0,1},\mathrm{2,3}$ into Eq. (4), we get

  Eq. (8)

  $$\begin{gathered} S_0=P_0{B}_{\mathrm{0,3}}(t_0)+P_1{B}_{\mathrm{1,3}}(t_0)+P_2{B}_{\mathrm{2,3}}(t_0)+P_3{B}_{\mathrm{3,3}}(t_0), \\ {S}_1=P_0{B}_{\mathrm{0,3}}(t_1)+P_1{B}_{\mathrm{1,3}}(t_1)+P_2{B}_{\mathrm{2,3}}(t_1)+P_3{B}_{\mathrm{3,3}}(t_1), \\ {S}_2=P_0{B}_{\mathrm{0,3}}(t_2)+P_1{B}_{\mathrm{1,3}}(t_2)+P_2{B}_{\mathrm{2,3}}(t_2)+P_3{B}_{\mathrm{3,3}}(t_2), \\ {S}_3=P_0{B}_{\mathrm{0,3}}(t_3)+P_1{B}_{\mathrm{1,3}}(t_3)+P_2{B}_{\mathrm{2,3}}(t_3)+P_3{B}_{\mathrm{3,3}}(t_3). \end{gathered}$$

Solving Eq. (8) for $P_0$, $P_1$, $P_2$, and $P_3$, we have

Eq. (9)

$$\begin{gathered} P_0=S_0, \\ {P}_1=(18\times S_1-9\times S_2-5\times S_0+2\times S_3)/6, \\ {P}_2=(18\times S_2-9\times S_1-5\times S_3+2\times S_0)/6, \\ {P}_3=S_3. \end{gathered}$$

Figure 3 shows the flowchart for employing the third-order Bernstein basis functions to characterize the spectral signatures of mixed pixels in MODIS data by using their coefficients (PSUIs).

Fig. 3

Flowchart for characterizing the spectral signatures of mixed pixels in MODIS data by using cubic Bezier curves.

Here, a Terra MODIS image (date: 2001324, time: 03:10) of the Pearl River Delta region of China was taken as an illustrative example of decomposing mixed pixels. After preprocessing the MODIS image (e.g., geometric correction and cloud masking⁴¹), the PSUIs were obtained using Eq. (9) for the mixed pixels (Fig. 4). Figure 4(a) presents the pseudocolor image derived from channels 7, 2, and 1 of the MODIS data. Figure 4(b) shows the distribution of the normalized difference water index (NDWI),⁴² which is used to evaluate the water distribution information in remote sensing applications,^43,44 whereas the normalized difference vegetation index (NDVI) is usually used to evaluate green coverage and vegetation growth [Fig. 4(c)]. Figure 4(d) shows the distribution of the normalized difference soil index (NDSI),⁴⁵ which is used to enhance soil information. The PSUIs, namely $P_0$, $P_1$, $P_2$, and $P_3$, are shown in Figs. 4(e)–4(h), respectively. The index $P_0$, which mainly reflects the distribution information for the $B_{\mathrm{0,3}}(t)$ function, can be used to identify the distribution of water, as NDWI does. The correlation coefficient between $P_0$ and NDWI is 0.98. The index $P_2$, which reflects the distribution information for the $B_{2,3}(t)$ function, may be used to estimate vegetation growth and to evaluate green coverage, as NDVI does. The correlation coefficient between $P_2$ and NDVI is 0.94. The index $P_3$, which reflects the distribution information for the $B_{\mathrm{3,3}}(t)$ function, can be applied to estimate the distribution of bare soil or outcropped areas, as NDSI does. The correlation coefficient between $P_3$ and NDSI is 0.98. The index $P_1$, as the coefficient of the $B_{\mathrm{1,3}}(t)$ function, can be correlated well with sediment-laden water, and it may have a potential application in estimating the sediment content of water. Figure 4 reveals that the $B_{\mathrm{0,3}}(t)$, $B_{\mathrm{1,3}}(t)$, $B_{2,3}(t)$, and $B_{\mathrm{3,3}}(t)$ functions can reflect information about water body, sediment-laden water, vegetation, and bare soil, respectively, by their coefficients ($P_0$, $P_1$, $P_2$, and $P_3$). Thus, the third-order Bernstein basis functions can be employed to characterize the spectral curves of mixed pixels in MODIS data with physical meaning, which is superior to principal components.

Fig. 4

(a) Pseudocolor image derived from MODIS data of the Pearl River Delta region of China (date: 2001324, time: 03:10. RGB: bands 7, 2, and 1); (b) NDWI from the MODIS data; (c) NDVI from the MODIS data; and (d) NDSI from the MODIS data. PSUIs derived from the MODIS data: (e) $P_0$, (f) $P_1$, (g) $P_2$, and (h) $P_3$.

2.3.

Abundance Calculation Based on the Calibration Model

The PSUIs $P_0$, $P_1$, $P_2$, and $P_3$, which are derived from a MODIS image by adopting Eq. (9), indicate the spectral signals from water body, sediment-laden water, vegetation, and bare soil, respectively. Because the PSUIs represent the relative proportions of ground features in each mixed pixel in the MODIS image, they need to be calibrated by means of the reference abundance values of ground features from high spatial resolution remote sensing images (e.g., Landsat ETM+ or QuickBird image) using the FCLS method, which creates a calibration model for calculating the abundances of the components of every mixed pixel in MODIS images. Here, a Landsat ETM+/OLI image is taken as an illustrative example. Because it is difficult to distinguish the sediment-laden water from a water body when classifying an ETM+/OLI image, the sediment-laden water and water body are classified as the same type (water body). Moreover, water body, vegetation, and bare soil are three basic categories of ground features on the earth’s surface,⁴⁶ which means that $P_0$, $P_2$, and $P_3$ contain most of the spectral information of each pixel. Thus, $P_0$, $P_2$, and $P_3$ are used for the calibration model. The steps for calibrating the PSUIs can be considered as follows:

• (1) Acquire the ETM+/OLI image (path/row: 122/044, date: 2001324) corresponding to the date and region of the MODIS image (date: 2001324, time: 03:10; Pearl River Delta region of China) from the USGS Global Visualization Viewer (GloVis), and classify it as a water body, vegetation, or bare soil using the maximum likelihood classification (MLC) method. Many researchers have used the MLC method for Landsat image classification and obtained satisfactory performance,^47–53 and the classification accuracy produced by the MLC method has been found to be comparable to other classification methods, such as support vector machine.^48,49,52

• (2) Collect a series of quasiground feature samples (i.e., water body, vegetation, or bare soil are the main ones in each sample) from the MODIS image. Here, a uniform sampling cell of $3\times 3\text{pixels}$ ($3\times 3\mathrm{km}$) was used for collecting these samples in order to reduce the projection error, and then the corresponding average values of $P_0$, $P_2$, and $P_3$ were respectively calculated for each sampling cell.

• (3) According to the latitudes and longitudes of the four corners of each sampling cell in the MODIS image, project the boundaries of the samples onto the ETM+/OLI image and the ETM+/OLI classification image (see Fig. 5). Then, respectively calculate the percentages of water body, vegetation, and bare soil pixels accounting for the total pixels in each projection scope in the ETM+/OLI classification image, which are taken as the reference abundances that are used to calibrate the PSUIs. The calibration model for the PSUIs can be expressed as

  Eq. (10)

  $$\begin{gathered} Y_{\mathrm{w}}=a_{10}+a_{11}{P}_0+a_{12}{P}_2+a_{13}{P}_3, \\ {Y}_{\mathrm{v}}=a_{20}+a_{21}{P}_0+a_{22}{P}_2+a_{23}{P}_3, \\ {Y}_{\mathrm{s}}=a_{30}+a_{31}{P}_0+a_{32}{P}_2+a_{33}{P}_3, \end{gathered}$$
  where $Y_{\mathrm{w}}$, $Y_{\mathrm{v}}$, and $Y_{\mathrm{s}}$ denote the abundances of water body, vegetation, and bare soil, respectively, which are obtained from the ETM+/OLI classification image; $P_0$, $P_2$, and $P_3$ represent, respectively, the PSUIs; and $a_{ij}$ ($i=1$, 2, 3, $j=0$, 1, 2, 3) are the fitting coefficients.

• (4) In the illustrative example, we took 189 samples from the MODIS and ETM+ classification images (Fig. 6). We substituted the abundance values of water body, vegetation, and bare soil obtained from the ETM+/OLI classification image and the PSUIs ($P_0$, $P_2$, and $P_3$) obtained from the MODIS image for the samples into Eq. (10), and then obtained the fitting coefficient $a_{ij}$ using a least squares method. The calibration model for each pixel in the MODIS image can be expressed as

  Eq. (11)

  $$\begin{gathered} Y_{\mathrm{w}}=0.5377+1.4790\times P_0-0.4161\times P_2-1.2738\times P_3, \\ {Y}_{\mathrm{v}}=1.6038-2.6723\times P_0+1.0573\times P_2-3.2340\times P_3, \\ {Y}_{\mathrm{s}}=-1.1416+1.1934\times P_0-0.6411\times P_2+4.5079\times P_3, \end{gathered}$$
  where $Y_{\mathrm{w}}$, $Y_{\mathrm{v}}$, and $Y_{\mathrm{s}}$ denote the abundances of water body, vegetation, and bare soil, respectively, and $P_0$, $P_2$, and $P_3$ are the PSUIs.

  The test results for the calibration model are shown in Table 1. All of the multiple correlation coefficients are larger than 0.97, indicating that there are significant linear correlations between the PSUIs ($P_0$, $P_2$, and $P_3$) that were derived from the MODIS data and the reference abundances of water body, vegetation, and bare soil that were obtained from the ETM+ classification image. The test results for the calibration model show that all the observed values for the $F$-test are evidently larger than the critical $F$-test value at the 99% confidence level [$F_{0.01}$ (3,185)]. Thus, there is a marked regression relationship between the PSUIs and the abundances of water body, vegetation, and bare soil, which ensures performance accuracy. The significance test for each PSUI shows that all the significance probabilities are larger than 99.00% and that each index has a significant effect on the abundances. Thus, the calibration model is acceptable.

• (5) Calculate the abundance of each component within each pixel in the MODIS image by substituting the PSUIs for every mixed pixel into the calibration model. If a calculated abundance fraction is less than 0, it is set to 0. The sum of the abundance fractions of different ground features within each pixel must be 1. If not, $Y_{\mathrm{w}}$, $Y_{\mathrm{v}}$, and $Y_{\mathrm{s}}$ should be normalized using the following expressions:

  Eq. (12)

  $$\begin{gathered} Y=Y_{\mathrm{w}}+Y_{\mathrm{v}}+Y_{\mathrm{s}}, \\ {Y}_{\mathrm{w}}^{\prime }=Y_{\mathrm{w}}/Y, \\ {Y}_{\mathrm{v}}^{\prime }=Y_{\mathrm{v}}/Y, \\ {Y}_{\mathrm{s}}^{\prime }=Y_{\mathrm{s}}/Y. \end{gathered}$$

• (6) Evaluate the accuracy of the component abundances obtained by decomposing the mixed pixels in MODIS images. In this accuracy evaluation, the error is defined by the difference between the calculated component abundance and the reference component abundance, where to ensure the objectivity of the accuracy evaluation, the reference abundance used to evaluate the accuracy of the calculated component abundances and the abundances employed to calibrate the PSUIs should be taken from the ETM+/OLI classification images at different times or in different scenes.

Fig. 5

Sampling cells of $3\times 3\text{pixels}$ of ground feature mixtures in a MODIS image, ETM+ image, and ETM+ classification image with the MLC method (date: 2001324). Black quadrangles denote the boundaries of the sampling areas. (a) Water body, (b) water body, (c) bare soil, and (d) vegetation.

Fig. 6

(a) Distribution of sampling cells in the MODIS image. Red squares represent the locations of 189 sampling cells. The positions of these sampling cells are projected onto (b) the ETM+ image and (c) the ETM+ classification image (black squares).

Table 1

Test results for the calibration model.

Calculated model	MCC	F-test	F0.01 (3185)	Significance probability (%)
				P0	P2	P3
$Y_{\mathrm{w}}$	0.979	1412.9	3.9	99.88	99.88	99.88
$Y_{\mathrm{v}}$	0.971	1004.8	3.9	99.82	99.74	99.74
$Y_s$	0.977	1303.4	3.9	99.87	99.81	99.56

Note: MCC: multiple correlation coefficients. MCC=ESSTSS, where the total sum of squares (TSS) is calculated by TSS=∑i=1n(Yi−Y¯)2, the explained sum of squares (ESS) is calculated by ESS=∑i=1n(Y^i−Y¯)2, Yi is the reference abundance of the i’th sample, Y^i is the estimated abundance of the i’th sample, Y¯ is the mean of the reference abundances of the samples, and n is the number of training samples (n=189).F-test: The observed value of F-test at 99% confidence level. F-test=ESS/pRSS/(n−p−1), where the residual sum of squares (RSS) is calculated by RSS=∑i=1n(Yi−Y^i)2, p is the number of explanatory variables (p=3), n-p-1 is the number of degrees of freedom, ESS and n are explained above.F0.01(3.185): the critical F-test value at 99% confidence level with a numerator degree of freedom of 3 and a denominator degree of freedom of 185. F0.01(3.185) is obtained from the table for critical values of the F distribution.54

The flowchart for the proposed approach to decomposing mixed pixels in MODIS images is shown in Fig. 7.

Fig. 7

Schematic description of the approach to decomposing mixed pixels in MODIS image data.

From the above, we can see that although Eq. (4) is still a linear mixture model as is the LSMM, this method using the third-order Bernstein basis functions is different from the LSMM-based approaches in that it does not need to resort to extracting endmember spectra. For the sake of convenience, hereinafter, using PSUIs for decomposing mixed pixels in the MODIS images is called the PSUI method.

3.

Experiments

3.2.

Application of the Calibration Model

In this section, we present the application of the calibration model [Eq. (11)] to MODIS images (see Table 2) in different areas or at different times to test the robustness and performance of the PSUI method (E1 to E5). The abundance maps of water body, vegetation, and bare soil for the MODIS images were then obtained (see Fig. 8).

Fig. 8

MODIS images (RGB: bands 7, 2, and 1) (first column) and their abundance maps of water body ($Y_w$) (second column), vegetation ($Y_v$) (third column), and bare soil ($Y_s$) (fourth column) in five experiments (E1 to E5). (a)–(d) E1, (e)–(h) E2, (i)–(l) E3, (m)–(p) E4, and (q)–(t) E5.

Five sets of sampling grids of $3\times 3\text{pixels}$, which were randomly collected from the MODIS images, were taken as test samples to evaluate the accuracy of the calculated abundances in these experiments. The accuracy evaluation results of these five experiments (Table 3) demonstrate good accuracy for decomposing mixed pixels in MODIS images in different areas at different times by using the calibration model. Therefore, the calibration model can be used for MODIS data in different areas or at different times, which means that there is no need to build a calibration model for every MODIS image.

Table 3

Accuracy evaluation of the abundances calculated by using the calibration model.

Experiment name	Ground feature	ME (%)	MAE (%)	P-10% (%)	P-20% (%)	RMSE
E1	Water body	$-2.1$	3.4	90.0	97.3	0.06
	Vegetation	0.9	8.8	62.7	95.7	0.11
	Bare soil	1.2	9.1	63.3	93.0	0.11
E2	Water body	$-1.3$	4.4	88.6	95.7	0.08
	Vegetation	4.0	7.4	71.4	91.0	0.11
	Bare soil	$-2.6$	6.5	75.7	93.3	0.09
E3	Water body	$-0.8$	2.4	94.8	97.4	0.03
	Vegetation	$-2.2$	3.4	85.7	98.3	0.04
	Bare soil	3.0	4.0	85.2	96.5	0.04
E4	Water body	3.6	7.6	91.0	95.7	0.10
	Vegetation	1.7	6.0	80.5	95.2	0.09
	Bare soil	$-5.3$	8.0	68.1	97.1	0.10
E5	Water body	0.8	2.8	96.0	98.4	0.05
	Vegetation	1.5	7.1	70.4	94.8	0.10
	Bare soil	$-2.3$	7.1	69.6	96.4	0.10

Note: ME, mean error; MAE, mean absolute error. Error = calculated abundance – reference abundance, where reference abundance was derived from the ETM+/OLI classification image corresponding to the date and region of the MODIS image from which test samples were collected.P-10% or P-20%: percentage of the samples with error less than 10% or 20% accounting for total samples.RMSE: root-mean-square error.

3.3.

Comparison with Conventional Methods

To examine the effectiveness of the PSUI method, we compared it with the PPI, N-FINDR, SMACC, and VCA methods using the same MODIS image, against the abundance values of the ground features derived from the ETM+/OLI classification image from the same day as the MODIS image. The PPI, N-FINDR, SMACC, and VCA methods are widely applied for endmember extraction due to their light computational burden and clear conceptual meaning.¹⁵ Detailed descriptions of these four methods can be found in the literature.^15,20–25

In this experiment (E6), the MODIS image (date: 2001356, time: 03:10) and ETM+ image (path/row: 122/044, date: 2001356) used for the method comparison are taken from the same area but not at the same time as those used for calibrating the PSUIs (date: 2001324).

The 295 sampling points of $3\times 3\text{pixels}$ that were randomly collected from the MODIS image (date: 2001356) were taken as test samples. As shown in Table 4, the mean error (ME), mean absolute error (MAE), root-mean-square error (RMSE), and root-mean-square abundance angle distance (rmsAAD) obtained by using the PSUI method are obviously smaller than those obtained by using the PPI, N-FINDR, SMACC, and VCA methods. Furthermore, the errors derived from the PSUI method are distributed around 0% and are centralized [see Fig. 9(a)], whereas those derived from the PPI, N-FINDR, SMACC, and VCA methods exhibit a more disperse distribution [see Figs. 9(b)–9(e)]. The accuracy evaluation results demonstrate that the PSUI method outperforms the PPI, N-FINDR, SMACC, and VCA methods.

Table 4

Accuracy comparison of the PSUI, PPI, N-FINDR, SMACC, and VCA methods.

Method	Ground feature	ME (%)	MAE (%)	P-10% (%)	P-20% (%)	RMSE	rmsAAD	Running time (s)
PSUI	Water body	1.6	5.9	81.4	98.3	0.08	0.22	0.21
	Vegetation	2.9	9.1	64.7	90.2	0.12
	Bare soil	$-\mathbf{4.5}$	9.4	64.4	88.1	0.13
PPI	Water body	$-9.6$	10.4	73.2	78.3	0.19	0.57	8.19
	Vegetation	19.7	22.0	39.7	50.5	0.29
	Bare soil	$-10.1$	25.3	29.5	44.4	0.31
N-FINDR	Water body	15.3	20.0	12.5	51.5	0.22	0.38	5.53
	Vegetation	$-8.1$	12.3	54.9	73.2	0.17
	Bare soil	$-7.3$	14.8	40.7	64.7	0.18
SMACC	Water body	0.6	9.2	64.4	81.0	0.13	0.30	—
	Vegetation	$-5.5$	13.7	48.1	69.5	0.18
	Bare soil	4.8	12.7	45.1	77.3	0.16
VCA	Water body	21.0	24.2	7.1	34.6	0.26	0.45	6.42
	Vegetation	$-13.9$	14.7	45.4	69.2	0.20
	Bare soil	$-7.1$	13.7	43.4	75.9	0.18

Note: ME, mean error; MAE, mean absolute error. Error = calculated abundance – reference abundance, where reference abundance was derived from the ETM+/OLI classification image corresponding to the date and region of the MODIS image from which test samples were collected.P-10% or P-20%: percentage of the samples with error less than 10% or 20% accounting for total samples.RMSE: root-mean-square error.rmsAAD: root-mean-square abundance angle distance. It demonstrates an overall measure of performance of each method. rmsAAD=[1N∑i=1N(AADai2)]1/2, where AADai=cos−1(aiTai^‖ai‖‖ai^‖) measures the similarity between reference abundances (ai) and calculated ones (ai^) of sampling grids; N is the number of sampling grids (N=295).The best results of the four algorithms are in bold font in the table.

Fig. 9

Error distributions of water body, vegetation, and bare soil abundances obtained from (a) the PSUI method, (b) the PPI method, (c) the N-FINDR method, (d) the SMACC method, and (e) the VCA method.

In the comparison experiment, the PSUI method and four conventional unmixing methods were performed with an Intel Core i7-8550U CPU running at 1.80 GHz with 8.0 GB RAM. The PPI, N-FINDR, and VCA methods were performed in MATLAB R2017b, and the running time of these methods was 8.19, 5.53 and 6.42 s, respectively. The running time for building and applying a calibration model for the PSUI method was 0.49 and 0.21 s. The SMACC method was performed in the Environment for Visualizing Images software (ENVI 5.1), which had the result that the running time of the method could not be obtained. The running times in Table 4 show that the PSUI method took less time than the PPI, N-FINDR, and VCA methods.

4.

Discussion

The existing methods of decomposing mixed pixels, based on either LSMM or NLSMM, are mainly based on pixel spectral information that is characterized by a single spectral curve composed of discrete data points and require extracting endmember spectra.^{1,15,16,18,20} The procedures adopted by the methods, such as the PPI²¹ and the SMACC,²³ have been quite successful when pure pixels are present in the original image data. However, it is very difficult to find pure pixels containing only one ground object in MODIS images with low spatial resolution. Many authors have argued that there are no pure pixels in remote sensing images with low spatial resolution.^17,27 Miao and Qi³¹ and Plaza et al.¹⁷ suggested that a trend in the hyperspectral imaging community was to design endmember identification algorithms that do not assume the presence of pure pixels to ensure the endmember accuracy and unmixing accuracy.

The PSUI method proposed herein provides a solution that is different from previous work on the effective decomposition of mixed pixels. This method does not need to resort to extracting endmember spectra from MODIS data. It was tested on five sets of MODIS and ETM+/OLI images, and satisfying unmixing results were obtained (see Fig. 8 and Table 3). The calibration model can be applied to MODIS data in different areas or at different times with high accuracy. The PSUI method was also compared with other methods using the same MODIS data, such as the PPI, N-FINDR, SMACC, and VCA, and the experimental results (Table 4) showed that the accuracy of the PSUI method was obviously higher than that of the PPI, N-FINDR, SMACC, or VCA methods.

In the PSUI method, the PSUIs quantify the relative proportions of spectrally distinct signals from several ground features in each mixed pixel of MODIS data, thus the indexes need to be calibrated with the abundance values of the ground features from a high spatial resolution remote sensing image such as Landsat ETM+ image. One might say that since the PSUIs need to be calibrated with the ETM+/OLI classification images, it would be more convenient to use the results from the ETM+ images directly. However, the low temporal resolution of the 16-day revisit cycle of Landsat ETM+ has long limited its use in many fields, such as studying global biophysical processes, understanding changes in the terrestrial carbon cycle, or mapping the quality and abundance of wildlife habitats.^55,56 MODIS visits the globe once or twice per day with coarse resolution of 250 to 1000 m. In addition, the calibration model is applicable for MODIS data in different areas or at different times, which means that there is no need to build a calibration model for every MODIS image. One of the advantages of the PSUI method is that it combines MODIS data of high temporal resolution with Landsat ETM+ data of high spatial resolution, which may be the reason the new method is superior to the PPI, N-FINDR, SMACC, and VCA methods in terms of decomposition accuracy for mixed pixels in MODIS images.

As we all know, there are 15 reflective spectral channels valid on land in a MODIS image, and these are distributed over a wavelength range of 405 to 2155 nm. These 15 reflective spectral channels can reflect the key spectral characteristics of ground features, such as the locations and intensities of absorption and reflection bands, which are obviously demonstrated in a spectral curve. Three very different ground features (i.e., water body, vegetation, and bare soil) having spectral curves that are easily distinguishable based on their peak locations are involved in the unmixing process. Thus, a good unmixing process for the PSUI method can be performed. However, the PPI, N-FINDR, SMACC, and VCA methods were originally proposed for hyperspectral data,^21–24 and thus would not be expected to perform for multispectral data with limited spectral resolution as well as for hyperspectral data. Furthermore, the endmembers in these conventional methods are specific components, i.e., specific types of mineral or vegetation.^21–24 There may be several specific types of vegetation and bare soil in a MODIS image. However, in the method comparison experiment, mixed pixels in MODIS data were decomposed into three general categories of water body, vegetation, and bare soil. Thus, the performance of the conventional methods may be affected.

For the PSUI method, the training samples used to establish the calibration model were derived from a MODIS image and an ETM+ classification image from the same day and in the same area, which were in almost the same atmospheric conditions. Furthermore, the unmixing accuracies of MODIS images without atmospheric correction were good, whether the MODIS images were the same as that used for the calibration model or not (see Table 3). Thus, atmospheric correction was not necessary for the PSUI method, which could save time and reduce workload for time series analysis with MODIS imagery. To examine the effectiveness of the PSUI method, it was compared with the PPI, N-FINDR, SMACC, and VCA methods using the same MODIS image without atmospheric correction. The conventional methods did not perform so well in this comparison experiment because they all required atmospheric correction.^21–24

The PSUI method, which is based on third-order Bernstein basis functions and does not resort to extracting endmember spectra, has been shown to be effective in decomposing mixed pixels in MODIS data. However, it should be noted that this study was the first attempt to decompose mixed pixels by characterizing the spectral curves of the mixed pixels in MODIS data with a set of Bernstein basis functions. There are still some limitations for the PSUI method. First, the PSUI method is now only suitable for decomposition into three general components (water body, vegetation, and bare soil) in images acquired by a coarse resolution multispectral sensor (e.g., MODIS). It would not be able to decompose mixed pixels into specific vegetation or soil types. Future studies should be carried out to apply the PSUI method to much more complicated ground feature situations. There are two situations: (1) if some of the ground features have very similar spectral signatures, spatiotemporal information as well as spectral information from MODIS data should be utilized comprehensively; or (2) if the high reflectance of each ground feature appears at different wavelengths, Bernstein basis functions of a higher order should be utilized. Second, the calibration model, without atmospheric correction, might work only at low aerosol optical depth (AOD), as the shape of the reflectance spectra at the top of the atmosphere would be highly dependent on the AOD. The impact of absorption and scattering of atmospheric aerosol on reflectance data varies with wavelength, which would change the shape of the spectral reflectance curves and should be corrected by an atmospheric correction algorithm (e.g., the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) algorithm⁵⁷). A new calibration model should be built and applied based on MODIS data with atmospheric correction if AOD is high.

5.

Conclusions

In this paper, the PSUI method, which provides a solution that is different from previous work on the decomposition of mixed pixels, was proposed. This method does not need to resort to extracting endmember spectra from MODIS data, which provides a new way of decomposing mixed pixels to assure the unmixing accuracy. In the PSUI method, the spectral integral area that is in the range enclosed by the spectral reflectance curves of ground features and the $x$-axis (in Cartesian coordinates) and a set of third-order Bernstein basis functions are applied to characterize the spectral curves of mixed pixels in a MODIS image, and the derived PSUIs (i.e., the coefficients of the basis functions) are used for representing the spectral characteristics of the mixed pixels. Then the PSUIs are calibrated with the abundance values of the ground features from a high spatial resolution remote sensing image such as Landsat ETM+ image, which creates a calibration model for calculating the abundances of the components of every mixed pixel in MODIS images. The calibration model is applicable for MODIS images in different areas or at different times, which has been proved by the experimental results using five sets of MODIS and Landsat EMT+/OLI images. The PSUI method was compared with four conventional methods, i.e., the PPI, N-FINDR, SMACC, and VCA. And the comparison results show that the PSUI method outperforms the other four methods for decomposing mixed pixels in MODIS data. Although the PSUI method performs well for decomposing mixed pixels in MODIS images with low AOD into three general categories of water body, vegetation, and bare soil, further study is needed to apply the PSUI method to MODIS images with much more complicated ground feature situations or high AOD.