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.

#### 2.3.2. 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 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. Cw, Cm and Soil brightness parameter (Rsoil) 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 (SL) 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, Cw, Cm, Rsoil, θs were known inputs. Then the mean values of Ns, LID at different growth stages were the final results. Finally, we go<sup>t</sup> the parameter set of Cw, Cm, Rsoil, θs, Ns, SL and LID at each growth stage, respectively (Table 3).

#### 2.3.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.

#### 2.3.4. 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–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:

$$J = \sum\_{i=1}^{N} \frac{\sqrt{\left(V \text{Im}\_i - VIs\_i\right)^2}}{V \text{Im}\_i} \tag{1}$$

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*.


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

## 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.

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 (R2) 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).

## **3. Results**

#### *3.1. 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.


**Table 5.** Correlations between agronomic variables (*n* = 384).

\*\* Model significance at the 0.01 probability level (*p* < 0.01).

#### *3.2. 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 NDLMA 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 NDLMA, 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.


**Table 6.** Correlations between agronomic variables and vegetation indices (*n* = 384).

\* Model significance at the 0.05 probability level (*p* < 0.05); \*\* Model significance at the 0.01 probability level (*p* < 0.01); # Model with no significance.

Correlations of LNC and CND to the fifteen vegetation indices were also analyzed (Table 6). The results showed that all vegetation indices, expect NDLMA, were identified as significantly correlated with LNC and CND, respectively. In particular, CIred edge, GNDVI, MCARI/MTVI2, mND705, ND705, SPVI, NDVIcanste, and NDRE showed relatively higher correlations with LNC (*r* ≥ 0.70) than the

others, while CIred edge, GNDVI, MSR, ND705, WDRVI, SPVI, NDVIcanste, 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.

#### *3.3. 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 ge<sup>t</sup> 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.

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

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<sup>2</sup> 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.69 and 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<sup>2</sup> 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<sup>2</sup> = 0.75 and RMSE = 0.73) was superior to the LAI estimation with the 1st-par setting (R<sup>2</sup> = 0.67 and RMSE = 0.74) and the 2nd-par setting (R<sup>2</sup> = 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<sup>2</sup> = 0.80 and RMSE = 0.69) appeared stable compared with LAI estimation with the 1st-par setting (R<sup>2</sup> = 0.81 and RMSE = 0.64) and the 2nd-par setting (R<sup>2</sup> = 0.76 and RMSE = 0.93). These results indicated that the 3rd-par setting to retrieve LAI had the best estimation accuracy.


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

