Assimilation of D-InSAR snow depth data by an ensemble Kalman filter

Snow depth mirrors regional climate change and is a vital parameter for medium- and long-term numerical climate prediction, numerical simulation of land-surface hydrological process, and water resource assessment. However, the quality of the available snow depth products retrieved from remote sensing is inevitably affected by cloud and mountain shadow, and the spatiotemporal resolution of the snow depth data cannot meet the need of hydrological research and decision-making assistance. Therefore, a method to enhance the accuracy of snow depth data is urgently required. In the present study, three kinds of snow depth data which included the D-InSAR data retrieved from the remote sensing images of Sentinel-1 synthetic aperture radar, the automatically measured data using ultrasonic snow depth detectors, and the manually measured data were assimilated based on ensemble Kalman filter. The assimilated snow depth data were spatiotemporally consecutive and integrated. Under the constraint of the measured data, the accuracy of the assimilated snow depth data was higher and met the need of subsequent research. The development of ultrasonic snow depth detector and the application of D-InSAR technology in snow depth inversion had greatly alleviated the insufficiency of snow depth data in types and quantity. At the same time, the assimilation of multi-source snow depth data by ensemble Kalman filter also provides high-precision data to support remote sensing hydrological research, water resource assessment, and snow disaster prevention and control program.


