InSAR-Coupled ERT Method and System for Inversion of Hydrothermal Conditions within Frozen Soil

Abstract

The present invention discloses a method and system for inverting the hydrothermal characteristics inside frozen soil by coupling InSAR and ERT, which relates to the technical field of inverting the hydrothermal characteristics inside frozen soil. The method includes: separating the freeze-thaw deformation data according to the InSAR surface deformation data; determining the depth of the dynamic freeze-thaw interface according to multi-period ERT data and inversely calculating the first water content; correcting the thickness of the active layer in the frost heaving theoretical model according to the depth of the dynamic freeze-thaw interface to obtain a corrected frost heaving theoretical model; inversely calculating the water content based on the corrected frost heaving theoretical model according to the freeze-thaw deformation data to obtain the second water content; optimizing the hydrothermal coupling numerical model by using the first water content; inputting the InSAR surface deformation data and the surface layer temperature into the optimized dynamic hydrothermal coupling numerical model to output the temperature and the third water content at different depths and different times inside the frozen soil in the study area. Finally, the high-precision dynamic inversion of the hydrothermal characteristics inside the frozen soil at the regional scale is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of inversion of hydrothermal characteristics inside frozen soil, and particularly to a method and system for inverting hydrothermal characteristics inside frozen soil by coupling InSAR and ERT. Background Art

[0002] The hydrothermal characteristics inside frozen soil are the key factors for revealing the degradation process of frozen soil under climate change and predicting the evolution of frozen soil disasters. High-precision inversion of the hydrothermal characteristics inside frozen soil is of great significance. Currently, the main inversion methods include on-site measurement method, model simulation method and satellite remote sensing inversion method: The on-site measurement method obtains single-point vertical hydrothermal profile data through equipment such as borehole temperature measurement and time domain reflectometer (TDR). Although the single-point detection accuracy can reach centimeter-level resolution, to obtain regional data, it is necessary to densely deploy and consume a large amount of resource costs; The model simulation method is mainly based on the law of conservation of matter and energy, mainly using land process models to simulate the energy cycle process of the atmosphere-vegetation-soil. Although it can obtain hydrothermal product data above the kilometer level, the spatial resolution is low and the data update cycle is long, which cannot meet the application requirements; Satellite remote sensing inversion methods, such as thermal infrared remote sensing and passive microwave remote sensing, can achieve large-scale surface water and heat monitoring at the spatial scale, but their detection depth is limited to the surface (<10 cm), and the hydrothermal inversion of the active layer of frozen soil (1-3 m inside) is completely ineffective. Therefore, the method for obtaining the hydrothermal characteristics inside frozen soil at the regional scale has not been effectively solved yet.

[0003] In view of this, the present application is specifically proposed. Summary of the Invention

[0004] The technical problem to be solved by the present invention is that the existing methods for inverting hydrothermal characteristics inside frozen soil have problems such as low inversion accuracy and poor universality. The purpose of the present invention is to provide a method and system for inverting hydrothermal characteristics inside frozen soil by coupling InSAR and ERT. By monitoring the surface deformation of InSAR and combining the frost heaving theory to obtain the water content, and obtaining the resistivity vertical distribution function and the dynamic value of the freeze-thaw interface depth through ERT resistivity tomography, a hydrothermal coupling numerical model considering time-varying characteristics is constructed, and finally, high-precision dynamic inversion of the hydrothermal characteristics inside frozen soil at the regional scale is realized, providing key data and solutions for frozen soil engineering safety, frozen soil disaster warning, etc.

[0005] Among them, InSAR (Interferometry Synthetic Aperture Radar) refers to synthetic aperture radar interferometry, and ERT (Electrical Resistance Tomography) refers to electrical resistance tomography technology.

[0006] The present invention is realized through the following technical solutions:

[0007] In a first aspect, the present invention provides a method for inverting the hydrothermal state inside frozen soil by coupling InSAR and ERT, the method comprising:

[0008] Obtaining InSAR surface deformation data, ERT data, and surface layer temperature of the same study area;

[0009] According to the InSAR surface deformation data, based on an improved InSAR frozen soil deformation separation model, separating the freeze-thaw deformation data;

[0010] According to multi-period ERT data, determining the depth of the dynamic freeze-thaw interface and inverting the first water content;

[0011] According to the depth of the dynamic freeze-thaw interface, correcting the thickness of the active layer in the frost heave theoretical model to obtain a corrected frost heave theoretical model;

[0012] According to the freeze-thaw deformation data, based on the corrected frost heave theoretical model, performing water content inversion to obtain the second water content;

[0013] Taking the second water content and the surface layer temperature as initial conditions, taking the depth of the dynamic freeze-thaw interface as a boundary condition, and inputting them into a hydrothermal coupling numerical model; and optimizing the hydrothermal coupling numerical model using the first water content to obtain an optimized dynamic hydrothermal coupling numerical model;

[0014] Inputting the InSAR surface deformation data and the surface layer temperature into the optimized dynamic hydrothermal coupling numerical model, and outputting the temperature and the third water content at different depths and different times inside the frozen soil in the study area.

[0015] Furthermore, the InSAR surface deformation data is surface deformation data obtained by processing SAR image data and DEM data based on InSAR time series interferometry technology;

[0016] The ERT data is obtained by arranging a resistivity tomography imaging system at typical points in the study area to obtain multi-period resistivity profile data ρ(z,t) covering the whole year;

[0017] The surface layer temperature is time series data Ts(t) of the surface temperature obtained through remote sensing product data, meteorological data, or ground monitoring stations.

[0018] Furthermore, separating the freeze-thaw deformation data based on the improved InSAR frozen soil deformation separation model includes:

[0019] According to multi-period SAR image data and DEM data, after image registration and terrain phase removal, selecting an interference pair by preset time baseline and spatial baseline thresholds to generate a differential interference image;

[0020] According to the differential interference image, selecting high coherence points as the target point set;

[0021] Using the minimum cost flow algorithm, perform phase unwrapping on the target point set to obtain the unwrapped phase;

[0022] According to the unwrapped phase, combine the internal hydrothermal factors to simulate the seasonal phase and improve the traditional frozen soil deformation separation model to obtain an improved InSAR frozen soil deformation separation model;

[0023] Based on the improved InSAR frozen soil deformation separation model, separate the seasonal deformation of the frozen soil and the deformation of the slope itself; and convert the seasonal deformation of the frozen soil into freeze-thaw deformation data.

[0024] Furthermore, the freeze-thaw deformation data is the surface freeze-thaw deformation amount d(t) that changes with time in the study area, where d(t)=[d(t 1 ),d(t 2 ),…,d(t N )], t 1 ~t N is the time series, and d(t N ) is the surface freeze-thaw deformation at the t N th time.

[0025] Furthermore, according to multi-period ERT data, determine the dynamic freeze-thaw interface depth and invert the first water content, including:

[0026] According to multi-period ERT data, determine the dynamic freeze-thaw interface depth at different times by identifying the resistivity gradient extreme values;

[0027] And in combination with Archie's law, invert the first water content θ ERT (z,t) at different depths and different times.

[0028] Furthermore, the second water content is the average soil water content corresponding to the freezing depth or the melting depth.

[0029] Furthermore, based on the modified frost heave theoretical model, perform water content inversion to obtain the second water content, including:

[0030] The water content inversion formula during the freezing period is:

[0031] ;

[0032] The water content inversion formula during the melting period is:

[0033] ;

[0034] In the formula, is the average soil water content at the time t when the freezing depth or the melting depth changes , that is, the second water content; is the residual water content after soil freezing; is the temperature correction coefficient; is the soil freezing temperature; is the surface temperature; is the change in freeze-thaw deformation, equal to d ( t ) - d( t -1), d ( t ) is the surface freeze-thaw deformation at time t , d ( t -1) is the freeze-thaw deformation at the previous time t of time t -1; is the ice resistivity; is the resistivity of water; is the change in the freezing or thawing depth of the soil layer, equal to z f ( t ) - z f ( t -1), z f ( t ) is the freeze-thaw interface depth at time t , z f ( t -1) is the freeze-thaw interface depth at the previous time t of time t -1, provided by ERT data.

[0035] Furthermore, the optimization steps of the hydrothermal coupling numerical model are as follows:

[0036] Based on the prior parameters, generate N parameter sets ; is the parameter vector to be optimized, i is a positive integer, is the effective soil thermal conductivity parameter, is the soil volume heat capacity, is the water permeability, is the residual unfrozen water content, is the expansion coefficient; T is the transpose of the vector;

[0037] For each parameter set , obtain the predicted water content at depth z at time t through the hydrothermal coupling numerical model ;

[0038] Arrange the first water content obtained by inversion according to time tGenerate the observation vector , θ ERT ( z m , t ) is the time inverted by ERT t The first water content at depth z m at;

[0039] Construct an objective function based on the predicted water content value and the observation vector Calculate the global residual sum of squares, N t is the observation time, N z is the maximum freezing or thawing thickness; θ ERT ( z , t ) is the time inverted by ERT t The first water content at depth z;

[0040] According to the global residual sum of squares, obtain the minimum residual through gradient calculation to update the hydrothermal characteristic parameters, and obtain an optimized dynamic hydrothermal coupling numerical model.

[0041] The above technical solution corrects the predicted water content value based on the observation vector and optimizes the hydrothermal coupling numerical model.

[0042] In a second aspect, the present invention further provides a system for inverting the internal hydrothermal state of frozen soil by coupling InSAR and ERT. The system includes:

[0043] An acquisition unit for acquiring InSAR surface deformation data, ERT data, and surface layer temperature of the same study area;

[0044] A separation unit for separating freeze-thaw deformation data based on an improved InSAR frozen soil deformation separation model according to the InSAR surface deformation data;

[0045] A determination unit for determining the depth of the dynamic freeze-thaw interface based on multi-period ERT data and inversely calculating the first water content;

[0046] A frost heave model correction unit for correcting the active layer thickness in the frost heave theoretical model according to the depth of the dynamic freeze-thaw interface to obtain a corrected frost heave theoretical model;

[0047] A water content inversion unit for inversely calculating the second water content based on the corrected frost heave theoretical model according to the freeze-thaw deformation data;

[0048] A hydrothermal model optimization unit is used to input a coupled hydrothermal numerical model with the second water content and the surface layer temperature as the initial conditions and the depth of the dynamic freeze-thaw interface as the boundary conditions; and optimize the coupled hydrothermal numerical model using the first water content to obtain an optimized dynamic coupled hydrothermal numerical model.

[0049] A frozen soil internal hydrothermal inversion unit is used to input InSAR surface deformation data and surface layer temperature into the optimized dynamic coupled hydrothermal numerical model and output the temperature and the third water content at different depths and different times inside the frozen soil in the study area.

[0050] Furthermore, a determination unit includes:

[0051] A determination subunit is used to determine the depth of the dynamic freeze-thaw interface at different times by identifying the extreme values of the resistivity gradient based on multi-period ERT data.

[0052] An inversion subunit is used to inversely obtain the first water content at different depths and different times in combination with Archie's law.

