The RTLSR model is the core algorithm used for generating the MODIS BRDF/Albedo products. The principle for retrieving BRDF parameters is as follows: for each pixel, multi-angle valid observations acquired within a 16-day moving window are used to fit a linear combination of three kernel functions—an isotropic kernel, the Ross_Thick volumetric scattering kernel, and the Li_Sparse-Reciprocal geometric-optical kernel—using the least squares method. This fitting process enables the estimation of directional reflectance parameters for the pixel (Li and Strahler, 1986; Roujean et al., 1992; Schaaf and Strahler, 1993). The mathematical formulation of this process is expressed as: R T L S R θ s , θ v , φ r = f i s o K i s o + f v o l K R o s s _ T h i c k θ s , θ v , φ r + f g e o K L S R θ s , θ v , φ r (1)

In Equation 1, K i s o represents the isotropic kernel, typically set as a constant 1; K R o s s − T h i c k is the Ross_Thick kernel; K L S R is the Li_Sparse-Reciprocal kernel; θ s denotes the solar zenith angle; θ ν is the sensor zenith angle; and φ r is the relative azimuth angle between the sun and sensor.

The Ross_Thick kernel is a semi-empirical kernel used in BRDF models to describe the volumetric scattering component of surface reflectance (Roujean et al., 1992). It characterizes the volumetric scattering effect caused by multiple scattering within dense vegetation canopies. Its mathematical expression is as follows: K R o s s _ T h i c k = π / 2 − ξ cos ⁡ ξ + sin ⁡ ξ cos ⁡ θ s + cos ⁡ ϑ v − π 4 (2) cos ⁡ ξ = cos ⁡ θ s ⁡ cos ⁡ θ v + sin ⁡ θ s ⁡ sin ⁡ θ v ⁡ cos ⁡ φ r (3)

In Equations 2, 3, ξ is the phase angle, which represents the single-scattering solution of the classical radiative transfer equation for a horizontally homogeneous vegetation canopy with a uniform leaf angle distribution and equal leaf reflectance and transmittance.

The Li_Sparse Reciprocal kernel is another key component in semi-empirical kernel-driven BRDF models, used to characterize the geometric-optical scattering effects of sparse vegetation or surface structures (Li et al., 1992; Wanner et al., 1995). It is derived from the proportion of illuminated to shaded areas in a scene consisting of randomly placed ellipsoidal tree crowns, where the tree crown shape is defined by the crown height h and the ratio of the vertical to horizontal crown b/r. Its mathematical expression is as follows: K L S R = O θ s , θ v , ϕ r − sec ⁡ θ s ′ − sec ⁡ θ v ′ + 1 2 1 + cos ⁡ ξ ′ sec ⁡ θ s ′ ⁡ sec ⁡ θ v ′ O = 1 π t − sin ⁡ t ⁡ cos ⁡ t sec ⁡ θ s ′ + sec ⁡ θ v ′ cos ⁡ t = h b D 2 + tan ⁡ θ s ′ ⁡ tan ⁡ θ v ′ ⁡ sin ⁡ ϕ r 2 D = tan 2 ⁡ θ s ′ + tan 2 ⁡ θ v ′ − 2 ⁡ tan ⁡ θ s ′ ⁡ tan ⁡ θ v ′ ⁡ cos ⁡ ϕ r cos ⁡ ξ ′ = cos ⁡ θ s ′ θ v ′ + sin ⁡ θ s ′ ⁡ sin ⁡ θ v ′ ⁡ cos ⁡ ϕ r θ s ′ = tan − 1 b r tan ⁡ θ s ′ θ v ′ = tan − 1 b r tan ⁡ θ v ′ (4)

In Equation 4, the MODIS BRDF Albedo algorithm, the tree crown center height is assumed to be 2b, and the horizontal crown radius is b. The angles θ s ′ and θ ν ′ are transformation parameters used to describe the tree crown as a spherical shape, set as θ s ′ = θ s , θ ν ′ = θ ν , respectively.

