Spectrum- and RGB-D-Based Image Fusion for the Prediction of Nitrogen Accumulation in Wheat

Abstract

The accurate estimation of nitrogen accumulation is of great significance to nitrogen fertilizer management in wheat production. To overcome the shortcomings of spectral technology, which ignores the anisotropy of canopy structure when predicting the nitrogen accumulation in wheat, resulting in low accuracy and unstable prediction results, we propose a method for predicting wheat nitrogen accumulation based on the fusion of spectral and canopy structure features. After depth images are repaired using a hole-filling algorithm, RGB images and depth images are fused through IHS transformation, and textural features of the fused images are then extracted in order to express the three-dimensional structural information of the canopy. The fused images contain depth information of the canopy, which breaks through the limitation of extracting canopy structure features from a two-dimensional image. By comparing the experimental results of multiple regression analyses and BP neural networks, we found that the characteristics of the canopy structure effectively compensated for the model prediction of nitrogen accumulation based only on spectral characteristics. Our prediction model displayed better accuracy and stability, with prediction accuracy values (R²) based on BP neural network for the leaf layer nitrogen accumulation (LNA) and shoot nitrogen accumulation (SNA) during a full growth period of 0.74 and 0.73, respectively, and corresponding relative root mean square errors (RRMSEs) of 40.13% and 35.73%.

1. Introduction

Nitrogen is a key element required for the growth and development of wheat, and thus exerts an important influence on its yield and quality. At present, during wheat production, in order to pursue high crop yields, problems exist such as excessive nitrogen fertilizer application and unscientific fertilization methods, which not only cause environmental pollution, but also reduce profits [1,2,3]. Therefore, the precise management of wheat fertilization and the realization of scientific nitrogen application not only have the production significance of saving nitrogen and increasing efficiency, but also the ecological significance of reducing farmland pollution.

A quick and effective grasp of the status of nitrogen nutrition indicators during the key growth period of wheat is a prerequisite for achieving accurate fertilizer management. With the continuous development of spectral analysis and image sensing, the timeliness of the non-destructive testing of nitrogen nutrition indices is steadily improving, and the testing method is becoming increasingly convenient [4,5,6,7]. The premise of spectral analysis technology is primarily to simplify and idealize the crop canopy as an isotropic and uniform Lambertian surface. In previous research, there have been many investigations on the real-time estimation of nitrogen nutrition indices using multi-spectral and hyper-spectral remote sensing equipment [5,8,9,10,11], reaching conclusions including that the sensitive bands of crop nitrogen nutrition status are mainly concentrated in the visible and near-infrared bands, with the “red edge” area being the most prominent [12,13]. Existing research on the interaction mechanism of the vegetation canopy and radiation, however, indicates that the assumption that the vegetation canopy should be regarded as a Lambertian surface is quite different from the actual situation. When the influence of the anisotropic characteristics of canopy structure on the spectrum is ignored, the detection results are not stable enough, and the accuracy is relatively poor [14,15]. Therefore, it is necessary to combine spectral characteristics and canopy structure characteristics in order to detect the nitrogen nutrition status of crops, because the combination of the two characteristics can have a more positive impact on the predicted results. At present, the extraction of canopy structure features mostly relies on digital image processing technology, and the differences in color and morphological characteristics of the crop canopy are expressed by analyzing image features [16,17,18]. Since only visible light band information is used, the R, G, and B values of the image are easily affected by lighting conditions at the time the image is captured, and, since it is difficult for a two-dimensional digital image to reflect the three-dimensional morphological information of the crop canopy, the accuracy of constructing a nitrogen nutrition index estimation model based on two-dimensional digital images is not ideal [19].

In recent years, depth information has gradually been applied to the three-dimensional reconstruction of plants and the monitoring of crop growth, which has the advantage of quickly reconstructing the three-dimensional structure of crops [20]. Compared with RGB and grayscale images, depth images contain more accurate geometric information of the features. Consumer-grade RGB-D cameras provide an inexpensive solution for the rapid acquisition of 3D crop information. Existing studies have shown that RGB-D cameras provide a convenient way to quickly obtain crop canopy structure information in a field environment [21,22]. The applications of RGB-D cameras include extracting crop canopy structure parameters through three-bit reconstruction [23,24], and using height differences to remove weed information [25].

