Spatial Variability Analysis of Within-Field Winter Wheat Nitrogen and Grain Quality Using Canopy Fluorescence Sensor Measurements

Wheat grain protein content (GPC) is a key component when evaluating wheat nutrition. It is also important to determine wheat GPC before harvest for agricultural and food process enterprises in order to optimize the wheat grading process. Wheat GPC across a field is spatially variable due to the inherent variability of soil properties and position in the landscape. The objectives of this field study were: (i) to assess the spatial and temporal variability of wheat nitrogen (N) attributes related to the grain quality of winter wheat production through canopy fluorescence sensor measurements; and (ii) to examine the influence of spatial variability of soil N and moisture across different growth stages on the wheat grain quality. A geostatistical approach was used to analyze data collected from 110 georeferenced locations. In particular, Ordinary Kriging Analysis (OKA) was used to produce maps of wheat GPC, GPC yield, and wheat canopy fluorescence parameters, including simple florescence ratio and Nitrogen Balance Indices (NBI). Soil Nitrate-Nitrogen (NO3-N) content and soil Time Domain Reflectometry (TDR) value in the study field were also interpolated through the OKA method. The fluorescence parameter maps, soil NO3-N and soil TDR maps obtained from the OKA output were compared with the wheat GPC and GPC yield maps in order to assess their relationships. The results of this study indicate that the NBI spatial variability map in the late stage of wheat growth can be used to distinguish areas that produce higher GPC.


Introduction
Wheat crops are more important in northern China than in other parts of China.The genetic background and, to a large extent, environmental factors, such as N availability, water and temperature conditions, determine the protein concentration in wheat [1][2][3][4].Therefore, real-time monitoring of wheat plant N status and pre-harvest prediction of grain yield or protein yield or both could assist producers in improving N management strategies as well as enabling yield or quality maps to be generated [5].Wheat GPC and GPC yield are important factors in the evaluation of wheat nutrition.Advanced knowledge of grain protein of a standing crop may also provide opportunities to manipulate inputs to optimize protein outputs.However, none of these pre-harvest strategies can be achieved unless there is a reliable technology to quantitatively forecast GPC of crops before they are fully ripened [6].
Remote sensing is a more feasible alternative to laboratory-based N analysis.It provides site-specific, non-destructive, economical, large-area estimation of the N status of a crop.This technique is effective in monitoring N status by estimating leaf chlorophyll (Chl) concentration which is primarily determined by N availability [7].In fact, nitrogen is a structural element of chlorophyll and protein molecules, and thereby affects formation of chloroplasts and accumulation of chlorophyll in them [8][9][10].Numerous studies have assessed the effects of N availability on canopy spectral reflectance measurements.Stone et al. demonstrated a high correlation between plant N spectral index (PNSI), the reciprocal of Normalized Difference Vegetation Index (NDVI), and total N uptake of wheat [11].The feasibility of field evaluation for crop N status using canopy reflectance spectra has also been tested [12,13].These reports showed that remote sensors were reliable in estimating plant N status.Wuest and Cassman suggest that the early season N environment has a significant influence on N partitioning at maturity [14].The ability to determine the N status of wheat and relate it to the N accumulation in the grain raises the possibility of predicting protein levels in wheat grain using remote sensing data.Hansen et al. tried to use early, repeated remote sensing multispectral data to predict grain yield and quality of winter wheat and spring barley [15].The N content of winter wheat at the anthesis stage was found to be indicative of the final protein content of the grain and the correlation coefficient between leaf N concentration at anthesis and GPC was 0.726 (n = 26) [16].
Light energy is absorbed by chlorophyll, carotenoids and other pigment molecules present in the photosynthetic antenna molecules in the thylakoid membranes of green plants.Chlorophyll fluorescence (ChlF) has been used for decades to elucidate the organization, function, and acclimation of the photosynthetic apparatus at the subcellular and leaf levels [17,18].Recent developments in sun-induced chlorophyll fluorescence suggest great potential in remote sensing to evaluate the effects of drought stress on chaparral communities [19,20].Fluorescence sensing methods have also been used to monitor crop physiology for years, and may offer advantages for N status estimation compared to reflectance-based methods.Canopy fluorescence emissions are associated with the biophysical attributes of crop growth that could potentially assist in the site-specific management of variable rate N fertilization application [21][22][23].Fluorescence spectroscopy, as an optical measurement alternative to reflectance, can be used to estimate the chlorophyll content using the ratio of red to far-red chlorophyll fluorescence bands [24,25].This ratio is reported to correlate to the N supply in winter barley and winter wheat [26], sugar beet [27], and winter oilseed rape [28].The recently developed Multiplex 3 portable sensor (FORCE-A, Orsay, Paris, France), based on the chlorophyll fluorescence mechanism, allows for the simultaneous detection of both chlorophyll and epidermal flavonol compounds [29].It has recently been applied to the detection of N deficiency in wheat [30] and turfgrass [31].
The majority of the reported studies limit their focus to crop N monitoring.There are few reports in the literature examining the application of analysis to spatial and temporal wheat N and wheat GPC, GPC yield.Moreover, spatial variability of soil properties may affect wheat N status, yield and grain quality within a crop.The magnitude and structure of this field variability suggest the suitability of site-specific management [32,33], which aims to increase both profitability and environmental protection by reducing the risk of pollution from chemical inputs applied in excessive amounts [34][35][36][37].Given the complexity of the interactions among the factors that affect grain yield and quality, the objectives of this research were: (i) to assess the spatial and temporal variability of wheat N attributes related to the grain quality of winter wheat production using canopy fluorescence; and (ii) to examine the influence of spatial variability of soil N and moisture content at different growth stages on wheat grain quality.

