A Levenberg–Marquardt Backpropagation Neural Network for Predicting Forest Growing Stock Based on the Least-Squares Equation Fitting Parameters

Abstract

Traditional field surveys are expensive, time-consuming, laborious, and difficult to perform, especially in mountainous and dense forests, which imposes a burden on forest management personnel and researchers. This study focuses on predicting forest growing stock, one of the most significant parameters of a forest resource assessment. First, three schemes were designed—Scheme 1, based on the study samples with mixed tree species; Scheme 2, based on the study samples divided into dominant tree species groups; and Scheme 3, based on the study samples divided by dominant tree species groups—the evaluation factors are fitted by least-squares equations, and the non-significant fitted-factors are removed. Second, an overall evaluation indicator system with 17 factors was established. Third, remote sensing images of Landsat Thematic Mapper, digital elevation model, and the inventory for forest management planning and design were integrated in the same database. Lastly, a backpropagation neural network based on the Levenberg–Marquardt algorithm was used to predict the forest growing stock. The results showed that the group estimation precision exceeded 90%, which is the highest standard of total sampling precision of inventory for forest management planning and design in China. The prediction results for distinguishing dominant tree species were better than for mixed dominant tree species. The results also showed that the performance metrics for prediction could be improved by least-squares equation fitting and significance filtering of the evaluation factors.

1. Introduction

To promote the management and long-term development of forest resources, managers and researchers need to retain the latest information about forest resources and track the spatial changes of forest landscapes [1]. Forest inventory can represent the forest extent, determine species composition, and track changes from the past to the present. Traditional ground surveys effectively provide objective and reliable information for monitoring and managing forest resources [2,3]. However, traditional ground-based field measurements are expensive, time-consuming, labor-intensive, and hard to implement, particularly in mountainous and forest areas. This lay a burden on forest management personnel and researchers [4,5].

In China, two main types of surveys of large-area forest resources are available: national forest inventory (NFI), which is repeated every five years, and inventory for forest management planning and design, which is conducted at 10-year intervals. However, for both ecological monitoring and forest product utilization, the long survey cycles do not satisfy the social needs [6].

Since the 1970s, satellite images, such as Landsat Thematic Mapper (TM) and Satellite Pour l’Observation de la Terre (SPOT), have been routinely employed to estimate continuous forest parameters. By combining remote sensing image with the actual plots data, managers and researchers can estimate forest variables [7,8,9].

Scholars have been trying to find more efficient and credible prediction models [10,11,12] to partially replace the traditional and expensive measuring techniques.

Due to their adaptability and flexibility, artificial neural networks (ANNs) have become a substitutable and effective method to simulate nonlinear and complex ecosystems like forests. With robust data structures and highly interrelated relationships, ANN models have become popular. Being only slightly affected by data quality and bias, ANNs can also learn more complicated models and data tendencies [2,13]. ANNs therefore have advantages that can be used to solve difficult problems in forest resource management [13,14,15,16].

Traditional BP algorithms are sensitive to convergence into a local minimum, which may lead to instability of the training or sub-optimal results [17]. The Levenberg-Marquardt (LM) algorithm provides is less sensitive to local converges and therefore, it provides a better learning training approach for the back-propagation network. It chooses a balanced tradeoff between the training speed and the stability [18].

As one of the fundamental indicators of forest inventory [19], forest growing stock helps to determine the productive capacity of an area identified as forest available for wood production. The spatiotemporal dynamics of forest growing stock of various tree species play a key role in the balance between a continuous supply of timber and the sustainable development of the whole ecosystem [20]. By applying biomass expansion factors, the forest growing stock can be transformed to estimate the aboveground and belowground biomass [21]. Moreover, the increase or decrease in forest growing stock is essential to calculate the carbon storage of the forestry department [22]. During the estimation or prediction of forest growing stock, data from permanent plots (the basic units of NFI, set up by systematic sampling methods at the intersection of kilometer networks referring to the topographic maps with 1:50,000 map scale, usually with an area of 0.0667 ha) are frequently chosen to validate various models [9,23,24,25,26]. However, in China, the number of subcompartments (the basic unit of inventory for forest management planning and design, divided by the terrain boundaries including ridge line, valleys, roads, etc. or forest ownership boundaries, with the maximum area of 15 ha in South China and 25 ha in other parts of China) is much larger than the number of permanent plots. For example, in Zhejiang province in China (Figure 1), the number of permanent forest plots was 4252 in 2004 [6]. In contrast, the number of subcompartments was 39,377 in the city of Longquan (Figure 1) in Zhejiang province in 2007. To reduce the cost of investigating, estimation of forest growing stock based on subcompartments is more valuable than methods based on permanent forest plots.

