Retrieving rice (Oryza sativa L.) net photosynthetic rate from UAV multispectral images based on machine learning methods

Photosynthesis is the key physiological activity in the process of crop growth and plays an irreplaceable role in carbon assimilation and yield formation. This study extracted rice (Oryza sativa L.) canopy reflectance based on the UAV multispectral images and analyzed the correlation between 25 vegetation indices (VIs), three textural indices (TIs), and net photosynthetic rate (Pn) at different growth stages. Linear regression (LR), support vector regression (SVR), gradient boosting decision tree (GBDT), random forest (RF), and multilayer perceptron neural network (MLP) models were employed for Pn estimation, and the modeling accuracy was compared under the input condition of VIs, VIs combined with TIs, and fusion of VIs and TIs with plant height (PH) and SPAD. The results showed that VIs and TIs generally had the relatively best correlation with Pn at the jointing–booting stage and the number of VIs with significant correlation (p< 0.05) was the largest. Therefore, the employed models could achieve the highest overall accuracy [coefficient of determination (R 2) of 0.383–0.938]. However, as the growth stage progressed, the correlation gradually weakened and resulted in accuracy decrease (R 2 of 0.258–0.928 and 0.125–0.863 at the heading–flowering and ripening stages, respectively). Among the tested models, GBDT and RF models could attain the best performance based on only VIs input (with R 2 ranging from 0.863 to 0.938 and from 0.815 to 0.872, respectively). Furthermore, the fusion input of VIs, TIs with PH, and SPAD could more effectively improve the model accuracy (R 2 increased by 0.049–0.249, 0.063–0.470, and 0.113–0.471, respectively, for three growth stages) compared with the input combination of VIs and TIs (R 2 increased by 0.015–0.090, 0.001–0.139, and 0.023–0.114). Therefore, the GBDT and RF model with fused input could be highly recommended for rice Pn estimation and the methods could also provide reference for Pn monitoring and further yield prediction at field scale.


