Integrating cellular automata, artificial neural network, and fuzzy set theory to simulate threatened orchards: application to Maragheh, Iran

Abstract

Urbanization processes challenge the growth of orchards in many cities in Iran. In Maragheh, orchards are crucial ecological, economical, and tourist sources. To explore orchards threatened by urban expansion, this study first aims to develop a new model by coupling cellular automata (CA) and artificial neural network with fuzzy set theory (CA–ANN–Fuzzy). While fuzzy set theory captures the uncertainty associated with transition rules, the ANN considers spatial and temporal nonlinearities of the driving forces underlying the urban growth processes. Second, the CA–ANN–Fuzzy model is compared with two existing approaches, namely a basic CA and a CA coupled with an ANN (CA–ANN). Third, we quantify the amount of orchard loss during the last three decades as well as for the upcoming years up to 2025. Results show that CA–ANN–Fuzzy with 83% kappa coefficient performs significantly better than conventional CA (with 51% kappa coefficient) and CA–ANN (with 79% kappa coefficient) models in simulating orchard loss. The historical data shows a considerable loss of 26% during the last three decades, while the CA–ANN–Fuzzy simulation reveals a considerable future loss of 7% of Maragheh’s orchards in 2025 due to urbanization. These areas require special attention and must be protected by the local government and decision-makers.

Keywords:

• threatened orchards
• urbanization
• cellular automata
• artificial neural network
• fuzzy set theory
• urban planning

1. Introduction

One of the key ecological pressures impacting on the environment and ecosystem is land-use change (LUC) due to rapid urbanization (Czamanski et al. 2008; Grimm et al. 2008; Haase and Nuissl 2010; Shafizadeh Moghadam and Helbich 2015). Big cities in Iran face tremendous urban expansion, resulting in degradation of natural areas in favor of built-up areas (Jarnagin 2004; Tayyebi, Pijanowski, and Tayyebi 2011a). The urbanization level¹ increased from 31% in 1956 to 70% in 2011 (Iran Statistical Center 2013). Rural–urban migration and rapid population growth are two essential causes of this urban expansion. In addition, macroeconomic policies in the real estate sector interfere with these urban growth processes. From 1990 to 2008, house prices increased by approximately 330% (Meidani, Zabihi, and Ashena 2011). As the prices of orchard products increased slightly compared to business or construction land prices (Khiabani 2004), it is monetarily more beneficial for farmers to sell their property for business and construction purposes. Therefore, due to limited unoccupied lands for urbanization, several fast-growing Iranian cities urbanized toward orchard lands (Tayyebi, Pijanowski, and Pekin 2011b; Zadeh 2008), causing an excessive shrinkage of orchard lands. While knowledge about the rate of orchard conversion is of significance for policy-makers and urban planners (Su et al. 2012; Zadeh 2008), no study has quantified the past and future orchard loss in Iran.

In this regard, LUC models (Verburg et al. 2004; Estoque and Murayama 2014) are valuable tools to detect ecological areas that are threatened due to urban growth (Shafizadeh Moghadam and Helbich 2015; Jokar Arsanjani, Helbich, and Mousivand 2014). Even though several LUC models have been proposed in the literature, including, for example, statistical approaches (Shafizadeh Moghadam and Helbich 2015) and geometrical approaches (Tayyebi, Pijanowski, and Pekin 2011b; Jokar Arsanjani, Helbich, and Mousivand 2014), among these models, cellular automata (CA; Santé et al. 2010) and linguistic approaches (Mantelas et al. 2012) such as fuzzy logic-based approaches (Grekousis, Manetos, and Photis 2013) and ANNs (Tayyebi, Tayyebi, and Khanna 2014; Basse et al. 2014) have received considerable attention due to their dynamic nature, being explicitly spatial with a highly flexible structure.

CA is widely used to simulate urban dynamics (Batty, Xie, and Sun 1999). CA calibration obtains a set of transition rules. However, extensive search strategies are required to explore appropriate combinations of transition rules during model calibration; otherwise, the process results in inaccurate predictions. The determination of transition rules is even more complex when additional spatial variables are considered. As a consequence, the calibration of CA models is challenging due to the determination of optimal weights for a large number of involved factors (Dietzel and Clarke 2006). Therefore, computer search algorithms such as ANNs (Li and Yeh 2002) are valuable to be coupled with CA to significantly reduce the model calibration time. For example, Li and Yeh (2001) used an ANN embedded in a CA model to simulate urbanization in China. Similarly, Almeida et al. (2008) combined ANN with CA for modeling intra-urban land-use dynamics. Finally, Shafizadeh Moghadam and Helbich (2013) utilized Markov chains and CA to simulate spatiotemporal urbanization processes in Mumbai, India. In the present study, we made use of this body of knowledge and integrated a CA model with ANN to enhance the model calibration time.