[0053] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0054] 1. The method and system for inverting the internal hydrothermal properties of frozen soil by coupling InSAR and ERT of the present invention obtain the water content through InSAR surface deformation monitoring and in combination with the frost heaving theory, obtain the resistivity vertical distribution function and the dynamic value of the freeze-thaw interface depth through ERT resistivity tomography, construct a coupled hydrothermal numerical model considering time-varying characteristics, and finally realize the high-precision dynamic inversion of the internal hydrothermal properties of frozen soil at the regional scale, with good universality, providing key data and solutions for the safety of frozen soil engineering, early warning of frozen soil disasters, etc.

[0055] 2. The present invention uses high-resolution SAR images such as Sentinel-1, LT-1, and GF-3, and obtains seasonal surface deformation through time-series InSAR time-series interferometry technology, and then inversely obtains the water content with high spatio-temporal resolution (≤12 days, ≤30 m), replacing traditional meteorological station interpolation and traditional remote sensing data.

[0056] 3. Traditional remote sensing means can only inversely obtain the water content on the surface (<10 cm). The present invention constructs a coupled hydrothermal model considering time-varying characteristics by combining InSAR and ERT data, realizes the dynamic calibration of the freeze-thaw interface depth z f (t), real-time tracks the freezing and melting processes of frozen soil, and dynamically inversely obtains the temperature and water content parameters at different depths inside the frozen soil.

