Identifying rice stress on a regional scale from multi-temporal satellite images using a Bayesian method

Identifying rice stress on a regional scale from multi-temporal satellite images using a Bayesian method

Meiling Liu, Tiejun Wang, Andrew K. Skidmore, Xiangnan Liu, Mengmeng Li PII: S0269-7491(18)33981-2

DOI: https://doi.org/10.1016/j.envpol.2019.01.024 Reference: ENPO 12063

To appear in: Environmental Pollution

Received Date: 31 August 2018 Revised Date: 5 November 2018 Accepted Date: 8 January 2019

Please cite this article as: Liu, M., Wang, T., Skidmore, A.K., Liu, X., Li, M., Identifying rice stress on a regional scale from multi-temporal satellite images using a Bayesian method, Environmental Pollution (2019), doi: https://doi.org/10.1016/j.envpol.2019.01.024.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

M

AN

US

CR

IP

T

AC

CE

PT

ED

M

AN

US

CR

IP

T

AC

CE

PT

ED

Identifying rice stress on a regional scale from multi-temporal satellite images 1

using a Bayesian method 2

Meiling Liu1,*, Tiejun Wang2, Andrew K. Skidmore 2,3, Xiangnan Liu1, Mengmeng Li2 3

1

School of Information Engineering, China University of Geosciences , Beijing 100083, 4

China 5

2

Faculty of Geo-Information Science and Earth Observation (ITC), University of 6

Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands 7

3

Department of Environmental Science, Macquarie University, NSW 2109, Australia 8

*Corresponding author. liuml@cugb.edu.cn, Tel: +86 10 82321060, Fax: +86 10 9

82322095 10

Abstract: 11

Crops are prone to various types of stress, such as caused by heavy metals, drought and 12

pest/disease, during their life cycle. Heavy metal stress in crops poses a serious threat to crop 13

quality and human health. However, differentiating between heavy metal and non-heavy 14

metal stress presents a great challenge, since responses to environmental stress in crops are 15

complex and uncertain, with different stressors possibly triggering similar canopy reflectance 16

responses. This study aims to infer the occurrence probability of heavy metal stress (i.e., Cd 17

stress) on a regional scale by integrating satellite-derived vegetation index and 18

spatio-temporal characteristics of different stressors with a Bayesian method. The study area 19

is located in the Hunan Province, China. Seven scenes of Sentinel-2 satellite images from 20

2016 and 2017 were collected, as well as Cd concentrations in the soil. First, the probability 21

of rice being stressed was screened using the normalized difference red-edge index (NDRE) 22

M

AN

US

CR

IP

T

AC

CE

PT

ED

at all the growth stages of rice. Further, the stressed rice was used as input, along with the 23

coefficients of spatio-temporal variation (CSTV) derived from NDRE, for a Bayesian method 24

to infer rice exposed to Cd pollution. The results demonstrated that NDRE was a sensitive 25

indicator for assessing stress levels in rice crops. The CSTV with a threshold of 2.7 26

successfully detected rice under Cd as well as abrupt stress on a regional scale. A high map 27

accuracy for Cd induced stress in rice was achieved with an accuracy of 81.57%. This study 28

suggests that vegetation index obtained from satellite images can assist in capturing crop 29

stress, and that the used Bayesian method can be very useful for distinguishing a specific 30

stressor in crops by incorporating temporal-spatial characteristic of different stressors in 31

crops into satellite-derived vegetation index . 32

Key words: Heavy metal stress; Sentinel-2 images; Bayesian method; Coefficients of 33 spatio-temporal variation. 34 35 Acknowledgments 36

This research was financially supported by the National Natural Science Foundation of China 37

(No.41601473 and No.41371407) and the China Scholarship Council (No.201706405037) 38

and co-funded by the Faculty of Geo-Information Science and Earth Observation (ITC), 39

University of Twente, the Netherlands. 40

41

1. Introduction 42

Rice (Oryza sativa L.) is among the most important crops and provides a staple food 43

source for nearly half of the world’s population (FAO, 2008). Rice production has made a 44

M

AN

US

CR

IP

T

AC

CE

PT

ED

significant impact on the lives of people worldwide. However, rice often grows in 45

unfavourable environmental conditions, exposed to, for instance, heavy metals in the soil, 46

drought, pests, diseases, and mismanagement, thus considerably reducing production 47

(Lichtenthaler, 2010). Excessive accumulation of heavy metals in agricultural soils may not 48

only influence crop yield, but also affect food quality and safety in a negative way, posing a 49

serious health threat to both animals and humans (Srivastava et al., 2017). As living standards 50

improve, the demand for high quality rice grains is increasing across the world (Zhang et al., 51

2005). Consequently, accurate detection of heavy metal stress levels in rice is essential for 52

food safety and living standard improvement. 53

Remote sensing has proved a valid method for detecting plant stress. For example, there 54

are studies examining plant response to heavy metal-contaminated soil in abandoned mines 55

and different field conditions (Ren et al., 2009; Chi et al., 2009; Liu et al., 2010, 2011a, 56

2011b; Wang et al., 2012; Liu et al., 2012; Wang et al., 2018), and other research 57

endeavouring to establish a relationship between pest/disease, nitrogen stress and reflectance 58

characteristic（Lemaire et al., 2008;Gilles et al.,2008; Sankaran et al., 2011; Zhang et al., 59

2012 ;Yuan et al., 2017). It is noted that to date experiments involving stress in crops were 60

performed with knowledge about the environmental conditions, such as heavy metal 61

contamination or pest/ disease. The fact that most studies have been conducted on a local 62

scale rather than on a large scale, forms one of the primary stumbling blocks for applying 63

remote sensing to precision detection of crop stress, primarily because crops growing in 64

ecosystems on larger scales may be associated with various stressors, including heavy metals, 65

diseases, insects, and drought. 66

M

AN

US

CR

IP

T

AC

CE

PT

ED

It is recognized that a multitude of stressors with different modes of action may trigger 67

similar physiological responses. For example, stress may induce: (i) changes in canopy 68

architecture and their internal structure, due both to a decrease of leaf area production and to 69

a destroyed cellular structure; (ii) pigments contents decrease due to a limitation of the 70

synthesis of chloroplasts; (iii) canopy temperature increase due to lower transpiration fluxes 71

(Chen and Kao, 1995; Larsson et al., 1998; Phadikar et al., 2012; Yuan et al., 2017). As a 72

result, the reflectance spectra of crops may exhibit a similar variation in the visible light and 73

near-infrared bands in a certain time period when crops are exposed to heavy metals, pests, 74

disease or nutrient stress (Yang et al., 2007; Liu et al., 2010;Mee et al., 2017). Therefore, it is 75

challenging to accurately identify stress factors in crops over large areas using a single-date 76

remote sensing imagery (West et al., 2003; Zhang et al., 2012). 77

However, different stressors in crops do cause variation in spatio-temporal 78

characteristics. Crops are regularly affected by pests, disease and nutrient stress, while crops 79

growing adjacent to industrial and mining areas may be subject to continuous stress exposure 80

to heavy metals. Heavy metals in the soil are viewed as nonpoint source pollution with high 81