Introduction
Photosynthesis is one of the most crucial parts of the global carbon and energy cycle (Reichstein et al., 2013;A Ivlev, 2017). The crop photosynthesis activities assimilate carbon dioxide (CO 2 ) and water (H 2 O) by using light energy to form organic matter and, therefore, are a key determinant of food production and security (Reich and Amundson, 1985;Long et al., 2006). Net photosynthetic rate (Pn) is the value of the total photosynthetic rate minus the respiration rate, which directly refers to the organic matter accumulated. Although researchers have gradually deepened the understanding of photosynthesis based on cell-scale gas exchange, current methods and equipment developed based on these theories are still mainly focused on the leaf level, which is time-consuming and has a poor representation (Stinziano et al., 2019). It is scientific to use large canopy photosynthesis and transpiration measurement system (CAPTS) (Song et al., 2016) to observe photosynthesis at the canopy scale, but the investment is too expensive to be popularized in regional-scale monitoring.
The mobile high-throughput phenotyping platforms (HTPPs) (Deery et al., 2014;Li et al., 2014) with RGB, fluorescence, hyperspectral, thermal, 3D laser, and computed tomography (CT) imaging sensors provide a non-destructive method for rapid crop phenotypic acquisition. In particular, a high-spectral-resolution spectroradiometer (Aguirre-Gomez et al., 2001;Meroni and Colombo, 2006) (most Fieldspec 4 or 4pro, Analytical Spectral Devices, ASD, Boulder, CO, USA) is the most physical and effective equipment for photosynthesis monitoring on the ground. The sensitive band reflectance or vegetation indices (VIs), generally including 2 or more band reflectance, was commonly used to establish a linear or nonlinear relationship with crop physiological and biochemical parameters. Qiu et al. (2015) comprehensively analyzed the correlation between main photosynthetic, fluorescence parameters and hyperspectral data in ear position leaves of maize and found that Dl699 had the best correlation with Pn. Sun et al. (2016) introduced wavelet analysis (WA) to select the sensitive bands of hyperspectral for estimating Pn of winter wheat on the leaf scale and found that the models based on WA were more accurate than the VIs method. Fu et al. (2019) constructed a stacking framework for retrieving the maximum carboxylation rate of Rubisco (V c,max ) and the maximum electron transport rate supporting RuBP regeneration (J max ) in the photosynthesis parameters of tobacco based on canopy hyperspectral reflectance, which further improve model accuracy compared with the basic models. Based on the advantages of ground platform on high-resolution continuous spectrum and texture features, the above research could provide practical and accurate estimation of photosynthetic parameters. However, the photosynthetic monitoring in the actual production field could hardly be represented due to environmental factors and the use of various equipment requires expertise.
As a new near-ground remote sensing approach, unmanned aerial vehicles (UAVs) (de Castro et al., 2021) can flexibly provide higherresolution and bigger-scale images by carrying different sensors (e.g., multispectral, hyperspectral, and thermal infrared cameras). It has already been widely used in the inversion of physiological and biochemical parameters such as plant height (PH) (Che et al., 2020), leaf area index (LAI) (Chen et al., 2022b), nutrient states (Xu et al., 2021), and aboveground biomass . Equipped with hyperspectral imaging (HSI) sensors, Liu and Peng. (2020) employed eight chlorophyll-related VIs for estimating maximum Pn, and proposed a model based on chlorophyll index (CI) and photosynthetically active radiation (PAR) for different rice varieties. However, the water vapor in the field (especially the paddy field) might have a great influence on the hyperspectral data. Additionally, high price and tedious data processing process (e.g., noise processing, dimension reduction, and spectral unmixing) have prevented the commercial application of this method. Therefore, multispectral sensors that can characterize key points (usually including blue, green, red, red edge, and near-infrared bands) in crop canopy spectra features are more commonly used for practical application on UAV. Chen et al. (2018) established linear inversion models of photosynthetic parameters at different time points in the cotton bud stage based on UAV six-band multispectral images. However, crossgrowth stage comparison and regression model selection could be done more comprehensively. Based on the relationship between VIs constructed from UAV multispectral image and photosynthetic parameters, Zhang et al. (2020a) explored the inversion method of diurnal variation of photosynthesis in rice canopy combined with the light response curve model and provided a method with physical basis for gross primary productivity (GPP) inversion, while the scale effect between 100-m UAV multispectral data and PAR monitoring data from a single point on the ground should be further discussed. On the other hand, the image obtained by UAV remote sensing has a higher resolution than satellite remote sensing; thus, it has more detailed texture features that can better reflect the difference in the set window size. Therefore, textural indices (TIs) are commonly introduced with VIs to improve the model accuracy. According to previous studies, TIs have a good correlation with aboveground biomass (Sarker and Nichol, 2011;Liu et al., 2019) and thus also have a good relationship with the accumulated amount of canopy elements (Pimstein et al., 2011;Lu et al., 2019;Zhang et al., 2021) (e.g., nitrogen, potassium, and chlorophyll). Zheng et al. (2019) have found that the normalized difference texture index (NDTI) is in good relationship with rice biomass and the fusion of NDTI with VIs improved the accuracy of biomass estimation. Similarly, Lu et al. (2021) and Zheng et al. (2020) demonstrated that the fusion of TIs and multispectral VIs could effectively improve the estimation of potassium accumulation and nitrogen accumulation in rice. Since the accumulated organic matter of photosynthesis can directly affect the basic growth indicators of rice such as plant height, tiller number, and leaf area index, TIs could also have untapped potential in Pn estimation.
For models employed in the inversion studies, linear regression or nonlinear regression were commonly used to construct inverse functions with definite expressions, but the accuracy is relatively low and poor in portability (Wan et al., 2021). Machine learning methods have been widely used in the regression and classification issues and have been proven to be fast, accurate, and good at generalization. WA (Bruce et al., 2002), partial least square regression (PLSR) , and least absolute shrinkage and selection operator (LASSO) (Yang and Bao, 2017) are usually used in HSI studies to reduce the high-dimension hyperspectral data to a few important components that are sensitive to the target parameters.
Artificial neural network (ANN) (Liang et al., 2015), kernel-based support vector machine regression (SVR) (Liang et al., 2016), and random forest (RF) (Cheng et al., 2022) regression methods have been most wildly employed to explore and fit the nonlinear relationship between reflectance or VIs and inversion objects. Other machine learning methods based on ANN, kernel function, and tree also have great potential in this issue. Yang et al. (2022) built a Bayesian neural network (BNN) model to predict potential maximum quantum yield (Fv/Fm) and two other chlorophyll fluorometer parameters of grape by quantifying the HSI response indices of photosynthetic pigments and water status parameters. Yuan et al. (2022) simulated the maximum carboxylation rate at 25°C (V m25 ) of crops over time based on the convolutional neural network (CNN) model combining flux and satellite remote sensing data to further improve the estimation accuracy of GPP. Based on the leaf phenotype data, Zhang et al. (2020b) established poplar Pn estimating models using the extreme gradient boosting model (XGBoost).  have also proven the good performance of machine learning models based on the rich feature input of the HSI data for photosynthetic parameter estimation. However, there is a lack of understanding of machine learning methods for photosynthesis parameters estimating with less reflectance features based on the multispectral data and less research on rice.
In this study, multispectral images of rice canopy were acquired by UAV, and the responses of multispectral reflectance features together with Pn, PH, and SPAD to the different nitrogen or leakage treatments were analyzed at different growth stages. The correlation between Pn, VIs, and TIs extracted from the multispectral reflectance was compared, and the VIs with relatively significant correlations were employed as input of the five machine learning models. Model performance comparison under different input combinations was performed, and the improvement of fusing TIs and basal growth index PH and SPAD was further analyzed. The final purpose is to explore an economical and accurate method at field scale for the estimation of Pn and photosynthesis stress detection during the whole growth season of rice.

