Surface Soil Moisture Estimation from Time Series of RADARSAT Constellation Mission Compact Polarimetric Data for the Identification of Water-Saturated Areas

Abstract

Soil moisture is one of the main factors affecting microwave radar backscatter from the ground. While there are other factors that affect backscatter levels (for instance, surface roughness, vegetation, and incident angle), relative variations in soil moisture can be estimated using space-based, medium resolution, multi-temporal synthetic aperture radar (SAR). Understanding the distribution and identification of water-saturated areas using SAR soil moisture can be important for wetland mapping. The SAR soil moisture retrieval algorithm provides a relative assessment and requires calibration over wet and dry periods. In this work, relative soil moisture indicators are derived from a time series of the RADARSAT Constellation Mission (RCM) SAR compact polarimetric (CP) data over reclaimed areas of an oil sands mine in Alberta, Canada. An evaluation of the soil moisture product is performed using in situ measurements showing agreement from June to September. The surface scattering component of m-chi CP decomposition and the RL SAR products demonstrated a good agreement with the field data (low RMSE values and a perfect alignment with field-identified wetlands).

1. Introduction

Prior to the development of oil sands, the site of Suncor’s Base Mine is estimated to have been composed of over 50% wetlands, with a high proportion of these being swamps. Oil sands mining operators are required to reclaim their leases to a suitable landscape once they have finished mining. The reclamation process is challenging and it requires significant resources (and high costs for the operator). Wetlands are areas where a water table is at, near, or just above the surface and where soils are water-saturated for a sufficient length of time such that excess water and the resulting low soil oxygen levels are principal determinants of vegetation and soil development [1]. Identifying unplanned wetlands, here called opportunistic wetlands, using remote sensing leads to cost savings compared to conventional methods that may include artificially establishing wetlands.

The information about the spatial and temporal variations of global surface water extents, comprising all surface waters (open water, wetlands, or rice paddies) at 25 km resolution, was provided by the Global Inundation Estimate from Multiple Satellites (GIEMS) product, which relies on microwave radiometer data [2]. Monthly surface water fraction (SWF) maps at a 10 times higher spatial resolution were generated using Cyclone Global Navigation Satellite System (CYGNSS) data [3]. A time series analysis of Sentinel-1 SAR soil moisture products showed [4], through time, the wetting and drying of newly reclaimed landscapes. Surface soil moisture (SSM) products [4] were generated at a relatively high pixel spacing (ranging from 10 m to 100 m), enabling the characterization of water-saturated areas at the field scale. This method was used with the intent of identifying the ideal locations for creating swamps and other wetlands. Furthermore, it uncovers the opportunity to investigate the reclaimed landscape for opportunistic wetland development and delineation; that is, wet areas that are naturally suited to becoming wetlands and require only a small amount of additional reclamation activities or are forming naturally on the landscaping.

Deriving soil moisture from SAR using the Envisat Advanced Synthetic Aperture Radar (ASAR) multi-polarization data at a high spatial resolution was previously demonstrated [5] to be an efficient method suitable for watershed-scale hydrological studies. The retrieval of sub-kilometer relative SSM was shown to be valuable for assessing changes and trends in soil moisture across various spatial resolutions (as fine as 100 m) [6] extending the method to both the interfield and intrafield scales. High-resolution soil moisture data can be extracted using dual-polarimetric Sentinel-1 data [7] from across the vegetated soil surface. The results for soil moisture content obtained using an artificial neural network from generated using Sentinel-1 and ALOS-2 data demonstrated a high accuracy at 30 m resolution over low-vegetation soils (crop, grass, and shrub) for both SARs [8]. The accuracy of the ALOS-2 soil moisture content retrieval method over high-vegetation soils (forest) was higher than that of Sentinel-1. A combination of incoherent and coherent change detection method to retrieve SSM from Sentinel-1 data demonstrated a good accuracy over bare or sparsely vegetated soils [9].

The termination of the Sentinel-1B operation in December 2021 resulted in a significant reduction in the quantity of Sentinel-1 datasets in time series acquired within the 2022 season over the Base Mine area. This reduced number of time series became insufficient for achieving the optimal performance of the change detection-based algorithm [6,10,11,12,13,14,15] used in our previous study [4]. A significant contribution came from the RADARSAT Constellation Mission (RCM) SAR satellites, which acquire substantial amount of compact polarimetric (CP) datasets.

