Refugee Camp Monitoring and Environmental Change Assessment of Kutupalong, Bangladesh, Based on Radar Imagery of Sentinel-1 and ALOS-2

Abstract

Approximately one million refugees of the Rohingya minority population in Myanmar crossed the border to Bangladesh on 25 August 2017, seeking shelter from systematic oppression and persecution. This led to a dramatic expansion of the Kutupalong refugee camp within a couple of months and a decrease of vegetation in the surrounding forests. As many humanitarian organizations demand frameworks for camp monitoring and environmental impact analysis, this study suggests a workflow based on spaceborne radar imagery to measure the expansion of settlements and the decrease of forests. Eleven image pairs of Sentinel-1 and ALOS-2, as well as a digital elevation model, were used for a supervised land cover classification. These were trained on automatically-derived reference areas retrieved from multispectral images to reduce required user input and increase transferability. Results show an overall decrease of vegetation of 1500 hectares, of which 20% were used to expand the camp and 80% were deforested, which matches findings from other studies of this case. The time-series analysis reduced the impact of seasonal variations on the results, and accuracies between 88% and 95% were achieved. The most important input variables for the classification were vegetation indices based on synthetic aperture radar (SAR) backscatter intensity, but topographic parameters also played a role.

1. Introduction

Remote sensing has become increasingly popular in the field of humanitarian action because it is an independent and reliable source of information and allows both a quick response to emergencies and the monitoring of gradual changes [1]. This is of particular importance when observations in the field area are not possible because of limited budget, legal restrictions, or security reasons [2,3,4]. Observation of specific places from space is not only crucial for decision making involving cases of humanitarian response to natural disasters or emergencies [5,6], but it also helps to develop a better general understanding of an area and how trends and temporal dynamics have shaped the spatial patterns of the present [7,8]. This temporal aspect is often more valuable than the level of spatial detail of the images alone, as it puts the information into a context and increases the interpretability of the retrieved information [9].

Institutionalized in the International Charter on Space and Major Disasters, a global collaboration of space agencies and disaster response units provides satellite images and value-added products to relief organizations to assist their operations [10]. However, other institutions and platforms also make operational use of satellite images for humanitarian assistance, such as UNOSAT (United Nations Institute for Training and Research (UNITAR) Operational Satellite Applications Programme, [11]), the United Nations Platform for Space-Based Information for Disaster Management and Emergency Response (UN-SPIDER, [12]), ReliefWeb [13], and the Copernicus Emergency Management Service (EMS) of the European Space Agency (ESA, [14]). However, these operational routines rely on scientific innovation and applied research. Many studies address the use of satellite imagery in humanitarian contexts. They range from short-term information needs, such as the retrieval of information on displaced persons and households based on very high-resolution (VHR) optical data [15,16,17,18], regarding the assessment of natural resources (e.g., firewood, groundwater) [19,20], developing long-term observations of environmental impacts [21,22,23,24,25], plus the documentation of human rights violations [26,27]. Especially in long-term applications, the availability of historical images is essential for humanitarian aid, and monitoring dynamics over long periods not only helps reconstruct ecosystem dynamics in an area, but also allows the development of disaster preparedness strategies, early warning systems, and risk adaptation measures for future events [28,29,30].

This study uses remote sensing to analyze the land cover changes related to forcibly displaced people in the Cox’s Bazar District of Bangladesh, close to the border of Myanmar (the historical background and a map of the area are provided in Section 2.1). The large amounts of people arriving since August 2017 (see Section 2.1) caused the massive expansion of the refugee camp of Kutupalong in 2017 and a considerable loss of forest areas. This case was investigated and documented in various research articles and reports. The most comprehensive and holistic was provided by the United Nations Development Program (UNDP) in 2018 [31]. It outlines the historical development of the humanitarian crisis and its effects on natural resources, ecosystems, and human health for the whole Cox’s Bazar District. Satellite images of Pléiades were used to qualitatively visualize the landscape changes around camp Kutupalong between 2002, 2009, and 2017. It stated that 1500 of the 26,000 hectares of the Ukhia, Whykong, and Teknaf forests were cleared because of the expansion of the camp. Of this impacted area, 53% were natural forests and 47% were plantations. These numbers are linked to biomass estimations and different forms of fuelwood consumptions based on findings from a study by the Food and Agriculture Organization of the United Nations (FAO), which made use of field measurements and a time-series of eight SPOT images between 1994 and 2016 [32]. Because fuelwood collection and bamboo use have high social relevance in this area, these points were specifically considered in environmental management and monitoring systems suggested in the UNDP report. Although little is mentioned on the methods of environmental mapping, it pledges for a monitoring of forest areas, growing stock, fragmented areas, areas affected by fire and soil erosion, and population density by means of remote sensing.

A work structure for an environmental impact assessment for the area was proposed by Honeth et al. [33] who suggest a classification system to quantify the magnitude of the impact for the infrastructures, food and water security, health, social issues, and land use sectors. They also underline the demand for energy as a severe threat to the forest areas and assume a minimum of 600 hectares cleared per month. Without monitoring and resource management, they expect a drastic loss of livelihood and an increased vulnerability to natural hazards.

Both reports highlight the demand for a consistent and continuous monitoring of forest areas, as field surveys are expensive and are often influenced by the subjectivity of the data collectors. To tackle this issue, remote sensing-based studies of this area were conducted, which are introduced in the following section. Imtiaz [34] employed products of Landsat 8 to quantify the loss of forest areas in the Teknaf sub-district based on the thresholding of the Normalized Difference Vegetation Index (NDVI). Images acquired in December of four consecutive years (2014–2017) were used to minimize the effect of seasonal variations regarding vegetation cover. The study showed a total loss of vegetation cover of 1284 hectares for the entire Teknaf sub-district, and a decrease of 102 hectares in the Teknaf Wildlife Sanctuary for the period of the incident (2016 to 2017), while the changes between 2014 and 2016 were widely negligible. By using the NDVI as an indicator for vegetation cover and a fixed threshold of 0.1 over which a pixel was considered vegetation, the study design is transparent and reproducible. However, no map on the spatial differentiation of changes is included and no accuracy assessment was conducted to estimate the quality of the results.

Labib et al. [35] compared images of Landsat 7 from March 2017 and December 2017 to investigate the decrease of forest area as a consequence of the extension of the camp. They derived brightness and greenness products and conducted a change vector analysis (CVA) between the two dates to extract the magnitude and direction of change between the three classes of water, forest, and no-forest. They calculated an area of 573 hectares for the transition between forest and camp, accounting for around 11% of the study area. The spatial implications of the changes are presented in a map, but the impact of seasonal variations on the result cannot be determined. Again, the methods are described in a transparent way, but no accuracy assessment is conducted.

A more comprehensive study, provided by Hassan et al. [23], compared two Sentinel-2 images of December 2016 (pre-influx) and December 2017 (post-influx) in a post-classification change detection approach. They classified the whole Teknaf sub-district into forest, camp, non-forest, and water bodies based on a random forest classifier to conduct an independent accuracy assessment based on field validation campaigns and a very high-resolution (VHR) QuickBird image, resulting in overall accuracies of 95% for both dates. Their reference data were collected during a field campaign in December 2017 to increase the independence of the validation. The study reported an increase of the camp area of Kutupalong–Balukhali of 1219 hectares and a forest loss of 2283 hectares for the investigated period. However, persistent cloud cover prevented them from using images during the rainy season to put their findings into a phenological context. Nonetheless, they highlight the importance of seasonal variations based on a comparison of NDVI between February 2017 and December 2017 along two buffered transects in the study area. The study concludes with a discussion on the environmental and socio-economic implications of the observed landscape changes.

