Estimation of DBH at Forest Stand Level Based on Multi-Parameters and Generalized Regression Neural Network

The diameter at breast height (DBH) is an important factor used to estimate important forestry indices like forest growing stock, basal area, biomass, and carbon stock. The traditional DBH ground surveys are time-consuming, labor-intensive, and expensive. To reduce the traditional ground surveys, this study focused on the prediction of unknown DBH in forest stands using existing measured data. As a comparison, the tree age was first used as the only independent variable in establishing 13 kinds of empirical models to fit the relationship between the age and DBH of the forest subcompartments and predict DBH growth. Second, the initial independent variables were extended to 19 parameters, including 8 ecological and biological factors and 11 remote sensing factors. By introducing the Spearman correlation analysis, the independent variable parameters were dimension-reduced to satisfy very significant conditions (p ≤ 0.01) and a relatively large correlation coefficient (r ≥ 0.1). Finally, the remaining independent variables were involved in the modeling and prediction of DBH using a multivariate linear regression (MLR) model and generalized regression neural network (GRNN) model. The (root-mean-squared errors) RMSEs of MLR and GRNN were 1.9976 cm and 1.9655 cm, respectively, and the R2 were 0.6459 and 0.6574 respectively, which were much better than the values for the 13 traditional empirical age–DBH models. The use of comprehensive factors is beneficial to improving the prediction accuracy of both the MLR and GRNN models. Regardless of whether remote sensing image factors were included, the experimental results produced by GRNN were better than MLR. By synthetically introducing ecological, biological, and remote sensing factors, GRNN produced the best results with 1.4688 cm in mean absolute error (MAE), 13.78% in MAPE, 1.9655 cm for the RMSE, 0.6574 for the R2, and 0.0810 for the Theil’s inequality coefficient (TIC), respectively. For modeling and prediction based on more complex tree species and a wider range of samples, GRNN is a desirable model with strong generalizability.


Introduction
In forest monitoring, the diameter at breast height (DBH) is the diameter of a cross-section of a tree trunk 1.3 m above the ground. The DBH is an important factor used to estimate important forestry indices like forest growing stock, basal area, biomass, and carbon stock. The traditional DBH ground surveys effectively provide objective and reliable monitoring and managing information for

Research Data
The research data included the inventory data for forest management planning and design, the digital elevation model (DEM), and remote sensing images from Landsat-8 Thematic Mapper (TM). The 2017 inventory data were supplied by the Forestry Bureau of Longquan. The TM and DEM data were both derived from the Geospatial Data Cloud [43].
Six factors were derived from the inventory data for forest management planning and design, including DBH, soil depth, humus thickness, canopy density, number of stems per hectare, and age of tree. Three factors, elevation, slope, and aspect, were obtained from DEM and 11 factors were derived from the remote sensing images, including B2-B7, NDVI, RVI, DVI, EVI, and RI.
The inventory data, containing 58,910 records of forest subcompartments, were obtained from the Forestry Bureau of Longquan in 2017. According to the rules for the technical operation of the inventory for forest management planning and design in Zhejiang province (2014) [44], the original inventory data acquisition procedure is as follows: Within the scope of the subcompartments, some standard plots or angle gauge plots are laid out using random, mechanical, or other sampling methods, and the corresponding forest factors are measured in the plots by professionals, and the subcompartment factors are then calculated. Therefore, whether modeling or predicting, all the

