SBAS-InSAR earth surface deformation inversion method based on mixed regularization improvement

Abstract

The invention discloses an improved SBAS-InSAR earth surface deformation inversion method based on hybrid regularization, and relates to the field of remote sensing image processing of synthetic aperture radar. The invention aims to solve the problem that the inversion result is unstable in the existing terrain inversion method. The method comprises the following steps: acquiring an SAR image, and preprocessing the SAR image to obtain a preprocessed SAR image; threshold values of a time baseline and a space baseline are set, interference pairs are obtained by using the preprocessed SAR image and the threshold values of the time baseline and the space baseline, and an original interferogram is generated; interference phases in the original interferogram are removed, phase unwrapping processing is carried out, and real differential interferometric phases are obtained; and constructing a design matrix, constructing a relationship between a phase change rate vector and an observation vector by using the design matrix and a real differential interference phase, constructing an objective function by using a hybrid regularization method and the relationship between the phase change rate vector and the observation vector, and solving the objective function to obtain a deformation rate, thereby completing surface deformation inversion. The method is used for surface deformation inversion.

Description

Technical Field

[0001] The present invention relates to the field of synthetic aperture radar remote sensing image processing, in particular to an improved SBAS-InSAR surface deformation inversion method based on hybrid regularization. Background Art

[0002] Surface deformation is an environmental geological phenomenon caused by compression of the Earth's crust, resulting in regional changes in surface elevation. This can cause permanent damage to environmental resources and human livelihoods. InSAR technology has been widely used in surface deformation monitoring. The acquisition of three-dimensional ground information effectively supports geological disaster monitoring and prevention. It also plays a vital role in predicting urban land subsidence and mining collapse.

[0003] SBAS-InSAR technology, proposed by Berardion et al. in 2002, improves the overall coherence of the interferometer pair by shortening both the temporal and spatial baselines, resulting in more accurate target point results. It plays a crucial role in acquiring ground deformation information. The introduction and development of SBAS-InSAR technology has further promoted the widespread application of remote sensing in surface monitoring and has become a focus of research on surface deformation in recent years. Using singular value decomposition and SBAS-InSAR, Berardion et al. used 44 SAR images to obtain ground deformation data for a city and a volcanic crater from 1992 to 2000. These data were consistent with GPS monitoring results, demonstrating the effectiveness of the SBAS-InSAR algorithm. In 2004, Riccardo Lanari improved Berardion's method, accounting for the influence of different phases on the monitoring results. He estimated the low-pass deformation phase from multi-look processed data, incorporated DEM error phase and atmospheric delay phase, and used this improved method to monitor surface deformation in the same area, with the results consistent with those from leveling. In 2008, Casu et al. used SBAS-InSAR technology to monitor surface deformation in a region, obtaining data on the evolution of land subsidence and its average rate from 1992 to 2000. In 2011, Pepe et al. first used SBAS-InSAR technology to monitor deformation of a volcanic cluster using ENVISAT ASAR imagery in both StripMap and ScanSAR modes. This combined processing of StripMap and ScanSAR imagery for the first time improved the temporal resolution of SBAS monitoring and laid the foundation for the combined processing of multi-angle and multi-mode data. In 2014, Calò et al. proposed the P-SBAS-InSAR technique based on a two-layer parallel approach. Using X- and C-band SAR imagery, they used the short baseline InSAR method to monitor landslide subsidence in the Ivancich region of Italy from 1992 to 2010. In 2023, Ma Chuang and others used SBAS-InSAR technology to monitor landslides in Baige, Polo Township, Jiangda County, Tibet Autonomous Region, based on 23 C-band Sentinel-1A ascending-orbit SAR images. They also used a velocity inverse model to predict the time of the disaster. The traditional short-baseline SBAS-InSAR method mentioned above can effectively utilize redundant interferometer pairs to obtain interferometric phases, overcoming the shortcomings of conventional D-InSAR technology, such as severe decorrelation and low deformation detection accuracy. It is highly adaptable to atmospheric, topographic, and surface changes. However, when these methods are faced with an insufficient number of interferometer pairs that meet the requirements, severe noise interference, or overlapping time baselines, the least squares solution will be unstable due to the rank deficiency of the design matrix, which in turn makes the surface deformation inversion results unstable. Summary of the Invention

[0004] The purpose of the present invention is to solve the problem of unstable inversion results in existing terrain inversion methods, and propose an improved SBAS-InSAR surface deformation inversion method based on hybrid regularization.

[0005] The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization is as follows:

[0006] Step 1: Acquire SAR images and preprocess the SAR images to obtain preprocessed SAR images;

[0007] Step 2: Set the thresholds of the time baseline and the space baseline, obtain the interferometric pairs using the pre-processed SAR image and the thresholds of the time baseline and the space baseline, and generate the original interferogram;

[0008] Step 3: Remove the interference phase in the original interferogram and perform phase unwrapping to obtain the true differential interferometry phase;

