SOIL INFILTRATION BASED ON BP NEURAL NETWORK AND GREY RELATIONAL ANALYSIS

Soil infiltration is a key link of the natural water cycle process. Studies on soil permeability are conducive for water resources assessment and estimation, runoff regulation and management, soil erosion modeling, nonpoint and point source pollution of farmland, among other aspects. The unequal influence of rainfall duration, rainfall intensity, antecedent soil moisture, vegetation cover, vegetation type, and slope gradient on soil cumulative infiltration was studied under simulated rainfall and different underlying surfaces. We established a six factor-model of soil cumulative infiltration by the improved back propagation (BP)-based artificial neural network algorithm with a momentum term and self-adjusting learning rate. Compared to the multiple nonlinear regression method, the stability and accuracy of the improved BP algorithm was better. Based on the improved BP model, the sensitive index of these six factors on soil cumulative infiltration was investigated. Secondly, the grey relational analysis method was used to individually study grey correlations among these six factors and soil cumulative infiltration. The results of the two methods were very similar. Rainfall duration was the most influential factor, followed by vegetation cover, vegetation type, rainfall intensity and antecedent soil moisture. The effect of slope gradient on soil cumulative infiltration was not significant.


INTRODUCTION
Soil infiltration is a key process of a natural water cycle.It is a movement of water vertically downwards into the soil, but also includes lateral and upward movement in the soil (Huang, 2000).The investigation of soil infiltration is very significant in theoretical and practical aspects and useful in many fields, e.g., for water resources assessment and estimation, runoff regulation and management, soil erosion modeling, nonpoint and point source pollution of farmland, etc. Numerous scholars have published a series of relevant papers on soil infiltration.Some of them carried out detailed studies on a physical model of soil infiltration (Green & Ampt, 1911;Philip, 1957a,b;Smith, 1972;Smith & Parlange, 1978;Lei et al., 1998).These physical models can describe the non-linear process of infiltration with some accuracy, but their parameters were complicated to solve and inconvenient to use.Of course, simplified infiltration formulas were proposed (Kostiakov, 1932;Horton, 1940;Fang et al., 1958;Jiang & Huang, 1986;Wang et al., 2003), which are easy to use and to simulate simple infiltration processes.However, the parameters in the simplified infiltration formulas have no actual physical meaning.In addition, some factors affect soil infiltration considerably, such as soil properties, underlying substrate, rainfall regime, soil moisture, land use patterns, etc. Branson et al. (1972) studied the effect of vegetation cover on soil infiltration and found that the effect was positive.Previous studies reported that there is no relation between soil infiltration and gradients (Singer & Blackard, 1982;Mah et al., 1992;Cerdà & García-Fayos, 1997).Assouline & Mualem (1997) modeled sealing due to rainfall as a function of rainfall intensity, the second moment of the drop-size density distribution, the maximal drop diameter, the compaction limit, and the initial shear strength of the soil, which depend on the initial soil bulk density and water content.Abu-Awwad (1997) studied the influence of surface crusts on water infiltration and redistribution, and found that sand columns could significantly increase the amount of moisture stored in the soil and reduce the amount of runoff.Some studies on infiltration, investigating mainly the effect of a single factor on soil infiltration have been carried out, but research on the important sequencing of multiple impact factors is relatively rare.
In this paper, we selected six factors: rainfall duration (Rd), rainfall intensity (Ri), antecedent soil moisture (Aw), vegetation cover (Vc), vegetation type (Vt) and slope gradient (Sg) as the main influencing factors of soil infiltration, and the BP neural network model of Ci was established, aside from the default factor used to test the sensitivity of each factor on cumulative infiltration (Ci).As the vegetation-related variables were not numerical, they were quantified first (Table 1).Bare land was assigned the arbitrary Vt value 0, while slopes with a legume (alfafa) and a graminaceous cover (wheat and ryegrass) were assigned Vt values of 1, 2 and 3, respectively.Grey relational analysis was used to calculate the grey correlation among six factors and Ci.Thus, the result of the BP default factor test was not only verified by the result of the grey relational analysis, but represented a selfvalidation of grey relational analysis as well.