All mentioned studies applied classical pixel-based techniques, which are often considered insufficient for the analysis of complex ecosystems because the spatial resolution of the sensor dictates the scale at which phenomena are described [36]. This was already declared by Woodcock and Strahler [37], who described how scale impacts the results of digital image analysis and that every environment has to be analyzed in a suitable spatial resolution. It also becomes evident in these studies that the NDVI is not always a suitable proxy for vegetation cover for various reasons: Firstly, indices and their relationship with biophysical parameters were observed to be highly scale-dependent [38,39]. This can especially become a problem when images of multiple sensors with different spatial and spectral resolutions are compared within a multi-temporal approach [40]. Secondly, there is an impact of soil reflectance which distorts the measure of green vegetation as anticipated by the NDVI [41]. For example, Todd et al. reported that targets with moist soil background showed higher NDVI values than dry targets with the same vegetation properties [42]. Especially in areas with pronounced seasonality where the green canopy cover changes during the year, such as Bangladesh in this case, measures that are more independent from spatial resolution and humidity have to be found.

Besides these spectral limitations, all studies name cloud cover as one main constraint in terms of data availability. The impacts of cloud cover and seasonal land cover variations on forest change detection are also widely discussed in various other studies [8,43,44,45]. For this reason, synthetic aperture radar (SAR) becomes an alternative for monitoring the environment, as it employs microwaves which are largely independent from cloud cover or daylight, and thus valid images that can be utilized for multi-temporal approaches are potentially available over the whole year and at more regular intervals. The potential of SAR for humanitarian action was initially acknowledged by Wegmüller et al. [46], who mentioned the quick mapping of land cover and hazards as useful applications in this domain. Besides the penetration of cloud cover, they name acquisition costs and the increased information content as the main advantages of radar images. Amitrano et al. [47] also discuss the capabilities of multi-temporal radar analyses for humanitarian action and present short case studies on land cover mapping, crop monitoring, seasonal water body mapping, and the derivation of digital elevation models with radar interferometry. The first environmental impact assessments were provided by Braun et al., who used time-series of radar data to assess the long-term changes of ecosystems around the refugee camps Dagahaley in Kenya [48] and Djabal in Chad [49]. Their studies underline the increased information content of microwave images for the detection of ecosystem changes and discuss how using multi-temporal data reduces the bias of seasonal variabilities which often superimpose gradual long-term changes [50,51]. While most SAR sensors make use of either horizontal or vertical waves to measure the backscatter intensity of a surface, polarimetric radar systems send and retrieve signals both horizontally and vertically. This allows for the ability to determine the contribution of different scattering mechanisms to the measured signal, which mainly consists of surface, dihedral, and volume scattering [52]. This form of signal decomposition allows for the ability to distinguish different types of land cover within a single radar image [53]; however, only a few polarimetric SAR satellites are operational at the moment [54].

This study provides a new method of change detection in the context of anthropogenic influence of displaced people as the main factor and addresses the following challenges which were identified above: (1) to develop a framework based on radar images which is independent from cloud cover, (2) to maintain objectivity at the selection of reference areas even when field surveys are not possible because of limited budget and/or security concerns, (3) to perform a multi-scale approach to reduce the scale-dependency of the generated change detection, and (4) to make use of polarimetric information for better classification accuracies via the combination of two SAR sensors of different wavelength and polarization. It furthermore reveals which wavelengths, polarimetric features, and spatial resolutions are most suitable to describe landscape changes at the local to regional scale.

The findings of this work are relevant with respect to the Forest Policy and the 20 year Forestry Sector Master Plan approved by the government of Bangladesh in the 1990s that aim at a national forest cover rate of 20% [55]. This work can furthermore assist larger monitoring frameworks on the impact of refugee camps on their environments, such as the Social and Environmental Screening Procedure developed by the United Nations Development Programme (UNDP, [56]). It was conducted in the study area by UNDP and the Ministry of Environment and Forest (MoEF) of Bangladesh to assess both the potential and probability of risks threatening the local population and their social, economic, and environmental implications in the Cox’s Bazar district [31]. Monitoring land cover change is furthermore relevant in accordance with the Bangladesh Environmental Conservation Act 1995 (ECA, [57]) and the Wildlife Preservation and Security Act 2012, both which protect the Teknaf wildlife sanctuary that is part of the study area [58]. On an international level, the Climate Resilient Ecosystems and Livelihoods project (CREL [59]) aims to improve the livelihood and environment in Bangladesh by improved governance of natural resources and social forest management [60]. All mentioned legal regulations and research programs underline the demand for a reliable and consistent monitoring framework of forest resources in the study area.

Currently available datasets on forest cover provide global maps, but at annual rates with no consideration of seasonality, such as the ALOS forest/non-forest maps (provided by the Japan Aerospace Exploration Agency (JAXA) [61]). Others provide more detailed classifications, but at low spatial resolutions, such as land cover products derived from MODIS (provided by NASA, 500 m [62]) or PROBA-V (provided by CCI/ESA [63]). None of them are sufficient to capture the dynamics of forest covers within the investigated period.

2. Materials and Methods

2.1. Study Area and Data

The study area lies in the Cox’s Bazar district in south-eastern Bangladesh and covers the area around Kutupalong which lies close to the border to Myanmar (Figure 1). Kutupalong became the world’s largest refugee camp in late 2017 when the Rohingya minorities fled from the neighboring Rakhine state, Myanmar, as a result of religious and ethnic oppression, which culminated in violent crackdowns and systematic killings beginning on 25 August 2017 [64]. Within a couple of weeks, the camp—which has existed since 1991—grew dramatically from around 100,000 to 440,000 inhabitants, and another 170,000 people settled down in spontaneous makeshift dwellings in the surrounding areas [31]. In 2018, the Kutupalong–Balukhali expansion (KBE) provided shelter for 550,000 Rohingya refugees [35]. These developments led to a massive decrease of forest cover and agricultural lands, which were converted to areas for housing, sanitation, and infrastructure facilities. However, this also endangers ecologically valuable areas, such as the Teknaf Wildlife Sanctuary (TWS), the proposed Inani National Park, and the Himchari National Park, which has already suffered from severe deforestation [31]. As outlined in the previous section, the aim of this study is to monitor these developments in a transparent and reproducible way.

To increase the potential of the radar images’ information, missions of Sentinel-1 (S1, ESA) and ALOS-2 (A2, JAXA) are combined in this study. This allows the measuring of backscatter intensity in all four polarizations (horizontal–horizontal (HH), horizontal–vertical (HV), vertical–horizontal (VH) vertical–vertical (VV)). Moreover, the different wavelengths of the sensors can be exploited more effectively—the C-band of Sentinel-1 operates at a comparably short wavelength (5.5 cm) which is strongly sensitive to the roughness of canopies and moisture of un-vegetated surfaces [65]. In turn, the large wavelength of ALOS-2 (22.9 cm) is able to penetrate surfaces to a certain extent and can therefore be used to describe structures under canopies, such as stems or branches [66]. This was already confirmed by Inoue et al., who report that the C-band is a more suitable canopy and grassland descriptor, while the L-band has more informative value regarding the biomass of vegetation [67]. Consequently, the combination of both sensors allows for the ability to analyze backscatter mechanisms of both surfaces and covered structures. This complementary use was already successfully demonstrated for biomass estimation [68,69], deforestation monitoring [70,71], the mapping of flooded vegetation [72], and rice or palm oil cultivation [73,74]. However, it has to be denoted that, while Sentinel-1 data can be freely accessed within the Copernicus Programme, ALOS-2 data is commercially distributed by the Japan Aerospace Exploration Agency (JAXA). This aspect will be discussed in Section 3.3.

As indicated in the studies presented in the previous section, the study area shows high phenological variations over the year. To be able to take these into consideration and prevent the deduction of false conclusions, SAR image pairs of S1 and A2 for eleven dates between 2016 and March 2019 (

Table 1. Synthetic aperture radar (SAR) data used in this study.

Month/Year	Season	Sentinel-1	ALOS-2	∆t 1	Orbit
Mar 2016	dry	21.03.2016	14.03.2016	7	descending
Jul 2016	rainy	29.07.2016	22.07.2016	7	ascending
Oct 2016	dry	09.10.2016	30.09.2016	9	ascending
Feb 2017	dry	18.02.2017	17.02.2017	1	ascending
Jun 2017	rainy	08.06.2017	05.06.2017	3	descending
Jul 2017	rainy	24.07.2017	21.07.2017	3	ascending
Feb 2018	dry	01.02.2018	02.02.2018	1	ascending
Jun 2018	rainy	03.06.2018	04.06.2018	1	descending
Jul 2018	rainy	19.07.2018	20.07.2018	1	ascending
Jan 2019	dry	15.01.2019	18.01.2019	3	ascending
Mar 2019	dry	04.03.2019	01.03.2019	3	ascending

¹ temporal baseline between both dates in days.

) were selected. As shown in Table 1, most of the pairs had a temporal baseline of three days or less. Furthermore, only images with the same flight direction (ascending or descending orbit) were selected to minimize the impact of topography within one investigated date.