Field Experimental Site
This study was carried out at the National Experimental Station for Precision Agriculture, which is located in the Changping District of Beijing (40 • 10.6 N, 116 • 26.3 E) (Figure 1).This site has a mean annual rainfall of 508 mm and a mean annual temperature of 13 • C. The winter wheat (Triticum aestivum L.) was sown on 27 September 2011, with a 15 cm row spacing and 300 kg/ha sowing density.The wheat cultivar was Jingdong 22 (Beijing Academy of Agricultural and Forestry Sciences, or BAAFS, Beijing, China), which is one of the main winter wheat varieties grown in Northern China.Fertilizer was applied on 29 September 2011 with 337.5 kg/ha of DAP (diammonium phosphate) and 150 kg/ha of urea.Supplementary fertilizer was applied on 5 April 2012 with 262.5 kg/ha of urea.

Field Experimental Site
This study was carried out at the National Experimental Station for Precision Agriculture, which is located in the Changping District of Beijing (40°10.6′N,116°26.3′E)(Figure 1).This site has a mean annual rainfall of 508 mm and a mean annual temperature of 13 °C.The winter wheat (Triticum aestivum L.) was sown on 27 September 2011, with a 15 cm row spacing and 300 kg/ha sowing density.The wheat cultivar was Jingdong 22 (Beijing Academy of Agricultural and Forestry Sciences, or BAAFS, Beijing, China), which is one of the main winter wheat varieties grown in Northern China.Fertilizer was applied on 29 September 2011 with 337.5 kg/ha of DAP (diammonium phosphate) and 150 kg/ha of urea.Supplementary fertilizer was applied on 5 April 2012 with 262.5 kg/ha of urea.

Data Collection
A square area of 1.1 ha was selected as the experimental area (Figure 1).This area was divided to 110 small plots with a size of 10 m × 10 m.In each plot, four 1 m 2 area distributed in the square were selected for canopy fluorescence spectral measurements, physiological and biochemical analyses.Measurements were performed at five growth stages: the wheat tillering stage (12 April 2012), jointing stage (27 April 2012), heading stage (10 May 2012), milking stage (24 May 2012) and ripening stage (6 June 2012).Winter wheat yield data and GPC data were collected during the wheat harvest season (16 June 2012).The Feekes scale is a widely used numerical scale that describes the growth and development of cereal crops.In this study, the wheat tillering stage (Feekes 2-and 3) begins with the emergence of lateral shoots (tillers) from the axils of the true leaves at the base of the main stem of the plant.Jointing or stem elongation is the next phase of growth (Feekes 4-9).The leaves of overwintering (dormant) wheat are generally short and lie rather flat.The heading stage begins when the tip of the spike (head) can be seen emerging from the flag leaf sheath (Feekes 10.1), and emergence continues until the head is completely emerged (Feekes 10.5).Milking of grain filling stage follows anthesis and refers to the period during which the kernel matures or ripens (Feekes 10.54 to 11.4) [38].

Data Collection
A square area of 1.1 ha was selected as the experimental area (Figure 1).This area was divided to 110 small plots with a size of 10 m × 10 m.In each plot, four 1 m 2 area distributed in the square were selected for canopy fluorescence spectral measurements, physiological and biochemical analyses.Measurements were performed at five growth stages: the wheat tillering stage (12 April 2012), jointing stage (27 April 2012), heading stage (10 May 2012), milking stage (24 May 2012) and ripening stage (6 June 2012).Winter wheat yield data and GPC data were collected during the wheat harvest season (16 June 2012).The Feekes scale is a widely used numerical scale that describes the growth and development of cereal crops.In this study, the wheat tillering stage (Feekes 2-and 3) begins with the emergence of lateral shoots (tillers) from the axils of the true leaves at the base of the main stem of the plant.Jointing or stem elongation is the next phase of growth (Feekes 4-9).The leaves of overwintering (dormant) wheat are generally short and lie rather flat.The heading stage begins when the tip of the spike (head) can be seen emerging from the flag leaf sheath (Feekes 10.1), and emergence continues until the head is completely emerged (Feekes 10.5).Milking of grain filling stage follows anthesis and refers to the period during which the kernel matures or ripens (Feekes 10.54 to 11.4) [38].