Site and experiment description
The experiment was carried out by the Soil and This study was conducted in a slope-adjustable soil box with a variable-intensity rainfall simulator.The experimental set-up and soil preparation were as follows: The whole experimental set-up consisted of four parts: (1) a needle plate with about 650 needles to produce raindrops that can be adjusted easily.Two types of needles were used: 7# to simulate low rainfall intensity (0.5-1.0 mm min -1 ), and 8# for high rainfall intensity (1.0-2.0 mm min -1 ); (2) a water tank (1.7 × 1.2 × 0.25 m 3 ) to ensure a certain depth of water and raindrops; (3) a constant pressure watering device to ensure a constant water depth in the water tank; (4) a vibration motor, to produce a uniform distribution of raindrops across the soil boxes.Rainfall area and height were 1.5 × 1.0 m 2 and 1.2 m, respectively.Rainfall intensities ranged from 0.3 mm min -1 to 2.5 mm min -1 , mean drop velocity was 4.78±0.25 m s -1 , and mean rainfall kinetic energy per unit and time was 0.2193 ± 0.12 J m -2 s -1 .
According to the soil box size used in laboratory experiments in many previous studies (Poulenard et al., 2001;Dunjó et al., 2004;Molina et al., 2007;Vahabi & Nikkami, 2008;Kato et al., 2009;Nassif, 1975;Poesen et al., 1994;Fohrer et al., 1999;Adekalu et al., 2007), we used soil boxes measuring 1.2 ×0.8 ×0.45 m (47.24 in×31.50in×17.72 in).They were mounted on four wheels to facilitate movement and a jack allowed an adjustment of the slope angle from 0° to 30°.Moreover, the bottom of the boxes was perforated with many evenly distributed holes, allowing vertical drainage, and there was a guiding gutter along the bottom of the soil box to drain the slope runoff.
The clay content of the soil in Linghou, Wuquan Town, which is located on the third floodplain and terrace of the Weihe River, is very high; while the soil in the Water-saving Exhibition Garden of Yangling, which is located on the first floodplain and terrace of the Weihe River, is rich in sand particles.Both soils are quite different from those generally found in the study area on the Loess Plateau.To obtain a soil texture similar to that commonly found on the Loess Plateau, the soils from these two sites were mixed at equal weight proportion, gently crushed before passing through a 10 mm sieve and finally mixed thoroughly to minimize the difference between treatments.Each soil box was filled with 40 cm of soil in four 10-cm layers, at 1.3 g cm -3 bulk density.Additionally, each soil layer was raked lightly before packing the next layer to diminish the discontinuity.Several physical properties of the final experimental soil are shown in table 2.
Two pasture grass and one crop species were tested as vegetation cover: ryegrass (Lolium perenne L.) and alfalfa (Medicago sativa L.) as pasture grass, which are common perennial components of the vegetation on the Chinese Loess Plateau.Spring wheat (Triticum aestivum L.) was chosen for being a crop that is widely used in northern Shaanxi.The planting density of ryegrass and alfalfa was around 2 × 3 cm, and 4 × 3 cm for spring wheat.Pasture grass and wheat were sown in early April 2009 and mid-March 2009, respectively, by broadcast sowing.The experiment was initiated in mid-May 2009; ryegrass and alfalfa in the different soil boxes were cut two to five times throughout the experimental period.The treatments consisted of: three soil boxes with vegetation cover (ryegrass, alfalfa and spring wheat), and one control box of bare soil.
Four surface gradients (5 o , 10 o , 15 o , and 20 o ) were examined, representing the range of the majority of slopes on the Loess Plateau (Wu et al., 2001).Based on the typical rainfall regime in the Northern Shaanxi Loess Plateau area (Wang et al., 1998), five rainfall intensities were employed (0.5, 0.75, 1.0, 1.5, and 2.0 mm min -1 , i.e., 1.18, 1.77, 2.36, 3.54, and 4.72 in h -1 ).Each experimental treatment was replicated twice.The rainfall duration in the experiments ranged from 40 to 90 min (on average approximately 60 min), and each event was ended when the runoff reached a stable state (i.e. when the runoff amount of two consecutive time intervals was identical).The average rainfall amount was 57.1 ± 11.6 mm.Per treatment, 20 rainfall events were simulated, and the interval between two rainfall events was around 15 days.