Introduction
Snow is a pivotal storage component in the hydrological cycle (Zhou et al. 2019) and a valuable water resource in Xinjiang Uygur Autonomous region (hereafter "Xinjiang"), China. River runoff replenishment, agriculture and animal husbandry production, and economic development in Xinjiang are dependent on snowmelt. Therefore, the amount of snowfall directly affects the socio-economic development of Xinjiang. Snow depth can respond to snow accumulation, to a certain extent. At the same time, as an indispensable input parameter of medium-and long-term numerical climate predictions, the simulation of land surface hydrological process, water resources assessments and early warning of snow disaster, and the accuracy of snow depth affect the quality of the output results (Armstrong and Brun 2010). Therefore, an accurate grasp of snow depth information is of great practical significance for evaluating water resource reserves and preventing snowmelt floods and avalanche disasters. Xinjiang is a typical arid and semiarid region and thus is very sensitive to climate change. In recent years, with climate change, extreme snowfall events had frequently occurred in Xinjiang, and the amount of snowfall had increased greatly. Moreover, the first snow day had slowly advanced, and the final snow day had obviously delayed. The snow day number, snow depth, and snow index had increased significantly on seasonal and interannual scales (Tang et al. 2013). However, the coarse precision of snow depth data used in the current numerical predictions cannot characterize the change of snowfall that occurs under the background of climate change or provide more flexible auxiliary decision-making information to deal with the threat increased snow disasters. In this context, it is mandatory to focus on how to improve the temporal and spatial resolution Responsible Editor: Amjad Kallel * Chengzhi Li xdlichengzhi@xju.edu.cn 1 of snow depth data in order to allocate precious snow resources more scientifically. Present methods acquiring snow depth data include field measurement and remote sensing inversion. Manual field measurement can even achieve the acquisition of millimeter snow depth data. And the underlying surface and meteorological conditions of the measured position can be obtained during field measurement and are helpful for analyzing the concrete reasons for the variation, distribution, and subregions of snow depth on a regional scale. However, the hourly activities of surveyors (Yue et al. 2016) can destroy the snow surface, break the original energy balance of the observation spot, expedite the melting speed of snow, alter the snow layer structure, and cause the measured snow depth to be lower than the true value. In addition, manual measurement is timeconsuming and labor-intensive. More importantly, the harsh geographical environment prevents investigators from measuring the depth of snow in a large scale. Remote sensing technology can extend the limited point-like representative information obtained by traditional measurement methods to a regional scale, which is more in line with objective needs. With the rapid development of technology, the new methods to accurately monitor many Earth system elements has been widely used in the monitoring and quantitative inversion of surface elements, such as global snow cover, soil moisture, and vegetation (Shi et al. 2012). According to the image type, the methods to obtain snow depth by remote sensing can be divided into optical image inversion, passive microwave inversion, and active microwave inversion. Optical image inversion is based on the principle that the surface reflectivity will increase as the snow depth adds, but will no longer aggrandize when the surface is completely covered by snow. Subsequently, a regression model is established for the inversion. The snow depth information in a large scale can be obtained by using this method, while the inversion errors are large because the relationship between surface reflectivity and snow depth is not very explicit (Liang et al. 2015). To a large extent, optical image is vulnerable to the influence of atmosphere, cloud, and mountain shadow. Passive microwave inversion of snow depth is based mainly on the principle that the microwave brightness temperature difference between the 18 and 37 GHz bands can effectively reflect the snow depth information (Chang et al. 1987;Che et al. 2008;Jiang et al. 2014a, b). Compared with visible and near infrared sensor, microwave sensor has all-weather observation ability and strong penetration ability, which can overcome the influence of atmosphere on imaging and the confusion of cloud on snow. However, this method is not effective to detect light snow (Kelly 2009) and usually underestimates the depth of snow. In addition, the resolution is low, generally between 0.5 and 25 Km, which cannot even reflect the difference of snow depth in adjacent snow units accurately and delicately. Some researchers use active microwave images to retrieve snow depth based on D-InSAR (Differential Interferometric Synthetic Aperture Radar) technology. This method eliminates the phases that are not related to snow from the total phases and then converts them to heights to obtain the snow depth (Li et al. 2014). Based on active microwave imaging, the accuracy of snow depth data has indeed been raised to the meter, although in areas with complex terrain, the accuracy will be affected by image overlay and shadows in some degrees. The evolution and application of snow remote sensing (Rees 2006) have resulted in the fusion of multiscale and multi-source data, which has improved the monitoring capacity of snow depth in terms of the spatial and temporal resolution (Yi et al. 2016) and has overcame the difficulties of artificial snow depth measurement. As the important parameters for snow hydrology and snow disaster models, the snow depth data retrieved by remote sensing can be used to interpret the snow melting process. Both the measured single point datum and the surface snow depth data retrieved from remote sensing images are instantaneous. However, the variation of snow depth needs succession both in time and space. These data have to be integrated to preserve the superiority of the measured and remote sensing-retrieved snow depth data. The ensemble Kalman filter assimilation method can achieve this objective.
In 1960, Hungarian mathematician RE (Rudolf Emil) Kerman introduced the recursive (cyclic) method to solve the problem of discrete data linear filtering (Welch and Bishop 1995), and the Kalman filter assimilation method followed. It was first applied in the domain of atmosphere and ocean research (Ghil and Malanotte-Rizzoli 1991). Epstein innovated the ensemble Kalman filter (Epstein 1969) based on the Kalman filter, and it combines the Kalman filter with ensemble prediction to form a fourdimensional assimilation method for sequential data without a linear model. The ensemble Kalman filtering method has been increasingly applied to the field of geoscience. The potential for data assimilation in snow cover research has been confirmed in many studies (Bergeron et al. 2016). Previous studies assimilated snow cover (Xu and Shu 2016), snow water equivalent, snow temperature (Granlund et al. 2009), snow depth, and other snow parameters. For example, to improve the quality of snow depth data, Zhiguang Tang assimilated the snow depth results of MODIS (moderate-resolution imaging spectroradiometer) passive microwave data with the snow depth data obtained by a snow model. The results showed that the assimilated snow depth was closer to the "true value" (Tang et al. 2016). Tao Che established a snow depth data assimilation system using numerical simulation, passive microwave radiation inversion, and comprehensive observation data, which solved the problem of a single snow depth data source (Che 2006). De Lannoy proposed four different methods to assimilate snow water equivalent data. After assimilating the snow water equivalent data with a 25 km spatial resolution retrieved by the AMSRE (Advanced Microwave Scanning Radiometer for the Earth Observing System), the 1 km spatial resolution SWE (snow water equivalent) data was obtained, which improved the quality of snow remote sensing products (De Lannoy et al. 2012). Considering the crucial character of snow depth data, improving the temporal and spatial accuracy of snow depth data will be tremendously beneficial to hydrological and climatic research. Although the aforementioned exploration can effectively improve the quality of snow depth data, the types of data sources used for assimilation and the resolution after assimilation still need to be supplemented and improved. In this paper, we consider snow depth data obtained from multiple sources, including automatic detection by ultrasonic snow depth detectors, manual measurements, and all-weather, multiangle, high-precision synthetic aperture radar (Sentinel-1) inversion, which is not affected by weather conditions. The ensemble Kalman filter method is used to assimilate the snow depth data from the three data sources to improve the temporal and spatial accuracy of snow depth information.

