Surface deformation inversion method based on splicing InSAR deformation field in overlapping area

Abstract

The present application provides a surface deformation inversion method based on overlapping area splicing InSAR deformation field, first, obtaining multiple frames of LOS images, image frame information and GNSS observation data within the scope of the study area; then classifying the LOS images to obtain the same track data set, and then screening the overlapping area to obtain the same track overlapping area and the corresponding GNSS station data; constructing the observation equation to solve and eliminate the error of the same track data, and splicing the same track data; performing spatial interpolation and projection on the GNSS three-dimensional deformation observation data to obtain the GNSS three-dimensional deformation projection data; combining the GNSS three-dimensional deformation projection data to perform weighted leveling of the error plane and execute the adjacent track data splicing; performing surface three-dimensional deformation inversion to obtain surface three-dimensional deformation data and calculate the strain rate data, and obtaining the evaluation result of surface three-dimensional deformation inversion. The method improves the precision of surface three-dimensional deformation inversion.

Description

Technical Field

[0001] This invention relates to the fields of geophysics and remote sensing image processing technology, specifically to a method for surface deformation inversion based on InSAR deformation fields stitched together from overlapping regions. Background Technology

[0002] Geological hazards are complex processes resulting from the combined effects of geodynamics and environmental factors, involving various deformation events, including but not limited to earthquakes, volcanic activity, landslides, ground subsidence, debris flows, land settlement, and freeze-thaw deformation caused by seasonal cycles. Regardless of the type of hazard, their development is typically accompanied by long-term or short-term strain accumulation, energy release, and stress redistribution. These processes exhibit significant non-uniformity and multi-scale coupling across temporal and spatial scales, with different mechanisms interacting through stress transfer, topographic deformation, and changes in the subsurface medium. To achieve effective risk assessment, early warning response, and post-disaster recovery planning, continuous, high-resolution reconstruction and monitoring of three-dimensional surface deformation over a wide area are necessary. This reconstruction can reveal the spatiotemporal evolution characteristics of different hazard processes, help infer the geometry of potential faults and landslide zones, the mechanical behavior of soil and rock masses, and the infiltration pathways of subsurface fluids, and explore driving mechanisms, ultimately supporting scenario-based hazard simulation and engineering decision-making.

[0003] Interferometric Synthetic Aperture Radar (InSAR) observation technology, with its advantages of high spatial resolution, wide coverage, and relative insensitivity to weather conditions, has become a key means of monitoring surface deformation. Traditional InSAR observations can only provide deformation information along the radar line of sight (LOS). To obtain a three-dimensional deformation field, it is usually necessary to combine multi-orbit, multi-view observation data, or introduce other observation methods, such as Global Navigation Satellite System (GNSS), pixel offset tracking (POT), Multiple Aperture Interferometry (MAI), or surface deformation models for constraint solving.

[0004] However, the problem becomes more complex in deformation inversion over large regions and long time scales. Obtaining extensive 3D deformation information requires integrating observational data from numerous different geographical regions, orbits, and time windows. Because the geometric relationships between different observation frames are not entirely consistent, unstable step phenomena often occur in overlapping regions. For example, observations from overlapping regions of frames from the same orbit may contain discrepancies, while irregular steps may appear in overlapping regions of frames from different orbits. These steps directly affect subsequent 3D inversion results, reducing inversion stability and accuracy, and impacting its spatiotemporal consistency.

[0005] The key to inverting large-area surface deformation using InSAR observation data lies in stitching different observation frames into a complete dataset. The core of the stitching operation is to unify the observation frames of different observation frames under a fixed reference frame, thereby reducing the geometric projection error caused by frame differences and suppressing the influence of long wavelengths. Currently, there are two main strategies for inverting large-area surface three-dimensional deformation: (1) Before calculating the deformation of each radar line of sight (LOS direction), the original observation images are processed using a long orbit window. This can obtain a series of continuous observation strips, each representing the LOS direction deformation information of an orbit. Subsequently, the data of the ascending orbit and the descending orbit are processed separately, and their data in the overlapping area are used, combined with the satellite's incident angle and azimuth angle, to solve for the large-area surface three-dimensional deformation. (2) Multiple frames of LOS direction deformation original observation frames that have been standardized within the study area are stitched together to cover the entire study area. Then, based on the data types of the ascending and descending orbits, the LOS-oriented deformation images are stitched together, and GNSS data or multi-source satellite data are used to invert the three-dimensional surface deformation of a large area.

[0006] Existing technologies have solved the problem of large-scale surface deformation inversion to a certain extent, but challenges such as resource constraints, accuracy issues, and process complexity still exist, which limit their application to a larger scale, higher accuracy, and greater efficiency. In general, there are three main shortcomings: (1) Computational resources and economic costs: Processing large-scale long-orbit InSAR datasets requires a lot of computing time and storage space, which not only increases hardware costs but also makes it difficult to carry out research work with limited resources. (2) Data consistency and accuracy issues: When using multiple images to cover the research area, there are challenges in ensuring image consistency and accuracy between frames. In the processing of multiple orbits and multiple observation frames, due to the influence of factors such as frame differences and atmospheric disturbances, the deviations and inconsistencies in the overlapping areas are often difficult to completely eliminate, which directly affects the reliability of the inversion results. (3) Complexity of the stitching process: Stitching images requires processing both along-orbit and cross-orbit data at the same time, which means that different geometric parameters need to be coordinated; secondly, the fusion of multi-source satellite data increases the complexity of the process, and geometric corrections between different sensors need to be considered; in addition, how to achieve a balance between geometric constraints and external reference data to optimize the stitching results is also a technical problem. Summary of the Invention

[0007] The purpose of this invention is to provide a method for surface deformation inversion based on InSAR deformation field stitching in overlapping regions. This method enables efficient stitching and unification of LOS images from large-area, multi-orbit InSAR observation data under limited computational and storage conditions. Furthermore, it allows for the inversion of surface 3D deformation and the calculation of strain rate data based on the stitched results combined with spatially interpolated GNSS 3D deformation observation data. This method significantly reduces computational and storage costs by combining GNSS 3D deformation observation data from both inside and outside orbits, employing a step-by-step solution and balancing strategy. Simultaneously, it improves the consistency of the deformation field in the stitched LOS image and the accuracy of subsequent inversion.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows.

[0009] The surface deformation inversion method based on overlapping region mosaicking InSAR deformation field includes the following steps:

[0010] Step S1: Acquire multiple frames of LOS images within the study area, the corresponding image frame information, and GNSS observation data. The GNSS observation data includes observation data, observation uncertainty, GNSS three-dimensional deformation observation data, and GNSS station data. The image frame information includes the coordinates, orbit type, and orbit number of the LOS images.

[0011] Step S2: Classify LOS images according to orbit type and orbit number. Group LOS images with the same orbit type and orbit number into a set of co-orbit data to obtain a co-orbit data set. Then, based on the co-orbit data set, filter out the overlapping areas of two geographically adjacent LOS images in the co-orbit data to obtain the overlapping areas of the same orbit. Finally, obtain GNSS station data and GNSS three-dimensional deformation observation data in the overlapping areas of the same orbit.

[0012] Step S3: Based on the GNSS station data and GNSS 3D deformation observation data in the overlapping area of ​​the same track, construct the observation equation, use the observation equation to solve the error plane between two geographically adjacent LOS images in the same track data, eliminate the error between adjacent LOS images and perform same track stitching to obtain a long strip LOS image stitched along the track direction.

[0013] Step S4: Spatial interpolation is performed on the GNSS three-dimensional deformation observation data, and then the spatially interpolated GNSS three-dimensional deformation observation data is projected onto the orbital plane where each corresponding LOS image is located to obtain GNSS three-dimensional deformation projection data.

[0014] Step S5: Obtain the strip-shaped LOS images that have been stitched along the track direction for two geographically adjacent tracks to form adjacent track data. Combine the data with GNSS 3D deformation projection data to solve the error plane between the two strip-shaped LOS images of adjacent tracks. Perform weighted balancing of the error plane on each strip-shaped LOS image to obtain the corrected strip-shaped LOS image. Then, stitch the images together with the adjacent tracks to obtain the stitched LOS image covering the study area.

[0015] Step S6: Combine the stitched LOS imagery covering the study area with the spatially interpolated GNSS 3D deformation observation data to construct the surface 3D deformation inversion equation, perform surface 3D deformation inversion, and obtain surface 3D deformation data.

[0016] Step S7: Calculate strain rate data based on the three-dimensional deformation data of the earth's surface to obtain the evaluation results of the three-dimensional deformation inversion of the earth's surface.

[0017] Furthermore, in step S1, the GNSS three-dimensional deformation observation data includes GNSS east-west deformation data, GNSS north-south deformation data, and GNSS vertical deformation data; the GNSS station data includes the number of GNSS stations and GNSS station location information; and the orbit type includes ascending orbit type and descending orbit type.

[0018] Furthermore, step S3 specifically includes the following steps:

[0019] Step S31, define the threshold for the number of GNSS stations as M. GNSS ;

[0020] Step S32: Based on the GNSS station data in the overlapping area of ​​two geographically adjacent LOS images in the same orbit, determine whether the number of GNSS stations in the overlapping area is greater than the threshold M. GNSS ;

[0021] If the number of GNSS stations within the overlapping area of ​​two geographically adjacent LOS images in the same track data is less than the threshold M GNSS Then, the observation equation is directly constructed based on the coordinates and observation data of the two LOS images;

[0022] If the number of GNSS stations within the overlapping area of ​​two geographically adjacent LOS images in the same track data is greater than or equal to the threshold M, then... GNSS Then, the coordinates and observation data of the two LOS images and the GNSS three-dimensional deformation observation data are combined to construct the observation equation;

[0023] The error plane between two geographically adjacent LOS images in the same track data is solved by using the observation equation, thus eliminating the error between adjacent LOS images.

[0024] Step S33: Stitch together the LOS images of the same track after error elimination to obtain a long strip LOS image stitched along the track direction.

[0025] Furthermore, in step S32, the observation equation directly constructed based on the coordinates and observation data of the two LOS images is as follows:

[0026] ,

[0027] In the formula, This represents the error plane between two geographically adjacent LOS images in the same orbital data. and These represent the observation data of two geographically adjacent LOS images in the same orbital data; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. , and It is the model factor of the error plane between two geographically adjacent LOS images in the same track data;

