Evaluating Di ﬀ erent Non-Destructive Estimation Methods for Winter Wheat ( Triticum aestivum L.) Nitrogen Status Based on Canopy Spectrum

: Compared to conventional laboratory testing methods, crop nitrogen estimation methods based on canopy spectral characteristics have advantages in terms of timeliness, cost, and practicality. A variety of rapid and non-destructive estimation methods based on the canopy spectrum have been developed on the scale of space, sky, and ground. In order to understand the di ﬀ erences in estimation accuracy and applicability of these methods, as well as for the convenience of users to select the suitable technology, models for estimation of nitrogen status of winter wheat were developed and compared for three methods: drone equipped with a multispectral camera, soil plant analysis development (SPAD) chlorophyll meter, and smartphone photography. Based on the correlations between observed nitrogen status in winter wheat and related vegetation indices, green normalized di ﬀ erence vegetation index (GNDVI) and visible atmospherically resistant index (VARI) were selected as the sensitive vegetation indices for the drone equipped with a multispectral camera and smartphone photography methods, respectively. The correlation coe ﬃ cients between GNDVI, SPAD, and VARI were 0.92 ** and 0.89 **, and that between SPAD and VARI was 0.90 **, which indicated that three vegetation indices for these three estimation methods were signiﬁcantly related to each other. The determination coe ﬃ cients of the 0–90 cm soil nitrate nitrogen content estimation models for the drone equipped with a multispectral camera, SPAD, and smartphone photography methods were 0.63, 0.54, and 0.81, respectively. In the estimation accuracy evaluation, the method of smartphone photography had the smallest root mean square error (RMSE = 9.80 mg / kg). The accuracy of the smartphone photography method was slightly higher than the other two methods. Due to the limitations of these models, it was found that the crop nitrogen estimation methods based on canopy spectrum were not suitable for the crops under severe phosphate deﬁciency. In addition, in estimation of soil nitrate nitrogen content, there were saturation responses in the estimation indicators of the three methods. In order to introduce these three methods in the precise management of nitrogen fertilizer, it is necessary to further improve their estimation models.

complementarity. Therefore, the objective of this research was to compare the three methods of nondestructive estimation of winter wheat N status, namely, a drone equipped with a multispectral camera (hereinafter referred to as UAVMC), SPAD, and smartphone photography (hereinafter referred to as PHONEP) methods. By performing this inter-comparison, it is expected to clarify the applicability and limitations of these three methods in the precision management of N fertilizer for winter wheat and provide a basis for the improvement for these methods and simplify the selection of estimation methods for the producer.

Field Experiments
The fertilizer level experiment was conducted at the Luancheng Agro-Ecosystem Experimental Station (37.8897°N, 114.6935°E) of the Chinese Academy of Sciences. This area has loam soil and is situated in the temperate semi-humid monsoon climate. The main crop rotation system is winter wheat and maize (Zea mays L.), and all their straws after harvest are returned to the field. The experiment has been established since 1997. There were 4 levels of N fertilizer applications: 0, 200, 400, and 600 kg(N)·ha −1 ·a −1 (marked as N0, N1, N2, and N3, respectively), 3 levels of phosphorus (P): 0, 32.5, and 64 kg(P)·ha −1 ·a −1 (marked as P0, P1, and P2, respectively), and 2 levels of potassium (K): 0 and 50 kg(K)·ha −1 ·a −1 (marked as K0 and K2, respectively). 1/4 of the total amount of N fertilizer was applied as winter wheat base fertilizer, and 1/4 was applied as a top dressing at the jointing stage of winter wheat, with the other 1/2 being applied as a top dressing for maize. 1/2 of the total K fertilizer was applied as the base fertilizer for winter wheat, and the other 1/2 was applied as the top dressing for maize. All P fertilizers were applied as base fertilizer for winter wheat. A total of 16 treatments were selected from the orthogonal incomplete design. A randomized 48 plots (8 × 14 m 2 each) with 3 replications of each treatment were arranged (as shown in Figure 1 and Table 1). The variety of the winter wheat used was "Kelong199".      1  2  3  4  5  6  7  8  9  10  11  12   24  23  22  21  20  19  18  17  16  15  14  13   25  26  27  28  29  30  31  32  33  34  35  36   48  47  46  45  44  43  42  41  40  39 38 37

