Comparison of artificial neural network and multiple linear regressions efficiency for predicting soil salinity in Yazd -Ardakan plain, central Iran

Soil salinity is one of the main concern in land degradation in the arid and semiarid regions (Tajgardan et al 2010).The extent of primary salt-affected soils in worldwide is about 955 M ha, while secondary salinization affects some 77 M ha with 58% in irrigated lands (Metternicht. & Zinck 2003). Accumulation of dispersive cations such as dissolved sodium in the soil solution and at the exchange phases (Ca+, Mg++, K+ ) affects soil physical properties such as soil structure, aggregate stability, infiltration rate and erosivity (Navarro‐Pedreño et al 2007), leading to decline in soil quality. The main effects of soil degradation are seen as declining productivity and quality, which adversely affect quality of life for those who are dependent on land as a source of livelihood. Direct measurement of Electrical Conductivity (EC) using traditional methods is time consuming and costly whereas application of remote sensing technology in mapping and monitoring degraded lands, especially in salt-affected soil in addition to having great potential in temporal and spatial scale, are timesaving, accurate and cost effective (Khan et al 2005).
There are different methods to model soil mapping using spectral data. For example multivariate techniques and artificial neural networks (ANN) as calibration tools in chemo metrics, can be applied to model the relationship between various soil chemical and physical  properties with spectral data (Were et al 2015, Yang et al 2003).
Multiple linear regressions (MLR) are frequently used to study soil properties (Zornoza et al 2007, Kalkhajeh et al 2012, Adams et al 2004) especially soil salinity (Arshad et al 2013, Tajgardan et a 2010, Lake et al 2009). MLR is generally helpful for observational studies where treatments are not imposed and the multiple explanatory variables tend not to be orthogonal with themselves.
ANN models have been effectively applied in many areas ranging from economics to science. ANNs are being widely used in soil science research. Using ANN models, some valuable information about locations is possible to collect where there are no available soil maps, by combining soil mapping data from other areas with landscape features to be representative for the spatial variation of soils. Hence ANNs are complex computer programs that able to model complicated functional relationships and applied for soil mapping problems using a set of variables related to soil types as training data in order to extend rules to unmapped area (Freire et al 2013). ANN as a strategy for digital mapping of soils associated with a certain predictive method can be found in some studies. Arruda et al., (2016) showed that digital techniques have the potential to match the natural soil distribution in the landscape. Bagheri Bodaghabadi et al., (2015) reported that the integration of data-mining methods and digital terrain analysis could provide more robust results for predicting soil units. Elarabi and Abdelgalil., (2014) reported that Artificial Neural Networks (ANN) have the acceptable ability to predict the soil classification and soil parameters in Sudan. Behrens et al., (2005) reported that Artificial Neural Networks (ANN) and digital terrain analysis are time- saving and cost effective and provides remarkable results.
This study was conducted to (i) evaluate the potential of ANN in predicting soil salinity and compare the results with those obtained by MLR using ASTER data, and (ii) to determine the importance of ASTER band ratios in explaining soil salinity by sensitivity analysis. The study was done in Yazd-Ardakan plain, Yazd province, located in central Iran.

2. Methods and Materials  
2. 1.  The Study area

The study area was Ardakan plain located in Yazd Province(central Iran) (Fig.1), with area of about 300000 ha extending from east of 53°   24´ - 54°   56´ longitude and a north of  31°   13´- 32°   36´ latitude. Elevation is 1500 m above the mean sea level. Most of the population and also farmlands and industrial towns are concentrated in this area.
Minimum and maximum annual rainfall are 0-14/6 mm (August-March), and daily temperatures are 6.38- 32.6 °C (Jan-Jul). Most of the area is on a plateau (Pediment and Alluvial Plain) characterized by gently rolling and undulating topography. The major land uses within the study area include agricultural uses and rangelands.
 The geological material of the area is composed of quaternary sediments and Neogene marl (that greatly influence soil salinity) and argillaceous sedimentary rocks (shale, slate, limestone) with some salt dome and salt and gypseous domes. Marl components have surrounded a relatively extensive area and include plentiful plaster and channel passing through the area that able to adversely affect the land resources especially Eocene  components (containing plaster) and salty domes  (even though their area is small). Based on the boundary of study area, the factors affecting on salinity of the border area and the less important components remained outside the boundary.
 
Figure.1 Location of the study site(Yazd-Ardakan plain) in Yazd province, central Iran.

