High-Resolution Inversion Method for the Snow Water Equivalent Based on the GF-3 Satellite and Optimized EQeau Model

: High-resolution snow water equivalent studies are important for obtaining a clear picture of the potential of water resources in arid areas, and SAR-based sensors can achieve meter-level snow water equivalent inversion. The advanced C-band SAR satellite Gaofen-3 (GF-3) can now achieve meter-level observations of the same area within one day and has great potential for the inversion of the snow water equivalent. The EQeau model is an empirical method for snow water equivalent inversion using C-band SAR satellites, but the model has major accuracy problems. In this paper, the EQeau model is improved by using classiﬁcation of underlying surface types and polarization decomposition, and the inversion of the snow water equivalent was also completed using the new data source GF-3 input model. The results found that: (1) the classiﬁcation of underlying surface types can signiﬁcantly improve the ﬁt between the snow thermal resistance and the backscattering coefﬁcient ratio; (2) the accuracy of the snow density extracted by the GF-3 satellite using the Singh– Cloude Three-Component Hybrid (S3H) decomposition is better than IDW spatial interpolation, and the overall RMSE can reach 0.005 g/cm 3 ; (3) the accuracy of the optimized EQeau model is signiﬁcantly improved, and the overall MRE is reduced from 27.4% to 10.3%. Compared with the original model, the optimized model is superior both in terms of veriﬁcation accuracy and image detail. In the future, with the combination of advanced technologies such as the Internet of Things (IoT), long, gapless, all-weather, and high-resolution snow water equivalent inversion can be achieved, which is conducive to the realization of all-weather monitoring of the regional snow water equivalent. Through the ﬁtting optimization of the snow thermal resistance and the spatial inversion optimization of snow density, the two input parameters of the EQeau model were optimized separately. The thermal resistance ﬁtting values and snow density of different ground snow types were put into Formula (6) to achieve the inversion of the SWE in the study area as a whole. The results are in


Introduction
Snow cover, as the most widely distributed component of the cryosphere and the most significant component corresponding to seasonal changes, can be used as one of the important indicators of climate change [1] and is also an important part of the earth's energy balance system [2]. The snow water equivalent (SWE) is the product of snow depth and density, and it is used to represent the amount of water stored in snow [3]. The successful estimation of the SWE can improve the accuracy of snowmelt runoff prediction and water supply management [4,5], which is of great significance for improving the utilization efficiency of water resources in arid regions and flood control prediction.
At present, many studies on the SWE are still conducted at a large scale, and the inversion mainly adopts field measurement, the long-term observation of ground stations, and the regional observation of satellite remote sensing [6][7][8]. The estimation of the SWE is the main problem to be solved in the hydrological research of basins [7]. Currently, there are three approaches to estimating SWE, namely, using field measurement data to calculate SWE [9][10][11][12], using optical remote sensing products and SWE reconstruction technology to achieve high-resolution SWE inversion [10,13,14], and using passive microwave remote sensing data to calculate SWE [15][16][17]. These three methods all have their own disadvantages and are not suitable for watershed high resolution SWE calculation.
Compared with the coarse resolution of passive microwave sensors, synthetic aperture radar (SAR) retains the advantages of all-day operation. At the same time, it does not have the disadvantage of optical remote sensing of being easily affected by the external environment, such as weather changes. Meanwhile, with its unique sensitivity to the dielectric and geometric features of the object, and its ability to provide potential scattering of medium properties and physical properties, SAR is one of the most promising methods for retrieving the physical properties of snow, such as SWE, snow depth, etc. [18]. The retrieval of the SWE using SAR mainly includes three methods. The first is the inversion algorithm based on the theoretical model [19,20], which was developed based on the backscattering theoretical model of snow cover and the underlying surface and is represented by the multifrequency and dual-polarization data inversion algorithm developed by Shi and Dozier [21]. Based on the theoretical model of the backscattering of snow, the attenuation coefficient and dielectric constant of snow were obtained by separating the contribution of snow depth and snow density to the backscattering coefficient, and then the snow depth and snow density were inverted to obtain the SWE. The derivation process of this kind of inversion algorithm is rigorous and the physical basis is clear. However, this kind of algorithm is too complex and requires the establishment of a complete snow parameter model simulation database, which has low computational efficiency and high data requirements [22]. Multi-frequency and multi-polarization data need to be input to solve the algorithm, and the high cost is not conducive to the popularization and application of the algorithm. The second method uses InSAR to achieve SWE inversion [23,24]. InSAR uses the SAR data of repeated transit to obtain the phase difference of images before and after snow cover. Moreover, there is a definite geometric relationship between the phase difference and the SWE, according to which the direct inversion of the SWE can be realized [25,26]. However, the method requires high-quality data, which limits the application. The third method is the SWE inversion algorithm based on the thermal resistance principle [27], which is called the EQeau model in subsequent developments. This method was first proposed by Bernier M. et al. in 1998 as an inversion algorithm for the shallow snow area in the Eaton River Basin of Quebec, Canada. The principle of this method is to establish a function of the semiempirical expression between the snow thermal resistance and the radar backscattering coefficient ratio to achieve SWE inversion at the watershed scale. This method was later applied to SWE inversion in Lebanon [28], the Manas River Basin in Xinjiang, China [29], and northeast China [30], and achieved good inversion results. The calculation formula of this method is relatively simple and has good universality and accuracy. Moreover, the method does not require SAR data with strict orbit requirements, and only C-band SAR combined with field measurement data can be used to achieve the inversion of the snow water equivalent. Combined with the analysis of model uncertainties by Chokmani et al., there are three main limitations affecting the application of the EQeau model at present: the low accuracy of the snow density generated by the spatial interpolation method [31], the poor correlation between the fitted snow thermal resistance and the backscattering coefficient ratio [31], and the main c-band SAR used, such as Sentinel-1 and RadarSat-2, lacking the temporal resolution required for snow product research [32][33][34]. Gaofen-3 (GF-3) is China's first C-band synthetic aperture radar (C-SAR) imaging satellite, which has up to twelve imaging modes and is one of the SAR satellites with the most imaging modes in the world [35]. The revisit period of the GF-3 satellite is three days. With the launch of the third satellite of GF-3 on 7 April 2022, the total revisit period can be reduced to one day, during which high precision inversion with a one-meter resolution and a one-day revisit time can be achieved. Compared with Sentinel-1 and RadarSat-2 satellites, it is more valuable to study snow parameters.
Based on the EQeau model, this study combines three factors that limit the application of the model and makes improvements to achieve the optimization of model accuracy. The Singh-Cloude Three-Component Hybrid (S3H) decomposition method is introduced to replace the spatial interpolation method to improve the accuracy of the snow density in the input model. The classification of the underlying surface type is added to improve the fitting accuracy of the snow thermal resistance and backscattering coefficient ratio. Combined with the GF-3 SAR imagery, we explore the method's potential when used together with the EQeau model. Our optimized model has the advantages of simple operation and high inversion accuracy, which may be useful for the future inversion of a snow water equivalent at the water-shed scale. Finally, our model provides a reliable scientific basis for hydrological decision-making at the watershed scale.