[0028] The specific steps for eliminating errors between adjacent LOS images by using the observation equation to solve for the error plane between two geographically adjacent LOS images in the same orbital data are as follows:

[0029] The model factors were obtained by solving the observation equation. , and Next, the two LOS images are defined as a single entity, and then the formula is used. To calculate the observation data used to eliminate the error plane between two LOS images, where, Represents the overall X-axis coordinates of the two LOS images. This represents the overall Y-axis coordinate of the two LOS images. This represents the error plane between two geographically adjacent LOS images in the same orbital data, expressed in coordinates. The observed data below;

[0030] Then, the error plane between two geographically adjacent LOS images in the same orbital data is evenly distributed to eliminate the error between the two LOS images; for the two LOS images, the observed data after eliminating the error plane are respectively... and .

[0031] Furthermore, in step S32, the observation equation constructed by combining the coordinates and observation data of the two LOS images with the GNSS three-dimensional deformation observation data is as follows:

[0032] ,

[0033] In the formula, and These represent the observation data of two geographically adjacent LOS images in the same orbital data; This represents the east-west deformation data of GNSS in 3D deformation observation data. This represents the north-south oriented deformation data in GNSS three-dimensional deformation observation data. This represents the GNSS vertical deformation data in GNSS three-dimensional deformation observation data; and These represent the satellite observation azimuth angles corresponding to two geographically adjacent LOS images in the same orbital data; and These represent the satellite observation incident angles corresponding to two geographically adjacent LOS images in the same orbital data; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. The east-west deformation of the t-th pixel within the overlapping region of two LOS images is represented by this value. The north-south deformation of the t-th pixel within the overlapping region of two LOS images represents the deformation of the Lt-th pixel. The vertical deformation of the t-th pixel within the overlapping region of two LOS images; , , and , , These represent the model factors of the error plane between two geographically adjacent LOS images and the GNSS 3D deformation observation data plane in the same track data;

[0034] The specific steps for eliminating errors between adjacent LOS images by using the observation equation to solve for the error plane between two geographically adjacent LOS images in the same orbital data are as follows:

[0035] The model factors were obtained by solving the observation equation. , , and , , Then, using the GNSS 3D deformation observation data plane as the reference plane for error elimination, the error planes between the two LOS images and the GNSS 3D deformation observation data plane are eliminated respectively, thereby eliminating the error planes between the two LOS images; specifically:

[0036] Using formula and formula Calculate the observation data for the error plane used to eliminate the gap between the two LOS imagery planes and the GNSS 3D deformation observation data plane; where, and These represent the X-axis and Y-axis coordinates of one of the LOS images, respectively. This indicates the error plane between the corresponding LOS imagery and the GNSS 3D deformation observation data plane in coordinates. The observed data below; and These represent the X-axis and Y-axis coordinates of another LOS image, respectively. This indicates the error plane between the corresponding LOS imagery and the GNSS 3D deformation observation data plane in coordinates. The observed data below;

[0037] For the two LOS images, the observed data after eliminating the error plane are respectively and .

[0038] Furthermore, step S4 specifically includes the following steps:

[0039] Step S41: Based on the uncertainty of the observed values, the GNSS three-dimensional deformation observation data within the study area is screened and then spatially interpolated to obtain spatially interpolated GNSS three-dimensional deformation observation data, including spatially interpolated GNSS east-west deformation data. Spatially interpolated GNSS north-south deformation data GNSS vertical deformation data after spatial interpolation ;

[0040] Step S42: Then, the spatially interpolated GNSS 3D deformation observation data is projected onto the orbital plane of each corresponding LOS image to obtain GNSS 3D deformation projection data. The projection formula is expressed as follows:

[0041] ,

[0042] In the formula, The data represents the spatially interpolated GNSS 3D deformation observation data projected onto the orbital plane of the LOS image, i.e., GNSS 3D deformation projection data; and These correspond to the azimuth and angle of incidence of the satellite's observation data during flight, respectively. This represents the spatially interpolated GNSS east-west deformation data. This represents the spatially interpolated GNSS north-south deformation data. This represents the spatially interpolated GNSS vertical deformation data.

[0043] Furthermore, step S5 specifically includes the following steps:

[0044] Step S51: Obtain the strip-shaped LOS images stitched along the track direction corresponding to two geographically adjacent tracks to form adjacent track data. Combine this with GNSS 3D deformation projection data to solve the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformation projection data plane. Determine the model factor of the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformation projection data plane. The calculation formula is as follows:

[0045] ,

[0046] In the formula, This represents the error plane between two long strip-shaped LOS images of adjacent orbits. This represents the error plane between one of the elongated LOS imagery and the GNSS 3D deformation projection data plane. This represents the error plane between another elongated LOS image and the GNSS 3D deformation projection data plane; and The observation data represent two strip-shaped LOS images of adjacent orbits, respectively; and These represent the GNSS 3D deformation projection data of two geographically adjacent tracks in the adjacent track data, respectively, and correspond to the GNSS 3D deformation projection data of two long strip-shaped LOS images of the adjacent tracks; and The X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of two long strip-shaped LOS images representing adjacent tracks; , , and , , The model factors represent the error planes between two strip-shaped LOS images of adjacent orbits and the GNSS three-dimensional deformation projection data plane, respectively.

[0047] The model factors were obtained by solving. , , and , , Then, using the GNSS 3D deformation projection data plane as the reference plane for error elimination, the error planes between the two strip-shaped LOS images and the GNSS 3D deformation projection data plane are eliminated respectively, thereby eliminating the error planes between the two strip-shaped LOS images; specifically:

[0048] Using formula and formula Calculate the observation data for the error plane used to eliminate the gap between the two strip-shaped LOS images and the GNSS 3D deformation projection data plane; where, and These represent the X-axis and Y-axis coordinates of one of the elongated LOS images, respectively. This indicates the error plane between the corresponding elongated LOS image and the GNSS 3D deformation projection data plane in coordinates. The observed data below; and These represent the X-axis and Y-axis coordinates of another elongated LOS image, respectively. This indicates the error plane between the corresponding elongated LOS image and the GNSS 3D deformation projection data plane in coordinates. The observed data below;

[0049] For two strip-shaped LOS images, after eliminating the error plane between the strip-shaped LOS images and the GNSS 3D deformation projection data plane, the resulting observation data are as follows: and ; and These represent the observation data of the two strip-shaped LOS images obtained after eliminating the error plane between the strip-shaped LOS image and the GNSS three-dimensional deformation projection data plane;

[0050] Step S52: Using the overlapping region data between two strip-shaped LOS images of adjacent tracks, solve for the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks, and determine the model factor of the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks. The calculation formula is as follows:

[0051] ,

[0052] In the formula, and The observation data represent two strip-shaped LOS images of adjacent orbits, respectively; and These represent the GNSS 3D deformation projection data of two geographically adjacent tracks in the adjacent track data, corresponding to the GNSS 3D deformation projection data of two long strip-shaped LOS images; and The X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of two long strip-shaped LOS images of adjacent tracks are respectively represented. , and The model factor representing the error plane between two adjacent LOS images caused by orbital errors;

[0053] The model factor of the error plane caused by orbital errors between two long strip-shaped LOS images of adjacent orbits is obtained by solving the problem. , and Next, the two elongated LOS images are defined as a single entity, and then the formula is used. To calculate the observed data for the error plane between two long strip LOS images of adjacent orbits, which is used to eliminate the error caused by orbital errors, where, The X-axis coordinates of the two elongated LOS images of adjacent orbits are represented. The Y-axis coordinate represents the overall coordinates of two long strip-shaped LOS images of adjacent orbits. The error plane representing the coordinate system between two long strip-shaped LOS images of adjacent orbits. The observation data below is used to eliminate the error plane caused by orbital errors between two long strip LOS images of adjacent orbits;

[0054] The error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is solved separately. The model factor of the error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is determined. The observed data used to eliminate the error plane caused by orbital error between the strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit are calculated separately.

[0055] The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its west is defined as the westward correction value of the strip-shaped LOS image, denoted as . The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its east of its orbit is defined as the eastward correction value of the strip-shaped LOS image, denoted as . ;

[0056] Then, based on the size of the area difference between the overlapping areas on the east and west sides of each strip-shaped LOS image track, weighted balancing of the error plane is performed on each strip-shaped LOS image. The specific steps are as follows:

[0057] Define the area of ​​the overlapping region on the eastern side of the elongated LOS image as S. 东 The area of ​​the overlapping region on the west side of the elongated LOS image is S. 西 ;

[0058] If the area difference between the overlapping regions on the east and west sides of the elongated LOS image is S 东 / S 西 ≤0.8 or S 东 / S 西 If the error plane correction is ≥1.2, then no error plane correction is performed on the east and west sides;

[0059] If the area difference between the overlapping regions on the east and west sides of the elongated LOS image is 0.8... 东 / S 西 <1.2, the correction value for weighted balancing of the error plane for each strip-shaped LOS image is defined as: The corrected observation data for the elongated LOS image are: ,in, This represents the observation data of the long strip-shaped LOS image after eliminating the error plane caused by overlapping regions;

[0060] ​Step S53: After weighted balancing of each strip-shaped LOS image with an error plane, a corrected strip-shaped LOS image is obtained; then, the corrected strip-shaped LOS images are stitched together according to the orbit type to obtain a stitched LOS image covering the study area, including a stitched ascending orbit LOS image covering the study area and a stitched descending orbit LOS image covering the study area.

[0061] Furthermore, in step S6, the three-dimensional deformation inversion equation of the Earth's surface is as follows:

[0062] ,

[0063] In the formula, The up-track LOS image representing the stitched-together coverage area of ​​the study region. The down-orbiting LOS image representing the stitched-together coverage area of ​​the study region. Represents spatially interpolated GNSS north-south deformation data; The satellite azimuth angle representing the orbital ascent type data. The satellite azimuth angle representing the down-orbit type data. Indicates the satellite incident angle for ascending orbit type data. Indicates the satellite incident angle for data of the down-orbit type; , and This represents the three-dimensional deformation rate field of the Earth's surface that needs to be inverted, where... Represents the vertical deformation rate field. Represents the north-south deformation rate field. The east-west deformation rate field represents the value that needs to be solved in the three-dimensional deformation inversion equation of the Earth's surface. Among these, the east-west deformation rate field... and vertical deformation rate field It is based on the stitched LOS imagery covering the research area. and LOS imagery The obtained north-south deformation rate field It is based on the spatially interpolated GNSS north-south deformation data The three sets of deformation rate field data obtained from the solution together constitute the three-dimensional deformation data of the Earth's surface.

