A method for identifying a water and heat control active area of an expansive soil slope

Abstract

The application discloses a kind of swelling soil side slope water heat control active area identification method, comprising: obtaining the vertical deformation time series of each pixel of target area and the rainfall time series of corresponding time period;Process the above two groups of time series, extract the periodic response characteristics representing the water-driven deformation process;Reflect the water-heat transfer model of the moisture diffusion mechanism and the depth integral deformation mechanism of swelling soil, to establish the physical mapping relationship between the theoretical prediction response characteristics and the normalized active area depth;Periodic response characteristics are substituted into water-heat transfer model to solve by inversion matching, and the total depth of each pixel is calculated by combining the moisture diffusion coefficient calibrated in advance.The present application effectively breaks through the theoretical upper limit of depth inversion of traditional single-layer time-domain model, realizes the accurate physical quantitative identification of active area depth under complex water-heat boundary.

Description

Technical Field

[0001] This invention belongs to the field of geotechnical engineering monitoring data processing, and in particular, it is a method for identifying the hydrothermal control activity zone of expansive soil slopes. Background Technology

[0002] Expansive soils, under the hydrothermal cycle of natural rainfall infiltration and evaporation drying, exhibit significant periodic volumetric expansion and contraction strain. Water-driven elastic deformation is concentrated within a specific geological depth range, known as the hydrothermal-controlled active zone. Accurately quantifying the depth of this active zone is a crucial physical and geometric boundary condition for constructing a coupled seepage-mechanism numerical model of slopes, determining the anchorage length of anti-slide piles, and laying out underground drainage networks.

[0003] Current methods for identifying the depth of active zones at the regional scale primarily utilize interferometric synthetic aperture radar (InSAR) technology to extract surface deformation sequences. The specific implementation typically involves: using polynomial fitting to eliminate long-term plastic trends in the deformation sequence, and calculating the effective rainfall sequence using a formula with a pre-set empirical attenuation coefficient; performing a cross-correlation residual search on the deformation sequence and the effective rainfall sequence in the time domain to obtain a single time-delay parameter of the rainfall response; directly equating this time delay to the water propagation time, substituting it into the analytical solution formula for a one-dimensional, single-layer, semi-infinite spatial diffusion front, and then calculating the depth of the active zone.

[0004] However, the aforementioned existing technologies suffer from theoretical limitations in inversion depth and a lack of decoupling in heterogeneous media conduction in their physical mapping mechanisms and computational frameworks. Specifically, existing schemes mistakenly equate the equivalent lag of the total deformation of the surface depth integral with the physical time of the water front propagating to the bottom. Since the phase delay of the depth integral deformation converges to a fixed quarter-pi constant as it diffuses into deeper layers, existing single-layer diffusion models have an insurmountable theoretical flaw in the upper limit of inversion depth, which may lead to numerical underestimation when dealing with deep active zones. In addition, existing technologies ignore the objective effects of preferential flow through surface shrinkage fissures in expansive soils and the filtering effect of surface water infiltration retention. The single mixed time delay extracted in the time domain alone cannot decouple the two-layer heterogeneous hydraulic response of the fissure network (in-phase) and the intact matrix (phase shift), and lacks a constraint mechanism for physical infiltration depth across the frequency domain, making the inversion equation prone to ill-conditioned solutions in complex porous media. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying the hydrothermal control activity zone of expansive soil slopes, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solution: A method for identifying hydrothermal control activity zones on expansive soil slopes, comprising:

[0007] Obtain the time series of vertical deformation of each pixel within the target slope area and the time series of rainfall within the corresponding time period;

[0008] Processing time series of vertical deformation and rainfall, extracting periodic response features characterizing the moisture-driven deformation process;

[0009] A hydrothermal transfer model reflecting the moisture diffusion mechanism and depth integral deformation mechanism of expansive soil was constructed. The hydrothermal transfer model established a physical mapping relationship between the theoretically predicted response characteristics and the normalized active zone depth.

[0010] The periodic response characteristics are substituted into the water and heat transfer model for inversion matching to obtain the normalized active zone depth. The total active zone depth of each pixel is then calculated by combining the pre-calibrated water diffusion coefficient.

[0011] Beneficial effects: (1) By constructing a fracture-enhanced double-layer hydrothermal transfer model, this invention effectively characterizes the heterogeneous transport mechanism where preferential flow from surface fractures and diffusion from deep intact matrix coexist in expansive soil, avoiding the mismatch between the traditional homogeneous model and the actual water migration process; (2) By establishing a depth integral deformation response transfer function, this invention breaks through the upper limit constraint of phase delay in the traditional single-layer diffusion front model in the deep active zone scenario, overcoming the problem that the depth of the active zone is easily underestimated; (3) By utilizing the multi-frequency periodic response characteristics to construct joint inversion constraints, this invention improves the well-posedness of the model parameter solution and the stability of the results, enabling the regional scale identification of the depth of the active zone of expansive soil slopes without relying on high-density borehole monitoring. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of a method for identifying hydrothermal control activity zones on expansive soil slopes, as provided in this application embodiment.

[0013] Figure 2 A flowchart illustrating the steps for extracting and characterizing the periodic response features of a water-driven deformation process, as provided in an embodiment of this application.

[0014] Figure 3 A flowchart illustrating the steps for obtaining the normalized fracture depth and the normalized thickness of the intact matrix layer through inversion, as provided in the embodiments of this application.

[0015] Figure 4 A flowchart illustrating the steps of reliable pixel screening using a preset threshold, as provided in an embodiment of this application. Detailed Implementation

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

[0017] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0018] like Figure 1 As shown, a method for identifying hydrothermal control activity zones on expansive soil slopes includes the following steps:

[0019] Obtain the time series of vertical deformation of each pixel within the target slope area and the time series of rainfall within the corresponding time period.

[0020] In this embodiment, the target slope area refers to the distribution area of ​​expansive soil with potential risks of swelling, shrinkage, deformation, or instability. A pixel is the basic spatial unit for remote sensing observation. The vertical deformation time series represents the sequence of uplift or subsidence of the Earth's surface over time in the direction of gravity; this time series can be obtained using satellite radar observation technology. The rainfall time series consists of meteorological data aligned with the deformation observation time within the same spatial range. Because expansive soils are water-sensitive, rainfall infiltration and subsequent evaporation are the core boundary conditions driving soil volume expansion and contraction, while vertical surface deformation is its macroscopic mechanical response. Therefore, obtaining both the vertical deformation time series and the rainfall time series—two time- and spatially aligned sequences—is the data foundation for establishing causal response correlations.

[0021] In some alternative implementations, the rainfall time series can be obtained by mapping daily rainfall data measured at surrounding meteorological stations to pixel locations using a spatial interpolation algorithm, or directly extracted from multi-source satellite remote sensing rainfall products.

[0022] By processing time series of vertical deformation and rainfall, periodic response features characterizing the moisture-driven deformation process are extracted.

[0023] Specifically, periodic response characteristics refer to the rhythmic changes in deformation signals caused by seasonal rainfall, which can be expressed as the phase lag or amplitude evolution of deformation relative to rainfall over time. Processing vertical deformation and rainfall time series involves using signal analysis to separate and purify the elastic expansion and contraction components directly caused by water turnover. The original total deformation of expansive soil slopes typically combines long-term irreversible plastic slippage, high-frequency environmental noise, and hydrothermal-driven periodic elastic expansion and contraction. Extracting periodic response characteristics that solely characterize water-driven processes effectively isolates geological interference from non-hydrological factors, allowing subsequent model inversion calculations to focus on water transport patterns within the soil, thus improving the signal-to-noise ratio for parameter identification.

[0024] Optionally, the extraction of periodic response features can be achieved through cross-correlation search in the time domain or harmonic decomposition in the frequency domain, to adapt to processing requirements under different data signal-to-noise ratio conditions. In specific implementation, when the target area has a stable annual or semi-annual rainfall rhythm and the periodic signal-to-noise ratio of the deformation time series is high, the frequency domain harmonic decomposition method is preferred for extracting periodic response features. Its input is the aligned vertical deformation time series and rainfall time series, and the output is the deformation phase, deformation amplitude, rainfall phase, and rainfall amplitude at at least one preset frequency. When the periodic rhythm of the target area is not significant, the frequency domain high-order harmonic signal is weak, or only a short time series can be obtained, the time domain cross-correlation or minimum residual search method can be used to identify the rainfall response lag time. Its input is the elastic deformation sequence and rainfall driving sequence after trend term removal, and the output is a single time lag parameter.

[0025] A hydrothermal transfer model reflecting the moisture diffusion mechanism and depth integral deformation mechanism of expansive soil was constructed. The hydrothermal transfer model established a physical mapping relationship between the theoretically predicted response characteristics and the normalized active zone depth.

[0026] Alternatively, based on the water diffusion mechanism and depth integral deformation mechanism of expansive soil, a hydrothermal transfer model is constructed to reflect the propagation of water to deeper layers in the active zone and the induction of depth integral deformation. The hydrothermal transfer model establishes a physical mapping relationship between theoretically predicted response characteristics and normalized active zone depth.

[0027] In other words, a hydrothermal transfer model reflecting the water diffusion mechanism and depth integral deformation mechanism of expansive soil is constructed. The hydrothermal transfer model establishes a physical mapping relationship between the theoretically predicted response characteristics and the normalized active zone depth based on the diffusion analytical solution and depth integral deformation response under periodic boundary conditions.

[0028] In this embodiment, the physical mapping relationship is not an empirical fitting relationship, but is derived jointly from the analytical solution of water diffusion under periodic boundary conditions and the depth integral deformation relationship. For example, the amplitude attenuation and phase lag expressions of water disturbance at different depths are obtained from the diffusion equation. The volumetric strain caused by the disturbance is integrated along the depth direction of the active zone to form a complex domain response expression of the total surface deformation. The theoretical phase delay is obtained by extracting the argument of the complex response, and a one-to-one correspondence between the theoretical phase delay and the normalized active zone depth is established.

[0029] The hydrothermal transfer model is a mathematical and physical representation system for the heterogeneous water migration process of expansive soil slopes. It uses surface moisture changes as the driving boundary and the fracture-matrix composite medium developed along the depth direction as the transport medium. For expansive soil slopes, rainfall infiltration does not diffuse uniformly in homogeneous soil, but rather forms rapid infiltration channels along surface drying shrinkage fractures, and continues to diffuse into the surrounding intact matrix layer at the bottom of the fractures. Therefore, the water diffusion mechanism is an equivalent diffusion mechanism under the coupling effect of rapid fracture transport and slow matrix diffusion. Its diffusion process can still be physically solved within a one-dimensional vertical framework, but the corresponding diffusion parameters are equivalent diffusion parameters reflecting the combined effect of fracture conduction and matrix diffusion. The depth integral deformation mechanism shows that the total surface deformation observed by satellite is not a deformation of a single shallow depth, but rather the cumulative integral of the volumetric strain of all soil layers from the surface downwards to the bottom of the water-affected active zone. The normalized active zone depth characterizes the relative proportional relationship between the periodic water infiltration characteristic scale and the actual physical depth of the active zone.

[0030] By constructing a hydrothermal transfer model, empirically dependent statistical fitting relationships can be transformed into mechanistic expressions with well-defined physical boundary conditions. This model not only describes the dynamic attenuation of water as it propagates deeper into the soil, but also explicitly demonstrates the accelerating effect of surface fissures on the infiltration process and the coupling transfer effect between fissures and the intact matrix. This establishes a physical correspondence between the periodic deformation response observed at the surface and the depth of the underground active zone, consistent with the actual engineering mechanisms of expansive soil.

[0031] As an example, the hydrothermal transfer model can be set as a homogeneous single-layer semi-infinite spatial physical structure, or, in order to more accurately simulate the drying shrinkage cracks in shallow expansive soils in nature, it can be set as a double-layer or multi-layer physical spatial structure containing different hydraulic conduction characteristics.

[0032] It should be noted that under the boundary conditions of periodic rainfall, the moisture disturbance in expansive soil is not instantaneously transmitted to the entire depth of the active zone, but gradually attenuates with increasing depth, resulting in a phase delay. Simultaneously, the deformation observed at the surface is not a local response at a single depth, but rather the total response after integrating the volumetric strain of all depth layers within the active zone along the depth direction. Therefore, extracting the phase delay of surface deformation relative to rainfall input essentially involves extracting the overall frequency response characteristics of the fracture-matrix transport system to periodic moisture input, and this frequency response characteristic has an analytically resolvable physical mapping relationship with the depth of the active zone.

