InSAR (Interferometric Synthetic Aperture Radar) deformation inversion method for lifting rail offset correction

Abstract

The invention discloses an InSAR (Interferometric Synthetic Aperture Radar) deformation inversion method for lifting rail offset correction, and mainly relates to the technical field of surface deformation monitoring. In order to solve the problem that the deformation monitoring quantity of the existing multi-dimensional small baseline set formed aperture radar interference is obviously lower than the actual deformation quantity, a time sequence deformation observation equation is established by using multi-orbit synthetic aperture radar data, an ascending-orbit small baseline set InSAR settlement monitoring point cloud and a descending-orbit small baseline set InSAR settlement monitoring point cloud are calculated respectively, point cloud registration is carried out, and the settlement monitoring point cloud of the ascending-orbit small baseline set InSAR settlement monitoring point cloud is calculated to obtain the settlement monitoring point cloud of the descending-orbit small baseline set InSAR settlement monitoring point cloud. Calculating the offset of the monitoring result of the lifting rail in the vertical, east-west and south-north directions; calculating to obtain an ascending orbit incidence angle correction number and a descending orbit incidence angle correction number of the synthetic aperture radar image by using the offsets of the ascending orbit monitoring result in the vertical direction, the east-west direction and the south-north direction; and finally, constructing a multi-dimensional small baseline set InSAR deformation inversion model containing lifting rail offset correction, solving to obtain a corrected deformation monitoring quantity, and improving the deformation monitoring precision.

Description

Technical Field

[0001] The present invention relates to the field of monitoring technologies, mainly to the field of surface three-dimensional deformation monitoring technologies, and specifically refers to an InSAR deformation inversion method for ascending / descending orbit offset correction. Background Art

[0002] The multi-dimensional deformation time series estimation technology that fuses multi-orbit SAR data has been widely applied in disasters such as earthquakes, volcanoes, landslides, and mining area subsidence, and can effectively detect surface three-dimensional vertical, east-west, and north-south deformation information, well making up for the defect of the acquisition of north-south deformation information due to the single-orbit flight mode. The Interferometric Synthetic Aperture Radar (InSAR) technology, with its advantages of all-weather, all-day, high resolution, and continuous spatial coverage, shows great potential in geological disaster monitoring, surface deformation research, etc. However, due to the problem of line-of-sight ambiguity in the monitoring process of multi-dimensional small baseline set synthetic aperture radar interferometry, the deformation monitoring quantity of multi-dimensional small baseline set synthetic aperture radar interferometry is significantly lower than the actual deformation quantity. There is an urgent need for an InSAR deformation inversion method that can solve the problems that the settlement information of the study area monitored by single-orbit SAR data shifts towards the line of sight and towards the proximal end, and the deformation quantity of multi-dimensional small baseline set synthetic aperture radar interferometry is significantly lower than the actual deformation quantity. Summary of the Invention

[0003] Aiming at the problems that the settlement information of the study area monitored by existing single-orbit SAR data shifts towards the line of sight and towards the proximal end, and the deformation quantity of multi-dimensional small baseline set synthetic aperture radar interferometry is significantly lower than the actual deformation quantity, the present invention provides an InSAR deformation inversion method for ascending / descending orbit offset correction, which performs point cloud registration on the ascending orbit small baseline set InSAR settlement monitoring point cloud and the descending orbit small baseline set InSAR settlement monitoring point cloud to calculate the offset between the two point clouds, calculates the ascending / descending orbit incident angle correction number according to the offset, and further constructs a multi-dimensional small baseline set InSAR deformation inversion model including ascending / descending orbit offset correction, improving the accuracy and accuracy of the settlement change result of the monitoring area.

[0004] An InSAR deformation inversion method for ascending / descending orbit offset correction includes:

[0005] S1. Obtain synthetic aperture radar (SAR) image data of the same area at different times, including ascending orbit SAR image data and descending orbit SAR image data, perform interferometric processing on the ascending orbit SAR image data and the descending orbit SAR image data to obtain interferometric pairs, and establish a time series deformation observation equation;

[0006] S2. Using the small baseline subset (SBAS) inversion algorithm, separately solve the ascending-track SAR image data and descending-track SAR image data after interferometric processing to obtain the ascending and descending-track monitoring results, including the ascending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and the descending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud;

[0007] S3. Perform point cloud registration on the ascending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and the descending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud, and calculate the offsets of the ascending and descending-track monitoring results in the vertical, east-west, and north-south directions;

[0008] S4. Based on the offsets of the ascending and descending-track monitoring results in the vertical, east-west, and north-south directions, calculate the ascending-track incident angle correction number and the descending-track incident angle correction number of the synthetic aperture radar (SAR) image;

[0009] S5. Based on the ascending-track incident angle correction number and the descending-track incident angle correction number, construct a multi-dimensional small baseline subset InSAR deformation inversion model including ascending and descending-track offset corrections, and use the model to solve for the deformation monitoring quantity after ascending and descending-track offset corrections.

[0010] S1 includes:

[0011] S1.1, Obtain N + 1 ascending-track synthetic aperture radar (SAR) image data covering the same area at times t0, t1,..., t N to generate M a interferometric pairs, calculate the phase difference between each interferometric pair, and obtain the linear equation of the ascending-track average deformation rate phase: a

[0012]

[0013] where K a is the coefficient matrix of the ascending-track SAR image data, V a is the estimated value of the ascending-track average deformation rate phase, ε a is the error coefficient of the ascending-track differential interferometric deformation phase matrix, is the ascending-track differential interferometric deformation phase matrix;

[0014] S1.2, Obtain N + 1 descending-track synthetic aperture radar (SAR) image data covering the same area at times t0, t1,..., t N to generate M d interferometric pairs, calculate the phase difference between each interferometric pair, and obtain the linear equation of the descending-track average deformation rate phase: d

