Estimation of Plant Height and Biomass of Rice Using Unmanned Aerial Vehicle

Abstract

Plant height and biomass are important indicators of rice yield. Here we combined measured plant physiological traits with a crop growth model driven by unmanned aerial vehicle spectral data to quantify the changes in rice plant height and biomass under different irrigation and fertilizer treatments. The study included two treatments: I—water availability factor (i.e., three drought objects, optimal, and excess water); and II—two levels of deep percolation and five nitrogen fertilization doses. The introduced model is extreme learning machine (ELM), back propagation neural network (BPNN), and particle swarm optimization-ELM (PSO-ELM), respectively. The results showed that: (1) Proper water level regulation (3~5 cm) significantly increased the accumulation of spike biomass, which was about 6% higher compared to that under flooded conditions. (2) For plant height inversion, the ELM model was optimal with a mean coefficient of determination of 0.78, a mean root mean square error of 0.26 cm, and a mean performance deviation rate of 2.08. For biomass inversion, the PSO-ELM model was optimal with a mean coefficient of determination of 0.88, a mean root mean square error of 3.8 g, and a mean performance deviation rate of 3.29. This study provided the possible opportunity for large-scale estimations of rice yield under environmental disturbances.

1. Introduction

Rice is one of the most important food crops in China and the world, providing about 19% of the energy for everyone on the planet and feeding nearly half of the world’s population [1,2]. Currently, rice covers about 25.8% of China’s area, nitrogen is the most sensitive factor in the nutrient elements affecting the growth, development, and yield of rice, and reasonable nitrogen application is the key measure for high yield and high quality of rice. At present, the utilization rate of nitrogen fertilizer in China’s farmland is far below the level of developed countries, only 30% to 35% [3], a large number of unused nitrogen fertilizers has led to a series of problems, such as rising soil fertility, the yield effect of fertilizers being gradually reduced, and the amount of nitrogen fertilizer residue per unit area has increased [4,5]. This has become one of the important causes of water environment and air pollution in China, which poses a threat to the sustainable development of China’s agriculture as well as to the survival of the human environment [6,7]. The latest research shows that crop production in the Yangtze River Basin flows a total of about 6 Tg of nitrogen into the water system, of which oilseed rape accounts for about 7.95%, and rice is 22.31%, which shows that rice is one of the main contributors to the nitrogen load of the water system [8].

It is of great significance to study the effects of different irrigation levels and nitrogen fertilizer application on the growth status of rice, to grasp the growth status of rice in real-time, and to ensure the adequate supply of rice, in order to maintain national food security, realize the green high yield, and high efficiency of crops and the sustainable development of agriculture.

Aboveground biomass (AGB) is an important indicator of crop growth [9]. Accurate measurement of AGB is a prerequisite for reasonable yield estimation during the whole rice crop period [10], while the weight of the spike (SB) in AGB directly reflects the weight of each grain and the percentage of grains in the spike. Therefore, proper monitoring of rice spike biomass is essential for precision agriculture. Similarly, plant height (PH) greatly influences the yield potential of rice and is a common indicator of rice growth status [11].

Currently, a common method for measuring rice spike biomass is direct measurement in the field. Although this method is straightforward and reliable, it requires separate destructive sampling of plants at different stages, which results in high crop losses [12]. At the same time, the process of acquiring and recording crop height is monotonous, labor-intensive, and time-consuming when the study area is large, and the accuracy of the results depends on the number and spatial distribution of the measurement samples.

With the rapid development of remote sensing technology, the combination of remote sensing technology with agriculture is becoming more and more common due to its short detection period, strong presentability, and the ability to observe large areas simultaneously [13]. Guo et al. [14] proposed a new method based on Random Volume over Ground (RVoG) model for rice height inversion in combination with PolInSAR; Sun et al. [15] introduced the time series Sentinel-1A data for the full growth period of sugarcane and used four empirical regression models to invert the height of sugarcane plants; and Wu et al. [16] developed an inversion model for leaf area index (LAI) of winter wheat based on parameter sensitivity analysis and model calibration using RADARST-2 data as the image source, which was feasible and applicable. Zhang et al. [17] extracted 27 polarization parameters from Radarsat-2 polarization synthetic aperture radar and studied the best inversion model for plant height, leaf area index, and biomass inversion based on sensitivity analysis. However, most of the above studies are based on satellite images such as Landsat, which have a low ground resolution and cannot guarantee the accuracy of crop growth estimation and are greatly affected by cloud cover and return time [18]. In addition, the large scale of satellite remote sensing does not meet the practical needs of small plots or research parks in southern China [19].

In recent years, unmanned aerial vehicles (UAV) have become an important tool for agricultural remote sensing due to their unique advantages such as fast, real-time, and high resolution in vegetation indices based on multispectral images [20,21,22]. Han et al. [23] used a UAV to collect crop and spectral data at different growth stages of wheat and constructed inverse models of aboveground biomass (AGB) and plant nitrogen accumulation (PNA) of wheat by band reflectance and vegetation index, and achieved good inversion results. Zhang et al. [24] used a combination of UVA and field-measured 1m² sample models to calculate a canopy height metric for each sample, enabling remote inversion of above-ground biomass at the plot scale. Prakriti Sharma et al. [25] used UAV to acquire multispectral data of oats before and after tasseling at different planting locations and used four models to construct the oat biomass estimation model, which explored the relationship between model prediction accuracy and planting area. However, the current research on inversion for vegetation indices ignores the fact that vegetation indices are the result of linear or non-linear combinations of different bands with very strong linear correlations between them, and that simply introducing various models may affect the inversion results and lead to distorted results.

