The research was carried out in the Shihezi, China (86.03°E, 44.18°N, a.s.l. 429 m) ( Figure 1 ). This region has a typical temperate continental climate. The average annual precipitation was only 125.9–207.7 mm, and the annual accumulated temperature (≥10°C) was 3,570°C–3,729°C. In addition, there was a large day–night temperature difference. The soil type was loam. Soil organic matter, alkaline hydrolyzable nitrogen, available phosphorus, and available potassium content were 20.10 g/kg, 60.92 mg/kg, 17.83 mg/kg, and 142 mg/kg, respectively. The previous crop was cotton.

In late April 2021, cotton seeds (variety Xinluzao 45) were sown, and drip irrigation and plastic film mulching were adopted [six rows and three irrigation tapes were under the mulching of one film (2.35 m in width)].

The randomized complete block design was adopted in this experiment, including five N treatments: 0 kg·ha⁻¹ [N0 (CK)], 110 kg·ha⁻¹ (N1), 220 kg·ha⁻¹ (N2), 330 kg·ha⁻¹ (N3, conventional application rate), and 440 kg·ha⁻¹ (N4). Each treatment had three replicates (plots). Each plot was 70.5 m² (7.05 m × 10 m), and the plot spacing was 4.7 m. Urea (N, 46%) of 345 kg·ha⁻¹ and potassium dihydrogen phosphate (P₂O₅, 52%; K₂O, 34%) of 240 kg·ha⁻¹ were applied. The first fertilization was conducted on June 12, and fertilizers were applied every 8–9 days thereafter (nine times totally).

The canopy spectrum of cotton plants without pests and diseases were collected at 12:00–14:00 on sunny days at the budding stage (June 21), flowering stage (July 10), initial bolling stage (July 25), and full bolling stage (August 18) using an ASD-FieldSpec Pro FRTM spectrometer. During acquisition, the probe of the sensor was downward, the vertical height was 50 cm, the sensor angle was 25°, and the scanning time was 0.2 s. Three sampling points were selected in each plot, and spectral acquisition was performed five times for each point. Finally, the average was calculated, which was the spectrum of the sample point. Spectral acquisition were conducted eight times totally during the whole growth period, and a total of 360 data were obtained. Previous studies have shown that the spectrum at 1,800 nm–2,500 nm is greatly affected by soil background and air moisture (Chandrasekar et al., 2022). Therefore, in this study, only spectral changes at 350 nm–1,800 nm were analyzed.

After spectral acquisition, three cotton plants were selected from each plot. Then, leaves were separated, placed in a crisper box, and sent to the lab for LAI determination using a LI-3100 leaf area meter. Eight samplings were conducted during the whole growth period, and a total of 360 cotton plant samples were collected. The calculation method of LAI is shown in Equation 1:

where K0 is the representative of cotton leaf area in a plot, N is the cotton plant number of one plot, and K is plot area.

