Detailed Investigation of Spectral Vegetation Indices for Fine Field-Scale Phenotyping

Abstract

Spectral vegetation indices (VIs) are a well-known and widely used method for crop state estimation. These technologies have great importance for plant state monitoring, especially for agriculture. The main aim is to assess the performance level of the selected VIs calculated from space-borne multispectral imagery and point-based field spectroscopy in application to crop state estimation. The results obtained indicate that space-borne VIs react on phenology. This feature makes it an appropriate data source for monitoring crop development, crop water needs and yield prediction. Field spectrometer VIs were sensitive for estimating pigment concentration and photosynthesis rate. Yet, a hypersensitivity of field spectral measures might lead to a very high variability of the calculated values. The results obtained in the second part of the presented study were reported on crop state estimated by 17 VIs known as sensitive to plant drought. An alternative approach for identification early stress by VIs proposed in this study is Principal Component Analysis (PCA). The results show that PCA has identified the degree of similarity of the different states and together with reference stress states from the control plot clearly estimated stress in the actual irrigated field, which was hard to detect by VIs values only.

1. Introduction

The 2030 Agenda represents an agreement between all 193 UN Member States to introduce a set of common strategies to achieve 17 goals (the Sustainable Development Goals, or SDGs) and 169 targets before the year 2030. SDGs are a collection of global goals to attain a better future. Sustainable agriculture is at the heart of this agenda. This goal is responsible for ensuring food production systems and implementing resilient agricultural methods, influencing the increase in production and productivity, assisting in maintaining ecosystems, adjusting to climate change and extreme weather [1, 2, 3, 4]. Thus, simultaneously taking into consideration improvement of land and soil quality. One of the ways that can help achieve sustainable agriculture is with Precision Agriculture (PA), a method to accurately apply the right treatment, in the right place, at the right time.

Timely detected crop stress allows rapid correspondence and adaptation of planned agricultural activities and preventing negative effects on the yield. Special attention is paid to water stress due to its effect on plant growth and yields [5]. Moving towards PA that stands for concept of managing crop fields considering spatial variation and local field requirements, involves data collection to characterize field spatial variability, mapping, decision-making, and management practice implementation. The growth number of precision agriculture applications has influenced the development of remote sensing technology owing to its ability to conduct higher spatial, spectral and temporal resolutions capabilities and cost-effectiveness. Remote sensing at visible and near-infrared wavelengths (vis–NIR) has been used to formulate many spectral indices for estimating crop properties (e.g. [6]). In irrigation management, crop water state can be estimated by water content or water potential in soil and plants [7]. These parameters can be measured directly and indirectly. Direct methods are gravimetric soil moisture, relative water content [8], energy status of soil water and plant water potential. Indirect estimation crop-water state performed by microwave or radar techniques [9], soil moisture balance calculations [10], and air-plant temperature differences [11]. Indirect methods have many advantages: they do not damage plants and soil structure; do not require laboratory conditions, expensive equipment, complex measurement protocols, and specific technical knowledge. Regardless to the measurements type, the prevalent method for crop stress detection in practice is a comparison of actual crop parameters with reference values for the normal state [12]. Visual agronomic field inspection for assessment wilting, morphometric changes (plant organ shrinkage), and growth rate [13] to estimate crop-water state is still the widespread approach in practice because of its traditional origins, low cost and ease application. The essential drawback of this visual inspection method is ability to detect stress at the obvious visual stage when a crop is already undergoing significant damage. In general any manually assessed plant trait in the traditional way is time-consuming, laborious; introducing errors and sometimes destructive [14, 15]. From 2000s, the era of non-destructive plant-phenotyping platforms, based on image-based techniques, has begun [16, 17, 18].

Spectroscopy can improve current agronomic inspection methods. The physical properties of plants directly influence its reflectivity at different spectral ranges [19]. When crop stress is clearly pronounced in visual spectrum range, it can be detected by traditional field inspection. Before this, stress has already been caused by crop properties and has influenced its reflectivity. Spectral tools provide measurements of plant reflectance with higher sensitivity and in a wider range than the human eye’s capability. Therefore, spectroscopy allows for detecting stress at earlier stages than traditional visual methods [20, 21].

Simple Vegetation Indices (VI) have significantly improved the ability and sensibility of the detection of green vegetation [22]. Vegetation indices are widely used to estimate crop growth status and crop parameters, such as biomass, yield, photosynthesis, Leaf Area Index, etc. [23, 24, 25, 26], and many studies have analyzed the potential of using spectral reflectance indices in wheat starting from the early years of the century (e.g., [27]) and up to the present day (e.g., [28]). High-resolution VIs may detect changes of wheat crop status and it might help to improve crop monitoring [29], nitrogen management [30], and crop yield estimation [31, 32]. Furthermore, it is indeed known that yield prediction while using VIs in wheat can be accurate two months prior to harvest, because yield estimates stabilize and especially during the flowering period, significant correlations between UAV-VIs and yield components were found [33, 34].

VIs estimate plant state by calculating ratios or more complicated mathematical models of reflectance measurements using different spectral wavelengths. Its development began with discovering a strong linear correlation between plant green biomass and the ratio of 2 spectral bands obtained by satellite imagery [35]. The concept of applying spectral data to assess plant parameters was introduced in the 1970s. The VIs era began with the 1972 launch of the first ERTS satellite (Landsat 1) with its MultiSpectral Scanner. The first index was the Normalized Difference Vegetation Index (NDVI) developed as quantitative measurement of vegetation conditions by calculating the ratio between visible and infrared (VIS/IR) spectral bands. Further on, NDVI was applied to assess plant health and estimate other physical properties presenting sufficiently good linear correlation with plant height [36] and asymptotic relationships with Green LAI [37], and has shown its ability to indicate different phenological stages [38]. Further high correlations of NDVI with crop biomass were discovered, Leaf-Area Index (LAI), absorbed photosynthetically active radiation, and canopy photosynthetic capacity [39]. However, NDVI does not have sufficient sensitivity to all crop features [40]. The next developed VIs were sensitive to crop photosynthesis and plant-water state [41, 42], there were several indices sensitive to crop pigmentation; e.g. chlorophyll, carotenoid, anthocyanin [43, 44].

As NDVI includes chlorophyll absorption band, it finds application in estimating chlorophyll and used as health state parameter as well [45, 46]. Likewise, NDVI’s response to physical crop characteristics and chlorophyll content makes it suitable tool to predict yield and detect N deficiency [47], estimate actual evaporation rates [48], assess fraction of Absorbed Photosynthetically Active Radiation [49], indicate soil salinity [50] and calculating crop coefficients for irrigation needs [51]. In addition, there were developed many modified VIs applied under detailed/dedicated conditions for estimating physical parameters and other specific tasks [24, 25, 26, 44]. Other approaches were associated with the ability for collecting spectral data at high spectral resolutions (1 nm). These approaches allowed obtaining changes in plant pigmentation, nutrient content, and chemical properties [32, 52, 53].

Nowadays, there are several dozen VIs and models for different spectral data types, which can be applied for estimating crop state, determining stress, salinity, diseases, and hazardous substances. Irrigation management is one of the most promising directions where VIs could be applied. For this task, there were already developed VIs that are directly sensitive to water absorption; e.g., Water Band index [54] and indirectly by estimating changes in pigmentation and photosynthesis rates [51]. It is important mentioning soil moisture indices such as Temperature Vegetation Dryness Index (TVDI) [55], on account of its potential applicability for estimation water state of agricultural fields. Besides estimation soil moisture TVDI provides information about groundwater depression cone which improves remote monitoring and allows to reduce in-situ measurement [56].

Despite the great importance and potential benefits of these studies, they do not consider the informativeness of indices without co-core physical measurements described above, such as water content, soil water balance and etc. Most of the indices based on water absorption bands require high-resolution spectral data in the short-wave infrared range: Normalized Difference Water Index (NDWI), Moisture Stress Index, Normalized Difference Infrared Index, and Normalized Multiband Drought Index [57, 58, 59, 60]. Spectral data in near-infrared range 900–970 nm provides an opportunity to estimate water content by Water Band Index [61]. Due to the high cost of obtaining required data, the above indices were not widely applied in agricultural management. Consequently, a number of studies were conducted to find a method of applying the current VIs from visible and near-infrared spectral ranges for water stress detection.

Several studies were devoted to discovering a direct correlation between widespread VIs (NDVI, Simple Ratio, Photochemical Reflectance Index, and Structure Insensitive Pigment Index) and leaf water-content, but did not find a significant correlation [62, 63]. Nevertheless, NDVI and Soil Adjusted Vegetation Index (SAVI) found application in irrigation management, because they were proved to be a strong approach for estimating crop coefficients (Kc) and predicting crop evapotranspiration [64, 65, 66, 67, 68]. Further investigation showed that combination between VIs and physical parameters might increase accuracy (R2 = 0.5–0.7) of estimated crop-water state, than single VIs analysis [69, 70]. Accordingly, considering water deficit by VIs, canopy temperature, air temperature, stomatal conductance, or stem water potential in one application might contribute to better detection. The PRI alone showed good correlation with plant water-content (R2 = 0.8) with no additional physical measurements in the model [71, 72]. However, this method also has a disadvantage, as a good PRI correlation can be found only at specific plant stages and when Photosynthetically Active Radiation was above 700 μmol/m2/s.

