The experiment was conducted at the Jiangning Campus of Hohai University, Nanjing City, Jiangsu Province in China (31°54’57” N, 118°46’37” E). A total of 22 plots were set in this study, with a length × width of 2.5 m × 2.0 m and a depth of 2.0 m, which were cultivated with a rice–wheat rotation for many years. Rice cultivar (Nanjing-9108) was transplanted on 4 July with a spacing of 20 cm × 15 cm and harvested on 25 October 2021 under the controlled irrigation and drainage scheme. In order to obtain various spectral characteristics and photosynthetic characteristics parameters of rice canopy at different stages, five nitrogen fertilizer levels (N1–N5: 0, 150, 225, 300 and 375 kg/ha total pure nitrogen) and two infiltration levels (W1 and W2: 3 and 5 mm/day) were applied. The above two-factor complete experimental scheme was used for randomized design within the 22 plots. N fertilizers were employed as the base fertilizer (5 July), tiller initiation fertilizer (14 July), and spikelet-developing fertilizer (14 August) with the proportion of 40%, 30%, and 30% of total pure nitrogen, respectively. Phosphate (P) and potassium (K) fertilizers were applied once as the base fertilizer. All plots were well managed with practices commonly adopted by local farmers. The basic properties of the test soils are listed in Table 1 (Chen et al., 2022a) and the location of the experimental area and the arrangement of experimental treatments are shown in Figure 1 .

A DJI Innovation’s Phantom4-M (P4M) was employed as the phenotyping platform in this study. It is equipped with a multispectral camera with six CMOS, including one color sensor for visible light imaging (RGB) and five monochrome sensors for multispectral imaging. Each sensor has an effective pixel of 2.08 million, a lens field angle of 62.7°, and a focal length of 5.74 mm. Specific parameters of the sensor are shown in Table 2 . The UAV-based multispectral image data were obtained under clear and cloudless weather conditions (10:00–14:00) at each rice growth stage. The UAV flew at an altitude of 15 m, with a heading overlap of 85% and a sideway overlap of 75%.

The multispectral original images of five bands acquired by each UAV flight sortie were exported into the PIE-UAV software (Piesat Information Technology Co., Ltd., China) to correct and splice into field orthophoto. The production steps of the orthophoto were as follows: (1) image matching: match the original images with 40,000 key and tie point limits by geographical location matching method; (2) image aligning: import ground control point (GCP) information and align the images with high adjustment accuracy, 0.05 pixel GCP measurement accuracy, and 0.5 pixel connection point matching accuracy; (3) DEM building: generate DEM data using a resolution of 1 GSD; (4) tessellation building: generate the tessellation line based on the Voronoi Geometry method; (5) orthophoto correction: correct the orthophoto with automatically calculated image resolution and the mosaic line mask method; (6) color balancing: set the number of pyramid layers to 3 for color homogenization of mosaic images; and (7) image mosaicking: resample the orthophoto by the cubic convolution method and export into Geo-Tiff format. The final size and resolution of the orthophoto were 9,012 × 5,126 pixels and 7.25 mm/pixel, respectively. ENVI 5.3 was used to perform layer stacking on Tiff images of each band to obtain five-band multispectral images, and the digital numbers (DNs) were transformed into reflectance by radiometric correction.

Simultaneous field measurements were conducted within the same day of the UAV multispectral image data acquisition, including rice PH, SPAD, and photosynthetic parameters Pn, Tr, and Gs. The PH values were measured with a soft ruler from the soil ground to the leaf tip (cm). The SPAD values were measured by the chlorophyll meter model (SPAD-502, Spectrum Technologies, Inc., NE, USA) and averaged from the measurements at the tip, middle, and base of each leaf. The photosynthetic parameters were measured by the portable photosynthesis system (LI-6800, LI-COR Inc., NE, USA) at 10:00–11:30 a.m. The measured leaf position was the middle of the latest fully unfolded leaf at the jointing–booting stage, and the middle of the panicle leaf at heading–flowering and ripening stages. Each parameter was averaged from three representative plants within a 30 cm × 30 cm quadrat, and three quadrats were measured for each plot. Thus, 60 groups of field data were obtained for each growth stage and 180 groups in the total growth season. The details of ground measurements and UAV flights are listed in Table 3 .

It should be noted that the weather conditions in early September 2021 were mainly cloudy and rainy; the measurement on 3 September was the only relatively ideal condition. Therefore, the ground SPAD and Pn measurement at the heading–flowering stage and photography of UAV were affected to a certain extent.