This study focuses on predicting forest growing stock. First, three models were designed according to whether the study samples were based on mixed tree species or dominant tree species groups and whether the evaluation factors could be fitted by the least-squares equation and filtering by setting the significance p ≤ 0.05. Second, we established an overall evaluation factor system with 17 factors. Third, we integrated the remote sensing images of Landsat TM, digital elevation model (DEM), and the inventory for forest management planning and design into a database. Last, the Levenberg–Marquardt backpropagation (LM-BP) model was used to predict the forest growing stock.

2. Materials

2.1. Study Area

The study area, Longquan city, located in the southwest of Zhejiang province, has a total area of 3059 km², which extends from 118°42′ E to 119°25′ E, and 27°42′ N to 28°20′ N. Longquan includes 22 basic forestry management departments that contain 19 townships, two forest farms, and one scenic area (Figure 1). There are 265,667 ha of forest resources, with 14.56 million cubic meters of forest growing stock and 84.2% forest coverage [27].

2.2. Evaluation Factors

Forest growth is comprehensively influenced by multiple environments that primarily contain three categories: climate, topography, and soil.

Due to the costly data, not all of the environmental factors are suitable for estimating forest growth [28], which has resulted in many researchers trying to choose some factors for their experiments and accordingly make some harvest predictions [29,30,31,32].

To balance the costs and the richness of information, we chose comprehensive evaluation factors (

Table 1. Comprehensive evaluation factors for forest growth.

Number	Evaluation Factor	Number	Evaluation Factor
1	Tree Age	2	Slope
3	Canopy Density	4	Soil Depth
5	A-layer Depth of Soil	6	Aspect
7	Elevation	8	Curvature
9	Solar Radiation Index	10	Topographic Humidity Index
11	NDVI	12	Band 1
13	Band 2	14	Band 3
15	Band 4	16	Band 5
17	Band 7

Note: NDVI is normalized differential vegetation index; bands 1 through 5, and band 7, are the bands of the images from Landsat Thematic Mapper.

), which had been applied in a previous paper [3].

2.3. Research Data

The research data included the administrative map Landsat TM, DEM, and inventory data for forest management planning and design.

The administrative map of the study area was provided by the Forestry Bureau of Longquan.

The first version data of the global digital elevation model (GDEM) were generated from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (a Japanese sensor which is one remote sensory device on board the Terra satellite launched into Earth orbit by NASA in 1999 and has been collecting data since February 2000) in 2009, with 30-m resolution, IMG data type, Universal Transverse Mercator (UTM) projection, and World Geodetic System 1984 (WGS84). The data set was supplied by the International Scientific & Technical Data Mirror Site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn).

The inventory data for forest management planning and design, containing 39,377 forest subcompartments with 83,078 subplots, were provided by the Forestry Bureau of Longquan in 2007. In order to eliminate erroneous and incorrect data, nonforest, zero-volume, and zero-canopy density subplot samples were removed by preliminary processing. Therefore, the survey data of forest resources planning finally included 28,707 forest subcompartments that consisted of 38,898 subplots. In the data, four dominant tree species (Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.), Masson pine (Pinus massoniana Lamb.), Taiwan pine (Pinus taiwanensis Hayata), and hard broadleaves accounted for the vast majority of the dominant tree species (

Table 2. Groups divided by the dominant tree species and the size of subplots.

Group Number	Dominant Tree Species	Number of Subplot	Total Number of Subplot	Proportion (%)
Group 1	Chinese fir	20,296	38,898	52.18%
Group 2	Masson pine	5989		15.40%
Group 3	Taiwan pine	2418		6.22%
Group 4	Hard broadleaves	9582		24.63%
	The other dominant tree species	613		1.58%

). We separated the research data into four groups (Table 2): Group 1, Chinese fir; Group 2, Masson pine; Group 3, Taiwan pine; and Group 4, hard broadleaves. The remaining 613 subplots, which consisted of 26 dominant tree species that included soft broadleaves, fruit trees, Chinese cryptomeria, etc., accounted for a very small proportion (1.58%) of the total data. Thus, the 613 subplots were removed because the small size and variety of dominant tree species might have caused the experiment results to be unstable and unreliable.

3. Methods

3.1. Study Scheme Design

