Estimating the Leaf Nitrogen Content with a New Feature Extracted from the Ultra-High Spectral and Spatial Resolution Images in Wheat

Abstract

Real-time and accurate monitoring of nitrogen content in crops is crucial for precision agriculture. Proximal sensing is the most common technique for monitoring crop traits, but it is often influenced by soil background and shadow effects. However, few studies have investigated the classification of different components of crop canopy, and the performance of spectral and textural indices from different components on estimating leaf nitrogen content (LNC) of wheat remains unexplored. This study aims to investigate a new feature extracted from near-ground hyperspectral imaging data to estimate precisely the LNC of wheat. In field experiments conducted over two years, we collected hyperspectral images at different rates of nitrogen and planting densities for several varieties of wheat throughout the growing season. We used traditional methods of classification (one unsupervised and one supervised method), spectral analysis (SA), textural analysis (TA), and integrated spectral and textural analysis (S-TA) to classify the images obtained as those of soil, panicles, sunlit leaves (SL), and shadowed leaves (SHL). The results show that the S-TA can provide a reasonable compromise between accuracy and efficiency (overall accuracy = 97.8%, Kappa coefficient = 0.971, and run time = 14 min), so the comparative results from S-TA were used to generate four target objects: the whole image (WI), all leaves (AL), SL, and SHL. Then, those objects were used to determine the relationships between the LNC and three types of indices: spectral indices (SIs), textural indices (TIs), and spectral and textural indices (STIs). All AL-derived indices achieved more stable relationships with the LNC than the WI-, SL-, and SHL-derived indices, and the AL-derived STI was the best index for estimating the LNC in terms of both calibration (${\mathrm{R}}_{\mathrm{c}}^2$ = 0.78, relative root mean-squared error (RRMSE_c) = 13.5%) and validation (${\mathrm{R}}_{\mathrm{v}}^2$ = 0.83, RRMSE_v = 10.9%). It suggests that extracting the spectral and textural features of all leaves from near-ground hyperspectral images can precisely estimate the LNC of wheat throughout the growing season. The workflow is promising for the LNC estimation of other crops and could be helpful for precision agriculture.

1. Introduction

Nitrogen is a crucial nutrient in crop growth that helps optimize fertilization management [1,2]. The reasonable application of nitrogen as fertilizer can not only improve the efficiency of its use as well as the yield and quality of the crop, but it can also reduce resource waste and environmental pollution [3]. Leaf nitrogen content (LNC) is an important indicator of the application of nitrogen as fertilizer in the early stages of crop growth [4,5], and it is intimately related to the quality of the final grain in later stages of growth [6,7]. Therefore, the accurate and timely quantification of LNC is a prerequisite for efficient agricultural management and beneficial for minimizing the environmental impact of farming [4,8].

Proximal sensing has emerged as an efficient technique in assessing precision agricultural. In particular, hyperspectral remote sensing has proved to be a promising approach in diagnosing the nitrogen content of wheat [9,10]. Hyperspectral remote sensing can be grouped into two categories: non-imaging and imaging spectrometry. A large amount of non-imaging spectral data have been used to monitor the nitrogen content of wheat [8,11]. However, the canopy spectra from non-imaging spectrometers inevitably contain a mixture of reflectance signals of crop organs and the soil as background to fields of wheat, especially in the early stages of growth (e.g., the greening and jointing stages) or in low-density wheat fields [10]. In fact, most of the studies based on non-imaging spectrometers focused on estimating nitrogen status in a dense canopy, such as after jointing stage or in high-density fields [8,9,12]. To reduce the impact of the background on the estimation of the nitrogen content of crops, many vegetation indices (VIs) have been developed [3,10], but the effect of the background has proved to be resilient to them [13]. Near-ground imaging spectrometers that combine images and spectra can provide data at very high spatial and spectral resolutions to improve the monitoring of crop growth [14]. Many studies have shown that near-ground imaging spectrometers are promising for assessing the nitrogen content of such crops as wheat [14], sugarcane [15], beet [16], and rice [17,18].

To estimate the nitrogen content of crops based on near-ground hyperspectral systems, previous studies have mainly focused on observations of well-illuminated leaves (i.e., sunlit leaves) in the field [14,16]. The crop canopy is not evenly illuminated by solar radiation [19]. Typically, leaves under direct sunlight constitute sunlit leaves (SL), and the blocking of direct sunlight onto lower leaves results in shaded leaves (SHL) [20]. Moreover, portions of the sunlit leaves (SL) and shaded leaves (SHL) change with the structure of the canopy of crop types, the solar angle, stages of crop growth, and stress-related conditions [21]. Ignoring signals of SHL might lead to the incorrect estimation of the nitrogen content of crops [18]. For example, the lower leaves (usually SHL) first respond to the early stage of nitrogen deficiency, but changes in the upper leaves (usually SL) are not prominent at this time, which leads to a deviation in the monitoring of the nitrogen content of crops [22]. Few studies have focused on determining the sunlit and shaded components of the canopy. Banerjee et al. [21] used thermal and visible cameras to identify the SL, SHL, and sunlit and shaded soils of wheat by using the support vector machine (SVM) before the heading stage, but they did not consider wheat panicles. Panicles gradually emerge from the sheath in canopies at advanced stages of reproduction [23]. Given that the reflectance signals of the canopy become more complex with the coexistence of leaves and panicles in case of wheat, uncertainties in the estimation of nitrogen content from the reflectance spectra of the canopy increase [18]. Zhou et al. [20] used a near-ground imaging spectrometer to separate the SL, SHL, and sunlit and shaded panicles by a decision tree. The results showed the potential of hyperspectral imagery for classifying the components of the canopy [20]. However, the method they used is empirical, and their findings might not be feasible for classifying wheat canopies and monitoring the nitrogen content of wheat owing to substantial differences in the structures of the canopies of rice and wheat.