To quantify the decrease of forest areas during the investigated period, a supervised classification was chosen (see Section 2.2). However, this required the input of training areas whose collection was either not possible in the field or biased by visual interpretation of existing images (see challenges outlined in Section 1). For this reason, optical data of Landsat ETM+ (L7), Landsat OLI (L8), and Sentinel-2 (S2) were selected on the basis of their acquisition date and amount of cloud cover (

Table 2. Optical data used in this study. L7 = Landsat ETM+, L8 = Landsat OLI, S2 = Sentinel-2.

Month/Year	Sensor	Date	Coverage 1	Cloud Cover 2
Mar 2016	L7	04.03.2016	63.0%	0.0%
Jul 2016	L8	27.07.2016	84.1%	68.7%
Oct 2016	L8	22.10.2016	100%	0.0%
Feb 2017	S2	18.02.2017	100%	0.00%
Jun 2017	L8	03.06.2017	100%	23.7%
Jul 2017	L7	13.07.2017	46.9%	69.2%
Feb 2018	S2	13.02.2018	100%	0.0%
Jun 2018	L8	30.05.2019	88.0%	13.0%
Jul 2018	S2	18.07.2019	94.5%	3.2%
Jan 2019	S2	19.01.2019	100%	0.0%
Mar 2019	S2	10.03.2019	100%	0.0%

¹ terrestrial part of the study area which is covered by valid image pixels; ² part of the valid image pixels which is covered by clouds.

). Even if some images had gaps because of the L7 SLC-off incident since 2003 [75], were partly covered by clouds, or did not fully cover the study area, they were still usable to retrieve training areas as described in the following section.

2.2. Data Processing

2.2.1. Overall Study Design

The overall workflow of the study is displayed in Figure 2. The single components will be explained in more detail in the subsequent sections, but the general approach is shortly outlined as follows. Spectral indices were derived on the optical images to serve as an input for a semi-automated pre-classification of the four target classes for each date (Section 2.2.2). The following classes were selected for the study area: (1) Vegetation—around 40% of the study area is covered by vegetation, mainly consisting of mixed tropical evergreen and semi-evergreen forest dominated by tall Garjan (Dipterocarpus spp.) trees [76]. Peaks of the hills are predominantly covered by sun grass (Imperata cylindrica), herbs, shrubs, and brush woods [55], but their share in the study area is low because these elevations are only reached further in the north and east (as indicated in Figure 1). (2) Open lands (abbreviated as “open” in the following sections)—includes all areas outside the forests which are only sparsely covered by grassy plants, barren lands, or degraded forests affected by shifting cultivation as identified by Giri and Shrestha [77]. (3) Water—the course and floodplains of the Naf River in the southeast of the study area. This also includes the frequently waterlogged marshes in this valley. (4) Informal settlements (abbreviated as “camp” in the following sections)—the study area only contains the refugee camp of Kutupalong and its extension areas. This class represents all areas that are covered by solid humanitarian facilities and simple dwellings of lighter construction materials.

Radar images of S1 and A2 were combined to compute a large feature space for the supervised classification (Section 2.2.4). As many of the target classes are also dependent from topographic conditions, a digital terrain analysis was conducted to provide information on landforms at several scales (Section 2.2.3). The training areas were then applied to the input features to generate a land cover prediction map for all eleven dates to identify the changes between 2016 and 2019 (Section 3.1). An accuracy assessment was performed based on validation samples which were independently extracted from the pre-classification (Section 3.2). Finally, the importance of all input features was evaluated to estimate their contribution to the resulting land cover classification (Section 3.3).

2.2.2. Processing of Optical Images and Retrieval of Training Areas

As shown in Table 2, some of the optical images closest to the acquisition dates of the SAR image pairs did not cover the full study area and unavoidably contained cloud cover. They were, however, still usable to derive training areas for the latter classification (see the next section). In the first step, cloud masking was conducted on the images based on IdePix (version 2.2, [78]) for Sentinel-2 data and the fmask algorithm (version 4.0, [79]) for Landsat products. In the second step, radiometric calibration of the images to retrieve bottom of atmosphere (BOA) surface reflectance using the calibration constants provided for Landsat by Chander et al. [80] and the sen2cor processor (version 2.5.5 [81]) for Sentinel-2. All image products were harmonized to a common spatial resolution of 10 m using bilinear resampling. Compared to nearest neighbor resampling, this technique alters the original pixel values and leads to a smoother-appearing image that can lead to problems in subsequent pattern recognition techniques [82]. As the optical images are used in this study to derive radiometric indices for the identification of small representative areas of extensive land cover classes, the smoothing might even contribute to larger homogenous areas.

The following radiometric indices were computed as a feature space for the automated extraction of training areas for all investigated dates:

• Normalized Difference Vegetation Index (NDVI [83], Equation (1))—This is the most frequently used vegetation index and is based on the red and near infra-red (nIR) band. It highlights the vitality of vegetation based on their chlorophyll content and was used in this study to distinguish vegetated from open or degraded areas.
• Enhanced Vegetation Index (EVI [84], Equation (2))—This index also addresses the identification of vegetation quality, but it was additionally selected for this study because it is said to have a higher sensitivity in high biomass areas. It is based on red, blue, and nIR, as well as additional constants—C1 and C2 are used to estimate the aerosol resistance, while G is an empirical gain factor, and L allows for the adjustment of the impact of soil.
• Modified Normalized Difference Water Index (MNDWI [85], Equation (3))—This index is based on a ratio between the reflectance in the green band and the middle infra-red (mIR) band to highlight areas that are covered by water. In this study, it was used to delineate the training areas for water.
• Built-up Index (BUI [86], Equation (4))—This rather complex index makes use of the red band and the two short-wave infra-red bands (swIR1, swIR2) and was used to highlight built-up surfaces for the delineation of the training areas of the refugee camp.
• Index-based Built-up Index (IBI [87], Equation (5))—This index consists of the three previously mentioned indices and is said to be more sensitive to built-up areas of different densities, which is important in this setting as parts of the refugee camp consist of different dwellings and construction materials. [31].
• Normalized Difference Build-up Index (NDBI [88], Equation (6))—This rather simple index uses the near and short-wave infra-red bands (nIR, swIR) which also highlight urban areas. It was additionally selected to mitigate the different center wavelengths and bandwidths of Landsat ETM+, Landsat Oli, and Sentinel-2.

NDVI = (nIR − red)/(nIR + red)

(1)

EVI = G × ((nIR − red)/(nIR + C1 × red − C2 × blue + L))

(2)

MNDWI = (green − mIR)/(green + mIR)

(3)

BUI = 2× ((red × swIR2) − (swIR1 × swIR1))/(red + swIR1) × (swIR1 + swIR2)

(4)

IBI = ((2 × NDBI) − (NDVI + MNDWI))/((2 × NDBI) + (NDVI + MNDWI))

(5)

NDBI = (swIR1 − nIR)/(swIR1 + nIR)

(6)

To minimize the effect of subjectivity and manual input in the study design, 25 points were digitized for each class, based on visual interpretation of each cloud-masked image. These 100 points were then used as input for a Classification and Regression Tree (CART [89]) to retrieve a coarse pre-classification of the four classes in each image. The trees were pruned to a maximum depth of six levels, wherein the minimum number of instances per leaf was two, and subsets containing less than five members were not split any further. This resulted in training accuracies between 87.3% and 96.9%. As shown by the trees for the single dates, which are displayed in Appendix A, the use of radiometric indices allowed for nodes to be achieved with high class purity within three levels for most of the dates, and six levels at maximum based on a small set of input rasters. The areas resulting from these classifications were then used as training inputs for the supervised classification (Section 2.2.4). Their sizes per analyzed date are summarized in