A limitation of conventional CA is the determination of complex relationships between the driving forces and their weights. ANNs allow modeling complex nonlinear relationships (Pijanowski et al. 2002, 2009) and weight estimation of individual driving forces in CA simulation (Lio and Yeh 2001; Tayyebi et al. 2014b). For example, a comparison among ANN, classification tree, and multivariate adaptive regression splines by Tayyebi et al. (2014b) found that ANN outperformed the other models in the short run (e.g. 5 years) and for coarse-resolution data (e.g. 1 km); however, similar accuracies are achieved by all models in the long run (e.g. 20 years) and for fine-resolution data (e.g. 30 m). A follow-up study (Tayyebi and Pijanowski 2014) reconfigured these three data mining models to concurrently simulate multiple land-use classes (agriculture, forest, and urban areas) in Wisconsin and Michigan, United States.

However, most of the recent LUC studies have mainly focused on the development of complex models while ignoring uncertainty in data and models (Tang and Pannell 2009). As highlighted by Tayyebi, Tayyebi, and Khanna (2014), uncertainty in LUC modeling originating from data is greater than other uncertainty dimensions (e.g. model parameter and model outcome). To take uncertainty into account, Wu (1996) proposed transition rules based on heuristics and fuzzy set theory to capture uncoordinated land development. Others, such as Pijanowski et al. (2007), coupled a back-cast ANN with groundwater models to quantify the contribution of land uses to the groundwater signal arriving at streams. Results showed that uncertainty in models and back-cast ANN affected the modeling outcome as well as their reliability for natural resource and land-use planning. Therefore, LUC models may not exactly represent the relationship between the dependent and independent drivers due to inherent uncertainties observed in existing data (Yeh and Li 2006; Christman et al. 2015). Based on these findings, we also integrated CA and ANN models with fuzzy set theory to capture the uncertainty associated with data as well as transition rules.

Based on these identified methodological research gaps and the importance for policy-makers to protect orchards from ongoing urbanization in Maragheh, Iran, the prime objectives of this paper are as follows: first, the development of a CA with ANN and utilizing fuzzy set theory for urban analysis, virtually not yet employed in empirical studies so far, second, to compare the power of the CA–ANN–Fuzzy model with two existing approaches, namely basic CA and CA–ANN, and third, to explore the extent of threatened orchards during the last three decades as well as for the upcoming year 2025 by predicting threatened locations. The remainder of the paper is structured as follows. While Section 2 introduces the model integration, Section 3 describes the study area and data. Next, Section 4 discusses the empirical results, and finally Section 5 highlights major conclusions.

2. Methods

This section introduces the existing models, including CA and CA-ANN, as well as the newly developed CA–ANN–Fuzzy model. This section also presents how three models were calibrated using historical data and details description of the questionnaire used to identify significant drivers of urbanization.

2.1. Conventional CA

CA encompasses the following components (Clarke, Hoppen, and Gaydos 1997; Arun and Katiyar 2013): spatial arrangement of cells having a specific state, neighborhood associations, and transition rules (Barredo et al. 2003). Time in CA is considered to be discrete while space is frequently represented using two-dimensional grid cells (Shafizadeh Moghadam and Helbich 2013). Each cell has a finite number of states, representing different land-use types. A state of a cell depends on the landscape transitions in the immediate neighborhood. At any point in time, each cell has a probability of being converted form one land use type to another, where the probability is given by means of a neighborhood function defined by a priori-determined transition rules that involve both adjacent cells as well as several contributing driving forces (Qi et al. 2004; He et al. 2008).

2.2. Artificial neural network

ANNs are machine-learning techniques capable of capturing nonlinear associations underlying land-use transformations by formulating the relationships between input and output variables (Pijanowski et al. 2002; Grekousis and Photis 2014). The multilayer perceptron, as used in this paper, typically includes an input layer, a hidden layer, and an output layer. These three layers are connected to each other in a feed-forward manner (Pijanowski et al. 2009; Tayyebi, Perry, and Tayyebi 2014a, 2014b). Neurons in the hidden layer assign the relationships between the input and output nodes and pass the result through a nonlinear activation function to the output neuron or adjacent neuron of the corresponding hidden layer.

