Frozen soil active layer thickness inversion method, medium and equipment

Abstract

The invention discloses a frozen soil active layer thickness inversion method, medium and equipment, and relates to the field of radar remote sensing image processing, and the frozen soil active layer thickness inversion method mainly comprises the steps: obtaining deformation information according to single-view complex orbit descending data, digital elevation model data and atmospheric delay data; according to the deformation information, obtaining to-be-selected long-term deformation and to-be-selected seasonal deformation of the permafrost by using multiple decomposition algorithms; evaluating the to-be-selected long-term deformation and the to-be-selected seasonal deformation, and determining long-term deformation and seasonal deformation; and obtaining the thickness of the active layer of the research area according to the seasonal deformation and the layered soil water content data. By implementing the frozen soil active layer thickness inversion method, the medium and the equipment provided by the invention, the accuracy of frozen soil active layer thickness inversion can be improved.

Description

Technical Field

[0001] The invention relates to the field of radar remote sensing image processing, and more specifically to a method, medium and equipment for inverting the thickness of a frozen soil active layer. Background Art

[0002] As one of the sources of the Yellow River, the permafrost in the source area of ​​the Yellow River plays a vital role in the storage and release of water resources. The melting of the permafrost layer will affect the infiltration of precipitation and the formation of surface runoff, directly affecting the water volume and quality of the Yellow River. Monitoring the deformation of permafrost can understand the dynamic trend of water resources, thereby providing data support for water resource management and regulation. At the same time, through the deformation of permafrost, the impact of climate change on the region can be more accurately assessed, and a scientific basis can be provided for responding to climate change. Therefore, it is crucial to monitor the deformation of permafrost over a large area and long time series in the source basin of the Yellow River. Traditional permafrost deformation monitoring mainly uses technical means such as field geophysical survey, leveling and GNSS (Global Navigation Satellite System) site measurement. These methods can obtain high-precision measurement data at a point scale, including permafrost deformation and ALT. However, permafrost areas are usually remote, with harsh environments and inconvenient transportation. The installation and maintenance of traditional monitoring facilities require a lot of manpower and material resources, and large-scale monitoring research is impossible. The development of remote sensing technology provides a new research method for permafrost environment research. Among them, Synthetic Aperture Radar Interferometry (InSAR) is a remote sensing technology based on radar wave phase information. It has the advantages of all-weather monitoring, large-scale coverage, and high precision. It provides a more effective way to monitor slow surface deformation. InSAR technology is also widely used in monitoring applications such as deformation and infrastructure in permafrost areas. In permafrost areas, the surface environment usually changes dramatically due to climate change, resulting in InSAR technology usually facing the impact of decoherence and atmospheric delay.

[0003] In order to overcome these limitations, scholars have proposed time series InSAR technology, such as Permanent Scatterer InSAR (PS-InSAR), Small Baseline Subset InSAR (SBAS-InSAR) and Distributed Scatterer InSAR (DS-InSAR). Time series InSAR technology has gradually become an important means of monitoring and research related to permafrost areas on the Qinghai-Tibet Plateau, and has made many advances. However, when the decomposition algorithm of the existing technology separates long-term deformation and seasonal deformation and simulates with a combination of linear and periodic functions, the periodic function does not match the seasonal deformation well, and it is impossible to deeply analyze the changes in the frozen ground surface, accurately calculate the relevant parameters, and accurately invert the ALT. Summary of the invention

[0004] The purpose of the present invention is to provide a method, medium and equipment for inverting the thickness of a frozen soil active layer, which can improve the inversion accuracy of the thickness of a frozen soil active layer.

[0005] The present invention provides a method for inverting the thickness of a permafrost active layer, comprising the following steps: S1: acquiring single-view complex descending orbit data, digital elevation model data and atmospheric delay data covering a study area; S2: acquiring deformation information according to the single-view complex descending orbit data, digital elevation model data and atmospheric delay data; S3: acquiring candidate long-term deformation and candidate seasonal deformation of permafrost according to the deformation information by using a variety of decomposition algorithms; S4: evaluating the candidate long-term deformation and the candidate seasonal deformation to determine the long-term deformation and the seasonal deformation; S5: acquiring stratified soil moisture content data, and acquiring the thickness of the active layer in the study area according to the seasonal deformation and the stratified soil moisture content data.

[0006] Furthermore, the single-look complex descending orbit data of the above-mentioned permafrost active layer thickness inversion method are single-look complex descending orbit data of the Sentinel-1A satellite, the digital elevation model data are digital elevation model data of the Space Shuttle Radar Topography Mission, and the atmospheric delay data are atmospheric delay data of the General Atmospheric Correction Online Service.