Estimating plant state only by VIs consists of calculating values and deciphering results. Some of the indices provide a range of optimal values for plant health; e.g., Structure Insensitive Pigment Index, Carotenoid Reflectance Index, Modified Red Edge Simple Ratio, etc. [25, 42, 44]. Other VIs only provide interpretations of physical processes based on increasing or reducing values; e.g., Anthocyanin Reflectance Index, Chlorophyll Absorption Ratio Index, Triangular Vegetation Index, etc. [26, 44, 73]. Since agricultural management requires a common design for a large area, a series of measurements (pixels from spectral images or points by a portable field spectrometer) should be converted to one field parameter, which can lead to loss of essential information. Likewise, the analysis could be complicated since the variations in the VIs from date to date are confounded with crop phenology and, possibly, atmospheric conditions. Consequently, field state analysis by VIs should include calculating values, merging them into a single field parameter, assessment variety, interpretation value according to the vegetation index response, specific environment and plant conditions.

In addition to the complications described above, the fact that the same VIs can be assessed by these two fundamentally different sources of spectral data, namely, space-borne and airborne imagery, and point-based spectral measurements, raises questions regarding its universality and accuracy. Several studies confirm that models based on field spectroscopy need direct comparisons, modification, and validation for application of space-borne and airborne spectral imagery [74]. The imagery data acquires spectral data from the whole scene, including soil, plants, and atmospherics [75]. Field point-based spectral measures could minimize the atmosphere’s effect and reduce signal noise from the soil by simply excluding it from the observed scene (ground footprint). Thus, it is important to remember that the VIs originally developed for one type of spectral data might lose estimation accuracy once applied to another data source without any additional validation or modification.

Sensitivity to plant stress at early stages, estimation of crop parameters, non-destructive method of observations, and existing low-cost equipment for spectral data collection, make VIs the prospective approach for monitoring crop-water state and implementing the optimization of an irrigation schedule. Despite all the advantages and attractiveness, the absence of clearly interpreted VIs’ behavior, relatively low-correlation of individual index with real-field water state, limits spectroscopic application approach in irrigation planning.

The main aim of this study is to assess the performance level of the selected VIs in application to different spectral data sources on actual irrigated agricultural fields without co-core physical control measurements. There were three specific objectives:

1. To determine suitability of different spatial and spectral (e.g. spaceborne multispectral imagery and point-based field spectroscopy) data sources for early detecting stress in the irrigated crops;

2. Define VIs corresponding to early water stress in the irrigated crops.

3. The aim of this study is to identify water stress on actual irrigated field by VIs.

For this propose, three crops are investigated: cotton, tomato and chickpea. The third objective was tested based on chickpea crop, one under actual agricultural management and other field with limited water treatments were studied. Stress identification was implemented using several methods. The first method is interpretation obtained VIs values according description from the original papers and their temporal behavior. The second is validation the correlation between physical field parameters and VIs stated in existing studies. The last is novel approach for estimation early crop stresses by principal component analysis (PCA) proposed in this study. PCA is introduced for processing the matrix of VIs, and further identification uniqueness and similarity of crop states and estimation stress by introduction reference VIs values of crop suffered from dryness.

2. Materials and methods

2.1 Study area

This study was carried out in Kibbutz Hazorea, in northern Israel (32°38′42.0”N 35°07′17.6″E). The fields are located in the typical Mediterranean climate with mild winters and dry summers. The average daily mean annual temperature is approximately 18°C and relative humidity is about 68%. Our observations were conducted during the 2015 summer growing season, from March 1 through August 31. The average daily temperature in spring is 18°C and in summer reaches 25°C. The total precipitation during the study period was 111.1 mm of which 59.9 fell in two days (April 10–11).

The soil is classified as Vertisol according to the USDA Soil Taxonomy [76]. The soil characteristics were as follows: heavy-textured soil, bulk density 1.8–2.0 g cm − 3, pH value 8.0–8.6, organic matter 5–10%, clay 50–60%, and CaCO3 9%. Vertisol is a widespread soil type which is prominent on almost every continent: Africa, India, Australia, southwestern USA (Texas), Uruguay, Paraguay and Argentina [77]. The main feature of this soil type is its high level of clay content (between 50 and 60%). As a consequence, it has a high moisture-holding capacity. Vertisol is known as suitable soil for agriculture, but the level of clay content requires very careful irrigation management [78]. Dry Vertisol soils are conglomerated and cracked. The specificity of the soil in Kibbutz Hazorea required great effort by the farmers: drainage, land reclamation, and finding appropriate agricultural crops. Initially, attempts were made with dry farming. In early years, vegetable gardens and grains were the leading sources of income. Gradually the high watering cost of vegetables melons and fruit trees were changed to more profitable irrigated crops [79].

The following three types of crops were chosen for this study: tomatoes, cotton, and chickpeas. These crops belong to various families of plants: tomatoes are Nightshades, cotton is Malvaceae, and chickpea is Fabaceae. The crops have different sensitivity to water scarcity, the strategy of development, and stress resistance. Chickpea (Cicer arietinum L. ‘Yarden’) was planted on 1st of January 2015 with density 190 kg seeds per hectare. The cover area was nearly 12 hectares. The rows were located at a distance of 90 cm from each other and had East/West direction. The crop was irrigated by a sprinkler on mobile irrigation systems. The irrigation period was 3 weeks since 18th of May. The total water amount was 2400 cubic meters per hectare. Watering was carried out in equal parts daily. Determination of actual phenological stage and assessment plant state by field inspections were conducted based on guideline WATERpak [80]. Tomatoes (Solanum lycopersicum ‘4107’) was planted on open ground on 15th of March 2015 with density 25000 seedlings per hectare. The cover area was nearly 13 hectares. The crop was irrigated by surface drip. Cotton (Gossypium hirsutum L. ‘HA − 195’) was planted on 28th of March 2015 with density 15000 seedlings per hectare. The cover area was nearly 10 hectares. The crop was irrigated by a sprinkler on mobile irrigation systems. Both crops have distance between rows 90 cm and rows had North/South direction. Criteria for estimation actual phenological stage and for assessment crop state by farmers field examination in tomatoes and cotton was provided in the relevant guidelines developed by NaanDanJain Irrigation Ltd. The daily irrigation amount to the crops was calculated by multiplying evapotranspiration to crop coefficient [81]. The obtained value was corrected according to recommendations based on field inspections.

Due to the reported impact of environmental characteristics on correlation strength between VIs and crop-water state, the study was conducted on actual irrigated field. To monitor a crop mainly suffering from water lack, a monitored area preferably should have risks of high temperatures and drought. Thus, crops grown on open fields in Israel are an appropriate object for the proposed study.

Chickpea was grown in open fields (Cicer arietinum L. ‘Yarden’) was chosen for investigation. The monitored chickpea was grown from January 25 until June 28, 2016, on open ground. Planting density was 190 kg seeds per hectare. The planted area was nearly 5.3 hectares. The crop was irrigated by surface drip. The plants have a 90 cm distance between rows, in a North/South configuration. Planting distance on the row was 10 cm. Limit of irrigation water amount for the whole growing season on the chickpea field was 120 mm and actual water consumption amounted 123.62 mm. The potential chickpea yield at the given seeding density declared by the seed producers is 0.55 t/ha (Nir Agricultural Works, Ltd). In practice, a yield of 0.5 t/ha is considered by farmers as good due to the low planting density proposed to prevent fungus development. In 2016, the chickpea yield from the monitored field was 0.47 t/ha.

Irrigation schedule was developed based on potential water consumption [10] using weekly weather forecast and crop coefficients provided by chickpea growth guide [80]. In addition, the chickpea field was equipped with a tensiometer, which was used to determine the metric water potential (soil moisture tension). The field feedback data obtained by tensiometer was applied to correcting and improving the weekly irrigation schedule. The irrigation starts on May 2 and continues for four weeks. The water was distributed unevenly: 26.5, 34, 15.5 and 46 mm per week according to the sequence of the weeks.

For the purposes of the study, control plot of chickpea was planted jointly with main chickpea field on the distance 5 meters. Plot size is 2x2 meter. All activities and treatments on the control plot match to the management on the main monitored field till flowering stage on April 25. Since this date the control plot was not irrigated.

2.2 Data collection

Field trips were conducted from January 2015 through August 2015. During these trips, the crops underwent agronomic inspection and spectral measurements. Crop phenological stages were determined according to crop guides [82]. Agronomic examination confirmed the state of the tomatoes’ and cotton’s normal development and health during the entire growth period. The actual tomato yield was 145 ton per hectare versus 160–170 ton per hectare maximum [82]. The actual cotton yield was 5.2 ton per hectare versus a potential 6.5–7 ton per hectare maximum [81].

The chickpeas developed under limited water treatments. Irrigation started in mid-May 2015 (at yield formation stage) and lasted two weeks. The agronomic inspection noted visual marks of plant dryness at the end of May 2015 [81]. Under-irrigation during vegetative and flowering stages leads to a slowdown in phenological development and a decrease in yield. The actual growth cycle of the chickpeas was 202 days compared to theoretical 150–170 days. The actual yield was 0.45 ton per hectare in contrast to 0.55 ton per hectare expected by farmers.

