A terrain soft constraint time sequence InSAR three-dimensional deformation joint solution method

Abstract

This application belongs to the technical fields of synthetic aperture radar interferometry, radar remote sensing deformation monitoring, mountain geological disaster monitoring, and three-dimensional surface deformation inversion. Specifically, it discloses a topographically constrained time-series InSAR three-dimensional deformation joint solution method, including: acquiring multi-temporal SAR images, digital elevation models, and satellite imaging geometric parameters; performing time-series InSAR processing to obtain the line-of-sight cumulative deformation time series and observation quality information; extracting topographic parameters; constructing topographic constraints and establishing topographic constraint applicability indices to determine constraint weights; establishing line-of-sight observation relationships and jointly constructing a three-dimensional deformation solution model with topographic constraints in a soft-constraint form; determining line-of-sight observation weights based on observation quality information, using the constraint weights as topographic constraint weights, constructing a joint weight matrix, and solving the model to obtain the three-dimensional deformation components. This application can improve the stability and reliability of three-dimensional deformation inversion in complex mountainous areas.

Description

Technical Field

[0001] This application belongs to the technical fields of synthetic aperture radar interferometry, radar remote sensing deformation monitoring, mountain geological disaster monitoring and three-dimensional surface deformation inversion, and more specifically, relates to a joint solution method for three-dimensional deformation of soft-constrained temporal InSAR. Background Technology

[0002] Temporal InSAR technology can acquire time series of surface deformation along the satellite line of sight using multi-temporal SAR imagery. It boasts advantages such as all-weather, all-time coverage, large scale, and high precision, and has been widely applied in monitoring landslides, glaciers, freeze-thaw slopes, mining subsidence, urban surface subsidence, and geological hazards. However, InSAR observations essentially only obtain a one-dimensional projection of surface deformation along the satellite line of sight, and cannot directly obtain the three-dimensional deformation components in the east-west, north-south, and vertical directions. Existing methods typically estimate east-west and vertical deformation through joint decomposition of ascending and descending InSAR data. However, due to the observation geometry limitations of near-polar orbit SAR satellites, the satellite line of sight is less sensitive to the north-south deformation component, making it difficult to stably recover the north-south deformation.

[0003] In complex mountainous regions such as high mountains and canyons, landslides, glaciers, and freeze-thaw slopes, surface movements often occur along the slope, the main sliding direction, or the glacier flow direction, with the north-south component potentially accounting for a large proportion. Ignoring north-south deformation in the 3D decomposition, or simply assuming zero north-south deformation, can lead to significant biases in the estimation of east-west and vertical deformation. Existing studies have utilized topographic constraint methods such as surface parallel flow models, slope aspect constraint models, and surface tangent plane constraint models to assist in InSAR 3D deformation inversion. However, existing methods still have the following shortcomings: they mostly treat terrain constraints as fixed hard constraints, which are difficult to adapt to situations such as local non-slope aspect movement, shear deformation, and non-rigid movement of landslides in complex mountainous areas; they lack pixel-level judgment on the applicability of terrain constraints, which can easily lead to misapplication of constraints in flat areas, areas with abrupt changes in slope aspect, low coherence areas, or areas with high terrain noise; they lack unified weight control for factors such as InSAR observation quality, slope aspect stability, slope magnitude, and terrain curvature; they do not adequately evaluate the ill-conditioning, residuals, and reliability of the three-dimensional solution model; and most methods focus on single-period displacement or average rate decomposition, lacking a three-dimensional time series joint solution framework oriented towards time-series InSAR cumulative deformation.

[0004] Therefore, improving the stability and reliability of time series inversion of three-dimensional deformation in complex mountainous areas is an urgent problem to be solved. Summary of the Invention

[0005] To address the shortcomings of existing technologies, the purpose of this application is to provide a joint solution method for three-dimensional deformation of terrain soft-constrained time-series InSAR, which can improve the stability and reliability of time-series inversion of three-dimensional deformation in complex mountainous areas.

[0006] To achieve the above objectives, in a first aspect, this application provides a method for joint calculation of three-dimensional deformation in terrain-constrained temporal InSAR, comprising the following steps: S10: Acquire multi-temporal SAR images, digital elevation model data, and satellite imaging geometric parameters of the study area; S20, perform time-series InSAR processing on the multi-temporal SAR images to obtain the satellite line-of-sight cumulative deformation time series and InSAR observation quality information at each observation time; S30, extract at least one topographic parameter from the slope, aspect, surface normal, topographic curvature, or main direction of motion of the study area based on the digital elevation model data; S40, construct terrain constraints based on the terrain parameters, the terrain constraints including at least one of surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints; construct terrain constraint applicability indices based on the terrain parameters and the InSAR observation quality information, and determine the constraint weights of the terrain constraints based on the terrain constraint applicability indices; S50, establish the line-of-sight observation relationship between the cumulative deformation along the satellite line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component based on the satellite imaging geometric parameters, and establish time continuity constraints for the three-dimensional deformation components or three-dimensional deformation rate at continuous observation times; combine the line-of-sight observation relationship, the time continuity constraints, and the terrain constraints in the form of soft constraints that allow the existence of constraint residuals to construct a time-series-oriented three-dimensional deformation joint solution model; S60, determine the line-of-sight observation weight based on the InSAR observation quality information, use the constraint weight as the terrain constraint weight, and combine it with the time continuity constraint to determine the time continuity constraint weight; construct a joint weight matrix including the line-of-sight observation weight, terrain constraint weight, and time continuity constraint weight, and solve the time-series-oriented three-dimensional deformation joint solution model to obtain the east-west deformation component, north-south deformation component, and vertical deformation component at each observation time.