In general, crop simulation models and statistical (machine) learning models are two common modeling approaches for crop yield and growth estimation. Crop simulation models incorporate crop growth information (such as plant physiological characteristics, nutrients, and environmental factors) to predict crop productivity [26]. The predictions are superior, but they require a large amount of actual crop data as input variables. In addition, environmental differences lead to poor generalization and generalizability of such models. Machine learning-based regression methods are not based directly on the physiological mechanisms of the plant, but rather on the prediction of various physiological indicators of the crop based on the predicted response information set (spectral data, measured value, etc.). In recent research, a number of scholars have proposed various machine algorithm models and demonstrated the superior performance of machine learning algorithms for crop indicator inversion [27,28,29]. Huang [30] proposed the Extreme Learning Machine (ELM), which has been shown to have extremely strong regression capabilities. Compared to the commonly used artificial neural network models and support vector machine models, ELM has better learning speed and generalization performance [31]. However, since the ELM model is based on a random allocation of the number of hidden layers, more hidden neurons are required under the same conditions than other tuning algorithms, increasing the complexity of the model [32]. In contrast, Particle Swarm Optimization (PSO) automatically learns the number of hidden layers from the training samples and can be used for optimization purposes [33].

Therefore, to address the problems of labor and yield loss caused by direct measurement in the field, the large reentry period and scale of satellite remote sensing, in this study, we set up rice experiments under different irrigation and nitrogen application rates to study their effects on rice plant height and biomass, and constructed a rice growth model based on multispectral images of rice in the experimental area obtained by UAV during different fertility periods. The objectives of this study are: (1) to determine appropriate irrigation and nitrogen fertilizer input thresholds for rice; (2) to develop ELM, BPNN, and PSO-ELM models for rice plant height and biomass during the plucking and gestation period, the spike and flowering period and the milking period, and to evaluate their accuracy; and (3) to develop a low-cost, non-destructive modeling method for estimating rice plant height and biomass by comparing the accuracy of different machine learning algorithms in the inversion of rice plant height and biomass, so as to meet the needs of precision agriculture development. This study provides a theoretical basis for large-scale and non-destructive monitoring of rice physiological indicators.

2. Materials and Methods

2.1. Study Area

The study area is located at the Jiangsu Engineering Research Centre for Efficient Use of Agricultural Soil and Water Resources and Carbon Sequestration and Emission Reduction (118°47′ E, 31°54′ N), as seen in Figure 1. The study area is a subtropical monsoon climate zone with an annual average temperature of 17.2 °C, a maximum temperature of 38.3 °C and a minimum temperature of −3.1 °C, an annual average humidity of 74% and a minimum humidity of 11%, a frost-free period of 259 d, an annual average sunshine duration of 1917.3 h, rain and heat at the same time, and an annual precipitation of 1106.5 mm.

2.2. Experimental Design

The rice variety used in this study was South Japonica 9108, which was transplanted on 5 July 2021, tillered on 13 July, plucked on 13 August, spiked and flowered on 26 August, milked on 8 September, and harvested on 26 October.

There were two trials with a total of 13 treatment groups (3 irrigation treatments and 10 nitrogen treatments). Three irrigation treatments were used in Trial 1: drought stress (W1), regular irrigation (W2), and flood stress (W3), while 5 different nitrogen treatments (N1~N5) and 2 leaching treatments (S1, S2) were used in Trial 2: S1N1, S1N2, S1N3, S1N4, S1N5, S2N1, S2N2, S2N3, S2N4 and S2N5. All agronomic treatment measures in this experiment were the same or similar to local agronomic measures, except for different gradients set under water and nitrogen regulation conditions to increase the diversity of sample data.

Table 1. Experiment design.

Test Treatment	Test Facilities	Process Number	Repeat	Water Level	Nitrogen Application kg/hm2	Deep Percolation mm/d	Measurement Indicators
IN	Potted plants	W1	8	80~85% Field Moisture Capacity	225	3	Plant height Biomass
		W2	8	3~5 cm	225	3
		W3	8	15~20 cm	225	3
NN	Soil pit method	S1N1	6	3~5 cm	0	3	Plant height
		S1N2	6	3~5 cm	150	3
		S1N3	6	3~5 cm	225	3
		S1N4	6	3~5 cm	300	3
		S1N5	6	3~5 cm	375	3
		S2N1	6	3~5 cm	0	5
		S2N2	6	3~5 cm	150	5
		S2N3	6	3~5 cm	225	5
		S2N4	6	3~5 cm	300	5
		S2N5	6	3~5 cm	375	5

Note: IN and NN denote irrigation treatment and nitrogen application treatment, respectively.

shows the treatment details of this study (irrigation treatments and nitrogen treatments are independent of each other).

2.3. Data Collection and Processing