ANNs learn through training processes, which find the relations between the input variables and the output variables (Pijanowski et al. 2002). During the training, the initial model parameters are modified repeatedly until model outcomes represent the observed output as accurately as possible (Liu 2009). The most popular learning algorithm used in geospatial simulation is backpropagation (Hagenauer and Helbich 2012). This learning algorithm is essentially a search procedure that attempts to minimize an error function (Wang and Li 2011). ANNs start by randomly assigning the weightings and calculating the mean squared error (MSE), after which the cycle continues until a terminating criterion is met. The MSE is measured by calculating the difference between the output from the network and the actual output presented in the training phase (Pijanowski et al. 2014). If the level of MSE is not attained, the error is distributed back to the neurons in the hidden layer, thus allowing them to update the weightings and mitigate the error (Wang and Li 2011). A suitability map is then produced using an activation function that applies a neuron which receives the weighted output from connected neurons of the preceding layers.

2.3. Integrating CA, ANN, and fuzzy set theory

CA–ANN–Fuzzy is a multi-model approach developed in this study. This approach comprises the following four steps: (1) fuzzy membership functions are applied to quantify the uncertainty associated with the driving forces; (2) an ANN captures the complex relationships between the normalized driving forces that encourage urbanization to encroach on the orchards and determines CA transition weights; (3) the normalized driving forces and their optimized weights (obtained from ANN) input to a fuzzy CA model to identify orchard conversion potential for built-up areas; and 4) to estimate the amount of threatened orchards through mean and standard deviation normalization. An overview of the methodological workflow is given in Figure 1.

Figure 1. Data preparation and modeling steps for simulating and predicting orchard loss.

2.3.1. Integrating fuzzy set theory into CA

As land-use transformations are complicated due to uncertainties (i.e. human decision-making behavior), fuzzy set theory is an attractive methodological extension to incorporate fuzziness in the development of CA transition rules and LUC monitoring (Santé et al. 2010). In a fuzzy CA, the transition of a cell state to another one is defined by means of a set of fuzzy logic-constrained rules, which are defined in a nondeterministic and probabilistic manner (Wu 1999). A fuzzy set is characterized by a fuzzy membership grade that scales from 0.0 to 1.0, where a value of 0 refers to nonmembership and a value of 1 represents complete membership (Zadeh 1992). Eastman (2009) argues that fuzzy sets provide a standardized measure and avoid the selection of a priori unknown Boolean constraints or cut-off values. The following fuzzy membership function has been applied:

(1)

where is the state of cell having a value between 0 and 1, is the variable for which the membership function is defined, and as well as are the lower (0) and higher (1) values for a cell.

In accordance with fuzzy set theory (Zadeh 1992), LUC is a continuous process where the state of a cell is defined by a membership function. Fuzzy CA uses the LUC probability between two cells (ij) as shown in Equation (2) (Li et al. 2011, 2013):

(2)

where is the state of cell at time , is the conversion possibility at time t for cell , is land-use conversion force of a cell at time , is the force weight of cell (calculated with ANN in this study: see Section 2.2), is the stochastic disturbance factor, is a uniform random variable within the range [0,1] in accordance with fuzzy membership function, and is the dispersion parameter to control the size of the stochastic perturbation.

The result of fuzzy CA is a nondeterministic output that shows the degree of cell conversion. However, it is required that the cells having high potential for conversion are shown as real values, such as threatened sites. To achieve this, a mean and standard deviation normalization are used to normalize fuzzy CA as shown in Equation (3):

(3)

where is the output value of or the fuzzy CA that ranges between −1 and 1 for cell and is the standard deviation for cells of at time . When all values are transformed in this fashion, the threatened sites are determined as given in the following Equation (4):

(4)

where 1 refers to threatened sites and 0 to unchanged areas.

2.3.2. Integration of fuzzy CA and ANN

The objective of CA calibration is to find the best combination of transition rules to model LUC (Batty, Xie, and Sun 1999). Transition rules have an important role in calibrating CAs, which might result in different simulation outputs (Yeh and Li 2001). It is challenging to define the best combinations of transition rules when there are many variables, because LUC patterns are complex (Luo, Yu, and Xin 2008). The variations are due to many possible ways of defining the transition rules. Land-change modelers usually integrate CA with other data mining tools to perform CA calibration. ANNs are appropriate algorithms for the calibration of CA parameters (Li and Yeh 2002) and usually have a higher overall accuracy than other models (e.g. logistic regression) for calibrating LUC models (Lin et al. 2011). It has been shown that ANNs as mathematical models are robust in calibrating CA models (Li et al. 2011; Wang and Li 2011; Li et al. 2013; Lee, Wu, and Lei 2013). Thus, in this study, we used ANN with a multilayer perceptron (see Section 2.2 for more details) to calibrate the fuzzy CA model.