The Multiplex Sensor and Nitrogen Balance Index (NBI)
Multiplex 3 (FORCE-A, Orsay, France) is a hand-held multi-parameter fluorescence sensor that operates with light-emitting diodes (LEDs) for excitation and filtered photodiodes for fluorescence detection, as described in the literature [29].This sensor detects a signal emitted by plant fluorescent pigments (fluorophores) after excitation.The four excitation channels are UV (375 nm), blue (450 nm), green (510 nm) and red (630 nm), with three detection channels being yellow, red and far-red.The Chl fluorescence in the red (RF) band at 680-690 nm and in the far-red (FRF) band at 730-780 nm was acquired sequentially at all the excitation wavelengths.A special diaphragm was mounted at the front of the sensor to illuminate a 4-cm-diameter surface (12.6 cm 2 ) at a 10-cm distance from the source and detectors.Each measurement consisted of a sequence of 250 flashes of two colors (UV and R).Different combinations of the RF and FRF fluorescence signals at the various excitation bands could be used as indices of different compounds, such as flavonoids, anthocyanins and Chl.The Chl fluorescence signals RF_R and FRF_R excited with red (R) light, FRF_UV excited with ultraviolet (UV) radiation, and FRF_G excited with green (G) light, were used to calculate the Chl indices SFR_G, SFR_R and two Nitrogen Balance Indices NBI_G and NBI_R (Table 1).In this study, Multiplex measurements were repeated 20 times for each plot and the average values were calculated to represent the fluorescence spectral indices of a particular plot.Multiplex measurements were performed over the five wheat growth stages: tillering, jointing, heading, milking and ripening.

Wheat Canopy Nitrogen Density (CND)
After wheat canopy fluorescence sensor measurements were completed, samples for determining leaf area index (LAI), specific leaf weight (SLW, g•m −2 ) and leaf N content (LNC, %) were collected on the same day.The plants in the FOV of the Multiplex were cut with a pair of scissors, placed in a plastic bag and transported to the laboratory for analysis.Each sample had all green leaves separated from the stems.LAI was determined by using a dry weight method [39].Leaf segments of about 0.04 m 2 cut from the central part of about 20 leaves were selected as reference leaves for LAI calculation.Both reference leaves and the remaining leaves were oven-dried at 70 • C to constant weight and weighed.LAI was calculated using the Equation ( 1): where S r (m 2 ) is the area of the reference leaves, W t (g) is the total dry weight of the sampled leaves, S l is the sampled land area (m 2 ), and the W r (g) is the dry weight of the reference leaves.SLW was calculated from the dried weight and area of the reference leaves.After LAI determination, leaf samples were ground and passed through a 40-mesh screen.LNC (%) was determined using the Kjeldahl method [40] with a B-339 Distillation Unit.Canopy N density (CND, g•m −2 ) was defined as total amount of N present in the canopy per unit area in this study.CND was calculated as an index of LNC, SLW and LAI using the Equation ( 2) [41].

Wheat Grain Yield and GPC
Wheat GPC samples for each plot were collected manually during the harvest season.Plants from four 0.3 m × 1 m area distributed in the square area were collected for each plot.These four samples were then mixed in laboratory and air dried.The wheat grain was threshed by a thresher and weighted for yield.Then the wheat grain samples were ground to flour and the GPC of winter wheat flour was determined with the NIR Instalab-610 (Newport Scientific Pty Ltd., Warriewood, New South Wales, Australia).All the GPC were normalized on a 14% wet weight basis.GPC yield was calculated using the equation: Soil NO 3 -N samples from 0 to 20 cm and 20 to 40 cm depth were collected in 110 locations on the 10 m × 10 m grid.The first set of soil samples was taken on 22 March 2012 at the wheat tillering stage.The second set of soil samples were taken on 19 June 2012, just prior to the wheat harvest.The location of each soil sample was recorded by a DGPS receiver with a station-based differential signal (DGPS, Trimble 5700 RTK, Trimble Inc., Sunnyvale, CA, USA).The soil samples were taken manually using a stainless steel probe.For each plot, a composite soil sample comprising five random samples from within a 5 m radius were taken, then placed into its designated plastic bags, mixed, and transported back to the laboratory.Soil NO 3 -N was analyzed using a rapid flow analyzer (Alpkem Corporation, Methodology A 303-S170, Clackamas, OR, USA), and the measured concentrations were expressed in mg•kg −1 dry soil weight.

Soil moisture sampling
Soil moisture is a critical and potentially highly variable component of the soil environment.Time Domain Reflectometry (TDR) is a proven technology for quickly and accurately determining volumetric water content (VWC) in soil.The VWC in the soil represents the fraction of the total volume of soil that is occupied by the water contained in the soil.TDR technique was initially applied to soil moisture determination by relating the apparent soil bulk permittivity (ε) to its VWC, θ (m − 3) [42,43].The underlying principal of TDR involves measuring the travel time of an electromagnetic wave along a waveguide.The speed of the wave in soil is dependent on the bulk permittivity (ε) of the soil matrix.The fact that water (ε = 80) has a much greater dielectric constant than air (ε = 1) or soil solids (ε = 3-7) is exploited to determine the VWC of the soil.The VWC measured by TDR is an average over the length of the waveguide.In this study, the spatial distribution of soil volumetric water content in the surface layer at depth of 20 cm was detected using a manual Field Scout TM TDR-300 Soil Moisture Meter (Spectrum Technologies, Inc., Aurora, IL, USA) at the 10 m × 10 m grid spacing across the five wheat growth stages in the wheat field.