Overview of the Study Area
The study area ( 32 29.168 E) is located in the middle reaches of the Kelan River Basin in Altay region, Xinjiang, China, as shown in Figure 1. Xinjiang, located in the hinterland of Eurasia and northwest China, is a semi-arid region that is highly sensitive to climate change. The region is characterized by a continental climate and extensive internal flow basins, which are entirely dependent on water from the surrounding mountains [36]. More than 80% of the surface runoff in Xinjiang comes from mountainous areas, and snow and ice meltwater account for more than 45% [37]. The northern Xinjiang region is one of the areas with frequent snowmelt flood disasters in Xinjiang. Against the backdrop of a warmer and wetter climate in northwest China, the hydrological processes of the basin have an obvious response to climate change. The main result is the increase in snow cover, leading to the advance of the flood season [38,39]. Therefore, the study of the physical parameters for snow cover in Xinjiang is not only a key indicator for understanding climate change and improving the assessment and utilization efficiency of regional water resources, but also provides more accurate suggestions for flood prevention and disaster prevention.
The water source of the region mainly depends on the supply of snowmelt runoff in the basin [36]. Kelan River is a tributary of the Irtysh River with a total length of 265 km and is mainly supplied by snow meltwater, rainfall water, and groundwater. Among them, snowmelt runoff contributes the most to the river, accounting for about 45% on average. The Kelan River basin has good snow research conditions. The snowfall period in the high mountain area is in early September, while in the plains area, it is in early and middle September, and it melts in the April of the following year [37]. At the same time, the Altay Meteorological Station located along the foothills has abundant snow measurement data, which can provide support for high-resolution snow research. The terrain of the study area gradually flattens from north to south. The northern region is mountainous, above 800 m. The central region is relatively flat, with an altitude of 700-800 m. The lowest altitude is about 600 m in the southwest, and the terrain in the northeast is relatively fragmented. According to the statistics on land use data (Figure 1), the types of land use in the study area in the order from most to least used are grasslands (48.86%), barren lands (41.79%), cropland (7.00%), and bodies of water (2.35%). The distribution among the categories is relatively regular. The water source of the region mainly depends on the supply of snowmelt runoff in the basin [36]. Kelan River is a tributary of the Irtysh River with a total length of 265 km and is mainly supplied by snow meltwater, rainfall water, and groundwater. Among them, snowmelt runoff contributes the most to the river, accounting for about 45% on average. The Kelan River basin has good snow research conditions. The snowfall period in the high mountain area is in early September, while in the plains area, it is in early and middle September, and it melts in the April of the following year [37]. At the same time, the Altay Meteorological Station located along the foothills has abundant snow measurement data, which can provide support for high-resolution snow research. The terrain of the study area gradually flattens from north to south. The northern region is mountainous, above 800 m. The central region is relatively flat, with an altitude of 700-800 m. The lowest altitude is about 600 m in the southwest, and the terrain in the northeast is relatively fragmented. According to the statistics on land use data (Figure 1), the types of land use in the study area in the order from most to least used are grasslands (48.86%), barren lands (41.79%), cropland (7.00%), and bodies of water (2.35%). The distribution among the categories is relatively regular.