Pearson’s correlation coefficients (r values) between the above 25 VIs and Pn at three growth stages are listed in Table 6 . The code for r value calculation and significance analysis was written in Python3.2 with the scipy package (1.20.3). Generally, the VIs had a better relationship with Pn at the jointing–booting stage, but the r value became worse as the growth stage advanced; however, it might be that the number of samples was relatively small and therefore no VIs passed the highly significant correlation test (p< 0.01). To be specific, CI_green, CVI, MNVI, NLI, OSAVI, RDVI, and RVI2 achieved a significant positive correlation (p< 0.01) with Pn at the jointing-boosting stage, with r value ranging from 0.3330 to 0.3893. NIRv, DVI, EVI, MSAVI, MSR, NDVI, RVI1, SIPI, and TCARI also had a satisfactory r value (p< 0.05) with absolute value between 0.2753 and 0.3262. At the heading–flowering stage, only CI_green, RVI2, and SIPI showed a significant relationship and NDVI, NLI, TCARI, and VDVI had a relatively higher r value. When the crop proceeded to the ripening stage, no VIs could achieve a satisfactory r value with Pn. The VIs employed for Pn estimation were thus selected based on the r value at different stages, and the selected VIs are shown in bold in Table 6 .

The NDTI, DTI, and RDTI were calculated using any two texture features from NO. 1 to NO. 40 and the correlation thermal map in Figure 4 was thus drawn according to the r value between TIs and Pn. Because the order of the two selected texture features was different, the correlation r value in the figure presented positive and negative axis symmetry. Taking the thermal map at the jointing–booting stage as an example, it could be found that there were obviously deeper red and blue lines in the DTI and RDTI figure, indicating that the DTI and RDTI composed of any of the NO. 1 (MEAN1) and NO. 25 (MEAN4) features with the other feature had a better relationship with Pn (r value mainly approximately 0.35 and 0.39, respectively). Scattered hot spots with high r value could be seen in the RDTI figure, but no dominant texture feature could be found. For each TI, the feature combination with the highest correlation with Pn was selected and the results at different growth stages are listed in Table 7 .

The accuracy comparison results between LR, SVR, GBDT, RF, and MLP models based on the selected VIs are listed in Table 8 and all criteria indices were calculated by the average of 10-fold cross-validation results. At the jointing–booting stage, most VIs showed good correlation with Pn and a total of 16 VIs were selected for modeling; therefore, the models achieved the relatively highest accuracy compared with those at other growth stages. Specifically, GBDT achieved the highest average accuracy (with an MSE of 0.253 μmol m⁻² s⁻¹, an MAE of 0.414 μmol m⁻² s⁻¹, an EVS of 0.938, and an R ² of 0.938), while SVR models attained the worst performance (with an MSE of 2.512 μmol m⁻² s⁻¹, an MAE of 1.155 μmol m⁻² s⁻¹, an EVS of 0.390, and an R ² of 0.383). RF, MLP, and LR models ranked second, third, and fourth, respectively. Figure 5A also shows that the Pn estimated value of GBDT was the closest to the measured Pn value in each validation set. The Pn estimated values of LR and SVR models were almost concentrated in the range of 24–27 μmol m⁻² s⁻¹; thus, the estimation accuracy was not satisfactorily compared with GBDT, RF, and MLP models, where the Pn values were relatively lower and higher (inside the blue and red dashed circle). The correlation between VIs and Pn gradually weakened and inputted VIs thus decreased in number as the growth stage progressed; therefore, the model estimation accuracy decreased to a certain extent without ranking change. In detail, the LR model suffered the biggest loss in estimation accuracy with an R ² value decreasing to 0.296 at the heading–flowering stage and then to 0.125 at the ripening stage. However, GBDT and RF models still showed good performance with R ² values of 0.928 and 0.869 at the heading–flowering stage and 0.863 and 0.815 at the ripening stage, respectively. Figures 5B, C also demonstrated that the GBDT and RF models could better describe the relationship between VIs and Pn in the value full range at different growth stages. Although the performance of the MLP model ranked third, the performance was not ideal when compared with GBDT and RF models for both lower and higher PN value estimation. In conclusion, the accuracy of modes for estimating Pn value at the jointing–booting stage was relatively best and the GBDT model could be highly recommended for Pn estimation during the rice whole growth season.

In order to further improve the Pn estimation accuracy, TIs (NDTI, DTI, and RDTI) and basal growth index (PH and SPAD) were introduced based on the input VIs. The accuracy comparison result of different models under different input combinations is shown in Table 9 . After adding the TIs inputs for Pn estimation, all models performed higher accuracy at the jointing–booting stage with the MSE of LR, SVR, GBDT, RF, and MLP decreasing by 0.090, 0.370, 0.126, 0.113, and 0.057 μmol m⁻² s⁻¹, respectively, and the R ² increasing by 0.022, 0.090, 0.031, 0.028, and 0.015, respectively. As the basal growth index further increased, the final R ² of the employed models increased to 0.792, 0.565, 0.987, 0.943, and 0.822, respectively. It is possible that the basal growth index improved the model accuracy slightly more than TIs because the difference in PH could reflect the stress of crops to a certain extent, and SPAD is directly related to chlorophyll content, which directly affects photosynthesis. The same improvement effect of model accuracy could also be found at the heading–flowering and ripening stages, and although the improvement of GBDT and RF models was relatively small (with R ² increasing to 0.062 and 0.031 for VIs + TIs input and 0.113 and 0.132 for VIS + TIs + PH and SPAD input), they were still the top two models among the employed models. The greatest improvement could be found in the LR model, with an R ² increase of 0.139 and 0.470 at the heading–flowering stage and 0.114 and 0.471 at the ripening stage under the VIs + TIs and VIS + TIs + PH and SPAD input combination, respectively. The accuracy of the SVR model also had been improved greatly, but it still is the lowest among the five models. In conclusion, both TIs and the basal growth index could obviously improve the model accuracy for Pn estimation and the PH and SPAD had a better effect compared with the TIs in this study, which significantly improved the model performance at the heading–flowering and ripening stages, especially for LR and MLP models.

