Local neural-network-weighted models for occurrence and number of down wood in natural forest ecosystem

The natural forest ecosystem has been affected by wind storms for years, which have caused several down wood (DW) and dramatically modified the fabric and size. Therefore, it is very important to explain the forest system by quantifying the spatial relationship between DW and environmental parameters. However, the spatial non-stationary characteristics caused by the terrain and stand environmental changes with distinct gradients may lead to an incomplete description of DW, the local neural-network-weighted models of geographically neural-network-weighted (GNNWR) models are introduced here. To verify the validity of models, our DW and environmental factors were applied to investigate of occurrence of DW and number of DW to establish the generalized linear (logistic and Poisson) models, geographically weighted regression (GWLR and GWPR) models and GNNWR (GNNWLR and GNNWPR) models. The results show that the GNNWR models show great advantages in the model-fitting performance, prediction performance, and the spatial Moran’s I of model residuals. In addition, GNNWR models can combine the geographic information system technology for accurately expressing the spatial distribution of DW relevant information to provide the key technology that can be used as the basis for human decision-making and management planning.


Materials
Study area. Liangshui National Nature Reserve is affiliated with the Northeast Forestry University, and the study area is located in Yichun City of northeast China (Fig. 1). It belongs to the southeast section of the Xiaoxing'an mountains range, the eastern slope of the Daridailing branch, with a gross area is 12,133 ha, the core area of 6394 ha is the study area. There are many primitive Pinus koraiensis preserved in China and secondary birch and broad-leaved forests covering different succession stages of mixed broad-leaved primitive Pinus koraiensis forests. Other major tree species includes Spruce, Fir, Juglans mandshurica, Betula platyphylla, etc.
The data comes from the 2019-2020 Forest Resources Planning and Design Survey. There are 32 compartments and 464 sub-compartments; 31 compartments and 443 sub-compartments belong to the arbor forest land. We have also carried out down deadwood volume (including down wood, standing deadwood, and others) surveys on the plots 36 . In order to further manage the area, according to the actual situation, we investigated the occurrence and number of valuable down wood through corner gauge points to sub-compartments in this research. Meanwhile, the tree species, number of living trees (NLT), forest stand mean height (H), and diameter at breast height (DBH) in sub-compartments were summarized by using an angle gauge. Moreover, terrain factors, such as DEM (m), slope (°), aspect and stand factors, mean age of living trees (years), canopy, vegetation coverage (%), etc., are also recorded here. www.nature.com/scientificreports/ Variable selection. This study aims to utilize terrain and stand variables to introduce independent and dependent variable models. The independent variables are ODW (binary) and the NDW (count). The final retention was determined by using stepwise regression (significance level α = 0.05) 37,38 and following five dependent variables: the number of living trees (NLT), canopy, forest stand mean height (H), slope, and forest stand mean DBH (DBH). Table 1 describes statistical variables related to statistical techniques. In order to study the non-stationarity in each direction, the variation of each variable along the longitude and latitude is calculated in Fig. 2. Considering the NDW, the incidence of the DW under different longitudes, latitudes, and the affect DW variables under any latitudes and longitudes, the following results are observed: the probability of more than half, the NDW is observed commonly in 0-25 n/ha, with an increase in the longitude the phenomenon of the DW appears for more than 50 n/ha, and DW more frequently in middle latitudes in the region, but cannot clear the trend of change. Moreover, it can be seen from the figure that the variation trends of the longitude and latitude of the other five independent variables are also inconsistent. Further, the longitude of the slope gradually increases and along the latitude of slope gradually decreases and then increases, its value is the lowest in the mid-latitude region. In conclusion, the variation of each variable significantly changes along the longitude and latitude, which further indicates that the stand environment shows high spatial non-stationarity.