Research Data
The research data included the inventory data for forest management planning and design, the digital elevation model (DEM), and remote sensing images from Landsat-8 Thematic Mapper (TM). The 2017 inventory data were supplied by the Forestry Bureau of Longquan. The TM and DEM data were both derived from the Geospatial Data Cloud [43].
Six factors were derived from the inventory data for forest management planning and design, including DBH, soil depth, humus thickness, canopy density, number of stems per hectare, and age of tree. Three factors, elevation, slope, and aspect, were obtained from DEM and 11 factors were derived from the remote sensing images, including B2-B7, NDVI, RVI, DVI, EVI, and RI.
The inventory data, containing 58,910 records of forest subcompartments, were obtained from the Forestry Bureau of Longquan in 2017. According to the rules for the technical operation of the inventory for forest management planning and design in Zhejiang province (2014) [44], the original inventory data acquisition procedure is as follows: Within the scope of the subcompartments, some standard plots or angle gauge plots are laid out using random, mechanical, or other sampling methods, and the corresponding forest factors are measured in the plots by professionals, and the subcompartment factors are then calculated. Therefore, whether modeling or predicting, all the data-including DBH, Forests 2019, 10, 778 5 of 18 soil depth, humus thickness, canopy density, number of stems per hectare, and age of tree-are presented as an average value at the forest subcompartment level. The unsuitable data were removed via preliminary processing, including records for zero-volume subcompartments which were located in non-forest plots, or those with subcompartment DBHs less than 5 cm, which do not meet the standard for measurement [44]. Therefore, the final inventory data included 35,763 records of forest subcompartments. These were randomly separated into two groups, with one containing 25,000 records for modeling, and the remaining 10,763 records used as testing data for prediction (Table 1). The DEM data were originally obtained from the first version data of the global DEM (GDEM) generated by the advanced spaceborne thermal emission and reflection radiometer, including four images in 2009, with 30 m resolution, IMG data type, World Geodetic System 1984 (WGS84), and the universal transverse mercator (UTM) projection. The data were obtained from the International Scientific and Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences [43].
The remote sensing images, which correspond fully with the timeline of the inventory data, were generated by Landsat-8 TM from Landsat satellites in 2017. There are two sensors in each Landsat satellite that collect data from discrete portions of the electromagnetic spectrum, known as bands, according to their wavelength. The Landsat-8 satellite's images, with the UTM/WGS 84 coordinate system, include 11 bands representing portions of the electromagnetic spectrum, consisting of visible bands and invisible bands to the human eye. Here, the visible bands include 5 bands, i.e., band 1, called the coastal band, with wavelengths of 0.433-0.453 µm; band 2, called the blue band, with wavelengths of 0.450-0.515 µm; band 3, called the green band, with wavelengths of 0.525-0.600 µm; band 4, called the red band, with wavelengths of 0.630-0.680 µm; and band 8, called the pan band, with wavelengths of 0.500-0.680 µm. The invisible bands include 6 bands, i.e., band 5, called the near infrared (NIR) band, with wavelengths from 0.845 to 0.885 µm; band 6, called the SWIR-1 band, with wavelengths of 1.560-1.660 µm; band 7, called the SWIR-2 band, with wavelengths of 2.100-2.300 µm; band 9, called the cirrus band, with wavelengths of 1.360-1.390 µm; band 10, called the TIRS-1 band, with wavelengths of 10.6-11.2 µm; and band 11, called the TIRS-2 band, with wavelengths of 11.5-12.5 µm. The spatial resolution of the bands is mostly 30 m except for band 8, which contains panchromatic data with a resolution of 15 m, and bands 10 and 11 provide data collected at a 100 m resolution.

Data Preprocessing and Integration
The geometric correction, radiometric calibration, and atmospheric correction of Landsat-8 TM images were handled by ENVI 5.3 (Exelis Visual Information Solutions, Boulder, CO, USA). Afterwards, VI was calculated from the processed images by the band math tool in ENVI 5.3. In ArcGIS 10.2 (Environmental Systems Research Institute, Redlands, CA, USA), the DEM data and remote sensing data were extracted based on each subcompartment. For modeling and prediction, all the preprocessed data, including the inventory data for forest management planning and design, DEM data, and remote sensing images from TM, were integrated into the same relational database.

Empirical Model
The traditional models used to predict DBH of the forest subcompartment are usually based on the relationship between DBH and the independent variables of tree age, height, and crown width. However, the costs of data-acquisition of the height and crown width of the trees are quite high in practical measurement. In contrast, the tree-age of the subcompartment refers to the average tree-age of the dominant species which is usually stable in a relative short period. In this way, as long as the initial values of the tree age are obtained successfully, only 1 needs to be added for each year after the measurement. Thus, the cost of data-acquisition is relatively low. Therefore, the 13 empirical models based on the age-DBH model that use tree age as the only independent variable [45][46][47][48] were selected to predict DBH in this paper ( Table 2).