In addition to spectral information, the spatial information (e.g., the form of texture) in hyperspectral imagery is worth considering for image classification. Image texture gives us information about the spatial arrangement of colors or intensities in an image or a selected region of it [24]. The applications of texture extraction to remote sensing image classifications can be traced back to the 1970s [25,26]. Since then, texture analysis has been widely used in monospectral and multispectral imagery based on data from satellites and unmanned aerial vehicles (UAVs) [27,28]. Although few studies have examined the performance of textural analysis for hyperspectral image classification, it has significant potential for use in the hyperspectral domain. For example, Ashoori et al. [29] evaluated the usefulness of textural measures for classifying crop types by using Earth Observing-1 (EO-1) Hyperion data, and they found that the overall accuracy improved by up to 7%. The results of images from a Visible/Infrared Imaging Spectrometer (AVIRIS) also yielded better classification based on textural analysis [30,31]. Zhang et al. [32] combined textural features with spectral indices from a pushbroom hyperspectral imager (PHI) to significantly improve the accuracy of crop classification and reduce edge effects. However, on the one hand, the spatial resolutions of hyperspectral sensors used in previous studies ranged from 1.2 m (PHI) to 30 m (EO-1 Hyperion); one the other hand, previous studies on crop classification have focused on identifying specific vegetation rather than components of the crop canopy. Therefore, the performance of near-ground hyperspectral images in terms of texture analysis in the context of the classification of the components of the crop canopy (e.g., SL and SHL) is promising but remains unexplored.

Based on image classification, crop growth monitoring can be further improved by selecting optimal spectral or/and textural indices, such as for the estimations of biomass [33,34], leaf area index (LAI) [35], and nitrogen content [28]. Dube and Mutanga [36] found that textural metrics could be used to characterize patterns of the structure of the canopy to enhance the sensitivity of Landsat-8 OLI data, especially under dense canopy coverage. Moreover, textural metrics can capture variations within and between rice canopies during early growth due to the consumption of nitrogen for vegetative development and its top dressing as fertilizer [28,37,38]. The application of the texture analysis was proposed to detect and quantify leaf nitrogen deficiency in maize [39,40]. However, little is known about the performance of texture analysis for near-ground hyperspectral images for estimating the LNC of wheat. Moreover, spectral indices (i.e., VIs) can be affected by the structure of the canopy of vegetation, canopy composition, and background and shadowing [41,42], and they often lose their sensitivity during stages of reproductive growth [43,44]. Therefore, making full use of near-ground hyperspectral data for monitoring the nitrogen content of crops by combining spectral and textural analysis offers considerable promise. However, the potential of spectral and textural indices from different objects (e.g., SL and SHL) for LNC estimation remains unclear.

To address the above-mentioned gaps in research, we investigated a new feature extracted from near-ground hyperspectral imaging data to estimate precisely the LNC of wheat. The objectives of this study are as follows: (1) to classify near-ground hyperspectral images into soil and wheat canopies (i.e., SL, SHL, and panicles); (2) to analyze the spectral variations in whole-image, all-leaves, SL, and SHL pixels under different growth conditions; and (3) to assess the potential of spectral and textural indices for estimating the LNC of wheat.

2. Materials and Methods

2.1. Experimental Design and LNC Measurements

Two field experiments were conducted over two growing seasons of winter wheat (2012–2013 and 2013–2014) in an experimental station of the National Engineering and Technology Center for Information Agriculture (NETCIA) located in Rugao city (120°45′ E, 32°16′ N), Jiangsu province, China. The annual regional precipitation and average annual temperature of this area were around 927.53 mm and 16.59 °C, respectively. The two experiments were designed using two varieties of wheat, two planting densities, and three rates of nitrogen application (

Table 1. The design of the two field trials.

Trials	Year	Variety	Planting Density (plants/ha)	Nitrogen Rate (kg/ha)	Sowing Date	Sampling Stage
Exp. 1	2012–2013	V1: Yangmai 18 V2: Shengxuan 6	D1: 150 D2: 300	N0: 0 N1: 150 N2: 300	5 November 2012	Greening Jointing Booting Heading
Exp. 2	2013–2014	V1: Yangmai 18 V2: Shengxuan 6	D1: 150 D2: 300	N0: 0 N1: 150 N2: 300	28 October 2013	Greening Jointing Booting Heading

). We arranged the plots randomly with three replications, each with an area of 35 m². Nitrogen fertilizers at the three different rates were applied, 50% as basal fertilizer on the sowing day and 50% at the jointing stage. The other 50% of the basal fertilizer included 120 kg/ha of P₂O₅ and 135 kg/ha of K₂O. All other agronomic management was undertaken according to local wheat production practices.