[0057] 4. The soil hydrothermal parameters input into the traditional hydrothermal model (parameters such as the effective thermal conductivity of the soil and the volumetric heat capacity of the soil etc. often rely on empirical values. Through in-situ observation combined with ERT, the present invention can obtain high-resolution (<0.1m) vertical water content, optimize the parameters of the traditional hydrothermal coupling model, and further improve the accuracy of the hydrothermal model;

[0058] 5. By integrating multi-source data such as InSAR, ERT, and in-situ observation, and combining time series analysis and the hydrothermal coupling numerical model, the present invention can provide technical solutions and effective data for inverting the internal hydrothermal field at the regional level. Brief Description of the Drawings

[0059] The drawings described herein are used to provide a further understanding of the embodiments of the present invention, form a part of this application, and do not limit the embodiments of the present invention. In the drawings:

[0060] Figure 1 is the flow chart of the method for inverting the internal hydrothermal of frozen soil by coupling InSAR and ERT of the present invention;

[0061] Figure 2 is the detailed flow chart of the method for inverting the internal hydrothermal of frozen soil by coupling InSAR and ERT of the present invention;

[0062] Figure 3 is the monitoring result diagram of the regional surface deformation rate of the present invention;

[0063] Figure 4 is the result diagram of the water content inverted by ERT of the present invention;

[0064] Figure 5 is the inversion diagram of regional deformation - water content of the present invention;

[0065] Figure 6 (a) is the temperature distribution diagram of the present invention in the region at a depth of 0.3m;

[0066] Figure 6 (b) is the temperature distribution diagram of the present invention in the region at a depth of 1m;

[0067] Figure 6 (c) is the temperature distribution diagram of the present invention in the region at a depth of 2m;

[0068] Figure 6 (d) is the temperature distribution diagram of the present invention in the region at a depth of 2.5m;

[0069] Figure 7 (a) is the water content distribution diagram of the present invention in the region at a depth of 0.3m;

[0070] Figure 7 (b) is the water content distribution diagram of the present invention in the region at a depth of 1m;

[0071] Figure 7 (c) is the water content distribution diagram of the present invention in the region at a depth of 2m;

[0072] Figure 7 (d) is the water content distribution diagram of the present invention in the region at a depth of 2.5m;

[0073] Figure 8 This is the structural block diagram of the hydrothermal inversion system for frozen soil by coupling InSAR and ERT in the present invention. Specific implementation manners

[0074] To make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below in conjunction with embodiments and drawings. The illustrative implementation manners and descriptions of the present invention are only used to explain the present invention and shall not be construed as a limitation to the present invention.

[0075] The innovation point of the present invention is the coupling of InSAR and ERT, which realizes the high-precision dynamic inversion of the hydrothermal characteristics inside frozen soil at the regional scale, and provides key data and solutions for the safety of frozen soil engineering, early warning of frozen soil disasters, etc.

[0076] First, by utilizing the high spatio-temporal resolution observation characteristics of InSAR time-series interferometry technology, combined with the frost heave theory and the high vertical resolution of ERT, the blind area of existing remote sensing technology for detecting the hydrothermal inside frozen soil is broken through, and the high-precision inversion of the internal hydrothermal field in multiple dimensions of time (≤12 days), space (≤30 m), and vertical direction (≤3 m) is realized. Specifically, the combined use of InSAR and ERT: InSAR can provide surface deformation data with high spatio-temporal resolution, while ERT can provide relatively accurate resistivity data over a long period for determining the depth of the dynamic freeze-thaw interface and the internal soil water content. Through the combination of the frost heave theory, the vertical distribution characteristics of the soil water content retrieved by remote sensing can be well compensated, providing multi-dimensional information on the internal water content of frozen soil.

[0077] Second, by dynamically updating the depth of the freeze-thaw interface through ERT, the problems of low accuracy and poor universality of regional hydrothermal inversion caused by factors such as unclear internal structure, inaccurate freeze-thaw process parameters, and uncertain freeze-thaw interface in traditional hydrothermal coupling numerical models are solved, thus solving the problem of low traditional accuracy. Specifically, the combination of ERT and the hydrothermal coupling numerical model: introducing ERT to obtain the time-varying function of the freeze-thaw interface depth, correcting the traditional hydrothermal coupling theory, solving the uncertainty of hydrothermal inversion parameters, and improving the accuracy of internal water content inversion.

[0078] The present invention processes multiple SAR images to obtain InSAR surface deformation data, separates the freeze-thaw deformation data (i.e., the surface freeze-thaw deformation amount) d(t), and combines the frost heave theory to invert the water content corresponding to the seasonal deformation as the initial condition of the hydrothermal coupling inside frozen soil; uses time-series resistivity tomography ERT to obtain the resistivity vertical distribution function ρ(z,t), and obtains the dynamic value z f (t) of the freeze-thaw interface depth of frozen soil as the boundary condition for internal hydrothermal inversion, constructs a hydrothermal coupling numerical model considering the time-varying characteristics of frozen soil, and takes the first water content θ measured during the ERT observation period ERT(z, t) is the verification data for the parameters of the hydrothermal coupling numerical model. Through the optimized model, by inputting the surface layer temperature and InSAR surface deformation data, the temperature T(z, t) and the third water content θ(z, t) of the hydrothermal distribution inside the regional frozen soil are inverted.

[0079] Example 1

[0080] As Figure 1 and Figure 2 shown Figure 1 is the flow chart of the method for inverting the hydrothermal inside the frozen soil by coupling InSAR and ERT of the present invention; Figure 2 is the detailed flow chart of the method for inverting the hydrothermal inside the frozen soil by coupling InSAR and ERT of the present invention; The method for inverting the hydrothermal inside the frozen soil by coupling InSAR and ERT of the present invention includes:

[0081] S1, obtaining InSAR surface deformation data, ERT data, and surface layer temperature of the same research area;

[0082] Specifically, the InSAR surface deformation data is the surface deformation data obtained by processing SAR image data and DEM data based on the InSAR time series interferometry technology;

[0083] The ERT data is the resistivity profile data ρ(z, t) covering multiple periods throughout the year obtained by deploying a resistivity tomography system at typical points in the research area;

[0084] The surface layer temperature is the time series data Ts(t) of the surface temperature obtained through remote sensing product data, meteorological data, or ground monitoring stations.

[0085] In this embodiment, Sentinel-1 SAR image data and SRTM 30 m resolution DEM data covering the research area from 2019 to 2020 are obtained. The time span needs to cover the entire freeze-thaw cycle, and image registration and topographic phase removal are performed to eliminate the error caused by uneven terrain. Precise registration of the master and slave images is completed through the weighted least squares method, and the registration accuracy is evaluated using the standard deviation of the offset of the homologous point offset, which is generally less than 0.001 pixel.

[0086] A resistivity tomography system is deployed on the same profile in the research area to obtain resistivity profile data covering multiple periods throughout the year (it is necessary to ensure that the measured data covers the key periods in the freeze-thaw cycle, the initial freezing period, the deepest freezing interface period, the initial melting period, and the maximum melting period), that is, the ERT data ρ(z, t);

[0087] The surface layer soil temperature Ts(t) is obtained by the soil temperature and humidity sensors deployed in the research area.

[0088] S2. Based on the InSAR surface deformation data and an improved InSAR frozen soil deformation separation model, separate the freeze-thaw deformation data;

[0089] Specifically, step S2 includes:

[0090] S21. According to multi-temporal SAR image data and DEM data, after image registration and topographic phase removal, select an interference pair composed of preset time baseline and spatial baseline thresholds to generate a differential interference image;

[0091] S22. According to the differential interference image, select high-coherence points as the target point set;

[0092] S23. Use the minimum cost flow algorithm to perform phase unwrapping on the target point set to obtain the unwrapped phase;

[0093] S24. According to the unwrapped phase, combine the internal hydrothermal factors to simulate the seasonal phase and improve the traditional frozen soil deformation separation model to obtain an improved InSAR frozen soil deformation separation model;

[0094] S25. Based on the improved InSAR frozen soil deformation separation model, separate the seasonal deformation of the frozen soil and the deformation of the slope itself; and convert the seasonal deformation of the frozen soil into freeze-thaw deformation data. The freeze-thaw deformation data is the surface freeze-thaw deformation amount d(t) that changes with time in the study area, that is, the long-time series surface freeze-thaw deformation data set in the study area; where d(t)=[d(t 1 ), d(t 2 ), …, d(t N )], t 1 ~t N is the time series, and d(t N ) is the surface freeze-thaw deformation at the t N th time.

[0095] In this embodiment, in step S21, an appropriate time baseline (≤24 days) and spatial baseline (≤200 m) threshold are selected to form an interference pair to generate a differential interference map; in step S22, by analyzing features such as amplitude stability, phase stability, and temporal correlation, an optimal threshold is set, and stable high-coherence points are selected as the target point set for subsequent phase extraction. In step S23, the minimum cost flow algorithm is used to perform unwrapping on the target point set to restore the true phase, and the formula is:

[0096] ;

[0097] where Δ φ is the unwrapped deformation phase, φ top is the phase related to elevation, φ atm is the phase related to atmospheric disturbance, φ seasonal is the seasonal deformation of the frozen soil, φdef Self-deformation under gravity drive, φ noise is other noise.

[0098] The separation of frozen soil deformation in step S24 includes:

[0099] ① Seasonal deformation separation and noise removal;

[0100] Considering the seasonal settlement and uplift of frozen soil, the frozen soil deformation model can be expressed as: ;

[0101] where λ is the radar wavelength, Δh is the height difference, φ atm-h and φ atm-tur : the atmospheric component of height and turbulent path delay caused by complex mountain effects, which is related to the change of the propagation path when the radar signal passes through the atmosphere and is removed by spatial and temporal filtering; φ top (Δh, t) is the elevation phase related to Δh at time t; φ atm-h (Δh, t) is the atmospheric phase related to Δh at time t; Δφ atm-tur (t) is the atmospheric component of turbulent path delay at time t; φ noise Removed by spatial filtering, φ def (t) is the slope deformation component at a specific time t, A is the phase constant, N is the number of internal environment factors; βn is the coefficient related to the internal environment factor (n = 1, 2, 3,...), Bn(t) is the phase within the time interval t in different interferograms; through measured data, it is determined that the parameters affecting frozen soil deformation are soil temperature and humidity, so N = 2, and the phase after noise removal can be expressed as:

[0102] ;

[0103] where, β 1 is the coefficient related to internal environment factor 1, B1(t) is the phase of internal environment factor 1 at time t; β 2 is the coefficient related to internal environment factor 2, B2(t) is the phase of internal environment factor 2 at time t;

[0104] ② Phase conversion and solution;

[0105] Substitute the unknown parameter matrix into the regression equation and solve it from the following matrix equations:

[0106] ;

[0107] ;

[0108] Among them, B is a matrix that defines the small baseline interferogram combination. The unique solution can be obtained by using the singular value decomposition (SVD) method. The inversion accuracy is evaluated by the phase residual to obtain φ. seasonal The seasonal deformation of frozen soil and φ def is the self-deformation driven by gravity. Finally, φ seasonal phase information is converted into the freeze-thaw deformation amount, and the surface freeze-thaw deformation amount d(t) that changes with time in the target area can be obtained. Figure 3 The figure shows the monitoring results of the surface deformation rate in the area from February to April 2020. The unit of the deformation rate is mm / d (millimeters per day), and "-" indicates the settlement amount.

[0109] S3. According to multi-period ERT data, determine the depth of the dynamic freeze-thaw interface and invert the first water content.

[0110] Specifically, according to multi-period ERT data, determining the depth of the dynamic freeze-thaw interface and inverting the first water content includes:

[0111] According to multi-period ERT data, by identifying the resistivity gradient extreme values, determine the depth of the dynamic freeze-thaw interface at different times.

[0112] And combined with Archie's law, invert the first water content θ ERT (z,t) at different depths and different times.

[0113] In this embodiment, step S3 is specifically as follows:

[0114] (1) Dynamic calibration of the depth of the dynamic freeze-thaw interface:

[0115] For each period of ERT data ρ(z,t k ), calculate the depth z according to the following, time t k downward vertical resistivity gradient , z is the depth, and Δz is the depth step.

[0116] ;

[0117] Among them, ρ(z,t k ) is the resistivity at depth z at time t k , ρ(z+Δz,t k ) is the resistivity at depth z+Δz at time t k ; ρ(z-Δz,t k ) is the resistivity at depth z-Δz at time t k ;

[0118] Time t kDepth of the lower freeze-thaw interface is the point of maximum resistivity gradient, and the dynamic value of the freeze-thaw interface depth varying with time is obtained according to the following formula.

[0119] ;

[0120] When ( T f is the soil freezing temperature), then the corrected freeze-thaw interface depth is:

[0121] ;

[0122] Where Z f is the freeze-thaw interface depth; is the time t k depth z f temperature at; is the temperature gradient; T is the temperature; z is the depth; is the change in temperature with depth; is the corrected freeze-thaw interface depth.

[0123] (2) Water content inversion based on Archie model

[0124] Based on Archie's law and combined with the dynamic freeze-thaw interface depth, a resistivity-water content relationship within frozen soil is constructed:

[0125] In the unfrozen region, the soil is mainly composed of liquid water and has good conductivity. The resistivity-water content relationship is:

[0126] ;

[0127] Where ϕ(z) is the soil porosity at depth z, z is the depth, ρ w is the resistivity of water, is the empirical coefficient reflecting soil particle characteristics, m is the pore index, s is the soil saturation, is the resistivity at depth z at time t.

[0128] In the frozen region, the relationship between resistivity and unfrozen water content is:

[0129] ;

[0130] Where ρ i is the ice resistivity, is the residual water content after soil freezing, and η is the unfrozen water influence coefficient. is the resistivity at depth z at time t.

[0131] (3)Parameter optimization

[0132] Layering soil temperature and humidity sensors along the profile in the typical monitoring area to obtain the measured data T sensor (z, t) and θ sensor (z, t) of soil temperature and humidity, and determine the resistivity-water content inversion parameters m, s through the least squares method, etc., and generate the resistivity inversion water content distribution map of the typical area.

[0133] As Figure 4 shown is the ERT inversion water content result map of a certain period during the thawing of frozen soil. In this embodiment, , m = 1, ϕ = 0.5 - 0.6, s = 3. Obtain the current freeze-thaw interface position z f As Figure 4 shown by the red dashed line in.

[0134] S4, according to the dynamic freeze-thaw interface depth, correct the active layer thickness in the frost heaving theoretical model to obtain the corrected frost heaving theoretical model;

[0135] In this embodiment, considering that the InSAR seasonal deformation is affected by the internal structure of frozen soil and the freeze-thaw interface depth, introduce the dynamic freeze-thaw interface depth observed by ERT to correct the active layer thickness in the frost heaving theoretical model.

[0136] S5, according to the freeze-thaw deformation data, perform water content inversion based on the corrected frost heaving theoretical model to obtain the second water content;

[0137] Specifically, the second water content is the average soil water content corresponding to the freezing depth or thawing depth.

[0138] In this embodiment, the surface freeze-thaw deformation amount d(t) is mainly caused by the freezing expansion or thawing settlement of water in the soil, and the deformation is related to the freezing or thawing thickness and water content. Considering that the heat conduction processes during the frost heaving and thawing processes of frozen soil are inconsistent, the surface temperature during the freezing period is unidirectional conduction from top to bottom, and during the thawing period, it is two-way conduction from top and bottom to the middle. Introduce the freeze-thaw interface depth z f (t) observed in-situ by ERT, and it can be obtained that during the freezing / thawing period, at time t the change in freezing or thawing thickness Δz f (t) of the average soil water content , and the relational expression is:

[0139] The water content inversion formula during the freezing period is:

[0140] ;

[0141] The water content inversion formula for the melting period is as follows:

[0142] ;

[0143] In the formula, is the time t The change in the freezing depth or melting depth of the average soil water content, that is, the second water content; is the residual water content after soil freezing; is the temperature correction coefficient; is the soil freezing temperature; is the surface temperature; is the change in freeze-thaw deformation, equal to d ( t ) - d( t - 1), d ( t ) is the surface freeze-thaw deformation at time t , d ( t - 1) is the freeze-thaw deformation at the previous time t of time t - 1; is the ice resistivity; is the resistivity of water; is the change in the freezing depth or melting depth of the soil layer, equal to z f ( t ) - z f ( t - 1), z f ( t ) is the freeze-thaw interface depth at time t , z f ( t - 1) is the freeze-thaw interface depth at the previous time t of time t - 1, provided by the ERT data; is the water-ice density ratio, with a value of 1.09.

[0144] According to the regional InSAR surface deformation data obtained previously, the regional deformation amount - water content is obtained as Figure 5 shown. The water content here is the second water content.

[0145] S6. Using the second water content and the surface layer temperature as the initial conditions, and the dynamic freeze-thaw interface depth as the boundary conditions, input into the hydrothermal coupling numerical model; and optimize the hydrothermal coupling numerical model using the first water content to obtain the optimized dynamic hydrothermal coupling numerical model; the initial water content is the water content observed from the InSAR surface deformation data;

[0146] Specifically, the temperature gradient change of frozen soil follows heat conduction, and the moisture migration follows Darcy's law. By combining ERT to obtain the dynamic freeze-thaw interface, a hydrothermal coupling numerical model considering time-varying characteristics is constructed. Using the initial water content and surface layer temperature obtained from InSAR deformation observations as the initial conditions, and the depth of the dynamic freeze-thaw interface observed by ERT as the boundary conditions, input into the hydrothermal coupling numerical model, and using the first water content θ ERT (z,t) inverted by ERT to optimize the model parameters (thermal conductivity k, moisture diffusion coefficient D, expansion coefficient α, etc.), and the optimal hydrothermal numerical solution is obtained through iterative calculation. Specifically:

[0147] 1) Establish a hydrothermal coupling numerical model considering time-varying characteristics