Study area
The target study area ( Fig. 1) for snow depth assimilation is located in the hinterland of Eurasia, in the Juntang Lake Basin in Hutubi County, Xinjiang, China (43°53′ 62″ N, 86°29′ 00″ E). It is an independent and complete semiclosed small watershed. The water flow in the basin flows through the Piedmont hills, into the plain, and ultimately exits into the Hongshan Reservoir. The overall shape of the study area is a shape of "V," and the area is approximately 841 km 2 . Scarce rainfall and high evaporation have led to a dry climate here. The annual mean temperature is 7.3°C; the multiyear daily minimum and maximum temperature are -27.9°C and 31.7°C, respectively. And the difference between the multiyear hottest and coldest monthly mean temperature is approximately 30°C . The main types of land use are natural pasture, woodland, shrub land, mining land, bare land, scenic spots, reservoir land, dryland, sandy land, highway land, rural roads and facilities, and agricultural land. Mountain meadows, plain meadows, and lowland meadows are the main vegetation types covering the basin. Snowfall in the study area belongs to dry and cold snow formed under the action of continental climate (Hu 2004), and it is mainly distributed in the low mountains, longitudinal valleys, and piedmont folds between 1000 m and 2000 m a.s.l. The first snowfall often occurs at the end of November, and the amount of snow increases rapidly in January, reaches its peak in February, and begins to melt and freeze in mid-March, and the seasonal snow completely disappears at the end of April. The maximum snow storage in 2016 is 79 × 10 4 m 3 , and the maximum snow depth is 60.51 cm. In recent years, snowmelt floods have frequently occurred in gullies. Many of the slopes, valleys, roads, and residential areas located downstream have been eroded and silted by sand, gravel, and soil carried by floods.

Data
The materials used in the snow depth data assimilation research include the data retrieved by Sentinel-1A image, the automatic measurement of the ultrasonic snow depth detector, and the manual field measurement.

Sentinel-1A image acquisition
Sentinel-1 is a dual constellation C-band satellite designed for the Global Environment and Safety Monitoring system project (the "Copernicus" Project), and it was launched in French Guiana in 2014. A single Sentinel-1 satellite maps the world every 12 days, and the repetition period of the two satellite constellations can be shortened to six days. There are four operating modes of the Sentinel-1 radar: IW (interferometric wide), WV (wave), SM (stripmap), and EW (extra-wide) mode. At present, only the IW mode image can be obtained in the study area. The IW mode is the default mode for radar to scan the land area, and the ground resolution is 5 × 20 m, with a width of 250 km. In this model, the surface is scanned progressively between the three subbands, and the radar beam is scanned repeatedly three times in one subband, which leads to strong imaging radiation performance. The study area is completely located in the IW1 band. The IW mode supports dual polarization (HH-HV and VV-VH) operation, and only VV-VH polarized images are available in the study area. Each mode of Sentinel-1 produces SAR-0 level, SLC-1 (single look complex) level, GRD-1 (ground range detected) level, and OCN-2 (ocean) level products. However, only the SLC-1 image contains phase information that is extremely important for inverting snow depth. Therefore, a set of SLC-1 image pairs was selected to obtain snow depth data.