Because spectral characteristics are affected by canopy structure during the prediction of nitrogen accumulation in wheat, the model prediction accuracy is not high. When extracting canopy structure parameters, RGB images cannot reflect the three-dimensional structural characteristics of the canopy, while depth images can better reflect the three-dimensional structural information of crops [23]. In order to consider both RGB images and depth images in the construction of canopy structure features, we propose a pixel-level fusion algorithm of RGB images and depth images based on the correlation between the corresponding pixels of the two image types. By constructing characteristic parameters that combine the color and structural information of the wheat canopy, we compensate for the previous model limitations when predicting the leaf layer nitrogen accumulation (LNA) and shoot nitrogen accumulation (SNA) of wheat with a single spectral feature, and establish a predictive model based on multi-dimensional features.

2. Materials and Methods

2.1. Wheat Experiment Design

From December 2017 to April 2018, a wheat experiment (Figure 1) was conducted at the Demonstration Base of the National Engineering and Technology Center for Information Agriculture in Rugao City (32°16′N, 120°46′E), Nantong City, Jiangsu Province, China. The wheat variety selected for the experiment was Shengxuan 6, and the experiment was divided into three nitrogen treatments, namely N0 (0 kg N/ha), N1 (180 kg N/ha), and N2 (360 kg N/ha). The planting densities had row spacings of 20 cm, 35 cm, and 50 cm. Three replicates were set for each group of experiments, for a total of 27 plots, with the area of each plot being 5 × 6 m².

2.2. Obtaining Wheat Agronomic Parameters

Based on the key fertilization periods of wheat, samples were taken at the jointing, booting, and heading stages. Due to the relatively long duration of the jointing stage, two samples were taken during this stage. From each plot, we randomly selected 15 wheat plants for sample analysis, and then separated them indoors according to their organs (stems, leaves, and ears). After the samples were fixed at 105 °C for 30 min, they were baked to a constant weight at 80 °C and then weighed. The individual organ weights, i.e., the leaf dry weight ($LDW$), stem dry weight ($S_{tem}DW$), and spike dry weight ($S_{pike}DW$), were measured. After each sample was crushed, the Kjeldahl method was used to determine the leaf nitrogen content ($LNC$), stem nitrogen content ($S_{tem}NC$), and spike nitrogen content ($S_{pike}NC$). We then calculated the leaf layer nitrogen accumulation ($LNA$) and shoot nitrogen accumulation ($SNA$) according to the following equations.

$$LNA(g/m^2)=LNC(\%)\times LDW(g/m^2)$$

(1)

$$SNA(g/m^2)=LNA(g/m^2)+S_{tem}NC(\%)\times S_{tem}DW(g/m^2)+S_{pike}NC(\%)\times S_{pike}DW(g/m^2)$$

(2)

2.3. Wheat Spectrum and Image Data Acquisition

The ASD FieldSpec HandHeld 2 was used to obtain the spectral data of the crop. The measurement position was 1 m from the top of the canopy. Each plot was measured three times at different locations, with the average value representing the spectral information of the plot. We used the Intel^® RealSense™ D415 Depth Camera to acquire wheat field images. As shown in Figure 2b, the camera has two infrared stereo cameras, one infrared emitter, and one color sensor. The infrared stereo camera was used to generate depth images, and the color sensor was used to generate RGB images, both with a resolution of 1280 × 720 pixels. The camera was positioned 70 cm from the crop canopy, and we selected three different locations in each plot for image acquisition. The actual image acquisition device in the field is shown in Figure 2.

2.4. RGB-D Image Preprocessing

2.4.1. Depth Image Restoration Method

Due to lighting conditions, infrared reflection characteristics of the surface material of the measured object, and area occlusion, depth images have the problem of missing data (appearing as holes), resulting in decreased accuracy of the target object information, which affects the subsequent feature extraction. As shown in Figure 3, we divide the occurrence of holes into two categories for discussion: local holes and large-area holes. Local holes appear within the leaf, and are mainly affected by the surrounding environmental conditions, such as light and vibration. The area of this type of hole is usually small, and the depth information is similar to the surrounding depth information. Large-area holes appear at the junctions of different areas. The reason for their occurrence is related to the application range of the sensor. They mainly appear at the transition of depth information on the edge of the leaf. Based on the analysis of different types of holes, we designed a hole-filling algorithm (HF), consisting of the following steps. First, perform column traversal on the depth information of a crop row image, find the positions of all the hole points (line), and calculate the difference H of the non-zero depth information before and after the point (line). By comparing H with an artificially set threshold T, the types of holes are distinguished. When H < T, it is a local hole, and the mean value of the non-zero depth information before and after the hole point (line) is interpolated. When H ≥ T, it is a large-area hole, and the non-zero depth information before and after the hole point (line) is separately fitted and interpolated. Since there were horizontal crop rows in the center of the images we examined, a curve was more suitable for fitting. We used a quadratic polynomial for fitting. HF is expressed as Equation (3), where $Hb$ and $Hf$ represent the non-zero depth information before and after the hole point (line), respectively; $x$ represents the column coordinates corresponding to the hole point; and $d(x)$ is the depth information to be filled.