For modeling and prediction, the remaining subcompartments were divided into two parts in which there were 11 basic forestry management departments for modeling and another 11 for predicting. Specifically, the modeling area included Badu town, Baoxi township, Chatian town, Jianchi town, Lanju township, Pingnan town, Shangyang town, city forest farm, Tashi township, Xiaomei town, and Zhuyang township (Figure 1); the remaining 11 basic forestry management departments were used for predicting (Figure 1). We then designed three schemes: Scheme 1 involved modeling and predicting for all mixed dominant tree species; Scheme 2 involved modeling and predicting for the four groups of dominant tree species, i.e., Chinese fir, Masson pine, Taiwan pine, and hard broadleaves; and Scheme 3 involved modeling and predicting for the four groups of dominant tree species based on the evaluation factors fitted by the least-squares equation.

3.2. Data Integration

For modeling and prediction, we integrated all the required forest resource survey data into the same relational database. The soil depth, A-layer depth of soil, tree age, canopy density, and forest growing stock data were obtained from the inventory for forest management planning and design. The elevation, slope, aspect, surface curvature, solar radiation index, and topographic humidity index data were first obtained from DEM. In addition, the normalized difference vegetation index (NDVI) and the characteristic values of the bands (Bands 1 to 5 and 7) were initially obtained from Landsat TM remote sensing images. The average unit volume (m³/ha) of the forest resources was the only predicted parameter.

3.3. Evaluation Factors Fitting

In Scheme 3, we pre-fitted the evaluation factors and the average volume per unit by the least-squares equation. The three evaluation factors—soil depth, A-layer depth of soil, and canopy density—did not need to be fit in advance because their increase in values usually indicates the increase in forest growing stock.

The detailed operating procedure was as follows:

Step 1: The normalized values for each evaluation factor of each sample were calculated with the following formula:

$$y_{i,j}=\frac{x_{i,j}-\min \left(x_i\right)}{\max \left(x_i\right)-\min \left(x_i\right)}$$

(1)

where x_i,j is the original value of the i^th evaluation factor of the j^th sample, max(x_i) is the maximum of all values for the i^th evaluation factor, min(x_i) is the minimum of all values for the i^th evaluation factor, and y_i,j is the normalized value of x_i,j.