\* Model significance at the 0.05 probability level (*p* < 0.05). \*\* Model significance at the 0.01 probability level (*p* < 0.01). # Model with no significance.

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 <sup>g</sup>·m<sup>−</sup><sup>2</sup> 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 <sup>g</sup>·m<sup>−</sup><sup>2</sup> 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<sup>2</sup> and RMSE values of 0.45 and 20.80 <sup>μ</sup>g·cm<sup>−</sup><sup>2</sup> for the 2nd-par setting, and 0.59 and 17.43 <sup>μ</sup>g·cm<sup>−</sup><sup>2</sup> for the 3rd-par setting, respectively. However, the R<sup>2</sup> and RMSE values for the 1st-par setting were 0.30 and 58.49 <sup>μ</sup>g·cm<sup>−</sup>2, respectively. The improvements at different growth stages were also significant compared with the 1st-par setting, with an increase in R<sup>2</sup> by 0–0.13 and 0.14–0.28, and a decrease in RMSE by 21.78–51.47 <sup>μ</sup>g·cm<sup>−</sup><sup>2</sup> and 27.08–52.78 <sup>μ</sup>g·cm<sup>−</sup><sup>2</sup> 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<sup>2</sup> = 0.46 and RMSE = 21.18 <sup>μ</sup>g·cm<sup>−</sup>2), followed by the 2nd-par setting (R<sup>2</sup> = 0.39 and RMSE = 24.15 <sup>μ</sup>g·cm<sup>−</sup>2), and finally the 1st-par setting ( R<sup>2</sup> = 0.34 and RMSE = 62.86 <sup>μ</sup>g·cm<sup>−</sup>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.

#### *3.4. 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.

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

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<sup>2</sup> 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<sup>2</sup> = 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<sup>2</sup> 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 ge<sup>t</sup> CND on the basis of LAI and LND (Figure 4d–f). The R<sup>2</sup> 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 <sup>g</sup>·m<sup>−</sup>2, respectively, and their values at all growth stages were 0.66 and 2.45 <sup>g</sup>·m<sup>−</sup><sup>2</sup> (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 <sup>g</sup>·m<sup>−</sup>2, respectively, and their values at all growth stages were 0.67 and 2.08 <sup>g</sup>·m<sup>−</sup><sup>2</sup> (Table 7). CND estimations for the 3rd-par setting showed a good performance, with a higher R<sup>2</sup> values and lower RMSE values across different stages, and the R<sup>2</sup> and RMSE values across different growth stages were 0.75 and 1.32 <sup>g</sup>·m<sup>−</sup>2, respectively. Then, the validation dataset was used to validate the model stability, and the model inversion with the 3rd-par setting (R<sup>2</sup> = 0.82 and RMSE = 0.95 <sup>g</sup>·m<sup>−</sup>2) also performed better than the inversions with the 1st-par setting (R<sup>2</sup> = 0.76 and RMSE = 1.47 <sup>g</sup>·m<sup>−</sup>2) and the 2nd-par setting (R<sup>2</sup> = 0.82 and RMSE = 1.03 <sup>g</sup>·m<sup>−</sup>2). The results of this study sugges<sup>t</sup> that the 3rd-par setting in the N-PROSAIL model inversion is the best choice.

#### *3.5. 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<sup>2</sup> 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<sup>2</sup> values of 0.58, 0.58, and 0.57, respectively (Figure 3a–c). Compared with LNC estimation using the N-PROSAIL model, the R<sup>2</sup> 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.

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


**Table 8.** Relationships between LNC and vegetation indices (*n* = 384).

#: 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.

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<sup>2</sup> 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<sup>2</sup> values for CND estimation were greater than 0.75 (Table 8), and three regression models by mND705, ND705 and NDVIcanste 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 <sup>g</sup>·m<sup>−</sup>2) performed better than the vegetation indices methods (RMSE = 1.26 − 1.78 <sup>g</sup>·m<sup>−</sup>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 grea<sup>t</sup> potential for LNC and CND estimation in winter wheat.

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

## **4. 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 agreemen<sup>t</sup> 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 <sup>g</sup>·m<sup>−</sup>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<sup>2</sup> values of 0.51 and 0.82 and lower RMSE values of 0.94% and 1.03 <sup>g</sup>·m<sup>−</sup><sup>2</sup> 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 Cw 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 <sup>g</sup>·m<sup>−</sup>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.

CIred edge, mND705 and NDVIcanste 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 NDVIcanste 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 managemen<sup>t</sup> strategies.

## **5. 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<sup>2</sup> = 0.80 and RMSE = 0.69) and LND (R<sup>2</sup> = 0.46 and RMSE = 21.18 <sup>μ</sup>g·cm<sup>−</sup>2) were consistent with the measured LAI and LND, respectively.

(2) LNC and CND were accurately estimated using the N-PROSAIL model, with R<sup>2</sup> and RMSE values of 0.75 and 0.38%, and 0.82 and 0.95 <sup>g</sup>·m<sup>−</sup>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.78 <sup>g</sup>·m<sup>−</sup>2). The N-PROSAIL model shows a grea<sup>t</sup> potential to estimate canopy N status at leaf and canopy scales in winter wheat.

**Author Contributions:** G.Y., Z.L. (Zhenhai Li) and C.Z. conceived and designed the experiments; Z.L. (Zhenhai Li) and H.Y. performed the experiments; Z.L. (Zhenhai Li) and X.J. analyzed the data; Z.L. (Zhenhai Li), Z.L. (Zhenhong Li) and C.Z. discussed and drafted the manuscript. Z.L. (Zhenhong Li), J.D. and B.C. revised the manuscript and edited English language. All authors read and approved the final version.

**Funding:** This study was supported by the National Natural Science Foundation of China (Grant No. 61661136003, 41471285, 41601369), the UK Science and Technology Facilities Council through the PAFiC project (Ref: ST/N006801/1), the Open Research Fund of Key Laboratory of Digital Earth Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences (No. 2016LDE008), and the National Key Technologies of Research and Development Program (2016YFD0300602-04).

**Acknowledgments:** We are grateful to Weiguo Li, Hong Chang and Ling Kong for data collection.

**Conflicts of Interest:** The authors declare no conflict of interest.