$$d(x)=\{\begin{array}{c}\frac{Hb+Hf}{2} H<T\\ a{x}^2+bx+c H\ge T\end{array}.$$

(3)

2.4.2. Wheat Canopy Segmentation Method Based on RGB-D Images

There are non-crop canopy features such as soil, straw, and weeds in canopy images of wheat during its early and mid-term growth. In order to eliminate the interference of non-canopy information in the images during feature extraction, the canopy portion of the images needs to be segmented and extracted. Common vegetation indices such as ExG and ExGR are greatly affected by environmental light when performing canopy segmentation and cannot remove green information such as weeds. When the light is strong, the color of the reflective leaf is closer to the color of the soil background, and is thus easily removed by mistake, resulting in poor segmentation. In order to avoid the influence of ambient light on the segmentation results, we propose segmenting the soil background in hue, saturation, lightness (HSI), and space. The value range of the H channel in the HSI color model is $[0,2\pi ]$, and the H values corresponding to yellow, green, and cyan are $\frac{\pi }{3}$, $\frac{2\pi }{3}$, and $\pi$, respectively. Due to different nitrogen treatments of wheat, the colors of the leaves varied, so the H channel value in the interval $(0,\frac{4\pi }{3})$ was retained as crop row information. In addition, considering the influence of field weed information on the extraction of crop row characteristics—since the wheat closes the rows after the jointing stage—weeds no longer affect the subsequent extraction of wheat canopy characteristics, and were only removed during the jointing stage via height information obtained from the depth images.

2.4.3. RGB-D Image Fusion Method for the Wheat Canopy

When most studies use RGB images to obtain canopy structure information, they extract coverage and other information through color differences, which cannot express the characteristics of wheat canopy structure in a three-dimensional space [26]. When using depth information for crop growth parameter inversion, researchers generally extract canopy height information to predict above-ground biomass and other agronomic parameters, and no further methods for extracting canopy structure features based on three-dimensional space have been proposed [27]. In order to simultaneously fuse RGB information and depth information to express the characteristics of canopy structure, we propose a pixel-level RGB-D image fusion method, which uses the connection between the corresponding pixels of RGB images and depth images to complement the two kinds of information in the modeling process.

Since the imaging origins of RGB images and depth images are inconsistent, depth information and RGB image information cannot correspond one-to-one. For the subsequent feature extraction, the first step is to align two images. The basic idea of image alignment is the conversion of different coordinate systems. By converting the coordinate system of a depth image to the world coordinate system and then to the coordinate system of an RGB image, the alignment of the two images is realized.

The basic idea of the RGB-D image fusion method is to use the physical relationship between intensity and phase. In the natural environment, color information in an RGB image indicates the reflection intensity of the object, and the intensity information is affected by different lighting conditions. The principle of acquiring depth information is that the sensor collects different phase information for different depth areas of the subject, and then converts these structural variations into depth information through the arithmetic unit. Therefore, the inverse square law can be used to describe the relationship between intensity and phase [28]:

$$d\approx I_R/{\varphi }^2,$$

(4)

where $I_R$ represents the intensity information of each point in the RGB image, corresponding to the grayscale value after the RGB image is converted to a grayscale image, and $d$ represents the depth information corresponding to $I_R$. From the inverse square law, the phase information $\varphi$ of each pixel in the image can be calculated. When performing image fusion, we use the idea of IHS transform to replace the I component in the HSI model with the phase information ($\varphi$) and depth information ($d$) for representation, thus obtaining a fused image based on HSP (where P is phase information) and HSD (where D is depth information) models. Since the value range of the I channel is $[0,1]$, in order to simplify the calculation, the depth information is normalized. The normalization formula is as follows:

$$\overline{d}=(d-d_{\min })/(d_{\max }-d_{\min }),$$

(5)

where $d$ is the actual depth information, $d_{\max }$ and $d_{\min }$ are the maximum and minimum depth information, respectively, and $\overline{d}$ represents the normalized results.

2.5. Wheat Canopy Feature Extraction

2.5.1. Spectral Feature Extraction

Previous studies have shown that the spectral reflectance of crops at 730 nm and 815 nm performs well in assessing leaf area index, pigment status, and nitrogen nutrition status [29,30]. Therefore, we selected the above two bands to construct the normalized difference red edge (NDRE) index, and the ratio vegetation index (RVI). The calculation formulas of vegetation indices are as follows:

$$NDRE=(R815-R730)/(R815+R730),$$

(6)

$$RVI=R815/R730,$$

(7)

where $R730$ and $R815$ represent the reflectance of wheat in the 730 and 815 nm bands, respectively.

2.5.2. Feature Extraction of Fused Images

In order to represent the information of the wheat canopy structure more accurately, a single row of crops was intercepted for feature extraction, and the crop row image block size intercepted in all periods was 150 × 500 pixels. After image fusion, the color co-occurrence matrix was used to extract texture features. The color co-occurrence matrix usually contains a large amount of information. In order to reduce the dimensionality of the feature space as much as possible to decrease the amount of calculation, while still retaining the description of the image texture information, we selected five feature statistics—entropy (ENT), inverse difference moment (IDM), contrast (CON), correlation (COR), and angular second moment (ASM)—to represent texture [31].

The specific calculation process of texture features is as follows: First, the RGB image is converted to color space, and its corresponding value in the HSD and HSP space models is calculated. In this study, the H, S, D, and P components were non-uniformly quantized into eight parts in order to reduce information redundancy and make the calculation of the co-occurrence matrix (CCM) easier. When calculating the CCM, set the four components as $C_1$, $C_2$, $C_3$, and $C_4$, set $m=C_k$, $n=C_{k^{\prime }}$ as the two components in the four-component combination space $(k,k^{\prime }\in \{1,2,3\})$, and use the color co-occurrence matrix $CC{M}_{mn}$ to represent the different components, $C_k$ and $C_{k^{\prime }}$, of the pixels in the image; specifically, the measurement of the interaction between $m$ and $n$ space. From this, $CC{M}_{HH}$, $CC{M}_{SS}$, $CC{M}_{PP}$, $CC{M}_{DD}$, $CC{M}_{HS}$, $CC{M}_{HP}$, $CC{M}_{SP}$, $CC{M}_{HD}$, and $CC{M}_{SD}$ can be calculated. Since the symbiosis matrices $CC{M}_{HH}$, $CC{M}_{SS}$, $CC{M}_{PP}$, and $CC{M}_{DD}$ only have elements on the diagonal, these elements have little effect on distinguishing various textures, and thus the remaining four matrices are used for calculation. For any pixel in the image, assuming that the value of its $k$th component is $i$, i.e., $m=i$, and the value of the $k^{\prime }$th component is $j$, i.e., $n=j$, then the element $CC{M}_{m.n}(i,j)$ in the matrix is used to represent the number of occurrences of such pixels in the image. The calculation formula of the color co-occurrence matrix is as follows:

$$CC{M}_{m,n}(i,j)={\displaystyle \sum _{\Delta x}{\displaystyle \sum _{\Delta y}\{\begin{array}{ll}1& m(x+\Delta x,y+\Delta y)=i \& n(x+\Delta x,y+\Delta y)=j)\\ 0& \mathrm{others}\end{array}}},$$

(8)

where the $i$th row and $j$th column elements in the color co-occurrence matrix are represented by $p(i,j)$.

2.6. Evaluation Index of the Fused Image

Two objective indicators—entropy and cross-entropy—were used to evaluate the quality of the fused images. Before image evaluation, each image in the HSI space is converted to an RGB image, and the grayscale operation is performed. The entropy of an image (E) is an important index for measuring the richness of image information. The larger the entropy, the more information the image has. Cross-entropy (C) can be used to determine the information difference between the grayscale distribution of two images and is a key index for evaluating the difference between two images. The smaller the cross-entropy, the smaller the difference between the fused image and the source image, and the more information the fusion method extracts from the source image. The calculation formulas of entropy and cross-entropy, respectively, are:

$$E={\displaystyle \sum _{i=0}^{L-1}{p}_i{\log }_2}{p}_i,$$

(9)

$$C={\displaystyle \sum _{i=0}^{L-1}{p}_i{\log }_2}\frac{p_i}{q_i},$$

(10)

where $L$ represents the grayscale level of the image, which is 256, and $p_i$ and $q_i$ represent the ratio of pixels with grayscale value to the total pixels in the source image and fused image, respectively.

2.7. Model Establishment and Evaluation Methods

2.7.1. Evaluation Method of the Hole Repair Algorithm

In order to evaluate the effectiveness and accuracy of the hole repair algorithm, we artificially created holes in existing depth images, including local holes and large-area holes at the edges, and evaluated the algorithm by calculating the true value and the RMSE of the filling value calculated by the algorithm. The calculation formula of the RMSE is as follows, in which $M_i$ represents the repair value of each hole point, $T_i$ represents the true value corresponding to the hole point, and n represents the number of hole points:

$$RMS{E}_{hole}=\sqrt{\frac{1}{n}({\displaystyle \sum _{i=1}^n{(M_i-T_i)}^2})}.$$

(11)

2.7.2. Establishment and Evaluation Method of the Wheat Nitrogen Accumulation Model

Due to high feature dimensionality, we selected two methods—multiple regression analysis and the BP neural network—to separately construct the model. In order to reduce the phenomenon of overfitting and improve the generalization ability of the model, the prediction result of the BP neural network model was the mean value of five trials, in which 50% of the samples were randomly selected each time as a training set (14 for a single growth period, and 54 for a full growth period), 20% as the validation set (five for a single growth period and 22 for a full growth period), and 30% as a test set (eight for a single growth period and 32 for a full growth period). The multiple linear regression model is as follows:

$$y={\displaystyle \sum _{i=1}^{num}{a}_i{x}_i+b},$$

(12)

where $num$ is the number of different features, $x_i$ represents the values of different features, $y$ represents the values of LNA and SNA, and $a_i$ and $b$ represent the parameters of the model. The structure of the BP neural network is presented in Figure 4. There are 10 neurons in the hidden layer, and the training function is Bayesian regularization. The model was evaluated by calculating R² and RRMSE. The calculation formulas of the RMSE and RRMSE are as follows:

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{(y_{pi}-y_r{}_i)}^2}},$$