Table A1. Size of the training areas retrieved from the CART in square kilometers.

Month/Year	Open	Vegetation	Water	Camp
Mar 2016	10.62	12.57	7.50	0.25
Jul 2016	2.81	4.46	2.10	0.29
Oct 2016	4.38	27.45	2.00	0.33
Feb 2017	19.70	39.12	9.25	0.57
Jun 2017	10.08	8.00	2.74	0.71
Jul 2017	7.44	3.93	4.47	1.03
Feb 2018	25.55	26.50	12.95	16.29
Jun 2018	11.26	8.94	1.71	3.10
Jul 2018	16.68	13.80	1.65	1.42
Jan 2019	37.28	19.43	8.22	18.12
Mar 2019	16.87	22.30	8.49	23.50

(Appendix B).

2.2.3. Processing of Topographic Input Features

As demonstrated by Chust et al., the integration of topographic features can significantly increase the quality of land cover classifications [90]. This was also proven in other studies and can be explained by the fact that the presence or absence of certain classes of land cover or land-use are strongly determined by elevation, slope, or aspect of surfaces [91,92,93]. For example, forest cover ends at certain elevations, water bodies only occur on flat terrain, and windward slopes are often characterized by different vegetation communities and soil types instead of downwind sides [94,95,96]. Furthermore, both optical and radar images underlie topographically-induced radiometric distortions in terms of illumination and orientation [97,98]. As shown in Figure 1, the study area shows pronounced topographical patterns, such as the ridges of over 100 m elevation ranging in the north-south direction within a narrow zone from the coast. To be able to describe the land cover with respect to these topographic conditions and to enhance the feature space of the supervised classification, a digital elevation model (DEM) at a spatial resolution of 30 m of the Shuttle Radar Topography Mission (SRTM [99]) was used to derive the following parameters: elevation, slope, northness, eastness, plan curvature, and profile curvature by means of digital terrain analysis [100], and a classification of morphometric features (peaks, ridges, passes, channels, pits, and planes [101]). In the first step, the data were resampled to a spatial resolution of 10 m using bilinear resampling to match the resolution of the final classification (see the next section). As these parameters are highly scale-dependent, this analysis was performed for window sizes of 1 × 1, 5 × 5, 9 × 9, 15 × 15, 25 × 25, and 49 × 49 pixels, leading to a total number of 42 topographic features.

2.2.4. Processing of Radar Input Features and Image Classification

Radar data of both S1 and A2 were radiometrically calibrated to Beta0 by applying calibration coefficients from the metadata [102,103] to retrieve radar brightness. Using the digital elevation model described in Section 2.2.3, terrain-flattened Gamma0 was derived, as introduced by Small et al. [98], to minimize the impact of the incidence angle on backscatter intensity. A refined Lee filter [104] with a window size of 5 × 5 pixels was applied to all images to reduce speckle effects. This adaptive speckle filter produces smoother images in both high and low backscatter areas and provides a good compromise between image quality and computing time [105]. For each date, image products of S1 and A2 were co-registered to a stack, based on their orbit information and geocoded using range Doppler terrain correction [106], so that data for each date consisted of all four polarizations (HH and HV from ALOS-2; VH and VV from Sentinel-1) at a common spatial resolution of 10 m.

To enhance patterns emerging at several spatial scales, the pre-processed SAR stacks were filtered by median filters with windows sizes of 5 × 5, 9 × 9, 15 × 15, 25 × 25, and 49 × 49 pixels, leading to a total number of 24 radar bands per analyzed date. As all four polarizations were available, the following radar measures were additionally derived for each date: Radar Vegetation Index (RVI [107], Equation (7)), Radar Forest Degradation Index (RDFI [108], Equation (8)), Canopy Structure Index (CSI [109], Equation (9)), average like-polarized magnitude (LK [109], Equation (10)), average cross-polarized magnitude (CS [109], Equation (11)), Volume Scattering Index (VSI [109], Equation (12)), standard deviation (StD, Equation (13)) and three combined indices (CI1, CI2, and CI3, Equation (14)), which are leaned towards the ratio of different polarizations as defined by the Pauli scattering elements [52]. Technically, all these indices can only be derived from fully-polarized radar data taken by one sensor within a single acquisition because they are sensitive to the orientation of objects on the measured surface [53]. This physical consistency is not given in our study, as the images of Sentinel-1 (VV and VH) and ALOS-2 (HH and HV) were not acquired at the same time, nor do they use the same frequency. Their interpretation regarding the initial use, namely the inherent physical characteristics of surfaces, should be made with respective care. Examples of the indices used in this study and how they are affected by the median filter products of different window sizes are given in Figure 3. It is nicely demonstrated how Figure 3A (based on median filter over 5 × 5 pixels) shows a higher level of detail, while Figure 3B (25 × 25 pixels) results in more homogenous areas favoring the classification of large patches of land cover. Figure 3C highlights the forest areas, represented in green, while the variations in vegetation density are well characterized by the CSI (Figure 3D).

RVI = (8 × HV)/(HH +VV + 2 × HV)

(7)

RFDI = (HH − HV)/(HH + HV)

(8)

CSI = VV/(VV + HH)

(9)

LK = (HH + VV)/2

(10)

CS = (HV + VH)/2

(11)

VSI = CS/(CS + LK)

(12)

StD = std(HH, HV, VH, VV)

(13)

CI1 = HH + VV; CI2 = HV + VH; CI3 = HH − VV

(14)

The feature space for the supervised classification consisted of 126 features per investigated date (84 radar features and 42 constant terrain features). They represent different types of landforms and land cover at multiple spatial scales, but since they are derived from a small amount of base data (four polarizations and one DEM), they are correlated with each other. To identify the features that are most suitable for the description of the study area, a random forest (RF) classifier was selected. As a method from the machine learning domain, it originates from simple classification trees, but additionally introduces a random component by repeatedly selecting random subsets from both the training samples and the input features within a defined number of independent iterations [110]. For each investigated date, 150 decision trees were calculated on the basis of one third of the available training pixels per class (Section 2.2.2), and twelve feature rasters (square root of the total 126 available features, potentially consisting of both radar and topography features). Accordingly, each decision tree used a different subset of training pixels and features to predict the distribution of the target classes (vegetation, open, water, camp) by a CART-like hierarchical ruleset. This led to a collection of 150 uncorrelated decision trees that were then used to determine the final classification of each pixel on the basis of the class majority of all iterations (Section 3.1). This had several advantages: Firstly, redundancies in the input features, as they were used in this approach, had no negative effect on the classification [111]. This is especially valuable in cases where the relevance of these features is not known in the beginning [112]. Additionally, smaller errors in the pre-classification derived from the spectral indices (Section 2.2.2) will not have a large impact by falsely training the classifier [110]. As a result of the permutation of the input features, their contribution to the final classification can be estimated, which not only identifies which features are most suitable to describe the study area, but also at which spatial scale (Section 3.3).

3. Results and Discussion

3.1. Land Cover Change

The land cover maps for all eleven dates are shown in Figure 4. Most strikingly, the area of the camp and the expansion of the camp after July 2017 is dominant throughout all images. The camp area increased from 36 hectares in early 2016 to over 1850 hectares at the beginning of 2019. The highest increase is observed between July 2017 and February 2018, where the camp area increased tenfold within only a few months (according to the reports, the largest immigration was between August and November 2017 [64]). Consistent with the findings of Hassan et al. [23], who reported an increase of 1219 hectares between December 2016 and December 2017, the findings of this approach seem plausible at first sight (camp growth of 1313 hectares between February 2017 and February 2018).