[0064] Furthermore, step S7 specifically includes the following steps:

[0065] Step S71: Calculate the horizontal strain rate tensor based on the three-dimensional surface deformation data obtained in step S6. The calculation formula is as follows:

[0066] ,

[0067] In the formula, Represents the horizontal strain rate tensor. , , , This represents the velocity gradient in the horizontal direction, where... This represents the rate of change of the east-west velocity with respect to the X-axis coordinate. This represents the rate of change of the north-south velocity with respect to the Y-axis coordinate. = , , Both are quantities with the same meaning, measuring the average effect of the change in east-west velocity in the north-south direction and the change in north-south velocity in the east-west direction, reflecting the shear deformation of the object; x represents the distance between two adjacent pixels in the east-west direction in the horizontal direction, and y represents the distance between two adjacent pixels in the north-south direction in the horizontal direction. This represents the partial derivative of the east-west deformation rate with respect to the X-axis. This represents the partial derivative of the east-west deformation rate with respect to the Y-axis. This represents the partial derivative of the north-south deformation rate with respect to the X-axis. The partial derivative of the north-south deformation rate with respect to the Y-axis is given.

[0068] Step S72: Based on the velocity gradients in each horizontal direction of the horizontal strain rate tensor obtained in step S71, calculate the horizontal expansion rate, shear rate, and the second-order invariants of the horizontal strain rate tensor. The calculation formulas are as follows:

[0069] ,

[0070] In the formula, Represents the horizontal expansion rate. Represents the shear rate. The second-order invariant representing the horizontal strain rate tensor;

[0071] Second-order invariant data of horizontal expansion rate, shear rate, and horizontal strain rate tensors can be used for crustal deformation analysis and geological engineering assessment. A larger positive value for the horizontal expansion rate indicates stronger stretching and a significant surface expansion trend, corresponding to normal fault zones and crustal extension zones. A larger negative value for the horizontal expansion rate indicates stronger compression and a significant surface shortening trend, corresponding to thrust fault zones and plate compression zones. A larger shear rate value indicates stronger shearing and more intense relative slippage of landmasses, corresponding to strike-slip fault zones and landmass boundary slip zones. A larger second-order invariant value for the horizontal strain rate tensor indicates greater overall deformation intensity and more intense crustal activity, corresponding to fault activity zones and high-risk earthquake zones.

[0072] The method of the present invention has the following technical effects:

[0073] This invention uses overlapping regions as correction units, combining GNSS 3D deformation projection data, Kriging interpolation, and layered (same-track priority followed by cross-track) error plane modeling and design weighting and balancing strategies. This achieves consistent stitching of multi-track LOS imagery with limited computing resources, thus stably and efficiently supporting the inversion of 3D surface deformation over large areas. Specifically, it includes:

[0074] (1) Improved consistency of deformation rate field of cross-frame and cross-track LOS images: By solving and correcting the error plane of adjacent frames in the same track, and then using overlapping area data and GNSS three-dimensional deformation projection data to constrain the error plane between adjacent tracks, the step and system deviation between frames and tracks are significantly reduced, resulting in a more continuous and smooth long track and cross-track LOS image set.

[0075] (2) Improved accuracy and stability of three-dimensional deformation inversion: By binding InSAR observation data, i.e. LOS image, to the reference frame of GNSS three-dimensional deformation observation data and using GNSS three-dimensional deformation projection data as an absolute constraint, the errors caused by baseline drift and inconsistency of reference surface are reduced, thus making the three-dimensional deformation rate field constrained and inverted by the joint ascending and descending orbit LOS image and GNSS three-dimensional deformation projection data more reliable and numerically stable.

[0076] (3) The adaptive fusion strategy reduces the risk of weak constraints: the solution method is selected according to the number of GNSS stations in the overlapping area (the observation consistency constraint is used when the GNSS station distribution is sparse, and the GNSS-InSAR joint constraint is used when the GNSS station distribution is sufficient), so as to avoid overfitting or solution instability in the case of sparse GNSS station distribution, and to balance robustness and accuracy.

[0077] (4) Reduced the impact of observation errors on subsequent geological interpretation: By controlling and balancing the inter-orbit error plane as a whole and using weight allocation in the overlapping area, the observation errors caused by orbital geometry, time baseline or atmospheric residuals can be effectively suppressed, thereby reducing the pseudo signals of strain rate and deformation rate fields and improving the reliability of geological interpretation.

[0078] (5) Significantly reduced processing costs under resource constraints: effectively reduced step and inconsistency in overlapping areas, and improved the intrinsic consistency of deformation rate field of LOS image after stitching; achieved a balance between geometric constraints and external references, and provided a continuous, reliable and accurate input dataset for large-scale regional inter-seismic surface 3D deformation inversion, thereby improving the accuracy of subsequent surface 3D deformation inversion; supported the fusion of multi-track and multi-source data, and adapted to the inversion needs of large-scale regions. Attached Figure Description

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

[0080] Figure 1 This is an overall flowchart of the surface deformation inversion method based on overlapping region stitching InSAR deformation field in this embodiment of the invention;

[0081] Figure 2 This is a spatial distribution map of the LOS imagery in an embodiment of the present invention;

[0082] Figure 3 This is a schematic diagram illustrating the principle of data stitching of LOS images in an embodiment of the present invention;

[0083] Figure 4 This refers to the unstitched LOS image provided in this embodiment of the invention;

[0084] Figure 5 This refers to the unstitched LOS image provided in this embodiment of the invention;

[0085] Figure 6 This is a result image of the Kriging interpolation GNSS three-dimensional deformation observation data provided in an embodiment of the present invention;

[0086] Figure 7 This refers to the stitched LOS image covering the research area provided in this embodiment of the invention;

[0087] Figure 8 This refers to the stitched LOS image covering the research area provided in this embodiment of the invention;

[0088] Figure 9 This is an example of an embodiment of the invention providing east-west and vertical deformation rate field maps derived from LOS image data;

[0089] Figure 10 This is a horizontal expansion rate distribution diagram provided in an embodiment of the present invention;

[0090] Figure 11 This is a shear rate distribution diagram provided in an embodiment of the present invention;

[0091] Figure 12 This is a second-order invariant plot of the horizontal strain rate tensor provided in this embodiment of the invention;

[0092] Figure 13 This is a comparison chart of the root mean square difference between the splicing results of the method of this invention and the traditional method of splicing only along the track. Detailed Implementation

[0093] To better understand the above-described objects, features, and advantages of the present invention, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the invention; however, the invention may be practiced in other ways different from those described herein, and therefore, the invention is not limited to the specific embodiments disclosed below.

[0094] To facilitate understanding of this invention, some technical terms specific to this art will be introduced first:

[0095] (1) The preliminary product obtained by processing satellite observation images through InSAR technology is called the LOS orientation deformation image of the InSAR image observation frame, or simply "LOS image".

[0096] (2) During the process of satellite observation of the Earth, the satellite's flight orbit can be divided into ascending orbit and descending orbit due to different flight directions. Among them, the orbit of the satellite flying from south to north and far away from the Earth's surface is called "ascending orbit", and the orbit of the satellite flying from north to south and close to the Earth's surface is called "descending orbit".

[0097] (3) To observe the entire Earth's surface, satellites need observation data from multiple geographically adjacent orbits to cover the entire Earth, as each ascent or descent orbit covers a relatively small area. If LOS images are in the same ascent or descent orbit, these LOS images are called "co-orbit data," while LOS image data from geographically adjacent orbits are called "adjacent orbit data."

[0098] (4) GNSS observation stations use GNSS observation equipment to receive product data processed from navigation satellite data, including GNSS east-west deformation data, GNSS north-south deformation data and GNSS vertical deformation data, collectively referred to as "GNSS three-dimensional deformation observation data".

[0099] Example:

[0100] like Figures 1-3 As shown, the surface deformation inversion method based on overlapping region stitched InSAR deformation field of the present invention includes the following steps:

[0101] Step S1, Obtain observation data:

[0102] Acquire multiple frames of LOS imagery, corresponding image frame information, and GNSS observation data within the study area. The GNSS observation data includes observation data, observation uncertainty, GNSS 3D deformation observation data, and GNSS station data. The image frame information includes the coordinates, orbit type, and orbit number of each LOS image.

[0103] The GNSS three-dimensional deformation observation data includes GNSS east-west deformation data, GNSS north-south deformation data, and GNSS vertical deformation data; the GNSS station data includes the number of GNSS stations (i.e., the number of GNSS observation stations) and GNSS station location information (i.e., the location information of GNSS observation stations); the orbit types include ascending orbit type and descending orbit type.

[0104] The observation data in this embodiment of the invention can be obtained from publicly available datasets such as the China Earthquake Science Data Sharing Center and the Nevada Geodetic Laboratory website.

[0105] Step S2, data preprocessing before splicing:

[0106] LOS images are classified according to orbit type and orbit number. LOS images with the same orbit type and orbit number are grouped into a co-orbit data set. Then, based on the co-orbit data set, overlapping areas between two geographically adjacent LOS images within the co-orbit data are selected to obtain the co-orbit overlapping areas. GNSS station data and GNSS 3D deformation observation data within these overlapping areas are then acquired. The specific steps are as follows:

[0107] Step S21: Divide the LOS images to be stitched into ascending-track type data and descending-track type data according to the ascending-track type and descending-track type (i.e., divide the LOS images into ascending-track LOS images and descending-track LOS images); in the subsequent data stitching process, ascending-track type data is stitched with ascending-track type data, and descending-track type data is stitched with descending-track type data; in this embodiment, the unstitched ascending-track LOS images are as follows: Figure 4 As shown, the unstitched LOS imagery is as follows: Figure 5 As shown;

[0108] Step S22: There are multiple tracks in both the ascending and descending track data. The ascending and descending track data are classified according to the track number. LOS images with the same track type and track number are grouped into a set of same track data to obtain the same track data set, which is used to prepare for subsequent same track data stitching.