We measured the LNC at critical stages of growth (from greening to heading) of winter wheat (Table 1). Note that an unknown LNC variation within plot was ignored given the same treatment. We collected 6–10 samples of wheat plants from each plot and picked their leaves, and the growth status of the selected samples was similar. All leaves were oven-dried at 105 °C for 30 min and then at 80 °C until they had attained a constant weight. The dried leaf samples were ground to pass through a 1 mm screen and stored in plastic bags for chemical analysis. The LNC was determined based on the micro-Keldjahl method [45] using the SEAL AutoAnalyzer 3 HR (SEAL Analytical, Ltd., German).

2.2. Near-Ground Hyperspectral Imagery

2.2.1. Data Acquisition

The hyperspectral images were acquired by using a near-ground hyperspectral imaging system (Figure 1). The system was featured an electric forklift, a guide rail, a guide rail trolley, a hyperspectral imaging spectrometer, and an external laptop (Figure 1a). The hyperspectral imaging spectrometer (ImSpector V10E-PS, SPECIM, Oulu, Finland) was a pushbroom scanning sensor with a field of view of 42.8°. It contained 520 bands in the visible and near-infrared (NIR) regions (360–1025 nm) at a spectral resolution of 2.8 nm. The sensor was mounted on the guide rail, which could be lifted to a maximum height of 3 m from the ground. We set a constant height of 1 m between the lens and the top of the canopy during the growing season. A barium sulfate (BaSO4) panel was placed on a tripod as a white reference (Figure 1a) for radiometric calibration. Hyperspectral images were collected at noon in clear weather under stable light conditions in each sampling stage (Table 1).

2.2.2. Data Preprocessing

The preprocessing of the hyperspectral data consisted of five steps: (1) light intensity correction, (2) dark current correction of the sensor, (3) noise correction, (4) radiometric correction, and (5) selection of sample images. First, the raw images were processed by FODIS (fiber optic downwelling irradiance sensor) correction. FODIS recorded the light intensity of each scan based on the manner of imaging of the spectrometer V10E-PS. The light intensity of raw images was corrected by the software SpecView (SPECIM, Oulu, Finland). Second, dark current correction was used to correct errors in the system. The systematic error is the digital number (DN) value of the hyperspectral image scanned using a covered lens. This error could be directly subtracted in SpecView. Third, the hyperspectral data were smoothed using the minimum noise fraction (MNF) transform in ENVI 5.1 (EXELIS, Boulder, CO, USA). Fourth, the final values of relative reflectance were converted from the DN values, which were converted in turn into radiance:

$$Rad=a\times DN+b,$$

(1)

where Rad, a, and b are the radiance (W·cm⁻²·sr⁻¹·nm⁻¹), gain, and bias of the sensor, respectively. The values of a and b were automatically acquired by SpecView during hyperspectral imaging. The relative reflectance values were obtained by:

$${Ref}_{target}=\frac{{Rad}_{target}}{Ra{d}_{panel}}\times Re{f}_{panel},$$

(2)

where Rad_target/Ref_target and Rad_panel/Ref_panel are the radiance/reflectance values of the target images and the calibration panel (Figure 1a), respectively. Fifth, we selected only images in which the wheat was growing evenly and not overexposed. To avoid BRDF (bi-directional reflectance distribution function) effects along the edges, we cropped these images at the central area (i.e., the middle half of the original images).

2.3. Image Classification

We used different methods to classify near-ground hyperspectral images into four classes, i.e., soil, sunlit leaf (SL), shaded leaf (SHL), and panicle. For each class, we selected thousands of the regions of interests (ROIs) during all growth stages of both Exp. 1 and Exp. 2 using ENVI 5.1 (EXELIS, Boulder, CO, USA). The ROIs of Exp. 1 were used as the training set, and the ROIs of Exp. 2 were used as the testing set to calculate the classification accuracy. Accordingly, we could select the best method for image classification (Figure 2a).

2.3.1. Spectral Analysis (SA)

An advantage of hyperspectral data is the high spectral resolution. Based on the spectral curves of different classes, we could find an optimal spectral combination and corresponding threshold for separating SL, SHL, soil, and panicles [20]. The average reflectance of each class is shown in Figure 3. The four classes had prominent differences in spectral characteristics in the green (540–570 nm), red (670–700 nm), and NIR (760–900 nm) regions. In general, trends of the curves of wheat canopy (SL, SHL, and panicle) were similar. They featured a reflection peak in the green region and a reflection valley in the red region. On the contrary, the soil had no significant reflection peak or valley. In the green region, the reflectance of SL was the highest at 553 nm, while that of SHL and soil were the lowest. In the red region, soil had the highest reflectance and SHL had the lowest reflectance at 680 nm. In the NIR region, the reflectance of SL was much higher than that of the other classes, and these differences were significant. According to these characteristics, the reflectance at 553 nm (R₅₅₃) was subtracted from that at 680 nm (R₆₈₀), so that soil with non-positive values could be separated from the wheat canopy with positive values. To increase the differences among different classes, the resulting values were multiplied by the reflectance at 770 nm (R₇₇₀). Thus, we constructed a spectral classification index (SCI) for image classification:

$${\mathrm{SCI}=(\mathrm{R}}_{553}{-\mathrm{R}}_{680})\times {\mathrm{R}}_{770}\times 100.$$

(3)

2.3.2. Textural Analysis (TA)