Aerial Photography of UAVMC and the Reference VIs
All of the experimental campaigns were conducted during the standing stage of the winter wheat on 12 April 2018 and 4 April 2019, before topdressing in the jointing stage. The weather conditions were all suitable for acquiring aerial photography (cloud-free without strong winds). Time was selected from Remote Sens. 2020, 12, 95 4 of 16 10:00 a.m. to 2:00 p.m. with a relatively high solar altitude. Parrot Sequoia was used as a multispectral camera in this study, which had two sensors: a multispectral sensor and a sunlight sensor, with a total weight of 105 g. Spectral band parameters of the Parrot Sequoia are shown in Table 2. The sensors were mounted on an unmanned aerial vehicle DJI Phantom 4 quadrotor drone to acquire aerial photography.
In order to cover all the experimental area with enough overlaps, the aerial routes were programmed by the drone ground station software (DJI Ground Station Pro, an iPad app). The flight altitude was set to 50 m, and the frontal and side overlaps were all set to 80%. Before each flight, the multispectral images of a calibration target provided by the manufacturer were recorded for the radiometric calibration in post-processing. After uploading the aerial route to the drone, the drone completed the programmed flight and the spectrum camera took the images automatically. Pix4D software was used for the post-processing of the multispectral images to obtain an entire image of the experimental area and its reflectance data for the 4 bands. Figure 2 shows the pseudo color multispectral image of the experimental area.

Aerial Photography of UAVMC and the Reference VIs
All of the experimental campaigns were conducted during the standing stage of the winter wheat on 12 April 2018 and 4 April 2019, before topdressing in the jointing stage. The weather conditions were all suitable for acquiring aerial photography (cloud-free without strong winds). Time was selected from 10:00 a.m. to 2:00 p.m. with a relatively high solar altitude. Parrot Sequoia was used as a multispectral camera in this study, which had two sensors: a multispectral sensor and a sunlight sensor, with a total weight of 105 g. Spectral band parameters of the Parrot Sequoia are shown in Table 2. The sensors were mounted on an unmanned aerial vehicle DJI Phantom 4 quadrotor drone to acquire aerial photography.
In order to cover all the experimental area with enough overlaps, the aerial routes were programmed by the drone ground station software (DJI Ground Station Pro, an iPad app). The flight altitude was set to 50 m, and the frontal and side overlaps were all set to 80%. Before each flight, the multispectral images of a calibration target provided by the manufacturer were recorded for the radiometric calibration in post-processing. After uploading the aerial route to the drone, the drone completed the programmed flight and the spectrum camera took the images automatically. Pix4D software was used for the post-processing of the multispectral images to obtain an entire image of the experimental area and its reflectance data for the 4 bands. Figure 2 shows the pseudo color multispectral image of the experimental area.  Referring to the related research in remote sensing monitoring of CNS, 11 commonly used VIs (see Table 3) that can be retrieved using the multispectral image were selected as reference VIs. The most sensitive VI will be selected from them to establish the estimation models for N status for winter wheat.  Referring to the related research in remote sensing monitoring of CNS, 11 commonly used VIs (see Table 3) that can be retrieved using the multispectral image were selected as reference VIs. The most sensitive VI will be selected from them to establish the estimation models for N status for winter wheat.

Name of VI Abbreviation Equation Reference
Normalized vegetation index Renormalized difference vegetation index RDVI RNDVI = (R NIR − R RED )/ (R NIR + R RED ) [34] Ratio vegetation index Where R RED , R GRE , and R NIR are reflectances of the red, green, and near-infrared bands, respectively.

Taking Photos of Winter Wheat with a Smartphone and the Reference Color-Based VIs
Studies have shown that there were no significant differences in the color parameters of the canopy photos obtained by different mobile phones [22]. Based on the previous research on the changes in canopy color parameters caused by different photographing methods, winter wheat was photographed by a smartphone at a height of 1 meter above the canopy and a 60 • angle to the ground [22]. Using automatic exposure mode, photos were taken at three different locations in each plot and saved as JPG format. Figure 3 shows the winter wheat canopy photos from plots with different N fertilizer levels. Where RRED, RGRE, and RNIR are reflectances of the red, green, and near-infrared bands, respectively.

