Generating a Baseline Map of Surface Fuel Loading Using Stratified Random Sampling Inventory Data through Cokriging and Multiple Linear Regression Methods

Abstract

Surface fuel loading is a key factor in controlling wildfires and planning sustainable forest management. Spatially explicit maps of surface fuel loading can highlight the risks of a forest fire. Geospatial information is critical in enabling careful use of deliberate fire setting and also helps to minimize the possibility of heat conduction over forest lands. In contrast to lidar sensing and/or optical sensing based methods, an approach of integrating in-situ fuel inventory data, geospatial interpolation techniques, and multiple linear regression methods provides an alternative approach to surface fuel load estimation and mapping over mountainous forests. Using a stratified random sampling based inventory and cokriging analysis, surface fuel loading data of 120 plots distributed over four kinds of fuel types were collected in order to develop a total surface fuel loading model (lntSFL-BioTopo model) and a fine surface fuel model (lnfSFL-BioTopo model) for generating tSFL and fSFL maps. Results showed that the combination of topographic parameters such as slope, aspect, and their cross products and the fuel types such as pine stand, non-pine conifer stand, broadleaf stand, and conifer–broadleaf mixed stand was able to appropriately describe the changes in surface fuel loads over a forest with diverse terrain morphology. Based on a cross-validation method, the estimation of tSFL and fSFL of the study site had an RMSE of 3.476 tons/ha and 3.384 tons/ha, respectively. In contrast to the average loading of all inventory plots, the estimation for tSFL and fSFL had a relative error of 38% (PRMSE). The reciprocal of estimation bias of both SFL-BioTopo models tended to be an exponential growth function of the amount of surface fuel load, indicating that the estimation accuracy of the proposed method is likely to be improved with further study. In the regression modeling, a natural logarithm transformation of the surface fuel loading prevented the outcome of negative estimates and thus improved the estimation. Based on the results, this paper defined a minimum sampling unit (MSU) as the area for collecting surface fuels for interpolation using a cokriging model. Allocating the MSUs at the boundary and center of a plot improved surface fuel load prediction and mapping.

1. Introduction

Wildfires are recognized as one of the major disturbances in terrestrial forest ecosystems. Fire can significantly change forest attributes and destroy the habitat of wildlife while at the same time, it can create another important habitat. Naturally occurring wildfires caused by lightning or extreme climate events (a long dry season or drought and high temperature) are generally inevitable. Fires are also frequently used to clear forest for large-scale plantation establishment [1] and small-scale agricultural uses [2,3,4]. Although rapid warming has recently resulted in more wildfires worldwide, human activities have been recognized as the major causes of wildfire [5]. Most human-caused wildfires can therefore be prevented by using fires responsibly and taking preventative measures [6].

During a long, dry heat-wave period, an original wildfire can occur naturally and may become uncontrolled when the weather conditions, topography, and fuels, the drivers of fire behavior, are suitable for fire spread or propagation [7,8]. From the viewpoint of ecosystem succession, a burnt forestland will recover over a long-term period of succession and consequently be suitable for the development of new communities of vegetation and wildlife [9]. However, large-scale and high-intensity fires are likely to occur simultaneously over a landscape and therefore cause profound societal impacts [10]. The most recent example of such extreme natural disturbance would be the 2019–2020 Black Summer megafires in Australia [11]. In contrast to prescribed fires, the biomass burning events release a significant amount of particular matter into the air [12]. Uncontrolled forest fires can be recurrent and cause a significant loss of aboveground biomass stocks. Observations of forest fires in northern Brazil showed the tropical rain forests lost more than 60% of the biomass stocks in 3–7 years after the last major fire in a series of recurrent fires; even for large trees with diameter > 50 cm, the biomass loss was about 54% in areas that burnt three times [13]. Similarly, the carbon loss of tropical forest in northern Australia caused by surface fires was equivalent to 46% of the annual net primary production (NPP) of the forest [14]. Although smaller carbon-stock reductions of about 6% of annual NPP were reported for temperate forests versus 11% for natural vegetation landscapes, wildfire-caused biomass consumption has become a significant source of carbon emissions globally [15]. Management of fuel loads therefore becomes an important issue for the safe use of fire, wildfire prevention, and REDD achievement conservation [16].

Fuel is the combustible biomass found in forests and can be divided into fine fuels such as leaves, grasses, and small twigs, and larger fuels such as shrubs, branches on the ground, downed trees, and logs [17]. From the standpoint of vertical dimension, the fuel can be classified as ground fuels, aerial fuels (trees, snags, and ladder fuels), and canopy fuels (green leaves and branches of crowns). Accordingly, a fire occurring in the humus layers, moving slowly, which can probably smolder for a long time, is called a ground fire; a fire burning only surface litter and duff is a surface fire, while a fire that burns trees over their entire height to the top is called a crown fire. Human-caused wildfire generally begins as a ground/surface fire at a small scale but occasionally it can become a crown fire and cause huge damage to the forest. Therefore, a map of fuel load distribution can highlight the prevalence of wildfire risk and guide people in more careful use of fire in order to minimize the possibility of heat transfer via conduction/radiation/convection in forests.

When considering the possibility of ignition of a forest fire, identifying fuels that are most likely to burn is related to the water content and the amount of surface fuels in a forest. In other words, the fuel loadings should include the quantities of duff, litter, fine-woody debris, and coarse woody debris or logs (fallen dead woods) distributed over the ground surface of the forest because these are the primary factors for predicting fire effects from on-site fuels [18]. The fuel-loading models (FLMs) designed by Sikkink et al. [18] emphasized the importance of fuel composition for ground and surface fires. When fires become uncontrollable, the amount of crown materials will be included in the fuels available for crown fires and the trunk may or may not be burnt in crown fires. Therefore, many studies have been conducted to explore methods for deriving the distribution of fuel or biomass (ton/ha) using a variety of data such as forest inventory data and airborne lidar scanning (ALS) data for aboveground biomass [19,20,21,22] and canopy biomass [23,24,25]. According to the remote sensing-based IPCC method, a canopy fuel map can be derived using an ALS canopy height model by segmenting every single tree, using, for example, mathematical morphology-based watershed segmentation [26,27,28,29], Multilevel Morphological Active Contour (MMAC) [30] or Multilevel Slicing And Coding (MSAC) techniques [31]. Moreover, the latest development of lidar sensing enables precise inventories of surface fuel and canopy fuel using mobile terrestrial lidar instruments [32,33,34,35]. To overcome the high cost of high-density point cloud ALS data, alternate methods for surface fuel loading (SFL) estimation can be based on mathematical/empirical models using inventory data such as vegetation/species maps and related environmental factors [18,36], satellite full-waveform lidar data [37,38], and photon lidar data [39]. More recently, the approach has been extended to integrate optical sensing images such as infrared orthophotos and QuickBird images with ALS data to improve fuel mapping accuracy [40,41,42,43]. The major strength of lidar technology in fuel estimation is the ability to retrieve fuel heights and discriminate between fuel types [36]. It is impossible to measure surface fuel loads directly from lidar data because lidar pulses can rarely penetrate the dense litter and duff layers on the bare earth surface and travel back to the sensor.