In addition to high spectral resolution, hyperspectral images also have high spatial resolution. We could see the spatial details of the soil, wheat leaves, and panicles in the images acquired. For image classification, we focused on the external structures of these features. Mathematical morphology (MM), a theory and technique for the analysis and processing of geometrical structures, forms the foundation of morphological image processing [46]. It has been widely used in remotely sensed image classification [47,48]. Here, we used MM to extract different image classes. First, edge features were extracted from the processed hyperspectral imagery. Second, morphological filtering was used to remove noise. Third, the basic operators of binary morphology were applied, including dilation, erosion, and opening and closing operators. All processes were performed in MATLAB (V2016a, MathWorks, Natick, MA, USA), and all parameters of morphological operations functions were set as default. The binary images were finally obtained.

2.3.3. Comparison and Evaluation

Based on the intersection of the optimal results from the spectral analysis (SA) and the textural analysis (TA), we combined the two methods as the integrated spectral and textural analysis (S-TA) to separate different classes. In addition, we employed unsupervised (Iterative Self-Organizing Data Analysis Technique Algorithm, ISODATA) and supervised (maximum likelihood, MLE) classification methods as the benchmarks, since both of ISODATA and MLE are classical methods for remotely sensed image classification [49,50,51,52]. The overall accuracy and kappa coefficient were used to evaluate classification accuracy. The overall accuracy (OA) was the ratio of correctly classified samples to the total number of samples [53]:

$$\mathrm{OA}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^r{x}_i},$$

(4)

where N is the total number of samples, r is the number of classification categories, and x_ii is the number of discriminant samples for samples in the ith category. A confusion matrix is often used to describe virtually the performance of a classification method. All correct predictions are located in the diagonal of the matrix. The kappa coefficient (K) is a more robust measure than simple percent agreement calculation, as it considers the possibility of the agreement occurring by chance. The calculation of K is based on the confusion matrix [53]:

$$\mathrm{K}=\frac{N{\displaystyle {\sum }_{i=1}^r{x}_{ii}}-{\displaystyle {\sum }_{i=1}^r(x_{i+}\cdot x_{+i})}}{N^2-{\displaystyle {\sum }_{i=1}^r(x_{i+}\cdot x_{+i})}},$$

(5)

where x_i+ and x_+i are the sums of the ith row and the ith column of the confusion matrix, respectively, and K ranges from −1 to 1. A larger value of K indicates a better accuracy of image classification. In addition, we used the median run time of different classification methods to evaluate their efficiency by repeating three times for each data set (i.e., Exp. 1 and Exp. 2).

2.4. LNC Estimation

Based on the results of image classification, different classes were extracted from the hyperspectral images. We investigated four target objects for LNC estimation: the whole image (WI), the image of all leaves excluding soil and panicles (AL), the image of sunlit leaf only (SL), and the image of shaded leaf only (SHL) (Figure 2b).

2.4.1. Selection of Characteristic Wavelengths

Owing to the large redundancy in the hyperspectral images, it was necessary to extract the characteristic wavelengths. Note that different feature selection methods would lead to different results of characteristic wavelengths. In this study, we used the successive projections’ algorithm (SPA) to select the characteristic wavelengths for the LNC. The SPA is a forward selection method and a deterministic search technique [54] that uses simple operations in a vector space to minimize variable collinearity. The results of variable selection of the SPA are reproducible, rendering it more robust to the selection of the validation set. It has been widely used for the selection of characteristic wavelengths from hyperspectral data [55]. We applied the SPA in MATLAB (V2016a, MathWorks, Natick, MA, USA) and set the number of characteristics in the range of 5–30. It began with one wavelength and incorporated a new one at each iteration until a specified number N of wavelengths had been reached. The final number of characteristics was determined by root mean-squared error cross-validation (RMSECV). The interested reader can refer to [54] for steps of the SPA, and the given version of the graphical user interface for MATLAB can be downloaded from http://www.ele.ita.br/~kawakami/spa/ (accessed on 5 January 2021).

2.4.2. Calculation of Textural Features