Ultrasonic snow depth data
Brazenec showed the capacity and superiority of ultrasonic technology in detecting snow depth (Brazenec et al. 2004). The ultrasonic snow depth detector was self-manufactured by our team based on the principle of acoustic ranging. It is composed of an ultrasonic ranging sensor (with a temperature compensation module), a microprocessor storage module, and a screen display module. A linear correlation occurs between the air temperature and the sound velocity in the temperature compensation module (Si and Lu 2008). The snow depth is obtained according to the relationship between the sound velocity and the distance.
The ultrasonic probe can be automatically measured in the range of 2-450 cm and is able to capture a minimum angle of 15°. The error is 0.2 cm. The probe emits a set of short sound square waves of 40 kHz every 30 min. The initial emission time is t 0 ; the distance between the ultrasonic probe and the air-snow cover interface is h; the time of ultrasonic return probe is t 1 ; and the duration of the high-level signal is Δt = t 1 -t 0 .
Assuming that the sound velocity is v, then The sound velocity v is a function of the air temperature T; Then, the depth of snow can be expressed as follows: In Formula (3), SD is the depth of snow, and H is the vertical height of the ultrasonic probe. The snow depth is the result of the dynamic change between the accumulation of new snow and the loss of old snow. Ultrasonic snow depth detector can well record the initial and changing snow depths of the entire snow season and snowmelt period and can obtain the centimeter-level series snow depth data set with a halfhour aging without disturbing the energy balance (Ryan et al. 2008). A total of 18 ultrasonic snow depth detectors had been installed in the basin (Fig. 1, blue circle). To reduce the influence of wind blowing, a windshield was set around the detector. After one year's test, the detector had operated normally, and the collection data were continuous and reliable, which showed that it had high cold and moisture resistance. Comparing the data recorded by the detector through conventional deviation correction processing with the data measured manually, the two types of data were more consistent (Yao and Ma 2015). Summing up, the application of ultrasonic snow depth detectors can supplement the types of snow depth data and improve the temporal and spatial resolution of the data.

Measured snow depth data
Only three artificially measured locations were selected in the basin (Table 1), which coincide with the stations 11, 12, and 2 of the ultrasonic snow depth detector. From January 4 to February 1, 2016, the characteristics of the snow depth, density, and moisture content were measured in strict accordance with the International Cryosphere Science Association measurement standard. It is worth noting that these properties can only be measured when more than half of the surface is covered by snow. The measurement is repeated three times at the same location every 2 h, and we take the average value as the final result. The measured snow depth is mainly used to verify the accuracy of the other snow depth data.

Theoretical considerations
The method employed involved:i) using D-InSAR technology to obtain the relative phase of Sentinel-1A SLC interference image pair, unwrapping the phases and eliminate the phases except snow, and then, the snow depth data is obtained by converting the phase of snow to height;ii) the Kalman filter method is used to assimilate the snow depth data automatically recorded by the ultrasonic sounder and the inverted snow depth data to obtain a snow depth value closer to the true value. In this chapter, the related theoretical support will be elaborated, and the visualization process is shown in Fig. 2.
Theory of remote sensing inversion of snow depth D-InSAR is an extension of the frontier InSAR (Interferometric Synthetic Aperture Radar) technology, which can obtain high-precision ground elevation information (Maghsoudi et al. 2013) and has great application potential in the field of geophysical parameter inversion, such as earthquake (Hu et al. 2012) and landslide deformation (Squarzoni et al. 2020). In recent years, it has been extended to the field of snow depth detection. In a nutshell, the principle of D-InSAR retrieving snow depth is that when the same antenna repeatedly transmits pulses to the same surface, it will produce two phase data of the original surface and the snow surface covered by fresh snow . The snow depth can be measured by converting the phase generated by the snowfall event into height. The detailed process is as follows:

i)
Choose the interference image pair. Whether the transit time of the image overlaps with the time period of the manual measurement and the snow depth data recorded by ultrasonic, and whether there is a snowfall event (usually the main image is the image after the snowfall event, the auxiliary image is the opposite), the atmospheric condition of the day and orbital position are the decisive conditions that limit the combination of image pairs (Di 2011). After screening, only the images on January 7, 12, and 23 and February 29, 2016, meet the demand.

