Attribution of Disturbance Agents to Forest Change Using a Landsat Time Series in Tropical Seasonal Forests in the Bago Mountains, Myanmar

In 2016, in response to forest loss, the Myanmar government banned logging operations for 1 year throughout the entire country and for 10 years in specific regions. However, it is unclear whether this measure will effectively reduce forest loss, because disturbance agents other than logging may have substantial effects on forest loss. In this study, we investigated an approach to attribute disturbance agents to forest loss, and we characterized the attribution of disturbance agents, as well as the areas affected by these agents, in tropical seasonal forests in the Bago Mountains, Myanmar. A trajectory-based analysis using a Landsat time series was performed to detect change pixels. After the aggregation process that grouped adjacent change pixels in the same year as patches, a change attribution was implemented using the spectral, geometric, and topographic information of each patch via random forest modeling. The attributed agents of change include “logging”, “plantation”, “shifting cultivation”, “urban expansion”, “water invasion”, “recovery”, “other change”, and “no change”. The overall accuracy of the attribution model at the patch and area levels was 84.7% and 96.0%, respectively. The estimated disturbance area from the attribution model accounted for 10.0% of the study area. The largest disturbance agent was found to be logging (59.8%), followed by water invasion (14.6%). This approach quantifies disturbance agents at both spatial and temporal scales in tropical seasonal forests, where limited information is available for forest management, thereby providing crucial information for assessing forest conditions in such environments.


