^{1}

^{1}

^{1}

^{*}

^{1}

^{2}

In natural forest ecosystems, there is often abundant down dead wood (DDW) due to wind disasters, which greatly changes the size and structure of forests. Accurately determining the DDW volume (DDWV) is crucial for sustaining forest management, predicting the dynamic changes in forest resources and assessing the risks of natural disasters or disturbances. However, existing models cannot accurately express the significant spatial nonstationarity or complexity in their spatial relationships. To this end, we established a geographically weighted deep neural network (GWDNN) model that constructs a spatially weighted neural network (SWNN) through geographic location data and builds a neural network through stand factors and remote sensing factors to improve the interpretability of the spatial model of DDWV. To verify the effectiveness of this method, using 2019 data from Liangshui National Nature Reserve, we compared model fit, predictive ability and residual spatial autocorrelation among the GWDNN model and four other spatial models: an ordinary least squares (OLS) model, a linear mixed model (LMM), a geographically weighted regression (GWR) model and a deep neural network (DNN) model. The experimental results show that the GWDNN model is far superior to the other four models according to various indicators; the coefficient of determination R^{2}, root mean square error (RMSE), mean absolute error (MAE), Moran’s I and Z-statistic values of the GWDNN model were 0.95, 1.05, 0.77, -0.01 and -0.06, respectively. In addition, compared with the other models, the GWDNN model can more accurately depict local spatial variations and details of the DDWV in Liangshui National Nature Reserve.