ii)
Calculate the space baseline. The baseline is the vector displacement between the positions of the two images (b in Fig. 3), used to evaluate the quality of the interference image pair.
Too short a baseline will reduce the accuracy of snow depth inversion. On the contrary, it will lead to loss of coherence, and the snow depth cannot be described without phase information. The suitable baseline length is 10-1000 m. The baseline can be expressed as follows: In Formula (4), B is the baseline, λ is the wavelength, R is the distance from the range, and R r is the distance from the azimuth.
By calculating the baseline length of the standby image (Table 2), the images of February 29 and January 23, 2016, were identified as the main image and auxiliary image to retrieve the snow depth, respectively.

iii)
Generate interference phase. In Fig. 3, SAR1 and SAR2 are the spatial positions 1 and 2 of the satellite antenna, θ is the incident angle of the satellite beam, P is the snow surface point, and R and R + ΔR are the oblique distances from SAR1 and SAR2 to the ground point P, respectively. Then there are: In the above formula, φ is the phase, where φ scatt1 andφ scatt2 are the random scattering phases of SAR1 and SAR2 respectively. Assuming φ scatt1 = φ scatt2 , the relative phase is When the snow surface changes from P to P', there will be a relative scattering displacement d projected in the oblique moment, that is, iv) Phase unwrapping. The phase change takes 2π as the period, and the phase exceeding 2π will restart and circulate. Phase unwrapping makes the phase correspond to the linearly changing terrain information and solves the problem of 2π ambiguity. The minimum cost flow (Pepe et al. 2011) uses a square grid to mask all pixels in the image whose coherence is less than the threshold. When there is a large area of low coherence or other factors limiting the growth that make unwrapping difficult, this algorithm can achieve better results than the area growth method.

v)
Generating snow phase. The total phase (φ total ) in the study area includes the φ atm caused by atmospheric delay, the φ orb yielded by orbit error, the φ noise due to speckle noise, the φ flat derived from the leveling effect, the geomorphological phase (φ topo ), and the snow phase (φ snow ).
Among them, the expression of φ snow is Under poor weather conditions, the phase collected by the antenna will contain the φ atm . However, in the corresponding time of the image pair, if the weather is clear and cloudless, therefore, the φ atm can be ignored. The route information in the original data is updated by using the POD (precise orbit determination) restituted orbit file provided by ESA (European Space Agency), and the φ orb introduced by the deviation orbit information is removed. Multiple scattering sources are observed in a single pixel in SAR spatial random coherent superposition, and noise will be generated (Herrera et al. 2007). The speckle not only affects the image quality but also makes interferometry difficult. The noise of the image can be reduced by filtering and multiview processing. In this paper, the Goldstein filtering method (Jiang et al. 2014) was used to filter both the real part and the imaginary part of the phase, and the spatial resolution information of the radar image was guaranteed to the maximum extent, while the texture attribute of the radar image was preserved by double filtering. On this basis, a small window (3 × 3) multiview Boxcar filter (Baier et al. 2018) was applied to remove the small noise. By refining the linear information of the image to the baseline, the ground fuzzy height represented by 2π in the phase cycle was corrected, and the φ flat was effectively removed. In summary, the nonsnow phase in the image is removed, and then the snow depth can be obtained by converting the φ snow information into height. Therefore, H snow is the snow depth. γ is the coherence coefficient between the interference image pair.

