Estimating Snow Depth Using Multi-Source Data Fusion Based on the D-InSAR Method and 3 DVAR Fusion Algorithm

Snow depth is a general input variable in many models of agriculture, hydrology, climate, and ecology. However, there are some uncertainties in the retrieval of snow depth by remote sensing. Errors occurred in snow depth evaluation under the D-InSAR methods will affect the accuracy of snow depth inversion to a certain extent. This study proposes a scheme to estimate spatial snow depth that combines remote sensing with site observation. On the one hand, this scheme adopts the Sentinel-1 C-band of the European Space Agency (ESA), making use of the two-pass method of differential interferometry for inversion of spatial snow depth. On the other hand, the 3DVAR (three dimensional variational) fusion algorithm is used to integrate actual snow depth data of virtual stations and real-world observation stations into the snow depth inversion results. Thus, the accuracy of snow inversion will be improved. This scheme is applied in the study area of Bayanbulak Basin, which is located in the central hinterland of Tianshan Mountains in Xinjiang, China. Observation data from stations in different altitudes are selected to test the fusion method. According to the results, most of the obtained snow depth values using interferometry are lower than the observed ones. However, after the fusion using the 3DVAR algorithm, the snow depth accuracy is slightly higher than it was in the inversion results (R2 = 0.31 vs. R2 = 0.50, RMSE = 2.51 cm vs. RMSE = 1.96 cm; R2 = 0.27 vs. R2 = 0.46, RMSE = 4.04 cm vs. RMSE = 3.65 cm). When compared with the inversion results, the relative error (RE) improved by 6.97% and 3.59%, respectively. This study shows that the scheme can effectively improve the accuracy of regional snow depth estimation. Therefore, its future application is of great potential.


Introduction
Seasonal snowpack in high latitude and altitude regions in the northern and southern hemispheres is particularly vulnerable to climate change.The high albedo and low thermal conductivity of snow have strong effects on the surface energy budget and the Earth's radiative balance [1].Therefore, the existence and change of snowpack exert great significance for the energy exchange between land, atmosphere, and water, as well as carbon recycling [2,3], which influences the utilization and management of water resources directly or indirectly, including farming, forestry, and animal husbandry.Moreover, extensive and continuous existence of snowpack may cause various natural disasters, and snow depth is an important indicator of such disasters [4,5].Therefore, snow depth monitoring is an important aspect of researching near-surface processes in cold regions, providing theoretical and scientific bases for natural resource management and strategies [6].
Due to the broad application of satellite remote sensing data in ice and snow research since the 1970s, researchers worldwide have utilized different remote sensing data to conduct studies on snow depth inversion to determine the spatial distribution of snow depth.There will be access to snow depth monitoring in areas with sparse observation stations [7][8][9][10][11][12][13][14][15][16][17].However, the relationship between the snow depth and the brightness value of the optical and near infrared bands is not clear yet [7].It is not suitable for the algorithm combining different bands of visible remote sensing data with the snow depth data that is observed by the ground stations for statistical regressive analysis to obtain snow depth.Microwave remote sensing can penetrate the clouds to obtain information about snow depth, which is irreplaceable in the remote sensing of snow [8].There have been many important academic achievements in the study of using remote sensing of passive/active microwaves to estimate the snow depth/snow water equivalent [9][10][11][12].Generally, the spatial resolution obtained from the remote sensing of passive microwave is rather coarse, and it is on the level of kilometres, which cannot satisfy the requirement of snow hydrological process research at a watershed scale.The remote sensing of active microwaves has fine spatial resolution and a higher sensibility to the snow parameters.Therefore, it is more applicable in a snow study at a watershed scale [13][14][15].However, the mechanism of active microwave of remote sensing is complicated in the forward-model-based inversion process that occurs in propagation among the snowpack, the earth surface below the snow cover and surface coverings (such as vegetation).When it is applied in mountain areas where the surface condition is rather complex, the back-scattering information received by the SAR sensor is hard to distinguish.Although it is available under certain conditions, the algorithm is not widely applied.Its global applicability remains to be checked [16,17].
SAR differential interferometry can obtain altitude information about topography and surface features with a fine resolution and monitoring of ground deformation at the millimetre level [18][19][20].Such inversion methods can also be applied to measure snow depth.The essence of D-InSAR is to use the correlation between snow depth and phase information to conduct snow depth inversion in SAR image pairs.To study the regional snow depth distribution and changes, researchers have also utilized SAR differential interferometry and selected active microwave data, such as the L-band in the Tandem X-band in TerraSAR.These studies have found out that SAR data from different bands have different sensibilities to the snow depth.Decorrelation of InSAR is influenced by snow conditions, such as snow humidity, melting, and refreezing [21][22][23][24][25].Moreover, InSAR's decoherence is influenced by the spatial and temporal baselines.The decoherence of the C-band becomes obvious when the timeframe is longer than six months.However, the decoherence of the L-band only becomes clear when the timeframe is longer than two years [26].In snow remote sensing research, the C-band of SAR is a major tool for map drawing and snow areas monitoring, due to its broad data sources and short revisiting period [27].Seasonal snow is usually unstable and exhibits a rapid change rate.The C-band has the advantage of forming image pairs with a short temporal baseline, which is conducive to the continuous monitoring of snow depth.
In spatial snow depth monitoring, remote sensing as a primary technology is able to reflect the spatial distribution of snow depth rapidly and comprehensively.However, the instantaneous characteristic of remote sensing is easily affected by the atmospheric conditions, ground heterology, and inversion algorithm.Thus, the results are not accurate enough.A data fusion algorithm can effectively combine the remote sensing inversion results with observations, so the accuracy of snow depth inversion can be improved so as to obtain a more accurate snow depth [28].The data fusion methods for integrating the two types of data include filtering algorithms and variation algorithms.Variation fusion is a modern method of data fusion that includes three-and four-dimensional variation.At present, the major weather forecasting centers around the world universally adopt these technologies, solving the variation minimization problem to obtain the optimal analysis solution and attain more reliable forecasting results.In recent years, researchers worldwide have worked hard to generate a better application of the three-dimensional variation fusion algorithm and to expand the applicable field of this technology [29,30].In related studies, researchers have used three dimensional variational (3DVAR) to assimilate observation information (temperature, humidity, intensity of pressure, wind speed, etc.) and satellite-borne detected data to improve the accuracy of background field generated by reanalysis data in the field of weather research and forecast (WRF) model to increase the reliability of meteorological forecasting results [31,32].
Therefore, this study proposes a scheme to estimate spatial snow depth in order to obtain the optimal estimation of the snow depth with a higher accuracy than shown in previous methods, which combines remote sensing with in situ data.This scheme carries out snow depth inversion by selecting the ESA's Sentinel-1 microwave remote sensing data (C-band) with the InSAR two-pass method and combining data from ground observation stations to conduct data fusion in the study area of Bayanbulak Basin in the central segment of Tianshan Mountains, Xinjiang, China, thereby it would provide foundation for future assessment and water resource management under climate changes.