Methods of Analysis
For each selected parameter at the different wheat growth stages, a statistical analysis was performed to determine the correlations of the four fluorescence spectral indices with wheat CND, GPC, and GPC yield.Experimental semi-variograms for these variables were analyzed based on the sample data and then the best fitted semi-variogram models were selected after evaluating the R 2 -values.Variowin 2.2 was used to compute the variograms and ArcGIS 9.0 (Esri, Redlands, CA, USA) was used for Kriging analysis and the generation of spatial variability maps.

Semi-Variogram Modeling
Geostatistics provide quantitative descriptions of natural variable distributions in space and time [44].Based on the regionalized variable theory, geostatistic methods assume that variables in an area exhibit both random and spatially structured properties [45].The experimental semi-variogram is a graphical representation of the mean square variability between neighboring points of distance (h) as shown in Figure 2 and Equation (4).For observations Z i , i = 1, . . . . . ., k at locations x i , . . . . . .x k , the empirical semi-variogram is defined as: where N(h) denotes the set of pairs of observation i, j such that x i − x j = h, and |N(h)| is the number of pairs in the set (generally, an "approximate distance" h is used, implemented using a certain tolerance).

Semi-Variogram Modeling
Geostatistics provide quantitative descriptions of natural variable distributions in space and time [44].Based on the regionalized variable theory, geostatistic methods assume that variables in an area exhibit both random and spatially structured properties [45].The experimental semi-variogram is a graphical representation of the mean square variability between neighboring points of distance (ℎ) as shown in Figure 2 and Equation (4).For observations Z  ,  = 1, … … ,  at locations x  , … … x  , the empirical semi-variogram is defined as: where (ℎ) denotes the set of pairs of observation ,  such that |   −   | = ℎ, and |(ℎ)| is the number of pairs in the set (generally, an "approximate distance" ℎ is used, implemented using a certain tolerance).The parameters nugget, sill and range (Figure 2) are often used to describe semi-variograms.Nugget is the height of the jump of the semi-variogram at the discontinuity at the origin.Sill is the limit of the variogram tending to infinity lag distances.Range is the distance in which the difference of the variogram from the sill becomes negligible.In models with a fixed sill, it is the distance at which this is first reached.For models with an asymptotic sill, it is conventionally taken to be the distance when the semi-variance first reaches 95% of the sill.
The empirical variogram cannot be computed at every lag distance ℎ and due to variation in the estimation, it is not ensured that it is a valid variogram, as defined above.However, some geostatistical methods such as kriging need valid semi-variograms.In applied geostatistics, the empirical variograms are thus often approximated by model function ensuring validity [47].Some important models include the spherical variogram model: the exponential variogram model: and Gaussian variogram model: The parameters nugget, sill and range (Figure 2) are often used to describe semi-variograms.Nugget is the height of the jump of the semi-variogram at the discontinuity at the origin.Sill is the limit of the variogram tending to infinity lag distances.Range is the distance in which the difference of the variogram from the sill becomes negligible.In models with a fixed sill, it is the distance at which this is first reached.For models with an asymptotic sill, it is conventionally taken to be the distance when the semi-variance first reaches 95% of the sill.
The empirical variogram cannot be computed at every lag distance h and due to variation in the estimation, it is not ensured that it is a valid variogram, as defined above.However, some geostatistical methods such as kriging need valid semi-variograms.In applied geostatistics, the empirical variograms are thus often approximated by model function ensuring validity [47].Some important models include the spherical variogram model: the exponential variogram model: and Gaussian variogram model: where h represents lag distance, a represents (practical) range, and c represents sill.

Ordinary Kriging Analysis
Ordinary Kriging Analysis (OKA) was used to generate the spatial distribution for wheat fluorescence spectral indices, wheat GPC, GPC yield data, and soil attribute data.Ordinary Kriging assumes that the mean of the process is constant and unknown within the spatial domain.A linear combination of available sample values is used for Ordinary Kriging estimation.In this study, the variables of interest at the un-sampled locations were estimated by the Ordinary Kriging method, which provides the best linear unbiased estimate of regionalized variable at an un-sampled location.

Relationship between Wheat CND, Fluorescence Spectral Indices and Wheat GPC
The formation of the protein content in grain is physically dependent on plant N accumulation and its translocation to grain during the grain filling stage.Wang et al. suggested that the leaf nitrogen concentration of winter wheat at anthesis stage was positively related to the final protein content of the grain [4].Fluorescence sensing methods have been used to monitor crop physiology for many years, and they may offer an alternative method in which to assess the N status diagnosis.
Significant correlations between wheat CND at tilling, jointing, heading, milking and ripening stages and GPC were found (Table 2).Table 2 also presents correlation coefficients between CND and the four fluorescence VIs at the five growth stages.The correlation coefficients of CND with NBI_G and NBI_R were consistently positive during the growth stages, and reached a significant level at the jointing, heading and ripening stages (Table 2).For SFR_G and SFR_R, negative correlations with CND at the tillering and jointing stages were observed, while significant positive correlations with CND were found at the heading, milking and ripening stages.The correlation coefficients between GPC, GPC yield, and the four fluorescence VIs at the five growth stages are shown in Table 3.Similar to the relationships between CND and fluorescence VIs, GPC was also negatively correlated with SFR_G and SFR_R at the tillering and jointing stages, and positively correlated at the heading, milking and ripening stages.As shown in Table 3, there existed significant positive correlations between GPC and the two NBI indices at all five stages.The relationships between the NBI and GPC yield were also significantly positive at all stages.
For the fluorescence spectral parameters, NBI_G showed the closest relationship with GPC at the ripening stage (r = 0.83 **), followed by NBI_G at the milking stage (r = 0.78 **) and the heading stage (r = 0.71 **).NBI_R was also significantly correlated with GPC at the jointing (r = 0.67 **) and tillering (r = 0.56 **) stages.Compared to the two SFR indices, the two NBI indices had stronger correlations with GPC.Both NBI_G and NBI_R were significantly correlated with GPC at all five wheat growth stages, indicating that the NBI can be used for monitoring wheat GPC status over time.