[0148] The water content in frozen soil is not statically distributed. Water undergoes phase changes (liquid water ↔ ice) during the freezing / melting process, accompanied by the release / absorption of latent heat. Moreover, during the freezing / melting process, the depth position of the freeze-thaw interface z f is not fixed. Based on the above, an improved hydrothermal coupling numerical model is established as follows:

[0149] Heat conduction equation: ;

[0150] Moisture migration equation: ;

[0151] The latent heat flux of phase change is: ;

[0152] In the formula, is the effective thermal conductivity of the soil, is the volumetric heat capacity of the soil, is the latent heat of water-ice phase change, θ i is the ice volume content, K w is the moisture permeability, which is related to soil type and temperature; and are the thermal conductivities of the frozen zone and the unfrozen zone respectively, is the temperature of the unfrozen zone, is the water source term, including factors such as precipitation and evaporation; is to solve the partial differential, is the gradient solution; is the temperature at depth z at time t; is the ice volume content at depth z at time t; is the water content at depth z at time t; z f (t) is the depth of the freeze-thaw interface at time t.

[0153] 2) Determine the boundary conditions of the equation

[0154] The upper boundary is the surface layer (z = 0), satisfying ; is the temperature of the surface layer at time t; is the temperature of the surface layer obtained by monitoring with sensors / weather stations, etc.;