2.2.  Soil sampling and laboratory analyses
Selection of the best band combination is usually the first step to compress remote sensed data by securing class separation (Metternicht. & Zinck 2003) In this study, VNIR and SWIR bands of ASTER were used. Table 1 contains the description of VNIR and SWIR bands. To assign sampling points, the best color composite representing the change of soil and geology was selected and the best ASTER three-band combination was 7-4-2 in RGB, respectively. Then, the sampling points were selected based on the spectral changes to find suitable special pattern of sampling points within the study area (Fig 2)

 
Table1  The description of VNIR and SWIR bands of Aster data
 
Subsystem    Band No    Spestral range
(µm)    Spatial resolution(m)    Quantization levels
VNIR    1    0.52-0.6    15    8 bits
    2    0.63-0.69        
    3N    0.78-0.86        
    3B    0.78-0.86        
SWIR    4    1.60-1.70    30    8 bits
    5    2.145-2.185        
    6    2.185-2.225        
    7    2.235-2.285        
    8    2.295-2.365        
    9    1.360-2.430        
TIR    10    8.125-8.475    90    12 bits
    11    8.475-8.825        
    12    8.925-9.275        
    13    10.25-10.95        
    14    10.95-11.65        
 

 
 
Fig.2 the color composite (7-4-2) for soil sampling (Aster Image) in Ardakan plain.
 

A total of 133 (0-10 Cm) soil samples were collected. At each point, five samples with about 10 m apart were collected and made into one composite sample to reduce micro variability. Soil samples were dried and passed through a 2-mm sieve. Electrical Conductivity (EC) was used to measure saturated extract.

2.3.  Processing of the Aster image
The VNIR and SWIR bands were geo-referenced with topography maps using image to map registration method and rectified to the UTM zone 40 cartographic projections. Nearest neighbor re-sampling model was applied with the absolute spatial accuracy (30 m, 15 m). The remote satellite images were processed by ILWIS software. In addition to main bands, some appropriate processing operations include first component selected from principal component analysis of VNIR and SWIR bands and standard principal component analysis of VNIR and SWIR data (Ren. & Abdelsalam 2006); NDVI (Normalized Difference Vegetation Index) was obtained from VNIR bands (Rouse Jr 1974) and selected ratio arithmetic bands were applied to determine and quantify salinity (Table 2).The spectral numbers (DN) of each sample point on main and processed bands were extracted using extract function in IDRISI software. These spectral numbers were used for statistical analysis.
 
Table 2 Main and processed bands
Bands    (bands made procedure)    No.
VNIR1............SWRI6    Main bands of ASTER except TIR bands 1    1
PCA1vnir    First component of PCA of VNIR bands    2
PCA1swir    First component of PCA of SWIR bands    3
PCA1vnir&swir    First component of PCA of VNIR&SWIR bands    4
NDVI    (vnir 3-vnir 2)/(vnir 3+vnir 2) 5    5
SUB48    (swir1/swir5)    6
SUM48    (swir1+swir5)    7
Sum19    (vnir1+swir6)    8
NDVI48     (swir1-swir5)/ (swir1/swir5)    
 

2.4  MLR and ANN modeling
The MLR is a method to model the linear relationship between dependent and one or more independent variables. The dependent and independent variables are called the predicted and predictors respectively. MLR is based on least squares: the model is calculated so that the sum-of-squares for difference between observed and predicted values is minimized (Mata 2011).
In present study, all sampled points were located in the original images and synthetic images were extracted. The results showed that all spectral values were normally distributed. Pearson linear correlation was used to determine the relationship between EC and spectral values in all bands, with a confidence level of 95%. Regression equation for those variables showed higher significant correlation with digital numbers. Regression analysis was performed using the following equation (1):
                                               (1)
where RSInd is the index derived from the remote sensing data and f is the regression function (Douaoui et al 2006). Spatial distribution of ECe was used to calculate this equation.
In this study, multilayer perceptron (MLP) with back-propagation learning rule was used. The MLP network [also termed back-propagation (BP) network](Fig. 3) is the most popular network to monitor engineering problems in the case of nonlinear mapping, (Haykin. & Network 2004).
Multi-layer Perceptron (MLP) is a supervised learning algorithm that learns a function   by training a dataset, where   is the number of dimensions for input and   is the number of dimensions for output. Given a set of features   and a target  , it can learn a non-linear function approximate for either classification or regression.

 

 
Figure. 3 Structure and flow chart for artificial neural network :One hidden layer MLP.
 