toxicity, persistency. They can impede crop growth during the whole growth cycle. Heavy 82

metal-induced stress in rice crops can be characterized as spatio-temporally stable and 83

continuous (Scudiero et al., 2015, Tian et al., 2017). Crops stressed by pest and/or disease are 84

not evenly distributed across a region, as they also depend on management as well as physical 85

barriers, such as mountains and rivers (Phadikar et al., 2012; Donatelli et al., 2017). Thus, 86

pest and disease stress tends to occur sporadically across an area and is typically spatially and 87

temporally transient across all stages of growth (Liu et al., 2012; Scudiero et al., 2015). 88

M

AN

US

CR

IP

T

AC

CE

PT

ED

Nutrient stress affects rice in a more persistent and stable pattern during the different growth 89

stages than stress caused by pests and diseases does (Mee et al., 2017). However, there is 90

clear evidence that areas affected by nutrient stress vary from year to year, as well as show 91

spatial heterogeneity, as rice growers adjust basic fertilization rates using the fertilization 92

status of the previous year, thus trying to avoid rice to be stressed by a nutrient surplus or 93

deficiency in consecutive years (Lemaire et al., 2008). In other words, heavy metal stress in 94

crops may be defined as relatively stable stress across space and over long time scales (e.g., a 95

decade), as well as persistent throughout all crop growth stages. This in contrast to 96

pests/disease and nutrient stress in crops, which may be defined as abrupt stress that is 97

typically transient in space and time. 98

To accurately detect rice stress on a regional scale in agro-ecosystems, spatio-temporal 99

characteristics should be considered as important prior information for discriminating heavy 100

metal stress from pests/disease and nutrient stress. One particular technique that holds 101

considerable promise in combining prior information and observation data is the Bayes' 102

method (Pearl, 1988). Compared with frequentist methods, Bayesian methods can make 103

better use of available information, and so typically produce good results. Bayes' theorem 104

allows us to combine a prior belief about the probability of an event (e.g., temporal-spatial 105

characteristic of different stressors in crop) with observation data (e.g., spectral data), 106

resulting in a new and more robust posterior probability distribution of an event. Bayesian 107

methods are also effective for quantifying uncertainty by calculating probability (Jensen et al., 108

1990; Hosack et al., 2008). Bayesian methods have been used extensively as tools for 109

understanding complex relationships in ecology and environmental sciences (Borsuk et al., 110

M

AN

US

CR

IP

T

AC

CE

PT

ED

2004; Holloman et al., 2004; Marcot et al., 2006; Dorner et al., 2007; Park and Stenstrom, 111

2008; Dlamini, 2010; Liedloff and Smith, 2010; Voie et al., 2010; Aguilera et al., 2011). 112

More recently, Bayesian methods have been shown to be useful in detecting crop disease 113

severity (Schikora, et al., 2010; Bauer et al., 2011; Hernandez-Rabadan et al., 2014; Raza et 114

al., 2014). However, the capability of Bayesian methods in mapping and monitoring different 115

stressors in crops has remained unexplored. 116

The objectives of this research are to (i) identify stressed rice and unstressed rice on a 117

regional scale using multi-temporal satellite-derived vegetation index, and (ii) map the 118

probability distribution of heavy metal stress in rice on the basis of temporal-spatial 119

characteristics of different stressors using a Bayesian method. 120

2. Study area and data 121

2.1. Study area 122

The study area is located in Hunan province, China,covering an area of 1175 km2 from 123

27°35'N to 28°N and 113°1'E to 113°17'E (Fig. 1), and it has one of the highest agricultural 124

yields in the province. The rice (type: Boyou 9083) used in the single-rice cropping region in 125

Hunan Province, is normally transplanted in early June and harvested in late September. The 126

predominant soil type in the study area is red soil (Orthic Acrisol, FAO-UNESCO system). 127

The top 20 cm of soil boast abundant fertility with figures for mean soil organic matter, total 128

nitrogen, total phosphorus, and total potassium of 30.43 g/kg, 1.79 g/kg, 0.82 g/kg, and 0.56 129

g/kg, respectively (Hu et al., 2013). The climate is subtropical monsoon with mean annual 130

temperatures of approximately 16 to 18 °C and a mean annual rainfall of approximately 1500 131

mm. Rice cultivation in this region can rely on sufficient irrigation water. Therefore, the rice 132

M

AN

US

CR

IP

T

AC

CE

PT

ED

is hardly affected by water stress, based on meteorological data collected from 2016 and 2017 133

(Table S1). 134

The study area is an important industrial zone. There are many factories producing 135

pollution, such as printing plants, paper mills, and clothing factories. The rice is irrigated with 136

water from the Xiang Jiang River, which contains industrial wastewater discharges (Zhao et 137

al., 2015; Wan et al., 2017). The concentration of Cd in paddy soils is reported to generally be 138

higher than Level III Soil Environmental Quality Standards in China (the standard threshold 139

value for Cd concentration in soil is 1 mg·kg-1), which is harmful to normal plant growth and 140

causes potential health risks to consumers (Zhang et al., 2009; Norse and Ju, 2015; Chen et 141

al., 2016). This Cd, which cannot be artificially controlled, is the main stressor, continuously 142

interfering with the normal growth of rice during the entire growth period. Additionally, the 143

rice planted in the study area is subjected to improper amounts of fertilizers and farm 144

chemicals or mismanaged, causing extra stress (e.g., pest, disease and nutrient stress), 145

according to field investigations and acquired pest and disease information obtained from the 146

database provided by the Plant Protection and Inspection Service, Hunan Province, China 147