[0033] The periodic response characteristics are substituted into the water and heat transfer model for inversion matching to obtain the normalized active zone depth. The total active zone depth of each pixel is then calculated by combining the pre-calibrated water diffusion coefficient.

[0034] Alternatively, the periodic response characteristics can be substituted into the hydrothermal transfer model for inversion matching, and the normalized fracture depth and the normalized thickness of the intact matrix layer can be obtained by joint inversion. The total depth of the active area of ​​each pixel can then be recovered by combining the pre-calibrated water diffusion coefficient.

[0035] For example, the extracted periodic response features are used as the input of the observation end, and the theoretical phase delay and / or theoretical amplitude ratio output by the water and heat transfer model at the corresponding frequency are used as the output of the model end. By constructing the residual function between the observation end and the model end for inversion matching, the parameter set that minimizes the residual function is searched to obtain the normalized active area depth. Then, combined with the pre-calibrated water diffusion coefficient and the scale conversion relationship defined by the model, the total active area depth of each pixel is calculated.

[0036] Specifically, inversion matching is a numerical optimization process, which searches for a set of solutions in a multidimensional parameter space such that the theoretically predicted response characteristics of the positive output of the hydrothermal transfer model infinitely approximate or even equal the periodic response characteristics actually extracted. The pre-calibrated moisture diffusion coefficient represents the inherent physical property of water conduction in the soil of the target area. After obtaining the normalized active zone depth by solving the matching equation, algebraic operations are performed on it with the feature variables containing the moisture diffusion coefficient, according to the physical definition of the model, to reconstruct the total depth of the active zone with actual spatial length dimensions. It should be noted that the moisture diffusion coefficient is an equivalent diffusion parameter that comprehensively reflects the conductivity of surface fissures and the diffusion of the underlying intact matrix in expansive soil, and is not a single diffusion parameter under the traditional homogeneous medium assumption.

[0037] This achieves a closed-loop transformation from dimensionless mathematical solutions to actual engineering physical quantities. By combining the pre-calibrated moisture diffusion coefficient for a specific region for conversion, it ensures that the total depth of the active area obtained by inversion truly reflects the local soil hydraulic differences, providing direct geometric boundary parameters for slope stability evaluation and drainage protection engineering design.

[0038] Furthermore, the process of obtaining the normalized active region depth can be performed using iterative optimization algorithms such as grid search or gradient descent, and the physical rationality of the calculation results can be screened according to preset quality control rules after the inversion output.

[0039] According to another aspect of this application, a method for identifying hydrothermal control activity zones on expansive soil slopes may further include: acquiring the vertical deformation time series of each pixel within the target slope area and the rainfall time series within the corresponding time period; performing frequency domain processing on the vertical deformation time series and the rainfall time series to extract the phase delay characterizing the water-driven deformation process at at least two preset frequencies; constructing a fracture-enhanced two-layer model, wherein the fracture-enhanced two-layer model is sequentially divided into a fracture layer and a intact matrix layer along the depth direction of the active zone, and a zero-flux boundary condition is applied at the bottom of the active zone, wherein the water response in the fracture layer is in phase with the surface driving signal, and the water in the intact matrix layer propagates through diffusion. A phase delay is generated. Based on the fracture-enhanced two-layer model, a complex numerical depth integral deformation response transfer function is established, thereby establishing a physical mapping relationship between the theoretically predicted phase delay and the normalized active zone depth parameters. The normalized active zone depth parameters include the normalized fracture depth and the normalized thickness of the intact matrix layer. The phase delay at at least two preset frequencies is substituted into the depth integral deformation response transfer function for joint inversion to obtain the normalized fracture depth and the normalized thickness of the intact matrix layer for each pixel. The characteristic depth of water infiltration is calculated according to the pre-calibrated water diffusion coefficient, and the total active zone depth of each pixel is recovered based on the normalized fracture depth and the normalized thickness of the intact matrix layer. This embodiment can realize regional-scale physical quantitative identification of the active zone depth of expansive soil slopes.

[0040] In one possible implementation, the vertical deformation time series of each pixel within the target slope area is obtained, including:

[0041] Obtain the cumulative deformation sequence of the synthetic aperture radar line of sight and the corresponding radar incident angle for each pixel.

[0042] In this embodiment, the synthetic aperture radar line-of-sight cumulative deformation sequence refers to the set of relative surface displacements continuously measured by the radar sensor at multiple observation time points along the line-of-sight direction of electromagnetic wave transmission and reception. The radar incident angle is the angle between the radar line-of-sight direction and the surface normal direction at the location of the target pixel.

[0043] To obtain a synthetic aperture radar (SAR) line-of-sight cumulative deformation sequence with engineering application accuracy, small baseline set interferometry (SLA) or permanent scatterer interferometry (PSI) is typically used to perform time-series processing on multi-temporal single-look complex image sets covering the target slope area. During differential interferometry and phase unwrapping, quality control parameters need to be set. Specifically, the coherence threshold is usually set between 0.7 and 0.8 to remove unreliable pixels with dense vegetation cover or in radar shadow areas. The spatial baseline threshold can be set to less than 200 meters, and the temporal baseline threshold can be set to less than 60 days. By strictly controlling the spatiotemporal decorrelation effect and assisting in the removal of terrain phase using a digital elevation model, speckle noise in the interferometric spectrum can be suppressed, enabling the final output sequence to capture millimeter-level minute surface displacement signals.

[0044] Based on the assumption that the periodic expansion and contraction of expansive soil is mainly in the vertical direction, the cumulative deformation sequence of the synthetic aperture radar line of sight is geometrically projected and transformed using the radar incident angle to calculate the vertical deformation time series that eliminates the influence of horizontal deformation.

[0045] Specifically, for gently sloping or compacted expansive soil slopes, the volumetric strain under seasonal wet-drying alternation mainly manifests as vertical uplift and settlement along the gravity direction, while the horizontal slip component is usually negligible. Geometric projection transformation refers to using spatial trigonometric transformation relationships to decompose and map the tilted deformation vector in the one-dimensional line-of-sight direction to the vertical direction.

[0046] For example, the specific calculation for geometric projection transformation using the radar incident angle can be implemented using the following formula:

[0047] d v (x, t) k )=d LOS (x, t) k ) / cosθ inc (x);

[0048] Where, d v (x, t) k Let ) be the pixel x at observation time t. k The cumulative deformation in the vertical direction, d LOS (x, t) k θ represents the corresponding acquired synthetic aperture radar line-of-sight cumulative deformation observation value. inc (x) represents the radar incident angle, and cos represents the cosine function operation.

[0049] The subsequent hydrothermal transfer physics model was built under one-dimensional vertical diffusion boundary conditions. The original line-of-sight observations were calibrated to vertical deformation, eliminating the systematic dimensional shifts caused by satellite orbital skew and radar side-view imaging geometry, and establishing the correspondence between surface deformation amplitude and subsurface infiltration depth in a unified one-dimensional vertical physical coordinate system.

[0050] In some alternative implementations, a meteorological input processing scheme is provided for acquiring rainfall time series data, using net infiltration as an alternative to raw rainfall data. For expansive soil slopes with exposed surfaces or significant evaporation, directly using raw rainfall as the system's water-driven boundary can lead to an overestimation of the effective water absorbed by the soil. In this case, the net infiltration time series can be calculated as the input entity for the rainfall time series.

[0051] Specifically, daily meteorological observation data of the study area are obtained, and basic indicators such as temperature, relative humidity, and wind speed are extracted. The potential evapotranspiration of the target area is calculated using a pre-set potential evapotranspiration model. The formula for calculating net infiltration based on potential evapotranspiration is as follows:

[0052] I net (t)=P(t)-η*ET0(t);

[0053] Among them, I net P(t) represents the net infiltration at a certain moment, P(t) represents the original rainfall at the same moment, ET0(t) represents the potential evapotranspiration, and η represents the surface evaporation reduction factor characterizing the surface vegetation or bare soil state. The preset surface evaporation reduction factor η is usually assigned a value of 0.6 to 0.8 for bare expansive soil slopes.

[0054] By calculating net infiltration, the water flux consumed by surface evaporation can be deducted in advance at the input of the driving signal. When the net infiltration value is less than zero, it objectively reflects the dryness of the surface due to evaporation. Since the evaporation process and the rainfall process have natural antiphase interference, using net infiltration to replace the original rainfall can naturally eliminate some environmental phase interference before subsequent feature extraction, thereby improving the quality of the input source of the inversion model.

[0055] like Figure 2 As shown, in an exemplary embodiment, the periodic response features include phase delay and amplitude ratio at at least two preset frequencies; extracting periodic response features characterizing the moisture-driven deformation process includes:

[0056] Frequency domain decomposition was performed on the vertical deformation time series and the rainfall time series to extract the deformation phase, deformation amplitude, rainfall phase and rainfall amplitude at each preset frequency.

[0057] In this embodiment, frequency domain decomposition refers to using mathematical transformation algorithms to convert a discrete observation sequence arranged in chronological order into a spectral energy map distributed by frequency. Specific transformation operations can be implemented by calling a discrete Fourier transform algorithm or a harmonic least squares fitting regression engine. The preset frequency refers to a reference working frequency artificially set according to the climate and meteorological patterns of the expansive soil region. For example, the preset frequency can be set as an annual cycle frequency and a semi-annual cycle frequency. The time span corresponding to the annual cycle frequency is one standard Earth year, and its angular frequency value is 2π. The time span corresponding to the semi-annual cycle frequency is half an Earth year, and its angular frequency value is 4π.

[0058] Furthermore, after performing frequency domain decomposition on the two sets of time series, the complex domain characteristic coefficients of each series at the annual and semi-annual cycle frequency points are extracted. The deformation phase and rainfall phase correspond to the argument of the complex domain characteristic coefficients, representing the initial time shift of the waveform of that frequency component. The deformation amplitude and rainfall amplitude correspond to the modulus of the complex domain characteristic coefficients, representing the absolute half-amplitude of the waveform of that frequency component.

[0059] In some optional implementations, considering that there may be significant temperature fluctuations in the target slope area, and that temperature fluctuations can indirectly affect soil moisture content by controlling the surface evaporation rate, a temperature effect subtraction step can be added. Specifically: Obtain the average air temperature sequence within the same time window as the vertical deformation time series. Perform the same frequency domain decomposition on the air temperature sequence to extract the temperature phase. Using a pre-configured multiple linear physical regression equation, separate and remove the deformation offset caused by the air temperature sequence from the total response of the vertical deformation time series. Perform frequency domain decomposition using the updated residual deformation sequence to cut off the phase interference of non-rainfall meteorological factors.

[0060] Preferably, before performing frequency domain decomposition, the vertical deformation time series and the rainfall time series are first resampled with a unified time step. Specifically, using a preset time resolution as the unified sampling interval, the vertical deformation observations with unequal time intervals are mapped onto a unified time axis through linear interpolation, spline interpolation, or by keeping the most recent observation; the rainfall time series is accumulated or resampled along the same time axis to form a synchronous discrete sequence that can be used for frequency domain decomposition.

[0061] The corresponding phase delay is calculated based on the deformation phase and rainfall phase at the same preset frequency.

[0062] In this embodiment, phase delay characterizes the degree of response lag on the time axis caused by the internal pore water pressure transmission and solid skeleton volumetric strain accumulation processes of the soil after receiving surface water infiltration boundary conditions. Within the frequency domain analysis framework, the degree of response lag is directly mapped to the angle difference.

[0063] Specifically, at the set annual cycle frequency node, the annual cycle deformation phase is subtracted from the annual cycle rainfall phase to obtain the preliminary angle difference. If the preliminary angle difference is negative or exceeds the standard circumferential range, a modulo 2π modulo operation is performed to forcibly calibrate the difference value and truncate it to the radian range of 0 to 2π. The above subtraction and truncation operations are performed simultaneously at the semi-annual cycle frequency node. Finally, independent and normalized phase delay data at two specific frequency scales are output.

[0064] The corresponding amplitude ratio is calculated based on the deformation amplitude and rainfall amplitude at the same preset frequency.

[0065] Specifically, the amplitude ratio is the proportionality coefficient between the output response intensity and the input driving intensity of the dynamic system, characterizing the transfer gain or attenuation characteristics of hydrothermal expansion energy at that frequency. For example, the calculation process involves dividing the deformation amplitude at the annual cycle frequency node by the rainfall amplitude at the corresponding node using a scalar division operation. This division calculation process is repeated for the semi-annual cycle frequency. Frequency domain decomposition technology is used to directly lock the phase delay and amplitude ratio at a specific frequency, avoiding the influence of high-frequency, high-order clutter generated by the time-domain superposition of complex meteorological systems.

