Establishing the Downscaling and Spatiotemporal Scale Conversion Models of NDVI Based on Fractal Methodology

Abstract

Scale effect is a crucial scientific problem in quantitative remote sensing (RS), and scholars attempt to solve it with scale conversion models, which can characterize the numerical relationship of RS land surface parameters at different resolutions (scales). As a significant land surface parameter, scale conversion of normalized difference vegetation index (NDVI) has been studied for a long time. Therefore, taking NDVI as an example, the development of scaling research is described and analyzed in the paper, and based on fractal theory, the development trends are discussed for land surface parameters in quantitative remote sensing. These are our conclusions: it will be the new trend to establish downscaling models based on fractal theory for land surface parameters in quantitative remote sensing; additionally, it still is the hotspot to establish spatiotemporal scale conversion models for land surface parameters in quantitative remote sensing in the future, and addressed on that, the multi-fractal scaling methodology is proposed, and its availability is analyzed in the paper, which presents significant potential.

1. Introduction

The scale problem is one of the important and fundamental problems of quantitative remote sensing [1, 2, 3]. Scholars have studied the scale effects of different remote sensing (RS) land surface parameters. The study of scale effect is conducive to the synergistic use of RS data of different spatial and temporal resolutions (scales) to solve the problem that “massive” RS images cannot be fully utilized and has important application potential and scientific research value [3]. In view of the spatiotemporal characteristics of the ground objects, the RS land surface parameters not only have spatial scale effects but also have temporal scale effects. Scholars have conducted extensive and in-depth research on the scale effect of land surface parameters, which includes the mechanism, manifestation, effect analysis, and solution of scale effects. The author has previously discussed it in detail [4]. Based on the above research aspects, scale conversion as a solution to scale effects has received attention. The scale conversion model can characterize the numerical or physical relationship of RS land surface parameter images at different resolutions (scales) and can quantitatively describe scale effects. This paper will also focus on the research progress of the spatial down-scaling and the spatiotemporal scaling.

2. Downscaling of NDVI based on fractal IFS

2.1 Review of downscaling of RS land surface parameters

Liang [1] has reviewed several current downscaling methods, including linear decomposition methods and nonlinear statistical decomposition methods, methods for generating continuous regions, normalized difference vegetation index (NDVI) time series decomposition, multi-resolution data fusion, the statistical downscaling method of global climate model products (GCM), etc. Further, Gao et al. [5], Zhu et al. [6], and Huang et al. [7, 8] have done systematic and effective work in the spatiotemporal fusion downscaling of land surface reflectance, which has become a research hot topic. The spectral-spatial feature fusion by Wang et al. [9, 10, 11, 12] and Shi and Wang [13] also achieved good results for subpixel mapping. These studies, however, scarcely considered the scale conversion process from the perspective of dynamics, which studies of surface parameter downscaling based on the fractal iterated function system (IFS) have paid attention to.

As a fractal branch of mathematics, because of its complete and rigorous theoretical system, it can systematically study the performance, nature, and causes of multi-scale characteristics of natural phenomena. In the fractal geometry theory system, in addition to the familiar fractal phenomenon description and fractal measurement, the internal causes or dynamic processes of mathematical fractals (interaction, feedback, and iteration, represented by IFS-iteration function system) and the physical causes of statistical fractals (such as critical or abrupt changes) are also important research contents of fractal geometry, and fractal geometry has become a part of nonlinear dynamics research [14]. Although the current research on fractal dynamics has just started, there are still many problems waiting to be solved, but its potential value and significance in dynamics research cannot be denied.

In quantitative remote sensing research, fractal methods are mostly used in the mapping of surface morphology (spatial structure) such as active radar imagery and snow and ocean imagery [15], but it also has important applications in scale conversion research and is further deepened and expanded. The use of fractals for surface parameter scale conversion modeling usually contains two important research components:

1. The performance of fractal features, that is, fractal metrics, and also the fractal dimension of the research object. For example, Zhang et al. [16, 17] used the information dimension method to describe the fractal dimension of leaf area index (LAI) scale conversion. Luan et al. [18, 19] and Wu et al. [20] used the similar dimension method to measure the fractal dimension of NDVI and LAI scale-up conversion, respectively.

2. The intrinsic nature of the fractal phenomenon, that is, the dynamics produced, which is the combined effect of multifactor surface effects.