[0009] Step 4: Construct a design matrix, use the design matrix and the true differential interferometer phase to construct the relationship between the phase change rate vector and the observation vector, and use the hybrid regularization method and the relationship between the phase change rate vector and the observation vector to construct the objective function. Solve the objective function to obtain the optimal phase change rate vector, thereby completing the surface deformation inversion.

[0010] Furthermore, the step 1 of acquiring the SAR image and preprocessing the SAR image to obtain the preprocessed SAR image is as follows:

[0011] Step 1: Acquire SAR images, convert the SAR images into single-look complex format, and then eliminate the phase deviation caused by orbit error;

[0012] Step 1 and 2: Crop the SAR image according to the study area to obtain the preprocessed SAR image.

[0013] Furthermore, in step 2, the thresholds of the time baseline and the space baseline are set, and the interferogram is obtained by using the pre-processed SAR image and the thresholds of the time baseline and the space baseline, and generating the original interferogram, specifically:

[0014] Step 21: At imaging time t0~t N Acquire N+1 SAR images for the same observation area, select one of the SAR images as the super master image, register the other SAR images with the super master image, and screen and separate the registered SAR images to obtain M interferometric image pairs.

[0015] Step 22: Calculate the phase Φ of each SAR image in the interferometric image pair and the differential interferometric phase ΔΦ of each interferometric image pair, thereby generating an original interferogram.

[0016] Furthermore, the step 21 of screening and sorting the registered SAR images is as follows:

[0017] First, all pre-processed SAR images are paired into interferometric image pairs;

[0018] Then, the thresholds of the time baseline and the space baseline are set, and the interference image pairs that are smaller than the thresholds of the time baseline and the space baseline are considered to meet the requirements;

[0019] The number of interference image pairs that meet the requirements is M, and M meets the following conditions:

[0020]

[0021] Among them, N+1 is the total number of SAR images in the current observation area.

[0022] Furthermore, in step 3, the interference phase in the original interference pattern is removed and phase unwrapping is performed to obtain the true differential interference phase, specifically:

[0023] Firstly, the differential interferometric phase of the original interferogram is obtained by using the dichotomy method, and then the external DEM data is introduced to remove the interference of the terrain phase.

[0024] Then, the flat-ground effect removal method based on frequency shift estimation is used to remove the flat-ground phase in the original interferogram.

[0025] Finally, the phase unwrapping method based on the minimum cost flow method of the network model is used to unwrap the phase of the original interferogram, and finally the true differential interferometry phase is obtained.

[0026] Furthermore, the design matrix is ​​constructed in step 4, and the relationship between the phase change rate vector and the observation vector is constructed using the design matrix and the true differential interferometry phase. The objective function is constructed using the hybrid regularization method and the relationship between the phase change rate vector and the observation vector. The objective function is solved to obtain the deformation rate, thereby completing the surface deformation inversion, specifically:

[0027] Step 41: Construct the design matrix A using the real differential interferometer phase;

[0028] Step 42: Use the design matrix and the true differential interferometer phase to construct the relationship between the phase change rate vector and the observation vector;

[0029] Step 4.3: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, and use the L curve to obtain the optimized regularization parameter λ opt, construct the objective function using the regularization parameter;

[0030] Step 44: Solve the objective function and obtain the optimal phase change rate to complete the surface deformation inversion.

[0031] Furthermore, the design matrix A is constructed using the real differential interference phase in step 41, specifically:

[0032] First, one SAR image in each interferometric image pair obtained in step 2 is used as the primary image, and the other SAR image is used as the auxiliary image;

[0033] Then, the design matrix A with M rows and N+1 columns is constructed using the main image and the auxiliary image;

[0034] The positions of "1" and "-1" in the row vector of the design matrix A represent the positions of interference pairs in an original interference pattern, and the column vector represents the SAR image at a certain moment;

[0035] In the matrix A, 1 represents the primary image, -1 represents the secondary image, and other positions are set to 0.

[0036] Furthermore, the relationship between the phase change rate vector and the observation vector is constructed by using the design matrix and the true differential interference phase in step 42, specifically:

[0037] Step 421: Use the true differential interferometry phase to construct the phase vector and observation vector, specifically:

[0038] Φ T =[φ(t0),φ(t1),…,φ(t N )]

[0039] ΔΦ T =[δφ1,δφ2,…,δφ M ]

[0040] Among them, φ(t N ) is the time t N The phase value on T is the phase vector, ΔΦ T is the observation vector, δφ M is the true differential interferometric phase of the Mth interferometric pair;

[0041] Step 422: Establish the relationship between the phase vector and the observation vector, specifically:

[0042] ΔΦ=AΦ

[0043]

[0044] IS=[IS1,IS2,…,IS M ]

[0045] IM=[IM1,IM2,…,IM M ]

[0046] Where ΔΦ j is the jth element in ΔΦ, Is the main image IM j The corresponding moment, Auxiliary image IS j The corresponding moment, Is the main image IM j The phase, It is auxiliary image IM j phase, IS is the auxiliary image label sequence, IM is the main image label sequence, j∈[1,M], IM j IS is the original label of the jth main image after sorting by time, j is the original label of the jth auxiliary image after sorting by time, IM j >IS j ;

