A method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy process

Abstract

The application discloses a kind of star-borne GNSS-R lake ice thickness inversion model construction methods of coupling freeze-thaw phase change energy process, method steps include: obtaining star-borne full polarization SMAP-R data, Tianmu-1 reflection measurement data, Fengyun-3 reflection measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical image, CIS ice map product and lake boundary vector data;Multi-source GNSS-R data and auxiliary data are preprocessed, and lake ice physics and dielectric parameter library are constructed;Lake ice coverage range and phenological date are obtained using the lake ice coverage range dynamic extraction and phenological event identification method based on response mechanism;Lake ice thermal growth / melting model is constructed, and the daily variation and seasonal evolution process of ice thickness are simulated;Key parameter sensitivity analysis and decoupling are carried out for ice thickness inversion;Ice thickness coupling inversion and precision verification are carried out.The application couples thermal growth / melting model and electromagnetic wave response model for the first time, constructs lake ice thickness physical inversion framework, and provides technical support for cold region water resource management.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary research direction of spaceborne GNSS reflection measurement and cryosphere science, and mainly involves a method for constructing a spaceborne GNSS-R lake ice thickness inversion model that couples the freeze-thaw phase change energy process. Background Technology

[0002] Of the world's 117 million lakes, more than half are intermittently or permanently covered by lake ice. Lake ice is not only an important component of the cryosphere but also an indicator of regional and even global climate. The seasonal occurrence of lake ice significantly alters the energy and heat exchange between lakes and the atmosphere, exacerbating climate change. Against the backdrop of climate change, the extent, phenological characteristics, and thickness of lake ice cover are changing rapidly, posing significant challenges to human activities and the lake ecosystems upon which they depend. Systematic monitoring of the seasonal and interannual variations in lake ice cover, phenology, thickness, and freeze-thaw states is crucial for understanding the mechanisms underlying climate change.

[0003] Conventional lake ice monitoring methods are limited by various factors (cost, spatial resolution, lighting / darkness conditions, etc.) and cannot be applied to small- to medium-sized lake ice scenarios with finer temporal and spatial resolutions. The emergence of spaceborne GNSS reflectance measurement technology has brought new developments to this challenge. Its high precision, all-weather, all-time high-resolution, and low-cost approach offers unique advantages for global lake ice monitoring and can fill gaps in global lake ice observation data. It is worth noting that while most studies on sea ice scenarios use a three-layer (air-ice-water) reflectance structure for sea ice thickness modeling, there are no reports of using spaceborne GNSS-R, especially multi-polarization spaceborne GNSS-R, for lake ice thickness estimation. However, the physical properties of lake ice (such as snow cover, ice porosity, and the complexity of the ice-water interface) differ significantly from those of sea ice, and directly applying a three-layer model may lead to insufficient inversion accuracy. Furthermore, existing models for sea ice thickness inversion suffer from severe ill-conditioning (multi-parameter coupling leads to non-unique solutions) and are heavily reliant on auxiliary products. Although thermodynamic models considering the energy balance at the ice surface (upper boundary condition) and electromagnetic coupling have been proposed to address this limitation, they lack an evolutionary description of the ice layer's thermal state. Specifically, they neglect the instantaneous balance of various energy fluxes at the ice surface (upper and / or lower boundary) at a specific moment, as well as the evolution of internal ice temperature with time and depth. This results in an imperfect physical constraint mechanism for coupling thermodynamic processes and electromagnetic responses, hindering high-precision inversion. Therefore, it is urgent to construct a lake ice thickness growth model coupled with freeze-thaw phase transition processes and establish a multi-polarization spaceborne GNSS-R lake ice thickness inversion paradigm driven by additional thermodynamic processes to achieve high-precision inversion of lake ice thickness with physical constraints.

[0004] In view of this, this invention proposes a method for constructing a spaceborne GNSS-R lake ice thickness inversion model that couples the freeze-thaw phase change energy process. For the first time, it couples the thermodynamic growth / ablation model with the electromagnetic wave response model to construct a physical inversion framework for lake ice thickness that is "driven by thermodynamic processes and constrained by electromagnetic observations". This method breaks through the ill-conditioned bottleneck of single-parameter inversion and the limitations of isolated parameter analysis, and provides technical support for the study of climate-cryosphere interaction and cold region water resource management. Summary of the Invention

[0005] To achieve the above-mentioned technical objectives and effects, the present invention provides the following technical solution:

[0006] This invention aims to propose a method for constructing a spaceborne GNSS-R lake ice thickness inversion model that couples the freeze-thaw phase transition energy process. This method constructs a lake ice thermodynamic growth / ablation model, establishes an electromagnetic-thermodynamic coupled inversion framework, performs dynamic response simulation of the phase transition interface, reveals the response mechanism of multipolar GNSS-R signals to the lake ice freeze-thaw process and its layered structure, and provides physical and theoretical support for subsequent parameter inversion.

[0007] This invention provides the following technical solution:

[0008] Step S1: Acquire spaceborne fully polarized SMAP-R data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data.

[0009] Step S2 involves preprocessing the multi-source GNSS-R data and auxiliary data, and constructing a library of lake ice physical and dielectric parameters.

[0010] Step S3: Use the response mechanism-based dynamic extraction and phenological event identification method to obtain the lake ice coverage and phenological dates.

[0011] Step S4: Construct a lake ice thermodynamic growth / ablation model to simulate the diurnal variation and seasonal evolution of ice thickness.

[0012] Step S5 involves performing sensitivity analysis and decoupling of key parameters for ice thickness inversion.

[0013] Step S6: Perform ice thickness coupling inversion and accuracy verification.

[0014] Preferably, the data preparation in step S1 includes the following:

[0015] Data preparation: Download the following data from the official website: SMAP-R full polarization satellite data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data.

[0016] Preferably, step S2, which involves preprocessing multi-source GNSS-R data and auxiliary data and constructing a parameter library, includes the following sub-steps:

[0017] S2.1 preprocesses multi-source GNSS-R data, including data cleaning, removing outlier DDMs, and removing land / coastline contamination data; feature extraction, extracting features such as SNR, reflectivity, peak DDM, and scattering cross section from DDMs of different polarization GNSS systems; geometric registration, establishing the lake-observation correspondence by matching the coordinates of the specular reflection point with the lake vector data; and incident angle normalization: establishing the variation law of SNR under different incident angles to provide input for subsequent models.

[0018] The parameter library in S2.2 mainly consists of two parts: basic parameters, including ice density, ice thermal conductivity, ice specific heat capacity, snow thermal conductivity, latent heat of melting, and ice surface albedo; and dielectric property parameters, including ice dielectric constant, snow dielectric constant, lake water dielectric constant, ice-water interface reflection constant, and ice layer attenuation coefficient.

[0019] According to the reflectivity formula, the total reflectivity is substituted into the formula to calculate the snow roughness:

[0020]

[0021]

[0022] The reflectance coefficient of the air-snow interface is... The reflectance of the snow-lake ice interface is The reflectance of the lake ice-lake water interface is In lossy media, wavenumber It is a complex number. , The complex wave number in the snow layer. is the complex wave number in the lake ice layer. The standard deviation of the snow surface roughness characterizes the snow scattering attenuation effect.

[0023] S.2.3 Based on daily true-color composite products from MODIS / Terra satellites, combined with Sentinel-2 optical imagery (under cloudless conditions) and CIS ice map products, label data for lake ice / open water surfaces is generated through visual interpretation.

[0024] Preferably, step S3, which uses a response mechanism-based method for dynamic extraction of lake ice coverage and identification of phenological events to obtain the lake ice coverage and phenological dates, includes the following sub-steps:

[0025] S3.1 Based on the advantages of GNSS-R data, the Random Forest (RF) method is used to separate water body and ice area data from a large amount of GNSS-R observation data as samples for subsequent fine classification.

[0026] S3.2 employs a text-encoded, multimodal GNSS-R lake ice fine classification model to perform refined identification of multiple states, including thin ice, thick ice, snow-covered ice, and mixed ice-water states. It defines a new continuous observation value for lake ice freeze-thaw cycles, providing a new quantitative indicator for the lake ice evolution process.

[0027]

[0028] In the formula, The seasonal scaling factor is based on the time of each grid cell. Sequence of lake surface reflectance Compared with seasonal reference freeze status and thawed state The relationship between time. The measured value, The definition is as follows:

[0029]