Ideally, BRDF modeling of the surface should be based on multi-angle satellite observations of the same target pixel at the same moment in time. However, in practical applications, the RTLSR model approximates this ideal observation condition by utilizing multi-temporal observations of the same pixel within a given time window from the MODIS sensor. This approach has two main limitations: first, it assumes that the structural and optical properties of the land surface target remain unchanged within the inversion time period; second, to obtain sufficient bidirectional observations for inverting BRDF kernel coefficients, multiple high-quality observations must be accumulated within a relatively short period (Vermote et al., 2009). Vermote’s study found that the normalized BRDF kernel coefficients can be linearly related to NDVI (Vermote et al., 2009), and later, Franch et al. validated the empirical relationship between spectral reflectance and NDVI (Franch et al., 2019). Based on this, Gao et al. further integrated soil moisture (SM) and NDVI into the RTLSR model, proposing the RTLSR_MP model.The formula is expressed as follows: R T L S R _ M P θ s , θ v , φ r , N D V I , S M = f i s o K i s o + f v o l K R o s s _ T h i c k θ s , θ v , φ r + f g e o K L S R θ s , θ v , φ r (5) f v o l = a 0 S M + a 1 N D V I f g e o = a 2 S M + a 3 N D V I (6)

In Equation 5, the formula, f v o l and f g e o represent the volume scattering kernel coefficient and geometric-optical kernel coefficient respectively, both accounting for NDVI and soil moisture (SM). The parameters a0, a1, a2 and a3 are model coefficients to be estimated.

Compared with the approach of Franch et al. (2014), who used NDVI and the normalized forms of f v o l and f g e o to linearly model and characterize BRDF shape, Gao’s method treats NDVI and SM as variables and employs an empirical model to estimate the dynamic changes of f v o l and f g e o . This approach not only captures the variation of BRDF shape with changes in surface structure and physical properties but also reflects the amplitude characteristics of surface bidirectional reflectance over the entire inversion period. Here, f i s o , treated as a constant over the inversion period, only characterizes the surface reflectance under the nadir solar-viewing geometry at a specific time during the cycle. Its value is expected to approximate the median integral value of the surface nadir reflectance variation over the entire annual inversion period.However, due to this adjustment, the f i s o loses its original physical meaning—namely, the surface bidirectional reflectance under stable conditions at nadir view and illumination. In practice, the f i s o derived from the RTLSR_MP model is treated as a constant, which serves as a compensation and adjustment for the magnitude of reflectance variation trends represented by the two kernel functions.

For the improved method, the RTLSR_MP model is used to fill daily surface reflectance gaps. The specific workflow (see Figure 1) is as follows.

1. Data Quality Control: Using quality control bands from MODIS observation data (e.g., the “Internal Snow Mask” bitmask in MOD09GA data), preprocess the surface reflectance bands into daily cloud-free, snow-free surface reflectance data and input these into the RTLSR_MP model.

For comparative analysis, the snow pixel removal method based on the NDVI threshold (NDVI <0) and the improved method incorporating the “Internal Snow Mask” snow mask are denoted as RTLSR_MP(NDVI) and RTLSR_MP*(NDVI), respectively. Correspondingly, the two improved models driven by EVI and LAI are denoted as RTLSR_MP*(EVI) and RTLSR_MP*(LAI).

The study area covers parts of East Asia and South Asia, spanning longitudes 73.05°E to 136.00°E and latitudes 15.05°N to 50.00°N, as shown in Figure 2. This region experiences significant seasonal climate variations, characterized by cold and dry winters and hot, humid summers, strongly influenced by the East Asian monsoon. The land cover types are diverse and complex, mainly including deciduous broadleaf forests, evergreen broadleaf forests, mixed forests, grasslands, sparse savannas, croplands, and areas of seasonal snow cover, with pronounced phenological changes and distinct snow seasons. Due to the frequent changes in land surface conditions and the significant influence of clouds, snow, and phenology on remote sensing observations, this region provides favorable experimental conditions for evaluating the adaptability of BRDF inversion models and optimization strategies.

To enhance the comparative analysis of the experimental results, in addition to evaluating the overall model inversion performance across the entire study area, this paper also selects representative sample plots of major vegetation types and two typical seasonal snow-covered regions within the area as key objects of analysis. The basic information of the seven representative vegetation sample plots is shown in Table 1.