[0066] In another exemplary embodiment, for data scenarios involving unimodal rainfall climate zones or extremely weak semi-annual cycle signals, the periodic response characteristic is the rainfall response lag time; the periodic response characteristics characterizing the water-driven deformation process are extracted, specifically through the following methods:

[0067] By removing the long-term plastic deformation component from the vertical deformation time series using a preset fitting function, a hydrothermal-driven elastic deformation time series is obtained.

[0068] In this embodiment, the surface deformation of the natural slope is a mixture of long-term gravitational creep slip and hydrothermal expansion and contraction cycles. The fitting function uses a cubic polynomial or other appropriate order polynomial equation to construct the time variable for baseline regression. After determining the coefficients of each term in the polynomial using the least squares inversion method, the corresponding long-term plastic deformation trend line is generated. The original vertical deformation time series is subtracted point by point from the long-term plastic deformation trend line to output a hydrothermal driven elastic deformation time series with the trend term eliminated.

[0069] The effective rainfall time series is obtained by calculating the rainfall time series using a preset attenuation formula.

[0070] Specifically, due to the accompanying surface runoff and capillary retention effects during rainfall infiltration, the contribution of historical rainfall to the current soil moisture content decreases non-linearly. The preset attenuation formula employs a cumulative algorithm incorporating an exponential attenuation factor for discrete calculation. The calculation formula can be:

[0071] P e (t)=∑[P(ti)*e(-λ*i) ];

[0072] Among them, P e P(t) represents the effective rainfall at time t, P(ti) represents the actual daily rainfall recorded on the i-th day prior to the time step, λ is the pre-configured rainfall attenuation coefficient, and e is the base of the natural logarithm. ∑ represents the summation recursively from i equal to 0 towards the historical time period with a step size of days, and the final number of days for the summation is limited by a preset upper limit of an empirical constant.

[0073] It should be noted that the rainfall attenuation coefficient λ can be determined by exponentially fitting historical rainfall-moisture content response data based on the soil water holding capacity and evaporation intensity of the target area. Those skilled in the art can select an appropriate value based on actual hydrological conditions. In one optional embodiment, the value of λ is between 0.01 and 0.10 days. -1 .

[0074] Based on a pre-configured time-domain linear hydrothermal driven model, a minimum residual search is performed on the time series of hydrothermal driven elastic deformation and the time series of effective rainfall to identify the rainfall response lag time.

[0075] In this embodiment, the time-domain linear hydrothermal driven model treats deformation as a lagged linear mapping of effective rainfall. A minimum residual search process is established using a daily iterative test procedure. A point-by-point scan is performed within a one-dimensional parameter space from day 0 to a preset upper limit of the search, shifting the entire effective rainfall time series historically by the given number of test days each time. The sum of squared residuals between the shifted series and the hydrothermal driven elastic deformation time series is calculated. The test day that causes the sum of squared residuals to reach its global minimum is identified and confirmed as the rainfall response lag time.

[0076] Furthermore, based on the principle of unified physical dimensions, a forced dimensional balancing operation is performed before outputting the rainfall response lag time for subsequent physical diffusion depth calculations. The calculation engine divides the extracted rainfall response lag time, which is in days, by the calendar constant 365.25, converting its numerical scale to a physical quantity in calendar years. This time dimension unification conversion mechanism eliminates the potential for computational crashes caused by substituting daily-scale time delays into annual-scale angular frequency partial differential diffusion models.

[0077] In this embodiment, the identified rainfall response lag time is considered as the equivalent propagation time under a single-layer equivalent diffusion framework and is used to approximate the depth of the active area. The rainfall response lag time reflects the comprehensive delayed response of surface deformation to periodic water input, rather than the arrival time of a single water front in a strict sense. Therefore, the resulting depth should be understood as a preliminary estimate under the assumption of single-layer equivalent diffusion.

[0078] According to one aspect of this application, the hydrothermal transfer model is a fracture-enhanced two-layer model; the fracture-enhanced two-layer model is divided into a fracture layer and a intact matrix layer along the depth direction of the active zone, and a zero-flux boundary condition is applied at the bottom of the active zone; wherein, the water response in the fracture layer is set to be in phase with the surface driving signal, and the water in the intact matrix layer is propagated through diffusion flow and generates a phase delay.

[0079] In this embodiment, the fractured layer is used to characterize the preferential infiltration channels formed by the drying shrinkage fractures on the surface of expansive soil, while the intact matrix layer is used to characterize the water migration zone below the fracture bottom, dominated by matrix diffusion. A physical framework for the one-dimensional boundary value problem is constructed through the vertical layering abstraction of the geological engineering profile. The fracture-enhanced two-layer model is a direct mathematical model of the natural property of expansive soil surface development of drying shrinkage fractures under alternating wet and dry conditions. The two-layer structure is adopted because expansive soil easily forms a network of drying shrinkage fractures on the surface under long-term alternating wet and dry conditions. The hydraulic conductivity of the fractured zone is significantly higher than that of the underlying intact soil matrix, resulting in a hierarchical characteristic of rapid surface preferential flow followed by slow diffusion in the deep matrix during rainfall infiltration. This hierarchical transport characteristic is the key hydrothermal response mechanism that distinguishes expansive soil slopes from homogeneous soil slopes, and it is also the physical basis for establishing an active zone depth inversion model.

[0080] Along a vertical axis from the surface to the subsurface, the fracture-enhanced two-layer model is divided into two heterogeneous control units. The upper layer is defined as the fracture layer, which is filled with interconnected desiccation macropore channels, forming a preferential flow network for rainfall infiltration. Within the fracture layer, the velocity of water conduction is much greater than the capillary diffusion velocity in conventional soil pores. Preferably, a quantitative criterion is set: when the ratio of the equivalent water diffusion coefficient of the fracture layer to the diffusion coefficient of the lower intact layer exceeds a preset constant, for example, more than two orders of magnitude, or when the hydraulic conduction characteristic time of the fracture layer is less than a preset quantile value of the periodic rainfall frequency, the hydrodynamic process of the fracture layer is considered to be in a quasi-steady state. Accordingly, it is assumed that the changes in matrix suction at each elevation node within the fracture layer and the surface rainfall boundary driving signal are not lagging behind each other on the time axis, i.e., they are in phase.

[0081] The region extending downwards from the bottom boundary of the fractured layer to the ultimate depth of the active zone is defined as the intact matrix layer. The intact matrix layer remains undamaged by desiccation cracks, preserving the original dense structure of the expansive soil. Moisture migrates slowly within the intact matrix layer according to Darcy's law and a one-dimensional diffusion partial differential equation controlled by pore water pressure. Due to diffusion resistance, the arrival time of the moisture disturbance wavefront exhibits an analytically calculable phase delay for each specific depth the wavefront advances downwards.

[0082] Furthermore, the bottom of the active zone refers to the maximum vertical depth extreme point where seasonal hydrothermal alternation has a substantial impact on soil moisture content. When assembling the model boundary conditions, the computational engine sets a zero-flux lower boundary condition with zero partial derivatives at the elevation coordinates of the bottom of the active zone. This boundary condition cuts off the downward transmission path of hydrological signals, indicating that the deep bedrock or stable soil layers below the active zone are in a constant moisture content state and no longer experience volume expansion or contraction caused by meteorological alternation.

[0083] In a preferred embodiment, the fracture-enhanced two-layer model establishes a physical mapping relationship through a depth integral deformation response transfer function containing complex values, the analytical expression of which is:

[0084] Φ bi (ξ2, r)=r+(1-r)*tanh[(1+i)ξ2] / [(1+i)ξ2];

[0085] Where, Φ bi ξ represents the theoretically predicted response characteristics; r represents the normalized fracture depth; ξ² represents the normalized thickness of the intact matrix layer; i represents the imaginary unit; tanh represents the hyperbolic tangent function.

[0086] The normalized active zone depth is composed of the normalized fracture depth to be inverted and the normalized thickness of the intact matrix layer.

[0087] It should be noted that Φ bi This refers to the complex numerical depth integral deformation response transfer function corresponding to the fracture-enhanced two-layer model, whose modulus |Φ bi | represents the theoretically predicted amplitude response, whose argument arg(Φ) bi ) indicates the theoretically predicted phase delay.

[0088] In this embodiment, a closed-form analytical solution for the physical boundary value problem in the periodic frequency domain is presented. By integrating the local volumetric strain of all infinitesimal layers from the bottom up along the depth coordinate axis, a total transfer matrix incorporating complex variable function properties is obtained. Specifically, the imaginary unit *i* is introduced because wave attenuation and phase lag are naturally coupled in the frequency domain. The tanh function is the necessary functional form derived from the partial differential equation in the complex plane after applying the zero-flux bottom boundary condition. The parameter *r* is the absolute ratio of the actual fracture depth boundary to the total depth boundary of the active zone, and its value is limited to the dimensionless open interval between 0 and 1. When *r* approaches 1, it indicates that the fracture penetrates the entire active zone, and the deformation of the entire soil layer is synchronized with rainfall; when *r* approaches 0, it indicates that the soil is homogeneous and without fractures. The parameter *ξ2* is the absolute ratio of the actual intact matrix layer absolute thickness to the characteristic depth of water infiltration corresponding to the current driving frequency, reflecting whether the intact soil layer appears extremely thick or thin relative to a specific rainfall cycle.

[0089] The fracture-enhanced two-layer model locks the normalized fracture depth *r* and normalized thickness *ξ2* within its computational framework. After obtaining the optimal solutions for the normalized fracture depth *r* and normalized thickness *ξ2* through subsequent joint inversion, the overall normalized active zone depth variable covering the entire profile can be obtained through simple algebraic transformations and assembly reconstruction, based on the definition of geophysical processes. Therefore, the active zone depth obtained through inversion is not the apparent diffusion depth under the traditional homogeneous diffusion assumption, but rather the equivalent hydrothermal influence depth after comprehensively considering the combined effects of rapid infiltration through expansive soil fractures and slow diffusion of the matrix, which better aligns with the engineering implications of the actual hydrothermal control range of expansive soil slopes.

[0090] In a further embodiment, the hydrothermal transfer model establishes a physical mapping relationship based on the depth integral deformation response transfer function, including:

[0091] Extract the real and imaginary parts of the depth integral deformation response transfer function in the complex plane;

[0092] Calculate the ratio of the negative value of the imaginary part to the real part, and extract the argument using the arctangent function to obtain the theoretical phase delay predicted by the model;

[0093] By establishing an equality constraint between the theoretical phase delay and the observed phase delay in the periodic response characteristics, the inversion matching of the normalized active region depth is achieved.

[0094] Alternatively, by establishing an equality constraint between the theoretical phase delay and the observed phase delay in the periodic response characteristics, the inversion matching of the normalized fracture depth and the normalized thickness of the intact matrix layer can be achieved.

[0095] This embodiment describes the transformation of the transfer function output, which has an abstract complex form, into a physical scalar that can be directly compared with satellite-measured deformation data. During the extraction operation, the computational engine utilizes the complex hyperbolic tangent identity and the principle of conjugate operation to transform the output term Φ... bi The explicit stripping and recombination process results in a pure real scalar field and an orthogonal scalar field carrying an imaginary unit factor. Upon entering the argument extraction module, the computational engine negates the entire stripped imaginary numerical field. A scalar division operation is then invoked, using the negated imaginary value as the numerator and the real value as the denominator to obtain the tangent ratio. A monotonically increasing arctangent mapping operation is applied to the tangent ratio. The calculated arctangent output radian value is the pure physical quantity without any complex features, i.e., the theoretical phase delay. The arctangent operation establishes a rigorous analytical channel from the partial differential diffusion space vector of moisture in the porous medium to the scalar time delay of surface deformation. In the actual parameter matching stage, this channel, along with the measured phase delay at a specific frequency extracted by InSAR, is substituted into both sides of the equation, forcing the equation to converge and be solved.

[0096] It should be noted that, in the frequency domain representation, the total surface deformation can be written as a complex periodic response, with its real part corresponding to the component in phase with the rainfall input and its imaginary part corresponding to the lagged component. Since the deeper the active zone, the longer it takes for moisture disturbance to propagate and accumulate in the deeper soil layers, the amplitude of the total response after depth integration systematically changes with the increase of the normalized active zone depth. Therefore, the theoretical phase delay is not simply a result of data post-processing, but rather a frequency domain projection of the groundwater migration depth structure onto the surface deformation. Thus, by matching the observed phase delay with the theoretical phase delay, the depth of the active zone can be inverted.