[0007] Furthermore, step S2 of the above-mentioned method for inverting the thickness of the permafrost active layer specifically includes: S21: according to the single-view complex descending orbit data, using the small baseline subset interferometry measurement method, according to the preset time baseline threshold and the preset space baseline threshold, obtaining an interference pair; S22: using the digital elevation model data to perform terrain phase elimination on the interference pair to obtain an interference fringe pattern; S23: using the atmospheric delay data to process and interpret the interference fringe pattern to obtain a new interference pattern; S24: using the minimum cost flow algorithm to perform phase unwrapping on the new interference pattern to obtain an unwrapped pattern; according to the interference pair, the interference fringe pattern, the new interference pattern and the unwrapped pattern, deformation information is obtained.

[0008] Furthermore, the preset time baseline threshold of the above-mentioned permafrost active layer thickness inversion method is 60 days, and the preset space baseline threshold is 150 meters.

[0009] Furthermore, the various decomposition algorithms of the above-mentioned frozen soil active layer thickness inversion method include sine-cosine decomposition algorithm, X11 decomposition algorithm, STL decomposition algorithm, decomposition algorithm based on independent component analysis, BEAST decomposition algorithm and RobustSTL decomposition algorithm.

[0010] Furthermore, step S4 of the above-mentioned permafrost active layer thickness inversion method specifically includes: S41: obtaining a simulated multi-year average deformation curve, and obtaining morphological similarity based on the simulated multi-year average deformation curve and the selected seasonal deformation using a dynamic time warping algorithm; S42: obtaining a root mean square error based on the selected seasonal deformation and the deformation information; S43: obtaining an optimal decomposition method based on the multiple decomposition algorithms using the morphological similarity and the root mean square error; and determining long-term deformation and seasonal deformation using the optimal decomposition method.

[0011] Furthermore, step S41 of the above-mentioned frozen soil active layer thickness inversion method specifically includes: obtaining a simulated multi-year average deformation curve, and using a dynamic time warping algorithm to obtain a morphological similarity based on the simulated multi-year average deformation curve and the selected seasonal deformation, such as the formula: , in, Represents the position in the cumulative distance matrix The value at Represents the Euclidean distance between two data points.

[0012] Furthermore, step S5 of the above-mentioned frozen soil active layer thickness inversion method specifically includes: obtaining layered soil moisture data, and obtaining the active layer thickness of the study area according to the seasonal deformation and the layered soil moisture data, such as the formula: , , in, is the thickness of the active layer, is the soil moisture layer number, For the The melting depth of the active layer within the soil layer, For the Seasonal sedimentation in the soil layer, and are the densities of pure water and ice, respectively, For the i The average soil moisture content of the layer, is a function of soil moisture content, is the soil depth, a, b, c, d and e are the coefficients of the polynomial fit function.

[0013] The present invention also provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the steps of the above-mentioned frozen soil active layer thickness inversion method are implemented.

[0014] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps of the above-mentioned method for inverting the thickness of the frozen soil active layer are implemented.

[0015] The implementation of the frozen soil active layer thickness inversion method, medium and equipment provided by the present invention has the following beneficial effects: The present invention uses SBAS-InSAR technology to obtain the surface deformation data of the study area without long-term series, uses a variety of deformation decomposition algorithms for seasonal deformation and long-term deformation in the study area, calculates the applicability of six decomposition algorithms in decomposing long-term deformation and seasonal deformation respectively, and uses two evaluation methods to evaluate all decomposition algorithms, selects the best decomposition algorithm to decompose the deformation of the long-term series into long-term deformation and seasonal deformation, deeply analyzes the relationship between different types of deformation and environmental factors, reveals the influence of different environmental factors on long-term deformation and seasonal deformation, based on deformation data and a variety of deformation decomposition algorithms, performs spatial simulation of deformation rate of frozen soil area in the Yellow River Basin, uses layered soil moisture to establish a piecewise constant ALT inversion model, combines deformation information and multi-layer soil moisture to invert the thickness of the active layer, analyzes the changes of frozen soil surface in more depth, and is more accurate in calculating related parameters, and obtains more accurate ALT inversion results, which provides a more accurate method for ALT extraction in the source area of ​​the Yellow River, and provides a reference for in-depth understanding of the hydrological, ecological and surface stability characteristics of the source area of ​​the Yellow River. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] The present invention will be further described below with reference to the accompanying drawings and embodiments, in which: Figure 1 It is a flow chart of the inversion method of the thickness of the frozen soil active layer provided by the present invention; Figure 2 It is the overall algorithm flow chart of the frozen soil active layer thickness inversion method provided by the present invention; Figure 3 It is a basic information schematic diagram of the research area provided by the present invention; Among them, (a) is a schematic diagram of the location of the Yellow River source, (b) is a schematic diagram of the coverage of SAR images, and (c) is a schematic diagram of the elevation of the Yellow River source; Figure 4 It is a schematic diagram of the inversion result of ALT provided by the present invention; Figure 5 It is a structural block diagram of the computer device provided by the present invention. DETAILED DESCRIPTION

[0017] In order to have a clearer understanding of the technical features, purposes and effects of the present invention, specific embodiments of the present invention are now described in detail with reference to the accompanying drawings.