Step 2: The average volume per unit (z_i,j) was grouped by rounding the values of y_i,j × 100, except for the three evaluation factors of soil depth, A-layer depth of soil, and canopy density. Then, their least-squares fitting equations between z_i,j and rounded y_i,j × 100 were established (

Table 3. The least-squares fitting equations for the evaluation factors for Scheme 3.

Number	Dominant Tree Species	Evaluation Factor	Number of Group	Least-Squares Fitting Equation	Correlation Coefficient (R)	Significance (p)
1	Chinese fir	Tree Age	44	−0.0427890198 × x2 + 6.4914989075 × x + −10.8999069429	0.8754	0.0000
2		Slope	86	−0.0115959263 × x2 + 0.6038878174 × x + 82.1824866639	−0.5336	0.0000
3		Aspect	99	−0.0033003048 × x2 + 0.4839666164 × x + 72.0559152822	0.3246	0.0010
4		Elevation	85	0.0023254337 × x2 + −0.2572556056 × x + 89.0182097406	−0.0945	0.3896
5		Curvature	64	0.0065366953 × x2 + −0.5371817451 × x + 90.5614723394	−0.1003	0.4305
6		Solar Radiation Index	82	−0.0026200461 × x2 + 0.3467036840 × x + 73.3921706799	0.0918	0.4120
7		Topographic Humidity Index	91	0.0043653961 × x2 + −0.2512069504 × x + 88.2155051143	0.1632	0.1221
8		NDVI	98	0.0022565071 × x2 + −0.1541807892 × x + 86.1625587414	0.1509	0.1380
9		Band 1	35	−0.0118922876 × x2 + 1.5080606424 × x + 43.2161510234	0.4637	0.0050
10		Band 2	22	−0.0127550587 × x2 + 1.4717083307 × x + 48.6242224527	0.3389	0.1228
11		Band 3	32	−0.0090806959 × x2 + 0.7815481904 × x + 70.1196285883	0.0042	0.9818
12		Band 4	53	0.0020084727 × x2 + −0.0805199224 × x + 84.5622590844	0.2199	0.1136
13		Band 5	80	−0.0042838563 × x2 + 0.4542927464 × x + 74.8158027009	0.0884	0.4355
14		Band 7	41	0.0004892917 × x2 + −0.0031706472 × x + 84.1887249310	0.0940	0.5586
15	Masson pine	Tree Age	38	0.0331501479 × x2 + 2.8278931378 × x + 15.2328173465	0.9482	0.0000
16		Slope	78	−0.0148356773 × x2 + 0.9016394596 × x + 67.9666405506	−0.3331	0.0029
17		Aspect	94	−0.0032729759 × x2 + 0.2348500208 × x + 75.6303889833	−0.1668	0.1081
18		Elevation	70	−0.0206467958 × x2 + 1.1473753074 × x + 65.8295224974	−0.4560	0.0001
19		Curvature	47	−0.0065277438 × x2 + 0.3054885854 × x + 74.8842350529	−0.0234	0.8758
20		Solar Radiation Index	67	−0.0116390188 × x2 + 1.7344282575 × x + 16.7808352190	0.5074	0.0000
21		Topographic Humidity Index	79	0.0068459173 × x2 + −0.6174814023 × x + 82.0851684875	0.0086	0.9402
22		NDVI	90	0.0040823038 × x2 + −0.1607289597 × x + 72.4870953181	0.2413	0.0220
23		Band 1	30	0.0215449663 × x2 + −1.4478075547 × x + 83.4298332936	0.4954	0.0054
24		Band 2	19	−0.0107774671 × x2 + 1.1612851686 × x + 43.9676859378	0.1036	0.6731
25		Band 3	28	−0.0056810995 × x2 + 0.6748729435 × x + 57.9892594691	0.1594	0.4178
26		Band 4	49	0.0051865652 × x2 + −0.4677663146 × x + 83.9501363669	0.0103	0.9442
27		Band 5	72	0.0010562299 × x2 + −0.2748031246 × x + 86.6400382977	−0.1567	0.1888
28		Band 7	39	0.0121557341 × x2 + −0.9500200238 × x + 83.6278316688	0.1422	0.3880
29	Taiwan pine	Tree Age	46	0.0108691950 × x2 + 2.8041406197 × x + 0.7548379420	0.9153	0.0000
30		Slope	71	0.0007210346 × x2 + −0.1981166266 × x+80.5153065313	−0.0985	0.4140
31		Aspect	92	−0.0033269117 × x2 + 0.4194282364 × x + 62.0507244967	0.1330	0.2064
32		Elevation	67	−0.0170165295 × x2 + 2.0442044776 × x + 13.8462032002	0.1801	0.1447
33		Curvature	52	−0.0089471129 × x2 + 0.1602318306 × x + 74.2890687778	−0.2853	0.0404
34		Solar Radiation Index	77	0.0024392679 × x2 + −0.1965197674 × x + 72.5025712358	0.1011	0.3816
35		Topographic Humidity Index	33	−0.0162849784 × x2 + 1.1347361067 × x + 56.4756479999	0.0568	0.7537
36		NDVI	85	−0.0018733247 × x2 + 0.3472925064 × x + 60.4853464235	0.2063	0.0581
37		Band 1	15	0.0584268211 × x2 + −2.8925728373 × x + 102.4213645020	0.1372	0.6259
38		Band 2	11	0.0462181612 × x2 + −2.0249896794 × x + 87.1407802797	0.4339	0.1824
39		Band 3	17	0.0520556446 × x2 + −1.8273924977 × x + 82.6589015604	0.3448	0.1753
40		Band 4	45	−0.0026199728 × x2 + 0.2061151757 × x + 69.9711400898	−0.0155	0.9195
41		Band 5	63	0.0012177555 × x2 + −0.1554842686 × x + 74.2397400623	−0.0815	0.5255
42		Band 7	32	0.0142484577 × x2 + −0.6341009191 × x + 75.5528181515	0.3688	0.0378
43	Hard broadleaves	Tree Age	62	0.0163003073 × x2 + 1.1534510399 × x + 17.1575498742	0.9163	0.0000
44		Slope	85	0.0036011768 × x2 + −0.6254997771 × x + 74.8775043187	−0.3450	0.0012
45		Aspect	96	0.0039735998 × x2 + −0.3577728963 × x + 58.2501579702	0.0960	0.3520
46		Elevation	83	0.0180034120 × x2 + −1.3930520666 × x + 71.4995888064	0.1853	0.0935
47		Curvature	60	0.0286950221 × x2 + −1.5689644141 × x + 64.8735422419	0.1606	0.2203
48		Solar Radiation Index	87	0.0040370576 × x2 + −0.2975027144 × x + 55.5262852886	0.1829	0.0900
49		Topographic Humidity Index	62	0.0277373791 × x2 + −1.8961270968 × x + 72.8165360650	0.2434	0.0566
50		NDVI	93	−0.0027252454 × x2 + 0.2307572964 × x + 48.9924512551	−0.0468	0.6558
51		Band 1	29	0.0626190912 × x2 + −4.5932489190 × x + 123.0293545383	0.3439	0.0678
52		Band 2	18	0.0292226631 × x2 + −2.1099893402 × x + 81.3753529658	0.5453	0.0192
53		Band 3	28	0.0226705539 × x2 + −1.2930119073 × x + 66.0857518410	0.4885	0.0084
54		Band 4	52	0.0029985089 × x2 + −0.4153878341 × x + 62.4330118645	−0.2118	0.1317
55		Band 5	69	0.0133684136 × x2 + −0.9407788330 × x + 62.4965426573	0.1885	0.1210
56		Band 7	33	0.0297143535 × x2 + −1.7693736017 × x + 67.6223561210	0.3281	0.0623

Note: NDVI is normalized differential vegetation index; bands 1 through 5, and band 7, are the bands of the images from Landsat Thematic Mapper.

).

