Remote Sensing of Leaf and Canopy Nitrogen Status in Winter Wheat ( Triticum aestivum L.) Based on N-PROSAIL Model

: Plant nitrogen (N) information has widely been estimated through empirical techniques using hyperspectral data. However, the physical model inversion approach on N spectral response has seldom developed and remains a challenge. In this study, an N-PROSAIL model based on the N-based PROSPECT model and the SAIL model canopy model was constructed and used for retrieving crop N status both at leaf and canopy scales. The results show that the third parameter (3rd-par) retrieving strategy (leaf area index (LAI) and leaf N density (LND) optimized where other parameters in the N-PROSAIL model are set at different values at each growth stage) exhibited the highest accuracy for LAI and LND estimation, which resulted in R 2 and RMSE values of 0.80 and 0.69, and 0.46 and 21.18 µ g · cm − 2 , respectively. It also showed good results with R 2 and RMSE values of 0.75 and 0.38% for leaf N concentration (LNC) and 0.82 and 0.95 g · m − 2 for canopy N density (CND), respectively. The N-PROSAIL model retrieving method performed better than the vegetation index regression model (LNC: RMSE = 0.48 − 0.64%; CND: RMSE = 1.26 − 1.78 g · m − 2 ). This study indicates the potential of using the N-PROSAIL model for crop N diagnosis on leaf and canopy scales in wheat.