To enable spatiotemporally continuous MODIS surface reflectance inversion across RTLSR_MP variants, input driving data with robust spatiotemporal continuity are essential. For efficient construction of continuous vegetation parameters (NDVI/EVI/LAI), this study employs the HANTS-GEE scalable software package (Zhou et al., 2023) on Google Earth Engine. Using MOD13Q1, MYD13Q1, and MOD15A3H products, harmonic analysis was applied to generate daily 2020 NDVI/EVI/LAI datasets. Their reconstruction efficacy has been validated in prior studies (Zhou et al., 2021; Zhou et al., 2023). Soil moisture data were sourced from the GRNN-based product developed by Cui et al. (2019), which integrates multi-source remote sensing data (MOD11C1, MOD13C1, 30-m SRTM DEM, and ECV soil moisture).

To standardize data scales and minimize errors from heterogeneous coordinate systems, all input datasets underwent uniform preprocessing. Specifically, data were projected to the WGS84 geographic coordinate system and resampled to a consistent 0.01-degree spatial resolution. Details of the input data specifications are summarized in Table 2.

In this study, the Root Mean Square Error (RMSE) serves as the primary metric for evaluating model accuracy. The formula is as follows: R M S E = 1 n ∑ i = 1 n P i − O i 2 (7)

In Equation 7, n represents the number of samples, i.e., the number of pixels within the selected area. Pi denotes the surface reflectance value retrieved by the model, while Oi represents the actual valid observed surface reflectance value from MOD09GA. RMSE provides a measure of the overall error between the inversion results and the actual valid observations. A smaller RMSE value indicates that the inversion results are closer to the true observations and that the accuracy is higher. To ensure the stability of RMSE and avoid bias caused by insufficient valid observations, the RMSE for a given day is excluded if the proportion of valid MOD09GA observation pixels meeting the conditions of “cloud-free, shadow-free, snow-free, and good quality” is less than 2% of the total pixels.

Considering that different gap-filling methods may not be applicable to all pixels in practical applications, this study introduces the metric of “filling rate” to evaluate the reconstruction efficiency of each method. The filling rate is defined as the ratio of the number of valid pixels successfully simulated by the model to the total number of pixels in the study area. The calculation formula is as follows: F i l l i n g   r a t e = N v a l i d 0 − 1 N t o t a l × 100 % (8)

In Equation 8, N v a l i d 0 − 1 represents the number of valid pixels successfully filled by the method (i.e., pixels with reflectance values within the range of 0–1), and N t o t a l denotes the total number of pixels in the simulated study area.

To evaluate the impact of the enhanced snow masking strategy on surface reflectance retrieval accuracy, this study implemented both the baseline RTLSR_MP(NDVI) method and the improved RTLSR_MP*(NDVI) method incorporating MODIS Internal Snow Mask for reconstructing MODIS band 1-7 surface reflectance throughout 2020. Figure 3 displays the spatial distribution of RMSE results for visible bands (Band1, Band3, and Band4).

By comparing the results shown in Figure 3, we found that both the RTLSR_MP (NDVI) and RTLSR_MP* (NDVI) methods exhibited an overall “low-high-low” trend in RMSE values across the study area as latitude increased. The regions with relatively higher RMSE were mainly concentrated on the Qinghai-Tibet Plateau and mid-to-high latitude seasonally snow-covered regions. Compared to the RTLSR_MP (NDVI) method, the RTLSR_MP* (NDVI) method achieved lower RMSE values in the mid-to-high latitude seasonally snow-covered regions, while both methods showed similar accuracy in low latitude regions. This indicates that stricter data quality control can effectively improve the accuracy of surface reflectance retrieval in seasonally snow-covered regions.

After presenting the spatial distribution characteristics of RMSE for the two methods, this study further quantitatively compared the fitting accuracy of surface reflectance for MODIS bands 1–7 between the RTLSR_MP (NDVI) and RTLSR_MP* (NDVI) methods over the entire study area and within seasonally snow-covered regions (see Table 3). The results show that RTLSR_MP* (NDVI) outperformed RTLSR_MP (NDVI) at both scales.