(http://www.hnagri.gov.cn) (Table S2). 148

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig. 1. Location of the study area in China and the distribution of the seven study sites 151

within the study region 152

2.2. Ground data 153

Soil Cd concentrations were downloaded from both the Environmental Monitoring 154

Centre of China and the Environmental Protection Agency in the Hunan Province, China. We 155

also carried out soil sampling from 2014 to 2016 as part of the experiments. Topsoil samples 156

(depth 0–20 cm) were collected in this study area. At each of the seven study sites, the soil of 157

four subplots was mixed. The metal content in the soil was analyzed by the Chinese Academy 158

of Agricultural Sciences, Beijing, China. Soil samples were also analyzed for 159

diethylene-triamine-penta-acetic acid-extractable heavy metals using the method of Lindsay 160

and Norvell (1978) and the AAS to determine metal concentrations. 41 soil Cd concentrations 161

were used, based on data collected from the Environmental Monitoring Center of China and 162

our soil samples, to evaluate whether it would be effective to employ the proposed model for 163

M

AN

US

CR

IP

T

AC

CE

PT

ED

assessing Cd stress in rice on a regional scale. 164

To establish the stress assessment model, in 2017, we conducted a field experiment at 165

seven chosen sites (Fig.1 (b)). Rice at Site 1 was unstressed and growing healthily. Rice at 166

Site 2 was found to be affected by pests and disease, and showed common visual symptoms 167

such as a change in color and morphology of the rice, as caused by the yellow stem borer, the 168

rice leaf folder, the rice bug, or by leaf blast, brown spot and variation in “vertical spread and 169

horizontal spread”. Rice in Sites 3-7 was exposed to Cd contamination. Detailed information 170

on the seven sites is presented in Table1. 171

Table 1 The description of seven sites used to establish a stress assessment model in the 172 study area 173 Experiment area Central Geographical location Average Cd concentration (mg.kg-1) Description Rice pixel number Training Pixel number Testing Pixel number Site 1 27°40'3"N,113°10'3"E 0.86 Unstressed 1173 784 389

Site 2 27°37'16"N,113°14'49"E 0.88

Pests and disease stress

1062 708 354

Site 3 27°40'40"N,113°8'1"E 3.51 Cd stress 644 429 215 Site 4 27°47'0"N,113°9'59"E 1.40 Cd stress 246 164 82 Site 5 27°49'43"N,113°2'22"E 7.81 Cd stress 50 33 17 Site 6 27°58'59"N,113°2'28"E 5.09 Cd stress 891 594 297 Site 7 27°59'0"N,113°5'59"E 5.88 Cd stress 808 539 269 174

2.3. Sentinel-2 satellite image and processing 175

Sentinel-2 satellite images provide a versatile set of 13 spectral bands spanning from the 176

visible and near infrared to the shortwave infrared. The Sentinel-2 image can be freely 177

downloaded from the European Space Agency (ESA) Sentinels Scientific Data Hub 178

(https://scihub.copernicus.eu). We downloaded seven Sentinel-2 images (with processing 179

M

AN

US

CR

IP

T

AC

CE

PT

ED

level 1C) acquired during the rice growing seasons of 2016 and 2017, namely on 28 August, 180

2016, as well as on 17 July, 24 July, 6 August, 21 August, 15 September and 30 September, 181

2017. 182

The processing level 1C includes radiometric and geometric corrections with sub-pixel 183

accuracy (ESA, 2015). We applied atmospheric correction to the Sentinel-2 images using 184

Sen2cor version 2.2.1 within the Sentinel-2 Toolbox (S2TBX), on the Sentinel Application 185

Platform (SNAP) version 4.0.2. In this study, the two bands, B05 (centered at 705 nm) from 186

698-713 nm and B07 (centered at 783 nm) from 773-793 nm, were used with the same spatial 187

resolution of 20 m. To avoid misjudgments arising from an abrupt change in the consecutive 188

periods, cloud and cloud shadows for each image were masked. 189

3. Methods 190

To identifying rice corps under Cd stress on a regional scale, the following methods 191

were implemented (see Fig. 2). (i) A supervised classification method was used to map rice 192

fields. (ii) Vegetation indices (VIs) derived from multi-temporal satellite images were applied 193

to map the unstressed rice and stressed rice. (iii) A Bayesian method was employed to detect 194

rice crops under Cd stress based on the stressed rice and the coefficients of spatio-temporal 195

variation of VIs. 196

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig.2. Flow chart of identifying rice corps under Cd stress on a regional scale 198

3.1. Mapping rice fields 199

A single-date, cloud-free Sentinel-2 image acquired on 16 September 2017, was used to 200

map rice distribution. To obtain the spatial distribution of rice in this study area, we 201

classified five land cover types (i.e., rice, forest, water, urban or bare land, and other) using a 202

maximum likelihood classifier, for which the training datasets were obtained from both 203

visually interpreting Google Earth imagery and field surveys. Total training and testing 204

sample pixels for the five land cover types were 318413 and 17096, respectively. And the 205

selected rice training and testing datasets amounted to 85585 and 4246 sample pixels, 206

respectively. We assessed the classification accuracy of the paddy rice using both overall 207

M

AN

US

CR

IP

T

AC

CE

PT

ED

accuracy and kappa coefficient. 208

3.2. Mapping stressed and unstressed rice crops 209

Based on our previous study, normalized red edge differences (NDRE) in bands B05 and 210

B07 proved to be the most effective indicators for monitoring stress levels on a regional scale 211

(Liu et al., 2018). Here, multi-temporal NDRE was used as the sensitive vegetation index (VI) 212

to determine the stressed and unstressed rice (Fig.2 (b). The unstressed and stressed rice 213

pixels from Site 2 and all other experimental Sites were selected, respectively (Table 1). In 214

addition, the selected rice pixels were divided into training and testing datasets, amounting to 215

3251 and 1623 sample pixels, respectively. A support vector machine (SVM) classifier was 216

used to distinguish between unstressed and stressed rice. We applied a kernel method and a 217

set of labelled pixels from unstressed and stressed rice to train a SVM classifier with a radial 218

basis function kernel to produce a probabilistic output (Lin et al., 2007). To evaluate the 219

accuracy of the stressed rice classification, a testing dataset was used to calculate a confusion 220

matrix, with the probability threshold set to 0.5. 221

3.3. Mapping rice crops under Cd stress 222

3.3.1. Calculating spatio-temporal characteristics of different stressors 223

The coefficients of spatio-temporal variation (CSTV) are referred to as the degree of 224

spatial variation for VI (i.e., NDRE) over a given time (i.e., rice growth stage). CSTV was 225

calculated to determine the spatio-temporal characteristics of different stressors (Liu et al., 226

2018). A greater CSTV represents a local instability that is significantly different from values 227

in the spatial neighborhood over the given period. In this study, the CSTV was calculated 228

using a threshold value to differentiate 'stable stress' from 'abrupt stress'. 229

M

AN

US

CR

IP

T

AC

CE

PT

ED

CSTV is computed based on the period-to-period variability of a given index time series 230

(Q_i ), Q_i and CSTVare expressed as follows: 231 1, , 1, / , i m i m i m i m i Q =V₊ −V or V₊ V , (1) 232 i Q Q CSTV δ − = , (2) 233

where

V

V

and m+1, respectively. Q_i , Q , and

δ

grid-averaged VI and standard deviation of Q along a given index, respectively. 236

From Eq.2, the CSTV is normalized by a mean and a standard deviation of one (

δ

is mandatory to know how to calculate the mean and standard deviation value of CSTV on 238

the basis of a specific range of spatial scale (i.e., several pixels or the whole region). When 239

monitoring Cd stress or abrupt stress on a large-scale, we need to know the suitable spatial 240

scale to capture the characteristic of different stressors accurately. To quantify the impact of 241

spatial scale on VI, the uniform grid was divided up for the whole region, with the number of 242

grid cells ranging from the original value of (22) to (2n). Moran I is a method for measuring 243

spatial auto correlation and can be used to select the optimal scale (Overmars et al., 2003; 244

Espindola et al., 2006; Johnson and Xie, 2011). Aiming to find the most suitable spatial scale, 245

the Moran’s I of VI in each grid was extracted and compared from different spatial scales. 246

The calculation of Moran’s I for each grid is expressed as: 247 1 1 2 1 1 1 ( ) ( ) ( ( ) ) n n ij i j i j n n n i ij i i j n w X X X X I X X w = = = = = − − = −

∑ ∑

∑

∑ ∑

Where n is the total number of grids, X_iandX are the VI value of the i _j th and j th pixel, 249

M

AN

US

CR

IP

T

AC

CE

PT

ED

respectively, X is the mean VI value of each grid. Each weight w is a measure of the ij

250

spatial adjacency of X_iand X_j . If X_iand X_j are adjacent: wij =1, otherwise wij =0.

251

Further, the mean, minimum, maximum, standard deviation of the Moran’s VI in each grid 252

were calculated at the respective spatial scale (q= 22, 24, 26…and 2n). Some researches 253

indicated that the break points of the scale parameter can be used as the optimal spatial scale 254

(Dra˘gu t et al., 2010, 2014). In this study, the lowest standard deviation of Moran’s I for each 255

grid was selected as the most suitable spatial scale for identifying the different type of 256

stressor. 257

3.3.2. Building a Bayesian network 258

In this study, a Bayesian method was employed to map rice crops under Cd stress. A 259

Bayesian network (BN) can be viewed as a probabilistic graphical model that encodes joint 260

probability distribution over a set of random variables (Nielsen and Jensen, 2007). It is used 261

for knowledge representation and reasoning about a data domain (Cheng and Wang, 2010; 262

Thirumuruganathan and Huber, 2011; Skidmore, 1989). BNs are based on the Bayes theorem 263

on conditional probability between two events A and B. The probability of A given that B 264

occurs p(A/B) is given by: 265 ( / ) ( ) ( / ) ( ) p B A p A p A B p B × = (4) 266

A BN structures a set of data (usually a finite set of random samples containing various 267

variables) and analyses their relationships (Aguilera et al., 2011). In this study, the structure 268

of the BN is based upon expert knowledge of the identification of rice under Cd stress from 269

multi-temporal Sentinel-2 images. It integrates VI and the spatio-temporal characteristic of 270

rice under different stressors. Fig. 2(c) presents a schematic overview of the proposed 271

M

AN

US

CR

IP

T

AC

CE

PT

ED

Bayesian network model. We assume that we have a class variable Ci, i.e., Cd stress and 272

abrupt stress. The nodeS denotes stressed rice, and the node HM n

S denotes Cd stress in rice. 273

The nodes NDRE1, NDRE2,…NDREm denotes NDRE value in different growth stages.

274

Following Bayes' theorem, 275 ( ). ( ) ( ) ( ) HM HM n n HM n p Ci p S Ci p Ci S p S = (5) 276 Where ( HM) n

p S is the same for all possible states of states of Ci , p Ci and ( ) ( HM ) n

p S Ci

277

are the prior probability and condition probability, respectively. 278

The p Ci ( ) in Eq. (5) models the prior distribution of different stressors in rice, in this

279

step, detected abnormal variation (i.e., abrupt stress) is excluded from Cd stress according to 280

a sequence of time series points. We want to avoid extreme probabilities, which would imply 281

absolute certainty about either Cd stress or abrupt stress. Therefore, CSTVs greater and less 282

than the threshold value were set to 0.05 and 0.95, respectively, for the prior probability by 283

introducing a block weighting function. 284

The condition probability p S( nHM Ci) in Eq. (5) is inferred by the iterative Bayesian 285

method, which is often applied to detect change characteristics based on each newly added 286

observation. Details may be found in the literature (Skidmore, 1989; Skidmore and Ryan, 287

1990; Reiche et al., 2015, 2018). In other words, p S( nHM Ci)from the first growth stage to the 288

nth growth stage, previous observation (n-1), the current observation (n), as well as upcoming 289

observations (n+1), can be stated as follows: 290 1 1 ( _nHM _i) ( _nHM _n HM). ( _n HM _i) p S C = p S S ₋ p S₋ C (6) 291 With ( HM ) n i

p S C being the posterior, denoting the conditional probability of Cd stress in 292

rice at_{n S}( HM) _given HM n

S as the new evidence of the nth future observation that can have a 293

M

AN

US

CR

IP

T

AC

CE

PT

ED

p S refers to the total probability of the signal Sn. In other words, the

294

posterior probability of Cd stress in rice in the (n-1)th observation will be taken as the 295

conditional probability of Cd stress in rice in the nth observation. p S( _nHM S_n−1HM) is the prior 296

probability, using CSTV calculated based on field observation and our published paper (Liu 297

et al., 2018). 298

For the initial step (n =2), ( 1 )

HM i p S C is calculated as follows: 299 1 1 1 2 ( n HM i) ( HM i) ( , ,... m)

p S₋ C = p S C ∝p S NDRE NDRE NDRE (7)

300

In Eq. (7), p S NDRE NDRE( 1, 2,...NDREm) denotes the probability of stressed rice. Here, 301

the probability of stressed rice is calculated by the SVM classifier with a radial basis function 302

kernel (See part 3.2). 303

3.3.3. Assessing accuracy of the classified rice under Cd stress 304

The Cd in soil and the probability of Cd stress in rice were compared to assess the 305

accuracy for mapping rice crops under Cd stress. According to the Soil Environmental 306

Quality Standards (SEQSs) in China (GB 15618-1995), Level III is referred to as the 307

threshold value of soil heavy metal concentration for ensuring crops and forest production 308

and normal growth of plants. For Cd, the threshold value soil of Level III SEQSs in China is 309

1mg·kg-1, therefore, a safe level (<=1 mg·kg-1.) and a polluted level (>1mg·kg-1) for Cd in 310

paddy soil were classified. In this study, Cd stress map only uses pixels with a probability (p) 311

greater than the probability threshold (p>0.5). Of a total of 41 soil sampling points, 38 312

sampling points were used and three sampling points were removed, because of their 313

geolocation in a cloud area. This results in four groups, determined by two factors (Cd and p): 314

(p>0.5, Cd>1), (p>0.5, Cd<=1), (p<=0.5, Cd>1) and (p<=0.5, Cd<=1). The accuracy of Cd 315

M

AN

US

CR

IP

T

AC

CE

PT

ED

stress assessment on a regional scale is then given by: 316

(

)

(

)

4.1. Classification of rice fields 319

The spatial distribution of rice crops was obtained by applying a maximum likelihood 320

classification method (Fig.3). The overall accuracy was 96.99%, and the kappa coefficient 321

was 0.95. The user’s accuracy and the producers’ accuracy of rice were 94.85% and 96.30%, 322

respectively. According to the average area of each field in Southern China, and four 323

neighbors pixel retains rice type. After classification, we applied the four-neighbor operation 324

rule to eliminate single pixels of rice within a cluster of no-rice pixels. 325

To avoid high rates of false detection of abrupt stress caused by clouds and cloud 326

shadows, the clouds and cloud shadows in planted rice regions were screened for all periods. 327

These pixels covered by clouds or cloud shadows (approximate 10% of all rice pixels) were 328

excluded before extracting rice stress information in the next step. 329

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig.3. Spatial distribution of rice in the study area (green color denotes rice). 331

4.2. Classification of unstressed rice and stressed rice 332

The unstressed and stressed rice were obtained by the SVM classifier. Fig. 4 shows the 333

probability results of rice stress in this study. The accuracy of the stressed rice classification 334

was used to calculate based on a confusion matrix, with the probability threshold set to 0.5 335

(see Table 2). This Table 2 shows an overall accuracy of 96.05%, and a corresponding kappa 336

coefficient of 0.90. In Fig.4 can be seen that most of the rice area was affected by stress 337

factors, and it was estimated that the pixels with a probability greater than the probability 338

threshold (p > 0.5) account for 82.64% of all the rice pixels. 339

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig.4. The probability of rice under stress in the study area 341

Table 2 Confusion matrix of identification of rice under stress based on SVM 342 Classification Unstressed rice (p<=0.5) Stressed rice (p>0.5) User’s accuracy Unstressed rice(p<=0.5) 357 32 91.77% Stressed rice(p>0.5) 32 1202 97.41% Producer’s accuracy 91.77% 97.41% Overall accuracy 96.05% Kappa coefficient 0.90 343

4.3. Classification of rice crops under Cd stress 344

Based on the spatio-temporal characteristics of different stressors in rice, CSTV was 345

used as prior information to distinguish between Cd stress and abrupt stress. Since the 346

selection of the threshold of CSTV may influence the area percentage of different stressors in 347

M

AN

US

CR

IP

T

AC

CE

PT

ED

rice and the accuracy of Cd assessment, thresholds of CSTV ranging from 1 to 3 were used, 348

with intervals of 0.1, leading to a variation in areal percentage and the accuracy of Cd 349

assessment (Fig.5 (a) and (b)). Fig.5 (a) shows that the higher threshold of CSTV due to the 350

higher confidence level can result in a gradual decrease in the areal percentage of rice under 351

abrupt stress and thus increase the areal percentage of rice under Cd stress. In other words, 352

changing the threshold of CTSV can affect the view on the differences in distribution of rice 353

under different stressors. Eq. (8) was used to calculate the accuracies shown in Fig.5 (b). It is 354

evident that the accuracy gradually improved with an increasing CSTV threshold value from 355

1 to 2.7, and that the accuracy remained constant when the CSTV value ranged from 2.7 to 3. 356

Thus, CSTV =2.7 was selected as most suitable threshold for identifying the different 357

stressors in rice. It can be noted that there is no disparity in accuracy with each differing 358

CSTV value. Further analysis showed that greater variations in accuracy occurred with 359

threshold intervals of 0.5, namely between 1and 1.5, between 1.5 and 2, between 2 and 2.5 360

and between 2.5 and 3. 361

The iterative Bayesian method was used to calculate the probability of rice under Cd 362

stress using the CSTV with a regular threshold of 2.7 in each period and probability of 363

stressed rice. The spatial distribution of Cd stress in rice is shown in Fig.5(c). The probability 364

of Cd stress in rice being greater than 0.5 occurred in most of the region, across 68% of the 365

total area of rice pixels. 366

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig.5 Spatial distribution of rice under Cd stress 368

369

The probability of Cd stress was compared with in situ Cd concentrations in soil 370

according to accuracy assessment method presented in Section 3.3.3. Four clusters were 371

obtained by setting threshold lines at Cd=1 and p=0.5 (Fig.6), with the largest concentrated 372

cluster located at Cd>1 and p>0.5, while the other sampling points were scattered across the 373

diagram. Further statistics regarding accuracy are presented in Table 3. A satisfying 374

performance was obtained with an overall accuracy of 81.57%. It is concluded that the Cd 375

concentrations in the soil are above the Level III Soil Environmental Quality Standards in 376

China, resulting in a higher possibility of crops becoming stressed. This is in agreement with 377

the strategy of the Soil Environmental Quality Standards in China. 378

M

AN

US

CR

IP

T

AC

CE

PT

ED

Fig.6 Comparison of the probability of Cd stress in rice and Cd concentrations in soil 380

Table 3 Classification accuracy of Cd stress in rice on a regional scale 381

Cd in soil (mg/kg) Number Probability Classification Accuracy

>1 32 <=0.5 5 84.37% >0.5 27 <=1 6 <=0.5 4 66.67% >0.5 2 Overall accuracy 81.57% 382 5. Discussion 383

5.1. Spatio-temporal characteristics of different stressors 384

In this study, CSTV of NDRE was calculated to capture the spatio-temporal 385

characteristics of different stressors (Liu et al., 2018). The temporal scale was based on the 386

cloud-free remote images at all growth stage of rice; while the suitable spatial scale of the 387

characteristic of different stressors was captured using Moran’s I. Based on Eq. 3, Moran’s I 388

of NDRE was calculated for each grid on a cloud-free day (i.e., 16th in September, 2017). The 389

change in spatial scale for monitoring stress in rice is displayed in Fig.S1. With increasing 390

M

AN

US

CR

IP

T

AC

CE

PT

ED

grid numbers in this study area from 22 to 214, the mean and minimum value of NDRE 391

Moran’s I kept decreasing gradually for each grid, whereas the maximum value of NDRE 392

Moran’s I kept increasing. The variation in standard deviation of NDRE Moran’s I kept 393

fluctuating. The minimum value of the standard deviation of NDRE Moran’s I was obtained 394

with a grid number of 26. Thus, the 26 grid number forms the optimal spatial scale in this 395

study. 396

Based on Eq. 2, the CSTV of NDRE was calculated at a different time. The ratio of the 397

same (or similar) growth period in different years for CSTV of NDRE was regarded as 398

indicator for an inter-annual anomaly to detect nutrient stress in rice.In this study, the data 399

collected at similar growth stages from the Sentinel-2 images on 28 August 2016 and 21 400

August 2017 were compared. The inter-variation of CSTV spatio-temporal characteristics is 401

displayed in Fig. S2(a), with the CSTV evidently showing little inter-annual variability, as the 402

CSTV for the majority of the area (about 95%) ranges from 0 to 2 while the area where the 403

CSTV is greater than 2 covers only approximately 5% of all rice pixels. This is due to the fact 404

that CSTV induced by soil nutrition can vary spatially and from year to year, and thus it 405

makes the distribution of abnormal VI over the same growth period in different years remains 406

localized. Therefore, it was concluded that nutrient stress in rice crops is not a routine event 407

and it can be detected when introducing a spatio-temporal anomaly level. 408

The CSTV of minus NDRE was calculated at two consecutive growth stages as shown 409

in Fig. S2 (b), (c), (d), (e) and (f). It can be clearly seen that the most CSTV in the rice area at 410

the earlier growth stage was concentrated between 0 and 1(Fig.S2 (b) and (c)). Comparatively 411

higher CSTV values (above 2) were present in local areas when the rice was in its later 412

M

AN

US

CR

IP

T

AC

CE

PT

ED

growth stages (Fig. S2 (d), (e) and (f). Detailed statistics of the areal percentages for each 413

CSTV classification are shown in Table S3. The area where the CSTV exceeded 2 was wider 414

at a later growth stage of the rice than at the earlier growth stage. Regardless of the growth 415

stage of rice in this study area, the CSTV of NDRE only exhibited greater values in the local 416

space. 417

5.2. Performance of the Bayesian method 418

This paper poses a representation of Cd stress in crops based on a probabilistic approach. 419

The initial step is to use a Bayesian method to detect stress levels in rice expressed by their 420

satellite-derived vegetation index. In general, crop stress is defined as: “any environmental 421

factor potentially unfavorable to crops” (Levitt, 1980). Rice growing in real agroecosystems 422

faces uncertainties whether they are stressed or not, so in the SVM classifier, RBF was used 423

for SVM kernel function to produce a probabilistic output, and to calculate the corresponding 424

conditional probabilities of classifying stressed rice. A high accuracy was obtained based on 425

the probabilistic output of the SVM classification (Table 2). The SVM classifier proved to be 426

effective for classification, especially for estimating conditional probabilities for variables, 427

where the training dataset was limited (Chang and Lin, 2011; Li et al., 2016). 428

Our results confirmed that vegetation index alone is insufficient to accurately identify 429

Cd stress in rice from various environmental stressors. The further step is to estimate the 430

probability of Cd in the soil causing the rice to be stressed, using the previous 431

satellite-derived vegetation index as input variables, as well as being iteratively updated with 432

spatio-temporal characteristics. Bayesian methods can make good use of various types of 433

available information, including prior information and temporal-spatial characteristics of 434

M

AN

US

CR

IP

T

AC

CE

PT

ED

different stressors. Moreover, a major advantage of using the iterative Bayesian method to 435

estimate Cd stress is that it directly accounts for the vegetation index at current growth stage 436

as a function of the vegetation index at upcoming growth stage using CSTV. The CSTV 437

indicated that Cd stress in rice is not only dependent upon stress-induced rice growing at a 438

specific growth stage but also on the spatial distribution of the stressed rice. The iterative 439

Bayesian method proved to be capable of monitoring the dynamics of the plant (Reiche et al. 440

2015, 2018). 441

The Bayesian method proved successful for detecting Cd stress in rice, when the 442

procedure was implemented using the collaborative system that integrated two classification 443

tasks, namely stressed rice and Cd stress in rice. The Bayesian method provided a flexible 444

framework for specifying expert knowledge or prior information and fitting complex 445

temporal-spatial characteristics that would have been difficult to achieve with frequentist 446

approaches. 447

5.3. Detection of regional Cd stress in rice crops 448

To monitor macro-level heavy metal pollution in agro-ecosystems, accurate and 449

large-scale information is needed, such as provided by satellite images. Hyperion imagery 450

proved feasible for monitoring large-scale rice under heavy metal contamination by mapping 451

the spatial distribution of Cu, Cd, as were As and Pb in Eucalyptus leaves (Liu et al., 2010; 452

Khalili et al., 2015). However, multispectral images have not been employed to directly 453

monitor plants under heavy metal stress on a large scale. Sentinel-2 is the first optical Earth 454

observation mission of its kind to include three bands in the “red edge,” providing key 455

information on vegetation state. Red edge refers to the wavelength range from 680 to 750 nm, 456

M

AN

US

CR

IP

T

AC

CE

PT

ED

within which the spectral reflectance of vegetation experiences an abrupt rise due to the 457

absorption of red radiation by chlorophyll and the strong reflection of red and near infrared 458

radiation (Horler et al., 1983). Some studies have shown that heavy metals can successfully 459

be studied in vegetation using red edge parameters, such as red edge position and red edge 460

area (Collins et al., 1983; Li et al., 2010, Clevers et al., 2004; Chi et al., 2006). Compared 461

with the previous SPOT and Landsat images, Sentinel-2 images are thus more useful for 462

monitoring vegetation healthy state, and hold great promise for accurately assessing Cd stress 463

in crops due to the provided “red edge” information combined with refined spatial and 464

temporal resolutions (Liu et al., 2018). 465

In this study, the probability of Cd stress was verified on the basis of collecting and 466

analyzing 41 soil samples, because crop heavy metals concentration is positively correlated to 467

the soil heavy metals concentrations (Zhao et al., 2012; McBride, 2002). In addition, our 468

multi-year field investigation showed that some stress symptoms of rice were displayed in 469

this study, such as relative lower chlorophyll and higher temperature in stressed rice canopy 470

virus healthy rice (Liu et al., 2016, 2018), when soil heavy metal concentrations are higher 471

than Level III Soil Environmental Quality Standards in China (Zhao et al., 2015). This study 472

confirms that multi-temporal or multi-year reflectance data can be useful to identify 473

spatio-temporal evidence of different stressors on crops. When crops are exposed to abrupt 474

stressors, the CSTV has a higher value since the VI is higher or lower than normal for a given 475

time interval within the growing season and therefore different CSTV indicates that crops are 476

infected by abrupt or stable stress at some consecutive growth stages. This may be a large 477

area where Cd stress on crops occurs persistently in all growth stages every year. Pests or 478

M

AN

US

CR

IP

T

AC

CE

PT

ED

diseases, on the other hand, do not occur during all growth stages of rice (Bushnell, 1984; 479

Lee, 1983). Some researchers confirm that spectral and spatio-temporal evidence is designed 480

to distinguish different crop types (Mukashema et al., 2014), and retrieve land use 481

classification information (Walde et al., 2014; Zhang and Du, 2015; Li et al., 2016). More 482

factors should be incorporated to improve the accuracy of such classification. 483

Although we obtained a high overall classification accuracy using our approach, further 484

improvement can be achieved, keeping the following aspects in mind. (i) Changes in 485

spectrum reflectance across time may be caused by land use changes, cloud cover or abrupt 486

stressors in crops. Therefore, it is important to discard misjudgments arising from unwanted 487

factors, such as cloud-covered or non-rice pixels. (ii) Regularly, several stress factors are 488

simultaneously active in crops, such as the frequently occurring combination of heat, heavy 489

metal stress and pest or disease in certain periods. It is difficult to distinguish between an 490

abrupt stressor and multi-stressors. Further research should improve the classification 491

accuracy through the use of highly complex signal analytical methods and multi-year datasets 492

to delineate between stable stress signals and abrupt stress signals. (iii) The threshold of 493

CSTV varies with the intensity of abrupt stress factors, crop species, their growing 494

environment, temporal resolution and spatial scale. The selection of this CSTV threshold 495

value is therefore dependent on the monitoring requirements of the user, and is always 496

associated with uncertainty. 497

6. Conclusion 498

Using satellite-derived vegetation index alone is insufficient to accurately detect 499

stressors of crops. However, Cd stress in crops is characterized by stability in space and time, 500

M

AN

US

CR

IP

T

AC

CE

PT

ED

thus differing from other stressors that are typically more transient. In our study, vegetation 501

index and spatio-temporal characteristics are employed as important evidence of Cd stress in 502

crops and are combined to gain the occurrence probability of Cd stress in rice from a 503

Bayesian network model. Firstly, NDRE is selected as sensitive indicator for screening stress 504

levels in rice crops. Secondly, many attempts have been made to identify the optimal spatial 505

scale and the CSTV proved to be successful in detecting rice under Cd stress and abrupt 506

stress with a suitable threshold. Lastly, a Bayesian method is used to effectively detect the 507

probability of Cd stress occurring in rice by using rice under stress and spatio-temporal 508

characteristics of different stressors as prior information derived from CSTV. 509

In this study, the Bayesian method is shown to be useful and effective for monitoring 510

crop stress. The possibility of rice being under stress can first be determined. Prior 511

information on the spatio-temporal characteristics of different stressors can then help 512

differentiate Cd stress in rice further. This study infers that a Bayesian method is well suited 513

as a stress assessment method, by integrating satellite-derived vegetation index with 514

spatio-temporal characteristics of different stressors, thus marking a significant step forward 515

in plant stress diagnosis. 516

517

Reference 518

Aguilera, P., Fernández, A., Fernández, R., Rumí, R., Salmerón, A., 2011. Bayesian networks 519

in environmental modelling. Environ. Model. Softw. 26, 1376–1388. 520

Bauer, S., Korč, F., Förstner, W. 2011. The potential of automatic methods of classification 521

to identify leaf diseases from multispectral images. Precision Agric., 12(3), 361-377. 522

M

AN

US

CR

IP

T

AC

CE

PT

ED

Borsuk, M. E., Stow, C. A., Reckhow, K. H. 2004. A Bayesian network of eutrophication 523

models for synthesis, prediction, and uncertainty analysis. Ecol.Modell., 173(2), 524

219-239. 525

Bushnell, W. 1984. 15–Structural and Physiological Alterations in Susceptible Host Tissue. 526

Cereal Rusts., 477-507. 527

Chang, C.-C., Lin, C.-J., 2011. LIBSVM: a library for support vector machines. ACM Trans. 528

Intell. Syst. Technol. 2 (3), 27:1–27:27. 529

Chen, D., Guo, H., Li, R., Li, L., Pan, G., Chang, A., Joseph, S., 2016. Low uptake affinity 530

cultivars with biochar to tackle Cd-tainted rice - A field study over four rice seasons in 531

Hunan, China. Sci. Total Environ. 541, 1489-1498. 532

Chen, S., Kao, C., 1995. Cd induced changes in proline level and peroxidase-activity in roots 533

of rice seedlings. Plant Growth Regul. 17 (1), 67-71. 534

Cheng, H., Wang, R., 2010. Semantic modeling of natural scenes based on contextual 535

Bayesian networks. Pattern Recognit. 43, 4042–4054. 536

Chi, G., Chen, X., Shi, Y., Liu, X., 2009. Spectral response of rice (Oryza sativa L.) leaves to 537

Fe2+ stress. Sci China. C. Life Sci. 52(8), 747-753. 538

Clevers, J., Kooistra, L., Salas, E., 2004. Study of heavy metal contamination in river 539

floodplains using the red-edge position in spectroscopic data. Int. J. Remote Sens. 25 540

(19), 3883–3895. 541

Collins, W., Chang, S. H., Raines, G. L., Canney, F., Ashley, R., 1983. Airborne 542

biogeophysical mapping of hidden mineral deposits. Econ. Geol. 78(4), 737-749. 543

Dlamini, W. 2010. A Bayesian belief network analysis of factors influencing wildfire 544

M

AN

US

CR

IP

T

AC

CE

PT

ED

occurrence in Swaziland. Environ.Modell.Softw, 25(2), 199-208. 545

Donatelli, M., Magarey, R., Bregaglio, S., Willocquet, L., Jpm, W., Savary, S. 2017. 546

Modelling the impacts of pests and diseases on agricultural systems. Agric. Syst., 155, 547

213-244. 548

Dorner, S., Shi, J.,.Swayne, D. 2007. Multi-objective modelling and decision support using a 549

Bayesian network approximation to a non-point source pollution model. 550

Environ.Modell.Softw. 22(2), 211-222. 551

Drǎguţ, L., Csillik, O., Eisank, C., Tiede, D., 2014. Automated parameterisation for 552

multi-scale image segmentation on multiple layers. ISPRS J. Photogramm. Remote Sens. 553

88, 119–127. 554

Drǎguţ, L., Tiede, D., Levick, S.R., 2010. ESP: A tool to estimate scale parameter for 555

multiresolution image segmentation of remotely sensed data. Int. J. Geogr. Inf. Sci. 24, 556

859–871. 557

Espindola, G.M., Camara, G., Reis, I.A., Bins, L.S., Monteiro, A.M., 2006. Parameter 558

selection for region-growing image segmentation algorithms using spatial 559

autocorrelation. Int. J. Remote Sens. 27, 3035–3040. 560

Food and Agriculture Organization (FAO) (2008). FAO and Sustainable Intensification of 561

Rice Production for Food Security. Rome, Italy: FAO. 562

Gilles, L., Marie-Hélène, J., François, G. 2008. Diagnosis tool for plant and crop N status in 563

vegetative stage: Theory and practices for crop N management. Eur. J. Agron., 28(4), 564

614-624. 565

Hernández-Rabadán, D., Ramos-Quintana, F., Guerrero, J. 2014. Integrating SOMs and a 566

M

AN

US

CR

IP

T

AC

CE

PT

ED

Bayesian classifier for segmenting diseased plants in uncontrolled environments. Sci. 567

World J., 2014, 1-13. 568

Holloman, C., Bortnick, S., Morara, M., Strauss, W.,.Calder, C. 2004. A Bayesian 569

Hierarchical Approach for Relating P M.Sub.2.5 Exposure to Cardiovascular Mortality 570

in North Carolina. Environ.Health.Presp, 112(13), 1282-1288. 571

Horler, D., Dockray, M., Barber, J. 1983. The red edge of plant leaf reflectance. Int. J. 572

Remote Sens., 4(2), 273-288. 573

Hosack, G., Hayes, K., Dambacher, J. 2008. Assessing Model Structure Uncertainty through 574

an Analysis of System Feedback and Bayesian Networks. Ecol. Appl., 18(4), 575

1070-1082. 576

Hu,R., Wang,S., Lin,J., Zhu,X., 2013. Current Status of Soil Fertility of Croplands in Hunan 577

and Its Change in Recent Three Decades. Soils, 45, 585–590 (in Chinese). 578

Jensen, F., Lauritzen, S., Olesen, K. 1990. Bayesian Updating in Causal Probabilistic 579

Networks by Local Computations. Comp.Stat., 5(4), 269-282. 580

Johnson, B., Xie, Z., 2011. Unsupervised image segmentation evaluation and refinement 581

using a multi-scale approach. ISPRS J. Photogramm. Remote Sens. 66, 473–483. 582

Khalili, R., Anvari, S., Honarmand, M. 2015. Combination of biochemical and hyperspectral 583

remote sensing methods for detection of heavy metal pollutions in eucalyptus leaves 584

(case study: the city of bam). ISPRS Archives., XL-1(W5), 379-385. 585

Larsson, E., Bornman, J., Asp, H., 1998. Influence of UV-B radiation and Cd2+ on 586

chlorophyll fluorescence, growth and nutrient content in Brassica napus. J. Exp. Bot. 587

49(323), 1031-1039. 588

M

AN

US

CR

IP

T

AC

CE

PT

ED

Lee, F. 1983. Rice sheath blight : A major rice disease. Plant Dis., 67(7), 829-832. 589

Lemaire, G., Jeuffroy, M., Gastal, F. 2008. Diagnosis tool for plant and crop N status in 590

vegetative stage: Theory and practices for crop N management. Eur. J. Agron, 28(4), 591

614-624. 592

Levitt, J. 1980. Responses of plants to environmental stresses. Academic Press. 1(5), 593

3642-3645. 594

Li, M., Stein, A., Bijker, W., Zhan, Q. 2016. Urban land use extraction from Very High 595

Resolution remote sensing imagery using a Bayesian network. ISPRS-J. Photogramm. 596

Remote Sens., 122, 192-205. 597

Li, N., Lue, J., Altermann, W., 2010. Applications of spectral analysis to monitoring of heavy 598

metal-induced contamination in vegetation. Spectrosc. Spect. Anal. 30,2508–2511. 599

Lichtenthaler, H. K. 2010. The Stress Concept in Plants: An Introduction. Ann. New. York. 600

Acad. Sci, 851(1), 187-198. 601

Liedloff, A.,.Smith, C. 2010. Predicting a ‘tree change’ in Australia's tropical savannas: 602

Combining different types of models to understand complex ecosystem behaviour. 603

Ecol.Modell. 221(21), 2565-2575. 604

Lin, H. T., Lin, C. J., Weng, R.C., 2007. A note on platts probabilistic outputs for support 605

vector machines. Mach. Learn. 68 (3), 267–276. 606

Lindsay ,W.; Norvell, W., 1978. Development of a DTPA soil test for zinc, iron, manganese, 607

and copper. Soil Sci. Soc. Am. J. 42(3), 421-428. 608

Liu, M., Liu, X., Ding, W., Wu, L., 2011a. Monitoring stress levels on rice with heavy metal 609

pollution from hyperspectral reflectance data using wavelet-fractal analysis. Int. J. Appl. 610

M

AN

US

CR

IP

T

AC

CE

PT

ED

Earth Observ. Geoinform. 13(2), 246-255. 611

Liu, M., Liu, X., Li, M., Fang, M., Chi, W., 2010. Neural-network model for estimating leaf 612

chlorophyll concentration in rice under stress from heavy metals using four spectral 613

indices. Biosyst. Eng. 106(3), 223-233. 614

Liu, M., Liu, X., Wu, L., Duan, L.,.Zhong, B. 2011b. Wavelet-based detection of crop zinc 615

stress assessment using hyperspectral reflectance. Comput. Geosci., 37(9), 1254-1263. 616

Liu, M., Liu, X., Wu, M., Li, L., Xiu, L., 2011c. Integrating spectral indices with 617

environmental parameters for estimating heavy metal concentrations in rice using a 618

dynamic fuzzy neural-network model. Comput. Geosci. 37(10), 1642-1652. 619

Liu, M., Liu, X., Zhang, B., Ding, C., 2016. Regional heavy metal pollution in crops by 620

integrating physiological function variability with spatio-temporal stability using 621

multi-temporal thermal remote sensing. Int. J. Appl. Earth Observ. Geoinform. 51, 622

91-102. 623

Liu, M., Wang, T., Skidmore, A.K., Liu, X., 2018. Heavy metal-induced stress in rice crops 624

detected using multi-temporal Sentinel-2 satellite images. Sci. Total Environ. 637–638, 625

18–29. 626

Liu, Z., Cheng, J., Huang, W., Li, C., Xu, X., Ding, X., Zhou, B., 2012. Hyperspectral 627

discrimination and response characteristics of stressed rice leaves caused by rice leaf 628

folder. In: Li, D.L., Chen, Y.Y. (Eds.), Computer and Computing Technologies in 629

Agriculture V, Pt II. Springer, Berlin, Heidelberg, Germany, 528–537. 630

Marcot, B., Steventon, J., Sutherland, G., McCann, R. 2006. Guidelines for developing and 631

updating Bayesian belief networks applied to ecological modeling and conservation. 632

M

AN

US

CR

IP

T

AC

CE

PT

ED

Can J.Forest Res., 36(12), 3063-3074. 633

McBride M.B., 2002. Cadmium uptake by crops estimated from soil total Cd and pH. Soil Sci. 634

167, 62-67. 635

Mee, C., Balasundram, S., Hanif, A. 2017. Detecting and monitoring plant nutrient stress 636

using remote sensing approaches: A review. Asian J.Plant Sci., 16(1), 1-8. 637

Mukashema, A., Veldkamp, A., Vrieling, A. 2014. Automated high resolution mapping of 638

coffee in Rwanda using an expert Bayesian network. Int J. Appl. Earth. Obs., 33(12), 639

331-340. 640

Nielsen, T., Jensen, F., 2007. Bayesian Networks and Decision Graphs. Springer 641

Science.Business Media. 642

Norse, D., Ju, X., 2015. Environmental costs of China’s food security. Agric. Ecosyst. 643

Environ. 209, 5–14. 644

Overmars, K.P., De Koning, G.H.J., Veldkamp, A., 2003. Spatial autocorrelation in 645

multi-scale land use models. Ecol. Modell. 164, 257–270. 646

Park, M., Stenstrom, M. 2008. Classifying environmentally significant urban land uses with 647

satellite imagery. J Environ Manage, 86(1), 181-192. 648

Pearl, J., 1988. Probabilistic reasoning in intelligent systems. San Mateo, CA: Morgan 649

Kaufmann. 650

Phadikar, S. 2012. Classification of Rice Leaf Diseases Based on Morphological Changes. Int. 651

J. Infor. Electron. Eng., 2(3), 460-463. 652

Raza, S., Smith, H., Clarkson, G., Taylor, G., Thompson, A., Clarkson, J.,.Rajpoot, N. 2014. 653

Automatic detection of regions in spinach canopies responding to soil moisture deficit 654

M

AN

US

CR

IP

T

AC

CE

PT

ED

using combined visible and thermal imagery. Plos One, 9(6). 655

Reiche, J., de Bruin, S., Hoekman, D.H., Verbesselt, J., Herold, M., 2015. A Bayesian 656

approach to combine Landsat and ALOS PALSAR time series for near real-time 657

deforestation detection. Remote Sens. 7, 4973–4996. 658

Reiche, J., Hamunyela, E., Verbesselt, J., Hoekman, D., Herold, M., 2018. Improving 659

near-real time deforestation monitoring in tropical dry forests by combining dense 660

Sentinel-1 time series with Landsat and ALOS-2 PALSAR-2. Remote Sens. Environ. 661

204,147–161. 662

Ren, H., Zhuang, D., Singh, A., Pan, J., Qiu, D., Shi, R. 2009. Estimation of As and Cu 663

contamination in agricultural soils Around a mining area by reflectance spectroscopy: a 664

case study. Pedosphere, 19(6), 719-726. 665

Sankaran, S., Ehsani, R. 2011. Visible-near infrared spectroscopy based citrus greening 666

detection: Evaluation of spectral feature extraction techniques. Crop Prot., 30(11), 667

1508-1513. 668

Schikora, M., Schikora, A., Kogel, K., Koch, W., Cremers, D. 2010. Probabilistic 669

Classification of Disease symptoms caused by Salmonella on Arabidopsis Plants. GI 670

Jahrestagung., 2, 874-879. 671

Scudiero, E., Skaggs, T., Corwin, D., 2015. Regional-scale soil salinity assessment using 672

Landsat ETM + canopy reflectance. Remote Sens. Environ. 169, 335-343. 673

Skidmore, A. K., 1989 .An expert system classifies eucalypt forest types using thematic 674

mapper data and a digital terrain model. Photogramm. Eng. Remote Sens., 55(10), 675

1449-1464. 676