Resolving inherent constraints in eutrophication monitoring of small lakes using multi-source satellites and machine learning

Lakes, as critical freshwater resources, play essential roles in sustaining the global hydrological cycle and ecological balance. They contribute significantly to biodiversity conservation, water supply regulation, and local climate moderation^1,2,3. However, existing research has predominantly focused on large lakes, while the ecological value and localized regulatory functions of small lakes—defined as those with <100 km² surface area—have been largely overlooked^4,5. Limited water exchange capacity and a higher susceptibility to external nutrient inputs make small lakes more prone to eutrophication, accompanied by lower self-purification potential and heightened ecological vulnerability⁶. Targets 6.3 and 6.6 of the United Nations Sustainable Development Goals emphasize the need for improving freshwater quality and protecting aquatic ecosystems⁷. In parallel, China’s “10-Point Water Plan” advocates for refined and localized waterbody management⁸. Accordingly, developing high-frequency, high-accuracy remote sensing monitoring systems tailored for small lakes is crucial for effective evaluation and adaptive regulation in water environment governance.

The Trophic State Index (TSI) is a widely used indicator to assess lake trophic status⁹. However, it fails to fully reflect the nitrogen-phosphorus imbalance characteristics in Chinese lakes and reservoirs¹⁰. To address this, the TSI has been modified by incorporating total nitrogen and the permanganate index, thereby creating a more comprehensive evaluation index, the Trophic Level Index (TLI), which is better suited to the unique characteristics of water bodies in China¹¹. The TLI integrates transparency (Secchi depth, SD), chemical oxygen demand by manganese (COD_Mn), total phosphorus (TP), total nitrogen (TN), and chlorophyll-a (Chla), offering a more accurate representation of eutrophication levels in lakes¹².

Recent remote sensing monitoring of lake eutrophication mainly used four methods: empirical models, semi-analytical models, machine learning methods, and deep learning methods^13,14,15. Empirical models are structurally simple and cost-effective, widely used to rapidly estimate water quality parameters¹⁶. However, they are highly sensitive to changes in water optical properties and struggle to adapt to complex or highly heterogeneous lake environments¹⁷. Semi-analytical models enhance the physical interpretability of models by parameterizing water spectral reflectance mechanisms. However, they rely heavily on prior knowledge and parameter calibration, making them difficult to scale for regional applications¹⁸. Although deep learning methods have strong nonlinear fitting capabilities and excel in large-sample, high-dimensional remote sensing scenarios, they face challenges in small-lake remote sensing monitoring due to limited sample sizes, high overfitting risks, and a lack of interpretability¹⁹.

In comparison, machine learning methods effectively explore nonlinear feature relationships under medium and small sample conditions. They offer strong modeling flexibility and good interpretability, providing an important tool for remote sensing water quality inversion²⁰. In January 2025, Hollmann et al.²¹ published in Nature the TPFN (Tabular Prior-Data Fitted Network) model, combining Bayesian priors and neural network architecture, enabling efficient regression modeling under small sample conditions without hyperparameter adjustment, demonstrating strong cross-task generalization ability. However, its applicability in remote sensing with multi-source, multi-dimensional, and temporal data remains unverified. Conversely, eutrophication studies of small lakes based on remote sensing are still significantly constrained by insufficient data. Furthermore, most studies are limited by reliance on a single remote sensing platform, often beset by insufficient spatial resolution, infrequent coverage cycles, and inadequate temporal continuity. These studies rarely provide systematic comparative analyses spanning multiple remote sensing images and models, thus neglecting the relationship between model interpretability, imagery, and TLI^22,23,24,25. These limitations hinder the high-frequency and dynamic monitoring of eutrophication processes in small lakes. Especially for lakes with limited areas and significant environmental changes, obtaining stable and continuous remote sensing data becomes a key obstacle to modeling.

Consequently, this study proposed an scenario-aware modeling framework to address the challenges of monitoring eutrophication in small lakes. This framework combined multiple remote sensing imagery sources, including the Sentinel-2 multispectral instrument (S2), GF-1, HJ-2, and PlanetScope (PS). This method could improve the spatial and temporal continuity of monitoring, ensuring more comprehensive coverage of the study area. Multiple machine learning models were employed, including CatBoost regression (CBR), XGBoost regression (XGBR), TPFN regression (TPFNR), and linear regression (LR). A systematic comparison of these four models evaluated their accuracy and cross-sensor generalization capabilities under multi-source remote sensing imagery. Thus, the model’s features were quantitatively analysed using SHapley Additive exPlanations (SHAP) to develop a robust and universally applicable TLI inversion model.

