# Study the Spatial-Temporal Variation of Wheat Growth Under Different Site-Specific Nitrogen Fertilization Approaches

Bei Cui, Wenjiang Huang, Xiaoyu Song, Huichun Ye, Yingying Dong

► To cite this version:

Bei Cui, Wenjiang Huang, Xiaoyu Song, Huichun Ye, Yingying Dong. Study the Spatial-Temporal Variation of Wheat Growth Under Different Site-Specific Nitrogen Fertilization Approaches. 11th International Conference on Computer and Computing Technologies in Agriculture (CCTA), Aug 2017, Jilin, China. pp.316-332, 10.1007/978-3-030-06137-1_29. hal-02124224

HAL Id: hal-02124224

https://inria.hal.science/hal-02124224

Submitted on 9 May 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Study the Spatial-temporal Variation of Wheat

Growth under Different Site-specific Nitrogen

Fertilization Approaches

Bei Cui¹,², Wenjiang Huang¹,²(✉), Xiaoyu Song³, Huichun Ye¹,², and Yingying Dong¹,²

¹ Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences,
Beijing 100094, PR China

² Key Laboratory of Earth Observation, Hainan Province, Sanya 572029, PR China
{cuibei, huangwj, yehc, dongyy} @ radi.ac.cn

³ Beijing Research Center for Information Technology in Agriculture, Beijing Academy of Agriculture and Forestry Sciences, Beijing 100097, PR China

songxy@nercita.org.cn

**Abstract.** Many variable fertilization approaches based on 'real-time' crop N status were developed for making N fertilizer management in precision agriculture. Unfortunately, to date, only few papers reported the effect of variable fertilization algorithms on the spatial and temporal variability of crop parameters. Based on these problems, this study designed three different variable fertilization algorithms based on vegetation index (Y), SPAD (S) and crop growth model (Z), respectively, with uniform fertilization and no fertilization as controls. Results showed that wheat growth had strong spatial dependence, which become stronger after fertilization. Wheat yield also had strong spatial dependence. There were some similar spatial distribution between NDVIs, soil TN and yield, indicating that spatial variability of yield had strong relationship with crop growth status and soil TN content. The site-specific fertilization treatment based on crop growth model (Z) had the best adjustment capacity to promote crop growth and yield, and reduce their spatial variation, compared with other fertilization treatments.

**Keywords:** site-specific N fertilization · winter wheat · spatial-temporal variation

# 1. Introduction

Nitrogen fertilizer is the necessary nutrient factor to ensure and maintain crop healthy growth and high production. Farmers prefer to apply excessive nitrogen fertilizers in order to gain high economic return, which may lead to the decrease of crop growth and grain production because of the osmotic limitations of soil [1, 2, 3]. In addition, in the conventional farmer management practice, the same
amount of N fertilizer is usually applied at different locations across the field, which ignore the spatial variability of crop growth status and soil parameters within a field. Therefore, many negative impacts were generated, such as N leaching, soil erosion, groundwater pollution, and even poor economic benefits [4]. Scientists found that these negative impacts could be prevented by more efficient use of nitrogen fertilizer in farming [5, 6]. In recent years, more researchers are paying close attention to site-specific fertilizer nitrogen management of precision agriculture. Many encouraging results were also reported that variable rate fertilization can increase crop yield and/or reduce fertilizer inputs. Griepentrog and Kyhn [7] pointed out that increasing yield in the high yields area by using variable rate fertilization is not possible but has 36% reduction in nitrogen inputs. Delin et al. [8] found that crop yield can be improved and N fertilizer also can be saved through a 3-year field investigation of site specific fertilization treatments. Cui et al. [9] reported that the soil Nmin-based on N management strategy can save more than half-N fertilizer without grain yield losses, compared with farmer's N practices. Li et al. [10] found that, compared with farmer practices, the N inputs of the sensor-based N management strategy and soil Nmin-based N management strategies decreased by 81.99% and 76.34%, while N use efficiencies improved by 367.94 % and 289.31 %, respectively. Colaco and Molin [11] studied a site-specific fertilization strip treatment experiment over 6 years and found that significant N fertilizer input reduction up to 39.6% without no increase in yield in field. Diacono et al. [12] summarized previous studies and found that sensor-based N management can highly increase N use efficiency of up to 368% and save N fertilizer from 10% to about 80% compared with common farmer practices.

Although crop historical yield maps can provide useful information for making site-specific
fertilization decision, which have been reported in previous researches, some uncontrollable factors,
such as rainfall, temperature and so on, may lead to crop yield variation between different growing
seasons. Therefore, the measurement and understanding ‘real-time’ crop N status are needed before
making N fertilizer recommendations. Given the time-consuming and high expensive of laboratory
test method, most researchers in recent years used various diagnostic tools, such as spectrometer [13,
14, 15, 16], hand-held chlorophyll meters [17, 18], Dualex® (FORCE-A, Orsay, France) [19], to help
them to obtain “real-time” crop growth information quickly, especially nitrogen status information.
Most of the tools for measurement crop N status mentioned above have been used to make
site-specific N management decision. Ehlert et al. [20] simulated late fertilizer rates based on existing
plant mass in cereals, indirectly measured by a mechanical sensor (pendulum meter), which is
mounted at the front of a tractor. Godwin et al. [21] proposed two N application strategies based on the
historic yield and density of shoots, which was determined by near real-time NDVI data from airborne
digital photography. Li et al. [10] developed an active canopy optical sensor-based N application
strategy for winter wheat, which was based on the predictability of NDVI to estimated yield. Morris et
al. [22] proposed a method to calculate topdressing N application rates, which is determined utilizing
wheat canopy NDVI measured by a GreenSeeker Hand Held Optical Sensor and an algorithm
developed at Oklahoma State University. Welsh et al. [23] proposed two N fertilizer application
approaches with two different estimates of yield potential: the first one is produced from historic yield
data and the second is shoot density maps derived from airborne digital photographic images.
Thomason et al. [24] calculated the variable N fertilizer rate based on the Virginia Wheat Algorithm
(VWA), of which wheat canopy reflectance was obtained using a GreenSeeker 03RT 200 system.
Unfortunately, to date, only few papers reported about the potential differences of these variable fertilization approaches.