(13)

$$RRMSE(\%)=\frac{100}{\overline{y_r}}\sqrt{\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{(y_{pi}-y_r{}_i)}^2}},$$

(14)

where $y_{pi}$, $y_{ri}$, and $\overline{y_r}$ represent the predicted value, measured value, and mean value of the wheat LNA and SNA, respectively, and $N$ is the sample size.

3. Results and Analysis

3.1. Hole-Repair Results

Since the T value setting in the HF algorithm will directly affect the classification results of the holes, it will thus affect the accuracy of the subsequent hole repair. In order to determine the best difference threshold T for the hole-repair algorithm, we employed different T values for experimentation. The results are listed in

Table 1. Hole-repair algorithm results.

T (mm)	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{h}\mathit{o}\mathit{l}\mathit{e}}$
20	21.09
30	20.32
40	11.44
50	14.59

. We chose a T value of 40 to repair the entire depth image. Figure 5b shows the repair result of the hole information. Compared with the original image shown in Figure 5a, it is obvious that the repaired image has no black holes, and the repair quality is good. By comparing Figure 5c,d, i.e., before and after the repair, it can clearly be seen that the edge information (the transition of depth information) and the information of the local hole have been repaired, and the overall curve is smooth.

3.2. Results of Wheat Canopy Segmentation

