UAV-Based Biomass Estimation for Rice-Combining Spectral, TIN-Based Structural and Meteorological Features

Abstract

Accurate estimation of above ground biomass (AGB) is very important for crop growth monitoring. The objective of this study was to estimate rice biomass by utilizing structural and meteorological features with widely used spectral features. Structural features were derived from the triangulated irregular network (TIN), which was directly built from structure from motion (SfM) point clouds. Growing degree days (GDD) was used as the meteorological feature. Three models were used to estimate rice AGB, including the simple linear regression (SLR) model, simple exponential regression (SER) model, and machine learning model (random forest). Compared to models that do not use structural and meteorological features (NDRE, R² = 0.64, RMSE = 286.79 g/m², MAE = 236.49 g/m²), models that include such features obtained better estimation accuracy (NDRE*Hcv/GDD, R² = 0.86, RMSE = 178.37 g/m², MAE = 127.34 g/m²). This study suggests that the estimation accuracy of rice biomass can benefit from the utilization of structural and meteorological features.

1. Introduction

Rice is an important cereal grain that provides a stable source of food for more than half of the world’s population [1]. Accurate estimation of rice parameters, like above ground biomass (AGB), are very important for rice growth monitoring, the assessment of rice nutrition status, and the prediction of rice yield [2]. Traditional methods to estimate biomass are both labor-intensive and time-consuming, and are difficult to apply over large areas [3]. Remote sensing is a non-contact and non-destructive measuring method that can acquire both spectral and structural properties of the target at different spatial and temporal scales. These characteristics make remote sensing an optimal method for large area biomass estimation [3].

Currently, the spectrum used for the estimation of AGB are sourced from either ground-based spectral data or satellite-based spectral data [3,4,5,6]. Obtaining ground-based spectral data is labor-intensive and time-consuming, which can cause scalability issues. On the other hand, satellite-based data may not provide sufficient data resolution for applications in precision agriculture. LiDAR (light detection and ranging) is another technique and has been used to estimate biomass. The high spatial resolution of LiDAR can provide a precise crop surface model (CSM), which further allows the biomass to be estimated using plant height [7,8]. However, a high computational resource is required for processing LiDAR data, which would not be ideal for applications that cover a large area [9].

The emergence and development of unmanned aerial vehicles (UAV) and lightweight sensors has filled the existing gap between satellite and terrestrial remote sensing, which offers a promising tool for precision agriculture. For example, Jin [10] estimated wheat plant density using red-green-blue (RGB) images acquired from an UAV at low altitude (3 m–7 m); Pádua [11] monitored the growth of grapevines with multi-temporal UAV-based RGB images; Zhou [12] used two cameras with different spectral ranges mounted on an UAV to acquire multi-temporal images to predict rice yield; Marcial-Pablo [13] used RGB and multispectral UAV images to estimate vegetation fraction of maize.

Three-dimensional (3D) point clouds generated from a series of overlapped 2D images acquired by an UAV using the Structure from motion (SfM) algorithm offer new options for the acquisition of crop surfaces. SfM is a computer vision technology that generates 3D geometry by repetitive bundle adjustment of the multiple unordered overlapped images and image matching techniques, like the scale-invariant feature transform (SIFT) [14,15]. Crop surface models (CSMs) derived from 3D point clouds contain crop canopy vertical distribution information, which can be used for crop monitoring, e.g., plant height measurement [16], biomass estimation [17,18,19,20], and yield prediction [21]. In addition to CSMs, RGB images, and multi- and hyper-spectral images acquired for UAV were combined with CSMs to estimate biomass [22,23,24,25,26,27].

Although CSMs derived from 3D point clouds were used in the research mentioned above, they focused mostly only on the plant heights derived from gridded CSMs. Triangulated irregular network (TIN) [28] is a kind of digital elevation model (DEM) derived from point clouds. TIN represents the surface with a series of continuous, non-overlapping, irregular triangles using the Delaunay triangulation algorithm. TIN has been used for studying the catchment hydrologic response [29], surface flow routing [30], and crown volume extraction [31]. Compared with gridded DEM, TIN is able to show surface structure details more accurately and more efficiently without the interpolation process [32,33].

Temperature is one of the most important driving variables for the simulation of crop growth and development [34]. Thermal time, mainly based on growing degree days (GDD), was introduced as the meteorological information in this study. GDD is calculated by subtracting a base temperature, below which no significant crop development is expected, from the daily mean temperature. Previous studies used GDD to predict winter wheat yield [35], estimate the growth stages of rice [36], map early season large-area winter crops [37], and estimate biomass accumulation of maize kernels [38].

The main goal of our study is to evaluate the potential of biomass estimation using multitemporal spectral and structural information obtained with an UAV. More specifically, our study aims to: (1) Evaluate the performance of AGB estimations based on multispectral features and structure features; (2) establish and evaluate models based on a combination of multispectral and structure features; and (3) establish and evaluate biomass estimation models based on a combination of spectral, structural, and meteorological features.