Considering the main growth periods of rice, the UAV multispectral remote sensing images and ground data were collected simultaneously on 20 August (jointing and booting stage, JB), 2 September (heading and flowering stage, HF), and 13 September (milk ripening stage, MM) in 2021. Among them, 24 sampling points (plant height and biomass at different fertility stages) were set up for different irrigation treatments and 60 sampling points (plant height at different fertility stages) were set up for different nitrogen application treatments. The UAV flight time was concentrated at 11:00~14:00, with clear weather and little wind; the ground data collection time was concentrated at 09:00~14:00.

2.3.1. In Situ Data Collection

Plant height: Three representative rice plants were selected to determine the height of the above-ground plant, from the surface of the soil to the highest leaf tip before flowering and from the surface of the soil to the top of the spike (without awning) after flowering.

Biomass (dry matter mass): For different irrigation treatments, at HF and MM stages, three well-grown rice plants were selected for each treatment as the average biomass of that treatment. The samples under the above three treatments were brought back to the laboratory and different organs of rice (e.g., roots, stems, leaves, and spikes) were killed at 105 °C for 30 min and dried at 85 °C, and the dry matter mass of each organ was weighed.

2.3.2. Multispectral Data Pre-Processing and Analysis

While collecting rice physiological information, the quasi-simultaneous use of the DJI Elf 4 Pro Multispectral Edition (P4_Multispectral) drone was used to acquire multispectral raw data by flying over the test area using its Mica Sense Red Edge five-channel multispectral sensor. The flight plan was set to fly at an altitude of 15 m, with a heading overlap of 85%, a side overlap of 75%, and a ground resolution of 7.25 mm/pixel.

The multispectral images acquired by the UAV were stitched and processed using Agisoft Metashape Professional software 2.1.0, and the processed images were subjected to waveband fusion and radiometric correction in combination with ENVI 5.5.3. According to the GPS coordinates of the actual measurement point, a 5 × 5-pixel size image was taken at the center of the actual ground measurement point, and the average reflectance of the samples in the area was taken as the spectral reflectance of that sampling point to obtain the spectral reflectance data of different bands.

2.3.3. Vegetation Index Construction

Vegetation indices are used to effectively differentiate features by amplifying the variability between different feature classes through operations (including addition, subtraction, multiplication, and division) between spectra and are important parameters for characterizing crop growth information. In this study, 22 vegetation indices were selected as indicators for the construction of rice plant height and biomass models according to the inversion needs of rice plant height and biomass (

Table 2. The calculation formulae of the vegetation indices.

Spectral Indices	Calculation Formula	Reference
Normalized Difference Vegetation Index (NDVI)	NDVI = (B5 − B4)/(B5 + B4)	[34]
Ratio Vegetation Index (RVI)	RVI = B5/B4	[35]
Enhanced Vegetation Index (EVI)	EVI = 2.5 × (B5 − B4)/(B5 + 6 × B4 − 7.5 × B2 + 1)	[36]
Difference Vegetation Index (DVI)	DVI = B5 − B4	[37]
Renormalized Difference Vegetation Index (RDVI)	RDVI = (B5 − B4)/Sqrt(B5 − B4)	[38]
Green Normalized Difference Vegetation Index (GNDVI)	GNDVI = (B5 − B3)/(B5 + B3)	[39]
Optimized Soil-Adjusted Vegetation Index (OSAVI)	OSAVI = (B5 − B4)/(B5 + B4 + 0.16)	[40]
Atmospherically Resistant Vegetation Index (ARVI)	ARVI = (B5 − 2 × B4 + B2)/(B5 + 2 × B4 + B2)	[41]
Structure Insensitive Pigment Index (SIPI)	SIPI = (B5 − B2)/(B5 − B4)	[42]
Chlorophyll Index (CIgreen)	CIgreen = (B5/B2) − 1	[43]
Chlorophyll Index (CIred edge)	CIred edge = (B5/B4) − 1	[44]
Chlorophyll vegetation index (CVI)	CVI = B5 × B3/B22	[45]
Green Index (GI)	GI = B2/B3	[36]
Modified Chlorophyll Absorption Reflectance Vegetation Index (MCARI)	MCARI = (B4 − B3)−0.2 × (B4 − B2)/(B4/B3)	[46]
Modified Non-linear Vegetation Index (MNVI)	MNVI = 1.5 × (B52 − B3)/(B52 + B3 + 0.5)	[47]
Modified Simple Ratio Index (MSR)	MSR = B5/B3 − 1/Sqrt(B5/B3 + 1)	[48]
MERIS Terrestrial Chlorophyll Index (MTCI)	MTCI = (B5 − B4)/(B4 − B3)	[49]
Non-linear Vegetation Index (NLI)	NLI = (B52 − B3)/(B52 + B3)	[50]
Ratio Vegetation Index2 (RVI2)	RVI2 = RVI = B5/B2	[51]
Modified Soil-Adjusted Vegetation Index (MSAVI)	MSAVI = (2 × B4 + 1 − Sqrt(2 × (2 × B4 + 1)−8 × (b4 − b3)))/2	[51]
Second Modified Triangular Vegetation Index (MTVI)	MTVI = 1.5 × (1.2 × (NIR-Green) − 2.5 × (Red-Green))Sqrt((2 × NIR + 1)² − (6 × NIR-5 Sqrt (Red)) − 0.5)	[51]
Transformed Chlorophyll Absorption in Reflectance Index (TCARI)	TCARI = 3((B4 − B3) − 0.2 × B4/B3)	[52]

Note: B1–B5 represents blue (B1), green (B2), red (B3), red-rimmed (B4), and near-infrared (B5), respectively.

). The vegetation indices selected for this study included three main aspects: (1) eliminating soil background values: EVI, DVI, RDVI, OSAVI, MSR, and MSAVI; (2) monitoring crop growth: NDVI, RVI, GNDVI, ARVI, SIPI, NLI, RVI1,2; and (3) assessing the chlorophyll content of crops: CI-green, CI- red edge, CVI, GI, MCARI, MTCI, TCARI, MTVI. For example, NDVI is more sensitive to plant leaf surface canopy structure and is suitable for plant height monitoring. EVI is an improvement of NDVI that removes the influence of background values such as soil on plant spectra. CI-green and others can be more sensitive to changes in plant chlorophyll, and can effectively monitor crop growth at different fertility stages.