Taking Photos of Winter Wheat with a Smartphone and the Reference Color-Based VIs
Studies have shown that there were no significant differences in the color parameters of the canopy photos obtained by different mobile phones [22]. Based on the previous research on the changes in canopy color parameters caused by different photographing methods, winter wheat was photographed by a smartphone at a height of 1 meter above the canopy and a 60° angle to the ground [22]. Using automatic exposure mode, photos were taken at three different locations in each plot and saved as JPG format. Figure 3 shows the winter wheat canopy photos from plots with different N fertilizer levels. A variety of color-based VIs can be obtained by combining calculations using the Red (R), Green (G), and Blue (B) color channel values of the crop canopy photo. Referring to the related research on the estimation of CNS using canopy photo, 10 color-based VIs were selected as reference VIs for the establishment of estimation models (Table 4). A variety of color-based VIs can be obtained by combining calculations using the Red (R), Green (G), and Blue (B) color channel values of the crop canopy photo. Referring to the related research on the estimation of CNS using canopy photo, 10 color-based VIs were selected as reference VIs for the establishment of estimation models (Table 4).
Since the photos of winter wheat in the field were taken at different time points, the values of the R, G, and B channels were susceptible to the intensity of the sun's illumination, suggesting that the absolute values of the single-color channels of the photos were not suitable for the estimation of CNS. To eliminate the effects of solar illumination changes, color-based VIs are often calculated using the ratio of different color channels or standardized form [37,38]. Therefore, this study normalized the original calculation equations of EXG and GMR in Table 4 by dividing (2G − R − B) and (G − R) by (R + G + B), respectively, in order to eliminate the influence of solar illumination changes. Table 4. Color-based VIs used in this study.

Name of VI Abbreviation Equation Reference
The dark green color index DGCI The difference between green and red GMR Where, R, G, and B are the values of the red, green, and blue channels, respectively.

Measurements of Nitrogen (N) Status of Winter Wheat
Conventional laboratory testing methods were adopted to measure the N status of winter wheat on 12 April 2018 and 4 April 2019. The above-ground part of winter wheat was sampled from 1 m 2 in each plot, and total N (TN) was measured using the Kjeldahl method. Winter wheat root zone soil layers in each plot were collected at the depths of 0-30 cm, 30-60 cm, and 60-90 cm, respectively. These soil samples were extracted with 1 mol·L −1 KCl and the nitrate nitrogen content was measured using the ultraviolet (UV) spectrometry method. 30 winter wheat samples were randomly selected in each plot. The SPAD value of the first fully expanded leaf of each sample was measured in the field with a SPAD-502 chlorophyll meter, and the average value was recorded as the SPAD value for this plot.

Analytical Methods
Considering the wide application of soil testing and formula fertilization technology in precision agriculture, the estimation models for soil nitrate nitrogen content were developed for the methods of UAVMC, SPAD, and PHONEP in this study respectively, with the nitrate nitrogen content of winter wheat root zone soil layers being the inversion targets.
For the UAVMC method, the first step was to remove pixels of bare soil from the multispectral images to avoid the contamination of bare ground on the analysis. To achieve this goal, the NDVI of this area was calculated using the multispectral image according to the equation listed in Table 3. Winter wheat and bare soil were distinguished to obtain the mask files for winter wheat in 2018 and 2019 based on the threshold for bare soil NDVI. The 11 VIs of winter wheat listed in Table 3 were calculated using multispectral images and the mask files for winter wheat. Correspondingly, the average value of these VIs for each plot was obtained using the map for the 48 plots. Based on the correlation analysis between the 11 VIs, TN of winter wheat, and the soil nitrate nitrogen content of each layer (excluding the data for the P0 plots, as for the other two methods), the VI with the largest correlation coefficient was selected to establish the estimation models for the inversion of soil nitrate nitrogen content in root layers.
For the SPAD method, the correlations between SPAD values, TN of winter wheat, and the soil nitrate nitrogen content of each layer were analyzed. Then, the estimation models were established based on the regression equations of SPAD values and soil nitrate nitrogen content in root layers. For the method of PHONEP, the analysis procedures were similar to the UAVMC method. To remove the non-vegetation pixels in each photo, SAVI Green and GMR were calculated using the R, G, and B of the winter wheat canopy photos referring to the equations in Table 4. According to the threshold (SAVI Green > 0 and GMR > 0), the leaf mask for each photo was obtained [44,45]. Using the corresponding leaf mask, the average values of 10 VIs of winter wheat in Table 4 were calculated for each canopy photo. There were 3 photos for each plot, and the average value of the 10 VIs in the 3 photos was taken as the value for this plot. Similarly, based on the correlations between the 10 VIs, the TN of winter wheat, and the soil nitrate nitrogen content for each layer, the VI with the largest correlation coefficient was selected to establish estimation models of soil nitrate nitrogen content in root layers.
All the data from 2018 and half of the data from 2019 (data of plots numbered 1 to 24 as shown in Figure 1) were used for the development of estimation models. The remaining data (data of plots numbered 25 to 48 in 2019) were used to evaluate the estimation accuracy of the three methods by the parameters of determination coefficient (R 2 ) and root mean square error (RMSE). The RMSE was calculated from: where n was the number of plots, andŷ, y i were the estimated and measured values, respectively. The correlations between the measured CNS (TN and nitrate nitrogen content at different soil layers) and reference VIs were analyzed using SPSS software. The calculation and spatial statistics of the VIs using multispectral images and canopy photos were implemented by the ERDAS IMAGINE software.