Relationship between Wheat GPC, GPC Yield, Soil Nitrogen and TDR
The correlation coefficients of soil NO 3 -N with wheat GPC and GPC yield in the two stages as well as those of soil TDR with wheat GPC and GPC yield at the five growth stages are summarized in Table 4.At the wheat tillering (22 March 2012) and harvest stages (19 June 2012), soil NO 3 -N at depths of 0-20 cm and 20-40 cm all showed significant relationships with wheat GPC and GPC yield.The strength of the relationship between GPC and soil TDR varied across different growth stages of winter wheat.Soil TDR in both the milking and ripening stages was significantly related with GPC, indicating that the soil moisture in the late growth stages of wheat affect the GPC status.

Within-Field Spatial Variability for Wheat Grain Quality Parameters
The spatial structure of wheat GPC and GPC yield was evaluated using isotropic semi-variogram models.The theoretical models were fitted using the experimental semi-variogram points to quantify spatial patterns with a least squares algorithm using the spherical, Gaussian and exponential models.Table 5 shows the statistical results of the best fit semi-variogram model parameters for wheat GPC and GPC yield.
The ratio of nugget to total semi-variance expressed as a percentage was used to classify spatial dependence: a ratio of less than 25% indicating strong spatial dependence; between 25% and 75% moderate spatial dependence; and greater than 75% weak spatial dependence [48,49].Nugget to sill ratio values for GPC and GPC yield were less than 25%, indicating strong spatial dependence.
The Gaussian semi-variogram model was found to be the best fit model for GPC and GPC yield.The results suggest that the spatial variations of GPC and GPC yield were primarily affected by structural factors.The spatial maps of GPC and GPC yield, obtained using Ordinary Kriging, revealed some distinct spatial patterns and a spatial association between them.In particular, from the maps of the two wheat grain quality parameters (GPC and GPC yield), it can be seen that GPC and CPC yield had a similar distribution in this field.The maps showed higher GPC and GPC yield (>14% and >10 T/ha, respectively) along the eastern portions of the field.While lower GPC and GPC yield values were distributed in the northern and middle parts of the field.The semi-variogram range of GPC was larger than that of GPC yield (Table 5), and the GPC distribution also looked more stable than that of GPC yield (Figure 3).The spatial maps of GPC and GPC yield, obtained using Ordinary Kriging, revealed some distinct spatial patterns and a spatial association between them.In particular, from the maps of the two wheat grain quality parameters (GPC and GPC yield), it can be seen that GPC and CPC yield had a similar distribution in this field.The maps showed higher GPC and GPC yield (>14% and >10 T/ha, respectively) along the eastern portions of the field.While lower GPC and GPC yield values were distributed in the northern and middle parts of the field.The semi-variogram range of GPC was larger than that of GPC yield (Table 5), and the GPC distribution also looked more stable than that of GPC yield (Figure 3).

Within-Field Spatial Variability of Wheat NBI Parameters
Table 6 shows the best fit model parameters for NBI_G and NBI_R at the five stages of wheat growth.The range for NBI increased slowly from 18.20 m at the tilling stage to 31.00 m at the ripping stage.Nugget to sill ratio values for NBI_R and NBI_G were less than 25%, indicating strong spatial dependence.The Gaussian and spherical semi-variogram models were found to be the best fit models for NBI_G and NBI_R.4 shows the winter wheat canopy NBI_G and NBI_R maps generated by the Ordinary

Within-Field Spatial Variability of Wheat NBI Parameters
Table 6 shows the best fit model parameters for NBI_G and NBI_R at the five stages of wheat growth.The range for NBI increased slowly from 18.20 m at the tilling stage to 31.00 m at the ripping stage.Nugget to sill ratio values for NBI_R and NBI_G were less than 25%, indicating strong spatial dependence.The Gaussian and spherical semi-variogram models were found to be the best fit models for NBI_G and NBI_R.
Figure 4 shows the winter wheat canopy NBI_G and NBI_R maps generated by the Ordinary Kriging during the five growth stages within the study field.These maps illustrate that the NBI indices were steadily distributed through all wheat growing stages.The wheat NBI_G and NBI_R spatial distributions were different from one growth stage to another.The maps of wheat NBI_G and NBI_R showed some similar spatial characteristics.In particular, higher values of NBI occurred in some areas localized along the northeastern parts of the study field (Figure 4), while lower NBI values were observed in the northwestern area of the field.Comparing Figures 3a and 4, the spatial distribution of GPC was very similar to that of NBI_R at the ripening stage (Figure 4e).While the GPC yield map (Figure 3b) was similar to the NBI_R map at the milking stage (Figure 4d).