2.3.4. Dimensionality Reduction in Data

Vegetation indices are linear or non-linear combinations between different remote sensing spectral bands, and there is a high degree of linear correlation between the indices. In multiple covariance regression models, the presence of strong correlations between independent variables can lead to distorted modeling results. The correlation between the independent variables can be measured by the variance inflation factor (VIF_j), when VIF_j > 10, severe covariance is generally considered to exist between the independent variables. The VIF principle is as follows:

$$VI{F}_j=\frac{1}{1-R_j^2},j=1,\dots ,p$$

(1)

where R_j represents the negative correlation coefficient of the independent variable on the remaining independent variables for regression analysis, and p is the number of all independent variables.

Factor analysis (FN) is a generalization and development of principal component analysis. FN combines objects (variables or samples) with intricate relationships into a small number of factors, retaining the intrinsic links between the factors and the original variables while diluting the linear relationships between the factor variables. The formula for the principle of factor analysis is as follows:

$$X_1=a_{11}{F}_1+a_{12}{F}_2+\dots +a_{1m}{F}_m+{\epsilon }_1$$

(2)

$$X_2=a_{21}{F}_1+a_{22}{F}_2+\dots +a_{2m}{F}_m+{\epsilon }_2$$

(3)

$$X_p=a_{p1}{F}_1+a_{p2}{F}_2+\dots +a_{pm}{F}_m+{\epsilon }_p$$

(4)

where a_ij represents the factor loading; A = (a_ij) is the factor loading matrix; ε is the Special factor.

The factor score function is derived from the factor analysis model as follows:

$$F_j=a_{j1}{X}_1+a_{j2}{X}_2+\dots +a_j{}_p{X}_p\underset{}{{\displaystyle }}j=1,2,\dots ,m$$

(5)

2.4. Model Construction and Accuracy Evaluation

2.4.1. ELM

The Extreme Learning Machine (ELM) is an algorithm proposed by Huang et al. for solving single hidden layer feedforward neural networks. The main advantages of ELM are that the input and hidden layer weights are randomly generated rather than iteratively adjusted and the output layer weights are only calculated once by Moore Penrose [53]. As a consequence, ELM has a much shorter computation time than traditional gradient-based methods. Furthermore, ELM theory shows that almost any non-linear segmented continuous random hidden node can be used for ELM.

The principle of ELM is shown in Figure 2:

The above structure can be expressed by the formula:

$$T={\left[t_1,t_2,\dots ,t_Q\right]}_{m\times Q},t_j=\left[\begin{array}{l}{t}_{1j}\\ t_{2j}\\ \dots \\ t_{mj}\end{array}\right]={\left[\begin{array}{c}{\displaystyle \sum _{i=1}^l{\beta }_{i1}}g\left({\omega }_i{x}_j+b_i\right)\\ {\displaystyle \sum _{i=1}^l{\beta }_{i2}}g\left({\omega }_i{x}_j+b_i\right)\\ \dots \\ {\displaystyle \sum _{i=1}^l{\beta }_{im}}g\left({\omega }_i{x}_j+b_i\right)\end{array}\right]}_{m\times 1}\left(j=1,2,\dots ,Q\right)$$

(6)

Equation (8) can be rewritten as:

$$H\beta =T$$

(7)

where T is Full output; L is the number of hidden layers; g(x) is the activation functions; n is the number of input layers; w, β is the connection weights between the layers; H is the implied layer output matrix.

2.4.2. PSO-ELM

ELM reduces the training time and model structure by randomly generating the input and hidden layer weights, and only the number of hidden layer neurons and activation function type needs to be determined to build the ELM model. However, in model prediction, different parameter settings will produce different prediction results. When the deviation of the randomly generated input and hidden layer weights is 0, some of the hidden layer nodes will be invalid [54]. Particle Swarm Optimization (PSO), on the other hand, is an optimization algorithm that can be used to adjust the parameters of machine learning models and is widely used in agricultural modeling. Therefore, PSO is considered to optimize the ELM model to improve the model accuracy.

PSO is an optimization algorithm inspired by the social behavior of flocks of birds. Because the algorithm has fewer parameters to calculate, only the position and velocity of the swarm of particles need to be calculated, it has the advantage of being simple to operate and easy to perform mathematical analysis [33,55]. The velocity and position of a swarm of particles are calculated using Equation (8).