The VIs have been widely used to monitor vegetation information. Especially, VIs can distinguish non-vegetation information such as water bodies and enhance vegetation information ((Donovan et al., 2021; (Morsy et al., 2022). In this study, 20 VIs ( Table 1 ) were selected for correlation analysis.

To explore the interaction between VIs and the correlations between VIs and LAI, a correlation matrix between VIs and LAI was generated using Origin. Pearson correlation coefficient indicates the level of correlation (−1–1). The higher the absolute value, the closer the correlation (Bermudez-Edo et al., 2018).

In this study, the entire dataset (360) was randomly divided into a modeling set (300) and a validation set (60) (5: 1) to ensure the generality and stability of the model (Joseph, 2022). The results of statistical analysis of the entire dataset, the modeling set, and the validation set ( Table 2 ) showed that the modeling set had a larger data range than the validation set and that the mean, minimum, maximum, standard deviation, and coefficient of variation of the modeling and validation sets were very similar. This indicates that the division of modeling and validation sets are very uniform, which is conducive to the accurate evaluation of model performance (Wu, 2023).

Multiple VIs with high correlation with LAI were selected to construct ELM, RF, BP, MLR, and SVM models for LAI estimation using the Scikit-Learn package in Python. Then, the accuracy of the constructed models were compared, and the cross-validation was performed to select the optimal modeling method. Based on the optimal model, the 20 VIs ( Table 1 ) were randomly arranged and fused, and the optimal fusion was selected. Figure 2 shows the specific workflow. In this study, the change trend of the spectral reflectance of cotton canopy and LAI during the whole growth period were first analyzed, and then, 20 VIs related to LAI were selected based on previous studies. After that, the correlations between LAI and VIs were analyzed, and the VIs with a correlation coefficient with LAI of no less than 0.6 (Meng et al., 2021) and a correlation coefficient with other VIs of less than 0.9 were selected (Pasternak and Pawluszek-Filipiak, 2022). Then, the LAI estimation models based on machine learning (ELM, RF, BP, MLR, and SVM) were constructed and compared to select the optimal model. Finally, the VIs participating in the modeling were randomly arranged and fused (i.e., for index 1, 2, and 3, there were seven fusions including 1, 2, 3, 1 + 2, 1 + 3, 2 + 3, and 1 + 2 + 3), and the model was verified to select the best vegetation index fusion model.

In this study, the R², RMSE, nRMSE, and MAE were used to assess the model accuracy. The larger the R², the better the model fit; the smaller the RMSE, nRMSE, and MAE, the higher the model accuracy. The calculation of R², RMSE, nRMSE, and MAE was based on the research of Wang et al. (2022) and Jiang et al. (2023).

Cotton LAI increased first and then reduced ( Figure 3 ), and the influence of different N treatments on cotton LAI was great throughout the entire growth period. The more N was applied, the greater was the LAI. The dynamics of cotton LAI in different treatments were similar throughout the entire growth period. The changes of LAI in different treatments were small at the beginning. Then, the LAI increased, and the difference in LAI between different treatments also increased. The LAI peaked in the late flowering stage and then declined.

The spectral reflectance changed consistently in different treatments ( Figure 4 ). In the visible region, there was little difference. With the growth of N dose, the spectral reflectance declined (N0 > N1 > N2 > N3 > N4). In addition, the reflectance varied obviously among stages in the NIR region.

In the bud stage, the reflectance of the N0 and N1 treatments were low (N0< N1) and that of the N2, N3, and N4 treatments were high (N4 > N3 > N2) in the NIR region. In addition, there was a significant difference in the reflectance between different treatments throughout the entire growth period. The reflectance increased with the increase in N dose (N4 > N3 > N2 > N1 > N0).

The coefficients of correlation between VIs and LAI were 0.14–0.67 ( Figure 5 ). LAI had a negative correlation with MTVI and RES and a positive correlation with other VIs. Among them, MNDSI had the highest correlation with LAI (0.67), and TVI had the lowest correlation with LAI (0.14). The correlation between VIs varied greatly. The correlation between SAVI and MTVI was the lowest (0.033), while that between EVI and RDVI and between RES and REP were the highest (0.98). In addition, RES had a positive correlation with MTVI, but a negative correlation with other VIs. MTVI had a positive correlation with RDVI, SAVI, TVI, and EVA, and a negative correlation with other VIs.

The 20 VIs were used to construct cotton LAI prediction models based on a single VI ( Table 3 ), and it was found that the accuracy of the models were generally low. To further explore the advantages of multi-VIs fusion and the performance of the constructed models, considering the multicollinearity between VIs and the overfitting of the model, this study selected the VIs with a correlation coefficient with LAI of not less than 0.6 and a correlation coefficient with other VIs of less than 0.9 (results of Section 3.3) for modeling ( Table 4 ), including: GNDVI, MSR, MNDSI, RVI, SRI, GRVI, NVI, REP, and CIred-edge. The comparison of the ELM, RF, BP, LR, and SVM modeling showed that the RF models had the highest accuracy, with R², MAE, RMSE, and NRMSE being 0.87, 0.37, 0.57, and 0.16, respectively, followed by the ELM, BP, and LR models. SVM models had the lowest accuracy, with R², MAE, RMSE, and NRMSE being 0.51, 0.84, 1.11, and 0.32, respectively.

The selected nine VIs were randomly arranged and fused. When there were 1, 2, 3, 4, 5, 6, 7, 8, and 9 VIs for fusion, 1, 15, 15, 15, 15, 15, 15, 9, and 1 combinations were randomly selected, respectively. Only the optimal fusions under different number of VIs are shown here. It was found that the RF model had the highest accuracy when MNDSI, SRI, GRVI, REP, CIred-edge, MSR, and NVI were fused ( Table 5 ). The RF-7 model had the highest R² (0.90) and the smallest RMSE (0.50), NRMSE (0.14), and MAE (0.26). The accuracy of the models built by a single VI was the lowest, with R², MAE, RMSE, and NRMSE being 0.63, 0.82, 1.03, and 0.30, respectively.

The results of linear fitting of predicted LAI by RF models ( Table 5 ) and measured LAI ( Figure 6 ) showed that the slope of the fitted line of measured LAI and RF-7 was closest to 1, and the slope of the fitted lines of measured LAI and single VIs was the lowest ( Figure 6A ). When 1, 2, 3, 4, 5, 6, and 7 VIs were fused ( Figures 6A–G ), the model accuracy increased with the increase in the number of VIs fused, but the accuracy decreased when 8 and 9 VIs were fused ( Figures 6H, I ). The slopes of all curves were less than 1, indicating that there was no over-fitting.

Leaf area index reflects the energy flow and material cycle in plant canopy and is closely related to crop yield. Therefore, it is commonly used in the study of photosynthesis, respiration, and carbon accumulation of vegetation (Lu et al., 2020; Sun et al., 2023). In this study, the cotton LAI under different N treatments showed an unimodal variation trend. This reflects that N has an obvious regulatory effect on cotton LAI, which is related to the growth characteristics of cotton at different stages (Ninkuu et al., 2023). From bud to flowering stage, due to the increase in cotton branches and leaves, the cotton leaf area continues to increase, resulting in a rapid increase in LAI, and the difference in LAI between N treatments also increases. This is similar to the results of Wen (2016). Nitrogen is the main nutrient required for cotton growth, and sufficient N supply could promote leaf enlargement. The bud stage is a key period for cotton vegetative growth. With the increase in temperature, cotton growth accelerates. This ultimately leads to increased LAI. During the flowering stage, cotton enters into the reproductive growth stage, and a large amount of organic matter is transferred to reproductive organs. In the late flowering stage, LAI is the largest, and then, LAI gradually declines in the bolling stage. However, Kumar et al. (2022) found that LAI began to decrease after the boll-opening stage. On the one hand, it may be caused by the different temperature, sunshine time, and rainfall in different regions (high temperature and long sunshine time can accelerate cotton growth, but rainfall slows down cotton growth). On the other hand, the increasing competition between vegetative and reproductive growth in the flowering and bud stages leads to leaf shedding and decreased LAI (Li, 2022). Especially, in this study, the LAI declined fastest in the N4 treatment. High N dose leads to too vigorous vegetative growth. The lower leaves block each other, leading to poor ventilation and early leaf aging and shedding. In addition, in our previous report (Fan et al., 2023), the LAI of cotton changed most significantly under the N3 treatment, while the LAI changed the most under N4 treatment in this study. This may be due to the fact that the N4 in this study was 440 kg N ha⁻¹. However, in previous study, the N3 was 360 kg N ha⁻¹ and the N4 was 480 kg N ha⁻¹. Although the N application rate N4 is the largest in this study, it is lower than that in previous study. This indicates that the demand for N fertilizer in cotton is limited, and it is not that the more nitrogen is fertilized, the greater the LAI.

The uptake, transport, and assimilation of N show dynamic changes in cotton leaves, and the nutrient composition and content in leaves also change continuously, resulting in changes in leaf color, morphology, and spectral reflectance (Sun et al., 2022). In the visible region, there was no significant difference in canopy spectral reflectance at different stages, and the reflectance decreased with the growth of N dose. In the NIR region, the difference increased, and the reflectance increased with the growth of N dose. This is consistent with the results of Liu et al. (2019). However, in our previous study (Fan et al., 2023), the canopy spectral reflectance reached the maximum under N3 treatment. In this study, the canopy spectral reflectance was the largest under N4 treatment. This may be due to the fact that although the N4 in this study is the largest, it is still lower than the N4 in the previous study. Excessive or small nitrogen fertilizer application can lead to a decrease in LAI, thereby reducing the canopy spectral reflectance.

The spectral reflectance of visible region mainly reflects leaf pigment status. In crop photosynthesis, chlorophyll reflects green light and absorbs red and violet light, resulting in low canopy spectral reflectance (Li et al., 2019; Yang et al., 2022). The reflectance in the NIR region is mainly affected by optical properties and leaf structure. In addition, the differences in intercellular space, composition, and shape are also influencing factors (Richetti et al., 2019; Spafford et al., 2021). Thus, the reflectance of the NIR region gradually reduces over time.

Vegetation index can provide crop growth information. At present, a large number of remote sensing studies of cotton LAI based on a single vegetation index have been carried out (Yan et al., 2022). In this study, the correlation analysis of LAI and VIs was carried out. It was found that MNDSI had the highest correlation with LAI. This is consistent with the results of Gao et al. (2013). MNDSI can better reflect the status of cotton LAI. The spectral bands in MNDSI has a strong absorption on the spectra of cotton LAI (Li et al., 2018). Therefore, future studies can try to accurately predict cotton LAI changes through this index. In this study, the correlations between some VIs were very high, with a correlation coefficient of 0.98 (EVI and RDVI, RES and REP). This may be due to the fact that these VIs contain the same spectral bands, resulting in very high multicollinearity between the VIs, so it is necessary to filter the VIs before fusion.

It was found that the ML models based on 20 VIs had a low accuracy. This may be due to the high correlations between some VIs. Therefore, the VIs with a correlation coefficient not less than 0.6 were selected for modeling. It was found that this improved the prediction accuracy of all models (ELM, RF, BP, MLR, and SVM). This indicates that LAI has its unique spectral features and corresponding VIs. Among the models, the RF models had the highest accuracy. This is consistent with the results of Han et al. (2019) and Wang et al. (2021). RF is a tree-based ensemble learning. It higher accuracy in this study may be due to that the data are trained using randomly selected subsets at each node, and the best performing prediction variables are selected to split nodes (insensitive to noise), which can effectively solve problems such as overfitting and collinearity (Lu and He, 2019).

In this study, multiple VIs were randomly arranged and fused, and the optimal fusion was selected for RF modeling. It was found that the model accuracy constructed using VI fusion was higher than that of the model constructed using a single VI. This is similar to the results of Qi et al. (2020) and Rumora et al. (2021). This may be due to that each VI is a physical structure and has limited valid bands (only two to three bands). It was also found that the model accuracy based on the fusion of seven VIs (RF-7) had the highest accuracy. This indicates that the fusion of the seven VIs synthesizes the structural characteristics of each VI and contains more reflectance information, which ultimately improves the LAI prediction accuracy. However, when eight and nine VIs were fused, the model accuracy began to decrease. This may be due to that (1) different VIs have different sensitivities to the same parameter; the model accuracy increases when the highly sensitive VIs are fused for modeling, and vice versa. (2) The influences of canopy structure, cover, and soil background brightness on different VIs are also different (Qiao et al., 2022). (3) There is data redundancy between VIs, which affects the model accuracy. Therefore, appropriate VI fusion can eliminate the influence of environmental factors, realize the complementarity of spectral features, and increase vegetation information, which ultimately improves model accuracy and stability.

In addition, hyperspectral data are more informative than VIs. It is widely used in different remote sensing analysis and feature classification tasks because hyperspectral data covers a wide range of bands and contains all the spectral information of features. However, it faces the problem of information redundancy, which increases the complexity of data processing and analysis, and may reduce the interpretability of the model (Jiang et al., 2022). In addition, sensitive wavelengths are more susceptible to interference from non-target factors, such as soil background and light conditions (Meiyan et al., 2023). In contrast, the VI is usually based on the combination of specific wavelengths, which can effectively eliminate the interference of soil signals and other factors, and is specially used for the monitoring of vegetation parameters. It is well targeted and has a high interpretability for vegetation-related research (Cao et al., 2021).

At present, the vegetation canopy spectral acquisition instruments mainly include hyperspectral sensors and multispectral (RGB) sensors. The wide spectral range of hyperspectra, including visible light and near-infrared spectra, provides rich spectral information, which enables hyperspectral data to capture detailed spectral characteristics of features (Liu et al., 2023a), such as spectral reflectance, absorption peaks, and absorption valleys. These detailed spectral features make hyperspectral data have higher performance and application potential in object detection. However, the number of RGB bands is small and contains less information (Liu et al., 2022a), making it difficult to fully reflect the changes in the spectral features of vegetation parameters. Nevertheless, due to the processing complexity and high cost of hyperspectral data, RGB still has certain advantages in some scenarios, especially in vegetation monitoring and simple object classification tasks. For example, Liu et al. (2022b) improved the estimation accuracy of potato aboveground biomass (AGB) by extracting different texture features from RGB images.

In recent years, many scholars have preprocessed hyperspectral data and achieved good estimation results. For example, Liu et al. (2022c) used the Savitzky-Golay (SG) smoothing to smooth the spectrum collected by unmanned aerial vehicles (UAVs) and the random forest to extract the characteristic bands, and found that the PLSR model had the highest estimation accuracy for AGB. Liu et al. (2023b) extracted spectral features from UAV-collected RGB images and hyperspectral data through wavelet transform and found that the fusion of the features of RGB images and hyperspectral data could improve the estimation accuracy of AGB. Since hyperspectra can provide comprehensive spectral information, future research can transform the raw hyperspectral data to highlight spectral features, extract features by some methods to reduce data redundancy, and calculate the vegetation index based on the spectral features to construct a model. In addition, RGB images and hyperspectral data can be combined to explore the potential of multi-source remote sensing data fusion in predicting cotton LAI.

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.