[0007] As a further preferred embodiment, the time-series InSAR processing includes at least one of image registration, interferometric network construction, differential interferometric processing, terrain phase removal, phase filtering, phase unwrapping, atmospheric delay correction, reference point correction, and time-series deformation inversion; the terrain parameters include at least one of slope, aspect, surface normal vector, terrain curvature, local aspect dispersion, geomorphic unit boundary, main slip direction, glacier flow direction, and geological structure direction.

[0008] As a further preferred embodiment, when the terrain constraint is a surface parallel flow constraint, the surface parallel flow constraint is expressed in the form of a soft constraint as follows: ,in This represents the constraint function constructed based on the surface parallel flow assumption. , , These represent the east-west deformation component, the north-south deformation component, and the vertical deformation component, respectively. This represents the surface parallel flow constraint residual.

[0009] As a further preferred embodiment, the soft constraint form refers to: allowing the existence of terrain constraint residuals in the joint solution model, and weighting the terrain constraint residuals according to the terrain constraint applicability index and the InSAR observation quality.

[0010] As a further preferred option, the terrain constraint applicability index The terrain constraint applicability index is determined based on at least one of slope, aspect local dispersion, terrain curvature, InSAR coherence coefficient, or phase unwrapping quality; the calculation formula for the terrain constraint applicability index is as follows: ,in Indicates slope; This represents slope aspect stability, which is calculated from the slope aspect dispersion within a local window. Indicates the curvature of the terrain; Indicates the InSAR coherence coefficient; This indicates the phase unwrapping quality.

[0011] As a further preferred embodiment, the joint weight matrix is ​​expressed as: ,in The line-of-sight observation weight is determined based on at least one of the InSAR coherence coefficient, phase standard deviation, unwrapping residual, and atmospheric delay correction quality. The terrain constraint weight is determined based on at least one of the following: slope, aspect stability, terrain curvature, geomorphic unit consistency, and confidence of the main movement direction. The weights representing time continuity constraints are determined based on at least one of time interval, degree of deformation nonlinearity, and time-series residuals. The regularization weight is determined based on at least one of the following: model ill-conditionedness, noise level, spatial smoothing requirements, and temporal stability requirements.

[0012] As a further preferred method, the three-dimensional deformation joint solution model is solved using the regularized least squares method, and the solution formula is as follows:

[0013] in, This represents the three-dimensional deformation vector to be determined, in a single-stage solution mode. , , , These represent the east-west, north-south, and vertical deformation components, respectively; in the time-series joint solution mode... It is an augmented vector containing the three-dimensional deformation components at all observation times; Represents the joint weight matrix; Represents the joint design matrix; Represents the observation vector; Represents the regularization parameter; This represents a spatial smoothing, temporal smoothing, or prior constraint matrix.

[0014] As a further preferred embodiment, during or after solving the three-dimensional deformation joint solution model, an ill-conditioning diagnosis is performed. The ill-conditioning diagnosis includes at least one of the following: joint design matrix condition number, parameter covariance matrix, posterior unit weight variance, line-of-sight residual, terrain constraint residual, and temporal deformation jump degree. When the ill-conditioning diagnosis index exceeds a preset threshold, the terrain constraint weight is adjusted, a regularization term is added, the reliability level of the three-dimensional deformation result is reduced, or the corresponding region is marked as an unsuitable three-dimensional solution region.

[0015] As a further preferred option, time continuity constraints or piecewise linear constraints are introduced into the three-dimensional deformation components or three-dimensional deformation rates at continuous observation times to suppress the three-dimensional deformation anomalous jumps caused by single-period InSAR noise.

[0016] Secondly, this application provides a terrain-soft-constrained temporal InSAR three-dimensional deformation joint solution system for implementing the method described in any one of the above, comprising: The data acquisition module is used to acquire multi-temporal SAR images, digital elevation model data, and satellite imaging geometric parameters of the study area; The temporal InSAR processing module is used to perform temporal InSAR processing on multi-temporal SAR images to obtain the satellite line-of-sight cumulative deformation time series and InSAR observation quality information at each observation time. The line-of-sight observation modeling module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite's line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component, based on the satellite imaging geometric parameters. The terrain parameter extraction module is used to extract at least one terrain parameter from the slope, aspect, surface normal vector, terrain curvature, or main direction of motion of the study area based on the digital elevation model data. The terrain constraint construction module is used to construct terrain constraints based on the terrain parameters. The terrain constraints include at least one of surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints. Then, the module constructs a terrain constraint applicability index based on the terrain parameters and the InSAR observation quality information, and determines the constraint weight of the terrain constraint based on the terrain constraint applicability index. The model building module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component based on the satellite imaging geometric parameters, and to establish time continuity constraints on the three-dimensional deformation components or three-dimensional deformation rate at continuous observation times; the line-of-sight observation relationship, the time continuity constraints, and the terrain constraints are combined in the form of soft constraints that allow the existence of constraint residuals to construct a time-series-oriented three-dimensional deformation joint solution model. The joint solution module is used to determine the line-of-sight observation weights based on the InSAR observation quality information, use the constraint weights as terrain constraint weights, and combine them with the time continuity constraint to determine the time continuity constraint weights; construct a joint weight matrix containing the line-of-sight observation weights, terrain constraint weights, and time continuity constraint weights, and solve the time-series-oriented three-dimensional deformation joint solution model to obtain the east-west deformation components, north-south deformation components, and vertical deformation components at each observation time; The quality evaluation module is used to perform pathological diagnosis, residual evaluation, and reliability classification on the three-dimensional deformation solution results; The results output module is used to output three-dimensional deformation time series, three-dimensional deformation rate field, deformation results along slope, and reliability evaluation results.