[0030] In the formula, the frozen and thawed reference states are calculated for each pixel in winter ( ) and summer ( Reflectance during the period ( The 5% and 95% quantiles of the cumulative distribution function. Different freeze-thaw states are related to the relevant analysis time periods.

[0031] S3.3 Based on the water body and ice zone data separated in S3.1, quantitative calculation of lake-scale coverage is performed; the lake is divided into 0.05°×0.05° grids, and multiple time-series observations within each grid are aggregated to obtain the ice condition probability at the grid scale; a lake ice cover distribution map is generated using the Kriging interpolation method, and a regression model is trained to output the daily lake ice cover area and coverage ratio.

[0032] S3.4 Automatic identification of lake ice phenology is performed based on the continuous curve of freezing intensity and the daily lake ice coverage ratio sequence to obtain key phenological index parameters (freezing day).

[0033] Preferably, step S4, which involves constructing the thermodynamic growth / ablation model, includes the following sub-steps:

[0034] S4.1 Initialize the thermodynamic model. The freezing date identified in step S3 is used as the start time of the thermodynamic model. Initial conditions are set as follows: ice thickness on the freezing date and snow thickness taken from ERA5 reanalysis data. Model time step is set as 1 hour (thermodynamic process) or 1 day.

[0035] S4.2 Solving for internal heat conduction, using a one-dimensional heat conduction equation to simulate the temperature evolution inside the ice layer:

[0036]

[0037] in, Shortwave radiation that penetrates ice decreases exponentially with depth:

[0038]

[0039] In the formula, The extinction coefficient is denoted as .

[0040] S4.3 Determining the upper boundary conditions based on ice surface energy balance calculations: This invention adopts the snow layer thermodynamic model (steady-state energy balance equation) proposed in existing research as the upper boundary condition, expressed as:

[0041]

[0042] In the formula, net shortwave radiation is: ,in (Downward shortwave radiation); Penetrating radiation Based on the parameterization of ice transparency, it is usually taken as 10-20% of the shortwave radiation; downward longwave radiation (From ERA5); Upward long-wave radiation ,in The ice surface temperature is to be determined. The emissivity of the ice surface is (~0.97). Stefan-Boltzmann constant; sensible heat flux ,in air density, The specific heat of air at constant pressure. Heat transfer coefficient; latent heat flux ,in To sublimate latent heat, The water vapor exchange coefficient, and Specific humidity of air and ice surface, respectively; thermal conductivity of ice layer. (Calculated from the temperature field inside the ice layer), solving this equation yields the ice surface temperature. .

[0043] S4.4 Determining the lower boundary conditions based on the heat flux at the ice-water interface requires considering the influence of the water body's heat storage on the ice bottom:

[0044]

[0045] In the formula, The density of water; This is the specific heat capacity of water; The heat exchange coefficient is taken as ~0.001. The velocity of the water flow under the ice (take an empirical value of 0.01-0.05 m / s); The subglacial water temperature is taken as ~0.5-2°C for shallow lakes and ~4°C for deep lakes. The freezing point temperature is taken as 0°C for freshwater lakes.

[0046] S4.5 Calculate the rate of ice thickness change based on the energy balance at the ice-water interface:

[0047]

[0048] in, For the latent heat of fusion (3.34 × 10⁻⁶) 5 J / kg). When the right end is positive, the ice bottom is melting; when it is negative, the ice bottom is freezing. Ice surface melting is determined by the ice surface temperature: when... At this time, ice surface melting occurs, and the melting rate is determined by the remaining energy. Snow layer evolution: Snow thickness is updated based on precipitation data, taking into account the snow compaction and melting processes.

[0049] Preferably, the sensitivity analysis and decoupling of key parameters for ice thickness inversion described in step S5 includes the following sub-steps:

[0050] S5.1 identifies and classifies core parameter sets such as geometric parameters (ice thickness, snow thickness, roughness), physical property parameters (ice dielectric, snow dielectric, liquid water content), and thermodynamic parameters (heat exchange coefficient, sub-ice flow velocity).

[0051] S5.2 Based on a three-layer simplified reflectivity model, a preliminary single-factor sensitivity analysis was conducted to assess the influence trends of incident angle, ice surface temperature, and snow depth on multipolar reflectivity.

[0052] Based on the above, S5.3 uses the Sobol global sensitivity analysis method based on the four-layer medium reflectivity coupling model to calculate the first-order sensitivity index and total effect index of each parameter to the multipolar observation, quantitatively evaluate the main effect and interaction effect of the parameters, and identify the dominant parameters and negligible parameters.

[0053] S5.4 For the strongly coupled parameter combinations revealed by sensitivity analysis, a multi-level decoupling strategy is adopted:

[0054] 1) Decoupling of multipolarization information utilizes the sensitivity differences of different polarization channels to separate parameter contributions;

[0055] 2) Multi-angle observation constraints combined with response differences at different incident angles enhance parameter identifiability;

[0056] 3) Temporal physical constraints introduce ice thickness growth smoothness as a regularization term to suppress abnormal fluctuations;

[0057] 4) External data assistance uses ERA5 snow depth products as a strong constraint on snow thickness to reduce parameter coupling.

[0058] The final result is a parameter sensitivity map (quantifying the contribution weight of each parameter), a parameter optimization scheme (clarifying the parameters to be inverted, fixed, or input), and an inversion parameter configuration table (providing the prior distribution and constraint range), which provides a complete theoretical basis for parameter configuration and uncertainty quantification of the coupled inversion framework.