2. Materials and Methods

2.1. Study Area

The study was conducted on 150 m by 70 m experimental fields in the Hainan Hybrid Rice Experimental Base of the State Key Laboratory of Hybrid Rice (Wuhan University), located in Lingshui, Hainan, China (18°32′2″ N, 110°2′34″ E) (Figure 1). Lingshui County is located in the southeast of Hainan Island, and the climate of Lingshui County is a tropical monsoon climate. The lowest temperature is 11 °C in February and the highest is 35 °C in July. The relatively high temperature in winter compared with the mainland makes Lingshui an ideal place for crop winter breeding from November the previous year to May the next year.

The study field was flattened and separated into eight districts by field ridges before transplanting. Two hundred and fifty-five different hybrid rice species were planted in total. The rice seeds were sowed on December 10, 2017 and transplanted on January 8, 2018. To ensure healthy growth, fertilizers were applied by the breeding scientists based on the growth stage of the rice. Insecticides and herbicides were applied to control pests and weeds as well. Districts 5 and 6 were selected as our research area, which have 42 hybrid species in total. These 42 hybrid species were recommended by the breeding scientists as our research species because these species are some of the most typical hybrid species of China and they have a wide variety of height, biomass, yield, etc. For these 42 hybrid species plots, 12 rows were planted and each row had double lines. The distance between the rows was 33 cm and the distance within rows was 20 cm. The plant density was approximately 23 plants per square meter. Each hybrid specie plot was divided into a 7 by 7 m measuring area for plant height, spectral, and structural information extraction and a 2 by 7 m sampling area for biomass estimation. Altogether, 11 ground control points (GCPs) made of 40 cm by 40 cm polyvinyl chloride (PVC) board and 1 m long wooden poles were distributed evenly across the field. Their coordinates were measured using the Hi-Target iRTK2 system (Hi-Target Surveying Instrument Co., Ltd., Guangzhou, China). The estimated accuracy of the GCPs was 2 cm in horizontal coordinates and 2 cm in height. The detailed plot arrangement is shown in Figure 1.

2.2. Field Data

Destructive samplings were conducted six times for each plot during the whole growing period (

Table 1. Data acquisition time, corresponding growth stage, day after sowing (DaS), and growing degree days (GDD). All the dates were in 2018.

Growth Stage	Biomass Sampling	Multispectral UAV	RGB UAV	Field Plant Height	DaS (d)	GDD (°C·d)
Late Tilling	Feb 25th	Feb 26th	Feb 28th	Mar 1st	77	489.5
Early Jointing	Mar 9th	Mar 11th	Mar 10th	-	89	665.5
Late Jointing	Mar 19th	Mar 18th	Mar 19th	-	99	798.5
Heading	Mar 31st	Apr 1st	Mar 31st	Mar 31st	111	957.5
Filling	Apr 17th	Apr 15th	Apr 15th	-	128	1212.5
Milking	Apr 28th	-	Apr 28th	Apr 24th	139	1400.5

Note: ‘-’ means lack of data. The multispectral UAV platform broke down during the last flight, causing the lack of data of the last sampling.

). For each sampling, three hills of rice plants were sampled near the center of the sampling area for each plot. After sampling, the plants’ roots were clipped, and the remaining part of the plant was cleaned and separated into the stems, leaves, and ears (after heading stage). All the samples were dried for 30 min at 105 °C and later at 80 °C until their weights remained unchanged, and then were weighed to calculated dry AGB.

Field plant height was measured on March 1, 2018, March 31, 2018, and April 24, 2018 for 42 selected plots. For each plot, the average plant height was measured with a T-shaped ruler. Three plants near the center of the measurement area were randomly selected and the height was measured as the height from the ground to the highest part of the rice plant.

GDD of each sampling date was calculated using Equation from [39]. In Equation, D_i is the ith day from D₁, which is the day of transplanting. T_MAX is the daily maximum temperature, T_MIN is the daily minimum temperature, and T_BASE is the temperature below which the process of growth does not progress. In this equation, if [(T_MAX + T_MIN)/2] < T_BASE, then [(T_MAX + T_MIN)/2] = T_BASE. T_BASE varies among different crops. For paddy rice, T_BASE is 10 °C [40]. We acquired the daily temperature information of Lingshui County from the Hainan Meteorology Administration.

$$GD{D}_{D_i}={\displaystyle \sum }_{D_1}^{D_i}\left\{\left[\frac{\left(T_{MAX}+T_{MIN}\right)}{\mathbf{2}}\right]-T_{BASE}\right\}$$

(1)

2.3. UAV Remote Sensing Data