$$v_{ij}^{\left(t+1\right)}=\omega v_{ij}^t+c_{1}ran{d}_1\left(pbes{t}_{ij}{}^t-x_{ij}^t\right)+c_{2}ran{d}_2\left(gbes{t}_{ij}^t-x_{ij}^t\right)x_{ij}^{t+1}=x_{ij}^t+v_{ij}^{t+1}$$

(8)

where t is an evolutionary generation, v_ij and x_ij stand for the velocity and position of particle i on dimension j, respectively. w is the inertia weight and it is used to balance the global exploration and local exploitation. The rand represents the random function, c₁ is the personal learning factor and c₂ is the social learning factor.

PSO-ELM is realized by MATLAB programming. The PSO parameters selected are as follows: the size of the population is 50, the maximum iteration is 50, the acceleration coefficients c₁ = c₂ = 2.0, and inertia weight = 0.6.

The PSO-ELM algorithm flow is shown in Figure 3.

2.4.3. BPNN

Back propagation neural network (BPNN) is based on the structure and principle of the human brain nervous system, using the gradient descent method to gradually find out the partial derivatives of the target function to different neuron weights, and use them as weights to improve the model accuracy and make the model learning achieve the expected effect [56]. BPNN model includes three parts: input layer, hidden layer, and output layer, and its algorithm flow is shown in Figure 4.

2.4.4. Accuracy Evaluation

Coefficient of determination (R²) Root mean square error (RMSE), and performance deviation ratio (RPD) were applied to assess model performance. R² reflects the degree of fit between predicted and tested values, the closer the value is to 1, the better the model fits. RMSE reflects the deviation of the predicted value from the tested value, the smaller the value, the better the prediction of the model. The RPD is divided into Class A (>2), Class B (1.4 ≤ RPD ≤ 2), and Class C (≤1.4), with a better prediction when the RPD is >2. The coefficients are calculated as follows:

$$R^2=\frac{{\displaystyle \sum _{i=1}^n{\left(Xi-\overline{X}\right)}^2{\left(Yi-\overline{Y}\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(Xi-\overline{X}\right)}^2{\displaystyle \sum _{i=1}^n{\left(Yi-\overline{Y}\right)}^2}}}$$

(9)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(Yi-\overline{Y}\right)}^2}}{n}}$$

(10)

$$RPD=\frac{\sqrt{\frac{{\displaystyle \sum _{i=1}^n\overline{V}}{\left(Y_i\right)}^2}{n}}}{RMSE}$$

(11)

where X represents the measured value; Y represents the predicted value; $\overline{X}$, $\overline{Y}$ represents the sample mean value; n represents the number of the samples; and V represents the variance of the predicted value.

3. Results

3.1. Characteristics of Situ Data Changes

3.1.1. Statistical Analysis of Plant Height and Biomass

The distribution of plant height and biomass data over the whole growth period of rice is shown in

Table 3. Statistical analysis of data.

	Treatment	Number of Samples	Max	Min	Average	Standard Deviation
JB	JB-Hi	24	85.5	78.3	81.8	2.10
	JB-Hn	60	84.75	76.5	81.3	1.65
HF	HF-Hi	24	100	91	96.4	2.46
	HF-Hn	60	95.1	83.8	89.4	2.78
	HF-Bi	24	49.60	42.62	45.66	2.27
MM	MM-Hi	24	101.2	90.7	95.7	2.6
	MM-Hn	60	97.4	77.6	90.79	4.2
	MM-Bi	24	79.46	45.04	57.37	10.24

Note: Hi, Hn and Bi denote plant height under irrigation treatment, plant height under nitrogen application treatment and biomass under nitrogen application treatment, respectively.

. The model training and test sets were partitioned using the trian_test_split data set random partitioning method in the Python Scikit Learn library for 300 sets of sample data over the full fertility period, with 90% of the samples used for training and 10% for testing.

3.1.2. Characteristics of Plant Height Changes in Different Treatments

The variation in plant height of the 60 samples selected for the nitrogen application and 24 samples for the irrigation treatment is shown in Figure 5. Under different experimental treatments, the plant height showed an upward trend with the advancement of the growth period and entered the high-speed growth stage after the JB stage. It tended to be stable at the end of the HF stage because the rice in the JB stage—the HF stage transitioned from vegetative growth to reproductive growth, rice growth shifted to the rice ear and grain growth stage, and the rice plant height stopped growing.

From Figure 5a, it can be seen that the plant height size under the three irrigation treatments showed that W2 > W3 > W1. The plant height increased significantly during the transition from the JB stage to the HF stage, probably because rice plant height growth was mainly concentrated during the JB stage, when rice changed from nutritional to reproductive growth, and the accumulated photosynthetic products of rice were best when the field was kept in a submerged shallow water layer [57].