Variation of CNS in the Fertilizer Level Experiment
The N supply for winter wheat was routinely determined by the nitrate content in root zone soil from the jointing stage to the harvest of winter wheat (0-90 cm) [47]. This practice was usually used in order to determine the amount of N fertilizer needed. The level of TN for winter wheat could also reflect the N status of the crop and was closely related to the N supply capacity of the root zone soil. Therefore, the nitrate nitrogen content in the 0-90 cm soil layer and the TN of winter wheat were measured in this study. As shown in Figure 4, due to the water, fertilizer, and other managements that were identical every year in the variable fertilizer level experiment, the two parameters measured during the same period in 2018 and 2019 were very close to each other. The biggest inter-annual differences of the nitrate nitrogen content of 0-90 cm soil and the TN of winter wheat were 28.20 mg/kg and 0.52%, respectively.
The distributions of the two parameters were consistent with the setting of the N fertilizer levels. With the increase of N fertilizer application, both parameters showed an upward trend. It is worth noting that the two parameters for the plots where P was 0 (numbered 10 to 24 in Figure 1). The content of soil nitrate nitrogen showed a rapid increase, while the TN of winter wheat in plots 16 to 24 stayed at around 3.05%. The severe shortage of P fertilizer inhibited the absorption of N fertilizer by winter wheat, and the accumulation of N fertilizer in the soil made the content of nitrate nitrogen of these plots significantly higher than that of other plots. Therefore, the data from these P0 plots were excluded during analysis.

Estimation Models for the Method of UAVMC
The correlations between the 11 spectral VIs of winter wheat, TN of plant, and the nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm soil layers were analyzed, and all the correlation coefficients are listed in Table 5. All the 11 spectral VIs of winter wheat were significantly and positively correlated to the TN of their plants. Among these VIs, the correlation coefficient between GNDVI and TN of plants was the highest, reaching to 0.90 **. Although the correlation coefficient between RVI and TN of plants was the lowest, it still reached to 0.83 **. The N status of winter wheat was closely related to its spectral characteristics.
During the standing stage of winter wheat, 92.7% of its total root length was concentrated in the soil root layer of 0-40 cm [48]. Therefore, the N status of winter wheat should be closely related to the N supply capacity of the upper soil layer. The correlation coefficients between 11 spectral VIs and the 0-30 cm soil nitrate nitrogen content were the highest. With the increase of soil depth, the correlation coefficients decreased gradually. Among all the spectral VIs, the correlation coefficients between GNDVI and soil nitrate nitrogen content of the three layers were the highest. They were 0.52 **, 0.48 **, and 0.42 ** from the top to bottom layers. For the whole soil layer (0-90 cm), the correlation coefficient was 0.52 **.