The learning process was performed using the well-known BP algorithm. The standard BP algorithm is based on the delta learning rule (Rumelhart. & Clelland 1986). Two main processes were performed in a BP algorithm, with a forward pass and a backward pass. In the forward pass, an output pattern was presented to the network and its effect propagated through the network, layer by layer. For each neuron, the input value was calculated by following Equation (2).
                                        (2)
where     is the input value of  th neuron in th layer; / is the weighted connection between  th neuron in   th layer and  th neuron in the ( )th layer;  is the output of  th neuron in the ( )th layer; /m is the number of neurons in the ( )th layer.
In each neuron, the value calculated from Equation (2) was transferred by an active function.
The common function for this purpose is the sigmoid function, calculated by Equation (3).
                                      (3)
The output of each neuron was computed and propagated through the sequential layers. Then, the final output of the network was prepared to compare with the target output.  In addition to correlation coefficient, the root mean square error (RMSE) as an objective function was calculated by Equation (4) (Degroot, 1986)
                                                
                                          (4)
Where  is the jth element of the target output related to the  th pattern;  is the output of  th neuron related to the  th pattern; / / is the number of patterns;   is the number of neurons in the output layer.

After computing the objective function, the second step of the BP algorithm, i.e. the backward process, was performed by back propagation of the network error in the previous layers. Using the gradient-descent technique, the weights were adjusted to reduce the network error by calculating Equation (5) (Rumelhart. & Clelland 1986)
 
                                  (5)                    
where  is the weight increment at the  th iteration (Epoch);  is the learning rate;  is the momentum term ( ).
This process was continued until the allowable network error was obtained. To design the artificial neural network, the measured field data were used. The number of available data was 133. The datasets were organized: 79 datasets (60%) were used for the learning process, 27(20%) datasets were used for testing, and 27(20%) for verification, respectively. The data sets for learning, testing, and verification processes were selected randomly at different points on landscape to avoid bias in estimation. The number of neurons in input and output layers depends on the independent and dependent variables, respectively. The network was designed with 18 parameters as input pattern and the EC as the output parameter. A network was designed to estimate EC. The number of hidden layers, number of neurons in the hidden layers, the parameter , and the number of iterations were selected by calibration through several test runs, and trial and error (Marquardt Levenberg learning rule). Various functions were tested for MLP neural networks and the tansigmoid function presented the best results. In order to identify the most important terrain attributes and band ratio representing EC, sensitivity analysis was done using the StatSoft method. A sensitivity ratio was calculated via dividing the total network error. A ratio greater than 1.0 implied that the variable plays an important role for the variability in EC. The higher the ratio, the more important the variable(Miao et al 2006).

3  Results and Discussion
3.1  Statistical characterization of data
 Descriptive statistics for ECe in the selected area showed 33.20 and 35.61dS/m as mean and standard deviation respectively (Table 3). The coefficient of variation (CV) for EC was 0.65. To define variability, we followed the system suggested by Wilding (Wilding 1985). According to this classification, any property with CV more than 0.35 is classified as high variable, therefore, EC indicated high variability in the selected area.








 
 
Table3: Statistical description of the electrical conductivity (n=133)
 
Variables    mean    Median    min    max    SD    CV
EC    33/20    17    2    136    35/61    0/65
 