GF-3 Image Data
According to the requirements of C-band SAR images by the EQeau model in autumn and winter, the SAR images of GF-3 in autumn and winter were obtained as the input data of the model. The GF-3 satellite is equipped with a C-band multi-polarization synthetic aperture radar, which was launched on 10 August 2016. It is capable of all-weather operations and can penetrate clouds, vegetation, sand layers, ice, and snow [35]. With a maximum spatial resolution of one meter, GF-3 is currently the SAR satellite with the most imaging modes, including twelve imaging modes such as strip imaging, scanning imaging, and cluster imaging [40]. The third GF-3 satellite was launched on 7 April 2022. Currently, the three GF-3 satellites are evenly distributed on the same orbital plane, and each satellite takes 99 min to orbit the Earth. The maximum revisit period for global observations has been shortened from 84 hours per satellite to twelve hours per satellite, and the average revisit period has been shortened from twelve hours per satellite to 4.8 h per satellite. The same region can be revisited at least two times per day and five times per day on average, and high dynamic imaging can be carried out in the global region. Therefore, GF-3 has the advantages of all-weather detection, high spatial resolution, high temporal resolution, etc. At present, the satellite is rarely applied to snow research. Users can log onto the GF-3 image resources satellite data distribution platform at the Center for Inquiry and Order (http://36.112.130.153:7777/DSSPlatform/index.html (accessed on 28 February 2022)). In this paper, two 1-A-level GF-3 images are used to achieve high-resolution SWE inversion based on the GF-3 through model optimization. The image parameters used in this paper are shown in Table 1. The spatial resolution of the two images is 8 m. This paper uses the overlapping area of the two images as the research area to carry out research and analysis ( Figure 1).

Field Measurement Data
Data input, parameter optimization, and accuracy verification of the EQeau model requires access to information on field parameters such as snow depth, snow density, snow water equivalent, and underlying conditions for the same period as the images. The synchronous measurement experiment on snow cover was carried out in the study area from 25 February to 1 March 2022. During the measurement period, there was no rainfall or snowfall in the study area, the weather was stable with no significant temperature or strong wind, and the maximum temperature was kept below 0 • C. Therefore, it was considered that the snow surface condition had not changed significantly. The distribution of measurement points is shown in Figure 1.
The main measurement parameters are snow density, snow depth, and the SWE. The main instruments used are a snow bucket, ruler, electronic scale, and electronic thermometer. Snow sampling is accomplished by inserting the bucket vertically into the bottom of the snow and lifting it along the bottom. At this point, the total weight and depth of the snow can be recorded. The final data for each measurement site were obtained by finding the average of three suitable points within 20 m. The unit of the SWE is g/cm 2 or kg/m 2 . Thus far, the snow depth, snow density, and SWE values of all measurement points are obtained as the optimization input parameters a model's optimization input parameters and verification parameters the measurement points is shown in Table 2.

Other Data
In addition to GF-3 remote sensing image and field measured data, the digital elevation model (DEM) and land use coverage are also needed in this paper to achieve the further detailed division of the study area. The DEM used was ASTER GDEM 30 m data (obtained at https://www.gscloud.cn/ (accessed on 10 January 2022)). Therefore, Yang et al.'s annual land cover data set in China based on Landsat-8 were selected [41]. The resolution of this dataset is 30 m, and the data are obtained from remote sensing images with a long time series and high precision. The topography and land use in the study area are shown in Figure 1.

Research Design
In this paper, the EQeau SWE inversion model is improved based on an advanced GF-3 SAR satellite's fully polarized images, combined with field observations of the study area in the same period, the S3H decomposition method, and the land classification of the study area. The research route is shown in Figure 2, which refers to the work of Singh et al. [18] and mainly includes the steps described below.

EQeau Model
The EQeau model was first proposed in 1998 and was applied by Bernier et al. to invert the SWE model of the Eton River Basin in Canada [27]. The EQeau model uses snow thermal resistance instead of snow depth as the input parameter of the snow water equivalent via an empirical formula. Meanwhile, the model is based on the correlation between the snow thermal resistance and backscatter coefficient ratio and is used jointly with Cband SAR images to obtain snow thermal resistance values in the study area. Moreover, the snow water equivalent is calculated by combining the spatial interpolation results of the observed snow density values. The principle of this model is as follows: Firstly, the model expression is established according to the empirical formula. In general, the snow water equivalent is calculated by the product of the snow depth and the snow density.
where represents the snow water equivalent (kg/m 2 ), represents snow depth (m), and represents snow density (kg/m 3 ).
At the same time, snow thermal resistance is the ratio of snow depth to snow thermal conductivity, which is an index characterizing snow's ability to hinder heat conduction.

=
(2) where represents the thermal resistance of snow cover (m 2 ×k/w) and represents snow thermal conductivity. According to Raudviki 's research results [42], snow thermal conductivity is a function of snow density, as shown in Formula (3): Step 1: The GF-3 image was preprocessed to calculate the ratio of backscattering coefficients.
Step 2: Combined with radar parameters and surface parameters, the semi-empirical formula of the backscattering coefficient ratio and snow thermal resistance was obtained according to the results of the ground classification, and the model influence of the constrained backscattering coefficient ratio was realized.
Step 3: The S3H decomposition was used to invert the snow density in the study area to reduce the impact of the snow density error on the model's accuracy.
Step 4: The SWE of the study area was calculated according to the optimized semiempirical formula, and the EQeau model was improved.

EQeau Model
The EQeau model was first proposed in 1998 and was applied by Bernier et al. to invert the SWE model of the Eton River Basin in Canada [27]. The EQeau model uses snow thermal resistance instead of snow depth as the input parameter of the snow water equivalent via an empirical formula. Meanwhile, the model is based on the correlation between the snow thermal resistance and backscatter coefficient ratio and is used jointly with C-band SAR images to obtain snow thermal resistance values in the study area. Moreover, the snow water equivalent is calculated by combining the spatial interpolation results of the observed snow density values. The principle of this model is as follows: Firstly, the model expression is established according to the empirical formula. In general, the snow water equivalent is calculated by the product of the snow depth and the snow density.
where SWE represents the snow water equivalent (kg/m 2 ), D S represents snow depth (m), and ρ S represents snow density (kg/m 3 ). At the same time, snow thermal resistance is the ratio of snow depth to snow thermal conductivity, which is an index characterizing snow's ability to hinder heat conduction.
where R S represents the thermal resistance of snow cover (m 2 × k/w) and K S represents snow thermal conductivity. According to Raudviki 's research results [42], snow thermal conductivity is a function of snow density, as shown in Formula (3): where A = 2.83056 × 10 −6 , B = −9.09947 × 10 −5 , and C = 3.19739 × 10 −2 . Thus, the calculation formula of the SWE can be transformed into Formula (4).
The second step is to obtain the snow density in the study area. According to the EQeau model, researchers usually use the spatial interpolation method of the observation points to obtain the snow density in the study area. The spatial interpolation method is a spatial statistical method employed to obtain the overall values of the study area based on the values and spatial relationships of the observation points.
The third step is to explore the relationship between the backward scattering coefficient ratio and the thermal resistance of snow accumulation in the study area using parameter fitting. According to the conclusions of Bernier et al., parameter fitting can be carried out between the backscattering coefficient ratio of the C-band SAR image and the snow thermal resistance. It is usually expressed in logarithmic form, as shown in Formula (5).
where BR is the backscattering coefficient ratio. Therefore, by establishing an empirical formula, the snow thermal resistance is replaced by the backscattering coefficient ratio, as shown in Formula (6), to achieve the purpose of using SAR to extract the snow water equivalent.
Chokmani et al. analyzed the EQeau model in 2006 and found that the main factors affecting the model lie in two areas [31]: the results of the spatial interpolation of snow density and the results of the fitting of the backward scattering coefficient ratio to the snow thermal resistance. To improve the accuracy of the model, it is necessary to also improve the accuracy of these two components.

Classification of Underlying Surface Types
According to the existing research conclusions [27,43], the influence of the backscattering coefficient ratio on the model mainly comes from the soil roughness, soil moisture, radar polarization mode, and local incidence angle, and these factors make the fit of the snow thermal resistance to the backscattering coefficient ratio decrease and affects the model's accuracy. The backscattering coefficient ratio is inversely related to the soil roughness and soil moisture, and positively related to the local incidence angle. Additionally, the backscattering coefficient ratio is influenced by the radar polarization mode. When the local incidence angle is not fixed, the polarization mode has more influence. However, the copolarization mode can suppress the effect of the local incidence angle on the backscattering coefficient ratio in winter and autumn and is the model's optimal polarization mode.
Since the land use types in the study area are not the same, the topography is undulating and the specific values of soil moisture and soil roughness in the study area are missing, it is not possible to achieve quantitative constraints on the influence of these factors on the backscattering coefficient ratio. Therefore, it is necessary to optimize the fit between the snow thermal resistance and the backscattering coefficient ratio in winter and autumn by optimizing the model's parameters to suppress the influence of the polarization mode and other effects for the specific topography and land use type conditions in the study area. According to existing studies, it has been shown that there is a relationship between soil roughness, soil moisture, and underlying surface type. The distribution of surface roughness and soil moisture is more concentrated in the same underlying surface type [44][45][46][47]. Therefore, adding the type of underlying surface to the model as a constraint condition can reduce the influence of external factors on the backscattering coefficient to improve the fitting accuracy. As shown in Figure 1B, the main underlying surface types in the study area include cropland, grassland and barren. Meanwhile, there is a large difference in the amount of solar radiation received on the sunny and shady slopes, which leads to a large difference in soil moisture. Therefore, the digital elevation model was used to extract the shady slope and sunny slopes for further classification. According to the conclusion of Wang et al., using the satellite incidence angle as the dividing boundary can effectively reduce the influence of local incidence angle on the fitting [29]. Finally, according to the type of underlying surface in the study area, the incident angle, shady slopes, and sunny slopes are introduced into each category for further division. Because the cropland is located on sunny slopes, and the local incidence angle is smaller the than incidence angle, it is not further divided. The research area was finally divided into six categories; the specific categories are shown in Table 3.

Singh-Cloude Three-Component Hybrid Decomposition (S3H Decomposition)
In previous studies of SWE inversion based on the EQeau model, the input of snow density was usually the average value of the measured data or spatial interpolation, which could not effectively describe the variation in the snow density in the study area. Other methods using SAR images to extract snow density, such as the co-polarized phase difference (CPD) method [48], have higher requirements on SAR images, which require two overlapping images as a reference. Thus, the method is less applicable. Singh-Cloude threecomponent hybrid (S3H) decomposition is a fully polarized SAR image decomposition method proposed by Singh et al. in 2017. S3H decomposition is a modification of Cloude's three-component mixed decomposition scheme [49]. It uses a rotating 3 × 3 coherence matrix [T] to rotate the orientation angle of the solution around the line of sight [50] and implements an extended volume scattering model. The principle of this method is mainly to use the volume scattering parameter obtained after the polarization decomposition and the Fresnel transmission coefficient to establish an equation to estimate for the dielectric constant of the snowpack. The snow density is then found by the empirical relationship between the dielectric constant of the snowpack and snow density. The volume scattering parameters obtained through decomposition have been used in snow density inversion research in many places, and good results have been achieved [18,51].
In the existing S3H extended decomposition method, the spherical shape of snow particles is taken into account to accurately express the volume scattering model [18]. The results of S3H decomposition mainly include surface scattering, double projectile scattering, and the volume scattering power coefficient, while on the surface adjacent to snow cover, double projectile scattering has little influence and can be ignored [52]. S3H is decomposed as follows: where m s , m d and m v are the surface scattering, double-bounce scattering, and volume scattering power coefficient, respectively. The data we selected are GF-3 polarization data, which can be transformed into polarization coherence matrix T(θ). |γ| 2 is the generalized volume scattering parameter. This method applies the orthogonality principle of surface and dihedral angular components. The orthogonality condition is expressed as α d + α s = π 2 , ∅ d − ∅ s = ±π [52]. ∅ d and ∅ s are the scattering phases of surface and double bounce scattering, respectively. Through orthogonality conditions, the two are combined into ∅ in Formula (7).
Moreover, since a large number of randomly distributed snow particles can be used as Rayleigh scattering objects, the average volume scattering coherence matrix can be expressed as the integral of the volume scattering coherence matrix of randomly distributed spherical particles in a resolution unit at all possible angles θ: Among them, γ HH and γ VV are Fresnel transmission coefficients of the HH and VV polarization of SAR images respectively. According to Fresnel's transmission equation and Snell's law, Fresnel's transmission coefficient is expressed as follows: Therefore, when the double-bounce scattering is assumed to be negligible, that is, when m d in Equation (7) tends toward zero, |γ| 2 can be written as follows: where T 11 (θ) denotes the elements of the first row and first column of the T(θ) matrix, and so on. Therefore, |γ| 2 can become the intermediate variable connecting Formula (7), Formula (10), Formula (11), and Formula (12), thus establishing the connection between s and the decomposition matrix, where s is the dielectric constant of snow. The dielectric constant of snow can be obtained by solving the above equation. Combined with the above formula, the following formula can be obtained: The derivation of Equation (14) is that the dielectric constant of the snowpack is an unknown parameter. The dielectric constant of the snowpack derived from Equation (14) is related to the snow density, which is calculated using the existing empirical formula [53]. The relationship between snow density and the dielectric constant of the snowpack is: Thus far, the polarization matrix information contained in the fully polarized SAR image can be obtained by the polarization decomposition method, and the snow density distributed in the study area can be obtained by establishing the relationship between the dielectric constant of the snowpack and the polarization scattering matrix. The results in the formula are algebraic calculations of the matrix after the polarization decomposition, we implement the algebraic calculation of the matrix through Matlab software. The snow density obtained by inversion is compared with the observed values and conventional spatial interpolation methods in Section 4.2.

Optimization of Snow Thermal Resistance Fitting
The EQeau model is a semi-empirical function. As analyzed in Section 3.2, the specific model fitting coefficient will be affected by external conditions. According to the auxiliary data and measured data, the study area is divided into land classes as shown in Table 3, excluding the water areas that are not suitable for the model in the study area. Observation points are divided into different categories according to the land use type of the location for the simulation and verification of different categories. The fitting points in each category accounted for about 80% of the total observation points. To minimize the influence of manual intervention on the experiment, we use random functions to determine the fitting points and verification points in each category. The numbers of fitting points and verification points are shown in Table 3. According to the experience of the EQeau model and the distribution law of measured data in this paper, the logarithmic function was selected as the fitting function to realize the function fitting between the ratio of measured snow thermal resistance and the backscattering coefficient after land classification. The R 2 value was used to measure the fitting degree between the ratio of the snow thermal resistance and the backscattering coefficient, and the root mean square error was added to measure the deviation between the observed value and the fitted value.
After classifying the study area (Figure 3), the type of snow around the measured values of the thermal resistance and the backward scattering coefficient ratio increased significantly; the unclassified area's fitting degree, which was lower than 0.62, increased to a fitting degree of over 0.75 in the classified area, including the optimal fitting degree of the bare land types, reaching a function fitting accuracy of more than 90%. It can be seen from the table that the root mean square error between the backscattering coefficient ratio fitted by the snow thermal resistance and the GF-3 inversion value is in the range of 0.195 to 0.722. It means that the fitting degree between the snow thermal resistance and the backscattering coefficient ratio is better after the classification of the study area, which reduces the error caused by the fitting between the two and achieves the purpose of reducing the influence of the backscattering coefficient ratio on the model's accuracy.  By integrating the fitting relation between the ratio of the snow thermal resistance and the backscattering coefficients in winter and autumn obtained by fitting the above regions, the inverse function was calculated and put into Formula (6), the model coefficients were revised, and finally, a complete EQeau model with optimized parameters suitable for the terrain and underlying surface conditions in the research area was formed. The model coefficients are shown in Table 4.  By integrating the fitting relation between the ratio of the snow thermal resistance and the backscattering coefficients in winter and autumn obtained by fitting the above regions, the inverse function was calculated and put into Formula (6), the model coefficients were revised, and finally, a complete EQeau model with optimized parameters suitable for the terrain and underlying surface conditions in the research area was formed. The model coefficients are shown in Table 4.

S3H Decomposition's Snow Density Extraction Results
S3H decomposition was used to establish the relationship between the snow dielectric constant and volume scattering parameters, and the snow dielectric constant of the snowcovered area in the study area was obtained. The density of snow cover in the study area was calculated according to Formula (15). At the same time, the conventional method was used to carry out spatial interpolation by fitting the measured values of snow density at the observation points. The selected spatial interpolation method is inverse distance weight interpolation, and the test points are reserved for the snow density interpolation test. As can be seen from the comparison above (Figure 4), the snow density obtained by the two methods is in the range of 0.13 g/cm 3 and 0.25 g/cm 3 . The densities were classified into four classes in the range of 0.03 g/cm 3 . Superimposed with a topographic map ( Figure 5B), it was found that the snow density in the study area decreased gradually along the elevation and from south to north along the latitude, which accords with the law of the snow density measurement in the study area. At the same time, the snow density obtained by S3H decomposition is more refined, and the variation in the snow density at the spatial scale is more obvious, which can better reflect the change in the snow density in the study area.

S3H Decomposition's Snow Density Extraction Results
S3H decomposition was used to establish the relationship between the snow diel tric constant and volume scattering parameters, and the snow dielectric constant of t snow-covered area in the study area was obtained. The density of snow cover in the stu area was calculated according to Formula (15). At the same time, the conventional meth was used to carry out spatial interpolation by fitting the measured values of snow dens at the observation points. The selected spatial interpolation method is inverse distan weight interpolation, and the test points are reserved for the snow density interpolati test.
As can be seen from the comparison above (Figure 4), the snow density obtained the two methods is in the range of 0.13 g/cm 3 and 0.25 g/cm 3 . The densities were classifi into four classes in the range of 0.03 g/cm 3 . Superimposed with a topographic map (Figu 5B), it was found that the snow density in the study area decreased gradually along t elevation and from south to north along the latitude, which accords with the law of t snow density measurement in the study area. At the same time, the snow density obtain by S3H decomposition is more refined, and the variation in the snow density at the spat scale is more obvious, which can better reflect the change in the snow density in the stu area.  The snow density value of the verification sample was compared with the results of two types of snow density generation ( Figure 6). The RMSE of snow density and the validation sample obtained by S3H decomposition and the IDW spatial interpolation method were 0.005 g/cm 3 and 0.014 g/cm 3 , respectively. The accuracy of the snow density inversion using S3H decomposition was better. The snow density value of the verification sample was compared with the results two types of snow density generation ( Figure 6). The RMSE of snow density and the v idation sample obtained by S3H decomposition and the IDW spatial interpolation meth were 0.005 g/cm 3 and 0.014 g/cm 3 , respectively. The accuracy of the snow density inve sion using S3H decomposition was better. After fitting, it can be found that the snow density obtained by S3H decompositi is closer to the actual value, which may be because compared with the spatial interpo tion method, which relies on the measured snow density distributed in different areas cannot reflect the spatial variation in the snow density. S3H decomposition can better r The snow density value of the verification sample was compared with the result two types of snow density generation ( Figure 6). The RMSE of snow density and the v idation sample obtained by S3H decomposition and the IDW spatial interpolation meth were 0.005 g/cm 3 and 0.014 g/cm 3 , respectively. The accuracy of the snow density inv sion using S3H decomposition was better. After fitting, it can be found that the snow density obtained by S3H decomposit is closer to the actual value, which may be because compared with the spatial interpo tion method, which relies on the measured snow density distributed in different area cannot reflect the spatial variation in the snow density. S3H decomposition can better flect the spatial distribution of the snow density at a high-resolution scale through scattering characteristics of the snow body. Therefore, it can be considered that the spa snow density distribution obtained by S3H decomposition is more accurate, and the sn density map obtained by S3H decomposition inversion can replace the traditional spa interpolation method of snow density generation to reduce the model error more ef tively. After fitting, it can be found that the snow density obtained by S3H decomposition is closer to the actual value, which may be because compared with the spatial interpolation method, which relies on the measured snow density distributed in different areas, it cannot reflect the spatial variation in the snow density. S3H decomposition can better reflect the spatial distribution of the snow density at a high-resolution scale through the scattering characteristics of the snow body. Therefore, it can be considered that the spatial snow density distribution obtained by S3H decomposition is more accurate, and the snow density map obtained by S3H decomposition inversion can replace the traditional spatial interpolation method of snow density generation to reduce the model error more effectively.

Results of the Optimized EQeau Model's Extraction of Snow Water Equivalent
Through the fitting optimization of the snow thermal resistance and the spatial inversion optimization of snow density, the two input parameters of the EQeau model were optimized separately. The thermal resistance fitting values and snow density of different ground snow types were put into Formula (6) to achieve the inversion of the SWE in the study area as a whole. The results are shown in Figure 7.

Results of the Optimized EQeau Model's Extraction of Snow Water Equivalent
Through the fitting optimization of the snow thermal resistance and the spatial version optimization of snow density, the two input parameters of the EQeau model w optimized separately. The thermal resistance fitting values and snow density of differe ground snow types were put into Formula (6) to achieve the inversion of the SWE in t study area as a whole. The results are shown in Figure 7. After the optimization of the model, the overall snow weight in the study area distributed between 23.7 mm and 87.6 mm, with an average of 48.7 mm. The SWE of t optimized model clearly shows an increasing trend along with the increase in latitu which is consistent with the trend of the terrain increasing from south to north in the stu area ( Figure 5A). Among different land types, the SWE of bare land and cultivated la on the sunny slope is the lowest, and the main SWE distribution is 20-40 mm. The SW obtained from the inversion is relatively high, mainly distributed at 60-90 mm in t grassland with a higher altitude and northern latitude. The performance of the SWE also different under different constraints of the same underlying surface, which is main reflected in the fact that the SWE in shady slope areas is generally higher than in sun slope areas. This may be because sunshine, wind, and other conditions affect the accum lation of snow on the sunny slope, making the SWE value on the shady slope higher.
At the same time, compared with the results of the model without optimization, t estimation of the SWE between different land types is clearer (Figure 7). The optimiz model can clearly distinguish that the snowfall value of the cultivated land in the northe part of the study area is low, but not mixed with the northern grassland. It can be se that the SWE distribution without optimization is mainly concentrated at 20-60 mm. the same time, the difference in the SWE between different land types decreases and difficult to distinguish. This is mainly because the model is not optimized and the fitti degree between the ratio of the snow thermal resistance and the backscattering coefficie is insufficient, which leads to great errors in the simulation of the overall snow therm resistance value. On the whole, the optimized model can better realize SWE inversion the study area. After the optimization of the model, the overall snow weight in the study area is distributed between 23.7 mm and 87.6 mm, with an average of 48.7 mm. The SWE of the optimized model clearly shows an increasing trend along with the increase in latitude, which is consistent with the trend of the terrain increasing from south to north in the study area ( Figure 5A). Among different land types, the SWE of bare land and cultivated land on the sunny slope is the lowest, and the main SWE distribution is 20-40 mm. The SWE obtained from the inversion is relatively high, mainly distributed at 60-90 mm in the grassland with a higher altitude and northern latitude. The performance of the SWE is also different under different constraints of the same underlying surface, which is mainly reflected in the fact that the SWE in shady slope areas is generally higher than in sunny slope areas. This may be because sunshine, wind, and other conditions affect the accumulation of snow on the sunny slope, making the SWE value on the shady slope higher.
At the same time, compared with the results of the model without optimization, the estimation of the SWE between different land types is clearer (Figure 7). The optimized model can clearly distinguish that the snowfall value of the cultivated land in the northern part of the study area is low, but not mixed with the northern grassland. It can be seen that the SWE distribution without optimization is mainly concentrated at 20-60 mm. At the same time, the difference in the SWE between different land types decreases and is difficult to distinguish. This is mainly because the model is not optimized and the fitting degree between the ratio of the snow thermal resistance and the backscattering coefficient is insufficient, which leads to great errors in the simulation of the overall snow thermal resistance value. On the whole, the optimized model can better realize SWE inversion in the study area.

Evaluation of Model Accuracy
In the results section, we show the spatial distribution characteristics of snow water equivalents retrieved by the model before and after optimization. In order to verify whether the results obtained by our model are good, it is necessary to compare the model accuracy before and after optimization as well as the accuracy of the optimized model and other methods. We calculated the root mean square error (RMSE), mean absolute error (MAE), and mean relative error (MRE) to compare the original model with the model under different constraints to reflect the difference in accuracy before and after the model's optimization. The values involved in the accuracy evaluation are randomly selected validation sample points. The specific values are shown in Tables 5 and 6. As can be seen from the values in the table, the simulation accuracy of the optimized model is significantly improved in all the local classes. The RMSE distribution of the optimized local classes ranges from 3.76 mm to 7.09 mm. Compared with the model before optimization, the average RMSE decreases by about 6.29 mm, and the estimation accuracy of the SWE is significantly improved. However, the unoptimized EQeau model is not ideal in the study area, with an overall MRE of 27.4%. Although it can still roughly reflect the distribution of the SWE in the study area, the fitting degree between the overall snow thermal resistance and the ratio of the backscattering coefficient decreases due to the lack of constraints, resulting in a wider range of model errors. At the same time, compared with previous studies without snow density process optimization, the overall MRE is 19.5% [54]. In this paper, S3H decomposition was added to optimize snow density, which further improved the model's accuracy.
At the same time, we present the results of several authors using different methods to invert the snow water equivalent for accuracy comparison ( Table 7). The comparison with other inversions of snow water equivalents using SAR data shows that our results are better in terms of accuracy and have the potential to be applied to basin-scale inversions of snow water equivalents. Table 7. Comparison of snow water equivalent accuracy among different researchers.

Researcher-Time The Average RMSE (mm)
Patil-2019 [48] 160 Conde-2019 [55] 5.3 Singh-2020 [56] 42.95 Santi-2021 [57] 34.8 Based on the analysis of the results of different researchers above, it can be seen that the inversion of snow water equivalent by remote sensing means will more or less have errors. The error of the model is mainly caused by the error of the measurements and the difference in the mathematical equation of the model in different study areas.
The error of the measurement can be reduced by improving the measurement means, such as using more accurate measuring instruments to measure the snow parameters. The error of the model itself mainly comes from the differences caused by applying the model to different research areas. There are different environment, climate and other conditions in different study areas, which will affect the physical relationship of related snow parameters. Therefore, the application of the improved EQeau model in other regions still needs to consider the environmental conditions of the study area and make corresponding adjustments. In Section 5.3, we will analyze the current limitations of the model and the outlook for the future development of the model.

Effect of the Number of Fitted Points on the Accuracy of the Model
In this study, measurement data are distributed at all locations in the study area, and the number of observation points also meets the simulation requirements of the model. The fitting between the snow thermal resistance and backscattering coefficient is performed at a ratio of eight to two for the fitting and validation points; however, the model accuracy between the fitting and validation points at different ratios also needs to be discussed in this paper. To avoid overfitting and underfitting, the ratios of 6:4 and 7:3 are chosen for comparison with the results in Table 5, respectively. The values involved in the accuracy evaluation are the validation sample points selected using random selection. The specific values are shown in Table 8. As can be seen from the values in Table 8, the model accuracy will vary with the number of fitting points and verification points. The more fitting points there are, the higher the simulation accuracy of the model. To avoid the insufficient number of verification points caused by increasing the number of fitting points, this paper selects the fitting points and verification points to conduct the test at a ratio of eight to two.

Uncertainty and Outlook
In this paper, based on the GF-3 SAR data, the EQeau model was optimized to re-alize SWE inversion at the watershed scale, and the results met the needs for practical application. However, it is still difficult to invert snow parameters at the watershed scale due to the difficulty in obtaining field observation data and the limitations of semi-empirical models.
The EQeau model itself has certain limitations. The application of the model should satisfy its basic conditions, the model needs to be applied to shallow snow areas with soil temperatures below 0 • C in winter, and the data used should also be limited to C-band SAR images [27]. EQeau has been used many times around the world. In 2001, Turcotte et al. discussed the method of combining the model with a distributed hydrological model to improve the simulation accuracy of spring snowmelt runoff [58]. In 2006, Chokmani et al. analyzed the uncertainty of the EQeau model in terms of the model coefficients, backscatter coefficient ratio and snow density, and analyzed their influence on the application of the model. At present, this method has been applied in many places in Canada and achieved good results [59]. Corbane et al. used the model in 2005 to retrieve snow water equivalents over the Lebanese plateau [28]. The conclusion proved that the basic conditions of the model should be explored before using it. At the same time, the terrain conditions and field snow observation conditions have a certain impact on the accuracy of the model application. Wang et al. and Zhu et al. applied the model in Xinjiang and northeast China, respectively, and achieved good results [30,54]. Combined with the application of the model by the above scholars, it can be posited that the EQeau model can be used in regions that meet the basic conditions of the model. The application of high-resolution remote sensing images provides conditions for generating high-resolution digital elevation models and land use coverage data. In addition, the GF-3 image includes a full polarization mode, allowing S3H decomposition to also be performed. We believe that the optimized model can be used in locations that meet the preconditions of the EQeau model. However, it needs to be divided according to the topographic features and underlying surface types of specific locations.
With the use of GF-3 satellites, the single-day review of C-band images has become a reality, allowing regional SWE monitoring to occur on a finer time scale. Our study proves that the Gaofen-3 satellite combined with the optimized EQeau model has a good application accuracy. At the same time, because snow-covered areas are mostly mountainous areas in middle and high latitudes, complicated topography and other conditions greatly limit the acquisition of field observation data, as well as the number and spatial distribution of observation points at the watershed scale. At the same time, due to the efficiency of manual work, it is difficult to achieve the synchronous observation of surface data with satellite images. All these factors will influence the final result of the model. Therefore, it is necessary to change the acquisition method of field-measured data to achieve timely, rapid, and accurate parameter observation. At present, with the rise of the Internet of Things and other technologies, some researchers have applied Internet of Things research to regional environmental monitoring. For example, Fang et al. designed an environmental Internet of Things system for environmental parameter monitoring with climate change and ecological monitoring as the theme [60]. Combined with technologies such as the Internet of Things (IoT), this means that more accurate measurements of snow parameters can be obtained through smarter devices. Meanwhile, the application of IoT can realize the measurement of snow parameters at exactly the same time as the satellite imaging so as to reduce the model error caused by the measurement value.
Therefore, we believe that further work is to realize the observation of snow parameters and surface parameters at a finer time scale by combining GF-3 images and the Internet of Things. In addition, the relationship between parameters affecting SWE is further summarized based on observation data for different regions of the model to achieve SWE inversion that is applicable to a larger range of river basins.

Conclusions
This paper is based on the advanced C-band SAR imagery GF-3, and referring to the existing research, the EQeau model is optimized to retrieve the high-resolution SWE. The application of GF-3 to basin-scale SWE studies has not been carried out yet, and this paper explores the feasibility of its application to SWE studies by improving the EQeau model. The input parameters of the model are mainly the radar backscatter coefficient ratio and snow density, and the two parts are improved using different methods. Firstly, we combine the actual measurement data of each land class with the land cover of the study area to classify the land classes and achieve the optimization of the fitting of the snow thermal resistance and backscattering coefficient ratio under different land classes. At the same time, the relationship between the snow dielectric constant and the decomposed volume scattering constant is established using the S3H polarization decomposition method to realize the inversion of high-resolution snow density in the study area.
The fitting degree between the snow thermal resistance and the backscattering coefficient ratio can be improved by using land classification. Most of the optimized R 2 ranges from 0.8 to 0.95, proving that the fitting accuracy has been significantly improved. The volume scattering constant is obtained by S3H decomposition, and the dielectric constant is calculated through integration with Fresnel's formula. The snow density is calculated by the linear empirical formula of the snow permittivity and snow density. The snow density obtained by S3H decomposition can better reflect the distribution of the snow density in the research area, and the accuracy is higher than that obtained by the traditional spatial interpolation method. The EQeau model using the optimized method has greatly improved the accuracy of the SWE inversion results in the study area, and the total RMSE can reach 5.76 mm, which is consistent with the SWE observation at the watershed scale.
Our study proves that the EQeau model has a good effect when combined with an advanced C-band SAR GF-3 image. Compared with the results before optimization, the MRE decreased from 27.4% to 10.3%. In future studies, the Internet of Things and other advanced technologies should be combined to obtain rapid and timely watershed parameter information and achieve SWE observation in different research demonstration areas.