Introduction
Nitrogen (N) is a critical nutrient element for maintaining photosynthesis, enhancing production, and improving grain quality in crops, but the excess use of N fertilizer also results in a series of plant and environmental problems (e.g., vigorous growth, and eutrophication) [1,2].Precision farming, which considers the crop spatial N distribution, plays an important role in solving these problems, whereby the accurate crop N estimation by remote sensing technology has the potential to precisely manage N by supplying a crop's N requirement at the right place and right time [3,4].
Moreover, the application of artificial intelligence methods in crop N estimation has been reported in many recent studies and has been proved to be a better predictor than using only sensitive spectral features or vegetation indices [3,4,[18][19][20][21][22][23][24].Among these, Hansen et al., Ecarnot et al., and Li et al. showed that partial least square regression (PLSR) could accurately predict LNC in winter wheat and winter oilseed rape [18][19][20].Miphokasap et al. indicated that estimaing PNC by stepwise multiple linear regression (SMLR) performed a higher estimation than the model by vegetation indices [4].Zhang et al. demonstrated that an artificial neural network (ANN) improved the prediction of LNC with consistently higher R 2 values, and was better than that by SMLR [21].Xu et al. proposed that the optimal combination principle (OCP) method to monitor LNC in barley would exhibit better performance than the vegetation indices methods [22,23].Yao et al. and Li et al. compared the method of SMLR, PLSR, ANN, and support vector machines regression (SVR), to monitor LNC, and showed PLSR and SVR to be preferred choices for estimating LNC, and ANN was also recommended when sufficient sample size was available [1,24].
Much research on crop N estimation has been reported, and some results were satisfactory.However, a statistical relationship between spectral information and crop N status cannot be expected everywhere and every time, even for a particular sensor [25,26].The physical model inversion approach on N spectral response has rarely been developed, except for an N-based PROSPECT model (N-PROSPECT) which was extended from a PROSPECT model by replacing the specific absorption coefficients of chlorophyll in the model with equivalent N absorption coefficients, and which could accurately simulate and retrieve leaf N density (LND) at the leaf scale [27].Thus, crop N status both at leaf and canopy scales could be retrieved through integrating the N-PROSPECT model and the SAIL canopy model [28], to be defined as the N-PROSAIL model in this study.In addition, an ill-posed problem, also called inaccurate inversion, existed and is unavoidable in model parameters inversion, which could be attributed to: (1) multiple solutions in the process of the inversion and (2) uncertainties from measurements and model assumptions [25,26].Using prior information has been demonstrated as a very efficient solution to this problem [29,30].
To develop an N estimation model based on physical model and explore the suitable strategies for estimating leaf and canopy N status, the objectives of this research were: (1) to develop the N-PROSAIL model to simulate canopy reflectance responses to leaf N density (LND); (2) to assess N status both at leaf and canopy scales, i.e., LNC and canopy N density (CND, g•m −2 soil), in winter wheat using the N-PROSAIL model; (3) to reduce the ill-posed inversion and improve the accuracy of related N variables by setting prior parameters values in the N-PROSAIL model at different growth stages; and (4) to evaluate the performance of the N-PROSAIL model method by comparing it with LNC and CND estimated by vegetation index methods.

Experimental Design
The experiments were conducted over four growing seasons, 2012-2013, 2013-2014, 2014-2015, and 2015-2016, at the Xiaotangshan National Experimental Station for Precision Agriculture, (40.17 • N, 116.43 • E) in Beijing, China (Figure 1).The test variables included various cultivars, different N fertilization rates and irrigation amounts (Table 1).Experiment (Exp.) 1 was carried out in 2012-2013 as a completely randomized design with two replications of four wheat cultivars and four N fertilizer applications rates.Exp. 2 and 3 were conducted in 2013-2015 with an orthogonal experimental design with three replications of two wheat cultivars, four N fertilizer applications rates, and three irrigation amounts.Exp. 4 was designed in 2015-2016 as a completely randomized design with three replications of two wheat cultivars and four N fertilizer applications rates.Other management procedures such as pest management, weed control, and phosphate and potassium fertilizer followed local standard practices for winter wheat production.
and 2015-2016, at the Xiaotangshan National Experimental Station for Precision Agriculture, (40.17°N, 116.43°E) in Beijing, China (Figure 1).The test variables included various cultivars, different N fertilization rates and irrigation amounts (Table 1).Experiment (Exp.) 1 was carried out in 2012-2013 as a completely randomized design with two replications of four wheat cultivars and four N fertilizer applications rates.Exp. 2 and 3 were conducted in 2013-2015 with an orthogonal experimental design with three replications of two wheat cultivars, four N fertilizer applications rates, and three irrigation amounts.Exp. 4 was designed in 2015-2016 as a completely randomized design with three replications of two wheat cultivars and four N fertilizer applications rates.Other management procedures such as pest management, weed control, and phosphate and potassium fertilizer followed local standard practices for winter wheat production.

Canopy Spectral Data
Canopy hyperspectral reflectances were obtained at four key growth stages (Zadoks growth stage: 31, 47, 65, 75 of winter wheat) [31] (Table 2).Canopy reflectances were measured by an ASD FieldSpec Handheld spectrometer (Analytical Spectral Devices Inc., USA) with a spectral range of 350-2500 nm.Under clear sky conditions between 10:00 and 14:00 Beijing time, the spectrometer was held at a height of 1.0 m above the canopy to ensure the same corresponding at different growth stages.A 40 cm by 40 cm BaSO4 calibration panel served as a black, baseline reflectance.To reduce the possible effects of sky and field conditions, spectral measurements were taken at four sites in each plot and then averaged to represent the canopy reflectance of each plot.Vegetation and panel radiance measurements were taken as the average of 20 scans at an optimized integration time, with a dark current correction for each spectrometry measurement.A panel radiance measurement was taken before and after each vegetation measurement by two scans.Bare soil refelectance was also acquired at each growth period so that it could be used as the input of the N-PROSAIL model.[31] (Table 2).Canopy reflectances were measured by an ASD FieldSpec Handheld spectrometer (Analytical Spectral Devices Inc., USA) with a spectral range of 350-2500 nm.Under clear sky conditions between 10:00 and 14:00 Beijing time, the spectrometer was held at a height of 1.0 m above the canopy to ensure the same corresponding at different growth stages.A 40 cm by 40 cm BaSO 4 calibration panel served as a black, baseline reflectance.To reduce the possible effects of sky and field conditions, spectral measurements were taken at four sites in each plot and then averaged to represent the canopy reflectance of each plot.Vegetation and panel radiance measurements were taken as the average of 20 scans at an optimized integration time, with a dark current correction for each spectrometry measurement.A panel radiance measurement was taken before and after each vegetation measurement by two scans.Bare soil refelectance was also acquired at each growth period so that it could be used as the input of the N-PROSAIL model.The aboveground vegetation at spectral measurement positions was collected immediately by randomly cutting 0.25 m 2 in each plot, and the number of tillers counted.Then, plant samples of 20 representative wheat tillers were randomly selected from the collected cut plants.All green leaves were separated from the stems.A laser leaf area meter (CI-203, CID Bio-Science Inc., WA, USA) was used to measure the leaf area, and electronic scales (±0.01 g) and drying oven were used to get the leaf dry mass, wet matter weight and dry matter weight (C m , g•m −2 ) [26].Equivalent water thickness (C w , g•cm −2 or cm) was calculated based on wet matter weight and dry matter weight [32].A Carlo-Erba NA 1500 dry combustion analyzer (Carlo Erba, Milan, Italy) [33] was used to measure leaf N concentration (LNC, %).Finally, leaf N density (LND, µg•cm −2 leaves) was calculated as the Remote Sens. 2018, 10, 1463 5 of 18 LNC multiplied by the Cm, and canopy N density (CND, g•m −2 soil) was acquired by leaf area index (LAI) multiplied by LND [27,34].

Inversion Procedure of LNC and CND Estimation
The inversion procedure of LNC and CND estimation with the N-PROSAIL model is shown in Figure 2. Specific steps were as follows: (iii) 3rd-par: LAI and LND were optimized and Cm was also set as fixed values at each growth stage since it has low variation in specific crop cultivar at one growth stage [35,36].The other six parameters in the N-PROSAIL model were set as 2nd-par.
(2) Vegetation indices calculation: simulated and measured vegetation indices (see Section 2.3.3) were computed and the optimal vegetation indices correlative with LAI, LND and Cm were chosen to construct cost function in optimization method.
(3) SCE-UA optimization: the SCE-UA method (see Section 2.3.4) was used to optimize parameters in the N-PROSAIL model.Parameters of the SCE-UA method were set and parameters in the N-PROSAIL model were updated until the number of optimization iterations was more than the maximum number of trials (maxn) allowed before optimization was terminated.
(4) Validation: the end LAI, LND and Cm were considered as optimized values and they were used to calculate LNC and CND.Three strategies were validated by comparing measured and simulated values.In addition, the inversion performance by the N-PROSAIL model was compared with LNC and CND estimation modeled by vegetation indices.

The N-PROSAIL Model
The N-PROSAIL model is a combination of the N-PROSPECT leaf model [27] and SAILH canopy model [28].At the leaf scale, PROSPECT uses a leaf structure parameter and leaf biochemical contents to simulate directional-hemispherical reflectance and transmittance of various leaves [37].
The N-PROSPECT model was developed from the PROSPECT model by replacing the specific absorption coefficients of chlorophyll in the PROSPECT model with equivalent N absorption coefficients [27,37].LNC can be estimated according to LND and Cm.At canopy scale, SAILH considers canopy structures (LAI, leaf angle distribution (θs)), soil brightness, and other angle information to generate canopy reflectance [28].LAI can be retrieved according to the SAIL model.Then, CND can be calculated by multiplying LAI and LND.During the inversion of the N-PROSAIL (1) Prior parameters initialization: parameters in the N-PROSAIL model (see Section 2.3.2) were initialized and simulated references can be simulated based on these set parameters.Three strategies of parameter (par) setting were tested and are as follows: (i) 1st-par: LAI, LND and Cm were optimized and the other six parameters in the N-PROSAIL model (Table 3) were set as fixed values in the whole growth period; (ii) 2nd-par: LAI, LND and Cm were optimized and the other six parameters in the N-PROSAIL model (Table 3) were set as fixed values at each growth stage; (iii) 3rd-par: LAI and LND were optimized and Cm was also set as fixed values at each growth stage since it has low variation in specific crop cultivar at one growth stage [35,36].The other six parameters in the N-PROSAIL model were set as 2nd-par.
(2) Vegetation indices calculation: simulated and measured vegetation indices (see Section 2.3.3) were computed and the optimal vegetation indices correlative with LAI, LND and Cm were chosen to construct cost function in optimization method.
(3) SCE-UA optimization: the SCE-UA method (see Section 2.3.4) was used to optimize parameters in the N-PROSAIL model.Parameters of the SCE-UA method were set and parameters in the N-PROSAIL model were updated until the number of optimization iterations was more than the maximum number of trials (maxn) allowed before optimization was terminated.
(4) Validation: the end LAI, LND and Cm were considered as optimized values and they were used to calculate LNC and CND.Three strategies were validated by comparing measured and simulated values.In addition, the inversion performance by the N-PROSAIL model was compared with LNC and CND estimation modeled by vegetation indices.# The ranges of Cm at the 2nd-par setting according to calibration dataset; ## the ranges of Cm at the 3rd-par setting according to calibration dataset.

The N-PROSAIL Model
The N-PROSAIL model is a combination of the N-PROSPECT leaf model [27] and SAILH canopy model [28].At the leaf scale, PROSPECT uses a leaf structure parameter and leaf biochemical contents to simulate directional-hemispherical reflectance and transmittance of various leaves [37].
The N-PROSPECT model was developed from the PROSPECT model by replacing the specific absorption coefficients of chlorophyll in the PROSPECT model with equivalent N absorption coefficients [27,37].LNC can be estimated according to LND and C m .At canopy scale, SAILH considers canopy structures (LAI, leaf angle distribution (θs)), soil brightness, and other angle information to generate canopy reflectance [28].LAI can be retrieved according to the SAIL model.Then, CND can be calculated by multiplying LAI and LND.During the inversion of the N-PROSAIL model to retrieve crop parameters, nine parameters needed to be determined (Table 3).These parameters were determined at different growth stages.LAI and LND were the main retrieval parameters and they were given intervals based on calibration dataset.C w , C m and Soil brightness parameter (R soil ) values were obtained by averaging measured values at different growth stages.The solar zenith angle (θ s ) was calculated for the time in the experiment when the hyperspectral data was measured.Leaf structure parameter (Ns), hot spot parameter (S L ) and leaf inclination distribution (LID) values were firstly optimized by the SCE-UA method using the 2012-2013, 2013-2014 and 2015-2016 wheat experiment data, where parameters of LAI, LND, C w , C m , R soil , θs were known inputs.Then the mean values of Ns, LID at different growth stages were the final results.Finally, we got the parameter set of C w , C m , R soil , θ s , Ns, S L and LID at each growth stage, respectively (Table 3).

Selection of Spectral Index
Fifteen vegetation indices correlated with agronomic parameters in previous results were calculated using the equations listed in Table 4. Using these vegetation indices has the following two purposes.(1) Selecting the best vegetation index correlation with LAI, LND and Cm to establish the cost function, an equation to evaluate the consistency between simulated and measured target, in the N-PROSAIL model (see Section 2.3.4).( 2) Establishing regression equations of LNC and CND, which were used to validate the inversion performance using the N-PROSAIL model.

SCE-UA Algorithm for LNC and CND Estimation
The SCE-UA method (the abbreviation for the Shuffled Complex Evolution method developed at the University of Arizona) proposed by Duan et al. is a general purpose global optimization program [50].The method has the advantages of search efficiency at high-parameter dimensionality, convergence speed and computational efficiency, and global searching stability [51,52], and it has been proved to be a useful and effective optimization method in past studies [53][54][55].Duan et al. [51] and Wang et al. [54] give a detailed description of the steps of the SCE-UA method.There are many parameters in the SCE-UA method, but most of them were set as default in the method.The number of complexes in a sample population (ngs) and the maximum number of trials (maxn) were determined by the actual condition and they are 2 and 1000 in this study, respectively [54].The cost function used to compare simulated with measured vegetation indices in this study was selected as follows: where J is the value of cost function and N is the determined number of vegetation indices.VIm and VIs are the measured vegetation indices and simulated vegetation indices by the N-PROSAIL model, respectively.In this study, three vegetation indices respectively correlated with LAI, LND, and Cm, were selected into the cost function.In the process of iterative inversion, the minimum of J value (minJ) was given as 5% to avoid model overfitting.Thus, the terminal condition happens when the iterative number is larger than maxn or the J value is less than minJ.Pearson Correlations (r) between vegetation indices and agronomic variables (LAI, LND, LNC, and CND) were analyzed using Microsoft Office Excel (Microsoft Corporation, Washington, DC, USA).The determination coefficient (R 2 ) and root mean square error (RMSE) were used to test the general performance of different models in this study.All calculations were made using the MATLAB (v2007, MathWorks Inc., Natick, MA, USA), and all graphs were made using the R statistical software RStudio (v1.0.44,RStudio Inc., Boston, MA, USA).

Correlations among LAI, Cm, LND, LNC, and CND
Correlation coefficients between agronomic variables were analyzed using the calibration set (Table 5).The results showed highly significant differences (p-value < 0.01) between agronomic variables, but correlation (r) values showed high differences.CND had the highest correlation with LAI and LND, with r values of 0.90 and 0.84, respectively.LND calculated from LNC showed a strong correlation (r = 0.73) with LNC, while LAI also demonstrated a high correlation with LNC (r = 0.66) although the two variables were acquired separately.Cm exhibited negative correlations with LAI, LNC, and CND and a positive correlation with LND, with r values of −0.55, −0.19, −0.30, and 0.45, respectively.

Correlations between Agronomic Variables and Vegetation Indices
Fifteen vegetation indices correlated with agronomic values were analyzed (Table 6).The results showed that all spectral indices were highly significantly related to LAI (p-value < 0.01) except ND LMA which had a correlation significant at the 0.05 level.The fourteen spectral indices except DCNI had correlations greater than 0.68, and MSR had the highest correlation with LAI, with r value of 0.80.Correlation coefficients between Cm and the spectral indices showed that fifteen vegetation indices, except DCNI, indicated highly significant differences (p-values < 0.01), but the absolute r values were only from 0.14 to 0.35, which were lower than the r values between LAI and the corresponding spectral indices.All spectral indices were highly significantly related to LND (p-value < 0.01).The maximum and minimum correlations were with MCARI/MTVI2 and ND LMA , with r values of −0.56 and −0.19, respectively.According to the correlations analysis, MSR, GI, and MCARI/MTVI2 were first used to develop the cost function in the N-PROSAIL model.Correlations of LNC and CND to the fifteen vegetation indices were also analyzed (Table 6).The results showed that all vegetation indices, expect ND LMA , were identified as significantly correlated with LNC and CND, respectively.In particular, CI red edge , GNDVI, MCARI/MTVI2, mND705, ND705, SPVI, NDVI canste , and NDRE showed relatively higher correlations with LNC (r ≥ 0.70) than the others, while CI red edge , GNDVI, MSR, ND705, WDRVI, SPVI, NDVI canste , and NDRE exhibited higher correlations with CND (r ≥ 0.79) than the others.Therefore, these vegetation indices could be used to further establish regression models with the purpose of comparing and evaluating the performance of LNC and CND estimation using the N-PROSAIL model.

LAI, LND, and Cm Estimation Using the N-PROSAIL Model Inversion
The N-PROSAIL model using the SCE-UA method was first applied to retrieve LAI, Cm, and LND in this study.Three parameter settings, 1st-par setting, 2nd-par setting, and 3rd-par setting, were tried in this process in order to get the best estimation of LNC and CND.The retrieved results of each agronomic variable with each parameter setting are shown in Figure 3 and Table 7.The results showed that a high consistency between the measured LAI and simulated LAI by the N-PROSAIL model inversion with the three parameter setting (Figure 3a-c).At different growth stages, the R 2 and RMSE values between the simulated LAI and measured LAI for the 1st-par setting, the 2nd-par setting, and the 3rd-par setting ranged from 0.45-0.69and 0.56-0.93,0.45-0.68 and 0.61-1.46,and 0.59-0.70 and 0.55-0.88,respectively (Table 7).The 3rd-par setting exhibited a relatively higher R 2 and lower RMSE than the other two parameter settings.The relationships between the simulated and measured LAIs at all growth stages were analyzed together, and the results showed that LAI estimation by the N-PROSAIL model inversion with the 3rd-par setting (R 2 = 0.75 and RMSE = 0.73) was superior to the LAI estimation with the 1st-par setting (R 2 = 0.67 and RMSE = 0.74) and the 2nd-par setting (R 2 = 0.67 and RMSE = 1.08).The independent data of 2014-2015 was used to test the estimation performance, and LAI estimation with the 3rd-par setting (R 2 = 0.80 and RMSE = 0.69) appeared stable compared with LAI estimation with the 1st-par setting (R 2 = 0.81 and RMSE = 0.64) and the 2nd-par setting (R 2 = 0.76 and RMSE = 0.93).These results indicated that the 3rd-par setting to retrieve LAI had the best estimation accuracy.Cm estimations with the 1st-par setting and the 2nd-par setting in the N-PROSAIL model were achieved (Table 7 and Figure 3d-f).The results showed no significant difference between measured Cm and estimated Cm with the 2nd-par setting across the validation set and with the 3rd-par setting across the calibration and validation set.In the 1st-par setting, the interval of Cm at different growth period was not considered, and the limited range of Cm was set the same at all growth periods.Many estimation results were ranged on both sides of the interval, and there were relatively high deviations, with RMSE values of 23.79 and 27.93 g•m −2 for the calibration set and validation set, respectively (Figure 3d).The deviation between the measured Cm and estimated Cm was decreased to 16.12 and 18.50 g•m −2 for the calibration set and the validation set in the 2nd-par setting, where the interval of Cm values were based on the statistical results of the calibration set at different growth periods.However, many estimated values ranged on both sides of the interval as well (Figure 3e).Cm estimations with the 3rd-par setting in the N-PROSAIL model showed relatively lower RMSE than that of the 1st-par and 2nd-par setting, with RMSE values of 6.10 and 8.19 g m -2 for the calibration set and validation set, respectively.The results showed that Cm inversion using the optimizing algorithm of this study had a high deviation even given the limited ranges at different growth stages, and one Cm value at each growth stage had the lowest deviation in this study.
Finally, LND estimations with the three parameters setting in the N-PROSAIL model were compared (Table 7 and Figure 3g-i).The retrieval accuracy of LND was effectively improved after considering the optimizing parameters (the 2nd-par and 3rd-par settings), with R 2 and RMSE values of 0.45 and 20.80 µg•cm −2 for the 2nd-par setting, and 0.59 and 17.43 µg•cm −2 for the 3rd-par setting, respectively.However, the R 2 and RMSE values for the 1st-par setting were 0.30 and 58.49µg•cm −2 , respectively.The improvements at different growth stages were also significant compared with the 1st-par setting, with an increase in R 2 by 0-0.13 and 0.14-0.28,and a decrease in RMSE by 21.78-51.47µg•cm −2 and 27.08-52.78µg•cm −2 for the 2nd-par setting and the 3rd-par setting, respectively.The phenomenon of overestimation at LND estimation by the N-PROSAIL model with the 1st-par setting was clear (Figure 3g).LND estimation with the optimized parameters setting (the 2nd-par setting and the 3rd-par setting) reduced the problem of LND overestimation.According to the above comparison, LND estimation with the 3rd-par setting had a slightly better performance than that with the 2nd-par setting (Figure 3h,i).The validation results showed that the best performance was produced by the 3rd-par setting (R 2 = 0.46 and RMSE = 21.18µg•cm −2 ), followed by the 2nd-par setting (R 2 = 0.39 and RMSE = 24.15µg•cm −2 ), and finally the 1st-par setting (R 2 = 0.34 and RMSE = 62.86 µg•cm −2 ).Therefore, the LND estimation confirmed the operational potential of the N-PROSAIL model inversion with optimizing parameters, and the retrieval without considering Cm inversion had the best LND estimation.

LNC and CND Estimation Based on LAI, LND, and Cm
The following estimations were to acquire LNC calculated by LND and Cm, and CND calculated by LAI and LND, respectively.The retrieved results of each agronomic variable with default parameters and optimized parameters were also compared in Figure 4 and Table 7.

Comparison of the N-PROSAIL Model Method with the Vegetation Index Method
To evaluate the performance of LNC and CND estimation by the N-PROSAIL method with parameter optimization (estimations of LNC and CND with the 3rd-par setting were using in the following comparison), the estimation results using the N-PROSAIL model method were compared with the estimation results by the vegetation index method.Fourteen spectral indices with significant relationships were used to fit the regression models (Linear model, Power model, Exponential model, and logarithmic) of LNC, and the best regression model was selected as optimal regression models (Table 8).In these regression models, ten indices with R 2 values for LNC and vegetation indices were greater than 0.50, seven were greater than 0.55, and CIred edge, mND705 and NDVIcanste had the highest R 2 values of 0.58, 0.58, and 0.57, respectively (Figure 3a-c).Compared with LNC estimation using the N-PROSAIL model, the R 2 value between the estimated LNC and measured LNC was 0.62, which was better than all the regression models by vegetation indices.The validation set was used to test the model accuracy, the estimated values and measured values were compared based on the RMSE (Table 8 and Figure 4a-d).The regression model between LNC and MCARI/MTVI2 had the lowest RMSE value of 0.48% among the fourteen regression models with the RMSE values ranging from 0.48% to 0.65% (Figure 5).The RMSE value from measured LNC and predicted LNC by the N-PROSAIL model was 0.38%, which was lower than the values from all regression models by vegetation index (Tables 7 and 8).The results also showed that the consistency between the estimated values and the measured values of the N-PROSAIL method performed better than any vegetation index methods.The results showed that the N-PROSAIL model inversion could Relative to the N-PROSAIL model inversion by the 1st-par and the 2nd-par setting, LNC estimation using the 3rd-par setting performed relatively better, with R 2 and RMSE values of 0.62% and 0.48% for the calibration set.The inversion of the N-PROSAIL model by the 3rd-par setting at different growth stages was also estimated better than by the 1st-par and the 2nd-par setting (Table 7).Many LNC estimations were overestimated, which resulted from the overestimation of LND and Cm (Figure 3d,g and Figure 4a).The 2nd-par setting considering the parameter settings at different growth stages mitigated the overestimation of LNC, but LNC estimation did not performed well (R 2 = 0.14 and RMSE = 0.87%).The independent data of 2014-2015 was used to validate the model stability, and the inversion of the N-PROSAIL model by the 3rd-par setting exhibited a superior result for LNC estimation with R 2 and RMSE values of 0.75 and 0.38%.This study further showed that calibrating parameters at different growth stages is necessary, and the retrieval without considering Cm inversion achieved a satisfactory estimation for LNC.
Together, using the N-PROSAIL model, it was able to get CND on the basis of LAI and LND (Figure 4d-f).The R 2 and RMSE values for the 1st-par setting at different growth stages ranged from 0.36 to 0.74 and 1.66 to 3.19 g•m −2 , respectively, and their values at all growth stages were 0.66 and 2.45 g•m −2 (Table 7).CND estimations for the 2nd-par setting across different growth stages ranged from 0.20 to 0.73 and 1.24 to 2.57 g•m −2 , respectively, and their values at all growth stages were 0.67 and 2.08 g•m −2 (Table 7).CND estimations for the 3rd-par setting showed a good performance, with a higher R 2 values and lower RMSE values across different stages, and the R 2 and RMSE values across different growth stages were 0.75 and 1.32 g•m −2 , respectively.Then, the validation dataset was used to validate the model stability, and the model inversion with the 3rd-par setting (R 2 = 0.82 and RMSE = 0.95 g•m −2 ) also performed better than the inversions with the 1st-par setting (R 2 = 0.76 and RMSE = 1.47 g•m −2 ) and the 2nd-par setting (R 2 = 0.82 and RMSE = 1.03 g•m −2 ).The results of this study suggest that the 3rd-par setting in the N-PROSAIL model inversion is the best choice.

Comparison of the N-PROSAIL Model Method with the Vegetation Index Method
To evaluate the performance of LNC and CND estimation by the N-PROSAIL method with parameter optimization (estimations of LNC and CND with the 3rd-par setting were using in the following comparison), the estimation results using the N-PROSAIL model method were compared with the estimation results by the vegetation index method.Fourteen spectral indices with significant relationships were used to fit the regression models (Linear model, Power model, Exponential model, and logarithmic) of LNC, and the best regression model was selected as optimal regression models (Table 8).In these regression models, ten indices with R 2 values for LNC and vegetation indices were greater than 0.50, seven were greater than 0.55, and CI red edge , mND705 and NDVI canste had the highest R 2 values of 0.58, 0.58, and 0.57, respectively (Figure 3a-c).Compared with LNC estimation using the N-PROSAIL model, the R 2 value between the estimated LNC and measured LNC was 0.62, which was better than all the regression models by vegetation indices.The validation set was used to test the model accuracy, the estimated values and measured values were compared based on the RMSE (Table 8 and Figure 4a-d).The regression model between LNC and MCARI/MTVI2 had the lowest RMSE value of 0.48% among the fourteen regression models with the RMSE values ranging from 0.48% to 0.65% (Figure 5).The RMSE value from measured LNC and predicted LNC by the N-PROSAIL model was 0.38%, which was lower than the values from all regression models by vegetation index (Tables 7 and 8).The results also showed that the consistency between the estimated values and the measured values of the N-PROSAIL method performed better than any vegetation index methods.The results showed that the N-PROSAIL model inversion could be a good method to estimate LNC.Similar comparisons were obtained for CND estimation between the N-PROSAIL model and vegetation indices methods.A highly significant regression relationship between measured CND and estimated CND was demonstrated both for the N-PROSAIL model method and the vegetation indices methods.The R 2 value between the estimated CND by the N-PROSAIL model and measured CND was 0.75.For the regression model by vegetation index, ten regression models with R 2 values for CND estimation were greater than 0.75 (Table 8), and three regression models by mND705,  Similar comparisons were obtained for CND estimation between the N-PROSAIL model and vegetation indices methods.A highly significant regression relationship between measured CND and estimated CND was demonstrated both for the N-PROSAIL model method and the vegetation indices methods.The R 2 value between the estimated CND by the N-PROSAIL model and measured CND was 0.75.For the regression model by vegetation index, ten regression models with R 2 values for CND estimation were greater than 0.75 (Table 8), and three regression models by mND705, ND705 and NDVI canste were up to 0.80 (Figure 6d-f).The results of model validation from the validation set showed that the N-PROSAIL method (RMSE = 0.95 g•m −2 ) performed better than the vegetation indices methods (RMSE = 1.26 − 1.78 g•m −2 ) for estimation of CND (Table 8 and Figure 6).An advantage of stability using the N-PROSAIL model inversion was demonstrated.Furthermore, there are underestimations at high CND values using vegetation indices estimation, but the N-PROSAIL method could mitigate the phenomenon of underestimation.Overall, our results indicated that using the N-PROSAIL model with parameters optimizing has great potential for LNC and CND estimation in winter wheat.

Discussion
The N-PROSAIL model integrating the N-PROSPECT model and the SAILH model were developed to retrieve N status at leaf scale (LNC) and canopy scale (CND) [27,28].LAI, Cm, and LND were retrieved from the N-PROSAIL model, and LNC and CND were calculated according to these relationships.LND had more variability on LNC than Cm.LNC showed a strong correlation

Discussion
The N-PROSAIL model integrating the N-PROSPECT model and the SAILH model were developed to retrieve N status at leaf scale (LNC) and canopy scale (CND) [27,28].LAI, Cm, and LND were retrieved from the N-PROSAIL model, and LNC and CND were calculated according to these relationships.LND had more variability on LNC than Cm.LNC showed a strong correlation with LND (r = 0.73) and relative low correlation with Cm (r = −0.19).This is because Cm is mainly determined by the plant cultivars and growth stage, and has little change when certain other conditions change [35].But LND composed of two compartments, high N concentration in metabolic tissues and low N concentration in plant architecture, changes dynamically with plant growth [56].Therefore, the results indicated the effect of LND on LNC estimation to be more sensitive to the effect from Cm.At the canopy scale, CND strongly reflects the variability of LAI (r = 0.90) as LND was relatively stable (r = 0.52).At different growth stages, canopy information, e.g., LAI, varied significantly, especially the variation between two growth stages [22,30].This is explained by the lower coefficient of variation of LND compared to LAI (Table 5).This conclusion is in agreement with the study of Darvishzadeh [30], who found that canopy chlorophyll content was more related to LAI with an r value of 0.94.
According to the results of LAI, Cm, and LND estimations, the relationship between measured LAI and its estimation reflects more consistence and accuracy than LND and Cm (Figure 3).The main reasons are as follows: (1) LAI is the canopy characteristics and one of the variables most affected by canopy reflectance, while Cm and LND are variables at leaf level and their variations were lower than LAI; (2) the most correlated vegetation indices between these three variables were selected to build the cost function.LAI showed best correlation with MSR (r = 0.80) and was also highly correlated with MCARI/MTVI2 (r = −0.69)and GI (r = 0.79) (Table 6).The deviation was also relatively lower than with the other two variables.The study result is in line with the findings of previous studies by Feret et al. [57], Darvishzadeh et al. [30], and Li et al. [26].
In this study, three inversion strategies of estimating LNC and CND were tried in order to improve the estimation accuracy.In the 1st-par setting, LAI, Cm, and LND were considered for retrieval, and the other parameters in the N-PROSAIL model were set default values.The results showed that many Cm values were estimated at both sides of the interval and overestimations of LND were obvious.The estimations of these three variables were improved by prior parameters initialization, which is to limit the interval of the three values and assigning different values to the other parameters at different growth stages (the 2nd-par setting).The RMSE of Cm between estimated and measured values greatly decreased from 27.93 to 18.50 g•m −2 , and the overestimation of LND was eliminated (Figure 3 and Table 7).These resulted in the improvement of LNC and CND estimation, with higher R 2 values of 0.51 and 0.82 and lower RMSE values of 0.94% and 1.03 g•m −2 for LNC and CND, respectively.The results indicated the necessity of giving different values of non-optimizing parameters in the N-PROSAIL model at different growth stages.As shown in this study, the statistical values of C w and calibrating solar zenith angle varied as plant growth progressed, and N and LID also showed difference at different growth stages (Table 3).The uncertainties of model inversion were reduced through giving the different values to these parameters at different growth stages.In the 3rd-par setting, we attempted not to retrieve the Cm value and set the mean values of Cm at each growth stage as the model input.The 3rd-par setting showed more improvement than the 2nd-par setting.The lower RMSE, 8.19 g•m −2 , for Cm between estimated and measured values demonstrated the lower certainty than the above two parameters setting, and the inversion of LAI and LND were also improved (Table 7).It has two explanations for this: (1) one parameter decreased can reduce the ill-posed problem in model parameters inversion (Figure 3d-f) [29,30]; (2) Cm is an important parameter in the crop growth model and has one initial value for specific crop cultivars and changes as growth goes on [35,36].Jacquemoud indicated that Cm was first assumed a constant for one plant in the PROSPECT model construction [37].Thus, the 3rd-par setting considered as a default value at each growth stage is reasonable.The estimated performance for every variable with the 3rd-par setting also exhibits the highest accuracy.Therefore, using the 3rd-par setting in the N-PROSAIL model is a better strategy for estimating plant N status.
CI red edge , mND705 and NDVI canste were selected as the best three vegetation indices to estimate LNC.They all showed overestimation at low LNC and these samples were mainly measured at Z.S. 75.CND estimation by mND705, ND705, and NDVI canste demonstrated the same phenomenon, with CND estimation at a high value and at Z.S. 31 showed underestimation.Estimated LNC and CND using the N-PROSAIL model showed a higher accuracy than using vegetation indices.Two advantages of using the N-PROSAIL model can explain the situation.Firstly, estimating LAI and LND by the N-PROSAIL model are interactional.The two parameters are retrieved all at once and adjusted according to changes to each other, and are acquired relatively accurately in the end.So LNC and CND estimations were taken into account the results of LAI and LND.Secondly, different parameters, except LAI and LND in N-PROSAIL, were calibrated as different values at various growth stages, which greatly reduced the deviation of model.The improvement of this step was significant according to the compared results of 1st-par setting and 3rd-par setting (Figure 4).
The results showed the potential of a priori information (setting different parameters values at various growth stages) in using N-PROSAIL model for LNC and CND retrieval in winter wheat, which is also suitable to apply in other crops (rice, maize, cotton, etc.) by giving corresponding crop parameters values.However, it is still necessary to define this information as accurately as possible.Four critical growth stages in winter wheat were selected in this study and parameters in the N-PROSAIL model at different growth stages were set respectively.Further studies should focus on considering more growth stages or a high temporal resolution [58], and in this situation, the considered fixed parameters set across different growth stage will increase.The estimation accuracy will be influenced if the growth stage was determined inaccurately.Moreover, fixed parameters were determined in the study area.Further studies should verify whether these parameters change across different regions.Future research should also focus on validating the model using multi-platforms, such as unmanned aerial vehicle and satellite platform.The study should finally point that accurate plant N estimation is an important step towards precision N management, and more studies are needed to develop N-PROSAIL model-based in-season N recommendation algorithms and management strategies.

Conclusions
In this study, the N-PROSAIL model was established to retrieve winter wheat LNC and CND at different growth stages.The results suggested that: (1) The 3rd-par setting retrieval strategies with LAI and LND optimized and other parameters in the N-PROSAIL model fixed at each growth stage exhibited the highest accuracy.The retrieved LAI (R 2 = 0.80 and RMSE = 0.69) and LND (R 2 = 0.46 and RMSE = 21.18µg•cm −2 ) were consistent with the measured LAI and LND, respectively.
(2) LNC and CND were accurately estimated using the N-PROSAIL model, with R 2 and RMSE values of 0.75 and 0.38%, and 0.82 and 0.95 g•m −2 , respectively.
(3) Compared with vegetation indices regression model, the N-PROSAIL model results reduced the problems of overestimation at low LNC and underestimation at high CND, and showed better performance than any vegetation index regression model (LNC: RMSE = 0.48−0.64%;CND: 1.26−1.78g•m −2 ).The N-PROSAIL model shows a great potential to estimate canopy N status at leaf and canopy scales in winter wheat.

Figure 1 .
Figure 1.Study site and experiment plot.Figure 1. Study site and experiment plot.

Figure 1 .
Figure 1.Study site and experiment plot.Figure 1. Study site and experiment plot.

Figure 2 .
Figure 2. Flowchart of estimating leaf N concentration (LNC) and canopy N density (CND) with the N-PROSAIL model and the Shuffled Complex Evolution method developed at the University of Arizona (SCE-UA) method.

Figure 2 .
Figure 2. Flowchart of estimating leaf N concentration (LNC) and canopy N density (CND) with the N-PROSAIL model and the Shuffled Complex Evolution method developed at the University of Arizona (SCE-UA) method.

2. 3
.5.Statistical Analysis Data collected from 2012-2013, 2013-2014, and 2015-2016 (n = 384) were mainly used for analyzing the correlation between vegetation indices and agronomic variables, calibrating the parameters of N-PROSAIL model, and developing the regression models by vegetation indices.Data collected from 2014-2015 (n = 192) were used to validate the estimating performance by N-PROSAIL model and regression models.

Figure 3 .
Figure 3.Comparison of measured and estimated values of leaf area index (LAI) (a-c), Cm (d-f), and LND (g-i) based on the N-PROSAIL model in winter wheat across the calibration set.

Figure 3 .
Figure 3.Comparison of measured and estimated values of leaf area index (LAI) (a-c), Cm (d-f), and LND (g-i) based on the N-PROSAIL model in winter wheat across the calibration set.

Figure 4 .
Figure 4. Comparison of measured and estimated values of LNC (a-c) and CND (d-f) based on the N-PROSAIL model in winter wheat across the calibration set.

Figure 4 .
Figure 4. Comparison of measured and estimated values of LNC (a-c) and CND (d-f) based on the N-PROSAIL model in winter wheat across the calibration set.

Figure 5 .
Figure 5.The relationships between measured and estimated values of LNC (a-d) and CND (e-h) in winter wheat by using data from 2014-2015 (n = 192).

Figure 5 .
Figure 5.The relationships between measured and estimated values of LNC (a-d) and CND (e-h) in winter wheat by using data from 2014-2015 (n = 192).
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 19there are underestimations at high CND values using vegetation indices estimation, but the N-PROSAIL method could mitigate the phenomenon of underestimation.Overall, our results indicated that using the N-PROSAIL model with parameters optimizing has great potential for LNC and CND estimation in winter wheat.

Figure 6 .
Figure 6.The best three regression models by vegetation index for LNC (a-c) and CND (d-f) estimation in winter wheat across the calibration set (n = 384).

Figure 6 .
Figure 6.The best three regression models by vegetation index for LNC (a-c) and CND (d-f) estimation in winter wheat across the calibration set (n = 384).

Table 1 .
Summary of cultivar, soil characteristics and treatments for the four experiments.

Table 2 .
List of acquired experiment data in the four wheat experiments.

Table 3 .
Specific ranges and setting for input parameters in the N-PROSAIL model.

Table 4 .
Summary of vegetation indices studied for the N-PROSAIL inversion and N estimation.

Table 7 .
Comparison of different agronomic parameters estimation with different parameter setting.

Table 7 .
Comparison of different agronomic parameters estimation with different parameter setting.
#: Linear and nonlinear regression (power regression, exponential regression and logarithmic regression) were conducted and listed are the optimal regression model of each vegetation indices.x: vegetation index; y: LNC or CND.