The specific objectives of this study are: (1) to enhance eutrophication monitoring capabilities for small lakes by integrating multi-sensor data in high-frequency field observations; (2) to investigate how the accuracy and generalization capabilities of different machine learning models vary across different remote sensing imagery and temporal dimensions; and (3) to determine the seasonal and spatial distribution characteristics of Dongqian Lake using TLI. This study provides the technical basis for developing high-resolution eutrophication monitoring systems for small lakes and offers a scientific foundation for formulating water quality management and ecological conservation policies.

Model framework overview

To assess the effectiveness and performance differences among various regression methods for predicting the trophic status of small lakes using different sensors, a scenario-aware modeling framework was developed. Representative machine learning models were selected: CBR, XGBR, TPFNR, and LR. Constructing the input dataset for remote sensing features involved combining features derived from various images with field sampling data. To assess the predictive performance and generalization capability of the models, model evaluation metrics were applied. Ultimately, the model demonstrating the highest accuracy and the lowest error across a range of remote sensing images was chosen for the final TLI inversion. The scenario-aware modeling inversion framework is illustrated in Fig. 1.

Statistical profile of in-situ water quality

The overall mean values of water quality parameters for the entire year were: DO at 8.38 ± 1.38 mg·L⁻¹, SD at 49.02 ± 14.85 cm, Chla at 38.79 ± 29.83 mg·m⁻³, COD_Mn at 4.27 ± 1.21 mg·L⁻¹, the TN/TP at 22.13, and TLI at 37.36 ± 4.99, indicating that the lake was in a mesotrophic state (Table 2) with phosphorus limitation. Correlation analysis of the overall data for the year revealed that TLI had a strong significant correlation with Chla (r ≈ 0.85, p < 0.001) (Fig. 2); significant correlations with TP (r ≈ 0.75, p < 0.01) and TN (r ≈ 0.66, p < 0.01); and a significant correlation with COD_Mn (r ≈ 0.55, p < 0.05). In addition, TLI showed a significant negative correlation with SD (r ≈ −0.70, p < 0.01). However, TLI had no significant correlations with Tem (r ≈ 0.50), pH (r ≈ 0.35), and DO (r ≈ −0.30), all with p > 0.05. To analyze the distribution characteristics of water quality parameters across different seasons, we employed histograms, normal distribution fitting curves, and violin plots for visual representation. Furthermore, we applied the Shapiro-Wilk test to evaluate normality and to emphasize the statistical distribution differences of the parameters on a seasonal scale.

TLI seasonal variations and related parameters demonstrated significant differences, with the distribution of each parameter significantly deviating from normal distribution (p < 0.0001)(Fig. 3). In spring, DQ Lake was generally in a mesotrophic state (TN/TP > 22.6, indicating phosphorus-limitation), with a slow increase in Chla, a relatively low TLI, moderate COD_Mn, higher DO, and good water transparency (SD), indicating stable water quality. In summer, the lake reached its peak eutrophication level, with TN/TP ranging from 9 to 22.6 (dual nitrogen and phosphorus limitation), and some areas dropping below 9 (nitrogen limitation). Chla increased rapidly, and TLI reached a light eutrophic state, with COD_Mn peaking at its highest level of the year. DO showed significant fluctuations throughout the day, and water transparency noticeably decreased. In autumn, the lake began its recovery process after eutrophication, with TN/TP levels remaining between 9 and 22.6, indicating dual nitrogen and phosphorus limitation. Although Chla decreased, it remained at a high level, while TLI exhibited a gradual decline. COD_Mn levels stayed elevated, low oxygen conditions continued, and water transparency remained poor, which indicates a slow recovery process. During winter, the lake reached a state of stabilization, with the TN/TP ratio exceeding 22.6, suggesting phosphorus limitation. Chla hit its lowest levels of the year, and TLI decreased in parallel. COD_Mn reached its lowest level of the year, while DO saw a significant increase. Additionally, water transparency improved markedly, leading to enhanced water quality. These seasonal differences validated the overall trends observed in the annual correlation analysis, further clarifying the seasonal patterns of water quality parameter variations.