3.2   MLR  and ANN modeling
Pearson linear correlation between digital numbers of Aster data and measured ECe values (Table 4) demonstrated that EC has significant correlation (p


 
Table4: Correlation coefficients between spectral indices and measured EC values
 
Spectral bands
    Correlation coefficients (EC )
SWIR1    -0.171*
SWIR2    -0.162*
SWIR3    -0.125*
SWIR4    -0.170*
SWIR5    -0.142*
SWIR6    -0.161*
VNIR1    -0.302**
VNIR2    -0.271**
VNIR3 -    -0.256**
NDVI    -0.177*
PCA1vnir    -0.176*
PCA1swir    -0. 182*
PCA1vnir&swir    -0.189*
SUM48    -0.173*
SUB48    0. 165*
 
**Correlation is significant at the 0.01 level. * Correlation is significant at the 0.05 level

 

The data obtained by optimal parameters (Table 5) of the final selected ANN models could be used to predict EC. The model had 18 input nodes and one output node. The hidden-layer node was 37. Also, the optimum iteration learning rate was determined as 7000. The ANN-EC models resulted in R2 and RMSE values of 0.79 and 0.11, respectively. Normalized predicted data against normalized observed data for testing data set were plotted (Fig.4) and the coefficients of determination (R2) were determined.

 
Table 5. Summary of the best results and optimum parameters of the artificial neural network modeling for estimating EC
 
Component    ANN structure    Transfer function    Iteration    Number of hidden layer    Number of hidden neurons
EC    18-37-1    Tansig    7000    1    37

 
 

 
Figure 4. Scatter plot displaying relationship between the measured and estimated values of the EC for test patterns of developed artificial neural network models in the selected site in Yazd-Ardakan plain, central Iran
The relative importance of soil and terrain attributes using sensitivity analysis based on sensitivity coefficients of the ANN models for EC is shown in Fig. 5. The variables with high values contributed greatly to the variability in soil salinity. Sum19 index was identified
as the most important factor to reveal EC. Other important factors for predicting EC are SWIR6, SWIR5, VNIR2, SUM48, SWIR1, VNIR3 with relative coefficients of sensitivity as 1, 0.99, 0.77, 0.73, 0.60, 0.59 and 0.48, respectively. NDVI and SWIR2 showed less sensitivity to other terrain and bands. Finally the results for ANN and MLR are summarized in Table 6.


 
 
Figure 5.   Histogram displaying the results of sensitivity analysis, relative sensitivity coefficients of selected soil and terrain attributes for soil salinity in Yazd-Ardakan plain, central Iran

Table 6.  Statistical validation of results by regression analyses and ANN for predicting EC in selected site (Yazd-Ardakan plain district, center of Iran).
MLR    ANN
EC=11.54+120VNIR1-169.6 NDVI48+0.7SWIR1-0.69VNIR2    y = 0.511x + 0.283        Model
            Statistical parameter
0.23    0.79    R2
0.33    0.11    RMSE

 
A comparative evaluation of the results by the two methods demonstrates that ANN model provides better accuracy as indicated by higher correlation coefficient and lower RMSE than those obtained by the linear regression method. In regression analyses, SWIR1, VNIR1, VNIR2 and NDVI48 were the main factors for predicting  soil salinity, but in the case of ANN, the most important factors for predicting EC were SUM19, SWIR6,SWIR5, VNIR2, SUM48, SWIR1 and VNIR3. According to nonlinear relationship between dependent and predictor variables, the neural network provided better prediction than regression equations.
This study contains promising results to model and map soil salinity based on ANN and ASTER images. In this way, Shahabi et al (2017) used ordinary kriging (OK), artificial neural networks (ANN) and multiple linear regressions (MLR) in order to model and predict the soil salinity in Dashte-Tabriz. The results showed that the ANN had the lowest RMSE and highest R2. Pachepsky et al (1996) investigated the accuracy of ANN and the regression method using correlation coefficient and the root mean square error. They reported that the neural network is able to predict the easily measurable soil parameters with better accuracy and less error. Benefits of remote sensing for identifying salt -affected soils and mapping them were improved by many researchers. Tayebi et al (2010) found VNIR and SWIR datasets of ASTER as the sustainable method for mapping and detecting salt-affected areas. Similar results have been reported by other researchers (Tajgardan et al 2007, Allbed et al 2014, Tayebi et a 2010)

4.  Conclusion
The present study reveals the combination of the measured and the remote sensing data (the bands and indices) into a regression model and an ANN model and offers a fast and inexpensive method to map and model the spatial variation in soil salinity .This combination into regression model was able to explain 23% of the spatial variation in the soil whereas ANN demonstrate 79% of variations in salinity of the selected area and this means that the ANN is more efficient in modeling and mapping than MLR. The higher correlation coefficient and lower root mean square revealed by backpropagation ANN clearly proves the greater efficiency of ANN in modeling rather than stepwise regression. Based on sensitivity analysis in neural network, SUM 19 and SWIR6 with 18 bands and their ratios are able to indicate soil salinity in study area.
The facility of developed model and degree of precision makes it as a desirable tool to predict soil salinity. Thus, this model can be applied by the decision makers in Ardakan plain and similar regions to perform or amplify effective soil reclamation programs that diminish or prevent high level of soil salinity. Although this study demonstrates that modeling and mapping soil salinity can be assumed with acceptable accuracy based on ASTER images, further research is required to consider the possibility of hyperspectral data in mapping and modeling soil salinity and investigating whether it can increase the accuracy in modeling and the mapping process.