The basic idea of the traditional color vegetation index for canopy segmentation is to use the difference between color information of the crop row and soil background, although the segmentation effect will be affected by factors such as light conditions and weeds. We converted the RGB color model to an HSI color model and performed canopy segmentation using the definition of hue to reduce the impact of light changes on crop row extraction. In addition, the depth information was used to remove weeds, thus reducing the interference of the green background of the non-wheat canopy on the extraction of crop rows. By analyzing the heights of wheat and weeds at the jointing stage, the maximum height of wheat and the average height of weeds under different nitrogen treatments were counted. The results are shown in Figure 6. Taking into account the height difference of the leaves at different positions of the wheat, only the pixels below 10% of the maximum height were regarded as weed information and removed. Figure 7 shows the crop row information extracted by the ExG and ExGR vegetation indices. It can be seen that the reflection of the leaves on the right side of the image has a poor segmentation effect. Figure 7b shows the canopy segmentation resulting from simultaneously removing the soil background and weed information, which is better than the segmentation effect of the vegetation index that only depends on color features. As seen in the figure, the outline of the crop rows is clear, and the reflective part of the leaves does not affect the segmentation.

3.3. Fused Image Quality Analysis

Figure 8 shows the image fusion results under low light, strong light, and normal light.

Table 2. Entropy and cross-entropy of source image and fused images.

	Quality	Weak Light	Strong Light	Normal Light