[0155] Time t The melting depth Δ at z f ( t ) The corresponding average water content , obtained from step S5;

[0156] The depth position z of the freeze-thaw interface f At: ; is the temperature at the depth position z of the freeze-thaw interface at time t f ; is the soil freezing temperature.

[0157] 3) Optimization of the parameters of the hydrothermal coupling process

[0158] Specifically, the optimization steps of the hydrothermal coupling numerical model are as follows:

[0159] Based on the prior parameters (soil effective thermal conductivity , soil volumetric heat capacity , water permeability , residual unfrozen water content , expansion coefficient ), generate N parameter sets ; is the parameter vector to be optimized, i is a positive integer, T is the transpose of the vector;

[0160] For each parameter set , obtain the predicted water content value z at depth t at time through the hydrothermal coupling numerical model;

[0161] Generate the first water content obtained by inversion as the observation vector t θ , θ ERT ( z m , t ) is the depth z at time t obtained by ERT inversion mThe first water content at

[0162] Construct an objective function based on the predicted water content value and the observation vector Calculate the global residual sum of squares N t is the observation time N z is the maximum freezing or thawing thickness θ ERT ( z , t ) is the time inversed by ERT t The first water content at depth z

[0163] According to the global residual sum of squares, obtain the minimum residual by gradient calculation to update the hydrothermal characteristic parameters, and obtain an optimized dynamic hydrothermal coupling numerical model.

[0164] The above technical solution corrects the predicted water content value based on the observation vector and optimizes the hydrothermal coupling numerical model.

[0165] S7. Input the InSAR surface deformation data and the surface layer temperature into the optimized dynamic hydrothermal coupling numerical model, and output the temperature T(z,t) and the third water content θ(z,t) at different depths and different times inside the permafrost in the study area.