Estimation Models for the Method of UAVMC
The correlations between the 11 spectral VIs of winter wheat, TN of plant, and the nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm soil layers were analyzed, and all the correlation coefficients are listed in Table 5. All the 11 spectral VIs of winter wheat were significantly and positively correlated to the TN of their plants. Among these VIs, the correlation coefficient between GNDVI and TN of plants was the highest, reaching to 0.90 **. Although the correlation coefficient between RVI and TN of plants was the lowest, it still reached to 0.83 **. The N status of winter wheat was closely related to its spectral characteristics.
During the standing stage of winter wheat, 92.7% of its total root length was concentrated in the soil root layer of 0-40 cm [48]. Therefore, the N status of winter wheat should be closely related to the N supply capacity of the upper soil layer. The correlation coefficients between 11 spectral VIs and the 0-30 cm soil nitrate nitrogen content were the highest. With the increase of soil depth, the correlation coefficients decreased gradually. Among all the spectral VIs, the correlation coefficients between GNDVI and soil nitrate nitrogen content of the three layers were the highest. They were 0.52 **, 0.48 **, and 0.42 ** from the top to bottom layers. For the whole soil layer (0-90 cm), the correlation coefficient was 0.52 **. Correlation analysis showed that GNDVI was the most sensitive spectral VI with N status of winter wheat. Therefore, GNDVI was selected as the estimation indicator, and the estimation models were established for the inversion of soil nitrate nitrogen content in different root depths of winter wheat. The models were: Estimation model for 0-30 cm:

Estimation Models for the SPAD Method
The correlation coefficients between the SPAD value of winter wheat leaves and TN of plants, nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm soil layers are shown in Table 6. SPAD and all winter wheat N status parameters were significantly correlated. The correlation coefficients between the SPAD value and soil nitrate nitrogen content decreased gradually with the increase of soil depth. Table 6. Correlation coefficients between SPAD value and N status of plant and soil. Taking the SPAD value as the estimation indicator, the estimation models for soil nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm layers were:

Estimation Models for the PHONEP Method
The correlation coefficients between 10 color-based VIs of winter wheat, TN of plants, and the nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm soil layers are shown in Table 7. All the color-based VIs were significantly correlated with the TN of the plants. Among them, EXG, GLI, NGI, and NRI were negatively correlated, and the rest were positively correlated. Affected by the distribution characteristics of winter wheat roots, the correlation coefficients between color-based VIs and soil nitrate nitrogen content decreased with the increase of depth. Among the 10 color-based VIs, VARI had the best correlation with TN of plants, and its correlation coefficient was 0.91 **. The correlation coefficients between VARI and the soil nitrate content in 0-30 cm, 30-60 cm, and 0-90 cm soil layers were the largest. In the soil layers 60-90 cm, the correlation coefficient of VARI was slightly lower than that of DGCI. Considering that the correlation coefficients between VARI and TN of plant, 0-30 cm, 30-60 cm, and 0-90 cm soil nitrate nitrogen content were the largest, VARI was selected as the most sensitive color-based VI to establish the estimation model. Similarly, VARI was also selected by other researchers to develop the crop N estimation models [22,48]. Table 7. Correlation coefficients between color-based VIs and nitrogen status of plant and soil. Taking the VARI as the estimation indicator, the estimation models for soil nitrate nitrogen content in 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm layers were:

Validation
In view of the guiding role of root zone soil nitrate nitrogen content in the precise management of N fertilizer, the accuracy of the estimation models for the 0-90 cm soil nitrate nitrogen content by using the methods of UAVMC, SPAD, and PHONEP were evaluated. Figure 5 shows the linear correlations between nitrate nitrogen content in the 0-90 cm soil obtained by the laboratory testing method and that estimated by the three methods. The PHONEP method had the best estimation accuracy (R 2 = 0.93 and RMSE = 9.80 mg/kg). The SPAD method had the lowest estimation accuracy (R 2 = 0.61 and RMSE = 19.80 mg/kg). For the method of UAVMC, R 2 is 0.86 and RMSE is 12.40 mg/kg. As shown in Figure 5, all three methods had relatively high estimation accuracy in the low-value areas of soil nitrate nitrogen content. While in the high-value areas of soil nitrate nitrogen content, the estimated values of the three methods were significantly lower than the measured values.

Comparison of the Three Estimation Methods
In this study, GNDVI, SPAD, and VARI were selected as estimation indicators for the three methods to establish estimation models for soil nitrate nitrogen content. The correlation coefficients between GNDVI, SPAD, and VARI were 0.92 ** and 0.89 ** respectively, and that between SPAD and VARI was 0.90 **. The significant correlation between these three indicators showed a high similarity of the three methods in the estimation of soil nitrate nitrogen content.
For the estimation models, the determination coefficients of the estimation equations for the PHONEP method were 0.82, 0.71, 0.67, and 0.81 for the soil layers of 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm respectively, which were all slightly higher than the other two methods. The determination coefficients and RMSEs in the estimation accuracy evaluation indicated that the PHONEP method was the best, while the SPAD method was the worst method in terms of the estimation accuracy.