[0047] Step 423: Use the relationship between the phase vector and the observation vector to construct the relationship between the phase change rate vector and the observation vector, specifically:

[0048] BV=ΔΦ

[0049] B(j,k)=t k+1 -t k (IS j +1<k<IM j )

[0050] V T =[v1,v2,…,v N ]

[0051]

[0052] Among them, δφ j is the jth element in ΔΦ, k is [IS j +1,IM j ], B is an M×N matrix, B(j,k) is the value of the jth row and kth column in B, t k+1 is the k+1th moment, t k is the kth moment, k is the time label, v k is the kth phase change rate, l∈[1,N], V T is the phase change rate vector, v N is the Nth phase change rate, l is the phase change rate index, v l is the lth phase change rate, φ(t l) is the time t l The phase value on φ(t l-1 ) is the time t l-1 The phase value on .

[0053] Furthermore, in step 43, an L curve is constructed based on the relationship between the phase change rate vector and the observation vector, and the L curve is used to obtain the optimized regularization parameter λ opt , the objective function is constructed using the regularization parameter, specifically:

[0054] Step 431: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, specifically:

[0055] L(λ)=(log||BV λ -ΔΦ||2,log||V λ ||2)

[0056] Among them, V λ is the phase change rate vector for a given λ value, L(λ) is the L curve, and λ is the regularization parameter;

[0057] Step 432: Use the L curve to optimize λ and obtain the optimized regularization parameter λ opt , specifically:

[0058] First, obtain the curvature of the point on the L curve, specifically:

[0059]

[0060] ρ=log||BV-ΔΦ||2

[0061] η=log||V||2

[0062] Where ρ is the residual norm, η is the solution norm, ρ' is the first derivative of ρ, ρ" is the second derivative of ρ, η' is the first derivative of η, η" is the second derivative of η, and δ1 is the curvature of the point on the L curve;

[0063] Then, the λ value corresponding to the maximum point of curvature is used as the optimized regularization parameter λ opt :

[0064]

[0065] Where δ(λ) is the maximum value of the point curvature on the L curve;

[0066] Step 433: Set the regularization parameter λ = λ opt , and use the regularization parameter to construct the objective function:

[0067]

[0068] Here, ||·||2 is the two-norm.

[0069] Furthermore, the objective function in step 44 is solved to obtain the optimal phase change rate, specifically:

[0070] First, when the regularization matrix is ​​the identity matrix I, the phase change rate estimate is obtained:

[0071]

[0072] in, is the estimated rate of phase change;

[0073] Then, perform singular value decomposition on the matrix B to obtain the eigenvalue matrix S as follows:

[0074]

[0075] Among them, U is an M×N matrix, and each column of U is BB T The characteristic vector of , W is an N×M matrix, each column of W is B T The eigenvector of B, matrix U and matrix W are both orthogonal matrices, u i are the diagonal elements of the matrix U, w i is the diagonal element of the matrix W, i is the index of the diagonal element, S is a diagonal matrix, and the value on the diagonal of S is B T B corresponds to the square root of the eigenvalue;

[0076] Finally, the optimal phase change rate vector is obtained using the characteristic matrix S and the phase change rate estimate, specifically:

[0077]

[0078] in, is the optimal phase change rate vector.

[0079] The beneficial effects of the present invention are:

[0080] The present invention proposes a method for inverting surface deformation information. The present invention combines ridge regression regularization with truncated singular value decomposition (TSVD) and then optimizes parameters in combination with the L curve. The present invention eliminates unstable parts in the solution by truncating small singular values, and then adds a regularization term to stabilize the solution. The present invention uses the L curve for parameter selection to balance the weights of data fitting and regularization terms. The present invention can significantly improve the robustness of deformation inversion. When the traditional short baseline set InSAR method has an insufficient number of interferometric pairs that meet the requirements, severe noise interference, or overlapping time baselines, the present invention can ensure the stability of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 Flowchart of the present invention;

[0082] Figure 2 Optimize the parameter curve graph for L curve;

[0083] Figure 3 Comparison chart of estimation results of four methods under highly ill-conditioned problems;

[0084] Figure 4 This is the study area map for deformation inversion;

[0085] Figure 5(a) shows the average deformation rate before the earthquake;

[0086] Figure 5(b) shows the average deformation rate during an earthquake;

[0087] Figure 5(c) shows the average deformation rate diagram after the earthquake;

[0088] Figure 6 is the overall average deformation rate diagram;

[0089] Figure 7 This is the curve of the accumulated deformation of the surface after the earthquake;

[0090] Figure 8(a) shows the first observed surface cumulative deformation map after the earthquake;

[0091] Figure 8(b) shows the second observed cumulative surface deformation map after the earthquake;

[0092] Figure 8(c) shows the third observed cumulative surface deformation map after the earthquake;

[0093] Figure 8(d) shows the fourth observed surface cumulative deformation map after the earthquake;