Methods
GL models. In many practical cases, the response variable is essentially binary 39,40 . That is, there are two possibilities in the result, namely, which are assigned a value of 0 (does not happen) or 1 (does happen). Here, we set ODW (y = 1) as P, and no ODW (y = 0) as 1−P. Logistic regression between the occurrence probability of the DW and their respective variables was established as shown in Eq. (1) 20,41 : where P represents the probability of the ODW and β 0 ∼ β k represents the regression coefficient of the model. k is the number of independent variables. The Poisson model is usually attributed to count data. Here, the Poisson model can be utilized to predict the NDW, in line with the Poisson variable Y and the Poisson probability density distribution function f (Y , µ) = e −µ µ Y /Y ! , where is the mathematical expectation and variance of the random variable Y , namely, E(Y ) = and V ar(Y ) = ; a monotonous average link function of response variables as a linear model is obtained by inducing some changes as shown in Eq. (2) 21 : It is assumed that the observed values are independent of each other, where µ represents the NDW and β 0 ∼ β k represents the regression coefficient of the model. k is the number of independent variables. All models estimate β 0 ∼ β k by the maximum likelihood method. GWR models. However, the above GL (logistic and Poisson) models are global in nature. The data collected in different geographical locations shows completely different results in the actual forestry survey due to the interference of different geographical environments and stand factors 42 . In the GL models analysis, it is often assumed that the estimated value of model coefficients are independent of the geographic location of the collected data, which leads to the estimated results tending to an average value. Thus, it can be inferred that all sample locations are based on unbiased estimates 43 . Therefore, the GL models show certain limitations in their application 35 . At present, GWR models are commonly used for improving the non-stationarity of space 44,45 .
Among them, we use the geographically weighted logistic regression (GWLR) model 46,47 to forecast the ODW. Eq. (3) is expressed as follows: The geographically weighted Poisson regression (GWPR) model 48 to forecast the NDW. Eq. (4) is expressed as follows: represents the coefficient of GWPR at the position i . All models estimate β 0 ∼ β k by the maximum likelihood method. GNNWR models. However, the kernel function of GWR models is relatively simple, and it is not easy to accurately model the complex stand geographical environment. For this reason, we propose a geographically neural-network-weighted regression (GNNWR) model, which is similar to the GWR models and uses the form of neural networks for defining the spatial non-stationary relationship 35 . Here, we integrate GWR into GL models, and the GNNWLR and GNNWPR models are shown in Eqs. (5) and (6): are the weight estimated at β 0 ∼ β k by using the corresponding logistic regression by the maximum likelihood method.
are the weight estimated at β 0 ∼ β k by using the corresponding Poisson regression by the maximum likelihood method.
The spatial weight w(u i , v i ) is expressed as shown in Eq. (7): www.nature.com/scientificreports/ In the GWR models, the spatial weights include the Gaussian function, double square function, etc., but their structure is relatively simple, so it is difficult to capture the complex relationship between the spatial distance and the non-stationary weight. Here, the calculation of the SWNN can be deemed to a complex nonlinear problem between the spatial distance and the weight 33,34 . The SWNN is utilized to develop the non-stationary weight matrix, and kernel weight is determined as expressed in Eq. (8): where d S i1 , d S i2 , · · · , d S in is the distance between point i and other samples. The neural network is designed using the Keras deep learning framework, also known as the deep neural network 49,50 . By setting the initial learning rate, the information between adjacent layers is considered fully connected (FC) and passed to the hidden layers [51][52][53] and the dropout algorithm is required to be used in the iterative process of the training models 54 . Here, we also set the batch size as the sample number of each training, defined as LeakyReLU 55 a nonlinear activation function for each network layer. The expressions for each layer are as expressed in Eq. (9).
where l is the number of all the layers in the training, x l are features of input layers, w t l is the weight matrix, b 1 is the offset parameter vector, y l are features of output layers, and σ is the activation function.
When the training process encounters the triggering condition of early stop or the time of epoch reaches the set maximum, the training is stopped. Once the training process is completed, the prescient capacity of the model is estimated by utilizing a validation set. In this paper, Huber 56 was used as the training loss function of the model, and the mean absolute error (MAE) of the validation set was used as the over-fitting evaluation index. By setting the maximum epoch, the trend of the MAE of the training set and the validation set was analyzed to find the optimal model parameters under the optimal number of iterations. To this end, its overall design framework of GNNWR (GNNWLR and GNNWPR) is shown in Fig. 3. Hyper-parameter settings of GNNWLR and GNNWPR models are shown in Table 2.
Model assessment. The model's performance is assessed by utilizing the coefficient of determination (R 2 ), which is utilized to assess the variability of the estimates, the standard deviation of the root-mean-square error (RMSE), and MAE estimates of the prediction error. The accuracy (acc) of 0 and 1 prediction under different  www.nature.com/scientificreports/ logistic models was also analyzed. Moreover, the revised Akaike information criterion (AICc) is utilized to select the best bandwidth of the GWR models. The residual of the model is defined as the difference between the observed values and the predicted values, and the spatial autocorrelation of each model residuals was analyzed by using correlation graphs of global Moran's I coefficients across a lag distance 36,57-59 .