Model Name Equation
Linear model y = a + bx Logarithmic curve model Note: y denotes the DBH of the forest subcompartments, x is the age of tree; a, b, c, and d are the parameters to be estimated.

Correlation Analysis
The nonparametric test is an important part of statistical analysis methods that is used to infer the population distribution using the sample data in situations where the total variance is unknown or uncertain. As the total variance in advance was unable to be estimated, the nonparametric test was preferable for inferring the population distribution. The Spearman rank correlation analysis is a specific method used to implement the nonparametric test to estimate the correlation between two variables. Here, the observed values of two variables needed to be paired with the graded data originally or converted from the observed data of the continuous variables without considering the overall distribution pattern and sample size of the two variables. The formula for calculating the correlation coefficient of the Spearman rank is shown in Equation (1), and the calculation is based on the rank of the value, rather than the value itself. The calculation was conducted in SPSS 20 (International Business Machines Corporation, Armonk, NY, USA) using the double tail test.
where R i is the rank of the ith observed value for variable x; S i is the rank of the ith observed value for variable y; R and S are the average ranks of variables x and y, respectively; and n denotes the total number of observed values.

Multiple Linear Regression Model
The multiple linear regression (MLR) is a statistical method. The MLR can be used to predict the values of a response variable using several explanatory variables [49] and fitting the relationship between the independent variables and dependent variable. It has been represented in Equation (2), and was implemented in SPSS 20.
where y denotes the dependent variable, xi denotes the ith independent variable, β 0 represents a constant, β i denotes the slope coefficient for the ith explanatory variable, and ε represents the residual. The process for solving the parameters for Equation (2) was as follows. First, leaving aside the remote sensing parameters, the four factors (elevation, age of tree, canopy density, and the number of trees per hectare) that were involved in the MLR analysis were obtained from the inventory data for forest management planning and design or DEM data. The various MLR selection methods include the forward, backward, stepwise, entering, and deletion methods. This study chose the stepwise method, also called the stepwise regression model, to filter the independent variables to feed to the MLR model, which is represented in Equation (3). This was implemented in SPSS 20.
where X 1 is the tree age, X 2 is the number of trees per hectare, X 3 denotes elevation, and X 4 is canopy density. Subsequently, the seven remote sensing image factors (bands 2-7 and DVI) together with the previous four factors (elevation, age of tree, canopy density, and number of trees per hectare), for a total of 11 factors, were introduced into the MLR model to analyze the linear relationship among DBH and the parameters. This study also chose the stepwise regression method to fit the MLR model in SPSS 20, and the fitting equation was established as: where the variables from X 1 to X 4 denote tree age, number of trees per hectare, elevation, and canopy density, respectively; X 5 denotes DVI; X 6 to X 8 denote bands 2 through 4, respectively; and X 9 and X 10 are bands 6 and 7, respectively. Note that band 5 was removed after the execution of the stepwise regression model.