[0166] In this embodiment, using the optimized dynamic hydrothermal coupling numerical model, on the premise that the internal structure of the regional permafrost is consistent, input the InSAR surface deformation data and the surface layer temperature to obtain the temperature T(z,t) and the third water content θ(z,t) at different depths and different times inside the regional permafrost. As shown in FIGS. 6(a)-6(d) and FIGS. 7(a)-7(d), FIGS. 6(a)-6(d) and FIGS. 7(a)-7(d) are the hydrothermal content maps of the region obtained according to the foregoing steps, reflecting that this method can track the freezing and thawing processes of the permafrost in real time and dynamically invert the temperature and water content parameters at different depths inside the upper permafrost in the region, providing an effective solution for inverting the internal hydrothermal field at the regional level.

[0167] Embodiment 2

[0168] As Figure 8 shown, the difference between this embodiment and Embodiment 1 is that this embodiment provides a permafrost internal hydrothermal inversion system coupling InSAR and ERT, and this system corresponds one-to-one to the permafrost internal hydrothermal inversion method coupling InSAR and ERT in Embodiment 1; this system includes:

[0169] An acquisition unit for acquiring InSAR surface deformation data, ERT data, and surface layer temperature of the same study area;

[0170] A separation unit, configured to separate freeze-thaw deformation data based on an improved InSAR frozen soil deformation separation model according to InSAR surface deformation data;

[0171] A determination unit, configured to determine the depth of the dynamic freeze-thaw interface according to multi-period ERT data and invert the first water content;

[0172] A frost heave model correction unit, configured to correct the thickness of the active layer in the frost heave theoretical model according to the depth of the dynamic freeze-thaw interface to obtain a corrected frost heave theoretical model;

[0173] A water content inversion unit, configured to invert the water content based on the corrected frost heave theoretical model according to the freeze-thaw deformation data to obtain the second water content;

[0174] A hydrothermal model optimization unit, configured to input a hydrothermal coupling numerical model with the second water content and the surface layer temperature as initial conditions and the depth of the dynamic freeze-thaw interface as boundary conditions; and optimize the hydrothermal coupling numerical model by using the first water content to obtain an optimized dynamic hydrothermal coupling numerical model;

[0175] A frozen soil internal hydrothermal inversion unit, configured to input the InSAR surface deformation data and the surface layer temperature into the optimized dynamic hydrothermal coupling numerical model and output the temperature T(z,t) and the third water content θ(z,t) at different depths and different times inside the frozen soil in the study area.

[0176] In this embodiment, the determination unit includes:

[0177] A determination subunit, configured to determine the depth of the dynamic freeze-thaw interface at different times by identifying the extreme values of the resistivity gradient according to multi-period ERT data;

[0178] An inversion subunit, configured to invert and obtain the first water content θ ERT (z,t) at different depths and different times in combination with Archie's law.

[0179] Wherein, the execution processes of each unit are carried out according to the process steps of the InSAR coupled ERT frozen soil internal hydrothermal inversion method in Embodiment 1, and will not be elaborated one by one in this embodiment.

[0180] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0181] This application is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram, and the combination of flows and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0182] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, such that the instructions stored in the computer-readable memory generate a manufactured article including instruction means that implement the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0183] These computer program instructions can also be loaded onto a computer or other programmable data processing device, such that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, so that the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one flow Figure 1 one flow or multiple flows and / or blocks Figure 1 or multiple blocks.

[0184] The specific embodiments described above further elaborate on the objectives, technical solutions, and beneficial effects of the present invention. It should be understood that the above description is only for the specific embodiments of the present invention and is not used to limit the protection scope of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included in the protection scope of the present invention.

Claims

1. The method of hydrothermal inversion of frozen soil internal by InSAR coupled with ERT is characterized by: The method includes: Obtain InSAR surface deformation data, ERT data and surface layer temperature in the same study area; According to the InSAR surface deformation data, freeze-thaw deformation data are separated based on the improved InSAR frozen soil deformation separation model; Based on multiple ERT data, the dynamic freeze-thaw interface depth is determined and the first water content is inverted; According to the depth of the dynamic freeze-thaw interface, the thickness of the active layer in the frost heave theoretical model is corrected to obtain the corrected frost heave theoretical model; According to the freeze-thaw deformation data, the water content is inverted based on the modified frost heave theory model to obtain the second water content; The second water content and the surface temperature are used as initial conditions, and the dynamic freeze-thaw interface depth is used as the boundary condition to input the hydrothermal coupling numerical model; and the hydrothermal coupling numerical model is optimized using the first water content to obtain an optimized dynamic hydrothermal coupling numerical model; The InSAR surface deformation data and surface layer temperature are input into the optimized dynamic water-heat coupling numerical model, and the temperature and third water content of the frozen soil at different depths in the study area at different times are output. The optimization steps of the hydrothermal coupling numerical model are: Based on the prior parameters, generate N Parameter Set ; is the parameter vector to be optimized, i is a positive integer, is the soil effective thermal conductivity parameter, is the volume heat capacity of soil, is the water permeability, is the residual unfrozen water content, is the expansion coefficient; [ ] T is the transpose of a vector; For each parameter set , the depth is obtained by the hydrothermal coupling numerical model z time t Predicted moisture content ; The first moisture content of the inversion is calculated by time t Generates the observation vector , θ ERT ( z m , t ) is the ERT reverse performance time t Lower depth z m The first moisture content at According to the water content prediction value and the observation vector, the objective function is constructed. Calculate the global residual sum of squares, N t is the observation time, N z is the maximum freezing or thawing thickness; θ ERT ( z , t ) is the ERT reverse performance time t The first water content at the lower depth z; According to the global residual square sum, the minimum residual is obtained through gradient calculation to update the hydrothermal characteristic parameters and obtain an optimized dynamic hydrothermal coupling numerical model.

2. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 1 is characterized in that: The InSAR surface deformation data is surface deformation data obtained by processing SAR image data and DEM data based on InSAR time series interferometry technology; The ERT data is obtained by deploying a resistivity tomography system in the study area to obtain multi-period resistivity profile data covering the entire year; The surface layer temperature is time series data of the surface temperature obtained through remote sensing product data, meteorological data or ground monitoring stations.

3. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 1 is characterized in that: Based on the improved InSAR frozen soil deformation separation model, freeze-thaw deformation data are separated, including: According to the multi-period SAR image data and DEM data, after image registration and terrain phase removal, the preset time baseline and space baseline threshold are selected to form an interference pair to generate a differential interferometric image; According to the differential interference image, selecting high coherence points as a target point set; Using a minimum cost flow algorithm, phase unwrapping is performed on the target point set to obtain an unwrapped phase; According to the unwrapped phase, the traditional frozen soil deformation separation model is improved by combining the internal hydrothermal factor to simulate the seasonal phase, and the improved InSAR frozen soil deformation separation model is obtained. Based on the improved InSAR frozen soil deformation separation model, the seasonal deformation of frozen soil and the deformation of the slope itself are separated; and the seasonal deformation of frozen soil is converted into freeze-thaw deformation data.

4. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 3 is characterized in that: The freeze-thaw deformation data is the surface freeze-thaw deformation d(t) that changes with time in the study area, where d(t)=[d(t1),d(t2),…,d(t N )],t1~t N is a time series, d(t N ) is the tth N Surface freeze-thaw deformation over time.

5. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 1 is characterized in that: Based on multiple ERT data, the dynamic freeze-thaw interface depth is determined and the first water content is inverted, including: Based on multi-period ERT data, the depth of the dynamic freeze-thaw interface at different times is determined by identifying the extreme values ​​of resistivity gradient; Combined with Archie's law, the first water content at different depths at different times is inverted.

6. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 1 is characterized in that: The second moisture content is the average soil moisture content corresponding to the freezing depth or the thawing depth.

7. The method for hydrothermal inversion of frozen soil interior by InSAR coupled with ERT according to claim 6 is characterized in that: Based on the modified frost heave theory model, the water content is inverted to obtain the second water content, including: The inversion formula of water content during freezing period is: ； The inversion formula of water content during the melting period is: ； In the formula, For time t Changes in freezing depth or thawing depth The average soil moisture content, i.e. the second moisture content; It is the residual water content of soil after freezing; is the temperature correction factor; is the soil freezing temperature; is the surface temperature; is the freeze-thaw deformation change, equal to d ( t )-d( t -1), d ( t ) is the time t The surface freeze-thaw deformation under d ( t -1) is time t The time before t -1 freeze-thaw deformation; is the ice resistivity; is the resistivity of water; is the change in soil freezing depth or melting depth, equal to z f ( t )-z f ( t -1), z f ( t ) is the time t Depth of freeze-thaw interface below, z f ( t -1) is time t The time before t The freeze-thaw interface depth of -1 is provided by ERT data.

8. The InSAR coupled ERT frozen soil internal water and heat inversion system is characterized by: The system includes: The acquisition unit is used to obtain InSAR surface deformation data, ERT data and surface layer temperature in the same study area; A separation unit is used to separate freeze-thaw deformation data based on the InSAR surface deformation data and the improved InSAR frozen soil deformation separation model; A determination unit is used to determine the depth of the dynamic freeze-thaw interface and invert the first water content based on multiple ERT data; A frost heave model correction unit is used to correct the thickness of the active layer in the frost heave theoretical model according to the depth of the dynamic freeze-thaw interface to obtain a corrected frost heave theoretical model; A water content inversion unit is used to perform water content inversion based on freeze-thaw deformation data and a modified frost heave theory model to obtain a second water content; The hydrothermal model optimization unit is used to input the hydrothermal coupling numerical model with the second water content and the surface temperature as initial conditions and the dynamic freeze-thaw interface depth as boundary conditions; and optimize the hydrothermal coupling numerical model with the first water content to obtain an optimized dynamic hydrothermal coupling numerical model; The hydrothermal inversion unit inside the frozen soil is used to input the InSAR surface deformation data and the surface layer temperature into the optimized dynamic hydrothermal coupling numerical model, and output the temperature and the third water content at different depths in the frozen soil of the study area at different times; The optimization steps of the hydrothermal coupling numerical model are: Based on the prior parameters, generate N Parameter Set ; is the parameter vector to be optimized, i is a positive integer, is the soil effective thermal conductivity parameter, is the volume heat capacity of soil, is the water permeability, is the residual unfrozen water content, is the expansion coefficient; [ ] T is the transpose of a vector; For each parameter set , the depth is obtained by the hydrothermal coupling numerical model z time t Predicted moisture content ; The first moisture content of the inversion is calculated by time t Generates the observation vector , θ ERT ( z m , t ) is the ERT reverse performance time t Lower depth z m The first moisture content at According to the water content prediction value and the observation vector, the objective function is constructed. Calculate the global residual sum of squares, N t is the observation time, N z is the maximum freezing or thawing thickness; θ ERT ( z , t ) is the ERT reverse performance time t The first water content at the lower depth z; According to the global residual square sum, the minimum residual is obtained through gradient calculation to update the hydrothermal characteristic parameters and obtain an optimized dynamic hydrothermal coupling numerical model.

9. The InSAR-coupled ERT frozen soil internal hydrothermal inversion system according to claim 8, characterized in that: The determining unit comprises: Determine subunits for determining the depth of the dynamic freeze-thaw interface at different times by identifying the extreme values ​​of resistivity gradient based on multi-period ERT data; The inversion subunit is used to combine Archie's law to invert the first water content at different depths at different times.

Images