[0017] Compared with the prior art, this application has the following beneficial effects: This application provides a joint solution method for three-dimensional deformation in InSAR with soft-constraint terrain. By introducing terrain constraint applicability discrimination, soft-constraint model construction, adaptive weighting mechanism, and improved ill-condition diagnosis, it significantly improves the stability and reliability of time series inversion of three-dimensional deformation in complex mountainous areas. Specifically, it achieves the following: (1) Alleviating the weak sensitivity of the north-south direction and improving the stability of the model. By extracting terrain parameters and combining line-of-sight (LOS) observations to construct a solution model, the matrix condition number is greatly improved, which effectively solves the serious ill-conditioned problem of the model caused by the insensitivity of LOS observations to north-south deformation and improves the recovery accuracy of the north-south component.

[0018] (2) A soft constraint mechanism is adopted to avoid the bias caused by hard constraints. The topographic priors such as surface parallel flow are introduced in the form of "soft constraints", which allows the existence of topographic constraint residuals. This avoids the serious solution bias caused by the traditional hard constraints that force the surface movement to be limited to a single direction, and ensures that the inversion results are more objective.

[0019] (3) Adaptive evaluation and weighting to enhance terrain adaptability. A pixel-level applicability index is constructed by combining terrain curvature, slope and InSAR observation quality to avoid misuse of constraints in flat or low coherence areas; at the same time, a joint weight matrix is ​​dynamically generated for adaptive weighting, which greatly improves the inversion stability in complex geological and noisy environments.

[0020] (4) Introducing time continuity constraints to ensure time series reliability. A solution framework is specifically designed for time series cumulative deformation. Time continuity constraints are added at continuous observation times, which effectively suppresses abnormal jumps caused by single-period InSAR noise and ensures that the output time series results are smooth and of high quality.

[0021] (5) Implement pathological diagnosis and grading, and quantify data confidence. Introduce pathological diagnosis indicators such as joint design matrix condition number and residuals, adaptively adjust weights or regularization terms when exceeding thresholds, and perform reliability grading on the final three-dimensional deformation results. This provides highly reliable three-dimensional deformation data support for subsequent disaster identification, risk analysis and engineering early warning. Attached Figure Description

[0022] Figure 1 This is the overall flowchart of the method provided in this application; Figure 2 This is a schematic diagram of the InSAR line-of-sight observation geometry provided in this application; Figure 3 This is a schematic diagram of the complex mountainous terrain constraints provided in this application; Figure 4 This is a flowchart of the terrain constraint applicability evaluation provided in this application; Figure 5 This is a structural diagram of the joint solution model provided in this application; Figure 6 This is a system architecture block diagram provided in this application; Figure 7 This is a comparison diagram of the north-south (N-direction) deformation rate inversion results of the simulated landslide area provided in the embodiments of this application; Figure 8 This is a bar chart comparing the root mean square error (RMSE) of three-dimensional deformation in the east-west, north-south, and vertical directions of the three solution schemes provided in the embodiments of this application. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0024] To address the issues of insensitivity of InSAR observations to north-south deformation in high-altitude areas, the ill-conditioned nature of traditional 3D deformation decomposition models, the difficulty in stably recovering north-south deformation, and the overlapping of vertical and horizontal deformations, which lead to low stability and reliability in the time series inversion of 3D deformation in complex mountainous terrain, this application proposes a terrain-constrained time-series InSAR 3D deformation joint solution method for complex mountainous terrain. Based on the cumulative deformation of time-series InSAR, a terrain constraint applicability judgment is introduced. Terrain constraints, including surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints, are added to the joint solution model in the form of soft constraints. Combined with adaptive adjustment of constraint weights, ill-conditioning diagnosis, temporal continuity control, and reliability grading mechanisms, this method improves the stability and reliability of the time series solution for 3D deformation in complex mountainous terrain.

[0025] like Figure 1 As shown, this application provides a method for joint calculation of three-dimensional deformation using soft-constrained temporal InSAR for complex mountainous terrain, comprising the following steps: S1, data acquisition; S2, temporal InSAR processing; S3, acquisition of LOS cumulative deformation time series; S4, terrain parameter extraction; S5, terrain constraint applicability evaluation; S6, soft-constrained joint model construction; S7, weighted, regularized, or robust calculation; S8, ill-conditioning diagnosis; S9, three-dimensional deformation time series output.

[0026] S8 is used as an independent step to evaluate the stability, solvability of the north-south component, constraint effectiveness, and result uncertainty of the joint model, and to adjust the constraint weights, model structure, or regularization parameters in S5, S6, or S7 when the preset conditions are not met.

[0027] In this application, terrain constraints refer to constraints constructed from prior information about terrain or geomorphology, such as digital elevation models, slope, aspect, surface normal vector, terrain curvature, principal slip direction, glacier flow direction, or geological structure direction, used to limit the physical rationality of the three-dimensional deformation solution. Soft constraints refer to a constraint implementation method that allows the existence of constraint residuals in the joint solution model and assigns weights based on the reliability of the constraints. Accordingly, soft terrain constraints refer to incorporating terrain constraints into the three-dimensional deformation joint solution model in the form of soft constraints, enabling it to provide prior terrain information while allowing for a certain deviation of the actual deformation from the ideal aspect, principal movement direction, or surface tangent plane.