Figure 4, furthermore, shows the impact of seasonal variations regarding vegetated and open areas—the only image pair found of S1 and A2 at the end of the rainy season was for October 2016. As shown in the figure, the classified vegetation areas dominate the study area, while vegetation coverage showed more differentiated patterns than in all other maps. As later shown in Section 3.2, this was not caused by low classification accuracy, but rather by the strong presence of vegetation, composed not only of trees, but also bushes and shrubs that are fully developed at the end of the rainy season and therefore additionally contribute to volume scattering. The remaining land areas were correspondingly classified as open areas. Inland water bodies showed larger spatial variations because of wetlands, temporarily flooded oxbows, and creeks along the Naf River in the southeast of the study area. Seemingly all these landscape features show similar backscatter characteristics and strongly vary between the rainy and dry seasons. However, this is of no concern for the assessment of camp growth and vegetation dynamics because they are rarely covered by forest and also have no suitable building ground [33].

The strong impact of phenology is also underlined in Figure 5A. It demonstrates the variation of vegetation cover during the investigated period in hectares (green line, primary y-axis) with the monthly precipitation for this area as retrieved from the Global Precipitation Climatology Centre (GPCC [113], blue areas, secondary y-axis). Although the vegetation cover might show a slightly negative trend (green line), the dynamics resulting from seasonal rainfall availability superimposed these patterns at this spatial scale. This underlined the indication that seasonal variations, or more precisely speaking, the date of the used images, had a large impact on temporal analyses in this area, and that studies using single images for comparison (see Section 1) potentially underlie a massive phenological bias that can lead to false conclusions. The water bodies in the study area (Figure 5A, blue line) underlay the variations discussed above, but showed an overall stable development. Clearly visible is the growth of the camp after July 2017 as a result of the refugee influx starting in August 2017. Unfortunately, no further image pairs of S1 and A2 are available before February 2018 for a more detailed temporal delineation of the camp expansion, especially for the period between 25 August and 11 September 2017 when the largest influx was reported [64].

To highlight the camp growth, but also to show its impact on the direct surrounding environment, Figure 5B shows the progress of the classes vegetation, open, and camp for an area within a 2 km distance from the main development axis of the camp (delineated in Figure 3D as a white dashed line for the axis and a gray dotted line for the area). It includes all portions of the main camp and its main expansion area in their present extent, but also the Balukhali expansion in the south. The decrease of forest areas became more clearly visible, and even though vegetation cover in October 2016 was above-average compared to the other investigated scenes, a strong decrease remained. What is interesting is that forest areas had already decreased in favor of open areas in the beginning of 2017, accordingly, before the camp expansion happened in August 2017. This can be partially explained by a preceding migration movement from Myanmar into the Cox’s Bazar district in October 2016, which was denoted as a “significant influx of people” by the Inter Sector Coordination Group [64]. This movement resulted in 74,000 arrivals until February 2017, of which 31,207 were registered in the Kutupalong makeshift settlement (KMS) and reportedly “thousands of newly built huts have appeared in the hills adjacent to KMS“ [114]. Thus, it is possible that the deforestation of the area close to Kutupalong already started before the massive influx of Rohingyas in August 2017 as a result of increased demand for construction materials to provide shelter and accommodation for the arrivals. This indication is underlined by the drastic increase of open areas around the center of the camp. Temporally delayed, these open areas were then converted to densely-populated camp sites after July 2017.

The decrease of forest areas in the direct surroundings of the camp matches the 2283 hectares reported by Hassan et al. [23] when the dates of July 2016 and February 2018 are compared (decrease of 2100 hectares), excluding the data from October 2016 (extraordinarily high amount of vegetation). As expected, the camp area increased the most between July 2017 and February 2018 (+1314 hectares) when most of the Rohingya crossed the border (see Section 1). A small reduction of camp area is observed between June and July 2018, which can partially be credited to demolitions “damaging thousands of tents” by monsoon rains in early June [115]. Nonetheless, it is most likely that smaller misclassifications played a role as well. They will be discussed in the next section.

Comparing the first and the last image utilized in this study (both acquired in March) a total loss of 1562 hectares of vegetation can be attested, of which 20% are now camp areas and 80% are open areas. This testifies both to the impact of the camp, but also the overall decrease of vegetation during the investigated period due to increased resource demand [33]. The camp area itself increased by 1830 hectares, which is demonstrated in Figure 6. It shows the largest transitions of vegetation to open areas (bright yellow) and the camp site (bright red). Most of the decrease happened in the west and south of the initial camp area (dark red), while the most stable areas of vegetation and open land are along the ridge near the coast and on the side of Myanmar.

A different form of change visualization is presented in Figure 7. This approach was originally suggested by Langer [116,117], who investigated landscape changes around the Lukole refugee camp in Tanzania based on a time-series of 42 Landsat images. For this study, squares of the size of 300 m width and height were generated, and the dominant landcover was assessed for each square and date on the basis of the determination of each square’s mode. The square was then divided into nine sub-squares the size of 100 m to which the land cover classifications between July 2016 and January 2019 were assigned. Accordingly, the change of landscapes can be visualized within each square of 300 m, as demonstrated in Figure 7. The initial area of the camp is characterized by consistently red squares, while the expansion sites south of it are gradually colored red in the lower parts of the squares. Additionally, the areas of conversion from forest to camp are well identifiable in the map on the right side (pink tones). They make up about 630 hectares, which correlates well with the 572 hectares of forest decrease reported by Labib et al. [35], who estimated the conversion between March and December 2017. The trends in Figure 7 (right) were calculated based on the first three (July 2016 to February 2017) and last three (June 2018 to January 2019) images of each square—it shows areas where vegetation has permanently changed to camp (light pink) and or to open areas (violet). The latter could also be some indication for transitions to camp areas, especially in the southwest and south of Kutupalong where new settlements are emerging [115]. Most of the forested areas of the Teknaf Wildlife Sanctuary in the south of the study area are characterized by permanent vegetation cover (similar to the areas shown in Figure 6), but two sub-squares indicate the clearance of forests, probably as a result of the population increase, as reported by Imtiaz [34].

3.2. Accuracy Assessment

Because decisions are to be based on the findings of this study and other studies on the forest change, as outlined in Section 1, the accuracy of the generated output maps must be known [118]. It strongly determines how well conclusions can be drawn from the retrieved land cover dynamics and at which points there might be weaknesses. Assessing the accuracy of digital maps was based on independently collected validation points which were sampled according to specific rules. One way to determine the required number of samples N according to the expected accuracy and confidence level is given by Fitzpatrick-Lins [119] (Equation (15)) based on the binomial distribution formula:

$$\mathrm{N}=\frac{Z^2\times \left(p\ast q\right)}{E^2},$$

(15)

where Z is the standard normal deviate for a 95% confidence level (here: 2), p is the expected accuracy of the map (here: 90%), q is 100 − p (here: 10), and E is the allowable error (here: 2.5%). This led to a required number of 576 points per analyzed date. As the different land cover classes were not equally distributed over the study area, we applied a stratified random sampling so that the number of samples represented the spatial share of each class—250 points for vegetation and 100 points for water (in accordance with the 41% of forest cover and 19.75% of floodplain in the greater study area, as reported by Moslehuddin [55]). One hundred and fifty points were sampled within areas classified as open, as the second largest class and the remaining 76 points were sampled within the camp areas, which did not represent their spatial share in the resulting map, but it is granted that enough sampling points lay within this class to provide a robust accuracy assessment. The points were automatically generated within each of the 11 analyzed maps and manually assigned to a class on the basis of visual interpretation of the corresponding optical image (Table 2) using the AcATaMa plugin (version 19.6.20, [120]) in accordance with the good sampling practices published by Olofsson et al. [121]. A confusion matrix was then computed for each date (

Table A2. Accuracy assessment of the land cover classification.

2016–03	V	O	W	C	2016–07	V	O	W	C
V	237	11	1	1	V	240	9	0	1
O	33	109	8	0	O	18	131	1	0
W	1	8	91	0	W	2	7	91	0
C	2	3	0	69	C	1	3	0	72
2016–10	V	O	W	C	2017–02	V	O	W	C
V	248	2	0	0	V	249	1	0	0
O	15	124	11	0	O	12	135	3	0
W	0	5	95	0	W	0	3	97	0
C	2	2	0	72	C	6	3	0	67
2017–06	V	O	W	U	2017–07	V	O	W	C
V	237	11	2	0	V	241	6	3	0
O	17	132	1	0	O	24	125	1	0
W	1	4	95	0	W	4	11	85	0
C	3	4	0	67	C	2	4	2	68
2018–02	V	O	W	U	2018–06	V	O	W	C
V	231	16	3	0	V	235	15	0	0
O	26	124	0	0	O	9	136	2	3
W	3	0	97	0	W	0	2	98	0
C	3	2	1	70	C	1	1	0	74
2018–07	V	O	W	U	2019–01	V	O	W	C
V	236	14	0	0	V	234	13	1	2
O	13	137	0	0	O	9	134	7	0
W	1	7	91	1	W	1	6	93	0
C	3	8	0	65	C	0	0	0	76
2019–03	V	O	W	U
V	239	8	2	1
O	37	110	2	1
W	0	8	92	0
C	4	1	0	71

, Appendix C). On this basis, the following overall accuracies were calculated: Mar 2016: 88.2%, Jul 2016: 92.7%, Oct 2016: 93.6%, Feb 2017: 95.1%, Jun 2017: 92.5%, Jul 2017: 90.3%, Feb 2019: 91.0%, Jun 2018: 94.4%, Jul 2018: 92.0%, Jan 2019: 93,4%, Mar 2019: 89.9%.