Selecting remote sensing features for TLI modeling

To identify the key features suitable for TLI remote sensing modeling, this study analyzed the relationship between features from the S2, GF-1, HJ-2, and PS datasets and TLI, using significance (p < 0.05) and r as criteria (Fig. 4). After outlier removal based on IQR method, 279, 78, 124, and 68 valid sampling points were retained for S2, GF, HJ, and PS data, respectively. In the S2 dataset, NDWI (r = 0.47), NGRDI (r = 0.46), and BGR (r = 0.46) showed strong significant correlations with TLI (p < 0.01). NDCI and CI, showing significant correlations (r = 0.31, p < 0.05), were included in the feature selection process. In the GF-1 dataset, NDWI (r = 0.64, p < 0.01) made the greatest contribution, followed by BGR and NGRDI (r = 0.43, p < 0.01). Although NDTI and some bands (B1-B4) exhibited negative correlations, they were also selected due to their statistical significance. In the HJ-2 dataset, there was a significant correlation between BGR and NGRDI (r = 0.44, p < 0.01), and bands B1, B2, and B5 satisfied the selection criteria (r ≈ 0.31–0.32, p < 0.05). In the PS dataset, the features B8, NDTI, NDCI, CI, GCI, and EVI also met the selection criteria (r ≈ 0.28–0.31, p < 0.05), while other features fell short of the required threshold.

Performance assessment of models with multi-source data

Estimating TLI across different satellite sensors

Table 1 displays the predictive performance of the regression models CBR, XGBR, LR, and TPFNR across various remote sensing datasets. Overall, the CBR model outperformed other models in terms of both prediction accuracy and stability. In the GF-1 and HJ-2 datasets, the CBR model achieved R of 0.94 and 0.83, respectively, with NRMSE of 0.13 and 0.16 and biases close to zero, indicating high prediction accuracy and robustness. In the S2 dataset, the correlation coefficients of the CBR and TPFNR models were similar (R = 0.80 and 0.81, respectively). However, CBR had significantly lower errors (NRMSE = 0.14, Bias = 0.15) than TPFNR (NRMSE = 0.16, Bias = 0.26), showing better prediction stability and accuracy. The high spatial resolution of PS remote sensing imagery led to increased computational errors in the model, which caused a decrease in the prediction accuracy of the CBR model (R = 0.64, NRMSE = 0.19, Bias = −0.09). However, it still surpassed other models, demonstrating its superior adaptability to the data. The XGBR model demonstrated high fitting accuracy on the training set (with R values generally close to 1). However, its generalization performance on the test set was insufficient, especially in the HJ-2 dataset, where the bias was large (Bias = 0.66). In the PS dataset, the NRMSE (0.21) was also higher than the CBR model, suggesting the presence of overfitting. The TPFNR model performed well in some datasets (e.g., S2 and GF-1) but demanded a significantly higher computational cost than other models (with training times reaching 68.5 seconds). Additionally, its performance in the PS dataset dropped substantially (R = 0.50, NRMSE = 0.24), indicating limitations in model complexity and generalization ability. The LR model recorded consistently poor performance across all datasets, with the highest error observed in the PS dataset (NRMSE = 0.24) and the lowest prediction accuracy (R = 0.47), demonstrating its inadequate ability to capture complex nonlinear relationships in remote sensing data.

To further assess the generalization capability of the models, a comparison was conducted between the predicted values of the CBR, XGBR, TPFNR, and LR models and observed values (Fig. 5). The analysis reveals significant differences among the models across various sensor types and sampling time points. Of the models evaluated, the CBR model exhibited the best generalization performance, with its predicted values closely aligning with the 1:1 reference line. For both S2 and GF-1 datasets, CBR achieved high consistency with observations, with no notable overestimation or underestimation tendencies. In the HJ-2 dataset, slight underestimation was observed at higher TLI, while minor overestimation was predicted at lower TLI. In the PS dataset, prediction scatter increased, particularly underestimating at high TLI values. In contrast, XGBR and TPFNR models displayed value-dependent error characteristics, particularly under HJ-2 and PS datasets. Significant underestimation was noted at high values, while predictions at low values deviated substantially from the 1:1 line, indicating potential overfitting. The LR model consistently underperformed across all datasets, with pronounced underestimation at high values and highly scattered results at low values, particularly for GF-1 and PS, leading to considerably higher error levels than other models