[0028] Temporal InSAR processing of multi-temporal SAR images includes: image registration; interferometric network construction; differential interferometry; terrain phase removal; phase filtering; coherence coefficient calculation; phase unwrapping; atmospheric delay correction; reference point correction; and time-series deformation inversion.

[0029] Temporal InSAR methods can include SBAS-InSAR, PS-InSAR, DS-InSAR, multi-master image InSAR, or other temporal InSAR methods.

[0030] like Figure 2 As shown, the projection relationship between satellite line-of-sight deformation and three-dimensional deformation components is established based on the satellite incident angle, satellite azimuth angle, or line-of-sight unit vector.

[0031] For a given pixel on the Earth's surface, the three-dimensional deformation vector can be expressed as: ,in, , , These represent the east-west, north-south, and vertical deformation components, respectively.

[0032] like Figure 3 As shown, topographic parameters for complex mountainous areas are extracted based on digital elevation models, including at least one of the following: Slope; aspect; surface normal vector; topographic curvature; local aspect dispersion; geomorphic unit boundary; main slip direction; glacier flow direction; geological structure direction.

[0033] The main direction of motion can be determined by the slope aspect of the DEM, or by the landslide boundary, glacier streamlines, geological survey data, optical image matching results, or historical deformation direction estimation results.

[0034] To avoid misuse of terrain constraints in unsuitable areas, this application constructs a terrain constraint applicability index.

[0035] like Figure 4 As shown, the terrain constraint applicability index can be determined by at least one of the following factors: Slope magnitude; local slope aspect stability; local slope aspect dispersion; terrain curvature; geomorphic unit consistency; InSAR coherence coefficient; phase unwrapping quality; atmospheric correction quality; external validation data accuracy; confidence of main motion direction.

[0036] For example, terrain constraint applicability indices can be defined:

[0037] in, Indicates slope; This represents slope aspect stability, which can be calculated from the slope aspect dispersion within a local window. The smaller the local slope aspect dispersion, the higher the slope aspect stability. Represents the curvature of the terrain. Within the same geomorphic unit, the smaller the absolute value of the curvature, the better the local slope continuity and the higher the reliability of the terrain constraint. Indicates the InSAR coherence coefficient; This indicates the phase unwrapping quality.

[0038] Indicates the applicability index of terrain constraints. Indicates by Terrain constraint weights are obtained through threshold discrimination and weight mapping. Used to evaluate whether a particular pixel or raster cell is suitable for enabling terrain constraints. This is used to control the contribution of the terrain constraint in the joint solution model.

[0039] When the slope is steep, the slope aspect is stable, the coherence coefficient is high, and the unwrapping quality is good, increase the weight of the terrain constraint; when the area is flat, the slope aspect changes abruptly, the coherence is low, or the unwrapping quality is poor, decrease the weight of the terrain constraint, or do not enable the constraint.

[0040] like Figure 5 As shown, the constraints in this application are not limited to a single form and may include topographic constraints and auxiliary constraints. Topographic constraints include surface parallel flow constraints, slope aspect or main motion direction constraints, and surface tangent plane constraints, etc.; auxiliary constraints include spatial consistency constraints and temporal continuity constraints, etc. These are described in detail below: 1. Terrain Constraints (1) Soft constraint form of surface parallel flow constraint Surface parallel flow constraints are a preferred form of terrain constraint. When surface motion mainly occurs along the slope, aspect, or main motion direction, surface parallel flow constraints can be established based on the aspect, surface normal vector, or main motion direction extracted from the digital elevation model. Unlike hard constraints that strictly limit surface motion to the slope direction or main motion direction, this application incorporates surface parallel flow constraints as soft constraints into the three-dimensional deformation joint solution model, allowing the existence of surface parallel flow constraint residuals, and assigning corresponding weights based on terrain constraint applicability indicators, InSAR observation quality, and the reliability of the main motion direction.

[0041]

[0042] in, This represents the constraint function constructed based on the surface parallel flow assumption. This represents the surface parallel flow constraint residual. When... When it approaches zero, it indicates that the surface motion is closer to the ideal state of surface parallel flow; when When the value is not zero, it indicates that the actual deformation direction is allowed to deviate to a certain extent from the slope direction or the main motion direction.

[0043] (2) Slope aspect or main motion direction constraint When the direction of horizontal movement on the ground is approximately along the slope aspect or the main sliding direction, a constraint on the direction of horizontal movement can be established. This constraint can be determined based on the slope aspect angle, glacier flow direction angle, landslide main sliding direction angle, or geological structure direction angle.

[0044] (3) Surface tangent plane constraint When surface motion mainly occurs along the surface tangent plane, a tangent plane constraint can be established based on the surface normal vector, indicating that the three-dimensional deformation vector is approximately orthogonal to the surface normal vector.

[0045] 2. Auxiliary constraint terms (1) Spatial consistency constraints For adjacent pixels within the same landslide body, glacial basin, or geomorphic unit, a local neighborhood deformation direction consistency constraint can be introduced to suppress abnormal directional jumps caused by noise.

[0046] (2) Time continuity constraint For continuous observation times, temporal continuity constraints on three-dimensional deformation components or deformation rates can be introduced to reduce the impact of single-period InSAR noise on the three-dimensional decomposition results.

[0047] The line-of-sight relationship and terrain constraint relationship are unified into a joint model:

[0048] in, The observation vector includes the deformation observations along the line of sight and constraint terms; The joint design matrix includes LOS projection coefficients and terrain constraint coefficients; Let be the three-dimensional deformation components to be determined; This is the residual vector.

[0049] In this application, the terrain constraint term is preferably added to the joint model as a soft constraint term, and its residual is not forced to be zero, but is assigned corresponding weights according to the applicability of the terrain constraint and the observation quality.

[0050] When the joint model simultaneously incorporates time continuity constraints and regularization constraints, its joint weight matrix can be expressed as:

[0051] in, The weight of InSAR line-of-sight observations can be determined based on coherence coefficient, phase standard deviation, unwrapping quality, atmospheric correction quality, etc. The terrain constraint weight can be determined based on factors such as slope, aspect stability, terrain curvature, geomorphic unit consistency, and the credibility of the main motion direction. Terrain constraints include at least one of surface parallel flow constraints, aspect constraints, surface tangent plane constraints, or main motion direction constraints. The weights representing time continuity constraints can be determined based on time intervals, the degree of deformation nonlinearity, or time-series residuals. The regularization weight can be determined based on the ill-conditioned nature of the model, noise level, spatial smoothing requirements, or temporal stability requirements.

[0052] By setting adaptive weights, this application can dynamically coordinate the contributions of LOS observation terms, terrain constraints, temporal continuity constraints, and regularization terms under different observation qualities and terrain conditions, thereby improving the stability and reliability of 3D deformation inversion in complex mountainous areas.

[0053] The joint model can be solved using any of the following methods: Weighted least squares; regularized least squares; robust estimation; iterative reweighted least squares; Bayesian estimation; piecewise joint time series estimation.

[0054] The weighted least squares form can be written as:

[0055] Regularized least squares can be written as:

[0056] in, Represents the regularization parameter. This represents a spatial smoothing, temporal smoothing, or prior constraint matrix.

[0057] To address the pathological issues caused by the weak sensitivity to north-south deformation in high-altitude regions, this application performs pathological diagnosis during or after the model solution process.

[0058] Pathological diagnostic indicators include at least one of the following: Joint design matrix condition number; parametric covariance matrix; posterior unit weight variance; LOS residual; terrain constraint residual; degree of anomaly in 3D deformation direction; degree of temporal deformation jump.

[0059] When the ill-conditioned nature of the model exceeds a preset threshold, at least one of the following actions will be performed: Adjust terrain constraint weights; add regularization terms; reduce the 3D solution level of low-confidence regions; output only LOS results or 2D decomposition results; mark as low-confidence 3D deformation regions.

[0060] The final output can be divided into: High-reliability 3D deformation region; Medium-reliability 3D deformation region; Low-reliability 3D deformation region; Unsuitable for 3D solution.

[0061] In a general implementation, this application is carried out according to the following process: S1, acquire multi-temporal SAR images, orbit data, DEM and imaging geometry parameters; S2, perform time-series InSAR processing on SAR images to obtain observation quality information such as unwrapped phase, coherence coefficient, and unwrapping quality; S3, convert the unwrapped phase into a LOS cumulative deformation time series; S4 extracts parameters such as slope, aspect, surface normal vector, and terrain curvature from the DEM; S5, build And map to get ; S6, construct the LOS observation term and the terrain soft constraint term to form a joint model; S7 can be solved using weighted, regularized, or robust methods. S8 performs condition number, residual, covariance, and temporal continuity diagnosis; S9 outputs the E / N / U three-dimensional deformation time series and reliability classification results.

[0062] Based on the same inventive concept, such as Figure 6 As shown, this application also provides a terrain soft-constraint temporal InSAR three-dimensional deformation joint solution system for implementing the above method, including: The data acquisition module is used to acquire multi-temporal SAR images, digital elevation model data, and satellite imaging geometric parameters of the study area; The temporal InSAR processing module is used to perform temporal InSAR processing on multi-temporal SAR images to obtain the satellite line-of-sight cumulative deformation time series and InSAR observation quality information at each observation time. The line-of-sight observation modeling module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite's line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component, based on the satellite imaging geometric parameters. The terrain parameter extraction module is used to extract at least one terrain parameter from the slope, aspect, surface normal vector, terrain curvature, or main direction of motion of the study area based on the digital elevation model data. The terrain constraint construction module is used to construct terrain constraints based on terrain parameters. The terrain constraints include at least one of surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints. Then, the terrain constraint applicability index is constructed based on the terrain parameters and InSAR observation quality information, and the constraint weight of the terrain constraint is determined based on the terrain constraint applicability index. The model building module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component based on the satellite imaging geometric parameters, and to establish time continuity constraints on the three-dimensional deformation components or three-dimensional deformation rate at continuous observation times; the line-of-sight observation relationship, time continuity constraints, and terrain constraints are combined in the form of soft constraints that allow the existence of constraint residuals to construct a time-series-oriented three-dimensional deformation joint solution model. The joint solution module is used to determine the line-of-sight observation weights based on InSAR observation quality information, use the constraint weights as terrain constraint weights, and combine them with the time continuity constraint to determine the time continuity constraint weights; construct a joint weight matrix that includes the line-of-sight observation weights, terrain constraint weights, and time continuity constraint weights, and solve the time-series-oriented three-dimensional deformation joint solution model to obtain the east-west deformation components, north-south deformation components, and vertical deformation components at each observation time. The quality evaluation module is used to perform pathological diagnosis, residual evaluation, and reliability classification on the three-dimensional deformation solution results; The results output module is used to output three-dimensional deformation time series, three-dimensional deformation rate field, deformation results along slope, and reliability evaluation results.