The literature on soil moisture estimation from full and/or compact SAR polarimetric data have shown the benefits of employing polarimetric decompositions to represent the SAR signal in terms of different scattering models [16,17,18,19,20,21,22]. A framework based on the polarimetric decomposition and copula quantile regression of RADARSAT-2 data improved the reliability of SSM estimation over corn-covered (i.e., vegetated) areas [23]. Five different decomposition approaches using various models for the individual scattering contributions have been investigated [22] using L-band full polarimetric data acquired with the airborne E-SAR system. The surface and dihedral components were used to retrieve SSM values for areas covered by agricultural vegetation. The study carried out by Hornacek et al. [16] utilized full polarimetric and simulated CP SAR data acquired by the airborne L-band uninhabited aerial vehicle synthetic aperture radar (UAVSAR) to estimate soil permittivity over croplands with vegetation cover. The method considered the extended Bragg scattering model for the soil surface utilizing the scattering-type parameter [24] (a function of the soil permittivity and the Barakat degree of polarization). This scattering-type parameter takes into account the depolarization information in relation to surface roughness thus enhancing the expected value of the inversion accuracies.

Several studies [17,25,26,27] investigated the effectiveness of C-band CP data for extracting soil moisture. A methodology based on the m-α CP decomposition [28] together with a surface component inversion was developed to retrieve SSM values for bare and vegetation-covered soils from RISAT-1 data [17]. Two surface scattering models, i.e., the integral equation method (IEM) and Extended Bragg, were reformulated for CP soil moisture inversion from the decomposed surface component. Another methodology to deliver accurate soil moisture products [27] was developed using the calibrated IEM multi-polarization inversion approach applied to RCM CH and CV simulated from RADARSAT-2 data. The potential of the RCM CP for soil moisture retrieval over bare soil was investigated by applying machine learning to 25 CP features [26]. An encouraging performance was achieved when using four CP features: (i) surface component of m-χ decomposition, (ii) coherency parameter (mu), phase difference (delta), and (iii) the Shannon Entropy (SE). The findings of a project on the evaluation and development of applications using RCM CP data [25] included an investigation of a theoretical approach based on IEM as well as the examination of effectiveness of simulated CP data for soil moisture estimation. The preferred CP configurations were suggested to employ circular Right–Left (RL) polarization (maximized backscattered intensity for surface scattering media) along with m-χ decomposition.

The primary goal of our paper is to apply the change detection method, previously successfully used for Sentinel-1 time series [4,13], for SSM estimation over the reclaimed areas at Suncor’s Base Mine from RCM CP data products, including RL polarization and m-χ decomposition. The paper describes the study area and data (Section 2), RCM CP data processing methods, and the SSM estimation results and provides discussion and conclusions (respectively, in Section 3, Section 4, Section 5 and Section 6).

2. Study Area and Data

2.1. Study Area

The study area (Figure 1), Suncor’s Base Mine, is located north of Fort McMurray in Alberta, Canada. The Base Mine includes areas of active excavation, soil stockpiling, water bodies, mining operations, extraction facilities, tailings management, and transportation infrastructure. Ground cover consists of a mixture of grasses and shrubs. In some areas, the vegetation density is light with minimal canopy cover. The recently reclaimed areas have ground cover which consists of saplings, shrubs, and grasses with sparse low-lying vegetation cover. Green dots distributed across the Base Mine (Figure 1) mark the locations of in situ meteorological stations that measure soil moisture.

2.2. RCM CP Data

The collection of a large quantity of RCM data was possible because there are three RCM satellites actively collecting data over the site. In total, 67 RCM CP datasets were available for the 2022 season, covering the period from 29 April to 31 October. This time period was used as the Alberta temperature records confirmed that there would be freezing temperatures up to 24 April and the presence of frost interferes with soil moisture estimation. The RCM CP data list included 48 datasets acquired in ScanSAR Medium Resolution 30 m mode (SC30MCP), which we used in our work. The RCM SC30MCP data have four-looks (2 × 2) and a pixel spacing of about 8 m x 23 m (range × azimuth). The swath width of SC30MCP is 125 km (accessible 350 km).

2.3. Electro-Optical Satellite Data

Very high resolution (VHR) electro-optical WorldView-2 and WorldView-3 scenes were acquired to cover the entire Base Mine four times during 2022 season. WorldView-2 has one panchromatic and eight multispectral bands. Its GSD (ground sample distance) at nadir is 46 cm for the panchromatic band and 1.85 m for multispectral imagery. WorldView-3 provides GSDs of 31 cm for the panchromatic band and 1.24 m for multispectral bands. WorldView-2/3 images were pansharpened with an output GSD of 50 cm. An example of a WorldView-2/3 image mosaic, acquired on 28 July 2022, is shown in Figure 1. VHR electro-optical satellite images were used to characterize variations in vegetation and generate water body masks (shown as a blue color in Figure 1).