Within-Field Spatial Variability of Soil NO3-N
Many studies report that external factors such as fertilization, irrigation and meteorology also influence the wheat GPC [50,51].In this study, soil NO3-N content at depths of 0-20 cm and 20-40 cm were sampled for each plot in the tillering and harvest stages.Table 7 shows the best fit semi-variogram model parameters for soil NO3-N content at 20-and 40-cm depths during the two growth stages.The wheat NBI_G and NBI_R spatial distributions were different from one growth stage to another.The maps of wheat NBI_G and NBI_R showed some similar spatial characteristics.In particular, higher values of NBI occurred in some areas localized along the northeastern parts of the study field (Figure 4), while lower NBI values were observed in the northwestern area of the field.Comparing Figures 3a and 4, the spatial distribution of GPC was very similar to that of NBI_R at the ripening stage (Figure 4e).While the GPC yield map (Figure 3b) was similar to the NBI_R map at the milking stage (Figure 4d).

Within-Field Spatial Variability of Soil NO 3 -N
Many studies report that external factors such as fertilization, irrigation and meteorology also influence the wheat GPC [50,51].In this study, soil NO 3 -N content at depths of 0-20 cm and 20-40 cm were sampled for each plot in the tillering and harvest stages.Table 7 shows the best fit semi-variogram model parameters for soil NO 3 -N content at 20-and 40-cm depths during the two growth stages.Nugget to sill ratio values for soil NO 3 -N content at depths of 0-20 cm and 20-40 cm during the wheat tillering stage were between 25% and 75%, indicating a moderate spatial dependence.The spherical semi-variogram models were found to be the best fit for soil NO 3 -N content at depths of 0-20 cm and 20-40 cm at the wheat tillering stage and for soil NO 3 -N content at depth of 0-20 cm at the harvest stage, while an exponential semi-variogram model was found to be the best fit for soil NO 3 -N content at the depth of 20-40 cm at the tillering stage.
The soil maps of NO 3 -N content at the depth of 0-20 cm at the tillering (Figure 5a) and harvest stages (Figure 5c) were similar, and the maps of soil NO 3 -N content at the depth of 20-40 cm at the tillering (Figure 5b) and harvest stages (Figure 5d) were also similar.Higher soil NO 3 -N content values were distributed along the eastern part of the field, while the lower soil NO 3 -N content values were observed in the middle and western parts of the field.
Remote Sens. 2017, 9, 237 11 of 18 Nugget to sill ratio values for soil NO3-N content at depths of 0-20 cm and 20-40 cm during the wheat tillering stage were between 25% and 75%, indicating a moderate spatial dependence.The spherical semi-variogram models were found to be the best fit for soil NO3-N content at depths of 0-20 cm and 20-40 cm at the wheat tillering stage and for soil NO3-N content at depth of 0-20 cm at the harvest stage, while an exponential semi-variogram model was found to be the best fit for soil NO3-N content at the depth of 20-40 cm at the tillering stage.
The soil maps of NO3-N content at the depth of 0-20 cm at the tillering (Figure 5a) and harvest stages (Figure 5c) were similar, and the maps of soil NO3-N content at the depth of 20-40 cm at the tillering (Figure 5b) and harvest stages (Figure 5d) were also similar.Higher soil NO3-N content values were distributed along the eastern part of the field, while the lower soil NO3-N content values were observed in the middle and western parts of the field.

Within-Field Spatial Variability of Soil TDR
Soil water status is critical to plant growth.The Time Domain Reflectometry (TDR) technique was used to detect the soil moisture at a depth of 0-20 cm during the five wheat growth stages.Table 8 shows the best fit semi-variogram model parameters for soil TDR at different growth stages.

Within-Field Spatial Variability of Soil TDR
Soil water status is critical to plant growth.The Time Domain Reflectometry (TDR) technique was used to detect the soil moisture at a depth of 0-20 cm during the five wheat growth stages.Table 8 shows the best fit semi-variogram model parameters for soil TDR at different growth stages.It can be seen from Table 8 that the range for TDR varied from 23.40 m to 132.0 m at the five stages of wheat growth.Nugget to sill ratio values for TDR at the five growth stages were all less than 25%, indicating strong spatial dependence.The exponential semi-variogram models were found to be the best fit for TDR at the tilling and jointing stages.The Gaussian and spherical semi-variogram models were the best fit for TDR at the heading, milking and ripening stages.The maps of soil TDR varied among the growth stages.Higher soil TDR values were distributed in the southern areas of the field at the tillering (Figure 6a) and jointing stages (Figure 6b).Then higher TDR values area expand to the eastern and northern part of the study field at the heading and milking stage (Figure 6c,d).At wheat ripening stage, the highest TDR values occurred in the western part of the field (Figure 6e).For the soil TDR Kriging map at the heading stage, higher soil TDR values were distributed along the eastern part of the field, similar to the Kriging maps of GPC and GPC yield (Figure 3).It can be seen from Table 8 that the range for TDR varied from 23.40 m to 132.0 m at the five stages of wheat growth.Nugget to sill ratio values for TDR at the five growth stages were all less than 25%, indicating strong spatial dependence.The exponential semi-variogram models were found to be the best fit for TDR at the tilling and jointing stages.The Gaussian and spherical semi-variogram models were the best fit for TDR at the heading, milking and ripening stages.The maps of soil TDR varied among the growth stages.Higher soil TDR values were distributed in the southern areas of the field at the tillering (Figure 6a) and jointing stages (Figure 6b).Then higher TDR values area expand to the eastern and northern part of the study field at the heading and milking stage (Figure 6c, 6d).At wheat ripening stage, the highest TDR values occurred in the western part of the field (Figure 6e).For the soil TDR Kriging map at the heading stage, higher soil TDR values were distributed along the eastern part of the field, similar to the Kriging maps of GPC and GPC yield (Figure 3).