Results
Model assessment. A total of 355 samples were randomly selected from 443 samples of the study area for model fitting and 88 samples for model validation. The standardized independent variables are NLT, canopy, H, slope, and DBH. The accuracy results are shown in Table 3. It can be seen that for both the training set and the validation set, the GNNWR (GNNWLR and GNNWPR) models are superior than GWR (GWLR and GWPR) models, and the GWR models are better than GL (Logistic and Poisson) models.
In order to verify that the spatial effects of GNNWR models include spatial non-stationarity, comparing spatial effect processes of different models can reduce the ability to misleading significance testing and prediction models. Moran's I and Z-value are calculated in Table 4. In order to better compare the spatial relationship of residuals of different models, spatial correlation graphs of residuals of different models are drawn with an interval of 300 m. Their average distance is about 300 m. Its neighbor pairs are 302, 1318, 2150, 2631, 3293 and 3751, respectively (Fig. 4).
From Table 4, it can be seen that GL and GWPR models (Z-value > Z α/2 = 1.96) indicate that the independent assumption of model residuals is contrary to these two types of models, and these results indicate a similar clustering pattern. In general, the ability to eliminate the autocorrelation of spatial residuals can be gleaned from the absolute value of Z-value, and it follows the order: GNNWR models > GWR models > GL models. It indicates that GNNWR models can effectively eliminate spatial non-stationarity. Moran's I of GNNWR models are also   www.nature.com/scientificreports/ relatively stable at any step size and close to 0, which indicates that GNNWR models demonstrate a good ability to maintain spatial stability.

Model parameter analysis.
According to the geographical location, the two types of models, GWR and GNNWR models, belong to local models and produce different geographical model coefficients. The descriptive statistics of the GWR and GNNWR model coefficients are summarized in Table 5. The distribution of GWR and GNNWR coefficients can be used to the non-stationary relationship between dependent and independent variables of stand environmental parameters 33,34 . The GWR and GNNWR models' coefficients show positive and negative fluctuations, indicating that the influence of the stand environment on DW may show opposite effects in different locations. Furthermore, the variation of the GNNWR coefficient is more dramatic than that of the response GWR coefficients, which may be the reason why the fitting and prediction performance of GNNWR has dropped significantly.