[0094] Figure 8(e) shows the fifth observed cumulative surface deformation map after the earthquake;

[0095] Figure 8(f) shows the sixth observed cumulative surface deformation map after the earthquake;

[0096] Figure 8(g) shows the seventh observed cumulative surface deformation map after the earthquake;

[0097] Figure 8(h) shows the eighth observed cumulative surface deformation after the earthquake;

[0098] Figure 8(i) shows the surface cumulative deformation diagram for the ninth observation after the earthquake;

[0099] Figure 8(j) shows the cumulative surface deformation at the tenth observation after the earthquake;

[0100] Figure 8(k) shows the cumulative surface deformation map for the eleventh observation after the earthquake;

[0101] Figure 8(l) shows the cumulative surface deformation observed at the twelfth time after the earthquake;

[0102] Figure 8(m) shows the cumulative surface deformation observed at the twelfth time after the earthquake;

[0103] Figure 8(n) shows the cumulative surface deformation observed at the twelfth time after the earthquake;

[0104] Figure 8(o) shows the cumulative surface deformation observed at the twelfth time after the earthquake. DETAILED DESCRIPTION

[0105] Specific implementation method 1: Figure 1 As shown in FIG, the specific process of the SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement in this embodiment is as follows:

[0106] Step 1: Acquire SAR images and preprocess them to obtain preprocessed SAR images. Specifically:

[0107] Step 1: Acquire SAR images, convert the SAR images into single-view complex format, and then use precise orbit data to eliminate phase deviation caused by orbit error;

[0108] Step 1 and 2: Crop the SAR image according to the study area to obtain the preprocessed SAR image.

[0109] Step 2: Based on the short baseline diversity principle, the thresholds of the time baseline and the space baseline are set. The interferogram is obtained using the preprocessed SAR image and the thresholds of the time baseline and the space baseline, and the original interferogram is generated. Specifically:

[0110] Step 21: At imaging time t0~t N Acquire N+1 SAR images for the same observation area, select one of the SAR images as the super master image, register the other SAR images with the super master image, and screen and separate the registered SAR images to obtain M interferometric image pairs.

[0111] The screening and classification of the registered SAR images are specifically as follows:

[0112] First, all pre-processed SAR images are paired into interferometric image pairs;

[0113] Paired pairing means any SAR image forms an interferometric pair with all other SAR images;

[0114] Then, thresholds of the time baseline and the space baseline are set. Interference image pairs that are smaller than the thresholds of the time baseline and the space baseline are considered to meet the requirements (diversity results). The number of interference image pairs that meet the requirements is M, and M meets the following conditions:

[0115]

[0116] Among them, N+1 is the total number of SAR images in the current observation area;

[0117] Step 22: Calculate the phase Φ of each SAR image in the interferometric image pair and the differential interferometric phase ΔΦ of each interferometric image pair, thereby generating an original interferogram.

[0118] Step 3: Remove the interference phase in the original interferogram and perform phase unwrapping to obtain the true differential interferometry phase. Specifically:

[0119] Firstly, the differential interferometric phase of the original interferogram is obtained by using the dichotomy method, and the external DEM data is introduced to remove the interference of the terrain phase.

[0120] Then, the flat-ground effect removal method based on frequency shift estimation is used to remove the flat-ground phase in the original interferogram.

[0121] Finally, the phase unwrapping method based on the minimum cost flow method of the network model is used to unwrap the phase of the original interferogram, and finally the true differential interferometry phase is obtained;

[0122] After removing the above flat ground effect, removing terrain phase interference and phase unwrapping, the true differential interferometry phase is obtained.

[0123] Step 4: Construct a design matrix. Use the design matrix and the true differential interferometry phase to construct the relationship between the phase change rate vector and the observation vector. Use the improved hybrid regularization method and the relationship between the phase change rate vector and the observation vector to construct the objective function. Solve the objective function to obtain the optimal phase change rate vector, thereby completing the surface deformation inversion. Specifically:

[0124] Step 4.1: Construct the design matrix A using the true differential interferometry phase:

[0125] First, one SAR image in each interferometric image pair obtained in step 2 is used as the primary image, and the other SAR image is used as the auxiliary image;

[0126] Then, the design matrix A with M rows and N+1 columns is constructed using the main image and the auxiliary image. The positions of "1" and "-1" in the row vector of the design matrix A are used to represent the selection of interference pairs in an original interferogram. The column vector is used to represent the SAR image at a certain moment. In the matrix A, 1 represents the main image, -1 represents the auxiliary image, and other positions are set to 0.

[0127] For example:

[0128]

[0129] Among them, the -1 and 1 positions in each row correspond to an image pair;

[0130] The design matrix A in this step is an approximate correlation matrix, in which the specific values ​​are mainly based on the actual grouping of the interference image pairs, and is used to mark the positions and time intervals of the main and auxiliary images.

[0131] Step 42: Use the design matrix and the true differential interferometer phase to construct the relationship between the phase change rate vector and the observation vector, specifically:

[0132] Step 421: Use the true differential interferometry phase to construct the phase vector and observation vector, specifically:

[0133] A certain coordinate is between t0 and t N The vector formed by the phase in the time series is an unknown quantity and can be expressed as:

[0134] Φ T =[φ(t0),φ(t1),…,φ(t N )]

[0135] Among them, φ(t N ) is the time t N The phase value on T is the phase vector;

[0136] In this step, if the SAR images at a certain moment are completely deleted (not included in any diversity result) after the registered SAR images are screened in step 21, the phase value corresponding to this moment is assigned a random value, which does not affect the calculation result when calculating the matrix A.

[0137] Then the vector formed by the true differential interference phase of the same coordinate in the time series is the observation vector, which is expressed as:

[0138] ΔΦ T =[δφ1,δφ2,…,δφ M ]

[0139] Where ΔΦ T is the observation vector, δφ M is the true differential interferometric phase of the Mth interferometric pair;

[0140] Step 422: Establish the relationship between the phase vector and the observation vector, specifically:

[0141] First, the main images in the M image pairs are sorted in chronological order. The main image label sequence is:

[0142] IM=[IM1,IM2,…,IM M ]

[0143] Among them, IM j is the original label of the jth main image after time sorting, j∈[1,M], IM is the main image label sequence;

[0144] The original label refers to the label of the image before it is sorted by time;

[0145] Then, the auxiliary images in the M image pairs are sorted in chronological order, and the auxiliary image label sequence is:

[0146] IS=[IS1,IS2,…,IS M ]

[0147] Among them, IS j is the original label of the jth auxiliary image after sorting by time, j∈[1,M], IS is the auxiliary image label sequence;

[0148] Then, according to the auxiliary image at time satisfying the relationship IM j >IS j (j=1,2,…,M), the observation equation is written as:

[0149]

[0150] in, Is the main image IM j The phase, It is auxiliary image IM j The phase, Is the main image IM j The corresponding moment, Auxiliary image IS j The corresponding moment, ΔΦ j is the observation equation of the jth interferometric image pair;

[0151] Finally, since the observation equation contains M equations and N+1 unknowns, the design matrix A is used to establish the relationship between the phase vector and the observation vector:

[0152] ΔΦ=AΦ

[0153] Step 423: Use the relationship between the phase vector and the observation vector to construct the relationship between the phase change rate vector and the observation vector, specifically:

[0154] First, in order to make the obtained solution more meaningful, the problem of solving the deformation phase can be transformed into the problem of solving the deformation rate. The phase change rate vector is:

[0155]

[0156] Among them, V T is the phase change rate vector, v N is the Nth phase change rate, l is the phase change rate index, v l is the lth phase change rate, l∈[1,N], φ(t l ) is the time t l The phase value on φ(t l-1 ) is the time t l-1 Phase value on ;

[0157] Then, the differential interference phase of the j-th interference pair is obtained, specifically:

[0158]

[0159] Among them, k is the time label, v k is the kth phase change rate, t k is the kth moment;

[0160] Then, the differential interferometry phase of the j-th interferometer pair is the integral of the phase change rate of each period over the time interval of the primary and auxiliary images. The relationship between the phase change rate vector and the observation vector is constructed based on the differential interferometry phase of the j-th interferometer pair:

[0161] BV=ΔΦ

[0162] B(j,k)=t k+1 -t k (IS j +1<k<IM j )

[0163] Where B is an M×N matrix, the value of the jth row and kth column is B(j,k), and the values ​​of the remaining positions are 0, t k+1 is the k+1th moment, t k is the kth moment;

[0164] Step 4.3: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, and use the L curve to obtain the optimized regularization parameter λ opt , the objective function is constructed using the regularization parameter, specifically:

[0165] For the equation BV = ΔΦ, (the same applies to AΦ = ΔΦ), if the design matrix B is full rank, the least squares method can be used directly to solve it so that the norm of the error vector reaches the minimum value. The least squares objective function is specifically:

[0166]

[0167] Then its least squares solution can be expressed as:

[0168]

[0169] in, is the estimated deformation rate;

[0170] Since the regularization parameter λ is an empirical value, setting it too large will lead to underfitting and excessive smoothing of the results, while setting it too small will lead to overfitting, amplifying noise and affecting the inversion results. Therefore, the present invention uses the L-curve method to optimize the regularization parameter λ.

[0171] Step 431: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, specifically:

[0172] L(λ)=(log||BV λ -ΔΦ||2,log||V λ ||2)

[0173] Among them, V λ is the phase change rate vector for a given λ value, L(λ) is the L curve, and λ is the regularization parameter;

[0174] Step 432: Figure 2 As shown, the L curve is used to optimize λ to obtain the optimized regularization parameter λ opt , specifically:

[0175] First, obtain the curvature of the point on the L curve, specifically:

[0176]

[0177] ρ=log||BV-ΔΦ||2

[0178] η=log||V||2

[0179] Where ρ is the residual norm, η is the solution norm, ρ' is the first derivative of ρ, ρ" is the second derivative of ρ, η' is the first derivative of η, η" is the second derivative of η, and δ1 is the curvature of the point on the L curve;