2.4. Questionnaire survey to obtain driving forces

To identify the driving forces that encourage urbanization of the orchards in Maragheh, we distributed a questionnaire among farmers, housebuilders, real estate agents, natural resource conservation organization (NRCO), and city authorities who worked in the garden protection sector. The sample was selected according to the following criteria: age (50 years older and above), continually resided in Maragheh for the last 15–40 years and considered by the locals to have sufficient knowledge about the orchards’ destruction and their correlation with urbanization. The questionnaire focused on 31 socioeconomic and spatial factors. The respondents were clustered to 30 farmers, 15 housebuilders, 12 real estate agents, 13 NRCO employees, and 20 employees of Maragheh municipality.

We used personal interviews guided by the questionnaire, which were answered by ticking on a 5-point Likert scale. The interview-based survey provided an in-depth understanding of the local actors’ knowledge about the driving forces and details about the main causes for the mentioned forces. Three essential questions of survey were: (1) Which factors are important in the conversion of orchards to built-up areas? (2) What is the logic behind these factors? 3) What is the relative importance of factors in orchards conversion to urban areas? Questions 1 and 3 were covered by the questionnaire and we interviewed the respondents for question 2.

After the questionnaire survey was carried out, factor analysis was used in SPSS to extract important factors. Factor analysis is a statistical method to explain the correlations between the observations in terms of the underlying factors, which are not directly observable (Yu et al. 2003). It is often applied to identify a small number of latent variables that describe most of the variance embedded in a larger number of observed variables (Thompson 2004).

2.5. Model calibration

It is necessary to calibrate the LUC model with reliable historic data before using the model for prediction purposes (Pijanowski et al. 2014). To calibrate the three models (CA, CA–ANN, and CA–ANN–Fuzzy), training data was used, and the test set was used to evaluate the performance of the three models. We used Kappa statistic to determine how well the model outputs matched the reference change maps. Kappa calculates the percentage success of a model relative to chance between two categorical datasets (Doswell, Davies-Jones, and Keller 1990). Since many LUC model outputs are categorical, Kappa is useful for comparing model outcome and reference change map using Equation (5) (Monserud and Leemans 1992; Hagen-Zanker and Martens 2008).

(5)

where is the percent correct for modeling output, the expected correct due to chance, and the percent correct when modeling is perfect, which is 100 or 1. It is generally considered (Tayyebi, Pijanowski, and Pekin 2015) that Kappa values for map agreement are as follows: >80% is excellent; 60–80% is very good; 40–60% is good; 20–40% is poor and <20% is very poor.

3. Study area, data preparation, and model implementation

3.1. Study area

The study area is Maragheh city, one of the oldest cities in Iran. It is located on the northwest side of the country (Figure 2). The area of the city is approximately 18.66 km². In terms of population, it is a mid-sized city and situated in the northeast of Ourmieh Lake and west of the protected fossil site. Maragheh is well-known for its favorable climate and attractive orchard landscape. The Safi Chay River crosses the east of Maragheh and creates a cultivated plain that makes Maragheh famous as a “Garden City” in Iran. In addition to valuable ecological features, orchards are important for the following reasons: (1) products of the orchards are an essential part of Maragheh’s economy; (2) the orchards are main landscape features and a key symbol of urban identity; and (3) the orchards are a crucial tourism attraction in Maragheh.

Figure 2. Study area located in northwest of Iran. For full color versions of the figures in this paper, please see the online version.

3.2. Data preparation

Three Landsat TM and ETM+ images (30 × 30 m) in 1984, 1998, and 2012 were used to extract land-use maps. The selection of the years was based on data availability as well as the following characteristics: the reduction of sessional effects, regular time intervals between the images, and cloud-free images. The three Landsat images (1984, 1998, and 2012) were classified using fuzzy ARTMAP algorithm (Carpenter, Grossberg, and Reynolds 1991), since this method usually provides a higher accuracy compared to alternative classification methods (Tong and Pern 2000; Muchoney and Williamson 2001; Azari and Reveshty 2013). We used stratified random sampling with 500 points over the entire landscape for training sites. These 500 points represented seven land-use classes, including: built-up areas, orchards, water, dry farming, grassland, irrigated agriculture land, regolith, and wasteland (Figure 3; Universal Transverse Mercator WGS 1984 zone 39 N). Reference data for the accuracy assessment of the classification outputs originate from the following sources: field observation, topographic maps, and aerial photos. The Kappa statistics were 88%, 93%, and 90% for 1984, 1998, and 2012, respectively.