Across the entire study area, the RTLSR_MP* (NDVI) method achieved slightly lower RMSE values in all bands compared to RTLSR_MP (NDVI), with the average RMSE for bands Band1–Band7 decreasing from 0.0283 to 0.0278, representing a 1.77% reduction. This indicates that, even in the absence of explicit snow mask information, the quality control strategy based on NDVI thresholding for snow pixel removal still has some applicability under an overall acceptable error margin.

In the seasonally snow-covered regions, the advantage of the improved model was more pronounced. The average RMSE of RTLSR_MP* (NDVI) decreased from 0.0377 to 0.0363, a reduction of 3.71%. Notably, in the visible bands (Band1, Band3, and Band4), the accuracy improvements were even more significant: RMSE reductions of 3.5%, 4.6%, and 4.4% respectively across the full study area, and 5.8%, 7.1%, and 6.4% within the snow-covered regions. This demonstrates that the RTLSR_MP* (NDVI) model, constructed by integrating MOD09GA snow mask information, effectively enhances model stability and accuracy in areas strongly affected by snow dynamics.

However, these results also reflect limitations of the RTLSR_MP (NDVI) method when dealing with complex snow conditions. For example, when snow cover is thin or vegetation protrudes through the snow layer, vegetation’s near-infrared reflectance significantly raises the pixel’s NIR value, causing NDVI to increase and possibly approach zero or even become positive (Klein et al., 1998). Additionally, the coarse spatial resolution of MODIS imagery means that a single pixel may be a mixture of snow, vegetation, and bare soil, further causing shifts and uncertainty in NDVI values (Salomonson and Appel, 2004). Under such circumstances, relying solely on NDVI-based snow pixel removal for quality control is limited and may fail to accurately identify and exclude snow interference, thereby impacting model inversion accuracy.

To address the limited inversion accuracy of previous methods in regions with significant vegetation phenological changes, this study attempts to replace NDVI with EVI and LAI—vegetation parameters that better reflect phenological variations and are less affected by saturation effects—as the driving factors in the model. Table 4 presents the RMSE and average RMSE of the model for bands Band 1 through Band 7 under the three different vegetation parameter drivers.

The results show that among the three parameters, the RTLSR_MP* (EVI) method driven by EVI achieves the lowest average RMSE, which is approximately 3.43% lower than the average RMSE of the NDVI-driven model (0.0291). Conversely, the LAI-driven model exhibits a higher average RMSE than the NDVI-driven one, with an increase of about 3.1%. Examining performance across individual bands, the EVI-driven model shows slightly higher RMSE in the non-near-infrared bands compared to the NDVI-driven model, but in the near-infrared bands (Band2 and Band5), which are more sensitive to vegetation changes, the RMSE decreases by 15.8% and 6.9%, respectively, demonstrating a clear advantage. The LAI-driven model also shows RMSE reductions of 10.3% and 4.1% in bands Band 2 and Band 5, respectively. Although its overall RMSE is slightly higher than the NDVI-driven model, it similarly displays better adaptability in the near-infrared bands. These results indicate that in spectral bands with high vegetation sensitivity, EVI and LAI can provide more accurate driving information as vegetation parameters, while in bands with lower vegetation sensitivity, the NDVI-driven model still holds certain advantages.

It is noteworthy that although both EVI and LAI can mitigate the saturation effect of NDVI to some extent, the EVI-driven model performs better overall. This is mainly attributed to EVI’s higher sensitivity and signal-to-noise ratio in areas of dense vegetation. The calculation of EVI relies on blue, red, and near-infrared reflectance from MODIS’s high-quality atmospheric correction products, which ensures good consistency and stability (Didan, 2021). In contrast, although LAI has a clear biophysical meaning, its retrieval depends on complex processes involving multi-band radiative transfer models and lookup tables (LUTs), including multivariate inputs, structure type identification, and fitting discrimination (Myneni et al., 2021) This makes LAI more susceptible to surface heterogeneity and observation conditions, resulting in higher uncertainty. Moreover, the original LAI time series usually exhibit strong fluctuations before HANTS interpolation, with limited usable data points. Consequently, the interpolation may fail to accurately restore phenological trends (Zhou et al., 2023), thereby affecting its stability as a driving factor and the model’s accuracy.