[0059] Preferably, the verification of the lake ice coupling inversion method described in step S6 includes the following steps:

[0060] To comprehensively verify the accuracy and applicability of the coupled inversion method, a triple-configuration algorithm was employed for evaluation. This invention uses a triple-configuration dataset composed of Tianmu-1 / FY-3 (FY-3) L1 GNSS-R data for ice thickness inversion and other product lake ice thickness datasets (i.e., ERA5, MODIS, ESA Lakes_cci, and ground-based measurements). Triple configuration comparisons were conducted: thermal model only, direct GNSS-R inversion only, and coupled inversion (the method of this invention). Uncertainty quantification was achieved through a triple strategy: Monte Carlo simulation (input error propagation), multi-model ensemble (model structure error), and Bayesian inversion (parameter posterior distribution). Results verification was conducted from three dimensions: time series comparison, spatial distribution comparison, and stratified error statistics (by ice thickness interval and lake type). The reliability, robustness, and universality of the method were systematically evaluated.

[0061] Compared with the prior art, the beneficial effects of the present invention are:

[0062] 1. For the first time, a thermodynamic growth / ablation model was coupled with an electromagnetic wave response model to construct a physical inversion framework for lake ice thickness driven by thermodynamic processes and constrained by electromagnetic observations, providing technical support for overcoming the bottleneck of single-parameter inversion and the limitations of isolated parameter analysis. Attached Figure Description

[0063] To more clearly illustrate the technical solution of the present invention, the accompanying drawings used in the description of the specific embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0064] This accompanying drawing is a schematic diagram of functional and connection relationships, and is for illustrative purposes only, not restrictive. The shape, relative proportions, and arrangement of components can be adjusted for engineering purposes without departing from the claims.

[0065] Figure 1 This is a flowchart illustrating an embodiment of the present invention.

[0066] Figure 2 This is a schematic diagram of the lake ice thermodynamic growth / ablation model in an embodiment of the present invention.

[0067] Figure 3 This is a flowchart of the electromagnetic-thermodynamic coupling inversion framework in an embodiment of the present invention. Detailed Implementation

[0068] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0069] The invention will now be further described with reference to the accompanying drawings.

[0070] To verify the feasibility of the method proposed in this invention, spaceborne fully polarized SMAP-R data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, MODIS data, ERA5 data, Sentinel-2 optical images, CIS ice map products, and lake boundary vector data were collected and selected from typical lakes including Qinghai Lake, Namtso (large deep-water lake), Hulun Lake (shallow-water lake), Qinghai-Tibet Plateau lakes (high-altitude snow-covered), and the Finnish Lake Group (medium and small-sized lakes) to conduct ice thickness monitoring experiments.

[0071] Example 1

[0072] A method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase transition energy processes, as shown in the appendix. Figure 1 This includes the following steps:

[0073] Step S1: Acquire spaceborne fully polarized SMAP-R data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data.

[0074] Step S2 involves preprocessing the multi-source GNSS-R data and auxiliary data, and constructing a library of lake ice physical and dielectric parameters.

[0075] Step S3: Use the response mechanism-based dynamic extraction and phenological event identification method to obtain the lake ice coverage and phenological dates.

[0076] Step S4: Construct a lake ice thermodynamic growth / ablation model to simulate the diurnal variation and seasonal evolution of ice thickness.

[0077] Step S5 involves performing sensitivity analysis and decoupling of key parameters for ice thickness inversion.

[0078] Step S6: Perform ice thickness coupling inversion and accuracy verification.

[0079] Preferably, the data preparation in step S1 includes the following:

[0080] Data preparation: Download the following data from the official website: SMAP-R full polarization satellite data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data.

[0081] Preferably, step S2, which involves preprocessing multi-source GNSS-R data and auxiliary data and constructing a parameter library, includes the following sub-steps:

[0082] S2.1 preprocesses multi-source GNSS-R data, including data cleaning, removing outlier DDMs, and removing land / coastline contamination data; feature extraction, extracting features such as SNR, reflectivity, peak DDM, and scattering cross section from DDMs of different polarization GNSS systems; geometric registration, establishing the lake-observation correspondence by matching the coordinates of the specular reflection point with the lake vector data; and incident angle normalization: establishing the variation law of SNR under different incident angles to provide input for subsequent models.

[0083] The parameter library in S2.2 mainly consists of two parts: basic parameters, including ice density, ice thermal conductivity, ice specific heat capacity, snow thermal conductivity, latent heat of melting, and ice surface albedo; and dielectric property parameters, including ice dielectric constant, snow dielectric constant, lake water dielectric constant, ice-water interface reflection constant, and ice layer attenuation coefficient.