The gray-level co-occurrence matrix (GLCM), the most commonly used texture algorithm, was used to calculate textural features of the four objects (i.e., WI, AL, SL, and SHL). Eight textural features were used: mean (ME), variance (VA), correlation (CC), homogeneity (HO), contrast (CO), dissimilarity (DI), entropy (EN), and second moment (SM). These features reflected the spatial distribution of a textural image (

Table 2. Textural features extracted from gray-level co-occurrence matrix (GLCM) and their representations.

Textural Feature	Formula 1	Description	Example
Mean (ME)	${\mu }_i={\displaystyle {\sum }_{i,j=0}^{n-1}i\cdot P_{i,j}}$ ${\mu }_j={\displaystyle {\sum }_{i,j=0}^{n-1}j\cdot P_{i,j}}$	A measure of the degree of texture rules
Variance (VA)	${\sigma }_i^2={\displaystyle {\sum }_{i,j=0}^{n-1}{P}_{i,j}}{(i-{\mu }_i)}^2$ ${\sigma }_j^2={\displaystyle {\sum }_{i,j=0}^{n-1}{P}_{i,j}}{(j-{\mu }_j)}^2$	A measure of the dispersion of the values around the mean
Correlation (CC)	${\displaystyle {\sum }_{i,j=0}^{n-1}{P}_{i,j}}\left[\frac{(i-{\mu }_i)(j-{\mu }_j)}{{\sigma }_i{\sigma }_j}\right]$	A measure of how correlated a pixel is to its neighbor over the whole image
Homogeneity (HO)	${\displaystyle {\sum }_{i,j=0}^{n-1}\frac{P_{i,j}}{1+\left|i-j\right|}}$	A measure of the local uniformity of a pixel pair
Contrast (CO)	${\displaystyle {\sum }_{i,j=0}^{n-1}{(i-j)}^2{P}_{i,j}}$	A measure of the local variations presented in an image
Dissimilarity (DI)	${\displaystyle {\sum }_{i,j=0}^{n-1}{P}_{i,j}\left|i-j\right|}$	A measure of the mean of the gray-level difference distribution of an image
Entropy (EN)	${\displaystyle {\sum }_{i,j=0}^{n-1}{P}_{i,j}(-\ln P_{i,j})}$	A measure of randomness degree of an image texture
Second moment (SM)	${\displaystyle {\sum }_{i,j=0}^{n-1}i\cdot P_{i,j}{}^2}$	A measure of global homogeneity of an image

¹P_i,j is the probability value from GLCM.

). To avoid complexity of computation, we calculated the eight textural features of the four objects at the selected characteristic wavelengths for LNC estimation. The features were calculated in ENVI 5.1 (EXELIS, Boulder, CO, USA) with the smallest window size (3 × 3 pixels).

2.4.3. Spectral and Textural Indices

We used three typical forms of two-band combinations for both spectral and textural indices, including the normalized difference index (NDI), ratio index (RI), and difference index (DI):

$$\mathrm{NDI}=\frac{x_{\lambda 1}-x_{\lambda 2}}{x_{\lambda 1}+x_{\lambda 2}},$$

(6)

$$\mathrm{RI}=\frac{x_{\lambda 1}}{x_{\lambda 2}},$$

(7)

$$\mathrm{DI}=x_{\lambda 1}-x_{\lambda 2},$$

(8)

where x is reflectance at a characteristic wavelength λ for spectral indices, or it is the textural feature (Table 2) at a characteristic wavelength λ for textural indices.

2.4.4. Estimation and Validation

Based on the NDI, RI, and DI, we built linear models of the LNC for the four objects (i.e., WI, AL, SL, and SHL) during all stages of wheat growth. Each of the estimated models were calibrated based on data from Exp. 1 and were validated using data from Exp. 2. The accuracy of the LNC models was evaluated by the coefficient of determination (R²) and the relative root mean-squared error (RRMSE):

$${\mathrm{R}}^2=\frac{\left[{\displaystyle {\sum }_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}\right]}{{\displaystyle {\sum }_{i=1}^n{(x_i-\overline{x})}^2}{\displaystyle {\sum }_{i=1}^n{(y_i-\overline{y})}^2}},$$

(9)

$$\mathrm{RRMSE}(\%)=\frac{100}{\overline{x}}\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{(y_i-x_i)}^2}},$$

(10)

where x, y, and n are the estimated LNC value, LNC measurements, and the number of samples, respectively.

3. Results

3.1. Near-Ground Hyperspectral Image Classification

Figure 4 shows the SCI frequency of the four classes (i.e., soil, panicle, SL, and SHL) in different stages. The ranges of the SCI values of soil, SL, and SHL were clearly demarcated in all stages. According to their ranges for the target components (

Table 3. The maximum and minimum SCI values of the target components at different stages. The values in the last column (all stages) are the corresponding minimum or maximum values in a row.

Target Component	Threshold	Greening	Jointing	Booting	Heading	All Stages
Soil	Min	−2.663	−3.742	−2.438	−1.414	−3.742
	Max	0.072	0.272	0.108	0.188	0.272
SHL	Min	0.303	0.602	0.596	0.522	0.303
	Max	1.411	2.009	2.063	2.072	2.072
SL	Min	2.203	2.215	3.250	2.980	2.203
	Max	9.795	9.865	9.794	9.444	9.865
Panicle	Min	—	—	—	−3.886	−3.886
	Max	—	—	—	8.104	8.104

), the soil, SHL, and SL could be separated as soil < 0.29 < SHL < 2.1 < SL. However, panicles could not be distinguished from the other components by SCI (Figure 4d and Table 3). By contrast, textural analysis based on the MM was able to extract panicles in the heading stage (Figure 5).

To further analyze the accuracy of different classification methods for each target class, a confusion matrix in the heading stage was shown in

Table 4. The confusion matrix (%) of each class based on different classification methods in the heading stage.

Target Class	Method	SL	SHL	Panicle	Soil
SL	ISODATA	96.58	0	3.42	0
	MLE	98.49	0	1.51	0
	SA	91.06	0	8.94	0
	TA	80.43	17.54	2.03	0
	S-TA	99.34	0	0.66	0
SHL	ISODATA	0.40	27.57	34.07	37.97
	MLE	0	90.09	0	9.91
	SA	0	89.20	10.80	0
	TA	6.23	93.37	0.40	0
	S-TA	0	98.92	1.08	0
Panicle	ISODATA	52.35	35.81	11.85	0
	MLE	2.32	5.81	91.86	0
	SA	18.86	30.35	48.50	2.28
	TA	7.38	0.72	91.91	0
	S-TA	0	0	100	0