Additional chickpea tested field was observed during the growth period from February to June 2016. In total, there were 8 observations on the main field. Each campaign consisted of visual crop inspection and health state assessment [13], identifying the phenological stage, height measures, and spectral data-collection. Crop health state was assessed by leaf color, dryness, and plant wilting. The phenological stage was identified according to the description provided by crop growth guides [81]. Crop height was defined as the average of 5 random measures from ground to top crop leaves.

Field spectral reflectance was measured directly on the leaves using a portable field spectrometer (USB4000, Ocean Optics Inc., USA [83]) acquiring data across the VIS and NIR range from 350 to 1100 nm with a resolution of 0.5 nm and an accuracy of 1 nm. The spectrometer was calibrated according to protocol [84] against a white Spectralon plate (Labsphere Inc., North Sutton, NH). The detailed spectra were collected at a nadir view angle approximately 5–10 cm above the crop’s leaf-scale by a bare fiber optic with 25° field-of-view positions. All measurements were carried out between 11:00 to 13:00 o’clock on sunny days when cloud cover was less than 10%. For observation, upper leaves orientated approximately perpendicularly relative to the sun were chosen. During the field trip, about 10 spectra were collected for each crop with 10-spectra repetition for each measurement. The plants were chosen randomly on the same 10x10 meter area. Selecting the control area size is conditioned to ensure uniform water distribution to crops due to homogeneous landscape, which as a complex with other identical parameters (crop and soil type, weather conditions and fertilizers) provides the expectation of the close water treatment for the tested crops. The spectral data from the chickpea field were collected from 5 x 5 meter control patch with a grid layout, and nearly 50–60 spectral signatures were obtained.

The plant water concentration (PWC) was introduced for validation results of agronomic inspections and providing accurate detection water stress in the crops. The plant water concentration is estimated by reflectance Water Index [85]. The ground-based reflectance measurements required for this method were obtained during the field campaigns according to the Penauelas protocol and processed according to Eq.(1).

$$\mathbf{PWC}=\mathbf{684}\ast \left(\mathbf{R}\mathbf{900}/\mathbf{R}\mathbf{970}\right)-\mathbf{620}$$E1

The plant water concentration obtained from the tested samples is performed in Table 1. In addition to visibly health plants tested in this study, the plant suffered from the dryness stress was detected in the cotton field (on 22nd of May) by PWC and introduced to the analysis as a reference of drought.

	Chickpea	Tomatoes	Cotton
22.01	94.17%	—	—
15.03	81.51%	—	—
19.04	89.15%	—	—
22.05	80.45%	87.33%	85.08% (63.85% stressed)
07.06	77.60%	84.17%	83.29%
28.07	—	78.59%	75.93%
31.08	—	—	68.10%

Table 1.

Plant water concentration (PWC) by Penauelas method.

For remote field observation, RapidEye data was kindly provided by the RESA project 597 [86]. The data was presented at preprocessing level 3A [87]. Atmospheric correction was conducted with the generic processing chain CATENA developed at the German Aerospace Center (DLR) [88]. RapidEye images have a spatial resolution of 5x5 meter per pixel and 5 spectral bands (440-510 nm; 520-590 nm; 630–685; 690-730 nm; 760-850 nm). Images were obtained from February 2015 through July 2015, with an average frequency of twice a month. The spatial resolution covered about 5 pixels that matched the study area 10x10 meter chosen for field spectral measures.

Soil moisture tension was measured by tensiometer (“Mottes Tensiometers” Company). The tensiometer consists of plastic tubes filled with distilled water and porous ceramic caps. The tubes were installed in the tomato field at depths of 30, 60, and 90 cm. High tensiometer values means dry soil, while low values indicate high soil moisture. The tensiometer logger recorded data every 20 minutes from April 11 until May 29 (Figure 1).

Figure 1.

Tensiometers records on depth 30, 60 and 90 cm from the chickpea field in 2016.

Meteorological data for this study was provided by the Israel Meteorological Service from Haifa University meteorology station, which is situated 16 km from the observed fields.

Field data for the further analysis consists of weekly irrigation amount records, tensiometers date, measured crop height, and estimated phenological stage. “Irrigation week amount” is an amount of water applied to the field during the week before a field campaign including precipitation. Tensiometer data considered the maximum and minimum recorded values (amplitude) during the last week before a field campaign (“Tens_max” and “Tens_min” respectively), and the actual value at the moment of field spectral measurements “Tens_actual” (at a depth of 60 cm). In this study, plant development was divided into 6 phenological stages (“Stage”): emergence, early vegetative, late vegetative, flowering, fruit formation, and ripening. Field input data for analysis is presented in Table 2.

ID	Date	Stage	Irrigation week amount (mm)	Height (cm)	Tens_max (hPa)	Tens_min (hPa)	Tens_actual (hPa)
C1	19.02	emergence	0	2	—	—	—
C2	16.03	early vegetative/ late vegetative	15	20	—	—	—
C3	03.04	late vegetative	0.4	32	—	—	—
C4	18.04	flowering	21.1	40	19.1	6.8	12.2
C5	06.05	flowering/ fruit formation	26.5	61	20.6	3.4	6.8
C6	11.05	fruit formation	34	66	21.6	6.4	11.3
C7	21.05	fruit formation	15.5	68	44.1	3.9	34.3
C8	30.05	fruit formation/ ripening	54.1	72	32.8	3.9	20.6

Table 2.

Crop parameters and watering conditions on the chickpea field in 2016.

According to tensiometer records, at the week before observation the soil water content on the depth 60 cm turns to the lowest values. It corresponds to the irrigation schedule: at that week the chickpea got the lowest amount of water. Field inspection on May 21 noted a relatively high concentration of yellow leaves not typical for this phenological stage.

Chickpea control plot has no irrigation since April 25. During the month since May 1, the control plot was observed every five days. At this period there was three rainy days from May 23 to May 28 and 8.1 mm fell. The plot monitoring consists of visual inspection, spectral measurements and collection leaf samples. The spectral data was collected according the same protocol like on the main field with a grid layout. Nearly 20 spectral signatures were obtained from the plot at each date. Also, at every observation five representative leaf samples were collected and sent to the lab measurements. In the lab, the leaves were weighted and dried by the oven at 70°C for 24 hours. Leaf water content (LWC) was defined (Eq.(2)) as the difference between initial weight and weight after drying divided to initial weight and converted to the percent [89].

$$\mathbf{LWC}=\left(\mathbf{\text{initial leaf weight}}-\mathbf{\text{weight of dried leaf}}\right)/\left(\mathbf{\text{initial leaf weight}}\right)\times \mathbf{10}$$E2

The dried samples were ground into powder and mixed with 100% acetone at the rate of 0.1 g of dry matter per 10 ml of acetone. The mixtures were centrifuged for 30 min in glass tubes to make the extract fully transparent. The resulting extracts were immediately measured by spectrometer USB4000. Specific absorption coefficients of Chl a and Chl b reported by Lichtenthaler [90] were used to estimated chlorophyll concentration in studied leaves. The results of measured leaves’ parameters are presented in Table 3.

ID	Date	LWC	Chlorophyll content mg/g
S1	May 1	$95.25\pm 2.17\%$	$1.33\pm 0.25$
S2	May 6	$95.12\pm 2.09\%$	$1.01\pm 0.26$
S3	May 11	$95.51\pm 2.14\%$	$1.04\pm 0.39$
S4	May 16	$96.33\pm 2.35\%$	$0.83\pm 0.37$
S5	May 21	$86.27\pm 6.72\%$	$0.76\pm 0.21$
S6	May 26	$69.01\pm 13.41\%$	$0.68\pm 0.16$
S7	May 30	$25.48\pm 13.72\%$	$0.49\pm 0.06$

Table 3.

Leaf water content and chlorophyll concentration estimated in control stressed chickpea leaves.

2.3 Data analysis

VIs for this study was selected according to the following criteria: 1) the considered indices should have application in irrigation management; 2) different types of VIs original response should be presented in this study; 3) required spectral data for VI calculation should match the spectral ranges obtained from the spectrometer and space-borne images. The last one, the considered VI should be developed by different technologies: remote sensing images (satellite and aerial) and field spectrometer. This requirement allows considering the influence of data source on processing results.

Statistical analysis was used to bring all the data to a single format, convenient for interpretation and comparative analysis. This goal was achieved by estimation VI’s average mean and internal variability. Measures of each crop at the same time by the same sensor (spectrometer or satellite) were merged into one dataset. From space-borne images, values of 5 pixels coinciding with the location of ground-based measurements were extracted. For further comparative assessment of variability, the number of field spectral measures was reduced from 10 to 5 by random sampling to match the number of values from spaceborne images.

The heterogeneity of the data (space-borne images and spectral plots) complicates the choice of analytic methods due to their applicability for both data sources. Sophisticated methods such as semivariogram have already proved their effectiveness for geostatistical analysis of remote sensing imagery [91]. However, the specificity of field spectral measures limits the suitable methods for assessment variety of VIs. For the purpose of further comparative analysis, the common statistical method was chosen: a coefficient of variability was applied to the quantification measured variability (Eq.(3)).

$$\mathbf{Cv}=\mathbf{s}/\mid \mathbf{x}̄\mid$$E3