Analyzing the contribution of SHAP-based feature

To clarify the contributions of individual remote sensing features to model predictions, SHAP was utilized to evaluate feature importance across four models: CBR, XGBR, TPFNR, and LR. (Fig. 6). The results showed both consistent patterns and notable differences in feature importance, attributed to variations in the spectral characteristics of the data and the structures of the models.

The CBR model consistently relied on vegetation indices and composite spectral features across all sensor types. In the S2 dataset, the primary predictors of CBR were identified as NDWI, BGR, NGRDI, and NDVI. For the GF-1 dataset, CBR mainly depended on BGR, B3, and NGRDI. In the HJ-2 dataset, the important features included NGRDI, BGR, NDTI, and B1. In the PS dataset, NDVI and NDTI emerged as the most influential features for CBR. The XGBR model integrated diverse spectral features to effectively utilize multi-dimensional spectral information. Within the S2 dataset, shortwave bands (B7, B6, B2) and index-based features (NDVI, NDCI, NGRDI) played a significant role in the XGBR model. The GF-1-based XGBR model highlighted a combination of BGR, NDWI, and NDTI, while the HJ-2 XGBR model focused on B6, NDTI, and B1. For the PS dataset, NDVI and NDTI were key variables in the XGBR model. In contrast, the TPFNR model showed strong feature selectivity, heavily relying on mid- to near-infrared bands and specific indices. In the S2 dataset, the key features for the TPFNR model included B7, B6, NDCI, and B11. For the GF-1 dataset, NGRDI and BGR were emphasized in the TPFNR model. The preferred features for the TPFNR model in the HJ-2 dataset included NGRDI, BGR, NDTI, and B1. In the PS dataset, the TPFNR model prioritized GCI, NDCI, and CI.

The LR model showed limited feature utilization across all datasets, with predictive power largely concentrated on a small subset of bands. In the HJ-2 dataset, B3, B5, and B2 accounted for most of the LR model’s variance, underscoring the linear model’s restricted capacity to capture complex spectral patterns.

Assessing TLI spatial and temporal variability

The integration of multi-source imagery effectively compensates for the lack of effective observations from individual sensors on specific sampling dates, ensuring temporal continuity between sampling and image matching, and enabling near-seamless spatiotemporal monitoring. Preprocessed images from S2, GF-1, HJ-2, and PS were fed into the CBR model to generate spatial distribution maps of TLI values for DQ Lake during various periods in 2023 (Fig. 7). Overall, TLI exhibited distinct seasonal dynamics. In early January, TLI values across the lake ranged from 32 to 33.5 and remained relatively low in early spring. From late spring to mid-May, TLI values experienced a significant increase, initially rising from 32 to 39 and later reaching levels between 44 and 49. During the summer months of June and July, TLI values continued to exhibit an upward trend, attaining peak values ranging from 34 to 36.25 and 38 to 46, with the exception of a sharp decline on June 14, 2023. However, by late August, TLI rebounded to higher levels, fluctuating between 37 and 46.

At the spatial scale, TLI displayed a consistent pattern of spatial heterogeneity, characterized by a distinct gradient: the values were higher in nearshore areas and decreased toward the lake center. The northern bay and southwestern coastal areas have consistently been zones with high TLI values at different points in time. In contrast, the central and southern open-water areas exhibited persistently low TLI levels. Notably, TLI values reached a maximum of 44 to 46 in nearshore zones during the summer, with occasional high-value regions emerging in shallow or shoreline areas on specific dates.