[0109] Step S23: View the LOS images in the same track data set according to geographical location. Multiple LOS images under the same track number form a strip of the same color. Figure 2(a) shows multiple north-south bands of the same color, composed of multiple polygonal image frames, belonging to the same orbit number within the same orbital dataset. For example, 026A, 128A, and 055A correspond to a north-south band of the ascending orbit type composed of multiple frames of LOS images, while 106D, 033D, 135D, and 062D correspond to a north-south band of the descending orbit type composed of multiple frames of LOS images. There are overlapping areas between multiple frames of LOS images, including overlapping areas within the same orbit (i.e., the overlapping areas of two geographically adjacent LOS images in the same orbital data, such as...). Figure 2 (b) Figure 2 (as shown in (d)) and the overlapping area between adjacent tracks (i.e., the overlapping area between two geographically adjacent LOS images under adjacent track numbers, such as...) Figure 2 (c) Figure 2 As shown in (e)); based on the same track data set, the overlapping areas are filtered out, and the number of GNSS stations and GNSS three-dimensional deformation observation data of the corresponding overlapping areas are obtained, so as to prepare for the use of GNSS three-dimensional deformation observation data in the subsequent same track data stitching.

[0110] Step S3, Data splicing on the same track:

[0111] Based on GNSS station data and GNSS 3D deformation observation data in the overlapping area of ​​the same orbit, an observation equation is constructed. This equation is then used to solve for the error plane between two geographically adjacent LOS images in the same orbit data. This process eliminates errors between adjacent LOS images, and they are then stitched together along the same orbit to obtain a long strip-shaped LOS image stitched along the orbital direction. The specific steps are as follows:

[0112] Step S31, define the threshold for the number of GNSS stations as M. GNSS Threshold M GNSS The number of GNSS stations must be greater than or equal to 6 to ensure that the observation equation has a unique solution. In this embodiment, the threshold M for the number of GNSS stations is... GNSS There are 10;

[0113] Step S32: Based on the GNSS station data in the overlapping area of ​​two geographically adjacent LOS images in the same orbit, determine whether the number of GNSS stations in the overlapping area is greater than the threshold M. GNSS If the number of GNSS stations in the overlapping area of ​​two geographically adjacent LOS images in the same track data is less than the threshold M GNSS Then, the observation equation is directly constructed based on the coordinates and observation data of the two LOS images; if the number of GNSS stations in the overlapping area of ​​two geographically adjacent LOS images in the same orbit data is greater than or equal to the threshold M GNSSThen, the coordinates and observation data of the two LOS images are combined with the GNSS 3D deformation observation data to construct an observation equation; the observation equation is used to solve for the error plane between two geographically adjacent LOS images in the same orbital data, thus eliminating the error between adjacent LOS images; the specific steps are as follows:

[0114] Step S321: If the number of GNSS stations in the overlapping area of ​​two geographically adjacent LOS images (i.e., LOS image A11 and LOS image A12) in the same track data is less than 10, then the observation equation is directly constructed based on the coordinates and observation data of the two LOS images. The error plane between the two geographically adjacent LOS images in the same track data is solved using the observation equation, and the model factor of the error plane between the two geographically adjacent LOS images in the same track data is determined.

[0115] The observation equations directly constructed based on two LOS images and image frame information are as follows:

[0116] ,

[0117] The above observation equation, after transformation, becomes:

[0118] ,

[0119] In the formula, This represents the error plane between two geographically adjacent LOS images in the same track data (i.e., the error plane between LOS image A11 and LOS image A12). and These represent the observation data of two geographically adjacent LOS images in the same orbital data, namely... The observation data represents LOS image A11. Observational data representing LOS image A12; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. , and It is the model factor of the error plane between two geographically adjacent LOS images in the same track data;

[0120] The model factor was obtained by solving the above observation equation. , and Next, the two LOS images are defined as a single entity, and then the formula is used. To calculate the observation data used to eliminate the error plane between two LOS images, where, Represents the overall X-axis coordinates of the two LOS images. This represents the overall Y-axis coordinate of the two LOS images. This represents the error plane between two geographically adjacent LOS images in the same orbital data in the coordinate system. The observed data below;

[0121] Since the azimuth and incident angles of satellite observation data are the same within the overlapping region of the same orbit, and without the constraint of GNSS 3D deformation observation data, the observation data of two geographically adjacent LOS images in the overlapping region should be identical. Therefore, the simplest method to eliminate the error plane distribution between two LOS images is to distribute the error plane equally. That is, for LOS image A11, the observation data after eliminating the error plane is... For LOS image A12, the observed data after eliminating the error plane is... LOS imagery is a raster image data composed of neatly arranged pixels. Each pixel has an observation value (the observation value of LOS imagery is the pixel value of the LOS imagery). By using the observation value data corresponding to each pixel on each LOS imagery, the LOS image after eliminating the error plane can be obtained.

[0122] Step S322: If the number of GNSS stations in the overlapping area of ​​two geographically adjacent LOS images (i.e., LOS image A11 and LOS image A12) in the same track data is greater than or equal to 10, construct the observation equation by combining the coordinates and observation data of the two LOS images and the GNSS three-dimensional deformation observation data, solve the error plane between the two LOS images and the GNSS three-dimensional deformation observation data plane respectively, and determine the model factor of the error plane between the two LOS images and the GNSS three-dimensional deformation observation data plane.

[0123] Taking the overlapping area of ​​two geographically adjacent LOS images in the same orbital data as an example, the observation equation constructed by combining the coordinates and observation data of the two LOS images with the GNSS three-dimensional deformation observation data is as follows:

[0124] ,

[0125] The above observation equation, after transformation, becomes:

[0126] ,

[0127] In the formula, and These represent the observation data of two geographically adjacent LOS images in the same orbital data, namely... The observation data represents LOS image A11. Observational data representing LOS image A12; This represents the east-west deformation data of GNSS in 3D deformation observation data. This represents the north-south oriented deformation data in GNSS three-dimensional deformation observation data. This represents the GNSS vertical deformation data in GNSS three-dimensional deformation observation data; and These correspond to the azimuth and incident angles of the observation data during satellite flight, respectively. This indicates the satellite observation azimuth angle corresponding to LOS image A11. This indicates the satellite observation azimuth angle corresponding to LOS image A12. This indicates the satellite observation incident angle corresponding to LOS image A11. This indicates the satellite observation incident angle corresponding to LOS image A12; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. The east-west deformation of the t-th pixel within the overlapping region of two LOS images is represented by this value. The north-south deformation of the t-th pixel within the overlapping region of two LOS images represents the deformation of the Lt-th pixel. The vertical deformation of the t-th pixel within the overlapping region of the two LOS images is represented; where, , , and , , These represent the model factors of the error plane between two geographically adjacent LOS imageries and the GNSS 3D deformation observation data plane in the same orbital data, respectively. , and The model factor representing the error plane between LOS image A11 and GNSS three-dimensional deformation observation data plane; , and The model factor represents the error plane between the LOS image A12 and the GNSS three-dimensional deformation observation data plane.

[0128] The model factor was obtained by solving the above observation equation. , and and model factors , and Then, using the GNSS 3D deformation observation data plane as the reference plane for error elimination, the error planes between the two LOS images and the GNSS 3D deformation observation data plane are eliminated respectively, thereby eliminating the error plane between the two LOS images; specifically:

[0129] Using formula To calculate the observed data for the error plane used to eliminate the gap between the LOS image A11 and the GNSS 3D deformation observation data plane; using the formula To calculate the observed data for the error plane used to eliminate the gap between the LOS image A12 and the GNSS 3D deformation observation data plane;

[0130] in, This represents the X-axis coordinate of LOS image A11. This represents the Y-axis coordinate of LOS image A11. This indicates the error plane between the LOS image A11 and the GNSS 3D deformation observation data plane in coordinate system. The observed data below; among them, This represents the X-axis coordinate of LOS image A12. This represents the Y-axis coordinate of LOS image A12. This indicates the error plane between the LOS image A12 and the GNSS 3D deformation observation data plane in coordinate system. The observed data below;

[0131] Therefore, for LOS image A11, the observed data after eliminating the error plane is... For LOS image A12, the observed data after eliminating the error plane is... .

[0132] Step S33: Based on multiple geographically adjacent LOS images in the same track data, the error plane between multiple pairs of adjacent LOS images is calculated in the same way. That is, the LOS images on the north and south sides of the track will only be corrected once, while those in the middle of the track will undergo two corrections from the LOS images on the north and south sides. This can eliminate or reduce the inconsistency between adjacent LOS images to the greatest extent. The LOS images on the same track after the error is eliminated are stitched together to obtain a smoother, continuous, and long strip-shaped LOS image along the track direction.

[0133] Step S4, secondary processing of GNSS three-dimensional deformation observation data:

[0134] Spatial interpolation is performed on the GNSS 3D deformation observation data, and then the spatially interpolated GNSS 3D deformation observation data is projected onto the orbital plane of each corresponding LOS image to obtain GNSS 3D deformation projection data; the specific steps are as follows:

[0135] Step S41: Based on the uncertainty of the observed values, the GNSS three-dimensional deformation observation data within the study area are screened and then spatially interpolated to obtain the spatially interpolated GNSS three-dimensional deformation observation data (including the spatially interpolated GNSS east-west deformation data). Spatially interpolated GNSS north-south deformation data Spatially interpolated GNSS vertical deformation data This facilitates subsequent use as a constraint condition for adjacent track splicing. Simultaneously, the spatially interpolated GNSS north-south deformation data... It can be directly used to participate in the inversion of subsequent three-dimensional surface deformation and the calculation of strain rate.

