Snow Depth Inversion in Forest Areas from Sentinel-1 Data Based on Phase Deviation Correction

Abstract

At present, snow depth inversion based on active microwave remote sensing is concerned essentially with areas having a relatively simple underlying surface. The existence of forests reduces the sensitivity of microwaves to snow, which often makes the snow depth inversion results uncertain. This paper presents a snow depth estimation algorithm for forest areas by introducing a forest phase to characterize the effect of forests on backscattering electromagnetic wave. Firstly, the interferogram is generated with the differential interference of two-pass master-slave Synthetic Aperture Radar (SAR) images, and the real phase under snow cover condition is obtained by phase unwrapping. Secondly, the phase models for forest and non-forest areas are constructed. The effects of forest cover are modeled as forest phase in the forest phase model, which is estimated under the assumption of snow depth consistency on both sides of the boundaries between forest and non-forest areas. Finally, snow depth is estimated by the snow phase-depth model. The correctness of the proposed forest snow depth inversion algorithm was verified by taking the Jiagedaqi area of Greater Xing’an Mountains as the study area and sentinel-1 dual polarization images as the data source. Finally, the snow depth distribution of the study area was obtained with a spatial resolution of 30 m on 7 December 2020. The experimental results show that the snow depth values estimated in Jiagedaqi area are mainly between 40–120 cm, and the average snow depth value is 80.27 cm. Taking the snow depth value of 84.69 cm reckoned from hourly accumulated snowfall in Jiagedaqi as the reference snow depth, the results of the estimated snow depth are relatively consistent and well-founded. With the introduction of the forest phase, the average snow depth values estimated in the forest area increase by 5.98 cm, which reduces the underestimation of the snow depth in forest areas.

1. Introduction

Snow cover plays a highly critical role in the global water cycle and energy exchange. As one of the major snow-covered areas in China, the Great Xing’an Mountains contain abundant snow resources. Accurate snow depth data are important for research on hydrological process, climate change, and natural disaster prediction [1,2,3]. On-site observation data of snow cover in winter are often limited by many factors, such as local weather conditions, scale of human resources, logistical support, and, particularly, complex terrains. Furthermore, the observation data of meteorological stations are often limited to a certain spatial range, which cannot meet the needs of monitoring large-scale snow cover dynamic changes. Therefore, remote sensing technology plays an important role in monitoring snow hydrological processes [4,5,6]. The remote sensing spectral bands commonly used for snow detection are visible light, near-infrared, and microwave bands. Among these, visible and near-infrared electromagnetic waves cannot penetrate snow and objects above snow (such as vegetation and forests), and are only suitable for detecting snow cover on bare land. Given the ability of microwaves to penetrate objects above snow and the snow layer itself, they can be used to estimate important snow parameters, such as snow depth and Snow Water Equivalent (SWE) [7].

At present, the snow depth inversion for forest areas mainly relies on passive microwave remote sensing. For example, Foster et al. [8] proposed a snow depth inversion algorithm under forest coverage. By modifying the regression coefficient in [9] and introducing the forest coverage, this algorithm can estimate better snow depth in most cases, but the snow depth values in the forest area are underestimated. Derksen et al. [10] proposed a snow depth inversion algorithm using a linear mixed pixel decomposition model when studying the snowfall in the forest area of western Canada. This method mainly considers four different surface types, namely, bare land, coniferous forest, deciduous forest, and sparse forest, and the weight of the mixed pixel decomposition is the coverage area ratio of each type in a pixel. To estimate snow depth with forest coverage in the main snowfall areas in China using passive microwave Scanning Multichannel Microwave Radiometer (SMMR) and Special Sensor Microwave/Imager (SSM/I) data, Che et al. [11] proposed different regression coefficients to obtain China’s long-term series snow depth dataset. To overcome the error of snow depth estimated under deep snowpack conditions without sacrificing performance under relatively thin snowpack conditions, Yang et al. [12] developed an empirical snow depth inversion algorithm suitable for FY-3D satellite images. An evaluation using meteorological station observations in 2014 and 2015 demonstrated that the root mean square error and bias of snow depth values estimated from FY-3D images were 6.6 and 0.2 cm, respectively, and their algorithm has advantages over other similar algorithms. Although passive microwave remote sensing-based snow depth inversion for forest area has achieved many results, the low spatial resolution of passive microwave remote sensing images limits its extensive application to snow depth inversion on a finer scale.