As can be seen from Figure 5b, under S1 (3 mm/d leakage), plant height increased with N application during the same growth period, with the maximum plant height occurring in the N5 treatment (375 kg/hm²); under S2 (5 mm/d leakage), plant height increased and then decreased with N application during the same fertility period, with the maximum plant height occurring in the N2 treatment (150 kg/hm²). The maximum plant height in the S1 treatment was greater than that in the S2 treatment when comparing longitudinally the rice plant height under different leaching conditions. This was due to the reduced loss of N fertilizer under low leaching conditions, which improved the efficiency of N fertilizer utilization and thus promoted rice growth. The results of the analysis indicated that the appropriate application of N fertilizer could significantly promote the growth of rice plant height when the irrigation water quantity was certain, which was consistent with the results of Zhu et al. [58].

3.1.3. Variation in Biomass Composition between Irrigation Treatments

Grain is formed by the gradual accumulation of biomass during the filling-maturation period, so the size of the rice spike biomass can directly affect the final yield of rice [59]. Changes in rice biomass composition under different irrigation treatments are shown in Figure 6. The spike biomass under the three irrigation treatments was W2 > W1 > W3, with the largest spike biomass in the W2 treatment, and significant differences in spike biomass among treatments (p < 0.05). In terms of the proportion of stem and leaf biomass, the distribution coefficients of stem and leaf biomass in the W3 treatment were higher than those in W1 (control) and W2, indicating that rice may transfer more photosynthetic products to tissues such as stems and leaves when subjected to flooding stress. This is because the concentration of ethylene and gibberellin in rice increased during flooding stress [60], promoting the growth of leaves, stems, etc. to obtain fresh air. Wu and others [60] found that this type of nutrient allocation can persist for a long time, resulting in a lack of carbon assimilates during the critical period of crop yield formation.

3.2. Correlation Analysis of Spectral Index with Rice Plant Height and Biomass

In order to discuss the correlation between the spectral indices and the measured plant height and biomass on the ground, a correlation model between the spectral indices and plant height and biomass at different growth stages was established, and the spectral indices whose correlation reached a significant level at each stage were preferably selected, as shown in

Table 4. Height-spectral index correlation analysis.

Spectral index	Correlation Coefficient (JB-H)	Spectral Index	Correlation Coefficient (HF-H)	Spectral Index	Correlation Coefficient (MM-H)
Band3	−0.256 **	TCARI	0.694 ***	Band2	0.494 ***
Clred edge	0.328 ***	SIPI	−0.586 ***	Band3	0.442 ***
Clgreen	0.216 **	RVI2	−0.633 ***	Band4	0.537 ***
Band2	−0.254 **	MTVI	0.575 ***	Band5	0.502 ***
Band4	−0.301	MCARI	0.504 ***	DVI	0.407 ***
ARVI	0.225 **	Band4	0.691 ***	MACRI	0.503 ***
MSR	0.266 **	GNDVI	0.732 ***	MNVI	0.344 ***
TCARI	−0.469 ***	ARVI	0.512 ***	MSAVI	0.324 ***
RVI2	0.216 **	Clgreen	−0.633 ***	MTVI	0.369 ***
NDVI	0.238 **	Clred edge	−0.48 ***	NLI	0.346 ***
MCARI	−0.373 ***	EVI	0.454 ***	TCARI	0.46 ***

Note: ***, ** represent 1%, 5% level of significance, respectively.

and

Table 5. Biomass-spectrum index correlation analysis.

Spectral Index	Correlation Coefficient (HF-B)	Spectral Index	Correlation Coefficient (MM-B)
Band5	0.559 ***	Band3	0.389 *
TCARI	0.568 ***	RVI1	−0.424 **
Band2	0.351 *	TCARI	0.565 ***
Band4	0.616 ***	MCARI	0.649 ***
RDVI	0.434 **	MTCI	−0.55 ***
EVI	0.458 **	MSR	−0.392 *
MCARI	0.491	GI	−0.439 *
DVI	0.564	CIred edge	−0.519 ***
MSAVI	0.385 *	Band4	0.886 ***
MNVI	0.504 **	Band5	0.433 **

Note: ***, **, * represent 1%, 5%, 10% level of significance, respectively.

.

3.3. Factor Analysis Dimensionality Reduction

The VIF between spectral indexes for the six treatments in this experiment is shown in

Table 6. Spectral index VIF statistics.

	JB	HF	MM
Irrigation treatment	>189	>95	>75
Nitrogen application	>15	>10	>140

. The results showed that there was a serious problem of multicollinearity between the spectral indexes.

According to Table 6, it can be seen that there is a strong multicollinearity between the spectral indices within different fertility periods. In this paper, factor analysis was used to extract potential co-linearity factors to solve the co-linearity problem.

Factor analysis was conducted on the preferred spectral variables (Table 4 and Table 5). The Kaiser-Meyer-Olkin (KMO) was 0.6 to 0.8 and the Bartlett’s sphericity test had a significance p-value of 0.000 ***. Bartlett’s sphericity test and factor analysis were valid.

3.4. Inverse Results of Rice Plant Height and Biomass from Different Models

The sensitive bands and vegetation indices after factor analysis were used to build ELM models for plant height and biomass, BP neural network and PSO-ELM models, respectively, in different fertility periods, and the inversion results are shown in

Table 7. Plant height and biomass inversion model results.

Model	Inverse Indicator	Training Set		Testing Set		RPD
		Rc2	RMSEc	Rv2	RMSEv