Principle of the ensemble Kalman filter
Data assimilation is the correction of data deviations. The ensemble Kalman filter originates from dynamic prediction theory (Ahlberg and Gustafsson 2012) and identifies the true value closest to the current state using the measured value and the predicted value of the next state. Inheriting the advantages of the Kalman filter, the amount of calculations is reduced, and it can be applied to nonlinear systems. The principle of the ensemble Kalman filter is based on the Monte Carlo method (Kang 2015) and includes disturbing the background field and observation value, predicting the background error of the model, and averaging the error as the best estimate. The analyses of the observed data and covariance update are synthesized, and the mean value of each updated variable is taken as the best estimate of the variable (Mi et al. 2012). Finally, each initial condition is combined to form a new matrix according to its distribution characteristics (such as a normal distribution) to generate a representative state (Reichle et al. 2002).
A state variable is defined when the snow depth measurement is initialized: where x i (i = 1,n) is a state variable satisfying the Gaussian distribution and T is the size of the set. The snow depth forecast value for each state variable at time t + 1 is as follows: where x f i;t þ 1 is the snow depth prediction value of the state variable i at t + 1; x a i;t is the snow depth analysis value of the state variable i at t time; F is a nonlinear model operator; ε i, t is the Gaussian white noise (background field disturbance) with an expectation of 0 and a variance of Q t ; and Q t is the model error variance matrix.
K is the ensemble Kalman filter gain matrix at t + 1 time: where P is the prediction field error covariance matrix. H is the observation operator, and R is the observation error covariance matrix. The analytical value of each state variable at time t + 1 is as follows: where x a i is the analytical value of the i state variable at t + 1, Y t + 1 is the set of observations at the same time, δ i, k is the disturbance of observation error, R k is the covariance matrix of observation error, and x a tþ1 is the mean value of the set of t + 1 time analysis, that is, the optimal value of the current state. The set covariance P a tþ1 becomes the state error at t + 1 according to P a tþ1 and ,x a tþ1 the snow depth state set at t + 1 time is generated, and then Formula (14) is reiterated to update the state and time.

D-InSAR inversion of snow depth
After performing the baseline calculation, orbit correction, calibration, interference generation, filtering, phase unwrapping, and terrain correction, the snow phase is obtained from the total phase of the Juntang Lake Basin. The period of the phase change is 2π, and the phase of the interferometric image pair in the study area is generally between 0 and π. After unwrapping, a small amount of rupture is observed in the mountain gully, while the effect in other areas is stable. The RMSE (root mean square error) of the X direction pixel is 0.63, and the RMSE of the Y direction is 0.59. The maximum pixel error of X and Y is 1.87 and 1.27, respectively. The surface deformation displacement in the simulated study area is ± 0.014 m. Figure 4 shows that the snow depth range on February 29, 2016, was 19-49 cm and the main snow depth value on that day was 44 cm. The spatial variation of snow depth is consistent with the trend of real snow distribution on the surface. The snow depth in the direction of the northwest monsoon and arctic water vapor source is higher than that in other directions. The snow depth on the south slope is generally smaller than that on the north slope, although with the increase of snow accumulation, the difference becomes small.

Ultrasonic snow depth measurement results
The observation time of the ultrasonic detectors was 47 days (January 20-March 6, 2016). During this period, the snow depth data of 18 sites (Table 3) were obtained every half hour. The Kriging spatial interpolation method was used to estimate the snow depth data (Zhu et al. 2010) for unknown points by using adjacent values based on spatial and temporal correlations . The special dates for snowfall and the first positive temperature but no snowfall date were interpolated. The ultrasonic snow depth data of Sentinel-1A transit time were selected to represent the snow depth of the day. Figure 5 shows that the snow depth range obtained by interpolation on January 23, 2016, is 9-45 cm. Most of the snow depth values are concentrated in the range of 21-29 cm, especially in areas where ultrasonic stations are not installed. The snow depth range interpolated on February 25, 2016, is 7-43 cm, and most of the snow depth values are concentrated in the range of 27-31 cm. Melting of snow occurred on this day. Although snowfall occurred in the study area on February 3 and 10 (less snow), due to the effect of wind blowing, the snow depth in the study area did not increase but rather decreased slightly on February 25.

