Daily High-Resolution Land Surface Freeze/Thaw Detection Using Sentinel-1 and AMSR2 Data

Abstract

High-resolution surface freeze/thaw (F/T) information is valuable for hydrological, frost creep and gelifluction/solifluction, and climate prediction studies. Currently, large-scale, high-resolution F/T detection is restricted by low spatial resolution of passive microwave remote sensing sensors or low temporal resolution of synthetic aperture radar (SAR) data. In this study, we propose a new method for detecting daily land surface F/T state at 1 km spatial resolution by combining the Sentinel-1 radar and the Advanced Microwave Scanning Radiometer 2 (AMSR2) with leaf area index (LAI) data. A non-linear relationship is established between the 1 km F/T index from Sentinel-1 with 1 km F/T index from AMSR2 (FTI) and 1 km LAI data. The 1 km FTI is a disaggregation of the 25 km FTI obtained from AMSR2. This non-linear relationship is then applied to daily 1 km FTI and LAI data to predict the 1 km daily F/T index, based on which the F/T status is detected with grid-cell-based F/T thresholds. The overall accuracy of this daily 1 km F/T is more than 88.1% when evaluated with the in situ 5 cm soil temperature over China and Canada. This study is valuable for detecting daily, high-resolution F/T status and is helpful for studies related to disaster and climate prediction.

1. Introduction

Approximately 50 million square kilometers of the land surface experiences soil freezing and thawing [1]. Part of this surface exhibits spatially discontinuous soil freezing and thawing changes due to the heterogeneity of topography, land cover, subsurface properties, snow cover, and other factors. The spatial and temporal change information of land surface freeze/thaw (F/T) is essential for hydrological, climatic, topographic deformation, and ecosystem processes [2,3,4,5,6]. The land surface F/T that occurs in permafrost regions increases emissions of carbon and nitrogen [7,8,9,10] and leads to uncertainties in climate change predictions [11]. Permafrost degradation in high mountains could lead to disasters, such as rock and ice avalanches [12]. High-resolution mapping of permafrost covering high-altitude areas is important for infrastructure planning [13]. High-resolution monitoring of surface F/T states is important for various research areas, particularly disasters, such as frost creep and gelifluction/solifluction [14], spring floods, and climate prediction [15,16]. Remotely sensed satellite data are efficient for detecting land surface F/T states at global and regional scales; satellite microwave remote sensing at long wavelengths is well suited for detecting the soil F/T state because of the ability to penetrate into the soil medium [17,18] and the high contrast in permittivity between liquid water and ice [19]. Active and passive microwave data have been widely used to detect land surface F/T in recent years. Several global and regional microwave surface F/T products have been produced with active and passive microwave sensors, such as the Advanced Microwave Scanning Radiometer 2 (AMSR2) [20,21], Soil Moisture and Ocean Salinity (SMOS) [22], Soil Moisture Active Passive (SMAP) [23], the Aquarius L-band radiometer [24], and the advanced scatterometer (ASCAT) [25]. All these F/T products detect the land surface F/T state at spatial resolutions from 25 km to 36 km, which cannot accurately describe the surface F/T process within heterogeneous pixels and does not satisfy the requirements of many applications, e.g., frost creep and gelifluction/solifluction and flood disaster monitoring at regional scales. The scattering properties of synthetic aperture radar (SAR) data, such as scattering type (surface, double-bounce, and volume scattering) and entropy, in different land covers (such as forest or non-forest, snow-covered or snow-free), land F/T status had significant seasonal changes and could be used for mapping land cover types and surficial geocryological characteristics of the permafrost active layer [26,27]. The obviously seasonal variations of the backscattered signal has been used to detect land surface F/T. Sentinel-1 (S1) satellite data, which can provide free high-resolution SAR data since 2015 and have high data volume coverage over study areas in this study, have been successfully used to detect land surface F/T [28,29,30,31,32]. However, the temporal resolution of Sentinel-1 depends on the geographical location. For example, Sentinel-1 is available daily in northern Europe but only once every 12 days in northeastern China. A high-resolution F/T detection algorithm is required to detect daily high-resolution surface F/T states over some areas where S1 data are not available daily and meet the demand for many regional frost creep and gelifluction/solifluction and climate change studies.

Some studies have developed F/T detection algorithm in which high-resolution optical/thermal data were used to downscale passive microwave data. Kou et al. [33] detected the land surface F/T state at 0.05° with the discriminant function algorithm (DFA) [34] from brightness temperature (Tb), which was downscaled by blending the land temperature data using a Bayesian maximum entropy method. Zhao et al. [35] established a relationship between the passive microwave F/T index (FTI) and thermal infrared observations to monitor high spatial resolution surface F/T states at 0.05°. However, these studies mainly concentrated on the Tibetan Plateau area by combining passive microwave and thermal infrared data. Additionally, using these methods might reduce the accuracy of passive microwave soil F/T detection [35] due to the different depths of F/T represented by passive microwave and thermal infrared.

Currently, there are a few studies regarding the detection of the F/T state that combine active and passive microwave satellite data. Zhong et al. [36] estimated more accurately F/T onsets by using a machine learning approach of combination SMAP and ASCAT. This study revealed the possibility of detecting F/T state by combing active and passive microwave remote sensing. The combination of passive and active data provides a potential means for downscaling soil moisture estimates to a finer spatial resolution [37,38,39]. The soil dielectric constant is highly dependent on the liquid water content; the higher the amount of free liquid water is, the higher the dielectric constant. For low dielectric constant values (dry or frozen soil), the backscattered signal from the soil surface is at its minimum, and the soil surface emission is at its maximum. Many methods, such as the Tb-based downscaling method [39], soil moisture-based downscaling method [37], and change detection method [40], have successfully estimated soil moisture at high spatial resolution by combining active and passive data. Surface F/T detection relies on observing changes in soil dielectric properties rather than the absolute value.

In this work, we demonstrate the feasibility of high spatial resolution surface F/T state determination by a combination of active and passive microwave data over four test regions located in China and Canada. The objectives of the study are to detect the daily land surface F/T state at a 1 km spatial resolution and to propose a new approach for detecting the F/T state combining radar and passive remote sensing data. The advanced integral equation method (AIEM) microwave emission and scattering model incorporating a dielectric constant of F/T soil is used to study the behavior of the radar and passive signals from thawed to frozen. Finally, a new approach is demonstrated for the detection of the surface F/T state at high spatial and temporal resolutions using Sentinel-1 and AMSR2 data by incorporating leaf area index (LAI) data. In Section 2, the study area and data used in this study are briefly described. In Section 3, a detailed description of our new algorithm development is presented. In Section 4 and Section 5, the results and discussion are given, followed by the conclusion in Section 6.

2. Materials

The dataset used consists of in situ data, Sentinel-1 SAR data, AMSR2 Tb data and auxiliary data. The auxiliary data include the ERA5-land reanalysis dataset and LAI data. The in situ soil moisture and temperature network located in China and Canada (Figure 1) was used to validate the accuracies of the F/T detection results.

2.1. In Situ Data

