Estimating Plant Nitrogen Concentration of Rice through Fusing Vegetation Indices and Color Moments Derived from UAV-RGB Images

: Estimating plant nitrogen concentration (PNC) has been conducted using vegetation indices (VIs) from UAV-based imagery, but color features have been rarely considered as additional variables. In this study, the VIs and color moments (color feature) were calculated from UAV-based RGB images, then partial least square regression (PLSR) and random forest regression (RF) models were established to estimate PNC through fusing VIs and color moments. The results demonstrated that the fusion of VIs and color moments as inputs yielded higher accuracies of PNC estimation compared to VIs or color moments as input; the RF models based on the combination of VIs and color moments (R 2 ranging from 0.69 to 0.91 and NRMSE ranging from 0.07 to 0.13) showed similar performances to the PLSR models (R 2 ranging from 0.68 to 0.87 and NRMSE ranging from 0.10 to 0.29); Among the top ﬁve important variables in the RF models, there was at least one variable which belonged to the color moments in different datasets, indicating the signiﬁcant contribution of color moments in improving PNC estimation accuracy. This revealed the great potential of combination of RGB-VIs and color moments for the estimation of rice PNC.


Introduction
Rice (Oryza sativa L.) is one of the most important crops in the world, feeding more than half of the world's population [1]. Nitrogen (N) is the essential nutrient for rice growth, and is also an important limiting factor in soil productivity. Plant nitrogen concentration (PNC) has been commonly used as a crop N status indicator [2]. Timely and accurate assessment of PNC to detect N excess or deficiency is essential for farmers to improve rice production and N use efficiency [3,4]. Laboratory analysis is one of the important ways to obtain crop N nutrition status. However, it is time-consuming and laborious to carry out field investigations and collect representative samples, and usually the obtained results are delayed [5]. Compared with traditional laboratory analysis (e.g., micro Kjeldahl method), the non-destructive and timely methods or tools to detect crop N status have been significantly increased over recent decades [6][7][8]. Remote sensing with the advantages of fast and non-destructive characterizations has been proved useful for acquiring related information of crop nutrition status [9].
On regional or global scales, satellite remote sensing data was analyzed for retrieval N status [9,10]. Spectral wavebands in green, red, and red edge regions were used to calculate vegetation indices (VIs) to establish the relationship with crop N [7]. VIs and biophysical variables derived from Sentinel-2 satellite images were examined for estimation of plant N status [11]. Loozen et al. [12] used satellite-based VIs and environmental variables to map crop N in European forests. Yet today, the application of satellite-based VIs is still challenged for crop N status assessment because of the limited spatial resolution, the infrequency of satellite overpasses, and the risk of poor image data quality due to atmospheric conditions [13].
Unmanned aerial vehicle (UAV) platforms with the advantages of low cost and high spatio-temporal resolution of image data have become a promising approach in monitoring crop growth status [14,15]. Particularly, the light-weight and cost-efficient consumer-grade UAVs with RGB digital cameras can be simply and flexibly operated [16]. VIs from RGB images can distinguish crops from background soil and other interferences [17]. Moreover, color images contain a large amount of information about phenotypic traits such as biomass, plant height, and leaf senescence [18][19][20]. Several studies have proved the potential of VIs and color parameters derived from RGB images to obtain crop growth and nutrition status [17,21]. Jiang et al. [22] confirmed that the applicability of a new proposed true color vegetation index (TCVI) based on UAV-based RGB images could improve the accuracy of leaf N concentration (LNC) estimation for winter wheat under different conditions. The combination of UAV-based VIs and color parameters could yield reliable estimates of plant N density in winter wheat [23]. Rorie et al. [24] found that there was a close relationship of the dark green color index (DGCI) calculated from digital images with LNC in corn, which suggested that color image analysis was an appropriate tool for crop N nutrition status assessment at field scale. The application of VIs from RGB images to assess crop N status has been considered a promising alternative to expensive sensors such as hyperspectral and multispectral techniques; some studies have considered the fusion of texture and spectral information from UAV-based RGB images for monitoring crop growth status [4,25,26].
Besides the texture and spectral information, color feature is one of the most important features in RGB images. Among the color features, color moments are very simple and effective color features proposed by Stricker and Orengo [27]. The mathematical basis of this method is that any color distribution in the image can be represented by its moments. However, the contribution of color features for crop N estimation remains unclear.
In this study, VIs and color moments were used to estimate rice PNC based on partial least square regression (PLSR) and random forest (RF) regression algorithms [12,26,28], in which the objectives are (i) to examine the performance of VIs and color moments from UAV-based RGB images for PNC estimation in rice using algorithms of PLSR and RF; (ii) to explore the potential of improving rice PNC estimation accuracy by fusing VIs and color moments based on the above two methods.