ELM	JB-H	0.54	0.58	0.68	0.69	1.48
	HF-H	0.81	0.07	0.82	0.01	2.36
	MM-H	0.89	0.06	0.83	0.09	2.39
	HF-B	0.69	3.55	0.75	1.76	2.02
	MM-B	0.86	3.91	0.84	2.47	2.50
BPNN	JB-H	0.26	0.53	0.27	0.85	1.17
	HF-H	0.70	0.04	0.76	0.77	2.03
	MM-H	0.72	0.03	0.79	0.52	2.17
	HF-B	0.64	0.27	0.71	2.06	1.59
	MM-B	0.65	2.42	0.61	6.29	1.59
PSO-ELM	JB-H	0.43	1.30	0.47	1.95	1.37
	HF-H	0.74	2.12	0.83	1.99	2.46
	MM-H	0.75	2.17	0.79	2.83	2.19
	HF-B	0.82	1.08	0.82	1.08	2.33
	MM-B	0.97	1.61	0.94	6.53	4.24

.

In terms of overall model quality, the three inversion models of rice plant height and biomass based on multispectral remote sensing proposed in this paper all achieved better inversion results, with most of the R² above 0.6 and the RPD above 1.7, which have high accuracy and good stability. However, the BPNN model was less effective in the inversion of plant height at the nodulation stage: R_v² was 0.269 and RMSE_v was 0.95.

Among the three plant height models for the same growth period, the JB-H optimal inversion model is the ELM model with R_v² of 0.68, RMSE_v of 0.69, and RPD of 1.48; the HF-H optimal inversion model is the ELM model with R_v² of 0.82, RMSE of 0.01, and RPD of 2.36; the MM-H optimal inversion model is the ELM model with R_v² of 0.83, RMSE_v of 0.09 and RPD of 2.39. The best inversion model for MM-H is the ELM model with R_v² of 0.83, RMSE_v of 0.09, and RPD of 2.39. Among the biomass inversion models, the best inversion models for HF-B and MM-B are both PSO-ELM models with R_v² of 0.82 and 0.94, RMSE_v of 1.08 and 6.53, and RPD of 2.33 and 4.24, respectively.

Since the ELM model has a special mechanism for selecting the hidden layer node parameters, which requires a large number of hidden layer nodes to achieve the target accuracy and has poor generalization capability, the PSO algorithm automatically learns the number of hidden layer nodes from the training data by optimizing the input layer weights and hidden layer deviations of the ELM for the purpose of optimization. The PSO evolutionary curve was shown in Figure 7, which showed that the PSO-ELM algorithm was highly optimized for the biomass inversion model, with the R² and RPD of the HF-B and MM-B training sets improving by 0.061, 0.1024, 0.7363, and 1.742, respectively. However, the PSO-ELM was less optimized for the plant height inversion model, probably because the plant height model had 84 sample data in different growth periods compared to 24 sample data for biomass in the same period, implying that the particle swarm optimized PSO-ELM model reduced the complexity of the model and was more accurate when dealing with fewer sample data.

Of the three models introduced in this experiment, ELM and PSO-ELM had the highest inversion accuracies for rice plant height and biomass, respectively, as shown in Table 7. The two models trained from the training set data still had high accuracy (R² > 0.68) after the test set data were imported and there was not much difference between the accuracy of the training set and that of the validation set, which indicated that their stability and adaptability were high. At the same time, this experiment obtained a large range of rice physiological data through different treatments, which expanded the data diversity and could be adapted to rice inversion under different planting conditions.

3.5. Validation of Rice Plant Height and Biomass Accuracy for Different Models

The accuracy of the inverse models of plant height and biomass constructed for different fertilization stages of rice was verified separately, as shown in Figure 8. Except for the poor inversion of JB-H by BPNN, the inversion R² of all three groups of models was higher than 0.4 and could reach up to 0.9 in different growth periods, showing a better inversion effect. The BPNN model was the best for MM-H with an R² of 0.71, while the ELM model was the best for HF-H, MM-H, and MM-B with R² of 0.81, 0.83, and 0.85, respectively. In summary, ELM and PSO-ELM performed best in the inversion of rice plant height and biomass. The BP algorithm is a gradient descent algorithm that attempts to improve the performance of a neural network by reducing the total error by changing the weights of the neural network along the gradient direction. In conclusion, the BP model is a process of modifying the weights and faces difficulties in manually adjusting the control parameters (e.g., learning rate and learning period) or local minima, which results in a weaker inversion accuracy compared to ELM and PSO-ELM.

4. Discussion

Plant height and biomass are important parameters for crop growth monitoring and yield estimation [61,62]. Rapid and accurate acquisition of crop height and biomass information through UAV multispectral remote sensing is the main direction for the future development of smart agriculture and precision farming. In this paper, we developed ELM, BPNN, and PSO-ELM rice plant height and biomass inversion models from multispectral remote sensing data under different irrigation and nitrogen application conditions, and evaluated the accuracy of model estimation.

The results of the experiment showed that the percentage of dry matter distributed among the various organs of rice under flooding stress differed significantly (p > 0.05) compared to other irrigation treatments, with the percentage of biomass of leaves and stems increasing and the percentage of biomass of spikes decreasing [63], which may be the result of rice adapting to flooded conditions. Nitrogen fertilizer application has the effect of increasing rice plant height only within a certain dosage range, and excessive application of nitrogen will result in an insignificant increase or even a decrease in plant height.