where Cv is a coefficient of variation, s is the standard deviation, and x̄ is the average mean.

Cv values close to zero were excluded from the analysis. Additionally, the temporal variability of each VI was estimated. The temporal variability was defined as the coefficient of variation between average means of datasets during the whole season. The absolute average means of VI also was calculated from average means of datasets during the whole season. The temporal variability and absolute average mean are presented in the result section and named “Total”.

The suitability of application the coefficient of variation was examined by normality Shapiro–Wilk W test in the SPSS environment [92]. The VIs corresponding to the requirements of the normal distribution was rescaled (Eq. (4)).

$$\mathbf{X}\_\mathbf{new}=\left(\mathbf{X}-\mathbf{X}\_\mathbf{min}\right)/\left(\mathbf{X}\_\mathbf{max}-\mathbf{X}\_\mathbf{min}\right)$$E4

where Xnew is a rescaled VIs value, X is the original value, Xmin and Xmax are minimum and maximum permissible VIs values respectively. The minimum and maximum values from “Range of values” in Table 3 are used for rescaling.

The sensitivity limits of VI’s, spatial and temporal variables were assessed by statistical methods for outlining extreme values (outliers). For this examination, the following methods were used: histogram of distribution, percentiles (Tukey’s Hinges), tests for normality (Kolmogorov–Smirnov, Shapiro–Wilk, Normal Q-Q plot), outlier labeling by quarters with g = 1.5 [93] and stem-and-leaf plot [94] to examine VIs spatial and temporal variables on data collected with both portable spectrometer and satellite.

The obtained spectral range allows calculating indices responding to visible and NIR diapason and excludes indices related to SWIR. To find the stress manifestation on different levels (physical and pigment), VIs related to different plant parameters and responded to drought stress are considered: biomass (NDVI, RENDVI, MRNDVI, MRESR), canopy coverage (MCARI2, MTVI, and MTVI2), photosynthesis rate (PRI), chlorophyll (MCARI, TCARI), carotenoid (CRI, PSRI, and SIPI), anthocyanin (ARI), and leaf water content (WBI). Table 4 shows VIs chosen for this study.

VI	Original response	Application for stress detection	Range of values	Reference
NDVI	Green biomass	Drought stress and nutrient uptake	0 to 1	[35]
WBI	Relative water content	Drought stress	0.8 to 1.2	[61]
ARI1	Anthocyanin	Rust infection, drought stress, low-temperature stress	0 to 3	[22]
ARI2	Anthocyanin	Rust infection, drought stress, low-temperature stress	0 to 2
CRI1	Carotenoid	Drought stress, salinity level	1 to 20	[44]
CRI2	Carotenoid	Drought stress, salinity level	1 to 20
PSRI	Carotenoid/ Chlorophyll	Leaf senescence and fruit ripening, salinity stress and water use efficiency	-1 to 1	[43]
PRI	Photosynthesis rate	Water stress	-1 to 1	[41]
SIPI	Carotenoid/ Chlorophyll a	Leaf senescence and unhealth plant, rust infection and aphid	0 to 2	[42]
MCARI	Chlorophyll	Stress detection, nitrogen availability and water stress	0 to 0.7	[95]
MCARI2	LAI	Heavy metals stress, water and nitrogen stresses	0 to 1	[26]
MRENDVI	Chlorophyll	Water and stress	0.2 to 0.8	[24]
MRESR	Chlorophyll	Drought stress	2 to 8	[25]
MTVI	LAI	Drought and salinity stressed	0 to 1	[26]
MTVI2	LAI	Water and nitrogen stresses	0 to 1	[26]
RENDVI	Green biomass	Salinity stress, early stress, stem water potential	0.2 to 0.9	[25]
TCARI	Chlorophyll	Nitrogen stresses	0 to 0.4	[26]

Table 4.

The list of considered vegetation indices and their stress response.

The estimation of VIs applicability for water-stress detection was performed separately for each data set. This estimation is based on a comparative analysis between VI behavior in health and presumably early drought-stressed crop. Tomatoes and cotton were considered as healthy crops. The chickpea lacked water during the month since the last rain (60 mm) on April 11, 2015, until the first irrigation on May 15, 2015. The chickpea state during this period was assumed as water lack suffered, and the absence of visible signs is perceived as the initial stage of (early) stress.

2.4 Regression analysis

Linear regression analysis was applied for estimation VIs’ response to chickpea irrigation. The analysis is designed to detect linear correlation between physical parameters (stage, height, weekly irrigation amounts, and tensiometer data) and calculated VIs averages. Input chickpea states for this analysis are C4-C8 (Table 2) when tensiometers measurements were obtained. The strength of correlation was assessed by Spearman and Pearson coefficients and two-tailed significance, since among the parameters there are both: scale and ordinal. Correlation analyses were performed in the SPSS environment.

2.5 Principle component analysis

The proposed method for estimation crop water stress by VIs is principal component analysis (PCA) implemented in MatLab environment [96]. This approach is a multivariate technique that analyzes a matrix of numerous intercorrelated quantitative dependent variables. PCA extracts the dominant patterns to new orthogonal variables called principal components and plotting the variables in new multi-demotion, where each dimension presents estimated principle component [97, 98]. The purpose of applying this method is consideration VIs values combination as single characteristic for estimating a pattern in crop behavior under stress. The VIs values from control chickpea plot were used as reference crops stress to the identification pattern corresponded to crop-water lack. The analysis was conducted on 3 input datasets: VIs averages with their variability, VIs averages and variability separately. For this analysis VIs values were rescaled by Eq. (4).

The first test of PCA determines each crop state as a point in multidimensional space, where every dimension is an individual parameter estimated from VIs averages and variability and plotting it in 2-dimensional space by determining 2 principal components. In this test the similarity of crop state assessed by mutual arrangement of points on the plot. The second PCA test was implemented in the reverse: combinations of indices and variability were plotted by 2 principal components and crop states became variables and were presented on the biplot as vectors showing the contribution to these components. The manner of vector presentation calls a correlation circle: if vectors with the same length are close to each other, they are significantly positively correlated (R close to 1); if they are orthogonal, they are not correlated (R close to 0); if these vectors are on opposite sides, they are negatively correlated (R close to −1). When the variables are close to the center, some information is carried on other axes, and then any interpretation might be hazardous.

3. Results

The results of VI variability assessment showed strong differences between space-borne imagery and portable field spectrometer data. Consequently, this section is subdivided by the data sources. For each data source, the VIs and their variability, sensitivity limits of variability were considered, and comparative analyses between healthy and stressed crops were performed.

3.1 Space-borne imagery data-source

The chickpea field was observed by space-borne imagery 11 times from February 1 until July 15, 2015, with a frequency of every 2–3 weeks. The first three observations in February and March characterized relatively low values of NDVI: 0.239, 0.325 and 0.323. Since NDVI ranges from 0.05 on bare soil to 0.90 on dense vegetation, the observed low values corresponded to the low density of crop seedlings. On May 4, NDVI considerable increased to 0.717, which corresponds to crop growth. Then, NDVI declined to a value of 0.493 on April 20, without obvious reasons. All three observations in May recorded high NDVI values (0.861, 0.840 and 0.835). The last three observations during the summer coincided with the crop ripening and show an NDVI decline: 0.606, 0.328 and 0.329.

The majority of the examined indices (NDVI, EVI, Green Atmospherically Resistant Index (GARI), Green Normalized Difference Vegetation Index (GNDVI), LAIwp, LAIc, Pigment Specific Simple Ratio (PSSR), TVI, MCARI2, Modified Triangular Vegetation Index (MTVI), MTVI2 and Modified Chlorophyll Absorption Ratio Index (MCARI)) display similar behavior, presenting low values during the first three observations; a significant increase on April 4; an inexplicable decline for the next observation; high values in May, and a decline during the summer. Slight differences in the behavior of these indices could be detected only with the first three observations. Some of the indices gradually increased during this period; e.g., EVI, GARI, LAIwp, TVI, MCARI2, MTVI, and MTVI2. An additional group had approximately equal values for the second and third observation: NDVI, LAIc, PSSR, and MCARI. Only the GNDVI had equal values for the first and third observations and a higher value in the second. Anyhow, these differences are insignificant in comparison with changes during the season.

Two indices related to carotenoid and anthocyanins display unique behavior in chickpeas. The Plant Senescence Reflectance Index (PSRI) behaves in an inverse manner with respect to health indices. At the beginning of the season, it had high values: 0.442 and 0.474. In April, it declined to 0.192, and then rose to 0.221. During May, the PSRI turned to its lowest values: 0.111, 0.061, and 0.155. In summer, it gradually increased: 0.327, 0.436, and 0.532. Despite the PSRI, a behavior pattern or correlation with the crop’s physical development in Anthocyanin Reflectance Index (ARI) was not found. At first observation, ARI was 4.617, the value then rose to 6.249. At a third observation, it dropped to 2.854. In April, ARI behavior matched the health indices: on April 4, it rose to 7.956 and on April 20, it declined to 3.034. During May, the ARI consistently gave values of 8.824, 6.907, and 7.914. During the summer, ARI also behaved without a pattern and displayed values of 4.979, 3.715, and 6.122. During the whole observation period, ARI values were approximately in the same range.