Most previous studies of site-specific fertilization management in precision agriculture focused more efforts on assessing the potential of variable N fertilizer application decision by comparing some key agriculture or soil properties (crop yield, nitrogen use efficiency, nitrogen inputs and etc.) to conventional fertilization decision using traditional statistics method, without analyzing its effect to the spatial distribution variability of crop and soil parameters. Crop growth and yield have spatial variation within a field due to the inherent spatial variability of soil parameters [25]. Reducing the spatial variability of crop growth is an effective approach to achieve the goal of high and stable grain production in precision agriculture. Therefore, understanding the effectiveness of site-specific N fertilizer management to reduce spatial variability of crop growth and yield is very important to make and select the optimal N application approach. There should be an important evaluation criterion for variable rate fertilization decision, but few researches have tried to analyze and explore its potential on the crop spatial distribution variation.

In this study, we designed three different site-specific fertilization approaches based on different N measurement tools, with no fertilization and uniform fertilization as controls, using traditional statistics and geo-statistics methods to analyze the effectiveness of different site-specific N fertilizer application decisions. The objectives of this research were: (i) to compare and assess the potential of different variable N rate fertilizer managements to improve wheat growth and yield; and (ii) to understand the spatial- temporal variability rule of wheat growth under site specific N fertilization treatment using geo-statists method and to assess the potential of variable fertilization to reduce spatial variation of crop growth and yield.

# 2. Materials and Methods

## 2.1 Experimental site description

The experimental site used in this study is located at the National Experimental Station for Precision Agriculture, which is situated in the Changping district of Beijing (40°10.6′N, 116°26.3′E), during the 2005-2006 growing seasons. The site is warm temperate zone with semi-humid continental monsoon climate with annual average rainfall of 508 mm and annual average temperature of 13°C. The main type of soil on the site is alluvial soil. The cropping systems of this site is winter wheat-corn rotation annually.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Study area located at the National Experimental Station for Precision Agriculture in Changping, Beijing</figcaption></figure>

## 2.2 Experiment design and fertilization algorithms

The crop planted was winter wheat cultivar 'Jingdong8', which is one of the main winter wheat cultivars planted in the Northern China. The sowing was carried out on September 26, 2005 at a depth of 5 cm, with a 15 cm distance between the rows and 330 kg/ha sowing density. The experimental area has no basal-dressing treatment and urea, as the N fertilizer, was only applied on April 14, 2006. The experiments were conducted in a similar manner except the N fertilizer application rates. The crop was harvested at June 20, 2006.

The experiment design in this study used five different N fertilizer application algorithms. The experimental region was divided into five management zones and each zone was further divided into 20 small experimental strips by 3 m × 3 m. The specified N fertilizer application strategy was applied in the corresponding zone. To avoid the interaction between the strips with different N fertilizer rate, there are isolation strip wider than 1 m between adjacent strips. The spatial distribution of experimental strips was shown in Fig. 2. Compared to conventional treatments (uniform fertilization and no fertilization), three different site-specific fertilization methods have been proposed to use the in-season application of N rates, which were based on wheat canopy OSAVI value (Y), the normalized SPAD value (S) and the spectral index and crop growth model (Z), respectively.

<figure><img src="image_3.png" /><figcaption>Fig. 2. The distribution of plots</figcaption></figure>

### 2.2.1 Variable fertilization treatments

All the optimal regression equations used in the following three variable fertilization approaches were determined by the historical winter wheat experiment data during 2003-2004, 2004-2005 growing
season in the same field. In order to ensure the comparability of winter wheat yield data under different fertilization treatments, the N fertilization rate of each strip calculated by different fertilization algorithms has been slightly adjusted so that each fertilization treatment zone were applied the same amount of N fertilizer.

Y treatment: A canopy sensor-based N prediction algorithm for winter wheat was developed. The first step was to calculate the in-season estimated yield (INSEY), where the sum of OSAVIs at standing and jointing stages, measured by a portable field spectro-radiometer (FieldSpec-FR2500, ASD, USA), were divided by GDD value (the number of days that the daily average temperature stability above zero from standing stage to jointing stage). The second step was to predict the potential goal yield (PGY) using the INSEY yield prediction formula. Next, the amount of N absorbed by wheat at jointing stage (AN) and the total amount of N required during the whole growing season of winter wheat (TN) were calculated using OSAVI for jointing stage and PGY prediction formulas, respectively. Finally, the N fertilizer recommendation is accomplished by subtracting AN from TN, which was the topdressing N rate at wheat jointing stage in this study. Total 20 replicated treatment strips were designed, namely Y-1, Y-2, Y-3......Y-20 from right to left of experimental site as shown in Fig.2.