[0063] The present application will now be described in conjunction with specific embodiments.

[0064] It should be noted that the method described in this application can be applied to various complex topographic deformation scenarios, including but not limited to landslides in mountainous areas, glacier flows, freeze-thaw slopes, mining slopes, accumulation slopes, debris flow source areas, and rockfall areas. The controlling factors and movement patterns of surface deformation differ in different application scenarios; therefore, the form of topographic constraints, constraint weights, and reliability evaluation methods adopted can be adjusted according to the research object. For example, in mountainous landslide scenarios, aspect, main sliding direction, or surface parallel flow constraints can be emphasized; in glacier scenarios, constraints can be constructed by combining ice surface aspect, glacier streamlines, or main flow direction; in freeze-thaw slope or mining slope scenarios, a joint judgment can be made by combining aspect constraints, surface tangent plane constraints, and vertical settlement characteristics.

[0065] To fully illustrate the content of this application, Examples 1 to 3 below will use examples of mountain landslides, three-dimensional flow monitoring of glaciers, and deformation monitoring of freeze-thaw slopes or mining slopes to illustrate the specific implementation process of this application in different terrain scenarios. In addition, Example 4 will construct a simulation working condition with known true values ​​to conduct a quantitative comparative test and analysis of the method of this application and the existing technical solutions to verify the technical effect of this application.

[0066] Example 1: Three-dimensional deformation time series inversion of landslides in mountainous areas A high-altitude landslide area was selected as the study area. This area has a large slope and obvious topographic relief. The landslide movement mainly occurs along the slope or the main sliding direction, and the north-south deformation may account for a large proportion.

[0067] First, multi-temporal SAR images and DEM data covering the study area are acquired. Then, the SAR images are processed using time-series InSAR to obtain the line-of-sight cumulative deformation time series.

[0068] Then, the slope, aspect, local aspect dispersion, and topographic curvature are calculated using the DEM to construct a topographic constraint applicability index. For areas with large slopes, stable aspects, and high coherence, the surface parallel flow constraint weight is increased; for flat areas, areas with abrupt aspect changes, or low coherence areas, the constraint weight is decreased.

[0069] Subsequently, LOS observation relationships and soft constraints on surface parallel flow were established, and a joint solution model was constructed. Weighted least squares or regularized least squares methods were used to solve the three-dimensional deformation components in the east-west, north-south, and vertical directions.

[0070] Finally, the results are evaluated for condition number, residuals, and temporal continuity, and three-dimensional deformation time series with high, medium, and low reliability levels are output.

[0071] Example 2: Three-dimensional flow monitoring of glaciers High-altitude glaciers were selected as the research subject. Glacier movement is usually driven by gravity, and the main direction of movement is related to the ice surface aspect or glacier streamlines.

[0072] First, multi-temporal SAR images and DEM data of the glacier area are acquired, and then the line-of-sight deformation time series is obtained through time-series InSAR processing.

[0073] Then, the ice surface slope and aspect are extracted from the DEM, or the main flow direction is extracted by combining the glacier basin boundary and glacier streamlines. This main flow direction is used as the topographic constraint direction and added to the joint solution model as a soft constraint.

[0074] By combining InSAR line-of-sight observations with glacier flow direction constraints, the east-west, north-south, and vertical motion components of the glacier are obtained, and the three-dimensional motion time series and reliability evaluation results of the glacier are output.

[0075] Example 3: Monitoring Deformation of Freeze-Thaw Slopes or Mining Area Slopes For freeze-thaw slopes or mining slopes, surface deformation may simultaneously include subsidence, horizontal slippage, and local slope movement.

[0076] This application can extract slope, aspect, and surface normal vectors from the DEM, and construct constraint applicability indices by combining coherence coefficient and unwrapping quality. When the slope motion characteristics are obvious, surface tangent plane constraints or aspect constraints are introduced; when vertical settlement is determined to be dominant, the aspect constraint weight is reduced to avoid incorrect constraints.

[0077] By using adaptive constraint weights and reliability grading, this application can adapt to three-dimensional deformation inversion under complex motion modes.

[0078] Example 4: Comparison Experiment of Existing Technical Solutions Based on Simulated Working Conditions To further verify the effectiveness and superiority of the proposed method for joint calculation of three-dimensional deformation in complex mountainous terrain using time-series InSAR with soft constraints, this embodiment constructs a complex mountainous landslide simulation with known true values ​​and conducts a comparative experiment with existing technologies.

[0079] (1) Simulation conditions and experimental setup. The experiment simulated a complex mountain landslide with an average slope of about 30° and a main slope orientation of east. The true three-dimensional motion of the landslide was set to include not only the sliding along the slope (mainly manifested as east-west and vertical displacement), but also a north-south (N) sliding component of about 15 mm / yr due to local terrain compression. Based on the typical near-polar orbit flight parameters of the Sentinel-1 satellite (ascending azimuth angle of about -12°, descending azimuth angle of about -168°, and incident angle of about 39°), the true three-dimensional deformation was projected onto the line of sight (LOS) of the satellite, and random Gaussian noise with a signal-to-noise ratio of 15 dB was added to generate simulated observation data in the LOS direction of the ascending and descending orbits.