Observations of tomato fields were conducted from April 4 until July 16, 2015. In total, 9 surveys were conducted with a frequency of 2–3 weeks. NDVI gradually increased from the beginning of observation until June 4: 0.201, 0.209, 0.465, 0.648, and 0.810. This behavior is a response to crop growth and an increase in green biomass. On June 9, a small decline of NDVI (0.777) was detected. Nevertheless, on the June 22 observation, the value rose and started to decline as related to the crop’s ripening (0.789, 0.597, and 0.585).

Similar to the chickpea field, the majority of indices calculated for tomatoes displayed close behavior: NDVI, EVI, GAVI, LAIc, PSSR, TVI, MCARI2, MTVI, MTVI2, and MCARI. These indices showed increments since June 4, a recession on June 9, a small rise in value at the next observation, and further decline until the end of the growing session. The recession on June 9, pronounced with different strength in these VIs. In some of them the decline was significantly pronounced: EVI (0.751–0.635-0.738), LAIc (1.734–1.422-1.528), TVI (31.649–26.600-32.476), MTVI (0.769–0.641-0.791), MCARI (0.307–0.158-0.344), MCARI2 (0.696–0.604-0.689) and MTVI2 (0.696–0.604-0.689). In other indices, the differences were barely noticeable: NDVI (0.810–0.777-0.789), GAVI (0.634–0.562-0.584), and PSSR (814.036–619.157-678.157).

On April 9, some indices displayed behavior close to NDVI, but without a decline. The GNDVI and LAIwp rose from April through June, and then gradually declined during the summer. Except for an increment in value at the last observation, ARI also matched this group of parabolic behavior. The PSRI in tomatoes did not have a behavior pattern. During the season, 4 cycles of increment and decline were observed, with the values ranging from 0.163 to 0.327.

The cotton field is located next to the tomato field, thus, its observations were carried out on the same dates: 9 surveys with a 2–3 week frequency from April 4 until July 16, 2015. The study was conducted from cotton seedlings until flowering. This explains NDVI behavior: low values in April–May (0.205, 0.219, 0.212, and 0.266), an increase in June (0.498, 0.570 and 0.858), and consistently high values during July (0.785 and 0.871). As in the other crops, NDVI behavior is repeated in other indices: EVI, GARI, GNDVI, LAIwp, LAIc, PSSR, TVI, MCARI2, MTVI, and MTVI2. On May 6 and July 2, NDVI showed a slight decline. The same reaction on crop state was pronounced in other indices with various intensities. LAIc, PSSR, MCARI2, and MTVI2, repeat NDVI’s behavior and displayed a decline on both dates. LAIwp, TVI, and MTVI declined only on May 6. GARI and GNDVI declined only on July 2.

All three pigment indices developed under laboratory conditions showed unanticipated behavior in cotton. MCARI grew unevenly during the two months from April 4 until June 4: 0.011, 0.014, 0.016, 0.016, and 0.082. The next 3 measurements during the summer show a decrease: 0.054, 0.047, and 0.028. The last measurement unexpectedly increased to 0.035. The PSRI also increased at the beginning of the season: 0.273, 0.302, 0.398, and 0.499. Since summer, PSRI behavior had not displayed a trend: 0.3, 0.393, 0.116, 0.108, and 0.182. Also, ARI did not exhibit a behavior trend during the whole season. The index growth changed to a regression and then reversed: 3.012, 2.844, 3.419, 4.898, 4.739, 5.878, 9.463, 6.889, and 11.286.

Abnormal behavior can be observed in chickpea: on April 20, increased VIs that are linked to carotenoid and sharply decreased VIs that are linked to greenness and physical parameters. During May, the VI recovers and returns back to expected behavior. A similar situation was observed in tomatoes during June 2015 and in cotton during July 2015. Unlike chickpeas, the VI decline in tomatoes and cotton was slight and not pronounced in all the VIs. Among the declined indices in tomatoes and cotton, there were NDVI, GARI, LAIc, PSSR, MCARI, MCARI2, and MTVI2.

The range of VI variabilities calculated from satellite images is from 0.3 to 34%. It was not possible to calculate variability for GARI because of mixed (negative and positive) values in datasets: the mean average was close to zero. A strong correlation between values or phenological stages and variability was not found.

The normality tests show the nonnormal distribution of measurement variability from the satellite data source. The outlines were determined in two ways. Steam-and-Leaf Plot defined extreme values above 16.1%. Labeling by quarters defined the limits of the acceptable range from −5.52% to 15.64%. According to this outlining, LAIwp, LAIc, and GARI should be excluded because of their high variability. Despite MCARI2; MTVI, and MTVI2 also showed outline variability. They will await further consideration because the extreme variability occurs only once – at the late ripening stage of chickpeas.

The temporal variability of VI as calculated by satellite images ranges from 22.95% to 137.82%. In the study case, the high temporal variability indicates a significant change in VIs during the growth period. Statistical methods did not define outlines. The smallest temporal variability (less 40%) was found in GNDVI, ARI, and PSRI. The highest variability (more than 100%) was observed in GARI, LAIc, PSSR, and MCARI.

Pattern behavior of all considered VI (except ARI) in healthy crops is linked to phenological development (as has already been described earlier). Both healthy crops had a small decline in greenness and physical parameters, and an increase in carotenoid during the flowering period (in cotton) and yield formation (in tomatoes). The chickpea’s VIs were also linked to phenological development and had a behavior similar to healthy crops. An exception is an abnormal behavior during April 2015 at the flowering stage. Despite healthy crop behavior in the middle of the growth period, a decline in chickpea greenness VIs was evident and pronounced in all VIs.

3.2 Spectrometer data-source

The chickpea field was examined by spectrometer during the growth period 5 times from the end of January until the beginning of June with an average frequency of once a month. NDVI calculated from spectrometer shows high values and increments during the season: 0.671, 0.740, 0.743 and 0.763. In this growth trend, the decline to 0.654 on May 22 seems abnormal. Similarly, high values and increases during the season with an abnormal decline in the end of May were shown by LAIc, PSSR, MCARI2, and MTVI2. Most of the health indices duplicate NDVI behavior but show a slight decline at the last observation: EVI, GARI, GNDVI, LAIwp, TVI, and MTVI.

ARI, PSRI, and MCARI, developed in laboratory conditions, display a common behavior trend in chickpeas: high value at the beginning, a decline during March, and an increase during April and May, and then a decrease at the end of observation in June.

The tomato field was observed 3 times during the irrigation period: May 22, June 7, and July 28. Analogous to the previous crop, spectrometer measurements recorded high values for the indices during crop growth: NDVI (0.905 and 0.915), EVI (0.797 and 0.983), GARI (0.764 and 0.770), GNDVI (0.758 and 0.761), LAIwp (5.005 and 5.251), LAIc (2.802 and 2.936), PSSR (2501 and 4367), TVI (27.445 and 37.120), MCARI2 (0.845 and 0.913), MTVI (0.735 and 1.002), and MTVI2 (0.845 and 0.913).

All indices, except those developed in laboratory conditions, veer to the lowest values for the last measurements: NDVI (0.822), EVI (0.680), GARI (0.557), GNDVI (0.615), LAIwp (2.830), LAIc (1.849), PSSR (1794), TVI (24.368), MCARI2 (0.723), MTVI (0.682), MTVI2 (0.723). This could be explained by the fact that the tomatoes were at their ripening stage. Thus, the difference in indices behavior is when the values were higher: on vegetative stage at May 22, or June 7. The NDVI, LAIc, and PSSR rise during this period. The EVI, LAIwp, TVI, MCARI2, MTVI, and MTVI2 declined from the vegetative to the flowering stage. GARI and GNDVI remained constant. Laboratory indices behave differently from each other. MCARI increased during the whole period of observation: 0.200, 0.221, and 0.324. The ARI displayed an increase and then decrease (1.210, 2.246, and 2.135); where the first value is significantly lower than the next two. PSRI values were close to zero: −0.003, −0.002, 0.010.

During the irrigation period, the cotton field was examined 4 times, monthly from May through August. Additionally, during the first field trip, a few plants were found with visual signs of leaf dehydration. These plants were measured separately from healthy plants. Thus, consideration of indices for cotton includes behavior studies and comparative analyses.

NDVI in cotton showed growth from May until the end of July (0.757, 0.815, and 0.837) and a decrease at the fruit formation stage at the end of August (0.780). As in the previous analysis, most of the indices match the NDVI trend. GARI, GNDVI, LAIwp, LAIc, PSSR, MCARI2, and MTVI2 also increased from May to July and declined at the end of August. However, some of the health indices did not react to crop changes at the fruit formation and growth stage during the whole observation period: EVI, TVI, and MTVI.

In laboratory indices, there was no trend corresponding to crop development. The ARI reached its highest value (0.549) when first measured, and then during the growth period values remained approximately the same: 0.301, 0.326, and 0.323. The MCARI behaved in a wavelike manner: 0.170, 0.247, 0.141 and 0.204. PSRI values were too low for detection of changes: 0.006, −0.009, −0.003 and − 0.001.