Step 3: In Table 3, R denotes the correlation coefficient between each evaluation and the average volume per unit, and p denotes the significance of the least square fitting regression. There were some evaluation factors with lower values of correlation coefficient (R) and higher values of significance (p) in Table 3. Usually, P has three values, 0.001, 0.01, and 0.05. By letting p equal to 0.05 in this paper, the evaluation factor would be retained if its p ≤ 0.05. Therefore, for the dominant tree species Chinese fir, the remained factors were tree ages, slope, aspect, and band 1; for Masson pine, the remained factors were tree ages, slope, elevation, solar radiation index, NDVI, and band 1; for Taiwan pine, the remained factors were tree ages, curvature, and band 7; and for hard broadleaves, the remained factors were tree ages, slope, band 2, and band 3. Including the 3 retained-factors mentioned above, soil depth, A-layer depth of soil, and canopy density, there were eight evaluation factors remained in modelling and prediction for Chinese fir, ten evaluation factors remained for Masson pine, seven evaluation factors remained for Taiwan pine, and eight evaluation factors remained for hard broadleaves. To all remained evaluation factors, the fitted values of x_i,j were recalculated by the fitting equations (Table 3) and were renormalized by Equation (1).

3.4. Improved BP Neural Network Model Based on LM Algorithm

Although there are many feasible network models, the backpropagation neural network (BPNN) has outstanding assortment ability and pattern identification in practical applications [33]. The BPNN model minimizes the mean square error between the predicted values and the expected values by modifying the connection weight of network [34]. Usually, a BPNN comprises three types of successive layers: input layer, hidden layer, and output layer. Its operating processes can be described in two steps. First, the input signal is forward-propagated from the input layer, via the hide layer, to the output layer. During this process, the weight value and offset value of the network are maintained constant, and the status of each layer of neurons will only affect the next layer of neurons. The difference between the predicted output and expected output of the network is defined as the error signal. Second, the error signal is backpropagated from the output layer to the input layer in a layer-by-layer manner. During this process, the weight value of network is regulated by the error feedback to make the predicted output of network more closer to the expected one [34].

The LM algorithm is essentially an evolution of the gradient descent algorithm and Newton algorithm, which can considerably speed up the convergence rate by reducing the iteration process and then produce more accurate data. Compared with the disadvantages of traditional BPNNs, such as slow convergence speed and local minimum problems, the convergence rate of the LM algorithm is the fastest of all traditional or improved networks. The improved BPNNs with the LM algorithm have been shown to achieve excellent results, especially when applied to evaluation and prediction [3,35,36,37,38,39].

Therefore, taking the comprehensive evaluation factors (Table 1) as input variables and the forest growing stock as the output variable, the LM-BP was selected as the modeling and prediction model in this paper. As mentioned above, the values of comprehensive evaluation factors were normalized to inputs by Equation (1), both in Scheme 1 and Scheme 2. In Scheme 3, the inputs were preprocessed by three steps. First, the values of input variables were pre-fitted with the evaluation factors and the average volume per unit from the least-squares method, except for the three evaluation factors of soil depth, A-layer depth of soil, and canopy density. Second, the evaluation factor would be retained if its significance p ≤ 0.05 when fitted by the least square method. Last, the inputs were normalized by Equation (1).

The LM-BP computer program was coded by MATLAB R2011b, and its main steps were as follows:

Step 1: The parameters were set as follows:

• Max_Epochs = 1000
• Input_Num = 17
• Output_Num = 1
• Hidden_Neuron_Num = 2 × Input_Neuron_Num + Output_Neuron_Num
• TransferFcn = {‘tansig’ ‘purelin’}
• TrainFcn = ‘trainlm’
• PerformFcn = ‘mse’

where Max_Epochs is the maximum of training epochs, Hidden_Neuron_Num is the number of nodes in the hidden layer, Input_Neuron_Num is the number of nodes in the input layer, Output_Neuron_Num is number of nodes in the output layer, tansig is the transfer function that transfers values to the hidden layer from the input layer, purelin is the transfer function that transfers values to the output layer from the hidden layer, trainlm is the training function corresponding to the LM algorithm, and mse (mean square error) is the performing function.