For the purposes of this study, data from four areas with distributed near-surface soil moisture and a temperature observation network could be obtained. Three in situ soil moisture and temperature observation networks, including the Genhe watershed, Saihanba, and Naqu, are located in China (Figure 1c–e), and another in situ soil moisture and temperature observation network, i.e., Real-time In Situ Monitoring for Agriculture (RISMA) network, located in Canada (Figure 1g). Each study area is affected by seasonal soil F/T and located within one Sentinel-1 image.

Genhe watershed study area is located in the northeast of China (50.25~51.25°N, 120~122°E). The prominent land cover types in the Genhe watershed, which has a combination of a cold and humid temperate climate and a continental monsoon climate, are forest, grassland, and cultivated land. The contents of clay and sand in soil over Genhe watershed are 36–52% and 3–29%, respectively, and the elevation ranges from 500 to 1420 m. The Genhe watershed is widely covered with permafrost. A soil temperature and moisture observation network has been established over the Genhe watershed since October 2013 [41,42]. Saihanba study area is located in the north of China (41.75~42.75°N, 115.5~117.75°E). Saihanba, with semi-arid and semi-humid climates, is mainly covered with forest and grassland with elevations ranging from 700 to 2000 m. The contents of clay and sand in soil over Saihanba are 9% and 79%, respectively. Seasonal frozen soil is widely distributed here. A soil moisture and temperature observation network has been established to measure data at depths of 5 cm and 10 cm since August 2018 over Saihanba [41]. Naqu, with a semiarid climate, is located in the central Qinghai-Tibetan Plateau and is distributed with permafrost (31.25~32°N, 91.5~92.5°E). The contents of clay and sand in soil over Naqu are 20% and 66%, respectively. With elevations ranging from 4280 to 5900, Naqu is mainly covered with alpine meadows. The in situ soil temperature at a depth of 5 cm was obtained from a multiscale soil moisture and temperature monitoring network [43]. The RISMA soil moisture and temperature network used in this study is located in Manitoba, Canada, covered with cultivated land and forest [44]. The contents of clay and sand in soil over RISMA are 9–70% and 3–90%, respectively. This RISMA network was established since 2011. The soil sensors measured the real dielectric permittivity, soil moisture, and soil temperature using Stevens Hydra Probe sensors at surface (0–5 cm), 5 cm, 20 cm, 50 cm, and 100 cm depths. Three hydra probe sensors were installed at each depth. Soil temperature data at a depth of 5 cm with an hour interval are downloaded from https://ismn.geo.tuwien.ac.at/en/and used in this study. RISMA has a cold climate and elevation ranging from 130 to 720 m. The RISMA study area (49~50.25°N, −100~−97°E) in this study includes the Manitoba RISMA network.

Taking into account the available time periods of the various soil moisture and temperature observation networks and the auxiliary data, a total of 19, 16, 23, and 13 ground sites from the Genhe watershed from June 2017 to May 2019, from Saihanba from August 2018 to December 2019, and from Naqu and RISMA from January 2018 to December 2019 were used in this study. Short descriptions of the in situ sites from the four study areas are given in

Table 1. Coordinates and land cover type of the in situ sites.

Area	Site Name	Lat (°)	Long (°)	Land Cover	Site Name	Lat (°)	Long (°)	Land Cover
Genhe	Site 1	50.507	120.529	Grasslands	Site 17	50.451	120.987	Forest
	Site 2	50.451	120.711	Forest	Site 18	50.327	120.484	Grasslands
	Site 3	50.448	120.834	Grasslands	Site 19	50.329	120.696	Cultivated
	Site 5	50.413	120.547	Grasslands	Site 20	50.311	120.589	Cultivated
	Site 9	50.556	120.955	Grasslands	Site 24	50.309	120.927	Cultivated
	Site 11	50.301	120.836	Cultivated	Site 26	50.256	120.948	Cultivated
	Site 12	50.367	120.883	Cultivated	Site 27	50.529	120.499	Grasslands
	Site 14	50.511	120.581	Grasslands	Site 28	50.463	120.537	Grasslands
	Site 15	50.575	120.843	Forest	Site 29	50.341	120.977	Forest
	Site 16	50.492	120.926	Forest
Saihanba	A3	42.312	117.242	Grasslands	P8	42.311	117.233	Forest
	A5	42.309	117.236	Grasslands	P9	42.249	117.294	Grasslands
	A6	42.308	117.241	Grasslands	P10	42.255	117.359	Forest
	A7	42.305	117.231	Grasslands	P11	42.201	117.199	Forest
	A11	42.307	117.233	Forest	P12	42.236	117.236	Grasslands
	P2	42.351	117.207	Forest	P13	42.164	117.302	Forest
	P6	42.367	117.296	Forest	P15	42.135	117.242	Forest
	P7	42.261	117.131	Forest	P16	42.149	117.371	Forest
Naqu	BC07	31.274	92.109	Grasslands	MS3513	31.677	91.842	Grasslands
	BC08	31.332	92.041	Grasslands	MS3518	31.661	91.794	Grasslands
	C1	31.683	91.771	Grasslands	MS3523	31.639	91.754	Grasslands
	C3	31.614	91.774	Grasslands	MS3527	31.614	91.739	Grasslands
	C4	31.618	91.841	Grasslands	MS3533	31.586	91.793	Grasslands
	CD01	31.712	92.458	Grasslands	MS3545	31.573	91.912	Grasslands
	CD07	31.495	92.132	Grasslands	MS3603	31.259	91.799	Grasslands
	F4	31.698	91.773	Grasslands	MSNQRW	31.463	92.017	Grasslands
	F5	31.693	91.786	Grasslands	P1	31.782	91.729	Grasslands
	MS3475	31.946	91.721	Grasslands	P10	31.807	91.845	Grasslands
	MS3494	31.805	91.749	Grasslands	P11	31.815	91.795	Grasslands
	MS3501	31.754	91.782	Grasslands
RISMA	MB1	49.562	−98.019	Cultivated	MB8	49.752	−97.982	Cultivated
	MB2	49.492	−97.933	Cultivated	MB9	49.694	−98.024	Cultivated
	MB3	49.519	−97.956	Cultivated	MB10	49.975	−97.348	Cultivated
	MB4	49.636	−97.988	Cultivated	MB11	50.111	−97.573	Cultivated
	MB5	49.621	−97.957	Cultivated	MB12	50.189	−97.598	Cultivated
	MB6	49.678	−97.959	Cultivated	MB13	49.932	−99.387	Cultivated
	MB7	49.665	−98.007	Cultivated