S treatment: The N fertilization rate applied in each plot was simulated based on the normalized SPAD value (NDSPAD) of the top two leaves of wheat at jointing stage, which was obtained by a chlorophyll meter (SPAD-520, MINOLIA, Japan). The first step was to calculate the yield potential for winter wheat (PGY) using the NDSPAD-yield prediction formula and the amount of N absorbed by wheat (AN) at jointing stage using the SPAD-N prediction formula. The second step was to predict the total amount of N required by wheat during the whole growth stages (TN) using the PGY-TN prediction formula. The last step was to calculate the topdressing N rate by subtracting AN from TN. Total 20 replicated treatment strips were designed, namely S-1, S-2, S-3......S-20 from right to left of experimental site as shown in Fig.2.

Z treatment: The N fertilization rate applied in each plot was simulated based on CERES-Wheat model and wheat canopy OSAVI value obtained by a portable field spectro-radiometer (FieldSpec-FR2500, ASD, USA) at jointing stage. Unlike the sensor-based variable fertilization algorithms (Y treatment) described above, the potential yield (PGY) for wheat was predicted based on CERES-Wheat model rather than canopy OSAVI. The total amount of N requirement during the whole growth stages of wheat (TN) was calculated in terms of the PGY-TN prediction formula, and the nitrogen content that has been absorbed by wheat (AN) at jointing stage was calculated according to the relationship between OSAVI and PFNU. The topdressing N rate was finally calculated by subtracting AN from TN. Total 20 replicated treatment strips were designed, namely Z-1, Z-2, Z-3......Z-20 from right to left of experimental site as shown in Fig.2.

### 2.2.2 Uniform fertilization treatment

The total amount of N fertilizer applied in uniform fertilization zone were equal to that in the other three site-specific fertilization zones. Each small strip was applied the same amount of N fertilizer, which was calculated by dividing the total amount by 20. Total 20 replicated experimental strips were designed, namely CK-1, CK-2, CK-3......CK-20 from right to left of experimental site as shown in Fig.2.
### 2.3.3 No fertilization treatment

As a control, no fertilization treatment was conducted at the experimental site. There was no N fertilizer applied in this treatment zone. Total 20 replicated treatment strips were designed, namely W-1, W-2, W-3......W-20 from right to left of experimental site as shown in Fig.2.

## 2.3 Data collection

### 2.3.1 Winter wheat canopy hyperspectral data

A portable field spectro-radiometer (FieldSpec-FR2500, ASD, USA) was used to measure wheat canopy hyperspectral data on April 13, 21, 26 and 29, May 10 and 29, June 8, 2006. The spectral range of spectro-radiometer is from 350 to 2500 nm and the spectral resolution is 1.4 nm for the region 350 to 1050 nm and 2 nm for the region 1050 to 2500 nm. To ensure the accuracy of measurement, the wheat canopy radiances was collected under clear, blue-sky conditions between 10:00 h and 14:00 h (Beijing Local Time) and the probe at a distance of around 1.3 m above the wheat canopy surface and a field of view of 25° was kept in order to maintain the same viewing geometry. The measured radiance was converted into absolute reflectance through calibration based on a white reference panel by formula (1). Total 20 individual measurements were carried out on each strip and were recorded by average as the accurate spectral value of this strip.

$$R_{\text{target}} = \frac{DN_{\text{target}}}{DN_{\text{reference}}} \times R_{\text{reference}} \quad (1)$$

Where, $R_{target}$ is canopy spectral reflectance of the targeted winter wheat, $DN_{target}$ and $DN_{reference}$ are the radiance of winter wheat canopy and the white reference panel collected by ASD spectro-radiometer, $R_{reference}$ is the white reference panel spectral reflectance provided by its manufacturing company.

### 2.3.2 Vegetation index

In the previous studies, it was frequently reported that crop canopy NDVI index showed the closest relationship with plant leaf area index, biomass, and it has become an important indicator to reflect crop growth status[26, 27, 28]. Several studies have also reported that NDVI is able to discriminate crop response to N fertilizer rates [29]. So, NDVI index was chosen to indirectly evaluate wheat growth situation under different fertilization applications in this study. The NDVI [30] was calculated with the following equation:

$$NDVI = (R_{800} - R_{680}) / (R_{800} - R_{680}) \quad (2)$$

Where, NDVI is the normalized differential vegetation index calculated with winter wheat canopy reflectance data, $R_{680}$ and $R_{800}$ are the winter wheat canopy reflectance at the wavelength of 680 nm and 800 nm, respectively.

### 2.3.3 Winter wheat yield data

Wheat samples were collected manually on June 19, 2006 from a 2 m² area for each strip and a total of 100 plant samples were obtained in the study field. After air-dry, the sampling plants thresh, and weigh, the grain yield was recorded with an electronic scale (±0.01g).

### 2.3.4 Soil total nitrogen (TN) samples

Soil total nitrogen (TN) samples from 0 to 30cm and 30 to 60 cm depth were collected in total 100
experimental plots on June 14, 2006, just prior to the wheat harvest. For each plot, five representative
soil-sampling sites were randomly selected within a 5 m radius near the center of plot, avoiding
transition areas with other plots. The soil samples from five sampling site at depths of 0-30 cm and
30-60 cm in each plot was taken by soil auger and then comprised and mixed to be as the sample of
this plot. The soil TN content at each plot was calculated using semimicro-kjeldahl method in the
laboratory.