The above statistical results reflect the overall accuracy differences among the models driven by different vegetation parameters. To further explore how these differences distribute across vegetation types, this study presents the spatial distribution of RMSE for the RTLSR_MP model driven by the three parameters in the near-infrared band (Band2), and summarizes the overall accuracy differences across seven typical vegetation types in the study area, as shown in Figure 4.

From the overall spatial distribution, the RTLSR_MP* (EVI) method shows relatively lower inversion accuracy mainly concentrated in regions with dense vegetation cover, such as the Indian Peninsula, Southeast Asia, and southern and northeastern China. Both RTLSR_MP* (EVI) and RTLSR_MP* (LAI) methods exhibit superior inversion performance in these areas.

According to the accuracy statistics for different typical vegetation types, the RTLSR_MP* (EVI) and RTLSR_MP* (LAI) methods show notable advantages in vegetation types with significant phenological changes. Specifically, these two methods exhibit significantly lower RMSE values than RTLSR_MP* (NDVI) in mixed forests, deciduous broadleaf forests, and woody savannas. In evergreen broadleaf forests, savannas, and croplands, their RMSEs are slightly lower than that of the NDVI-driven method. For grasslands, the RMSE values of all three methods are relatively close.

To further validate the above statistical results, this study conducted a time series and spatial distribution analysis of surface reflectance reconstruction accuracy for the three models across seven typical vegetation sample regions (S1–S7).

As shown in Figure 5, in spectral bands with low vegetation sensitivity (e.g., Band1), the RMSE time series of the three methods exhibit relatively stable performance across all sample areas, with no significant differences. However, in spectral bands with high vegetation sensitivity (e.g., Band2), the NDVI-driven model shows significantly greater fluctuations in RMSE time series compared to the EVI- and LAI-driven models in sample areas such as S1 (savanna), S2 (woody savanna), S3 (mixed forest), and S4 (deciduous broadleaf forest). In contrast, in S5 (grassland), S6 (cropland), and S7 (evergreen broadleaf forest), the time series accuracy of the three methods is relatively consistent. Overall, these results are consistent with the preceding statistical analysis.

It should be noted that although regions S6 (cropland) and S7 (evergreen broadleaf forest) may also be affected by NDVI saturation, the three models still exhibit comparable reconstruction accuracy in these areas. This may be attributed to the following reasons: in cropland region, NDVI saturation typically occurs during the peak growing season, which is relatively short in duration and often coincides with periods of high cloud cover. As a result, the saturated observations contribute less weight to the inversion process. In the case of evergreen broadleaf forests, although the vegetation coverage is high, intra-annual variation is minimal (as shown in Table 1, the annual NDVI standard deviation is only 0.04, less than one-fifth of that in deciduous broadleaf forests), making the impact of NDVI saturation relatively weak in this vegetation type.

Figure 6 further illustrates that in typical vegetation sample regions such as S1–S4, the EVI- and LAI-driven models significantly outperform the NDVI-driven model in Band2, while in regions S5–S6, the differences among the three models are not significant (Note: The gaps observed in the RTLSR_MP*(LAI) results are caused by the failure of interpolation due to an insufficient number of valid data points in the LAI time series available for the HANTS algorithm). This further verifies, from a spatial perspective, the limitation of the NDVI-driven RTLSR_MP model in regions with significant vegetation phenological changes, where it is prone to saturation effects.

Through the work presented in Sections 4.2, 4.3, the improvements in seasonally snow-covered regions and regions with significant vegetation phenological changes have been preliminarily established. Building on this, we will further evaluate the differences between the RTLSR_MP*(EVI) method, the original RTLSR_MP(NDVI) method, and the MODIS surface reflectance derived from MCD43A1 BRDF parameters.

First, an example of reconstruction performance across different seasons is shown in Figure 7. It can be observed that in southern China and the South Asia region (including the South Asian subcontinent and Southeast Asia), cloud cover dominates the imagery throughout the year. The surface reflectance derived from the MCD43A1 Full inversions parameter (Schaaf, 2021), which is based on the RTLSR model, demonstrates certain gap-filling capability in high-latitude regions but exhibits a large number of invalid pixels in mid-to low-latitude areas. The alternative MCD43A1 Magnitude inversions parameter (Strugnell and Lucht, 2001) partially complements these gaps in mid-to low-latitude regions, yet substantial data voids remain in southern China and South Asia. In contrast, both the improved RTLSR_MP*(EVI) and the original RTLSR_MP(NDVI) methods show considerable data gaps in high-latitude regions during spring and winter but outperform the MCD43A1 products in gap-filling ability at mid-to low-latitudes. Notably, in summer and autumn, the reconstruction quality of RTLSR_MP*(EVI) and RTLSR_MP(NDVI) further improves, demonstrating superior surface reflectance retrieval in these seasons.