Generalized Regression Neural Network
The generalized regression neural network (GRNN) was first proposed by Specht in 1991 [50]. On the basis of the estimation of a probability distribution function, the GRNN can solve nonlinear approximative problems. The GRNN gradually converges to the optimal regression surface as the number of training samples increases. The advantages of the GRNN include its simplicity and high accuracy in modeling and prediction. The GRNN architecture includes four layers ( Figure 2): Input, hidden, summation, and output [51].
On the basis of the estimation of a probability distribution function, the GRNN can solve nonlinear approximative problems. The GRNN gradually converges to the optimal regression surface as the number of training samples increases. The advantages of the GRNN include its simplicity and high accuracy in modeling and prediction. The GRNN architecture includes four layers ( Figure 2): Input, hidden, summation, and output [51]. In the input layer, the amount of neurons is determined by the number of the input variables, and each neuron represents an input variable. The input layer does not process the signal but nonlinearly transmits the signal to the hidden layer.
In the hidden layer, the neurons accept the information from the input layer and express the relationship (Pi) between the input layer and hidden layer with a pattern of radial basic function as the activation function. The radial basic function is usually represented by the Gaussian function: where i = 1, 2, …, n, with n representing the number of training cases; Pi denotes the multivariate Gaussian function; σ denotes the smoothing parameter, which is also called the spread parameter; X represents the input variable; and Xi denotes a training sample for the ith neuron of the hidden layer.
There are two summations in the summation layer, Sw and Ss, which are defined by Equations (6) and (7), respectively. Sw is the weighted summation of the summation layer and Ss is the arithmetic summation of the summation layer, where wi represents the weight of the ith neuron of the hidden layer to connect to the summation layer. In the input layer, the amount of neurons is determined by the number of the input variables, and each neuron represents an input variable. The input layer does not process the signal but nonlinearly transmits the signal to the hidden layer.
In the hidden layer, the neurons accept the information from the input layer and express the relationship (P i ) between the input layer and hidden layer with a pattern of radial basic function as the activation function. The radial basic function is usually represented by the Gaussian function: where I = 1, 2, . . . , n, with n representing the number of training cases; P i denotes the multivariate Gaussian function; σ denotes the smoothing parameter, which is also called the spread parameter; X represents the input variable; and X i denotes a training sample for the ith neuron of the hidden layer. There are two summations in the summation layer, S w and S s , which are defined by Equations (6) and (7), respectively. S w is the weighted summation of the summation layer and S s is the arithmetic summation of the summation layer, where w i represents the weight of the ith neuron of the hidden layer to connect to the summation layer.
In the output layer, the output variable y is calculated using Equation (8). The dimension of the output matrix of the final training of the neural network depends on the number of neurons.
Forests 2019, 10, 778 9 of 18 The GRNN model was implemented in MATLAB 2017b (The MathWorks, Natick, MA, USA). Its network structure was determined by: where x denotes the input matrix for the independent variables, y denotes the output matrix for the dependent variables, and spread is the smoothing factor. The value of this factor was tested by cross-validation.

Model Performance Metrics
The following indices can be used to assess the performance of each model, where Y i represents the measured values, y i represents the predicted values for sample i, Y denotes the predicted mean, y denotes the observed mean, and n denotes the total number of samples.
The first statistical metric is the mean absolute error (MAE), expressed by The mean absolute percentage error (MAPE) is represented by The root-mean-squared error (RMSE), one of the most commonly used indicators for evaluating performance, is expressed by The coefficient of determination (R 2 ) is a statistical index used to describe the degree to which a change in the results can be explained by a change in the independent variables. R 2 ranges from 0 to 1; the closer R 2 is to 1, the better the model effect. The values of R 2 are calculated by: The Theil's inequality coefficient (TIC) is a measure of the difference between the measured and predicted values. When the value is between 0 and 1, and the closer the value is to 0, the better the prediction effect. The TIC represented by:

Modeling and Prediction Based on Empirical Models
The data from 25,000 samples were fed to the 13 empirical models. The results are illustrated in Table 3. The modeling results of the 13 empirical models showed that the R 2 values were between 0.438 and 0.592. The F-test was performed for the first 13 models, showing that the p-values (0.000) were less than 0.01. This indicates a very high level of significance for the models.
Subsequently, the tree age data from the remaining 10,763 samples were fed to the above 13 empirical models to predict DBHs of the forest subcompartments. The results are illustrated in Table 4 and Figure 3. The MAPEs ranged from 15.63 to 19.30%, which means that the accuracies of the prediction were above 80%. Compared to the modeling results, the fluctuation in R 2 obviously increased with a range of 0.2492-0.5092 in the prediction results. Even if the four models (linear, compound curve, growth curve, and exponential curve) produced the same R 2 of 0.455 in modeling, the prediction ability of the linear model, with an R 2 of 0.4374, was obviously better than the curve models (compound, growth curve, and exponential curve) with an R 2 of 0.2492. The Gompertz model provided the strongest generalizability with the best R 2 of 0.5092, even if the modeling R 2 was only 0.496.