Experiment Design
The field experiment ( Figure 1) was carried out at the Tangquan experimental station of the Nanjing institute of Soil Science, Chinese Academy of Sciences, Nanjing city, Jiangsu province, China (32 • 04 15 N, 118 • 28 21 E). The experimental station belongs to the subtropical monsoon climate zone. The average annual rainfall is 1102.2 mm. The average annual temperature is approximately 15.4 • C. The soil in the Tangquan experimental station was paddy soil with 22.26 g/kg organic matter, 1.31 g/kg total N, 15.41 mg/kg Olsen-P, and 146.4 mg/kg NH 4 OAc-K. In the experimental station, the rice cultivars were Wuyunjing 23 and Nanjing 5055, respectively in 2018 and 2019. Treatments consisted of five N fertilization rates with four replications in a total of 20 plots. All of the plots had dimensions of 4 m × 10 m. In Table 1, N0 represents no fertilizer was applied; N1 represents conventional fertilizer; N2, N3, and N4 represents N mixed by urea and coated urea with 30%, 40%, and 50% slow-release N respectively. In 2018, the rice plants were transplanted on 12 June and harvested on 22 October. In 2019, the rice plants were transplanted on 12 June and harvested on 12 November.

Determination of PNC
The ground destructive samplings were carried out at the tillering, jointing, and flowering stages ( Figure 1 and Table 2). One hill of rice was destructively and randomly sampled within the sampling region of each plot at each sampling time ( Figure 1). All of the plant samples were oven-dried at 105 ∘ C for 30 mins, followed by 80 ∘ C until constant weight and then ground for chemical analysis in the laboratory. PNC (%) was determined by using the micro-Kjeldahl method [29].  The ground destructive samplings were carried out at the tillering, jointing, and flowering stages ( Figure 1 and Table 2). One hill of rice was destructively and randomly sampled within the sampling region of each plot at each sampling time ( Figure 1). All of the plant samples were oven-dried at 105 • C for 30 mins, followed by 80 • C until constant weight and then ground for chemical analysis in the laboratory. PNC (%) was determined by using the micro-Kjeldahl method [29].

Image Acquisition
The acquisition of UAV-based images was implemented on the same dates as those of rice plant sampling ( Table 2). A consumer-grade camera mounted on a UAV (Phantom 4 Pro, SZ DJI Technology Co., Ltd., Shenzhen, China) was used in this study. The UAV had flight duration of about 30 mins, depending on the actual working conditions. The digital camera was equipped with a one-inch complementary metal-oxide semiconductor (CMOS) sensor, which had a spatial resolution of approximate 20 mega pixels. The UAV was flown over the rice fields at fight altitudes of 50 m above ground level (ground sampling distance: 1.36 cm). The flights were carried out in stable ambient light conditions between 11:00 am and 13:00 pm by setting an auto-exposure time. The UAV aviated automatically based on the preset flight waypoints, leading to approximately 80% forward overlap and 60% side overlap controlled by the Pix4Dcapture application.