Step two: Net creation.

• Net = newff (I, O, Hidden_Neuron_Num)

where I is the input matrix and O is the output vector.

Step three: Net training.

• [Net TR] = train (Net, I, O)
• Step 4: Prediction.
• y = sim (Net, I_test)

where Net is the trained net, I_test is the input matrix of predicting samples, and y is the prediction result.

3.5. Model Performance Metrics

The following equations were used to calculate the accuracy of each model. By letting t_i represent the measured values, y_i represent the predicted values at sample i, t denote the observed mean, and y denote the predicted mean, the first statistical measurement, i.e., the group absolute percentage error (GAPE), was defined with Equation (2):

$$\mathrm{GAPE}=\left|\frac{{{\displaystyle \sum }}_{i=1}^n\left(t_i-y_i\right)}{{{\displaystyle \sum }}_{i=1}^n{t}_i}\right|\times 100\%$$

(2)

where n is the number of samples in the test dataset.

The second criterion is the mean absolute error (MAE) as defined in Equation (3):

$$\mathrm{MAE}=\frac{1}{n}\left|{\displaystyle \sum }_{i=1}^n\left(t_i-y_i\right)\right|$$

(3)

Then, the mean absolute percentage error (MAPE) is defined by Equation (4):

$$\mathrm{MAPE}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\frac{\left|t_i-y_i\right|}{\left|t_i\right|}\times 100\%$$

(4)

The root mean squared error (RMSE) is one of the most frequently used factors for assessing the capability of ANNs, which can be obtained with Equation (5):

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(t_i-y_i\right)}^2}$$

(5)

The factor of agreement (IA) is a nondimensional factor of calculation that can be obtained with Equation (6):