Correlation Analysis
The significance level was used to compute the confidence level with a value between 0 and 1, usually denoted by symbol p. Due to the rule that the lower the significance level, the higher the confidence level, and the lower the probability of making mistakes, choosing an appropriate significance level is necessary to guarantee the proper credibility level. A two-tailed test was used to compute the significance level (p) here, and the correlation analysis was involved in the calculation of the correlation coefficient (r).
Since forest growth is affected by many environmental factors, such as topography and soil [52,53], the DBH of trees of the same age may vary considerably in different environments. From the inventory data for forest management planning and design and DEM, this study initially chose eight influencing factors-elevation, slope, slope aspect, thickness of soil layer, humus thickness, age of tree, canopy density, and number of trees per hectare-to analyze the correlation coefficients [54][55][56]. The results ( Table 5) show that the six factors (elevation, slope aspect, humus thickness, age of tree, canopy density, and number of trees per hectare) reached a very significant level (p ≤ 0.01), with a relatively large correlation coefficient (r ≥ 0.1) for four factors, namely, the elevation, age of the tree, canopy density, and the number of trees per hectare, which were included in further modeling and prediction by MLR and GRNN.  Note: * and ** indicate the two-tailed test was significant at p ≤ 0.05 and p ≤ 0.01, respectively.
From the remote sensing images and based on previous research, 11 influencing factors were chosen, including bands 2-7, NDVI, RVI, DVI, EVI, and RI, as the primitive factors.
After the correlation analysis, the results (Table 6) showed that except for RI, the remaining 10 factors reached a very significant level (p ≤ 0.01), and the correlation coefficient was relatively large (r ≥ 0.1) in the seven factors of bands 2-7 and DVI, which were included in further modeling and prediction by MLR and GRNN.

Multiple Linear Regression
The testing data of 10,763 records were used for the prediction of DBH using Equations (3) and (4). Compared to the results produced by the traditional empirical models (Table 4, Figure 3), better outcomes were obtained. The results are displayed in Table 7 and Figure 4.

Multiple Linear Regression
The testing data of 10,763 records were used for the prediction of DBH using Equations (3) and (4). Compared to the results produced by the traditional empirical models (Table 4, Figure 3), better outcomes were obtained. The results are displayed in Table 7 and

Modeling and Prediction Based on GRNN
As with the MLR method, the remote sensing parameters were first left aside, and the four factors of the elevation, age of tree, canopy density, and the number of trees per hectare were used in the GRNN. Subsequently, a total of 11 factors-the seven remote sensing image factors, including bands 2-7 and DVI together with elevation, age of tree, canopy density, and number of trees per hectare-were used as the input for the GRNN model to estimate DBH.
The experiments showed that the optimal spread of the smoothing factor is 0.1 ( Figure 5 and Table 5).

Modeling and Prediction Based on GRNN
As with the MLR method, the remote sensing parameters were first left aside, and the four factors of the elevation, age of tree, canopy density, and the number of trees per hectare were used in the GRNN. Subsequently, a total of 11 factors-the seven remote sensing image factors, including bands 2-7 and DVI together with elevation, age of tree, canopy density, and number of trees per hectare-were used as the input for the GRNN model to estimate DBH.
The experiments showed that the optimal spread of the smoothing factor is 0.1 ( Figure 5 and Table 5).

Discussion
Based on the research data that were divided into 25,000 training samples and 10,763 testing samples (Table 1), the 13 empirical age-DBH models (Table 2), in which the only independent variable was the tree age, were used to fit the relationship between the tree age and DBH of the forest Figure 5. Scatter plot of the prediction of DBH for forest subcompartments by GRNN (a) with four independent variables of the elevation, age of tree, canopy density, and number of trees per hectare; and (b) with 11 independent variables including the elevation, age of tree, canopy density, number of trees per hectare, and seven remote sensing image factors including bands 2-7 and DVI.