[0180] Then, the λ value corresponding to the maximum point of curvature is used as the optimized regularization parameter λ opt :

[0181]

[0182] Where δ(λ) is the maximum value of the point curvature on the L curve;

[0183] Step 433: Set the regularization parameter λ = λ opt , after introducing the ridge regression regularization method, the objective function is constructed using the regularization parameters:

[0184]

[0185] Among them, ||·||2 is the two-norm;

[0186] Step 4. Solve the objective function and obtain the optimal phase change rate to complete the surface deformation inversion. Specifically:

[0187] When the regularization matrix is ​​the identity matrix I, the solution can be obtained:

[0188]

[0189] in, is the estimated rate of phase change;

[0190] When λ = 0, the solution obtained is the traditional least squares solution;

[0191] For the TSVD method, the singular value decomposition is performed on the matrix B to obtain the eigenvalue matrix S, as follows:

[0192]

[0193] Among them, U is an M×N matrix, each column is BB T The characteristic vector of , W is an N×M matrix, each column is B T The eigenvector of B, matrix U and matrix W are both orthogonal matrices, σ i are the diagonal elements of the matrix S, u i are the diagonal elements of the matrix U, w i are the diagonal elements of the matrix W, i is the index of the diagonal element, S is a diagonal matrix and its diagonal values ​​are B T B corresponds to the square root of the eigenvalue. Since the rank of matrix B is N-L+1, the matrix S has only N-L+1 non-zero diagonal values, and L is a constant.

[0194] Set the eigenvalue threshold τ=εσ1,ε∈[10 -3 ,10 -1 ], σ1 is the largest eigenvalue, and the phase change rate can be expressed as:

[0195]

[0196] in, is the phase change rate vector obtained after singular value decomposition, ε is a constant, and k' is an intermediate variable;

[0197] Apply the ridge regression method on the low-dimensional space after TSVD truncation, that is, perform ridge regression on the retained singular values. The solution after combining ridge regression regularization and TSVD method can be expressed as:

[0198]

[0199] in, It is the phase change rate vector (optimal phase change rate vector) obtained after combining ridge regression regularization with the TSVD method, which characterizes the surface deformation rate and then obtains the surface accumulation variable, thereby completing the inversion of landmark deformation information.

[0200] Example: In order to verify the beneficial effects of the present invention, the present invention conducted the following experiments:

[0201] Example 1:

[0202] This embodiment simulates a multi-scale real deformation scene, and simulates highly pathological data by constructing a rank-deficient matrix B. At the same time, high noise interference is introduced into the observation data to verify the anti-interference effect of the present invention.

[0203] The length of the time series is set to 40, the number of interference patterns is set to 60, and the first 30 rows of the matrix B are composed of linearly correlated sequences, making the matrix rank-deficient, with a rank of r = 30.

[0204] The generated deformation signal consists of three parts: linear surface deformation, periodic seasonal deformation and high-frequency noise. For this scenario and the ill-conditioned matrix B, the least squares method, ridge regression method, TSVD method and hybrid regularization method are used to invert the deformation, and the final result is as follows: Figure 3 Table 1 shows the performance comparison of the four methods for deformation parameter inversion. The hybrid regularization method can effectively and significantly improve the robustness of deformation inversion.

[0205] Table 1

[0206]

[0207] Example 2:

[0208] This example selects 15 scenes of Sentinel-1 radar data with relatively uniform time series as the radar data source for time-series InSAR analysis to carry out surface deformation monitoring and regularity research. Radar SLC images are obtained by screening at an average time interval of 12 days. The acquired SAR images cover the period from August 26, 2024 to February 22, 2025. The imaging mode is IW mode, vertically polarized SAR images are selected, and the orbit direction is ascending. Figure 4 It is the study area for deformation inversion.

[0209] In the analysis of measured deformation data based on Sentinel-1 satellites, earthquake events have become a key node affecting the surface deformation pattern. At 9:05 on January 7, 2025, a magnitude 6.8 earthquake occurred in Dingri County, Shigatse City, Tibet Autonomous Region (28.50 degrees north latitude, 87.45 degrees east longitude), with a focal depth of 10 kilometers. In order to study the impact of earthquakes on surface deformation, this embodiment selects the study area as the vicinity of the epicenter of Dingri County, which is connected to the epicenter area through the Pengqu River. The terrain there is mainly mountainous, and the flat areas are part of cities and towns. After the 15 SAR images are freely combined, 105 pairs of interference images can be generated. The time and space baseline parameters of these SAR images are obtained, the time baseline threshold is set to 40 days, and the space baseline threshold is set to 200 meters. After screening, there are 36 pairs of interference images that meet the conditions. After differential interference processing, the deformation phase is extracted for subsequent inversion. This embodiment introduces the Copernicus 30-meter precision digital elevation model (Copernicus DEM) and adopts the dichotomy method to obtain the differential interferometry phase. After removing the flat ground phase and phase unwrapping, the real deformation phase is obtained, which is only related to the surface deformation. The design matrix is ​​constructed according to the time-space baseline combination method in units of days, and the surface deformation rate and deformation amount are inverted by the hybrid regularization method combined with the L-curve parameter optimization. The two SAR images on January 5 and January 17, 2025, have sudden changes in deformation amount and deformation rate due to the earthquake, so the average deformation rate before and after the earthquake is drawn respectively. Figure 5(a)-Figure 5(c) The overall deformation rate in the time series is shown as Figure 6 shown.