Original image	Entropy	6.6657	6.8190	6.7336
Fusion image1	Entropy	6.9854	7.1891	7.1824
	Cross-entropy	0.1318	0.0896	0.0854
Fusion image2	Entropy	6.7133	6.8487	6.8019
	Cross-entropy	0.1518	0.1147	0.1827

shows the quality evaluation results of fused images under different lighting conditions, in which fused image 1 represents the fused image based on phase information, and fused image 2 represents the fused image based on depth information. The calculation results reveal that the entropy of the fused image under different lighting conditions is greater than that of the original image, while the cross-entropy is smaller, indicating that the quality of the fused image has been improved and the information of the original image has been better retained.

3.4. Forecast Results of Wheat Nitrogen Accumulation

3.4.1. Feature Screening of Fused Images of the Wheat Canopy

A total of 25 features were extracted from the color co-occurrence matrix of two fused images. The feature data was linearly fitted with the wheat LNA and SNA, and feature screening was performed by comparing the R² magnitude and the results are shown in Figure 9. In order to select textural features that have strong applicability in different growth periods, four features of angular second moment and entropy derived from the CCM_HP and CCM_HD were selected in the early jointing, late jointing, and booting stages. During the heading stage, however, the canopy structure changes significantly due to the wheat heading. Four features were selected, including the contrast derived from the CCM_HS, the second moment and entropy of the CCM_HP symbiosis matrix angle, and the entropy calculated in the CCM_HD.

3.4.2. Forecast Results of Nitrogen Accumulation in Wheat during Different Growth Periods

Table 3. Comparison of LNA and SNA estimates derived from different features in the early jointing stage, in which F1 represents the feature set extracted based on the fused image, F2 represents the set based on the fused image features and the RVI, and F3 represents the set based on the fused image features and the NDRE.

LNA					SNA
Feature	Model	R2	RMSE	RRMSE	Feature	Model	R2	RMSE	RRMSE
NDRE	LR	0.37	1.03	45.47%	NDRE	LR	0.31	2.92	47.91%
RVI	LR	0.30	1.09	48.12%	RVI	LR	0.24	3.04	49.88%
F1	MR	0.63	0.85	37.53%	F1	MR	0.63	2.28	37.41%
F1	BP	0.69	0.93	41.06%	F1	BP	0.67	2.33	38.23%
F2	MR	0.65	0.84	37.08%	F2	MR	0.63	2.31	37.90%
F2	BP	0.72	0.94	41.50%	F2	BP	0.70	2.24	36.75%
F3	MR	0.66	0.83	36.64%	F3	MR	0.64	2.29	37.57%
F3	BP	0.70	0.89	39.29%	F3	BP	0.69	2.36	38.72%

shows the prediction ability of wheat LNA and SNA with different features during the early stage of jointing. By comparing the prediction results of different features, it was discerned that the prediction accuracy based on spectral features is significantly lower than that based on fusion features. This is because in the early stage of jointing, wheat plants are small, and soil background has a huge impact on the collection of canopy reflectivity, resulting in low prediction accuracy. Since the features extracted by image fusion are not affected by soil background, the prediction effect of fused image features is better than that of the spectral features.

As the wheat grows and forms joints, it gradually closes the rows, and the collection of canopy reflectance is no longer disturbed by soil background.

Table 4. Comparison of LNA and SNA estimates derived from different features in the late jointing stage.

LNA					SNA
Feature	Model	R2	RMSE	RRMSE	Feature	Model	R2	RMSE	RRMSE
NDRE	LR	0.68	0.97	30.30%	NDRE	LR	0.72	1.24	24.20%
RVI	LR	0.70	0.93	29.05%	RVI	LR	0.72	1.23	24.01%
F1	MR	0.58	1.17	36.54%	F1	MR	0.62	1.63	31.81%
F1	BP	0.69	1.07	33.42%	F1	BP	0.71	1.57	30.64%
F2	MR	0.76	0.91	28.42%	F2	MR	0.80	1.20	23.42%
F2	BP	0.86	1.00	31.23%	F2	BP	0.88	1.09	21.27%
F3	MR	0.75	0.93	29.05%	F3	MR	0.81	1.19	23.23%
F3	BP	0.81	0.92	28.73%	F3	BP	0.86	1.06	20.69%

