Assessing the Potential of the DART Model to Discrete Return LiDAR Simulation—Application to Fuel Type Mapping

: Fuel type is one of the key factors for analyzing the potential of ﬁre ignition and propagation in agricultural and forest environments. The increase of three-dimensional datasets provided by active sensors, such as LiDAR (Light Detection and Ranging), has improved the classiﬁcation of fuel types through empirical modelling. Empirical methods are site and sensor speciﬁc while Radiative Transfer Models (RTM) approaches provide broader universality. The aim of this work is to analyze the suitability of Discrete Anisotropic Radiative Transfer (DART) model to replicate low density small-footprint Airborne Laser Scanning (ALS) measurements and subsequent fuel type classiﬁcation. Field data measured in 104 plots are used as ground truth to simulate LiDAR response based on the sensor and ﬂight characteristics of low-density ALS data captured by the Spanish National Plan for Aerial Orthophotography (PNOA) in two different dates (2011 and 2016). The accuracy assessment of the DART simulations is performed using Spearman rank correlation coefﬁcients between the simulated metrics and the ALS-PNOA ones. The results show that 32% of the computed metrics overpassed a correlation value of 0.80 between simulated and ALS-PNOA metrics in 2011 and 28% in 2016. The highest correlations were related to high height percentiles, canopy variability metrics as for example standard deviation and Rumple diversity index, reaching correlation values over 0.94. Two metric selection approaches and Support Vector Machine classiﬁcation method with variants were compared to classify fuel types. The best-ﬁtted classiﬁcation model, trained with the DART simulated sample and validated with ALS-PNOA data, was obtained using Support Vector Machine method with radial kernel. The overall accuracy of the classiﬁcation after validation was 88% and 91% for the 2011 and 2016 years, respectively. The use of DART demonstrates its value for simulating generalizable 3D data for fuel type classiﬁcation providing relevant information for forest managers in ﬁre prevention and extinction.


Introduction
Fuel types are defined by Merrill and Alexander [1] as "an identifiable association of fuel elements of distinctive species, form, size, arrangement and continuity that will exhibit characteristic fire behavior under defined burning conditions." Fuel type mapping is crucial for forest management and fire risk assessment as the spatial distribution of fuel affects wildfire ignition and propagation.  In situ data measured from July to September 2014 in 104 field plots were used to adjust ground-truth to fuel type simulations. A stratified random sampling technique was applied to define the field plot location ensuring that it covers the range of terrain slopes and vegetation cover within the study area, which were estimated through ALS data. The Aleppo pine forests play an important role in the protection and recovery of forest in Mediterranean region characterized by nutrient-poor, gypsiferous soils, as it is the case of the site under study, since this species is practically the only one adapted to the adverse climatic and edaphic conditions of the area. The area presents a hilly topography, with altitudes ranging from about 400 m to 750 m a.s.l. The climate of the region is Mediterranean with continental features, characterized by irregular annual precipitation, cold winters and hot and dry summers [35].