Figure 3. Land-use maps of Maragheh in (a) 1984, (b) 1998, and (c) 2012.

The questionnaire focused on 31 socioeconomic and spatial factors. Factor analysis revealed that the following variables were relevant: (1) vegetation; (2) land prices; (3) distance to central business district (CBD); (4) slope; and (5) accessibility to transportation networks.

3.2.1. Vegetation

The normalized difference vegetation index (NDVI) based on Landsat data was computed to reveal the variety of vegetation characteristics that include quantity, productivity, biomass, etc. (Equation (5); Myneni et al. 1995):

(5)

where are the band 4 and band 3 of Landsat ETM⁺ data, respectively. The NDVI varies between –1 and 1. Areas containing dense vegetation usually show values close to 1, indicating a reduced potential for orchard conversion to housing. In contrast, areas containing sparse vegetation show NDVI values closer to −1, indicating that these areas have a high potential for the orchard conversion to urban areas (Myneni et al. 1995).

3.2.2. Land prices

Land prices in the neighborhoods as well as the rate of price changes in the last 10 years were extracted from the database of a real estate agency. The major purpose was to have general knowledge of land prices in the neighborhoods and the rate of price changes over the last 10 years. The questionnaire respondents were selected according to their willingness to provide information and whether they continually worked in real estate businesses in the last 10 years and had sufficient knowledge about land price changes. Land prices and their intensity changes in the last 10 years were mapped at local neighborhoods and rasterized at a 30 m resolution.

3.2.3. Distance to the central business district

To create CBD layer, buffer function was used in GIS to calculate the distance from the closest CBD (Table 1).

Table 1. Data sources and types.

Dataset	Source	Date	Resolution
Landsat images (TM 1984, ETM 1998, ETM+ 2012)	US Geological Survey	1984, 1998, 2012	30 m
NDVI	US Geological Survey	2012	30 m
Digital elevation model	Aster (NASA)	2012	30 m
Land price Map	Own survey	2012	30 m
Transportation network	Extracted from country road maps	2010	30 m
CBD	Extracted from land use map	2012	30 m

3.2.4. Slope

Slope values were extracted from Aster digital elevation map (DEM) dataset. The DEM has a 15 × 15 m spatial resolution.

3.2.5. Accessibility to main transport networks

To create accessibility to the main transport network layer, the buffer function was used in GIS to calculate the distance from the closest transport networks.

3.3. Model implementation

ANN was used to obtain optimized weights of CA parameters using the MATLAB software. Two groups of data were used to build ANN: (1) historical LUCs between 1984, 1998, and 2012; and (2) land-use conversion forces. Before training the ANN, three primary tasks were carried out: First, the data were converted to the same spatial resolution (the spatial resolution of 30 m because Landsat TM and ETM⁺ data were used). Second, as the data had different ranges, we normalized them into the range between 0 and 1 using a linear membership function in fuzzy set theory. Third, the data were divided into training and testing data (Tayyebi and Pijanowski 2014). We randomly selected 70% of the entire data as training data and the remaining data was reserved as testing data (Tayyebi, Tayyebi, and Khanna 2014).

The ANN consists of an input layer, one hidden layer, and an output layer (Tayyebi et al. 2013). In this study, we used five inputs (i.e. land-use types and spatial variables mentioned before) and one output (representing the orchards change and no-change in the historical LUC in Tables 2 and 3). We also combined water and protected areas in a new layer called “exclusionary zone” and masked these areas out of the simulation and further analysis. After obtaining the optimal weights for the three models in the training run, the whole data were used to simulate conversion of the orchards. IDRISI Selva version (Clark Labs, Clark University, Worcester, MA, USA) was utilized for CA simulation. We verified the predictive ability of the three models through the simulation of orchard conversion in 2012 based on their status in 1998 and comparison with reference land-use map.

Table 2. LUC between 1984 and 2012. The numbers of cells and percent values are given.

	Built-up	Orchard	Water	Dry farming	Wasteland	Grassland	Irrigated agriculture	Total	Gains and losses