The mathematical basis of fractal generation is IFS. Kim and Barros [21] first constructed the r function from the dynamic factors (soil sediment content, vegetation water content) of soil moisture scale conversion and then established the IFS to describe the soil moisture downscaling, and the conversion effect was good. The model can describe the dynamic process of soil moisture scale conversion, which has physical significance and demonstrates the advantages of downscaling surface parameters based on fractal IFS. In general, there is currently little research into the causes of fractal dynamics. In mathematics, the fractal IFS is a continuously iterative calculation based on the whole research object [14], and the RS land surface parameter image is created in units of local pixels. This ensures that the mathematical IFS vertical conversion factor (r function) is usually constant [21], while the vertical conversion factor of RS land surface parameters (such as soil moisture) is based on the physical elements of each pixel (such as sandy soil). The amount of space and the vegetation water content varies dynamically and temporally [21]. This is why the IFS function can describe the scale switching dynamics of surface parameters and why the model has certain physical meanings. The vertical conversion factor is used to describe the interscale conversion of surface parameter values and is the key to determining the IFS function. Different surface parameters have different values due to the spatial distribution and scale conversion factors (or dynamic factors), and the vertical conversion factor (r function) contains different types of variables and function forms. How to determine the r function is the difficulty in determining the IFS function, which is also an important reason why the latter is less frequently applied in descriptions of quantitative RS land surface parameter scale conversion. Therefore, the NDVI downscaling model based on the fractal IFS function can be considered to describe the dynamic process of scale conversion. This research covers a wide area and is of great significance. The following is a description of a preliminary implementation [22].

2.2 Methodology

How does one build an NDVI downscaling model based on the fractal IFS function? The following points need to be considered: first, how to identify the sensitive factors affecting the spatial distribution and scale effect of NDVI for NDVI; second, how to use this sensitive factor to establish the vertical scale conversion factor r function in the IFS and then determine the IFS function to achieve NDVI downscaling; and finally, how to evaluate the downscaled conversion results. The solution incorporating these considerations is described below.

2.2.1 Identify sensitive factors

According to the above description, water body is an important parameter affecting the spatial distribution and scale effect of NDVI; thus it can be determined that the pixel water parameter is one of the important dynamic factors of NDVI scale conversion. In addition, Wen et al. [23] gave a method for albedo conversion from small-scale to large-scale images and used the pixel topographical influencing factors to correct the converted results, which demonstrated that the method was effective for albedo scale conversion of rugged terrain. Considering the close relationship between the surface reflectivity and the surface albedo, and that the surface reflectance is the basic parameter for calculating NDVI, the topographic factor parameter can be determined as one of the important kinetic factors for NDVI scale conversion. Therefore, the important dynamic factors in NDVI spatial distribution and scale conversion are determined to be the pixel water parameters and topographic factors.

2.2.2 Determine the vertical conversion factor r function and establish the IFS function

Referring to Kim [21], IFS formula (1), horizontal transformation formula (2), and vertical transformation formula (3) for large-scale surface parameter pixel downscaling are obtained as follows. The IFS formula is calculated by pixel-by-pixel sliding. Get the full image downscaling results:

$${\left.{\mathit{IFS}}^{i,j}\right|}_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)=\left(p_n\left(x^i\right),q_m\left(,y^j\right),I_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)\right),$$E1

$$\{\begin{array}{c}{p}_n\left(x^i\right)=x_{n-1}^i+\alpha \left(x^i-x_0^i\right)\\ q_m\left(,y^i\right)=y_{m-1}^j+\alpha \left(y^j-y_0^j\right)\end{array},$$E2

$$I_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)=\left(e_{n,m}{x}^i+f_{n,m}{y}^j+g_{n,m}{x}^i{y}^j+r_1\left(x^i,y^j\right)s^{\mathit{ij}}+k_{n,m}\right)\times r_2\left(x^i,y^j\right),$$E3