Active microwave remote sensing has shown advantages in the field of snow depth inversion due to its high spatial resolution and high sensitivity to snow characteristic parameters. Currently, various snow depth inversion methods based on active microwave remote sensing have been proposed. For example, in snow depth inversion based on SAR interferometry, Shi et al. [13] demonstrated a technique to map snow cover with both the backscattering and coherence (with and without snow covered images); the study found that coherence measurements provide a much easier way to map snow-covered areas. Guneriussen et al. [14] confirmed that the microwave refraction in dry snow has a significant effect on the interference phase by analyzing the ERS-1/2 C band SAR interferometric images, and deduced the relationship between the change in SWE and the interference phase. In snow depth inversion based on Copolar Phase Difference (CPD), Leinss et al. [15] constructed a model of the relationship between snow microstructure and CPD based on TerraSAR-X VV and HH data, and determined an empirical model of VV and HH CPD and fresh snow depth. In snow depth inversion based on data fusion, Liu et al. [16] proposed a snow depth inversion scheme that fuses remote sensing data and site observation data. On the basis of using differential interferometry based on Sentinel-1 SAR images to estimate snow depth, this scheme further uses the Three-Dimensional Variational (3DVAR) fusion algorithm to integrate snow depth data from virtual stations and real-world observation stations into the snow depth inversion results. It was applied to the Bayinbuluke basin in the Tianshan Mountains, Xinjiang, China, and the results showed that the accuracy of snow depth estimation can be effectively improved. In the snow depth inversion based on interpolation, three kinds of snow depth data, namely, the D-InSAR data obtained from the remote sensing images of Sentinel-1 SAR, the automatically measured data using ultrasonic snow depth detectors, and the manually measured data, are assimilated based on an ensemble Kalman filter. Under the constraint of the measured data, the accuracy of the assimilated snow depth value is higher and can meet the needs of subsequent research [17]. In the snow depth inversion based on the Artificial Neural Network (ANN) method, Santi et al. [18] used ANN to estimate snow depth and SWE using X and Ku band airborne Snow SAR data. In the snow depth inversion based on the Machine Learning (ML) method, Zhu et al. [19] proposed a snow depth inversion method for SAR data based on parameter selection in combination with the Correlation Coefficient Method (CCM) and ML. This method was used to estimate snow depth in northeastern China, and the results were verified by combining the snow depth data from meteorological stations and on-site observation data; the study showed that the best snow depth estimation results can be obtained by combining the parameters selected by the CCM and XGBoost algorithms. Although the current snow depth inversion based on active microwave data has achieved some promising results, most of the related research focuses on the relatively simple case of the underlying surface. Using active microwave remote sensing for complex underlying surfaces, such as forest areas, snow depth inversion faces significant challenges and has opportunities for further development.

In order to study active microwave remote sensing-based snow depth inversion in forest areas, this paper presents a snow depth inversion algorithm based on the correction of snow phase deviation in forest areas by introducing the forest phase to describe forest cover effects. Firstly, two Sentinel-1 C band SAR images under snow and snow-free conditions in Jiagedaqi area of Greater Xing’an Mountains are selected for the two-pass differential interferometry to obtain the interferogram, and the real phase is obtained by phase unwrapping. Secondly, phase models for forest and non-forest areas are constructed, and a forest phase is introduced into the forest phase model to describe the impact of forest cover. The forest phase is estimated by the assumption of the consistency of snow depth on both sides of the boundary between forest and non-forest areas, and modeling the effects of snow cover in both areas as snow phases. Finally, according to the relationship between the snow phase and snow depth, the distribution of the snow depth is estimated. The results are compared with the reference snow depth, proving the effectiveness of the proposed method. The comparison of the snow depth before and after the introduction of the forest phase and the source of error are analyzed.

2. Study Area and Data Description

2.1. Overview of the Study Area