Producer’s accuracy (PA) refers to the commission error and indicates how reliable the classification of a pixel is, while user’s accuracy (UA) refers to the omission error and indicates the completeness of the determination of a class [122]. Cohen’s kappa additionally determines how likely the agreement between validation samples and classification occurred by chance and is often claimed to be a more robust measure than the simple percent agreement [123].

Table 3. Error measures of the land cover classification.

	2016–03			2016–07			2016–10			2017–02			2017–06			2017–07
	PA	UA	K	PA	UA	K	PA	UA	K	PA	UA	K	PA	UA	K	PA	UA	K
V	86.8	94.8	0.8	92.0	96.0	0.7	93.6	99.2	0.9	93.3	99.6	0.9	91.9	94.8	0.9	89.0	96.8	0.8
O	83.2	72.7	0.7	87.3	87.3	0.8	93.2	82.7	0.8	95.1	90.0	0.8	87.4	88.0	0.6	85.1	84.0	0.7
W	91.0	91.0	1.0	98.9	91.0	0.9	89.6	95.0	1.0	97.0	97.0	1.0	96.9	95.0	0.8	95.5	84.0	1.0
C	98.6	93.2	1.0	98.6	94.7	1.0	100.0	94.7	1.0	100.0	88.2	0.9	100.0	90.5	1.0	100.0	89.5	1.0
	2018–02			2018–06			2018–07			2019–01			2019–03
	PA	UA	K	PA	UA	K	PA	UA	K	PA	UA	K	PA	UA	K
V	87.9	92.8	0.8	95.9	94.4	0.9	93.7	94.4	0.9	95.9	93.6	0.8	85.2	96.8	0.8
O	87.3	82.7	0.8	88.3	90.7	0.8	82.1	92.0	0.7	87.0	89.3	0.7	87.5	74.7	0.9
W	98.0	98.0	1.0	98.0	98.0	0.9	100.0	92.0	1.0	94.0	94.0	0.9	100.0	92.0	1.0
C	100.0	92.1	0.9	97.4	97.4	1.0	100.0	84.2	1.0	97.4	100.0	1.0	100.0	94.7	1.0

V = vegetation, O = open, W = water, C = camp; PA = producer’s accuracy [%], UA = user’s accuracy [%], K = kappa.

shows that the high vegetation cover in October 2016 was not a result of over-estimation, as it can be attributed to the time of image acquisition which lies at the end of the rainy season with temperatures above 30 °C and precipitation of up to 1000 mm per month, leading to a pronounced vegetation cover [113]. Nonetheless, the PA of vegetation of 93.6 indicated a slight over-estimation of this date, but it was comparable to those of other dates, ranging from 85.2% to 95.9% in total for the vegetation class. Accordingly, the classified vegetation cover can be trusted throughout the study to a high degree, which was additionally confirmed by the kappa values between 0.8 and 0.9.

Throughout almost all investigated dates, open areas had the lowest producer’s accuracies, which led to the conclusion they were over-estimated in many cases. This was likely at the edges of large connected patches of vegetation, in areas where the vegetation cover is spatially fragmented, or areas which are generally characterized by lower vegetation, such as shrubs or bushes. The lowest overall accuracies observed in this study were achieved for 2016–03, 2017–07, and 2019–03, which were caused by falsely classified open areas (UA of 72.2%, 84.0%, and 74.7%, respectively) and under-estimated vegetation areas (PA of 86.8%, 89.0%, and 85.2%, respectively).

Lastly, the camp areas were classified with very high kappa values and low errors of commission and omission. This means the expansion of the camps and the accompanying decrease of vegetation not only corresponded to the values found in other studies (see Section 3.1) but they were also well described by SAR data in this approach, most probably due to an increase of dihedral scattering caused by double bounce effects at artificial structures.

One drawback to mention here is that the observed accuracy was a result of the study design, starting with the selection of 25 training samples for the pre-classification based on the radiometric indices. Although the fundamental idea of this study is based on the assumption that this number is enough for classifying an entire image based on radar and topographic features, it has to be noted that small changes in these training samples potentially had a large impact on the final results. This was partly because they defined the reference areas (as derived from radiometric indices) for the training of the actual land cover classification, leading to different rulesets within the classification algorithm. Moreover, their small number was sensitive to outliers which did not fully represent the respective class. As already indicated in Section 2.2.4, this risk was minimized by the massive permutation of input features within the random forest classifier, making this classifier robust to outliers and favoring features with a higher significance to distinguish the target classes [110]. Also, as shown by the independently computed overall accuracies, beginning with only 25 visually identified training points only can lead to very good results.

As for the stratified random sampling used in the accuracy assessment, it is suggested by Townshend and Justice that, instead of using single pixels, clusters of connected pixels (e.g., blocks of 3 × 3 pixels) should be analyzed as a more robust measure for the representation land cover classes and their spatial distribution of land cover classes [124]. While this form of sampling clearly reduces the costs for data collection at an early stage of the study [121], it increases the number of samples to be manually labelled nearly 10-fold in the later validation part.

To put the findings of this study into the context of the current scientific debate, we compared the land cover maps of this study to those created by Hassan et al. in 2018 [23]. Figure 8A shows the pre-influx agreement between the results from Hassan et al., based on optical images from December 2016), with the closest date analyzed in this study from February 2017. The overall statistical agreement was only 66.7% (pre-influx, Figure 8A) and 60.8% (post-influx, Figure 8B), which seems low at first sight. However, the comparison reveals some interesting points about the quality of this study that will be discussed as follows.

Figure 8 shows that the camp area was nearly the same in both studies, but larger parts, which were classified as vegetation in this study, were open areas in the study of Hassan et al. (Figure 8A bright blue areas). This can only partly be explained by an overestimation of vegetation areas in our study, as indicated by the producer’s accuracy of 93.3% (Table 3). More importantly, the temporal delay (here: 48 days) between both observations might already contain phenological effects caused by the rainy season gradually starting in February (Figure 5A), resulting in more areas considered as healthy vegetation. Lastly, both approaches were based on different input data—while the study of Hassan et al. used multi-spectral satellite data that mostly measured the amount of chlorophyll of a surface using infra-red information, this study made use of radar data, which is sensitive to backscatter mechanisms. That means that areas which had low chlorophyll values in the optical images might still be recognized as vegetation by the radar images because of their branch structures and crown surfaces.