2.4. Lidar Data

Additionally, a VHR LIDAR dataset acquired on 1 June 2022 was used to orthorectify SAR images and to calculate the local incidence angle (LIA) for each SAR image.

2.5. Field Data

Reference data for comparison with satellite-derived information comprised in situ observations of soil moisture from 18 meteorological stations (examples are shown in Figure 1). The amount of water in the soil (soil moisture) was monitored throughout the area of interest by measuring soil water content at various depths at the stations located throughout the mine area. Daily mean measurements [in %] from 1 May to 1 October 2022 were used for validation. Data from seven soil moisture monitoring stations were selected for comparison with the SAR SSM products. The meteorological stations recorded soil moisture at multiple depths. Only the 5 cm depth measurement was used for comparison. Daily precipitation data were also available from a weather station located at the Base Mine.

In 2017, an opportunistic wetland study was conducted to map the extent of wetlands [29]. The project involved field verifications of wetland delineations and classifications. These field-identified wetland locations can serve to verify water-saturated areas identified with higher values of satellite-derived SSM.

3. Methods

The RCM data processing workflow included the following seven steps:

(1)

Sigma naught (${\sigma }^0$) calibration;

(1)

Speckle filtering;

(2)

CP decomposition;

(3)

Geometric terrain correction and LIA calculation;

(4)

Water bodies masking;

(5)

Incidence angle normalization;

(6)

Relative SSM estimation using time series of CP products as input.

Steps 1–4 were implemented in the ESA SNAP software v.9.0 as a graph, enabling batch processing of all RCM CP datasets. Steps 5–6 were implemented using MATLAB R2023a allowing for flexibility in data analysis and visualization.

The RCM datasets acquired in SC30MCP mode were initially provided in multi-look complex format with CH and CV polarizations. The CP decompositions and certain other required methods are implemented in the ESA SNAP v.9.0 and CATALYST Professional v.2223 (former PCI Geomatica) software packages. There is also an option for users to order various RCM CP products (RL, RR, CP decompositions) generated with SAR Toolbox Version 3.4, which is implemented within the Earth Observation Data Management System (EODMS). In our analysis we use intensities of CH, CV, RL, and products of m-χ (m-chi) decomposition [30], i.e., surface (single-bounce (and Bragg) backscattering), double bounce, and volume (the randomly polarized constituent).

Considering that SC30MCP products are multi-looked (2 × 2), we investigated three different types of products regarding the window size of the averaging speckle filter: no filtering, 3 × 3, and 10 × 10. The 3 × 3 filter provides a better compromise between preserving spatial details and speckle noise reduction. Due to high topographical variations in the Base Mine area, the images of any point having elevation will have parallax in images acquired under different incidence angles [31]. This effect becomes even more significant when images are acquired from different (ascending and descending) orbits. Therefore, geometric terrain correction is critical to compensate for possible image distortions so that the geometric representation of the image will be as close as possible to the real world, especially for high resolution SSM estimation. The Range–Doppler terrain correction using the VHR LIDAR digital elevation model (DEM) was used to derive the precise geolocation information enabling co-registration (perfect alignment) of SAR images.

All RCM images were acquired under the nominal incidence angles in a range between 23.8– and 44.0 deg. Within this range, the CP calibration biases are negligible for SC30MCP mode. The LIA was calculated for each RCM scene. An example of LIA values ranging from 0 to 60 deg over an RCM scene fragment acquired on 2022–06-03 (under the nominal incidence angle of 26.5 deg) is shown in Figure 2 (top). In order to mitigate the dependence (Figure 2 (middle)) of ${\sigma }^0$ on LIA (${\theta }_{inc}$), the backscatter value ${\sigma }^o\left({\theta }_{inc},t\right)$ for each image (acquired at time $t$) can be normalized to a reference angle [13]. The normalization was performed for all CP data products but Figure 2 (bottom) shows only for surface component of m-chi decomposition.

According to RCM specifications, the nominal incidence angles for SC30MCP mode vary from 17.3 to 48.3 deg, depending on the beam mode. Therefore, for the reference incidence angle we selected the mean value of 32.8 deg. For an LIA (${\theta }_{inc}$) different from the reference angle, the normalized backscatter value (Figure 2, bottom) to an LIA of 32.8 deg was modeled using linear regression:

$${\sigma }^{\circ }\left({32.8}^{\circ },t\right)={\sigma }^{\circ }\left({\theta }_{inc},t\right)-{\beta }_r\left({\theta }_{inc}-{32.8}^{\circ }\right) \left[\mathrm{dB}\right],$$

(1)

where ${\beta }_r \left[d\mathrm{B}{/}^{\circ }\right]$ is the linear regression slope parameter, which characterizes the dependency of the backscatter on LIA. The linear regression slope parameter (${\beta }_r$) was determined as a best fit (in a least-squares sense) for the data. It was implemented using the MATLAB function “polyfit”. For the surface component of m-χ decomposition, ${\beta }_r=-0.0671$.

Then, the relative SSM for each RCM CP product can be found by application of the change detection method [4,11,12,13,32] with the following equation:

$$SSM\left(t\right)={\displaystyle \frac{{\sigma }^0\left({32.8}^{\circ },t\right)-{\sigma }_{\mathrm{dry}}^0\left({32.8}^{\circ }\right)}{{\sigma }_{\mathrm{wet}}^0\left({32.8}^{\circ }\right)-{\sigma }_{\mathrm{dry}}^0\left({32.8}^{\circ }\right)}},$$

(2)

where ${\sigma }^0\left({32.8}^{\circ },t\right)$ [dB]—backscatter coefficient at the reference angle and time t;

${\sigma }_{\mathrm{dry}}^0\left({32.8}^{\circ }\right)$ [dB]—backscatter coefficient observed under dry soil conditions (min value of ${\sigma }^0\left({32.8}^o,t\right)$ in SAR time series);

${\sigma }_{\mathrm{wet}}^0\left({32.8}^o\right)$—backscatter coefficient corresponding to wettest soil (maximum of ${\sigma }^0\left({32.8}^o,t\right)$ in SAR time series).

The relative SSM content ranges from zero in dry soil to unity (or 100%) in a completely saturated soil.

4. Results

4.1. Comparison of Satellite-Derived Estimates with In Situ Observations

The SAR SSM, calculated using Equation (2), is expressed as a relative value ranging from 0 to 1. In order to compare it with in situ soil moisture [in %], the SAR SSM can be rescaled to match its mean and variance to the reference data [13,33]. The comparison of different rescaling techniques [33] demonstrated that a simple rescaling procedure with a linear relation to match the range of variability in the corresponding soil moisture [34] is both effective and preferable. In our case, the linear rescaling SAR SSM values did not significantly change or expand the variance, indicating that the sensitivity of the SAR to soil moisture variations remained consistent. Figure 3 shows in situ soil moisture observed at one meteorological station, SSM derived from RL product, and precipitation measurements.

Table 1. RMSE values [in %] for SAR and in situ soil moisture.

CH	CV	m-χ Decomposition			RL
		Surface	Double Bounce	Volume
4.8	4.1	4.4	5.7	4.6	4.3

shows the Root Mean Square Error (RMSE) values for the evaluation of soil moisture between SAR-derived products and in situ measurements collected at the meteorological station. The lowest RMSE values for this station were achieved with CV, RL, and the surface component of m-χ decomposition. The scattering mechanisms remain relatively stable with soil moisture changes (i.e., χ remains constant), thus the change occurs primarily in the span and degree of polarization, where the three m-χ parameters are linearly linked. The difference in RMSE among the components of m-χ decomposition can be attributed to the signal-to-noise ratio of the intensities of those scattering mechanisms.

However, the challenge persists when attempting to compare in situ soil moisture at a station with spatially distributed SAR SSM products, especially in our case with a 30 m resolution. Therefore, it is important to compare the SAR products with the wetland areas identified in the field.

4.2. Comparison with Previously Reclaimed Wetlands

In order to obtain a single product from all the relative SSM results calculated for different dates, an average SSM can be generated.

An opportunistic wetland study was conducted to map wetland extents [29], involving field verifications of wetland delineations and classifications. Figure 4 presents all six RCM CP products, including CV, CH, m-χ decomposition (double bounce, volume, and surface components), and RL, allowing for a visual analysis of CP SAR SSM compared to field-delineated wetlands. By overlaying these CP SAR SSM products with field-delineated wetlands, the effectiveness of each product in capturing the spatial distribution and extent of water-saturated areas (with high soil moisture) can be evaluated. The example in Figure 4 shows that higher values of SSM derived from the surface component of the of m-χ decomposition and RL products have a good correspondence to the field-verified wetland locations, shown as black polygons. That means that utilizing these RCM CP products can help us identify water-saturated areas. The successful outcomes observed in surface scattering with RL polarization can be attributed to the phenomenon where a circular wave, when interacting with a surface, has its rotation direction inverted, following the backward scattering convention. The underestimation of SSM by the CH and CV-based products can be explained by their higher sensitivity to vegetation levels.