Built-up	6078	2438	363	2591	2569	4613	1901	20553	13853
	7.29	2.93	0.44	3.11	3.08	5.53	2.28	24.66	16.62
Orchard	114	12978	204	311	118	53	793	14571	−5028
	0.14	15.57	0.24	0.37	0.14	0.06	0.95	17.48	–6.03
Water	239	1087	171	145	203	92	193	2130	1283
	0.29	1.3	0.21	0.17	0.24	0.11	0.23	2.56	1.54
Dry farming	148	1947	70	18657	2981	2717	2316	28836	715
	0.18	2.34	0.08	22.38	3.58	3.26	2.78	34.6	0.86
Wasteland	77	154	11	2442	2694	3045	393	8816	−458
	0.09	0.18	0.01	2.93	3.23	3.65	0.47	10.58	–0.55
Grassland	26	150	10	3844	699	1744	591	7064	−5254
	0.03	0.18	0.01	4.61	0.84	2.09	0.71	8.48	–6.3
Irrigated agriculture	18	845	18	131	10	54	303	1379	−5111
	0.02	1.01	0.02	0.16	0.01	0.06	0.36	1.65	–6.14
Total	6700	19599	847	28121	9274	12318	6490	83349	0
	8.04	23.51	1.02	33.74	11.13	14.78	7.79	100

Table 3. LUC between 1998 and 2012. The numbers of cells and percent values are given.

	Built-up	Orchard	Water	Dry farming	Wasteland	Grassland	Irrigated agriculture	Total	Gains and losses
Built-up	11541	1543	325	2344	1886	1683	1231	20553	7203
	13.85	1.85	0.39	2.81	2.26	2.02	1.48	24.66	8.64
Orchard	244	13776	129	186	6	6	224	14571	−5721
	0.29	16.53	0.15	0.22	0.01	0.01	0.27	17.48	–6.87
Water	421	1055	209	197	31	54	163	2130	1349
	0.51	1.27	0.25	0.24	0.04	0.06	0.2	2.56	1.62
Dry farming	876	2593	80	19217	1295	2415	2360	28836	−945
	1.05	3.11	0.1	23.06	1.55	2.9	2.83	34.6	–1.13
Wasteland	126	111	6	2950	2701	2316	606	8816	2694
	0.15	0.13	0.01	3.54	3.24	2.78	0.73	10.58	3.23
Grassland	82	205	24	4719	197	888	949	7064	−313
	0.1	0.25	0.03	5.66	0.24	1.07	1.14	8.48	–0.37
Irrigated agriculture	60	1009	8	168	6	15	113	1379	−4267
	0.07	1.21	0.01	0.2	0.01	0.02	0.14	1.65	–5.12
Total	13350	20292	781	29781	6122	7377	5646	83349	0
	16.02	24.35	0.94	35.73	7.35	8.85	6.77	100

4. Results

4.1. Analysis of data from interviews

To select the potential driving forces for orchard conversion to built-up areas, questionnaire data were analyzed using factor analysis. The five driving forces described in Table 4 could explain 74.6% of the total variation.

Table 4. Factor analysis statistics for elicited forces.

Variables	Eigenvalue	Variance percent (%)
Vegetation	8.67	21.14
Land price in the neighborhoods	7.321	17.54
Slope	6.781	14.21
Accessibility to transport networks	6.374	12.14
CBD	5.21	9.57

To assess the main reasons stated by the respondents, to extract the driving forces, interviews were carried out with the selected sample. Vegetation factor explained 21.14% of variation in the orchard conversion to built-up areas. In Maragheh, the key reasons for orchard destruction were found to be monetary benefits. Land regulations prohibit orchard conversion for construction purposes; the orchards are kept dried by owners who can then receive the permit to build houses for their wasted lands. Therefore, the orchards with low-density vegetation have a high potential to be changed to built-up areas.

Land prices in the neighborhoods are an essential factor in tendency of landowners to sell their lands for construction. This factor explained 17.54% of variations. According to the land price survey, the north and northwest neighborhoods had the highest land prices and farmers were more likely to sell their orchards for built-up purposes in these areas. Several factors potentially led to urbanization and orchard conversion, including distance to CBD, slope, and accessibility to transportation networks. Slope variable explained 14.21% of variation. Steep slopes were found to be less attractive for urban expansion, and the cells showing a pronounced slope have a low potential for urban development compared to lower slopes. Therefore, the orchards which are located in the steep and irregular slopes have a less potential for urbanization.