The same applies to the post-influx comparison (Figure 8B), which also shows a considerable discrepancy between vegetated and open areas along their borders. Two further differences can be observed here: Firstly, the camp is larger in the western part in the map of Hassan et al. (Figure 8B dark blue). These areas concentrate at the western and southern borders of the main camp area (part of the Balukhali expansion area), but also of the smaller site in the south, and it is likely that they contain dwellings which are the most recently constructed. As resources of solid materials have become short since August 2017 [32], one explanation for their under-representation in this study could be that these buildings are no longer covered by corrugated metals, or at least to a smaller degree, as a consequence of material shortage and a shift to lighter construction materials, such as wood, leaves, or plastic (all used for covering the roofs [31]). Furthermore, less densely built parts of the camp were prone to under-estimations of this land cover class. This slight under-estimation was also indicated, but not fully explained by the user’s accuracies (92.1%) and kappa value (0.9) of the camp class from February 2018 because in the image used for training (13.02.2018) the classifier was closer to the radar image (01.02.2019) than the one used by Hassan et al. (15.12.2017). Nonetheless, the ambiguity in the backscatter intensity of solid and light dwellings can become a problem in studies dealing with humanitarian settings, as already discussed by Braun [125], however, more research is required to understand the differences between their radar signatures [126].

3.3. Feature Importance

Besides the land cover maps illustrated in Figure 4, the permutation importance (PI) and the Gini decrease (GD), as measured for the contribution of each feature (input raster), were calculated for the random forest classification. PI measured the mean decrease in the model accuracy by systematically comparing the effect of randomly shuffled input variables [127]. Accordingly, the higher a feature’s permutation importance, the more important it was for the respective class. GD made use of the Gini impurity, which determined how well a split separates the permuted samples at each node to estimate each variable’s importance on the basis of its occurrence throughout all trees [128]. As an index value, it had no upper boundary, but the higher the GD, the larger the importance of the variable for the classification. For an evaluation of the feature importance, the 500 features with the highest PI throughout all investigated dates were selected. Of the features, 82 were based on topographic parameters (Section 2.2.3) and 418 were SAR-based features (Section 2.2.4). This confirms the design of the study, which combined spectral information with topographic characteristics to predict the land cover classes, where the radar parameters played a significantly larger role. As shown in

Table 4. Frequency of the 500 most important features used for the classification. DEM = digital elevation model.

Feature	Feature Origin	Frequency
VH	SAR (S1)	85
VV	SAR (S1)	51
HH	SAR (A2)	43
Radar Vegetation Index (RVI)	SAR (S1 + A2)	40
Cross-polarized magnitude (CS)	SAR (S1 + A2)	38
Average like-polarized magnitude (LK)	SAR (S1 + A2)	26
Combined index (CI1)	SAR (S1 + A2)	23
HV	SAR (A2)	21
Standard deviation (StDev)	SAR (S1 + A2)	20
Radar Forest Degradation Index (RFDI)	SAR (S1 + A2)	20
Canopy Structure Index (CSI)	SAR (S1 + A2)	19
Volume Scattering Index (VSI)	SAR (S1 + A2)	18
Elevation	DEM	18
Aspect: northness	DEM	18
Morphometric feature	DEM	17
Slope	DEM	17
Combined index (CI3)	SAR (S1 + A2)	14
Profile curvature	DEM	6
Aspect: eastness	DEM	6

, the most frequent feature used throughout all the analyzed dates was VH polarization of Sentinel-1, followed by VV and HH. Accordingly, backscatter intensity of the SAR satellites already had a large informative value for the classification of land cover in the study area. However, indices such as the Radar Vegetation Index (RVI) and the Radar Forest Degradation Index (RFDI) strongly contributed to higher prediction accuracies. As expected, information of the cross-polarized channels (VH, HV, cross-polarized magnitude CS) played a larger role because of their sensitivity to volume scattering caused by vegetation. Table 4 also shows that all features based on topography were of secondary importance, most likely because they were stable throughout the whole period, while the seasonal variations could only be described by the changes of SAR backscatter. Still, the elevation and exposition of areas seemed to play a role in the study area, most probably for water bodies and permanent forest cover.

In the context of the features used for the classification and their contribution to the final result, it must also be discussed that the indices derived from combinations of both SAR sensors were potentially distorted by the relatively stronger level of backscatter from the L-band compared to the C-band. Accordingly, these indices might be dominated by the contribution of HH and HV (ALOS-2). In this study, the magnitude of input features probably diluted this effect, but the actual difference of both sensors should be investigated in further studies, probably raising the need for a proper mutual normalization of both sensors.

The values in

Table 5. The five most important features per class according to their permutation importance (PI) [%]. CI = combined index.

#	Vegetation	Open	Water	Camp
1	RFDI: 3.1	RFDI: 3.6	slope: 3.1	HH: 3.4
2	VSI: 2.2	VSI: 3.2	elevation: 1.8	RFDI: 2.7
3	elevation: 2.0	slope: 2.5	HH: 1.5	CI1: 2.4
4	CI2: 1.5	elevation: 2.2	StDev: 0.7	StDev: 1.8
5	RVI: 1.2	HH: 1.9	VH: 0.5	elevation: 0.9

show the permutation importance per class, namely how much worse (in percent of the achieved training accuracy) a classifier would perform if the named feature had random values. The table shows that the RFDI and the VSI were the most important variables to discriminate vegetation from open areas, but also that topographical features played a role here, as the forests are restricted to certain heights and are absent on steep slopes [55]. Also, as indicated above, the combined index CI2, potentially highlighting complex surfaces, such as crowns and other volumes, was among the five most important variables for the identification vegetation. The water class was mainly determined by areas with low elevation and small slopes, but also very low values in the HH polarization. Lastly, the camp was characterized by high HH values, but also by CI1, as was probably retrieved from the many horizontal structures. The high importance of the standard deviation among all four polarizations (StDev) again underlined the complex backscatter characteristics of settlements of light construction materials.

Aside from the thematic information content of the input data for the classification, all features were evaluated regarding the window size, which was used for their computing. This refers to the scale of the digital terrain analysis (Section 2.2.3) and the pixel size of the median filter applied to all SAR data (Section 2.2.4). Figure 9 shows that large window sizes were selected more frequently during the random forest classification compared with small ones. Of the 500 most important features as listed in Table 4, 124 were of the largest size of 49 pixels (corresponding to a spatial scale of 490 m). Generally, large window sizes have higher prediction values than smaller ones, most probably because they are more suitable to predict large consistent areas. This was already discovered by Braun et al. [48], who reported a dominance of large-scale features for land cover assessments around the refugee camp of Dagahaley in Kenya. As another possible explanation, products of smaller window sizes are often characterized by topographic variations or speckle effects, which do not represent the spatial distribution of the expected land cover. In our study, features with window sizes of five pixels slightly exceeded the linear trend. Considering the spatial resolution of the used raster, this corresponded to a size of 50 m, which apparently played a more important role in the study area. Consequently, it might be sufficient for future approaches to derive filter products of window sizes of 5, 25, and 49 pixels to keep computing time low, while results of similar prediction quality can be expected.

Lastly, the importance of specific features was evaluated by calculating the average PI and GD for each combination of feature and window size throughout all 11 investigated dates.

Table 6. Features with the highest average permutation importance throughout all classifications.

Feature and Window Size	PI [%]	Frequency
Elevation 49	21.1	3
Radar Forest Degradation Index (RFDI) 49	14.7	6
Slope 25	13.2	2
Aspect: northness 49	10.2	2
Slope 49	9.5	1
Volume Scattering Index (VSI) 49	8.3	9
Standard deviation (StDev) 49	7.1	8
Slope 15	7.0	3
HH 5	6.6	1
Radar Forest Degradation Index (RFDI) 25	6.5	3

and

Table 7. Features with the highest average Gini decrease throughout all classifications. GD = Gini decrease.

Feature and Size	GD Score	Frequency
Elevation 49	1808.1	3
Radar Forest Degradation Index (RFDI) 49	764.3	6
HH 15	585.0	4
Radar Forest Degradation Index (RFDI) 25	580.6	4
Standard deviation (StDev) 49	436.9	8
Radar Vegetation Index (RVI) 49	433.7	8
Aspect: northness 49	424.0	2
VSI Volume Scattering Index (VSI) 25	376.2	6
Combined index 2 (CI2) 9	296.1	3
Volume Scattering Index (VSI) 49	292.9	9