Results of snow depth assimilation
The ensemble Kalman filter is a sequential assimilation (Evensen 2003). The snow depth data with different data sources, spatiotemporal resolutions, and errors are superimposed on the numerical dynamic model. According to the Monte Carlo method, an optimal solution is found between the numerical model solution and the snow depth data, and it is used to update the initial snow depth field. A continuous loop such as this can update the field and obtain the optimal solution (Roth et al. 2015). Then, the error is corrected, the snow depth prediction value of the numerical dynamic model close to the observed snow depth value is maximized, and the process of snow depth assimilation is completed (Sun 2012). To avoid computing singularities (Zhang et al. 2011), the data were assimilated into columns and columns of the same (2823, 2823) square matrix. Figure 6 shows that the snow depth data obtained by ensemble Kalman filter assimilation are consistent with the measured values at the mean and maximum. Therefore, the optimal estimation of the optimal value is realized. The R 2 value of the assimilated snow depth data is 0.815.

Accuracy verification
To verify the accuracy of snow depth inversion by the D-InSAR method and the effect of the ensemble Kalman filter assimilation of snow depth data, the measured data, the ultrasonic snow depth data, the D-InSAR inversion snow depth data, and the assimilated snow depth data are compared (Fig. 7), and the main results are as follows: (i) Manual measurements are performed only at three stations that coincide with ultrasonic detection stations 11, 10, and 2. The snow depth data (Table 3) obtained by the ultrasonic detector are verified using the manual measurement results (Table 1) as the standard, and all of these data are slightly larger than the manually measured results, although the error does not exceed 3 cm. The minimum error of snow depth data obtained by ultrasonic detection is 3.7%, and the maximum error is 11.8%. Therefore, the accuracy of other snow depth data can be verified using the ultrasonic snow depth. (ii) Although the ultrasonic snow depth data are the input data for the Kriging interpolation, the values of the stations with snow depths less than 40 cm after interpolation are underestimated, and most of the snow depth values greater than 40 cm are overestimated. The overall average accuracy is 82.74%. (iii) Most of the inversion values are lower than the snow depth values obtained by the ultrasonic detector. The accuracy of the data is low, with an average of 78.23%. (iv) Although some of the snow depth data assimilated by the ensemble Kalman filter are underestimated and some are overestimated, these snow depth data are closest to the true value, with an average accuracy of 92.03%.