5. Discussion

Previous studies have demonstrated the suitability of multi-polarization data for SAR-based soil moisture monitoring across different scales. The literature highlights the benefits of using polarimetric decompositions and different scattering models to enhance soil moisture retrieval accuracy. The focus of our work was on the exploration and evaluation of different RCM CP data and their decomposition products for SSM retrieval using the change detection method [6,10,11,12,13,14,15] previously demonstrated with Sentinel-1 for the same application [4]. A comparison of the results for the SAR SSM retrieved from RCM CP data products with Sentinel-1 SSM [35] demonstrated that the m-χ surface component was in better agreement with field data. The important finding was the demonstration of the suitability of RL polarization along with surface component of m-χ decomposition for SSM retrieval.

The objective of our study was not to extract absolute soil moisture. Instead, it aimed to explore the effectiveness of using time series of SAR imagery to identify trends in relative soil moisture, as indicated by the soil moisture index. A persistent challenge in our and other studies on soil moisture estimation is the comparison of in situ measurements (using a contact-based sensor) with spatially distributed SAR SSM products, particularly given the 30 m resolution of the SAR data. Such a spatial resolution can complicate direct comparisons due to the heterogeneous nature of soil moisture distribution over larger areas. Consequently, it was an important finding that the SAR-derived products complemented the field-identified wetland areas to facilitate a better validation of soil moisture estimates. The correspondence between higher values of SSM, derived from the surface component of m-χ decomposition and RL products, and field-verified wetland locations meant that we could successfully identify water-saturated areas. The usage of these products highlights the effectiveness of RCM CP data in environmental monitoring involving soil moisture. This SAR-based technique is especially suitable for non-vegetated or sparsely vegetated flat areas which are typical for newly reclaimed oil sands mining landscapes. Compared to traditional field-based methods, the presented approach is more efficient and cost-effective.

Future studies will include verification with SSM measured at different depths across multiple stations, taking into account major contributing factors such as SAR speckle noise, surface roughness, vegetation, topography, and soil characteristics. Analysis of soil moisture within different soil layers will provide a possibility for better understanding of soil moisture dynamics and impacts of different variables. Further studies can also investigate combining/fusion (e.g., weighted sum) of the CP SAR SSM products to identify the best approach for addressing some of these contributing factors (surface roughness, vegetation, and soil characteristics) and tackling the challenges associated with polarimetric SAR sensitivity to soil moisture. Such fusion could lead to a more robust method for interpreting CP SAR data, improving the accuracy of soil moisture estimation for various environmental conditions.

The penetration capability of L-band signals is much deeper, providing soil moisture information for deeper layers, than the C-band signals. For this reason, the use of L-band data for soil moisture retrieval will be explored with available SAR missions (e.g., upcoming NISAR). Additionally, a comparison of C-band and L-band SAR SSM products will allow for a more comprehensive evaluation of soil moisture retrieval accuracy and its dependence on the landcover types, soil depth, and different environmental conditions.

6. Conclusions

The study confirmed that the surface component of the m-χ decomposition, CV, and RL products derived from RCM CP data demonstrate a good performance in SSM estimation, as validated against field measurements (at 5 cm depth) obtained from a meteorological station. To tackle the challenge of comparing single-point measurements at a station with spatially distributed RCM SSM products with a 30 m resolution, we utilized the field-verified wetlands to identify areas with higher SSM. In this scenario, the higher values of SSM derived from the surface component the of m-χ decomposition and RL products show better correspondence than CH and CV-based products, which underestimate SSM due to their higher sensitivity to vegetation levels.

A wide range of environmental conditions, characterized by variations in precipitation and temperature, were considered for validation during the snow-free season from May to October. Additionally, terrain properties were taken into account using various sources of information, including VHR lidar DEM, vegetation cover, and the proximity to water bodies, allowing for a comprehensive validation process for the diverse environmental and terrain conditions. Future work will extend this validation framework to include more variables, such as soil characteristics.

The novelty of our work was in comparing and demonstrating different RCM CP products for SSM retrieval and in demonstrating the suitability of RL polarization. Our SSM products can be used as an input in machine learning models for wetland classification.