[0097] In some alternative implementations, for geological environments with light shallow weathering or dominated by undisturbed soil and lacking a significant preferential flow network, a degraded alternative hydrothermal model construction scheme is provided. Specifically, the hydrothermal transfer model is a single-layer homogeneous semi-infinite space model. The single-layer homogeneous semi-infinite space model assumes that water spreads to deeper layers within the boundary according to a one-dimensional diffusion law. The phase delay of the corresponding surface depth integral deformation exhibits a nonlinear monotonically increasing trend with the increase of the normalized active zone depth and converges to the upper limit of the constant.

[0098] In this embodiment, the fracture layer configuration parameters are not considered. The entire active region is assumed to be a homogeneous porous medium with only a single diffusion coefficient, and the zero flux cutoff at the bottom is abandoned in the mathematical derivation. Instead, the diffusion boundary is extended to the negative infinity direction of the coordinate axis to simplify the analytical expression and reduce computational overhead.

[0099] However, the single-layer homogeneous semi-infinite space model has inherent physical limitations. As shown in the diffusion differential-integral equation of the single-layer homogeneous semi-infinite space model, regardless of the magnitude of the normalized active region depth variable, the predicted phase lag will not exceed a defined upper limit in radians, which approximates a quarter of pi. When the normalized depth coefficient is increased beyond this limit, the output feedback curve of the phase delay will exhibit oscillations or saturation decay. Therefore, when the actual measured phase delay extracted by the satellite observation system exceeds the upper limit of the mathematical extremum, the iterative algorithm within the single-layer homogeneous semi-infinite space model will face divergence failure or enter a multi-solution trap. This single-layer model branch can only be loaded and run when the boundary condition that the target project belongs to a shallow active region scenario is clearly determined.

[0100] In a further embodiment, the water and heat transfer model also includes a surface water conversion transfer function in series; the surface water conversion transfer function is a first-order low-pass filter containing the surface conversion time constant to be inverted, used to compensate for the physical lag in the process of actual rainfall being converted into effective surface water input.

[0101] Alternatively, a surface water transfer function is set in series before the water and heat transfer model to characterize the phase lag and amplitude attenuation in the process of converting the original rainfall input into effective surface water input; the surface water transfer function is a first-order low-pass filter, which contains the surface conversion time constant to be inverted.

[0102] In this embodiment, after the original rainfall in the actual meteorological environment falls onto the soil surface, it undergoes a series of meteorological-to-hydrological water resource transformation processes, including interception by surface vegetation, water filling of depressions, and overcoming the initial drying resistance of the surface soil. These transformation processes consume a specific amount of time, causing the effective water waveform that actually penetrates the surface interface and enters the soil pore network to participate in subsequent deep diffusion flow to attenuate and shift in time relative to the original rainfall waveform measured in the sky. If this preceding time-consuming process is ignored and model docking calculations are performed directly, the calculation engine will misallocate and superimpose the surface water retention time into the deep diffusion propagation time, thereby triggering a cumulative cascading error that overestimates the depth of the underlying active zone.

[0103] To correct boundary alignment errors, the computational engine forcibly incorporates an independent signal modification stage, namely the surface water transfer function (SFD), at the pre-data input of the main transfer model. The SFD is defined as a linear, time-invariant first-order low-pass filter operator structure at the system dynamics level. For example, the frequency domain network expression of the first-order low-pass filter is defined as:

[0104] H surface (ω)=1 / (1+i*ω*τ s );

[0105] Among them, H surface (ω) represents the complex frequency response result of the extracted surface water transformation transfer function, ω is the preset specific working angular frequency parameter corresponding to the extraction of periodic features, i is the imaginary unit used in complex operations, and τ s The introduced surface transformation time constant parameter is to be inverted.

[0106] The surface transformation time constant τ to be inverted s With a clear physical time dimension, such as the Gregorian year, it quantitatively characterizes the absorption and damping capacity of the target slope soil surface system for meteorological precipitation. For sparsely vegetated or highly porous weathered expansive soil slopes, the calculated time constant may fluctuate within several to more than ten days; for sites with dense surfaces or attached slope protection grid buffer layers, the constant value will broaden. The system opens the surface transformation time constant parameter to be inverted as part of the unknown inversion target matrix, and entrusts it to the subsequent joint optimization algorithm for adaptive convergence, thus eliminating the algorithmic uncertainty caused by subjectively given empirical reduction coefficients.

[0107] The overall phase delay of the theoretically predicted response characteristics is composed of the superposition of the phase delay of the surface water transformation transfer function and the phase delay of the fracture-enhanced two-layer model.

[0108] In this embodiment, the phase transfer superposition criterion under the cascaded architecture is standardized. When the multiphase medium transfer channels from the sky boundary into the surface and then into the deep bedrock are connected in series, for any selected single frequency point, the total time delay effect accumulated by the hydrological wave signal as it traverses the entire cascaded system is equal to the linear sum of the independent delay radians generated by each cascaded sub-module. For the introduced surface water conversion stage, the computational engine independently calls the complex phase extraction algorithm to obtain the first-level phase damping triggered by this module. The first-level phase delay is obtained by calculating the arctangent of the product of the angular frequency variable and the surface conversion time constant to be inverted.

[0109] Furthermore, the second-level physical diffusion phase delay calculated for the same frequency point by the fracture-enhanced two-layer model issued by the core engine is extracted. The calculation engine executes a scalar addition instruction to add the first-level infiltration phase delay scalar to the second-level deep diffusion phase delay scalar, synthesizing an array of theoretically predicted total phase delays containing full-path physical information. This synthesized array replaces the single delay array output by the original independent two-layer model and is encapsulated and imported into the subsequent target residual construction equation, where it is aligned end-to-end with the apparent total phase delay sequence intercepted by the multi-frequency satellite observation terminal.

[0110] In some optional implementations, regarding the control constraints on computing resources, if the capillary water absorption characteristics of the surface soil in the target area are extremely strong and there is no water accumulation barrier on the slope, as estimated in advance based on meteorological station data, resulting in the upper limit of the surface transformation time constant being constrained and assessed to be lower than the preset tolerance short-term threshold (e.g., less than five days), the computing system can trigger the rule bypass module. At this time, the system determines that the first-level phase delay injected by the front-end surface filtering link decays to a negligible high-order small perturbation range relative to the annual cycle span, and automatically bypasses the surface water transformation transfer function, directly backs back and outputs the theoretical calculation branch of the simple fracture-enhanced two-layer model, in order to reduce the optimization iteration overhead of the joint inversion process in the high-dimensional parameter space.

[0111] In one embodiment of this application, the periodic response features include phase delays extracted from at least two preset frequencies of vertical deformation time series and rainfall time series; substituting the periodic response features into a hydrothermal transfer model for inversion matching includes:

[0112] By using the phase delay at at least two preset frequencies as the observation end and the transfer function of the depth integral deformation response at different frequencies as the model end, a joint inversion equation set is constructed.

[0113] Specifically, after introducing the geometric parameters of the fracture layer, the theoretical transfer equation for a single frequency node transforms into an ill-conditioned underdetermined system containing two independent unknowns: the normalized fracture depth and the normalized thickness of the intact matrix layer. This means that a single equation cannot pinpoint the exact physical truth solution. To reconstruct the convergence constraints of the mathematical framework, the computational engine needs to inject more mutually orthogonal boundary information into the parameter system.

[0114] A joint inversion equation set is constructed by extracting the measured phase delay values ​​corresponding to at least two preset frequencies from the Fourier system in parallel, and loading them into the left side of the equation to form a multidimensional observation-side vector set. Correspondingly, multiple sets of depth integral deformation response transfer function calculation channels, equidistant from the number of extracted frequencies, are expanded in parallel on the right side of the equation to form the corresponding multidimensional model-side vector set. By forcibly constraining the equality of each sub-element on both sides of the equation, a simultaneous overdetermined equation structure is formed.

[0115] In a preferred implementation, when constructing the joint inversion equation set, cross-frequency constraints are constructed using the physical relationship between water penetration depths at different frequencies.

[0116] In this embodiment, if the unknowns in the constructed multidimensional equation system are allowed to vary freely and independently in different frequency channels, each additional observation frequency point will introduce a new batch of unknown derivative terms, making it impossible to close the solution loop. Cross-frequency constraints are a key mechanism for reducing the degree of freedom of parameters and improving the well-posedness of the equation system.

[0117] Since the rate at which water travels through the target porous matrix is ​​determined by the inherent pore network structure and physical conduction-diffusion rate of the specific soil layer, it does not vary or drift with changes in the driving frequency of the external boundary. Therefore, it can be determined that the absolute value of the physical diffusion coefficient of the medium remains constant across different frequency bands. Furthermore, for the same spatial observation pixel, the spatial distance from the bottom of the geological fracture interface to the surface, i.e., the total depth of the active zone, is also an objectively invariant physical quantity.

[0118] Set at least two preset frequencies, including an annual cycle frequency and a semi-annual cycle frequency. The expression for the cross-frequency constraint is:

[0119] ξ 2，2 =ξ 2，1 *2 1 / 2 ;

[0120] Where, ξ 2，1 ξ represents the normalized thickness of the intact matrix layer at the annual cycle frequency. 2，2 This represents the normalized thickness of the complete matrix layer at a six-month cycle frequency.

[0121] In the joint inversion equations, the normalized fracture depth is set to remain constant at different preset frequencies.

[0122] Specifically, according to the laws of partial differential diffusion theory, the specific depth at which the external environmental wave signal penetrates and attenuates to a preset boundary intensity is defined as the characteristic depth of water penetration. This characteristic physical quantity scale has an absolute inverse proportional relationship with the arithmetic square root of the excitation signal angular frequency parameter.

[0123] When the semi-annual cycle frequency is obtained and its angular frequency is quantized to be twice the annual cycle frequency, due to the constant physical diffusion coefficient law, the physical depth of water infiltration characteristics for the semi-annual high-frequency cycle will be forcibly compressed to twice the characteristic depth corresponding to the original annual cycle. 1 / 2 One-third. Meanwhile, the normalized thickness variable of the intact matrix layer is itself equal to the quotient of the actual physical thickness scalar of the intact layer divided by the corresponding water penetration characteristic depth scalar. Assuming the dividend, i.e., the physical thickness, remains constant, compressing and reducing the divisor will necessarily amplify the quotient by the same proportion.

[0124] From this, the expression for the cross-frequency constraint condition was deduced and assembled. Simultaneously, when registering variables in the background database, the computation engine pointed all the normalized fracture depth symbol pointers for the annual and semi-annual frequency channels to the same memory register address, thus performing a forced constant binding operation. Through physical locking without redundant degrees of freedom, the originally divergent multi-path equation system was reduced in dimension and converged into a well-posed equation system containing only two core unknown parameters.

[0125] Solve the joint inversion equations simultaneously to obtain the normalized fracture depth and the normalized thickness of the intact matrix layer.

[0126] like Figure 3 As shown, in a preferred embodiment, the inversion to obtain the normalized fracture depth and the normalized thickness of the intact matrix layer includes:

[0127] A uniform grid is constructed within a preset parameter space, and the sum of squares of the dual-frequency phase residuals between the phase predicted by the depth integral deformation response transfer function and the phase extracted from the observation end is calculated one by one at each grid point.

[0128] In this embodiment, a coarse mesh search is performed to prevent the subsequent gradient optimization module from getting trapped in local optima due to inappropriate random point placement. The upper and lower boundaries of the parameter space are defined manually based on geological knowledge. For example, the normalized fracture depth search coordinate axis boundary is specified to span the interval from 0.0 to 0.8, and the intact matrix layer thickness search coordinate axis boundary spans the interval from 0.1 to 10.0. The computation engine deploys a dense cross-test grid of points within this two-dimensional planar region according to a predetermined discrete scattering. At each discrete two-dimensional parameter coordinate intersection (r, ξ)... 2，1On the engine, the forward start model generates annual frequency prediction phase and semi-annual frequency prediction phase. The angular difference deviation components between these two predictions and the radar measured values ​​are calculated. Each of these deviation components is squared individually and then summed and integrated. The total cost function value corresponding to this discrete grid point is output and recorded, which is the sum of squares of the dual-frequency phase residuals.

[0129] The parameters corresponding to the grid point that minimizes the sum of squares of the dual-frequency phase residuals are selected to form the initial solution. Starting from the initial solution, the sum of squares of the dual-frequency phase residuals is nonlinearly optimized using the Levenberg-Marquardt algorithm to obtain the optimal normalized fracture depth and the normalized thickness of the intact matrix layer.