where ${\left.{\mathit{IFS}}^{i,j}\right|}_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)$ represents the surface parameter of the pixel at the $\left(i,j\right)$ location when the large-scale pixel of the surface parameter is downscaled to the small-scale image of the $n\times m$ dimension; $x^i$, $y^j$, and $s^{\mathit{ij}}$ correspond, respectively, to the x-direction coordinate $p_n\left(x^i\right)$, the y-direction coordinate $q_m\left(,y^j\right)$, and the surface parameter values $I_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)$ of the three-dimensional data of the pixel; $x_{n-1}^i$ and $x_0^i$ represent, respectively, the x-direction starting coordinate of the $\left(i,j\right)$ pixel in the $n\times m$ dimensional small-scale image and the x-direction starting coordinate of the large-scale pixel; $\alpha$ represents the downscaling ratio (small-scale/large-scale, which is less than or equal to 1); $e_{n,m}$, $f_{n,m}$, $g_{n,m}$, and $k_{n,m}$ are, respectively, functions of the x and y coordinates of the lower left corner and upper right corner of the large-scale pixel, the downscaled surface parameter data, and the vertical scale conversion surface function; and $r_1\left(x^i,y^j\right)$ and $r_2\left(x^i,y^j\right)$ represent, respectively, the two different vertical conversion factors in the vertical scale conversion surface function. Reference [21] should be consulted for the parameters or factors not represented in the formula, which will not be explained here. Generally, the $p_n\left(x^i\right)$ and $q_m\left(,y^j\right)$ coordinates of the $\left(i,j\right)$ pixel are obtained by dividing the large-scale pixel equally into 1/α parts, and the $I_{n,m}\left(x^i,y^j,s^{\mathit{ij}}\right)$ calculation is the key. In formula (3), $r_2\left(x^i,y^j\right)$ is the same as the $r_1\left(x^i,y^j\right)$ function, but their argument coefficients are different.

For NDVI, $g_{n,\mathrm{m}}$, $e_{n,\mathrm{m}}$, $f_{n,\mathrm{m}}$, and $k_{n,\mathrm{m}}$ represent the functions of the$\left(n,m\right)$ pixel of the downscaled NDVI image. Based on the special downscaled NDVI 3D values of the four corner pixels, $g_{n,\mathrm{m}}$, $e_{n,\mathrm{m}}$, $f_{n,\mathrm{m}}$, and $k_{n,\mathrm{m}}$ can be calculated as formulas (4)–(11):

$$g_{n,\mathrm{m}}=\frac{s_{n-1,m-1}-s_{n-1,m}-s_{n,m-1}+s_{n,m}-R_1^g}{x_0{y}_0-x_N{y}_0-x_0{y}_M+x_N{y}_M},$$E4

$$e_{n,\mathrm{m}}=\frac{s_{n-1,m-1}-s_{n,m-1}-g_{n,\mathrm{m}}\left(x_0{y}_0-x_N{y}_0\right)-R_1^e}{x_0-x_N},$$E5

$$f_{n,\mathrm{m}}=\frac{s_{n-1,m-1}-s_{n-1,m}-g_{n,\mathrm{m}}\left(x_0{y}_0-x_0{y}_M\right)-R_1^f}{y_0-y_M},$$E6

$$k_{n,\mathrm{m}}=s_{n,m}-e_{n,\mathrm{m}}{x}_N-f_{n,\mathrm{m}}{y}_M-g_{n,\mathrm{m}}{x}_N{y}_M-R_1^k.$$E7

Furthermore,

$$R_1^g=r_1\left(x_{n-1}^i,y_{m-1}^j\right)s_{0,0}^{\mathit{ij}}-r_1\left(x_n^i,y_{m-1}^j\right)s_{N,0}^{\mathit{ij}}-r_1\left(x_{n-1}^i,y_m^j\right)s_{0,M}^{\mathit{ij}}+r_1\left(x_n^i,y_m^j\right)s_{N,M}^{\mathit{ij}},$$E8

$$R_1^e=r_1\left(x_{n-1}^i,y_{m-1}^j\right)s_{0,0}^{\mathit{ij}}-r_1\left(x_n^i,y_{m-1}^j\right)s_{N,0}^{\mathit{ij}},$$E9

$$R_1^f=r_1\left(x_{n-1}^i,y_{m-1}^j\right)s_{0,0}^{\mathit{ij}}-r_1\left(x_{n-1}^i,y_m^j\right)s_{0,M}^{\mathit{ij}},$$E10

$$R_1^k=r_1\left(x_n^i,y_m^j\right)s_{N,M}^{\mathit{ij}}.$$E11

Therefore, the calculation of the $r_1(x^i,y^j)$ function is significant, and $r_1\left(n,m\right)$ ($0\le r_1\le 1$) is used to adjust the NDVI surface roughness. The following treatment focuses on establishing the vertical transformation formula for NDVI, that is, the determination of the r function (containing $r_1\left(x^i,y^j\right)$ and $r_2\left(x^i,y^j\right)$).

Based on the above sensitivity factors, a vertical conversion factor r function can be constructed:

$$r=\gamma \times S_{\text{water}}+\beta \times s+\delta ,$$E12