The ground sample results showed that the SPAD and Pn generally increased with nitrogen application (N1–N3 levels), then decreased at the N4 level, and finally recovered at the N5 level under the same leakage treatment at the jointing–booting stage, which also indicated that proper N application could improve the photosynthesis, while excessive N application not only had a poor effect on photosynthesis, but also affected plant growth and increased risk of contamination during leakage and drainage. This phenomenon was similar to that found by Cisse et al. (2020). It might due to the fact that there was no significant difference in Rubisco activity and non-photochemical quenching (NPQ) between the N4–N5 and the N3 level; thus, excessive energy could not be dissipated by NPQ, leading to oxidative stress, resulting in a decrease in Pn when excessive nitrogen was applied. An opposite trend could be found for canopy multispectral reflectance as the N application increased. Band 1 to band 4 generally decreased when N application increased from the N1–N3 level, while they slightly increased at the N4 level then decreased again at the N5 level; however, band 5 reflectance consistently increased with N level, which was consistent with previous studies on other crops Qiu et al. (2015). Generally, crops with good growth have a lower reflectance and a steeper increase from the red to the NIR band; thus, the variation trend of the canopy reflectance was consistent with SPAD and PN, which also provides a theoretical basis for the inversion of photosynthetic characteristic parameters using VIs.

Based on the result and analysis in Section 3.1.1, the relationship between VIs and Pn decreased with the growth stage, which might be due to the influence of heading and flowering, although panicle photosynthesis is also an important part of crop canopy photosynthesis and contributes significantly to grain formation. However, due to the limitation that the photosynthetic measurement equipment used in this paper could only be used to measure leaves, the canopy photosynthesis was thus approximately the photosynthetic capacity of the included leaves. Therefore, after rice heading and flowering, the spectral reflectance of the canopy was affected to a certain extent and the correlation based on the above data decreased significantly. In order to improve the estimation accuracy of Pn at the heading–flowering stage, the image segmentation should be carried out first to remove the panicle reflectance image. Meanwhile, the method of canopy photosynthesis measurement and the effect of panicle on photosynthetic contribution and reflectance should be revised and improved in future studies.

According to the correlation analysis results of TIs during the whole growth season in Figure 4 , it could be concluded that most of the dominant texture features were extracted from band 1 to band 3, while the TIs constructed by the features extracted from band 4 and band 5 were not that satisfactory. Specifically, MEAN and COR texture features were more included in the optimal features for NDTI, DTI, and RDTI construction, because MEAN represents the average of moving windows containing targets and backgrounds, which can smooth the image and reduce the interference of background factors, while COR can reflect the local gray-level correlation in the image and distinguish the differences of image texture in all directions. In addition, window sizes of 6 × 6, 9 × 9, and 12 × 12 were also employed in the GLCM for textural feature extraction; however, the performance was not much different from that of the 3 × 3 window size, which might be due to the high resolution (4.5 mm/pixel) of each pixel in the image. Therefore, improving the resolution of visible images to extract higher-quality TIs might be an economical and practical approach to improve the accuracy of Pn estimation.

The model performance under the different input combinations of VIs, VIs + TIs, and VIs + TIs + PH and SPAD concludes that the fusion of VIs and TIs could effectively improve the accuracy of Pn estimation because the VIs contain the canopy reflectance features and are more sensitive to the nutrition variations, while the TIs could better reflect the difference in canopy structure. The results were also consistent with previous studies (Liu et al., 2019; Zhang et al., 2021) in that this fusion combination could improve the estimation accuracy of biomass, LAI, nitrogen nutrition, and potassium nutrition and accumulation. Furthermore, the addition of PH and SPAD brought higher accuracy, which could attribute its success to the high correlation between PH and N nutrition status, SPAD and chlorophyll content, brought by the obvious difference of 5 N treatment. To sum up, more different kinds of data could introduce more direct or indirect related features, and it is also suggested that stacking and blending ensemble learning methods (Wu et al., 2021) could be employed to combine the model ability of feature extraction and analysis based on different principles in future research to improve the model accuracy for Pn estimation, which is also the purpose and significance of developing agricultural big data and agricultural intelligent models.