Image Mosaic
The orthophotos of the experimental site were generated by the Pix4Dmapper Version 1.1.38 (https://www.pix4d.com/, accessed on 10 September in 2019), which provides an automated pipeline through steps of image alignment, matching, mosaicking, constructing dense point cloud, and finally generating orthoimages with the Geo-TIF format.

Calculation of VIs
Before calculation of VIs from the orthophotos at multi stages, an empirical line correction method was used to carry out radiometric calibration [30]. In addition, we used the optimal index method, proposed by Qiu et al. [31], to remove the background effects. Firstly, the VI which had the best correlation with PNC at each phenological stage was selected as the optimal vegetation index (OVI). Then, we used ArcGIS 10.1 (ESRI, Redlands, CA) to divide OVI into five levels based on the natural fault zone method. Finally, the three middle levels were selected as the canopy area. After background removal, ArcGIS 10.1 was used to draw the region of interest (ROI) of each plot in the non-sampling area ( Figure  1). In this study, most of the selected color indices were well-known and have been studied in leaf chlorophyll content, biomass, and N status estimation [13,17,26,32]. The R, G, and B channels from ortho-images were used for calculating twelve VIs (Table 3). Then, we used the "Zonal Statistics As Table" function in ArcGIS 10.1 to extract the mean reflectance value of ROI in each plot.

Calculation of Color Moments
Due to the information of color distribution being mainly concentrated in low order moments, only the first, second, and third order moments of color could be used to express the color distribution. Compared with color histograms which divide the color space into several small intervals and calculate the number of pixels in each interval, another advantage of this method was that it did not need color space quantization and had lower dimensions of feature vector. Accordingly, 9 components (3 color components, 2 lower order moments on each component) could cover the color moments of an image, which is very concise compared with other color features. The color moments were calculated as follows [39]: where p i,j represents the probability of the pixels of the gray value with j in i th color channel. N was the total number pixels in the ROI of each plot ( Figure 1). The entries µ i (1 ≤ i ≤ 3) represent the average color in each color channel. The entries σ i and s i represent the variance and skewness of each color channel, respectively. In this study, we used the HSV color space, and then the nine color features which are expressed in Table 3.

PLSR
PLSR is one of the most common multivariate regression algorithms applied to deal with data with collinear variables. PLSR was inherent in the multiple linear regression (MLR) and able to avoid multi-collinearity [40]. In this study, we used PLSR to examine Remote Sens. 2021, 13, 1620 6 of 18 the relationships of RGB-VIs and color moments with the rice PNC. The number of latent variables was a critical parameter in the PLSR algorithm, and was determined following the method described by Fassio, Cozzolino et al. [41]. The criterion of the selected latent variables was based on the minimum predicted residual sum of squares (PRESS) [42] with 10-fold cross validation in calibration datasets.

Random Forests
The algorithm of Random forests was developed by Breiman et al. [43] and shows a promising capability to avoid overfitting by sampling the predictor space randomly. It can construct non-linear relationships without the limitations of the assumptions of variable distributions and dependency [12]. In addition, RF could effectively evaluate the importance of independent variables and partly deal with multicollinearity among variables while having great tolerance to noises and outliers. In this study, we implemented the RF in the python 3.6 environment using the "RandomForestRegressor" function in "sklearn" package. Three tuning parameters (max_depth, min_samples_split, min_samples_leaf) were implemented in RF by using "GridSearchCV" module with 10-fold cross validation in calibration datasets. Other parameters were set to default values.

Statistical Analysis
In this study, we used the data in 2019 as the calibration set and the data in 2018 as the validation set ( Figure 2). Before establishment of estimation models, the correlations between VIs, color moments, and PNC were analyzed using Pearson's correlation coefficient.
Paired t-tests was conducted to analyze the variations between the measured and estimated PNC values in calibration and validation datasets by using "ttest_real" function from "scipy" package in python 3.6 environment. The means were compared at the 5% level of significance by the t-tests [44]. The fitness was assessed from a 1:1 line of the estimated and measured PNC values. The performance of models to estimate PNC was evaluated by comparing the differences in coefficient of determination (R 2 ) in Equation (4) and normalized root mean square error (NRMSE) in Equation (5). These statistical indicators were expressed as follow: where n was the number of observations, O denoted the average value of measured PNC, P i and O i were the estimated and observed values of PNC, respectively. Generally, the simulation was considered excellent when NRMSE was less than 10%, good if NRMSE was greater than 10% and less than 20%, fair if NRMSE was greater than 20% and less than 30%, and poor if NRMSE was greater than 30% [45].
where n was the number of observations, denoted the average value of measured PNC, and were the estimated and observed values of PNC, respectively. Generally, the simulation was considered excellent when NRMSE was less than 10%, good if NRMSE was greater than 10% and less than 20%, fair if NRMSE was greater than 20% and less than 30%, and poor if NRMSE was greater than 30% [45].

Descriptive Analysis of Measured PNC Data
Across growth stages, years, and applied N rates, the PNC of rice ranged from 0.7% to 3.5% (Table 4). Because of the N "dilution effect" [46], the average values of the observed PNC decreased from the tillering stage to the flowering stage in the calibration and validation sets. In the calibration set, the average PNC values decreased from 2.3% at the tillering stage and 1.7% at the jointing stage to 1.1% at the flowering stage. In the validation set, the average PNC values were close to those in the calibration set at the same growth stage which decreased from 2.5% at the tillering stage and 1.8% at the jointing stage to 1.0% at the flowering stage.

Descriptive Analysis of Measured PNC Data
Across growth stages, years, and applied N rates, the PNC of rice ranged from 0.7% to 3.5% (Table 4). Because of the N "dilution effect" [46], the average values of the observed PNC decreased from the tillering stage to the flowering stage in the calibration and validation sets. In the calibration set, the average PNC values decreased from 2.3% at the tillering stage and 1.7% at the jointing stage to 1.1% at the flowering stage. In the validation set, the average PNC values were close to those in the calibration set at the same growth stage which decreased from 2.5% at the tillering stage and 1.8% at the jointing stage to 1.0% at the flowering stage.

Correlations between RGB-VIs, Color Moments and PNC
To evaluate the performance of the 12 RGB-VIs and 9 color moments obtained from the UAV-based images, the Pearson correlation analysis between PNC and these variables was implemented across the growth stages (Figure 3). At the tillering stage, high correlations were found between RGB-VIs and rice PNC such as NRI (r = −0.77) and VARI (r = 0.77). Except for S_var (p > 0.05), other color moments were significantly correlated with PNC and the absolute values of r ranged from 0.58 to 0.83. At the jointing stage, very strong correlations were found between RGB-VIs and PNC such as G/R (r = 0.89) and NGRDI (r = 0.89). Color moments showed slightly lower correlations with PNC. The top three color moments were V, S, H and the absolute values of r ranged from 0.67 to 0.88. At the flowering stage, G/R revealed the highest correlation with PNC reaching 0.84 among all RGB-VIs. In addition, H had the highest correlation with PNC reaching 0.77 among color moments. three color moments were V, S, H and the absolute values of r ranged from 0.67 to 0.88. At the flowering stage, G/R revealed the highest correlation with PNC reaching 0.84 among all RGB-VIs. In addition, H had the highest correlation with PNC reaching 0.77 among color moments.
Among the VIs and color moments, some variables were constantly well correlated with PNC across the growth stages such as G/R, NGRDI, INT, ExR, H, and V with r values ranging from 0.74 to 0.89, 0.74 to 0.89, −0.7 to −0.87, −0.72 to −0.85, 0.67 to 0.77 and −0.64 to −0.88, respectively. However, for some other variables, such as NRI, NGI, VARI and H_ske (Figure 3), the correlations changed randomly even when they had a good correlation with PNC at the three individual stages.  (Figure 3), the correlations changed randomly even when they had a good correlation with PNC at the three individual stages. Figure 4 demonstrates the change process of the mean PRESS values with the increase of latent components by using the 10-fold cross-validation method in the calibration datasets. When the number of PLSR components increased, the mean PRESS values decreased first of all, then increased slowly, and remained relatively stable finally except for the flowering stage. When using all the RGB-VIs as input variables, the PLSR models achieved minimum mean PRESS values with 2, 2, 2, and 3 components, respectively at the tillering, jointing, flowering, and combined stages. While using all the color moments as inputs, the PLSR models attained the minimum mean PRESS values with 5, 4, 1, and 3 components, respectively at the tillering, jointing, flowering, and combined stages. When combining all the RGB-VIs and color moments as inputs, the minimum mean PRESS values Remote Sens. 2021, 13, 1620 9 of 18 were achieved by the PLSR models with 6, 2, 6, and 6 components at the tillering, jointing, flowering, and combined stages, respectively. of latent components by using the 10-fold cross-validation method in the calibration datasets. When the number of PLSR components increased, the mean PRESS values decreased first of all, then increased slowly, and remained relatively stable finally except for the flowering stage. When using all the RGB-VIs as input variables, the PLSR models achieved minimum mean PRESS values with 2, 2, 2, and 3 components, respectively at the tillering, jointing, flowering, and combined stages. While using all the color moments as inputs, the PLSR models attained the minimum mean PRESS values with 5, 4, 1, and 3 components, respectively at the tillering, jointing, flowering, and combined stages. When combining all the RGB-VIs and color moments as inputs, the minimum mean PRESS values were achieved by the PLSR models with 6, 2, 6, and 6 components at the tillering, jointing, flowering, and combined stages, respectively.  Table 5 shows the results from the PLSR models. In the calibration datasets, the PLSR models using color moments as inputs obtained slightly better results than those using RGB-VIs as inputs, with R 2 ranged from 0.72 to 0.89 and NRMSE ranged from 0.10 to 0.18 across the single stages and combined stages. Except for the jointing stage, the PLSR using all variables performed best, with R 2 values of 0.79, 0.80. and 0.83, respectively at the tillering, flowering, and combined stages. In validation datasets, the models using all variables showed higher accuracy and lower NRMSE values compared to the models using either RGB-VIs variables or color moment variables. The values of R 2 and NRMSE ranged from 0.68 to 0.87 and 0.10 to 0.29, respectively at the single and combined stages. In addition, the color moments only type of model showed lower R 2 values (0.32 and 0.33) compared to the model using only RGB-VIs variables (R 2 = 0.63 and 0.60), respectively at the tillering and flowering stages.  Table 5 shows the results from the PLSR models. In the calibration datasets, the PLSR models using color moments as inputs obtained slightly better results than those using RGB-VIs as inputs, with R 2 ranged from 0.72 to 0.89 and NRMSE ranged from 0.10 to 0.18 across the single stages and combined stages. Except for the jointing stage, the PLSR using all variables performed best, with R 2 values of 0.79, 0.80. and 0.83, respectively at the tillering, flowering, and combined stages. In validation datasets, the models using all variables showed higher accuracy and lower NRMSE values compared to the models using either RGB-VIs variables or color moment variables. The values of R 2 and NRMSE ranged from 0.68 to 0.87 and 0.10 to 0.29, respectively at the single and combined stages. In addition, the color moments only type of model showed lower R 2 values (0.32 and 0.33) compared to the model using only RGB-VIs variables (R 2 = 0.63 and 0.60), respectively at the tillering and flowering stages. Scatterplots of estimated and measured values for PNC across single and combined stages are presented in Figure 5. The black line is the 1:1 line used to observe the distribution of scatters. At the flowering stage, the PLSR model showed over-estimated PNC values at high level of the measured PNC values in the validation dataset when using either RGB-VIs or color moment variables (Figure 5c,g). The paired t test analysis showed significant difference (p < 0.05) between the estimated and measured PNC in the validation datasets at flowering stages for models using RGB-VIs only or color moments only under no N stress treatment. However, the overestimation of PNC values could be significantly eliminated (p > 0.05) by the PLSR model when considering all variables (Figure 5k). Generally, the distributions of scatters obtained from models using all variables as inputs were closer to the 1:1 line compared to those of scatters obtained from models using RGB-VIs only or color moments only as inputs.   Table  1; "Full N" represents no N stress treatment with the combination of N treatments ranging from N1 to N4 in Table 1.

RF Analysis
Before estimating the PNC by RF models, we used the method of grid search with 10-fold cross validation to obtain the optimal model parameter sets in the calibration datasets (Table 6). Then, the top five important variables for rice PNC estimation for all RF Figure 5. Scatterplots between measured and estimated PNC values (%N) based on PLSR models in the calibration and validation datasets across single and combined stages. The embedding histogram in each subplot represents the means ± standard of error in different N treatments in the calibration and validation datasets. Note: "N0" represents N0 in Table 1; "Full N" represents no N stress treatment with the combination of N treatments ranging from N1 to N4 in Table 1.

RF Analysis
Before estimating the PNC by RF models, we used the method of grid search with 10-fold cross validation to obtain the optimal model parameter sets in the calibration datasets (Table 6). Then, the top five important variables for rice PNC estimation for all RF models are presented in Table 7. For the RF models based on only RGB-VIs, the rank of important variables was inconsistent across the growth stages. At the tillering stage, the difference of the importance values was relatively small and the values of importance ranged from 0.1 to 0.16. At the flowering and combined stages, the most important variables were ExR and NRI, respectively. The values of importance (0.35 and 0.44) of these two variables showed a great influence for PNC estimation. For the RF models based on only color moments, the two most important variables had a dominant effect on the model performance at the jointing, flowering and combined stages while the importance values reached over 0.79. For the models including all variables as inputs, there was at least one out of the five most important variables which belonged to the color moments at the single and combined stages.   Table 8 displays the calibration and validation statistic of estimated PNC from the RF models. In calibration datasets, the RF models including only RGB-VIs performed well with R 2 ranging from 0.79 to 0.93 and NRMSE from 0.10 to 0.19. Only color moments as inputs could get similar results to those using only RGB-VIs as inputs, with R 2 ranging from 0.73 to 0.94 and NRMSE from 0.10 to 0.16 across the single and combined stages. The models including RGB-VIs and color moments yielded superior performance compared to the models using either RGB-VIs or color moments. The values of R 2 and NRMSE ranged from 0.83 to 0.95 and 0.08 to 0.15, respectively at the single and combined stages. In the validation datasets, all the models obtained slightly lower R 2 values compared to those in the calibration datasets. In contrast, the NRMSE values were still low and ranged from 0.07 to 0.11 across the single stages. In addition, the fusion of RGB-VIs and color moments improved the accuracy of the PNC estimation with R 2 ranging from 0.69 to 0.91 and NRMSE from 0.07 to 0.13 compared to those models using RGB-VIs only or color moments only as inputs. This indicated the RF models could effectively predict PNC in rice.  Figure 6 displays the scatterplots of the estimated and measured values for PNC across the single and combined stages. In the calibration datasets, the distributions of scatters are close to the 1:1 line. However, a slight overestimation of PNC values occurred in the low level of PNC as indicated by the deviation of the 1:1 line when using only color moments as inputs at the flowering and combined stages in the validation datasets (Figure 6g,h). There was significant difference (p < 0.05) between the estimated and measured PNC when using color moments as inputs for the RF model in the validation dataset at the flowering stage under different N conditions. When using all variables as inputs, the distributions of scatters obtained from the RF models were close to the 1:1 line in the calibration and validation datasets (Figure 6 i,j,k,l).  Figure 6. Scatterplots between measured and estimated PNC values (%N) based on RF models in the calibration and validation datasets across the growth stages. The embedding histogram in each subplot represents the means ± standard of error in different N treatments in the calibration and validation datasets. Note: "N0" represents N0 in Table 1; "Full N" represents no N stress treatment with the combination of N treatments ranging from N1 to N4 in Table 1.

Comparisons of the PLSR and RF Models
In this study, the performance of the PLSR and RF models was tested to assess rice PNC under varying N fertilizer application rates at different growth stages. These two kinds of models were multivariate regression methods and good at dealing with a lot of predicators which were cross-correlated [25,47]. Based on the results of Table 5 and Table  8, the RF regression models outperformed the PLSR regression models for PNC estimation in the same validation dataset. In similar studies, Maimaitijiang et al. [48] compared Figure 6. Scatterplots between measured and estimated PNC values (%N) based on RF models in the calibration and validation datasets across the growth stages. The embedding histogram in each subplot represents the means ± standard of error in different N treatments in the calibration and validation datasets. Note: "N0" represents N0 in Table 1; "Full N" represents no N stress treatment with the combination of N treatments ranging from N1 to N4 in Table 1.

Comparisons of the PLSR and RF Models
In this study, the performance of the PLSR and RF models was tested to assess rice PNC under varying N fertilizer application rates at different growth stages. These two kinds of models were multivariate regression methods and good at dealing with a lot of predicators which were cross-correlated [25,47]. Based on the results of Tables 5 and 8, the RF regression models outperformed the PLSR regression models for PNC estimation in the same validation dataset. In similar studies, Maimaitijiang et al. [48] compared the performance of PLSR, RF, extreme learning regression (ELR), and support vector regression (SVR) in estimating crop LNC, using satellite-based VIs and UAV-based canopy structure information as inputs. Their results showed that the RF model performed best for N estimation. Osco et al. [49] evaluated the performance of nine machine learning models for maize LNC predictions using UAV-based multispectral imagery. Their results indicated that the RF model performed better than other models for LNC. Liang et al. [50] also concluded that the RF model was preferred to predict LNC compared to the least square support vector model (LS-SVR). Our results were consistent with the aforementioned studies. Accordingly, the RF model based on decision tree algorithms is more robust and appropriate to assess crop N status compared to other machine learning algorithms.

Fusion of RGB-VIs and Color Moments for PNC Estimation
Many previous studies proposed new VIs to improve crop N status estimation in their publications [7,51]. However, some other studies also concluded that these VIs had poor and unstable performance in N status estimation [4,52,53]. In recent years, texture indices have been successfully used to monitor crop growth status for precision agriculture [4,25]. Compared with texture or shape information, it was found that it was easier to extract color moments [39]. In addition, color moments based on color distribution features could be matched more efficiently and robustly than color histograms [27]. In this study, the color moments were for the first time investigated for crop N estimation. Generally, the combination of color feature information (color moments) and RGB-VIs derived from UAVbased imagery improved the performance of PLSR and RF models for PNC estimation at the single and combined stages (Tables 5 and 8).
In the PLSR models, the number of latent components was the only tuning parameter. The PLSR models which used VIs or color moments as inputs, tended to obtain minimum Mean PRESS values when the numbers of components were between 2 and 5 (Figure 4a,b). When combing VIs and color moments as inputs, the optimal number of components was approximately six at the individual and combined stages (Figure 4c). This indicated that when using multiple types of variables as inputs, more valuable components would be screened out by the PLSR models to improve model accuracy compared to the model with a single type of variables as inputs. From the rank of the top five feature importance values in the RF models, there was at least one feature variable from color moments at the single and combined stages (Table 7). Additionally, the first three moments of H demonstrated a considerable contribution in improving model performance in RF models at the tillering stages when using all variables as inputs (Table 7). Thus, the models when adding color moments had the potential to improve the accuracy of crop N status estimation. Although, the fusion of color moments and RGB-VIs reduced the NRMSE values to a certain extent, the NRMSE values of RF and PLSR models still reached 0.15 to 0.17 and 0.13 to 0.29, respectively at the combined stages in the calibration and validation datasets. This was partly due to the changing variation of the morphological characteristics in bringing in several variables (NRI, NGI, or VARI) which changed their correlation at the different phenological stages though they were well correlated with PNC ( Figure 3). When these variables were included as inputs at the combined stages, the model resulted in larger uncertainty compared to models based on the single stage.

Implications of UAV-Based RGB Imagery for Crop Monitoring
In recent decades, various sensors mounted on UAVs have been widely used to monitor crop growth traits such as LiDAR [54], hyperspectral [55] and multispectral [26] devices, NIR [56][57][58], and RGB cameras [59,60]. Compared to other sensors, the UAV-RGB camera was cheaper and lighter which made it possible for a longer working time to apply to larger areas in precision agriculture. With the advantage of ultrahigh-ground-resolution, the UAV-RGB images could be used for precision classification of crops and estimation of crop biochemical parameters. In addition, several previous studies demonstrated that the Digital Surface Models (DSMs) derived from UAV-RGB imagery were efficient in extracting a vegetation canopy structure [61][62][63]. Thus, the ultra-high resolution images from UAVs make it more powerful for field-scale applications in precision agriculture. In this study, we used the data fusion of color moments and VIs from UAV-RGB imagery to estimate rice PNC across cultivars and growth stages. The performance of the models was comparable to previous studies [3,23].
Although the performances of RF models were better by fusion of color moments and RGB-VIs than those of PLSR models, the other promising methods such as ensemble models [64], Long Short-Term Memory (LSTM) [65,66], and convolutional neural network (CNN) [67,68] need to be further explored to improve the accuracy of crop N estimation. Furthermore, the imageries over multiple years in more locations with different soil characteristics and climate conditions are essential for practical application.

Conclusions
The incorporation of color moments for both PLSR and RF models yielded more accurate estimations of PNC in rice with NRMSE reducing by 9% to 58% compared to results of models based on RGB-VIs across single and combined stages. The RF models with VIs and color moments showed a similar performance to the PLSR models for PNC estimation. The accuracies of RF models based on color moments only (R 2 = 0.57-0.75 and NRMSE = 0.08-0.11) were comparable to results of RF models using VIs only (R 2 = 0.62-0.82 and NRMSE = 0.08-0.09) at single stages in validation datasets. The analysis of important variables in RF models showed that color moments played an important role in PNC estimation, which demonstrated that color moments significantly improved the estimation of rice PNC. Future work should be undertaken to further improve the accuracy of the model for monitoring crop growth status by combining the texture, color features, structure from motion (sfM), and VIs. Meanwhile, these features are easily and readily extracted from UAV-based images.