In practical remote sensing applications, inversion models generally require data processing from multiple satellite platforms, which exhibit significant discrepancies in sensor specifications, spatial resolution, spectral configurations, and imaging conditions. Therefore, ensuring consistent generalization performance across diverse sensor configurations is paramount for reliable water quality prediction. The results indicate that all models experience a decline in prediction performance when applied across different sensors. However, the extent of degradation and stability present significant variations among the models. Among the evaluated models, CBR demonstrates higher stability and generalization capability across sensors by partitioning the feature space into numerous regions and performing weighted aggregation of similar samples within the same region. The CBR model demonstrates the most consistent performance, exhibiting a reduction in the R ranging from 0.07 to 0.15. It is noteworthy that the model shows relatively high prediction accuracy when transitioning from higher-resolution data, such as GF-1, to lower-resolution data, such as HJ-2, and from S2 to PS. SHAP-based feature importance analysis confirms the model’s consistent reliance on vegetation-related spectral indices (e.g., NDWI, NGRDI, CI, and BGR) across all sensor datasets, indicating strong adaptability to heterogeneous data sources. Conversely, the XGBR model displays notable variations in cross-sensor prediction performance, particularly in scenarios characterized by substantial disparities in spatial resolution. This sensitivity is likely attributable to the model’s high dependence on input feature consistency. The TPFNR model demonstrates notable performance instability across sensors and incurred substantially higher computational costs due to its structural complexity. This finding indicates that the quality and consistency of the data significantly influence the model’s generalization. The LR model demonstrates the least effective cross-sensor generalization capacity. The linear structure of the model proves inadequate in capturing the nonlinear spectral variability inherent in remote sensing data, thereby limiting its predictive capability across different sensor platforms. Consequently, models such as CBR, which demonstrate superior generalization performance, should be prioritized in practical applications of remote sensing-based water quality monitoring. For models showing weaker generalization capacity (e.g., XGBR and TPFNR), it is imperative to exercise additional caution to mitigate the uncertainties introduced by sensor-specific differences.

The variations in model prediction performance can be attributed to various factors, including model structure, preprocessing strategies, image quality, and observed values. Meanwhile, the variations in spectral features of different sensors are primarily influenced by their physical characteristics:

1. (1)

   Differences in band availability and spectral response functions are key factors. For example, S2 and PS are associated with red edge bands centered around 705‒740 nm, whereas GF-1 and HJ-2 satellites lack dedicated red edge coverage. Consequently, red edge dependent indices, NDCI and CI, used in this study are standard configurations for S2 and PS data, but require alternatives, NDCI* for NDVI-type and CI* for NIR/R-1 on HJ-2.

2. (2)

   Due to differences in spatial resolution and shoreline proximity effects, subpixel mixing is more obvious in narrow waterways and along bright shorelines. As a result, indices constructed from red/green bands (e.g., NGRDI, BGR) exhibit increased sensitivity to these non-water signals, leading to reduced robustness and greater susceptibility to positive bias near shorelines.

3. (3)

   Different sensors utilize distinct atmospheric correction workflows. S2 adopt the 6S model, GF-1 and HJ-2 use the FLAASH model, and PS uses the DSF model, respectively. Residual errors persist in the visible–near-infrared (VIS/NIR) bands, due to different assumptions and inputs among the methods, such as incomplete removal of haze path radiation or effects from shoreline proximity.

Figure 7 shows the spatial distribution of the TLI consistently register higher values in the nearshore zone. This phenomenon may be attributed partly to human and natural factors, such as tributary nutrient inputs, enclosed bays, and extensive aquatic grass beds. However, some higher values may be induced by optical interference factors, including shallow-water reflectance, shoreline shading effects, and signals from aquatic vegetation that enhance surface reflectance. Although high-frequency sampling can enhance reliability, some factors may reduce the accuracy, such as nearshore substrate heterogeneity, anthropogenic disturbance, and hydrodynamic influence. Consequently, during the image pre-processing stage, we first excluded nearshore areas by applying a water body mask. After inversion, high-value nearshore regions were designated as “hotspots” and analyzed for their rationality based on spatial patterns and temporal consistency.

Despite the comprehensive evaluation of model performance and cross-sensor generalization using multi-source remote sensing data, some limitations have remained. They include the relatively short temporal span of available data, the limited spatial coverage of in-situ sampling, and the lack of fully integrated multi-source strategies. To further advance eutrophication monitoring in small lakes, future studies should focus on extending the temporal span of remote sensing and in-situ datasets, refining data preprocessing and sampling protocols, and promoting the integrated use of domestic and international satellite resources, thereby improving model accuracy, robustness, and generalizability.

Study area