Comparison of cotton health and stress on May 22 shows a reaction of all the indices on visually detected leaf dehydration. Health indices developed in the field conditions have values in stressed crops that were 1.5–2 times lower than in healthy ones: NDVI (0.537–0.757), EVI 0.545–0.916), GARI (0.350–0.587), GNDVI (0.464–0.601), LAIwp (1.111–3.044), LAIc (0.542–1.207), PSSR (306.455–528.177), TVI (35.787–24.459), MCARI2 (0.488–0.746), MTVI (0.640–0.945), and MTVI2 (0.488–0.746). Indices developed under laboratory conditions in contrast with field conditions, have higher values in stressed crops: ARI (2.823–0.549), MCARI (0.170–0.305) and PSRI (0.151–0.006).

Variability in spectrometer measurements is also rather different from satellite images. The common range of variance in spectrometers’ VIs is from 1.9% to 127.94%. PSRI for all crops and ARI for chickpeas have negative values. Therefore, estimating their variability is not possible.

The normality tests show the non-normal distribution of variability from the spectrometer. The outliers were determined by two methods. Steam-and-Leaf Plot defined extreme values after 50.28%. Labeling by quarters defined the limits of the acceptable range from −18.30% to 52.26%. Outlining variables belong to GARI (1 value), LAIwp (1 value), LAIc (3 values), PSSR (3 values), and MCARI (3 values). All these indices had extreme variability in chickpeas on May 22, 2015.The temporal variability of VI, calculated by the spectrometer, ranging from 3.84% to 46.01%. Both statistical methods, Stem-and-Leaf Plot and labeling by quarters determined only one extreme value – 46.01%. According to outlining, the optimal range of temporal variability is from 3.84% to 33.29%, and PSSR had an extreme value.

VIs calculated by field spectrometry could be grouped by behavior trends. The first group of variability showed a gradual increase in cotton and tomatoes during yield formation and ripening (Figure 2a). NDVI, LAIc, PSSR, MCARI2 (and MTVI2 – had the same values as MCARI2) belong to this group. In chickpeas, these indices displayed common behavior as well, but differ from cotton and tomatoes: high value at the establishment, lowest values at vegetative stage, a slight increase at flowering and yield formation, and an abnormally high value on May 22, 2015. In MCARI2 and MTVI2, the increase was expressed slightly for this date. The next group includes LAIwp, GARI, and GNDVI (Figure 2b). These indices behaved in tomatoes in a manner similar to the previous group. In cotton and chickpeas, it also matched behavior in the previous group, except for relatively high values at the late vegetative stage. The third group consists of EVI, TVI, and MTVI (Figure 2c). These VIs have a common unique type of behavior for each crop type. MCARI did not match any group (Figure 2d).

Figure 2.

Variability trends of VIs by spectrometer in different crop types, a is the first behavior group (NDVI, LAIc, PSSR, MCARI2), b is the second behavior group (GARI, GNDVI, LAIwp), c is the third behavior group (EVI, TVI, MTVI),d is indices not falling into the trends (MCARI). The x-axis labels: e – establishment, v – vegetative stage, f – flowering, y – yield formation, r – ripening.

Summarizing the above, VIs calculated by a spectrometer had no common behavior pattern for all crops. For chickpeas, all VIs related to chlorophyll, green biomass, and LAI have comparable behavior during growth season: increasing since establishment until flowering, and declining after yield formation until ripening. What is unusual is that the values of the indices for May 2015 were lower than those for June 2015. In tomatoes and cotton, VIs related to health can be divided into two groups according to their behavior. The first group includes NDVI, GARI, GNDVI, LAIpw, LAIc, PSSR, MCARI2, and MTVI2. These indices have a direct correlation with crop development. The second group consists of indices from EVI, TVI, and MTVI. These indices show abnormal behavior for healthy crops: in tomatoes, the decrease begins at the late vegetative stage, while cotton has low VI values during the growth period and then begins to rise after the flowering stage. Indices developed in laboratory conditions (ARI, MCARI, and PSRI) show unique behavior for each crop and there is no possibility to distinguish any trends.

Behavior patterns for healthy crops (cotton and tomatoes) could be clearly defined in VIs related to phenological development: NDVI, GARI, GNDVI, LAIpw, LAIc, PSSR, MCARI2, and MTVI2. VIs grew until flowering and decreased during yield formation and ripening. In chickpeas, the same VI behavior was observed except at season end. At yield formation during May 2015, this decrease was pronounced in the VIs related to phenology. Regardless of expected continuing decline at yield formation in June 2015, all the indices increased. Thus, the decline in May 2015 during the drought period could be considered as a reaction to stress. In other indices, behavior patterns for healthy crops were not defined and consequently cannot reveal abnormal behavior in chickpea stress.

3.3 Validation

The first results show that health crops have an approximately common value of VIs obtained by the field spectroscopy. In the same time, results give the right to assume that crops under stress that situated in the same area have a high spatial variability of indices. Anyhow, a limited dataset of field spectral data for comparative analysis is not suitable for the approval of variability as a stress indicator. Thus, to verify the hypothesis about indication stress by variability, all field spectral measures obtained in the study was merged in “health” and “stress” groups. Control “stress” group consists of 18 point spectral measurements from the stressed cotton on 22.05.15. Control “health” groups are presented by 40 values of cotton (10 measures X 4 dates) and 20 (10 measures X 2 dates) of tomatoes field spectral measurements collected during the vegetative and fruit formation stages. Tested groups consist of chickpea field spectral measures. The “health” tested group consists of 20 (10 measures X 2 dates) spectral measures from the vegetative and flowering stage in chickpea. The proposed “stress” group is presented by 10 spectrums on 22.05.15 when the variability in the first results was the highest. For this test, indices with higher (more than 40%) variability detected in the first analysis were applied: GARI, GNDVI, LAIwp, LAIc, and PSSR.

The test shows (Figure 3), that variability of GARI, GNDVI, LAIwp and LAIc in the group of “health” cops merged from several dates is lower than in the “stress” dataset obtained on the one date. Also, the results display that level of variety of the control “stress” cotton and tasted “stress” chickpea are approximately the same and rather far from the “health” groups (especially in GARI, LAIwp, and LAIc). Only PSSR shows high variability in all groups. This fact was already noticed at the first stage of the analysis. The absence of a relationship with the stress in crops and high variability of the values indicates the non-occurrence of the use of the PSSR algorithm for field spectroscopy.

Figure 3.

Variability of GARI, GNDVI, LAIwp, LAIc and PSSR obtained by the spectrometer in “health” and “stress” crop groups order as follow: “Health” cotton, “Health” tomato, “Health” chickpea, “Stress” cotton, “Stress” chickpea, from dark gray to light gray.

3.4 Early water stress detection by VIs

The first stage of analysis is calculation VIs averages and consideration their temporal behavior (Figure 4). The W value of the Shapiro–Wilk test is significant for all considered VIs. The VIs informativeness assessment for stress estimation based on their response. All graphs are presented in the range of allowable values from the original papers or in the range of optimal values for green vegetation where it was possible (CRIs, PRI, PSRI, SIPI).

Figure 4.

VIs’ changes over time on the chickpea field in 2016: a) anthocyanin indices; b) carotenoid indices; c) chlorophyll indices; d) canopy cover indices; e) NDVI group; f) MRESR; g) WBI; h) SIPI; i) PRI; j) PSRI.

Interpretation of VIs responded to anthocyanin and carotenoids (Figure 2a,b) is ambiguous. On one side, in the stress conditions, the concentration of anthocyanin and carotenoids should increase. This principle lays in interpretation of ARI and CRI behavior from original papers: weakening vegetation contains higher concentration of these pigments, thus increase indices values indicates stress. However, the content of these pigments also relates to chlorophyll level [99]. Subsequently, decline of anthocyanin and carotenoids can be consequences of chlorophyll reduction, which is also a sign of stress. This contradiction complicates the interpretation of Anthocyanin and Carotenoid Reflectance Indices. The highest values of ARIs (Figure 2a) were obtained at the beginning and at the end of season. The lowest values of this index correspond to the chickpea vegetative period, after which the anthocyanin concentration began to grow smoothly and slowed down on May 21 only, when chickpea get least watering. The high values of ARI at the beginning and end of observation corresponds to chickpea extreme stages: emergence and ripening. Detected sensitivity of ARI to critical stages coincides with the results of previous studies of this pigment [100, 101]. The decline on May 21 is abnormal and can be reaction on stress, but the existing information is not enough for an unambiguous interpretation. Like ARIs, CRIs have smooth graph with high values at the edges. But unlike anthocyanin, the level of carotenoids on May 21 increase.

Indices developed for estimating green LAI and biomass, chlorophyll content and health (Figure 4b–f) have common temporal trend: the lowest values at emergence, increase and highest values on April 3, and notable decline on May 21. At other dates these indices have stable values and did not manifest any crop changes. Even deviation of these indices is almost inessential concerning the declared common range for green vegetation because of data obtaining on leaf scale [102], it is enough to identify chickpea at late vegetative stage as most stable and health and detect signs of weakness on May 21. However, the absence of a reference range for leaf scale spectral measurements does not allow to assess the significance of weakness and define it as stress. PRI (Figure 4g) related to photosynthesis rate match the trend of health indices behavior, but with less intensity of expression. On the smooth graph of PRI, the most pronounced is decline on May 21. Based on these graphs it can be concluded that chickpea experienced weakness and decreased activity at the week with the least watering.