Surface fuel inventories involve the complicated task of the collection of live biomass such as grasses, forbs, and small woody plants and dead fuels such as duffs, debris, and fallen dead wood over a large area of forest floor. Lidar technology is considered to be an efficient method of gathering information of detailed biomaterials distributed on the land surface. However, there still are difficulties in the determination of surface fuel loads over a wide range of forests using this technology. In practice, the determination of surface fuel loads is a process related to the collection and analysis of diverse surface fuels, fuel bed depth, and bulk density in the field. Fuel depth varies widely, for example, from 0.3 m for grasses to 1.8 m for shrub fields and 0.06–0.30 m for timber litter in forests [44]. Therefore, the surface fuels lying on the fuel bed 0.3 m above the earth surface are generally collected to account for the bulk density [45]. In addition, the ground is generally covered with surface mass due to the weathering effect and accumulation of fine materials. In spite of the ability of laser pulses to penetrate canopy gaps and reach the ground, collecting sufficient numbers of point clouds that lie on the surface fuel and the earth surface remains difficult and therefore, characterizing surface fuel though lidar point clouds is still challenging.

Surface fuel loads and bulk density are primarily subject to forest structures related to species composition, phenology, and canopy height [41,46,47,48,49,50]. Distribution of surface fuel loads is therefore a consequence of the interaction of multiple factors such as forest type or overstorey/understorey species, topographic relief, and climate. A traditional forest fuel inventory is capable of collecting the loads and bed depth of surface fuel and can further differentiate and measure a variety of fuel sizes, for example, duff mass, litter mass, fine-woody debris, coarse-woody debris, and fallen dead wood (FDW). In practice, a field inventory of surface fuel loads within the whole area of sample plots has to collect and measure the masses of the litter layer and duff layer. The work is time consuming and labor intensive and significantly disturbs the surface stability and seed bank when the plot covers a large area, for example, 10 by 10 square meters or even larger. The distribution of fuel masses is most likely spatially dependent in a relatively local space. This kind of spatial autocorrelation can be described via geostatistics which incorporates the spatial coordinates of sample data to interpolate values for locations where samples were not taken. Geostatistical methods such as kriging and cokriging (refer to Section 2.2 for the details) capture observed spatial dependence among regional points. In addition, the distribution of forests is typically a consequence of natural processes in which terrain morphology such as the slope, aspect, and elevation may influence the amount of solar radiance, temperature, humidity, and water availability on the slope, and these further affect the weathering and accumulation of surface fuel mass on a slope. Thus, the integration of traditional inventory and geostatistical methods could be helpful for mapping fuel distribution. A geospatially explicit continuous map of surface fuel loads can be a fuel baseline for the use of fire protection scenario generation and forest management. Therefore, the objective of this study was to propose an algorithm for generating surface fuel load maps through the integration of forest types, topographic variables, and in-situ inventory mass data using geostatistical analysis and multiple linear regression methods. The surface fuel load was presented as two types: fSFL (the fine- and coarse surface fuel loads) and tSFL (the total surface fuel load, i.e., the sum of fSFL and the FDW mass). With detailed composition complexity of fSFL and tSFL distribution over the study site, the uncertainty in mapping the surface fuel load was further examined, and a strategy to generate the surface fuel load was suggested.

2. Materials and Methods

2.1. Study Site and Data Acquisition