[0018] Figure 1 The schematic diagram of the inversion method of the thickness of the frozen soil active layer in this embodiment is shown. In this embodiment, the inversion method of the thickness of the frozen soil active layer includes the following steps: S1: Acquire single-view complex descending orbit data, digital elevation model data and atmospheric delay data covering the study area; Specifically, the single look complex down-orbit data is the single look complex (SLC) down-orbit data of the Sentinel-1A satellite, the digital elevation model data is the Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) data, and the atmospheric delay data is the Generic Atmospheric Correction Online Service (GACOS) atmospheric delay data; S2: Obtaining deformation information according to the single-view complex descending orbit data, digital elevation model data and atmospheric delay data; Specifically, step S2 includes: S21: obtaining an interference pair according to the single-view complex descending orbit data by using a small baseline subset interferometry method according to a preset time baseline threshold and a preset space baseline threshold; Specifically, the preset time baseline threshold is 60 days, and the preset space baseline threshold is 150 meters; S22: performing terrain phase elimination on the interference pair using the digital elevation model data to obtain an interference fringe pattern; Specifically, in step S22, the Space Shuttle Radar Topography Mission Digital Elevation Model data is used as an external DEM to eliminate the terrain phase of the interferogram, and a multi-look operation of 8 pixels in the range and 2 pixels in the azimuth is performed during the interferometric processing to suppress noise and improve the signal-to-noise ratio of the interferogram; S23: Processing and interpreting the interference fringe pattern using the atmospheric delay data to obtain a new interference pattern; Specifically, in step S23, after obtaining the interference fringe pattern, the atmospheric delay data of the General Atmospheric Correction Online Service GACOS is used to perform data processing and interpretation to eliminate or reduce the phase change caused by the atmospheric delay in the source area of ​​the Yellow River and reduce the influence of the atmospheric delay; S24: performing phase unwrapping on the new interference pattern using a minimum cost flow algorithm to obtain an unwrapping pattern; and obtaining deformation information according to the interference pair, the interference fringe pattern, the new interference pattern and the unwrapping pattern; Specifically, in step S24, the minimum cost flow (MCF) algorithm of the Delaunay triangulation network is used to perform phase unwrapping; after the interference processing is completed, the quality of the interference pair is judged according to the coherence coefficient diagram, the interference diagram and the unwrapping diagram, and the interference pairs with low coherence and poor unwrapping quality are removed to avoid affecting the inversion effect; S3: Based on the deformation information, a variety of decomposition algorithms are used to obtain the long-term deformation and seasonal deformation of permafrost; Specifically, the multiple decomposition algorithms include a sine-cosine decomposition algorithm, an X11 decomposition algorithm, an STL decomposition algorithm, a decomposition algorithm based on independent component analysis, a BEAST decomposition algorithm, and a RobustSTL decomposition algorithm; Specifically, in step S3, six decomposition algorithms are applied to decompose the deformation, which are as follows: (1) Sin-cosine decomposition algorithm:

[0019] In the formula, d is the cumulative deformation; v is the linear deformation rate; t To obtain the time accumulation of the image moment relative to the reference time; T is the seasonal deformation cycle, usually taken as 1 year; a 1 , a 2 , b are the parameters to be determined for the model; (2) X11 decomposition algorithm: First, the original data is smoothed by Henderson weighted moving average to obtain the trend term; then, the Musgrave asymmetric moving average is used to fill the endpoint values ​​that cannot be accurately estimated by Henderson weighted moving average to capture long-term, medium-term and short-term seasonal changes; next, the seasonal factor is estimated for the adjusted observations to obtain the adjusted trend term; finally, the trend factor is multiplied by the seasonal factor to obtain the final trend term. The whole process is iterated until convergence is reached; (3) STL decomposition algorithm: The STL decomposition algorithm consists of two processes: an inner loop and an outer loop. The inner loop divides the time series into periodic subsequences by subtracting the estimated value of the current trend from the original sequence; each periodic subsequence is smoothed by Locally Estimated Scatterplot Smoothing (LOESS) and then processed by a low-pass filter to obtain the seasonal component; finally, the seasonal component is subtracted from the original data to obtain the remaining sequence; the remaining sequence is smoothed again with LOESS to obtain the final trend; LOESS smoothing curve g(x) The calculation method is as follows: Assuming x i For the observation value, choose an integer q As the closest x The number of values, then when q ≤ n Time Based x i arrive x The distance is x i Set a neighborhood weight, the expression is:

[0020]

[0021] In the formula, W represents a weight function, The neighborhood weight of and x The distance between The longer, the weight The smaller; finally, at point Find a weight value The fitted value of g ( x ); The main function of the outer loop is to adjust the robust weights. If there are outliers in the sequence, the remainder will be larger, which can reduce the impact; (4) Decomposition algorithm based on independent component analysis: Independent Component Analysis (ICA) is a computational signal separation method that can effectively decompose mixed signals into a linear combination of a set of statistically independent signal sources. The algorithm first performs zero mean and whitening on the original deformation values ​​to generate an orthogonal matrix. Then, the mixed signal is linearly transformed and represented as an uncorrelated variable with a variance of 1 by whitening. The mixed matrix after zero mean and whitening preprocessing is represented by W, and the problem becomes:

[0022] In the formula, is the orthogonal adjustment matrix; then, the source matrix is ​​estimated by the following formula:

[0023] In the formula ， is the unmixing matrix W ; In order to obtain the unmixing matrix W, ICA must select an objective function to measure independence and find a suitable optimization algorithm; different objective functions directly affect the consistency and robustness of the ICA algorithm, and the selected optimization algorithm determines the convergence speed and stability of the algorithm; currently, the most widely used ICA algorithm is the fast fixed point algorithm (Fast ICA), which has a cubic linear convergence speed and is significantly better than other algorithms; therefore, using FastICA iterative calculation, by maximizing the spatial non-Gaussian source, the source matrix S, the mixing matrix W and the unmixing matrix U are derived; (5) BEAST decomposition algorithm: BEAST decomposition is performed through two sets of parameters, one for seasonality and the other for trend; seasonality and trend are modeled by piecewise harmonic models and piecewise linear functions, respectively, and the mutation points are implicitly encoded in the seasonal and trend signals; the seasonal signal S(t) is modeled as a linear combination of a set of harmonic terms, and the specific expression is as follows:

[0024] In the formula, P is the period of the seasonal signal; p is the number of segments of the seasonal signal, L k is the order of the harmonic model; and is the model coefficient, which represents the amplitude of the harmonic term; trend signal T(t) Modeling is done through piecewise linear function, the expression is:

[0025] In the formula, m is the number of trend change points, a j and b j are the model coefficients; (6) Robust STL decomposition algorithm: The algorithm first uses a bilateral filter to denoise the input time series data, and then robustly extracts the trend by solving the least absolute deviation (LAD) regression with sparse regularization; since the seasonal component is assumed to change slowly, a seasonal difference operation is performed on the denoised signal to reduce the seasonal effect; then a non-local seasonal filter is used to obtain the seasonal signal component; unlike the traditional method of estimating seasonality by considering only subsequences, this method uses K neighborhoods centered on various points; the weights assigned to these neighborhoods are based on temporal distance and seasonal similarity, making the algorithm more robust to outliers and seasonal changes; finally, an adjustment is made to ensure that the sum of all seasonal components within a cycle is zero; this adjustment is achieved by removing the mean, which corresponds to the estimation of the trend point; the algorithm iteratively improves the estimates of trend, seasonality, and residual components to improve accuracy; S4: Evaluate the candidate long-term deformation and the candidate seasonal deformation to determine the long-term deformation and the seasonal deformation; Specifically, step S4 includes: S41: obtaining a simulated multi-year average deformation curve, and obtaining morphological similarity using a dynamic time warping algorithm according to the simulated multi-year average deformation curve and the selected seasonal deformation; Specifically, step S41 includes: obtaining a simulated multi-year average deformation curve, and obtaining a morphological similarity using a dynamic time warping algorithm according to the simulated multi-year average deformation curve and the selected seasonal deformation, such as the formula:

[0026] in, Represents the position in the cumulative distance matrix The value at Represents the Euclidean distance between two data points; Specifically, in step S41, the dynamic time warping (DTW) algorithm is used to determine the degree of seasonal deformation after decomposition by different algorithms and its morphology, taking the simulated multi-year average deformation curve as the standard; S42: obtaining a root mean square error according to the selected seasonal deformation and the deformation information; Specifically, in step S42, the root mean square error RMSE is calculated using the residual value between the deformation decomposed by each algorithm and the original InSAR deformation. The larger the calculated value of RMSE is, the worse the performance of the corresponding method is. S43: obtaining an optimal decomposition method according to the multiple decomposition algorithms, using the morphological similarity and the root mean square error; and determining long-term deformation and seasonal deformation using the optimal decomposition method; S5: Obtaining layered soil moisture data, and obtaining the thickness of the active layer in the study area according to the seasonal deformation and the layered soil moisture data; Specifically, step S5 includes: obtaining stratified soil moisture data, and obtaining the active layer thickness of the study area according to the seasonal deformation and the stratified soil moisture data, such as the formula: , , in, is the thickness of the active layer, is the soil moisture layer number, For the The melting depth of the active layer within the soil layer, For the Seasonal sedimentation in the soil layer, and are the densities of pure water and ice, respectively, For the i The average soil moisture content of the layer, is a function of soil moisture content, is the soil depth, a, b, c, d and e are the coefficients of the polynomial fitting function; Specifically, based on seasonal deformation, the active layer thickness ALT of the study area was inverted based on the soil moisture content at different depths fitted by a polynomial function, and the accuracy of the method was verified using measured data; the active layer thickness ALT inversion model based on soil moisture content assumes that the seasonal surface settlement is caused by the volume change when the porous ice in the active layer is transformed into water; the relationship between frozen soil deformation and ALT can be expressed as:

[0027] In the formula, is the soil moisture content; and are the densities of pure water and ice respectively; It is the seasonal deformation of permafrost; is the active layer thickness ALT; The above formula is transformed by integration and the multi-layer soil moisture data is introduced to obtain the ALT expression:

[0028] In the formula, is the data soil moisture layer number; It is the melting depth of the active layer within a certain soil layer; is the seasonal sedimentation in a certain soil layer; No. i Mean soil moisture content in the layer; This embodiment introduces a polynomial on the basis of the ALT inversion model based on multi-layer soil moisture content. This polynomial function is intended to establish the relationship between underground soil moisture content and soil depth. By fitting the ERA5-Land layered soil moisture data of the climate reanalysis data, a smoother and continuous quartic polynomial is used to replace the original piecewise constant, which not only optimizes the continuity of the data, but also more accurately reflects the change of soil moisture content, so that the trend of soil moisture content changing with depth can be captured more accurately, so as to invert a more accurate ALT. The soil moisture polynomial fitting function established for the original data is shown as follows:

[0029] In the formula, is a function of soil moisture content; h is the soil depth; a, b, c, d and e is the coefficient of the function to be determined; Collect the measured ALT data in the study area, extract the inverted ALT at the same position as the measured point from the inversion results, and compare it with the measured ALT at the point to confirm the accuracy of the inversion.

[0030] In some embodiments, the above-mentioned method for inverting the thickness of the frozen soil active layer can also be implemented in the following manner.

[0031] Please refer to Figure 2 The present invention provides a frozen soil active layer thickness inversion method considering deformation and multi-layer soil moisture using Sentinel-1A images covering several consecutive years of the study area and the corresponding GACOS atmospheric delay data. Figure 3 .

[0032] In this embodiment, the inversion method of the thickness of the frozen soil active layer includes the following steps: Step 1: Obtain the Sentinel-1A Single Look Complex (SLC) descending orbit data and SRTM DEM data covering the study area, and download the GACOS atmospheric delay data covering the study area; Specifically, the area where the active layer thickness needs to be calculated is determined, and the Sentinel-1ASLC image data covering the study area and the matching GACOS atmospheric delay data are downloaded; Step 2: Perform interferometric pair combination, interferometric processing, track refinement and re-flattening, deformation inversion, and geocoding on the Sentinel-1A images obtained in step 1, and then obtain the final deformation information; Specifically, Step 2.1: Set the temporal baseline threshold of the SBAS-InSAR process to 60 days and the spatial baseline threshold to 150 meters, and generate interferometric pairs that meet the above thresholds; Step 2.2: Use SRTM data as external DEM to eliminate the terrain phase of the interferogram; perform multi-look operations of 8 pixels in the range and 2 pixels in the azimuth during interferometric processing to suppress noise and improve the signal-to-noise ratio of the interferogram; Step 2.3: After obtaining the interference fringe pattern, use GACOS data to process and interpret the data to eliminate or reduce the phase change caused by atmospheric delay in the source region of the Yellow River and reduce the impact of atmospheric delay; Step 2.4: Use the minimum cost flow (MCF) algorithm of the Delaunay triangulation network for phase unwrapping; after the interference processing is completed, the quality of the interference pair is judged according to the coherence coefficient diagram, interference diagram and unwrapping diagram, and the interference pairs with low coherence and poor unwrapping quality are removed to avoid affecting the inversion effect; Step 3: Six different decomposition algorithms, including Sin-Cos, X11, STL, BEAST, FastICA, and RobustSTL, are used to extract the long-term and seasonal deformations of permafrost; Specifically, Step 3.1: Six decomposition algorithms are applied to decompose the deformation. The mathematical models or solution steps of some algorithms are as follows: (1) Sin-cosine decomposition algorithm:

[0033] Where d is the cumulative deformation; v is the linear deformation rate; t is the cumulative time relative to the reference time when the image is acquired; T is the seasonal deformation period, which is usually 1 year; a1, a2, and b are the parameters to be determined for the model; (2) X11 decomposition algorithm: First, the original data is smoothed by Henderson weighted moving average to obtain the trend term; then, the Musgrave asymmetric moving average is used to fill the endpoint values ​​that cannot be accurately estimated by Henderson weighted moving average to capture long-term, medium-term and short-term seasonal changes; next, the seasonal factor is estimated for the adjusted observations to obtain the adjusted trend term; finally, the trend factor is multiplied by the seasonal factor to obtain the final trend term; the whole process is iterated until convergence is reached; (3) STL decomposition algorithm: The STL decomposition algorithm consists of two processes: an inner loop and an outer loop. The inner loop divides the time series into periodic subsequences by subtracting the estimated value of the current trend from the original sequence; each periodic subsequence is smoothed by LOESS and then processed by a low-pass filter to obtain the seasonal component; finally, the seasonal component is subtracted from the original data to obtain the remaining sequence; the remaining sequence is smoothed again with LOESS to obtain the final trend; LOESS smoothing curve g(x) The calculation method is as follows; assuming x i For the observation value, choose an integer q As the closest x The number of values, then when q ≤ n Time Based x i arrive x The distance is x i Set a neighborhood weight, the expression is:

[0034]

[0035] In the formula, W represents a weight function, The neighborhood weight, and x The distance between The longer, the weight The smaller; finally, at point Find a weight value The fitted value of g(x) ; The main function of the outer loop is to adjust the robust weights. If there are outliers in the sequence, the remainder will be larger, which can reduce the impact; (4) Decomposition algorithm based on independent component analysis: ICA (Independent Component Analysis) is a computational signal separation method that can effectively decompose mixed signals into a linear combination of a set of statistically independent signal sources. The algorithm first performs zero mean and whitening on the original deformation values ​​to generate an orthogonal matrix. Then, the mixed signal is linearly transformed and represented as an uncorrelated variable with a variance of 1 by whitening. The mixed matrix after zero mean and whitening preprocessing is represented by W, and the problem becomes:

[0036] In the formula, is the orthogonal adjustment matrix; then, the source matrix is ​​estimated by the following formula:

[0037] In the formula ， is the unmixing matrix W ; In order to obtain the unmixing matrix W, ICA must select an objective function to measure independence and find a suitable optimization algorithm; different objective functions directly affect the consistency and robustness of the ICA algorithm, and the selected optimization algorithm determines the convergence speed and stability of the algorithm; currently, the most widely used ICA algorithm is the fast fixed point algorithm (FastICA), which has a cubic linear convergence speed and is significantly better than other algorithms; therefore, using FastICA iterative calculation, by maximizing the spatial non-Gaussian source, the source matrix S, the mixing matrix W and the unmixing matrix U are derived; (5) BEAST algorithm: BEAST decomposition is performed through two sets of parameters, one for seasonality and the other for trend; seasonality and trend are modeled by piecewise harmonic models and piecewise linear functions, respectively, and the mutation points are implicitly encoded in the seasonal and trend signals; the seasonal signal S(t) is modeled as a linear combination of a set of harmonic terms, and the specific expression is as follows:

[0038] In the formula, P is the period of the seasonal signal; p is the number of segments of the seasonal signal, L k is the order of the harmonic model; and is the model coefficient, which represents the amplitude of the harmonic term; trend signal T(t) Modeling is done through piecewise linear function, the expression is:

[0039] In the formula, m is the number of trend change points, aj and b j are the model coefficients; (6) RobustSTL decomposition algorithm First, the input time series data is denoised using a bilateral filter, and then the trend is robustly extracted by solving the least absolute deviation (LAD) regression with sparse regularization; since the seasonal component is assumed to change slowly, a seasonal difference operation is performed on the denoised signal to reduce the seasonal effect; then a non-local seasonal filter is used to obtain the seasonal signal component; unlike the traditional method of estimating seasonality by considering only subsequences, this method uses K neighborhoods centered on various points; the weights assigned to these neighborhoods are based on temporal distance and seasonal similarity, making the algorithm more robust to outliers and seasonal changes; finally, an adjustment is made to ensure that the sum of all seasonal components within a cycle is zero; this adjustment is achieved by removing the mean, which corresponds to the estimation of the trend point; the algorithm iteratively improves the estimates of trend, seasonality, and residual components to improve accuracy; Step 4: Combining the two evaluation indicators, select the optimal algorithm to extract the long-term deformation and seasonal deformation of the frozen ground surface in the study area; Specifically, Step 4.1: Analyze the shape of the decomposed seasonal deformation curve to evaluate the performance of the algorithm; display the multi-year average of the same date of the time series seasonal deformation curve into one cycle, and observe the shape of the seasonal deformation curve obtained by different algorithms; Step 4.2: Using the simulated multi-year average deformation curve as the standard, the Dynamic Time Warping (DTW) algorithm is used to determine the degree of seasonal deformation after decomposition by different algorithms and its morphology; the calculation formula is as follows:

[0040] in D(i,j) Represents the position in the cumulative distance matrix (i,j) The value at d(x i ,y j ) Represents the Euclidean distance between two data points; Step 4.3: Calculate RMSE using the residual value between the deformation decomposed by each algorithm and the original InSAR deformation. The larger the calculated value of RMSE, the worse the performance of the corresponding method. Step 4.4: Combine the evaluation results of the two methods and select the optimal decomposition method; Step 5: Nine environmental factors including elevation, slope, and aspect were selected, and the geographic detector model was used to analyze the influence of the driving factors of long-term deformation and seasonal deformation; Specifically, Step 5.1: Geographic detectors are a set of statistical methods designed to detect spatial heterogeneity and reveal the driving forces behind it. The main idea is to divide the study space into multiple sub-regions according to the variables, and then compare the spatial variance within each sub-region and between different sub-regions to evaluate the determining power of potential explanatory variables. Its components include: (1) Differentiation and factor detection: Differentiation and factor detection is based on the q statistic to quantify the explanatory power of a factor X on the spatial differentiation of deformation. Its expression is:

[0041] in, q is the explanatory power of a factor on deformation, and its value is between 0 and 1. q The larger the value, the stronger the explanatory power of the factor on the deformation, and vice versa; h is the stratification of the factor, that is, the number of categories for classification or division; Nh and N Represents categories h and the number of samples across the region; Represents categories h and variance across the region; (2) Interaction detection: Interaction detection is used to calculate the interaction between any two factors, that is, to evaluate whether the explanatory power of factors X1 and X2 on the dependent variable (deformation) is enhanced or weakened when they act together. First, the q statistics of factors X1 and X2 on the deformation are calculated and marked as q(X1) and q(X2). Then, the interaction q statistic of factors X1 and X2 is calculated and marked as q(X1∩X2). By comparing the interaction q statistics of two factors and the q statistics of the two factors, the interaction relationship can be divided into five types, as shown in Table 1: Table 1: Comparison of interaction relationships

[0042] Step 6: Based on seasonal deformation, the ALT of the study area is inverted based on the soil moisture content at different depths fitted by the polynomial function, and the accuracy of the method is verified using measured data; Specifically, Step 6.1: The ALT inversion model based on soil moisture assumes that the seasonal surface settlement is caused by the volume change when the active layer pore ice is transformed into water; the relationship between frozen soil deformation and ALT can be expressed as:

[0043] In the formula, Mv is the soil moisture content; are the densities of pure water and ice respectively; It is the seasonal deformation of permafrost; H is ALT; the above formula is transformed by integration and the multi-layer soil moisture data is introduced to obtain the ALT expression:

[0044] In the formula, H is ALT; N is the data soil moisture layer number; It is the melting depth of the active layer within a certain soil layer; is the seasonal sedimentation in a certain soil layer; No. i Mean soil moisture content in the layer; This method introduces a polynomial on the basis of the ALT inversion model based on multi-layer soil moisture content; this polynomial function aims to establish the relationship between underground soil moisture content and soil depth; by fitting the ERA5-Land layered soil moisture content data, a smoother and continuous quartic polynomial is used to replace the original piecewise constant; this not only optimizes the continuity of the data, but also more accurately reflects the changes in soil moisture content; this improvement makes it possible to more accurately capture the trend of soil moisture content changing with depth, so as to invert a more accurate ALT; the soil moisture polynomial fitting function established for the original data is shown as follows:

[0045] In the formula, is a function of soil moisture content; h is the soil depth; a, b, c, d and e is the coefficient of the function to be determined; Step 6.2: Collect the measured ALT data in the study area, extract the inverted ALT at the same location as the measured point from the inversion results, and compare it with the measured ALT at the same location to confirm the accuracy of the inversion. For the ALT inversion results based on polynomial fitting and soil moisture content, please refer to Figure 4 .

[0046] This embodiment provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the above-mentioned frozen soil active layer thickness inversion method are implemented. Among them, the storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), a random access memory (RAM), a flash memory (Flash Memory), a hard disk (HDD) or a solid-state drive (SSD), etc.; the storage medium can also include a combination of the above-mentioned types of memory.

[0047] This embodiment provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps of the above-mentioned method for inverting the thickness of the frozen soil active layer are implemented.

[0048] like Figure 5 As shown, the computer device may include: at least one processor 121, such as a CPU (Central Processing Unit), at least one communication interface 123, a memory 124, and at least one communication bus 122. Among them, the communication bus 122 is used to realize the connection and communication between these components. Among them, the communication interface 123 may include a display screen (Display), a keyboard (Keyboard), and the optional communication interface 123 may also include a standard wired interface and a wireless interface. The memory 124 may be a high-speed RAM memory (Random Access Memory, volatile random access memory), or a non-volatile memory (non-volatile memory), such as at least one disk memory. The memory 124 may also be at least one storage device located away from the aforementioned processor 121. Among them, the memory 124 stores an application program, and the processor 121 calls the program code stored in the memory 124 to perform any of the above method steps. Among them, the communication bus 122 may be a peripheral component interconnect (PCI) bus or an extended industry standard architecture (EISA) bus, etc. The communication bus 122 can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 5Only one line is used to represent it, but it does not mean that there is only one bus or one type of bus. Among them, the memory 124 may include a volatile memory (English: volatile memory), such as a random access memory (English: random-access memory, abbreviated: RAM); the memory may also include a non-volatile memory (English: non-volatile memory), such as a flash memory (English: flash memory), a hard disk drive (English: hard disk drive, abbreviated: HDD) or a solid-state drive (English: solid-state drive, abbreviated: SSD); the memory 124 may also include a combination of the above-mentioned types of memory. Among them, the processor 121 may be a central processing unit (English: central processing unit, abbreviated: CPU), a network processor (English: network processor, abbreviated: NP) or a combination of CPU and NP. Among them, the processor 121 may further include a hardware chip. The above-mentioned hardware chip may be an application-specific integrated circuit (English: application-specific integrated circuit, abbreviated: ASIC), a programmable logic device (English: programmable logic device, abbreviated: PLD) or a combination thereof. The PLD can be a complex programmable logic device (CPLD), a field-programmable gate array (FPGA), a generic array logic (GAL), or any combination thereof. Optionally, the memory 124 is also used to store program instructions. The processor 121 can call the program instructions to implement the inversion method for the thickness of the frozen soil active layer in this embodiment.

