Modeling of trees failure under windstorm in harvested Hyrcanian forests using machine learning techniques
Site selection
Hyrcanian temperate forests, dominated by old broadleaf trees, are located in the north of Iran adjacent to the Caspian Sea. We selected the Neka Zalemroud forest in Mazandaran province as the study area for this research (36° 26′ 09″ to 36° 30′ 47″ N latitude and 53° 20′ 34″ to 53° 31′ 51″ E longitude) as the boundaries of sampling area have been illustrated in Fig. 1 by the authors. This forest has been covered with 2533 hectares of old broadleaf trees such as Fagus orientalis, Carpinus betulus, Quercus castanafolia, Acer velutinum, Acer cappadocicum, Parrotia persica and some other species. This region is affected by permanent winds and the maximum wind speed is in the range of 10 to > 30 m/s. The windstorms with more than 100 km/h cause damage forest trees, including uprooting or stem breakage. The windstorms cause many tree failures annually, which is in result of the tree harvesting and gap creation in forest stands. Hence, we aimed to identify the uprooted or stem broken trees using the sample plots data after the windstorm with 100 km/h (4 h) on 21 March 2018.
Figure 1
The location of study area and sample plots (QGIS 3.12.0, https://www.qgis.org/en/site/).
Full size image
Methods
Data from field measurement of 600 sample plots in the study area were summarized before and after the windstorm in March 2018. ARC MAP 9.3 software was used to spatially locate the sample plots on the forest map. Also, the slope of plots was defined by topography map (1:50,000) in this software. The permanent sample plots were created by Department of Natural Resources of Mazandaran Province (in the structure of the forest management plan) as a part of government-sponsored “permanent sample plot” protocol to monitor the changes in the number of trees and species48. In this protocol, the centers of these sample plots are marked with a metal rod and the geographical coordinates are recorded. The attributes of all trees in the sample plots are recorded consisting of species, diameter at the breast height, tree’s height, number of trees, distance to the center of the plot, crown diameter and land, soil and climate characteristics of the plot location. On average, in each plot, there are about 10 to 20 trees with a diameter of more than 7 cm (marketable size), and using the recorded data of the trees, it is possible to identify each tree in subsequent monitoring. As the protocol dictates, preliminary data were recorded on the stage of creating plots by experts from the Provincial Department of Natural Resources in 201649. After the windstorm on 21 March 2018, these plots were investigated by the authors to find damaged and undamaged trees. The area of each sample plot was 1000 m2 and in the shape of the circle with 17.84 m radius. These plots were investigated to record harmed trees which were uprooted or stem broken. Sample plots are clustered into susceptible (with harmed trees) and unsusceptible (without harmed trees) plots. We recorded the forest stand attributes, as well as tree attributes, which influence the likelihood of trees to be damaged in windstorms. Therefore, we investigate the damaged and resistant trees in the plots within 1000 m2 area. The stand variables were recorded in the plot area such as the mean of the trees’ height in the plot. Some plots (306 plots) were determined where the damaged trees were observed and some others (294 plots) were determined where the target tree in the middle of the plot was undisturbed. In a literature review, we found that stand and tree characteristics influence wind-susceptibility of trees (such as tree diameter and height, spread crown area, rooting system, soil type, stand density, topography, etc.1,3,19,20,21). Therefore, 15 variables in 600 sample plots were recorded that are divided into two categories: 1. Forest stand variables: Plot Slope (PS) (%), Soil Depth (SD) in plot area (cm), trees Mean Diameter at the Breast Height (MDBH) at plot (cm), trees Mean Height (MH) in plot (m), trees Density (De) in plot (Number of trees), trees Diversity (Di) in plot (Number of tree species), Number of Thicker trees in diameter (than the target tree) (NTh), and Number of Taller trees (than the target tree) (NTa).
2.Tree variables (bigger than marketable size): Tree Area (TA) (occupied area by the tree) (m2), Tree Diameter at the Breast Height (TDBH) (cm), Tree Height (TH) (m), Tree Crown Diameter (TCD) (m), Mean Distance from Neighbor trees (MeDN) (m), Minimum Distance from Neighbor trees (MiDN) (m), and Maximum Distance from Neighbor trees (MaDN) (m). Tree heights and some more data were recorded in permanent sample plots before the windstorm. Therefore, we used the recorded data at the stage of creating plots, for damaged trees. In fact, 10 to 20 trees are recorded in each plot and using the recorded data of the trees (species, diameter, distance to the center of the plot, etc.) it is possible to identify target trees (damaged or undamaged) in subsequent monitoring. There are some other factors that influence wind-susceptibility of trees such as forest edge, tree diversity, and land form. In this research, we were looking for the impact of forest plan activities on tree failure; therefore, we neutralized forest edge effects by selecting sample plots inside the forest stands which are far from the forest edges. In land form variables, land slope was considered in stand variables, but altitude and geographical aspect of the hill, as well as tree diversity, were omitted because of limited variation in the samples.
Tree Failure Model (TFM) was developed by recording 15 variables of the trees and forest stands in 600 selected sample plots. In fact, we designed three TFMs with three modeling techniques to achieve the most accurate one based on model accuracy assessment. The damaged trees in the plots were identified by two features: (1) The tree was uprooted by wind forces along the windstorm. (2) The tree was uprooted or stem broken. Leaning trees (under windstorm force) were also counted as uprooted trees. The damaged tree was chosen as the target tree in susceptible plots (plots with a damaged tree). Also, the central tree was selected as the target tree in the wind-unsusceptible plots (plots without damaged tree). Indeed, the output of TFM will be in two classes of damaged trees (1) and stable trees (0). Hence, the response or output of the model will be a discrete class {0,1}. The accuracy of the model is assessed by confusion matrices which detect the number of accurate and false classifications of sample trees.
The new mathematical modeling approaches and machine learning techniques are needed to cover limitation in data collection and forest inventories. Indeed, ANNs are over 50 years old, but just not often applied yet. History of ANN development represents this fact that learning algorithms and structure of neural networks are developed every year. Therefore, we developed artificial intelligence techniques in natural phenomena modeling11,17 namely MLP, RBFNN, and SVM.
Machine learning techniques rely on a specific concept that is “a set of weak learners develop a single strong learner” (by Freund and Schapire22, Breiman23 and Breiman et al.24). As it is known, a weak learner is a classifier correlating slightly with the real (target) classification; while a strong learner is well correlated with real (target) classification. Machine learning algorithms are trying to combine weak classification rules in a one strong classification rule. Based on this approach, we used 15 variables (even if we believe that these variables are not related to tree susceptibility in the wind) and tested some algorithms and variable weights to make weak learners or rules and combining them into one strong rule in the structure of three machine learning techniques.
Multi-layer perceptron (MLP) neural network
ANNs use different methods, such as feed-forward, backward, recurrent and other, to teach the network for output prediction. MLP is a multi-layer form of Feed-forward neural networks without any cycle or loop. In Feed-forward neural networks, the information analysis is performed in one direction from the input layer, through the hidden layer to the output layer47. In this learning method, the errors of the network propagate from the output layer to inputs to revise the weights of input variables. MLP is a multi-layer Artificial Neural Network (ANN) model with self-learning mechanism which uses samples for classification. Indeed, MLP has been using some interconnected processing elements that are called PEs (Processing Elements). MLP learns by using samples and transfer functions which are applied between neurons and hidden layers in a computer program17. In the training process, each PE receives signals periodically from other PEs and sends the new signal to other processors. Considering inputs, MLP adjusts the weights of neurons continuously, and the learning process is completed.
We used some activation functions (such as logarithmic sigmoid, hyperbolic tangent, and linear transfer functions) which determine the relation between inputs and outputs and these functions were tested to achieve the best performance of MLP. Back Propagation (BP) method propagates the error of outputs to the input layer where the first random weights have been assigned. The weights of the network inputs will be justified until the best performance of the network is reached; and after that, the learning process will be completed7,11. Errors between Ynet (MLP output) and Y (real class of tree failure) are decreased by BP when the weight of neurons or Processing Elements (PEs) (w) and input variables (x) come to the best performance, and the output of jth PE on the kth layer (PEkj) will be achieved by Eq. (1):
$$net_{j}^{k} = mathop sum limits_{i = 0}^{n} w_{ji} x_{ji}$$
(1)
Transfer functions are used in the structure of network, and neuron output value is determined by (Eq. 2).
$$Y_{net} = smallint net_{j}$$
(2)
Finally, weights of t samples will be adjusted by delta rule which has been summarized in Eq. (3).
$$w_{ji}^{t} = w_{ji}^{t – 1} + Delta w_{ji}^{t}$$
(3)
By using the ANN function in MATLAB R2013b, 360 uniformly distributed random samples (60% of 600 samples) were defined as training data set. 120 evenly distributed random samples (20% of 600 samples) were defined as validation data set, and 120 samples (20% of 600 samples) were determined as test data set. All data were normalized to the interval of 0 to 1 using the Min–Max technique by mapminmax function in MATLAB R2013b (Refer to Demuth and Beale25 for MATLAB codes for MLP neural network development and related preprocessing algorithms).
Radial basis function neural network (RBFNN)
Radial basis function neural network is architecturally similar to the MLP with different activation function in the hidden layer. RBFNNs have been used in function approximation and classification in researches of the last decade13,26,27,28. RBFNN uses samples in two data sets of training and test. The radial function is applied in each neuron of the hidden layer; and the number of neurons depends on input matrix of variables. Considering two classes of trees failure (0 and 1) in this research, we have two output layers in the structure of RBFNN. Gaussian function is the most frequently used function in the hidden layers of RBFNN13,27. The Gaussian function can find the center of circular classifiers successfully. The Gaussian function regulates the centre of mentioned circular classifiers by Eq. (4).
$$R_{j} left( x right) = expleft( {frac{{||x – a_{j} ||^{2} }}{{2sigma^{2} }}} right)$$
(4)
In Eq. (4), input variables are structured in the matrix “x”, radial basis function has been defined as Rj(x), centre of RBF function is presented as aj, and we have a positive real number as “ϭ”. The outputs of network will be calculated by an output function Eq. (5).
$$y_{k} = mathop sum limits_{j = 1}^{m} w_{jk} R_{j} left( x right) + b_{j}$$
(5)
In Eq. (5), the number of calculation nodes in the structure of hidden layers (j), the number of neurons (m), the weights of neurons (wik), and a bias value (bj) have been used to calculate output (yk).
Neuron weights (wjk) are updated continuously to decrease output errors until network training process comes to end. Network performance is calculated when the number of neurons and the weights of neuron or layers are fixed13. (Refer to Demuth and Beale25 for MATLAB codes for RBF neural network development and related preprocessing algorithms).
Support vector machine (SVM)
SVM is one of the machine learning techniques that requires quite a lot of data for training, but this method also provides more accurate results than other methods when the volume of training data is limited29,30. Therefore, SVM has been used for modeling in this paper to deal with this issue with the collected data along forest inventory.
As a classifier technique, SVM aims to determine the largest margin in decision boundaries that could separate classes of decision31. SVM is looking for the largest margin in the boundaries of classification when the uncertainties in the decision are expected27. This method of prediction minimizes the probability of over-fitting in classes limits of tree failure.
We have two datasets of training and test in the structure of SVM. The values of target are structured in a n-dimensional matrix so it is possible to find the most accurate boundaries and margins. Equation (6) is the SVM model and equation parameters define: y(x) = SVM output, α_i = a multiplier, K = kernel function, and b = threshold parameter.
$$yleft( x right) = mathop sum limits_{i = 1}^{n} alpha_{i} Kleft( {x_{i} ,x_{j} } right) + b$$
(6)
Then we defined Gaussian Radial Basis Function (RBF) in Eq. (7), in the context of non-linear SVM. As we know, RBF is the most popular function in the context of SVM with remarkable ability to control generalization of SVM classifier.
$$Kleft( {x_{i} ,x_{j} } right) = {text{exp}}left( { – gamma| x_{i} – x_{j}|^{2} } right)$$
(7)
The parameters of Eq. (7) are: xi and xj = samples and γ = kernel parameter.
Finally, primal problem in Eq. (8) should be minimized to achieve the most accurate SVM for tree failure prediction.
$$frac{1}{2}|w|^{2} + Cmathop sum limits_{i = 1}^{n} xi_{i}$$
(8)
In Eq. (8), the parameters are: 1/2||w||2 = the margin, Σξi = training errors and C = the tuning parameter.
SVM uses samples in two data sets of training and test. The classes of the target will be summarized in a n-dimensional matrix to determine the nearest classification boundaries and margins. (Refer to Demuth and Beale25 for MATLAB codes for SVM neural network development and related preprocessing algorithms).
Model selection
MLP, RBFNN, and SVM models were run on the training dataset with 15 tree variables as inputs and tree failure classes in 600 selected trees as output. To evaluate prediction accuracy of the MLP, RBFNN, and SVM models, we used the confusion matrix to determine the percentage of accuracy in target tree classification. Also, the number of trees with accurate and failed classification will be detected11.
Sensitivity analysis
Sensitivity analysis was designed to prioritize the most accurate model variables with respect to the significance of variables in output. Sensitivity analysis defines the usefulness of variables in model predictions. In sensitivity analysis, we changed each variable in the range of standard deviation with 50 steps while the other variables were fixed at the value of the average. Then, the standard deviation of outputs for each variable changes was measured as model sensitivity for that variable. Variables with high value in the outputs standard deviation are the most important variables with more influence on model outputs. The trend of model output changes with changing the most significant variables, in the range of standard deviation (50 steps) was illustrated in some figures to find out the way that model outputs are changing with variable changes (negatively or positively) (Refer to Kalantary et al.27 and Jahani et al.13.
Environmental decision support system (EDSS) tool
Finally a user friendly GUI (Graphical User Interface) tool was designed as an EDSS for susceptible tree identification in windstorms. It is applicable for forest managers who are looking for hazardous trees to plan for tree protection and increase the forest stands resistant against windstorms. ANN models use a huge matrix of weights; so model execution should be in the mathematical software (in this research MATLAB R2013b). Users, who are not familiar with the software, need a simple tool to run the model on new samples and get the results of prediction. To design EDSS tool, we developed a GUI extension in MATLAB R2013b software. With this tool, users enter the values of trees and stands variables (based on forest inventory data in other target forests) and the susceptible trees will be identified only by pushing a button. The model will be run on the data and the model outputs for each tree will be appeared in a table (0 or 1). More