. Figure 1c–e, and Figure 1g show the land cover map and distribution of in situ stations of the Genhe watershed, Naqu, Saihanba, and RISMA. The GlobeLand30–2010 land cover map provided by the National Geomatics Center of China was used here, and the land cover map was produced by integrating pixel- and object-based methods with knowledge and over 10,000 Landsat satellite images with more than 80% overall accuracy [45] (http://www.webmap.cn/mapDataAction.do?method=globalLandCover, accessed on 5 December 2021)

2.2. Satellite Data

Sentinel-1 comprises a constellation of two polar-orbiting satellites (Sentinel-1A (S1A) and -1B (S1B)), performing the C-band (5.406 GHz) synthetic aperture radar (SAR) imaging. S1A and S1B were launched in April 2014 and April 2016, respectively. Sentinel-1 data products are available from https://scihub.copernicus.eu/dhus/#/home, accessed on 20 December 2021. Only S1B has acquired data over the Genhe watershed and RISMA, and only S1A has acquired data over Naqu. Both S1A and S1B acquired data over Saihanba, but to match the descending AMSR2 data with overpass time at 1:30 am, S1B data with overpass time at 22:20 pm were used. The revisit frequency of S1 data over the four study areas is as long as 12 days. A short description of the S1 data of the four study areas is given in

Table 2. Detailed information on S1 data over the four study areas. Number means the number of S1 SAR images during the time period.

Area	S1A/S1B	Time Period	Number	Overpass Time
Genhe watershed	S1B	June 2017–May 2019	55	22:10
Saihanba	S1B	August 2018–December 2019	33	22:20
Naqu	S1A	January 2018–December 2019	59	23:50
RISMA	S1B	January 2018–December 2019	57	00:23

.

The S1 SAR data used in this study are interferometric wide swath mode (IW) images in both VV and VH polarizations that were generated from the Level-1 ground range detected (GRD) product. First, the S1 SAR images were preprocessed. The preprocessing steps include multilooking processing, radiometric correction, filtering, and geometric correction using Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data. Then, the image pixel values were converted to backscatter coefficients (${\sigma }^0$). The backscatter coefficients obtained from different incidence angles $\theta$ (${\sigma }_{\theta }^0$) were normalized to 40° (${\sigma }_{40°}^0$) according to the square cosine correction equation [46,47]:

$${\sigma }_{40°}^0=\frac{{\sigma }_{\theta }^0{\cos }^2\left(40°\right)}{{\cos }^2\left(\theta \right)}$$

(1)

Finally, the spatial resolution of the preprocessed S1 data was resampled from 20 m to 1 km using the average sampling method.

The passive microwave data used in this study are AMSR2 L3 Tb (6.925 GHz and 36.5 GHz) with a spatial resolution of 0.25° provided by the Japan Aerospace Exploration Agency (JAXA) (https://gportal.jaxa.jp/gpr/information/product, accessed on 20 December 2021). The AMSR2 Tb data were selected from descending orbits in the same time period as the S1 data over different areas. Then, the spatial resolution of the Tb data was resampled to 1 km with a simple weighted method.

2.3. Auxiliary Data

ERA5-land reanalysis data from https://cds.climate.copernicus.eu/were used in this work for both the development and validation of F/T estimates. The ERA5-land dataset provides 0.1° hourly air (~2 m height), soil (layer 1: 0~7 cm) and hourly snow depth data in gridded format.

The Global Land Surface Satellite (GLASS) LAI product, with a spatial resolution of 500 m and a temporal resolution of 8 days, was generated from time series MODIS surface reflectance data using general regression neural networks [48,49]. To match the resolution of S1B data, the GLASS LAI data were resampled to a 1 km spatial resolution using the average sampling method. Considering the small change in vegetation, the LAI within 8 days was treated as the same value.

3. Methods

3.1. F/T Indices Derived from Active and Passive Microwave Observations

3.1.1. F/T Indices from Radar Observations

Algorithms for F/T detection from SAR are typically based on the threshold detection approach [28,29,30] and may include the additional use of models [31,32]. Threshold detection algorithms are based on examining the temporal change in the backscatter coefficient relative to reference thresholds acquired during frozen and thawed seasons. Threshold algorithms are often called seasonal threshold algorithms (STAs). Derksen et al. [23] identified the land surface F/T state from SMAP satellite radar data by examining the temporal changes in the seasonal scale index (SSI). SSI is a linear stretching of the time series radar backscatter and is linearly related to radar backscatter with an approximate range of 0~1. The SSI obtains a value of 0 when the backscatter is the same as the frozen reference and a value of 1 when the backscatter signal is the same as the thaw reference. In this study, we first implemented the STA for S1 data to detect the surface F/T state. The original SSI in Derksen et al. [23] is defined as follows:

$$\mathrm{SSI}=\left({\mathsf{\sigma }}_t-{\mathsf{\sigma }}_{fre}\right)/\left({\mathsf{\sigma }}_{thaw}-{\mathsf{\sigma }}_{fre}\right)$$

(2)

where ${\mathsf{\sigma }}_t$ is the backscatter coefficient at time t, ${\mathsf{\sigma }}_{fre}$ is the average of 10 minimum backscatter values obtained from winter (December, January and February), and ${\mathsf{\sigma }}_{thaw}$ is the average of 10 maximum backscatter values obtained from summer (June, July, and August). ${\mathsf{\sigma }}_{fre}$ and ${\mathsf{\sigma }}_{thaw}$ can be approximated as the minimum and maximum backscatter, respectively. Finally, the land surface F/T states are categorized by the following threshold value:

$$\mathrm{SSI}>\mathrm{T}, \mathrm{thaw}$$

(3)

$$\mathrm{SSI}<\mathrm{T}, \mathrm{frozen}$$

(4)

where T is the threshold. In previous work, T was set to 0.5 [23] or was optimized using in situ soil/air temperature data observed from measurement stations [29,30]. However, the feasibility of T with a fixed value or optimized value is limited by certain regions. In this study, we defined the pixel-based T with the distribution fit of SSI and ERA5-land air temperature, which is described in detail in Section 3.2.3. Figure 2 and Figure 3 show the time series of the boxplot of SSI calculated with S1 data (noted as SSI_S1) at VH and VV polarization over four study areas, respectively. The green line is the median of SSI_S1. In Figure 2 and Figure 3, the land surface is easy to classify as a thawed or frozen state when the SSI_S1 is close to 1 or 0, respectively. To compare the F/T detecting ability of SSI_S1 at VH and VV polarization, we calculated the difference in SSI_S1 between the minimum value in summer and the maximum value in winter. The larger the difference in SSI_S1 between summer and winter, the more sensitive it is to land surface F/T. In Figure 2, the difference in SSI_S1 at VH polarization is 0.83, 0.7, 0.4, and 0.74 over the Genhe watershed, Saihanba, Naqu, and RISMA, respectively. In Figure 3, the difference in SSI_S1 at VV polarization is 0.72, 0.55, 0.63, and 0.6 over the Genhe watershed, Saihanba, Naqu, and RISMA, respectively. Obviously, SSI_S1 at VH polarization is more suitable for F/T detection over the Genhe watershed, Saihanba, and RISMA, and SSI_S1 at VV polarization is more suitable for F/T detection over Naqu. In previous studies, Baghdadi et al. [29], Fayad et al. [30], and Rodionova [50] found that VH polarization is more suitable for F/T detection than VV polarization over different study areas, although VV polarization is more sensitive to the soil dielectric constant. We found that the land surface of the above three studies is mainly covered with forest, crop, shrub, and grassland, which is similar to the Genhe watershed, Saihanba and RISMA areas. The sensitivity of VH polarization to volume scattering from vegetation may cause a greater difference in SSI_S1 between summer and winter. Over the Naqu area, the effect of volume scattering from meadows is almost negligible; therefore, SSI_S1 at VV polarization is more sensitive to land surface F/T. Therefore, SSI_S1 at VH polarization is selected to monitor land surface F/T over the Genhe watershed, Saihanba and RISMA areas, and SSI_S1 at VV polarization is selected to classify land surface F/T over the Naqu area.

3.1.2. F/T Indices from Passive Microwave Observations

DFA [21,34], which is based on AMSR2 data, was chosen to detect land surface F/T states for passive microwave remote sensing data. Both the SMAP and AMSR2 F/T detection algorithms have F/T indices to identify surface F/T, while the F/T detection index of SMAP (i.e., normalized polarization ratio (NPR)) is sensitive to the change in soil liquid water content and that of AMSR2 (i.e., FTI) is sensitive to changes in soil temperature [51]. In this study, the passive microwave F/T detection index was mainly used to provide surface temperature information. DFA examines the 2-D spatial distribution of V polarization Tb of 36.5 GHz (Tb36.5 V) and the quasi-emissivity (Qe) from AMSR2 satellite data relative to signatures acquired during frozen and thawed states. These two factors can distinguish F/T by Fisher discriminant analysis. Tb36.5 V was selected because it is the most sensitive band to land surface temperature [52]. As the ratio of Tb was between any lower frequency and 36.5 V, Qe was mainly influenced by changes in the dielectric constant. Several studies [21,34,53] have used DFA to successfully monitor F/T states. To obtain the DFA over the Genhe watershed, based on the in situ 5 cm soil temperature data and AMSR2 C and Ka-band data, Fisher linear discriminant analysis was used to parameterize the DFA. The C band was chosen because of its similar wavelength to the Sentinel-1B. The FTI was computed as follows:

$$\mathrm{FTI}=a\times \mathrm{Tb}36.5\mathrm{V}+b\times \left(\mathrm{Tb}6.925\mathrm{H}/\mathrm{Tb}36.5\mathrm{V}\right)+\mathrm{c}$$

(5)

where FTI corresponds to the passive F/T detection index, and Tb6.925H corresponds to a Tb of 6.925 GHz at H-polarization, as H-polarization is more sensitive to soil liquid water content. The coefficients in Equation (5) are parameterized by areas with in situ 5 cm soil temperature and vary with study area and shown in

Table 3. The coefficients in Equation (5) over four study areas.

Area	a	b	c
Genhe watershed	−0.4552	−44.3366	160.5139
Saihanba	−0.3668	−35.4591	129.0621
Naqu	−0.279	40.1433	35.5872
RISMA	−0.2266	−3.3469	60.7911

. The coefficients of Genhe watershed and Saihanba are obviously different from Naqu and RISMA. Genhe watershed and Saihanba have similar location and closer coefficients. Finally, the F/T state is classified according to the following formula, where T is set to 0:

$$\mathrm{FTI} < \mathrm{T}, \mathrm{thaw}$$

(6)

$$\mathrm{FTI} > \mathrm{T}, \mathrm{frozen}$$

(7)

3.2. Development of the High-Resolution F/T Detection Algorithm

The FTI and SSI have obvious F/T season characteristics, as shown in Figure 4, and can be used to detect the surface F/T state. The seasonal variation in the FTI was consistent with the seasonal changes in soil temperature. Over bare land surface, the SSI changes during the F/T transition are mainly influenced by the dielectric constant. Typically, when the soil freezing process starts, the SSI decreases rapidly due to the drastic decrease in free liquid water within the top layer of the soil surface. As the soil freezing front progresses in the deeper layers and the soil temperature decreases, the SSI decrease slows down and tends to stabilize. Similarly, during the soil thawing process, the increase in SSI is typically rapid at the beginning and slows down as the thawing process continues. Both FTI and SSI have seasonal responses to land surface F/T, thus it is potential to detect surface F/T by analyzing the relationship between the FTI and SSI. The relationship between the dielectric constant and temperature typically follows a logistic curve during the soil surface F/T process [54,55,56]. In our study area, the soil dielectric constant is mainly determined by the soil water content and its phase. A near-linear relationship between soil liquid-water content and radar backscatter has been demonstrated and implemented in high-resolution soil moisture estimation studies [37,40,57]. As SSI is linearly related to backscatter, SSI is linearly related to the dielectric constant. Zhao et al. [35] developed a high-resolution near-surface F/T detection algorithm based on the linear relationship between FTI and thermal infrared temperature with an accuracy higher than 70% over the Qinghai–Tibet Plateau region. In a study validating the performance of the DFA parameterized by Wang et al. [21] in China, Wang et al. [51] also found that the FTI was linearly correlated with surface temperature. Naeimi et al. [25] found that the behavior of backscatter concerning soil temperature followed a logistic curve in most cases based on the ASCAT data, especially in high latitude areas with sparse vegetation. Without considering the influence of vegetation on backscatter, the behavior of the SSI from the soil (noted as SSI_soil) and the FTI is expected to follow a logistic curve, similar to the relationship between the dielectric constant and soil surface temperature. The relationship between SSI_soil and FTI is expressed as follows:

$${\mathrm{SSI}}_{\mathrm{soil}}=\frac{1}{1+a\times \exp \left(-b\times \mathrm{FTI}\right)}$$

(8)

where a and b are the model parameters.

3.2.1. AIEM Model Simulations of FTI and SSI

An analysis of the potential of combining S1 C-band data and AMSR2 data to detect land surface F/T states was carried out using simulations of Tb and backscatter. The radar backscattering coefficients and passive Tb of bare soils under frozen and unfrozen conditions were simulated using the AIEM [58], incorporated with the dielectric constant model of frozen and unfrozen soil (Zhang’s model) [55,59,60,61]. One of the factors used to determine microwave emissions from a rough surface is the dielectric constant of soil. Dobson et al. [59] developed a dielectric mixing model for thawed soil that has been widely applied in the field of microwave remote sensing. An extension of the semi-empirical dielectric constant model was developed by adding an ice item to respond to frozen conditions [55,60]. These two dielectric mixing models were used to simulate the dielectric constant of F/T soil. Mironov et al. [54] validated the Zhang’s model [55,60] with dielectric data for soils and showed that Zhang’s model could simulate the dielectric constant of frozen soil well, with root mean square errors (RMSEs) of 2.16 and 1.18 for the real and imaginary parts, respectively. AIEM was implemented to simulate the radar backscattering coefficient and Tb from the soil surface. As an improvement to the integral equation method (IEM), the AIEM reserves complementary components for the scattered fields and removes the simplifying assumption in the spectral representation of Green’s function [62]. The AIEM was used because it provides a better performance than the IEM [58]. The AIEM can better explain the surface microwave scattering and accurately simulate the ground radiometer data at AMSR-E frequencies [63,64]. The effects of frozen soil volume scattering on higher frequencies (Ku and Ka-band) rather than the C band have been observed and modeled in a few studies [65,66]. Although AIEM does not consider the volume scattering of frozen soil, it has been shown to successfully simulate surface emissions from frozen soil [34,53].

For simulations, the F/T soil dielectric constant was estimated by the dielectric mixing model. The dielectric constant estimates were set as input data for the AIEM to simulate the backscattering coefficients in VV and VH polarizations at an S1 frequency of 5.406 GHz and to calculate the passive microwave Tb at 6.925 GHz and 36.5 GHz. The input parameters for simulation of the radar backscattering coefficient and passive microwave Tb were randomly generated and are shown in

Table 4. Input parameters in the simulation of the radar backscattering coefficient and passive microwave Tb. RMS height indicates root mean square height. CL means correlation length.

Simulation	Frequency (GHz)	Theta (°)	Soil Moisture (m3/m3)	Soil Temperature (°C)	RMS Height (cm)	CL (cm)	Sand Content (%)	Clay Content (%)
Radar	5.406	40	0.02–0.44	−50~50	2	20	40	30
Passive	6.925	55	0.02–0.44	−50~50	2	20	40	30
	36.5

. We assumed that the roughness parameters and soil texture of a fixed pixel did not change over time, and the parameter was set to a constant value.

Figure 5 shows the behavior of SSI vs. FTI based on simulation data. To highlight the relationship between the SSI and FTI in the F/T transition, the simulation data from Table 4 with soil moisture higher than 0.2 m³/m³ were selected when obvious dielectric constant changes occurred during the F/T transitions. A logistic relationship is observed between the SSI and FTI with a coefficient of determination of 0.97 from the simulation data in Figure 5.

3.2.2. Effect of Vegetation on F/T Determination

The scatters between the satellite-based SSI_S1 and FTI over cultivated land, grassland, and forest in plots Figure 6a–c follow the logistic curve with a coefficient of determination from 0.88 to 0.91. Even though the relationship between the SSI_soil and FTI in Equation (8) was established for bare land, it also applies for vegetated surfaces. The correlation coefficient is slightly lower over grassland and forest. Additionally, vegetation cover affects the surface freeze onset. The freezing onset of most vegetation-covered surfaces occurs later than that of bare land due to the insulation effect of the vegetation layer [67]. Forests and grasslands are widely distributed in the study areas. Considering the non-negligible effect of vegetation, LAI data were introduced during algorithm development. According to the water cloud model (WCM) [68,69], the backscatter from vegetation can be expressed as follows:

$${\sigma }_{veg}^0=E\times V_1\times \cos \left(\theta \right)\times \left(1-\exp \left(-\frac{2\times F\times V_2}{\cos \left(\theta \right)}\right)\right)$$

(9)

where E and F are the parameters of the model and $\theta$ is the incidence angle of radar. V1 and V₂, which describe the canopy, have various descriptors in previous studies [69,70,71,72] and can be described as an exponential form of the LAI. Since SSI and backscatter are linearly correlated, the relationship between backscatter from vegetation and LAI in WCM is used to express the relationship between the SSI from vegetation (SSI_veg) and the LAI. According to Equation (9), the relationship between SSI_veg and LAI can be simplified as follows:

$${\mathrm{SSI}}_{\mathrm{veg}}=E\times {\mathrm{LAI}}^{L1}\times \left(1-\exp \left(-F\times {\mathrm{LAI}}^{L2}\right)\right)$$

(10)

where E, F, L1, and L2 are fitting coefficients. The influence of vegetation based on Equation (10) was added to Equation (8), and the SSI can be written as follows:

$$\mathrm{SSI}=\frac{a}{1+b\times \exp \left(-c\times \mathrm{FTI}\right)}+d+E\times {\mathrm{LAI}}^{L1}\times \left(1-\exp \left(-F\times {\mathrm{LAI}}^{L2}\right)\right)$$

(11)

where a, b, c, d, E, F, L1, and L2 are the function coefficients. Nonlinear fitting was performed pixel by pixel based on the time series of the SSI_S1, FTI and LAI data. In this method, the relationships were established at the 1 km scale. The regressed coefficients were applied to the daily 1 km FTI and LAI data to estimate the new time-continuous 1 km SSI data (SSI_new). Taking the three in situ sites in Figure 6 as an example, after comprehensively considering the soil information from the FTI and vegetation from the LAI, the coefficient of determination between the SSI_new and SSI_S1 increased to 0.92, 0.89, and 0.88 for site 12, site 3, and site 17, respectively. Figure 7 shows the relationship between SSI_S1 and SSI_new both considering and not considering the LAI at sites 12, 3, and 17. A small increase in the coefficient of determination can be seen when using LAI information. This improvement was mainly reflected when the value of SSI was relatively small.

3.2.3. Pixel-Based Threshold of SSI for Detecting F/T

To detect surface F/T at the 1 km grid from the SSI, a threshold (T) needs to be determined. Previous studies basically defined T as 0.5 [23,50] or optimized it using in situ soil/air temperature data [29,30] when detecting land surface F/T with radar data. However, the feasibility of T with a fixed value or optimized value is limited by certain regions. Considering the heterogeneity of the land surface (such as land cover and terrain), a pixel-based T is defined to monitor the land surface F/T status with SSI over different areas. The T is defined as follows:

$$\mathrm{T}=\{\begin{array}{c}\mathrm{SSI}\_\mathrm{T} \mathrm{T}\_\min \le \mathrm{SSI}\_\mathrm{T}\le \mathrm{T}\_\max \\ \mathrm{T}\_\mathrm{refer} \mathrm{SSI}\_\mathrm{T}<\mathrm{T}\_\mathrm{minor} \mathrm{SSI}\_\mathrm{T}>\mathrm{T}\_\max \end{array}$$

(12)

where SST_T is the threshold defined with SSI and T_min, T_max, and T_refer are the thresholds defined with ERA5-land air temperature. The ERA5-land air temperature data were used here to ensure that T is effective since the SSI_T could be very low or high and lead to low accuracy in some pixels with strong heterogeneity of land cover. To analyze the possible maximum and minimum frozen conditions, T_min and T_max are defined with time series ERA5-land daily minimum and maximum temperatures, respectively. We assume that the SSI_T should be within the range between T_min and T_max. If not, T_refer is defined with the time series ERA5-land daily 1:00 AM temperature, which corresponds to the overpass time of AMSR2. To find T_min, T_max, and T_refer, the F/T classification accuracy of SSI was validated for different threshold values (from min of SSI to max of SSI with an increment of 0.01) with ERA5-land air minimum, maximum, and 1:00 AM temperatures, respectively. The F/T status of ERA5-land air temperature is defined with a threshold of 0 °C. Then, T_min, T_max, and T_refer are found with threshold values when the accuracy is highest. However, the biases between the in situ soil temperature and ERA5-land hourly air temperature over the four study areas are −1.18 °C, −1.18 °C, −7.83 °C, and −2.1 °C over the Genhe watershed, Saihanba, Naqu, and RISMA, respectively. Considering the larger bias of ERA5-land hourly air temperature over the Naqu area, to avoid introducing more error, the SSI_T is the threshold over the Naqu area. The cold bias over the Naqu area was also reported by previous studies [73,74,75]. To define the SSI_T, we estimated the probability density distribution curve of the SSI during late autumn to early spring. The SSI_T is fitted as the SSI value that is symmetrical to the minimum values of SSI (SSI_min). The mean of the probability density distribution curve is the axis of symmetry. Figure 8 shows an example of determining SSI_T based on a histogram and the probability density distribution curve of SSI_new during late autumn to early spring. Finally, the land surface F/T status was detected by SSI_new based on Equations (3), (4), and (12). To compare the land surface F/T detection results with S1, Equation (12) was also applied to SSI_s1 to monitor the F/T status.

3.3. Evaluation Metrics

The F/T classification accuracy was evaluated using three indices [33] as follows:

$$\mathrm{F}\_\mathrm{right}=\mathrm{FF}/\left(\mathrm{FF}+\mathrm{FT}\right)*100\%$$

(13)

$$\mathrm{T}\_\mathrm{right}=\mathrm{TT}/\left(\mathrm{TT}+\mathrm{TF}\right)*100\%$$

(14)

$$\mathrm{Total}\_\mathrm{right}=\left(\mathrm{FF}+\mathrm{TT}\right)/\left(\mathrm{FF}+\mathrm{FT}+\mathrm{TF}+\mathrm{TT}\right)*100\%$$

(15)

where FF and TT represent the correctly classified numbers of frozen soil and thawed soil, respectively. FT and TF represent false detection flags of frozen soil and thawed soil, respectively. Equations (13)–(15) represent freeze accuracy, thaw accuracy, and overall accuracy, respectively. The in situ F/T state is defined by soil temperature data corresponding to the overpass time of AMSR2 with a threshold of 0 °C.

4. Results

4.1. F/T_new Classification Assessment

In this section, the 5 cm depth of in situ soil temperature was used to evaluate the performance of the new F/T algorithm. Figure 9 shows the classification accuracies of F/T_new at each site from June 2017 to May 2019 over the Genhe watershed, from August 2018 to December 2019 over Saihanba and from January 2018 to December 2019 over Naqu and RISMA. As shown in

Table 5. Classification accuracy (%) statistics for F/T_new over four study areas.

Area	F_Right	T_Right	Total_Right
Genhe watershed	93.64	93.19	93.29
Saihanba	93.85	89.08	90.79
Naqu	93.13	85.77	88.53
RISMA	85.24	90.2	88.1

, the total accuracies for all ground sites were 93.29%, 90.79%, 88.53%, and 88.1% over the Genhe watershed, Saihanba, Naqu, and RISMA, respectively. The mean total accuracy of all study areas was 90.18%, which is higher than the required measurement accuracy (90%) of soil F/T based on the essential climate variable (ECV) observation requirements from the Global Climate Observing System (GCOS) (https://gcos.wmo.int/en/essential-climate-variables/about/requirements, accessed on 8 July 2021). Furthermore, the overall accuracies over four study areas are higher than the baseline mission requirement of 80% accuracy of SMAP [23].

The accuracy at site 11 over the Genhe watershed area with farmland is lower than that at other sites with lower thaw accuracy. This is because the roughness of ground soil over the farmland decreased after harvest and plowing. The smoother surface caused by farming activities in the autumn would lead to a decrease in backscatter coefficients [76]. In this situation, the drop in backscatter coefficients before surface freezing may lead to false frozen flags being detected by F/T_new. The accuracies of CD07, MS3527, and P1 over the Naqu area are significantly lower than those at other sites. Figure 10 shows the time series of the SSI_S1, in situ soil moisture and ERA5-land snow depth at site MS3527 over the Naqu area. At site MS3527, SSI_S1 in winter is higher (approximately 0.3) in the case of deep snow (~50 cm in December 2018) due to the effect of snow volume scattering. In Figure 10a, soil moisture increases to 0.2 m³/m³ in May and June with the thawing of the land surface and then decreases to 0.02~0.15 m³/m³ in July and August due to drought. The SSI_S1 would decrease to −0.1 with decreasing soil moisture in summer. With the combined effect of deep snow in winter and drought in summer, the high values of SSI_S1 in winter (~0.3) were larger than the low values of SSI_S1 in summer (~−0.1), leading to false frozen flags in summer. The accuracy over RISMA with cultivated surfaces is lower than that over other study areas. Figure 11 gives the time series of the SSI_new, ERA5-land temperature, in situ soil temperature and the threshold of F/T classification at site MB3 over the RISMA area. The pink dotted line is the threshold of SSI_new. In Figure 11, SSI_new at site MB3 in winter (~0.4) is as high as SSI_new in spring, which causes false thaw flags in winter and lower frozen accuracy due to the effect of snow melting. The high SSI_new in winter usually corresponded to the increase in ERA5-land air temperature when snowmelt may occur.

4.2. Comparisons of the SSI_S1B and SSI_new Results

In this section, to see the performance of the new algorithm (SSI_new) in determining the surface F/T status, we presented an example of the spatial distribution of SSI_S1B, SSI_new, the retrieved F/T maps of S1B (F/T_S1B), the retrieved F/T maps of SSInew (F/T_new) and ERA5-land air temperature over the Genhe watershed. Then, time series of F/T detection results were presented over the study areas. Figure 12 provides examples of the SSI_S1B ((a), (f)), the SSI_new ((b), (g)), the F/T_S1B ((c), (h)), the F/T_new ((d), (i)), and the ERA5-land air temperature ((e), (j)) on 24 October 2018, 5 November 2018, respectively (i.e., the freezing period) over Genhe watershed. The spatial resolutions of SSI_S1B and SSI_new are 1 km, and the ERA5-land air temperature is 10 km. Although ERA5-land air temperature has a coarser spatial resolution than SSI_new, it can still provide temperature distribution information for reference. Most of the SSI_S1B and SSI_new values over the study area decreased from approximately 0.5 to nearly 0 during this period, which indicated freezing of the land surface. The images in Figure 12b,g are the SSI_new obtained from this study, which has similar spatial details to the SSI_S1B (Figure 12a,f), but the value is higher than that of the SSI_S1B. At the bottom of Figure 12f, the backscattering coefficients showed a distinct bright stripe, which may be caused by noise. In Figure 12g, the bright stripe was significantly improved by SSI_new. The images in Figure 12c,d and Figure 12h,i can be compared with the ERA5-land air temperature images shown in Figure 12e,j. The maps of F/T_S1B and F/T_new correspond well with the ERA5-land air temperature. The land surface was mainly unfrozen on 24 October 2018 and began to freeze on 5 November 2018. In Figure 12e, the air temperature in the southwestern part of the study area is approximately 2 °C, and the corresponding F/T_new states were mainly thawed; the temperature in the northeastern part of the study area is approximately 0 °C, and the corresponding F/T_new states were both thawed and frozen.

Figure 13, Figure 14 and Figure 15 show time series of F/T_new, F/T_S1B and air temperature over 13 days over the Genhe watershed, Saihanba and RISMA areas, 2018. During the 13 days of each area, F/T_new provided daily F/T detection maps, while F/T_S1B provided only two days of F/T detection maps. Figure 13 shows the results over the Genhe watershed from 24 October 2018 to 5 November 2018. On 24 October, the F/T_new estimated that the study area was partly frozen, and on 25 October, it was mainly thawed, both of which correspond well with temperature. During 26–29 October, the F/T_new maps showed increasing amounts of frozen soil as the air temperature dropped. From 30 October to 2 November, the air temperature increased, and the surface was again estimated to be thawed in the F/T_new maps. After 3 November, the air temperature suddenly dropped, and the surface soil was estimated to be mainly frozen in the F/T_new maps. Figure 14 shows the results over Saihanba from 10 November 2019 to 22 November 2019. During this period, the land surface was undergoing a freezing process and almost completely frozen from 15 November. Daily F/T_new maps captured the thawing on 13–14 November and corresponded well with ERA5-land air temperature. Figure 15 shows the results over RISMA from 5 March 2019 to 17 March 2019. The surface was predominantly frozen during this time period, but a brief thaw occurred from 13 March to 15 March. The advantage of F/T_new over F/T_S1B in terms of temporal resolution is evident in capturing the land surface F/T states during the F/T transition period (i.e., spring and autumn). As mentioned in Section 3.2.3, the ERA5-land air temperature has a significant bias in the Naqu area, so the comparative analysis between daily F/T_new and ERA5-land air temperature is not shown here.

5. Discussion

5.1. The Feasibility of the Method Developed in This Study

In this section, the feasibility of the method developed in this study was discussed in three stages: (1) analysis of the quality of the regression model in Equation (11) between Sentinel 1 radar, AMSR2 radiometer and GLASS LAI data; (2) exploration of the influence of the difference in the surface F/T state between the S1B and AMSR2 overpass times; and (3) discovery of the role of LAI data in the method proposed in this study.

The purpose of this paper was to estimate the land surface F/T state with a high spatial and temporal resolution and to investigate a new approach for detecting the F/T state at a high spatial resolution of 1 km by combining radar and passive remote sensing data. The high temporal and spatial surface F/T detection algorithm developed in this study is particularly useful in areas with heterogeneous land classes, soil types, and topographic distributions. Typically, for such areas, the soil F/T transitions also progress heterogeneously according to the target characteristics. In previous studies, the retrieval of high-resolution land parameters, such as the F/T state and soil moisture, was mostly based on a model established on the same scale by combining multisource remote sensing [35,38,39,57,77]. The other group of methods estimated high-resolution land parameters from coarse-scale data by introducing ancillary data to describe the spatial heterogeneity [78,79]. The method presented here is based on a nonlinear regression model between active microwave, passive microwave, and optical remote sensing data, including Sentinel 1 radar, AMSR2 radiometer, and GLASS LAI product. This study first combined these three categories of data to provide high-resolution F/T estimates both temporally and spatially. The quality of SSI_new obtained by the regression model in Equation (11) was evaluated with the RMSE between SSI_new and SSI_S1, which is shown in Equation (16):

$$\mathrm{RMSE}=\sqrt{\frac{1}{m}{\displaystyle \sum }_{i=1}^m{\left(f\left(y_i\right)-y_i\right)}^2}$$

(16)

where m is the number of data points used for regression, $f\left(y_i\right)$ is the SSI_new, and $y_i$ is the SSI_S1. Figure 16 provides the histogram and spatial distribution of RMSE over four study areas. The artificial surfaces and water bodies were masked. Except for the some mixed pixels covered with cultivated land and grassland or grassland and forest, the RMSE was lower than 0.2 over the Genhe watershed, Saihanba, and Naqu areas. Most of the RMSE in the RISMA area are below 0.25. The smaller the value of RMSE is, the smaller the difference between SSI_new and SSI_S1. Generally, smaller RMSE values obtained in Equation (16) correspond to a more reliable SSI_new and a higher probability of accurate monitoring of land surface F/T. Figure 17 provides two examples of the relationships among SSI_S1, FTI, LAI, and SSI_new from the Genhe watershed, one with low RMSE (a–c) and another with high RMSE (d–f). We found that even for the high RMSE example (Figure 17d–f), the regression model could capture the soil seasonal characteristics and soil state changes. In Figure 17d, the SSI_S1 increased during winter (the value of FTI was approximately 10) due to the effect of snow. Volume scattering and absorption of the microwave signal within the snow layer increases with deepening snow cover [25], which possibly leads to an increase in the SSI_S1. This finding also corresponded well with Figure 10. These results demonstrate the feasibility of the regression model proposed in this paper.

The high consistency of the surface F/T state between the S1 and AMSR2 overpass times would decide the basis for algorithm development. As shown in Table 2, the overpass time of S1B over the Genhe watershed (22:00) differs the most from that of AMSR2 (1:30). To combine S1B and AMSR2 to detect the surface F/T state, an analysis was performed to determine the difference in surface F/T indicated by the in situ 5 cm soil temperature between 10:00 PM and 1:30 AM over the Genhe watershed. Figure 18 shows the comparison of in situ 5 cm soil temperature between 10:00 PM and 1:30 AM over cultivated land, grassland, and forest in the Genhe watershed. The temperatures between the two times show little difference, especially at approximately 0 °C. The consistency of surface F/T flagged at 10:00 PM and 1:30 AM was 93.72%, 98.85%, and 99.04% for cultivated land, grassland, and forest, respectively. The consistency of cultivated land was lower due to a lack of vegetation cover, which can cause a small temperature difference between 10:00 PM and 1:30 AM.

The presented method uses the FTI and LAI to estimate the SSI_new and further detects surface F/T. The FTI and SSI, as microwave-based F/T indices, are suitable for monitoring the surface F/T status. LAI data were introduced into the F/T index-based logistic model to compensate for the effect of the vegetation layer on surface F/T status determination. The vegetation layer directly influences the observations, particularly the radar signal. The vegetation layer also thermally insulates the soil beneath it, affecting the soil freezing and thawing processes. To clarify the role of LAI in detecting the surface F/T state, the F/T index-based logistic model in Equation (8) without considering the LAI was executed to estimate the SSI_{new_noLAI}. Then, the surface F/T state was detected by the SSI_{new_noLAI} using the same method as in the SSI_new. The in situ soil temperature was used to evaluate the accuracies of F/T_new and F/T_{new_noLAI}, as shown in Figure 19. The x-axis and y-axis show the accuracies of F/T_new and F/T_{new_noLAI}, respectively.

Table 6. Total classification accuracy (%) statistics for F/T_new and F/T_{new_noLAI}.

	Area	Cultivated	Grassland	Forest
F/Tnew_noLAI	Genhe	88.66	88.23	86.16
	Saihanba	/	91.3	90.22
	Naqu	/	86.4	/
	RISMA	86.87	/	/
F/Tnew	Genhe	92.24	93.36	93.79
	Saihanba	/	91.53	90.34
	Naqu	/	88.53	/
	RISMA	88.1	/	/

provides the classification accuracy statistics for F/T_{new_noLAI} and F/T_new over different land cover types. Figure 20 shows the frozen days of F/T_new and F/T_{new_noLAI} from June 2017 to May 2018 over the Genhe watershed area as an example. The two indices have similar spatial distributions, but F/T_new, which introduces high-resolution vegetation information, has more spatial details. In interlaced area of cultivated and grasslands in the southwest Genhe watershed area, frozen days of F/T_new in Figure 20a are more than that of F/T_{new_noLAI} in Figure 20b over cultivated area. In addition, in interlaced areas of grassland and forest, frozen days of F/T_new in Figure 20a are more than that of F/T_{new_noLAI} in Figure 20b over grasslands. The difference of frozen days between F/T_new and F/T_{new_noLAI} mainly occurs in the land cover which accounts for a small proportion of all land covers within passive microwave pixel. From the comparison at these four regions shown on Table 6, the overall accuracy of F/T_new is greater than F/T_{new_noLAI} over different land cover types, especially in Genhe watershed where the heterogeneity of land cover within passive microwave pixel scale is greater. The introduction of LAI data obviously improved the F/T detection accuracy in the case of heterogeneous land cover within the passive microwave pixel scale. This is because LAI can indicate high-resolution vegetation green-up date which is close to surface thawing date [80,81,82]. In summary, by introducing LAI data, the validation accuracies of the vegetation-covered surface improved, and the F/T maps had more spatial details.

5.2. F/T Classification Accuracy Comparisons between F/T_S1 and F/T_new

Figure 21 provides comparisons of the thawed (a), frozen (b), and total (c) accuracies between F/T_S1 and F/T_new over four study areas at a 1 km resolution.

Table 7. Classification accuracy (%) statistics for F/T_S1 and F/T_new.

	Area	FF	FT	TT	TF	F_Right	T_Right	Total_Right
F/TS1	Genhe	357	28	486	55	92.73	89.83	91.04
	Saihanba	221	16	229	62	93.25	78.69	85.23
	Naqu	558	38	555	206	93.62	72.93	82.02
	RISMA	236	60	365	80	79.73	82.02	81.11
F/Tnew	Genhe	368	17	505	36	95.58	93.35	94.28
	Saihanba	221	16	242	49	93.25	83.16	87.69
	Naqu	552	44	649	112	92.62	85.28	88.5
	RISMA	247	49	417	28	83.45	93.71	89.61

provides the classification accuracy statistics for F/T_S1 and F/T_new. The result comparisons were calculated based on all available S1 data for each study area, as shown in Table 2. The mean thaw, freezing, and total accuracies for almost all in situ sites of F/T_new were higher than those of F/T_S1 over the four study areas. The F/T_S1 estimates were retrieved using only S1 radar data, and the F/T_new estimates included the combined use of S1 and FTI data. Multifrequency passive microwave remote sensing is sensitive to both soil moisture and soil temperature and is suitable for monitoring surface F/T. Radar signals are susceptible to roughness and vegetation when monitoring surface F/T. From our results comparison, F/T_new, which contains passive information, has higher accuracies than F/T_S1. Additionally, anomalies caused by noise may exist in the S1 data (Figure 12b). SSI_new could reduce abnormal phenomena by considering time series data. Therefore, F/T_new detection had a better performance than F/T_S1.

5.3. F/T Classification Accuracy Comparisons between F/T_new and F/T_AMSR2

Figure 22 provides comparisons of the thaw (a), frozen (b), and total (c) accuracies between F/T detected from AMSR2 (F/T_AMSR2) and SSI_new at the site level.

Table 8. Classification accuracy (%) statistics for F/T_AMSR2 and F/T_new.

	Area	FF	FT	TT	TF	F_Right	T_Right	Total_Right
F/TAMSR2	Genhe	4529	725	6897	134	86.2	98.09	93.01
	Saihanba	2262	311	3819	392	87.91	90.69	89.64
	Naqu	4868	513	5958	1024	90.47	85.33	87.57
	RISMA	2891	795	5083	386	78.43	92.94	87.1
F/Tnew	Genhe	4888	366	6555	476	93.03	93.23	93.29
	Saihanba	2823	189	4353	539	93.73	88.98	90.79
	Naqu	6729	541	8136	1384	92.56	85.46	88.53
	RISMA	3273	570	5088	559	85.17	90.1	88.1

provides the classification accuracy statistics for F/T_AMSR2 and F/T_new. The total accuracy of F/T_new was higher than that of F/T_AMSR2, while F/T_AMSR2 had a higher thaw accuracy. The theoretical soil penetration depth of the C band is 1–5 cm [83]. Unlike the F/T detection algorithm of AMSR2, which could reflect the 5 cm depth F/T state because it was parameterized by the 5 cm depth soil temperature, the F/T_new detected F/T using the C band, which cannot reflect the soil F/T at the 5 cm depth. Even though the soil at the 5 cm depth was thawed, it could not be detected by F/T_new with the C band. Thus, F/T_new detected more frozen soil than F/T_AMSR2. Figure 23 shows the frozen days calculated from the F/T_new (a) at 1 km and F/T_AMSR2 (b) at 25 km from June 2017 to May 2018 over the Genhe watershed area as an example. The white color in Figure 23a indicates artificial surfaces and water bodies that were masked. Additionally, the spatial patterns of frozen days calculated from the F/T_new and F/T_AMSR2 data were similar over the study area, with the spatial trend of fewer frozen days in the southwestern part and more in the northeastern part of the study area. The frozen days calculated from the F/T_new obviously have more details than that of AMSR2.

6. Conclusions

This paper investigated a new approach for detecting daily land surface F/T at 1 km with AMSR2, S1 and LAI data over the Genhe watershed, the Saihanba and Naqu area located in China and the RISMA network located in Canada, respectively. The proposed approach is derived based on a logistic relationship between the SSI_S1 derived from Sentinel-1 SAR data and AMSR2 FTI. The logistic relation was confirmed by AIEM simulations and observational data. Since vegetation affects microwave signals and the freezing time of the land surface, LAI data were introduced to express the influence of vegetation and provide spatial patterns at the kilometer scale for F/T status determination. The daily F/T state at 1 km could be estimated from the SSI_new obtained from FTI and LAI data. The in situ 5 cm depth soil temperature data from four soil temperature measurement networks were used to evaluate the classification accuracies of F/T_new. The mean total accuracies obtained from in situ sites are 93.29%, 90.79%, 88.53%, and 88.1% over the Genhe watershed, Saihanba, Naqu and RISMA, respectively. By introducing LAI data, the classification accuracies of F/T_new over the vegetation-covered surface maximum increased by 7.6%. The results comparison between F/T_new and F/T_S1 at 1 km showed that this method could reduce some S1 data anomalies and obtain higher accuracies. Additionally, the accuracy comparison results between F/T_new at 1 km and F/T_AMSR2 at 25 km showed that the high-resolution results could more accurately reflect the land surface state. In this study, daily 1 km land surface F/T could be detected by combining with Sentinel-1 SAR, AMSR2, and LAI data. The accuracy increased by 5.17% and 0.85%, reaching 90% and 90.18% when compared with the stand-alone S1 F/T results at 0.01° and the AMSR2 F/T results at 0.25° resolution. Despite this method may be affected by the change in roughness over cultivated areas caused by agricultural activity, leading to false detection flags for thaw surfaces. The overall accuracy over cultivated areas was still above 90%. This indicated that combined SAR and passive microwave remote sensing can monitor relatively high spatial resolution surface F/T at 1 km with high accurately than through passive satellite sensors. The accuracies have improved over vegetation-covered surfaces by considering the effect of vegetation.

The combination of S1 and AMSR2 has the potential to detect land surface F/T states with high spatial and temporal resolutions over a large scale, even at the global scale. In future work, we will apply a temperature-based filter to remove false detections during winter caused by the snow effect. In addition, the new approach developed in this work is established based on time series data for each pixel, and the algorithm coefficients obtained from the available Sentinel-1 data could be used to estimate the SSI_new during periods before Sentinel-1 was launched. Thus, the daily land surface F/T state at 1 km before the launch of Sentinel-1 could be estimated. It is expected to derive a long series of daily land surface F/T states at a spatial resolution of 1 km since 2002.