,

Table 5. Comparison of LNA and SNA estimates derived from different features in the booting stage.

LNA					SNA
Feature	Model	R2	RMSE	RRMSE	Feature	Model	R2	RMSE	RRMSE
NDRE	LR	0.69	2.00	35.19%	NDRE	LR	0.71	2.90	31.57%
RVI	LR	0.70	1.99	35.01%	RVI	LR	0.74	2.90	31.57%
F1	MR	0.71	2.08	36.60%	F1	MR	0.71	3.09	33.64%
F1	BP	0.76	1.70	29.91%	F1	BP	0.82	2.16	23.52%
F2	MR	0.76	1.93	33.96%	F2	MR	0.78	2.78	30.27%
F2	BP	0.81	1.35	23.75%	F2	BP	0.89	1.95	21.23%
F3	MR	0.76	1.92	33.78%	F3	MR	0.78	2.76	30.05%
F3	BP	0.80	1.99	35.01%	F3	BP	0.83	2.38	25.91%

and

Table 6. Comparison of LNA and SNA estimates derived from different features in the heading stage.

LNA					SNA
Feature	Model	R2	RMSE	RRMSE	Feature	Model	R2	RMSE	RRMSE
NDRE	LR	0.71	1.68	38.19%	NDRE	LR	0.77	2.87	32.08%
RVI	LR	0.74	1.61	36.60%	RVI	LR	0.79	2.71	30.29%
F1	MR	0.63	2.04	46.37%	F1	MR	0.79	2.62	29.28%
F1	BP	0.74	1.60	36.37%	F1	BP	0.78	2.49	27.83%
F2	MR	0.76	1.93	43.87%	F2	MR	0.84	2.48	27.72%
F2	BP	0.87	1.17	26.59%	F2	BP	0.88	2.36	26.38%
F3	MR	0.77	1.65	37.51%	F3	MR	0.82	2.79	31.18%
F3	BP	0.84	1.40	31.82%	F3	BP	0.87	2.18	24.36%

present the prediction results of different characteristics in the late jointing, booting, and heading stages, respectively. The prediction accuracy of wheat LNA and SNA based on spectral characteristics during these stages is higher and tends to be stable. Meanwhile, the prediction accuracy based on fused image features is close to the prediction accuracy during the early stage of jointing, indicating that the canopy structure features reflected by fused images are not affected by the growth period and are relatively consistent in predicting wheat LNA and SNA at different growth stages.

By comparing the prediction results of wheat LNA and SNA based on fused image features at different growth stages listed in Table 3, Table 4, Table 5 and Table 6 it was found that the prediction of SNA based on fused image features is better. This is because the depth information in the fused image not only contains the structure of the wheat leaf layer, but also reflects the height of the wheat stalk, so it can more accurately express the nitrogen accumulation of the aboveground wheat population.

Table 3, Table 4, Table 5, Table 6 and

Table 7. Comparison of LNA and SNA estimates derived from different features in all stages.

LNA					SNA
Feature	Model	R2	RMSE	RRMSE	Feature	Model	R2	RMSE	RRMSE
NDRE	LR	0.63	1.75	45.02%	NDRE	LR	0.56	3.16	43.42%
RVI	LR	0.66	1.68	43.22%	RVI	LR	0.59	3.05	41.91%
F1	MR	0.46	2.12	54.53%	F1	MR	0.54	3.29	45.21%
F1	BP	0.61	2.15	55.31%	F1	BP	0.66	2.83	38.89%
F2	MR	0.72	1.53	39.36%	F2	MR	0.70	2.66	36.55%
F2	BP	0.74	1.56	40.13%	F2	BP	0.73	2.60	35.73%
F3	MR	0.70	1.58	40.64%	F3	MR	0.68	2.72	37.38%
F3	BP	0.77	1.29	33.18%	F3	BP	0.71	2.67	36.69%