where $S_{\text{water}}$ represents the pixel water parameter; s represents the topographic information, taking into account the magnitude of the r function; the normalized difference water index (NDWI) and slope (calculated from the digital elevation model (DEM) image) represent, respectively, the water body effect and the topographic influence in the pixel; $\gamma$ and $\beta$ are the coefficients of the two parameters, respectively; and $\delta$ represents the adjustment constant. Two different orders of magnitude of r are calculated as follows:

$$r_1={\gamma }_1\times S_{\text{water}}+{\beta }_1\times s+{\delta }_1,$$E13

$$r_2={\gamma }_2\times S_{\text{water}}+{\beta }_2\times s+{\delta }_2.$$E14

For NDVI, the $\gamma$, $\beta$, and $\delta$ coefficients can be calculated by linear regression between the high-resolution NDVI image and its NDWI/slope images.

Following construction of the r function, formulas (1)–(3) can be solved in combination with other known conditions, and NDVI downscaling can be achieved.

2.2.3 Evaluation of downscaling results

In order to obtain more accurate downscaling results, if the resolution of the low resolution image is too different from the resolution of the target resolution image (such as downscaling from 250 m MODIS NDVI to 30 m NDVI), a hierarchical downscaling method will be adopted. First, the low-resolution surface parameter image is downscaled to an intermediate resolution image, and then the intermediate resolution image is further downscaled to the target resolution image, which can largely guarantee the accuracy of the result.

Referring to the study by Kim and Barros [21], the accuracy of the downscaled results can be evaluated using statistical indicators such as the maximum, minimum, variance, and standard deviation (compared to high-resolution NDVI images). Moreover, the histograms of the downscaled NDVI and true NDVI images were drawn and compared, and their correlation coefficient was calculated. With those indexes, the accuracy of the downscaled images and methodology could be validated.

2.3 Experiment and result analysis

2.3.1 Experiment

As the best indicator of the status of vegetation growth and vegetation coverage, the normalized difference vegetation index is widely used in the study of environmental (climate) changes, crop yield estimation, and other fields. Among existing vegetation index products, the moderate resolution imaging spectroradiometer (MODIS) vegetation index products are highly valued for their ease of use, ready availability, global coverage, and continuous phase. They have been widely used in studies of forest fires [24, 25], grassland vegetation growth [26, 27], drought [28, 29], land desertification [30], and other studies involving ecological environment monitoring. The maximum spatial resolution of MODIS vegetation index products, however, is only 250 m. The validation of this remote sensing land surface parameter is an important issue that cannot be avoided [31, 32, 33] and needs to be carried out by means of scale conversion. The most representative MODIS NDVI product, namely, MOD13 Q1, will be studied in this paper, which will also focus on establishing a downscaled model of NDVI and validating the MOD13 Q1 product based on it.

This is the experiment. A Landsat8 OLI NDVI image (Figure 1) was utilized to validate a MODIS NDVI image (MOD13 Q1, Figure 2) with nearest imaging time in Xiamen, China. Based on the downscaling formulas in Section 2.2, the MOD13 Q1 image of Xiamen was directly downscaled by ⅛ multiples, and the 30 m downscaled NDVI was obtained as Figure 3. The histograms of the original and processed NDVI images are drawn as Figure 4, and the statistics and correlation coefficients of the NDVI images are presented in Table 1. Based on these data, the downscaled results were evaluated and the MOD13 Q1 image was validated.

Figure 1.

30 m OLI NDVI image.

Figure 2.

240 m MOD13 Q1 image.

Figure 3.

30 m downscaled image of MOD13 Q.

Figure 4.

Histograms of original and processed NDVI images. (a) MOD13 Q1, (b) 30 m downscaled MOD13 Q1, (c) OLI NDVI, (d) difference image between Figures 3 and 1 (Figure 3−Figure 1).

NDVI images	Maximum	Minimum	Mean	Variance	Correlation coefficient
MOD13 Q1	0.999100	−0.561500	0.397400	0.001100	—
30 m downscaled MOD13 Q1	0.999950	−0.999793	0.326412	0.206495	0.937449
OLI NDVI	0.984436	−0.961210	0.299418	0.227313
Difference image (Figure 3−Figure 1)	1.672050	−1.088930	−0.029594	0.027916	—

Table 1.

Statistics of original and processed NDVI images.

Note: The correlation coefficient is the one between the downscaled MOD13 Q1 and the OLI NDVI.

2.3.2 Result analysis