[0130] For example, the extracted dual-frequency phase delay is used as the input observation, and a global grid search is performed within a preset parameter space to locate the approximate minimum region of the sum of squared residuals. Using the parameters of this approximate minimum region as the initial solution, a Levenberg-Marquardt nonlinear optimization algorithm is initiated. After several iterations, it outputs the optimal combination of normalized fracture depth and normalized thickness of the intact matrix layer that minimizes the sum of squared residuals of the dual-frequency phase delay. When the change in the objective function residual between the two consecutive outputs is less than a preset tolerance, convergence is determined, and the iteration terminates.

[0131] In some optional implementations, a contingency plan is also provided for the degraded multi-frequency joint solution caused by the lack of specific meteorological stomata data. Specifically, target execution areas such as North China are often subject to a single-peak dominant precipitation cycle, and the energy of the semi-annual harmonic frequency bands extracted by waveform decomposition is extremely thin. If the calculation engine detects that the signal-to-noise ratio value fed back by the input terminal pixel at the semi-annual frequency point is less than the set reliable interception threshold requirement, for example, below the indicator scale of 3.0, the preset protection interruption mechanism is triggered. At this time, the execution of the dual-frequency equation overdetermined assembly and grouping step is immediately abandoned, and the adaptive degraded reconstruction is rolled back and loaded into the single-layer homogeneous model single-frequency measurement execution branch. On this degraded measurement branch, only the annual cycle phase single observation index stream that meets the quality access conditions is extracted and loaded, and single-frequency band search and comparison measurement is carried out for the only unknown variable scale, so that the system operation framework continues to be robust and coherent without overflowing errors under extreme conditions of data scarcity.

[0132] like Figure 4 As shown, in one possible embodiment, before calculating the total active area depth of each pixel, a step of reliable pixel screening using a preset threshold is further included, specifically including:

[0133] The frequency domain signal-to-noise ratio, inversion fitting residual, and linear deformation rate extracted from the vertical deformation time series of each pixel are obtained.

[0134] In other words, the frequency domain signal-to-noise ratio of each pixel obtained during the frequency domain decomposition process, the inversion fitting residual generated during the inversion matching process, and the linear deformation rate extracted from the vertical deformation time series by linear regression are obtained.

[0135] In this embodiment, after the inversion engine completes the core joint numerical optimization and before the spatial physical parameter restoration projection is initiated, the intervention system control flow acts as a data filtering barrier. The reason for constructing this filtering barrier is that the synthetic aperture radar full-frame image matrix, which actually covers tens or even hundreds of square kilometers, contains non-expansive soil areas, such as hard bedrock outcrops, artificial impermeable paving structures in urban areas, as well as catastrophic areas that do not conform to the assumed physical deformation laws and isolated areas of pure digital noise caused by the temporal decorrelation of vegetation. If all pixels in the entire domain are directly allowed to penetrate to the subsequent product calculation terminal without discrimination, a large number of physically meaningless deep numerical or complex errors will be generated.

[0136] Specifically, the obtained frequency domain signal-to-noise ratio parameter characterizes the intensity of suppression of the proportion of meteorological driving periodic signal energy in the target radar pixel reflection signal relative to the disordered random background noise. The obtained inversion fitting residual is used to characterize the degree of residual deviation between the theoretical phase delay and the measured phase delay during the inversion matching process; the obtained linear deformation rate is used to characterize the long-term unidirectional variation trend of the vertical deformation time series along the time direction after removing periodic fluctuations.

[0137] Determine whether each pixel simultaneously meets the pre-configured lower threshold for signal-to-noise ratio, upper threshold for residual, and upper threshold for linear deformation rate.

[0138] In this embodiment, the judgment logic is set as a parallel Boolean logic AND operation barrier, requiring that the parameters in the three dimensions must simultaneously pass through the separately calibrated tolerance boundary conditions.

[0139] The first verification standard is a pre-configured lower limit threshold for the signal-to-noise ratio (SNR) for the six-month cycle signal frequency band. Since the higher-order fluctuation characteristics of the six-month cycle are often weak and easily submerged in speckle noise, setting a lower limit threshold for the SNR, typically a pure scalar ratio greater than or equal to 3.0, ensures that the phase delay data source extracted to drive the constraint equations is based on a real, existing periodic rhythm, rather than a digital spurious root pieced together by software algorithms that erroneously force harmonic decomposition and fitting of white noise.

[0140] The second verification criterion is the upper limit threshold for the residuals. When the calculated inversion fitting residuals are too large, it indicates that the volumetric lesion process occurring at the current pixel space has become incompatible with the expansion and contraction porous medium model forcibly constructed by the core algorithm in terms of its physical origin. The upper limit threshold for the residuals needs to be defined and adapted based on the baseline of the basic data array's fault tolerance index given in the initial time-series inversion step. Based on the InSAR deformation measurement accuracy and periodic signal intensity of the study area, the statistical distribution of the inversion residuals across the entire pixel area can be analyzed to select the residual value corresponding to a preset percentile as the upper limit threshold for the residuals.

[0141] The third verification line is the upper limit threshold of the linear deformation rate. Setting an upper limit interception checkpoint can isolate and distinguish irreversible continuous collapse movements such as landslide slip surfaces. The vertical undulation characteristics caused by purely hydrothermal alternation activity zones have extremely high closed reciprocating elasticity. If the pre-regression test detects that the pixel exhibits extremely high-level continuous unidirectional subsidence or rise, exceeding the pre-configured limit rate benchmark zone, for example, in absolute values ​​of tens of millimeters per year, the system determines that the dominant deformation force of the micro-plot is no longer controlled by the hydrothermal rhythm diffusion, but is dominated and taken over by the sliding of the deep gravity fault surface, and therefore it is excluded from the scope.

[0142] Optionally, the judgment criteria also include embedded verification constraints on the physical rationality based on the analytical results of the core inversion structure itself. That is, to detect whether the calculated optimal two-layer model geometric parameters themselves cross the axiom red line of nature. For example, to check whether the mathematical feedback value of the normalized fracture depth has been overflowed by the nonlinear derivation engine to a physical inconsistency range greater than or equal to 1.0 and less than or equal to zero, and at the same time to judge whether the thickness variable of its complete matrix layer has been extremely compressed below the preset minimum spatial limit mark, such as the thickness scale coefficient being less than the limit value of 0.1.

[0143] Cells that do not meet any threshold condition are removed to obtain a set of reliable cells that have passed the screening. The total depth of the active area is calculated only for the cells in the set of reliable cells.

[0144] Furthermore, it is determined whether the normalized fracture depth obtained by inversion is within the preset physical range and whether the normalized thickness of the intact matrix layer is greater than the preset lower threshold, and pixels that do not meet the physical rationality conditions are removed.

[0145] Specifically, once any single pixel triggers any of the aforementioned boundary warning breakpoint instructions, the data flow channel to that pixel's memory address array will be directly cut off, assigned a null value, and removed, terminating all subsequent multiply-accumulate extension calculations. After spatial traversal and cleaning, the index number list of all continuous or discrete patch data remaining and delineated on the system matrix bridge is integrated, compiled, and packaged, which is considered the reliable pixel set output that has passed the final security verification. The subsequent physical final length dimension restoration decoding equation, which incorporates the pre-calibrated moisture diffusion coefficient scalar constant, is only used for point-to-point precise calculation of the specific spatial nodes registered within the selected reliable pixel set. For pixels that fail the screening, null or invalid values ​​can be assigned to the output active area depth distribution map, or they can be removed using a mask, to avoid misusing pixels that do not meet the physical model assumptions or have insufficient signal quality for subsequent spatial analysis, engineering zoning, and statistical evaluation.

[0146] In some alternative implementations, the operational points for quality control and reliable pixel screening and interception processes can also be designed to be moved forward. After completing the initial preprocessing and unwrapping of the global temporal image to obtain the coarse-grained vertical deformation backbone sequence, an early truncation and cleaning task is immediately initiated before the periodic response features are extracted and analyzed. This involves calculating and statistically refining the linear deformation rate index and the ratio of high-frequency divergence noise interference for all pixels. By reducing the burden on nodes of the spatial screen in advance, the cost of the huge amount of computing power wasted in the subsequent joint iterative search in the high-dimensional complex domain nonlinear joint inversion model can be reduced.

[0147] According to one aspect of this application, the pre-calibrated moisture diffusion coefficient is obtained through the following pre-construction steps:

[0148] Obtain the measured active zone depth and measured fracture depth at a given calibration profile.

[0149] In this embodiment, the given calibration profile refers to a geological exploration point pre-established by manual or mechanical means through excavation and drilling within the target slope area. The calibration profile is preferably located in a position with representative expansive soil distribution, complete monitoring data, and fracture development characteristics consistent with the main conditions of the study area. To ensure the reliability of the calibration results, profiles with clear periodic deformation signals, insignificant long-term slippage, and where reference values ​​for the active zone depth and fracture depth can be obtained through on-site observation are preferred as calibration profiles. The measured active zone depth refers to the physical scale distance from the deepest boundary of significant seasonal fluctuations in the stratum's water content to the surface, determined by a sequence of water content sensor probes deployed year-round on-site or by periodic sampling and indoor drying comparison. The measured fracture depth refers to the depth of the bottom cutoff surface, characterized by a preferential water-conducting network of shrinkage-induced cracks, directly defined through on-site surveying using pit sketching and logging, water injection tests, or ground-penetrating radar scanning.

[0150] Obtaining the measured depth of active zones and measured fracture depths can break through the fictitious space generated by the preceding theoretical calculations, which only contain dimensionless structural parameters. Typically, such in-situ geophysical exploration operations are costly and damaging to the slope structure, making it impossible to deploy high-density grids on vast mountain slopes. Therefore, the obtained attributes of measured active zone depths and measured fracture depths naturally belong to extremely scarce and discrete point sampling samples.

[0151] Extract the periodic response features at the calibration profile.

[0152] Specifically, the extraction process involves mapping the latitude and longitude plane coordinates of a given calibration profile to the pixel matrix raster coordinate system generated by the synthetic aperture radar processing system using a geodetic coordinate system transformation function. This locates the specific remote sensing radar pixel patch covering the borehole's actual coordinates. Then, the dedicated response dataset generated after the spatiotemporal signal decomposition of this specific spatial pixel unit is extracted. This clearly reveals the inherent measured phase delay radian parameter set of this reference point under a preset working period, i.e., the periodic response feature scalar package.

[0153] The measured active zone depth, measured fracture depth, and periodic response characteristics at the calibration profile are substituted into the phase equation corresponding to the hydrothermal transfer model for inverse deduction, and the physical diffusion parameters corresponding to the calibration profile are calculated.

[0154] In this embodiment, the reverse inversion can logically reverse and interchange the known and unknown values ​​in the forward inversion process. In the conventional forward mode, the system assumes the known soil seepage resistance constant to find the unknown stratum interface boundary. In the offline calibration mode, the system substitutes the explicitly obtained stratum geometric thickness into the inverted structure of the transcendental differential equation to obtain the theoretical diffusivity value of the medium that makes the geometric field and the observation lag time self-consistent.

[0155] In a preferred implementation, the physical diffusion parameters corresponding to the calibration profile are calculated and solved analytically using the following formula:

[0156] α m =(ω1*H m，cal 2 ) / (2*[ξ 2，1 (x cal )] 2 );