summarize the 10 combinations of features and window sizes with the highest PI and GD scores. Interestingly, elevation 49 had the highest average PI and GD scores, but only occurred in three of the 11 analyzed dates. This indicated a strong correlation between some of the automatically selected samples, as described in Section 2.2.2, with topographical features. Interestingly, the Volume Scattering Index (VSI) and the Radar Forest Degradation Index (RFDI) were among the most important features (Table 4) and showed generally large PI and GD scores (Table 6 and Table 7), but CI2 was not among these lists. Its role for vegetation is indicated in Table 5, however, attesting a high significance for the vegetation class. We conclude that VSI and RFDI may be the more important variables for the entire study area as a whole, as they also provide an explanation for the spatial distribution of the other classes (open, water, and camp). Furthermore, the information content of the combined indices (CI) was not directly representative of the scattering mechanisms, as was predictable by a canonical polarimetric matrix derived from fully-polarimetric sensors, as discussed in Section 2.2.4. The high PI scores of slope (sizes 25 and 49) could result from their high importance for the classification of water bodies, which only occur in flat areas. Of all original SAR backscatter features, the co-polarization HH of ALOS-2 had the highest frequency and the highest average variable importance. Furthermore, it was the only variable of this type among the 10 with the highest average scores based on window sizes of 5 and 15 pixels. This confirms the idea that even if most of the important variables were derived from large window sizes, a combination of features derived at different scales would lead to the most promising result. Despite the many indications on variable importance, we did not observe a prevalence of one sensor for the final result. Moreover, the combined measures (e.g., RFDI, VSI) were the ones that systematically achieved the highest scores. For the discrimination of forests from grasslands, the contribution of the sensor’s different wavelengths would have become more evident. This was already demonstrated by Braun et al. [129], who combined short and long radar signals of active and passive sensors to derive regional maps of above-ground biomass in Senegal.

By the time this article was written, the only satellites capturing radar images at the L-band were ALOS-2 and SAOCOM-1A, which are commercially distributed by JAXA and CONAE respectively [130]. This means that, currently, no L-band data is freely available (besides archived data of ALOS PALSAR between 2007 and 2011). This can be an obstacle for humanitarian organizations that do not have a budget to purchase satellite data. The ALOS-2 products used in this study cost ¥240,000 (roughly 2000 Euro) per scene [131]. Accordingly, time-series analyses such as this one can become expensive, and the transferability of this approach is limited in this respect. However, many organizations from the humanitarian domain have established agreements with operators and distributors so to get data at a discount, or even for free [132], or have increased their budget for digital technologies [133,134]. Lastly, new satellites will emerge in the upcoming years that will provide images acquired at the L-band (NISAR [2021], TanDEM-L [2023]), and even P-band (BIOMASS [2012]), entirely for free [130]. Against this background, it is important that methods already exist that address the special challenges of sensor combinations within humanitarian settings.

To conclude, the random forest classifier effectively handled the partially-redundant feature space and selected the ones which were most suitable towards predicting the land cover classes on the basis of the automatically derived training samples. This makes this approach largely operational, as it is no longer dependent upon an extensive collection of training samples by experts, nor from the manual selection of input features for the classification.

4. Conclusions and Outlook

This study showed how radar data of different frequencies can be integrated for land cover assessments and the monitoring of environmental conditions in situations of humanitarian crisis. It complements existing studies that, because of the limited availability of optical data, only roughly estimate the decrease of forests and the growth of camps. This study also introduces a high level of automation, as user input is only needed at the beginning for the selection of a few representative reference samples. As demonstrated in this study, a land cover classification with accuracies above 90% can be achieved by the selection of only 25 sample points per class. This number was proven to be a good trade-off between scientific accuracy and transparency, but also applicability and reproducibility. Applying radiometric indices on optical data to construct simple CARTS was used to identify all areas which fulfil similar criteria and therefore were suitable for training the classifier, even if the optical data were obscured by clouds, gaps, or did not fully cover the area of interest. Especially in cases of emergencies, this brings enormous benefits, as the workflow can be transparently transferred and automatically applied to other data sources or different areas. The only decision to be made by the users is the selection of target classes and the identification of a small amount of representative training areas. It is yet to be tested if other settings, for example, those based on different sensor combinations or located within other ecosystems, in achieving the same level of accuracy as observed in this study.

The use of microwaves additionally increased the quality of the land cover assessments, as they allowed a higher number of input images, especially in the rainy seasons. This helps to consider phenologic variations of the images and reduces their potential bias impacting the results as the images are analyzed as parts of a time series instead of a simple pre- and post-event change detection based on a single image pair. The benefit of time-series analyses over simple binary approaches has been demonstrated in numerous studies and should therefore also be integrated into monitoring frameworks of the humanitarian domain [7,47,135,136,137]. Furthermore, radar images can increase the quality of land cover classifications because they measure physical surface characteristics instead of spectral reflectance. This makes them especially valuable in areas of high vegetation cover where higher vegetation and grasslands have similar spectral characteristics. Among other features, this was accomplished by the computation of combined indices (CI) based on polarizations of both sensors. Although they were not based on the complex information of fully-polarimetric imagery, their role for the classification of vegetation and camp areas was demonstrated by their high feature importance—vegetation was found to be characterized best by the Volume Scattering Index (VSI) and the Radar Forest Degradation Index (RFDI), and the combined index 2 (CI2) seemed to be more sensitive towards the volumes and complex structures of crowns. In turn, the combined index 1 (CI) element was found to play a large role for the determination of camp areas where horizontal objects caused higher direct backscatter intensity (Table 5).

The study also showed that topography played a role in the identification of classes, especially in cases where large areas were covered by homogenous classes, such as water bodies and large evergreen forest areas. Topographic features turned out to be of largest information content when retrieved at large window sizes of 25 or 49 pixels, which correspond to distances of 250 to 500 m. It appeared that these were spatial scales that represented the connection between land cover and landforms best. This was underlined by the fact that 16% of the highest features were DEM-based, and were also slope and aspect features. However, these factors were constant throughout all analyzed dates and therefore only contributed to certain areas, while the seasonal variations were highlighted by changes in backscatter intensities and its derivations.

All classified dates showed high classification accuracies, especially for water and camp areas. A comparison with a reference study [23] showed that the camp areas might have been slightly under-estimated at the edges of the expansion areas where the dwelling density is low. This, as well as the high impact of seasonal variations, has to be considered when land cover maps are evaluated for the monitoring of camp expansion and forest degradation. Furthermore, the benefit of a higher amount of input images became evident here. Although the temporal difference between optical images used for the selection of training areas and the radar images used for classification was kept small, this delay potentially introduces errors in the training step initially, which significantly reduces the quality of the later classification. In an ideal case, optical and radar images are acquired at the same date. With the increasing availability of both radar (e.g., the upcoming Radarsat constellation mission [138] plus many of the planned SmallSAT missions [139]) and multi-spectral sensors (e.g., Planet Scope [140]), this will be more likely in the near future. Nonetheless, a minimum number of image pairs over the investigated period is required, as well as a budget to cover the acquisition costs of commercial satellite data.

Of course, the high classification accuracies are also a product of a low number of target classes. In this study, a discrimination of vegetation, open, camp, and water was sufficient, but in ecosystems with more complex land cover or different seasonal dynamics, a higher number of target classes (e.g., grassland, shrubs, but also solid urban areas vs. camp dwellings of light construction materials) should be considered. The limits regarding the number of classes and usable image pairs over a longer period are to be tested in further studies.

In a humanitarian context, the integration of SAR data not only gives higher flexibility regarding the number of images and the definition of suitable acquisition dates, but also allows for the identification of landcover types which are related to characteristics of surfaces (size, form, roughness, material, and scattering mechanisms). This makes them a valuable source of information for future emergencies, given that enough studies exist that demonstrate their applicability and transferability. More research is needed to investigate the benefits of radar data in humanitarian settings.

The presented insights on the development of vegetation patterns and the growth of the Kutupalong refugee camp are of higher detail than in all previously published studies. In particular, the incorporation of data from all seasons helps to understand the phenological variations and puts the overall developments into context. The suggested approach can be automated to a very large degree so that it can feed new data into monitoring frameworks, such as those conducted by the UNDP [56], to ensure the protection of the Teknaf Wildlife Sanctuary and the livelihood of the people living in this area.