Study Area
The Bayanbulak Basin is located in the central segment of Tianshan Mountains.It is surrounded by mountains on three sides with semi-closed geographic characteristics.It possesses features of a high and cold basin.Bayanbulak Basin includes the Big and Small Yultuz Basins.The Small Yultuz Basin, with a total area of 6417.96km 2 (Figure 1), is taken as the study area in this study.The flat-bottom basin altitude is higher than 2400 m, but lower than 2600 m.It is surrounded by large mountains with altitudes in the range of 2458-4685 m.There are also certain areas with permanent snow and glaciers.Due to the basin's unique cold, high-altitude climate, as well as its distinct landform, various alpine steppes and meadow ecosystems are present (the landcover types in the study area is a simple underlying surface structure).There are large areas of wet meadow and lake, with rich aquatic plants and animals.This forms a distinct inland wetlands ecosystem that plays a significant role in water yield adjustment, water storage, and the maintenance of the regional water balance [33].
Remote Sens. 2017, 9, 1195 3 of 16 variational (3DVAR) to assimilate observation information (temperature, humidity, intensity of pressure, wind speed, etc.) and satellite-borne detected data to improve the accuracy of background field generated by reanalysis data in the field of weather research and forecast (WRF) model to increase the reliability of meteorological forecasting results [31,32].Therefore, this study proposes a scheme to estimate spatial snow depth in order to obtain the optimal estimation of the snow depth with a higher accuracy than shown in previous methods, which combines remote sensing with in situ data.This scheme carries out snow depth inversion by selecting the ESA's Sentinel-1 microwave remote sensing data (C-band) with the InSAR two-pass method and combining data from ground observation stations to conduct data fusion in the study area of Bayanbulak Basin in the central segment of Tianshan Mountains, Xinjiang, China, thereby it would provide foundation for future assessment and water resource management under climate changes.

Study Area
The Bayanbulak Basin is located in the central segment of Tianshan Mountains.It is surrounded by mountains on three sides with semi-closed geographic characteristics.It possesses features of a high and cold basin.Bayanbulak Basin includes the Big and Small Yultuz Basins.The Small Yultuz Basin, with a total area of 6417.96km 2 (Figure 1), is taken as the study area in this study.The flatbottom basin altitude is higher than 2400 m, but lower than 2600 m.It is surrounded by large mountains with altitudes in the range of 2458-4685 m.There are also certain areas with permanent snow and glaciers.Due to the basin's unique cold, high-altitude climate, as well as its distinct landform, various alpine steppes and meadow ecosystems are present (the landcover types in the study area is a simple underlying surface structure).There are large areas of wet meadow and lake, with rich aquatic plants and animals.This forms a distinct inland wetlands ecosystem that plays a significant role in water yield adjustment, water storage, and the maintenance of the regional water balance [33].
The high elevation/cold climate is dominant in the Bayanbulak Basin.In the last 52 years, the records indicate that the winter lasts up to seven months, and that the annual average temperature is −4.6 °C, with a minimum temperature of −48 °C.There are 139.3snow-covered days, and the maximum depth of frozen ground is 750 cm.Besides, an increasing trend in the annual precipitation and rainy days in the Bayanbulak Basin has occurred, with values of 9.5 mm/10 a and 3.2 day/10 a, respectively.However, the increase of annual precipitation in cold seasons is more than that in warm seasons [34].The high elevation/cold climate is dominant in the Bayanbulak Basin.In the last 52 years, the records indicate that the winter lasts up to seven months, and that the annual average temperature is −4.6 • C, with a minimum temperature of −48 • C.There are 139.3snow-covered days, and the maximum depth of frozen ground is 750 cm.Besides, an increasing trend in the annual precipitation and rainy days in the Bayanbulak Basin has occurred, with values of 9.5 mm/10 a and 3.2 day/10 a, respectively.However, the increase of annual precipitation in cold seasons is more than that in warm seasons [34].
The Kaidu River, which originates from the Bayanbulak Basin, is a river of great importance for the Xinjiang economy.It flows into the largest inland freshwater lake in China, Bosten Lake.As the source of the Tarim River, Bosten Lake is a crucial component in a north-to-south water diversion project.It is responsible for delivering emergency water to the downstream ecosystem of the Tarim River.Therefore, the Bayanbulak Basin is a vital water source in the oasis of Southern Xinjiang.