**2.4 Analyses method and software**

The conventional statistics for raw experimental data was conducted using SPSS 22 software for
Windows (SPSSInc., Chicago, IL, USA) to get mean, maximum, minimum, standard deviation (SD),
kurtosis, skewness and coefficient of variation (CV). MATLAB R2015a was used to get bar charts for
wheat canopy NDVI and grain yield with the plotting function.

The semi-variogram model was used to indicate spatial structure and the degree of spatial
dependence of wheat canopy NDVI, grain yield and soil parameters. The semi-variogram model is
described usually by five mainly parameters, i.e. nugget, partial sill, sill, range and Nugget/Sill.
Nugget means the spatial variation of variable led by random factors and partial sill means the spatial
variation led by structure factors such as soil texture, crop variety and etc. Sill is the sum of nugget and
partial sill, which means the total spatial variation of variable. The range means the max distance
between two locations with spatial dependence. The last but the most important parameter is the value
of Nugget/Sill ratio, which can directly indicate the degree of spatial dependence and the reason led to
spatial variation by random or structure factors. The best fit theoretical semi-variogram model for
different variables were simulated and the fitting accuracy were evaluated according to R² and
residual sums of squares (RSS) indicators. All above analysis were conducted using GS 7+ software.

Ordinary Kriging interpolation method was used to plot the spatial distribution maps of wheat canopy NDVI, grain yield and soil parameters. Ordinary kriging is a widely used geostatistical interpolation method, which assumes the mean of regional variable is an unknown constant. Ordinary Kriging interpolation was implemented by using the geostatistical Analyst Tools of ArcGIS (ArcGIS 10.2, ESRI, Redland, CA, USA).

**3. Results**

**3.1 Analysis wheat parameters under different fertilization treatments**

**3.1.1 Analysis wheat growth under different fertilization treatments**

The basic statistical information of winter wheat canopy normalized differential vegetation index (NDVI) value under different fertilization treatments at six key growth stages, were shown in Fig. 3. The average NDVI values of wheat canopy presented a consistent tendency among different fertilization treatments, which was firstly increased to the maximum value on May 10, 2006, then decreased gradually. In addition, the average values of NDVI varied with different fertilization treatments at each growing stage, for example, the sequence was S, Y, Z and CK treatment from big value to small on April 13, 2006; however, after fertilization, the order of the average value of NDVI was changed among the five fertilization treatments at each observation stages. The average NDVI values for fertilization treatments, including CK, Y, S and Z, were significantly higher than the values for no fertilization treatment (W), except NDVI for Y was a little lower than W on April 21, 2006. The
Z treatment was the highest value of NDVI, followed by CK treatment at each stage after fertilization. The above results indicate that N fertilizer application management can effectively promote crop growth and different nitrogen fertilizer management's methods made different capacities of adjustment crop growth. The fertilization method based on crop growth model and canopy spectral index (Z) has the best adjustment capacity to promote crop growth.

<figure><img src="image_4.png" /><figcaption>**Fig.3.** The average value and CV of winter wheat canopy NDVIs at five different fertilization treatments from jointing stage to ripening stage</figcaption></figure>

The changes of variation coefficients of NDVI values for different fertilization treatments at different stages are given in Fig. 3 (b). A higher variation coefficient of NDVI value for no N fertilizer application treatment (W) is observed from April 21 to June 8, compared to other fertilization treatments at corresponding observation stages, which indicates that fertilization can effectively reduce the spatial variability of crop growth. For the four fertilization treatments, the variation degree of NDVI for Y and S treatments were lower than that for CK treatments after fertilization. Furthermore, firstly the NDVI for Z treatments had a higher variation coefficient than the other fertilization treatments, and then decrease to the smallest variation coefficient among the four fertilization treatments. The variation coefficients of NDVI after fertilization were lower than that at the observation stage of April 13, except the period of 8 June.

### 3.1.2 Analysis of wheat yield under different fertilization treatments

The wheat yield and its coefficient of variation for different fertilization treatments were shown in Fig. 4. From Fig. 4 (a), it is observed that the total yield of wheat in the N fertilization zones, including Z, CK, Y and S zone were apparently higher than the yield with no fertilization (W). Among the three variable N fertilizer application zones, the zone applied fertilizer based on crop growth (Z) had the highest yield, followed by S, Y zone with an increase in yield up to 18.44%, 12.08% and 11.74%, respectively, from CK zone.
<figure><img src="image_5.png" /><figcaption>Fig.4. The average value and coefficient of variation of winter wheat yield</figcaption></figure>

The variation coefficients of yield in different fertilization treatment zones were given in Fig. 4 (b), which varied between 10.78% and 27.29%, indicating that the yield of all the treatments had a moderate variation degree. The yield in W zone had the biggest variation degree (27.27%), and the S zone had the smallest variation degree (10.78%). The coefficient of variation of yield for Z zone was second-smallest among these treatments. The degree of variation of yield for the CK zone was similar to the degree for the Y zone. The coefficient of variation of yield for S, Z and Y zone was decreased by 37.61%, 16.04% and 1.93% than CK zone, respectively. In short, Z fertilization approach can effectively improve the wheat yield and S approach had the biggest help to reduce the degree of variation.