ALS Data
The ALS data for simulating and validating the model is freely provided by the Spanish National Plan for Aerial Orthophotography (PNOA) through the National Center for Geographic Information (CNIG) (http://centrodedescargas.cnig.es). The ALS data were acquired in 2011 and 2016 with two slightly different acquisition specifications. The first coverage was captured in several surveys conducted between January and February 2011 with a Leica ALS60 sensor. The second campaign was conducted between September Remote Sens. 2021, 13, 342 4 of 20 and November 2016 with a Leica ALS80 sensor. Both sensors can record up to four returns per pulse and operate at a wavelength of 1064 nm. Data are delivered in 2 km × 2 km tiles of classified points in LAS binary file, format v. 1.2, with coordinate system in Universal Transversal Mercator (UTM) units, Zone 30, datum European Terrestrial Reference System 1989(ETRS 1989, (EPSG 25830). The flying height of the first and second ALS campaigns were around 3000 and 3150 m above ground level. Detailed information on the respective acquisition specifications are shown in Table 1.

Field Data
In situ data measured from July to September 2014 in 104 field plots were used to adjust ground-truth to fuel type simulations. A stratified random sampling technique was applied to define the field plot location ensuring that it covers the range of terrain slopes and vegetation cover within the study area, which were estimated through ALS data. The centroid of the circular plots (15 m radius) was positioned in the field using a Leica VIVA GS15 CS10 GNSS real-time kinematic Global Positioning System with an average accuracy of the planimetric coordinates of 0.33 m. The total tree height (h) and the green crown height were measured in all trees with a diameter at breast height (dbh) higher than 7.5 cm using a Vertex instrument for precise height measurement (Haglöf Sweden). Tree diameters were measured at breast height at the standard height of 1.3 m, using a Mantax Precision Blue diameter caliper (Haglöf Sweden). Additionally, the average height of the different shrub species and their coverage percentage at different height levels with respect to the plot surface were collected.
Fuel type was assigned using the Prometheus classification [36] that is based on the type, height and coverage percentage of the propagation elements [8]. This classification comprises seven categories ( Figure 2): one grass cover, three shrub covers with different mean heights (0 to 0.6 m, 0.6 to 2 m, 2 to 4 m) and three different tree covers (with no understory, with small understory, with understory connected to the base of the canopy).
Top-of-canopy reflectance measurements were acquired during the same field campaign using an Analytical Spectral Devices spectrometer (ASD FieldSpec 4 SR) in the 400−2500 nm spectral range (spectral resolution of 3-10 nm at Full Width at Half Maximum (FWHM) and a sampling interval of 1 nm). Reflectance was calibrated using a white Spectralon panel (Labsphere Inc., North Sutton, NH, USA) registered before every sample measurement. Official procedures of field spectrometry were applied to guarantee the quality of acquisitions [37,38] according to illumination conditions (clear days and close to the solar noon) and improvement of signal-to-noise ratio by subtraction of the dark current signal and spectrum average (25 measurements each). As a result, we obtained a total of 330 absolute reflectance spectra (with an average value of 5-10 different spectra for each species). tralon panel (Labsphere Inc., North Sutton, NH, USA) registered before every sample measurement. Official procedures of field spectrometry were applied to guarantee the quality of acquisitions [37,38] according to illumination conditions (clear days and close to the solar noon) and improvement of signal-to-noise ratio by subtraction of the dark current signal and spectrum average (25 measurements each). As a result, we obtained a total of 330 absolute reflectance spectra (with an average value of 5-10 different spectra for each species).

Simulation in DART Model
The scene generated in DART represents a plot with a flat surface of 30 × 30 m in order to resemble the field plot area (15 m radius). In DART, to constructing the scenes,

Simulation in DART Model
The scene generated in DART represents a plot with a flat surface of 30 × 30 m in order to resemble the field plot area (15 m radius). In DART, to constructing the scenes, different elements are provided such as plots with different characteristics (soil, soil+vegetation and vegetation) and trees. For more information on the scene components see Lamelas et al. [24] and Roberts et al. [25]. These elements are created using voxels that can be filled using turbid medium or facets. In our case a voxel of 0.5 m size was selected to represent the general scene, the plots were filled with turbid medium and the facets were used for the tree stratum. A voxel size of 1 m and the use of voxels instead of facets in trees were also tested obtaining worst results.
Different parameters are required to generate the plots and trees, related mainly to their size and structure, leaf area index (LAI), reflectance and transmittance values.
The simulation of the grasslands and shrublands was performed using plots and adjusted to the coverage percentage and height measured in the field in the case of those species with coverage percentages greater than 4% of the plot area, otherwise this surface was assigned to a species with similar characteristics and height. The trees were parameterized using the mean and standard deviation of the plot information. The shape of the Remote Sens. 2021, 13, 342 6 of 20 crown was selected from truncal to ellipsoidal according to the species to be represented, that is, truncal for genus Pinus and ellipsoidal for Quercus. It should be mentioned that the crown width variable had to be estimated from the diameter and height of the trees following the equations developed by Condés and Sterba [39].
The plots and the trees location follow a random spatial distribution since information on their location was not available. However, this was not considered a handicap since the methodology followed uses the vertical distribution in height of the returns and not their horizontal distribution.
The simulation requires the LAI values for the different species. In absence of field data, the LAI was estimated from two Sentinel 2-A top of canopy normalized reflectance data scenes (Level-2A images) using the Biophysical Processor tool integrated in the SNAP software [40]. The first image was captured on January 12th, 2016 and the second one on October 21st, 2016 ( Table 2). Pure pixels, characteristic of the different covers, were selected to extract the LAI value from the Level-2B Biophysical product. The time lapse in years from the simulation (2011) to the first available scene (2016) was not considered a handicap since the LAI value varies with a seasonal pattern and this was covered. In absence of laboratory reflectance and transmittance information, the values for some land covers were provided by Dr. Mariano García (University of Alcalá) and Dr. Olga Rosero (University of Zaragoza) (personal communication) (see Table 3). In the case of vegetation, reflectance was measured in laboratory conditions with a leaf clamp, while transmittance values were estimated using PROPECT and LIBERTY, adjusting them to the curve of the measured reflectance. Later, the values of these large groups were assigned to the different species located in the study area, taking as reference the reflectance measured in the field. In addition, in all scenes it is required to include the reflectance and transmittance of the terrain or bare soil. In this case, the reflectance value provided was measured with a contact probe [41].  Table 4 presents the ALS sensor parameters entered in the DART model, adjusted to the real capture. As some of the parameters were not provided by the sensor developer (e.g., area of LiDAR-PNOA sensor), they were suggested by Dr. Tiangang Yin, developer of DART, based on his expertise.

Processing of ALS-PNOA-Data and Simulated DART Point Clouds
The first processing step of ALS-PNOA data was noise removal. Then, ground points were classified using the multiscale curvature classification algorithm [42], implemented in the MCC 2.1 command-line tool, according to Montealegre et al. [43]. The Point-TIN-Raster interpolation method [44], implemented in ArcGIS 10.5 software (ESRI, Redlands, CA, USA), was applied to the ground points to produce a digital elevation model (DEM) with a 1-m grid size, following Montealegre et al. [45]. The ground elevation value of the DEM was subtracted from the ALS point height in order to obtain the normalized heights using FUSION LDV 3.30 open source software [46].
DART point clouds are stored in *.txt file format. The simulated point clouds were converted into LAS format using the "txt2las" tool available in LAStools software (https: //rapidlasso.com/lastools/). Normalization of the data was not required since the heights were referred to the ground level. Both point clouds, simulated by DART and acquired by PNOA, were clipped to fit the 15 m radius of the field plots.
Additionally, three structural diversity indices (DI) were computed. The Foliage Height Diversity Index or also called LiDAR height diversity index (LHDI) [47], Equation (1), which is an adaptation of the Shannon (H ) index, the LiDAR height evenness index (LHEI) proposed by Listopad et al. [47], Equation (2), that adapts the Pielou (J ) index and Rumple index [48] as a measure of roughness or structural heterogeneity, Equation (3).
where p is the proportion of returns at regular intervals of 0.5 m or at defined Prometheus classification height intervals (h). The first step to compute LHDI and LHEI was the calculation of return proportion at different height intervals using the "strata" switch within the "Cloud Metrics" command of FUSION. Thus, regular intervals of 0.5 m were selected according to Listopad et al. [47].
Rumple index is the ratio of three-dimensional canopy surface model (CSM) to ground area. Rumple was computed as the ratio between the sum of the three-dimensional area of triangles from CSM grid points to the two-dimensional area of the grid cell surface. The CSMs were created for each plot using a 1.5 m pixel and a 3 × 3 smoothing algorithm, considering point cloud density of ALS-PNOA. The surface area of each CSM 1.5 m pixel is computed by creating triangles that fit the centroid of the pixel and those of the neighboring ones. CSMs were created using the highest returns of each height range to account for canopy roughness. Rumple was computed for Prometheus classification height ranges (0.6, 2 and 4) to characterize heterogeneity within each stratum and for the overall forest canopy. The three diversity indexes were generated in R environment for both, simulated and ALS-PNOA data.

Accuracy Assessment of DART Simulations
The accuracy assessment of DART simulations was performed by comparing the simulated data with ALS-PNOA data in the 104 field plots. A series of statistics previously calculated, that commonly used in forestry that describes the canopy height (CHM), canopy height variability (CHVM), canopy density (CDM) and three diversity indexes were compared using the Spearman correlation coefficient. The Spearman correlation coefficient ranges from −1 to 1, values closest to 1 indicates highest positive correlation, values closest to −1 indicate highest negative correlations and values closest to 0 indicates null correlation [49].

Fuel Type Classification
As mentioned in Section 2.2.2. Field data, fuel type was assigned to the field plots using the Prometheus classification that is based on the type, height and coverage percentage of the fire propagation elements. These fuel types were assigned to the point clouds from both PNOA captures clipped to the field plot extension and the corresponding simulations.
The most explanatory LiDAR simulated metrics for fuel model discrimination, were selected using two selection methods according to Domingo et al. [50]: (i) Spearman rank correlation selection method; (ii) all subset selection, considering four different approaches: comprehensive, forward, backward and sequential replacement. Metric selection was performed independently for 2011 and 2016.
All subset selection determines the best variables of a group, without considering the rest of the variables [51]. Four searching techniques were tested: exhaustive, backward, forward and sequential replacement (seqrep). The maximum size of subsets was set to 6, while tests were performed between 4 to 6 subsets. Spearman's correlation and all subset selection were computed within R environment. R package "leaps" and specifically the "regsubsets" function was applied for all subset selection.
Fuel type classification was performed for 2011 and 2016 using the Support Vector Machine (SVM) artificial intelligence method [52] and including the most suitable LiDAR simulated metrics determined by the selection methods. SVM method allows multiclass classification assigning each class to the one with higher probability. In this sense, the "Cclassification" was selected using the "e1071" package in R environment. The classification was trained and parameterized using the simulated DART metrics, being validated with the metrics obtained from ALS-PNOA. An SVM is a supervised learning algorithm that allows pattern recognition and is based on the hypothesis that data are separable into classes in space, trying to find the optimal separation between classes through multidimensional hyperplanes. The data located in the hyperplanes are called support vectors, being these the most complex to classify, since there is less separability between classes. The SVM models with radial kernel and linear kernel were generated. The cost and gamma parameters were parameterized using the intervals 1-1000 and 0.01-1, respectively, in accordance with Domingo et al. [50]. The classification overall accuracy, confusion matrices, user's and producer's accuracy for the fitting and validation phases were evaluated to compare and, subsequently, determine the best classification model [53].  Table 6 (see Table 1 in Appendix A for all correlation values). The highest correlations are associated to high height percentiles and to CHVM metrics such as standard deviation and variance for both 2011 and 2016. Furthermore, several metrics from CDM and DI reach values over 0.90, as for example mean above_4 or Rumple. CHM related with lower heights, as for example low percentiles or minimum height, present lower correlation than CHM high height metrics for both 2011 and 2016 ALS-PNOA captures. A similar trend is observed for CDM. Metrics related to lower strata show lower correlation than those from higher strata. Diversity indices (DI) show high correlation values, while LHEI present lower values close to 0.75 for both years.   Table 6 shows the metrics with correlation coefficients between the forest fuel models and simulated point cloud metrics for both years, 2011 and 2016, that were subsequently selected for fuel model classification. The metrics show slight differences between years 2011 and 2016 and similar trends respect to the most suitable metrics. In this sense, high percentiles, CHVM as for example variance and Rumple diversity index show high correlations. CDM associated to low strata and high strata show higher correlation than intermediate height strata.  The simulated LiDAR metrics selected by all subsect selection approaches are presented in Table 7. Considering the maximum number of selected metrics, the selection methods included generally one or two metrics related to CHM, one metric associated with CHVM and one to two metrics that express the canopy density in 2011. The selection approaches did not include CHVM for the 2016 year, selecting only CHM and CD metrics. Furthermore, DI were not selected for any of the years and approaches.  Table 8 shows the overall accuracy of the two best classification models using the most explanatory metrics for fuel type classification that were determined by the Spearman's rank correlation. The model for both years, 2011 and 2016, included five simulated metrics: the 80 th percentile of return heights (P80), the coefficient of variation of the L moments (Elev. L. CV), the mean height of 0 to 0.6 strata returns (Mean 0_0.60), Rumple and LHDI. The best classification method was SVM with radial kernel, with an overall classification accuracy after validation of 0.88 for 2011 and 0.91 for 2016. The models were tuned with a cost value of 50 and a gamma value of 0.15. Lower accuracies were obtained with SVM with linear kernel that shows a decrease in overall accuracy of 0.19 for both 2011 and 2016. Table 8. Comparison between support vector machine with radial kernel (SVMr) and support vector machine with linear kernel (SVMl) classification method using overall accuracy for the selected metrics based on Spearman correlation. OA stands for overall accuracy.

Metrics
Year  Table 9 shows the performance of models computed using all subsect selection metrics and SVM with radial kernel (see Table 2 Appendix A for results using SVM with linear kernel). The selected metrics have overall accuracy values ranging from 0.72, when using the sequential replacement and exhaustive approaches, up to 0.74, when using the forward or backward selected metrics for 2011. Backward provide higher overall accuracy (0.84) than other all subsects selection approaches for 2016 year. However, none of the last mentioned models, computed using all subset selection metrics, improved the performance obtained using Spearman rank correlation coefficients as a selection method, being disregarded for subsequent fuel model classification. According to the results shown in Table 10 for year 2011, derived from the validation sample, there is confusion between types 1 and 2 as well as between types 3 and 2, which could be related with terrain complexity that may blurs height differentiation. Table 11 shows the results, derived from the validation sample, for year 2016. Misclassification is found between types 1 and 2 as well as between types 3 and 2, as was previously reported for 2011. Furthermore, confusion is found between types 5 and 6 and slight confusion between types 6 and 7. The mean user's classification accuracy was 88.5 % for 2011 and 93.8% for 2016, respectively. The mean producer's accuracy ranged from 89.3% in 2011 to 91.3% in 2016.

Discussion
The simulation of discrete-return low density ALS data provides essential information to support forest management at regional scales due to the widespread use of ALS data in operational forestry [25]. Field surveys have been traditionally carried out to derive ground truth data for training modes but these are challenging and expensive [13]. The simulation of 3D data using RTM could help reducing fieldwork and increasing training samples size. This is even more relevant when performing studies at regional scale and using nationwide ALS coverages of low point densities. The goodness of the results shows the potential of simulation with radiative transfer models in order to exploit more effectively an existing resource, such as ALS-PNOA data, in our case by fuel type classification. Furthermore, fuel type classifications are relevant to support preventive actions, manage forest fires and assist on fire modelling [54]. These classifications are especially relevant in forested areas recurrently affected by wildfires as the case of Mediterranean basin [19].
The simulation in DART requires a great amount of information, related to the sensor characteristics, the vegetation in the area to be simulated, its reflectance and transmittance values and so forth. In this sense, the program will make better simulations the more accurate the initial information entered. In this work, reflectance, transmittance and LAI values had to be estimated in absence of field or laboratory information. Concerning the reflectance and transmittance values, Lamelas et al. [24] used the vegetation and soil databases available in DART to include the optical properties of turbid medium instead of using specific values for the species present in the study area, as these authors do not compare or validate their simulations with real LiDAR captures. Roberts et al. [25] used reflectance values measured in the field with an spectrometer. We preferred to select laboratory values of species located in Mediterranean environments, adjusted to the species located in the study area and provided by experts, as the DART model is developed to work with laboratory values instead of field data. With respect to the LAI values Lamelas et al. used these parameters to simulate the fuel load of the different fuel types simulated, but, as mentioned before, they did not try to simulate real plots. Roberts et al. [25] imported 3D trees generated in external modelling software to simulate vegetation, directly creating DART voxels comprised of triangles and parallelograms with discrete reflectance and transmittance characteristics that do not require assigning a LAI value. Our objective was to analyze the accuracy of the DART model to replicate the low density small-footprint ALS data captured in the PNOA project and to assess the ability of simulations for model training to classify the fuel types of the study are, accordingly, in absence of field data, the LAI value to be assigned to the turbid medium was estimated applying the Biophysical Processor tool integrated in the SNAP software to two Sentinel 2-A top of canopy normalized reflectance data scenes (Level-2A images).
Our results demonstrate that DART simulations could replicate low-density discretereturn ALS data across 104 field plots and derive commonly used 3D metrics. The simulated metrics show an average correlation value of 0.55 in 2011 and 0.50 in 2016 with respect to ALS-PNOA. A correlation value of 0.8 was overpassed by 32% of the simulated metrics commonly used for forestry modelling in 2011 and 28% in 2016, respectively. High height percentiles, canopy variability metrics and Rumple diversity index reached correlation values over 0.94. Lower correlation values were found for low height percentiles and low strata density metrics. The high correlation values, especially for higher or maximum CHM, are in accordance with Robert et al. [25], that observed absolute errors around 0.5 m when replicating ALS measurements with DART.
The relationships between fuel types and ALS metrics was assessed by two selection methods: (i) Spearman rank correlation coefficients [55], showing good explanatory power as the most suitable metrics for modelling were selected by this method in accordance with Domingo et al. [19]; and by (ii) all subsets selection approach, that provided some variables as the mean height of 0 to 0.6 strata returns. The variables selected by the most accurate model are the 80th percentile of return heights (P80), the coefficient of variation of the L moments (Elev. L. CV), the mean height of 0 to 0.6 strata returns (Mean 0_0.60), Rumple and LHDI. Our results agree with Valvuena et al. [56], who concluded that the L-moments from the distribution of ALS returns can have a direct relationship to forest structural characteristics at the community level. The use of structural complexity indexes, such as Rumple and LHDI, in fuel classification is considered to generate parsimonious models in accordance with Domingo et al. [19] and Gelabert et al. [57], reducing the number of metrics used when applying the height bin approach [18,57]. The selection of a high percentile (P80) instead of a lower percentile selected in previous studies [19] may have been caused by the low correlation between simulation metrics and real ones in the low percentile's metrics. In this sense, although our results are satisfactory for our final objective, fuel types mapping, in future research it may be convenient to test the suitability of using synthetic trees with explicitly defined crown architectures [22] or specific laboratory reflectance and transmittance information.
The comparison between kernel approaches for SVM showed that the highest accuracy to classify Prometheus fuel types were produced using radial kernel as was previously reported by García et al. [14] and Domingo et al. [19]. The performance of the classification with an overall accuracy value of 88% for 2011 and 91% for 2016 shows similar results to the ones obtained by García et al. [14] (88% overall agreement) and Alonso-Benito et al. [58] (85% overall agreement) that classified Prometheus fuel types using LiDAR and high-resolution images. Higher classification performance than Huesca et al. [13] and Domingo et al. [19] was found, who classified Prometheus fuel types using low-density ALS-PNOA data from the first and second coverage, respectively. Their results highlight that, though the fusion of ALS with multispectral data and the use of higher point density increase fuel classification performance, the high structural complexity of vegetation in Mediterranean environments constitute a handicap in classification performance as reported in our previous work [19]. In this sense, the use of simulations can increase the number and variability in the sample to train the models, improving the final precision and reducing costs.
Furthermore, although overall accuracy is high there are some confusion in classification between fuel types 5, 6 and 7 that was also previously reported by García et al. [14] and Huesca et al. [13]. Similarly, there exits confusion between types 1 to 3, also reported by Huesca et al. [13], that could be related in our study area with steep slopes that generate a decrease in digital elevation model accuracy [42] and consequently blurs height differentiation.
In summary, the present study shows the utility of DART simulations to generate accurate 3D data for fuel type classification. This encourages to conduct more research to predict different forestry metrics in Mediterranean forests through simulation. Furthermore, it should be considered to analyze other methods and techniques to assess DART simulations accuracy as the ones proposed by Roberts et al. [25]. In addition, extra analysis might focus on increasing the sample size in predicting fuel types using simulated plots resulting from expert knowledge instead of plots previously measured in field. Finally, it would be interesting to analyze the sensitivity of the simulations to the use of field and laboratory information related to LAI, reflectance and transmittance values.

Conclusions
This study assessed the usefulness of DART model to simulate low density smallfootprint ALS data and subsequently train and classify Prometheus fuel types in Mediterranean environment. The high correlation of metrics such us High height percentiles, canopy variability metrics and Rumple diversity demonstrates that DART simulations could replicate low-density discrete-return ALS data. Spearman rank coefficient was the most powerful selection method to generate a representative and meaningful fuel types classification. The SVM with radial kernel method produced the most accurate fuel type classification model, which included five ALS metrics: the 80th percentile of return heights, the coefficient of variation of the L moments, the mean height of 0 to 0.6 strata returns, Rumple and LHDI diversity indexes. The classification was trained using DART simulated data and validated with ALS-PNOA data from 2011 and 2016 coverages. The best-fitted classification achieved an overall accuracy of 88% for 2011 and 91% for 2016 years, respectively. The results revealed that DART simulations provide suitable 3D data that can be used for model training, implying an important decrease of human and economic resources invested in field surveys conducted for fuel mapping and open a promising line of research to improve simulation by analyzing the difficulties encountered in this study.