By analyzing Figures 1–4 and Table 1, it is found that:

1. Compared with the real 30 m OLI NDVI image, the 30 m downscaled MOD13 Q1 image has smaller differences in maximum value, minimum value, mean value, and variance. The correlation coefficient between the two images is 0.93, which is highly correlated. The overall quality of the NDVI image obtained by downscaling the MOD13 Q1 image is considered to be good, indicating that the overall quality of MOD13 Q1 is good. Comparing Figure 4(a) and (b), there is a certain similarity between the distribution patterns of the two images, which indicates that the downscaled image retains the spatial distribution structure of the original image to a good degree, which proves to some extent that the original MOD13 Q1 image is of good quality. In addition, comparing Figure 4(b) and (c), it is found that the downscaled NDVI image has a higher proportion in the vicinity of the zero value (mainly artificial features) than the real image. In the range of 0.2–0.6, the difference is greater. The downscaled image generally has a higher proportion in this range of values, and the histogram is smoother, indicating that the image recognition of the NDVI difference is not high. Referring to the correlation between Figure 4(b) and (a), it is known that MOD13 Q1 also has these problems within the abovementioned range of values. Analysis of the original MOD13 Q1 image shows that it is a 16-day NDVI composite product, and each pixel takes the maximum value of NDVI within 16 days as the result of the product release. Therefore, the histogram is reasonable to a certain degree in the larger value area. At the same time, the histogram distribution of the difference image indicates that the pixel values are distributed in the range [−1, 1] and the distribution pattern is low on both sides and high in the middle (approximately a value of 0), which also indicates that the downscaled image and the real image are highly consistent.

2. Further, the analysis of Table 1 shows that the maximum value of the difference image exceeds the range of [−1, 1]. This may be due to a certain error which is caused by the MOD13 Q1 and OLI image during preprocessing process (atmospheric correction, geometric correction, etc.), which causes a large difference in pixel values between the MOD13 Q1 downscale image and the OLI NDVI image. However, the analysis of Table 1 shows that the mean and variance of the difference image are small, so the above abnormal situation only occupies a small space and does not affect the overall evaluation conclusion.

According to the above analysis, the overall quality of the MOD13 Q1 downscaled image is good, indicating that the overall quality of MOD13 Q1 is good. In the NDVI range of values from 0.2 to 0.6, MOD13 Q1 is overestimated, and its discrimination ability of NDVI difference is low, which should be taken into account in practical applications.

2.4 Discussion

Based on the fractal iterated function system, downscaling models of remote sensing land surface parameters can be established. The models can then be merged with more ancillary data, which relate to the scale effects of land surface parameters. Therefore, the models are of benefit for obtaining accurate downscaled results.

In summary, although the breadth and depth of the fractal IFS application in establishing RS land surface parameters downscaling models is still insufficient, the inherent physical meaning and the advantages of the dynamic process expression of this method confer great potential on it, which needs further investigation. It is expected to become a new universal method for quantitative downscaling of RS land surface parameters and lead to the discovery of new research methods.

3. Establishing spatiotemporal scale conversion models of RS land surface parameters based on multi-fractal theory and method

3.1 Review of establishing spatiotemporal scale conversion models of RS land surface parameters

The phase is an important feature of RS images. When the phase changes, the spectrum of the objects in the image changes accordingly. Then, the parameters calculated based on the spectral information will also change, such as surface reflectivity, NDVI, and so on. The temporal response of RS land surface parameters will be further reflected in the variation of its spatial scale conversion model, i.e., the phase characteristics of spatial scale effects.

In order to quantitatively characterize the phase characteristics of spatial scale effects, that is, to establish a spatiotemporal scale conversion model (also called a spatiotemporal scaling fusion model), scholars combined the advantages of higher temporal-resolution feature of low spatial-resolution images and higher spatial-resolution feature of medium spatial-resolution images, and carried out a series of studies on spatiotemporal fusion of remotely sensed surface parameters such as surface reflectance [34], land surface temperature [35, 36], vegetation indexes [37], leaf area index [38], and so on. And then Huang et al. [39] reviewed this and presented the systematic achievements in theory and application. From the theoretical basis of spatiotemporal scale conversion fusion (the spatial scale consistency of the time-phase variation model and the time consistency of the spatial downscaling model) to the type division of the spatiotemporal scale conversion fusion algorithms (the algorithms based on features’ components, the algorithms based on surface spatial information, the algorithms based on features’ temporal change, and the combination algorithms of the ones above), and then to the key problems and challenges encountered in existing research (the imaging geometry and radiation characteristics, differences between multi-source RS images, the complexity of subpixel unmixing models, the complexity of features’ temporal change models, etc.), and the possible development trend in future (improvement in the versatility and robustness of the algorithms), he made a detailed and in-depth explanation, so that we have a more comprehensive understanding of the development of spatiotemporal scale conversion fusion research. In fact, in addition to this method, the multi-fractal method has important potential to solve the above problems [21, 40]. The following is an example of NDVI analysis and how to establish a spatiotemporal scale conversion model (or spatiotemporal scaling fusion model) based on multi-fractal theory and method.