## 3.2 Analysis spatial variability of wheat canopy NDVI

To understand the changes of spatial structure and distribution of wheat growth status between before and after fertilization is very important to help to realize the effect of different N fertilizer application strategies on spatial variation of crop growth. In this study, geo-statistics method was used to analyze the spatial structure of NDVI in N fertilizer application zone, including CK, Y, S and Z zone, and ordinary kriging interpolation method was used to map the spatial distribution of wheat growth status. In order to ensure interpolation accuracy, all samples data of fertilization region were implemented by Geo-statistics analysis.

The precondition for geo-statistics analysis is that all data studied should be normally distributed. Therefore, all outliers should be firstly checked before the test of normal distribution using the Estimating Neighborhood Method [31]. If the data are still non-normality distribution after correction outliers, some efforts of data transform (log, square-root and etc.) should to conducted to ensure that the data to be normal distribution. Through the above-mentioned processing, the values of skewness and kurtosis of all data used in this study were very small as shown in Table 1, indicating that the NDVIs were normal distribution.

<table><caption>**Table 1** The coefficients of skewness and kurtosis for winter wheat NDVI</caption><thead><tr><th>Date</th><th>skewness</th><th>kurtosis</th><th>Date</th><th>skewness</th><th>kurtosis</th></tr></thead><tbody><tr><td>April 13</td><td>0.273</td><td>0.667</td><td>May 10</td><td>-0.477</td><td>-0.075</td></tr><tr><td>April 21</td><td>0.043</td><td>-0.335</td><td>May 29</td><td>-0.342</td><td>0.074</td></tr><tr><td>April 29</td><td>0.160</td><td>0.611</td><td>Jun 08</td><td>0.012</td><td>-0.823</td></tr></tbody></table>
The theoretical model parameters of best fit for wheat canopy NDVI at different growth stages were given in Table 2. The spherical model was found to be the best fit model for NDVI on April 21, April 26, May 10 and June 08 and the exponential model was found to be the best fit model for NDVI on April 13, April 29 and May 29. The range for NDVI varied from 14.10 m to 29.70 m from May 10 to May 29. Nugget/sill ratio values for NDVI at the six growth stages were all less than 25%, indicating that there were strong spatial dependence and the spatial variability was mainly resulted from structural factors such as wheat cultivar and soil characters, etc. In addition, it was also observed that there was obvious decrease of the nugget/sill ratio values during the period of April 21 to June 8 when compared with that at the period of April 13, indicating that the spatial dependence of wheat growth was enhanced further after fertilization.

<table><caption>**Table 2** Best fit theoretical semi-variogram model and parameters for NDVI</caption><thead><tr><th></th><th>Model</th><th>Range(m)</th><th>Nugget</th><th>Sill</th><th>Nugget/Sill</th><th>R<sup>2</sup></th><th>RSS</th></tr></thead><tbody><tr><td>NDVI (13/04)</td><td>Exponential</td><td>20.40</td><td>0.00033</td><td>0.00325</td><td>10.18%</td><td>0.81</td><td>1.423E-07</td></tr><tr><td>NDVI (21/04)</td><td>Spherical</td><td>16.90</td><td>0.00005</td><td>0.00327</td><td>1.53%</td><td>0.87</td><td>2.168E-07</td></tr><tr><td>NDVI (29/04)</td><td>Exponential</td><td>17.30</td><td>0.00006</td><td>0.00342</td><td>1.75%</td><td>0.94</td><td>5.296E-08</td></tr><tr><td>NDVI (10/05)</td><td>Spherical</td><td>14.10</td><td>0.00005</td><td>0.00209</td><td>2.15%</td><td>0.74</td><td>2.606E-08</td></tr><tr><td>NDVI (29/05)</td><td>Exponential</td><td>27.00</td><td>0.00011</td><td>0.00154</td><td>6.87%</td><td>0.98</td><td>4.698E-09</td></tr><tr><td>NDVI (08/06)</td><td>Spherical</td><td>17.00</td><td>0.00001</td><td>0.00580</td><td>0.17%</td><td>0.86</td><td>9.704E-07</td></tr></tbody></table>

Fig.5. visually illustrated the spatial distribution of wheat canopy NDVI on N fertilizer application zones during wheat growing season using the ordinary Kriging method. Z, CK, Y and S zone, in turn, were distributed from up to bottom on the each NDVI maps. To further understand the spatial variation, the NDVI values were divided into ten levels by geometric intervals method. Although the wheat canopy NDVI spatial distributions were different from one growth stage to another, there were still some similar spatial distribution characteristics between them. The higher NDVI value were mainly located along the southern boundary of the field, and lower NDVI were mainly located in the northern regions on April 13, 2006. However, compared with the spatial distribution of NDVI on April 13, 2006, it showed a contrary distribution tendency at the period of April 21 to June 8, which was the higher NDVI value were mainly distributed in the northern regions of field and the lower NDVI value were mainly distributed in southern regions. The phenomenon mentioned above indicated that wheat growth status had obvious change within the field and different N fertilizer application strategies were the main reason for these. Furthermore, higher NDVI value and lower NDVI value appeared some strip distribution characteristics in north-south direction during the period of wheat growing stages (Fig.5).
<figure><img src="image_6.png" /><figcaption>Fig. 5. The spatial distribution of wheat canopy NDVI at critical periods on N fertilizer application zone</figcaption></figure>

3.3 Analysis spatial variability of wheat yield