Most of the respondents mentioned that new settlements are more probable to be built near transport networks. Therefore, the orchards that have a small distance to transport links have a high potential of conversion to built-up areas. This variable explained 12.14% of variation. The orchards near the CBD have a higher destruction potential than the orchards that are further apart. This variable explained 9.57% of variation. Seventeen percent of the orchards that were converted to built-up areas lie in the first 400 m, 39% lie in the second 400 m, and 24% lie in the third 400 m. The amount of the orchard conversion to built-up areas reduced in the following rings. Figure 4 shows the CBD influence on orchard conversion to built-up areas in 1984–2012.

Figure 4. CBD influence on orchards conversion to built-up areas between 1984 and 2012.

4.2. Historical LUC

We used a contingency table (Christman et al. 2015) to summarize the temporal and spatial changes in Maragheh. Tables 4, 5a and b show the LUC in three time periods: 1984–1998, 1984–2012, and 1998–2012. Also, the tables show gains and losses of each land-use class. In the time period 1984–1998, 36.95% of the land uses changed to other land-use classes. This amount increased to 41.87% in 1998–2012. Water, dry farming, and wasteland were the classes that showed higher stability in 1984–2012. The proportion of built-up areas increased notably compared to other classes from 1984 to 2012 (16.62% increased). The proportion of built-up areas increased by 7.98% and 8.64% in the first (1984–1998) and the second periods (1998–2012), respectively. In contrast, orchards, irrigated agriculture, and grassland reduced, with a loss of more than 6% for each category in 1984–2012. Most losses of grassland occurred in the first time period (1984–1998), but orchards and irrigated losses intensified in second period (1998–2012). Between 1998 and 2012, more than 5721 cells (6.78%) of the orchards were converted to gain other categories, especially built-up areas (1543 cells) and dry farming (2593 cells).

Table 5. LUC between 1984 and 1998. The numbers of cells and percent values are given.

	Built-up	Orchard	Water	Dry farming	Wasteland	Grassland	Irrigated agriculture	Total	Gains and losses
Built-up	5967	1079	344	863	2234	1911	952	13350	6650
	7.16	1.29	0.41	1.04	2.68	2.29	1.14	16.02	7.98
Orchard	185	17429	265	312	176	56	1869	20292	693
	0.22	20.91	0.32	0.37	0.21	0.07	2.24	24.35	0.84
Water	128	254	181	17	77	1	123	781	−66
	0.15	0.3	0.22	0.02	0.09	0	0.15	0.94	–0.08
Dry farming	90	269	17	20925	2619	4233	1628	29781	1660
	0.11	0.32	0.02	25.11	3.14	5.08	1.95	35.73	1.99
Wasteland	164	98	4	1097	2900	1596	263	6122	−3152
	0.2	0.12	0	1.32	3.48	1.91	0.32	7.35	–3.78
Grassland	31	34	0	2912	567	3665	168	7377	−4941
	0.04	0.04	0	3.49	0.68	4.4	0.2	8.85	–5.93
Irrigated agriculture	135	436	36	1995	701	856	1487	5646	−844
	0.16	0.52	0.04	2.39	0.84	1.03	1.78	6.77	–1.02
Total	6700	19599	847	28121	9274	12318	6490	83349	0
	8.04	23.51	1.02	33.74	11.13	14.78	7.79	100

4.3. Simulation of orchard loss

Figure 5 shows the MSE across the training cycles of ANN. The MSE starts at around 0.05 and drops linearly thorough three cycles. We stopped training at 70 cycles, where the MSE was 0.015 at epoch 64. After factor weights were optimized by ANN, we obtained the neighborhood size through empirically testing the model accuracy. The CA, CA–ANN, and CA–ANN–Fuzzy were run in different neighborhood sizes with the same combination of input data. The comparison revealed that a 7 × 7 neighborhood size had the highest accuracy (Table 6). Kappa coefficients were 0.51, 0.79, and 0.83 for CA, CA–ANN, and CA–ANN–Fuzzy, respectively. The simulation results were acceptable for CA–ANN and CA–ANN–Fuzzy since overall Kappa statistics for these models were over 0.80. However, we used CA–ANN–Fuzzy with the best performance to simulate future orchard conversion.

Table 6. Accuracy for CA, CA–ANN, and CA–ANN–Fuzzy with different neighborhood sizes.

Model	Neighborhood size	Kappa	P-values
CA	9 (3 × 3)	0.39	<0.0001
	25 (5 × 5)	0.49	<0.0001
	49 (7 × 7)	0.51	<0.0001