Measurements
During each rainfall event, runoff and sediments were collected in a 1000 mL standard cylinder every 0.5-20 min; after settling of the turbid water, the clear water volume was regarded as the runoff amount.
Vegetation cover was measured by a photographic method.Firstly, a digital camera (Fuji J25) was used to take 3 to 5 JPG images to record the vegetation growth status.Care was taken to ensure that the Table 1.Quantification of the non-numerical vegetation type (Vt) variables camera lens was vertical to the slope surface and the lens always at the same distance from the object.Secondly, these photos were transformed into TIFF images with a lab color channel, using the software package Photoshop CS 3.0 TM .Finally, based on the color split and color gradation analysis by the software package Image-J TM , the vegetation cover was expressed as the ratio of the green to the total area of the photo.The accuracy and stability of this method was satisfactory (absolute error about 2 %) (Niu, 2009).
In addition, the result calculated by the photographic method was modified with a visual estimation method.
During the whole experiment, the vegetation cover was measured 45 times and the geometrical mean was around 61.0 %, with a maximum of 97.1 % and a minimum of 21.0 %.Generally, the vegetation cover with ryegrass was the densest of the three species.
Rainfall intensity was measured by a portable rain gauge.Furthermore, before each artificial rainfall event, rainfall intensity was calibrated three times at six measuring points and the rainfall uniformity coefficient was calculated by the following equation: where K represents uniformity coefficient, %; I e is the rainfall intensity at each measuring point, mm min -1 ; I m is the mean rainfall intensity of six measuring points (in mm min -1 ).The experiment could be run when K was not less than 80 %.
The soil water content was measured by the ovendrying method.Samples were initially taken from each soil layer, every 10 cm, by a soil drill, then dried to constant weight at 105 o C in an air-forced oven and finally weighed and the soil water content calculated.We used the following formula to determine Ci: where, Ra and Trv are the rainfall amount (mm), and the total slope runoff volume (mm), for each rainfall event, respectively.The software packages Matlab 7.10.0 and Origin 8.0 were used for the statistical analysis and to draw figures, respectively.

Introduction to BP neural networks and grey relational analysis 1. BP Neural networks and their development
Artificial neural networks appeared around 1940s, and have been widely used due to their high-dimensional and parallel distributed processing, for being self-adaptive, self-organizing, self-learning, and because of other excellent features (Chen, 2000).An integrated neural network should include three parts: an input layer for to receive input signals, an output layer for exporting and a hidden layer for converting and processing (Figure 1).
A neural network, using the back-propagation algorithm has outstanding characteristics, e.g., various types of function approximation (especially for nonlinear functions), pattern recognition, data classification, compression, etc (Song et al., 2008).The steps underlying the algorithm are: learning information is propagated forward, while error is propagated backward, and network weights and error are continuously adjusted for repeated training.Finally it was verified whether the output value approached the expected value as close as possible.The network training can be interrupted and resumed, until the output error is less than the specification error or the iteration loop reaches the desired maximum value (Han, 2006).
The major drawback of the fundamental BP algorithm is the low training speed.Therefore, in this paper, a momentum term was added to the basic BP algorithm.It reflected the previous experience of the network and had a damping effect on the subsequent adjustment, and effectively avoided the local extremum (Han, 2002;Feng et al., 2007).Using W and X to represent the weight matrix and input vector, respectively, the weight adjustment formula including the momentum term can be expressed as follows: (3)  where α ∈ (0, 1) is the momentum factor and α = 0.95 in this paper (Cong, 1998); δ is the curved surface of the error gradient; η is the learning rate; t n is the training time.
Moreover, the learning rate of the standard BP algorithm is a constant, which would increase iterative loops, as function calculation would skip the narrow gap when the learning rate is too high.
However, thel adaptive learning rate is selfadjusting, and consequently the decreased speed of the network error will ensure a better learning rate of the network to continue training work and maintain the convergence speed high and shorten the training time.The formula of the adaptive learning rate is as follows: where both β and θ are adjusted ratios of the learning rate, and β > 1, θ < 1, generally β = 1.05, θ = 0.7 (Yang, 2001); E RME is the network