Field Data Collection
Data of snow depth and snow density for this study were obtained from 12 temporal observation stations (set up by the Xinjiang Institute of Ecology and Geography of the Chinese Academy of Sciences; Figure 1) in the Bayanbulak Basin in the central segment of Tianshan Mountains, Xinjiang.The layout of observation station is evenly distributed as much as possible that taking factors like differences in terrain and meteorological into consideration.In addition, gradient observation is carried out in the north-south direction, forming a fan-shaped observation network.
The snow depth monitoring devices are based on the principle of ultrasonic ranging to obtain snow depth, automatically.The data acquisition frequency is twice per day, and the data acquisition times occur at 10:00 a.m. and 20:00 p.m. (UTC+8).Besides, snow density is calculated according to the snow depth that measured and the weight of the snowpack above the snow pillow.The data of snow depth and snow density for 18 December 2016 and 11 January 2017 were used for this study.All of these data passed a quality control (3σ principle).
Data of six observation stations (snow-depth site 1, 3, 4, 7, 9, and 11) are conducted with 3DVAR fusion algorithm to obtain optimal snow depth distribution in Bayanbulak Basin, while the data of the rest six stations are used as verification data for the analysis of snow depth fusion result.

Sentinel-1 SAR Images and Preprocessing
As the successor satellite to ERS-2 and Envisat, Sentinel-1 is an important component of Global Monitoring for Environment and Security (GMES), which is supported by the European Commission and European Space Agency (ESA).It consists of Sentinel-1A and Sentinel-1B.Sentinel-1A was launched on 3 April 2014, while Sentinel-1B was launched on 25 April 2016.Sentinel-1 can provide SAR C-band data, which can be repeatedly observed.The revisiting period of a single satellite is 12 days.When the two satellites form a network, the revisiting period is six days, with free data-sharing services.
Sentinel-1 carries a C-band SAR sensor with a frequency of 5.405 GHz.There are four types of acquisition modes: stripmap mode (SM, 80 km Swath, 5 m × 5 m spatial resolution), interferometric wide swath (IW, 250 km Swath, 5 m × 20 m spatial resolution), extra-wide swath (EW, 400 km Swath, 25 m × 100 m spatial resolution), and wave mode (WM, 20 km × 20 km, 5 m × 20 m spatial resolution).Level-1 also has two types of data products, namely Single Look Complex (SLC) and Ground Range Detected (GRD).The data mode employed in this study was IW SLC, with a wide swath of 250 km and a spatial resolution of 5 m × 20 m.
The objective of this study was to retrieve dry snow depth (the minimum temperature in Bayanbulak Basin on 18 December was −26 • C, and the maximum temperature was −17 • C; the minimum temperature in Bayanbulak Basin on 11 January was −38 • C, and the maximum temperature was −11 • C. Therefore, the snow was assumed to be dry snow under these conditions).Meanwhile, the temporal baseline and space baseline in the two-pass method were taken into consideration.What is more, the two images acquired by Sentinel-1B and the corresponding accurate orbital data AUX_POEORB (precise orbital ephemeris parameters that were generated 21 days later after acquiring the data; the precision was within 5 cm) were also included.Table 1 shows the details of the selected images.
The two images selected in this study were imaged by Sentinel-1B on the descent orbit from the right lateral view over the study area, and the images were side-inverted.The Sentinel Application Platform (SNAP, a common architecture for all Sentinel Toolboxes that was jointly developed by Brockmann Consult and Array Systems Computing) was adopted for preprocessing.The preprocessing steps are described below.
(1) Precise matching of the master and slave images.Orbital parameters were used to roughly determine the deviation between the master and slave images.For a precise matching of image pairs, the R-D model that was developed by Brown [35] and earth model equation were used.
After the calculation of the matched model, the real and imaginary parts of image pairs were resampled to calculate the interferogram.The bilinear interpolation method was adopted for this process.(2) Generation of the interferogram.After the resampling of the real and imaginary parts of the image pairs, the master image and slave image were subjected to conjugate multiplication.Then, the interferometric phase images were generated.(3) Removal of flat earth effect and terrain height effect.Since the flat earth effect and terrain height effect are contained in the InSAR interferogram, they needed to be removed before phase unwrapping to obtain the phase information concerning snow depth [36].(4) Removal of speckle noise.The Goldstein algorithm was used in this study for noise reduction in interferogram.Then, the Boxcar filter (small 3 × 3 window) was used to generate a further smoothing effect of this interferogram [37,38].(5) Phase unwrapping.Phase unwrapping is one of the most crucial steps in InSAR.This procedure directly determines the accuracy of snow depth inversion [39].Minimum cost flow (MCF) was adopted for the unwrapping (Figure 2).