Soil	ISODATA	0	0	0.21	99.79
	MLE	0	0.02	0	99.98
	SA	0	0	0	100
	TA	0	0	0	100
	S-TA	0	0	0	100

, since only the imagery in the heading stage contained all the four targets (especially panicles). For the target class SL, all methods classified part of it incorrectly as a panicle. The highest percentage of SL being incorrectly classified as panicles was 8.94% using the SA, and the lowest was 0.66% using the S-TA. In addition to being misclassified as panicles, 17.54% of SL was misclassified as SHL by the TA. For the target class SHL, ISODATA delivered the poorest performance, with the highest misclassification of 72.43%, including 37.97% as soil, 34.07% as panicle, and 0.40% as SL. On the contrary, the S-TA delivered the best performance because it recorded the lowest rate of misclassification of the SHL (1.08%). For the panicle as target class, most methods misclassified it as a wheat leaf (SL and SHL), but the S-TA correctly extracted all panicles. For soil as the target class, most methods were able to identify all soil pixels from the hyperspectral images, but the ISODATA and MLE misclassified 0.21% and 0.02% of the soils as panicles, respectively.

In general, the classification accuracy of the S-TA in all stages was the highest (OA = 97.8%, K = 0.971), which is followed by the MLE (OA = 94.5%, K = 0.942) (

Table 5. The classification accuracy and efficiency of different methods for all stages.

Method	Overall Accuracy (%)	Kappa Coefficient	Run Time (min)
ISODATA	71.6	0.632	85
MLE	94.5	0.942	100
SA	89.4	0.866	6
TA	89.6	0.868	12
S-TA	97.8	0.971	14

). However, the median run time of the MLE was more than six times longer than that of the S-TA. The unsupervised method ISODATA was the least efficient, with a value of OA of 71.6%, K = 0.632, and run time of 85 min.

3.2. Changes in Reflectance of Different Objects under Different Conditions

Based on the results of image classification, we compared the averaged changes in the reflectance of the four objects (i.e., WI, AL, SL, and SHL) under different conditions (Figure 6). In general, the reflectance of WI increased, and gaps among the reflectance of AL, SL, and SHL decreased when the nitrogen rate increased from N0 to N1 (Figure 6a–c) and planting density increased from D1 to D2 (Figure 6b,e). This is because the canopy coverage increased and the canopy spectra were less affected by the background of soil. Moreover, the spectra of different varieties of wheat were different (Figure 6b,d), which might have been related to the type of wheat plant. Yangmai 8 (V1) is erect and thus was more affected by the background of soil than the spread type (V2: Shengxuan 6). Regardless of the different conditions, the spectra of AL with the background of soil removed had a blue shift in the red edge region. This indicates that the soil as background would cause a red edge shift of the spectra of the wheat canopy.

3.3. LNC Estimation Using Spectral and Texture Indices of Different Objects

The characteristic wavelengths of different objects were extracted by the SPA and sorted by sensitivity to the LNC (

Table 6. The characteristic wavelengths of different objects using the successive projections algorithm (SPA).

Object	Characteristic Wavelength (nm) 1
WI	537.46, 479.26, 732.78, 400.49, 931.88, 945.34
AL	688.72, 758.28, 810.48, 931.88, 432.28, 664.39
SL	723.46, 457.03, 690.16, 400.49, 945.21
SHL	759.68, 696.42, 550.36, 840.38, 412.49, 933.21
WI	537.46, 479.26, 732.78, 400.49, 931.88, 945.34

¹ The characteristic wavelengths of different objects are sorted by their sensitivity to LNC.

). Based on Equations (6)–(8), we constructed spectral and textural indices of the four objects (i.e., WI, AL, SL, and SHL). By comparing all combinations of the characteristic wavelengths or textural features at these wavelengths, we selected the optimal indices (spectral index (SI), textural index (TI), and spectral and textural indices (STI)) for each object (

Table 7. The optimal spectral index (SI), textural index (TI), and spectral and textural index (STI) of the four objects, and the corresponding statistics for LNC estimations and validations.

Object	Index 1	${\mathbf{R}}_{\mathbf{c}}^2$	RRMSEc (%)	${\mathbf{R}}_{\mathbf{v}}^2$	RRMSEv (%)
WI	RSI (R945.34, R537.46)	0.64	16.3	0.54	17.8
	NDTI (T400.49ME, T945.34EN)	0.40	18.4	0.41	23.2
	NDSTI (R931.88, T931.88EN)	0.66	15.9	0.67	14.4
AL	DSI (R688.72, R432.38)	0.75	13.9	0.66	16.2
	RTI (T810.48HO, T931.88ME)	0.70	14.4	0.73	13.6
	NDSTI (R931.88, T432.38SM)	0.78	13.5	0.83	10.9
SL	NDSI (R945.21, R723.46)	0.50	15.2	0.49	18.0
	NDTI (T945.21DI, T690.16SM)	0.27	19.0	0.33	19.9
	NDSTI (R945.21, T945.21SM)	0.53	15.4	0.52	17.0
SHL	NDSI (R696.42, R412.49)	0.52	15.4	0.27	23.1
	NDTI (T933.21ME, T696.42SM)	0.34	17.9	0.39	21.8
	NDSTI (R933.21, T933.21SM)	0.49	16.1	0.46	17.5

¹R_λ is the reflectance at a characteristic wavelength λ (Table 5). T_λ is the textural feature (Table 2) at a characteristic wavelength λ (Table 5).

). In terms of their optimal composition, most indices were composed of one NIR wavelength and one visible wavelength. The most commonly selected textural feature was the second moment (SM), which was followed by the mean (ME). The optimal STIs of the four objects were the types of combination of the NDSTI with different spectral and textural features. For different objects, the accuracies of LNC calibration of SL and SHL were poor, with ${\mathrm{R}}_{\mathrm{c}}^2$ lower than 0.54. On the contrary, the LNC calibrations of AL were the best, with ${\mathrm{R}}_{\mathrm{c}}^2$ higher than 0.70. In terms of spectral and textural indices, the LNC~STI models were generally more accurate than the LNC~SI and LNC~TI models of each object. Although the accuracy of the LNC~SI models was not much lower than that of the LNC~STI models, the LNC~SI models were not robust. Conversely, LNC estimations based on the STI were more stable. For example, the accuracy of validation of the LNC~DSI model of AL decreased to values of ${\mathrm{R}}_{\mathrm{v}}^2$ of 0.66 and RRMSE of 13.9% (compared with an ${\mathrm{R}}_{\mathrm{c}}^2$ of 0.75 and RRMSE_c of 16.2%), while the ${\mathrm{R}}_{\mathrm{v}}^2$ of the LNC~NDSTI model was 0.83, and its RRMSE_v was 10.9%. Moreover, although the accuracy of the LNC~NDSI model during calibration was the highest for SHL with an ${\mathrm{R}}_{\mathrm{c}}^2$ of 0.52 and RRMSE_c of 15.4%, that in validation decreased to an ${\mathrm{R}}_{\mathrm{v}}^2$ of 0.27 and RRMSE_v of 23.1%. Overall, the best LNC estimations during both calibration and validation were obtained by the LNC~NDSTI (R_931.88, T_432.38SM) model of AL.