[0136] In this embodiment, spatial interpolation is performed on GNSS station data with observation uncertainty less than 0.7 mm / yr (i.e., 0.7 mm / year) for both east-west and north-south GNSS deformation data. For vertical GNSS deformation data, spatial interpolation is performed on all GNSS station data. Spatial interpolation methods include Kriging interpolation, inverse distance weighting, and Gaussian process regression interpolation. In this embodiment, Kriging interpolation is used. The results of spatial interpolation of the GNSS three-dimensional deformation observation data are as follows: Figure 6 As shown, Figure 6 (a) Figure 6 (b) Figure 6 (c) shows the spatially interpolated GNSS three-dimensional deformation observation data, in which... Figure 6 In the text, (a) represents the spatially interpolated GNSS east-west deformation data; Figure 6 (b) in the figure represents the spatially interpolated GNSS north-south deformation data; Figure 6 (c) in the text represents the spatially interpolated GNSS vertical deformation data; Figure 6 (d) in Figure 6 (e) in Figure 6 (f) in the figure shows the uncertainty of the three-dimensional deformation result after spatial interpolation, where, Figure 6 In this context, (d) represents the uncertainty of the GNSS east-west deformation result after spatial interpolation; Figure 6 In this context, (d) represents the uncertainty of the GNSS north-south deformation result after spatial interpolation; Figure 6 In this context, (d) represents the uncertainty of the vertical deformation result of the GNSS after spatial interpolation;

[0137] Step S42: Then, the spatially interpolated GNSS 3D deformation observation data is projected onto the orbital plane (LOS direction plane) of each corresponding LOS image, so that the spatial resolution of the spatially interpolated GNSS 3D deformation observation data is the same as that of the LOS image, thus obtaining GNSS 3D deformation projection data. The projection formula is expressed as follows:

[0138] ,

[0139] In the formula, This represents the data after spatially interpolated GNSS 3D deformation observation data is projected onto the orbital plane of the LOS image, i.e., GNSS 3D deformation projection data. and These correspond to the azimuth and angle of incidence of the satellite's observation data during flight, respectively. This represents the spatially interpolated GNSS three-dimensional deformation observation data, where, This represents the spatially interpolated GNSS east-west deformation data. This represents the spatially interpolated GNSS north-south deformation data. This represents the spatially interpolated GNSS vertical deformation data.

[0140] Step S5, adjacent track data stitching:

[0141] The process involves acquiring stitched LOS (Longest-of-Sight) images along the track direction for two geographically adjacent orbits to form adjacent-track data. Combined with GNSS 3D deformation projection data, the error plane between the two stitched LOS images of adjacent orbits is solved. Weighted balancing of the error plane is then performed on each stitched LOS image to obtain a corrected stitched LOS image, which is then stitched together with the adjacent orbits to obtain a stitched LOS image covering the study area. The specific steps are as follows:

[0142] Step S51: Obtain the strip-shaped LOS images (i.e., strip-shaped LOS image A1 and strip-shaped LOS image A2) that have been stitched along the track direction for two geographically adjacent tracks, form adjacent track data, and combine them with GNSS three-dimensional deformation projection data to solve the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS three-dimensional deformation projection data plane in the adjacent track data, and determine the model factor of the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS three-dimensional deformation projection data plane;

[0143] The error plane between two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformable projection data plane is determined by solving the adjacent track data. The calculation formula for the model factor of the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformable projection data plane is as follows:

[0144] ,

[0145] The above calculation formula, after being transformed, becomes:

[0146] ,

[0147] In the formula, The error plane representing the relationship between two strip-shaped LOS images (strip-shaped LOS image A1 and strip-shaped LOS image A2) on adjacent orbits. This represents the error plane between the elongated LOS image A1 and the GNSS 3D deformation projection data plane. This represents the error plane between the elongated LOS image A2 and the GNSS 3D deformation projection data plane. and The observed data represent two strip-shaped LOS images of adjacent orbits, respectively. The observation data representing the elongated LOS image A1. Observational data representing the elongated LOS image A2; and These represent the GNSS 3D deformation projection data of two geographically adjacent orbits in the adjacent orbit data, namely... GNSS 3D deformation projection data corresponding to the elongated LOS image A1 GNSS 3D deformation projection data corresponding to the elongated LOS image A2; and The x-axis and y-axis coordinates of the t-th pixel within the overlapping region of two elongated LOS images representing adjacent orbits; where... , , and , , The model factors represent the error planes between two strip-shaped LOS images of adjacent orbits and the GNSS 3D deformation projection data plane, respectively. , and The model factor representing the error plane between the strip-shaped LOS image A1 and the GNSS three-dimensional deformation projection data plane; , and The model factor represents the error plane between the strip-shaped LOS image A2 and the GNSS three-dimensional deformation projection data plane.

[0148] The model factors were obtained by solving the above calculation formula. , and and model factors , and Then, using the GNSS 3D deformation projection data plane as the reference plane for error elimination, the error planes between the two strip-shaped LOS images and the GNSS 3D deformation projection data plane are eliminated respectively, thereby eliminating the error planes between the two strip-shaped LOS images; specifically:

[0149] Using formula To calculate the observed data used to eliminate the error plane between the strip-shaped LOS image A1 and the GNSS 3D deformation projection data plane, the formula is used. To calculate the observation data for the error plane used to eliminate the gap between the strip-shaped LOS image A2 and the GNSS three-dimensional deformation projection data plane;

[0150] in, This represents the X-axis coordinate of the elongated LOS image A1. This represents the Y-axis coordinate of the elongated LOS image A1. The error plane representing the coordinate system between the elongated LOS image A1 and the GNSS 3D deformation projection data plane. The observed data below; This represents the X-axis coordinate of the elongated LOS image A2. This represents the Y-axis coordinate of the elongated LOS image A2. The error plane representing the coordinate system between the elongated LOS image A2 and the GNSS 3D deformation projection data plane. The observed data below.

[0151] Therefore, for the elongated LOS image A1, after eliminating the error plane between the elongated LOS image A1 and the GNSS 3D deformation projection data plane, the obtained observation data is: For the elongated LOS image A2, after eliminating the error plane between the elongated LOS image A2 and the GNSS 3D deformation projection data plane, the obtained observation data is: ; and These represent the observation data of the two strip-shaped LOS images obtained after eliminating the error plane between the strip-shaped LOS image and the GNSS three-dimensional deformation projection data plane.

[0152] Step S52: Since the azimuth and incident angles of the satellite flight observation data corresponding to different orbits are different, there is a fixed step plane between two strip-shaped LOS images of adjacent orbits. Based on the overlapping area data between multiple geographically adjacent strip-shaped LOS images that have been stitched along the orbit direction in the adjacent orbit data, the error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is solved, and the western correction value of each strip-shaped LOS image is calculated. and East Side Correction Value Then, based on the area difference of the overlapping regions on the east and west sides of each strip-shaped LOS image track, weighted balancing of the error plane is performed on each strip-shaped LOS image. This allows for a secondary overall balancing of each strip-shaped LOS image using a fixed step plane between two adjacent tracks. Specifically, the following steps are included:

[0153] Step S521: Using the overlapping region data between two strip-shaped LOS images of adjacent tracks, solve for the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks, and determine the model factor of the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks. The calculation formula is as follows:

[0154] ,

[0155] The above calculation formula, after being transformed, becomes:

[0156] ,

[0157] In the formula, and The observed data represent two strip-shaped LOS images of adjacent orbits, respectively. The observation data representing the elongated LOS image A1. Observational data representing the elongated LOS image A2; and These represent the GNSS 3D deformation projection data of two geographically adjacent orbits in the adjacent orbit data, namely... GNSS 3D deformation projection data corresponding to the elongated LOS image A1 GNSS 3D deformation projection data corresponding to the elongated LOS image A2; and The X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of two long strip-shaped LOS images of adjacent tracks are respectively represented. , and The model factor representing the error plane between two adjacent LOS images caused by orbital errors;

[0158] The model factor for the error plane caused by orbital error between two long strip-shaped LOS images of adjacent orbits is obtained by solving the above calculation formula. , and Next, the two elongated LOS images are defined as a single entity, and then the formula is used. To calculate the observed data for the error plane between two long strip LOS images of adjacent orbits, which is used to eliminate the error caused by orbital errors, where, The X-axis coordinates of the two elongated LOS images of adjacent orbits are represented. The Y-axis coordinate represents the overall coordinates of two long strip-shaped LOS images of adjacent orbits. This represents the error plane between two strip-shaped LOS images (strip-shaped LOS image A1 and strip-shaped LOS image A2) on adjacent orbits in coordinates. The observation data below is used to eliminate the error plane caused by orbital errors between two strip-shaped LOS images (strip-shaped LOS image A1 and strip-shaped LOS image A2) of adjacent orbits;

[0159] Based on the above principle, the error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is solved respectively, and the model factor of the error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is determined; and the observed data used to eliminate the error plane caused by orbital error between the strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit are calculated respectively.

[0160] The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its west is defined as the westward correction value of the strip-shaped LOS image, denoted as . The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its east of its orbit is defined as the eastward correction value of the strip-shaped LOS image, denoted as . ;

[0161] Taking the strip-shaped LOS image A2 as an example, its two adjacent strip-shaped LOS images on the east and west sides are strip-shaped LOS image A3 and strip-shaped LOS image A1, respectively. Then, the observed data used to eliminate the error plane caused by orbital errors between strip-shaped LOS image A2 and strip-shaped LOS image A1... The western correction value is the value used to correct the error between the long strip LOS image A2 and the long strip LOS image A3. The eastern correction value is the value used to correct the error between the long strip LOS image A2 and the long strip LOS image A3 caused by orbital errors.

[0162] Step S522: Then, based on the size of the area difference of the overlapping region on the east and west sides of each strip-shaped LOS image track (i.e., the size of the area difference of the overlapping region between a single strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its track), weighted balancing of the error plane is performed on each strip-shaped LOS image. The specific steps are as follows:

[0163] The area of ​​the overlapping region on the eastern side of a strip-shaped LOS image (i.e., the area of ​​the overlapping region between a single strip-shaped LOS image and a strip-shaped LOS image on the eastern side of its orbit) is defined as S. 东 The area of ​​the overlapping region on the west side of the strip-shaped LOS image (i.e., the area of ​​the overlapping region between a single strip-shaped LOS image and the strip-shaped LOS image on the west side of its orbit) is S. 西 ;

[0164] If the area of ​​the overlapping region on the east and west sides of a long strip-shaped LOS image differs significantly (S) 东 / S 西 ≤0.8 or S 东 / S 西 If the error plane is ≥1.2), then no error plane correction is performed on the east and west sides to avoid insufficient correction accuracy from increasing the error.