Through the outlier diagnosis and normal distribution test, the values of skewness and kurtosis for
yield were both very small, i.e. 0.111 and -0.055, respectively, indicating that the wheat yield data in
this study were normally distributed. Table 3 showed the best fit theoretical semi-variogram model
and its parameters for wheat yield. The Gaussian model was found to be the best fit model for wheat
yield. The range for yield was 12.82 m and the nugget/sill ratio value was 15.18 %, indicating strong spatial dependence and weak spatial variability.

<table><caption>**Table 3** Best fit theoretical semi-variogram model and parameters for yield</caption><thead><tr><th>Model</th><th>Range(m)</th><th>Nugget</th><th>Sill</th><th>Nugget/Sill</th><th>R<sup>2</sup></th><th>RSS</th></tr></thead><tbody><tr><td>Gaussian</td><td>12.82</td><td>90000</td><td>593000</td><td>15.18%</td><td>0.81</td><td>8.340E+09</td></tr></tbody></table>

The spatial distribution pattern of wheat yield on N fertilizer application zones was shown in Fig. 6. Z, CK, Y and S zone, in turn, were distributed from up to bottom on the yield maps. Higher value areas of wheat yield mainly occurred on the northern (Z zone), and secondary southern area (S zone) of the field, while lower yield areas are mainly distributed on the middle region of the field (Y zone and CK zone). Compared with the NDVI maps in Figure 5, the spatial distribution pattern of yield had some similar spatial distribution characteristics with them. Lower yield pattern were consist with the spatial distribution pattern of lower NDVI on the NDVI maps at different growing stages described above except June 8, 2006. Some higher yield area also has higher NDVI. Based on the above evident, a conclusion can be drawn that the areas with stable high or low crop growth status after fertilization usually have high or low yield on different N fertilizer application zones.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** The spatial distribution of wheat yield on N fertilizer application zones (kg/ha)</figcaption></figure>

### 3.4 Within-field Spatial variability of Soil Total Nitrogen

The inherent spatial variability of soil parameters is one of main reasons of the spatial variation of crop growth and yield. In this study, the total nitrogen of soil was assumed to be uniform distribution across the field before sowing with low degree of variation of 1.2 %. And during the growing period from sowing (September 26, 2005) to jointing stage (April 13, 2006), no fertilizer was applied in this experiment field. Prior to the wheat harvest, soil total nitrogen (TN) samples were collected in total 100 experimental plots on June 14, 2006. Table 4 showed the best fit theoretical semi-variogram model and parameters for soil total nitrogen (TN). The Spherical semi-variogram model was found to be the best fit model for soil TN both at depths of 0-30 cm and 30-60cm. Nugget/sill ratio values for soil TN at depths of 0-30 cm and 30-60 cm were all less than 25%, indicating strong spatial dependence.

<table><caption>**Table 4** Best fit theoretical semi-variogram model and parameters for soil total nitrogen (TN)</caption><thead><tr><th></th><th>Model</th><th>Range(m)</th><th>Nugget</th><th>Sill</th><th>Nugget/Sill</th><th>R<sup>2</sup></th><th>RSS</th></tr></thead><tbody><tr><td>TN<br>(0-30cm)</td><td>Spherical</td><td>25.80</td><td>0.000005</td><td>0.000043</td><td>11.63%</td><td>0.75</td><td>1.18E-10</td></tr><tr><td>TN<br>(30-60cm)</td><td>Spherical</td><td>23.50</td><td>0.000008</td><td>0.000043</td><td>18.60%</td><td>0.64</td><td>1.27E-10</td></tr></tbody></table>
The spatial distribution of soil total nitrogen (TN) at depth of 0-30 cm and 30-60 cm were shown in Fig. 7. The map of soil TN content at depth of 0-30 cm were similar with the map of TN content at depth of 30-60 cm in the study area. Soil TN existed obviously strip spatial distribution characteristics in north-south direction both at the depth of 0-30 cm and 30-60 cm. There were three main higher TN content regions and four main lower TN content regions on maps, which were similar to the wheat yield spatial distribution maps. Except no fertilization treatment, other fertilization treatments had little difference to the soil TN content at harvest stage.

<figure><img src="image_8.png" /><figcaption>Fig.7. The spatial distribution of soil total nitrogen (TN) at depth of 0-30 cm and 30-60 cm</figcaption></figure>

## 4. Discussion

In-field, wheat growth and yield show obvious spatial variation due to spatial variability of soil nutrient. The results shown in this study is in accordance with previously reported results with similar spatial distribution pattern of wheat yield (Fig. 6) with NDVIs (Fig. 5) and soil total nitrogen (Fig. 7). And Soil TN content at depth 0-30 cm was strongly related with wheat yield and canopy NDVIs, except NDVI at June 8, 2006 (Table 5). That is to say, although using different fertilization treatments, the spatial distribution of crop growth and yield were still deeply influenced by soil nitrogen nutrient in the field.