Effect of P Fertilizer Shortage on CNS Estimation
Considering the synergistic effects of N, P, and K in the nutrient intake of crops, the severe shortage of P fertilizer could inhibit the absorption of N by crops, and therefore, the P0 plots were excluded from the analyses. In order to analyze the impact of P fertilizer shortage on the diagnosis of CNS, taking VARI as an example, the estimation models with and without the data from P0 plots were compared.
In the fertilizer level experiment, plots with P0 contained four levels of N fertilizer, N0, N1, N2, and N3. As shown in Figure 6a, the points with P0 are concentrated at the lower left of the figure. This finding indicated that even if the N fertilizer was sufficient (such as N3), the TN of winter wheat plant was still generally low. If the points with P0 were removed, the correlation between VARI and TN of winter wheat plants was significantly improved (Figure 6b). The smaller VARI was correlated with the lower nitrate nitrogen content in the soil layer of 0-90 cm as shown in Equation (13). However, in Figure 6c, there were some points with low VARI and high soil nitrate nitrogen content in the 0-90 cm soil layer. This could be attributed to the shortage of P fertilizer. When this happened, the NUE of winter wheat was reduced, leading to the retention of applied N in the soil. If these P0 points were excluded, the determination coefficient of the fitted equation for VARI and the 0-90 cm

Comparison of the Three Estimation Methods
In this study, GNDVI, SPAD, and VARI were selected as estimation indicators for the three methods to establish estimation models for soil nitrate nitrogen content. The correlation coefficients between GNDVI, SPAD, and VARI were 0.92 ** and 0.89 ** respectively, and that between SPAD and VARI was 0.90 **. The significant correlation between these three indicators showed a high similarity of the three methods in the estimation of soil nitrate nitrogen content.
For the estimation models, the determination coefficients of the estimation equations for the PHONEP method were 0.82, 0.71, 0.67, and 0.81 for the soil layers of 0-30 cm, 30-60 cm, 60-90 cm, and 0-90 cm respectively, which were all slightly higher than the other two methods. The determination coefficients and RMSEs in the estimation accuracy evaluation indicated that the PHONEP method was the best, while the SPAD method was the worst method in terms of the estimation accuracy.

Effect of P Fertilizer Shortage on CNS Estimation
Considering the synergistic effects of N, P, and K in the nutrient intake of crops, the severe shortage of P fertilizer could inhibit the absorption of N by crops, and therefore, the P0 plots were excluded from the analyses. In order to analyze the impact of P fertilizer shortage on the diagnosis of CNS, taking VARI as an example, the estimation models with and without the data from P0 plots were compared.
In the fertilizer level experiment, plots with P0 contained four levels of N fertilizer, N0, N1, N2, and N3. As shown in Figure 6a, the points with P0 are concentrated at the lower left of the figure. This finding indicated that even if the N fertilizer was sufficient (such as N3), the TN of winter wheat plant was still generally low. If the points with P0 were removed, the correlation between VARI and TN of winter wheat plants was significantly improved (Figure 6b). The smaller VARI was correlated with the lower nitrate nitrogen content in the soil layer of 0-90 cm as shown in Equation (13). However, in Figure 6c, there were some points with low VARI and high soil nitrate nitrogen content in the 0-90 cm soil layer. This could be attributed to the shortage of P fertilizer. When this happened, the NUE of winter wheat was reduced, leading to the retention of applied N in the soil. If these P0 points were excluded, the determination coefficient of the fitted equation for VARI and the 0-90 cm soil nitrate nitrogen content was significantly increased (Figure 6d). Results showed that the estimation method based on the canopy spectrum of winter wheat was not suitable for the estimation of N status of crops under a severe P fertilizer shortage. In actual production, the amount of P fertilizer applied was generally not too low to affect the absorption of N fertilizer by crops. Nevertheless, we should pay attention to the existence of this phenomenon. Due to the fact that more than 80% of the K absorbed by crops was stored in its straw, these K could be quickly released into the soil by straw returning [49]. In this way, plots with K0 in the fertilizer level experiment could get enough K to support the growth of winter wheat. Therefore, data from plots with K0 were involved in the development of estimation models.