Dongqian (hereinafter referred to as “DQ”) Lake is a typical small coastal plain lake with shallow depths. Human activities have significantly impacted it, making it an ideal case study for investigating eutrophication in small lakes. Situated in the eastern part of Yinzhou District in Ningbo City, Zhejiang Province, China, it is the largest natural freshwater lake in the province (Fig. 8)²⁶. The elongated lake is 8.5 km long (north-south axis) and 4.5 km wide, with 22 km² water surface area and 45 km shoreline. With an average depth of 2.2 m, the water volume amounts to 33.9 million m³. DQ Lake is classified as a lagoon type. The region has a subtropical monsoon climate, with mild and humid conditions throughout the year, an average annual temperature of 15.4 °C²⁷. The lake is pivotal in providing water for drinking and agricultural irrigation for Ningbo, which is important for regional water resource security and ecosystem services.

In-situ water quality assessments

The water-quality dataset was provided by the Ningbo Ecological and Environmental Bureau, covering the sampling period from January to October 2023. A total of 26 fixed monitoring points were established around the lake (Fig. 8c). Monitoring was conducted 27 times, resulting in 702 samples. Sampling was performed four times per month from May to September, and once per month for the remaining months. Temperature (Tem), Dissolved Oxygen (DO), pH and Secchi Depth (SD) were measured in situ. Chemical Oxygen Demand-Manganese (COD_Mn), Total Phosphorus (TP), Total Nitrogen (TN) and Chlorophyll-a (Chla) were analyzed in the laboratory. Table 2 shows the descriptive statistics for the water quality parameters.

Satellite data acquisition and processing

This study’s satellite imagery primarily comprised S2 data retrieved via Google Earth Engine (GEE), supplemented by HJ-2 and GF-1 imagery obtained from the China Centre for Resources Satellite Data and Application (CRESDA), and PS imagery sourced from Planet Labs (Supplementary Table S1 for details)²⁸. The following processing steps were implemented:

1. (1)

   All imagery was matched within a ± 2-day window of the sampling date for each field survey. When multiple images were available, the one with the closest temporal synchronization, cloud-free condition, and highest spatial resolution was selected (Supplementary Table S2 for details).

2. (2)

   Reliable auxiliary data on aerosols and water vapor were available, and where the distance between land and water bodies was moderate, the 6S reflectance model was utilized with Sentinel-2 satellite data²⁹. MODTRAN effectively handled proximity effects and path radiation. In cases where bright coastlines or densely populated urban areas caused notable proximity effects and haze, the FLAASH model was applied for GF-1 and HJ-2 satellites³⁰. Dark spectral fitting (DSF) was proven to be highly effective in both inland and turbid waters without relying on external aerosol data. Therefore, the DSF albedo model was applied to PlanetScope satellite imagery^31,32.

3. (3)

   Water bodies were extracted using the Normalized Difference Water Index (NDWI)³³.

4. (4)

   To minimize interference from shoreline mixing, image data showing a water proportion less than 0.8 within a 3 × 3 pixel window were excluded.

Methodology for TLI estimation

Outlier detection was performed using the Interquartile Range (IQR) method, which can ensure the inversion accuracy of the model. Each variable’s first (Q1) and third (Q3) quartiles were calculated, and IQR was defined as IQR = Q3 − Q1³⁴. Observations falling outside the range of [Q1 − 1.5 × IQR, Q3 + 1.5 × IQR] were identified as outliers and removed³⁵.

TLI was constructed with Chla (mg m⁻³) as the core parameter, along with SD (cm), TP (mg L⁻¹), TN (mg L⁻¹), and COD_Mn (mg L⁻¹). TLI was calculated based on the national standard³⁶:

n the formula, TLI(j) represents the composite index for parameter j, with the associated weight Wj. The r_ij refers to the correlation coefficient between the reference Chla and each parameter j (r²_Chla = 1, r²_TP = 0.7056, r²_TN = 0.6724, r²_SD = 0.6889, r²_CODMn = 0.6889)³⁷. The final TLI value is calculated as the weighted average of the individual indices. The trophic state classification standard for Chinese lakes (reservoirs) is shown in Table 3³⁸:

In lake ecosystems, the TN/TP ratio is commonly utilized to identify the primary limiting nutrient factor, indicating the relative availability of nitrogen or phosphorus during phytoplankton growth³⁹. A TN/TP ratio below 9 suggests that the water may be nitrogen-limited, which can promote the growth of nitrogen-fixing cyanobacteria. On the other hand, when the TN/TP ratio exceeds 22.6, phosphorus limitation may occur, hindering phytoplankton growth. If the TN/TP ratio falls between these two thresholds, the water body may experience simultaneous nitrogen and phosphorus limitations⁴⁰. Analyzing the TN/TP ratio is crucial for understanding how changes in lake nutrient profiles relate to risks of eutrophication.

To enhance the accuracy of remote sensing data inversion models and to improve their sensitivity to spectral variations, this study developed a multidimensional spectral feature dataset. This dataset was derived from atmospheric-corrected surface reflectance products obtained from various sources of remote sensing imagery. It is important to note that differences exist in the band configurations across different remote sensing platforms (see Supplement Table S3 for details). The dataset mainly includes: (1) Water Optical Characteristics indicators, such as the Normalized Difference Turbidity Index (NDTI) and Normalized Suspended Sediment Index (NDSSI), which primarily reflect changes in water turbidity and suspended particulate matter^41,42. (2) Algal Biological Response indicators, including the Phytoplankton Algae Index (FAI), Chlorophyll Index (CI), Normalized Chlorophyll-a Index (NDCI), and Green Chlorophyll Index (GCI), which characterize the dynamic changes in algal pigment content and chlorophyll-a in water bodies^43,44,45. (3) Key Bands and Ratio Combinations, such as the Blue/Green Ratio (BGR) and various single-band reflectance indices (B1−B12), which help capture water absorption characteristics and spectral curve shape variations⁴⁶. During the feature selection process, the Pearson correlation coefficient (r) was utilized to assess the relationship between each variable and TLI. Variables that demonstrated insufficient significance were eliminated based on a two-tailed t-test (p < 0.05)⁴⁷.

Machine learning models for TLI retrieval. CBR is an improved gradient-boosting method that utilizes ordered target encoding and oblivious decision trees⁴⁸. It efficiently handles heterogeneous features without complex preprocessing, offering strong generalization capabilities and interpretability. XGBR offers advantages such as regularization control, missing value handling, and parallel computation, making it effective for modeling nonlinear relationships⁴⁹. Particularly well-suited for high-dimensional and noisy data, it is widely applied in remote sensing inversion and feature recognition tasks. TPFNR is based on a Prior-Function Network structure, leveraging meta-learning from millions of tabular tasks for model initialization. It does not require manual parameter tuning and can achieve rapid convergence and robust predictions in small-sample scenarios. The classic linear model LR has a simple modeling mechanism and low computational cost, making it suitable for modeling linear or approximately linear relationships⁵⁰.

The dataset was randomly split into training (70%) and testing (30%) sets (Rhodes et al. 2023). For CBR and XGBR, hyperparameter tuning was performed using grid search combined with k-fold cross-validation. The key adjusted parameters were learning rate, L2 regularization (CBR), and gamma (XGBR)⁵¹. TPFNR and LR were automated modeling methods requiring no manual parameter configuration. After model training, the SHapley Additive exPlanations (SHAP) method was used to explain feature importance in each model, enhancing model interpretability⁵². The TPFNR model was implemented using the Python TabPFN package. CBR and XGBR were executed via CatBoost (v1.2.3) and XGBoost (v2.0.3), respectively, while the LR model was supported by scikit-learn (v1.3.1). The model parameter configurations utilized for TLI prediction in multi-sensor experiments are detailed in Table 4.

To systematically evaluate the predictive performance of the models in TLI inversion, this study utilizes four key metrics: Pearson correlation coefficient (r), Normalized Root Mean Square Error (NRMSE), Bias, and Unbiased Mean Absolute Percentage Difference (UMAP)^53,54. These metrics comprehensively reflect model performance from the perspectives of correlation, margin of error, and systematic bias, ensuring the thoroughness and robustness of the evaluation⁵⁵. The following equations were used to calculate:

Where x_p represents the predicted values obtained from the model and x_t denotes the observed values. The Unbiased Difference (UD) is the residuals resulting from the LR fitting between the predicted and observed values. UD reflects how accurately the model’s predictions align with the observations.