[0080] (2) Comparison of Scheme Designs. Using the same LOS observation simulation data and DEM data, the following three solution schemes were set up for comparison under the same conditions: Scheme A (Prior Art 1: Traditional Unconstrained Method): Only the LOS observation geometry of the ascending and descending rails was used to establish simultaneous equations for three-dimensional decomposition and solution. Scheme B (Prior Art 2: Traditional Hard Constraint Method): The assumption of surface parallel flow (SPF) was introduced, and this terrain constraint was used as a fixed hard constraint (i.e., the deformation was strictly assumed to occur along the ideal slope, and residuals were not allowed to exist). Scheme C (Method of this Application): Terrain parameters were extracted and terrain constraint applicability indexes were calculated. A solution model was jointly constructed using soft constraints, allowing residuals to exist and adaptively allocating constraint weights. At the same time, ill-conditioning was performed to ensure the stability of the model solution.

[0081] (3) Performance data test and comparison results. The inversion results of the three schemes were compared with the set true values ​​of deformation, and the root mean square error (RMSE) of the deformation rate in each dimension was calculated. The performance test data of key indicators are shown in Table 1: Table 1

[0082] (4) Comparison and analysis of results, combining experimental data and figures ( Figure 7 , Figure 8 Analysis reveals the following: ① Regarding the insensitivity to north-south deformation: Due to the limitations of SAR satellite observation geometry, the model condition number of Scheme A (unconstrained) is as high as 135.2, resulting in a severely ill-conditioned solution and extremely high north-south (N-direction) deformation error (RMSE reaches 32.6 mm / yr), making it completely impossible to recover the true north-south sliding characteristics. ② Regarding the additional bias problem of traditional hard constraints: Although Scheme B significantly reduces the condition number and improves the accuracy in the north-south direction, the simulated landslide body has local deformation components that are not strictly parallel to the ideal slope. Forcibly intervening with topographic priors as "hard constraints" causes the propagation and transfer of errors, resulting in a severe degradation of accuracy in the east-west (E-direction) and vertical (U-direction) directions (E-direction RMSE deteriorates from 4.5 to 8.7 mm / yr). ③ The unexpected technical effect of this application: This application (Scheme C) effectively balances the weights of InSAR observation and topographic priors by introducing a "topographic constraint applicability evaluation" and a "soft constraint weighting" mechanism. The results show that this method not only effectively overcomes the ill-conditioned nature of the model, successfully reducing the north-south RMSE to 6.5 mm / yr (an improvement of approximately 80% in accuracy compared to the unconstrained scheme), but also avoids the additional biases caused by hard constraints, maintaining optimal accuracy in both the east-west and vertical directions. This simulation fully demonstrates the extremely high stability and reliability of this application in solving complex mountain three-dimensional deformation problems.

[0083] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for joint calculation of three-dimensional deformation in temporal InSAR with soft-constraint terrain, characterized in that, Includes the following steps: S10: Acquire multi-temporal SAR images, digital elevation model data, and satellite imaging geometric parameters of the study area; S20, perform time-series InSAR processing on the multi-temporal SAR images to obtain the satellite line-of-sight cumulative deformation time series and InSAR observation quality information at each observation time; S30, extract at least one topographic parameter from the slope, aspect, surface normal, topographic curvature, or main direction of motion of the study area based on the digital elevation model data; S40, construct terrain constraints based on the terrain parameters, the terrain constraints including at least one of surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints; construct terrain constraint applicability indices based on the terrain parameters and the InSAR observation quality information, and determine the constraint weights of the terrain constraints based on the terrain constraint applicability indices; S50, establish the line-of-sight observation relationship between the cumulative deformation along the satellite line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component based on the satellite imaging geometric parameters, and establish time continuity constraints for the three-dimensional deformation components or three-dimensional deformation rate at continuous observation times; combine the line-of-sight observation relationship, the time continuity constraints, and the terrain constraints in the form of soft constraints that allow the existence of constraint residuals to construct a time-series-oriented three-dimensional deformation joint solution model; S60, determine the line-of-sight observation weight based on the InSAR observation quality information, use the constraint weight as the terrain constraint weight, and combine it with the time continuity constraint to determine the time continuity constraint weight; construct a joint weight matrix including the line-of-sight observation weight, terrain constraint weight, and time continuity constraint weight, and solve the time-series-oriented three-dimensional deformation joint solution model to obtain the east-west deformation component, north-south deformation component, and vertical deformation component at each observation time.

2. The method for joint calculation of three-dimensional deformation in terrain-constrained temporal InSAR as described in claim 1, characterized in that, The time-series InSAR processing includes at least one of image registration, interferometric network construction, differential interferometric processing, terrain phase removal, phase filtering, phase unwrapping, atmospheric delay correction, reference point correction, and time-series deformation inversion; the terrain parameters include at least one of slope, aspect, surface normal vector, terrain curvature, local aspect dispersion, geomorphic unit boundary, main slip direction, glacier flow direction, and geological structure direction.

3. The method for joint calculation of three-dimensional deformation in temporal InSAR with terrain soft constraints as described in claim 1, characterized in that, When the terrain constraint is a surface parallel flow constraint, the surface parallel flow constraint is expressed in the form of a soft constraint as follows: ,in This represents the constraint function constructed based on the surface parallel flow assumption. , , These represent the east-west deformation component, the north-south deformation component, and the vertical deformation component, respectively. This represents the surface parallel flow constraint residual.

4. The method for joint calculation of three-dimensional deformation in temporal InSAR with soft-constraint terrain as described in claim 1, characterized in that, The soft constraint form refers to allowing the existence of terrain constraint residuals in the joint solution model, and weighting the terrain constraint residuals according to the terrain constraint applicability index and the InSAR observation quality.

5. The method for joint calculation of three-dimensional deformation in temporal InSAR with soft-constraint terrain as described in claim 1, characterized in that, The terrain constraint applicability index The terrain constraint applicability index is determined based on at least one of slope, aspect local dispersion, terrain curvature, InSAR coherence coefficient, or phase unwrapping quality; the calculation formula for the terrain constraint applicability index is as follows: ,in Indicates slope; This represents slope aspect stability, which is calculated from the slope aspect dispersion within a local window. Indicates the curvature of the terrain; Indicates the InSAR coherence coefficient; This indicates the phase unwrapping quality.

6. The method for joint calculation of three-dimensional deformation in terrain-constrained temporal InSAR as described in claim 1, characterized in that, The joint weight matrix is ​​expressed as follows: ,in The line-of-sight observation weight is determined based on at least one of the InSAR coherence coefficient, phase standard deviation, unwrapping residual, and atmospheric delay correction quality. The terrain constraint weight is determined based on at least one of the following: slope, aspect stability, terrain curvature, geomorphic unit consistency, and confidence of the main movement direction. The weights representing time continuity constraints are determined based on at least one of time interval, degree of deformation nonlinearity, and time-series residuals. The regularization weight is determined based on at least one of the following: model ill-conditionedness, noise level, spatial smoothing requirements, and temporal stability requirements.

7. The method for joint calculation of three-dimensional deformation in temporal InSAR with terrain soft constraints as described in claim 1, characterized in that, The three-dimensional deformation joint solution model is solved using the regularized least squares method, and the solution formula is as follows: in, This represents the three-dimensional deformation vector to be determined, in a single-stage solution mode. , , , These represent the east-west, north-south, and vertical deformation components, respectively; in the time-series joint solution mode... It is an augmented vector containing the three-dimensional deformation components at all observation times; Represents the joint weight matrix; Represents the joint design matrix; Represents the observation vector; Represents the regularization parameter; This represents a spatial smoothing, temporal smoothing, or prior constraint matrix.

8. The method for joint calculation of three-dimensional deformation in terrain-constrained temporal InSAR as described in claim 1, characterized in that, During or after solving the three-dimensional deformation joint solution model, an ill-conditioning diagnosis is performed. The ill-conditioning diagnosis includes at least one of the following: joint design matrix condition number, parameter covariance matrix, posterior unit weight variance, line-of-sight residual, terrain constraint residual, and temporal deformation jump degree. When the ill-conditioning diagnosis index exceeds a preset threshold, the terrain constraint weight is adjusted, a regularization term is added, the reliability level of the three-dimensional deformation result is reduced, or the corresponding region is marked as an unsuitable three-dimensional solution region.

9. The method for joint calculation of three-dimensional deformation in temporal InSAR with soft-constraint terrain as described in claim 1, characterized in that, Temporal continuity constraints or piecewise linear constraints are introduced into the three-dimensional deformation components or three-dimensional deformation rates at continuous observation times to suppress anomalous jumps in three-dimensional deformation caused by single-period InSAR noise.

10. A joint solution system for three-dimensional deformation of terrain-constrained temporal InSAR, characterized in that, To implement the method of any one of claims 1 to 9, comprising: The data acquisition module is used to acquire multi-temporal SAR images, digital elevation model data, and satellite imaging geometric parameters of the study area; The temporal InSAR processing module is used to perform temporal InSAR processing on multi-temporal SAR images to obtain the satellite line-of-sight cumulative deformation time series and InSAR observation quality information at each observation time. The line-of-sight observation modeling module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite's line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component, based on the satellite imaging geometric parameters. The terrain parameter extraction module is used to extract at least one terrain parameter from the slope, aspect, surface normal vector, terrain curvature, or main direction of motion of the study area based on the digital elevation model data. The terrain constraint construction module is used to construct terrain constraints based on the terrain parameters. The terrain constraints include at least one of surface parallel flow constraints, slope aspect constraints, surface tangent plane constraints, or main motion direction constraints. Then, the module constructs a terrain constraint applicability index based on the terrain parameters and the InSAR observation quality information, and determines the constraint weight of the terrain constraint based on the terrain constraint applicability index. The model building module is used to establish the line-of-sight observation relationship between the cumulative deformation along the satellite line of sight and the east-west deformation component, the north-south deformation component, and the vertical deformation component based on the satellite imaging geometric parameters, and to establish time continuity constraints on the three-dimensional deformation components or three-dimensional deformation rate at continuous observation times; the line-of-sight observation relationship, the time continuity constraints, and the terrain constraints are combined in the form of soft constraints that allow the existence of constraint residuals to construct a time-series-oriented three-dimensional deformation joint solution model. The joint solution module is used to determine the line-of-sight observation weights based on the InSAR observation quality information, use the constraint weights as terrain constraint weights, and combine them with the time continuity constraint to determine the time continuity constraint weights; construct a joint weight matrix containing the line-of-sight observation weights, terrain constraint weights, and time continuity constraint weights, and solve the time-series-oriented three-dimensional deformation joint solution model to obtain the east-west deformation components, north-south deformation components, and vertical deformation components at each observation time; The quality evaluation module is used to perform pathological diagnosis, residual evaluation, and reliability classification on the three-dimensional deformation solution results; The results output module is used to output three-dimensional deformation time series, three-dimensional deformation rate field, deformation results along slope, and reliability evaluation results.

Images