The Mini-MCA 12 multispectral camera (Tetracam Inc., Chatsworth, CA, USA) mounted on a DJI S1000 octocopter (SZ DJI Technology Co., Ltd., Shenzhen, China) was used to obtain multispectral images for the study area. The Mini-MCA had 12 bands with different band widths (

Table 2. Center wavelength and full width at half-maximum (FWHM) bandwidth of each spectral band of the Mini-MCA 12 multispectral camera.

Spectral Band	Center Wavelength (nm)	Bandwidth FWHM (nm)	Spectral Band	Center Wavelength (nm)	Bandwidth FWHM (nm)
Blue	490	10		700	10
	520	10	Rededge	720	10
Green	550	10	NIR	800	10
	570	10		850	10
Red	670	10		900	20
	680	10		950	40

). These bands were selected since they are commonly used for estimating vegetation parameters [41,42,43]. The cameras were co-registered in the laboratory before flight using a camera distortion correction model so that the corresponding pixels of each lens were spatially overlapped in the same focal plane [44,45]. A gimbal stable platform was equipped on the UAV to ensure a nadir view, which minimized the influence of UAV fluctuations during image acquisition. Each flight was carried out on a sunny day between 11:00 and 13:00 local time when the changes in light intensity were small. The flight height was 120 m and the ground sampling distance (GSD) was 6.5 cm/pixel so that the study area could be covered in a single shot.

RGB images were acquired using a DJI Phantom 4 Professional (SZ DJI Technology Co., Ltd., Shenzhen, China). The camera mounted on the UAV was the DJI FC6310. The camera has a field of view (FOV) of 84° and an image resolution of 5472 × 3648, approximately 20 million pixels. The flight missions were planned using DJI Ground Station Pro (v. 1.8.5.). The flight path remained unchanged for all flights. The camera’s direction was set perpendicular to the ground to acquire nadir images.

Due to the physiological characteristics of paddy rice, the rice field was constantly covered by water and the maximum water depth can be more than 5 cm in the early stages, which can cause sun glint in the RGB images (Figure 2), affecting the generation of point clouds and digital orthophoto maps (DOMs). Moreover, the hotspot effect [46,47,48] that exists in the vegetation canopy (Figure 2) has a strong influence on the quality of point clouds and DOMs as well [48]. Both the sun glint and the hotspot effect will appear when the sun elevation angle is greater than 90 minus 1/2 the FOV of the camera, so it is important to analyze the solar elevation angle before the planning of flights to avoid the influences caused by these two effects.

The solar elevation angle, α, was calculated using Equation:

$$\sin \alpha = \sin \delta \ast \sin \phi +\cos \delta \ast \cos \phi \ast \cos \omega$$

(2)

where δ is the solar declination, φ is the latitude, and ω is the solar hour angle. Here, δ is related to the flight date and ω is related to the flight date and time. For the detailed formula, please refer to [49].