Discussion
Based on the research data that were divided into 25,000 training samples and 10,763 testing samples (Table 1), the 13 empirical age-DBH models (Table 2), in which the only independent variable was the tree age, were used to fit the relationship between the tree age and DBH of the forest subcompartments (Table 3) and to predict the DBH ( Table 4). The modeling results of the 13 empirical models had R 2 values between 0.438 and 0.592. The MAPEs ranged from 15.63 to 19.3%, which means that the accuracies of the predictions were above 80%. The amplitude of R 2 prediction results obviously increased with a range of 0.2492-0.5092 compared to the modeling results. However, the Gompertz model provided the strongest generalizability with the best performance metrics of 1.7622 cm in MAE, 15.63% in MAPE, and 0.5092 in R 2 , which were even better than the modeling result. These findings are similar to those reported by Xu et al. and Lin et al. [57,58]. With the various size classes of DBH, Xu et al., taking 90 felled maple trees (Liquidambar formosana Hance) as the study object, and Lin et al., taking 90 felled Schima superba trees (Schima superba Gardn. et Champ.) as the study object, found that Gompertz was the optimal model, whether for DBH, tree height, or volume [57,58].
After the correlation analysis, the four factors from the inventory data for forest management planning and design and DEM were retained for further modeling and prediction by MLR and GRNN, namely, the elevation, age of the tree, canopy density, and the number of trees per hectare, as these satisfied p ≤ 0.01 and r ≥ 0.1. The order of absolute values of the correlation coefficients from high to low was the age of the tree, canopy density, the number of trees per hectare, and elevation. The R for the relationship between the number of trees per hectare and the DBH was negative, which indicated that competition among the trees affected the growth in DBH. This finding is similar to those reported by Luo et al. [53], with a negative correlation between the number of trees per hectare and DBH, and similar to Escalante et al. [54], in which forest stand simulations under fixed conditions indicated that the probabilities of diameter growth increase with site productivity increase, and decrease with an increase in the stand density index.
From the remote sensing images, the seven factors of bands 2-7 and DVI were also retained. The order of the absolute values of the correlation coefficients from high to low was B3, B2, B5, DVI, B6, B4, and B7.
The prediction results obtained by MLR and GRNN (Table 7) when using more environmental factors, such as topography and soil, were much better than those produced by the 13 empirical age-DBH models. This indicates that the environmental factors are beneficial for the estimation of DBH.
By adding the remote sensing image factors, the experimental results from MLR were further improved from 1.6006 cm in MAE, 14.96% in MAPE, 2.0707 cm in RMSE, 0.6193 in R 2 , and 0.0854 in TIC to 1.5300 cm in MAE, 14.38% in MAPE, 1.9976 cm in RMSE, 0.6459 in R 2 , and 0.0823 in TIC. Correspondingly, the experimental results from GRNN were improved from 1.4926 cm in MAE, 13.99% in MAPE, 1.9980 cm in RMSE, 0.6484 in R 2 , and 0.0826 in TIC to 1.4688 cm in MAE, 13.78% in MAPE, 1.9655 cm in RMSE, 0.6574 in R 2 , and 0.0810 in TIC. This indicates that the remote sensing image factors are beneficial to the estimation of DBH.
Regardless of whether remote sensing image factors were included, the experimental results from GRNN-with lower MAE, MAPE, RMSE, and higher R 2 -were better than those of MLR. In the forestry application of artificial neural networks (ANNs), He et al. used the ANN method to simulate and predict the DBH for different age groups and showed that the mean predictive accuracy of DBH for bamboo stands was satisfactory [16]. Vieira et al. estimated the growth of DBH and height of eucalyptus trees in 398 plots using artificial intelligence techniques (AIT), including ANN and an adaptive neuro-fuzzy inference system, and obtained good accuracy [59]. Although different from previous studies, this study still achieved an acceptable accuracy to predict DBH at the subcompartment's level with GRNN based on more diversified tree species and a wider range of samples totaling 35,763. Thus, the prediction results are more convincing.

Conclusions
The prediction results by MLR and GRNN when using more environmental factors, such as topography and soil, were much better than those produced using traditional empirical age-DBH models. The use of comprehensive factors, including ecological and biological factors, was beneficial for improving the prediction accuracy of both the MLR and GRNN models. When remote sensing factors were included, the performance metrics were further improved in the MLR and GRNN models, with MAE decreasing by 0.0238 and 0.0706 cm, MAPE decreasing by 0.21% and 0.58%, RMSE decreasing by 0.0325 and 0.0731 cm, R 2 increasing by 0.09 and 0.0266, TIC decreasing by 0.0031 and 0.0016, respectively. Regardless of whether remote sensing image factors were included, the experimental results produced using GRNN were better than those using MLR. For modeling and prediction based on more complex tree species and a wider range of samples, the GRNN model is more desirable due to its stronger generalizability.