[0210] The surface deformation before the earthquake manifested as linear deformation caused by seasonal settlement or tectonic activity; it showed a slow accumulation feature related to the meteorological line, with a small deformation rate, and an average deformation rate of -0.8mm / day to 0.8mm / day; during the earthquake, the ground showed violent uplift and settlement due to crustal movement, and the average deformation rate during the observation period could reach up to 20mm / day; after the main shock, due to the influence of multiple aftershocks, secondary displacement occurred, and the average surface rate in the observation area was -0.2mm / day to 1.4mm / day, which was higher than the deformation rate before the earthquake.

[0211] The cumulative deformation of the high coherence point time series after inversion is as follows Figure 7 As shown in Figure 2. The surface deformation before the earthquake was seasonal, and the overall deformation showed obvious topographic characteristics, with the cumulative deformation within 50 mm. The inversion results within the time interval of the earthquake showed obvious deformation changes, with the cumulative deformation reaching up to 276 mm. The cumulative surface deformation after the earthquake is shown in Figure 2. Figure 8(a)-Figure 8(o) As shown, the overall trend is slightly subsiding. The inversion data of the cumulative deformation show that the west side of the earthquake rupture zone is subsiding and the east side is uplifting, which is consistent with the actual situation.

Claims

1. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization is characterized by The specific process of the method is: Step 1: Acquire SAR images and preprocess the SAR images to obtain preprocessed SAR images; Step 2: Set the thresholds of the time baseline and the space baseline, obtain the interferometric pairs using the pre-processed SAR image and the thresholds of the time baseline and the space baseline, and generate the original interferogram; Step 3: Remove the interference phase in the original interferogram and perform phase unwrapping to obtain the true differential interferometry phase; Step 4: Construct a design matrix, use the design matrix and the true differential interferometer phase to construct the relationship between the phase change rate vector and the observation vector, and use the hybrid regularization method and the relationship between the phase change rate vector and the observation vector to construct the objective function. Solve the objective function to obtain the optimal phase change rate vector, thereby completing the surface deformation inversion.

2. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 1, characterized in that: The step 1 of acquiring a SAR image and preprocessing the SAR image to obtain a preprocessed SAR image is as follows: Step 1: Acquire SAR images, convert the SAR images into single-look complex format, and then eliminate the phase deviation caused by orbit error; Step 1 and 2: Crop the SAR image according to the study area to obtain the preprocessed SAR image.

3. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 2, characterized in that: In step 2, the thresholds of the time baseline and the space baseline are set, and the interferogram is obtained using the pre-processed SAR image and the thresholds of the time baseline and the space baseline, and the original interferogram is generated, specifically: Step 21: At imaging time t0~t N Acquire N+1 SAR images for the same observation area, select one of the SAR images as the super master image, register the other SAR images with the super master image, and screen and separate the registered SAR images to obtain M interferometric image pairs. Step 22: Calculate the phase Φ of each SAR image in the interferometric image pair and the differential interferometric phase ΔΦ of each interferometric image pair, thereby generating an original interferogram.

4. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 3, characterized in that: The step 21 of screening and dividing the registered SAR images is as follows: First, all pre-processed SAR images are paired into interferometric image pairs; Then, the thresholds of the time baseline and the space baseline are set, and the interference image pairs that are smaller than the thresholds of the time baseline and the space baseline are considered to meet the requirements; The number of interference image pairs that meet the requirements is M, and M meets the following conditions: Among them, N+1 is the total number of SAR images in the current observation area.

5. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 4, characterized in that: The step 3 removes the interference phase in the original interferogram and performs phase unwrapping processing to obtain the true differential interferogram phase, specifically: Firstly, the differential interferometric phase of the original interferogram is obtained by using the dichotomy method, and then the external DEM data is introduced to remove the interference of the terrain phase. Then, the flat-ground effect removal method based on frequency shift estimation is used to remove the flat-ground phase in the original interferogram. Finally, the phase unwrapping method based on the minimum cost flow method of the network model is used to unwrap the phase of the original interferogram, and finally the true differential interferometry phase is obtained.

6. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 5, characterized in that: In step 4, the design matrix is ​​constructed, and the relationship between the phase change rate vector and the observation vector is constructed using the design matrix and the true differential interferometry phase. The objective function is constructed using the hybrid regularization method and the relationship between the phase change rate vector and the observation vector. The objective function is solved to obtain the deformation rate, thereby completing the surface deformation inversion, specifically: Step 41: Construct the design matrix A using the real differential interferometer phase; Step 42: Use the design matrix and the true differential interferometer phase to construct the relationship between the phase change rate vector and the observation vector; Step 4.3: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, and use the L curve to obtain the optimized regularization parameter λ opt , construct the objective function using the regularization parameter; Step 44: Solve the objective function and obtain the optimal phase change rate to complete the surface deformation inversion.

7. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 6, characterized in that: The design matrix A is constructed using the real differential interferometry phase in step 41, specifically: First, one SAR image in each interferometric image pair obtained in step 2 is used as the primary image, and the other SAR image is used as the auxiliary image; Then, the design matrix A with M rows and N+1 columns is constructed using the main image and the auxiliary image; The positions of "1" and "-1" in the row vector of the design matrix A represent the positions of interference pairs in an original interference pattern, and the column vector represents the SAR image at a certain moment; In the matrix A, 1 represents the primary image, -1 represents the secondary image, and other positions are set to 0.

8. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 7, characterized in that: The relationship between the phase change rate vector and the observation vector is constructed by using the design matrix and the true differential interferometry phase in step 42, specifically: Step 421: Use the true differential interferometry phase to construct the phase vector and observation vector, specifically: F T =[φ(t0),φ(t1),…,φ(t N )] DF T =[δφ1,δφ2,…,δφ M ] Among them, φ(t N ) is the time t N The phase value on T is the phase vector, ΔΦ T is the observation vector, δφ M is the true differential interferometric phase of the Mth interferometric pair; Step 422: Establish the relationship between the phase vector and the observation vector, specifically: ΔΦ=AΦ IS=[IS1,IS2,…,IS M ] IM=[IM1,IM2,…,IM M ] Where ΔΦ j is the jth element in ΔΦ, Is the main image IM j The corresponding moment, Auxiliary image IS j The corresponding moment, Is the main image IM j The phase, Auxiliary Image IM j phase, IS is the auxiliary image label sequence, IM is the main image label sequence, j∈[1,M], IM j IS is the original label of the jth main image after sorting by time, j is the original label of the jth auxiliary image after sorting by time, IM j >IS j ; Step 423: Use the relationship between the phase vector and the observation vector to construct the relationship between the phase change rate vector and the observation vector, specifically: Among them, δφ j is the jth element in ΔΦ, k is [IS j +1,IM j ], B is an M×N matrix, B(j,k) is the value of the jth row and kth column in B, t k+1 is the k+1th moment, t k is the kth moment, k is the time label, v k is the kth phase change rate, l∈[1,N], V T is the phase change rate vector, v N is the Nth phase change rate, l is the phase change rate index, v l is the lth phase change rate, φ(t l ) is the time t l The phase value on φ(t l-1 ) is the time t l-1 The phase value on .

9. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 8, characterized in that: In step 43, the L curve is constructed according to the relationship between the phase change rate vector and the observation vector, and the L curve is used to obtain the optimized regularization parameter λ opt , the objective function is constructed using the regularization parameter, specifically: Step 431: Construct an L curve based on the relationship between the phase change rate vector and the observation vector, specifically: L(λ)=(log||BV λ -ΔΦ||2,log||V λ ||2) Among them, V λ is the phase change rate vector for a given λ value, L(λ) is the L curve, and λ is the regularization parameter; Step 432: Use the L curve to optimize λ and obtain the optimized regularization parameter λ opt , specifically: First, obtain the curvature of the point on the L curve, specifically: ρ=log||BV-ΔΦ||2 η=log||V||2 Where ρ is the residual norm, η is the solution norm, ρ' is the first derivative of ρ, ρ" is the second derivative of ρ, η' is the first derivative of η, η" is the second derivative of η, and δ1 is the curvature of the point on the L curve; Then, the λ value corresponding to the maximum point of curvature is used as the optimized regularization parameter λ opt : Where δ(λ) is the maximum value of the point curvature on the L curve; Step 433: Set the regularization parameter λ = λ opt , and use the regularization parameter to construct the objective function: Here, ||·||2 is the two-norm.

10. The improved SBAS-InSAR surface deformation inversion method based on hybrid regularization according to claim 9, characterized in that: The objective function in step 44 is solved to obtain the optimal phase change rate, specifically: First, when the regularization matrix is ​​the identity matrix I, the phase change rate estimate is obtained: in, is the estimated rate of phase change; Then, perform singular value decomposition on the matrix B to obtain the eigenvalue matrix S as follows: Among them, U is an M×N matrix, and each column of U is BB T The characteristic vector of , W is an N×M matrix, each column of W is B T The eigenvector of B, matrix U and matrix W are both orthogonal matrices, u i are the diagonal elements of the matrix U, w i is the diagonal element of the matrix W, i is the index of the diagonal element, S is a diagonal matrix, and the value on the diagonal of S is B T B corresponds to the square root of the eigenvalue; Finally, the optimal phase change rate vector is obtained using the characteristic matrix S and the phase change rate estimate, specifically: in, is the optimal phase change rate vector.

Images