Introduction
Tropical forests make important contributions to the global carbon cycle, biodiversity, and ecosystem processes because they are carbon-and species-rich ecosystems [1][2][3]. However, deforestation and forest degradation in tropical forests have been of major concern in the past few decades [4][5][6][7][8]. Especially, tropical forests in Southeast Asia face one of the highest rates of deforestation and forest degradation among tropical regions [9][10][11]. Thus, consistent and periodic monitoring of forest conditions is crucial for sustainable forest management.
Myanmar is one of the most forested countries on the mainland of Southeast Asia. However, the forests in Myanmar are being deforested and degraded rapidly. For example, according to the country report of the Global Forest Resources Assessment 2015 [12], the forest cover of Myanmar decreased from 51.5% to 42.9% between 2000 and 2015. The annual forest loss between 2010 and 2015 reached 546,000 ha, which was the third largest forest loss in the world [13]. In 2016, in response to the forest loss, the Myanmar government banned logging activities for one year throughout the entire country and for 10 years in specific regions. However, it is unclear whether the ban will effectively reduce forest loss, as previous studies have shown that disturbance agents other than logging, such as dam construction and shifting cultivation, play an important role in deforestation in Myanmar [14][15][16][17]. Thus, the attribution of disturbance agents and the quantification of areas affected by disturbance agents must be characterized in a consistent manner to obtain accurate information. However, to our knowledge, no study has attributed the areas affected by different disturbance agents in Myanmar.
While field measurements provide highly accurate information, they are costly, labor intensive, and limited in spatial scale. Thus, the information from field measurements is limited in the spatial and temporal dimensions. Satellite images, because of their immense spatial dimension and high temporal frequency, may reliably detect forest changes over large areas and provide insights into the agents underlying these changes [18][19][20]. Landsat data are especially well suited to detecting forest changes over decades and a variety of spatial extents [21,22]. Free access to the Landsat archives provides an opportunity to develop new methodologies that use dense Landsat time series for detecting forest disturbance and/or recovery [23][24][25][26][27][28][29][30]. Dense Landsat time series have preferable characteristics to capture accurate forest change information [29][30][31][32][33]. Furthermore, analyzing the spectral trajectory of dense Landsat time series enables major and minor forest changes to be detected reliably [24,[34][35][36].
While many studies have quantified pre-and post-disturbance events using dense Landsat time series [37][38][39][40], detecting disturbance agents using dense Landsat time series is a relatively new approach. Kennedy et al. [41] tested a method for forest change agent attribution using annual Landsat images in the Puget Sound region, USA. They showed that it was reasonably accurate for a variety of agents across different land cover types (84% overall accuracy). Similarly, Hermosilla et al. [42] performed a spectral trend analysis to attribute various forest change types over a vast area in Canada. They achieved an overall accuracy of 91.6% for object-level change attributions. These two studies focused on boreal or temperate forests; however, to date, there is a paucity of such studies in natural tropical forests. Because disturbance agents vary depending on watersheds [41], countries [43], or broad regions [44,45], the developed methods may not be applicable to other landscapes [43]. Here, we evaluated the applicability of a time series of annual Landsat images for attributing forest disturbances in a tropical seasonal forest.
Our main objectives were (1) to test and demonstrate the applicability of a trajectory-based approach using a time series of annual Landsat images to attribute a wide variety of disturbance agents and (2) to characterize the areas affected by these agents in tropical seasonal forests in the Bago Mountains, Myanmar. Assessing agents of change and their impacts will provide comprehensive understanding of forest change in the study area, which will benefit forest management.

Study Area
The study area is located in the Bago Mountains in central Myanmar ( Figure 1). The area covers approximately 2.31 million ha, which accounts for 3.4% of the total land area of Myanmar. The area is dominated by mixed deciduous forests, which comprise the most economically productive forest type in Myanmar. Thus, the study area is a major timber-extraction area in Myanmar, and selective logging has been conducted to produce timber products under the Myanmar Selection System since 1856 [15]. Under the Myanmar Selection System, economically important hardwood species larger than a predefined minimum stem size are selectively harvested based on a few rules, such as a 30-year felling cycle and an annual allowable cut for the timber volume based on forest inventories.
In addition to the selective logging activities (i.e., legal logging activity), several other factors, such as illegal logging, plantations, shifting cultivation, dam construction, agricultural expansion, and urban development, cause forest disturbances in the study area. While there is no detailed estimate of the illegal logging intensity in this region, several studies implied that illegal logging activities caused severe forest losses [14,15]. Shifting cultivation has also been conducted by indigenous people since the 19th century [16]. 1856 [15]. Under the Myanmar Selection System, economically important hardwood species larger  than a predefined minimum stem size are selectively harvested based on a few rules, such as a 30year felling cycle and an annual allowable cut for the timber volume based on forest inventories. In addition to the selective logging activities (i.e., legal logging activity), several other factors, such as illegal logging, plantations, shifting cultivation, dam construction, agricultural expansion, and urban development, cause forest disturbances in the study area. While there is no detailed estimate of the illegal logging intensity in this region, several studies implied that illegal logging activities caused severe forest losses [14,15]. Shifting cultivation has also been conducted by indigenous people since the 19th century [16].  [46] was used to define the country borders.

Remote Sensing Data
We selected a total of 126 Level 1 Terrain corrected Landsat TM, ETM+, and OLI images located on five paths/rows of the World Reference System 2 (i.e., 132/47, 132/48, 133/46, 133/47, and 133/48) from 2000 to 2014 (Table S1). Only images acquired from November to February were used to avoid leaf-off conditions. For each year of the period, one cloud-free and Scan Line Corrector (SLC)-on image was selected to construct a Landsat time series. In instances where there were no such images, up to three SLC-off images or partially cloud-contaminated images were selected to create gap-free images using a gap-filling procedure.

Topographic Data
A digital elevation model (DEM) was used to model the terrain in the study area. We used 30-m resolution data from the Shuttle Radar Topographic Mission (SRTM) obtained from the United States Geological Survey archives [47]. The original DEM was resampled and co-registered with the Landsat imagery of the study area. Elevation, slope, and aspect were extracted from the DEM and used to perform a topographic correction. These topographic variables were also used in the construction of the change attribution model.  [46] was used to define the country borders.

Remote Sensing Data
We selected a total of 126 Level 1 Terrain corrected Landsat TM, ETM+, and OLI images located on five paths/rows of the World Reference System 2 (i.e., 132/47, 132/48, 133/46, 133/47, and 133/48) from 2000 to 2014 (Table S1). Only images acquired from November to February were used to avoid leaf-off conditions. For each year of the period, one cloud-free and Scan Line Corrector (SLC)-on image was selected to construct a Landsat time series. In instances where there were no such images, up to three SLC-off images or partially cloud-contaminated images were selected to create gap-free images using a gap-filling procedure.

Topographic Data
A digital elevation model (DEM) was used to model the terrain in the study area. We used 30-m resolution data from the Shuttle Radar Topographic Mission (SRTM) obtained from the United States Geological Survey archives [47]. The original DEM was resampled and co-registered with the Landsat imagery of the study area. Elevation, slope, and aspect were extracted from the DEM and used to perform a topographic correction. These topographic variables were also used in the construction of the change attribution model.

Methods
The change detection, characterization, and attribution processes are summarized in the flowchart ( Figure 2). The processing workflow includes: (1) image pre-processing, which includes atmospheric correction, normalization, topographic correction, cloud masking, and gap filling; (2) forest change detection using the LandTrendr algorithm [24]; followed by (3) spatial aggregation of change pixels; and (4) change attribution to disturbance agents and recovery conditions using a random forest (RF) model [48].

Methods
The change detection, characterization, and attribution processes are summarized in the flowchart ( Figure 2). The processing workflow includes: (1) image pre-processing, which includes atmospheric correction, normalization, topographic correction, cloud masking, and gap filling; (2) forest change detection using the LandTrendr algorithm [24]; followed by (3) spatial aggregation of change pixels; and (4) change attribution to disturbance agents and recovery conditions using a random forest (RF) model [48].

Pre-Processing of the Landsat Time Series
Pre-processing of the Landsat time series consisted of atmospheric correction, normalization, topographic correction, cloud masking, and gap filling. A cloud-free image was selected as a reference image for each path/row, and it was atmospherically corrected to the surface reflectance image using the Landsat Ecosystem Disturbance Adaptive Processing System, which uses a 6S radiative transfer model to correct for atmospheric effects on a given date [49]. All other Landsat images were normalized to the reference images using the iteratively reweighted multivariate alteration detection (MAD) algorithm [50]. Then, the C-correction [51] was applied to all the images to remove topographic effects in the Landsat time series. Clouds and cloud shadows were automatically masked using the fmask algorithm with the default parameters [52]. Subsequently, data gaps due to SLC-off and clouds were filled with a weighted linear regression procedure [53] and the regression method proposed by [54] using auxiliary Landsat images in the same year, respectively.

Forest Change Detection
LandTrendr, which is a temporal segmentation and fitting algorithm, was used to conduct a trajectory-based change detection. The algorithm investigates the disturbance and growth history by fitting straight line segments capturing the overall shape of the trajectory to the spectral values of each pixel. A comprehensive explanation of the LandTrendr algorithm is provided in Kennedy et al. [24]. The LandTrendr algorithm is applied to a single spectral index or variable. In this study, Tasseled Cap Wetness (TCW), which is derived from a Tasseled Cap Transformation [55], was used as the spectral index because of its sensitivity to disturbance events [18,24] and water contents [56]

Pre-Processing of the Landsat Time Series
Pre-processing of the Landsat time series consisted of atmospheric correction, normalization, topographic correction, cloud masking, and gap filling. A cloud-free image was selected as a reference image for each path/row, and it was atmospherically corrected to the surface reflectance image using the Landsat Ecosystem Disturbance Adaptive Processing System, which uses a 6S radiative transfer model to correct for atmospheric effects on a given date [49]. All other Landsat images were normalized to the reference images using the iteratively reweighted multivariate alteration detection (MAD) algorithm [50]. Then, the C-correction [51] was applied to all the images to remove topographic effects in the Landsat time series. Clouds and cloud shadows were automatically masked using the fmask algorithm with the default parameters [52]. Subsequently, data gaps due to SLC-off and clouds were filled with a weighted linear regression procedure [53] and the regression method proposed by [54] using auxiliary Landsat images in the same year, respectively.

Forest Change Detection
LandTrendr, which is a temporal segmentation and fitting algorithm, was used to conduct a trajectory-based change detection. The algorithm investigates the disturbance and growth history by fitting straight line segments capturing the overall shape of the trajectory to the spectral values of each pixel. A comprehensive explanation of the LandTrendr algorithm is provided in Kennedy et al. [24]. The LandTrendr algorithm is applied to a single spectral index or variable. In this study, Tasseled Cap Wetness (TCW), which is derived from a Tasseled Cap Transformation [55], was used as the spectral index because of its sensitivity to disturbance events [18,24] and water contents [56] (see next paragraph). To implement the LandTrendr algorithm, a set of parameters and threshold values was used to divide the spectral trajectory into several segments and to label each segment as either disturbance or no-disturbance. For these parameters, we used the parameters tuned by our prior studies for tropical forests in Myanmar [57]. Threshold values of 0.05 and −0.05 were set to define forest change. Thus, only segments with TCW changes ≥0.05 or ≤−0.05 were assigned as forest change. The accuracy of forest change detection was investigated using 379 pixel-based samples that were distributed randomly over the study area. We determined the threshold value based on a study conducted in the same study area [58]. This study determined the relationship between the threshold value in LandTrendr and the producer's and user's accuracies of change class. The study indicated that the producer's accuracy for the change class was close to 100% until the threshold value reached 0.05. Then, the producer's accuracy decreased with increasing threshold values greater than 0.05. We determined the threshold with a focus on the producer's accuracy of the change class in order to reduce the misdetection of disturbances. As a result, the detected change areas might include a certain amount of false positive. Thus, in the change attribution process, we defined a "no change" class, which includes areas with change detection errors (see Section 3.4).
In this study, we assigned both positive and negative spectral changes as forest change (i.e., 0.05 and −0.05, respectively), while previous studies detected only positive or negative spectral changes as forest disturbance. This is because the TCW trajectory caused by water invasion has the opposite temporal evolution compared with the trajectories of other agents ( Figure 3). Higher TCW values are often found over forests, and lower TCW values are more likely in open stands or clearcut areas, compared with dense forests. Thus, generally, TCW decreases when there is a disturbance event that removes vegetation (Figure 3). In contrast, TCW increases when there is a disturbance caused by water invasion, as TCW is the indicator of water contents and its value is high in bodies of water [56]. (see next paragraph). To implement the LandTrendr algorithm, a set of parameters and threshold values was used to divide the spectral trajectory into several segments and to label each segment as either disturbance or no-disturbance. For these parameters, we used the parameters tuned by our prior studies for tropical forests in Myanmar [57]. Threshold values of 0.05 and −0.05 were set to define forest change. Thus, only segments with TCW changes ≥0.05 or ≤−0.05 were assigned as forest change. The accuracy of forest change detection was investigated using 379 pixel-based samples that were distributed randomly over the study area. We determined the threshold value based on a study conducted in the same study area [58]. This study determined the relationship between the threshold value in LandTrendr and the producer's and user's accuracies of change class. The study indicated that the producer's accuracy for the change class was close to 100% until the threshold value reached 0.05. Then, the producer's accuracy decreased with increasing threshold values greater than 0.05. We determined the threshold with a focus on the producer's accuracy of the change class in order to reduce the misdetection of disturbances. As a result, the detected change areas might include a certain amount of false positive. Thus, in the change attribution process, we defined a "no change" class, which includes areas with change detection errors (see Section 3.4).
In this study, we assigned both positive and negative spectral changes as forest change (i.e., 0.05 and −0.05, respectively), while previous studies detected only positive or negative spectral changes as forest disturbance. This is because the TCW trajectory caused by water invasion has the opposite temporal evolution compared with the trajectories of other agents ( Figure 3). Higher TCW values are often found over forests, and lower TCW values are more likely in open stands or clearcut areas, compared with dense forests. Thus, generally, TCW decreases when there is a disturbance event that removes vegetation (Figure 3). In contrast, TCW increases when there is a disturbance caused by water invasion, as TCW is the indicator of water contents and its value is high in bodies of water [56].

Spatial Aggregation of Change Pixels
Anthropogenic forest disturbances rarely occur at the single-pixel level, but rather over larger geographic areas (patches) [41]. In this study, we grouped adjacent, detected change pixels in the

Spatial Aggregation of Change Pixels
Anthropogenic forest disturbances rarely occur at the single-pixel level, but rather over larger geographic areas (patches) [41]. In this study, we grouped adjacent, detected change pixels in the same year as the patches, and we assumed that the detected, adjacent change pixels, which were categorized as a disturbance or recovery, in the same year were caused by the same disturbance event. To reduce the inclusion of erroneous detections, only disturbance patches greater than three pixels were retained. From this point onward, all disturbance events were assessed at the patch level, which maintained the change characteristics produced by LandTrendr at each pixel to calculate the patch-based metrics in the following step. Land cover conditions after changes were also considered using a similar patch-based approach.

Variable Extraction from the Patches
For each patch, we extracted trajectory-based variables to characterize disturbance agents. From the Landsat time series, we extracted variables, such as the duration of change events, the magnitude of change events, and pre-and post-events conditions. Duration was defined as the time over which the change occurred. Magnitude was calculated as the difference of the average spectral values of the patch before and after the change event. Pre-and post-event conditions were calculated as the average and standard deviation of the spectral values of the patch before and after the disturbance, respectively. The magnitude of the change events and the pre-and post-events conditions were computed from the Tasseled Cap Brightness (TCB), Tasseled Cap Greenness (TCG), and TCW derived from the Tasseled Cap Transformation [55] and Tasseled Cap Angle (TCA), which is based on the angle formed by the TCB and TCG [59]. Topographic conditions were characterized using the average elevation and average slope. The shape of each patch was also considered by calculating the area, perimeter, and compactness [60]. In total, we extracted 29 variables from each patch (Table 1).

Change Attribution
A change attribution was implemented by constructing a RF model with reference samples and applying the RF model to all the detected change patches. We randomly selected 1411 sample patches, and we manually interpreted the disturbance agents and recovery condition using high-spatial-resolution images from Google Earth™. The disturbance agents were defined as "logging", "plantation", "shifting cultivation", "urban expansion", and "water invasion", as shown in Table 2.
Regarding logging, we included areas affected by logging activity that did not show any transformation into other land use types. Because it was difficult to distinguish disturbances caused by selective logging activities (legal logging) from other logging activities, such as illegal logging, only using Google Earth, we did not separately attribute disturbances caused by selective logging from other logging disturbances. The attribution also considered "recovery", "other change", and "no change" patches. Recovery is attributed to areas showing vegetation recovery after a disturbance event. Because this study captured both positive and negative spectral change values as forest change, the changes in this study include not only forest disturbance, but also forest recovery. The other change areas include changes occurring outside forested areas (e.g., agricultural land). Because LandTrendr captures not only forest cover change, but also temporary variations such as forest health, we assigned such negative temporary variations as no change, along with areas that were not affected by any of the disturbance events considered in this study. The no change class also included areas with change detection errors. Table 2. Definition of disturbance agents in this study.

Agent Class Description
Logging Both short-and long-term stand replacing forest disturbances caused by human intervention such as logging, without transformation to other land use types.

Plantation
Short-term stand replacing forest disturbances followed by plantation development.
Shifting cultivation Short-term disturbances caused by shifting cultivation.

Urban expansion
Both short-and long-term stand replacing forest disturbances caused by city and agricultural land expansion, and road construction.

Water invasion
Water body invasion to forest caused by dam construction.

Recovery
Vegetation recovery occurred after disturbances.
Other change Changes in agricultural land, shrub and water body.

No change
No forest disturbances nor other changes.
We used 25% of the reference samples as test data, while 75% of the samples were used to construct RF models for agent attribution and recovery condition. The constructed RF models were evaluated using the separately kept test data, and confusion matrices and Kappa coefficients were produced to evaluate the accuracy at the patch and area levels. Then, the RF models were applied to all the disturbance patch maps in each year. Subsequently, the temporal trends of the disturbance area were investigated using a Mann-Kendall trend test [61,62]. The attribution accuracy was measured using two confusion matrices for the patch and area levels.
To run the RF model, we used the "random forest" package [63] in version 3.24 of R statistical package [64]. The R package "caret" [65] was also used for tuning the RF models with a 10-fold cross-validation.

Accuracy Assessments for Pixel-Based Change Detection
The pixel-based accuracy assessment for change detection showed that the overall accuracy was 71.8%. The producer's and user's accuracies for change were 91.9% and 35.8%, respectively (Table 3), while the producer's and user's accuracies for no change were 67.8% and 97.8%, respectively.

Accuracy Assessments for Patch-Based Change Attribution
The RF model for disturbance attribution achieved an overall accuracy of 84.7% and a Kappa coefficient of 0.81 at the patch level (Table 4). Water invasion showed the lowest producer's accuracy (73.3%) and the highest accuracy was achieved by urban expansion (93.3%). Shifting cultivation and urban expansion achieved the highest user's accuracy (100%), while no change showed the lowest user's accuracy (77.0%). The overall accuracy and Kappa coefficient at the area level were 96.0% and 0.95, respectively (Table 5). In general, compared with the patch-level accuracy, both the user's and producer's accuracies for logging, water invasion, and no change were higher at the area level. The producer's accuracy for plantation was higher at the area level, while its user's accuracy was lower. Figure S1 represents random forest variable importance in terms of the mean Gini decrease, which is the sum of the decrease in Gini impurity when a given variable is used, then averaged over the number of the trees in the forest [66,67]. The TCW magnitude (Mag_TCW) was the most important variable (83.8), followed by the TCA magnitude (Mag_TCA) (66.2) and the pre-disturbance mean TCA (Pre_mean_TCA) (60.7).

Forest Changes in the Bago Mountains
The disturbance area in the Bago Mountains from 2001 to 2013 was estimated by applying the RF model to attribute disturbance agents to all the patches in each year (Figure 4). The Mann-Kendall trend test showed that there was no significant increasing trend for the disturbance area (p > 0.05, Figure 4). In general, 10.0% of the study area experienced a forest disturbance in the past 15 years. Logging accounted for 59.8% of the total disturbance area, while plantation development, shifting cultivation, and urban expansion accounted for 8.4%, 10.4%, and 6.8%, respectively, of the disturbance area. Water invasion accounted for only 14.6% of the disturbance area. The plantation development and urban expansion areas increased significantly in recent years (p < 0.05), while shifting cultivation decreased significantly (p < 0.05, Figure 5). No significant trends were found for logging and water invasion. The spatial distributions of the disturbance agents were mapped by aggregating 26 townships in the study area ( Figure 6). Logging accounted for the largest proportions of disturbance in 23 townships. Plantation accounted for significant proportions of disturbance in the Bago (Pegu) and Hlegu townships, which are in the southern part of the study area, while urban expansion was highly distributed in Pyinmana, which is the northeastern part of the study area. Shifting cultivation and water invasion were found mainly in the eastern (the Yedashe, Thoungoo, and Pyu townships) and southeastern parts (the Pyu, Kyauktaga, and Bago (Pegu) townships), respectively, of the study area.       Forest recovery was observed in 9.1% of the study area, and it did not show any significant trend (Figure 7). In 2001, the compositions of former disturbances were not quantified, because the disturbance events occurred before 2001, which was outside our study period. The largest recoveries were recorded in logging areas (46.6%), followed by shifting cultivation (11.2%), plantation (4.2%), urbanization (0.6%), and water invasion (0.1%). These disturbance agents accounted for 62.7% of the total recovery. Forest recovery was observed in 9.1% of the study area, and it did not show any significant trend (Figure 7). In 2001, the compositions of former disturbances were not quantified, because the disturbance events occurred before 2001, which was outside our study period. The largest recoveries were recorded in logging areas (46.6%), followed by shifting cultivation (11.2%), plantation (4.2%), urbanization (0.6%), and water invasion (0.1%). These disturbance agents accounted for 62.7% of the total recovery.

Discussion
The RF models used for the attribution of disturbance agents produced substantial overall accuracies (84.7% at the patch level and 96.0% at the area level). The overall accuracy of the disturbance attribution is comparable to the results obtained in other environments by previous studies (84.0% for Kennedy et al. [41] and 91.4% for Hermosilla et al. [42]). Our study confirmed that the RF models achieved accurate attributions of the changes in tropical forests. The importance of the predictor variables showed that the trajectory-based variables, magnitude, pre-mean, and post-mean were the four most important variables. They were also selected from the spectral TCW, TCA, and TCG indices as the four most important variables. This indicates that the inclusion of disturbance information from various spectral indices is important for the attribution of disturbance agents. The results from this study indicate that the trajectory-based approach is suitable for providing accurate change attributions. It should be noted that, however, we used only TCW for forest change detection. Because using various spectral indices is important for the change attribution, the inclusion of various spectral indices may also increase the accuracy of forest change detection. Further study is required to develop ways to detect forest change using results of LandTrendr analyses on several indices simultaneously.
When the RF model was applied for the attribution of disturbance agents for the entire study area, 10% of the study area was shown to be affected by disturbance events during the past 15 years. At the national level, areas of deforestation and forest degradation are increasing in Myanmar [12]. While the study area accounted for only 3.4% of the total land area of Myanmar, we did not observe any significant increasing trend of deforestation and forest degradation, in contrast to the countrylevel trend. One possible reason for the lack of a trend could be that the study area has been a major timber-extraction area since 1856. Thus, the forests in this study area have been under high pressure, and, thus, the disturbance intensity in the study area has remained high during the past 15 years.
Most of the disturbances occurred because of a single disturbance agent (i.e., logging, which accounted for approximately 59.8% of the disturbed areas). In this study, because of the lack of spatially explicit reference samples for mapping detailed agents, the logging class includes both illegal and legal logging activities. Thus, we interpreted the proportion of the area affected by logging with some caution. Because previous studies showed that severe forest disturbances were not caused by illegal logging, but rather by legal logging activities [14,15], the results from our study may reveal

Discussion
The RF models used for the attribution of disturbance agents produced substantial overall accuracies (84.7% at the patch level and 96.0% at the area level). The overall accuracy of the disturbance attribution is comparable to the results obtained in other environments by previous studies (84.0% for Kennedy et al. [41] and 91.4% for Hermosilla et al. [42]). Our study confirmed that the RF models achieved accurate attributions of the changes in tropical forests. The importance of the predictor variables showed that the trajectory-based variables, magnitude, pre-mean, and post-mean were the four most important variables. They were also selected from the spectral TCW, TCA, and TCG indices as the four most important variables. This indicates that the inclusion of disturbance information from various spectral indices is important for the attribution of disturbance agents. The results from this study indicate that the trajectory-based approach is suitable for providing accurate change attributions. It should be noted that, however, we used only TCW for forest change detection. Because using various spectral indices is important for the change attribution, the inclusion of various spectral indices may also increase the accuracy of forest change detection. Further study is required to develop ways to detect forest change using results of LandTrendr analyses on several indices simultaneously. When the RF model was applied for the attribution of disturbance agents for the entire study area, 10% of the study area was shown to be affected by disturbance events during the past 15 years. At the national level, areas of deforestation and forest degradation are increasing in Myanmar [12]. While the study area accounted for only 3.4% of the total land area of Myanmar, we did not observe any significant increasing trend of deforestation and forest degradation, in contrast to the country-level trend. One possible reason for the lack of a trend could be that the study area has been a major timber-extraction area since 1856. Thus, the forests in this study area have been under high pressure, and, thus, the disturbance intensity in the study area has remained high during the past 15 years.
Most of the disturbances occurred because of a single disturbance agent (i.e., logging, which accounted for approximately 59.8% of the disturbed areas). In this study, because of the lack of spatially explicit reference samples for mapping detailed agents, the logging class includes both illegal and legal logging activities. Thus, we interpreted the proportion of the area affected by logging with some caution. Because previous studies showed that severe forest disturbances were not caused by illegal logging, but rather by legal logging activities [14,15], the results from our study may reveal the crucial state of illegal logging activity in Myanmar. However, this result may indicate that the influence of legal logging on the forest is more severe than previous studies expected. Further studies that distinctly quantify both illegal and legal logging disturbance areas are required for a more comprehensive assessment.
Disturbance agents other than logging (i.e., plantation, shifting cultivation, urban expansion, and water invasion) were concentrated in particular areas. For instance, disturbances caused by urban expansion were found predominantly in the northeastern parts of the study area ( Figure 6). These disturbance areas are located close to the new capital city of Myanmar, which was recently constructed in 2005. Disturbances caused by shifting cultivation were observed mainly in the eastern part of the study area. These results are also consistent with expectations, because shifting cultivation has been only permitted in specific regions [68]. The areas affected by plantation development increased recently ( Figure 5), and they were observed mainly in the southern part of the study area. Plantations in this region include both teak and rubber plantations, and especially rubber plantations have been expanding recently in Southeast Asia [69,70]. Thus, our study reflects the recent development of plantations in this region.
The disturbances caused by water invasion represented the second largest forest disturbance in the study area (14.6%). Such disturbances were observed in several townships in the eastern part of the study area. In this region, several dams were constructed during the past 15 years. Our results reflect the frequent and large coverage of the disturbances caused by the construction of these dams. No significant trend was found for the areas affected by water invasion, but large forest disturbances were caused by water invasion in 2001, 2011, and 2013 ( Figure 5). In these years, the total disturbance areas also showed a large increase (Figure 4), and the disturbances caused by water invasion made a substantial contribution to the total disturbance area. Thus, we conclude that the construction of dams has had a major impact on the forests in Myanmar. In 2013, Wang et al. [71] estimated that 48 hydropower projects were planned, under construction, or already operational in Myanmar. These dams may have large impacts on carbon emissions and biodiversity in Myanmar forests [71], as confirmed by our study. While previous studies suggested that logging has a substantial impact on forest disturbances in Myanmar [14,15], little attention has been paid to dam construction; thus, further attention should be paid to the impacts of dams.
We also found a large amount of forest recovery in the study area (9.1%), as disturbances caused by logging, plantation, and shifting cultivation were followed by forest recovery in the region (Figure 7). This result implies that most of the disturbances caused by logging were followed by recovery in this region. This result, however, does not consider the recovery of forest structure, such as canopy height, aboveground biomass, and species composition, after the disturbance events, as it only considers the spectral recovery of TCW. While TCW is highly correlated with forest structure, the recovery of TCW does not necessarily indicate the recovery of the forest structure to its previous condition. Thus, recovery assessments in this study should be interpreted cautiously.
The methodological framework described here showed reliable results for quantifying disturbance areas and disturbance agents, and it enabled us to provide both spatial and temporal characterizations of forest disturbances. The spectral information derived from the time series of annual Landsat images was used to extract a variety of change metrics, which effectively distinguish between the different states of forest change caused by different disturbance agents. By mapping forest disturbances in time and space, we can estimate forest conditions and forest changes from the local to national levels. Together with disturbance agents, forest change maps provide a greater understanding of past and future forest conditions, such as deforestation, forest degradation, and carbon emissions, because different disturbance agents cause different forest structures after disturbance events.
Future studies will include assessments of more detailed agents of forest disturbances, such as those caused exclusively by illegal logging. In 2016, the Myanmar government banned all logging activities, even those that were previously legal. Thus, by applying the methodology described in this study using the Landsat time series and including data from 2016, we can quantify the effect of illegal logging. Other remote-sensing products, such as imagery from unmanned aerial vehicles or light detection and ranging, may aid the collection of detailed information regarding disturbance agents. Combining Landsat time series with these remote-sensing products can provide more useful information for understanding forest conditions.

Conclusions
The pre-and post-disturbance information derived from forest change detection using a time series of annual Landsat images provides a comprehensive understanding of forest changes. Here, we presented the methodology for mapping disturbance areas, disturbance agents, and the following recovery conditions in tropical seasonal forests in Myanmar. The results of disturbance agent attribution using the RF model showed that patch-level and area-level models achieved overall accuracies of 84.7% and 96%, respectively. The estimated disturbance area from the disturbance attribution model accounted for 10% of the total study area. Disturbances caused by logging accounted for the largest area, followed by those caused by water invasion and shifting cultivation. The methodological framework described here can be used to map widely distributed disturbance agents at both the spatial and temporal scales, thus providing critical information for assessing forest conditions in tropical seasonal forests.
Supplementary Materials: The following are available online at www.mdpi.com/1999-4907/8/6/218/s1, Figure S1: The importance of variables for the random forest disturbance agent model, Table S1: Landsat images used in each World Reference System 2 path/row.