The Saturation Response of the Estimation Indices
Some VIs (such as NDVI) might reach saturation in high-density vegetation areas [50], and the SPAD method had this limitation when N supply was sufficient for the crop as well [51]. Taking the measurements in this study as an example (excluding P0 plots), the saturation response of the estimation VIs (GNDVI, SPAD, and VARI) of the three methods in inversion of the 0-90 cm soil nitrate nitrogen content was analyzed and is shown in Figure 7. The distributions of the three VIs were similar to each other. In the low-value region of soil nitrate nitrogen content in the 0-90 cm depth, all the three VIs and the soil nitrate nitrogen content showed the quadratic polynomial correlations. With the increase of soil nitrate nitrogen content, the three VIs rapidly increased and Due to the fact that more than 80% of the K absorbed by crops was stored in its straw, these K could be quickly released into the soil by straw returning [49]. In this way, plots with K0 in the fertilizer level experiment could get enough K to support the growth of winter wheat. Therefore, data from plots with K0 were involved in the development of estimation models.

The Saturation Response of the Estimation Indices
Some VIs (such as NDVI) might reach saturation in high-density vegetation areas [50], and the SPAD method had this limitation when N supply was sufficient for the crop as well [51]. Taking the measurements in this study as an example (excluding P0 plots), the saturation response of the estimation VIs (GNDVI, SPAD, and VARI) of the three methods in inversion of the 0-90 cm soil nitrate nitrogen content was analyzed and is shown in Figure 7. The distributions of the three VIs were similar to each other. In the low-value region of soil nitrate nitrogen content in the 0-90 cm depth, all the three VIs and the soil nitrate nitrogen content showed the quadratic polynomial correlations. With the increase of soil nitrate nitrogen content, the three VIs rapidly increased and reached their maximum values, i.e., their saturation points. After the saturation points, the three VIs did not significantly increase with the increase of nitrate nitrogen content in the soil. As shown in Figure 7a,b, the values of GNDVI and SPAD were hardly changed after their saturation points. Figure 7c showed that the value of soil nitrate nitrogen content at the saturation point of VARI was greater than that of GNDVI and SPAD. Moreover, after its saturation point, with the increase of soil nitrate nitrogen content, VARI tended to increase slowly. This indicated that the PHONEP method was relatively better than the other two methods. Therefore, it is necessary to pay attention to the saturation response of VIs for the estimation of CNS, or to improve the estimation accuracy by using the method of segmentation.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 16 nitrogen content, VARI tended to increase slowly. This indicated that the PHONEP method was relatively better than the other two methods. Therefore, it is necessary to pay attention to the saturation response of VIs for the estimation of CNS, or to improve the estimation accuracy by using the method of segmentation.

Conclusions
Fertilizers play irreplaceable roles in agricultural production, making it possible to provide more food from limited arable land for a growing population. However, excessive application for fertilizers and low use efficiency have increased the production cost and environmental pollution. It is crucial to implement rational use of chemical fertilizers under the premise of ensuring substantial grain yields. In order to achieve this goal, practical crop nutrient estimation methods and fertilization techniques are needed, especially for the N fertilizer that has the largest application amount. Compared to the laboratory testing method, rapid non-destructive estimation methods of CNS based on the canopy spectrum have more advantages in terms of timeliness, cost, practicality, and thus, have received extensive attention and research.
Based on the fertilizer level experiment conducted during 2018 and 2019, GNDVI and VARI were selected as the sensitive VIs for the methods UAVMC and PHONEP, respectively. The correlation coefficients between GNDVI and SPAD, VARI were 0.92 ** and 0.89 **, and the correlation coefficient between SPAD and VARI was 0.90 **. The estimation VIs for the three methods were significantly correlated with each other. The determination coefficients of the 0-90 cm soil nitrate nitrogen content estimation models for the UAVMC, SPAD, and PHONEP methods were 0.63, 0.54, and 0.81, respectively. Moreover, the PHONEP method had the smallest RMSE (9.8 mg/kg) in the estimation accuracy evaluation. Therefore, in terms of estimation accuracy, the PHONEP method was slightly higher than the other two methods. In addition, it was found that CNS estimation methods based on canopy spectrum were not suitable for crops under severe P deficiency due to the synergistic effects between different nutrients. At the same time, the estimation VIs of the three methods had a saturation response in the estimation of soil nitrate nitrogen content. The method of UAVMC is suitable for rapid CNS estimation for large-area crops, and SPAD and PHONEP are convenient for the random sampling and monitoring of crops. All three methods have good application prospects, but their estimation models still need to be further improved in order to achieve an efficient application in the precise management of N fertilizer.