Firstly, the correlation between spectral index, plant height and biomass was analyzed and the spectral index introduced into the model was preferred based on significant correlation (p < 0.01) and correlation coefficient (R > 0.4). However, the preferred vegetation indices such as TCARI and RDVI are linear combinations of R, G, B, RE, and NIR, with extremely strong covariance between the respective variables, and higher multicollinearity may lead to large variance in the least squares estimates of the beta coefficients in the regression equation and affect the final results [56,64]. This study extracted potential common factors through factor analysis and transformed the input variables from spectral indices to a total of 10 indicators from T1-T10, diluting the autocorrelation between factors, avoiding to a certain extent the interaction between independent variables and improving the stability of the model.

Rice undergoes three stages: the jointing and booting stage, the heading and flowering stage, and the milky ripening stage. Under the ELM algorithm, the plant height and biomass inversion accuracies increased as the rice gradually matured, and the highest inversion accuracy among the three stages was at the milky ripening stage. When rice is at the jointing and booting stage, the plant is shorter (this is why this paper did not study the biomass inversion at the jointing and booting stage of rice), and when extracting the region of interest (ROI), the leaf blades are thin and narrow and difficult to recognize, and the collected spectral data are easily affected by the soil and other influences, which leads to low accuracy.

Previous studies have shown that ELM usually shows high accuracy in inversion studies of crop growth and soil physical and chemical properties. The only parameter that needs to be defined for ELM is the number of nodes in the hidden layer. Unlike traditional learning algorithms for neural networks (e.g., BP algorithm), which may face difficulties in manually adjusting the control parameters, learning rate, and learning period. ELM is fully automated and does not require iterative adjustments, and theoretically does not require user intervention, resulting in a relatively higher accuracy of the model. The ELM, BPNN, and PSO-ELM models were compared and analyzed for plant height and biomass inversion at different growth stages, and it was found that the ELM model was the most suitable for rice plant height inversion, which was consistent with previous studies. For rice biomass inversion, the optimized PSO-ELM model inversion was better than the other two models, with R_v² above 0.81 and 2.33. However, the two inversion models were less effective in inversion of rice plant height at the jointing and booting stage, probably due to the low density of rice at the plucking stage, where the spectral information of the rice canopy easily overlapped with the soil, resulting in the inability of the multispectral camera to obtain complete information.

Comparing the PSO-optimized model with the ELM model, PSO-ELM improved the model accuracy in the rice biomass inversion model. The R² of HF-B and MM-B improved by 0.06 and 0.10, respectively, but in the strain height inversion model, PSO was less optimized compared to the ELM model, and the R² of JB-H and MM-H decreased by 0.22 and 0.03, respectively. This difference may be due to the random generation of hidden layers in the ELM model, which usually requires more hidden neurons than traditional neural networks to achieve matching performance, but a principle of the ELM model is that the number of hidden layer neurons is not greater than the number of input samples, which is less in the biomass inversion model and affects the accuracy of the model [53]. In contrast, PSO-ELM achieves optimization by automatically learning from the training data to obtain the number of hidden layer nodes.

As a method to estimate rice plant height and biomass by modeling, the data set is an important factor in determining the results of the model. It has been demonstrated that the performance of machine learning models is highly dependent on the amount of data. Zhu et al. [12] obtained 150 training samples and 50 test samples by planting Maize in the field and achieved relatively stable model results. In this experiment, due to the experimental site, rice was cultivated by planting pots. Compared with the field experiment, the experimental gradient and the number of samples obtained were small. In terms of data significance and model stability, small sample data may have certain defects. This is where the follow-up research needs to focus on improving.

5. Conclusions

In this study, we analyzed the changes in rice plant height and biomass under different irrigation and nitrogen application treatments, and collected multispectral remote sensing data during the whole growth period of rice using UAV multispectral technology, combined with three machine learning methods to construct an inverse model of rice plant height and biomass, and obtained the following conclusions.

Rice plant height was highest at all growth stages under different irrigation treatments W2 (water level regulation 3–5 cm); among the different N application treatments, rice plant height increased with increasing N application under low percolation conditions (3 mm/d). Under high percolation conditions (5 mm/d) rice plant height increased and then decreased with N fertilizer application at different growth stages, with the maximum plant height occurring at N2 (150 kg/hm²). Flooding stress affected the principle of biomass (dry matter mass) distribution in various tissues of rice, and the dry matter mass of the spike decreased with stress, eventually leading to a reduction in rice yield.

The results of this study showed that the best model for plant height inversion during the whole rice reproductive period was the ELM model, with R_v² of 0.684 and RPD of 1.48 in the inversion of plant height at the jointing and booting stage, R_v² of 0.82 and RPD of 2.36 in the inversion of plant height at the heading and flowering stage, and R_v² of 0.75 and RPD of 2.02 in the inversion of plant height at the milk ripening stage. The best model for biomass inversion was the PSO-ELM model, with an R_v² of 0.82 and RPD of 2.33 in the inversion of heading and flowering biomass and an R_v² of 0.94 and RPD of 4.24 in the inversion of milk ripening biomass.