Visualized analysis of model parameters. The visual distribution analysis of five predictive variables
(NLT, canopy, H, slope, and DBH) is shown in Fig. 5a-e and according to the compartments. For further study the spatial correlation between ODW and NDW, the distribution diagrams of the model coefficients of the five predictor variables of the GNNWLR (Fig. 5f-j) and GNNWPR (Fig. 5k-o) models are presented. For the convenience of analysis, we divided the study area into nine orientations ( (Fig. 5p). In terms of the NLT (Fig. 5a), there are fewer NLT in the N and S areas and more in the M area. Canopy (> 0.7) (Fig. 5b) is mainly distributed in the S and N areas; H (Fig. 5c) is higher in the N and SE. There are many rivers and roads in M area of Liangshui Nature Reserve. The slope in this area is slow (< 11°) (Fig. 5d). The DBH (Fig. 5e) of big living trees is > 35 cm and mainly clustered in N, W, and E areas. Therefore, it can be seen that in the N and SE areas, there are trees with less NLT and big DBH. The M area has medium-sized trees with a gentle slope. In the S area, there are more small trees with more NLT. Their slope is steep in the N, NE, E, and SE areas. The coefficients vary significantly in space and show several directional patterns, and the two symbols and sizes are often heterogeneous, which also explains the significant spatial non-stationarity of the DW space.

Discussion
As the study area of Liangshui National Nature Reserve, there is a large DW affected by the wind storms, especially most of the trees are old primitive mixed broadleaf-conifer forest. The trees that have passed the mature stage become more and more vulnerable in the face of wind as the H taller, and the canopy grows larger, the physiological aging and the ability to resist diseases and insect pests decrease, and most of them end their life cycle in the state of breakage or uprooting 60,61 and such disaster changed its native forest stand structure 62 , which www.nature.com/scientificreports/ is associated with the hardwood species changing into the broad-leaved mixed forest and shows the positive pioneer trees, such as oak, birch, and aspen 63 . This indicates the formation of the secondary forest. Not only that, but it also leads to the loss of a lot of valuable timber in this area 64 . Therefore, it shows the economic value to the region (not rotten wood), and buck and forest investigate the environmental factors by using the data of the ODW (binary) 19,20 and the NDW (count) 21 to establish the GL, GWR, and GNNWR models, respectively. To better understand spatial non-stationarity caused by DW, the terrain with obvious gradient change and stand environment are used in the study area. www.nature.com/scientificreports/ Here, we compare the three methods and find that the non-stationary matrix constructed by the SWNN can well solve the complex nonlinear problem between the weight space distance and the weight space and thus can better predict the spatial change of the DW 33,34 . This method also has good results to solve complex non-stationary and nonlinear problems in coastal environments [33][34][35] . It is found that the GNNWR model shows absolute advantages in both the training set and the validation set 35 . At the same time, we utilized Moran's I and Z-value to check the non-stationarity of model space and observed that GL and GWR models would have independent assumptions and inefficient model coefficient estimates. In addition, according to the spatial residual correlation graph drawn at every 300 m, the Moran's I of GNNWR models at any step size is relatively stable and close to 0, which indicates that GNNWR models show a very good ability to maintain spatial stability 65 . At this moment, drop lines (logistic regression) and point plot (Poisson regression) with the ground survey data, respectively, and draw its 1:1 linear fitting diagram. The results of GNNWR models are closest to the ground-truth data (Fig. 6).
The distribution of ODW and NDW are easily affected by the climate, stand, terrain, and tree factors 61 . For the investigation, variables in the Liangshui Nature Reserve, four stand variables (NLT, canopy, H, and DBH), and one terrain variable (slope) 66 were selected by stepwise regression. NLT represents the density of the stand, and the canopy represents the size of the canopy structure. The DBH and H represent the tree size [67][68][69] ; the slope represents the geographical characteristics of stands 66 .
In GNNWR models, the symbols and sizes of the model coefficients are often heterogeneous; that is, different stand environments show different or even opposite effects on the distribution of the DW 33,34 . Here, we analyze the 5 variables of GNNWLR and GNNWPR models and the 9 orientations of model coefficients(ODW and NDW), and we can judge whether they are positively or negatively correlated with the corresponding dependent variables from the positive and negative mean values of variable model coefficients. However, we can see from the range of the box that these are not the only relations (Fig. 7), which will vary according to different regions, which also shows the existence of non-stationarity in the spatial model 35 . In the E and W orientations, the ODW is high and the NDW is large 64 . In the N, NE and NW orientations, the ODW is high, but NDW not large. In the region of N and NE, NLT is less, but the trees size and the tree canopy are large and the slope is steep. However, in the NW area, trees are of medium size and more NLT. However, SE, S and SW orientations are not easy to have DW, and the NDW is small. The tree size is small, the canopy is big, NLT is medium, the slope is gentle. In the M orientation, there are certain ODW and NDW 63 .
Finally, the predicted ODW and NDW on the visualization distribution maps of GNNWR closest to the ground-truth results, combined with statistical graphics and GIS mapping ability 70 , can provide inverted wooden key visual information, and forest managers can distinguish areas that are affected by disaster in order to help the administrator to provide timber and forest terrain variables and the relationship between the detailed information. Therefore, GNNWLR models are used to evaluate the ODW in a given area or small class, and GNNWPR models are used to predict the NDW. In order to prevent wind damage in the future, stand density, structure, and species composition can be changed in local areas 63 , which can assist in decision making and management planning for rational afforestation and management activities in the natural forest ecosystem.

Conclusion
Logistic regression (GL, GWLR, and GNNWLR) models and Poisson regression (GL, GWPR, and GNNWPR) models were used to model terrain variables and stand variables, respectively, for predicting the ODW and NDW. The analysis indicated that GNNWR models offer greater advantages than GWR and GL models in model-fitting