[0084] According to the reflectivity formula, the total reflectivity is substituted into the formula to calculate the snow roughness:

[0085]

[0086]

[0087] The reflectance of the air-snow interface is... The reflectance of the snow-lake ice interface is The reflectance of the lake ice-lake water interface is In lossy media, wavenumber It is a complex number. , The complex wave number in the snow layer. is the complex wave number in the lake ice layer. The standard deviation of the snow surface roughness characterizes the snow scattering attenuation effect.

[0088] S.2.3 Based on daily true-color composite products from MODIS / Terra satellites, combined with Sentinel-2 optical imagery (under cloudless conditions) and CIS ice map products, label data for lake ice / open water surfaces is generated through visual interpretation.

[0089] Preferably, step S3, which uses a response mechanism-based method for dynamic extraction of lake ice coverage and identification of phenological events to obtain the lake ice coverage and phenological dates, includes the following sub-steps:

[0090] S3.1 Based on the advantages of GNSS-R data, the Random Forest (RF) method is used to separate water body and ice area data from a large amount of GNSS-R observation data as samples for subsequent fine classification.

[0091] S3.2 employs a text-encoded, multimodal GNSS-R lake ice fine classification model to perform refined identification of multiple states, including thin ice, thick ice, snow-covered ice, and mixed ice-water states. It defines a new continuous observation value for lake ice freeze-thaw cycles, providing a new quantitative indicator for the lake ice evolution process.

[0092]

[0093] In the formula, The seasonal scaling factor is based on the time of each grid cell. Sequence of lake surface reflectance Compared with seasonal reference freeze status and thawed state The relationship between time. The measured value, The definition is as follows:

[0094]

[0095] In the formula, the frozen and thawed reference states are calculated for each pixel in winter ( ) and summer ( Reflectance during the period ( The 5% and 95% quantiles of the cumulative distribution function. Different freeze-thaw states are related to the relevant analysis time periods.

[0096] S3.3 Based on the water body and ice zone data separated in S3.1, quantitative calculation of lake-scale coverage is performed; the lake is divided into 0.05°×0.05° grids, and multiple time-series observations within each grid are aggregated to obtain the ice condition probability at the grid scale; a lake ice cover distribution map is generated using the Kriging interpolation method, and a regression model is trained to output the daily lake ice cover area and coverage ratio.

[0097] S3.4 Automatic identification of lake ice phenology is performed based on the continuous curve of freezing intensity and the daily lake ice coverage ratio sequence to obtain key phenological index parameters (freezing day).

[0098] Preferably, step S4, which involves constructing the thermodynamic growth / ablation model, includes the following sub-steps:

[0099] S4.1 Initialize the thermodynamic model. The freezing date identified in step S3 is used as the start time of the thermodynamic model. Initial conditions are set as follows: ice thickness on the freezing date and snow thickness taken from ERA5 reanalysis data. Model time step is set as 1 hour (thermodynamic process) or 1 day.

[0100] S4.2 Solving for internal heat conduction, using a one-dimensional heat conduction equation to simulate the temperature evolution inside the ice layer:

[0101]

[0102] in, Shortwave radiation that penetrates ice decreases exponentially with depth:

[0103]

[0104] In the formula, The extinction coefficient is denoted as .

[0105] S4.3 Determining the upper boundary conditions based on ice surface energy balance calculations: This invention adopts the snow layer thermodynamic model (steady-state energy balance equation) proposed in existing research as the upper boundary condition, expressed as:

[0106]

[0107] In the formula, net shortwave radiation is: ,in (Downward shortwave radiation); penetrating radiation Based on the parameterization of ice transparency, it is usually taken as 10-20% of the shortwave radiation; downward longwave radiation (From ERA5); Upward long-wave radiation ,in The ice surface temperature is to be determined. The emissivity of the ice surface is (~0.97). Stefan-Boltzmann constant; sensible heat flux ,in air density, The specific heat of air at constant pressure. Heat transfer coefficient; latent heat flux ,in To sublimate latent heat, The water vapor exchange coefficient, and Specific humidity of air and ice surface, respectively; thermal conductivity of ice layer. (Calculated from the temperature field inside the ice layer), solving this equation yields the ice surface temperature. .

[0108] S4.4 Determining the lower boundary conditions based on the heat flux at the ice-water interface requires considering the influence of the water body's heat storage on the ice bottom:

[0109]

[0110] In the formula, The density of water; This is the specific heat capacity of water; The heat exchange coefficient is taken as ~0.001. The velocity of the water flow under the ice (take an empirical value of 0.01-0.05 m / s); The subglacial water temperature is taken as ~0.5-2°C for shallow lakes and ~4°C for deep lakes. The freezing point temperature is taken as 0°C for freshwater lakes.

[0111] S4.5 Calculate the rate of ice thickness change based on the energy balance at the ice-water interface:

[0112]

[0113] in, For the latent heat of fusion (3.34 × 10⁻⁶) 5 J / kg). When the right end is positive, the ice bottom is melting; when it is negative, the ice bottom is freezing. Ice surface melting is determined by the ice surface temperature: when... At this time, ice surface melting occurs, and the melting rate is determined by the remaining energy. Snow layer evolution: Snow thickness is updated based on precipitation data, taking into account the snow compaction and melting processes.

[0114] Preferably, the sensitivity analysis and decoupling of key parameters for ice thickness inversion described in step S5 includes the following sub-steps:

[0115] S5.1 identifies and classifies core parameter sets such as geometric parameters (ice thickness, snow thickness, roughness), physical property parameters (ice dielectric, snow dielectric, liquid water content), and thermodynamic parameters (heat exchange coefficient, sub-ice flow velocity).

[0116] S5.2 Based on a three-layer simplified reflectivity model, a preliminary single-factor sensitivity analysis was conducted to assess the influence trends of incident angle, ice surface temperature, and snow depth on multipolar reflectivity.

[0117] Based on the above, S5.3 uses the Sobol global sensitivity analysis method based on the four-layer medium reflectivity coupling model to calculate the first-order sensitivity index and total effect index of each parameter to the multipolar observation, quantitatively evaluate the main effect and interaction effect of the parameters, and identify the dominant parameters and negligible parameters.

[0118] S5.4 For the strongly coupled parameter combinations revealed by sensitivity analysis, a multi-level decoupling strategy is adopted:

[0119] 1) Decoupling of multipolarization information utilizes the sensitivity differences of different polarization channels to separate parameter contributions;