[0165] If the area difference between the overlapping regions on the east and west sides of the elongated LOS image is not significant (0.8), 东 / S 西 <1.2, S 东 S represents the area of ​​the overlapping region on the eastern side of the elongated LOS image. 西 (where the area of ​​the overlapping region on the west side of the strip-shaped LOS image is defined as the area), and the correction value for weighted balancing of the error plane for each strip-shaped LOS image is defined as follows: The corrected observation data for the elongated LOS image are: ,in, This represents the observation data of the long strip-shaped LOS image obtained in step S51 after eliminating the error plane caused by the overlapping region;

[0166] ​Step S53: After weighted balancing of each strip-shaped LOS image with an error plane, a corrected strip-shaped LOS image is obtained. Then, the corrected strip-shaped LOS images are stitched together to obtain a stitched LOS image covering the study area, including a stitched ascending orbit LOS image covering the study area and a stitched descending orbit LOS image covering the study area. In this embodiment, the stitched ascending orbit LOS image covering the study area is as follows: Figure 7 As shown, the stitched LOS image covering the study area is as follows: Figure 8 As shown.

[0167] Step S6, 3D deformation inversion of the Earth's surface:

[0168] By combining the stitched LOS images of both ascending and descending orbits covering the study area, and the spatially interpolated GNSS north-south deformation data obtained in step S4, a three-dimensional surface deformation inversion equation is constructed. The three-dimensional surface deformation inversion is then performed, and the three-dimensional surface deformation rate field is solved based on the equation to obtain the three-dimensional surface deformation data. The three-dimensional surface deformation inversion equation is as follows:

[0169] ,

[0170] The above three-dimensional deformation inversion equation for the Earth's surface, after deformation, becomes:

[0171] ,

[0172] In the formula, The up-track LOS image representing the stitched-together coverage area of ​​the study region. The down-orbiting LOS image representing the stitched-together coverage area of ​​the study region. Represents spatially interpolated GNSS north-south deformation data; and These correspond to the azimuth and incident angles of the observation data during satellite flight, respectively. The satellite azimuth angle representing the orbital ascent type data. The satellite azimuth angle representing the down-orbit type data. Indicates the satellite incident angle for ascending orbit type data. Indicates the satellite incident angle for data of the down-orbit type; , and This represents the three-dimensional deformation rate field of the Earth's surface that needs to be inverted, where... Represents the vertical deformation rate field. Represents the north-south deformation rate field. The east-west deformation rate field represents the value that needs to be solved in the three-dimensional deformation inversion equation of the Earth's surface. Among these, the east-west deformation rate field... and vertical deformation rate field It is based on the stitched LOS imagery covering the research area. and LOS imagery The obtained north-south deformation rate field It is based on the spatially interpolated GNSS north-south deformation data The three sets of deformation rate field data obtained from the solution together constitute the three-dimensional deformation data of the Earth's surface.