<table><caption>Table 5 The relationship between NDVIs, yield and soil</caption><thead><tr><th></th><th>NDVI (13/04)</th><th>NDVI (21/04)</th><th>NDVI (29/04)</th><th>NDVI (10/05)</th><th>NDVI (29/05)</th><th>NDVI (08/06)</th><th>TN1</th><th>TN2</th><th>Yield</th></tr></thead><tbody><tr><td>NDVI (13/04)</td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NDVI (21/04)</td><td>0.691 <sup>**</sup></td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NDVI (29/04)</td><td>0.580 <sup>**</sup></td><td>0.757 <sup>**</sup></td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NDVI (10/05)</td><td>0.493 <sup>**</sup></td><td>0.571 <sup>**</sup></td><td>0.567 <sup>**</sup></td><td>1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NDVI (29/05)</td><td>0.316 <sup>**</sup></td><td>0.490 <sup>**</sup></td><td>0.544 <sup>**</sup></td><td>0.657 <sup>**</sup></td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>NDVI (08/06)</td><td>-0.197</td><td>-0.122</td><td>0.015</td><td>0.086</td><td>0.394 <sup>**</sup></td><td>1</td><td></td><td></td><td></td></tr><tr><td>TN (0-30cm)</td><td>0.415 <sup>**</sup></td><td>0.407 <sup>**</sup></td><td>0.421 <sup>**</sup></td><td>0.495 <sup>**</sup></td><td>0.348 <sup>**</sup></td><td>-0.092</td><td>1</td><td></td><td></td></tr></tbody></table>
<table><tr><td>TN (30-60cm)</td><td>0.201</td><td>0.215</td><td>0.207</td><td>0.278</td><td>0.186</td><td>-0.087</td><td>0.660</td><td>1</td></tr><tr><td>Yield</td><td>0.262</td><td>0.439</td><td>0.726</td><td>0.689</td><td>0.680</td><td>0.317</td><td>0.313</td><td>0.310</td></tr></table>

Note: * Significant at the 0.05 level; ** Significant at the 0.01 level

## 5. Conclusion

This study compared and assessed the potential of three different site-specific N fertilization decisions to improve crop growth status and grain production, with uniform fertilization and no fertilization as controls, by using traditional and geo-statistical method. The results of this study were as followings.

The fertilization strategy based on crop growth model (Z) had the best capacity to promote crop growth and grain yield, compared to other fertilization treatments (CK, Y and S), but S treatment had the biggest capacity to reduce the degree of variation in yield. Wheat canopy NDVIs, grain yield and soil TN on fertilization area had strong spatial dependence and the spatial variability of them were mainly resulted from structural factors instead of random factors. The spatial dependence of wheat growth was enhanced after fertilization. The spatial distribution of yield had some similar spatial distribution characteristics with NDVIs and soil TN on fertilization area. And the areas with stable high or low crop growth status after fertilization usually had high or low yield on different N fertilizer application zones.

**Funding Information.** This study was supported by National Key R&D Program of China (2016YFD0300601), National Natural Science Foundation of China (41501468, 41601466), the Agricultural Science and Technology Innovation of Sanya (2015KJ04), the Natural Science Foundation of Hainan Province, China (20164179, 2016CXTD015), the Technology Research, Development and Promotion Program of Hainan Province, China (ZDXM2015102), the Hainan Provincial Department of Science and Technology under Grant (ZDKJ2016021).

## References

1. Wang, Z., Li, S.: Effects of nitrogen and phosphorous fertilization on plant growth and nitrate accumulation in vegetables. J. Plant Nutr. 27, 539-556 (2004).

2. Albornoz, F., Lieth, J.H.: Over fertilization limits lettuce productivity because of osmotic stress. Chilean J. Agric. Res. 75, 284-290 (2015b).

3. Francisco, A.: Crop responses to nitrogen over fertilization: A review. Sci. Hortic. 205, 79-83 (2016).

4. Montemurro, F.: Different nitrogen fertilization sources, soil tillage, and crop rotations in winter wheat: effect on yield, quality, and nitrogen utilization. J. Plant Nutr. 32, 1-18 (2009).

5. Meyer-Aurich, A., Weersink, A., Gandorfer, M., & Wagner, P.: Optimal site-specific fertilization and harvesting strategies with respect to crop yield and quality response to nitrogen. Agr. Syst. 103, 478-485 (2010b).

6. Arne, M. R., Henning, K.: Predicting the site specific soil N supply under winter wheat in Germany. Nutr. Cycl. Agroecosys. 1-11 (2017).

7. Griepentrog, H. W. and Kyhn, M.: Strategies for site specific fertilization in a highly productive agricultural region. The 5th International Conference on Precision Agriculture, Minneapolis, USA, July 2000.
8. Delin, S., Lindén, B., Berglund, K.: Yield and protein response to fertilizer nitrogen in different parts of a cereal field: potential of site-specific fertilization. Eur. J. Agron. 22, 325-336 (2005).

9. Cui, Z.L., Zhang, F.S., Chen, X.P., Miao, Y.X., Li, J.L., Shi, L.W., Xu, J.F., Ye, Y.L., Liu, C.S., Yang, Z.P., Zhang, Q., Huang, S.M., Bao, D.J.: On-farm evaluation of an in-season nitrogen management strategy based on soil Nmin test. Field Crops Res. 105, 48-55 (2008).

10. Li, F., Miao, Y., Zhang, F., Cui, Z., Li, R., Chen, X., Zhang, H., Schroder, J., Raun, W.R., Jia, L.: In-season optical sensing improves nitrogen-use efficiency for winter wheat. Soil Sci. Soc. Am. J.73, 1566-1574 (2009).

11. Colac, A. F., Molin, J. P.: Variable rate fertilization in citrus: a long term study. Precis. Agric.18, 169-191 (2017).