3.2 Spatiotemporal scale conversion models of NDVI based on multi-fractal theory and method

As the best indicator of vegetation growth status and vegetation coverage, NDVI has typical phenological characteristics. This means that in the same area where the surface cover type is unchanged, the physiological characteristics and external forms of the plant can change significantly in different growth stages, and this change will be directly reflected in the changes in image spectrum and NDVI of the surface. Furthermore, the NDVI spatial scale conversion model based on RS images of different growth periods (i.e., different phases) will also change. How to effectively reflect the influence of the phase characteristics of RS images on the construction of this model and then construct a more universal NDVI spatial scale conversion model that can be integrated with surface phenological features, namely, NDVI spatiotemporal scale conversion model? This issue has important research value. Kim and Barros [21] proposed the idea of multi-fractal method for multi-temporal remote sensing soil moisture spatial down-scaling model to describe the phase characteristics of soil moisture spatial down-scaling, but did not do specific research.

Referring to the existing knowledge, the specific method of establishing the NDVI spatiotemporal scale conversion model is given here: first, analyze the surface condition of the study area, determine the type of the main cover of the study area, and based on its phenological knowledge, select enough low and medium-high spatial resolution images finely corresponding to important “nodes” of vegetation throughout the growing season; secondly, the NDVI spatial downscaling models for different growing stages “nodes” are constructed based on the down-scaling methods such as fractal IFS; third, according to multi-fractal theory and method, using the time phase as a factor in the fractal dimension calculation method, the models corresponding to each growth stages are “fused” to obtain a unified and full growth period NDVI scale conversion model (i.e., NDVI spatiotemporal scale conversion model). At this time, the time phase (i.e., different growth stages) has been embodied as a parameter in the model. This model is more universal than the downscaling model based on the single phase image. To obtain a medium to high spatial resolution NDVI image of a certain phase during vegetation growth, the corresponding phase and the low spatial resolution NDVI image of the phase are brought into the model calculation. Of course, this method requires the research object to have a more significant phase or time periodicity, and the established spatiotemporal scale conversion model is more accurate.

Besides, there is another method of multi-fractal modeling of NDVI spatiotemporal scaling. The implementation idea is similar with Section 2.2, while the r function changes. The r functional parameters may need to be recalibrated when the spatial distribution of vegetation cover changes obviously with time (e.g., sowing stage, heading stage, maturity stage, etc.). Therefore, r function will be merged with temporal parameters of NDVI distribution, such as LAI. And the multi-fractal model of NDVI spatiotemporal scaling should be a function of NDVI to capture temporal changes in relation to ancillary data such as LAI.

Although the multi-fractal theory and method has advantages in constructing a spatiotemporal scale conversion model of RS land surface parameters, the theory and implementation of this method are more complicated, and few research cases are currently seen. However, this method is expected to become a new method for the construction of spatiotemporal scale conversion model of RS land surface parameters, which is worthy of further study.

4. Conclusions

Taking normalized difference vegetation index (NDVI) as an example, the establishment of scaling models based on fractal theory was described and analyzed in the paper. It was concluded that fractal iterated function system was an effective methodology to establish downscaling models for remote sensing land surface parameters such as NDVI and multi-fractal modeling may be a novel methodology to establish spatiotemporal scale conversion models for land surface parameters such as NDVI in the future.

Acknowledgments

This work was supported by the National Natural Science Foundation of China “Coupling of NDVI’s up-scaling and down-scaling fusing with ground objects classification” (No. 41601350), the Natural Science Foundation of Fujian Province, China “Research on NDVI’s scaling fusing with ground objects classification” (No. 2017J05069), and the “Scientific Research Climbing Plan” Project from Xiamen University of Technology “Spatial Distribution Estimation and Dynamic Monitoring of Soil Organic Matter Based on Multi-source and Heterogeneous Data” (No. XPDKT19010).