Precision Agric (2009) 10:471–487 DOI 10.1007/s11119-009-9108-2 The delineation of agricultural management zones with high resolution remotely sensed data Xiaoyu Song Æ Jihua Wang Æ Wenjiang Huang Æ Liangyun Liu Æ Guangjian Yan Æ Ruiliang Pu Published online: 11 February 2009
Springer Science+Business Media, LLC 2009 Abstract Remote sensing (RS) techniques have been widely considered to be a promising source of information for land management decisions. The objective of this study was to develop and compare different methods of delineating management zones (MZs) in a field of winter wheat. Soil and yield samples were collected, and five main crop nutrients were analyzed: total nitrogen (TN), nitrate nitrogen (NN), available phosphorus (AP), extractable potassium (EP) and organic matter (OM). At the wheat heading stage, a scene of Quickbird imagery was acquired and processed, and the optimized soil-adjusted vegetation index (OSAVI) was determined. A fuzzy k-means clustering algorithm was used to define MZs, along with fuzzy performance index (FPI), and modified partition entropy (MPE) for determining the optimal number of clusters. The results showed that the optimal number of MZs for the present study area was three. The MZs were delineated in three ways; based on soil and yield data, crop RS information and the combination of soil, yield and RS information. The evaluation of each set of MZs showed that the three methods of delineating zones can all decrease the variance of the crop nutrients, wheat spectral parameters and yield within the different zones. Considering the consistent relationship between the crop nutrients, wheat yield and the wheat spectral parameters, satellite remote sensing shows promise as a tool for assessing the variation in soil properties and yield in X. Song (&)
J. Wang
W. Huang
L. Liu National Engineering Research Center for Information Technology in Agriculture, Beijing 100097, China e-mail: Songxy@nercita.org.cn; song.xiao.yu@163.com X. Song
G. Yan College of Geography/Research Center for Remote Sensing and GIS, Beijing Normal University, Beijing 100875, China X. Song
J. Wang
W. Huang
L. Liu Key Laboratory for Information Technologies in Agriculture, The Ministry of Agriculture, Beijing, China R. Pu Department of Geography, University of South Florida, Tampa, FL, USA 123 472 Precision Agric (2009) 10:471–487 arable fields. The results of this study suggest that management zone delineation using RS data was reliable and feasible. Keywords Quickbird imagery
Management zone (MZ)
Optimized soil-adjusted vegetation index (OSAVI)
Winter wheat Introduction Recent advances in technology for variable-rate applications in precision agriculture, combined with advances in global positioning systems (GPS), geographic information systems (GIS) and remote sensing (RS), have provided powerful analytical tools for farm management. This has been termed ‘precision crop management’ (PCM) and is defined as an information- and technology-based agricultural management system to identify, analyze and manage soil spatial and temporal variation within fields for optimum profitability, sustainability and protection of the environment. Variable-rate technology (VRT), probably the best-developed part of the PCM system (Searcy 1995), applies production inputs at rates appropriate to the soil and plant conditions within fields. Variable-rate technology has been shown to be successful for several agricultural products, such as herbicides (Mortensen et al. 1995), fertilizers (Schueller 1992), insecticides (Weisz et al. 1996) and seeds. A popular basis for applying VRT is the sitespecific management zone (MZ) method (Khosla et al. 2002; Koch et al. 2004). High quality information about the variation in soil and crop properties is the basis of effective management. There are several approaches for the delineation of site-specific MZs (Rodrigo and Oscar 2007). The first is based on soil information (Nolan et al. 2000; Fraisse et al. 2001). Sampling the soil on a grid to determine the chemical and physical properties is used widely for characterizing spatial variation and zone delineation (Franzen et al. 2000; Flowers et al. 2005). The second approach is based on yield maps from several seasons (Taylor et al. 2001; Blackmore 2000; Diker et al. 2004). Soil zones based on topography have also been used for determining management zones. More practical approaches, such as the farmer’s knowledge, have also been considered in delineating MZs. For example, Fleming et al. (2004) used a farmer’s experience in combination with aerial photographs to define areas of similar soil properties. Precision agriculture (PA) was introduced to China in the early 1990s and it has developed rapidly in recent years. For example, Hu et al. (1999) found that the scale of variation for soil nitrate- nitrogen (NN) was less 10 m, whereas that for ammonia–nitrogen, OM, AP and TN varied from 22.4 to 58.5 m. Gao et al. (2002) studied site-specific fertilizer management based on the spatial variation in crop nutrients and yield. Xue et al. (2004) applied fertilizer nitrogen for winter wheat at a variable-rate based on an 18 9 18 m2 grid-sampling of alkali-soluble soil nitrogen and a yield map. Li et al. (2005) delineated management zones for precision agriculture using yield data over 4 years at the National Experimental Station for Precision Agriculture of China. The research by Zhang et al. (2006) showed that variable-rate fertilizer application on maize resulted in a higher yield and smaller coefficient of variation (CV) for the crop nutrients than conventional uniform fertilizer application. In principle, these methods seem logical, but there are limitations. Soil sampling on a grid can provide basic information for delineating management zones, but the sampling density should be related to the degree of variation at the site. It should be intensive enough 123 Precision Agric (2009) 10:471–487 473 to be able to create sound management zones. For this reason, soil sampling is labor intensive, time consuming and sometimes temporally variable, which means that the soil needs to be sampled every growing season to determine crop nutrient levels in fields where there is variable-rate fertilizer application (Gotway et al. 1996; Fleming et al. 2001; Nolan et al. 2000; Koch and Khosla 2003). Yield maps produced from GPS-equipped grain yield monitors also provide information about the spatial variation of a field. Taylor et al. (2001) discussed several methods for using yield monitor data to develop yield goal maps. However, their yield prediction results were inconsistent across site-years, and they concluded that several years of data would be needed to develop appropriate yield goal maps for a given location. Remotely sensed images acquired by aircraft and satellite-based sensors have the potential to provide a synoptic view of an entire field, not just the sample sites, within the time and space requirements of PCM applications. Therefore, RS techniques have been considered as a promising source of information for land management decisions (Johannsen and Barney 1981). Images from sensors have a unique role for monitoring seasonally variable crop and soil conditions, and for time-specific and time-critical crop management. One of the strengths of RS is the ability to obtain spatially intensive information rapidly over large areas, thus creating the potential to supply information about spatial variation of N requirements, which is less expensive and more convenient than currently used sources. Considerable efforts have been devoted to detecting crop N stress using remote sensing. In general, different levels of nitrogen stress can be detected easily from aerial images (Blackmer et al. 1996; Ashcroft et al. 1990; Blackmer and White 1998). Nitrogen stress increases canopy reflectance over all visible wavelengths because of a shortage of chlorophylII and other light-absorbing pigments. Nitrogen stress may also decrease canopy reflectance in near-infrared wavelengths (Beatty et al. 2000; McMurtrey et al. 1994). Indices combining information from visible and near-infrared regions may maximize sensitivity to N stress. Wright et al. (2002) used Quickbird imagery and aerial photography provided by Environmental Mapping and Remote Sensing (EMARS) to manage wheat grain protein. Song et al. (2004) used three airborne Push-room Hyperspectral Image (PHI) sensor images to investigate the effect of soil nitrogen supplies and variable-rate fertilizer application on winter wheat growth. Recently, some attempts have been made to develop and evaluate decision algorithms for N management based on RS. Scharf and Lory (2002) developed a decision algorithm for side-dressing corn based on aerial photographs. This algorithm is based on green-light intensity from unfertilized corn relative to well-fertilized corn; the larger the difference in color, the more side-dressed N is recommended for the previously unfertilized corn. Liang et al. (2005) studied the algorithm for variable-rate nitrogen application based on the reflected spectrum from a wheat canopy, and the vegetation index OSAVI was used to determine the amount of nitrogen fertilizer recommended for variable-rate management. The results indicated that, compared to traditional uniform application, variable-rate fertilizer application reduced the variation in wheat yield, ear numbers and dry biomass, but it did not increase crop yield and grain protein content significantly. The objective of this study was to use multi-spectral and high spatial resolution RS data (Quickbird image) to generate within-field MZs. Three different MZ delineation methods that based on soil nutrient data, RS information, and soil, spectral and yield data were evaluated in this study. 123 474 Precision Agric (2009) 10:471–487 Materials and methods A field experiment at the National Experimental Station for Precision Agriculture of China was designed during the 2005/2006 winter wheat growing season. Soil samples were taken, and five crop nutrients, TN, NN, AP, EP and OM were determined using standard laboratory procedures (Zhang 2005a, b). Meanwhile, one scene of Quickbird imagery during the wheat’s heading stage was acquired and processed. The spectral parameter of the optimized soil-adjusted vegetation index (OSAVI) was extracted from the imagery. Study site This study was done at the National Experimental Station for Precision Agriculture, which is located in the Changping district of Beijing (4010.60 N, 11626.30 E). This site has a mean annual rainfall of 508 mm and a mean annual temperature of 13C. The winter wheat was sown on September 29, 2005 with a 15 cm row spacing. The wheat cultivar was Jingdong 8, which is one of the main winter wheat varieties in Northern China. The base fertilizer was applied on September 27, 2005 and supplementary fertilizer was applied on April 22, 2006. A square of 90 9 90 m2 in the wheat field was selected as the soil sampling area (Fig. 1). Data collection Soil sampling Soil sampling was carried out after the study area had been subdivided into 81 small plots of 100 m2. The first set of soil samples was taken on April 4, 2006 just before fertilizer application at the wheat elongation stage. The second set of soil samples was taken on June 16, 2006 just before the wheat harvest. The location of each soil sample was recorded by a GPS receiver with a station-based differential signal (DGPS, Trimble 5700 RTK). 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 and a depth of 0–30 cm was taken; there were 81 soil samples in total. Winter wheat yield data On June 15, 2006 geo-referenced wheat yield samples were collected and threshed manually to determine the yield. For each plot, three samples, which were randomly distributed on the plots, were collected and combined for a single plot yield. The sampling area for each wheat plot was 3 m2. All wheat grain was expressed on a dry matter basis after the samples were oven-dried at 60C. Quickbird imagery The Quickbird imagery (at wheat heading stage) covering the entire study site was recorded on May 2, 2006. Quickbird is a commercial satellite that records high-resolution imagery globally. Its onboard sensor has a multispectral scanner equipped with three visible (VIS), one near infrared (NIR) and one panchromatic (PAN) band. Visible and NIR bands acquire data with 2.4-m resolution, and PAN data are acquired at a 60- and 70-cm 123 Precision Agric (2009) 10:471–487 475 Fig. 1 Quickbird remotely sensed image of the study site spatial resolution. Data were recorded on days with \20% cloud cover, and as close to solar noon as possible. Geometric and atmospheric corrections were applied to the Quickbird imagery. Relatively coarse geometric correction was carried out by choosing ground control points (GCP) and the precise geometric correction was completed using DGPS GCPs. Atmospheric correction was done using the ‘empirical line method’ (ELM) (Farrand et al. 1994) based on ground spectral measurements taken with an ASD FieldSpec Pro spectrometer (Analytical Spectral Devices, Boulder, CO, USA). A cement field in the National Experimental Station for Precision Agriculture was selected as the bright object, and a reservoir with pure water near the station was selected as the dark object. After the geometric and atmospheric corrections, a pan-sharpening process was applied to the Quickbird’s four multi-spectral bands by using the pan band. Consequently, the Quickbird image has four 0.6 m nominal spatial resolution bands. 123 476 Precision Agric (2009) 10:471–487 Extraction of vegetation indices Soil reflectance is often greater than crop reflectance in the visible wavelengths and can, therefore, interfere with remote estimates of crop reflectance. The optimized soil-adjusted vegetation index (OSAVI) was developed by Rondeaux et al. (1996) using bidirectional reflectance in the near-infrared (NIR) and red (R) bands: OSAVI ¼ ð1 þ 0:16Þðqnir  qr Þ=ðqnir þ qr þ 0:16Þ; ð1Þ where qnir and qr are the reflectance of near infrared and red bands, respectively. The soil adjustment coefficient (0.16) was selected as the optimal value to minimize variation arising from the soil background. The residual variation in OSAVI caused by the soil was spread evenly across the full range (0–1) of crop groundcover, making this index particularly suitable for agricultural applications. The OSAVI index is closely related to many vegetation properties, such as leaf area index, vegetation cover, vegetation biomass and crop growth, therefore, it is often used to monitor crop growth and to predict crop yield. In the present study, OSAVI was calculated to reflect wheat cover and growth. Methods of analysis Statistical analyses, including the Kappa coefficient, were performed to make a systematic comparison of RS-, soil and yield- and RS, soil and yield-based methods for generating MZs. (Ahn and Mezzich 1989; Erik 1996). Ordinary kriging and fuzzy k-means clustering were also used in this study. Conventional statistics were computed with Excel 2003. Variowin 2.2 (Pannatier 1996) was used to compute the variograms and Arcgis9.0 was used for kriging (ESRI 2001). The Quickbird image was analyzed and displayed with ENVI 4. 1 (ITT 2003), and Fuzme3.5 was used for the fuzzy k-means clustering algorithm in this study (Minasny 1999). Geostatistical analysis The geostatistical tools of variography and kriging were used in this study (Edward and Srivastava 1989; Zhang 2005a, b). To characterize the spatial variation of the soil properties and wheat yield data, experimental variograms were computed on the selected soil variables and yield data using Variowin 2.2. Isotropic models were fitted to the experimental variograms to determine the type of spatial structure. The parameters of the fitted models were then used with the data for ordinary kriging in Arcgis9.0. The values of the selected soil properties and yield were predicted on a 0.6 m grid at un-sampled locations and the predictions were used to produce contour maps. Fuzzy k-means clustering algorithm One approach to fuzzy classification is fuzzy k-means (de Gruijter and McBratney 1988). It partitions a set of data {x1, x2, … xn} into k clusters to minimize the squared distances between a point and the class centroid within the cluster. The objective function is: J¼ c n X X i¼1 subject to the following conditions: 123 k¼1 m/ik d2 ðxi ck Þ; ð2Þ Precision Agric (2009) 10:471–487 477 k X mik ¼ 1 i ¼ 1; 2; . . .; n; k¼1 n X mik [ 0 k ¼ 1; 2; . . .; c; i¼1 mik 2 f0; 1g i ¼ 1; 2; . . .; n; k ¼ 1; 2; . . .; c; ð3Þ where n is the number of data, c is the number of classes, mik is the fuzzy membership matrix, / is fuzziness weighting exponent, ck is the vector representing the centroid of class k, xi is the vector representing an individual datum i, and d2 ðxi ck Þis the squared distance between xi and ck according to a chosen definition of distance. The latter, denoted by dik2 / for simplicity, is the fuzzy exponent and it ranges from 1 to . It determines the degree of fuzziness of the final solution, which is the degree of overlap between groups. With / ¼ 1; the solution is a hard partition. As / approaches infinity, the solution approaches its greatest degree of fuzziness. Minimization of the objective function, J, provides the solution for the membership function (Bezdek 1981): 2=ð/1Þ d mik ¼ Pc ik j¼1 i ¼ 1; 2; . . .; n; k ¼ 1; 2; . . .; c 2=ð/1Þ ð4Þ dij Pn ck ¼ Pi¼1 n m/ik xi i¼1 m/ik k ¼ 1; 2; . . .; c: ð5Þ To determine the optimal number of classes, the fuzziness performance index (FPI) and modified partition entropy (MPE) (Roubens 1982) were used. The FPI is defined as: FPI ¼ 1  cF  1 ; c1 ð6Þ where c is the number of classes and F is the partition coefficient. F¼ c n X 1X ðmik Þ2 : n i¼1 k¼1 ð7Þ The modified partition entropy (MPE) is given by: MPE ¼ H ; log c ð8Þ where H is the entropy function: H¼ n X c 1X mik  logðmik Þ: n i¼1 k¼1 ð9Þ The FPI estimates the degree of fuzziness generated by a specified number of classes. Values of FPI and MPE may range from 0 to 1. Values approaching 0 indicate distinct classes with little membership sharing, whereas values near 1 indicate non-distinct classes with a large degree of membership sharing. The MPE estimates the degree of disorganization created by a specified number of classes. The optimal number of clusters for each index occurs when the index is at the minimum, representing the least membership sharing 123 478 Precision Agric (2009) 10:471–487 and the greatest amount of organization as a result of the clustering process (Fridgen et al. 2004). Methods of generating management zones Based on the soil, yield and OSAVI data, MZs were generated by the fuzzy k-means algorithm. Three sets of data were used to generate the MZs and the results were compared in this study. Firstly, the kriged predictions for soil OM, AP, EK and wheat yield were selected as the input data for clustering and to generate the MZ maps. Secondly, RS data in the form of OSAVI were used to generate an MZ map. Lastly, all soil, yield and RS data were combined to generate an MZ map. Results and discussion Summary statistics for soil, wheat yield and RS data The soil sample data collected on April 4 and June 16, 2006 were analyzed statistically (Table 1) before the geostatistical analysis. In accordance with the nutrient classification standard of China’s second soil survey (Guo et al. 2004), soil AP content at the study site was at a high level(150–200 mg kg-1); OM was at a moderate level (0.10–0.15%); EK was at a low level (5–10 mg kg-1) and NN was also at a low level. Wheat yield can be affected by many factors, such as crop nutrient status, fertilizer regime, irrigation, weather, etc. In this study, the relations among the wheat yield, crop nutrients and wheat spectral information were analyzed; the results are given in Table 2. The results show that the wheat’s Quickbird OSAVI have large correlation coefficients with wheat yield and nutrient status. Wheat yield also has strong positive correlations with soil AP and TN on April 4, as well as soil OM and TN on June 16. Soil EK has weak correlations with yield and OSAVI. The correlation of 0.94 between OM and TN measured on June 16 is strong. The coefficient between June OM and April TN is 0.71 and 0.72 for Table 1 Results of statistical analysis of soil properties and wheat yield Soil nutrients Minimum Maximum Mean Standard deviation CV (%) April-AP (mg kg-1) 2.78 15.82 8.13 3.12 38.43 April-NN (mg kg-1) 10.03 44.46 20.79 7.30 35.09 0.09 0.13 0.11 0.01 9.69 102.80 237.20 148.02 23.06 15.58 June-AP (mg kg-1) 3.21 16.51 8.39 2.72 32.41 June-NN (mg kg-1) 0.96 15.13 3.93 2.27 57.70 June-TN (%) 0.09 0.13 0.10 0.01 9.00 86.42 157.13 112.95 12.97 11.49 April-TN (%) April-EK (mg kg-1) June-EK (mg kg-1) June-OM (%) Yield (kg ha-1) 1.43 2937.5 2.27 5986.7 1.81 4591.9 0.18 9.97 722.54 15.73 Note: n = 81, n number of soil samples, TN soil total nitrogen, NN soil nitrate nitrogen, AP soil available phosphorus, EK soil extractable potassium, OM soil organic matter, CV coefficient of variation 123 Precision Agric (2009) 10:471–487 479 Table 2 Correlation coefficients between crop nutrients and Quickbird spectral data, OSAVI and wheat yield data Crop nutrients Wheat yield (kg ha-1) Band1 Band2 Band3 Band4 OSAVI April-AP (mg kg-1) 0.43 -0.44 -0.45 -0.45 0.56 0.51 April-TN (%) 0.30 -0.23 -0.22 -0.26 0.36 0.31 -0.31 April-EK (mg kg-1) -0.02 0.31 0.35 0.30 -0.31 June-AP (mg
kg-1) 0.24 -0.39 -0.44 -0.41 0.47 0.45 June-TN (%) 0.37 -0.39 -0.41 -0.41 0.51 0.47 June-OM (%) 0.41 -0.42 -0.44 -0.44 0.54 0.50 Wheat yield (kg ha-1) 1.00 -0.55 -0.64 -0.70 0.75 0.76 Note: TN soil total nitrogen, NN soil nitrate nitrogen, AP soil available phosphorus, EK soil extractable potassium, OM soil organic matter June TN and April TN. Based on the strong correlation with yield, the soil AP and EK on April 4 and OM on June 16 were selected as variables for kriging. Kriging of crop nutrients The spatial structure of soil AP and EK on April 4 and OM on June 16 was evaluated by isotropic variogram models. The experimental variograms and the fitted models for the three soil variables and wheat yield are shown in Fig. 2. The model types and their parameters are given in Table 3. The degree of spatial dependence was determined by the nugget: sill ratio, which represented the percentage of unexplained variance (nugget) to the total variation, expressed by nugget ? spatially correlated variance (the sill). A nugget: sill ratio \25% indicates strong spatial dependence, whereas one [75% indicates weak spatial dependence (Cambardella et al. 1994). In all cases the soil properties show moderate to Fig. 2 Experimental variogram (symbols) and fitted model (solid line) of soil: (a) OM, (b) AP, (c) EK and (d) wheat yield 123 480 Precision Agric (2009) 10:471–487 Table 3 Isotropic variogram model parameters for crop nutrients and wheat yield Nugget (c0) Spatially dependent c0/(c0 ? c1) Range (m) RSME component (c1) Crop nutrients Model June-OM (%) Spherical 0.0068 0.0312 0.1789 38.38 April-AP (mg kg-1) Spherical 3.117 1.947 0.6155 43.15 0.9644 April-EK (mg kg-1) Spherical 1.0050 196.1 360.1 0.3526 22.67 1.0070 Wheat yield (kg ha-1) Spherical 5654.6 4324.0 0.5667 45.14 1.0040 Note: AP soil available phosphorus, EK soil extractable potassium, OM soil organic matter strong spatial dependence (Table 3). The ranges of spatial dependence vary between 17 and 41 m (Table 3). Based on the models and parameters given in Table 3, ordinary kriging was applied to the three soil properties AP, AK, OM and wheat yield. The prediction grid was 0.6 9 0.6 m2, which was the same as the resolution of Quickbird OSAVI. This enabled the plots to be divided into several classes to determine homogenous zones. The contour maps obtained for the four variables are shown in Fig. 3. It can be seen from Fig. 3 that large values of OM and AP are in the middle-northern part of the study area, and small values occur in the middle. The large values of EK are in the (a) OSAVI (a) (c) (b) Soil OM (%) 0.31-0.40 1.49 1.61 0.40-0.45 1.61 1.72 0.45-0.50 1.72 1.84 0.50-0.55 1.84 1.96 0.55-0.60 1.96 2.08 0.60-0.71 2.08 2.20 (b) (e) (d) (c) Soil AP (mg·kg -1 ) (d) Soil EK (mg·kg -1 ) (e) Wheat Yield (kg·ha 4.28-5.02 115.01-129.16 3727.11-3994.78 5.02-5.76 129.16-143.31 3994.78-4262.34 5.76-6.50 143.31-157.46 4262.34-4529.90 6.50-7.24 157.46-171.61 4529.90-4797.47 7.24-7.98 171.61-185.76 4797.47-5065.02 7.98-8.73 185.76-199.91 5065.02-5332.58 -1 ) Fig. 3 (a) Map of OSAVI, and kriged maps of: (b) soil organic matter (OM), (c) soil available P (AP), (d) extractable K (EK) and (e) wheat yield 123 Precision Agric (2009) 10:471–487 481 middle and small values are in the south of the study area. The OSAVI image is similar to the map of wheat yield, with large values in the NW and SW of the study site, and small values in the middle. The soil OM and AP also have a similar pattern of variation to those of yield and OSAVI. The results of the correlation analysis in Table 2 show that OSAVI has moderate to strong correlations with the interpolated values of soil AP, OM, EK and wheat yield; the correlation coefficients are 0.62, 0.53, -0.47 and 0.72, respectively. As an important index to reflect crop cover and growth, OSAVI was used together with the kriged values of OM, AP, EK and wheat yield for fuzzy k-means cluster analysis to identify areas that are homogenous in terms of landscape and soil conditions for site-specific management. Determining the optimal number of MZs The OSAVI and kriged values of soil AP, EK, OM and wheat yield (Fig. 3) were analyzed with the fuzzy k-means algorithm to assign each pixel of 0.6 9 0.6 m2 to an appropriate class based on the minimization of the generalized least squares error function of each pixel. We considered six clusters to be the maximum number of practical management zones. To determine the optimal number of classes we increased the number of classes one at a time from two to six. The functions FPI and MPE were used to determine the optimum number of clusters. The results from the two indices are plotted in Fig. 4. The optimal number of clusters for each index occurs when the index is at the minimum. Figure 4 shows that FPI and MPE had the same change in trend with an increase in cluster number, and the minimum FPI and MPE were obtained with 3 clusters for our study area. This indicates that the sum of squares for members within a cluster is minimized while the sum of squares between members of different clusters is maximized. The final decision of how many clusters to use for creating management zones when the two cluster validity indices are dissimilar may require additional verification. For example, when developing productivity zones, the optimal number of clusters might be determined by comparing the within-zone variation in yield with the increase in the number of clusters (Fridgen et al. 2004). 0.55 0.60 0.55 0.50 0.50 0.40 MPE FPI 0.45 0.45 0.40 0.35 0.35 0.30 0.30 2 3 4 5 6 Number of clusters Soil and yield based MPE OSAVI based MPE Soil, yield and OSAVI based MPE Soil and yield based FPI OSAVI based FPI Soil, yield and OSAVI based FPI Fig. 4 Graphs of fuzzy performance index (FPI) and modified partition entropy (MPE) plotted against number of clusters based on the results of fuzzy k-means classification using different combinations of data (soil, yield and OSAVI) 123 482 Precision Agric (2009) 10:471–487 Management zone delineation and evaluation Based on the optimal number of classes, three MZ maps were generated using different data; soil and yield, RS data and RS, soil and yield. Figure 5 shows the resulting maps. The Kappa coefficient was then used to compare the homogeneity of the zones in the three different MZ maps. The results show that zones based on soil and yield (Fig 5a), and soil, yield and RS (Fig. 5b) data are the most similar; the similarity between Fig. 5a and c is the greatest with a Kappa coefficient 0.91. The Kappa coefficient based on a comparison of the zones in Fig. 5a and b was 0.16. The results indicate a poor agreement between the zones in Fig. 5 a and b, and for Fig. 5b and c. Within each class, crop and yield parameters could be quantified, analyzed and compared over other zones. The statistical analyses indicate significant differences between the crop nutrients and yield in each zone of three maps. The MZ 3 had the largest nutrient (a) (c) (b) Management zone 1 Management zone 2 Management zone 3 (a) MZ map based on soil and yield data; (b) MZ map based on RS data; (c) MZ map based on soil, yield and RS data Fig. 5 Maps of management zones generated by fuzzy k-means classification using different combinations of data (soil, yield and OSAVI) Table 4 Coefficients of variance for MZs based on different methods of generating the fuzzy classes by fuzzy k-means Method of MZ generation Soil and yield based MZ RS (OSAVI) based MZ Soil, yield and RS based MZ MZ OM (%) AP (mg kg-1) EK (mg kg-1) OSAVI Mean CV (%) Mean Mean CV (%) Mean CV (%) Mean 1 1.69 4.67 5.17 7.88 162.1 6.50 0.62 10.83 4192.6 6.29 2 1.80 3.78 6.31 6.71 140.8 7.99 0.71 6.79 4616.3 4.94 3 1.98 6.38 7.91 5.23 147.0 5.41 0.74 6.36 4926.7 3.37 1 1.69 4.04 5.43 16.06 159.8 6.61 0.53 8.95 4016.9 6.33 2 1.77 7.22 6.06 15.58 150.7 9.94 0.66 4.41 4447.1 6.52 3 1.86 6.53 6.96 12.36 141.8 6.74 0.74 7.82 4760.0 10.26 1 1.69 4.80 5.21 9.21 162.3 6.44 0.61 10.69 4185.6 6.16 2 1.80 3.82 6.27 6.90 140.7 7.96 0.71 5.51 4617.8 4.87 3 1.97 6.40 7.88 5.46 146.8 5.42 0.74 6.13 4922.4 3.36 CV (%) Yield (kg ha-1) CV (%) Note: CV coefficient of variation, MZ management zone, MZ 1 management zone1, MZ 2 management zone2, MZ 3 Management zone3, AP soil available phosphorus, EK soil extractable potassium, OM soil organic matter 123 483 CV (%) Precision Agric (2009) 10:471–487 CV (%) MZ 1 CV (%) MZ 2 MZ 3 CV for Whole study site CV for OSAVI based MZ CV for Soil and yield based MZ CV for Soil, yield and OSAVI based MZ Fig. 6 Histograms of CV for the whole study site and for the three MZs based on different combinations of data (soil, yield and OSAVI) 123 484 Precision Agric (2009) 10:471–487 status and potential crop productivity, whereas MZ 1 had the lowest. Table 4 shows the results of analyzing the variation in soil properties, yield and OSAVI in all three zones for the three different maps (Fig. 5). Figure 6 shows the histograms of the CVs for all variables used in clustering for the three MZ maps. From this figure, we can see that the CV for soil AP decreases dramatically in all MZs for the three maps. Although it increases for soil EK in MZ 2 of Fig. 5b together with OSAVI in MZ 1 of Fig. 5a and c, the CV for wheat yield and soil OM decrease to different degrees in all MZs. To evaluate the three methods of delineating MZs, the CVs for soil, wheat yield and RS data for the whole study site (before interpolation) were compared with those for each zone and in the three MZ maps (after interpolation). The results of the percentage decrease in CV for the whole site and for the zones are given in Table 5. This table shows that the CVs of crop OM and EP decrease considerably in the three MZs where EK increases by 4.94% in MZ 2 and decreases by 30.2% and 28.8% for MZ 1 and MZ 3, respectively, of the RSbased map (Fig. 5b). The CVs for OSAVI drops most in MZs of the RS-based map (Fig. 5b). It increases in MZ 1 and decreases in other two zones of the soil-yield based and soil-yield-RS based maps (Fig. 5a, c). The results also show that the CVs for wheat yield in the three zones for all maps decrease, especially for the soil and yield-based map and soil, yield and RS-based map. They decrease by 18.97%, 36.46%, 56.60% for the soil and yieldbased map (Fig. 5a) and by 20.69%, 37.25%, 56.78% for the soil, yield and RS-based map (Fig. 5c). In the RS based map yield CVs decrease by 18.49%, 16.11% and 31.43% for the three zones, respectively. Discussion and application of the method By contrast with traditional VRT decisions based on soil sampling and yield data, remotely sensed imagery can provide within-season monitoring to aid management of the growing crop. Our study proposed optimal procedures for delineating MZs with Quickbird imagery on wheat at the heading stage. In addition, the method is reliable compared with two other methods to delineate MZs based on soil and yield, and soil, yield and RS data. The RS method for MZ delineation described in this research is based on one field and this needs to Table 5 Percentage difference between the CV for soil, RS and wheat yield data for the whole site before interpolation and after classification of the interpolated data Method of MZ generation Soil and yield based MZ RS based MZ Soil, yield and RS based MZ MZ OM (%) AP (mg kg-1) EK (mg kg-1) OSAVI Yield (kg ha-1) 1 38.52 66.42 31.31 -3.35 18.97 2 50.31 71.40 15.61 35.14 36.46 3 16.10 77.69 42.83 39.31 56.60 1 21.97 31.55 30.20 14.57 18.49 2 5.02 33.57 -4.94 57.92 16.11 3 11.65 47.33 28.81 56.25 31.43 1 36.82 60.74 31.94 -2.00 20.69 2 49.75 70.57 15.88 47.42 37.25 3 15.87 76.69 42.69 41.48 56.78 Note: CV coefficient of variation, MZ management zone, MZ 1 management zone1, MZ 2 management zone2, MZ 3 management zone3, AP soil available phosphorus, EK soil extractable potassium, OM soil organic matter 123 Precision Agric (2009) 10:471–487 485 be extended by further intensive and in-depth research of a larger area based on RS imagery with a different spatial resolution. In addition, there should be further studies to include other environmental features that influence crop growth, such as irrigation water, plant diseases and insect pests. The application and use of this method for farm managers may require help from professional personnel for RS data processing and analysis. Therefore, it is necessary to develop simple software for users to popularize this method in the future. Conclusion The spatial variation in soil, yield and RS data were well represented using geostatistics. The fuzzy k-means clustering approach was then used to classify the kriged soil, yield and RS values to create homogeneous management zones. For the experimental area, Quickbird imagery confirmed that the variation in wheat growth estimated by the vegetation index OSAVI was similar to the variation in the soil and yield data. Statistical analyses showed that the Quickbird OSAVI data can not only reflect the spatial variation in wheat growth during the early growing stage, but also the spatial variation in soil properties and yield. Therefore, it could be used to provide a valuable diagnostic to guide the generation of MZs. Satellite remote sensing shows promise as a tool for assessing the variation in soil properties and yield over arable fields. The results of this study suggest that management zone delineation using RS data was reliable and feasible. Acknowledgements This work was financially supported by 863 Program (2006AA10A302, 2006AA12Z138, 2006AA210108) from the National High-Tech R&D Program of China, the 973 Program from the State Basic Research Development Program of China (No. 2007CB714406), National Natural Science Foundation of China (40701120), and Beijing Municipal Natural Science Foundation (4092017). References Ahn, C. W., & Mezzich, J. E. (1989). PROPOV-K: a FORTRAN program for computing a kappa coefficient using a proportional overlap procedure. Computers and Biomedical Research, 22, 415–423. doi: 10.1016/0010-4809(89)90035-9. Ashcroft, P. M., Catt, J. A., Curran, P. J., Munden, J., & Webster, R. (1990). The relation between reflected radiation and yield on the broadbalk winter wheat experiment. International Journal of Remote Sensing, 10, 1821–1836. Beatty, M. K., Johannsen, C. J., & Ross, K. (2000). In situ detection of leaf chlorophyll content and leaf nitrogen content in Zea mays L. with remote sensing. In 2000 annual meeting abstracts, Minneapolis, MN. ASA, CSSA, and SSSA, Madison, WI, pp. 281. Bezdek, J. C. (1981). Pattern recognition with fuzzy objective function algorithms. New York: Plenum Press. Blackmer, A. M., & White, S. E. (1998). Using precision farming technologies to improve management of soil and fertilizer nitrogen. Australian Journal of Agricultural Research, 49, 555–564. Blackmer, T. M., Schepers, J. S., Varvel, G. E., & Meyer, G. E. (1996). Analysis of aerial photography for nitrogen stress within corn fields. Agronomy Journal, 88, 729–733. Blackmore, S. (2000). The interpretation of trends from multiple yield maps. Computers and Electronics in Agriculture, 26, 37–51. doi:10.1016/S0168-1699(99)00075-7. Cambardella, C. A., Moorman, T. B., Novak, J. M., Parkin, T. B., Kalen, D. L., Tuco, R. F., et al. (1994). Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal, 58, 1501–1511. de Gruijter, J. J., & McBratney, A. B. (1988). A modified fuzzy k means for predictive classification. In H. H. Bock (Ed.), Classification and related methods of data analysis (pp. 97–104). Amsterdam: Elsevier Science. 123 486 Precision Agric (2009) 10:471–487 Diker, K., Heermann, D. F., & Brodahl, M. K. (2004). Frequency analysis of yield for delineating yield response zones. Precision Agriculture, 5, 435–444. doi:10.1007/s11119-004-5318-9. Edward, H. I., & Srivastava, R. M. (1989). Applied Geostatistics (pp. 141–163). Oxford, NY: Oxford University Press. ENVI 4.1. (2003). ENVI user’s guide. Boulder: ITT Industries, Inc. Erik, N. (1996). Use of the weighted Kappa coefficient in classification error assessment of thematic maps. International Journal of Geographical Information Science, 10(5), 591–603. doi:10.1080/ 02693799608902099. ESRI. (2001). Using ArcGIS geostatistical analyst. Redlands, CA: ESRI, Inc. Farrand, W. H., Singer, R. B., & Merenyi, E. (1994). Retrieval of apparent surface reflectance from AVIRIS data - a comparison of empirical line, radiative transfer, and spectral mixture methods. Remote Sensing of Environment, 47, 311–321. doi:10.1016/0034-4257(94)90099-X. Fleming, K. L., Heermann, D. F., & Westfall, D. G. (2004). Evaluating soil color with farmer input and apparent soil electrical conductivity for management zone delineation. Agronomy Journal, 96, 1581– 1587. Fleming, K. L., Westfall, D. G., & Wausch, W. C. (2001). Evaluating management zone technology and grid soil sampling for variable rate nitrogen application. In precision agriculture [CD-ROM]. Proceedings of the 5th international conference, Minneapolis, MN, ASA, CSSA, and SSSA, Madison, WI, 16–19. Flowers, M., Weisz, R., & White, J. G. (2005). Yield-based management zones and grid sampling strategies: Describing soil test and nutrient variability. Agronomy Journal, 97, 968–982. Fraisse, C. W., Sudduth, K. A., & Kitchen, N. R. (2001). Delineation of site-specific management zones by unsupervised classification of topographic attributes and soil electrical conductivity. Transactions of the ASAE, 44, 155–166. Franzen, D. W., Halvorson, A. D., & Hoffman,V. L. (2000). Management zones for soil N and P levels in the Northern Great Plains. In P. C. Robert et al. (Eds.), Proceedings of the 5th international conference on precision agriculture, Bloomington, MN [CD-ROM]. ASA, CSSA, and SSSA, Madison, WI. Fridgen, J. J., Kitchen, N. R., Sudduth, K. A., Drummond, S. T., Wiebold, W. J., & Fraisse, C. W. (2004). Management zone analyst (MZA): Software for subfield management zone delineation. Agronomy Journal, 96, 100–108. Gao, X. Z., Hu, K. L., Guo, Y., Li, B. G., Ma, Y. T., Du, S., et al. (2002). Spatial variability of soil nutrients, crop yield and site-specific fertilizer management. Scientia Agricultura Sinica, 35, 660–666 (in Chinese). Gotway, C. A., Fergusong, R. B., & Hergert, G. W. (1996). The effects of mapping scale on variable-rate fertilizer recommendations for corn. In P. C. Robert, R. H. Rust, & W. E. Larson (Eds.), Precision agriculture, Proceedings of the 3rd international conference (pp. 321–330), Minneapolis, MN, ASA, CSSA, and SSSA, Madison, WI, 23–26. Guo, Y. F., Cui, K. Z., Zhang, L. J., Wu, Z. R., Fen, X. L., Zhen, B. F., et al. (2004). Dynamic study of soil furrow layer nutrient on south of Tianjin. Science and Technology of Tianjin Agricultural and Forestry, 177, 31–33 (in Chinese). Hu, K. L., Li, B. G., Lin, Q. M., Li, G. T., & Chen, D. L. (1999). Spatial variability of soil nutrient in wheat field. Transactions of the Chinese Society Agricultural Engineering, 15(3), 33–38 (in Chinese). Johannsen, C. J., & Barney, T. W. (1981). Remote sensing applications for resource management. Journal of Soil and Water Conservation, 36, 128–131. Khosla, R., Fleming, K., Delagado, J. A., Shaver, T. M., & West-fall, D. G. (2002). Use of site-specific management zones to improve nitrogen management for precision agriculture. Journal of Soil Water Conservation, 57(6), 513–518. Koch, B., & Khosla, R. (2003). The role of precision agriculture in cropping systems. Journal of Crop Production, 8, 361–381. doi:10.1300/J144v09n01_02. Koch, B., Khosla, R., Frasier, W. M., Westfall, D. G., & Inman, D. (2004). Economic feasibility of variablerate nitrogen application utilizing site-specific management zones. Agronomy Journal, 96, 1572–1580. Li, X., Pan, Y. C., Zhao, C. J., Wang, J. H., Bao, Y. S., & Wang, J. D. (2005). Delineation and scale effect of precision agriculture management zones using yield monitor data over four years. Scientia Agricultura Sinica, 38, 1825–1833 (in Chinese). Liang, H. X., Zhao, C. J., Huang, W. J., Liu, L. Y., Wang, J. H., & Ma, Y. H. (2005).Variable-rate nitrogen application algorithm based on canopy reflected spectrum and its influence on wheat. In Multispectral and hyperspectral remote sensing instruments and applications II, Proceedings of international society for optical engineering, (SPIE), 2005 (Vol. 5655, pp. 522–530). Bellingham, WA. McMurtrey, III, J. E., Chappelle, E. W., Kim, M. S., Meisinger, J. J., & Corp, L. A. (1994). Distinguishing nitrogen fertilization levels in field corn (Zea mays L.) with actively induced fluorescence and passive reflectance measurements. Remote Sensing of the Environment, 47, 36–44. 123 Precision Agric (2009) 10:471–487 487 Minasny, B. (1999). A computer program for fuzzy k-means analysis. Updated July 1, 2002, from http://www.usyd.edu.au/su/agric/acpa/fkme/FuzME.html. Mortensen, D. A., Johnson, G. A., Wyse, D. Y., & Martin, A. R. (1995). Managing spatially variable weed populations. In P. C. Robert, R. H. Rust, & W. E. Larson (Eds.), Proceedings of the 2nd international conference on site specific management for agricultural systems (pp. 397–415). Minneapolis, MN. Madison, WI: American Society for Agronomy. Nolan, S. C, Goddard, T. W, Lohstraeter, G., & Coen, G. M. (2000). Assessing management units on rolling topography. In P. C. Robert, R. H. Rust, & W. E. Larson (Eds.), Proceedings of the 5th international conference on precision and other resource management. CD-ROM. ASA/CSSA/SSSA, Madison, WI, USA. Pannatier, Y. (1996). Variowin Software for spatial data analysis in 2D. New York: Springer. Rodrigo, A. O., & Oscar, A. S. (2007). Determination of management zones in corn (Zea mays L.) based on soil fertility. Computers and Electronics in Agriculture, 58, 49–59. doi:10.1016/j.compag.2006.12.011. Rondeaux, G., Steven, M., & Baret, F. (1996). Optimization of soil adjusted vegetation indices. Remote Sensing of Environment, 55, 95–107. doi:10.1016/0034-4257(95)00186-7. Roubens, M. (1982). Fuzzy clustering algorithms and their cluster validity. European Journal of Operational Research, 10(3), 294–301. Scharf, P. C., & Lory, J. A. (2002). Calibrating corn color from aerial photographs to predict sidedress N need. Agronomy Journal, 94, 397–404. Schueller, J. K. (1992). A review and integrating analysis of spatially-variable control of crop production. Fertilizer Resolution, 33, 1–34. doi:10.1007/BF01058007. Searcy, S. W. (1995). Engineering systems for site-specific management: Opportunities and limitations. In Proceedings of site-specific management for agriculture system (pp. 603–608). Minneapolis, MN, ASA-CSSA-SSSA, Madison,WI. Song, X. Y., Wang, J. H., Xue, X. Z., Liu, L. L., Chen, L. P., & Zhao, C. J. (2004). Assessment of the influence of soil nitrogen supplies and variable fertilization on winter wheat growth condition using airborne hyperspectral image. Transactions of the Chinese Society of Agricultural Engineering, 20(4), 45–49 (in Chinese). Taylor, R. K., Kluitenberg, G. J., Schrock, M. D., Zhang, N., Schmidt, J. P., & Havlin, J. L. (2001). Using yield monitor data to determine spatial crop production potential. Transactions of the ASAE, 44, 1409– 1414. Weisz, R., Fleischer, S.J., & Smilowitz, Z. (1996). Site-specific integrated pest management for high-value crops: impact on potato pest management. Journal of Economic Entomology, 89, 501D–509. Wright, D. L., Rasmussen, V. P., Baker, D. J. (2002). Using remote sensing to manage wheat grain protein. NASA SSC Report, ARC-USU-001-02. Affiliated Research Center Final Reports [CD-ROM]. Earth Science Applications Directorate, National Aeronautics and Space Administration, John C. Stennis Space Center, Mississippi. Xue, X. Z., Chen, L. P., Sun, Z. G., & Zhao, C. J. (2004). Results of variable-rate nitrogen fertilization of winter wheat based on soil fertility and yield map. Transactions of the Chinese Society of Agricultural Engineering, 20(3), 59–62 (in Chinese). Zhang, R. D. (2005a). Theory and application of spatial variability. Beijing: Science Press (in Chinese). Zhang, X. F. (2005b). Chemical analysis for agricultural (pp. 158–169). Beijing: Chemical Industry Press (in Chinese). Zhang, S. H., Ma, C. L., Li, W., & Xu, Y. (2006). Experimental study on the influence of variable rate fertilization on Maize yield and soil nutrients. Transactions of the Chinese Society of Agricultural Engineering, 22(8), 64–67 (in Chinese). 123