Processing of the Moderate Resolution Imaging Spectrometer (MODIS) Snow Cover Data and Virtual Site Arrangement
The Moderate Resolution Imaging Spectrometer (MODIS) snow cover data can be freely downloaded from the United States National Snow and Ice Data Center (NSIDC).It includes the eight-day composite snow cover data MOD10A2 and MYD10A2 of the Terra satellite (covering the area in the morning) and Aqua satellite (covering the area in the afternoon) at 500-m resolution.Based on the longitude and latitude of the typical area in Bayanbulak Basin (Figure 1), the image with the tracking number h24v04 was sufficient to cover the entire study area.MOD10A2 and MYD10A2 on 18 December 2016 and 11 January 2017 were selected for this study, and the MODIS reprojection tool (MRT) was used for data extraction, reprojection, and cropping based on the extent of the study area.
There is a specific cloud detection algorithm in MOD10A2 and MYD10A2 [40][41][42].Due to the limitations of the cloud removing algorithm, the SCA data in the study was gathered from adjacent pixels in a relatively short contiguous time to reflect the largest range of snow cover within eight days, which compensates the low sampling density.In addition, the establishment of virtual sites in the study area can reduce the smoothness of the interpolation results from cokriging [43].Figure 3a,b show the snow cover for the study area on 18 December 2016 and 11 January 2017, after eliminating clouds.
(4) Removal of speckle noise.The Goldstein algorithm was used in this study for noise reduction in interferogram.Then, the Boxcar filter (small 3 × 3 window) was used to generate a further smoothing effect of this interferogram [37,38].(5) Phase unwrapping.Phase unwrapping is one of the most crucial steps in InSAR.This procedure directly determines the accuracy of snow depth inversion [39].Minimum cost flow (MCF) was adopted for the unwrapping (Figure 2).

Processing of the Moderate Resolution Imaging Spectrometer (MODIS) Snow Cover Data and Virtual Site Arrangement
The Moderate Resolution Imaging Spectrometer (MODIS) snow cover data can be freely downloaded from the United States National Snow and Ice Data Center (NSIDC).It includes the eight-day composite snow cover data MOD10A2 and MYD10A2 of the Terra satellite (covering the area in the morning) and Aqua satellite (covering the area in the afternoon) at 500-m resolution.Based on the longitude and latitude of the typical area in Bayanbulak Basin (Figure 1), the image with the tracking number h24v04 was sufficient to cover the entire study area.MOD10A2 and MYD10A2 on 18 December 2016 and 11 January 2017 were selected for this study, and the MODIS reprojection tool (MRT) was used for data extraction, reprojection, and cropping based on the extent of the study area.
There is a specific cloud detection algorithm in MOD10A2 and MYD10A2 [40][41][42].Due to the limitations of the cloud removing algorithm, the SCA data in the study was gathered from adjacent pixels in a relatively short contiguous time to reflect the largest range of snow cover within eight days, which compensates the low sampling density.In addition, the establishment of virtual sites in the study area can reduce the smoothness of the interpolation results from cokriging [43].Figure 3a,b show the snow cover for the study area on 18 December 2016 and 11 January 2017, after eliminating clouds.
We analyze the snow-free pixels in the rectified MODIS snow cover data and set up virtual observation stations, using the following rules [43]: (1) There is no snow in the eight neighbouring pixels of the virtual station.
(2) The distance between the observation stations could be no less than the smallest distance between the existing stations after adding the virtual stations to achieve an even distribution and density of observation stations.