Study area
The experiment was conducted at the Jiangning Campus of Hohai University, Nanjing City, Jiangsu Province in China (31°54'57" N, 118°4 6'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.

UAV based multispectral data acquisition and processing
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.

Field data collection
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.

Vegetation index and textural index calculation 2.4.1 Vegetation index calculation
VI is established by the linear or nonlinear combination of different spectral band reflectances, which is a common method to retrieve physiological and biochemical indicators of crops (Zeng et al., 2022). A set of 25 commonly used VIs were employed in this study to investigate the relationship between VIs and rice photosynthetic parameters. Threshold processing was firstly performed on the stacked multispectral image to eliminate the influence of water on the reflectance. The canopy reflectance of each band within the 30 cm × 30 cm region of interest (ROI) was then averaged to calculate the VIs of each plot. The involved VIs and formulas are listed in Table 4.

Textural index calculation
Gray-level cooccurrence matrix (GLCM) (Haralick et al., 1973) was applied in this study to extract eight texture features from each band in the stacked image, including mean (MEAN), variance (VAR), homogeneity (HOM), contrast (CON), dissimilarity (DIS), entropy (ENT), second moment (SEC), and correlation (COR), and a total of 40 texture features (with a 3 × 3 pixel window size) were obtained. Study area and experiment treatments. W represents the leakage treatments (including W1: 3mm/day and W2: 5mm/day); N represents the nitrogen treatments (including N1-N5: 0, 150, 225, 300 and 375 kg/ha total pure nitrogen, respectively); GCP is abbreviation of gourd control points for geometric correction; Ground measurements in each sample point were averaged from 3 representative plants.
The same ROI size with VIs was used to extract texture features and the average value was taken. The extracted GLCM texture features were numbered in the order of MEAN, VAR, HOMO, CON, DIS, ENT, SEC, and COR of each band (in the order of band 1 to band 5) from 1 to 40. The normalized difference textural index (NDTI), difference textural index (DTI), and renormalized difference textural index (RDTI) were selected to construct TI involving two different texture features. The TI formulas are as follows: where T 1 and T 2 represented two random different texture features.

Modeling and validation 2.5.1 Machine learning regression methods
Linear regression (LR), support vector regression (SVR), gradient boosting decision tree (GBDT), random forest (RF), and multilayer perceptron neural network (MLP) were employed in this study for Pn estimation. The gridsearch tuning results for the hyperparameters of each model are listed in Table 5, where the unmentioned hyperparameters were the default values.
(1) Linear regression: LR is a traditional algorithm based on classical statistics, which is the most commonly used model in the spectral inversion research because of its simple construction form and strong interpretation. Combined with the correlation analysis, the relationship between variables and target parameters can be directly reflected. In this study, the LR model with the ordinary least squares method was used for Pn multiple regression.
(2) Support vector regression: SVR is an important application branch of support vector machine (SVM), which seeks the optimal hyperplane by minimizing the total deviation of all sample points from the hyperplane (Cortes and Vapnik, 1995). Unlike ordinary least squares, the SVR model sets a threshold ϵ around the regression line such that all data points within ϵ are not penalized for their errors.
Kernel function, gamma, and C are crucial parameters in the SVR model and have been tuned through the gridsearch method in the sklearn package.
(3) Gradient boosting decision tree: GBDT is an iterative decision tree algorithm with a "boosting" ensemble learning method (Friedman, 2001;Wu et al., 2020). The basic learners [usually classification and regression tree (CART)] in the GBDT model have strong dependencies between each other and are trained by progressive iterations based on the residuals. The results of all basic learners are added together as the final output, which grant GBDT great advantages in overfitting and computational cost fields and reduce bias at the same time.
(4) Random forest: RF is one of the most popular tree algorithms proposed by Breiman (2001) based on the bagging idea of ensemble learning. RF applies the "bootstrap" method to retrieve samples to train the N basic learners (usually CART) in parallel without dependence. The final output of the RF model is derived by combining results of the basic models with the "voting" method, which makes the RF model insensitive to outlier variable values.
(5) Multilayer perceptron neural network: MLP is generally composed of a fully connected input layer, a hidden layer, and an output layer, in which the hidden layer can be multiple (Khoshhal and Mokarram, 2012). As the most basic form of feed-forward neural network, the MLP model has been widely applied in the analysis of various complex problems and is also the foundation of CNN, deep neural network (DNN), and other complex neural networks. A typical three-layer MLP model was used in this study and parameters were well tuned.

Model validation and evaluation
Sixty groups of multispectral data and field measured data in each growth stage were divided into training and validation sets by the 10fold cross-validation method. Each time, 90% and 10% of the data were employed as training and validation sets, respectively; this process continues 10 times until all the samples have been predicted once and only once. The model final performance was averaged from the evaluation criteria in the cross-validation. To comprehensively evaluate the model performance of Pn estimation, the mean square error (MSE), mean absolute error (MAE), explained variance score (EVS), and coefficient of determination (R 2 ) were considered in this study as the evaluation criteria. All model code

Vegetation index Abbreviations and formula
NIRv (Zeng et al., 2019 Chlorophyll Index Green (Gitelson et al., 2005) Chlorophyll Index Red Edge (Gitelson et al., 2005) Chlorophyll Vegetation Index (Gitelson et al., 2003) Greenness Index (Huete et al., 2002) GI Green Normalized Difference Vegetation (Smith et al., 1995 Modified Chlorophyll Absorption in Reflectance Index (Gitelson et al., 1996) Modified Soil Adjusted Vegetation Index (Gong et al., 2003) Modified Simple Ratio (Goel and Qin, 1994) Modified Triangular Vegetation Index (Dash and Curran, 2004) Nonlinear Vegetation Index (Haboudane et al., 2004) Normalized Difference Vegetation Index (Tucker et al., 1979) Optimization of Soil-Adjusted Vegetation Index (Rondeaux et al., 1996) Renormalized Difference Vegetation Index (Roujean and Breon, 1995 Ratio Vegetation Index 1 (Birth and McVey, 1968 Transformed Chlorophyll Absorption in Reflectance Index (Haboudane et al., 2002) Triangular Vegetation Index (Broge and Leblanc, 2001) Visible Difference Vegetation Index  VDVI where y i ,ŷ i and yepresent the measured value, the mean measured value, and the estimated value, respectively. n represents the number of the results. VAR represents the variance of the results. MSE and MAE are in the same unit with the measured value, ranging from 0 (optimum value) to +∞ (worst value). EVS and R 2 are dimensionless, ranging from 0 (worst value) to 1 (optimum value).

Results
3.1 Rice photosynthetic traits and canopy multispectral feature response to different treatments 3.1.1 Rice plant height, SPAD, and net photosynthetic rate The PH, SPAD, and Pn of tested rice at the jointing-booting, heading-flowering, and ripening stage are shown in Figure 2, respectively. PH increased obviously with the increase of nitrogen (N) application and advancement of growth stage, while it decreased slightly with ear filling at the ripening stage (Figure 2A). When PH reached the highest at the heading-flowering stage, the average rice PH with N2-N5 level under W1 (low leakage) treatment was 10.11%, 14.51%, 18.53%, and 21.23% higher than that under the N1 level, respectively, and 18.48%, 15.59%, 15.47%, and 15.59% higher than the N1 level, respectively, for W2 treatment. PH under W2 (high leakage) treatment was significantly higher than those under W1 treatment at N1 and N2 levels in the jointing-booting and heading-flowering stage, respectively, but the difference was not obvious under other N applications.
As shown in Figure 2B, SPAD generally showed a trend of initially increasing (N1-N3) then decreasing (N3-N4) and finally ending up with a small increase (N4-N5) with the increase of nitrogen application under the same leakage conditions. The maximum value of SPAD in each growth stage almost appeared at the N3 level, while the lowest value was found at the N1 level without N application. The SPAD value reached the maximum at the headingflowering stage, and the average SPAD values of N2-N5 levels were 12.27%, 8.30%, 1.62%, and 10.29% higher than the N1 level for W1 leakage treatment and 1.44%, 3.60%, 0.36%, and 1.08% higher than the N1 level for W2 leakage treatment, respectively. For the same N application level, W1 leakage treatment could increase the SPAD value of N1-N5 levels by −0.18%, 10.48%, 4.35%, 1.08%, and 8.91%, respectively, compared with W2.
It can be seen from Figure 2C that the Pn of rice decreased with the advancement of the rice growth stage. Pn under different N treatments showed a generally similar change trend with SPAD, increasing with the increase of N application at N1 to N3 levels and reaching the maximum at N3, decreasing at N4, and reverting at the N5 level [the Pn increase from N4 to N5 is not significant (p > 0.05)]. At the jointing-booting stage when photosynthesis was most vigorous, the average Pn values under N2-N5 levels with W1 treatment were 9.78%, 21.74%, 8.69%, and 13.04% higher than the N1 level, respectively, and 16.36%, 22.87%, 9.72%, and 17.16% higher than the N1 level under the conditions of W2 treatment, respectively, which indicated that excessive application of N fertilizer might inhibit photosynthesis. Under the same N application level, the Pn with W1 treatment was slightly higher than that with W2 treatment, indicating that low leakage intensity could promote leaf photosynthesis to a certain extent. Figure 3 illustrates the rice canopy multispectral reflectance of blue (band 1), green (band 2), red (band 3), red edge (band 4), and near infrared (band 5) with different treatments at three growth stages. In general, the average reflectance value of all five bands decreased as the growth stage progressed. As one of the most representative features in the crop spectral curve, the band 5 reflectance value ranged from 0.337 to 0.465 at the jointing-booting stage, while it decreased slightly to 0.327-0.449 at the heading- flowering stage, and finally dropped to approximately 0.030 at the ripening stage. Band 2 was the most intuitive band visible to the naked eye that could represent the nutrient statue and growth stage of crops, which reached a maximum of 0.160-0.221 at the jointing-booting stage, decreased to 0.139-0.174 with heading and flowering, and finally decreased to 0.027-0.038 with the yellowing of leaves and ears at the ripening stage. For the same N application level at the jointing-booting stage, the average reflectance values of band 1, band 2, and band 3 under the W1 leakage treatment were generally lower than those under the W2 leakage treatment, except that W1N3 had higher band 1 and band 3 values than W2N3. However, the average reflectance values of band 4 and band 5 showed opposite trends; specifically, W1 leakage treatment could increase the reflectance of band 4 and band 5 compared with W2 treatment and the N3 level improved the most. It could also be found from Figure 3A that lower leakage treatment had a steeper increase from band 3 to band 4, which indicated the better growth status. Under the same leakage treatment, the average reflectance values of band 1 to band 4 generally presented a trend of initially decreasing (N1-N3), then increasing (N3-N4), and finally ending up with a small decrease (N4-N5), while the reflectance value of band 5 increased with the N application level, but the law was not PH, SPAD and Pn of rice response to different treatments at different stages. (A) PH represents plant height (cm); (B) SPAD is relative chlorophyll content; (C) Pn represents the net photosynthetic rate (umol m-2s-1). W represents the leakage treatments (including W1: 3mm/day and W2: 5mm/day); N represents the nitrogen treatments (including N1-N5: 0, 75, 150, 225 and 300 kg/ha total pure nitrogen, respectively). The above ground measurements were conducted at the same time with UAV flight at 3 growth stages.

Rice canopy multispectral reflectance characteristics
obvious. Moreover, the average reflectance value of these five bands did not show a certain rule in the heading-flowering stage and the ripening stage, which might be due to the influence of rice heading and flowering, here leaves and ear yellowing in the growth process on the multispectral characteristics of the canopy.

Correlation analysis between VIs, TIs, and Pn
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 jointingboosting 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.

Estimation rice Pn from VIs at different growth stages
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 crossvalidation 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 mmol m −2 s −1 , an MAE of 0.414 mmol m −2 s −1 , an EVS of 0.938, and an R 2 of 0.938), while SVR models attained the worst performance (with an MSE of 2.512 mmol m −2 s −1 , an MAE of 1.155 mmol m −2 s −1 , an EVS of 0.390, and an R 2 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 mmol m −2 s −1 ; thus,

FIGURE 3
Average value of canopy band reflectance response to different treatments at different stages. (A-C) are the reflectance at jointing-booting, headingflowering and ripening stage, respectively. Band1-band5 represent blue (450±16 nm), green (560±16 nm), red (650±16 nm), red edge (730±16 nm), nearinfrared (840±26 nm) reflectance, respectively. W represents the leakage treatments (including W1: 3mm/day and W2: 5mm/day); N represents the nitrogen treatments (including N1-N5: 0, 75, 150, 225 and 300 kg/ha total pure nitrogen, respectively). 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 2 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 2 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.

Estimation rice Pn from fused VIs, TIs, and basal growth index
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 mmol m −2 s −1 , respectively, and the R 2 increasing by 0.022, 0.090, 0.031, 0.028, and 0.015, respectively. As the basal growth index further increased, the final R 2 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 2 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 2 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 headingflowering and ripening stages, especially for LR and MLP models.

Relationship between rice growth and canopy multispectral feature
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 nonphotochemical 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.

Limitations and suggestions on Pn estimation using VIs and TIs in this study
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 Correlation coefficient between Pn and TIs with different textural features combination. (A-C) are the correlation coefficient values at jointing-booting, heading-flowering and ripening stage, respectively. NDTI, DTI, RDTI represent normalized difference textural index, difference textural index and renormalized difference textural index, respectively. X-axis and Y-axis legends are the texture features in order of NO.1-40. The coloration in the thermal map is based on the correlation (r value) between TIs and Pn. 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.

Influence of input combination on Pn estimation performance
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 Comparison of estimated and measured Pn values of different models at certain stages. (A-C) are the comparison of estimated and measured Pn values at jointing-booting, heading-flowering and ripening stage, respectively. LR, SVR, GBDT, RF and MLP represent linear regression, support vector regression, gradient boosting decision tree, random forest and multilayer perceptron neural networks model, respectively. Pn represents the net photosynthetic rate (umol m-2s-1). The dashed blue and red circle in each figure are used to compare the fitting between the estimated and measured value when Pn value is low and high, respectively.
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  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.

Conclusion
This paper studied and revealed the responses of canopy multispectral band reflectance and rice net photosynthetic rate (Pn) to different nitrogen applications and leakage treatments through different growth stages under controlled irrigation and drainage schemes. The relationship between VIs, TIs, and Pn based on the UAV multispectral image was comprehensively analyzed and focused on. The performance of LR, SVR, GBDT, RF, and MLP models for Pn estimation under different input combinations was evaluated and compared at the jointing-booting, heading-flowering and ripening stages. The results indicated that the selected VIs and TIs had a relatively better correlation relationship with Pn at the jointingbooting stage, while only a moderate correlation at the heading-flowering stage and an unsatisfactory correlation at the ripening stage could be found. Therefore, the employed models generally had a better performance during the jointing-booting stage and the accuracy decreased as the growth stage progressed. Among the five used models, GBDT and RF models achieved the highest and most stable accuracy in the whole growth season and could be highly recommended for Pn estimation in the paddy field. Meanwhile, the fusion of VIs with TIs and basal growth index could significantly improve the model accuracy, and the plant height (PH) and SPAD had a better effect on performance improvement compared with NDTI, DTI, and RDTI employed in this study. The techniques and results presented in this paper could be valuable for rice field-scale photosynthetic monitoring, which could assist further stress detection and yield prediction.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
Conceptualization, TW and WZ. Methodology, TW and MC. Software, TW. Validation, TW, ZW and MC. Formal analysis, SW. Data curation, LQ. Writing-original draft preparation, TW. Writingreview and editing, TW, and LQ. Visualization, SW. Supervision, GS and XJ. Project administration, XJ. Funding acquisition, XJ. All authors contributed to the article and approved the submitted version.

Funding
This study is financially supported by the Jiangsu Provincial Key Research and Development Program (BE2022390) and the Postgraduate Research and Practice Innovation Program of Jiangsu Province (KYCX21_0537).