[0173] Figure 9 This embodiment displays east-west and vertical deformation rate field maps retrieved from LOS image data; these deformation rate field maps visually demonstrate the east-west and vertical deformation rate field maps of the study area (northern edge of the Tibetan Plateau). ) and vertical ( The three-dimensional crustal movement characteristics on the surface.

[0174] Figure 9 The east-west deformation rate field shown in (a) exhibits a highly regular north-south zoning characteristic: the southern study area (the interior of the Tibetan Plateau) shows a significant eastward movement, corresponding to... Figure 9 The red area in (a) shows a rate of approximately 8–10 mm / yr (8–10 mm / year), while the northern study area (the stable Alashan block) maintains a lower eastward rate, corresponding to… Figure 9 The blue region in (a) shows a velocity of approximately 2 mm / yr. This strong velocity gradient is concentrated at the main boundaries of the Haiyuan Fault and the Qilian Mountain Fault Zone, quantitatively reflecting the dynamic process of the eastward "escape" of materials from the Tibetan Plateau under the compression of the Indian Plate.

[0175] Figure 9 The vertical deformation rate field shown in (b) reveals a more complex signal of local tectonic evolution and surface adjustment, with rate values ​​ranging from -5 to 1 mm / yr. Experiments have observed a clear uplift trend in the uplifted areas of mountains such as the Qilian Mountains, corresponding to… Figure 9 The orange-red area in (b) corresponds to the tectonic setting of regional crustal shortening and outward thrust thickening; while at some basin edges or fault intersections, there are obvious subsidence centers, corresponding to Figure 9 The dark blue area in (b) reflects the complex strain distribution in this region, or the influence of local non-tectonic factors (such as groundwater utilization and permafrost changes). By integrating the horizontal and vertical deformation field results, this data not only precisely defines the slip distribution of the plateau boundary faults, but also provides crucial boundary constraints for establishing a dynamic model of the "compression-uplift-escape" dynamics of the northeastern margin of the Tibetan Plateau through the coupling characteristics of the vertical and horizontal directions.

[0176] Step S7, Calculation of strain rate data:

[0177] Based on three-dimensional surface deformation data, strain rate data is calculated. This strain rate data includes second-order invariants of horizontal expansion rate, shear rate, and horizontal strain rate tensor, yielding the evaluation results of the three-dimensional surface deformation inversion. The specific steps include:

[0178] Step S71: Calculate the horizontal strain rate tensor based on the three-dimensional surface deformation data obtained in step S6. The calculation formula is as follows:

[0179] ,

[0180] In the formula, Represents the horizontal strain rate tensor. , , , This represents the velocity gradient in the horizontal direction, where... This represents the rate of change of the east-west velocity with respect to the X-axis coordinate. This represents the rate of change of the north-south velocity with respect to the Y-axis coordinate. = , , Both are quantities with the same meaning, measuring the average effect of the change in east-west velocity in the north-south direction and the change in north-south velocity in the east-west direction, reflecting the shear deformation of the object; x represents the distance between two adjacent pixels in the east-west direction in the horizontal direction, and y represents the distance between two adjacent pixels in the north-south direction in the horizontal direction. Taking the partial derivative of the east-west deformation rate with respect to the X-axis is a mathematical operation. This represents the partial derivative of the east-west deformation rate with respect to the Y-axis. This represents the partial derivative of the north-south deformation rate with respect to the X-axis. The partial derivative of the north-south deformation rate with respect to the Y-axis is given.

[0181] Step S72: Based on the velocity gradients in each horizontal direction of the horizontal strain rate tensor obtained in step S71, calculate the horizontal expansion rate, shear rate, and the second-order invariants of the horizontal strain rate tensor. The calculation formulas are as follows:

[0182] ,

[0183] In the formula, Represents the horizontal expansion rate. Represents the shear rate. The second-order invariant representing the horizontal strain rate tensor;

[0184] Second-order invariant data of horizontal expansion rate, shear rate, and horizontal strain rate tensors can be used for crustal deformation analysis and geological engineering assessment. A larger positive value for the horizontal expansion rate indicates stronger stretching and a significant surface expansion trend, corresponding to normal fault zones and crustal extension zones. A larger negative value for the horizontal expansion rate indicates stronger compression and a significant surface shortening trend, corresponding to thrust fault zones and plate compression zones. A larger shear rate value indicates stronger shearing and more intense relative slippage of landmasses, corresponding to strike-slip fault zones and landmass boundary slip zones. A larger second-order invariant value for the horizontal strain rate tensor indicates greater overall deformation intensity and more intense crustal activity, corresponding to fault activity zones and high-risk earthquake zones.

[0185] The horizontal expansion rate distribution diagram calculated in this embodiment is as follows: Figure 10 As shown in the figure, the shear rate distribution diagram is as follows: Figure 11 As shown, the second-order invariant plot of the horizontal strain rate tensor is as follows: Figure 12 As shown.

[0186] Figure 10 The horizontal expansion rate distribution map shown clearly demonstrates the strong spatial heterogeneity and tectonic control characteristics of crustal deformation in the study area (northeastern margin of the Tibetan Plateau). Figure 10 The numerical distribution of horizontal expansion rate in the figure shows that the horizontal expansion rate in this region exhibits a clear pattern of "banded accumulation and blocky distribution." Specifically, the high strain rate region (i.e., the dark blue area) has a horizontal expansion rate of approximately 60 nst / yr (i.e., 60 nanostrain / year), highly concentrated at the edges of the main active fault zones. In particular, the high-value bands trending near-east-west and northwest across the middle of the figure accurately delineate the activity traces of the Haiyuan Fault, Qilian Mountain Fault, and East Kunlun Fault, reflecting that these fault zones, acting as strong deformation boundaries, absorbed most of the strain energy from the northeastward compression of the plateau. Meanwhile, Figure 10 The large red areas (regions with horizontal expansion rates close to 0) correspond to relatively rigid tectonic units within the study area (the edges of the Alashan Block and the Ordos Block, as well as the interior of the plateau), indicating lower deformation rates and stronger overall stability within the blocks. This characteristic of "high-strain zones surrounding low-strain blocks" strongly supports the crustal movement model in this region, suggesting that deformation mainly occurs at block boundaries rather than within the blocks. Further analysis reveals that the convergence and turning points of multiple high-strain zones near 104°E to 106°E indicate strain intensification caused by the obstruction of surrounding stable blocks during the outward expansion of plateau material. This provides intuitive and precise physical evidence for identifying regional seismic hazard zones and understanding the dynamic evolution of the northeastern margin of the Tibetan Plateau.

[0187] Figure 11The shear rate distribution map shown further reveals the dynamic characteristics of crustal deformation in the study area (northern margin of the Tibetan Plateau), which exhibits strong spatial localization and strike-slip tectonic control. Experimental results show that the spatial distribution of shear rates is highly consistent with the main fault zones in the region, with values ​​fluctuating dramatically between -50 and 50 nst / yr (nanostrain / year). Figure 11 The most striking feature is the continuous high-value band of deep red across the central part (approximately 37°N), which accurately depicts the geometry of the Haiyuan Fault Zone. Its significant negative values ​​directly reflect the strong left-lateral strike-slip nature and high strain accumulation rate of this fault zone. Meanwhile, the deep red band in the south (approximately 99°–102°E, 33°–35°N) corresponds to the eastern extension of the East Kunlun Fault Zone, while the alternating blue and red areas distributed in the Qilian Mountains (approximately 102°–106°E, 33°–35°N) reveal the complex shear strain field and its moderating role in absorbing block rotation and compression. In contrast, Figure 11 The white and light-colored regions (regions with shear rates close to 0) are mainly distributed in the southern part of the Alashan Block and the western edge of the Ordos Block, indicating extremely low shear deformation within these rigid blocks. This pattern of "strong shear zones surrounding weak shear blocks" strongly supports the fault-controlled model of crustal movement on the northeastern edge of the Tibetan Plateau, suggesting that tectonic displacement within the region is mainly released through strike-slip movement of boundary faults. These data not only quantify the activity intensity of major faults but also provide core evidence for determining the shear strain background of high-risk earthquake zones.

[0188] Figure 12 The second-order invariant plot of the horizontal strain rate tensor clearly demonstrates the spatial localization and tectonic control characteristics of the total crustal deformation intensity in the study area (northern margin of the Tibetan Plateau). As a comprehensive indicator of the intensity of deformation, the second-order invariant of the horizontal strain rate tensor exhibits a clear pattern of "linear concentration and block stability." In the high-value region (i.e., the dark blue area), the peak value of the second-order invariant of the horizontal strain rate tensor approaches 80 nst / yr (nanostrain / year), highly defining the geometric outline of the main active fault zones in the region. Among these, the transverse... Figure 12 The Haiyuan fault zone in the middle, the Qilian Mountain fault system in the north, and the East Kunlun fault zone in the south are the most prominent. These high-value bands visually depict the tectonic boundaries where strain energy accumulation is most concentrated, reflecting that crustal deformation under plate compression is mainly absorbed and regulated by these large fault zones. Meanwhile, Figure 12The large red background (regions where the second-order invariant of the horizontal strain rate tensor is close to 0) corresponds to relatively stable tectonic units within the Alashan Block, the Ordos Block, and the Tibetan Plateau, indicating extremely low deformation rates within these rigid blocks. This characteristic of "high-strain intensity zones surrounding low-strain rigid blocks" strongly supports the conclusion that deformation in this region follows a micro-block movement model, meaning that the total deformation intensity is not diffusely distributed within the region but exhibits a strong boundary effect. Furthermore, the strain intensification and arc-shaped transition of the second-order invariant at the eastern boundary (approximately 104°–106° E) due to the obstruction of stable blocks quantitatively reveals the stress field evolution characteristics during the northeastward "escape" of plateau material, providing core physical evidence for identifying long-term seismic hazard zones and understanding the dynamic mechanism of the outward expansion of the Tibetan Plateau.

[0189] The method of this invention can significantly reduce processing costs under resource-constrained conditions; effectively reduce step and inconsistency in overlapping areas, and improve the intrinsic consistency of the deformation rate field of the stitched LOS image; achieve a balance between geometric constraints and external references, thereby improving the accuracy of subsequent three-dimensional deformation inversion; and support the fusion of multi-track and multi-source data to adapt to the inversion needs of large areas.

[0190] This invention employs a combination of along-track and cross-track stitching, significantly improving stitching accuracy compared to traditional along-track stitching methods. In this embodiment, the LOS image stitched using this method is compared to that stitched using traditional along-track stitching. The root mean square error (RMSE) is calculated based on the projection of GNSS 3D deformation observation data onto the LOS direction. The comparison of the RMSE results is shown in the figure below. Figure 13As shown, the horizontal axis represents the InSAR comparison clusters for different orbits, with blue representing the results of along-orbit stitching only, and orange representing the results of both along-orbit and cross-orbit stitching. The vertical axis represents the comparison with GNSS data RMSE values. Compared with the traditional along-orbit stitching method, the along-orbit stitching plus cross-orbit stitching method of this invention significantly improves the accuracy of surface deformation monitoring; specifically, the RMSE value corresponding to orbit 055A decreased from 3.33 mm / yr (mm / year) to 3.20 mm / yr (mm / year), a reduction of 3.93%; the RMSE value corresponding to orbit 164D improved from 3.29 mm / yr to 3.16 mm / yr, a reduction of 3.80%. Notably, orbit 135D showed the most significant improvement, with its RMSE value dropping sharply from 2.45 mm / yr to 1.70 mm / yr, a reduction of 30.31%. The RMSE value for track 106D decreased from 2.22 mm / yr to 2.03 mm / yr, a reduction of 8.23%; the RMSE value for track 062D decreased from 1.80 mm / yr to 1.68 mm / yr, a reduction of 6.49%; and the RMSE value for track 128A decreased from 1.66 mm / yr to 1.59 mm / yr, a reduction of 3.97%. The decrease in RMSE value indicates improved stitching accuracy. The magnitude of error reduction is highly correlated with the initial error value; the abnormal performance of track 135D suggests that its effectiveness in correcting high-error datasets is not very robust.

[0191] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for surface deformation inversion based on InSAR deformation fields stitched together in overlapping regions, characterized in that, Includes the following steps: Step S1: Acquire multiple frames of LOS images within the study area, the corresponding image frame information, and GNSS observation data. The GNSS observation data includes observation data, observation uncertainty, GNSS three-dimensional deformation observation data, and GNSS station data. The image frame information includes the coordinates, orbit type, and orbit number of the LOS images. Step S2: Classify LOS images according to orbit type and orbit number. Group LOS images with the same orbit type and orbit number into a set of co-orbit data to obtain a co-orbit data set. Then, based on the co-orbit data set, filter out the overlapping areas of two geographically adjacent LOS images in the co-orbit data to obtain the overlapping areas of the same orbit. Finally, obtain GNSS station data and GNSS three-dimensional deformation observation data in the overlapping areas of the same orbit. Step S3: Based on the GNSS station data and GNSS 3D deformation observation data in the overlapping area of ​​the same track, construct the observation equation, use the observation equation to solve the error plane between two geographically adjacent LOS images in the same track data, eliminate the error between adjacent LOS images and perform same track stitching to obtain a long strip LOS image stitched along the track direction. Step S4 involves spatially interpolating the GNSS 3D deformation observation data within the study area, and then projecting the spatially interpolated GNSS 3D deformation observation data onto the orbital plane of each corresponding LOS image to obtain GNSS 3D deformation projection data; specifically, this includes the following steps: Step S41: Based on the uncertainty of the observed values, the GNSS three-dimensional deformation observation data within the study area are screened and then spatially interpolated to obtain the spatially interpolated GNSS three-dimensional deformation observation data, including the spatially interpolated GNSS east-west deformation data. Spatially interpolated GNSS north-south deformation data GNSS vertical deformation data after spatial interpolation ; Step S42: Project the spatially interpolated GNSS 3D deformation observation data onto the orbital plane of each corresponding LOS image to obtain GNSS 3D deformation projection data. The projection formula is expressed as follows: ， In the formula, The data represents the spatially interpolated GNSS 3D deformation observation data projected onto the orbital plane of the LOS image, i.e., GNSS 3D deformation projection data; and These correspond to the azimuth and incident angles of the observation data during satellite flight, respectively. This represents the spatially interpolated GNSS east-west deformation data. This represents the spatially interpolated GNSS north-south deformation data. Represents the spatially interpolated GNSS vertical deformation data; Step S5: Obtain the strip-shaped LOS images that have been stitched along the track direction for two geographically adjacent tracks to form adjacent track data. Combine the data with GNSS 3D deformation projection data to solve the error plane between the two strip-shaped LOS images of adjacent tracks. Perform weighted balancing of the error plane on each strip-shaped LOS image to obtain the corrected strip-shaped LOS image. Then, stitch the images together with the adjacent tracks to obtain the stitched LOS image covering the study area. Step S6: Combine the stitched LOS imagery covering the study area with the spatially interpolated GNSS 3D deformation observation data to construct the surface 3D deformation inversion equation, perform surface 3D deformation inversion, and obtain surface 3D deformation data. Step S7: Calculate strain rate data based on the three-dimensional deformation data of the earth's surface to obtain the evaluation results of the three-dimensional deformation inversion of the earth's surface.

2. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 1, characterized in that, In step S1, the GNSS three-dimensional deformation observation data includes GNSS east-west deformation data, GNSS north-south deformation data, and GNSS vertical deformation data; GNSS site data includes the number of GNSS sites and their location information; Track types include ascending track and descending track.

3. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 2, characterized in that, Step S3 specifically includes the following steps: Step S31, define the threshold for the number of GNSS stations as M. GNSS ; Step S32: Based on the GNSS station data in the overlapping area of ​​two geographically adjacent LOS images in the same orbit, determine whether the number of GNSS stations in the overlapping area is greater than M. GNSS ; If the number of GNSS stations is less than M within the overlapping area of ​​two geographically adjacent LOS images in the same track data, then... GNSS Then, the observation equation is directly constructed based on the coordinates and observation data of the two LOS images; If the number of GNSS stations is greater than or equal to M within the overlapping area of ​​two geographically adjacent LOS images in the same track data, then... GNSS Then, the coordinates and observation data of the two LOS images and the GNSS three-dimensional deformation observation data are combined to construct the observation equation; The error plane between two geographically adjacent LOS images in the same track data is solved by using the observation equation, thus eliminating the error between adjacent LOS images. Step S33: Stitch together the LOS images of the same track after error elimination to obtain a long strip LOS image stitched along the track direction.

4. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 3, characterized in that, In step S32, the observation equation is directly constructed based on the coordinates and observation data of the two LOS images as follows: ， In the formula, This represents the error plane between two geographically adjacent LOS images in the same orbital data. and These represent the observation data of two geographically adjacent LOS images in the same orbital data; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. , and It is the model factor of the error plane between two geographically adjacent LOS images in the same track data; The specific steps for eliminating errors between adjacent LOS images are as follows: The error plane between two geographically adjacent LOS images in the same orbital data is solved using the observation equation. The model factors were obtained by solving the observation equation. , and Next, the two LOS images are defined as a single entity, and then the formula is used. To calculate the observation data used to eliminate the error plane between two LOS images, where, Represents the overall X-axis coordinates of the two LOS images. This represents the overall Y-axis coordinate of the two LOS images. This represents the error plane between two geographically adjacent LOS images in the same orbital data, expressed in coordinates. The observed data below; Then, the error plane between two geographically adjacent LOS images in the same orbital data is evenly distributed to eliminate the error between the two LOS images; for the two LOS images, the observed data after eliminating the error plane are respectively... and .

5. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 3, characterized in that, In step S32, the following observation equation is constructed by combining the coordinates and observation data of the two LOS images with the GNSS three-dimensional deformation observation data: ， In the formula, and These represent the observation data of two geographically adjacent LOS images in the same orbital data; This represents the east-west deformation data of GNSS in 3D deformation observation data. This represents the north-south deformation data in GNSS three-dimensional deformation observation data. This represents the GNSS vertical deformation data in GNSS three-dimensional deformation observation data; and These represent the satellite observation azimuth angles corresponding to two geographically adjacent LOS images in the same orbital data; and These represent the satellite observation incident angles corresponding to two geographically adjacent LOS images in the same orbital data; and These represent the X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of the two LOS images, respectively. The east-west deformation of the t-th pixel within the overlapping region of two LOS images is represented by this value. The north-south deformation of the t-th pixel within the overlapping region of two LOS images represents the deformation of the Lt-th pixel. The vertical deformation of the t-th pixel within the overlapping region of two LOS images; , , and , , These represent the model factors of the error plane between two geographically adjacent LOS images and the GNSS 3D deformation observation data plane in the same track data; The specific steps for eliminating errors between adjacent LOS images are as follows: The error plane between two geographically adjacent LOS images in the same orbital data is solved using the observation equation. The model factors were obtained by solving the observation equation. , , and , , Then, using the GNSS 3D deformation observation data plane as the reference plane for error elimination, the error planes between the two LOS images and the GNSS 3D deformation observation data plane are eliminated respectively, thereby eliminating the error plane between the two LOS images; specifically: Using formula and formula Calculate the observation data for the error plane used to eliminate the gap between the two LOS imagery planes and the GNSS 3D deformation observation data plane; where, and These represent the X-axis and Y-axis coordinates of one of the LOS images, respectively. This indicates the error plane between the corresponding LOS imagery and the GNSS 3D deformation observation data plane in coordinates. The observed data below; and These represent the X-axis and Y-axis coordinates of another LOS image, respectively. This indicates the error plane between the corresponding LOS imagery and the GNSS 3D deformation observation data plane in coordinates. The observed data below; For the two LOS images, the observed data after eliminating the error plane are respectively and .

6. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 5, characterized in that, Step S5 specifically includes the following steps: Step S51: Obtain the strip-shaped LOS images stitched along the track direction corresponding to two geographically adjacent tracks to form adjacent track data. Combine this with GNSS 3D deformation projection data to solve the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformation projection data plane. Determine the model factor of the error plane between the two strip-shaped LOS images of adjacent tracks and the GNSS 3D deformation projection data plane. The calculation formula is as follows: ， In the formula, This represents the error plane between two long strip-shaped LOS images of adjacent orbits. This represents the error plane between one of the elongated LOS imagery and the GNSS 3D deformation projection data plane. This represents the error plane between another elongated LOS image and the GNSS 3D deformation projection data plane; and The observation data represent two strip-shaped LOS images of adjacent orbits, respectively; and These represent the GNSS 3D deformation projection data of two geographically adjacent tracks in the adjacent track data, respectively, and correspond to the GNSS 3D deformation projection data of two long strip-shaped LOS images of the adjacent tracks; and The X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of two long strip-shaped LOS images representing adjacent tracks; , , and , , The model factors represent the error planes between two strip-shaped LOS images of adjacent orbits and the GNSS three-dimensional deformation projection data plane, respectively. Solving for the model factors , , and , , Then, using the GNSS 3D deformation projection data plane as the reference plane for error elimination, the error planes between the two strip-shaped LOS images and the GNSS 3D deformation projection data plane are eliminated respectively, thereby eliminating the error planes between the two strip-shaped LOS images; specifically: Using formula and formula Calculate the observation data for the error plane used to eliminate the gap between the two strip-shaped LOS images and the GNSS 3D deformation projection data plane; where, and These represent the X-axis and Y-axis coordinates of one of the elongated LOS images, respectively. This indicates the error plane between the corresponding elongated LOS image and the GNSS 3D deformation projection data plane in coordinates. The observed data below; and These represent the X-axis and Y-axis coordinates of another elongated LOS image, respectively. This indicates the error plane between the corresponding elongated LOS image and the GNSS 3D deformation projection data plane in coordinates. The observed data below; For two strip-shaped LOS images, after eliminating the error plane between the strip-shaped LOS images and the GNSS 3D deformation projection data plane, the resulting observation data are as follows: and ; and These represent the observation data of the two strip-shaped LOS images obtained after eliminating the error plane between the strip-shaped LOS image and the GNSS three-dimensional deformation projection data plane; Step S52: Using the overlapping region data between two strip-shaped LOS images of adjacent tracks, solve for the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks, and determine the model factor of the error plane caused by track errors between the two strip-shaped LOS images of adjacent tracks. The calculation formula is as follows: ， In the formula, and The observation data represent two strip-shaped LOS images of adjacent orbits, respectively; and These represent the GNSS 3D deformation projection data on two geographically adjacent tracks in the adjacent track data, corresponding to the GNSS 3D deformation projection data of two long strip-shaped LOS images; and The X-axis and Y-axis coordinates of the t-th pixel within the overlapping region of two long strip-shaped LOS images of adjacent tracks are respectively represented. , and The model factor representing the error plane between two adjacent LOS images caused by orbital errors; The model factor of the error plane caused by orbital errors between two long strip-shaped LOS images of adjacent orbits is obtained by solving the problem. , and Next, the two elongated LOS images are defined as a single entity, and then the formula is used. To calculate the observed data for the error plane between two long strip LOS images of adjacent orbits, which is used to eliminate the error caused by orbital errors, where, The X-axis coordinates of the two elongated LOS images of adjacent orbits are represented. The Y-axis coordinate of the two elongated LOS images of adjacent orbits is represented. The error plane representing the coordinate system between two long strip-shaped LOS images of adjacent orbits. The observation data below is used to eliminate the error plane caused by orbital errors between two long strip LOS images of adjacent orbits; The error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is solved separately. The model factor of the error plane caused by orbital error between each strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit is determined. The observed data used to eliminate the error plane caused by orbital error between the strip-shaped LOS image and the two strip-shaped LOS images on the east and west sides of its orbit are calculated separately. The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its west is defined as the westward correction value of the strip-shaped LOS image, denoted as . The observed data used to eliminate the error plane caused by orbital errors between the strip-shaped LOS image and the strip-shaped LOS image to its east of its orbit is defined as the eastward correction value of the strip-shaped LOS image, denoted as . ; Then, based on the size of the area difference between the overlapping areas on the east and west sides of each strip-shaped LOS image track, weighted balancing of the error plane is performed on each strip-shaped LOS image. The specific steps are as follows: Define the area of ​​the overlapping region on the eastern side of the elongated LOS image as S. 东 The area of ​​the overlapping region on the west side of the elongated LOS image is S. 西 ; If the area difference between the overlapping regions on the east and west sides of the elongated LOS image is S 东 / S 西 ≤0.8 or S 东 / S 西 If the error plane correction is ≥1.2, then no error plane correction is performed on the east and west sides; If the area difference between the overlapping regions on the east and west sides of the elongated LOS image is 0.8... 东 / S 西 <1.2, the correction value for weighted balancing of the error plane for each strip-shaped LOS image is defined as: The corrected observation data for the elongated LOS image are: ,in, This represents the observation data of the long strip-shaped LOS image after eliminating the error plane caused by overlapping regions;​ Step S53: After weighted balancing of each strip-shaped LOS image with an error plane, a corrected strip-shaped LOS image is obtained; then, the corrected strip-shaped LOS images are stitched together according to the orbit type to obtain a stitched LOS image covering the study area, including a stitched ascending orbit LOS image covering the study area and a stitched descending orbit LOS image covering the study area.

7. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 6, characterized in that, In step S6, the three-dimensional deformation inversion equation of the Earth's surface is as follows: ， In the formula, The up-track LOS image representing the stitched-together coverage area of ​​the study region. The down-orbiting LOS image representing the stitched-together coverage area of ​​the study region. Represents spatially interpolated GNSS north-south deformation data; The satellite azimuth angle representing the orbital ascent type data. The satellite azimuth angle representing the down-orbit type data. Indicates the satellite incident angle for ascending orbit type data. Indicates the satellite incident angle for data of the down-orbit type; , and This represents the three-dimensional deformation rate field of the Earth's surface that needs to be inverted, where... Represents the vertical deformation rate field. Represents the north-south deformation rate field. The east-west deformation rate field represents the value that needs to be solved in the three-dimensional deformation inversion equation of the Earth's surface. Among these, the east-west deformation rate field... and vertical deformation rate field It is based on the stitched LOS imagery covering the research area. and LOS imagery The obtained north-south deformation rate field It is based on the spatially interpolated GNSS north-south deformation data The three sets of deformation rate field data obtained from the solution together constitute the three-dimensional deformation data of the Earth's surface.

8. The surface deformation inversion method based on overlapping region stitched InSAR deformation field according to claim 7, characterized in that, Step S7 specifically includes the following steps: Step S71: Calculate the horizontal strain rate tensor based on the three-dimensional surface deformation data obtained in step S6. The calculation formula is as follows: ， In the formula, Represents the horizontal strain rate tensor. , , , This represents the velocity gradient in the horizontal direction, where... This represents the rate of change of the east-west velocity with respect to the X-axis coordinate. This represents the rate of change of the north-south velocity with respect to the Y-axis coordinate. = , is the average effect of the change of east-west velocity in the north-south direction and the change of north-south velocity in the east-west direction, reflecting the shear deformation of the object; x represents the distance between two adjacent pixels in the east-west direction in the horizontal direction, and y represents the distance between two adjacent pixels in the north-south direction in the horizontal direction. This represents the partial derivative of the east-west deformation rate with respect to the X-axis. This represents the partial derivative of the east-west deformation rate with respect to the Y-axis. This represents the partial derivative of the north-south deformation rate with respect to the X-axis. The partial derivative of the north-south deformation rate with respect to the Y-axis is given. Step S72: Based on the velocity gradients in each horizontal direction of the horizontal strain rate tensor obtained in step S71, calculate the horizontal expansion rate, shear rate, and the second-order invariants of the horizontal strain rate tensor. The calculation formulas are as follows: ， In the formula, Represents the horizontal expansion rate. Represents the shear rate. The second-order invariant representing the horizontal strain rate tensor; Second-order invariant data of horizontal expansion rate, shear rate, and horizontal strain rate tensors can be used for crustal deformation analysis and geological engineering assessment.

Images