As shown in Figure 8, we further quantified the time series of the fill rates in 2020 for MODIS surface reflectance reconstructed by the RTLSR_MP*(EVI) method, the original RTLSR_MP(NDVI) method, and the MCD43A1 products. The results indicate that throughout 2020, the daily filling rates of RTLSR_MP*(EVI) and RTLSR_MP(NDVI) were generally higher than those of the MCD43A1 Full inversions and MCD43A1 Magnitude inversions parameters. At the beginning and end of the year, due to the influence of snow, all four models showed relatively low filling rates (Cui et al., 2019; Schaaf, 2021). During the middle of the year, between day 125 and day 315, the filling rates of RTLSR_MP and RTLSR_MP* ranged between 70% and 80%, whereas the filling rate of the MCD43A1 Full inversions parameter significantly declined, reaching its minimum of around 40% near day 200. The filling rate of the MCD43A1 Magnitude inversions parameter fluctuated between 20% and 40% for nearly the entire year.

The MCD43A1 Full inversions parameter is derived from the RTLSR model applied over a 16-day period. Consequently, in high-latitude regions with relatively low cloud cover, sufficient valid observations (≥7) are usually available within the 16-day window to successfully drive the inversion model. However, in mid-to low-latitude regions with heavier cloud cover or during summer and autumn, the number of valid observations within 16 days often falls below this threshold (<7), resulting in missing data for the MCD43A1 Full inversions product. In such cases, valid observations fewer than 7 days can be utilized to construct the MCD43A1 Magnitude inversions parameter, which dynamically adjusts based on a prior database to partially compensate for data gaps. Nevertheless, when no valid observations are available within the 16-day period, the model cannot perform any form of inversion, leading to complete data loss for the corresponding pixels during that time interval.

Compared with the MCD43A1-based gap-filling methods, the RTLSR_MP(NDVI) and RTLSR_MP*(EVI) approaches incorporate multi-source surface parameters with stronger spatiotemporal continuity, which to some extent reduces the RTLSR model’s dependence on the number of valid observations. This relaxation of the inversion requirements allows the model to function effectively even in persistently cloudy regions or periods with fewer observations, provided the inversion window is relatively long, thereby improving the spatial coverage (Gao et al., 2020).

However, the effectiveness of this approach strongly depends on the spatiotemporal completeness and quality of the external driving parameters. When the input parameters contain substantial gaps or poor-quality data, model inversion may still fail. For example, due to widespread missing soil moisture (SM) data in high-latitude regions during winter, both RTLSR_MP*(EVI) and RTLSR_MP(NDVI) failed to perform inversions in northern parts of the study area. Similarly, in arid and desert regions of northwestern China, vegetation parameters (such as EVI or NDVI) still exhibit long periods of missing data even after interpolation using the HANTS algorithm, making inversion in these regions difficult.

Finally, to evaluate reconstruction accuracy, we used the MCD43A1 Full inversions parameter as the reference for comparing the performance of RTLSR_MP*(EVI) and RTLSR_MP(NDVI) due to its relatively high quality, consistent input data, and model similarity. The results are presented in Table 5.

The results indicate that the improved RTLSR_MP*(EVI) method achieves higher accuracy, with RMSE values across all spectral bands more closely aligned with those of the MCD43A1 Full inversions parameter. Specifically, the RTLSR_MP*(EVI) method yields an average RMSE that is 0.012 lower than that of the original RTLSR_MP(NDVI) method, representing an improvement of approximately 4.5%. The enhancements are particularly notable in vegetation-sensitive bands Band2 and Band5, where the RMSE is reduced by 0.055 (approximately 14.3%) and 0.017 (approximately 6.3%), respectively.