Grey relational analysis
Grey relational analysis is an important kind of multivariate analysis, based on the Grey system theory, which assumes that a random process is a grey quantity variable in an area of a certain amplitude and time zone, as proposed by Deng (1989).It can measure the correlation between series and belongs to the category of data-analysis or geometric methods.Grey relational analysis is based on geometrical mathematics, in accordance with the principles of normality, symmetry, entirety, and proximity.The purpose of grey relational analysis is to search for primary relationships among the factors and to determine the relevant ones that significantly influence certain objectives (Fu et al., 2011).Grey relational analysis is defined as a quantitative analysis to detect a trend in the varietal system, and the more similar the developing trends, the greater is the relational extent.The way to compare the relational extent among factors is called relational coefficient or relational grade method.The measure of relevance between two systems or two factors is known as the relational grade, which describes the trend of relationship between an objective sequence (collection of measurements) and a reference sequence in the system.If both tend towards coincidence, the relational grade is scored high, and vice versa.Grey relational analysis is widely used, since it enriches the traditional mathematical statistical methods and compensates their deficiencies (Liu et al., 1999).
The objective sequence Y 0 (p)=(Y 1 (1),Y 1 (2),…,Y 1 (M)), and the reference sequence X q (p)=(X 1 (1),X 1 (2),…, X 1 (M)), q=1, 2, 3,…, N, have the sequence length N. The research target is therefore to determine the closeness, i.e., to determine the relational grade between X q (p) and Y 0 (p).Similarly to the common statistical analysis methods, grey analysis firstly requires an appropriate normalization of raw data to remove anomalies associated with different measurement units and scales.The raw data can be transformed into dimensionless forms by standard processing, using formula (5): X q (p)=(x q (p)-avx q )/σ q (5) where avx q and σ q are the average value and standard deviation of the sequence x q (p), respectively.Standard processing eliminated not only the influence of the dimension and the order of magnitude, but also diminished the variation degree of the index.
According to the grey relational space definition, the relational grade between x q and y 0 was defined by the following formula: where p=1, 2, 3,…, N. On average, the numerical values of the relational grades ranged from 0 to 1.A marked influence is indicated by r > 0.9, a relatively marked influence by r > 0.8, a detectable influence by r > 0.7, and a negligible influence by r < 0.6.

BP Neural network model of Ci
The BP model of Ci with six input parameters was established using the neural network toolbox function in the Matlab 7.10.0program, and then trained.The number of neurons in the input and output layers are easily determined (6 and 1, respectively).Generally, one hidden layer is enough for the input-output mapping.However, in order to describe the complicated nonlinear soil water movement process accurately, two hidden layers were used, an arrangement that seems to improve the learning convergence speed and mapping relations.After several adjustments and trainings for the hidden layers, it was found that 15 neurons with logarithmic sigmoid transfer-function in the first hidden layer and one neuron with linear transferfunction in the second hidden layer were most efficient and effective.In other words, the network structure was designed (6×15×1×1) for the learning scheme.
Wang Juan et al.Ci = 0.410×Sg -0.058 × (Vt+1/10000) -0.048 × (Vc+1/10000) 0.086 ×Ri 0.983 ×Rd 1.145 ×Aw -0.016 R 2 = 0.9550, F=717.2 > F 0.01 (4, 43) = 3.803 (11) The relative errors of the BP model were insignificantly small, ranging from -10% to 10% and the average of the absolute value of the relative error was only 3.99 %.However, the MSE and average relative error of the multiple nonlinear regression equation were 0.1969 and 24.65 %, respectively.Its fluctuation of relative error was much greater than that of the BP model.This results from figure 2 evidently indicated that the improved BP model could more accurately simulate and map Ci than the multiple nonlinear regression equation.Obviously, the BP model has more advantages than the regression equation to calculate Ci and exactly represent a process of soil water movement.

Sensitivity analysis for factors: based on the BP model
It is difficult to establish a relationship between soil infiltration and impact factors based on a large number of influencing factors.However, it is necessary to find the main factor influencing soil infiltration.Based on the above-established BP model, the default factor test of the BP model was described in this paper.Through a default factor test for each input factor, the sensitivity of each input factor was determined in relation to the output factor, according to the ratio of the test error of each factor by that of the full-factor model.
where R i is the sensitivity index of the i th factor; RMSE i is the model testing error of the i th default input factor; RMSE is the model testing error of the full input factor.If RMSE i >RMSE j , the i th input factor is more sensitive and important than the j th input factor.Besides, if R i >1, the i th input factor is the sensitive factor.All 75 groups of measured data constituted the sample for the sensitivity analysis, and then the default factor test for input layer factors was performed one by one, a 5-factor model of Ci was established using the improved BP neural network algorithm (results in Table 3).The sensitivity extent of six factors on Ci was: Rd>Vc>Vt>Ri> Aw>Sg.The sensitivity index of Rd was higher than 4.3, and indicated that Rd was the most important factor.Vc index followed the Rd index, and the lowest index was Sg (only 1.2), which revealed that Sg was not significant for Ci , compared to the other five factors.