Coarse woody debris (CWD) is an important material in forest ecosystems (

At present, commonly used spatial effect models can be divided into two categories: parametric models and nonparametric models (

With the development of computers, nonparametric models, including machine learning methods such as k-nearest neighbor (KNN), artificial neural network (ANN -

By realizing accurate modeling of the DDWV with stand factors and remote sensing factors, this study demonstrated the potential of GWDNNs to describe spatial relationships. This method establishes a two-layer DNN structure, in which one layer of the neural network is used to establish the spatially weighted matrix, and the other layer is used to estimate a parameter for each variable and realizes precise deconstruction and efficient calculation of spatial nonstationarity. To investigate whether the GWDNN model can improve the prediction of DDWV, the Liangshui National Nature Reserve in the city of Yichun (China) was considered as a study area. The coefficient of determination R^{2}, root mean square error (RMSE), mean absolute error (MAE), Moran’s I and Z-statistic values of five models (OLS model, LMM, GWR, DNN and GWDNN models) were compared to assess model fit, prediction accuracy and residual spatial autocorrelation. The advantages of the GWDNN model in predicting spatial relationships were verified.

The study area of Liangshui National Nature Reserve is located in the city of Yichun, Heilongjiang Province, on the east slope of the Dalidailing branch south of the southern Xiaoxing’an Mountains and is surrounded by 6 forest farms of the Dailing Forestry Experimental Bureau. The geographical location is 47° 07′ 19″-47° 14′ 40″ N, 128° 48′ 08″-128° 55′ 45″ E (

The data were collected during an investigation for the forest resource planning and design of the Liangshui reserve in 2019, which is called the second-class survey. For the purpose of forest resource management, 31 compartments and 464 subcompartments were defined. In addition, systematic sampling was conducted every 1 km. There were 130 fixed sample plots. Each plot was round-shaped with an area of 0.06 ha. The 130 sample plots comprised 65 coniferous and broad-leaved mixed forest plots, 31 broad-leaved mixed forest plots, and 34 coniferous relatively pure forests. The main variables recorded in 130 fixed plots were geographic location (latitude and longitude), elevation (from digital elevation model - DEM, m), slope (SLOPE, deg), aspect, average age of the dominant tree species (AGE, years), canopy, average diameter at breast height of the dominant tree species (DBH, cm), total standing forest stock (TSFS, m^{3}) and DDWV (m^{3}). For the determination of DDWV, we recorded tree species, quantity and DBH of the DDW in each sample plot and calculated DDWV through the unitary volume table for Liangshui National Nature Reserve.

The experimental data included the fixed sample plot data and visible light (R, G and B bands) remote sensing images, which were derived from the Chinese Academy of Forestry (CAF)’s LiCHy Hyperspectra (AISA Eagle II) airborne observation system (^{-1}. The image rate can be up to 160 Hz, the focal length is 18.1 mm, and the field of view (FOV) is 37.7″, with a spectral resolution of 9.6 nm and a spatial resolution of 0.52 m. In the imaging process of remote sensing images, due to changes in flight attitude, height and speed, the image pixels are distorted, offset, compressed or stretched relative to the actual position of the ground target. Therefore, we used the ENVI RPC model to perform geometric corrections for these deformations. Furthermore, for atmospheric correction, the ENVI FLAASH module was used to eliminate the absorption and scattering of electromagnetic waves from sunlight and ground objects by the atmosphere (

where

The stepwise regression method was used to screen the seven stand factors and seven remote sensing factors. Then, to examine the models for multicollinearity, variance inflation factors (VIFs) were calculated (

The traditional linear regression model is based on the OLS method to establish a linear relationship, as shown in

where

The LMM can be used to incorporate fixed effects and a single compartment random effect, with 31 levels (

where

The GWR model is an extension of the traditional linear regression model and can establish local parameter estimates (

where _{i}, _{i}) denotes the spatial coordinates of location _{0}(_{i, }_{i}) is the intercept, and β_{k}(_{i, }_{i}) is the coefficient of

where _{i, }_{i}) is an

while for the bisquare function (

where _{ij} is the Euclidean distance. The model parameter estimation and prediction accuracies largely depend on the bandwidth choice. We used the corrected Akaike information criterion (AICc) to select the GWR model with the optimal bandwidth (

With the development of computer technology, DNNs have a wide range of applications in various fields and can reliably determine the complex relationship between features (

where _{l} is an _{l}^{T} is the weight matrix, _{l} is the offset parameter vector, _{l} is the output mapping of the layer, and

In the DNN model, 4 layers of the neural network are designed to fit the complex relationships between the independent and dependent variables. Among them, the five selected variables conform to the input layer, the 2 layers in the middle are the hidden layers, and the prediction results are the output layer. The hidden layer uses a fully connected layer. The activation function is the rectified linear unit (ReLU) function, which is computationally efficient, has good sparsity and does not require unsupervised pretraining. The parameter initialization strategy enables the neural network model to learn useful information during the training process. However, the neural network model is prone to fall into local optima and is easily overfit because the initial model is complex and the amount of noise in the training set is too large. It is necessary to use the dropout algorithm to discard neurons in the hidden layer with a certain probability in the iterative training process, but their weights will be retained. The network participates in the next training iteration, and the weights are updated, which is equivalent to using multiple models to train the same dataset. Finally, the weights are assigned to the neurons to optimize the network structure.

DNNs rarely consider the spatial weight problem independently and rarely consider or directly add location coordinates into the network as variables (

where _{i}, _{i}) is a 6×6 diagonal matrix representing the geographical weights for observation

In the GWR model, spatial weighting is achieved with a function such as the Gaussian function or double square function, but its structure is simple, and it is difficult to accurately estimate nonstationarity due to the difficulty in choosing the correct kernel function for complex relationships. To address this problem, using the superior fitting ability of neural networks,

where [_{i1}_{i2}_{ij}] are the distances from point

The entire training processes of the DNN and GWDNN are shown in

The OLS model and LMM were implemented with the “nlme” library in R, and the GWR model was implemented using the “mgwr” library in Python. The DNN and GWDNN models were implemented using the “TensorFlow” and “Keras” libraries in Python.

Based on ground-truthed DDWV samples, a total of 104 samples were randomly selected from among the 130 samples from Liangshui National Nature Reserve for model training, and 26 samples were used for model verification. The R^{2}, RMSE and MAE values were used to assess model fit and prediction ability.

Spatial effects include spatial autocorrelation and nonstationarity. Ignoring the spatial effects of a model will misleadingly reduce the significance of a test and the model’s predictive ability (

where _{i} and _{j} are the values of the sample locations _{ij}(

The global Moran’s

This study used Z-statistics (

where (

The dependent variable and the five standardized independent variables SLOPE, AGE, TSFS, VDVI and VEG were used to establish the five models. The regression coefficients are significant at the

The setting of the parameters in the DNN and GWDNN models is particularly important in the training process. The DNN is fully connected with a total of 4 layers. The network layer includes 1 input layer, 2 hidden layers and 1 output layer. The numbers of neurons are set to

With reference to the ground-truthed DDWV data, the results of the five different models were compared. According to the general kriging interpolation in ArcGIS® v. 10.4 (ESRI, West Redlands, CA, USA) geostatistical analysis, the spatial interpolation process is similar to the weighted moving average (

The overall trend in the measured DDWV and the spatial distributions from the five models are consistent. The DDWV values ranged from 2.00 to 7.00 m^{3}. To facilitate comparison, the values were divided into seven levels by equal intervals of 1 m^{3} (^{3}, while low DDWV regions are located near the river and road in the middle and western parts of Liangshui National Nature Reserve, with values ranging from 0.00 to 2.00 m^{3}. The DNN model predicted poor levels in this region. The map resulting from the GWDNN model exhibited more explicit spatial variations than the other maps.

On the whole, the GWDNN model results are closer to the measured data and can describe the spatial distribution of DDWV in more detail than the other model results. Such maps can guide resource allocation for forest management. The forest DDWV mapping results of these models were insufficient for evaluation in this study, as we were limited by sample size. Future verification work should be conducted; conventionally, such verification work is performed by using independent sample sets or acknowledged high-accuracy results such as airborne data, especially unmanned aerial vehicle LiDAR data (

Currently, in the already existing models for mortality and CWD distribution, the predictors are usually stand factors, which have been confirmed to have great potential in previous studies (

The model coefficients of traditional parametric models include a substantial amount of information (

The fitting accuracy and prediction accuracy were compared among the five models based on the R^{2}, RMSE, MAE, Moran’s

The potential of the different models was demonstrated during the modeling process: establishment of the OLS model is simple and fast, and the correlations between different variables can be found intuitively. The GWR model considers the model results under different bandwidths. The advantage of the GWR model over the LMM is the possibility of determining the spatial location of every parameter without additional measurements (

DDWV was predicted by using parametric and nonparametric models of stand and remote sensing factors. The nonparametric models were found to be far superior to the parametric models in terms of fitting and verification accuracy, and the Moran

In summary, the results of this study suggest that the feasibility of the GWDNN model as a spatial model should be further investigated and promoted in the future. DNNs are popular algorithms that can mine data better than other networks. However, due to the lack of a clear connection between network parameters and approximate mathematical functions, DNNs are often called “black boxes”.

Y.S. and W.J. designed the research; Y.S. Y.C. and K.X. performed the experiments; Y.S. and Z.A. conducted the analysis and wrote the paper. All authors reviewed the manuscript.

This research was funded by the National Natural Science Foundation of China, grant no. 31870622, and the Special Fund Project for Basic Research in Central Universities, grant no. 2572019CP08. The authors are grateful to the First Investigation and Planning Institute of Forestry and Grassland in Heilongjiang Province and Liangshui National Nature Reserve for providing the second-class survey data in 2019.

The study area in Yichun (northeastern China).

The proposed GWDNN model for point

Flow chart of the DNN and GWDNN training process.

The spatial correlation between residuals for the five models.

Spatial distribution of (a) the ground-truth DDWV (down dead wood volume, m^{3}), and (b-f) the results of the five predicting models.

Descriptive statistics of the dependent variable, stand variables and remote sensing variables.

Variables | Name | Mean | SD | Min | Q1 | Q3 | Max |
---|---|---|---|---|---|---|---|

Dependentvariable | Down dead wood volume (DDWV, m3) | 3.44 | 4.72 | 0.01 | 0.73 | 3.9 | 25.63 |

Standvariables | Elevation (DEM, m) | 401.06 | 76.22 | 270 | 344 | 430 | 638 |

Slope (SLOPE, deg) | 8.58 | 5.21 | 1 | 6 | 10 | 36 | |

Aspect | 5.45 | 2.31 | 1 | 4 | 7 | 9 | |

Average age of dominant tree species (AGE, years) | 98.08 | 47.57 | 30 | 55 | 150 | 170 | |

Canopy | 0.53 | 0.14 | 0.3 | 0.4 | 0.6 | 0.9 | |

Average diameter at breast height od dominant tree species (DBH, cm) | 40.36 | 18.36 | 11.6 | 26 | 54 | 91 | |

Total standing forest stock (TSFS, m³) | 14.94 | 8.17 | 0.89 | 9.4 | 18.45 | 45.97 | |

Remote sensing variables | Visible-band difference vegetation index (VDVI) | 0.33 | 0.03 | 0.27 | 0.31 | 0.36 | 0.41 |

Excess green index (EXG) | 0.5 | 0.05 | 0.39 | 0.47 | 0.54 | 0.64 | |

Red:green ratio index (RGRI) | 0.58 | 0.05 | 0.44 | 0.54 | 0.61 | 0.7 | |

Excess green minus excess red index (EXGR) | 0.6 | 0.09 | 0.41 | 0.55 | 0.67 | 0.86 | |

Vegetation index (VEG) | 1.95 | 0.01 | 1.67 | 1.88 | 2.02 | 2.29 | |

Color index of vegetation (CIVE) | 18.39 | 0.02 | 18.33 | 18.38 | 18.4 | 18.44 | |

Normalized green-red difference index (NGRDI) | 0.27 | 0.04 | 0.18 | 0.25 | 0.3 | 0.39 |

OLS model, LMM, GWR model and GWDNN model results.

Model | Estimate | Intercept | SLOPE | AGE | TSFS | VDVI | VEG |
---|---|---|---|---|---|---|---|

OLS | Coefficient | 3.58 | 1.39 | 2.65 | -2.04 | 1.34 | -1.79 |

SE | 0.31 | 0.33 | 0.42 | 0.41 | 0.34 | 0.37 | |

P-value | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | <0.001 | |

VIF | - | 1.19 | 1.81 | 1.67 | 1.23 | 1.33 | |

LMM | Coefficient | 3.54 | 1.29 | 2.55 | -2.08 | 1.36 | -1.68 |

SE | 0.37 | 0.39 | 0.44 | 0.41 | 0.36 | 0.39 | |

P-value | <0.001 | 0.002 | <0.001 | <0.001 | <0.001 | <0.001 | |

GWR | Mean_Coefficient | 3.28 | 1.16 | 2.46 | -2 | 1.29 | -1.58 |

Min_Coefficient | 2.6 | 0.14 | 1.95 | -3.66 | 0.36 | -2.37 | |

Median_Coefficient | 3.15 | 1.32 | 2.24 | -1.82 | 1.13 | -1.57 | |

Max_Coefficient | 4.43 | 1.81 | 3.74 | -1.19 | 2.43 | -0.78 | |

GWDNN | Mean_Coefficient | 2.32 | 0.05 | 1.67 | -1.08 | 0.27 | -1.05 |

Min_Coefficient | -2.2 | -0.14 | -0.25 | -9.26 | -0.16 | -2.85 | |

Median_Coefficient | 2.18 | 0.04 | 1 | -0.42 | 0.32 | -0.9 | |

Max_Coefficient | 6.98 | 0.35 | 10.98 | 10.36 | 0.94 | -0.29 |

DNN model and GWDNN model hyperparameter settings.

Hyperparameter | DNNmodel | GWDNNmodel |
---|---|---|

Input | n×5 | n×130 (A),n×6 (BC) |

Hidden1 | n×32 | n×128, n×64 |

Hidden2 | n×32 | n×128, n×64 |

Hidden3 | - | n×128 |

Output | n×1 | n×1 |

LR | 0.001 | 0.001 |

Epochmaximum | 2000 | 2000 |

Dropout | 0.5 | 0.4 |

Batchsize | 10 | 5 |

Epoch stop | 1175 | 792 |

Simple cross-validation model evaluation of the training set and validation set.

Models | Training set | Validation set | ||||
---|---|---|---|---|---|---|

R^{2} |
RMSE | MAE | R^{2} |
RMSE | MAE | |

OLS | 0.57 | 3.08 | 2.26 | 0.43 | 3.60 | 2.57 |

LMM | 0.71 | 2.52 | 1.81 | 0.46 | 3.51 | 2.69 |

GWR | 0.69 | 2.59 | 1.91 | 0.52 | 3.31 | 2.34 |

DNN | 0.84 | 1.86 | 1.33 | 0.83 | 1.96 | 1.40 |

GWDNN | 0.95 | 1.05 | 0.77 | 0.85 | 1.82 | 1.39 |

Global Moran’s

Model | Moran’s I | Z-statistic |
---|---|---|

OLS | 0.24 | 2.29 |

LMM | 0.06 | 0.8 |

GWR | 0.05 | 0.64 |

DNN | 0.17 | 1.95 |

GWDNN | -0.01 | -0.06 |