[0015] ​​

[0016] Among them, K d is the coefficient matrix of the descending orbit SAR image data, V d is the estimated value of the average deformation rate phase of the descending orbit, and ε d is the error coefficient of the differential interferometric deformation phase matrix of the descending orbit, and is the differential interferometric deformation phase matrix of the descending orbit.

[0017] S2 includes:

[0018] S2.1. Using the small baseline subset (SBAS) inversion algorithm, the ascending orbit SAR image data after interference processing is solved to obtain the ascending orbit small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the rank of the coefficient matrix K a of the ascending orbit SAR image data is N a1 , and M a ≥N a1 , the solution of the estimated value V a of the average deformation rate phase of the ascending orbit based on the least squares method is:

[0019]

[0020] When the coefficient matrix K a of the ascending orbit SAR image data is not full rank, perform singular value decomposition on K a to obtain:

[0021] K a =U a ∑ a X a t ;

[0022] Among them, ∑ a represents a diagonal matrix of size [M a ×N a , and the diagonal elements are the eigenvalues λ a of K a K a T . T represents the transpose of the matrix, () -1 represents the inverse matrix, and U a , X a are orthogonal matrices composed of the eigenvectors obtained from K a K a T and K a T a K a ;

[0023] Based on the coefficient matrix K a of the ascending orbit SAR image data after singular value decomposition, the estimated value V a of the average deformation rate phase of the ascending orbit is calculated as:

[0024]

[0025] Integrate the ascending orbit average deformation rate phase estimate V a to obtain the settlement amount M of the ascending orbit monitoring result in the vertical direction A :

[0026]

[0027] where θ a is the incident angle of the ascending orbit SAR image, and is the deformation amount in the line-of-sight direction of the ascending orbit SAR image.

[0028] S2 also includes:

[0029] S2.2. Use the small baseline subset SBAS inversion algorithm to solve the descending orbit SAR image data after interference processing to obtain the descending orbit small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the rank of the coefficient matrix K d of the descending orbit SAR image data is N d1 and M d ≥N d1 , the solution of the descending orbit average deformation rate phase estimate V d based on the least squares method is:

[0030]

[0031] When the coefficient matrix K d of the descending orbit SAR image data is not full rank, perform singular value decomposition on K d to obtain:

[0032] K d =U d ∑ d X d T ;

[0033] where ∑ d represents a diagonal matrix of size [M d ×N d , the diagonal elements are the eigenvalues λ d K d T of K d , T represents the transpose of the matrix, () -1 represents the inverse matrix, and U d , X d are the left singular vectors and right singular vectors of K d K d T , K d T Kd The orthogonal matrix composed of the obtained eigenvectors;

[0034] Based on the coefficient matrix K of the descending-track SAR image data after singular value decomposition d , calculate the descending-track average deformation rate phase estimate V d as:

[0035]

[0036] Integrate the descending-track average deformation rate phase estimate V d to obtain the settlement amount M in the vertical direction of the descending-track monitoring result D :

[0037]

[0038] where θ d is the incident angle of the descending-track SAR image, and is the deformation amount in the line-of-sight direction of the descending-track SAR image.

[0039] In S3, let the ascending-track small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud be the source point cloud P a , and the descending-track small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud be the reference point cloud P d . Through rigid translation transformation of the ascending-track and descending-track settlement monitoring point clouds, point cloud registration is achieved, and we get:

[0040] P a =P b +Δ p ;

[0041] where Δ p represents the translation matrix, and Δ p =[Δ x ,Δ y ,Δ z , and Δ x , Δ y and Δ z respectively represent the translation amounts of the source point cloud P a in the x, y, and z coordinate axes;

[0042] Calculate the centroids N a of the source point cloud P d and the reference point cloud P a 、N d , and we get:

[0043]

[0044] where n is the number of corresponding point pairs, They respectively represent the corresponding points of the $i$-th pair in the source point cloud and the reference point cloud. The translation matrix $\Delta$ is obtained by solving the centroids of the source point cloud and the reference point cloud. p , and the expression is:

[0045] $\Delta$ p $= \overline{N}$ a $- \overline{N'}$ d ;

[0046] Based on the solved rigid transformation translation matrix $\Delta$ p , rough alignment of the point cloud is performed. The positional relationships of the corresponding points in the source point cloud and the reference point cloud are compared to identify the overlapping and non-overlapping regions of the source point cloud and the reference point cloud. The part of the source point cloud outside the overlapping region is cropped, and the translation transformation parameters are calculated.

[0047] According to the solved rigid transformation translation matrix, an objective function is constructed to obtain:

[0048]

[0049] where are the translation transformation parameters of the rigid transformation, is the data point in the source point cloud for the $i$-th pair of correspondences, $k$ is the number of iterations. The minimum value of the objective function is taken to calculate the translation transformation parameters of the rigid transformation. According to the translation transformation parameters of the rigid transformation, a rigid transformation is performed to obtain a new source point cloud ($P'$ a ). k+1 It is:

[0050]

[0051] The mean square error after the rigid transformation is calculated, and its expression is:

[0052]

[0053] where $d$ k , $d$ k+1 respectively represent the mean square errors of the $k$-th and $(k + 1)$-th iterations. A threshold $\varepsilon$ is set. When $d$ k+1 $- d$ k $< \varepsilon$, or when the maximum number of iterations $k$ max is satisfied, the iteration stops. Otherwise, the iteration process is repeated, and the translation transformation parameters obtained from the last iteration are used as the transformation matrix $\Delta$ p $= [\Delta$ x , $\Delta$ y , $\Delta$ z . Based on the transformation matrix, the spatial position of the source point cloud is moved into the space of the reference point cloud to complete the registration of the source point cloud and the reference point cloud. At the same time, the offsets of the ascending and descending orbit monitoring results in the east-west, north-south, and vertical directions are obtained, that is, the translation amounts of the source point cloud $P$ a on the $x$, $y$, and $z$ coordinate axes.

[0054] In S4, the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ of the SAR image is:

[0055]

[0056] Among them, R is the distance from the SAR satellite radar to the ground monitoring point, θ is the incident angle before correction of the SAR image, that is, the incident angle of the SAR image, θ′ is the incident angle of the image after correction of the SAR image, that is, the incident angle correction of the SAR image. According to the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ of the SAR image, based on the offsets in the vertical, east-west, and north-south directions of the ascending and descending orbit monitoring results, the ascending orbit incident angle correction number θ′ of each corresponding point is calculated. a :

[0057]

[0058] Among them, (θ′ a ) i represents the ascending orbit incident angle correction number of the i-th pair of corresponding points;

[0059] According to the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ, based on the offsets in the vertical, east-west, and north-south directions of the ascending and descending orbit monitoring results, the descending orbit incident angle correction number θ′ of each corresponding point is calculated. d :

[0060]

[0061] Among them, (θ′ d ) i represents the descending orbit incident angle correction number of the i-th pair of corresponding points.

[0062] In S5, based on the ascending orbit incident angle correction number and the descending orbit incident angle offset correction number, a multi-dimensional small baseline set InSAR deformation inversion model including ascending and descending orbit offset corrections is constructed. The relationship between the line-of-sight deformation amount and the vertical, east-west, and north-south deformation amounts of the ground surface is:

[0063]

[0064] Among them, D los is the line-of-sight deformation amount, M p is the vertical deformation amount, M E is the east-west deformation amount, M N is the north-south deformation amount, θ′ is the incident angle of the image after SAR image offset correction, β is the sensor azimuth angle of the SAR satellite. According to the differential interferometric phase diagrams of the ascending and descending orbits, a linear equation of the average phase deformation rate at each interpolation time node is obtained:

[0065] CV ENP =μφ;

[0066] Where C is the coefficient matrix of the average phase deformation rate, V ENP =[V E , V N , V P T is the deformation rate matrix in the east-west, north-south, and vertical directions at each moment to be solved, V E is the rate vector of the average phase deformation rate in the east-west direction, V N is the rate vector of the average phase deformation rate in the north-south direction, V P is the rate vector of the average phase deformation rate in the vertical direction, μφ = [μφ a μφ d T is the deformation phase vector, V is the deformation phase coefficient, φ a is the true deformation phase vector after unwrapping the ascending-track interferogram, φ d is the true deformation phase vector after unwrapping the descending-track interferogram.

[0067] By solving the linear equation of the average phase deformation rate, the surface deformation rates in the east-west, north-south, and vertical directions are obtained. Assuming that the single-look complex (SLC) image after fusing the ascending and descending tracks has T interpolated time nodes, the number of ascending-track differential interferograms and the number of descending-track differential interferograms for the solution are N A and N D , respectively. Then, the linear equation of the average phase deformation rate is represented in matrix form as:

[0068]

[0069] Where, is the deformation phase vector corresponding to the j a th ascending-track differential interferogram; is the deformation phase vector corresponding to the j d th descending-track differential interferogram; V represents the average deformation rate; is the average deformation rate in the east-west direction corresponding to the tth time node; is the average deformation rate in the north-south direction corresponding to the tth time node; is the average deformation rate in the vertical direction corresponding to the tth time node; is the deformation amount in the east-west direction of the ascending-track settlement point cloud; is the deformation amount in the north-south direction of the ascending-track settlement point cloud; is the deformation amount in the vertical direction of the ascending-track settlement point cloud; is the deformation amount in the east-west direction of the descending-track settlement point cloud; ​​is the deformation amount of the descending orbit settlement point cloud in the north-south direction; is the deformation amount of the descending orbit settlement point cloud in the vertical direction; Δt1, Δt2, and Δt3 represent time intervals;

[0070] W = [W E W N W P = [-cosβsinθ′ sinβsinθ′ cosθ′];

[0071] where W is the deformation amount in the line-of-sight direction of the SAR image, W E is the deformation amount in the east-west direction of the SAR image, W N is the deformation amount in the north-south direction of the SAR image, W P is the deformation amount in the vertical direction of the SAR image.

[0072] Calculate the deformation parameters using the singular value decomposition method, then calculate the average deformation rate V and perform integral processing on the average deformation rate V to calculate the cumulative settlement amounts M in the east-west, north-south, and vertical directions with the offset correction value added E M N M P :

[0073]

[0074] where N SAR is the total number of SAR images, represents the average deformation rate in the east-west direction of the i SAR th SAR image, represents the average deformation rate in the north-south direction of the i SAR th SAR image, represents the average deformation rate in the vertical direction of the i SAR th SAR image, represents the time interval corresponding to the i SAr th SAR image.

[0075] Compared with the prior art, the present invention has the following beneficial effects:

[0076] The present invention calculates the offset between the point cloud settlement results of the ascending orbit SAR image and the descending orbit SAR image through point cloud registration, calculates the incident angle of the ascending orbit SAR image and the incident angle correction value of the descending orbit SAR image based on the offset between the point cloud settlement results of the ascending orbit SAR image and the descending orbit SAR image, constructs a multi-dimensional small baseline set InSAR observation equation containing the ascending and descending orbit offset correction and performs the solution to obtain a more accurate and reliable monitoring area settlement change result. Description of the Drawings

[0077] Figure 1 Flowchart of the multi-dimensional small baseline set InSAR deformation inversion method including ascending / descending orbit offset correction provided by an embodiment of the present invention;

[0078] Figure 2 Comparison chart of ascending orbit SAR image settlement monitoring results and true values provided by an embodiment of the present invention;

[0079] Figure 3 Comparison chart of descending orbit SAR image settlement monitoring results and true values provided by an embodiment of the present invention;

[0080] Figure 4 Comparison chart of multi-dimensional small baseline set InSAR settlement monitoring results before correction and true values provided by an embodiment of the present invention;

[0081] Figure 5 Comparison chart of multi-dimensional small baseline set InSAR settlement monitoring results after correction and true values provided by an embodiment of the present invention. Detailed implementation manners

[0082] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions in the present invention will be clearly and completely described below. Apparently, the described embodiments are some but not all of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present invention without creative efforts shall fall within the scope of protection of the present invention.

[0083] As Figure 1 shown, an InSAR deformation inversion method for ascending / descending orbit offset correction includes:

[0084] S1. Obtain synthetic aperture radar (SAR) image data of the same area at different times, including ascending orbit SAR image data and descending orbit SAR image data, perform interferometric processing on the ascending orbit SAR image data and the descending orbit SAR image data to obtain interferometric pairs, and establish a time series deformation observation equation;

[0085] S1.1, obtain N N +1 ascending orbit synthetic aperture radar (SAR) image data covering the same area at times t0, t1,..., t a , generate M a interferometric pairs, calculate the phase difference between each interferometric pair, and obtain a linear equation for the ascending orbit average deformation rate phase:

[0086]

[0087] wherein, K a is the coefficient matrix of the ascending orbit SAR image data, V a is the estimated value of the ascending orbit average deformation rate phase, εa is the error coefficient of the ascending orbit differential interferometric deformation phase matrix, is the ascending orbit differential interferometric deformation phase matrix;

[0088] S1.2, Obtain N N +1 descending orbit synthetic aperture radar (SAR) image data covering the same area at times t0, t1,..., t d Generate M d interferometric pairs, calculate the phase difference between each interferometric pair, and obtain the linear equation of the descending orbit average deformation rate phase:

[0089]

[0090] where, K d is the coefficient matrix of the descending orbit SAR image data, V d is the estimated value of the descending orbit average deformation rate phase, ε d is the error coefficient of the descending orbit differential interferometric deformation phase matrix, is the descending orbit differential interferometric deformation phase matrix.

[0091] S2. Use the small baseline subset (SBAS) inversion algorithm to separately solve the ascending orbit SAR image data and the descending orbit SAR image data after interferometric processing, and obtain the ascending and descending orbit monitoring results, including the ascending orbit small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and the descending orbit small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud;

[0092] S2.1, Use the small baseline subset (SBAS) inversion algorithm to solve the ascending orbit SAR image data after interferometric processing, and obtain the ascending orbit small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the rank of the coefficient matrix K a of the ascending orbit SAR image data is N a1 , and M a ≥N a1 , based on the least squares method, the solution of the estimated value V a of the ascending orbit average deformation rate phase is:

[0093]

[0094] When the coefficient matrix K a of the ascending orbit SAR image data is non-full rank, perform singular value decomposition on K a to obtain:

[0095] K a = U a Σ a X a T ;

[0096] where, Σa represents a diagonal matrix of size [M a ×N a , with diagonal elements being K a K a T 's eigenvalue λ a , T represents the transpose of the matrix, and () -1 represents the inverse matrix, and U a , X a are respectively the orthogonal matrices composed of the eigenvectors obtained from K a K a T , K a T K a ;

[0097] Based on the coefficient matrix K a of the ascending-track SAR image data after singular value decomposition, the ascending-track average deformation rate phase estimate V a is calculated as:

[0098]

[0099] Integrating the ascending-track average deformation rate phase estimate V a yields the settlement amount M A in the vertical direction of the ascending-track monitoring result:

[0100]

[0101] where θ a is the incident angle of the ascending-track SAR image, and

[0102] is the deformation amount in the line-of-sight direction of the ascending-track SAR image. d S2.2. Using the small baseline subset SBAS inversion algorithm, the descending-track SAR image data after interference processing is solved to obtain the descending-track small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the rank of the coefficient matrix K d1 of the descending-track SAR image data is N d , and M d1 ≥N d , the solution of the descending-track average deformation rate phase estimate V

[0103]

[0104] When the coefficient matrix K d of the descending-track SAR image data is non-full rank, performing singular value decomposition on K d yields:

[0105] Kd = U d ∑ d X d T ;

[0106] where ∑ d represents a diagonal matrix of size [M d × N d , with diagonal elements being the eigenvalues K d λ d T of K d , T represents the transpose of the matrix, and () -1 represents the inverse matrix. U d , X d are respectively orthogonal matrices composed of the eigenvectors obtained for K d K d T , K d T K d ;

[0107] Based on the coefficient matrix K d of the descending-orbit SAR image data after singular value decomposition, the descending-orbit mean deformation rate phase estimate V d is calculated as:

[0108]

[0109] Integrating the descending-orbit mean deformation rate phase estimate V d yields the settlement amount K D in the vertical direction of the descending-orbit monitoring result:

[0110]

[0111] where θ d is the incident angle of the descending-orbit SAR image, and is the deformation amount in the line-of-sight direction of the descending-orbit SAR image.

[0112] S3. Perform point cloud registration on the ascending-orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and the descending-orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud, and calculate the offsets of the ascending and descending orbit monitoring results in the vertical, east-west, and north-south directions;

[0113] Let the ascending-orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud be the source point cloud P a , and the descending-orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud be the reference point cloud P d . By performing a rigid translation transformation on the ascending and descending orbit settlement monitoring point clouds, point cloud registration is achieved, resulting in:

[0114] P a = P b + Δ [ ;

[0115] where, Δ [ represents the translation matrix, and Δ p = [Δ x , Δ y , Δ z , where Δ x , Δ y and Δ z respectively represent the translation amounts of the source point cloud P a along the x, y, and z coordinate axes;

[0116] Calculate the centroids N a and N d of the source point cloud P a and the reference point cloud P d , respectively, to obtain:

[0117]

[0118] where n is the number of corresponding point pairs, respectively represent the i-th pair of corresponding points in the source point cloud and the reference point cloud, and the translation matrix Δ p is obtained by solving the centroids of the source point cloud and the reference point cloud. The expression is:

[0119] Δ p = N a - N d ;

[0120] Based on the solved rigid transformation translation matrix Δ p , perform rough alignment of the point cloud, compare the positional relationships of the corresponding points in the source point cloud and the reference point cloud, identify the overlapping and non-overlapping regions of the source point cloud and the reference point cloud, crop the part of the source point cloud outside the overlapping region, and calculate the translation transformation parameters.

[0121] According to the solved rigid transformation translation matrix, construct an objective function to obtain:

[0122]

[0123] where, is the translation transformation parameter of the rigid change, is the data point in the source point cloud for the i-th pair of correspondences, k is the number of iterations, take the minimum value of the objective function, calculate the translation transformation parameter of the rigid change, and perform a rigid transformation according to the translation transformation parameter of the rigid change to obtain a new source point cloud (P' a ) k+1 as:

[0124]

[0125] Calculate the mean square error after the rigid change, and its expression is:

[0126]

[0127] where d k , d k+1 respectively represent the mean square errors after the k-th and (K + 1)-th iterations. Set a threshold ε. When d k+1 - d k < ε, or when the maximum number of iterations k max is satisfied, stop the iteration. Otherwise, repeat the iteration process. Use the translation transformation parameters obtained in the last iteration as the transformation matrix Δ p = [Δ x , Δ y , Δ z . Based on the transformation matrix, move the spatial position of the source point cloud to the space of the reference point cloud to complete the registration of the source point cloud and the reference point cloud, and at the same time obtain the offsets of the ascending and descending orbit monitoring results in the east-west, north-south, and vertical directions, that is, the translation amounts of the source point cloud P a on the x, y, and z coordinate axes.

[0128] S4. Based on the offsets of the ascending and descending orbit monitoring results in the vertical, east-west, and north-south directions, calculate the ascending orbit incidence angle correction number and the descending orbit incidence angle correction number of the synthetic aperture radar (SAR) image;

[0129] The relationship between the distance from the SAR satellite radar to the ground monitoring point and the incidence angle θ of the SAR image is:

[0130]

[0131] where R is the distance from the SAR satellite radar to the ground monitoring point, θ is the incidence angle before the correction of the SAR image, that is, the incidence angle of the SAR image, and θ′ is the incidence angle of the SAR image after the correction, that is, the corrected incidence angle of the SAR image. According to the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incidence angle θ of the SAR image, based on the offsets of the ascending and descending orbit monitoring results in the vertical, east-west, and north-south directions, calculate the ascending orbit incidence angle correction number θ′ of each corresponding point:

[0132]

[0133] where (θ′ a ) i represents the ascending orbit incidence angle correction number of the i-th pair of corresponding points;

[0134] According to the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ, and based on the offsets in the vertical, east-west, and north-south directions from the ascending and descending orbit monitoring results, the descending orbit incident angle correction θ′ for each corresponding point is calculated. d :

[0135]

[0136] Among them, (θ′ d ) i represents the descending orbit incident angle correction for the i-th pair of corresponding points.

[0137] S5. Based on the ascending orbit incident angle correction and the descending orbit incident angle correction, construct a multi-dimensional small baseline set InSAR deformation inversion model that includes ascending and descending orbit offset corrections. Use the model to solve for the deformation monitoring quantity after ascending and descending orbit offset corrections. The relationships between the line-of-sight deformation quantity and the vertical, east-west, and north-south deformation quantities on the ground surface are as follows:

[0138]

[0139] Among them, D los is the line-of-sight deformation quantity, M p is the vertical deformation quantity, M E is the east-west deformation quantity, M N is the north-south deformation quantity, θ′ is the image incident angle after SAR image offset correction, β is the sensor azimuth angle of the SAR satellite. The linear equation for the average phase deformation rate at each interpolation time node is obtained from the ascending and descending orbit differential interferometric phase diagrams:

[0140] CV ENP = μφ;

[0141] Among them, C is the coefficient matrix of the average phase deformation rate, V ENP = [V E , V N , V P T is the deformation rate matrix in the east-west, north-south, and vertical directions at each moment to be solved, V E is the rate vector of the average phase deformation rate in the east-west direction, V N is the rate vector of the average phase deformation rate in the north-south direction, V P is the rate vector of the average phase deformation rate in the vertical direction, μφ = [μφ a μφ d T is the deformation phase vector, μ is the deformation phase coefficient, φ a is the true deformation phase vector after unwrapping the ascending orbit interferogram, φ d is the true deformation phase vector after unwrapping the descending orbit interferogram.

[0142] By solving the linear equation of the average phase deformation rate, the surface deformation rates in the east-west, north-south, and vertical directions are obtained. Assuming that the single-look complex (SLC) image after fusing ascending and descending orbits has obtained T interpolated time nodes, the number of ascending-orbit differential interferograms and descending-orbit differential interferograms for solving are N A and N D , respectively. Then, the linear equation of the average phase deformation rate is represented in matrix form as:

[0143]

[0144] where is the deformation phase vector corresponding to the j a -th ascending-orbit differential interferogram; is the deformation phase vector corresponding to the j d -th descending-orbit differential interferogram; V represents the average deformation rate; is the average deformation rate in the east-west direction corresponding to the t-th time node; is the average deformation rate in the north-south direction corresponding to the t-th time node; is the average deformation rate in the vertical direction corresponding to the t-th time node; is the deformation amount of the ascending-orbit settlement point cloud in the east-west direction; is the deformation amount of the ascending-orbit settlement point cloud in the north-south direction; is the deformation amount of the ascending-orbit settlement point cloud in the vertical direction; is the deformation amount of the descending-orbit settlement point cloud in the east-west direction; is the deformation amount of the descending-orbit settlement point cloud in the north-south direction; is the deformation amount of the descending-orbit settlement point cloud in the vertical direction; Δt1, Δt2, and Δt3 represent time intervals;

[0145] W = [W E W N W P = [-cosβsinθ′ sinβsinθ′ cosθ′];

[0146] where W is the deformation amount in the line-of-sight direction of the SAR image, W E is the deformation amount in the east-west direction of the SAR image, W N is the deformation amount in the north-south direction of the SAR image, and W P is the deformation amount in the vertical direction of the SAR image

[0147] Calculate the deformation parameters using the singular value decomposition method, including the surface deformation rates in the east-west, north-south, and vertical directions, the deformation amounts in the east-west, north-south, and vertical directions of the ascending-track settlement point cloud and the descending-track settlement point cloud. Based on the average phase deformation rate equation, calculate the average deformation rate V and perform integral processing on the average deformation rate V to calculate the cumulative settlement amounts M in the east-west, north-south, and vertical directions after adding the offset correction value E , M N , M P :

[0148]

[0149] Among them, N SAR is the total number of SAR images, represents the average deformation rate in the east-west direction of the i SAR -th SAR image, represents the average deformation rate in the north-south direction of the i SAR -th SAR image, represents the average deformation rate in the vertical direction of the i SAR -th SAR image, represents the time interval corresponding to the i SAR -th SAR image.

[0150] As Figure 2 and Figure 3 shown, due to the combined influence of factors such as the satellite observation angle, surface cover, and atmospheric interference under different orbits, there is a line-of-sight ambiguity problem. When the horizontal deformation in the study area is large, the difference between the vertical deformation amount value monitored by the single-orbit InSAR technology and the actual surface vertical deformation amount value will increase and there will be an offset of the deformation center. The ascending and descending-track settlement monitoring ranges shift towards the proximal end of their respective line-of-sight directions, that is, relative to the true value of the point cloud settlement, the monitoring result of the ascending-track SAR image point cloud settlement shifts to the left, i.e., westward, and the monitoring result of the descending-track SAR image point cloud settlement shifts to the right, i.e., eastward, and there are large errors between their monitoring results and the true value. In addition, the monitoring result of the ascending-track SAR image point cloud settlement shows a relatively smooth change trend at some positions, while the monitoring result of the descending-track SAR image point cloud settlement is more complex and variable.

[0151] As Figure 4 and Figure 5 shown, Figure 4Among them, the multi-dimensional small baseline set InSAR monitoring technology before correction can reflect the surface settlement to a certain extent. The deformation monitoring area is located between the ascending and descending orbit small baseline set InSAR results and basically covers the same position. The monitoring results basically coincide with the true values, and can improve the situation where the deformation areas obtained by the ascending and descending orbit small baseline set InSAR methods are offset. However, the multi-dimensional small baseline set InSAR monitoring uses SAR images from different perspectives of multiple orbits. Due to the problem of line-of-sight ambiguity, there is a problem that the deformation monitoring result quantity of the multi-dimensional small baseline set InSAR is significantly lower than the actual deformation quantity during the monitoring process in some areas. Figure 5 It is the monitoring result of the multi-dimensional small baseline set InSAR deformation inversion model including ascending and descending orbit offset correction constructed based on the embodiments of the present invention. By performing point cloud registration on the ascending orbit small baseline set InSAR settlement monitoring point cloud and the descending orbit small baseline set InSAR settlement monitoring point cloud to calculate the offset between the two point clouds, calculating the ascending and descending orbit incident angle correction numbers according to the offset, and then constructing a multi-dimensional small baseline set InSAR deformation inversion model including ascending and descending orbit offset correction, and using the model to solve for the corrected deformation monitoring quantity. While improving the situation where the deformation areas obtained by the ascending and descending orbit small baseline set InSAR methods are offset, it solves the problem that the deformation monitoring result quantity of the multi-dimensional small baseline set InSAR is significantly lower than the actual deformation quantity, improves the accuracy and accuracy of the settlement change results in the monitoring area, and makes the deformation monitoring quantity approximately close to the true value.

[0152] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. An InSAR deformation inversion method for correcting the elevation track offset, characterized in that Including: S1. Obtain synthetic aperture radar (SAR) image data at the same location but different times, including ascending-track SAR image data and descending-track SAR image data. Perform interferometric processing on the ascending-track SAR image data and the descending-track SAR image data to obtain an interferometric pair, and establish a time-series deformation observation equation. S2. Use the small baseline subset (SBAS) inversion algorithm to separately solve the interferometrically processed ascending-track SAR image data and descending-track SAR image data to obtain ascending and descending track monitoring results, including an ascending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and a descending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. S3. Perform point cloud registration on the ascending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud and the descending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud, and calculate the offsets of the ascending and descending track monitoring results in the vertical, east-west, and north-south directions. S4. Based on the offsets of the ascending and descending track monitoring results in the vertical, east-west, and north-south directions, calculate the ascending-track incident angle correction number and the descending-track incident angle correction number of the synthetic aperture radar (SAR) image. S5. Based on the ascending-track incident angle correction number and the descending-track incident angle correction number, construct a multi-dimensional small baseline subset InSAR deformation inversion model that includes ascending and descending track offset corrections, and use the model to solve for the deformation measurement quantity after ascending and descending track offset corrections.

2. The InSAR deformation inversion method for correcting the elevation track offset according to claim 1, wherein S1 includes: S1.1, obtain t0, t1,..., t N The N a +1 ascending-track synthetic aperture radar (SAR) image data covering the same area in terms of time, generate M a interferometric pairs, calculate the phase difference between each interferometric pair, and obtain the linear equation of the ascending-track average deformation rate phase: Among them, K a is the coefficient matrix of the ascending-track SAR image data, V a is the ascending-track average deformation rate phase estimation value, ε a is the error coefficient of the ascending-track differential interferometric deformation phase matrix, is the ascending-track differential interferometric deformation phase matrix; S1.2, obtain t0, t1,..., t N N d +1 descending orbit synthetic aperture radar (SAR) image data covering the same area, generate M d interferometric pairs, calculate the phase difference between each interferometric pair, and obtain the linear equation of the descending orbit mean deformation rate phase: Among them, K d is the coefficient matrix of the descending orbit SAR image data, V d is the phase estimation of the descending orbit average deformation rate, ε d is the error coefficient of the descending orbit differential interferometric deformation phase matrix, is the descending orbit differential interferometric deformation phase matrix.

3. An InSAR deformation inversion method for correcting the offset of a lifting track according to claim 2, characterized in that S2 Including: S2.1, Using the small baseline subset (SBAS) inversion algorithm, solve the ascending-track SAR image data after interference processing to obtain the ascending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the coefficient matrix K a has a rank of N a1 , and M a ≥N a1 , the solution of the ascending-track average deformation rate phase estimate V a obtained based on the least squares method is: When the coefficient matrix K of the ascending orbit SAR image data a is rank-deficient, for K a perform singular value decomposition to obtain: K a = U a ∑ a X a T ; Among them, ∑ a represents a diagonal matrix of size [M a ×N a , with diagonal elements being the eigenvalues λ a K a T of K, T represents the transpose of a matrix, () a represents the inverse matrix, and U -1 , X a are orthogonal matrices composed of the eigenvectors obtained from K a respectively; a K a T K a T K a The orthogonal matrix formed by the obtained eigenvectors; Coefficient matrix K of up-track SAR image data after singular value decomposition a , the phase estimation V of the average deformation rate of the up-track is calculated a as follows: Integrate the ascending orbit average deformation rate phase estimate V a to obtain the settlement amount M of the ascending orbit monitoring result in the vertical direction A : Among them, θ a is the incident angle of the ascending-track SAR image, and is the deformation amount in the line-of-sight direction of the ascending-track SAR image.

4. The InSAR deformation inversion method for correcting the offset of the lifting track according to claim 3, characterized in that, S2 also includes: S2.2, using the small baseline subset (SBAS) inversion algorithm, solve the descending-track SAR image data after interference processing to obtain the descending-track small baseline subset synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud. When the coefficient matrix K d has a rank of N d1 , and M d ≥N d1 , the solution of the descending-track average deformation rate phase estimate V d obtained based on the least squares method is: When the coefficient matrix K of the descending orbit SAR image data d is not of full rank, perform a singular value decomposition on K d to obtain: K d = U d ∑ d X d T ; Among them, ∑ d represents a diagonal matrix of size [M d ×N d , with diagonal elements being the eigenvalues λ d K d T of K d , T represents the transpose of the matrix, () -1 represents the inverse matrix, and U d , X d are respectively orthogonal matrices composed of the eigenvectors obtained from K d K d T , K d T K d ; Coefficient matrix K of descending orbit SAR image data after singular value decomposition d , the descending orbit average deformation rate phase estimate V d is calculated as follows: Integrate the descending orbit average deformation rate phase estimate V d to obtain the settlement amount M of the descending orbit monitoring result in the vertical direction D : Among them, θ d is the incident angle of the descending orbit SAR image, and is the deformation amount in the line-of-sight direction of the descending orbit SAR image.

5. An InSAR deformation inversion method for correcting the offset of a lifting track according to claim 4, characterized in that In S3, set the ascending orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud as the source point cloud P a , and set the descending orbit small baseline set synthetic aperture radar interferometry (InSAR) settlement monitoring point cloud as the reference point cloud P d , and perform a rigid translation transformation on the ascending and descending orbit settlement monitoring point clouds to achieve point cloud registration, obtaining: P a = P b + Δ p ; Among them, Δ p represents the translation matrix, and Δ p = [Δ x , Δ y , Δ z , where Δ x , Δ y and Δ z respectively represent the translation amounts of the source point cloud P a along the x, y, and z coordinate axes; Calculate the centroids of the source point cloud P a and the reference point cloud P d respectively, and obtain the centroids N a 、N d : where n is the number of corresponding point pairs, which respectively represent the o-th corresponding points in the source point cloud and the reference point cloud. The translation matrix Δ is obtained by solving the centroids of the source point cloud and the reference point cloud p , and the expression is: Δ p = N a - N d ; Translation matrix Δ of rigid transformation based on solution p Perform rough alignment of the point cloud, compare the positional relationships of corresponding points in the source point cloud and the reference point cloud, identify the overlapping and non-overlapping regions between the source point cloud and the reference point cloud, crop the part of the source point cloud outside the overlapping region, and calculate the translation transformation parameters.

6. The InSAR deformation inversion method for correcting the offset of the lifting track according to claim 5, characterized in that, According to the solved rigid transformation translation matrix, construct an objective function to obtain: Among them, is the translation transformation parameter with rigid change, is the data point in the source point cloud corresponding to the i-th pair, k is the number of iterations, taking the minimum value of the objective function, calculating the translation transformation parameter with rigid change, and performing a rigid transformation according to the translation transformation parameter with rigid change to obtain a new source point cloud (P′ a ) k+1 is: Calculate the mean square error after rigid change, and its expression is: where d k and d k+1 respectively represent the mean square errors of the k-th and (k + 1)-th iterations. A threshold ε is set. When d k+1 - d k < ε, or when the maximum number of iterations k max is reached, the iteration is stopped. Otherwise, the iteration process is repeated, and the translation transformation parameters obtained from the last iteration are used as the transformation matrix Δ p = [Δ x , Δ y , Δ z . Based on the transformation matrix, the spatial position of the source point cloud is moved into the space of the reference point cloud to complete the registration of the source point cloud and the reference point cloud. At the same time, the offsets of the ascending and descending orbit monitoring results in the east-west, north-south, and vertical directions are obtained, that is, the translation amounts of the source point cloud P a on the x, y, and z coordinate axes.

7. An InSAR deformation inversion method for correcting the offset of a lifting track according to claim 6, characterized in that, In S4, the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ of the SAR image is: Among them, R is the distance from the SAR satellite radar to the ground monitoring point, θ is the incident angle before the SAR image correction, that is, the incident angle of the SAR image, and θ′ is the incident angle of the image after the SAR image correction, that is, the SAR image correction incident angle. According to the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ of the SAR image, and based on the offset in the vertical, east-west, and north-south directions of the ascending and descending orbit monitoring results, the ascending orbit incident angle correction number θ′ of each corresponding point is calculated. a : Among them, (θ′ a ) i represents the ascending orbit incidence angle correction of the corresponding points of the i-th pair; Based on the relationship between the distance from the SAR satellite radar to the ground monitoring point and the incident angle θ, and the offsets in the vertical, east-west, and north-south directions based on the ascending and descending orbit monitoring results, the descending orbit incident angle correction θ′ for each corresponding point is calculated. d : Among them, (θ′ d ) i represents the correction of the descending orbit incidence angle of the corresponding points of the i-th pair.

8. An InSAR deformation inversion method for correcting elevation track offset according to claim 7, characterized in that, In S5, based on the ascending-track incident angle correction number and the descending-track incident angle offset correction number, construct a multi-dimensional small baseline subset InSAR deformation inversion model that includes ascending and descending track offset corrections. The relationship between the line-of-sight deformation quantity and the vertical, east-west, and north-south deformation quantities on the ground surface is: Among them, D los is the line-of-sight deformation variable, M p is the vertical deformation variable, M E is the east-west deformation variable, M N is the north-south deformation variable, θ′ is the image incidence angle after SAR image offset correction, β is the sensor azimuth angle of the SAR satellite. The linear equation of the average phase deformation rate at each interpolation time node is obtained from the differential interferometric phase diagrams of ascending and descending orbits: CV ENP =μφ; Among them, C is the coefficient matrix of the average phase deformation rate, and V ENP =[V E , V N , V P T is the deformation rate matrix in the east-west, north-south, and vertical directions at each moment to be solved, V E is the rate vector of the average phase deformation rate in the east-west direction, V N is the rate vector of the average phase deformation rate in the north-south direction, V P is the rate vector of the average phase deformation rate in the vertical direction, μφ = [μφ a μφ d T is the deformation phase vector, μ is the deformation phase coefficient, and φ a is the true deformation phase vector after unwrapping the ascending orbit interferogram, and φ d is the true deformation phase vector after unwrapping the descending orbit interferogram.​​ 9. The InSAR deformation inversion method for correcting the offset of the lifting track according to claim 8, characterized in that By solving the linear equation of the average phase deformation rate, the surface deformation rates in the east-west, north-south, and vertical directions are obtained. Assuming that the single-look complex (SLC) images after the combination of ascending and descending orbits have obtained T interpolated time nodes, the number of ascending-orbit differential interferograms and the number of descending-orbit differential interferograms for the solution are N A and N D , respectively. Then, the linear equation of the average phase deformation rate is expressed in matrix form as: Among them, is the deformation phase vector corresponding to the j a th ascending orbit differential interferogram; is the deformation phase vector corresponding to the j d th descending orbit differential interferogram; V represents the average deformation rate; is the average deformation rate in the east-west direction corresponding to the t-th time node; is the average deformation rate in the north-south direction corresponding to the t-th time node; is the average deformation rate in the vertical direction corresponding to the t-th time node; is the deformation amount of the ascending orbit settlement point cloud in the east-west direction; is the deformation amount of the ascending orbit settlement point cloud in the north-south direction; is the deformation amount of the ascending orbit settlement point cloud in the vertical direction; is the deformation amount of the descending orbit settlement point cloud in the east-west direction; is the deformation amount of the descending orbit settlement point cloud in the north-south direction; is the deformation amount of the descending orbit settlement point cloud in the vertical direction; Δt1, Δt2, Δt3 represent time intervals; W = [W E W N W P = [-cosβsinθ′ sinβsinθ′ cosθ′]; Among them, \(W\) is the deformation in the line-of-sight direction of the SAR image, \(W\) E is the deformation in the east-west direction of the SAR image, \(W\) N is the deformation in the north-south direction of the SAR image, \(W\) P is the vertical deformation of the SAR image.

10. An InSAR deformation inversion method for correcting the offset of a lifting track according to claim 9, characterized in that Calculate the deformation parameters using the singular value decomposition method, and then solve for the average deformation rate V. Integrate the average deformation rate V to calculate the cumulative settlement amounts M in the east-west, north-south, and vertical directions with the offset correction value added E , M N , M P : Among them, N SAR is the total number of SAR images, represents the average deformation rate in the east-west direction of the i SAR -th SAR image, represents the average deformation rate in the north-south direction of the i SAR -th SAR image, represents the average deformation rate in the vertical direction of the i SAR -th SAR image, represents the time interval corresponding to the i SAR -th SAR image.

Images