12. Mariangela, D., Pietro, R., Francesco, M.: Precision nitrogen management of wheat. A review. Agronomy for Sustainable Development, 33 (1), 219-241(2013).

13. Clevers, J. G. P. W., Kooistra, L.: Using hyperspectral remote sensing data for retrieving canopy chlorophyll and nitrogen content. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(2), 574-583 (2012).

14. Clevers, J. G. P. W., Gitelson, A. A.: Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3. International Journal of Applied Earth Observation and Geo-information, 23, 344-351 (2013).

15. Li, F., Mistele, B., Hu, Y. C., Yue, X. L., Yue, S. C., Miao, Y. X., Chen, X. P., Cui, Z. L., Meng, Q. F. Schmidhalter, U.: Remotely estimating aerial N status of phenologically differing winter wheat cultivars grown in contrasting climatic and geographic zones in China and Germany. Field Crops Research, 138, 21-32 (2012).

16. Wang, W., Yao, X., Yao, X. F., Tian, Y. C., Liu, X. J., Wei, J. N., Cao, X., Zhu, Y.: Estimating leaf nitrogen concentration with three-band vegetation indices in rice and wheat. Field Crops Research, 129: 90-98 (2012).

17. Ladha, J. K., Pathak, H., Krupnik, T. J., Six, J., Kesse, C. V.: Efficiency of fertilizer nitrogen in cereal production: retrospects and prospects. Adv. Agron. 87, 85-156 (2005).

18. Montemurro, F., Maiorana, M., Ferri, D., Convertini, G.: Nitrogen indicators, uptake and utilization efficiency in a maize and barley rotation cropped at different levels and sources of N fertilization. Field Crop Res. 99, 114-124 (2006).

19. Cartelat, A., Cerovic, Z. G., Goulas, Y., Meyer, S., Lelarge, C., Prioul, J. L., Barbottin, A., Jeuffroy, M. H., Gate, P., Agati, G., Moya, I.: Optically assessed contents of leaf polyphenolics and chlorophyll as indicators of nitrogen deficiency in wheat (Triticum aestivumL.). Field Crop Res. 91, 35-49 (2005).

20. Ehlert, D., Schmerler, J., Voelker U.: Variable rate nitrogen fertilization of winter wheat based on a crop density sensor. Precis. Agric. 5, 263-273 (2004).

21. Godwin, R. J., Richards, T. E., Wood, G. A., Welsh, J. P., Knight, S. M.: An economic analysis of the potential for precision farming in UK cereal production. Biosyst Eng. 84, 533-545 (2003b).

22. Morris, K. B., Martin, K. L., Freeman, K. W., Teal, R. K., Girma, K., Arnall, D. B., Hodgen, P. J., Mosali, J., Raun, W. R., Solie J. B.: Mid-season recovery from nitrogen stress in winter wheat. J. Plant Nutr. 29, 727-745 (2006).

23. Singh, B., Sharma, R. K., Kaur, J., Jat, M. L., Martin, K. L., Singh, Y., Singh, V., Chandna, P., Choudhary, O. P., Gupta, R. K., Thind, H. S., Singh, J. U. H. S., Khurana, H. S., Kumar, A., Uppal, R. K., Vashistha, M., Raun, W. R., Gupta, R.: Assessment of the nitrogen management strategy using an optical sensor for irrigated
wheat. Agron. Sustain. Dev. 31(3), 589-603 (2011).

24. Thomason, W. E., Phillips, S. B., Davis, P. H., Warren, J. G., Alley, M. M., Reiter, M. S.: Variable nitrogen rate determination from plant spectral reflectance in soft red winter wheat. Precis. Agric. 12, 666-681(2011).

25. Gao, X. Z., Hu, K. L., Guo, Y., Li, B. G., Ma, Y. T., Du, S., Wang, Y. H.: Spatial variability of soil nutrients and crop yield and site-specific fertilizer management. Scientia Agriculture Sinica. 35(6):660-666 (2002). (in Chinese with English Abstract)

26. Fensholt, R., Proud, S.R. Evaluation of earth observation based global long term vegetation trends-comparing GIMMS and MODIS global NDVI time series. Remote Sens. Environ. 119, 131-147 (2012).

27. Wen, L., Saintilan, N., Yang, X., Hunter, S., Mawer, D.: MODIS NDVI based metrics improve habitat suitability modelling in fragmented patchy floodplains. Remote Sens. Appl.: Soc. Environ. 1, 85-97 (2015).

28. Gerardo, E. S., Christian, G. P.H., Ingo, J. H., Amanda, D. R., Pablo, M. V.: Tree senescence as a direct measure of habitat quality: Linking red-edge Vegetation Indices to space use by Magellanic woodpeckers. Remote Sens. Environ. 193, 1-10 (2017).

29. Bruno, B., Costanza, F., Davide, C., Urs, S.: Variable rate nitrogen fertilizer response in wheat using remote sensing. Precis. Agric. 1-15 (2015).

30. Rouse, J.W., Haas, R.H., Schell, J.A., Deering, D.W.: Monitoring vegetation systems in the Great Plains with ERTS. Proceedings of the 3rd ERTS Symposium. U.S. Government Printing Office, Washington D.C., 1309-1317 (1974).

31. Ding, X.D.: Comparisons among the methods of handling outliers. Computing techniques for geophysical and geochemical exploration 18(1), 71-77 (1996). (in Chinese with English Abstract)