[0120] 2) Multi-angle observation constraints combined with response differences at different incident angles enhance parameter identifiability;

[0121] 3) Temporal physical constraints introduce ice thickness growth smoothness as a regularization term to suppress abnormal fluctuations;

[0122] 4) External data assistance uses ERA5 snow depth products as a strong constraint on snow thickness to reduce parameter coupling.

[0123] The final result is a parameter sensitivity map (quantifying the contribution weight of each parameter), a parameter optimization scheme (clarifying the parameters to be inverted, fixed, or input), and an inversion parameter configuration table (providing the prior distribution and constraint range), which provides a complete theoretical basis for parameter configuration and uncertainty quantification of the coupled inversion framework.

[0124] Preferably, the verification of the lake ice coupling inversion method described in step S6 includes the following steps:

[0125] To comprehensively verify the accuracy and applicability of the coupled inversion method, a triple-configuration algorithm was employed for evaluation. This invention uses a triple-configuration dataset composed of Tianmu-1 / FY-3 (FY-3) L1 GNSS-R data for ice thickness inversion and other product lake ice thickness datasets (i.e., ERA5, MODIS, ESA Lakes_cci, and ground-based measurements). Triple configuration comparisons were conducted: thermal model only, direct GNSS-R inversion only, and coupled inversion (the method of this invention). Uncertainty quantification was achieved through a triple strategy: Monte Carlo simulation (input error propagation), multi-model ensemble (model structure error), and Bayesian inversion (parameter posterior distribution). Results verification was conducted from three dimensions: time series comparison, spatial distribution comparison, and stratified error statistics (by ice thickness interval and lake type). The reliability, robustness, and universality of the method were systematically evaluated.

[0126] Appendix Figure 3 The invention provides a detailed schematic diagram of the construction of the electromagnetic-thermodynamic coupled inversion framework. The proposed solution solves the ill-conditioned bottleneck of single-parameter inversion and the limitations of isolated parameter analysis, providing a solid physical foundation and theoretical methodological innovation for multi-parameter collaborative inversion of lake ice.

[0127] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0128] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase transition energy processes, characterized in that, Includes the following steps: Step S1: Acquire spaceborne fully polarized SMAP-R data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data. Step S2 involves preprocessing the multi-source GNSS-R data and auxiliary data, and constructing a library of lake ice physical and dielectric parameters. Step S3: Use the response mechanism-based dynamic extraction and phenological event identification method to obtain the lake ice coverage and phenological dates. Step S4: Construct a lake ice thermodynamic growth / ablation model to simulate the diurnal variation and seasonal evolution of ice thickness. Step S5 involves performing sensitivity analysis and decoupling of key parameters for ice thickness inversion. Step S6: Perform ice thickness coupling inversion and accuracy verification.

2. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy processes according to claim 1, characterized in that: Step S1 includes the following: Data preparation: Download the following data from the official website: SMAP-R full polarization satellite data, Tianmu-1 reflectance measurement data, Fengyun-3 reflectance measurement data, ERA5-Land data, MODIS data, Sentinel-2 optical imagery, CIS ice map products, and lake boundary vector data.

3. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase transition energy processes according to claim 1, characterized in that: Step S2 includes the following: S2.1 preprocesses multi-source GNSS-R data, including data cleaning, removal of outlier DDMs, and removal of land / coastline contamination data; Feature extraction: Extracting features such as SNR, reflectivity, peak DDM, and scattering cross section from DDMs of different polarization GNSS systems; Geometric registration: Based on the matching of the coordinates of the mirror reflection point with the lake vector data, a lake-observation correspondence is established. Incident angle normalization: The variation law of SNR under different incident angles is established to provide input for subsequent models. The parameter library in S2.2 mainly consists of two parts: basic parameters, including ice density, ice thermal conductivity, ice specific heat capacity, snow thermal conductivity, latent heat of melting, and ice surface albedo; and dielectric property parameters, including ice dielectric constant, snow dielectric constant, lake water dielectric constant, ice-water interface reflection constant, and ice layer attenuation coefficient. According to the reflectivity formula, the total reflectivity is substituted into the formula to calculate the snow roughness: ； ； The reflectance coefficient of the air-snow interface is... The reflectance of the snow-lake ice interface is The reflectance of the lake ice-lake water interface is In lossy media, wavenumber It is a complex number. , The complex wave number in the snow layer. is the complex wave number in the lake ice layer. The standard deviation of the snow surface roughness characterizes the snow scattering attenuation effect. S.2.3 Based on daily true-color composite products from MODIS / Terra satellites, combined with Sentinel-2 optical imagery (under cloudless conditions) and CIS ice map products, label data for lake ice / open water surfaces is generated through visual interpretation.

4. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy processes according to claim 1, characterized in that: Step S3 includes the following: S3.1 Based on the advantages of GNSS-R data, the Random Forest (RF) method is used to separate water body and ice area data from a large amount of GNSS-R observation data as samples for subsequent fine classification. S3.2 employs a text-encoded, multimodal GNSS-R lake ice fine classification model for refined identification of various states, including thin ice, thick ice, snow-covered ice, and ice-water mixture. It defines a new continuous observation value for lake ice freeze-thaw cycles, providing a new quantitative indicator for the lake ice evolution process. ； In the formula, The seasonal scaling factor is based on the time of each grid cell. Sequence of lake surface reflectance Compared with seasonal reference freeze status and thawed state The relationship between time. The measured value, The definition is as follows: ； In the formula, the frozen and thawed reference states are calculated for each pixel in winter ( ) and summer ( Reflectance during the period ( The 5% and 95% quantiles of the cumulative distribution function. Different freeze-thaw states are related to the relevant analysis time periods. S3.3 Based on the water body and ice zone data separated in S3.1, quantitative calculation of lake-scale coverage is performed; the lake is divided into 0.05°×0.05° grids, and multiple time-series observations within each grid are aggregated to obtain the ice condition probability at the grid scale; a lake ice cover distribution map is generated using the Kriging interpolation method, and a regression model is trained to output the daily lake ice cover area and coverage ratio. S3.4 Automatic identification of lake ice phenology is performed based on the continuous curve of freezing intensity and the daily lake ice coverage ratio sequence to obtain key phenological index parameters (freezing day).

5. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy processes according to claim 1, characterized in that: Step S4: Consider the instantaneous equilibrium relationship of various energy fluxes on the ice surface (or the upper boundary of the ice layer) at a specific moment, including the following: S4.1 Initialize the thermodynamic model. The freezing date identified in step S3 is used as the start time for the thermodynamic model. Initial conditions are set as follows: ice thickness on the freezing date... , the snow is thick Data taken from ERA5 reanalysis; model time step set: 1 hour (thermodynamic process) or 1 day. S4.2 Solving for internal heat conduction, using a one-dimensional heat conduction equation to simulate the temperature evolution inside the ice layer: ； in, Shortwave radiation that penetrates ice decreases exponentially with depth: ； In the formula, The extinction coefficient is denoted as . S4.3 Determining the upper boundary conditions based on ice surface energy balance calculations: This invention adopts the snow layer thermodynamic model (steady-state energy balance equation) proposed in existing research as the upper boundary condition, expressed as: ； In the formula, net shortwave radiation is: ,in (Downward shortwave radiation); penetrating radiation Based on the parameterization of ice transparency, it is usually taken as 10-20% of the shortwave radiation; downward longwave radiation (From ERA5); Upward long-wave radiation ,in The ice surface temperature is to be determined. The emissivity of the ice surface is (~0.97). Stefan-Boltzmann constant; sensible heat flux ,in air density, The specific heat of air at constant pressure. Heat transfer coefficient; latent heat flux ,in To sublimate latent heat, The water vapor exchange coefficient, and Specific humidity of air and ice surface, respectively; thermal conductivity of ice layer. (Calculated from the temperature field inside the ice layer), solving this equation yields the ice surface temperature. . S4.4 Determining the lower boundary conditions based on the heat flux at the ice-water interface requires considering the influence of the water body's heat storage on the ice bottom: ； In the formula, The density of water; This is the specific heat capacity of water; The heat exchange coefficient is taken as ~0.

001. The velocity of the water flow under the ice (take an empirical value of 0.01-0.05 m / s); The subglacial water temperature is taken as ~0.5-2°C for shallow lakes and ~4°C for deep lakes. The freezing point temperature is taken as 0°C for freshwater lakes. S4.5 Calculate the rate of ice thickness change based on the energy balance at the ice-water interface: ； in, For the latent heat of fusion (3.34 × 10⁻⁶) 5 J / kg). When the right end is positive, the ice bottom is melting; when it is negative, the ice bottom is freezing. Ice surface melting is determined by the ice surface temperature: when... At this time, ice surface melting occurs, and the melting rate is determined by the remaining energy. Snow layer evolution: Snow thickness is updated based on precipitation data, taking into account the snow compaction and melting processes.

6. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy processes according to claim 1, characterized in that: Step S5: This invention conducts key parameter analysis and decoupling, employing a progressive sensitivity analysis strategy, including the following: S5.1 identifies and classifies core parameter sets such as geometric parameters (ice thickness, snow thickness, roughness), physical property parameters (ice dielectric, snow dielectric, liquid water content), and thermodynamic parameters (heat exchange coefficient, sub-ice flow velocity). S5.2 Based on a three-layer simplified reflectivity model, a preliminary single-factor sensitivity analysis was conducted to assess the influence trends of incident angle, ice surface temperature, and snow depth on multipolar reflectivity. Based on the above, S5.3 uses the Sobol global sensitivity analysis method based on the four-layer medium reflectivity coupling model to calculate the first-order sensitivity index and total effect index of each parameter to the multipolar observation, quantitatively evaluate the main effect and interaction effect of the parameters, and identify the dominant parameters and negligible parameters. S5.4 For the strongly coupled parameter combinations revealed by sensitivity analysis, a multi-level decoupling strategy is adopted: 1) Decoupling of multipolarization information utilizes the sensitivity differences of different polarization channels to separate parameter contributions; 2) Multi-angle observation constraints combined with response differences at different incident angles enhance parameter identifiability; 3) Temporal physical constraints introduce ice thickness growth smoothness as a regularization term to suppress abnormal fluctuations; 4) External data assistance uses ERA5 snow depth products as a strong constraint on snow thickness to reduce parameter coupling. The final result is a parameter sensitivity map (quantifying the contribution weight of each parameter), a parameter optimization scheme (clarifying the parameters to be inverted, fixed, or input), and an inversion parameter configuration table (providing the prior distribution and constraint range), which provides a complete theoretical basis for parameter configuration and uncertainty quantification of the coupled inversion framework.

7. The method for constructing a spaceborne GNSS-R lake ice thickness inversion model coupled with freeze-thaw phase change energy processes according to claim 1, characterized in that: Step S6 includes the following: To comprehensively verify the accuracy and applicability of the coupled inversion method, a triple-configuration algorithm was employed for evaluation. This invention uses a triple-configuration dataset composed of Tianmu-1 / FY-3 (FY-3) L1 GNSS-R data for ice thickness inversion and other product lake ice thickness datasets (i.e., ERA5, MODIS, ESA Lakes_cci, and ground-based measurements). Triple configuration comparisons were conducted: thermal model only, direct GNSS-R inversion only, and coupled inversion (the method of this invention). Uncertainty quantification was achieved through a triple strategy: Monte Carlo simulation (input error propagation), multi-model ensemble (model structure error), and Bayesian inversion (parameter posterior distribution). Results verification was conducted from three dimensions: time series comparison, spatial distribution comparison, and stratified error statistics (by ice thickness interval and lake type). The reliability, robustness, and universality of the method were systematically evaluated.

Images