Discussion
The snow depth data of inversion and Kriging interpolation are assimilated by using ensemble Kalman filter method. Because the underlying surface in winter becomes simpler, and the resolution of the source data image is high, therefore, the image pair used for inversion presents strong coherence, which leads to inversion results that are closer to the measured results on the ground. The inversion results show that the snow depth has the characteristics of a transverse "strip"  distribution, the variation of snow depth in each strip is slight, and the value at the intersection between adjacent strips is between the ranges. The snow depth decreases gradually from west to east. The above phenomena are associated with the relationship between the distribution of snow and the altitude. The snow depth is thick at high altitude and thinner at relatively low altitudes. In addition, the vertical band distribution of vegetation intercepts part of the solar and surface radiation and affects the distribution of snow. In high-altitude and mountainous areas, the coherence of image pairs is strong, and the retrieved snow depth results are closer to the measured values. However, due to the interference of artificial buildings and human activities (such as clearing road snow), areas located closer to urban residents present a greater deviation between the inversion results and the measured values. Topographic conditions also impact the distribution of snow depth. The depth of snow cover on high mountain windward slopes is significantly higher than that in other areas because the warm and humid air flow entering the study area is forced to rise by the terrain, the water vapor condenses, and snow falls after reaching these conditions, which results in a large difference in the snow depth of the windward slope and the leeward break. The spatial resolution of snow depth data based on Sentinel-1 imagery is 5 × 20 m. Compared with existing remote sensing snow depth data, the spatial resolution is improved, although certain errors remain. Because the Sentinel-1 C band is very sensitive to subtle changes in snow properties, changes in the snow layer structure (Fig. 8) directly affect the radar echo. In the early stage of snowfall, the particles on the surface of the snow are small, and the snow only has a density of 0.043 g/cm 3 , although the porosity and permeability of the snow are very large. With the effects of wind, temperature, and other meteorological factors on snow, the particle size of snow gradually develops from the surface to a fine grain snow, medium grain snow, and then coarse grain snow structure. With the continuous action of these factors, the bottom layer gradually appears to become cemented deep frost and polymerized deep frost. New snow will continue to fall on the surface, and over time, the surface will begin to appear as frost, with the snow density reaching 0.120 g/cm 3 . By the end of February, under increasing temperature, the melting and freezing of snow continue to occur in the snow layer. In addition, fine, medium, and coarse snow are no longer observed in the snow layer, which consists mainly of cemented deep frost and aggregated deep frost until the snow completely melts and disappears. These changes in snow structure affect the scattering of ground objects and radar echoes, thus resulting in the overestimation or underestimation of the snow depth. In addition, it is assumed that the snow is in an ideal state (dry snow) when the snow depth is inversed based on the D-InSAR method. However, in practice, SAR is sensitive to the internal state of snow; therefore, the retrieved snow depth is overestimated in the seasonal snow distribution area when snow melts and freezes repeatedly.
The snow depth data obtained by the ultrasonic snow depth detector have high accuracy and high time resolution. An ultrasonic snow depth detector is superior to manual measurement in both the quantity of data and acquisition method. Based on the ultrasonic snow depth data of 18 stations, the snow depth data of unknown points were obtained by Kriging interpolation, and surface snow depth data covering the study area were generated. The results show that the snow depth of each unknown point is controlled by the adjacent known highprecision data and the interpolation results are basically consistent with the geographic distribution trend of the actual snow depth in the study area. Because an ultrasonic snow depth detector station was not placed at high altitude, the Fig. 7 Comparison of snow depth data maximum snow depth is mainly distributed in the middle and south, and the lowest snow depth is mainly concentrated in the northeast of the study area. The interpolated snow depth range is in the ultrasonic data range, although the maximum value is smaller than the ultrasonic snow depth value because the ordinary Kriging method smoothed the ultrasonic snow depth data to some extent, which led to the underestimation of deep snow depths and overestimation of shallow snow depths. Generally, the snow depth range is within the measured value range, the interpolation results are natural, and there are no abnormal values.
The snow depth data obtained by the data assimilation method integrate the advantages of high spatial resolution surface data and high temporal resolution point data; thus, the snow depth value is closer to the measured value. Compared with only retrieved snow depth data, the introduction of ultrasonic snow depth data weakens the absolute influence of coherence on the accuracy of snow depth data, preserves the data accuracy of the regions with strong coherence, and corrects the regions with low coherence. The assimilation data can determine the average value of the interpolated snow depth data in the highaltitude region using the retrieved snow depth data and equalize the distribution of the snow depth value. In general, the assimilated snow depth is closer to the real value on the ground.

Conclusion
In this paper, the D-InSAR method was used to retrieve Sentinel-1A snow depth data, and the Kriging interpolation method was used to obtain surface ultrasonic snow depth data. The ensemble Kalman filter was used to assimilate two types of homologous data. The main results are as follows: i) The distribution of snow depth obtained by inversion is consistent with the trend of the measured values. At high altitudes and with little or no human activity, the snow depth is close to the measured value. In the areas with urban residents and along roads and buildings, the data deviate from the measured values.

ii)
The snow depth data obtained by interpolation are generally consistent with the geographic distribution trend of the actual snow depth in the study area. Because of the "smoothing effect" of interpolation, the range is smaller than that of the ultrasonic data. The snow depth value is within the range, and no abnormal values are observed; therefore, the interpolation results are reliable. iii) The snow depth data obtained by assimilation combine the above two types of snow depth data and are restricted by the accuracy of the ultrasonic data and measured data; however, they are the closest to the measured value. The feasibility of the ensemble Kalman filter Fig. 8 Typical snow structure of the study area assimilation method in improving the spatial and temporal resolution of snow depth data is verified.
This research extends the sources of snow depth data and improves the spatial and temporal resolution of snow depth data. The results can improve the quality of input parameters for hydrological models and meteorological numerical predictions.