show the prediction accuracy of different features in different models. In all growth periods, the prediction results of the BP neural network were found to be better than those of the multiple regression model. There is an incomplete linear correlation between the fused image features and LNA and SNA, and since the multiple regression model cannot accurately describe the nonlinear correlation, the result is a certain deviation. The prediction results of all growth periods reveal that the combination of RVI and fused image features performs best, with R² > 0.7 in the BP neural network. In addition, in the multiple regression model, the integration of the fused image features also significantly improved the prediction accuracy of the spectral features for LNA and SNA. This proves that the fused image features expressing canopy structure parameters can compensate for the LNA and SNA prediction models based on spectral features to a certain extent, and effectively improve the prediction accuracy of the models.

4. Discussion

4.1. Prediction Model Performance of Wheat Nitrogen Accumulation Based on Spectral Characteristics

In previous studies, multispectral remote sensing technology has been used to estimate the amount of nitrogen accumulation in wheat [32,33]. It was found that in the range f 550–750 nm for winter wheat, the spectral reflectance of the canopy changes with the change of leaf nitrogen content [34]. The spectral information in the red light portion of the spectrum is often used as an index for predicting nitrogen content [35], which proves the potential of spectral features to predict nitrogen accumulation in wheat. However, the spectral reflectance of the crop canopy exhibits obvious anisotropic characteristics due to the differences in canopy structure. The nitrogen accumulation prediction model based on spectral features ignores the influence of canopy structure on spectral information, resulting in low prediction accuracy and insufficient stability of prediction models [7,13,34]. In our study, depth images and RGB images were combined to extract wheat canopy texture features with depth information, and wheat canopy structure information was expressed from a three-dimensional perspective. The experimental results revealed that when predicting the nitrogen accumulation in wheat, the three-dimensional structural characteristics of the wheat canopy compensated for the shortcomings of the nitrogen accumulation prediction model, effectively improving the prediction accuracy of the model.

4.2. Expression Method of Wheat Canopy Structure Based on RGB-D Fused Images

Image fusion is generally divided into three levels: pixel-level, feature-level, and decision-level. At present, in most studies, feature-level fusion methods have often been used when predicting crop canopy structure parameters with RGB-D data [36,37], i.e., extracting features from RGB images and depth images separately, and predicting parameters by combining features. This method, however, often ignores the relationship between the corresponding pixels in the RGB images and depth images. The pixel-level fusion method can retain as much detailed information of the original image as possible. After fusing depth images with RGB images at the pixel level in our study, the extracted fused image features contained the depth information of crop rows, which could break through the limitations of two-dimensional images and express canopy structure information of wheat more accurately. The extracted fused image features cannot quantitatively describe canopy structure information, however. In subsequent research, we will continue to explore the extraction of more detailed canopy structure features from point cloud information, such as leaf inclination.

4.3. Generalization Ability of the Wheat Nitrogen Accumulation Prediction Model

The generalization ability of a model refers to its ability to adapt to fresh samples. The generalization ability directly determines the stability and versatility of the model. We utilized the independent experimental data of a wheat variety under three nitrogen treatments and two density treatments to construct and verify the model, and to improve the generalization ability of the model through cross-validation. When applying the model to other varieties, however, the model parameters still need to be corrected. This is because the plant height, leaf inclination angle, tiller number, and leaf development process of different varieties are quite different. In the follow-up study, different plant types of wheat varieties will be tested in order to further improve the generalization ability of the model among different varieties of wheat.

5. Conclusions

We proposed a pixel-level RGB-D image fusion method, which expresses three-dimensional structural information of the canopy by extracting textural features of fused images, thereby compensating the model for predicting nitrogen accumulation based only on spectral characteristics. We evaluated and analyzed the compensation effect by using two models—multiple regression analysis and the BP neural network. Our work is summarized as follows:

(1)

We proposed a hole-filling algorithm with depth information. When T = 40 mm, the repair accuracy is the best (RMSE = 11.44), effectively avoiding the loss of features due to information loss.

(2)

We used IHS transformation to fuse RGB images and depth images. The textural features of the fused images contain depth information, which can break through the limitation of extracting canopy structure features from a two-dimensional image and express the three-dimensional structure information of wheat canopy.

(3)

Our test results revealed that in different models, the combination of the characteristics of the fused image with the spectral characteristics can predict nitrogen accumulation in wheat more accurately. In the prediction of a full growth period, the best prediction accuracy values (R²) of the combined features for LNA and SNA were 0.74 and 0.73, respectively, with corresponding RRMSEs of 40.13% and 35.73%, indicating good predictive ability.