Discussion
Nitrogen (N) is one of the essential nutrients for growth, yield, and quality in wheat.Numerous field studies have been conducted to assess the ability to detect nitrogen stress and then apply different amounts of nitrogen fertilizer to different parts of a field according to the crop requirements

Discussion
Nitrogen (N) is one of the essential nutrients for growth, yield, and quality in wheat.Numerous field studies have been conducted to assess the ability to detect nitrogen stress and then apply different amounts of nitrogen fertilizer to different parts of a field according to the crop requirements and the soil condition [52].Schepers and Holland related crop vigor to yield and suggested that the variation in the soil contribute to the yield variation [53].Engel et al. showed that it was possible to relate nitrogen management to protein in wheat during natural water stress conditions at the grain-filling stage [54].With this information, it was possible to maintain the protein levels close to the critical values for high-quality wheat grain.Reyns et al. evaluated wheat yield, grain moisture and GPC across a field in Belgium and observed that they were not uniform across the field [55].Stewart et al. evaluated yield and protein content in durum wheat and found areas of the field with lower available soil water capacity leading to water stress during grain-fill where there was little soil organic matter and the soil texture was coarse; these effects reduced grain yield, but increased protein levels in the grain [56].In a related study, Pettersson et al. were able to develop a relationship for barley (Hordeum vulgare L.) yield and grain protein using relectance values at the early vegetative growth stage, an index for elevated daily maximum temperatures during grain-filling and normalized electrical conductivity (ECa) of the soil.These indices relate the early-season vigor with potential weather stress on the plant during grain-filling and a surrogate for water availability from the soil [57].
Non-destructive, rapid diagnosis of plant N status is important for efficient wheat growth and the management of N nutrition.The fluorescence indices, SFR (index of chlorophyll content) and NBI (index of N sufficiency) measured with the Multiplex 3 sensor, were strongly related with crop canopy N density (CND) for most of the measurements made throughout the winter wheat growth stages.These results from this study are in agreement with recent reports which describe the potential of fluorescence sensing methods for the monitoring of crop N status of other crop species [58][59][60], indicating that fluorescence indices are also suitable to determine the N status and grain quality in winter wheat.
This study is one of the first to relate the fluorescence indices to in-field wheat N and grain quality.The differences in the spatial distributions of the wheat N status from one growth stage to another were probably due to variable soil nutrition, water content and management factors influencing crop growth.Table 9 shows the correlation analysis results relating the soil NO 3 -N and soil TDR maps with the NBI maps at the five winter wheat growth stages.It revealed that the wheat canopy NBI values were influenced by soil NO 3 -N content, particularly at wheat tillering stage.The correlation coefficients of TDR with NBI_R and NBI_G were positive at the tillering, jointing and heading stages Negative correlations between TDR and NBI_R at the milking and ripening stages were observed in Table 9.The number of samples is more than 1000.
Table 10 shows the correlation analysis results relating the soil NO 3 -N, soil TDR maps and NBI maps with wheat GPC and GPC yield maps at the five winter wheat growth stages.The correlation coefficients of NBI with GPC and GPC yield were consistently positive during the wheat growth stages.Positive correlations between tillering and harvest stage Soil NO 3 -N and GPC, GPC yield were observed at all wheat growth stages in Table 10.For TDR, negative correlations with GPC and GPC yield at the jointing, milking and ripening stages were observed, while positive correlations with GPC and GPC yield were found at the tillering and heading, stages.It revealed that the water stress at wheat grain-filling stage reduced grain yield, but interacted with soil nitrogen content to increase protein levels in the grain.These results agree with the observations of Stewart et al. [56].The most significant factors influencing crop growth in a dry environment such as the study site are soil physical properties, such as texture, bulk density and organic matter, which control the water-holding capacity.Lòpez-Bellido et al. suggested that a high amount of rainfall in the vegetative period was positively correlated to yield, due to the clay texture of vertisol that absorbs a large amount of water and retains it for a long period [61].The spatial patterns of winter wheat yield and quality parameters were possibly related to the changes in spatial variation of available soil water over different growth stages.These results agree with the observations of several other studies [62][63][64].In general, water deficit during the wheat growth period and around anthesis causes yield losses due to reduction in potential grain number per unit of land area [65].Drought stress and high temperatures during grain filling can reduce the mean kernel weight by decreasing daily rates of translocation of carbohydrate reserves from the vegetative organs of the plant to the grain [66].The N nutrition is largely considered as the primary factor influencing both the production of high yields and the quality of the grain, as it influences protein concentration [67].The observed spatial GPC and GPC yield patterns in this study illustrate the significant influence that soil nutrient and water patterns have over the five winter wheat growth stages.There is evidence that combining canopy fluorescence measurements with soil maps and agronomic assessments will provide new insights into improved wheat grain protein management practices.