The forest in the area of the Dajiaxi Working Circle, managed by the Taiwan Forestry Bureau, was selected for this study. The forest located at 121°07′42″–121°27′03″E and 24°08′09″–24°26′31″N in north central Taiwan (Figure 1) has been frequently disturbed by fires during the last few decades. According to the official records revealed on the website of the Forestry Bureau (https://forecast.forest.gov.tw/Forecast/# accessed on 10 February 2021), a total of 2171 fires occurred in the 37 working circles of national forests from 1963 to 2019, averaging 37 ± 39 fires per year, with a minimum and maximum of three and 252 fires, respectively, while averaging 55 ± 74 per working circle in which the minimum frequency was 2 and the maximum was 287. Historical records of the Dajiaxi Working Circle revealed that larger fires occurred frequently during the period from 1963 to 1990. The annual occurrence frequency ranged from 0 to 20, and the burnt area was on average 526 ± 717 ha per year. Due to the implementation of a fire-fighting approach that incorporates the Incident Command System and Government Flying Service, the fire frequency and burnt area, after 1990, was significantly reduced to 0 to 4 and 39 ± 84 ha per year, respectively. Most of the fires occurred in the dry season from winter to the early spring; during this period, fire is most likely caused by careless use of fire. The prevalence of forest fire in the Dajiaxi Working Circle was around 10% and ranked 3rd among the 37 forest management units. The Dajiaxi national forest is therefore officially considered as a hotspot area for fire.

The altitude of the Island of Taiwan ranges from 0 to 3950 m. The forests on the island vary dramatically along with the changes in altitude, temperature, and latitude. With respect to the altitudinal variation, the forest is divided into the foothill zone (tropical forest), submontane zone (subtropical forest), montane zone (warm-temperate and temperate forests), upper montane (cool-temperate forest), subalpine zone (cold-temperate forest), and alpine zone (subarctic forest) [51]. The elevation of this site ranges between 1115 and 3885 m across the subtropical–temperate–subarctic forest zones, and the temperate forest dominates this area. The forest types include conifers, mixed pine–conifer–broadleaf (hereafter mixed), and broadleaf forests. Specifically, this site has a lot of Taiwan red pine (Pinus taiwanensis) plantations which were originally managed for wood production and thus account for a major part of the coniferous forest even though it has been suffering from a high risk of forest fire for decades. The pine is therefore listed together with the conifer, mixed, and broadleaf forests as one of the forest types for fuel load inventory and SFL modeling. Correspondingly, each of the forest types was sequentially encoded as 1: Taiwan red pine, 2: conifer, 3: mixed, and 4: broadleaf in this study. The map of forest types was generated using high-resolution ortho-photos obtained from the Taiwan Forestry Bureau.

The airborne lidar scanner data were acquired on 14 December 2018 via Strong Engineering Consulting Company using a P68C-TC plane. A small-footprint, full-waveform lidar system (Riegl LMS-Q780) mounted on the aircraft provided high-accuracy point cloud data. Lidar data were collected at an operating flight altitude of 3400–4000 m (or 1970–2370 m above ground level) with a laser pulse repetition rate of 240–270 KHz. The resulting lidar dataset with ground and canopy point cloud density around 2.5 and 15 points per square meters, respectively, was used to produce a 1.0-m cell resolution of a rasterized digital elevation model (DEM) and digital surface model (DSM) using a linear interpolation technique. Both DSM and DEM were used to produce CHM data for aboveground biomass mapping in a previous study, while only DEM data were used to derive topographic parameters such as elevation, slope, and aspect for this study. The range of the DEM data, classified degree slope (CS), and classified degree aspect (CA), as well as the reciprocal of degree slope (RDS) are shown in Figure 2. Forest type and topographical variables are abbreviated as BioTopo variables hereafter.

In implementing airborne laser scanning, aerial photographs with 80% endlap and 50% sidelap were acquired using a PhaseONE iXA 180 camera. The geometric distortion of the aerial photography had a value of 0.0169°, 0.0209°, and 0.0188° for the yaw, pitch, and roll error, respectively, which accounted for an overall error of 0.0328°. The aerial photographs were used to generate a 0.2-m cell resolution orthophoto whose x-, y-, and z-coordinates had RMSE values of 5.05, 2.45, and 1.36 cm, respectively, accounting for an overall RMSE of 5.78 cm, based on the RTK-based ground control points. With the high-resolution DEM and orthophotos of the study site, a series of route planning was implemented in advance for in-situ inventory.

2.2. Surface Fuel Load Inventory Using Stratified Random Sampling and Cokriging Analysis

This study designed a three-level sampling strategy (Figure 3) to collect sufficient SFL sample data. The sampling process was first implemented using forest type as the strata and basemap-ID as the sampling unit. A series of random numbers indicating a sampling unit was selected in which a level-1 plot with an area of 20 m by 20 m was drawn for the pine, conifer, mixed, and broadleaf forestsover the study site (Figure 4). Second, a level-1 plot was divided into four subplots (level-2) with a size of 10 m by 10 m. Third, every subplot was evenly subdivided into 1-m-gridded microplots (the size is identical to topographical variables derived from ALS DEM data), and the particular microplots (level-3) located at the four corners and the center of each level-2 subplot whose surface fuels including dull, litter, fine-woody debris, coarse-woody debris, and fallen dead wood were investigated. Fourth, an ordinary cokriging method (Equation (1)) was applied to derive a distribution map of fSFL over the level-1 plot in the 1-m gridded cells, and finally, the fSFL of a level-2 subplot was determined by aggregating the values of all 1-m cells within the area of the corresponding subplot. After that, the tSFL of a subplot was determined by summing up the fSFL and the observed FDW mass within the corresponding area.

The in-situ field inventory was carried out between February and July in 2019. Virtual base station real-time kinematic (VBS-RTK) technology was applied to locate the inventory plots in the forest. In addition, a GeoSLAM terrestrial mobile lidar scanner ZEB Horizon was used to capture 3D data for the plots for another aboveground biomass study (Figure 5). Live and dead surface fuels, including 1-hr, 10-hr, 100-hr, and 1000-hr timelag fuels, were collected separately [18,36,52]. For each fuel category, the fresh (or wet) weight and absolute dry weight were measured in-situ and in the laboratory, respectively. For each category, a fuel sample of 500 g was sent to be oven-dried at a temperature of 105 °C for measuring the absolute dry weight [53]. With the dry­–fresh weight ratio of each subsample, the value of surface fuel loads of every plot was determined.

Under normal circumstances, surface fuel distribution is mainly governed by the location of vegetation/trees, which are distributed randomly/systematically over the forestland area in natural/planted forests, and is most likely related to topography, particularly the slope and aspect. This is because litterfall as well as fallen dead logs naturally move downslope due to gravity and collect at a position where the movement is blocked by topography or objects. In addition, when the fallen dead wood decomposes/decays, the rolling debris will accumulate more at the bottom of slopes over a limited local space. In other words, the amount of fine surface fuel or dead biomass can gradually change in terms of the physical characteristics over a slope, and fallen dead wood can occasionally interrupt the continuity in stands [54]. The possibility of discontinuity/continuity of surface fuel distribution increases with an increase in the area of interest of factors such as powerful typhoons or tropical cyclones that frequently bring high winds and rainfall, causing biomass to move to lowland areas, inland lakes or the ocean. A plot-based forest inventory is generally area limited, and attributes of the trees and the ground surface within a plot tend to be homogeneous; thus, the surface fuel properties of points in a plot are supposed to be spatially dependent. Because cokriging analysis conducts interpolation according to a semivariogram model, the technique can determine the spatial dependence among points [55,56]. In cokriging analysis, the fSFL in the level-3 1-m cells within a level-1 plot can therefore be presented as a function of the observed amount of surface fuel (the primary variable), with the slope degree, degree aspect, and fuel bed depth the three secondary variables. The cokriging system containing one primary and three secondary variables is defined as:

$$z_0^{*}={\displaystyle \sum }_{i=1}^n{\alpha }_i{z}_i+{\displaystyle \sum }_{j=1}^m{\beta }_j{x}_j+{\displaystyle \sum }_{k=1}^p{\gamma }_k{v}_k+{\displaystyle \sum }_{l=1}^q{\tau }_l{u}_l$$

(1)

where $z_0^{*}$ is the estimate of Z at location 0; z₁, z₂, …, z_n are the primary data at n locations; x₁, x₂, …, x_m, v₁, v₂, …, v_p, and u₁, u₂, …, u_l are the secondary data at m, p, and q locations, respectively. α₁, α₂, …, α_n, β₁, β₂, …, β_m, γ₁, γ₂, …, γ_p, and $l_1, l_2, \dots , l_q$ are cokriging weights to be determined.

In the cokriging analysis, a variety of semivariogram models such as exponential, circular, spherical, and Gaussian models was tested and evaluated using the cross-validation method. Only the semivariogram model fitted with a smaller RMSE was used to generate the fSFL map of a plot. The microplot fSFLs over a level-1 plot area were further used to aggregate the level-2 subplot-based fSFL values and tSFL.

2.3. Modeling of Surface Fuel Loads Using Multiple Linear Regression

In this study, the regression coefficients (B) of a multiple linear regression (Equation (2)) were estimated by the method of ordinary least squares using Equation (3).

$$Y=XB+\in$$

(2)

$$\widehat{\mathbf{B}}={\left(X^{\prime }X\right)}^{-1}{X}^{\prime }Y$$

(3)

In Equation (3), X is the matrix of independent variables including CS, CA, RDS, and NFT (the normalized forest type, which is determined as the ratio of FT to the maximum value of FT) and Y is the dependent vector lnSFL. X’X = R is the correlation matrix of independent variables as each of them is standardized by its own mean and standard deviation. In this study, 120 level-2 plots were collected and randomly divided into training and assessing sub-datasets. Based on leave-20%-of-the-plots-out (5-fold) cross-validation [38], all of the sample plots were evenly grouped into five assessing datasets. As a result, the average and standard deviation of performance measures were determined. In order to explore if estimation bias was related to fuel type, an additional evaluation was implemented based on leave-1-fuel-type-out cross-validation.

In the regression analysis, two surface fuel load models were generated, that is, the fSFL-BioTopo model and tSFL-BioTopo model in which fSFL and tSFL represent the fine surface fuel loads and total surface fuel load, respectively, and BioTopo is associated with the biological and topographic variables. The fSFL and tSFL models derived from the training dataset were further applied to derive a distribution map of the fSFL and tSFL for the Dajiaxi National Forest. Accuracy of the fSFL and tSFL maps was evaluated by a cross-validation method via the root-mean-square error (RMSE) and the percentage root-mean-square error (PRMSE). The RMSE is a scale-dependent accuracy measure and is presented on the same scale as the surface fuel load; in contrast, the PRMSE is scale-independent and measures the accuracy as an error percentage relative to the average of observations [22].

3. Results

3.1. The Derived fSFL Semivariogram Models and Their Performance in Estimating Level-1 Plot Surface Fuel Loads

Based on the rule of the smallest estimation bias in the cross-validation of a cokriging method in the level-3 microplot, a semivariogram model with the smallest RMSE was applied to generate a 1-m cell SFL map for the corresponding level-1 plot. Detailed information of the fitted models, prediction maps, RMSEs, and the percentage of RMSE related to the mean average of fSFL observations (PRMSEs) in the inventory data is shown in

Table 1. The ordinary cokriging method-derived maps of 1-m cell fine surface fuel loads (fSFL) and for the inventory plots of forest types.

Forest Type	Plot 1	Plot 2	Plot 3	Plot 4	Plot 5	Plot 6	Plot 7	Plot 8
Pine
Model *	Gaus (1)	Exp (2)	Exp (3)	Exp (3)	Exp (4)	Spher (1)	Exp (4)	Exp (4)
RMSE	0.3036	0.3158	0.1604	0.1999	0.1762	0.1809	0.2667	0.0902
PRMSE	26.31	28.71	15.31	20.91	16.94	19.54	27.61	4.83
Conifer
Model *	Exp (4)	Exp (4)	Exp (4)	Exp (1)	Exp (1)	Exp (1)	Exp (2)	Exp (4)
RMSE	0.2428	0.2594	0.2023	0.1687	0.3012	0.1423	0.2750	0.2387
PRMSE	32.81	31.18	18.46	29.91	43.15	22.80	39.51	35.42
Mixed								x
Model *	Exp (3)	Exp (5)	Exp (2)	Cir (4)	Exp (4)	Exp (3)	Exp (6)	x
RMSE	0.6082	0.7723	0.2239	0.2777	0.4112	0.2295	0.1988	x
PRMSE	39.14	50.95	36.58	35.06	31.29	30.68	48.88	x
Broadleaf								x
Model *	Exp (2)	Exp (4)	Exp (1)	Exp (5)	Exp (4)	Exp (4)	Exp (2)	x
RMSE	0.2229	0.1215	0.0934	0.2685	0.1332	0.1712	0.3411	x
PRMSE	24.44	17.21	21.13	45.05	23.13	29.12	43.61	x

*: Exp, Spher, Cir, and Gaus represent the exponential, spherical, circular, and Gaussian semivariogram models, respectively. The codes (1), (2), (3), (4), (5), and (6) after the type of semivariogram model indicate a combination of secondary variables used in deriving that semivariogram model. Correspondingly, the codes represent slope, slope–aspect, slope–fuel bed depth, fuel bed depth, slope–aspect–fuel bed depth, and aspect–fuel bed depth, respectively.

. As can be seen, the best prediction of fSFL for every single plot of the forest types was mostly achieved by an exponential semivariogram model. Out of 30 level-1 plots with the tested secondary variables, the variable slope was the most frequently used as the supplementary data to describe the spatial change in fSFL. Table 1 also shows 17 of the 30 models whose fSFL was predicted via the slope or simultaneously via the aspect and/or fuel bed depth. Fuel bed depth was another frequently selected variable which was used alone to account for the distribution of the amount of surface fuel in 11 models and another five models when combined with other secondary variables. In contrast, the aspect appeared to be additional supplementary data when accompanied with the slope.

Recall the bias of the fSFL prediction map: the best accuracy had an RMSE of 0.0902 kg/100 m² or an equivalent error rate of PRMSE = 4.83% for a pine forest while the lowest accuracy had an RMSE of 0.7723 kg/100 m² and an PRMSE of 50.95% for a mixed forest. On average, the method of integrating 3-level stratified random sampling and cokriging analysis was able to derive the amount of surface fuels at an RMSE of 0.2533 ± 0.1390 kg/100 m² or a PRMSE of 29.66 ± 10.64%. In addition, the SFL distribution within the area of every single plot showed quite a different pattern among plots of the pine, conifer, mixed, and broadleaf forests. The natural variation of fSFL in a variety of forest types and topographic features revealed the spatial heterogeneity of fSFL distribution in forests, indicating that the method proposed makes sense for gathering plot-based surface fuel loads with a low cost of labor and time for forest inventories.

3.2. The Level-2 Subplot-Based fSFL-BioTopo Models and Their Performance in Generating the fSFL Map of the Whole Forest

The amount of fine surface fuels of the 120 level-2 subplots in an area of 10 × 10 m² is shown in

Table 2. The aggregated amount of fine surface fuel loads (fSFL) in level-2 subplots.

Forest Type	Subplot ¶	Plot 1	Plot 2	Plot 3	Plot 4	Plot 5	Plot 6	Plot 7	Plot 8	AVG	STD
Pine (n1 = 32)	LL	149.75	132.98	75.89	161.72	121.82	63.54	92.59	79.22	107.32	24.56
	LR	122.12	120.37	125.91	126.11	139.38	101.47	87.41	83.19
	UL	139.70	120.52	73.71	119.99	82.82	84.83	86.04	109.47
	UR	120.34	110.68	128.56	78.36	97.3	112.78	76.40	109.42
Conifer (n2 = 32)	LL	89.93	90.90	113.00	91.88	76.28	55.86	66.12	60.31	81.34	22.90
	LR	95.32	76.19	135.27	46.10	76.17	40.48	82.05	48.35
	UL	93.78	85.37	110.91	68.92	99.70	82.24	90.00	92.86
	UR	85.60	56.18	138.06	49.73	67.62	68.18	84.53	84.99
Mixed (n3 = 28)	LL	86.04	223.01	104.74	87.04	225.84	73.07	36.48	x	99.58	56.11
	LR	80.86	227.40	48.80	94.39	174.59	89.64	41.86	x
	UL	87.53	138.60	63.41	80.39	111.58	73.27	56.27	x
	UR	91.69	180.47	39.15	87.23	104.04	46.44	34.30	x
Broadleaf (n4 = 28)	LL	83.50	67.42	54.91	71.69	56.90	65.14	92.40	x	73.54	16.28
	LR	107.28	74.87	53.79	64.40	62.95	58.54	104.31	x
	UL	74.98	92.06	49.62	85.58	63.65	73.2	84.48	x
	UR	100.47	77.08	44.76	68.04	70.53	64.07	92.55	x

^¶: The abbreviations LL/LR/UL/UR indicate the level-2 subplot on the lower left/lower right/upper left/upper right of a level-1 plot listed in Table 1. The value of each entry has a unit of kg/pixel. A pixel has an area of 100 m². AVG and STD represent the respective average and standard deviation of the values for a forest type.

(hereafter a pixel) and was aggregated from the cokriging-derived level-1 SFL map. As can be seen, the pine stands had fSFL values that ranged from 73.71 kg/pixel to 149.75 kg/pixel, with an average of 107.32 ± 24.56 kg/pixel greater than 99.58 ± 56.11 kg/pixel, 81.34 ± 22.90 kg/pixel, and 73.54 ± 16.28 kg/pixel of the mixed stands, conifer stands, and broadleaf stands, respectively. The broadleaf stands on average had obviously smaller fSFL while the mixed stands, whose fSFL values displayed a significant and dramatic change indicated by their standard deviation, was almost three times larger than those of the pine, conifer, and broadleaf stands.

The ANOVA test showed that the derived fSFL-BioTopo model using the fSFL data of level-2 subplots displayed an R-squared value of 0.162 (F = 3.096, p < 0.005, n = 120). In Equation (4), the dependent variable is the natural log-transformed fSFL (lnfSFL, kg/pixel), and the independent variables are the original or first-order topographical variables (NFT, AC, and SC), their second-order interaction product (NFTxAC, NFTxSC, and ACxSC) and third-order interaction product (NFTxACxSC). Based on the cross-validation test, this model was able to achieve an accuracy of RMSE = 34.10 kg/pixel and PRMSE = 37.59%. In contrast, when the 6-class SC was replaced by four classes (that is, 1: <5°, 2: 5–10°, 3: 10–20°, 4: >20°) and the 8-class AC was regrouped as four classes (1: North, 2: East, 3: South, 4: West), the alternative lnfSFL-BioTopo model (Equation (5)) displayed an R-squared value of 0.173 (F = 8.063, p < 0.001, n = 120) and had an RMSE = 33.07 kg/pixel and PRMSE = 38.03%. The RMSE and PRMSE differences between Equations (4) and (5) were 1.03 kg/pixel and 0.44%, respectively, which accounts for a relative change rate of 3% and 1% in the RMSE and PRMSE. It was therefore concluded that the performance of the two models was almost identical.

$$\begin{gathered} \mathrm{lnfSFL}=3.785174+0.819635*\mathrm{NFT}+0.154689*\mathrm{AC}+0.416842*\mathrm{SC}-0.230969*\mathrm{NFTxAC}-0.671540*\mathrm{NFTxSC} \\ -0.072628*\mathrm{SCxAC}+0.127985*\mathrm{NFTxACxSC} \end{gathered}$$

(4)

$$\begin{gathered} \mathrm{lnfSFL}=3.911091+1.209008*\mathrm{NFT}+0.076095*\mathrm{AC}+0.375311*\mathrm{SC}-0.237776*\mathrm{NFTxAC}-0.842360*\mathrm{NFTxSC} \\ -0.040685*\mathrm{SCxAC}+0.130257*\mathrm{NFTxACxSC} \end{gathered}$$

(5)

3.3. The Level-2 Subplot-Based tSFL-BioTopo Model for Total Surface Fuel Loading Estimation

The detailed information of total surface fuel loads for the 120 subplots is shown in

Table 3. The aggregated amount of total surface fuel loads (tSFL) in level-2 subplots.

Forest Type	Subplot ¶	Plot 1	Plot 2	Plot 3	Plot 4	Plot 5	Plot 6	Plot 7	Plot 8	AVG	STD
Pine (n1 = 32)	LL	149.75	145.32	75.89	161.72	125.38	70.93	92.59	109.61	114.03	24.99
	LR	122.12	132.71	125.91	126.11	142.93	108.86	87.41	113.58
	UL	139.70	132.86	73.71	119.99	86.38	92.22	86.04	139.86
	UR	120.34	123.02	128.56	78.36	100.86	120.17	76.40	139.81
Conifer (n2 = 32)	LL	89.93	90.90	113.00	91.88	77.19	55.86	66.12	114.54	88.23	26.70
	LR	95.32	76.19	135.27	46.10	77.07	40.48	82.05	102.59
	UL	93.78	85.37	110.91	68.92	100.61	82.24	90.00	147.10
	UR	85.60	56.18	138.06	49.73	68.53	68.18	84.53	139.22
Mixed (n3 = 28)	LL	90.18	224.96	104.74	87.04	226.88	83.38	37.88	x	102.27	55.71
	LR	85.00	229.36	48.80	94.39	175.62	99.94	43.26	x
	UL	91.67	140.55	63.41	80.39	112.61	83.58	57.67	x
	UR	95.83	182.43	39.15	87.23	105.08	56.75	35.70	x
Broadleaf (n4 = 28)	LL	83.50	67.42	54.91	72.18	56.90	88.72	92.40	x	76.98	16.63
	LR	107.28	74.87	53.79	64.89	62.95	82.12	104.31	x
	UL	74.98	92.06	49.62	86.06	63.65	96.78	84.48	x
	UR	100.47	77.08	44.76	68.53	70.53	87.65	92.55	x

^¶: The same as in Table 2. The italics is used to highigh the partiular plots.

in which the italic numbers indicate the fallen dead wood mass of that particular subplot. On average, the largest amount of FDW mass within the subplot area was found in the conifer stand, with an average of 27.57 ± 26.66 kg/pixel, followed in descending order by pine stands (13.42 ± 10.28 kg/pixel), broadleaf stands (12.04 ± 11.55 kg/pixel), and mixed forest stands (3.77 ± 3.44 kg/pixel). In contrast to the fSFL, the increasing amount of FDW mass in pine, conifer, mixed, and broadleaf forest stands was around 6.25%, 8.47%, 2.70%, and 4.68%, respectively; this indicates that the prevalence rate of FDW in the conifer and pine stands was significantly higher than that in the mixed and broadleaf stands.

The inclusion of FDW mass did not change the relationship of the total surface fuel loads among the forest types, that is, the pine stand had the highest amount of tSFL, followed by the mixed, conifer, and broadleaf stands. The significant variation in tSFL in the mixed stands remained significantly larger than that in the other forest types. The R-squared value of the derived lntSFL-BioTopo model as shown in Equation (6) was 0.144 (F = 2.701, p < 0.013, n = 120). The performance of this model in predicting tSFL of the whole forest had an RMSE of 35.02 kg/pixel, corresponding to a PRMSE of 36.57%. Similarly, Equation (7) shows the alternative lntSFL-BioTopo model for tSFL estimation using NFT, 4-classes SC, 4-classes AC, and their interaction product variables. The R-squared value of this model was 0.168 (F = 7.836, p < 0.001, n = 120), and RMSE and PRMSE were 33.81 kg/pixel and 37.85%, respectively. The performance measures of the two models were also quite close, with a difference in RMSE and PRMSE of 1.21 kg/pixel and 1.28%, respectively. In contrast to Equation (6), the relative change in the two measures of Equation (7) was 3% and 4%.

$$\begin{gathered} \mathrm{lntSFL}=4.481948+0.171296*\mathrm{NFT}+0.082297*\mathrm{AC}+0.153252*\mathrm{SC}-0.170669*\mathrm{NFTxAC}-0.364650*\mathrm{NFTxSC} \\ -0.041167*\mathrm{SCxAC}+0.092501*\mathrm{NFTxACxSC} \end{gathered}$$

(6)

$$\begin{gathered} \mathrm{lntSFL}=4.433072+0.635037*\mathrm{NFT}-0.012423*\mathrm{AC}+0.163312*\mathrm{SC}-0.145670*\mathrm{NFTxAC}-0.541648*\mathrm{NFTxSC} \\ -0.000879*\mathrm{SCxAC}+0.077361*\mathrm{NFTxACxSC} \end{gathered}$$

(7)

For the 5-fold cross-validation, the respective average RMSE and PRMSE were 34.57 ± 8.76 kg/pixel and 37.58 ± 6.11% for the fSFL and 35.36 ± 9.09 kg/pixel and 36.38 ± 5.98% for the tSFL. However, for the leave-1-fuel-type-out cross-validation, a significant difference in estimation performance among the four trials was observed. When the pine, conifer, mixed, or broadleaf plots were sequentially excluded from modeling and then used for validation, the PRMSE for the fSFL estimation was 36%, 34%, 62%, and 109%, respectively, but 32%, 35%, 62%, and 102% for the tSFL estimation (

Table 4. Summary of the performance of fSFL and fSFL models based on cross-validation.

Model ¶	fSFL Model R2	RMSE (kg/m2)	PRMSE (%)	tSFL Model R2	RMSE (kg/m2)	PRMSE (%)
Equation (4)/Equation (6)	0.162 (F = 3.096, p < 0.005, n = 120)	34.10	37.59	0.144 (F = 2.701, p < 0.013, n = 120)	35.02	36.57
Equation (5)/Equation (7)	0.173 (F = 8.063, p < 0.001, n = 120)	33.07	38.03	0.168 (F = 7.836, p < 0.001, n = 120)	33.81	37.85
DeGroup 1	0.154 (F = 2.295, p = 0.034, n = 96)	23.28	28.56	0.128 (F = 1.844, p = 0.089, n = 96)	24.10	28.84
DeGroup 2	0.167 (F = 2.526, p = 0.020, n = 96)	25.82	33.75	0.164 (F = 2.469, p = 0.023, n = 96)	25.73	31.44
DeGroup 3	0.182 (F = 2.801, p = 0.011, n = 96)	46.70	46.30	0.145 (F = 2.128, p = 0.049, n = 96)	47.96	45.59
DeGroup 4	0.136 (F = 1.986, p = 0.066, n = 96)	39.25	37.91	0.120 (F = 1.713, p = 0.116, n = 96)	39.05	36.11
DeGroup 5	0.193 (F = 3.016, p = 0.007, n = 96)	37.82	41.36	0.188 (F = 2.908, p = 0.009, n = 96)	39.94	39.90
DePine	0.106 (F = 1.352, p = 0.237, n = 88)	38.47	35.84	0.051 (F = 0.609, p = 0.747, n = 88)	37.01	32.45
DeConifer	0.232 (F = 3.455, p = 0.003, n = 88)	27.82	34.20	0.240 (F = 3.607, p = 0.002, n = 88)	30.70	34.79
DeMixed	0.257 (F = 4.160, p = 0.001, n = 92)	61.56	61.82	0.246 (F = 5.231, p < 0.001, n = 92)	63.02	61.62
DeBroadleaf	0.342 (F = 6.242, p < 0.001, n = 92)	80.11	108.94	0.258 (F = 4.167, p = 0.001, n = 92)	78.89	102.48

^¶: The specific models (Equations (4) and (6)) used for generating surface fuel load maps of the whole study site. Performance measures of the models were determined based on all of the plots. DeGroup 1–5 represents the five evaluations of leave-20%-of-the-plots-out cross-validation; the respective average RMSE and PRMSE were 34.57 ± 8.76 kg/pixel and 37.58 ± 6.11% for the fSFL and 35.36 ± 9.09 kg/pixel and 36.38 ± 5.98% for the tSFL. DePine, DeConifer, DeMixed, and DeBroadleaf represent the four evaluations of leave-1-fuel-type-out cross-validation; the respective average RMSE and PRMSE were 51.99 ± 20.31 kg/pixel and 60.20 ± 30.20% for the fSFL and 52.40 ± 19.51 kg/pixel and 57.83 ± 28.21% for the tSFL. Both cross-validations included six slope classes and eight aspect classes.

). Obviously, the mixed and broadleaf surface fuel loads tended to be significantly over-estimated if they were not included in modeling. Both fSFL and tSFL of the broadleaf stands appeared to be significantly over-estimated because the surface fuel load in the stands was lower but under-estimated for the mixed stands due to the evidently diverse surface fuel load in the stands. In other words, a smaller PRMSE only occurred in the estimation when the plots of pine or conifer were not used for modeling, and the prediction bias appeared to be related to fuel type. Collecting appropriate numbers of plots for each fuel type for deriving a general model or having a sufficient number of plots for deriving fuel type-specific models should be able to prevent the significant over-estimation and/or under-estimation problem.

4. Discussion

4.1. The Uncertainty of Surface Fuel Loading Estimation in fSFL and tSFL Models

Because the relative error in the two fSFL-BioTopo models (Equations (4) and (5)) was small, and that of the two tSFL-BioTopo models (Equations (6) and (7)) was almost identical, estimation performance of the paired models for both fSFL and tSFL can be considered equal. The uncertainty of surface fuel models is therefore discussed based primarily on the estimation of Equations (4) and (6).

The predicted values of surface fuel loading over the whole area of the study site are shown in Figure 6 and are summarized in

Table 5. A summary of surface fuel loadings with respect to forest types in Dajiaxi National Forest.

Models	Forest Types	Areas (ha)	Minimum (ton/ha)	Maximum (ton/ha)	Average (ton/ha)	STD (ton/ha)	Total (tons)
lnfSFL-BioTopo	Pine	13,070	1.42	18.44	10.67	1.72	139,445.59
(Equation (4))	Conifer	14,039	1.04	13.33	9.29	1.10	130,433.75
	Mixed	7001	1.02	12.90	8.22	1.53	57,555.40
	Broadleaf	7280	0.66	13.96	7.18	2.40	52,283.57
	Sum	41,390	0.66	18.44	9.17	2.08	379,718.31
lntSFL-BioTopo	Pine	13,070	1.28	11.95	9.61	1.01	125,665.94
(Equation (6))	Conifer	14,039	1.03	11.62	8.81	1.03	123,659.56
	Mixed	7001	1.06	12.81	8.40	1.49	58,835.28
	Broadleaf	7280	0.76	14.32	7.71	2.25	56,147.02
	Sum	41,390	0.76	14.32	8.80	1.55	364,307.80
Difference	Pine	13,070	−6.50	1.34	−1.05	1.04	−13,779.65
(fSFL–tSFL)	Conifer	14,039	−3.34	0.83	−0.48	0.45	−6774.19
	Mix	7001	−2.52	0.79	0.18	0.26	1279.88
	Broadleaf	7280	−2.19	1.10	0.53	0.32	3863.45
	Sum	41,390	−6.50	1.34	−0.37	0.89	−15,410.51

. Based on the prediction maps of the whole forest, the fSFL mass of the pine stands ranged from 1.42 to 18.44 ton/ha and averaged 10.67 ± 1.72 ton/ha. Taking into account the forest areas, there was approximately 379,718.31 tons of fine surface fuel within the pine stands. This is the largest amount of fine fuel mass among the forest types over the whole forest. In contrast, the fine fuel mass of the conifer, mixed, and broadleaf stands averaged 9.29 ± 1.10 ton/ha, 8.22 ± 1.53 ton/ha, and 7.18 ± 2.40 ton/ha, respectively, and resulted in a total mass of 130,433.75 tons, 57,555.40 tons, and 52,283.57 tons. The relative amount of fine surface fuel mass among the four forest types derived from the lnfSFL-BioTopo model is quite similar to that observed in the subplot values. A similar trend appeared in the total surface fuel map derived from the lntSFL-BioTopo model. In addition, the amount of tSFL of the forest types showed the same sequence of pine > conifer > mixed > broadleaf, with an average value of 9.61 ± 1.01 ton/ha, 8.81 ± 1.03 ton/ha, 8.40 ± 1.49 ton/ha, and 7.71 ± 2.25 ton/ha, respectively. The results show that the lnfSFL-BioTopo and lntSFL-BioTopo models are capable of performing fSFL and tSFL estimations, and the estimates are generally consistent with the sampling inventory results. This reveals that the distribution map of surface fuel loading generated by each of the models is able to provide a baseline for accounting for the accumulation of surface fuel loads over time.

As noted by comparing the descriptive statistics of the two models in Table 5, the mixed and broadleaf fSFL estimate was generally smaller than the tSFL estimate by 0.18 ton/ha and 0.53 ton/ha while the pine and conifer stands’ fSFL was generally greater than the tSFL by an amount of 1.05 ton/ha and 0.48 ton/ha. Mathematically, the situation of fSFL > tSFL in a forest stand should not happen according to their definitions as described in Section 2.2. Recall the descriptive statistics of inventory data shown in Table 2 and Table 3: the fallen dead wood mass in the pine and conifer stands was almost three times higher than that in the mixed and broadleaf stands. In view of the range of surface fuel loading estimates of the models, the tSFL of pine and conifer stands was apparently smaller than the fSFL, with values of 10.66 vs. 17.02 and 10.60 vs. 12.29. In contrast to the mixed and broadleaf stands (11.75 vs. 11.87 and 13.56 vs. 13.30 for the range of fSFL and tSFL estimates), a significant uncertainty occurred in the estimation of the lntSFL-BioTopo model, and the source of uncertainty should be the inclusion of fallen dead wood mass. This kind of estimation uncertainty was also observed in the estimation bias of the subplots (Figure 7) in which an extra amount of bias in the estimates of tSFL was highlighted by an arrow with respect to those corresponding hollow bars.

4.2. The Dependency of Estimation Bias on the Amount of Surface Fuel Loads

In Section 3, the cross-validation tests showed that fSFL and tSFL were predicted with a similar accuracy of RMSE (33.84 kg/pixel vs. 34.76 kg/pixel) and PRMSE (37.29% vs. 36.28%), indicating that the lnfSFL-BioTopo and lntSFL-BioTopo models have almost the same ability to predict surface fuel mass in the forest. The similarity was revealed through the bias of the inventory subplots as shown in the bar chart in Figure 7. However, the difference in fSFL and tSFL in a one-hectare areal basis over the whole forest of the study site showed more than 50% of the areas whose fSFL was larger than tSFL. This is particularly evident in the pine and coniferous forests (Figure 8).

Although the estimation appeared to be bias-independent and randomly based on the scatter plot of bias vs. estimates in Figure 9a,b, the prediction bias still revealed a linear dependency on the observed value. This is evident in Figure 9c,d, where the original variable was converted to the deviation from the mean of observed values, i.e., fSFL-fSFL_avg or tSFL-tSFL_avg. The prediction bias can be presented as a negative linear function of the transformed variable, i.e., the bias is most likely to be compensated by the fuel mass itself, and the bias-adjustment value or compensatory value y_com can be determined by the linear models shown in Figure 9c,d. The compensated fSFL and tSFL estimates can be retrieved by subtracting the y_com from the original estimates of fSFL and tSFL using Equations (4) and (6), respectively. Accordingly, the compensated estimates of fSFL and tSFL through the bias-adjustment models in Figure 9c,d were significantly improved to 11.64 kg/pixel and 12.84% and 11.37 kg/pixel and 11.87% for RMSE and PRMSE, respectively. The bar chart shown in Figure 10a,b demonstrates the improvement of prediction bias for the original estimation and the compensated estimation of fSFL and tSFL with respect to each of the subplots.

4.3. A Possible Strategy for Improving Surface fuel Load Mapping

The prediction bias of the lnfSFL-BioTopo model and the lntSFL-BioTopo model (Equations (4) and (6)) can also be presented as a nonlinear function of the observed value of surface fuel loads. Estimates derived from the regression models can be over- or under-estimated and correspondingly generate a positive or negative prediction bias, determined as $\widehat{y}-y$. Each bias can be adjusted to positive by introducing an additive component, c, to compensate for the negative values without changing the relationship between the bias and the observed values. Assuming that the compensated offset value c is ≥ the maximum observed value of tSFL, the reciprocal transformation of “fSFL_bias + c” and “tSFL_bias + c” can be presented as an exponential growth function of the observed value of fSFL or tSFL, respectively (Figure 11). The transformed bias is helpful to diagnose the prediction bias behavior with respect to the original scale of the fSFL and tSFL observed values. As shown in Figure 11a,b, the R-squared value of the two exponential growth models was 0.8872 and 0.8713 when a constant value of 250 was assigned to the compensative offset.

In the multiple linear regression models (Equations (4) and (6)), a larger surface fuel load tended to be underestimated while a smaller fuel load was most likely overestimated, revealing that a smaller range of surface fuel load estimates was made by the lnfSFL-BioTopo model and the lntSFL-BioTopo model. Similarly, based on the reciprocal transformation, a smaller surface fuel load was generally found to have a smaller value of transformed bias and vice versa. The exponential growth function in Figure 11a,b shows that the greater the surface fuel load, the more significant the bias in the estimation. This is most likely induced by a shortage of samples of larger fuel loads in deriving the multiple linear regression models because the number of samples with respect to the diverse amounts of surface fuel loads is generally proportional to the size of the corresponding populations.

Figure 12 provides a generalized logistic distribution of the tSFL and fSFL observed values. The probability density function of this distribution is abbreviated as GL(x; α, β, γ for the shape, scale, and location parameters. Specifically, the distributions for tSFL and fSFL are GL(tSFL; 0.1804, 18.57, 90.04) and GL(fSFL; 0.2109, 17.38, 84.35), respectively. Figure 12a,b reveals a right skewed distribution, indicating the rareness of larger surface fuel load samples. In this study, the standard deviation and mean of the inventory data were 36.63 kg/pixel and 95.76 kg/pixel for the fSFL in the forests and 35.80 kg/pixel and 90.70 kg/pixel for the tSFL, which resulted in a coefficient of variation (CV) of around 0.38–0.39 for the forests. Statistically, the CV is used to examine the extent of data variability in relation to the arithmetic mean of a variable. The evidence of both the generalized logistic distributions of samples and CV suggests that increasing the number of inventory samples covering a wide range of a variable, particularly with a sufficient number of observations for diverse fuel loads, would be expected to upgrade a model’s estimation performance.

4.4. An Examination of the Appropriateness of the Cokriging-Based Surface Fuel Mapping Method

An additional independent inventory of surface fuel loads was carried out in an area of 10 × 10 square meters for examining the appropriateness of the sampling scheme in deriving the plot-based SFL map. The SFL of the whole plot was 100collected. As shown in

Table 6. A comparison of surface fuel load mapping using different numbers of samples within a spatial scale of a 10-m size plot.

Template No.	1	2	3	4	5	6
Locational template and the number of samples (NS) for deriving model ¶
	NS = 5@1 m	NS = 5@2 m	NS = 9@2 m	NS = 13@2 m	NS = 25@2 m	NS = 100@1 m
SFL Prediction map
Semivariogram model	Exponential	Exponential	Gaussian	Gaussian	Gaussian	Exponential
Secondary variables	slope, aspect, fuel bed depth	slope, aspect, fuel bed depth	slope, aspect, fuel bed depth	slope, aspect, fuel bed depth	slope, aspect, fuel bed depth	slope, aspect, fuel bed depth
RMSE (kg/m2)	0.59	0.65	0.80	0.69	0.40	0.05
PRMSE (%)	26.58	29.71	39.05	34.34	17.54	2.30

^¶ The black and gray boxes show the training and testing samples for cokriging model derivation and accuracy evaluation, respectively. The size value comes after the symbol “@” presenting the area of a microplot.

, a few samples of 2-m size microplots located in the black cells of the locational template (template no. 2–5, with 5, 9, 13, and 25 samples), were used to derive a cokriging semivariogram model for generating the SFL prediction map. The RMSE of the four models was 0.65, 0.80, 0.69, and 0.40 kg/m², which is equivalent to a PRMSE of 29.71%, 39.05%, 34.34%, and 17.54%, respectively, indicating the best accuracy was achieved by a geospatial cokriging model that used all samples from the whole plot area at the scale of a 2-m size sampling scheme. The difference in PRMSE between the fourth case and the first case was approximately 12%, indicating the proposed method to collect surface fuel loads of inventory plots is a viable approach.

In contrast, templates 1 and 6 show alternative sampling schemes at a 1-m scale. The cokriging model derived using template number 6 had a prediction accuracy of RMSE 0.05 kg/m² and PRMSE 2.30%; the predicted values were quite close to the observed values. Compared with template 2, the smaller values of both RMSE and PRMSE achieved by template 1 also reveal the appropriateness of the proposed method in reducing labor costs while retaining accuracy for generating surface fuel load maps. The 1-m size cell is therefore defined as the minimum sampling unit (MSU).

Table 7. A comparison of appropriateness of extended prediction based on the predefined cokriging surface fuel load map ^¶.

7	8	9	10	11	12	13	14
0.60/26.85%	0.59/26.30%	0.59/26.12%	0.59/26.35%	0.72/31.80%	0.79/35.42%	0.79/35.42%	0.89/39.92%
5 (CenLL)	5 (CenLR)	5 (CenUL)	5 (CenUR)	5 (LL)	5 (LR)	5 (UL)	5 (UR)

^¶ The black and gray boxes show the training and testing samples for cokriging model derivation and accuracy evaluation, respectively. The surface fuel load map for the particular area bounded by the training samples was first generated through the cokriging semivariogram model and then extended to the extent defined by templates no. 7–14. The xx/xx% values below each map specify the RMSE/PRMSE that was derived from the corresponding template.

further demonstrates the extrapolation map of surface fuel loads extended from a limited predefined map of the five sample microplots. The RMSE and PRMSE did not show significant differences among the four centralized sampling schemes (denoted as CenLL, CenLR, CenUL, CenUR), but the two measures showed an evident increase in the four edged sampling schemes. The results indicate that the extrapolation is not appropriate for deriving surface fuel loads of inventory plots, and the MSUs should be distributed along the boundary as well as in the center of a plot.

In published articles, much research demonstrates that the performance in the estimation of surface fuel loads using lidar data and optical images varies dramatically. From the perspective of PRMSE, the integrated approach of diverse remote sensing data was around 20–38% for a dense coniferous forest located in central Greece [33] and 37–98% for a bark beetle-affected forest in eastern Grand County in north-central Colorado [57]. A better accuracy of PRMSE around 5–47% was achieved for an upland oak-dominated forest in Kentucky using small footprint full-waveform lidar data [58]. According to Franke et al. [59], a PRMSE ranging from 21 to 41% was achieved when measuring diverse coarse woody debris in a forest savanna in the Brazilian Cerrado using only multi-temporal Landsat OLI images. In contrast, the 37% of PRMSE achieved in this study is quite close to the moderate accuracy revealed in previous studies, reflecting the potential of applying stratified random sampling to forest types, topographical variables, and inventory data in generating baseline information for surface fuel loading. The classification of forest types was determined based on the biological conditions of the study site. This system is quite similar to the fuel types of coniferous, deciduous, mixed wood, slash, and open grassland, as defined in the Canadian Fire Behavior Prediction System (FBPS) [36]. Various research has demonstrated the feasibility of integrating ALS and optical images to map the fuel types (alternatively fuel models), such as the ones defined by the Northern Forest Fire Laboratory (NFFL) [36,41,49,50,60,61]. The proposed algorithm for mapping the surface fuel load is therefore most likely able to substitute the fuel types of the FBPS, NFFL, and NFDRS classification systems to moderately improve the mapping performance for forests with undulating terrain morphology in mountainous area.

5. Conclusions

Surface fuel load estimation and mapping are of particular importance, not merely for the origin of fire and better prediction of fire spread and intensity, but also for understanding the source of soil nutrients. Decomposed vegetative mass can easily enter into nutrient cycling and support the needs of vegetation growth and forest development. An in-situ field inventory is the direct method to collect real data of surface fuel load in a forest, but the method is only implementable over limited or small areas in the view of costs of labor, time efficiency, and finance. In general, the surface fuel load over a forest area is a result of the accumulation of litterfall and dried and short-lived vegetation mass as well as fallen dead wood, and wildfire consumes surface fuel. A fire-behavior model can be used to chart the post-fire fuel dynamics when the dynamics of fuel accumulation are established based on the historical records of fire regimes and initial surface fuel load [62]. Empirical models with predictors derived from lidar data and/or satellite images provide alternative methods of estimation at a certain accuracy or uncertainty. In fact, surface fuel masses generally form a dense cover on the ground surface. It is impossible for lidar pulses to penetrate the fuel bed and reach the bare ground. Consequently, measuring the surface fuel load directly with lidar point cloud data from the air and even the ground is obviously quite challenging.

The amount of surface fuel mass will change over time due to diverse influences induced by the interaction of biological, physical, and climatological factors; therefore, the information revealed through a map of surface fuel distribution is most likely valid or practical only for a limited period. Frequent updates of fuel maps are required for efficient management of forest fuel in order to control fire risks. In view of the economic cost of fuel mapping, a method of deriving accurate surface fuel load maps is needed that will be both easier to implement and at a lower cost. Considering the complexity of undulating terrain morphology and inaccessibility of vehicles in mountainous areas, the proposed method for estimation and mapping of surface fuel loads using topographic variables and classified fuel models (forest types) is highly appropriate to meet this need. To implement the 3-level stratified random sampling based approach for surface fuel load mapping, the user should apply a fuel type (also fuel models) classification as needed and then carry out inventories to collect data for generating a map of the surface fuel load.

For deriving a reliable prediction of surface fuel loads of an inventory plot, it is recommended that the minimum sampling unit for collecting surface fuel should include the four corners and the central position of an inventory plot in order to allow the cokriging method to achieve an accurate prediction. In the proposed method, the orthogonal decomposition design of the plots was mainly for the convenience of linking with raster-based remote sensing data. The size of a plot area can be determined based on the pixel size of the fuel type map and the topographic map. The plot design can be flexible in the geometries and sizes that allow compatibility with a variety of forest inventory systems, for example, the forest inventory and analysis (FIA) plot design of the USDA that establishes three additional plots next to a core plot at a fixed distance and in three directions [63]. The plot design of the Indonesia National Forest Inventory is a systematic cluster design. It is composed of clusters, temporary sample plots, and permanent sample plots. The fuel load data collected via the Cluster-TSP-PSP plot design [64,65] can directly adopt the proposed approach.

Biological variables are the primary leading factors of the surface fuel load. Some of the factors are likely to change over time due to growth, competition among trees, and disturbances. To address the temporal-related changes, a spatiotemporal dynamic model of biological mass transition would be a critical solution for better surface fuel load management.