PSRI (Figure 4h) shows carotenoid/chlorophyll ratio and originally proposed for estimation canopy senescence. Increase of the index value is interpreted as stress. In this study the index behaves opposite to indices of green LAI (Figure 4d–f), except May 21 when did not increased as expected. Thus, PSRI does not confirm stress estimated by previous group of health indices. Last two graphs on the Figure (Figure 4i,j) presents relative water content by WBI and carotenoid/chlorophyll ratio estimated by SIPI. Both graphs represent straight lines that do not reveal any significant changes during the time. On closer examination it can be seen that SIPI is similar to PSRI but with less expressed values.

As a result of the VIs temporal behavior analysis, one can assume the chickpea stress on May 21. This assumption is confirmed only by some of the examined indices and it is difficult to determine the degree of stress significance. The assumption of stress was made based on VIs temporal changes and information of irrigation amount. The absolute values of the indices on May 21 belongs to range of optimal values for green vegetation.

3.5 VIs variability

Variability is a complicating factor in interpreting VIs behavior and identifying stress. Since at one observation there was a measured series of spectral signatures, the VIs can vary throughout the field. Thus, consideration of VIs as the average of the dataset, as at the previous stage of analysis, should be combined with a VIs’ variability assessment. The high level of variety (more than 50%) makes the average of VIs uninformative for estimating field crop state. As was mentioned above, the results of the Shapiro–Wilk test is significant for all applied VIs that allow to consider their variability. The results of calculation variability are presented in the Table 5 and marked according to intensity.

Table 5.

The variability (Cv) of the considered VIs on the tomato field in 2016:

  – less than 10%;

  – from 10% to 25%;

  – from 25% to 50%;

  – more than 50%.

The highest level of variability was detected in indices related to anthocyanin and carotenoid. Also, MCARI has variability more than 50% since chickpea flowering. Other indices were stable and has variability near or les 25%. WBI and SIPI that did show temporal deviation have the least variability from 0.4 to 5.7%.

To date, the role of the variability of leaf-scale spectral measurements for crop state estimation has been poorly studied, but already found that high variability indicates different stresses in leaves [103, 104]. In the next section of PCA, the VIs averages will be studied commonly and separately with their variability to identify its effect on crop state estimation.

3.6 Regression analysis

Considered in this study VIs did not reveal linear correlation with irrigation amount and crop height. CRIs, PRI, MRENDVI, TCARI show strong relationship with maximal drought obtained by tensiometers (R2 = 0.75–0.95) and the same indices together with MRESR and RENDVI fit the actual tensiometer data with accuracy R2 = 0.81–0.96). At the same time no relation stronger than R2 = 0.59 was found with minimum values obtained by tensiometers. Also, several indices reacted to the phenological stage. Nevertheless, these relationships were weak (R2 = 0.5–0.7). In this case of study, only 6 of 17 considered indices reacts on soil water content on chickpea field. These results highlight, the complexity of application VIs as estimator of physical parameters. Such approach requires validation for each specific crop and field to identify the most suitable and informative VIs. Also, incomplete correlation effects on crop state estimation accuracy that can reduces the possibility of detecting water stress at early stages.

3.7 Principal component analysis

The advantage of implementation PCA is the ability to find unique behavior patterns by combination of indices, despite their relatively close values in the irrigated chickpea that limits such approaches like interpretation temporal behavior and regression analysis. The field parameters were excluded from this analysis and crop water stress is identified by spectral data only. The detection of stress is based on comparison of the similarity considered and reference stressed crops. Despite the close values of the indices, their combination turns to the unique signature that describes the singularity of crop state at each campaign.

The variability and regression analysis results were used to define the input VIs for PCA: the indices with variability more than 50% were excluded (ARIs, CRIs and MCARI) as non-informative. The first runs of PCA were failed due to significant differences in crop states at specific development stages (Figure 5). The initial chickpea stage C1 (Table 2), and late stress stages from the control site S6 and S7 (Table 3) were plotted by PCA far from other states and make deviation of other states insignificant and hard identifiable. The second reason of closeness the plotted states in the first results of PCA is slight VIs deviation because of too wide range of values used for normalization.

Figure 5.

Principal component analysis for the chickpea states from the field and control plot: a) first PCA test for VIs; b) second PCA test for VIs.

Based on the experience of the first results, new configurations for PCA have been developed. The extreme chickpea states (C1, C2, S6, S7) when watering is inappropriate were excluded from the analysis. To increase the sensitivity to changes in crop states, the considered VIs were normalized (Eq. (4)) by minimum and maximum values obtained in these states. After these changes, the analysis became more accurate and allowed to distinguish slight variation in chickpea (Figure 5). In the first PCA test, the crop states are plotting according to estimated 2 principal components makes possible to estimate their similarity based on their closeness in space (Figure 6a,c,e). In the second PCA test, the crop state is characterized by vectors (Figure 6b,d,f). The similar crop states have the same length and angle of the vector in the second PCA test. In the all performed tests, the first and second principal components together explained 70–80% of the variance.

Figure 6.

Principal component analysis for the chickpea states C3-C8 and S1-S5: a) first PCA test for VIs; b) second PCA test for VIs; c) first PCA test for VIs’ variability; d) second PCA test for VIs’ variability; e) first PCA test for VIs and their variability; f) second PCA test for VIs and their variability.

In the first tests, C3, C8 and S5 are most distant from the plot center and can be called extreme states. C3 is situated at the opposite part of plot from C8 and S5 and has opposite value of the first principal component. States C8 and S5 are located on the same plot side by the first principle component, but at a considerable distance from each other. While states C7 and S4 have almost same location on the plots. States C4-C6 and S1-S3 lay in the middle of the plot. Among the first PCA tests, variability and VIs with their variability plotted the chickpea states such way that it can be seen the dynamics of state changes. At the beginning of the experiment the state of the chickpea on the control plot is close to the chickpea on the field (Figure 6c,e). During first weeks of the experiment, when stress was not significant, the states C5, C6, S2 and S3 still has similar location, but they still located on the different sides of the second principle component axis.

The second PCA tests were less informative. The vectors of chickpea states were distributed in several groups and have different length that limits interpretation and identification correlations. The most graphic is the plot of second PCA test for VIs and their variability (Figure 6f). All chickpea states are distributed between two opposite conditions: active growth at vegetative stage C3 and suffered from drought control chickpea S5. Stressed and ripening states C7, C8, S4, S5 have negative values of the first component, while other states have positive.

The results of PCA tests indicates similarity of chickpea state on May 21 with crop suffered from drought. The assumption of stress in chickpea field on May 21 corresponds to the low soil moisture at week before the spectral measurement. The implemented PCA test indicates the early water stress in chickpea by spectral data only, while other considered in this study traditional methods of VIs analysis could not uniquely detect lack of water in chickpea field. The obtained results confirm already discovered PCA’s applicability for precise analysis in agricultural management [105] like estimation nitrogen concentration [106] or detecting diseases [107].

4. Discussion

The discussion section presents a recommendation for applying VI’s for determining water stress in irrigated crops. Differences in the variance of intensity, temporal changes, and correlation to phenological development were mainly dependent on the data source. VIs from satellite data sources are strongly related to phenological stages and significantly changes with time. This correlation is caused by the common measuring of vegetation and soil background. Since crop growth and leaf development, the noise from soil is reducing and greenness rises. Thus, satellite imagery collects not only spectral plant characteristics but also records physical development. On the one hand, it is a good approach to estimate crop physical parameters: LAI and biomass. On the other hand, it reduces sensitivity to changes in crop pigmentation.

In the satellite data source, the variability was estimated by few image pixels (due to the field size). Therefore, the variance in measurements could be interpreted as differences in crop characteristics on adjacent 5 x 5 m areas. All indices showed low variability (less 20%). This can be explained by crops conditions that should not be significant different across the field. Nevertheless, the crop characteristics cannot be absolutely the same and variability shows up in satellite’s sensitivity to the differences. The lowest variability was observed in indices developed for satellite data sources. The indices developed by field and laboratory spectroscopy, and modified for airborne detection provide high variance regardless the proposed response (physical or pigmentation). Variability in TVI with its modifications reacts on crop type. In chickpeas, these indices have lower variance than in tomatoes and cotton.

To unequivocally estimate water stress in chickpeas by VIs using satellite images is complicated. The abnormal VI behavior for April 2015, in chickpeas, could be considered a reaction to water stress. There were already 9 days of drought since the last intensive precipitation (60 mm). The recommended level of irrigation at this growth stage is nearly 5 mm per day [80]. It is logical to assume that water stress would result and reach a maximum at the beginning of the irrigation period, in May 2015. However, VIs calculated from May 1, 2015, toMay 16, 2015, did not detect stress. Moreover, all VIs testified to the maximally a healthy and strong chickpea state on May 1, 2015, and May 6, 2015. Thus, abnormal VI behavior on April 20, 2015, is not related to water crop stress.

Measurements by field spectrometer have typically high variability and slightly pronounced changes over time. Despite satellite data, direct leaf measurements reduced soil background noise and minimized the influence of physical plant parameters. It can explain weekly correlation with phenological development and lower temporal variability than in VIs calculated by satellite images. Measurements on leaf-scale by high-resolution spectral equipment provided pigmentation sensitivity.