Given the flight date, flight time, and the field’s longitude and latitude, the solar elevation angle was calculated. The open source solar library, Solar Position Algorithm (SPA), was used here (https://rredc.nrel.gov/solar/codesandalgorithms/spa/). After the calculation and field test, we set the flight time of the RGB camera between 9 a.m. and 10 a.m., when the solar angle would not cause the sun glint and hotspot effects and the solar illumination was ensured.

Due to the changeable costal weather, the flights were conducted one or two days after the ground sampling. Unfortunately, the DJI S1000 UAV platform broke down during the last flight, which led to the lack of multispectral images for the last sampling.

2.4. Data Processing

2.4.1. Radiometric Processing

The image digital numbers (DN) were transformed to surface reflectance using the empirical line method (ELM) [50]. Six calibration blankets (nominal reflectance was 0.03, 0.12, 0.24, 0.36, 0.56, and 0.8, respectively) made with highly durable woven polyester fabric were placed in the cameras’ FOV as radiometric reference targets. Since a linear relationship was assumed between the surface reflectance and DNs, the canopy surface can be calculated as [51,52]:

$$\begin{gathered} {\rho }_{\lambda }=D{N}_{\lambda }\times G_{\lambda }+B_{\lambda } \\ \left(\lambda =490, 520, 550, 570, 670, 680, 700, 720, 800, 850, 900, \mathrm{and} 950 \mathrm{nm}\right) \end{gathered}$$

(3)

$$\left(\begin{array}{c}0.03\\ 0.12\\ 0.24\\ 0.36\\ 0.56\\ 0.8\end{array}\right)=\left(\begin{array}{c}D{N}_{0.03}\\ D{N}_{0.12}\\ D{N}_{0.24}\\ D{N}_{0.36}\\ D{N}_{0.56}\\ D{N}_{0.8}\end{array}\right)\times G_{\lambda }+B_{\lambda }$$

(4)

where ρ_λ and DN_λ are the surface reflectance and digital number, respectively, of a given pixel at the wavelength, λ; G_λ and B_λ are the gains and bias of the camera at different wavelengths, respectively. For each wavelength, G_λ and B_λ are calculated using the least-square method from ρ and DN values of six calibration targets (Equation).

For each of the 42 plots, a region-of-interest was selected (around 4900 pixels) in the measuring area. The plot-level reflectance was calculated as the average reflectance over the pixels in the region-of-interest. Vegetation indices (VIs) were then retrieved from the plot-level reflectance (

Table 3. Vegetation indices used in this study.

Vegetation Indices	Formula	References
NDVI	$\left({\rho }_{800}-{\rho }_{670}\right)/\left({\rho }_{800}+{\rho }_{670}\right)$	[53]
GNDVI	$\left({\rho }_{800}-{\rho }_{550}\right)/\left({\rho }_{800}+{\rho }_{550}\right)$	[54]
NDRE	$\left({\rho }_{800}-{\rho }_{720}\right)/\left({\rho }_{800}+{\rho }_{720}\right)$	[55]
SR	${\rho }_{800}/{\rho }_{700}$	[56]
TVI	$0.5\ast \left[120\ast \left({\rho }_{800}-{\rho }_{550}\right)-200\ast \left({\rho }_{670}-{\rho }_{550}\right)\right]$	[57]
EVI2	$2.5\ast \left({\rho }_{800}-{\rho }_{670}\right)/\left({\rho }_{800}+6\ast {\rho }_{670}-7.5\ast {\rho }_{490}+1\right)$	[58]
RDVI	$\left({\rho }_{800}-{\rho }_{670}\right)/\sqrt{\left({\rho }_{800}+{\rho }_{670}\right)}$	[59]
CIrededge	${\rho }_{800}/{\rho }_{720}-1$	[60]
CIgreen	${\rho }_{800}/{\rho }_{550}-1$	[60]
LCI	$\left({\rho }_{850}-{\rho }_{710}\right)/\left({\rho }_{850}-{\rho }_{680}\right)$	[61]
MCARI	$\left[\left({\rho }_{700}-{\rho }_{670}\right)-0.2\ast \left({\rho }_{700}-{\rho }_{550}\right)\right]/\left({\rho }_{700}/{\rho }_{670}\right)$	[62]
Band Reflectance	Reflectance of 12 bands of MCA

). The radiometric processing was done using the ENVI software (Exelis Visual Information Solutions, Boulder, Colorado, USA).

2.4.2. Point Cloud Processing

After every flight, the SfM algorithm was employed to generate point clouds to acquire TIN-based structural features for each plot.

Point cloud data used in this study were produced using the UAV RGB images and the Pix4D mapper software (Pix4D SA, Lausanne, Switzerland). First, the images captured from the same flight were imported into the software and all the GCPs were manually marked. The images were then aligned, mosaicked, and geo-referenced by feature matching [14] and structure from motion (SfM) [15] algorithms embedded in the software. Then, the SfM algorithm reconstructed the 3D scene and elements of the interior and exterior orientation of each image, using geo-reference information, camera coordinates, and GCP coordinates. After this, a dense point cloud was generated from a huge amount of common features in the images. The digital surface model (DSM) was then derived using the dense point cloud. For each flight, 4 GCP points were randomly selected as check points to evaluate the accuracy of the DSMs. The accuracy of the DSMs is shown in

Table 4. The accuracy of digital surface models (DSMs) of each flight.

Flight Date	Error X (cm)	Error Y (cm)	Error Z (cm)
Jan 8th	2.17	1.50	1.63
Feb 28th	2.28	1.49	1.10
Mar 10th	2.43	2.69	1.99
Mar 19th	2.65	3.20	2.23
Mar 31st	3.06	2.55	1.74
Apr 15th	3.23	2.36	1.67
Apr 28th	3.48	2.60	2.00

Error X, Y, Z represent the root mean square value (RMSE) of the X, Y, Z coordinates.

. The DSMs were not used in this study except for the DSM generated from the extra flight on January 8, 2018, which was used as the base DEM for the calculation of the TIN-based structural features.

For a 3D point cloud, P(x,y,z), the true point heights were derived by subtracting the corresponding height of the base DEM from original point heights. The corresponding height of the DEM, H_DEM(x,y), was calculated using the bilinear interpolation method with the nearest four points of the DEM. The true point height, H_t(x,y), can be obtained as in Equation. H_p(x,y) is the z value of 3D point cloud, P(x,y,z).

$$H_t\left(x,y\right)=H_p\left(x,y\right)-H_{DEM}\left(x,y\right)$$

(5)

Once the true point height was derived, the original point height was replaced by the true point height, resulting in a true point cloud representing the canopy surface. Point clouds of the 42 plots were manually extracted for further analysis. Each plot point cloud was filtered using the statistical outlier removal (SOR) algorithm to remove possible outliers. After that, TINs were created using the Delaunay triangulation algorithm. For each plot, a total of 17 TIN-derived 3D features were retrieved from each plot-based TIN (

Table 5. 3D features derived from triangulated irregular network (TIN).

3D Features	Description
Hmax	Maximum rice height
Hmean	Mean rice height
Hstd	Standard deviation of rice height
Hcv	Coefficient of variation of rice height
Hen	Entropy of rice height
Hre	Relief of rice height
Hp50	50th percentile of rice height
Hp75	75th percentile of rice height
Hp90	90th percentile of rice height
Hp95	95th percentile of rice height
Hp99	99th percentile of rice height
Amean	Mean area of TIN triangles
APmean	Mean area of TIN triangles’ projection on the horizontal plane
Aratio	Amean/APmean
Smean	Mean slope of TIN triangles
Sstd	Standard deviation of slope of TIN triangles
CSmean	Cosine of slope of TIN triangles

).

2.5. Statistical Methods

The simple linear regression (SLR) model and simple exponential regression (SER) model were used to estimate AGB. The SLR model was selected since it is the most common model for crop parameter estimation. Previous studies indicated that AGB has an exponential relationship with vegetable indices and indices related with plant height [2,17,27], so the SER model was selected as well. Univariate and multivariate models were used in the SLR and SER models. The data used for the estimation include the band reflectance, vegetation indices (VIs), TIN-based structural features (TINs), and meteorological features (GDD). For the multivariate models, two combination forms were considered: Multispectral and TIN-based structural features; multispectral, TIN-based structural features, and meteorological features (GDD).

Because SLR and SER cannot solve multiple correlation problems, only a few indices were used in SLR and SER. Machine learning algorithms, such as random forest (RF), are known be efficient algorithms for dealing with multiple correlation problems. RF is a data analysis and statistical method proposed by Breiman and Culter for classification and regression [63]. RF has good tolerance to outliers and noise, makes full use of input data and avoids over fitting when the number of regression tree increases.

For each measurement, we acquired 42 datasets, with 210 datasets in total. The models were tested using the 5-fold cross validation. The model performance was evaluated using the coefficient of determination (R²), root mean square value (RMSE), and mean absolute error (MAE). In this study, MATLAB R2018a (MathWorks, Inc., Natick, Massachusetts, USA) and Python language was used to process point clouds and perform the data analysis. A flowchart of the study is provided in Figure 3.

3. Results

3.1. Plant Height Estimation

For plant height estimation, we evaluated the correlationship between Hmean, Hp99, Hp95, Hp90, Hp75, Hp50 (Table 5), and the manually measured plant height for the entire 126 field measured plant heights. All the mentioned height features had strong linear correlationships (R² > 0.85) with the manually measured plant height. The mean UAV plant height had the strongest correlation, with R² = 0.9058 and RMSE = 0.2658 m (Figure 4), although the mean UAV plant height was 0.264 m smaller than the manually measured height. The coefficient of variation (CV) for the mean UAV plant height was 17.53% and for the manually measured plant height was 12.10%.

3.2. Univariate Biomass Estimation

Correlation coefficients between features and AGB or ln (AGB) (natural log-transformed AGB) were calculated to select proper features for further analysis and to decide which model is best for AGB estimation; here, AGB for the SLR model and ln (AGB) for the SER model. Features that had higher correlation coefficients with either AGB or ln (AGB) were selected to build univariate and multivariate AGB estimation models.

3.2.1. Multispectral Biomass Estimation

Most of the specific bands and vegetation indices had good correlations with AGB except for the near infrared bands (ρ800, ρ850, ρ900, and ρ950), TVI, EVI2, and RDVI (Figure 5). Specific bands were positively related to AGB while most of the VIs except for MCARI were negatively related to AGB. The difference between the correlation coefficients of SLR and SER was small for correlated specific bands (Δr < 0.02), which indicated that both models could be used for AGB estimation for specific bands. For correlated VIs, the correlation coefficients of SLR were higher than SER. Therefore, SLR was selected to build models for specific bands and VIs. Variables with correlation coefficients higher than 0.6 or lower than –0.6 (ρ520, ρ570, ρ670, ρ680, ρ700, ρ720, NDVI, GNDVI, NDRE, SR, CIrededge, CIgreen, LCI, and MCARI) (Figure 5) were selected to build univariate multispectral AGB estimation models. Their performance was moderate. NDRE had the best estimation accuracy (R² = 0.64, RMSE = 286.79 g/m², MAE = 236.49 g/m²) while CIrededge and ρ700 also had better estimation accuracy (R² = 0.63, RMSE = 288.95 g/m², MAE = 243.77 g/m², R² = 0.61, RMSE = 297.33 g/m², MAE = 231.38 g/m²). The statistical parameters of the multispectral regression models are shown in

Table 6. Statistical parameters of regression models of multispectral features and above ground biomass (AGB) (SLR). The best model is bolded.

Model	Variables	R-Squared (R2)	RMSE (g/m2)	MAE (g/m2)
SLR	ρ520	0.54	322.14	257.26
	ρ570	0.55	318.23	243.49
	ρ670	0.37	375.99	293.50
	ρ680	0.34	385.73	302.29
	ρ700	0.61	297.33	231.38
	ρ720	0.46	346.81	277.22
	NDVI	0.41	394.27	282.82
	GNDVI	0.52	329.06	273.81
	NDRE	0.64	286.79	236.49
	SR	0.58	307.31	247.25
	CIrededge	0.63	288.95	243.77
	CIgreen	0.39	369.98	304.66
	LCI	0.54	322.55	260.72
	MCARI	0.42	360.48	281.88

.

3.2.2. Structural Biomass Estimation

Most TIN-based structural features were significantly correlated with AGB except Hre, Astd, Acv, and Sstd (Figure 6). Variables representing plant height had positive relationships with AGB while calculated plant height indices had negative relationships with AGB. Among other TIN-based structural features, Amean, APmean, Astd, Sstd, Scv, and CSmean had positive relationships with AGB while Acv, Aratio, and Smean had negative relationships with AGB. For the H-related variables (Hmax, Hmean, etc.), the correlation coefficients of SER were higher or slightly lower (Δr < 0.01) than SLR. For the other TIN-based structural features, the correlation coefficients of SLR were higher than SER except for Acv, which had the lowest r of −0.190 (SLR) and −0.279 (SER) and a Δr of 0.089. Therefore, the models were built separately.

Since Hmean, Hp50, Hp75, Hp90, and Hp99 are all variables representing plant height, only the variable with the highest correlation coefficient (Hmean) was used as a plant height feature for AGB estimation and further analysis. Variables with high correlation coefficients (Hmax, Hmean, Hstd, Hcv, Hen, Amean, APmean, Aratio, Smean, Scv, and CSmean) (Figure 6) were selected to build univariate structural AGB estimation models. Hcv had the best estimation accuracy among the structural features (R² = 0.77, RMSE = 224.69 g/m², MAE = 175.31 g/m²). Hstd and Hen also had better estimation accuracy than the other structural features, but their performance was lower than Hcv (R² = 0.56, RMSE = 311.82 g/m², MAE = 250.10 g/m², R² = 0.55, RMSE = 315.47 g/m², MAE = 250.64 g/m²). The statistical parameters of the regression models are shown in

Table 7. Statistical parameters of regression models of TIN-based structural features and AGB (SLR) or ln (AGB) (SER). The best model is bolded.

Models	Variables	R-Squared (R2)	RMSE (g/m2)	MAE (g/m2)
SER	Hmax	0.37	376.08	299.60
	Hmean	0.53	321.91	251.89
	Hstd	0.56	311.82	250.10
	Hcv	0.77	224.69	175.31
	Hen	0.55	315.47	250.64
SLR	Amean	0.37	375.46	316.94
	APmean	0.51	330.77	269.51
	Aratio	0.50	337.84	288.29
	Smean	0.50	336.83	276.74
	Scv	0.45	352.26	291.23
	CSmean	0.50	333.55	275.51

.

3.3. Multivariate Biomass Estimation

We combined the structural and meteorological features with multispectral features. The combination of structural and meteorological features with multispectral features improved the estimation accuracy of AGB and the best model form varied between different combinations of variables.

3.3.1. Multispectral Features Combining Structural Features

Combining multispectral features with TIN-based structural features improved the estimation accuracy (

Table 8. Model type, combination form, and statistical parameters of the best-combined regression models of each selected multispectral features and structure features.

Models	Combination Form	R-squared (R2)	RMSE (g/m2)	MAE (g/m2)
SLR	ρ520/Hcv	0.73	247.40	190.49
SLR	ρ570/Hcv	0.76	233.47	173.11
SLR	ρ670×APmean	0.55	317.56	246.02
SLR	ρ680×APmean	0.55	317.09	242.37
SLR	ρ700/Hcv	0.78	221.67	169.61
SLR	ρ720/Hcv	0.79	214.28	160.42
SER	NDVI×Hcv	0.79	213.84	167.28
SER	GNDVI×Hcv	0.82	201.08	158.36
SER	NDRE×Hcv	0.83	192.78	156.48
SER	SR/APmean	0.81	205.61	163.73
SER	CIrededge×Hcv	0.83	196.28	160.38
SER	CIgreen×Hcv	0.78	223.13	182.60
SER	LCI×Hcv	0.82	200.01	159.71
SLR	MCARI×APmean	0.64	284.26	229.07

). However, the models and combination forms were different based on different features. SLR models had the best performance for specific bands and MCARI (R² = 0.55–0.79, RMSE = 317.56–214.28 g/m², MAE = 246.02–169.61 g/m²) while SER models suited other VIs (R² = 0.77–0.83, RMSE = 223.13–192.78 g/m², MAE = 182.60–156.48 g/m²). The best model was NDRE×Hcv with R² = 0.83, RMSE = 192.78 g/m², and MAE = 156.48 g/m². The best structural features selected for the combination were Hcv and APmean, but the combination form varied. Those features that were positively related to AGB were divided by Hcv or multiplied by APmean and those that were negatively related to AGB were multiplied by Hcv or divided by APmean.

3.3.2. Multispectral, Structural, and Meteorological Features

For the combination of three types of features, we selected some of the best multivariate combinations of multispectral-structural features and added meteorological features into the combinations. The results are shown in

Table 9. Statistical parameters of the best-combined regression models of multispectral, structural, and meteorological features.

Models	Variables	R-Squared (R2)	RMSE (g/m2)	MAE (g/m2)
SLR	ρ570/Hcv×GDD	0.82	200.27	156.37
SLR	ρ700/Hcv×GDD	0.82	202.44	163.16
SLR	ρ720/Hcv×GDD	0.85	179.80	141.83
SER	NDVI×Hcv/GDD	0.84	189.01	136.29
SER	GNDVI×Hcv/GDD	0.85	183.38	131.70
SER	NDRE×Hcv/GDD	0.86	178.37	127.34
SER	SR/APmean/GDD	0.81	204.94	153.01
SER	CIrededge×Hcv/GDD	0.85	181.45	130.00
SER	CIgreen×Hcv/GDD	0.83	198.38	145.12
SER	LCI×Hcv/GDD	0.86	178.49	127.53

. The addition of the meteorological feature improved the estimation accuracy with R² higher than 0.8, RMSE lower than 205 g/m², and MAE lower than 165 g/m² for every combination. The combination model forms were the same with the multispectral-structural models, which indicated that adding the meteorological feature had no influence on the combination model forms. Figure 7 shows the improvement of the estimation accuracy of the multivariate models.

3.4. Improvement of Estimation Accuracy Using Random Forest

The RF method was used to combine multispectral, structural, and meteorological features to estimate rice AGB. All the multispectral and structural features were used in the RF method. The statistical parameters are shown in

Table 10. Statistical parameters of the estimation of AGB using random forest (RF).

Model	Variables	R2	RMSE (g/m2)	MAE (g/m2)
RFa	Multispectral	0.88	158.40	116.69
RFb	3D	0.84	173.33	120.43
RFc	Multispectral, 3D	0.89	150.78	111.40
RFd	Multispectral, 3D, meteorological	0.92	126.28	92.18

. The estimation with only structural features yielded the lowest result (R² = 0.84, RMSE = 173.33 g/m², MAE = 120.43 g/m²). The model using only multispectral features was better (R² = 0.88, RMSE = 158.40 g/m², MAE = 116.69 g/m²). The combination of the multispectral, structural, and meteorological features was the best estimation model (R² = 0.92, RMSE = 126.28 g/m², MAE = 92.18 g/m²).

4. Discussion

The results of this study showed that both multispectral features and TIN-based structural features are capable of estimating rice AGB. The combination of these two types of features improves the estimation accuracy and the addition of meteorological information (GDD) further improves the performance.

When a single spectral variable was used for AGB estimation, we found that NDRE, CIrededge, and ρ700 had better performance among all the remotely sensed features. NDRE and CIrededge are red-edge-based vegetation indices and ρ700 is included in the red-edge spectrum. Previous study suggested that red-edge-based indices have a relationship with biomass. Mutanga and Kanke [64,65] reported that vegetation indices based on red-edge bands have strong relationships with AGB compared with red-based indices. Cheng [66] suggested that CIrededge has a significant relationship with leaf biomass and weaker relationship with total biomass. In our study, we found that not only red-edge-based vegetation indices, but also red-edge band reflectance have good relationships with biomass; this phenomenon might be caused by the fact that the red-edge spectrum are sensitive to chlorophyll absorption [60]. Compared with NDRE, NDVI had a weaker relationship with biomass (R² = 0.41) (Table 6). We believe this difference was caused by the saturation in the red spectral band after canopy closure [64], which, in this case, appeared since the first flight.

Among the TIN-based structural features, Hcv was the best predictor for AGB estimation. Previous study [17,22] suggested that plant height is a good predictor for biomass. However, in our study, we found that plant height was moderately correlated with biomass (R² = 0.53). This moderate correlation may be caused by the variation of the stem diameter, canopy thickness, leaf width, and leaf thickness among different hybrid species, which influence the biomass of the stem and leaves and cannot be reflected by crop height. Hcv is the coefficient of variance of the plant height, which reflects the roughness of the crop canopy within a plot. Higher Hcv values indicate a more uneven crop canopy surface [67]. The results in Table 7 are in line with the findings in [23,68,69], which suggested that Hcv is an important parameter for estimating biomass.

It is clear that height-related structural features are often exponentially related to biomass as reported in [17,69] and confirmed in this study. Interestingly, however, our results show that area-related and angle-related structural features are linearly related to biomass. To the best of the authors’ knowledge, such a relationship has not been reported before, and further investigation is required to confirm this relationship.

Aasen [70] suggested that the combination of multiple features may improve the estimation accuracy of crop parameters, which was confirmed in this study. Combining multispectral features with TIN-based structural features improved the estimation accuracy of biomass. This is probably explained by the fact that TIN-based structural features are structural features and they contain information about the crop canopy distribution both horizontally and vertically. The combination takes advantage of both features, i.e., the biochemical information and biophysical information. In this study, the multivariate models with Hcv and APmean produced the best accuracies while models with Hmean had a moderate performance. The Hmean did not appear to be the most influential factor for the model, which implies that plant height might not be the best feature for multivariate biomass estimation. Moreover, TIN was generated directly from dense point clouds created by the SfM algorithm while CSMs required the interpolation procedure of dense point clouds, which needs extra time and may cause a loss of detailed information of the canopy structure. Therefore, we believe the TIN-based structural features, like Hcv and APmean, are more appropriate for biomass estimation. The combination changed the estimation forms of VIs while the estimation forms of the band reflectance were unchanged, indicating that VIs are more easily influenced than band reflectance when combined with structural features. It should be noted that the combination with structural features significantly reduced the difference in the estimation accuracy between NDRE and NDVI (Table 8), indicating that the combination with structural features can relieve the saturation of NDVI after the canopy closure. A possible explanation is that NDVI is a vegetation index based on the reflectance of the red band (670 nm) and near-infrared (NIR) band (800 nm). When canopy closure appears (canopy cover reaches 100%), the absorption of red light reaches its maximum value [71]. At the same time, NIR reflectance still increases with biomass accumulation because the addition of leaves results in multiple scattering [72]. The imbalance between the reflectance of the red band and NIR band results in an almost unchanged NDVI ratio while biomass keeps increasing. Structural features used in this study, like Hcv and APmean, describe the vertical and horizontal distribution of rice leaves. By combining structural features with NDVI, the lack of distribution information of leaves was compensated, thus improving the estimation accuracy of NDVI. However, this hypothesis needs to be tested in a future study.

Traditionally, temporal variables, like GDD, were used to predict biomass [73]. In this study, we used GDD as one of the model features to support biomass estimation. The results show an improved estimation accuracy, which is in line with the results in [73]. This suggests that it is worthwhile to consider more features in the estimation model that are not directly related to the properties of biomass.

We applied the random forest algorithm to different combinations of features to verify our results and the algorithm performed better than the univariate and multivariate estimation models. The random forest algorithm makes full use of all input data and tolerates outliers and noise, hence we chose it as the verification algorithm. The best model was RFd, which contains all spectral, structural and meteorological features, while the model containing only structural features (RFb) performed the worst. The results confirmed that combining different types of features improves estimation accuracy.

Previous studies used RGB indices/multispectral indices in combination with RGB SfM height/multispectral SfM height to estimate biomass. Either the spectral resolution or the image resolution was low, reducing estimation accuracy. To solve the problem of resolution, Tilly [2] used TLS and a field spectrometer to acquire plant height and vegetation indices. However, both instruments are ground-based sensors and cannot be applied over a large area. Furthermore, they are expensive and require extra training to operate. We used a combination of the MCA multispectral camera (12 bands) and DJI FC6310 RGB camera (20 M pixels) to acquire high spectral resolution and high spatial resolution images. Furthermore, we used TIN instead of raster DEM to ensure that the high spatial resolution was preserved. Additional attention needs to be paid to avoid hotspot and sun glint effects (Figure 2) since both introduce unexpected error [48]. Our future work will evaluate the transferability of this multivariate biomass estimation model to varying crops and climate conditions.

5. Conclusions

In this study, we developed an approach to estimate rice AGB using multispectral, structural, and meteorological features. We proposed TIN-based structural features instead of the commonly used plant height for more detailed information about the rice canopy structure. Multispectral, structural features and their combinations were used to estimate rice AGB. TIN-based structural features performed better than plant height for univariate and multivariate AGB estimation, indicating that more detailed structural information should be considered for AGB estimation. The saturation problem of NDVI for AGB estimation was relieved after combining with structural features. Both the simple regression result (SLR and SER) and machine learning result (RF) demonstrated that the combination of multispectral and structural features had better AGB estimation accuracy than the use of multispectral or structural features alone, and the addition of the meteorological feature to the combination further improved the estimation accuracy.