Snow Depth Inversion Algorithm
Viewed as an apparent and refractive medium, the scattering effect of snow can be ignored because the size of ice particles in the snow is less than 1 mm and the microwave length is longer than 1.0 cm.Since the refractive index of the snow-soil boundary exhibits a great difference (microwave's refractive index is εr = 1.0-1.8 in dry snow, εr = 2.5 in dry soil, and εr ≈ 40.0 in wet soil), the refractive effect in the snow-soil boundary could be ignored [23].We analyze the snow-free pixels in the rectified MODIS snow cover data and set up virtual observation stations, using the following rules [43]: (1) There is no snow in the eight neighbouring pixels of the virtual station.
(2) The distance between the observation stations could be no less than the smallest distance between the existing stations after adding the virtual stations to achieve an even distribution and density of observation stations.

Snow Depth Inversion Algorithm
Viewed as an apparent and refractive medium, the scattering effect of snow can be ignored because the size of ice particles in the snow is less than 1 mm and the microwave length is longer than 1.0 cm.Since the refractive index of the snow-soil boundary exhibits a great difference (microwave's refractive index is ε r = 1.0-1.8 in dry snow, ε r = 2.5 in dry soil, and ε r ≈ 40.0 in wet soil), the refractive effect in the snow-soil boundary could be ignored [23].
The underlying surface of study area is simple, and there is no forest canopy to intercept the snow.Thus, the slight phase caused by underlying floor is also ignored.The total phase in the study area includes [44]: In this equation, ϕ atm is the phase that is caused by the atmosphere (mainly caused by the delay of troposphere and ionised layer); ϕ orb is the phase caused by orbital deviation; ϕ noise is the phase caused by speckle noise (including thermal noise, sampling, and matching errors, etc.); ϕ flat is the phase that is caused by the flat earth effect; ϕ topo is the phase caused by the terrain height effect; and, ϕ snow is the phased caused by snow depth.
According to the geometrical relationship between snow phase and the travel path of the radar beam, the calculation equation for phase ϕ snow is [25,44]: Here, λ is the microwave length (cm) of the C-band of Sentinel-1; θ i is the satellite incidence angle ( • ); d s is the snow depth (cm); and, ε s is the dielectric constant of snow (ε 0 ).Moreover, ε s is closely related to the snow density, and the correlation equation is as follows [25]: where ρ is the snow density (g/cm 3 ).Thus, ϕ snow can be calculated by selecting master and slave images for interferometry and the removal of other phase information (this study ignores the atmosphere delay effect of the troposphere and ionised layer on the C-band).

The Cokriging Model
As an interpolation technique, cokriging allows for multiple variables correlated with explanatory variables, which optimize the estimated result [43].Target variables that are difficult to observe can be locally estimated by auxiliary variables that are easy to observe and control through establishing the cross-covariance function: Here, X k (x 0 ) is the selected auxiliary variable; p indicates the number of selected auxiliary variables; β is the regression coefficient; and, ε(x 0 ) is the residual.The spatial variation of snow depth is complicated (it can vary in different directions or with different scales).Therefore, the following auxiliary variables are available to choose from [45]: Northness is defined as the cosine of the aspect and has values ranging from −1 to 1.A value of −1 represents a slope that is facing directly south and a value of 1 indicates a slope facing directly north.
Eastness (Figure 3b) is defined as the sine of the aspect with values ranging from −1 to 1. −1 represents an east-facing slope and 1 indicates a west-facing slope.
Topographic ruggedness index (TRI) is defined as the standard deviation of the elevation in a 270 m × 270 m (3 × 3 pixels) window.High TRI values are indicative of areas with a large fluctuation in elevation.
Additionally, auxiliary variables such as height and slope are available.

The 3DVAR Land Data Assimilation System Scheme
This study introduces the 3DVAR method to estimate snow depth while minimizing the deviation between the inversion and actual condition.The method consists of three parts, as described below.

Cost Function
The basic concept of the 3DVAR fusion is to find the minimum of a quadratic function (called the cost function) that presents the difference of analysis from the observation and background.In this data fusion scheme, the cost function is taken as where x and x b are analysis and background vectors of length of N, while y is an observation vector with length M. B is an N × N background error covariance matrix, while O is an M × M observation error covariance matrix.R is an observation operator, which maps analysis variables into observation variables and observation positions.Under the condition that observation variables coincide with model variables, R is a simple linear interpolation operator.Otherwise, R is a complex nonlinear operator.

The Background and Observation Error Model
In essence, the two parts of the cost function is a weighted linear combination with an analysis of background and observation.The inverses of B and O (B −1 and O −1 ) are regarded as weightings of the background and observation that contributed in the analysis.Clearly, the result of fusion can be attributed to the background error covariance B and observation error covariance O [46].
The snow depth inversion result determined by D-InSAR reflects the spatial distribution of snow depth (acts as background), while the actual points (including virtual and observation stations) provide spot information about snow depth.In order to acquire background error covariance, errors between inversion and actuality (including both virtual and observation stations) are obtained point by point.Then, cokriging is used to obtain the background error.The background error covariance is obtained accordingly.
By conducting spatial interpolation for all of the actual points, the observation field in the 3DVAR fusion model is obtained.Moreover, errors between the observation field and actuality would be considered point by point.Finally, cokriging is used to obtain the observation error.The observation error covariance is obtained accordingly.

The Principle of 3DVAR and the Cost Function Minimization Algorithm
For the solution of the minimal value of formula (5), there are some algorithms (conjugate gradient methods, quasi-newton method, etc.) that can be used to obtain an optimal solution [46].This study adopts the method of optimal linear unbiased estimation: the gradient of cost function for the control variables is 0, and the analytical solution is as follows:

Inversion Results for Snow Depth
From Equations ( 2) and (3), it is clear that the regional snow depth distribution can be obtained from the snow density ρ and incidence angle θ i .In this study area, years of observations demonstrated that: the snow density varied between 0.15 g /cm 3 and 0.25 g/cm 3 from accumulation period to the deterioration of snow before the melting period [25,47].Therefore, this study assumes that the snow density is 0.18 g/cm 3 , and the local incident angle replaces the central incidence angle in Equation (2) (Figure 4).Decorrelation sources, such as the thermal noise, matching errors, and perpendicular baseline, will cause decorrelation.Errors are inevitable in the process of phase removal that includes the elimination of the orbit effect, flat earth effect, and terrain effect [48].Thus, correction of the unwrapping phase is necessary.Since the snow depth phase of the snow-free area is the minimal value of the phase in interferogram, it is set as the base to retrieve snow depth.The inversion of snow depth with a spatial resolution of 13.89 m can be obtained (Figure 5(a1,a2)).
It can be seen from Figure 5(a1,a2) that the snow depth distribution were 0.00-36.00cm and 0.00-32.39cm, respectively.In Figure 5(a1), snow depth gradually increased from marginal mountain to centre mountain (the snow depth of the northeast basin is deeper than that of the other regions).The snow depth of the southern slopes is generally thinner than that of the northern slope.While in Figure 5(a2), the snow depth of the northwest basin is deeper than that of the other regions.Decorrelation sources, such as the thermal noise, matching errors, and perpendicular baseline, will cause decorrelation.Errors are inevitable in the process of phase removal that includes the elimination of the orbit effect, flat earth effect, and terrain effect [48].Thus, correction of the unwrapping phase is necessary.Since the snow depth phase of the snow-free area is the minimal value of the phase in interferogram, it is set as the base to retrieve snow depth.The inversion of snow depth with a spatial resolution of 13.89 m can be obtained (Figure 5(a1,a2)).
It can be seen from Figure 5(a1,a2) that the snow depth distribution were 0.00-36.00cm and 0.00-32.39cm, respectively.In Figure 5(a1), snow depth gradually increased from marginal mountain to centre mountain (the snow depth of the northeast basin is deeper than that of the other regions).The snow depth of the southern slopes is generally thinner than that of the northern slope.While in Figure 5(a2), the snow depth of the northwest basin is deeper than that of the other regions.

Background Error Covariance and Observed Error Covariance in the 3DVAR Fusion Algorithm
Since the snow depth data from the actual points were more reliable than the inversion results, actual snow depth was adopted to correct the inversion results.
depth with a spatial resolution of 13.89 m can be obtained (Figure 5(a1,a2)).
It can be seen from Figure 5(a1,a2) that the snow depth distribution were 0.00-36.00cm and 0.00-32.39cm, respectively.In Figure 5(a1), snow depth gradually increased from marginal mountain to centre mountain (the snow depth of the northeast basin is deeper than that of the other regions).The snow depth of the southern slopes is generally thinner than that of the northern slope.While in Figure 5(a2), the snow depth of the northwest basin is deeper than that of the other regions.Based on the above analysis, the terrain ruggedness index was used as the covariate.While the snow depth inversion errors in actual points were employed as the master variable, the semi-variance function was employed to analyze the variation and structure of the spatial distribution of snow depth inversion error in the study area.Then, the circular model was chosen as the semi-variance model to obtain the background error concerning snow depth.Thus, the background error covariance was obtained on this foundation, respectively.

Establishment of the Observation Error Covariance in the 3DVAR Fusion Algorithm
Obtained by the cokriging, the observation error covariance in the 3DVAR fusion model is based on the observation error of snow depth in the study area.Therefore, height-the factor that is significantly correlated with the snow depth (r 1 = 0.64, p 1 = 0.003; r 2 = 0.49, p 2 = 0.023)-was used as covariate.To obtain the spatial distribution of snow depth, the J-Bessel model was selected for semi-variance in the cokriging method.On this basis, the point-by-point observation errors between the observation field and actuality were used as the master variable, and height-the factor that was significantly correlated with snow depth inversion error (r 1 = 0.64, p 1 = 0.003; r 2 = 0.57, p 2 = 0.011)-was used as the covariate.The spherical model was selected as the cokriging of the semi-variance model to determine the distribution of the observation error in the study area.The observation error covariance was obtained on this foundation, respectively.

Estimation of Snow Depth Based on the 3DVAR Fusion Algorithm
Based on 3DVAR fusion algorithm in ( 5) and ( 6), Figure 5(b1,b2) illustrates the estimated results of the spatial distribution of snow depth.The snow depth was in the range of 0.00-34.81cm and 0.00-29.51cm, respectively.
The snow depth distribution between 3DVAR fusion and inversion were compared.In the estimated results of 3DVAR, the area with a snow depth ranging from 0.00 to 5.00 cm increased by 0.34% and 6.76%.What is more, the area with a snow depth in the range of 5.00-10.00cm increased by 0.88% and 13.21%, respectively, which indicated that the values of low snow depth areas increased after the 3DVAR fusion method was applied in correction.

Evaluation of Accuracy
Figure 6 illustrates the comparison between the actual snow depth obtained from the observation stations, inversion, and 3DVAR fusion snow depth.It is evident that most of snow depth obtained by the two-pass method of D-InSAR method were lower than the actual data.However, the inversion result changed after the 3DVAR fusion: Figure 6a is the one of a 3DVAR fusion result in this study.The snow depth values for observation points 2, 8, and 12 were higher than those from the 3DVAR fusion.The corrected snow depth values for observation point 10 were lower than the inversion data.Moreover, the value of R 2 with 0.31 increased to 0.50.In addition, Figure 6b is the other one in this study, and the snow depth values for observation points 2, 5, 10, and 12 were higher than those from the 3DVAR fusion, the corrected snow depth values for observation points 5, 6, 10, and 12 lower than the inversion data.Moreover, the value of R 2 with 0.27 increased to 0.46.To illustrate the effectiveness of the 3DVAR fusion method for integrating actuality and inversion to form the optimized snow depth estimation, the 3DVAR fusion results of the snow depth at 6 observation points were tested using the following statistical methods: relative error (RE), the mean estimation error (MEE), root mean square error (RMSE), and mean absolute estimation error (MAEE).The results are shown in Table 2.It can be seen from Table 2 that the 3DVAR fusion results of snow depth were slightly higher To illustrate the effectiveness of the 3DVAR fusion method for integrating actuality and inversion to form the optimized snow depth estimation, the 3DVAR fusion results of the snow depth at 6 observation points were tested using the following statistical methods: relative error (RE), the mean estimation error (MEE), root mean square error (RMSE), and mean absolute estimation error (MAEE).The results are shown in Table 2.
It can be seen from Table 2 that the 3DVAR fusion results of snow depth were slightly higher than the inversion results (R 2 = 0.31 vs. R 2 = 0.50, RMSE = 2.51 cm vs. RMSE = 1.96 cm; R 2 = 0.27 vs. R 2 = 0.46, RMSE = 4.04 cm vs. RMSE = 3.65 cm).In addition, when compared with the inversion results without fusion, the accuracy tests of 3DVAR fusion, including RE, MAEE, and RMSE, improve with the minimum of 3.59%, 12.22%, and 0.97%.It can be concluded from these accuracy tests that 3DVAR fusion can be used to improve the estimation accuracy of the spatial distribution of regional snow depth.The inversion results of the spatial distribution of regional snow depth in this study were acquired using the two-pass method of D-InSAR.Correlation analysis was conducted between the snow depth inversion error of observation points and the height, TRI, eastness, northness, and slope.It is found that snow depth inversion error was significantly correlated with TRI only (r 1 = 0.82, p 1 = 0.00; r 2 = 0.73, p 2 = 0.01).It shows that while TRI increases, the radar echoes from earth surface belowed the snow cover decrease correspondingly.Therefore, this reduction renders a positive bias in the snow phase, which causes underestimated snow depth.
There are numerous studies of applying active microwave remote sensing into snow depth inversion, and researchers have used the VV/HH of the X-band, L-band and C-band via polarization phase differences to conduct snow depth inversion [15,[49][50][51].However, the snow depth in these studies is measured by phase differences in radar time series.Since several parameters are needed in the snow depth inversion, the principles in these studies are complicated.The two-pass method in interferometry takes advantage of the geometrical relationship between the snow phase and the travel path of the radar beams to accomplish snow depth inversion, which is advantageous to theoretical operation.
In future studies, InSAR of different polarization modes, such as VH, HH, and HV, could be conducted.Meanwhile, simultaneous observations on snow properties, such as the snow density, structure of the snow layer, and air humidity should also be launched.Crucial terrain parameters, including slope and terrain roughness, can be combined to optimize models of the snow phase.Thereby the accuracy of snow depth inversion is improved.

Limitations of D-InSAR Method
Although there are many advantages of the two-pass method (e.g., high spatial resolution), it also has limitations.One of the major limitations is that the coherence between corresponding pixels of SAR images is significantly decreased due to the decorrelation of temporal and spatial data.As results, observation points 7, 9, 10, and 11 in the study area showed a low coherence (γ < 0.25; Figure 7), and the inversion error at these observation points was higher than 35%, indicating that the accuracy of snow depth inversion is also limited by the coherence.
In this study, the D-InSAR method assumed that the snowpack was in dry conditions as only dry snow can be viewed as a transparent media layer.In actual, SAR is extremely sensitive to the internal conditions of the snow layer.After repeated melting and freezing, the distribution of snow layers and the snow density cannot be evenly separated.Figure 8 shows the observed snow density in six observation stations (snow-depth site 2, 5, 6, 8, 10, and 12).It can be seen that the snow density of station whose actual snow depth higher than inversion ones is less than 0.18 g/cm 3 .While, the snow density of station whose actual snow depth lower than inversion ones is more than 0.18 g/cm 3 .This study relied on the statistical relationship between the dielectric constant and snow density, which will cause errors.It is difficult to obtain real-time information about the snow density, dielectric coefficient, and snow layers in the study area.
including slope and terrain roughness, can be combined to optimize models of the snow phase.Thereby the accuracy of snow depth inversion is improved.

Limitations of D-InSAR Method
Although there are many advantages of the two-pass method (e.g., high spatial resolution), it also has limitations.One of the major limitations is that the coherence between corresponding pixels of SAR images is significantly decreased due to the decorrelation of temporal and spatial data.As results, observation points 7, 9, 10, and 11 in the study area showed a low coherence (γ < 0.25; Figure 7), and the inversion error at these observation points was higher than 35%, indicating that the accuracy of snow depth inversion is also limited by the coherence.In this study, the D-InSAR method assumed that the snowpack was in dry conditions as only dry snow can be viewed as a transparent media layer.In actual, SAR is extremely sensitive to the internal conditions of the snow layer.After repeated melting and freezing, the distribution of snow layers and the snow density cannot be evenly separated.Figure 8 shows the observed snow density in six observation stations (snow-depth site 2, 5, 6, 8, 10, and 12).It can be seen that the snow density of station whose actual snow depth higher than inversion ones is less than 0.18 g/cm 3 .While, the snow density of station whose actual snow depth lower than inversion ones is more than 0.18 g/cm 3 .This study relied on the statistical relationship between the dielectric constant and snow density, which will cause errors.It is difficult to obtain real-time information about the snow density, dielectric coefficient, and snow layers in the study area.It should be noted that the total phase that was measured in this study included the atmosphere phase (ϕatm), as well as the noise phase, which causes errors in the inversion results.Moreover, errors were caused by water vapor content, and detailed observation about precipitable water vapor in the atmosphere during the same period is not sufficient.The inversion did not delete the corresponding atmosphere phase and noise phase from the remained phase, which may have been one of the reasons for the snow depth inversion errors.

Limitations Related to the Establishment of Background Error and Observation Error in the 3DVAR Fusion Algorithm
Background error is a key factor in data fusion.Thus, the accuracy of the background error should be quantified.However, the actual snow depth distribution is unknown, and we cannot obtain the real background error.We can only make use of geo-statistical methods to estimate the background error based on certain assumptions.Generally, the selection of a statistical approach for background error is based on an observation network and the prior knowledge of errors.Cokriging is an easy approach for the calculation of the background error, however, it relies greatly on the quality of semi-variance simulation.Thus, there are some limitations in it.
This study applied spatial interpolation of the samples of snow depth to obtain the spatial distribution of snow depth to gain the spatial distribution of snow depth observation errors using cokriging (as a representative interpolation algorithm).Such an approach makes good use of the related structures and provides error estimation of the interpolation results.However, it fails to fully It should be noted that the total phase that was measured in this study included the atmosphere phase (ϕ atm ), as well as the noise phase, which causes errors in the inversion results.Moreover, errors were caused by water vapor content, and detailed observation about precipitable water vapor in the atmosphere during the same period is not sufficient.The inversion did not delete the corresponding atmosphere phase and noise phase from the remained phase, which may have been one of the reasons for the snow depth inversion errors.

Limitations Related to the Establishment of Background Error and Observation Error in the 3DVAR Fusion Algorithm
Background error is a key factor in data fusion.Thus, the accuracy of the background error should be quantified.However, the actual snow depth distribution is unknown, and we cannot obtain the real background error.We can only make use of geo-statistical methods to estimate the background error based on certain assumptions.Generally, the selection of a statistical approach for background error is based on an observation network and the prior knowledge of errors.Cokriging is an easy approach for the calculation of the background error, however, it relies greatly on the quality of semi-variance simulation.Thus, there are some limitations in it.
This study applied spatial interpolation of the samples of snow depth to obtain the spatial distribution of snow depth to gain the spatial distribution of snow depth observation errors using cokriging (as a representative interpolation algorithm).Such an approach makes good use of the related structures and provides error estimation of the interpolation results.However, it fails to fully consider the environmental factors that are closely correlated with snow depth.What is more, the accuracy of interpolation is influenced by the sampling density.Therefore, the interpolation accuracy may be further improved if the factors concerning observation error can be fully analyzed more suitable environmental factors are selected as covariate to analyze the spatial variation of snow depth's observation error.

Limitations of the Data Fusion Methods
The 3DVAR fusion method that is used in this study can obtain the minimal vector satisfying the cost function by introducing observation error covariance and background error covariance, that is, by using observation data to correct the snow depth inversion result.In contrast to the least-squares method, the analysis solution is acquired by a global minimum of the cost function for the background and observation fields.Therefore, its constraint features avoid unstable results from snow depth correction or abnormal values.
Although this method was easy to execute when compared with sequential fusion methods (such as the ensemble Kalman filter), the background error needed to be based on geostatistical analysis and did not change along with time sequence (i.e., it lacked time sequence updating of error covariance).Therefore, the accuracy of estimation may be further improved if we can combine the updating structure of error covariance in the Kalman filter and the statistics of snow depth estimation to optimize and adjust the error covariance.

Conclusions
In this study, we proposed a snow depth spatial distribution estimation scheme that used a 3DVAR fusion method that is derived from a numerical meteorological model to combine remote sensing inversion and observation.This study employed the Sentinel-1 SAR data to form interferometry image pairs for snow depth inversion using the two-pass method.On this basis, we also made use of the 3DVAR fusion method to integrate actuality to correct snow inversion errors.The experiment was carried out in the Bayanbulak Basin in central segment of Tianshan Mountains.The major conclusions are as follows: (1) After correct preprocessing of Sentinel-1 data, the SAR differential interferometry technology of the two-pass method can be utilized for snow depth inversion.However, the variance of certain factors, such as coherence, may influence the accuracy of snow depth inversion.The inversion in this study, which was based on D-InSAR, is a face-to-face process.Therefore, it is more difficult to guarantee a high coherence for a whole region than it is for limited points.It is possible to make use of the SAR differential interferometry technology of the two-pass method to select limited points in areas of high coherence to acquire snow depth.(2) The study introduced the background and observation error covariance to obtain the minimal analysis vector that would satisfy the cost function, that is, by making use of observation data to correct snow depth inversion.The study indicated that after the 3DVAR fusion, every accuracy indicator of snow depth inversion had been improved.

Figure 1 .
Figure 1.Sketch map of the study area of Bayanbulak.Figure 1. Sketch map of the study area of Bayanbulak.

Figure 1 .
Figure 1.Sketch map of the study area of Bayanbulak.Figure 1. Sketch map of the study area of Bayanbulak.

Figure 3 .
Figure 3. Map of snow cover distribution in a typical region of the Bayanbulak Basin on 18 December 2016 (a) and 11 January 2017 (b).

Figure 3 .
Figure 3. Map of snow cover distribution in a typical region of the Bayanbulak Basin on 18 December 2016 (a) and 11 January 2017 (b).

Figure 4 .
Figure 4. Images of the local incidence angle on 18 December 2016 (a) and 11 January 2017 (b).

Figure 4 .
Figure 4. Images of the local incidence angle on 18 December 2016 (a) and 11 January 2017 (b).
Remote Sens. 2017,9, 1195   11 of 16than those from the 3DVAR fusion, the corrected snow depth values for observation points 5, 6, 10, and 12 lower than the inversion data.Moreover, the value of R 2 with 0.27 increased to 0.46.

Figure 6 .
Figure 6.Comparison between actual snow depth, inversion, and the 3DVAR fusion snow depth on 18 December 2016 (a) and 11 January 2017 (b).

Figure 6 .
Figure 6.Comparison between actual snow depth, inversion, and the 3DVAR fusion snow depth on 18 December 2016 (a) and 11 January 2017 (b).

Figure 8 .
Figure 8.The verified snow density in observation stations on 18 December 2016 (a) and 11 January 2017 (b).

Table 1 .
List of Sentinel-1B SAR Images Acquired Over the Study Area.

Table 1 .
List of Sentinel-1B SAR Images Acquired Over the Study Area.

Table 2 .
Comparison of Precision Evaluation Indices between Inversion and 3DVAR Fusion.

Table 2 .
Comparison of Precision Evaluation Indices between Inversion and 3DVAR Fusion.