Based on the overall optimal indices of the NDSTI (R_931.88, T_432.38SM) of AL, Figure 7 shows the validation of LNC estimations under different conditions and in different stages. In general, the scatters were distributed along the 1:1 line. The LNC was accurately estimated, regardless of different rates of nitrogen (Figure 7a), planting densities (Figure 7b), wheat varieties (Figure 7c), and stages of growth (Figure 7d).

4. Discussion

4.1. Assessment of Classification Methods

We compared the performance of five classification methods: ISODATA, MLE, SA, TA, and S-TA. ISODATA is an unsupervised method of classification that could not distinguish different target components, and its overall classification accuracy and kappa coefficient were the lowest (Table 5). Compared with ISODATA, the overall accuracy and kappa coefficient of supervised classification (MLE) were 22.9% and 0.31 higher, respectively (Table 5). This shows that selecting appropriate and accurate training samples is crucial to the accuracy of classification. However, different kinds of image data require different training samples, and the use of supervised classification requires sophisticated hardware for data processing. In this study, the median run time of the MLE was the longest (100 min): 1.2 times that of ISODATA. The fastest classification method was the SA. It was constructed as a decision tree of the SCI based on the spectral characteristics of SL, SHL, panicles, and soil (Figure 2). Our results show that the differences in SCI among soil, SL, and SHL were clear (Figure 4), but panicles were easily misclassified as the other components by the SA (Table 4). This suggests that the SCI could not be used to extract wheat panicles, even though it was the optimal index according to our results. We might find a spectral index to distinguish leaves from panicles, but the limitations of the SA can affect the results of classification. A decision tree classification has a boundary effect, such that some pixels with values close to the threshold are misclassified [20].

Given that the spectral classifications (i.e., ISODATA, MLE, SA) are strongly dependent on differences in the spectral values of pixels [21], we investigated the TA based on the structure of spatial information of the images. Compared with the SA, the TA could extract panicles more accurately in the heading stage (Figure 5 and Table 4) because the spectral characteristics of wheat panicles are similar to those of leaves (Figure 2), but their textural characteristics are different. This suggests that the TA can improve classification accuracy by mitigating spectral confusion among spectrally similar classes [33]. However, the classification accuracy of the TA did not improve for all growth stages compared with that of the SA, and its overall accuracy and kappa coefficient were similar to those of the SA (Table 5). When we used only the TA, SL and SHL could not be distinguished (Table 4), owing to similar textural features of the leaves, but SA could help in this case. By combining the SA and TA, the S-TA boasts the complementary advantages of spectral and textural information. Moreover, compared with traditional classification methods, the S-TA was able to classify different components of winter wheat with higher accuracy, better robustness, and in shorter time (Table 5).

4.2. Soil Background and Shadow Effects in LNC Estimation

Soil as background affects the spectra of the crop canopy owing to the combined effects of planting density, canopy structure, and growth stage [56], which in turn affect LNC estimation [10]. The experimental design satisfactorily represented an environment in which the crop and soil coexisted in different conditions of planting density, wheat variety, and nitrogen rate during the greening, jointing, booting, and heading stages. Our results showed that the soil background would affect the spectral response of winter wheat regardless of different conditions (Figure 6). Furthermore, the correlations between LNC and the canopy indices were significantly affected by the soil background (Table 7), which is consistent with previous studies [10]. When excluding the soil background from WI, R² increased and RRMSE decreased for AL in both LNC calibration and validation (Table 7). Therefore, the soil background should be reduced to improve the LNC estimation of wheat.