The study area of this paper is Jiagedaqi area, which is located in Oroqen Autonomous Banner of Nei Mongol Autonomous region, but administratively belongs to Greater Xing’an Mountains of Heilongjiang province, as shown in Figure 1a; the figure is the standard base map of China, from the Ministry of Natural Resources Mapping and Geographic Information website (https://zwfw.ch.mnr.gov.cn/index, accessed on 3 August 2021). Figure 1b is the image map of Jiagedaqi, the geographical coordinates are 123°45′~124°26′E and 50°09′~50°35′N, with a total area of 1359 km². The forest area is 902.07 km², accounting for 66.38% of the total area, and the local vegetation is mainly Larix gmelinii forest, mixed with Pine sylvestris forest, Betula platyphylla forest, etc. The Gan River, the only surface water resource in the whole Jiagedaqi area, is a branch of the Nenjiang River, and passes through the area from west to east, with a length of 55 km. Figure 1c shows the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data with a spatial resolution of 30 m, which can be downloaded from http://dwtkns.com/srtm30m/, accessed on 14 March 2021. The terrain is high in the northwest and low in the southeast, with an average elevation of 472 m, belonging to a low mountain hilly area. The northwest wind prevails in winter, with the windward slope in the northwest and the leeward slope in the southeast. The northwest slopes have a deeper snow distribution than the southeast slopes.

The Jiagedaqi area has a cold temperate continental monsoon climate, with a cold and long winter and abundant snowfall. The average annual temperature is −1.5 °C, the extreme minimum temperature is −45.4 °C, and the frost-free period is 85–130 days. The first snowfall in 2020 occurred on 21 October. At the beginning of snow cover, the snow particle size is small; however, due to the effect of wind, temperature, and other factors, the snow particle size gradually develops from the surface layer to the deep layer into a fine–medium–coarse particle size structure. As a result of the increase in temperature, the snow constantly melts and freezes, and the snow then becomes frost until it completely disappears. The research on snow depth provides important scientific theories for guiding local forestry production and coping with climate change (note: the data of the study area are from the official government website of Jiagedaqi area, http://www.jgdq.gov.cn/about/, accessed on 12 January 2021).

2.2. Data Description

The C-band Sentinel-1 dual polarization Interferometric Wide (IW) swath Level-1 Single Look Complex (SLC) data product is selected for the experiment, and is provided by the European Space Agency (ESA) through the Copernicus Programme (https://scihub. copernicus.eu/dhus/#/home, accessed on 4 February 2021). Sentinel-1B IW data have two polarizational states, VV and VH, with a temporal resolution of 12 days, and a spatial resolution of 5 × 20 m [20]. The Sentinel-1 dataset is very advantageous in snow depth interferometry due to its high temporal resolution, high spatial resolution, and free availability. Sentinel-1 precision orbit files are generated 21 days after acquiring the Sentinel-1B data. For this paper, 9 August 2020 (snow-free)–7 December 2020 (snow) image pairs were selected as the raw data for interferometric snow depth inversion, taking into account the spatial and temporal baselines, and the image pair coherences in two-pass differential interferometry.

Table 1. Sentinel-1 SAR image information in the study area.

Master Image	Slave Image	Time Baseline/day	Vertical Baseline/m	Polarization Mode
9 August 2020	7 December 2020	120	12.82	VV

lists the experimental data information.

Global land cover classification data (GlobeLand30 V2020) with a spatial resolution of 30 m, released by the Ministry of Natural Resources of China (http://www.globallandcover.com/, accessed on 5 October 2021), is also used in this study. The data accuracy evaluation is led by the Aerospace Information Research Institute of the Chinese Academy of Sciences. Based on the landscape shape index sampling model, a full set of data points are arranged, with a total of more than 230,000 samples; the overall accuracy is 85.72% and Kappa coefficient is 0.82.

3. Methodology

In the proposed snow depth inversion method based on the snow phase, firstly, the interferogram is generated by the two-pass differential interference of the master and slave images, and the real phase is obtained by phase unwrapping. Secondly, according to the scattering mechanism of each scattering target in the SAR signal propagation path, the phase models for forest and non-forest areas are constructed, respectively. In the phase model for forest areas, the forest phase is introduced to describe the influence of forest cover. Finally, the snow depth is estimated according to the relationship between the snow phase and snow depth.

3.1. Interferogram Generation and Phase Unwrapping

To generate an interferogram using two-pass differential interference, it is necessary to use SAR images in both snow-free and snow periods in the same area as the master and slave interferometric images, respectively. That is, let Z^(M) = {Z_i^(M), i = 1, …, N} and Z^(S) = {Z_i^(S), i = 1, …, N} for the master (snow-free) and slave (snow) SAR images, respectively, where i is the pixel index, N is the number of pixels, Z_i^(M) = A_i^(M)exp(jφ_i^(M)) and Z_i^(S) = A_i^(S)exp(jφ_i^(S)) are complex signals of pixel i in the master and slave images, respectively, A_i^(M) and A_i^(S) are the signal modulus values, respectively, and φ_i^(M) and φ_i^(S) are signal phase values, which are related to the signal path and random phases caused by different scattering characteristics. Figure 2 shows the master–slave SAR image interference geometry, where R is the distance from the ground point P to the master satellite antenna S₁, R+ΔR is the distance from the ground point P to the slave satellite antenna S₂, and ΔR is the slant distance of the two paths.

Accordingly, φ_i^(M) and φ_i^(S) can be calculated as follows, respectively:

$${\phi }_i{}^{\left(M\right)}=\frac{4\pi }{\lambda }R+\arg \left(A_i{}^{\left(M\right)}\right)$$

(1)

$${\phi }_i{}^{\left(S\right)}=\frac{4\pi }{\lambda }\left(R+\mathsf{\Delta }R\right)+\arg \left(A_i{}^{\left(S\right)}\right)$$

(2)

where λ is the wavelength, and arg (A_i^(M)) and arg (A_i^(S)) are the random phases of pixel i in the master and slave images, respectively. Assume that the random phases caused by different scattering properties in two SAR images are equal, that is:

$$\arg \left(A_i{}^{\left(M\right)}\right)=\arg \left(A_i{}^{\left(S\right)}\right)$$

(3)

The interferogram can be obtained by multiplying the conjugates of the master and slave image complex signals [21], that is:

$$Z_i{}^{\left(M\right)}{Z}_i{}^{*\left(S\right)}=\left|A_i{}^{\left(M\right)}{A}_i{}^{\left(S\right)}\right|\exp \left\{j\left({\phi }_i{}^{\left(M\right)}-{\phi }_i{}^{\left(S\right)}\right)\right\}=\left|A_i{}^{\left(M\right)}{A}_i{}^{\left(S\right)}\right|\exp \left(-\frac{4\pi }{\lambda }\mathsf{\Delta }R\right)$$

(4)

where * represents complex conjugation. The phase difference in interferogram is the phase principal value ∆φ_i, which ranges between [−π, π]:

$$\mathsf{\Delta }{\phi }_i=-\frac{4\pi }{\lambda }\mathsf{\Delta }R$$

(5)

In order to obtain the real phase value, an integer multiple of 2π should be added to or subtracted from the phase principal value ∆φ_i. Phase unwrapping is the process of restoring the phase principal value to the real phase φ_i, that is:

$${\phi }_i=\mathsf{\Delta }{\phi }_i+2\pi f\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}f=0,\pm 1,\pm 2,\cdot \cdot \cdot$$

(6)

3.2. Snow Phase Model

According to GlobeLand30 V2020 30 m resolution land cover classification data, Jiagedaqi area can be divided into forest and non-forest areas. The ground objects corresponding to non-forest areas are mainly cultivated land, grassland, and artificial surfaces. The classification results are shown in Figure 3. The phase components of forest and non-forest areas are different, so phase models are constructed, respectively.

3.2.1. Snow Phase Model in Non-Forest Area

The backscattering of radar waves penetrating the snow layer and reaching the ground includes three parts: the scattering of radar waves at the air-snow interface, the snow–ground interface, and scattering in the interior of the snow layer [22]. The division between dry and wet snow is related to the temperature of the snowpack. At temperatures below 0 °C, the snowpack is assumed to remain dry, while above 0 °C the snow is considered to be wet [23]. Figure 4a,b shows the microwave scattering mechanisms in wet and dry snow states, respectively. Wet snow is a mixture of ice crystals, air, and water [24]. Due to the high absorption of microwaves in the wet snow layer, the absorption coefficient increases with the increase in liquid water content, so the backscattering mainly comes from the air–snow interface, while the scattering at the snow-ground interface can generally be negligible [25,26,27]. Dry snow is a mixture of ice crystals and air [24]; due to the fact that the difference in dielectric constant between air and dry snow is small enough, most microwaves incident on the snow surface can be transmitted through the interface. In addition, the microwave attenuation in the dry snow layer is very small, and the backscattering received by the radar mainly comes from the snow–ground interface [28].

In fact, snow usually undergoes qualitative changes due to atmospheric and surface temperature, wind, and multiple snowfalls, whereas microwaves are very sensitive to the internal state of snow, and snow layers can cause microwaves to refract multiple times internally, resulting in phase delay errors. In addition, when the snow cover is thick, significant volume scattering can also occur, which affects the interference phase; however, the influence mechanism is complicated. The snow cover generally lasts from November to April of the following year in the study area [29]. From the beginning of snow cover to December in 2020, the temperatures were below 0 °C, and only in November were the temperatures occasionally higher than 0 °C. This paper proposes a simplified model; the influence of volume scattering is not considered, the snow layer is assumed to be dry snow, and radar waves can penetrate the snow layer to reach the ground surface.

The total phase φ_i^t1 under snow cover conditions in non-forest area can be modeled as:

$${\phi }_i{}^{t1}=\left({\phi }_{atm}+{\phi }_{noise}+{\phi }_{orb}\right)+\left({\phi }_i{}^{flat}+{\phi }_i{}^{topo}+{\phi }_i{}^{snow1}\right)$$

(7)

where φ_atm is the phase caused by the atmospheric action, and the factors affecting the atmospheric phase are quite complex, which can be ignored in clear weather; φ_noise is the phase of the system noise. Due to the influences of coherent speckle, thermal noise, and data processing noise in the SAR imaging process, noises inevitably appear in the interferogram. They can be greatly reduced by filtering the SAR image before differential interferometry and on the interferogram before phase unwrapping; φ_orb is the phase caused by the orbit error. The precise orbit file provides accurate satellite position and velocity information, which can be used to estimate the error; φ_i^flat is the phase caused by the flat-earth effect at pixel i. Even in the interferogram of the flat ground environment, the phases of different pixels change periodically in the ground range direction; this phenomenon is called the flat-earth effect [30]. The phase can be estimated by orbit parameters and metadata information; φ_i^topo is the phase difference due to changes in the distance between the satellite and target for topography at pixel i [31], which can be removed by subtracting the topographic interferogram simulated by external DEM data from the interferogram. The topographic interferogram simulated by external DEM data is divided in two steps; in first step, the DEM is radar-coded to the coordinate systems of the master image. Per DEM point, the master coordinate (real valued) and the computed reference phase are saved to a file. Then, the reference phase is interpolated to the integer grid of master coordinates. A linear interpolation based on a Delaunay triangulation is used. The DEM should be 10–15% bigger in its extent (wider and higher) than the interferogram. φ_i^snow1 is the phase caused by the snow layer in the non-forest area at pixel i, which can be obtained from Equation (7):

$${\phi }_i{}^{snow1}={\phi }_i{}^{t1}-\left({\phi }_{atm}+{\phi }_{noise}+{\phi }_{orb}\right)-\left({\phi }_i{}^{flat}+{\phi }_i{}^{topo}\right)$$

(8)

3.2.2. Snow Phase Model in Forest Area

Figure 5 shows the microwave scattering process for snow cover in forest area. The existence of the forest increases the backscattering of forest, while also attenuating the backscattering process of snow and reducing the sensitivity of microwave to snow scattering.

The forest cover effect is modeled as the forest phase φ_i^forest. When modeling the total phase φ_i^t2 obtained by interfering and unwrapping the two SAR images in forest areas with and without the presence of snow cover, in addition to considering the difference in snow phase between forest and non-forest areas, by assuming that the other phases are the same as the above, the total phase φ_i^t2 can be modeled as:

$${\phi }_i{}^{t2}=\left({\phi }_{atm}+{\phi }_{noise}+{\phi }_{orb}\right)+\left({\phi }_i{}^{flat}+{\phi }_i{}^{topo}+{\phi }_i{}^{snow2}+{\phi }_i{}^{forest}\right)$$

(9)

where φ_i^snow2 is the snow phase in a forest area. From Equation (9), the relationship between the snow phase φ_i^snow2 and the forest phase φ_i^forest can be obtained as:

$${\phi }_i{}^{snow2}={\phi }_i{}^{t2}-\left({\phi }_{atm}+{\phi }_{noise}+{\phi }_{orb}\right)-\left({\phi }_i{}^{flat}+{\phi }_i{}^{topo}+{\phi }_i{}^{forest}\right)$$

(10)

Since it is difficult to model the forest phase φ_i^forest, this paper proposes its estimation method. Assume that at the boundary of forest and non-forest areas, adjacent non-forest and forest pixels have the same snow phase, that is, φ_x^snow1 = φ_y^snow2, where x and y represent non-forest and forest pixels on the sides of the boundary, respectively. It is also considered that the phase difference caused by the terrain and flat-earth effect in the neighborhood pixels can be negligible, and the forest phase close to the boundary is deduced as:

$${\phi }_y{}^{forest}={\phi }_y{}^{t2}-{\phi }_x{}^{t1}$$

(11)

Take the mean φ_y^forest as an estimate value of the forest phase, that is:

$${\phi }_i{}^{forest}=Mean\left({\phi }_y{}^{forest}\right)=Mean\left({\phi }_y{}^{t2}\right)-Mean\left({\phi }_x{}^{t1}\right)$$

(12)

where Mean(·) is the mean operation.

3.3. Snow Depth Inversion

Figure 6 shows the difference in the propagation paths of microwaves with and without snow.

When there is snow cover, the microwave penetrates the snow to reach the surface by (∆R_a + ∆R_r), and when there is no snow cover, the distance is ∆R_s, and the difference between the two distances is:

$$\mathsf{\Delta }R=\mathsf{\Delta }{R}_s-\left(\mathsf{\Delta }{R}_a+\mathsf{\Delta }{R}_r\right)=\mathsf{\Delta }{R}_1-\mathsf{\Delta }{R}_r$$

(13)

where φ^(snow) = {φ_i^(snow), i = 1, …, N} is the snow phase map, and φ_i^(snow) is the snow phase value at pixel i:

$${\phi }_i{}^{\left(snow\right)}=\left\{\begin{array}{l}{\phi }_i{}^{snow1}=-\frac{4\pi }{\lambda }\left(\mathsf{\Delta }{R}_1-\mathsf{\Delta }{R}_r\right)+2\pi k \hspace{0.33em}k=\pm 0,\pm 1,\pm 2,\cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot forest area\hspace{0.33em}\\ {\phi }_i{}^{snow2}=-\frac{4\pi }{\lambda }\left(\mathsf{\Delta }{R}_1-\Delta R_r\right)-{\phi }_i{}^{forest}+2\pi k \hspace{0.33em}k=\pm 0,\pm 1,\pm 2,\hspace{0.33em}\cdot \cdot \cdot non-forest area\end{array}\right.\hspace{0.33em}$$

(14)

According to Figure 6, the geometric relationship of microwave propagation path can be deduced as:

$$\mathsf{\Delta }{R}_1=\mathsf{\Delta }{R}_r\cos \left({\theta }_i-{\theta }_s\right)$$

(15)

$$\mathsf{\Delta }{R}_r=\frac{d_s}{\cos {\theta }_s}$$

(16)

The snow permittivity ε_s and the refractive index n satisfy the relationship ε_s = n², and the snow density ρ satisfies the relationship ε_s = 1 + 1.6ρ + 1.86ρ³, the refractive index n = sinθ_i/sinθ_s, where θ_i is the incident angle and θ_s is angle of refraction. d^(s) = {d_i^(s), i = 1, …, N} is the snow depth map, and d_i^(s) is the snow depth value at pixel i. According to the above relationship, the snow depth value can be calculated as:

$$d_i{}^{\left(s\right)}=-\frac{\lambda }{4\pi }\frac{{\phi }_i{}^{\left(snow\right)}}{\cos {\theta }_i-\sqrt{\left(1+1.6\rho +1.86{\rho }^3\right)-\sin {\theta }_i{}^2}}$$

(17)

4. Results and Analysis

In order to obtain the local incident angle and the phase information after phase unwrapping, this experiment uses the SAR image processing software SNAP (http://step.esa.int/main/download/snap-download/, accessed on 14 January 2021) developed by ESA to process Sentinel-1 data with two-pass differential interferometry. In the experiment, the master and slave images are set as SAR images obtained during snow-free and snow-covered periods in Jiagedaqi area, respectively. Figure 7 shows the flow chart for obtaining initial data.

Orbit corrections can effectively avoid systematic errors caused by orbit factors. For image registration, the resampling method is based on bilinear interpolation, and by enhancing the spectral diversity, offset correction of the range and azimuth directions is performed on the slave band. The interferogram is obtained by multiplying the conjugate of the complex values of the master and slave images, as shown in Figure 8a; the topographic phase is removed from the interferogram for phase flattening. Using Goldstein phase filtering can greatly reduce the residual error in subsequent processing and improve the phase unwrapping accuracy. In the Linux environment, the phase unwrapping after filtering was performed with the statistical cost network flow algorithm Snaphu; the statistical cost mode is Differential Interferometry (DEFO), and the network flow generation method is Minimum Cost Flow (MCF). The phase unwrapping results are shown in Figure 8b. If the phase unwrapping is performed in the Windows environment, Snaphu needs to be invoked through PowerShell. For range Doppler terrain correction of the image used in this paper, the pixel spacing was set to 30 m, the reference system was WGS84, the map projection was UTM 51N and the resampling method was based on bilinear interpolation. The obtained local incident angle information is shown in Figure 8c.

4.1. Snow Depth Inversion Results

The phase unwrapping map is divided into forest areas and non-forest areas, and the total phase averages of forest areas and non-forest areas on both sides of the boundary are 27.567 and 29.784, respectively. According to Equation (12), the phase caused by the forest effect is −2.217, and, according to Equation (10), the phase diagram of the forest phase removal in the forest area is obtained, that is, the snow phase diagram. A box-plot can be used to visually identify outliers in the data, where a number beyond the interval of mean ±3 std is considered as an outlier. It is used to perform outlier analysis of the snow phase in forest and non-forest areas, respectively. As shown in Figure 9a, there is no abnormal value in the snow phase in non-forest areas, while a snow phase in forest areas larger than 54.8 is recognized as a large abnormal value, and less than 4.38 as a small abnormal value. The locations of the outliers are found on the land cover classification map, as shown in Figure 9b. The outliers are mostly distributed near the edge of the forest area close to the residential area. The number of outlier phases in the forest area accounts for 0.02% of the total number of snow phases in the forest area, which is relatively small. In this experiment, artificial correction was applied, and the large outliers were uniformly set as 54.8, and the small outliers were uniformly set as 4.38, the outlier-corrected box-plot is shown in Figure 10a. After removing the forest phase and correcting the outlier phase, the snow phase can finally be obtained, as shown in Figure 10b.

In this paper, the dry snow density ρ = 0.18 g/cm³ is taken as the empirical value, the radar wavelength λ = 5.5466 cm, the local incident angle is used to replace the satellite central incident angle [32], and the snow depth is estimated according to the relationship between the snow phase and snow depth, that is, Equation (17). The obtained spatial distribution of snow depth with a spatial resolution of 30 m is shown in Figure 11. From the histogram of the snow depth frequency distribution (Figure 12), it can be seen that there are abnormal peaks in the snow depth range of 100–110 cm. These abnormal values are concentrated in the red circle area shown in Figure 11, corresponding to the urban area of Jiagedaqi. Figure 12 also shows the histogram of the snow depth frequency distribution after excluding the urban area, and the peak phenomenon is significantly removed. It can be seen that the urban area has a greater impact on the snow depth inversion results, while the snow distribution state in the urban area can be changed by human activities, such as cleaning snow. The proposed method is mainly aimed at snow depth inversion in a natural state, so human living areas should be excluded as much as possible.

At present, the snow weather forecast mainly relates to snowfall, and lacks prediction and estimation of snow depth. In 2009, the China Meteorological Administration used snow observation data to analyze the relationship between snowfall and snow depth, and concluded that the ratio of snow depth to snowfall is above 1 cm·mm⁻¹ in northeast China [33]. We also found relevant records on the Network Political platform of Shuangyashan (http://www.syswlwz.gov.cn/, accessed on 8 March 2022), where the snowfall and snow depth are usually converted according to a ratio of 1:(10–15). Of course, this ratio is not fixed, and changes with the state of snow. Generally, the ratio of dry snow is greater than that of wet snow; in other words, the lower the water content, the closer the ratio to 1:15. Since the snow in the study area is dry snow, in this paper, we set the ratio between the snowfall and the snow depth as 1:15. It can be seen from Figure 12 that the snow depth is mainly concentrated between 40 and 120 cm, and the average snow depth is 80.27 cm. The hourly snowfall data from Jiagedaqi meteorological station are provided by Xihe Energy Big Data platform (https://www.xihe-energy.com/, accessed on 9 March 2022). The accumulated snowfall was 56.46 mm of water equivalent, as shown in Figure 13; the relationship between snow depth and snowfall of 1:15 was converted, and the snow depth on 7 December 2020 was 84.69 cm. Taking this result as the reference snow depth, it can be seen that the average estimated snow depth is relatively close to it.

4.2. Analysis of Snow Depth Inversion Results

As shown in Figure 12, the estimated snow depth ranges from 0 to 160 cm, and is mainly concentrated from 40 to 120 cm. From the perspective of spatial distribution, the northern area of Jiagedaqi has deeper snow cover than the southern area; these values are about 110 and 60 cm, respectively. The northwest wind prevails in the Jiagedaqi area in winter due to the difference in air pressure between the high pressure of Siberia and the low pressure of northeast China, and the windward slope is in the northwest. The northwest has a deeper snow distribution than the southeast, which is more consistent with the spatial distribution of snow depth in the inversion of this paper. The analysis of snow depth inversion results is described from the comparative analysis of snow depth inversion before and after the introduction of the forest phase in forest areas and the source of snow depth inversion errors.

4.2.1. Comparative Analysis of Snow Depth before and after the Introduction of Forest Phase

Figure 14a shows the snow depth difference map obtained by subtracting the snow depth values for the forest area estimated with and without the forest phase (that is, the traditional algorithm without considering the influence of forest is used to retrieve snow depth), and Figure 14b is the frequency distribution histogram of the difference. From Figure 14b, it can be seen that the differences are positive and mainly distributed between 4 and 8 cm, with an average value of 5.98 cm. The introduction of the forest phase increases the estimated snow depth values in the forest area. From the comparison of the average snow depth of the whole study area, the average snow depth without the introduction of the forest phase is 77.13 cm, the average snow depth with the introduction of the forest phase is 80.27 cm, and the reference average snow depth value is 84.19 cm. The introduction of the forest phase provides a new method for snow depth inversion in forest area.

4.2.2. Incoherence Error Analysis

As an important parameter to measure the complex correlation between the amplitude and phase information of two interference signals, coherence is closely related to the accuracy of interferometry. The coherence coefficient of the two images is defined as:

$$r=\frac{\left(\left|y_1{y}_2{}^{*}\right|\right)}{\sqrt{\left(y_1{y}_1{}^{*}\right)\left(y_2{y}_2{}^{*}\right)}}$$

(18)

where y₁ and y₂ are random signals, r is the coherence coefficient, and the value range is [0, 1], which is used to measure the degree of complex signal from complete incoherence to complete coherence [34]. There are many factors that affect coherence, including change in land cover type, atmospheric effect, temporal baseline, spatial baseline, and noise introduced in data processing.

(1)

Incoherence caused by changes in land cover

The change in the land cover type is the main factor causing the spatial incoherence of SAR images. The areas with obvious changes in the coherence coefficient were selected for analysis. The location is 124°03′~124°11′E and 50°23′~50°27′N. As shown in Figure 15a, the better coherence values are mainly located on the center of the cropped image, and the coherence around the central area is relatively low. It can be seen from the land cover classification map of the corresponding areas (Figure 15b) that the center area is a relatively flat man-made surface area, whereas the surrounding areas are mostly cultivated land, forest areas, and grasslands. The boundary of the coherence change areas is more consistent with the boundary of the change in the ground objects; except for the difference in the coherence caused by the characteristics of the ground objects, this does not rule out the deterioration in the coherence caused by the influence of human activities. The coherence coefficient is an important indicator to measure the quality of the interferogram, and low coherence affects the accuracy of phase unwrapping, which in turn affects the snow depth inversion results.

(2)

Incoherence caused by atmospheric effects

Figure 16 shows the average high temperature, average low temperature, extreme high temperature, and extreme low temperature from January to December in Jiagedaqi area in 2020; the temperature data were taken from the weather network (https://tianqi.com/, accessed on 3 January 2022), which records the temperature data of each city from the weather forecast information of the day. The lowest temperature appeared in January, February, and December, and the highest temperature appeared in June to August. In the experiment, images from August and December were selected. The temperature difference between the two months was greater than 30 °C, and the atmospheric state was different during satellite imaging, so the atmospheric delay error could not be avoided in the interference results. In this experiment, the atmospheric phase was not removed, which may cause the estimated snow depth values in some areas to be larger than the real values.

4.2.3. Error Analysis of Input Parameter Selection

According to Equation (17) of the relationship between the snow phase and snow depth, snow density ρ and incidence angle θ_i are the main input parameters; dry snow density ρ = 0.18 g/cm³ was selected as the empirical value, and the local incident angle was used to replace the satellite center incident angle. All of these factors can introduce error into snow depth inversion.

(1)

Error caused by snow density ρ

The form of snow usually changes with time, resulting in corresponding changes in snow density. Under the condition of conservation of mass, changes in snow density lead to changes in snow depth (snow compaction increases and depth decreases) [35]. Because it is affected by factors such as wind speed, altitude, and temperature, the snow density in the study area is often unevenly distributed, and using the same snow density as an input parameter is obviously insufficient. Accurate snow density is very important for snow depth inversion accuracy, and it is necessary to measure the snow density of different terrains and different underlying surfaces.

(2)

The local incident angle replacing the satellite central incident angle

Using repeated orbital interferometry to estimate snow depth, the satellite incident angle is an important parameter; however, for areas with large ground fluctuations, there is a deviation between the satellite incident angle and the local incident angle. In this experiment, the local incident angle was used to replace the satellite central incident angle to weaken the bias effects. The estimated snow depth value is relatively sensitive to local incident angle changes. Figure 17a shows the local incident angle histogram corresponding to the study area, and its values are focused between 20° and 50°. Figure 17b simulates the relationship between the local incident angle and the snow depth when the snow phase is 30 and the snow density is 0.18 g/cm³. In general, as the local incidence angle increases, the snow depth value estimated decreases.

5. Conclusions

Based on the C-band Sentinel-1 SAR data covering Jiagedaqi area of Great Xing’an Mountains, this paper uses the two-pass differential interferometry method to generate interferograms, and uses Snaphu to perform phase unwrapping in the Linux environment to obtain real phase information. A phase model is constructed based on the scattering mechanism of the target, and the forest phase is introduced to correct the snow phase deviation in forest areas, so as to less the effect of the forest on the snow phase. Then, snow depth inversion is carried out according to the relationship between snow phase and snow depth. Finally, the snow depth distribution of the study area is obtained with a spatial resolution of 30 m on 7 December 2020. The experimental data are used to analyze the error of the estimated snow depth results. The main conclusions are as follows.

The snow depth in Jiagedaqi area is mainly concentrated between 40 and 120 cm, and the average snow depth is 80.27 cm. According to the hourly snowfall observation data provided by the Xihe Energy Big Data platform, the reference snow depth value is 84.69 cm, which is converted according to the ratio between the snowfall and snow depth of 1:15, and the average snow depth is relatively close to this value. With the introduction of the forest phase, the average snow depth in forest area increases by 5.98 cm, thus demonstrating a new method for snow depth inversion in forest areas. In the future, snow depth measurements can be made directly in forest areas to further validate the effectiveness of the proposed methodology.

Snow depth inversion is limited by the coherence of interference image pairs due to factors such as changes in land cover, atmospheric effects, and snow conditions. Human activity areas, such as residential areas, may introduce interference phase due to artificial surfaces, and human activities, such as snow clearing, change the snow depth and snow distribution, resulting in large errors in the estimated snow depth results. The proposed algorithm may be not suitable for snow depth inversion in these areas. Because a synchronous observation experiment on the ground was not carried out to obtain the accurate snow content required for parameters such as snow density and snow permittivity, relevant experiments can be carried out to optimize model parameters to reduce snow depth estimation errors. In addition, the observation data of meteorological stations in forest areas are often fewer and lack representativeness, which limits the development of snow depth inversion research in forest areas.

The interaction between forest and snow cover is very complex, and the spatial distribution of the forest, canopy density, branches and leaves, tree species, and canopy structure characteristics can cause errors. For example, compared with Larix gmelinii, the evergreen coniferous forest Pinus sylvestris has higher canopy closure and less exposed ground surface in the forest area. In addition Larix gmelinii loses its leaves in winter; although the branches and trunks of the forest canopy remain, the exposed ground surface is relatively greater. Thus, the phase values contributed by the forest differ within a single pixel. The study in this paper uses the same forest phase value for the calculation, which is obviously insufficient. In future research, the classification of different types of forest area can be refined, and experiments can be carried out for different species of forest area to obtain a more accurate value of the forest phase. Jiagedaqi is divided into forest and non-forest areas according to land cover classification data, and mixed pixels containing forest and non-forest areas were not considered. Future research can be conducted for the mixed pixel problem, such as by adding a mixed pixel decomposition model, etc. In addition, one pixel on both sides of the forest and non-forest borders can be excluded, and divided into three classes of forest, non-forest, and excluded areas.

The snow depth inversion in this paper is based on the assumption that the snow cover is in a homogeneous and dry state, and radar waves can penetrate the snow layer to reach the ground surface. Strictly speaking, the media that exist are inhomogeneous and even layered in nature. That is, there is no ideal homogeneous media in nature; only at a specific incident frequency, incident angle, or media state, can some media be treated as homogeneous media. Therefore, the assumption will have an impact on the experimental results. When selecting an image date for the snow cover period, the snow state should be as close to the assumed state as possible. Alternatively, at the beginning of the experiment, the identification of wet and dry snow can be performed first. Then, the wet snow and dry snow depth inversions can be performed separately using different methods.