CA–ANN	9 (3 × 3)	0.61	<0.0001
	25 (5 × 5)	0.77	<0.0001
	49 (7 × 7)	0.79	<0.0001
CA–ANN–Fuzzy	9 (3 × 3)	0.63	<0.0001
	25 (5 × 5)	0.79	<0.0001
	49 (7 × 7)	0.82	<0.0001

Figure 5. Training run of ANN MSE and MSE value across cycles.

After model verification, we simulated urban developments for the year 2025 to estimate the quantity of potential orchard loss due to urbanization. We used CA–ANN–Fuzzy to simulate urbanization in 2025 since this model had the best performance compared to other models. This simulation is based on spatial variables, land-use map for 2012, data combination based on ANN results, and a 7 × 7 neighborhood size. The first output of the model is a suitability map showing the locations that are more likely to be converted from orchards to urban areas. A pixel value closer to 1 indicates a higher change probability to urban areas while values closer to 0 are more likely to stay as orchards. Figure 6 shows orchard conversion potential obtained by fuzzy CA–ANN modeling.

Figure 6. Suitability map of orchard conversion to built-up areas obtained by fuzzy CA–ANN.

We specifically identified the threatened sites or the orchards that are more endangered and require protection by authorities. Since CA models lack the ability to determine the actual amount of change, mean and standard deviation normalization were utilized to convert fuzzy CA output to simulated change. According to the CA–ANN–Fuzzy model, 0.963 km² or 7.34% of Maragheh orchards will be destructed as a consequence of urbanization processes in 2025 (Figure 7). Results also showed that the threatened sites are located along the edge of urban and suburban areas (in less than 1 km from the edge of the built-up area).

Figure 7. Predicting and comparing the threatened orchard in 2025 using (a) CA, (b) CA–ANN, and (c) CA–ANN–Fuzzy.

5. Discussion and conclusions

Orchards of Maragheh play a significant role from ecological, economical, and tourist perspectives. In the recent years, however, the demand for built-up areas caused encroachment on the orchards and damaged ecological resources. This study attempted to determine the impacts of urbanization on orchard areas and identify areas at risk, to preserve the orchards. Assessment of the historical land-use data revealed a substantial decline in orchard areas. The historical data showed a considerable loss of 25.65% during the last three decades. The predicted map indicated that other 7.34% of Maragheh orchards will be replaced by urban constructions by 2025.

This study also compared the power of three land-use models (CA, CA–ANN, and CA–ANN–Fuzzy) for simulating orchard loss. This study showed the flexibility of CA to be integrated with other models, such as fuzzy set theory and ANN, to identify areas that are at risk and need to be protected. Results also showed that both CA–ANN and CA–ANN–Fuzzy performed better than conventional CA. This is due to the fact that conventional CA cannot detect nonlinear spatial and temporal patterns in land-use data, while CA–ANN and CA–ANN–Fuzzy both use ANN to identify nonlinear patterns in data. While both CA–ANN and CA–ANN–Fuzzy were able to simulate complex patterns of orchard loss, CA–ANN–Fuzzy performed slightly better than CA–ANN. This is due to the fact that fuzzy set theory captured the uncertainty associated with transition rules.

This study also demonstrated that incorporating questionnaires in spatial LUC models can provide better simulation results. The results of factor analysis conducted on survey data showed that human-induced forces such as house prices are highly influential and drive the urbanization to encroach on the orchards. Generally, more pronounced land prices increase the tendency to sell orchards for housing purposes. According to Azari and Reveshty (2013), there has been a moderate correlation (42%) between accessibility to main transport links and built-up areas in Maragheh during 1984–2011. Therefore, the orchards, being highly accessible, are facing high conversion potential for urban construction. Jadidyan (2012) argues that accessibility and building commercial areas are essential factors for the construction of built-up areas near the transportation network.

Results showed that the threatened orchards located in suburban areas are next to the existing built-up areas. Besides ecological, economical, and tourist issues, the suburban orchards function as a green belt. Suburban orchards also reduce air pollution caused by factories located in the surroundings of Maragheh. Therefore, these areas require more attention from authorities to ensure the ecological integrity of such landscapes. For future research, we suggest the use of high-resolution remote-sensing images to explore changes on a more detailed scale.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Mehdi Azari http://orcid.org/0000-0002-4156-6022

Notes

1. It reflects the ratio of urban population to total population.