Crop pigments are a suitable stress indicator. Stressors affecting plants primarily lead to impaired or perturbed metabolism and photosynthesis [108]. Chlorophyll level is the basic parameter of photosynthesis. Thus, an abnormal chlorophyll level points to stress. Apart from metabolism, pigments play a protective role against stress. The increase of carotenoid and anthocyanin indicate crop weakness [22, 44]. That is why pigments are considered to be the optimal approach for early, in vivo, detection of abiotic plant stress [86].

Particular attention should be paid to the variance of measurements during one field trip; as stress leads to changes in pigment production, the adjustment is not evenly implemented. Therefore, the increase in pigmentation variability could be a sign of stress. The highest values of VI variances was detected at emergence, ripening stage, and at water deficit in chickpeas on May 22, 2015. Emergence, ripening, and water lack could be considered as stress periods. At emergence stage, the crop is a weak, unformed plant that is practically non-resistant to external conditions. At ripening stage, the plant stops supplying nutrients to leaves and begins dying. Prolonged water deficiency leads to crop drought and a limitation of nutrients for leaves. The unstable state of the plant due to these stresses leads to uneven leaf nourishment resulting in pigmentation variability. Likewise, high variability of point-by-point field spectral measures can be an indicator of unstable and stressed crop state.

Almost all VIs, apart from the highest variability during emergency, ripening, and stress, additionally showed increases in variance at other times. These defined a few behavior trends in VIs variance. The presence of trends in variability may mean that common increase of VI variances was a reaction to some unstable or stressed crop state that was not noticeable during visual inspection. This assumption gives grounds to study variability as a parameter for stress estimation during early stages in irrigated crops. Hence, the high temporal frequency of field spectral measurement could clarify variability trend in different phenological stages under normal conditions that allow accurate and timely highlight fluctuations associated with stress.

VIs calculated from spectrometers were more sensitive to stress than VI from satellites. Most of VI calculated by spectrometer displayed an abnormal behavior in chickpeas on May 22, 2015, when there was an expected lack of water due to a long drought. The reaction to stress was estimated only in VIs related to phenological development. VIs with high sensitivity to pigmentation displayed complex behavior during the growth period when a behavior pattern could not be observed in healthy crops requiring further stress detection by comparative analysis.

Determining lack of water in irrigated crops is a complicated task. The goal is to define stress before water lack causes problems in physical development. Consequently, stress should be detected at an early stage when it only affects metabolism and related pigmentation. Thereafter, VIs linked to crop pigmentation are expected to be more applicable than VIs related to LAI and green biomass. However, this study has shown that method and scale of spectral data collection from the irrigated crops significantly influence VI correlation with various parameters.

VIs calculated from satellite images are strongly related to physical crop development regardless of proposed physical or pigmentation sensitivity. Type of inputs applicated to develop VIs (satellite and airborne images, field and lab spectroscopy) not effect on VI sensitivity as was expected. Likewise, satellite spectral images could be considered as a useful source for estimating crop physical parameters and other tasks related to it, such as crop classification. The low variance of VIs on the same date and high temporal variability allow for clearly estimating the crop state for the whole field and observe crop development. At the same time, these reasons lead to low-pigmentation sensitivity and make satellite images unsuitable for detecting water stress in irrigated crops at early stages.

Point-by-point field spectroscopy measurements of the irrigated crops provide results with high sensitivity to crop pigmentation. Measuring reflectance with an accuracy of 1–2 nm allows identifying and quantifying leaf pigments [109, 110]. However, successful early stress detection depends not only on the degree of sensors‘sensitivity to pigmentation but also the approach for estimating the crop state. Among the considered VIs, a few (NDVI, GARI, GNDVI, LAIpw, LAIc, PSSR, MCARI2, and MTVI2) were found exhibiting a correlation with phenological development. This correlation was weaker than in VIs calculated from satellite images. However, it was enough to define pattern behavior by crop development stages in healthy crops and detect abnormal behavior in drought-stricken chickpeas. The accurate estimating crop stress by VIs responded to the leaf pigmentation requires a comprehensive analysis that includes the specifics of the development strategy to each crop species and the metabolism features. Without these specific parameters, there could be identified only global trends of change in pigment composition that limits the opportunity to detect stress at the early stage.

Another feature of spectrometer measurement is hypersensitivity that leads to high variability of values on the relatively small and homogenous area. The variance of some VIs reached 127% High variability casts doubt on the measurements’ validity. Outlining extreme values defined the optimal range of variability in this study - from 18.30% to 52.26%. Such high variability rates in field spectrometer measurements can significantly reduce accuracy and clarity of crop state estimation. Since the variability can reach extreme values, this parameter should be tested when the agriculture field is monitored by VIs from point spectral measures. Despite to this shortcoming, the variability itself can provide information about plant state: the study has shown that variability has certain behavior trend; e.g., high variability rates were associated with periods of crop stress (GARI, LAIwp, and LAIc). Since the level of variability is responded to uneven changes in pigment production, the increase of Vis variations indicates alteration of behavior strategy: the transition of phenological stage or reaction to uncomfortable environmental conditions (i.e. stress). Thus, VIs variability could be a suitable approach to detecting stress at the early stages.

In practical agriculture, crops should be in optimal conditions, since stresses due to imperfections in the management can lead to lower yield. Spectral measurements on leaf scale allows for tracking numerous slight changes in the plant associated with both early stresses and natural processes characteristic of each phenological stage [111]. On the one hand, this level of sensitivity is necessary to detect early stress in irrigated crops [112, 113]. On the other hand, it requires an approach for isolating external stress from natural processes. Also, the downside of high-precision measurements is high variability of spectra, that doubts the results [114], which was confirmed in this study. However, the fact that variability of leaf scale spectral measurements is also sign of stress forces to leave this parameter in the analysis. For this purpose, it is necessary to consider a dataset of parameters (VIs averages and their variability) that characterize different plant processes. The advantage of PCA is ability to process commonly the VIs and their variability and results of this study proves that inclusion of variability improves accuracy of stress estimation by PCA together with reference stress states.

5. Conclusion

To date, there are many scientific works claiming the advantages of using spectral indices in agricultural needs, including the detection of stress. Unfortunately, the results of this study discover: the fact that the index showed a response to stress is not enough to clear stress detection by VIs in practice. The stress detection task is complicated by the fact that an agricultural field requires a unified management approach and hence one parameter value of crop state. It means that the dataset of obtained spectral data (point measures by field spectrometer or a number of pixels from spectral imagery) should be transformed into a single unified value, which will be used in further decision-making about field treatments. Exactly this requirement has been met at the first stage of the analysis: calculation of the VIs average for each dataset. Meanwhile, results of the variability estimation indicate that field spectroscopy could provide high spatial variability that depreciates the significance of the VIs average itself. The intensity of variation obtained in this study was in the range from 0.4–177% on the relatively small test area (5 x 5 meter).

The next stage of analysis was dedicated to identifying in practice response of VIs to field parameters related to irrigation. 6 from 17 considered VIs show strong correlation with soil moisture on the chickpea field, while other 11 VIs also proposed for water state estimation did not show significant correlation. Based on these results it can be argued, monitoring crop state by VIs based on their correlation with physical parameters need to be supported by validation from each specific crop and field to define the informative VIs in the specific case. Such requirement can significantly reduce the attractiveness of this approach for practical applications.

The novel method considered in this paper for early stress detection by VIs, was developed based on the experience that the joint consideration of parameters is more effective than a separate one, which was stated in the “Introduction.” The proposed consideration of common VIs and their variability was implemented through PCA. The results met expectations: the joint examination of the parameters turned out to be much more precise in determining crop state than the individual arguments in the first stages of analysis. Proposed in this study VIs and their variability processing by PCA in common with spectral library of reference crop states allow detect stress on early stages. An additional advantage of PCA for estimating crop state is the presentation of results as a grid of vectors, where extreme positions define opposite crop states (active growth and death), while other states within this grid can be estimated from the proximity to one or the other extreme state.

To conclude, field spectroscopy is a promising technology for detecting water stress in irrigated crops at an early stage. Nevertheless, existing remote spectral methods should be supported by additional physically measured ground-based data. VIs as a separate independent approach for crop monitoring is limited by the difficulty of unambiguous interpreting. Hypersensitivity of spectral measurements could provide early stress detection on the one hand, and on the other hand, it requires a modified approach for accurate analysis and interpretation. Determining crop stress by spectroscopy as an independent method should be implemented through common consideration of VI behavior patterns in healthy crops, the intensity of correlation with physical development, and rate of measurement variance on a particular date.

Acknowledgments

The authors want to thank the German Aerospace Center (DLR), RESA Science Team, Neustrelitz for the support by providing the satellite data of the RapidEye Science Archive (proposal no. 597) through funding from the German Federal Ministry for Economic Affairs and Energy. We also thank the staff of the DLR (Oberpfaffenhofen), namely Thomas Krauss and Peter Fischer, for the atmospheric correction using the generic processing chain CATENA.

This work was supported by the Israel Water Authority [Grant #41683, 2014], and by Israel Water Authority (grant number 41683); and the Israel Innovation Authority, Israel-Italy Joint Innovation Program for Industrial, Scientific and Technological Cooperation (grant number 60283).