[0157] Where, α m The physical diffusion parameters of the complete matrix layer are represented; ω1 represents the annual periodic angular frequency used when extracting periodic response features; H m，cal ξ represents the thickness of the intact matrix layer at the calibration profile, obtained by subtracting the measured active zone depth from the measured fracture depth; 2，1 (x cal The normalized thickness at the annual periodic frequency is obtained by substituting the periodic response characteristics at the calibration profile into the depth integral deformation response transfer function.

[0158] The final calibrated regional material flow conduction coefficient is proportional to the square of the actual span thickness of the target intact dense soil layer, and inversely proportional to the square of the normalized attenuation scaling factor determined by the feedback of surface remote sensing waveforms.

[0159] The physical diffusion parameters are used as pre-calibrated moisture diffusion coefficients required for inversion calculations for each pixel.

[0160] Specifically, the constant values ​​of individual physical diffusion parameters obtained through reverse engineering from a single or limited number of test benchmark boreholes are extracted and encapsulated. These encapsulated variables are then distributed to the main control depth inversion module program responsible for performing depth calculations across the full-frame image radar map matrix. During the decoding process—which involves calculating and reconstructing the specific physical meter-level depth from an abstract, dimensionless state for thousands of blind pixels on the map—a pre-calibrated moisture diffusion coefficient constant is unconditionally invoked and embedded as a unified physical multiplier benchmark.

[0161] In some optional implementations, a multi-point collaborative median smoothing strategy is provided to counteract the extremely heterogeneous spatial variability of mineral composition and compaction density in natural geological bodies. When the system obtains three or more sets of borehole data from given calibration profiles at different discrete locations across a large area of ​​slope through various means, the system drives a loop control program to execute and complete the entire set of reverse extraction steps for each independent borehole. The system summarizes and extracts an array list containing physical diffusion parameters of multiple individual calibration boreholes with numerical drift between them. Instead of using the conventional arithmetic mean, the array sorting function is called to rearrange the set sequence into a stepped array according to numerical size, and the median value at the center of the sorted array is extracted as the pre-calibrated regional generalized water diffusion coefficient constant that ultimately represents the entire large-scale analysis area with the greatest generalization and compatibility. By employing a nonlinear selection strategy, the destructive and pulling effects of extreme deviations in gross parameter errors caused by the presence of abnormally connected seepage funnels or extremely dense clay clumps in a particular isolated plot on the global parameter representative baseline system were eliminated.

[0162] In one embodiment of this application, the method further includes a step of verifying the self-consistency of the hydrothermal transfer model based on amplitude consistency, specifically including:

[0163] Obtain the pre-configured rainfall-suction conversion coefficient and the swelling-shrinkage coefficient of expansive soil.

[0164] In this embodiment, rainfall is a hydrological flux measured in millimeters. The direct physical cause driving the shrinkage and expansion of the expansive soil particle skeleton due to volumetric strain is the change in matrix suction at the water-air interface of the microscopic meniscus within the soil pores, measured in kilopascals (kPa). Under natural soil physics, rainfall in millimeters cannot be directly equated to the kilopascal suction scalar. To bridge the gap in magnitude and unit between these two heterogeneous physical fields, a pre-configured rainfall-suction conversion coefficient is required. A linear approximate equivalent conversion ratio is defined for the effective atmospheric water input per millimeter at a specific soil surface, which is converted into the corresponding average matrix suction change in kilopascals. The composite derivation dimension is set to kilopascals per millimeter.

[0165] Furthermore, the swelling and shrinkage coefficient of expansive soil is an empirical constant mechanical parameter reflecting the strength of the inherent hydrophilic swelling potential of the minerals within that specific soil type. It represents the absolute variable in length of the vertical surface deformation caused by a change in the unit matrix suction within the soil profile at a unit depth. Its physical dimensions and composite units are expressed in millimeters per kilopascal per meter, and it is typically obtained through standardized indoor uniaxial consolidation apparatus compression-swelling comparative geotechnical tests.

[0166] The theoretically predicted deformation amplitude is obtained by calculating the observed rainfall amplitude, rainfall-suction conversion coefficient, expansion and contraction coefficient, total depth of the active zone, and theoretically predicted amplitude response output by the hydrothermal transfer model, extracted from the rainfall time series.

[0167] Specifically, the observed rainfall amplitude extracted from the rainfall time series is the peak extremum parameter of the rainfall wave driving source energy obtained through Fourier band hull decomposition, which retains the original millimeter meteorological dimension. The theoretically predicted amplitude response output by the water and heat transfer model is the purely dimensionless proportional attenuation coefficient output after the system's complex two-layer transfer matrix operation. The total depth of the active area is the thickness scale scalar carrying the meter dimension output by the model's depth principal inversion task.

[0168] In a preferred embodiment, the theoretically predicted deformation amplitude is obtained, and dimensional matching and physical transformation are performed using the following product formula:

[0169] A model =β s *κ*A P *H*|Φ bi |;

[0170] Among them, A model Indicates the theoretically predicted deformation amplitude; β s Indicates the expansion / contraction coefficient; κ represents the rainfall-suction conversion coefficient used to convert rainfall units to suction units; A P Indicates the observed rainfall amplitude; H represents the total depth of the active area; |Φ bi | represents the theoretically predicted amplitude response obtained by taking the absolute value of the depth integral deformation response transfer function.

[0171] In this embodiment, dimensional consistency is achieved on both sides of the equation. During the parameter multiplication derivation, the kilopascal dimension of the denominator of the expansion / contraction coefficient and the kilopascal dimension of the numerator of the rainfall-suction conversion coefficient cancel each other out during the multiplication elimination cross-comparison; simultaneously, the meter-level dimension of the denominator of the expansion / contraction coefficient and the meter-level spatial length characteristic dimension of the total depth of the active area are also multiplied and canceled out; the millimeter dimension introduced in the denominator of the rainfall-suction conversion coefficient and the inherent millimeter dimension carried by the observed rainfall amplitude are multiplied and canceled out. After the full-chain dimensional multiplication cross-elimination and clearing, the final residual dimensional attribute of the complex aggregation factor cluster on the right side of the equation is millimeters. The derived dimensions and the expected ground undulation displacement amplitude physical length index to be generated on the left side of the equation achieve overlap and balance in a mathematical dimension sense.

[0172] The theoretically predicted deformation amplitude is compared with the observed deformation amplitude extracted from the vertical deformation time series, and the model self-consistency evaluation results are output.

[0173] Specifically, the comparison step is a model health status review bypass mechanism designed based on the dual-frequency cross-isolation verification principle. In all model building and splicing, as well as parameter nonlinear deep optimization search operation loops, the system only extracts and utilizes the phase delay characteristic parameters in each frequency component, while deliberately shelving and sealing the amplitude energy intensity characteristics of all waveforms without intervention or retrieval.

[0174] The observed deformation amplitude extracted from the vertical deformation time series is the measured value of the physical maximum offset scale of surface waves actually captured and recorded by the radar terminal. The system calculates the absolute difference between the theoretically predicted deformation amplitude and the observed deformation amplitude. If this difference is within the set tolerance confidence interval, the output model self-consistency evaluation result is certified, confirming that the deep geometric parameter network architecture derived by the system not only closely matches the time delay rhythm of wave propagation, but also derives the correct deformation peak extreme value in the energy conservation attenuation dimension. Otherwise, an alarm log is triggered, and an abnormal calibration code is thrown.

[0175] In some optional implementations, for field-data-poor mapping scenarios where it is difficult to obtain accurate absolute values ​​of the rainfall-suction conversion coefficient due to the lack of high-precision micro-tensimeter matrix arrays deployed in real-time on-site, a weakened, downgraded qualitative trend spatial comparison strategy is provided. In this case, the system abandons the calculation of rigid values ​​with clear absolute dimensions and millimeter precision. Instead, it treats the conversion coefficient components and expansion / contraction coefficient terms as a block of unknown common constant coefficients uniformly distributed across the entire field. The system then calculates the gradient divergence map of the spatial distribution of relative fluctuation amplitudes remaining at all pixel grid points after deducting this constant block. This map is then compared with the measured thermal relative weight distribution cloud map of the corresponding pixel model array extracted from actual radar scans using Spearman's spatial correlation topological fitting statistical scoring. This provides a rating feedback index to assess the reasonableness and reliability of the model's spatial relative trend distribution array distribution.

[0176] This embodiment introduces a cross-physical field connection transformation constant, which allows the signal amplitude characteristics that are not constrained and invoked in the main process of deep inversion to be included in the closed loop for verifying the rationality of the constructed water and heat transfer model.

[0177] According to one aspect of this application, a system for identifying hydrothermal control activity zones on expansive soil slopes includes:

[0178] The data acquisition module is used to acquire the vertical deformation time series of each pixel within the target slope area and the rainfall time series within the corresponding time period.

[0179] Specifically, the data acquisition module is manifested in the system's physical architecture as a cluster of data communication and cached input interfaces. It includes application programming interfaces (APIs) for establishing network connections with satellite ground receiving stations or cloud-based remote sensing databases, and a dedicated protocol stack for parsing meteorological reports shared by meteorological bureaus. The acquired data arrays are pushed into an independent memory heap area allocated to the data acquisition module for data cleaning, removing invalid frames or missing transition points. At the system hardware level, the data acquisition module relies on a front-end data server pool configured with a high-speed network interface controller and a large-capacity hard disk array to ensure the throughput and resident capacity of multi-source, massive spatiotemporal sequence data.

[0180] The feature extraction module is used to process vertical deformation time series and rainfall time series to extract periodic response features that characterize the moisture-driven deformation process.

[0181] In this embodiment, the feature extraction module is essentially a digital signal processing engine residing within the system's core processing unit. This engine incorporates a Fast Fourier Transform operator package and a matrix least squares fitting function library. Upon being awakened by the front-end data stream, the engine directly takes over the original waveform matrix in the memory heap, calls the underlying mathematical operation core to perform detrending operations and multiharmonic frequency decomposition, calculates the phase angle and magnitude at the specified target frequency, and uses built-in instructions to generate phase difference vectors and amplitude ratio vectors. The feature extraction module heavily relies on the computer's floating-point arithmetic unit and multi-threaded concurrent processing mechanism to handle the massive parallel computing overhead covering hundreds of thousands of pixels.

[0182] The model building module is used to construct a hydrothermal transfer model that reflects the water diffusion mechanism and depth integral deformation mechanism of expansive soil. The hydrothermal transfer model establishes a physical mapping relationship between the theoretically predicted response characteristics and the normalized active zone depth.

[0183] Specifically, the model building module is the internal logic and rule definition center of the system. It does not process specific numerical values, but rather is responsible for instantiating the physical architecture boundaries in memory. It embeds components supporting complex number operations and, by declaring constant constraints, setting the zero-flux boundary function stack for differential equations, and loading transcendental function pointers, solidifies and reconstructs the natural diffusion laws of the external environment into a set of pure machine code functions. In actual operation, all subsequent forward computations and derivations use the virtual mathematical space instance registered and mapped by the model building module as the host runtime environment.

[0184] The depth inversion module is used to substitute the periodic response features into the water and heat transfer model for inversion matching, solve for the normalized active area depth, and calculate the total active area depth of each pixel by combining the pre-calibrated water diffusion coefficient.

[0185] In this embodiment, the depth inversion module is the terminal execution carrier for the entire system's computing power consumption. It deploys a nonlinear optimizer, such as a Gauss-Newton or Levenberg-Marquardt optimization closed-loop program. Simultaneously, the depth inversion module extracts measured feature vectors from the feature extraction module and extracts virtual mapping function structures from the model building module, initiating thousands of iterative iterations under parameter tolerance constraints. Its hardware support environment is typically built on a multi-core CPU architecture; for large-scale parameter spaces, it can also call the unified computing device architecture interface of the graphics processor to perform tensor-level concurrent acceleration. After calculating the optimal dimensionless parameters, the depth inversion module reads and extracts the calibrated diffusion constant from the system's solid-state storage medium or configuration file, performs the final dimensional multiplication decoding conversion, and pushes the calculated spatial absolute depth distribution field map to the display rendering port or generates a geological engineering file.

[0186] According to another aspect of this application, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the steps of the method for identifying hydrothermal control activity zones of expansive soil slopes as described in any of the above embodiments.

[0187] In this embodiment, the computer-readable storage medium defines the physical solidified form of the system's underlying implementation. The computer-readable storage medium specifically includes, but is not limited to, non-volatile solid-state drives, static random access memory arrays, optical disc drives, or virtual object storage volumes distributed in cloud platform data centers. The computer program code set burned and resided on it includes all the linear formula streams, discrimination and verification barriers, parameter degradation bypasses, and data interaction communication logic disclosed in detail in the preceding embodiments. When the central processing unit of the electronic computing device equipped with this medium reads, compiles, and distributes these binary or high-level language instruction sets according to a predetermined timing sequence, the electronic computing device is physically reshaped and transformed into a customized analysis host with dedicated functions for quantitative mapping of expansive soil active areas.

[0188] In a preferred detailed embodiment, a method for identifying hydrothermal control activity zones on expansive soil slopes specifically includes the following process:

[0189] Obtain the time series of vertical deformation of each pixel within the target slope area and the time series of rainfall within the corresponding time period.

[0190] In this embodiment, expansive soil slope areas with stable seasonal rainfall patterns are preferred. The vertical deformation time series of each pixel is obtained using multi-temporal SAR imagery through time-series interferometry, while daily rainfall data corresponding to the SAR observation time window are also acquired. To ensure consistency in the frequency domain analysis, the vertical deformation time series and rainfall time series are aligned to a unified time step and resampled.

[0191] Frequency domain decomposition was performed on the time series of vertical deformation and rainfall to extract periodic response features at at least two preset frequencies.

[0192] In this embodiment, the preset frequencies preferably include an annual cycle frequency and a semi-annual cycle frequency. The deformation phase and amplitude of the vertical deformation time series under the annual and semi-annual cycles are extracted through frequency domain decomposition, and the rainfall phase and amplitude of the rainfall time series under the annual and semi-annual cycles are also extracted. The annual cycle phase delay and the semi-annual cycle phase delay are calculated based on the difference between the deformation phase and the rainfall phase at the same frequency; the annual cycle amplitude ratio and the semi-annual cycle amplitude ratio are calculated based on the ratio of the deformation amplitude to the rainfall amplitude at the same frequency. The annual cycle phase delay and the semi-annual cycle phase delay are used as the observation inputs for subsequent dual-frequency joint inversion.

[0193] Construct the depth integral deformation response transfer function corresponding to the fracture-enhanced two-layer model.

[0194] In this embodiment, the active zone is divided into a fractured layer and a intact matrix layer along the depth direction. The fractured layer characterizes the preferential infiltration channels formed by the drying shrinkage fractures on the surface of the expansive soil, while the intact matrix layer characterizes the diffusion-dominated water migration region below the bottom of the fractured layer. Based on the approximately in-phase response within the fractured layer, the diffusion phase shift response within the intact matrix layer, and the zero flux boundary condition at the bottom of the active zone, a fracture-enhanced bilayer transfer function Φ is constructed. bi The transfer function is a complex-valued function, with its argument corresponding to the theoretically predicted phase delay and its magnitude corresponding to the theoretically predicted amplitude response. A physical mapping relationship is established between the observed phase delay and the normalized fracture depth r and the normalized thickness ξ2 of the intact matrix layer through the transfer function.

[0195] A dual-frequency joint inversion equation set was constructed and the normalized fracture depth and the normalized thickness of the intact matrix layer were solved.

[0196] In this embodiment, the annual and semi-annual phase delays are used as the observation end, and the theoretical phase delay output by the fracture-enhanced bilayer transfer function at the corresponding frequency is used as the model end to construct a dual-frequency joint inversion equation set. Simultaneously, based on the relationship between the characteristic depth of water infiltration at different frequencies, cross-frequency constraints are introduced to ensure that the normalized thickness of the intact matrix layer at the semi-annual frequency and the normalized thickness of the intact matrix layer at the annual frequency satisfy a fixed proportional relationship, and the normalized fracture depth is set to remain consistent across different frequencies.

[0197] In an exemplary scenario, if the annual phase delay extracted from the observation end is 0.50 radians and the semi-annual phase delay is 0.75 radians, a uniform grid can be constructed within a preset parameter space. The sum of squared residuals for the dual-frequency phase is calculated for each grid point, and the parameter set with the smallest sum of squared residuals is selected as the initial solution. Starting from the initial solution, the Levenberg-Marquardt nonlinear optimization algorithm is used to continue searching for the optimal solution. After iteration, the normalized fracture depth r is approximately 0.21, and the normalized thickness ξ of the intact matrix layer under the annual cycle is obtained. 2，1 The optimal parameter set is approximately 0.98.

[0198] Physical diffusion parameters of the intact matrix layer are determined based on known geometric information from the calibration profile.

[0199] Specifically, one or more calibration profiles are selected. For each calibration profile, the measured value H of the total depth of its active area is obtained. cal and the measured value of crack depth H c，cal And calculate the thickness H of the intact matrix layer at the calibration profile. m，cal =H cal -H c，cal Simultaneously, the annual periodic phase lag of the pixel containing the calibration profile is extracted, and the corresponding ξ at that profile is obtained by inverting the transfer function. 2，1 Based on the thickness H of the intact matrix layer at the calibration profile. m，cal , annual periodic angular frequency ω1 and ξ 2，1 The physical diffusion parameter α of the complete matrix layer was obtained by reverse calculation. m .

[0200] In scenarios with multiple calibration profiles, multiple α values ​​can be calculated separately. m The values ​​are obtained by selecting the median or by least squares statistical estimation to obtain pre-calibrated diffusion parameters that are representative of the region, so as to reduce the bias caused by local abnormal profiles.

[0201] Restore the total depth of the activity area and form a spatial distribution result.

[0202] In this embodiment, dual-frequency joint inversion is performed on all pixels that pass the reliability screening to obtain the normalized fracture depth r and the normalized thickness ξ of the intact matrix layer for each pixel. 2，1 Combined with the pre-calibrated α m And the scale conversion relationship at the corresponding frequency, to restore the complete matrix layer thickness H of each pixel. m Based on the geometric relationship between the fracture layer thickness and the intact matrix layer thickness, the total depth H of the active region of each pixel was calculated. Among them, the intact matrix layer thickness H is [value missing] at the annual cycle frequency. m Normalized thickness ξ of the intact matrix layer 2，1The water penetration characteristic depth is obtained by multiplying it by the water penetration characteristic depth at the corresponding frequency; the water penetration characteristic depth is obtained by the pre-calibrated physical diffusion parameter α of the intact matrix layer. m The annual periodic angular frequency ω1 is used to determine this. Specifically, the process for recovering the total depth of the active region of each pixel is as follows: based on the pre-calibrated physical diffusion parameters α of the intact matrix layer... m Calculate the characteristic depth d of water infiltration using the annual cycle angular frequency ω1. p :

[0203] d p =sqrt(2·α m / ω1);

[0204] The normalized thickness ξ of the complete matrix layer obtained for each pixel during inversion is calculated. 2，1 With the characteristic depth d of water penetration p Multiplying these yields the physical thickness H of the complete matrix layer. m :

[0205] H m =ξ 2，1 ·d p ;

[0206] According to the definition of normalized fracture depth r, r = H c / H and the total depth of the activity area H=H c +H m Based on the geometric relationships, the total depth H of the active region of each pixel is calculated:

[0207] H=H m / (1−r);

[0208] Where H c The thickness of the fracture layer is given. This leads to a spatial distribution map of the total depth of the active zone within the study area.

[0209] For example, the total depth of the output active zone shows a certain spatial distribution difference, with a larger active zone depth corresponding to the slope toe area and a smaller active zone depth corresponding to the slope shoulder area. Compared with field monitoring, the inversion results of the two-layer model better reflect the actual engineering characteristics of the coexistence of preferential flow in surface fissures and diffusion in deep matrix in expansive soil than the rapid preliminary estimation results of single-layer equivalent diffusion.

[0210] It should be noted that the specific values ​​in the above embodiments are only used to illustrate the solution process and parameter update method of the present invention, and do not constitute a limitation on the scope of protection of the present invention.

[0211] Furthermore, when the signal-to-noise ratio of the signal corresponding to the semi-annual cycle frequency of the target area is lower than the preset threshold, the dual-frequency constraint features cannot be stably extracted, or the prior information related to the fracture layer parameters cannot be obtained, the system can automatically switch to the simplified implementation method to perform a rapid preliminary estimation of the active area depth.

[0212] In another simplified implementation, for scenarios where surface fissures are not significantly developed, the active zone is shallow, only a preliminary estimate of the active zone depth is needed quickly, or it is difficult to obtain the multi-frequency constraint information required by the aforementioned fissure-enhanced two-layer model, a time-domain response identification method based on the single-layer equivalent diffusion assumption can be used to make a preliminary estimate of the depth of the hydrothermal control active zone of expansive soil slopes, specifically including:

[0213] Multi-temporal synthetic aperture radar remote sensing images of the study area were acquired, and the slope deformation time series was obtained through time series interferometry.

[0214] The deformation signal is decomposed into components to extract the hydrothermal driven deformation signal;

[0215] A hydrothermal driven response model was established by combining rainfall and temperature time series to identify the deformation response lag time caused by rainfall infiltration.

[0216] Under the assumption of single-layer equivalent diffusion, an approximate mapping relationship between deformation response time delay and moisture disturbance propagation depth is established based on response lag time and soil moisture diffusion characteristics, and the depth of hydrothermal control zone of expansive soil slope is preliminarily estimated.

[0217] The spatial distribution of depth in the active area is obtained by calculating remote sensing pixels.

[0218] In this simplified embodiment, the depth of the hydrothermal control zone of expansive soil slopes is preliminarily estimated by utilizing the deformation response lag time caused by rainfall infiltration in remotely sensed deformation time series. The deformation component decomposition step includes extracting the long-term deformation trend using a trend fitting method and using the remaining deformation signal as the hydrothermal driven elastic deformation component. The hydrothermal driven response model establishes the deformation driving relationship through the rainfall infiltration process and the temperature change process. The deformation response lag time caused by rainfall infiltration and the deformation response lag time caused by temperature change are determined by a joint search method.

[0219] It should be noted that the diffusion model used in this simplified embodiment is a single-layer equivalent diffusion model. It does not explicitly distinguish the stratified coupling effect between preferential flow through surface fissures and diffusion into the deep intact matrix of expansive soil. Instead, it uniformly characterizes water transport capacity through equivalent diffusion parameters. Therefore, the results obtained in this embodiment are approximate estimates of the active zone depth, suitable for rapid screening, initial value setting, or comparison and verification with the results of the preferred embodiment.

[0220] Furthermore, calibrating the equivalent diffusion parameters using on-site monitoring data can improve the reliability of the initial depth estimation results for the active area. By performing pixel-by-pixel calculations on the remote sensing pixels of the study area, the spatial distribution of the active area depth can be obtained, and regional partitioning can be performed based on the depth size of the active area. The deformation response lag time is identified using time series analysis methods, including one or more of correlation analysis, spectral analysis, machine learning, or deep learning methods. The deformation response lag time is obtained by identifying the time response relationship between the remote sensing deformation time series and the rainfall time series.

[0221] According to another aspect of this application, in a complete implementation process, a method for identifying hydrothermal control activity zones of expansive soil slopes can be executed in the following order:

[0222] Step A1: Acquire multi-temporal SAR images and rainfall time series of the target area, and complete time alignment and necessary resampling processing;

[0223] Step A2: Obtain the time series of vertical deformation for each pixel through time-series InSAR processing;

[0224] Step A3: Detrend the vertical deformation time series to separate the hydrothermal driven elastic deformation signal;

[0225] Step A4: Extract periodic response features from the vertical deformation time series and rainfall time series, preferably including phase delay and amplitude ratio;

[0226] Step A5: Construct a water and heat transfer model and calculate the theoretically predicted response characteristics under a given set of parameters;

[0227] Step A6: Parameter inversion is achieved by minimizing the residual function to obtain the normalized fracture depth, the normalized thickness of the intact matrix layer, and / or the normalized active zone depth.

[0228] Step A7: Calculate the total depth of the active area by combining the pre-calibrated diffusion coefficient and scale conversion relationship;

[0229] Step A8: Perform reliable pixel screening on the results and output the active area depth distribution map and related evaluation results.

[0230] In summary, this application does not follow the approach of directly equating the surface deformation response lag with the water propagation time. Instead, it establishes a depth integral deformation response transfer function in the expansive soil fracture-matrix bilayer medium, converting the observed phase delay into an invertible active zone depth parameter. Furthermore, it utilizes multi-frequency periodic response characteristics to construct joint constraints, thereby achieving a stable solution for the active zone depth. Specifically, by constructing a hydrothermal transfer model under the coupled interaction of the expansive soil fracture layer and the intact matrix layer, the active zone depth identification problem is elevated from empirical statistical time delay matching to a physical inversion problem with clear boundary conditions and medium structure constraints. By establishing a depth integral deformation response transfer function, an analytical mapping relationship is established between the periodic deformation phase delay observed on the surface and the underground active zone depth, avoiding the misinterpretation of the total surface deformation as a local propagation response at a single depth. By introducing periodic response characteristics at at least two preset frequencies to construct a joint inversion equation set, a synergistic constraint is achieved on the normalized fracture depth and the normalized thickness of the intact matrix layer, improving the stability and physical interpretability of the active zone depth inversion results.

[0231] In a detailed embodiment, a certain expansive soil canal slope area is taken as the research object. Satellite remote sensing monitoring data and meteorological data are used to identify and analyze the hydrothermal control activity zone of the slope. The slope is a typical deep excavation expansive soil canal slope, with a slope height of about 20m, a slope angle of about 25° to 30°, and the soil type is weak to moderate expansive soil. The average annual rainfall in the area is about 800mm.

[0232] Multi-temporal synthetic aperture radar (SAR) imagery data of the study area was acquired, and remote sensing image sequences covering multiple seasonal cycles were selected. In this embodiment, Sentinel-1 satellite SAR imagery data was used: the time span was 2018-2023; the number of images was 120; the spatial resolution was approximately 20m; digital elevation model (DEM) data of the study area was also acquired to eliminate the influence of topographic phase. In addition, meteorological data of the study area were acquired, including time series of rainfall and temperature, and the acquisition time of the remote sensing images was aligned with the time of the meteorological data.

[0233] Multi-temporal SAR images were processed using time-series synthetic aperture radar interferometry (InSAR). Fine image registration was performed, and small baseline interferometric pairs were established. Interferograms were generated, and phase filtering, phase unwrapping, and terrain phase elimination were applied to obtain the deformation time series d(t) for each pixel.

[0234] Since slope deformation includes both long-term slip deformation and hydrothermal driven deformation, it is necessary to decompose the deformation signal. The total deformation is expressed as:

[0235] d(t)=d pla (t)+d ela (t);

[0236] Where d(t) is the deformation value of the pixel at time t; d pla (t) represents the long-term plastic deformation component; d ela (t) represents the hydrothermal driven elastic expansion and contraction component.

[0237] Long-term deformations are represented by polynomial functions:

[0238] d pla (t) = vt + at 2 +bt 3 ;

[0239] Where v is the long-term deformation rate, and a and b are higher-order deformation coefficients.

[0240] The long-term deformation component was obtained by fitting using the least squares method, and this component was then removed from the total deformation to obtain the hydrothermal driven deformation signal d. ela (t).

[0241] The swelling and shrinkage deformation of expansive soil is mainly affected by rainfall infiltration and evaporation. Therefore, a hydrothermal driven model is established:

[0242] d ela (t)= α1P e (t-τ1)+ α2T m (t-τ2)+c;

[0243] Where P e (t) represents the effective rainfall, T m (t) represents the average temperature, τ1 represents the rainfall response lag time, τ2 represents the temperature response lag time; α1, α2, and c are model parameters.

[0244] The formula for calculating effective rainfall is:

[0245] P e (t)=∑ i=0 k P(ti)e -λi ;

[0246] Where P(ti) is the rainfall on day ti, λ is the rainfall attenuation coefficient, i is the day index, k is the maximum number of backtracking days, and e is the natural constant.

[0247] The response lag times of rainfall and temperature were identified using a joint search algorithm. The lag time ranges were defined as: τ1∈[0, 30], τ2∈[0, 15]. The model residuals were calculated under different lag combinations.

[0248] SSE=∑(d obs -d model ) 2 ;

[0249] Where SSE is the sum of squared residuals, ∑ represents the summation range covering all observation times, and d obs d represents the measured deformation value. model These are the model's predicted values.

[0250] Select the lag combination (τ1*, τ2*) with the smallest residual, where τ1* is the rainfall response lag time with the smallest residual and τ2* is the temperature response lag time with the smallest residual.

[0251] Moisture changes in expansive soil can be approximated as a one-dimensional equivalent diffusion process. It should be noted that this equivalent diffusion process is a simplified representation of the combined effect of rapid transport through surface fissures and slow diffusion through the underlying matrix. It is suitable for scenarios where fissures are not significant or only a rapid initial estimation is required, and is not intended to replace the fissure-enhanced two-layer model described in the aforementioned implementation. Under the assumption of single-layer equivalent diffusion, the governing equation for moisture migration can be expressed as:

[0252] Ψγ / Ψt=α·(Ψ 2 γ / Ψz 2 );

[0253] Where Ψ is the partial derivative; γ is the volumetric water content; t represents time; z represents the vertical depth coordinate; and α represents the equivalent diffusion coefficient.

[0254] Based on the analytical solution of the diffusion equation, an approximate relationship between the propagation depth and propagation time of moisture disturbance can be obtained:

[0255] L = 2·sqrt(α·Δt);

[0256] Where L is the depth of moisture disturbance propagation, and Δt is the propagation time.

[0257] The rainfall response lag time obtained through time-domain identification is used as the equivalent propagation time, i.e.:

[0258] Δt=τ1*;

[0259] This yields a preliminary estimate of the depth of the active zone on the expansive soil slope:

[0260] H = 2·sqrt(α·τ1*);

[0261] Where H is the approximate value of the active zone depth calculated under the assumption of single-layer equivalent diffusion. Compared with the inversion results of the fracture-enhanced two-layer model, this result does not explicitly consider the layered coupling effect between the surface fracture layer and the deep intact matrix layer, and is therefore more suitable as a rapid estimate, initial parameter value, or comparative verification value.

[0262] Furthermore, due to differences in mineral composition, fissure development, and soil compaction in expansive soils from different regions, it is necessary to calibrate the equivalent diffusion coefficient in this embodiment. If a monitoring profile can obtain the measured depth z of the active area... ref and the corresponding response lag time τ ref The equivalent diffusion coefficient can then be obtained using the following formula: α=z ref 2 / 4τ ref It should be noted that the α obtained here is a comprehensive parameter under the assumption of single-layer equivalent diffusion. Its physical meaning is to provide a unified approximate characterization of the surface fracture conduction effect and the underlying matrix diffusion effect. Therefore, this parameter is mainly used for rapid estimation in this embodiment, and does not replace the physical diffusion parameter obtained by inversion based on the fracture-enhanced bilayer model in the aforementioned embodiments.

[0263] The approximate active area depth H(x, y) is calculated for all remote sensing pixels in the study area, generating a preliminary spatial distribution map of the active area depth. A reliability index is also established.

[0264] CI = 1 - σ z / z*;

[0265] Where CI is the credibility index, and σ z Let z be the standard deviation of the active zone depth, and z* be the average active zone depth.

[0266] When CI > 0.8, the active zone identification result is considered reliable. In this embodiment, the identified active zone depth distribution is basically consistent with the field monitoring results. The results show that this embodiment can achieve rapid preliminary identification of the active zone depth of expansive soil slopes without the need for extensive borehole monitoring. This embodiment does not rely on the fracture-enhanced two-layer model and multi-frequency joint constraints as its core, but instead uses a single-layer equivalent diffusion approximation to quickly estimate the hydrothermal control active zone of expansive soil slopes. It has low data requirements and a simple calculation process, making it suitable for use when fracture stratification information is insufficient, high-order frequency domain features are unstable, or when only a rapid preliminary estimate of the active zone depth is needed in engineering applications. In scenarios requiring higher accuracy in characterizing the preferential flow of surface fractures and the diffusion stratification effect of deep matrix in expansive soil, the fracture-enhanced two-layer model implementation is preferred for active zone depth inversion.

[0267] This invention replaces the traditional single-layer diffusion front formula by constructing a depth integral deformation response transfer function derived based on zero-flux boundaries and the complex domain. This breaks through the mathematical bottleneck of the phase delay approaching the upper limit of the constant in the original model and solves the problem of underestimation of the depth of deep active zones. A fracture-enhanced two-layer physical architecture is reconstructed, and a surface infiltration low-pass filter is connected in series at the data input end. This successfully decouples the heterogeneous hydraulic responses of surface water retention, shallow fracture preferential flow (in phase), and deep matrix slow diffusion flow (phase delay), eliminating the systematic error of mistaking surface evaporation and infiltration loss for deep physical thickness. Specific dual-frequency features are extracted using frequency domain decomposition, and cross-frequency hard constraints are constructed using the physical correlation between high and low frequency water infiltration depths. This transforms the ill-conditioned underdetermined system into a closed-loop optimal overdetermined equation set, achieving convergent solution of multidimensional parameters. Thus, a complete physical closed loop was realized, from periodic rainfall input - groundwater migration - depth integral deformation response - surface observation phase delay - active area depth inversion, avoiding the problem of unclear mechanism caused by using phase lag only as an empirical statistic.

[0268] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for identifying hydrothermal control activity zones on expansive soil slopes, characterized in that, include: Obtain the time series of vertical deformation of each pixel within the target slope area and the time series of rainfall within the corresponding time period; Processing time series of vertical deformation and rainfall, extracting periodic response features characterizing the moisture-driven deformation process; A hydrothermal transfer model reflecting the moisture diffusion mechanism and depth integral deformation mechanism of expansive soil was constructed. The hydrothermal transfer model established a physical mapping relationship between the theoretically predicted response characteristics and the normalized active zone depth. The periodic response characteristics are substituted into the water and heat transfer model for inversion matching to obtain the normalized active area depth. The total active area depth of each pixel is then calculated by combining the pre-calibrated water diffusion coefficient. The periodic response characteristics include phase delay and amplitude ratio at at least two preset frequencies; Extracting periodic response features characterizing water-driven deformation processes, including: Frequency domain decomposition was performed on the vertical deformation time series and the rainfall time series to extract the deformation phase, deformation amplitude, rainfall phase and rainfall amplitude at each preset frequency; The corresponding phase delay is calculated based on the deformation phase and rainfall phase at the same preset frequency; The corresponding amplitude ratio is calculated based on the deformation amplitude and rainfall amplitude at the same preset frequency; The hydrothermal transfer model is a fracture-enhanced two-layer model; The fracture-enhanced bilayer model is divided into a fracture layer and a intact matrix layer along the depth direction of the active region, and a zero-flux boundary condition is applied at the bottom of the active region. In this study, the moisture response within the fractured layer is set to be in phase with the surface driving signal, while the moisture within the intact matrix layer propagates through diffusion flow and generates a phase delay.

2. The method according to claim 1, characterized in that, The fracture-enhanced two-layer model establishes a physical mapping relationship through a depth integral deformation response transfer function containing complex values. Its analytical expression is as follows: Φ bi (ξ2, r) = r + (1 - r) * tanh[(1 + i)ξ2] / [(1 + i)ξ2]; where Φ bi represents the theoretically predicted response characteristic; r represents the normalized crack depth; ξ2represents the normalized thickness of the intact substrate layer; i represents the imaginary unit; tanh represents the hyperbolic tangent function; The normalized active zone depth is composed of the normalized fracture depth to be inverted and the normalized thickness of the intact matrix layer.

3. The method according to claim 2, characterized in that, The periodic response features include phase delays extracted from at least two preset frequencies of vertical deformation time series and rainfall time series; The periodic response characteristics are substituted into the water and heat transfer model for inversion matching, including: Using the phase delay at at least two preset frequencies as the observation end and the depth integral deformation response transfer function at different frequencies as the model end, a joint inversion equation set is constructed. Solve the joint inversion equations simultaneously to obtain the normalized fracture depth and the normalized thickness of the intact matrix layer.

4. The method according to claim 3, characterized in that, When constructing the joint inversion equation set, cross-frequency constraints are constructed using the physical relationship between water infiltration depths at different frequencies. Set at least two preset frequencies, including an annual cycle frequency and a semi-annual cycle frequency. The expression for the cross-frequency constraint is: ξ 2，2 =ξ 2，1 *2 1 / 2 ; Where, ξ 2，1 ξ represents the normalized thickness of the intact matrix layer at the annual cycle frequency. 2，2 This represents the normalized thickness of the complete matrix layer at a six-month cycle frequency. In the joint inversion equations, the normalized fracture depth is set to remain constant at different preset frequencies.

5. The method according to claim 4, characterized in that, The inversion yields the normalized fracture depth and the normalized thickness of the intact matrix layer, including: A uniform grid is constructed within a preset parameter space, and the sum of squares of the dual-frequency phase residuals between the phase predicted by the depth integral deformation response transfer function and the phase extracted from the observation end is calculated one by one at each grid point. The initial solution is formed by selecting the parameters corresponding to the grid point that minimizes the sum of squares of the dual-frequency phase residuals. Starting with the initial solution, the Levenberg-Marquardt algorithm is used to perform nonlinear optimization on the sum of squares of the dual-frequency phase residuals, and the optimal normalized fracture depth and normalized thickness of the intact matrix layer are obtained by inversion.

6. The method according to claim 1, characterized in that, Before calculating the total active area depth of each pixel, a reliable pixel screening step using a preset threshold is also included, specifically including: The frequency domain signal-to-noise ratio, inversion fitting residual, and linear deformation rate extracted from the vertical deformation time series of each pixel are obtained. Determine whether each pixel simultaneously meets the pre-configured lower threshold for signal-to-noise ratio, upper threshold for residual, and upper threshold for linear deformation rate; Cells that do not meet any threshold condition are removed to obtain a set of reliable cells that have passed the screening. The total depth of the active area is calculated only for the cells in the set of reliable cells.

7. The method according to claim 1, characterized in that, The pre-calibrated moisture diffusion coefficient was obtained through the following pre-construction steps: Obtain the measured active zone depth and measured fracture depth at a given calibration profile; Extract the periodic response features at the calibration profile; The measured active zone depth, measured fracture depth, and periodic response characteristics at the calibration profile are substituted into the phase equation of the hydrothermal transfer model for inverse deduction, and the corresponding physical diffusion parameters at the calibration profile are calculated. The physical diffusion parameters are used as pre-calibrated moisture diffusion coefficients required for inversion calculations for each pixel.

Images