[0049] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation modes, which are merely illustrative rather than restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the present invention and the claims, all of which are within the protection of the present invention.

Claims

1. A method for inverting the thickness of a frozen soil active layer, characterized in that: The following steps are involved: S1: Acquire single-view complex descending orbit data, digital elevation model data and atmospheric delay data covering the study area; S2: Obtaining deformation information according to the single-view complex descending orbit data, digital elevation model data and atmospheric delay data; S3: according to the deformation information, using a variety of decomposition algorithms, obtaining candidate long-term deformation and candidate seasonal deformation of permafrost; S4: Evaluate the candidate long-term deformation and the candidate seasonal deformation to determine the long-term deformation and the seasonal deformation; S5: Obtain stratified soil moisture data, and obtain the active layer thickness of the study area based on the seasonal deformation and the stratified soil moisture data.

2. The method for inverting the thickness of the frozen soil active layer according to claim 1 is characterized in that: The single-view complex descending orbit data is the single-view complex descending orbit data of the Sentinel-1A satellite, the digital elevation model data is the digital elevation model data of the Space Shuttle Radar Topography Mission, and the atmospheric delay data is the atmospheric delay data of the General Atmospheric Correction Online Service.

3. The method for inverting the thickness of the frozen soil active layer according to claim 1 is characterized in that: Step S2 specifically includes: S21: obtaining an interference pair according to the single-view complex descending orbit data by using a small baseline subset interferometry method according to a preset time baseline threshold and a preset space baseline threshold; S22: performing terrain phase elimination on the interference pair using the digital elevation model data to obtain an interference fringe pattern; S23: Processing and interpreting the interference fringe pattern using the atmospheric delay data to obtain a new interference pattern; S24: performing phase unwrapping on the new interference pattern using a minimum cost flow algorithm to obtain an unwrapping pattern; and obtaining deformation information according to the interference pair, the interference fringe pattern, the new interference pattern and the unwrapping pattern.

4. The method for inverting the thickness of the frozen soil active layer according to claim 3 is characterized in that: The preset time baseline threshold is 60 days, and the preset space baseline threshold is 150 meters.

5. The method for inverting the thickness of the frozen soil active layer according to claim 1 is characterized in that: The multiple decomposition algorithms include a sine-cosine decomposition algorithm, an X11 decomposition algorithm, an STL decomposition algorithm, a decomposition algorithm based on independent component analysis, a BEAST decomposition algorithm and a RobustSTL decomposition algorithm.

6. The method for inverting the thickness of the frozen soil active layer according to claim 1, characterized in that: Step S4 specifically includes: S41: obtaining a simulated multi-year average deformation curve, and obtaining morphological similarity using a dynamic time warping algorithm according to the simulated multi-year average deformation curve and the selected seasonal deformation; S42: obtaining a root mean square error according to the selected seasonal deformation and the deformation information; S43: According to the multiple decomposition algorithms, using the morphological similarity and the root mean square error, an optimal decomposition method is obtained; using the optimal decomposition method, long-term deformation and seasonal deformation are determined.

7. The method for inverting the thickness of the frozen soil active layer according to claim 6 is characterized in that: Step S41 specifically includes: obtaining a simulated multi-year average deformation curve, and obtaining a morphological similarity using a dynamic time warping algorithm according to the simulated multi-year average deformation curve and the selected seasonal deformation, such as the formula: ， in, Represents the position in the cumulative distance matrix The value at Represents the Euclidean distance between two data points.

8. The method for inverting the thickness of the frozen soil active layer according to claim 1, characterized in that: Step S5 specifically includes: obtaining layered soil moisture data, and obtaining the active layer thickness of the study area according to the seasonal deformation and the layered soil moisture data, such as the formula: ， ， in, is the thickness of the active layer, is the soil moisture layer number, For the The melting depth of the active layer within the soil layer, For the Seasonal sedimentation in the soil layer, and are the densities of pure water and ice, respectively, For the i The average soil moisture content of the layer, is a function of soil moisture content, is the soil depth, a, b, c, d and e are the coefficients of the polynomial fit function.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method for inverting the thickness of the frozen soil active layer as described in any one of claims 1 to 8 are implemented.

10. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the method for inverting the thickness of the active layer of frozen soil as described in any one of claims 1-8 are implemented.

Images