Conclusions
This study demonstrated that the fluorescence spectral indices SFR and NBI measured with the Multiplex 3 sensor can be used as indicators of the N status and grain quality characteristics for winter wheat.Ordinary kriging analysis enabled spatially explicit evaluation of relationships between wheat canopy N and GPC and GPC yield over different wheat growing stages.The results of this study have indicated that the NBI Kriging maps at the ripening and heading stages were very similar to the GPC map, while the NBI Kriging map at the heading stage was similar to the GPC yield map.There is evidence to suggest that the NBI spatial variation in the late growth stage of wheat can be used to distinguish the areas that produce grain with higher protein content, which will aid in maximizing the premium price of grain produced within the field.The study also highlighted the influence of soil NO 3 -N and soil moisture conditions on wheat growth and GPC in a spatially explicit manner.

Figure 1 .
Figure 1.Study area located at the National Experimental Station for Precision Agriculture in Changping, Beijing.

Figure 1 .
Figure 1.Study area located at the National Experimental Station for Precision Agriculture in Changping, Beijing.

Figure 3 .
Figure 3. Kriging maps for wheat grain quality: (a) wheat GPC map; and (b) wheat GPC Yield Map.

Figure 3 .
Figure 3. Kriging maps for wheat grain quality: (a) wheat GPC map; and (b) wheat GPC Yield Map.

Figure 4 .
Figure 4. Kriging maps for wheat NBI over different growth stages: (a) NBI_ G maps at tillering stage; (b) NBI_G maps at jointing stage; (c) NBI_R maps at heading stage; (d) NBI_R maps at milking stage; and (e) NBI_R maps at ripening stage.

Figure 4 .
Figure 4. Kriging maps for wheat NBI over different growth stages: (a) NBI_ G maps at tillering stage; (b) NBI_G maps at jointing stage; (c) NBI_R maps at heading stage; (d) NBI_R maps at milking stage; and (e) NBI_R maps at ripening stage.

Figure 5 .
Figure 5. Kriging maps for soil NO3-N content at different wheat growth stages: (a) soil NO3-N in depth of 0-20 cm in tillering stage; (b) soil NO3-N in depth of 20-40 cm in tillering stage; (c) soil NO3-N in depth of 0-20 cm in harvest stage; and (d) soil NO3-N in depth of 20-40 cm in tillering stage.

Figure 5 .
Figure 5. Kriging maps for soil NO 3 -N content at different wheat growth stages: (a) soil NO 3 -N in depth of 0-20 cm in tillering stage; (b) soil NO 3 -N in depth of 20-40 cm in tillering stage; (c) soil NO 3 -N in depth of 0-20 cm in harvest stage; and (d) soil NO 3 -N in depth of 20-40 cm in tillering stage.

Table 1 .
Four fluorescence spectral indices of Multiplex sensor.

Table 3 .
Pearson correlation coefficients (r) between GPC, GPC Yield and the winter wheat canopy fluorescence spectral indices.
Note:1GPC is grain protein content; 2 GPC yield is grain protein content yield.* and ** indicate significance at 0.05 and 0.01 levels, respectively.The number of samples is 106.

Table 4 .
Pearson correlation coefficients (r) between GPC, GPC Yield and soil NO 3 -N and soil TDR.
Note. * and ** indicate significance at 0.05 and 0.01 levels, respectively.The number of samples is 106.

Table 5 .
Best fit semi-variogram models and parameters for wheat GPC and GPC yield.

Table 5 .
Best fit semi-variogram models and parameters for wheat GPC and GPC yield.

Table 6 .
Best fit semi-variogram models and parameters for NBI_G, NBI_R at different growth stages.

Table 6 .
Best fit semi-variogram models and parameters for NBI_G, NBI_R at different growth stages.

Table 7 .
Best fit semi-variogram models and parameters for soil NO3-N content at 20-and 40-cm depths.

Table 7 .
Best fit semi-variogram models and parameters for soil NO 3 -N content at 20-and 40-cm depths.

Table 8 .
Best fit semi-variogram models and parameters for soil TDR during the wheat growth season.

Table 8 .
Best fit semi-variogram models and parameters for soil TDR during the wheat growth season.

Table 9 .
Pearson correlation coefficients (r) of soil NO 3 -N content and soil TDR maps with NBI maps at different winter wheat growth stages.

Table 10 .
Pearson correlation coefficients (r) of soil NO 3 -N content, soil TDR and NBI maps with GPC and GPC yield maps at different winter wheat growth stages.