Grey relational analysis for factors
Our second research objective was the objective sequence; knowing that Ci can accurately represent soil permeability in a rainfall event, the factors Sg, Vt, Vc, Ri, Rd and Aw were regarded as reference sequences with impact on the Ci behavior.The Of all 75 experimental data groups, 50 groups were randomly selected as the samples to train and build the BP network, while the other 25 data groups were used for testing.The first step was the standard data preprocessing, according to equation ( 8): Where X i is the i th input data after standard preprocessing; x i is the i th raw data prior to standard processing; x max and x min are the maximum and minimum data.This processing method decreases the sample distribution range, which is between [0, 1], and is more appropriate for BP neural network training (Zhang et al., 2007).
After 5000 training cycles for the BP network, the MSE (mean square error) of the test process was 0.1003, and the best-fitting equation between measured and simulated values was: y simulated = 0.9955 x measured , R 2 = 0.9872, F = 690.1,p = 0.000.
In addition, the same 50 data groups were used to fit a multiple nonlinear regression equation to Ci and six factors.The other 25 data groups were then used for testing.Considering the factual causation between six impact factors and Ci, a typical multiplication event was formed by Ci and impact factors, according to the basic principles of probability, which was expressed as: where a, b, c, d, e, f, and g are undetermined parameters.Due to several zero values of Vt or Vc, these cannot be used as base variables in the power function.So formula (9) was changed into formula (10), as follows: Ci = a×Sg b × (Vt+1/10000) c ×(Vc+1/10000) d ×Ri e ×Rd f ×Aw g (10) The specific forms of the multiple nonlinear regression equation for Ci were

^
objective and reference sequence length was 75 and the data were the same as those used in the sensitivity analysis.In a first step, the raw data were preprocessed by standard processing with equation (8); secondly, the grey relational coefficient was calculated by equation ( 6); then grey relational grades were calculated by equation ( 7); and at last, arranged in relational grade order (Table 4).
The grey relation of factors and Ci were: Rd>Vc>Vt>Ri> Aw>Sg, in agreement with the result of the default factor test of the BP neural network.Grey relation of Rd was highest (0.8169), indicating that Rd had a very important influence on cumulative infiltration.The grey relation of Vc, Vt and Ri ranged from 0.7 to 0.8, while that of Aw and Sg was relatively small (only about 0.7).

DISCUSSION AND CONCLUSIONS
The Ci can be expressed by the formula: , where i t is the instantaneous infiltration rate, and t 1 and t 2 are two time moments of infiltration.From the above formula, it was easy to find that time is the most important and direct factor affecting soil infiltration, and the infiltration amount increased with increasing time periods.However, Vc, Vt, Ri, Aw, Sg and other factors indirectly affect soil infiltration.The Vc is relevant for reducing the damage of kinetic energy of raindrops on the soil surface structure and to avoid crust formation, by maintaining the soil infiltration rate high for a long period (Huang et al., 2010).In addition, plant roots could form pore structures in the soil that differ according to the plant type, which can influence soil infiltration.A different Vt is also a significant factor for Ci (Zhao et al., 2001;Li & Shao, 2007).There is no consistent information in the literature about the effect of Sg on soil infiltration, and most authors agreed that the effect of Sg on Ci is not significant (Singer & Blackard, 1982;Mah et al., 1992;Cerdà & García-Fayos, 1997).The default factor test of the BP neural network and grey relational analysis could be used to determine the correlative extent of factors and Ci, such as Rd, Ri, Aw, Vc, Vt, and Sg.The ranking of results by the two methods was the same: Rd>Vc>Vt>Ri> Aw>Sg.This not only demonstrates the reliability and accuracy of the default factor test for BP neural networks, but also completes the selfvalidation of grey relational analysis.
This experiment was conducted using disturbed soil under specific conditions, which is a qualitative simulation study.Therefore, it may be difficult to transfer the results of this study to the field, due to the effect the specific experimental conditions will have on the results.Nevertheless, the results of this study could also be used for comparative purposes and provide some useful information for a better understanding of the different influence of impact factors on cumulative infiltration in a small scale slope.

Figure 1 .
Figure 1.Structured chart of neural network.

Figure 2 .
Figure 2. Comparison of Ci between a six factor BP model and multiple nonlinear regression equation.
Water Conservation Research Institute, of the Chinese Academy, at the field monitoring station of rainfall runoff regulation in Linghou, Wuquan, District of Yangling, Province Shaanxi.Yangling is located on the western Guanzhong Plain of Shaanxi, north of the Weihe River (107 o 59' ~ 108 o 08' E, 34 o 14' ~ 34 o 20' N; 418.0 ~ 540.1 m asl).