In spite of the soil background, shadow effects can influence the accuracy of LNC estimation. However, previous studies have often focused on the sunlit pixels of the canopy for monitoring crops while neglecting shadow effects [57,58,59]. Our results showed low correlations between the indices and the LNC in SL without using SHL pixels (Table 7), which is mainly because the multiple scatterings among SHL would lead to a longer optical length and weaker signal strength [60,61] that would enhance the apparent absorptive characteristics of the leaves [62,63,64]. Zhou et al. [20] also calculated the gain from the SHL pixels, but for the estimation of chlorophyll content in rice, and they obtained higher accuracy when using SHL than when using SL. However, our results showed that the accuracies of LNC estimations in SHL and SL were similar, and the highest was in AL (Table 7). This might have occurred because LNC represents information on the wheat population at the leaf level, but SL pixels represent only parts of wheat leaves. Therefore, we suggest using the features of all leaves (including SL and SHL) to estimate the LNC of wheat.

4.3. Further Improvement of the LNC Estimation

Assessing the predictive performance of different indices showed that STIs provided higher accuracy of estimation of the LNC of wheat. In addition, the LNC estimation was the most stable when applied to the AL object regardless of different conditions (nitrogen rate, planting density, and wheat variety) in different growth stages (Figure 7). The premise is to eliminate the confounding effects of other components on the reflectance spectra of the canopy [18]. The TA could exclude soil background and panicles, thereby improving the LNC estimation. In contrast, the contribution of textural index to the LNC estimation was not significant (Table 7). The mechanistic link between LNC and textural index should be further investigated. The S-TA can provide a good compromise between accuracy and efficiency in extracting the features of all wheat leaves from the canopy image. Therefore, the integration of spectral and textural information offers a promising approach to monitor the nitrogen content of plants to enable precision farming throughout the growing season of winter wheat.

In particular, near-ground hyperspectral imagery with ultra-high spatial and spectral resolutions can provide more details. Thus, we should investigate the effects of STI selection, spatial resolution, and the estimated models on crop monitoring. When selecting optimal texture features, a moving window and image bands of appropriate sizes should be considered [65]. As the window size does not have a significant impact on experimental crop fields, where crops are often planted in rows [33], we used the recommended window size (3 × 3). However, the effect of window size should be investigated in future work, especially for large and heterogeneous fields. Moreover, textural features are computed from the spectral images, and so the corresponding spectral distribution is represented in image space [66]. Since the results of characteristic wavelengths might be different (or sometimes contradictory) from different feature selection methods, the characteristic wavelengths for LNC from the SPA in this study should be further investigated. We should compare the SPA to other methods (e.g., principal component analysis, particle swarm algorithm, maximum noise fraction, etc.) and investigate their impact on LAI estimations in the future work. Previous studies have focused on applications of spectral and textural analyses using satellite data (mostly Landsat TM) [27,67] and images captured by unmanned aerial vehicles (UAVs) (mostly digital and multispectral cameras) [28,33]. As the spatial resolution increases, more spatial details can be observed, but ultra-high spatial imagery often causes noise [68,69]. Therefore, different methods of sensitivity analysis and texture algorithms should be investigated. The impact of spatial resolution on image classification and crop monitoring should also be further assessed, and the efficiency of application can be improved at the optimal spatial resolution [18]. In addition to STI selection and spatial resolution, integrating spectral and textural information can help extend the potential of use of the estimation models. Although the simple statistical regression models of the LNC and STI showed good robustness, the potential for the use of large amounts of data needs to be further explored. In particular, machine learning might benefit LNC estimation as textural features provide a new dimension and enlarge the training data. In this regard, more prediction algorithms will be investigated in our future work.

5. Conclusions

We used near-ground hyperspectral data to estimate the LNC of winter wheat through the entire growing season in two-year field experiments under different nitrogen rates, planting densities, and wheat varieties. To classify the soil, panicles, sunlit leaves (SL), and shadowed leaves (SHL), the combination of spectral and textural analysis (S-TA) provided a good compromise between accuracy and efficiency (overall accuracy = 97.8%, kappa coefficient = 0.971, and run time = 14 min) compared with unsupervised classification (ISODATA), supervised classification (MLE), spectral analysis (SA), and textural analysis (TA). Based on the S-TA, we classified four objects—the whole image (WI), all leaves (AL), SL, and SHL—and extracted three indices (SI, TI, and STI) from them. All AL-derived indices achieved more stable relationships with the LNC than the WI-, SL-, and SHL-derived indices due to the elimination of the soil background and the consideration of shadow effects. Therefore, we suggest using the features of all leaves (including SL and SHL) to estimate LNC. In particular, the AL-derived NDSTI (R_931.88, T_432.38SM) delivered the best performance in terms of LNC estimations in both calibration (${\mathrm{R}}_{\mathrm{c}}^2$ = 0.78, RRMSE_c = 13.5%) and validation (${\mathrm{R}}_{\mathrm{v}}^2$ = 0.83, RRMSE_v = 10.9%). It suggests that the integration of spectral and textural information offers a promising approach to both image classification and crop monitoring. In the future work, the effects of STI selection, spatial resolution, and the estimated models need to be further investigated.