$$\mathrm{IA}=\frac{{{\displaystyle \sum }}_{i=1}^n{\left(t_i-y_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(\left|t_i-\overline{t}\right|+\left|y_i-\overline{y}\right|\right)}^2}$$

(6)

The IA distinguishes between the expected and measured values, and it is restricted from 0 to 1. High value of IA means good consistency between values.

The correlation coefficient R² is defined with Equation (7):

$$R^2=\frac{{\left[{{\displaystyle \sum }}_{i=1}^n\left(t_i-\overline{t}\right)\left(y_i-\overline{y}\right)\right]}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(t_i-\overline{t}\right)}^2\times {{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(7)

4. Results

4.1. Modeling

During the modeling period, the optimal model was chosen from 10 experiments performed for each scheme. The results are presented in

Table 4. Detailed performance metrics of the modeling results for each scheme.

Scheme	Dominant Tree Species	Numbers of Sample	Average Measured Value (m3)	Average Predicted Value (m3)	GAPE (%)	MAPE (%)	MAE (m3)	RMSE (m3)	IA	R2
Scheme 1	Mixed	16112	76.8375	75.2040	2.1251	41.0680	20.4495	28.119	0.8773	0.6388
Scheme 2	Chinese fir	9529	85.0350	83.4120	1.9097	29.4761	18.0450	24.6525	0.9117	0.7212
	Masson pine	1872	79.1520	76.2300	3.6916	30.1500	18.1140	23.9415	0.9080	0.7107
	Taiwan pine	1457	70.5510	66.8235	5.2822	25.0961	14.9265	20.4195	0.9276	0.7605
	Hard broadleaves	3014	54.5385	54.3390	0.3673	44.5836	14.9295	21.5490	0.9181	0.7372
Scheme 3	Chinese fir	9529	85.0354	86.0979	1.2495	31.3399	18.3058	25.0376	0.9087	0.7116
	Masson pine	1872	79.1524	79.9650	1.0267	31.4837	18.4711	24.4194	0.9041	0.6949
	Taiwan pine	1457	70.5505	70.4016	0.2110	27.2093	15.4221	21.4346	0.9153	0.7261
	Hard broadleaves	3014	54.5391	54.2795	0.4760	45.3362	15.6780	23.0508	0.9049	0.6991

Note: GAPE is the group absolute percentage error defined by Equation (2); MAPE is the mean absolute percentage error defined by Equation (4); MAE is the mean absolute error defined by Equation (3); RMSE is the root mean squared error defined by Equation (5); IA is the factor of agreement defined by Equation (6); R² is the correlation coefficient defined by Equation (7).

.

In the three modeling schemes, the GAPE values were much less than 10%, meaning the accuracies of the group estimation were above 90%. The performance metrics MAE, RMSE, IA, and R² in Scheme 1 were all significantly inferior to those in Schemes 2 and 3, and the MAPE in Scheme 1 was generally inferior to that in Schemes 2 and 3. Through a pooled analysis of Scheme 1 to Schemes 2 and 3, the R² were significantly improved from 0.6388 to 0.7427 and 0.7272, respectively, and MAPE decreased from 41.0680% to 32.0223% and 33.6355%, respectively. Schemes 2 and 3 were much better than Scheme 1. Except GAPE, other performance metrics of Scheme 3 were slightly worse than those of Scheme 2 in the modelling results (

Table 5. The total performance metrics for each scheme for modeling and predicting.

Modeling or Predicting	Scheme Name	Number of Sample	Average Measured Value (m3)	Average Predicted Value (m3)	GAPE (%)	MAPE (%)	MAE (m3)	RMSE (m3)	IA	R2
Modeling	Scheme 1	16,112	76.8375	75.2040	2.1251	41.0680	20.4495	28.1190	0.8773	0.6388
	Scheme 2	15,872	77.2215	75.5205	2.2011	32.0223	17.1750	23.6415	0.9202	0.7427
	Scheme 3	15,872	77.2215	77.8916	0.8687	33.6355	17.5616	24.2850	0.9157	0.7272
Predicting	Scheme 1	22,786	72.8010	69.3690	4.7149	50.2675	24.1500	33.6555	0.8296	0.5279
	Scheme 2	22,413	73.2135	71.2950	2.6217	41.7724	20.9190	28.9950	0.8883	0.6460
	Scheme 3	22,413	73.2135	73.3415	0.1739	37.5268	19.5685	27.4908	0.9036	0.6823

Note: GAPE is the group absolute percentage error defined by Equation (2); MAPE is the mean absolute percentage error defined by Equation (4); MAE is the mean absolute error defined by Equation (3); RMSE is the root mean squared error defined by Equation (5); IA is the factor of agreement defined by Equation (6); R² is the correlation coefficient defined by Equation (7).

).

4.2. Predicting

The prediction results are presented in

Table 6. Detailed performance metrics of the prediction results for each scheme.

Scheme	Dominant Tree Species	Number of Sample	Average Measured Value (m3)	Average Predicted Value (m3)	GAPE (%)	MAPE (%)	MAE (m3)	RMSE (m3)	IA	R2
Scheme 1	Mixed	22,786	72.8010	69.3690	4.7149	50.2675	24.1500	33.6555	0.8296	0.5279
Scheme 2	Chinese fir	10,767	83.0250	80.9520	2.4971	31.5641	20.8275	28.6515	0.8862	0.6455
	Masson pine	4117	87.3510	80.2590	8.1192	37.8890	24.6330	34.0260	0.8657	0.5979
	Taiwan pine	961	82.4640	69.6090	15.5888	35.9952	25.5375	36.1290	0.8105	0.4979
	Hard broadleaves	6568	46.9170	50.0910	6.7675	61.7867	18.0675	24.6255	0.8604	0.5647
Scheme 3	Chinese fir	10,767	83.0250	83.6362	0.7361	32.8997	20.3668	28.2005	0.8897	0.6544
	Masson pine	4117	87.3510	84.8161	2.9019	34.4703	22.5600	31.4144	0.8942	0.6523
	Taiwan pine	961	82.4640	76.0733	7.7503	35.4529	22.6440	31.8299	0.8332	0.5601
	Hard broadleaves	6568	46.9170	48.8730	4.1704	47.3315	15.9346	22.5483	0.8893	0.6375

Note: GAPE is the group absolute percentage error defined by Equation (2); MAPE is the mean absolute percentage error defined by Equation (4); MAE is the mean absolute error defined by Equation (3); RMSE is the root mean squared error defined by Equation (5); IA is the factor of agreement defined by Equation (6); R² is the correlation coefficient defined by Equation (7).

and Figure 2. Table 6 shows the sample size and performance metrics of each scheme. Figure 2 shows the scatter plots between predicted values and measured values of forest growing stock.

During prediction, except the scheme for Taiwan pine as the dominant tree, the GAPE values of the other schemes were less than 10%, which means that the accuracies of the group estimation were above 90%. In Scheme 2, the GAPE for predicting Taiwan pine was as much as 15.5888%, so its group estimation precision was 84.4112%. However, in Scheme 3, the GAPE for predicting Taiwan pine has been improved greatly that in Scheme 2.

The pooled analysis for Schemes 2 and 3, achieved by merging the prediction results of the four dominant tree species into a whole of prediction result, and the total performance metrics were recalculated and are shown in Table 5. From Scheme 1 to Schemes 2 and 3, the R² were improved from 0.5279 to 0.6460 and 0.6823, respectively, and MAPE decreased from 50.2675% to 41.7724% and 37.5268%, respectively. The other performance metrics, including GAPE, MAE, RMSE, and IA, all improved to some extent from Scheme 1 to Schemes 2 and 3.

The GAPEs in all schemes were less than 5% (Table 5), which indicated that the total group estimation precisions were much greater than the highest standard of the total sampling precision of forest resource inventory in China (GB/T 26424-2010).

5. Discussion

The results of modeling indicated that the performance metrics for distinguishing dominant tree species were much better than that for determining mixed dominant tree species. However, except GAPE, other performance metrics of Scheme 3 were slightly worse than those of Scheme 2 in the modelling results. This probably occurred because we used the values fitted between the evaluation factors and the average volume per unit using the least-squares equations instead of the original values of the evaluation factors before modeling, and remove the non-significant evaluation factors when their P-values were greater than 0.05, meaning the fitted values were closer to their averages than their original values; meanwhile, some detail in the information could have also been lost.

The prediction results indicated that the performance metrics for distinguishing dominant tree species were also better than for determining mixed dominant tree species, which is similar to most research results [40,41]. This result was apparently different from the published forest growing stock simulation [8], which used the k-nearest-neighbor estimation method to estimate the forest growing stock based on Landsat TM images and stand-level field inventory data. The results of that study showed that the volumetric estimation errors of different tree species were significantly higher than those of the overall estimations.

Regardless of the modeling results, all the prediction performance metrics of Scheme 3 were better than those of Scheme 2, implying that the prediction results that were based on the least-squares equation fitting and significance p filtering for evaluation factors were more stable and reliable.

In Schemes 1 and 2, the predicted values were systematically lower than the measured values from the modeling results and the prediction results. This may be due to the modeling areas being chosen manually rather than randomly. However, in Scheme 3, for both modeling and prediction, the average predicted values were very close to the average measured values. This may be explained by the fitted factor values obtained using the least-squares equations being closer to their averages than their original values and removing the non-significant evaluation factors to reduce the random noise, meaning the predicted results would be closer to the expected values.

Despite high MAPE values, the group estimation accuracy from our study may be good enough for some purposes. For example, for forest management planning and design at a regional level, the estimated results should be a partial substitution for field investigation. An acceptable level of estimation error for total forest volume is about 20% when the inventory area is larger than 30 ha [8,42]. The total performance metrics indicate that all GAPE values were lower than 5%, meaning the overall accuracies were over 95%, which exceed the highest standard of total sampling precision (90%) of inventory for forest management planning and design in China. Especially in Scheme 3, both the modeling results and the predicting results had GAPE values lower than 2%, meaning the overall accuracies were over 98%. Therefore, the results should be used in forest management planning and design in situations where there is a lack of available measured data.

The prediction results could likely be further improved with further research. There are, however, two limitations to our study: overfitting and limited network layers. In our study, we found that higher modeling accuracy did not necessarily mean higher prediction accuracy, which might have been caused by overfitting in the modeling procedure. Due to the full connection network structure, the number of the layers was very limited, therefore, the prediction accuracy was also limited. The next study plan may depend on a new model, such as the deep-learning model.

6. Conclusions

Three conclusions were drawn from the data analyses in this study. First, the prediction results of all schemes showed that their average predicted values were close to the average measured values and all GAPE values were lower than 5% (Table 5), indicating that group estimation precision of all predicting schemes were much greater than 90%, which is the highest standard of the total sampling precision for volume of forest resource inventory in China. Second, the prediction results for distinguishing dominant tree species were better than for mixed dominant tree species. Lastly, the results showed that the performance metrics for prediction could be improved using least-squares equation fitting for the evaluation factors and by using the significance P filtering the evaluation factors.

We chose one part of an area in the research region as the modeling area and the other as the prediction area instead of randomly choosing modeling samples [3,43]. This procedure is more convenient in practice because the field measurements in the inventory for forest management planning and design are recorded in some contiguous areas while others are predicted. This research also reflects the characteristics of sufficient samples and reduced data costs during the prediction in which the canopy density is the only factor that should be remeasured by observation or a line-intercept method in the plot. Therefore, the costs of data acquisition are much lower than those of traditional models that are used to survey the diameter or height of each tree [3,44,45]. In addition, the pattern is more suitable for the prediction of volume for forest under the state of natural growth or with minimal human interference. This model offers the possibility to monitor forest biomass within a short period based on a conversion of volume to biomass.