Surface deformation inversion method based on hybrid regularization improved SBAS-InSAR

Abstract

The method for improving SBAS-InSAR surface deformation inversion based on mixed regularization relates to the field of remote sensing image processing of synthetic aperture radar.The present application is to solve the problem of unstable inversion results of the existing terrain inversion method.The present application comprises: obtaining SAR images, and pre-processing the SAR images to obtain pre-processed SAR images; setting the threshold of time baseline and space baseline, obtaining the interference pair by using the pre-processed SAR images and the threshold of time baseline and space baseline, and generating the original interference graph; removing the interference phase in the original interference graph and performing phase unwrapping processing to obtain the real differential interference phase; constructing a design matrix, constructing the relationship between the phase change rate vector and the observation vector by using the design matrix and the real differential interference phase, and constructing the objective function by using the mixed regularization method and the relationship between the phase change rate vector and the observation vector, and solving the objective function to obtain the deformation rate, so as to complete the surface deformation inversion.The present application is used for surface deformation inversion.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing of synthetic aperture radar, and particularly to a SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement. Background Technology

[0002] Surface deformation is an environmental geological phenomenon caused by the compression of rock masses on the Earth's surface, resulting in changes in regional surface elevation information. It can cause permanent damage to environmental resources and human production and life. InSAR technology has been widely used in surface deformation monitoring. The acquisition of three-dimensional ground information can effectively carry out geological disaster monitoring and prevention. In addition, it also plays an important role in the prediction of urban land subsidence and mining area collapse.

[0003] Proposed by Berardion et al. in 2002, SBAS-InSAR technology improves the overall coherence of interferometric pairs by shortening the temporal and spatial baselines, thus obtaining more accurate target point results. It plays a crucial role in acquiring ground deformation information. The proposal and development of SBAS-InSAR technology have further promoted the widespread application of remote sensing in the field of land surface monitoring, and in recent years it has become a focus of scientists' research on land surface deformation. Berardion et al. used singular value decomposition and SBAS-InSAR technology to obtain ground deformation results for a city and a volcano from 1992 to 2000 using 44 SAR images. The results were consistent with GPS monitoring results, proving the effectiveness of the SBAS-InSAR algorithm. In 2004, Riccardo Lanari improved Berardion's method, considering the influence of different phases on the monitoring results. He estimated the low-pass deformation phase in the multi-look processed data, introduced DEM error phase and atmospheric delay phase, and used the improved method to monitor land deformation in the same area. The results were consistent with the leveling measurements. In 2008, Casu et al. used SBAS-InSAR technology to monitor surface deformation in a certain region, obtaining the land subsidence evolution process and average subsidence rate from 1992 to 2000. In 2011, Pepe et al., for the first time, used SBAS-InSAR technology and ENVISAT ASAR stripmap and scanSAR images to monitor the deformation of a volcanic cluster. This was the first time that images from different working modes, such as Stripmap and ScanSAR, were jointly processed, improving the temporal resolution of SBAS technology monitoring and laying the foundation for the joint processing of multi-angle and multi-working mode data. In 2014, Calò et al. proposed the P-SBAS-InSAR technology based on a two-layer parallel method and used SAR image data in the X and C bands to obtain landslide subsidence monitoring results in the Ivancich region of Italy from 1992 to 2010 using the short baseline set InSAR method. In 2023, Ma Chuang et al. used SBAS-InSAR technology, based on 23 C-band Sentinel-1A ascending-orbit SAR images, to monitor a landslide in Baige Village, Boluo Township, Jiangda County, Tibet Autonomous Region, and predicted the time of the disaster using a velocity reciprocal model. The aforementioned traditional short baseline set SBAS-InSAR method can effectively utilize redundant interferometric pairs for interferometric phase acquisition, overcoming the shortcomings of conventional D-InSAR technology such as severe decorrelation and low deformation detection accuracy, and is highly adaptable to atmospheric, topographic, and surface changes. However, when faced with insufficient numbers of qualified interferometric pairs, severe noise interference, or overlapping time baselines, these methods may suffer from least-squares instability due to rank deficiency in the design matrix, leading to unstable surface deformation inversion results. Summary of the Invention

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

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

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

[0007] Step 2: Set the thresholds for the temporal and spatial baselines, obtain the interferometric pairs using the preprocessed SAR image and the thresholds for the temporal and spatial baselines, and generate the original interferogram;

[0008] Step 3: Remove the interfering phase from the original interferogram and perform phase unwrapping to obtain the true differential interferometric phase;

[0009] Step 4: Construct the design matrix. Use the design matrix and the true differential interferometric phase to construct the relationship between the phase change rate vector and the observation vector. Then, 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 acquisition of SAR images in step one, and the preprocessing of the SAR images to obtain preprocessed SAR images, specifically involves:

[0011] Step 1: Acquire SAR images, convert SAR images to single-look complex format, and then eliminate phase deviations caused by orbital errors;

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

[0013] Furthermore, in step two, setting thresholds for the time and spatial baselines, using the preprocessed SAR image and the thresholds for the time and spatial baselines to obtain interferometric pairs, and generating the original interferogram, specifically involves:

[0014] Step 2.1, during the imaging time t0~t N N+1 SAR images are acquired for the same observation area. One of the SAR images is selected as the super master image. The other SAR images are registered with the super master image. The registered SAR images are then filtered and divided to obtain M interferometric pairs.

[0015] Step 2: Calculate the phase Φ of each SAR image in the interferometric image pair and the differential interferometric phase ΔΦ of each interferometric image pair to generate the original interferogram.

[0016] Furthermore, the screening and diversity analysis of the registered SAR images in step two-one specifically involves:

[0017] First, all the preprocessed SAR images are paired up to form interferometric pairs;

[0018] Then, set thresholds for the time baseline and the spatial baseline. Interference pairs that are less than the thresholds for the time baseline and the spatial baseline are considered to be qualified interference pairs.

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

[0020]

[0021] N+1 represents the total number of SAR images in the current observation area.

[0022] Furthermore, in step three, removing the interference phase from the original interferogram and performing phase unwrapping to obtain the true differential interferometric phase specifically involves:

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

[0024] Then, a flat-ground phase in the original interferogram is removed using a frequency shift estimation-based method for removing flat-ground effects.

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

[0026] Further, in step four, the design matrix is ​​constructed, and the relationship between the phase change rate vector and the observation vector is established using the design matrix and the true differential interferometric phase. Then, a target function is constructed using a hybrid regularization method and the relationship between the phase change rate vector and the observation vector. Solving the target function yields the deformation rate, thus completing the surface deformation inversion. Specifically:

[0027] Step 41: Construct the design matrix A using the actual differential interferometric phase;

[0028] Step 4.2: Construct the relationship between the phase change rate vector and the observation vector using the design matrix and the true differential interferometric phase;

[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 λ. optThe objective function is constructed using regularization parameters;

[0030] Step 4: Solve for the objective function to obtain the optimal phase change rate, thereby completing the surface deformation inversion.

[0031] Furthermore, the construction of the design matrix A using the true differential interferometric phase in step four-one specifically involves:

[0032] First, take one SAR image from each interferometric image pair obtained in step two as the main image and the other SAR image as the auxiliary image.

[0033] Then, a 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 vectors of the design matrix A represent the positions of the interferometric pairs in an original interferogram, and the column vectors represent the SAR image at a certain moment.

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

[0036] Furthermore, the step four-two, which involves constructing the relationship between the phase change rate vector and the observation vector using the design matrix and the true differential interferometric phase, specifically involves:

[0037] Step 421: Construct the phase vector and observation vector using the true differential interferometric phase, specifically as follows:

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

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

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

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

[0042] ΔΦ=AΦ

[0043]

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

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

[0046] Where, ΔΦ j It is the j-th element in ΔΦ. It is the main image IM j At the corresponding time, It is an auxiliary image IS j At the corresponding time, It is the main image IM j phase, It is an auxiliary image IM j The phase, IS is the secondary image label sequence, IM is the primary image label sequence, j∈[1,M], IM j It is the original label of the j-th main image after sorting by time, IS j It is the original label of the j-th auxiliary image after sorting by time, IM j >IS j ;

[0047] Step 423: Construct the relationship between the phase change rate vector and the observation vector using the relationship between the phase vector and the observation vector, specifically as follows:

[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] Where, δφ j It is the j-th element in ΔΦ, and k is [IS j +1,IM j The values ​​in ] are given, where B is an M×N matrix, B(j,k) is the value in the j-th row and k-th column of B, and t k+1 It is the (k+1)th time, t k It is the k-th moment, where k is the time index, v k It is the k-th phase change rate, l∈[1,N], V T It is the phase change rate vector, v N It is the Nth phase change rate, l is the phase change rate label, v l It is the l-th phase change rate, φ(t) l) is time t l The phase value on, φ(t) l-1 ) is time t l-1 The phase value on.

[0053] Furthermore, in step four-three, an L-curve is constructed based on the relationship between the phase change rate vector and the observation vector, and the optimized regularization parameter λ is obtained using the L-curve. opt The objective function is constructed using regularization parameters, specifically as follows:

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

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

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

[0057] Step 432: Optimize λ using the L-curve to obtain the optimized regularization parameter λ. opt Specifically:

[0058] First, obtain the curvature of a 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 a point on the L curve.

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

[0064]

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

[0066] Step 433: Let the regularization parameter λ = λ opt Construct the objective function using regularization parameters:

[0067]

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

[0069] Furthermore, in step four, solving the objective function to obtain the optimal phase change rate specifically involves:

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

[0071]

[0072] in, It is an estimate of the rate of phase change;

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

[0074]

[0075] Where U is an M×N matrix, and each column of U is BB. T The eigenvectors of W are given by the matrix B, where W is an N×M matrix and each column of W is a B-value. T The eigenvectors of B, matrices U and W are both orthogonal matrices, u i w is the diagonal element of matrix U. i Let be the diagonal element of matrix W, i be the label of the diagonal element, S be a diagonal matrix, and the values ​​on the diagonal of S be 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 estimated phase change rate, specifically:

[0077]

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

[0079] The beneficial effects of this invention are as follows:

[0080] This invention proposes a method for inverting surface deformation information. It combines ridge regression regularization with truncated singular value decomposition (TSVD), and further optimizes parameters using L-curves. This invention eliminates unstable parts of the solution by truncating small singular values, and then adds a regularization term to stabilize the solution. By utilizing L-curves for parameter selection, this invention balances the weights of data fitting and the regularization term, significantly improving the robustness of deformation inversion. When traditional short-baseline InSAR methods suffer from insufficient numbers of qualified interferometric pairs, severe noise interference, or overlapping time baselines, this invention ensures stable inversion results. Attached Figure Description

[0081] Figure 1 This is a flowchart of the present invention;

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

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

[0084] Figure 4 This is a map showing the study area 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 the earthquake;

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

[0088] Figure 6 This is a diagram showing the overall average deformation rate.

[0089] Figure 7 This is a graph showing the cumulative deformation of the Earth's surface after the earthquake.

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

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

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

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

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

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

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

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

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

[0099] Figure 8(j) shows the cumulative surface deformation observed ten times after the earthquake;

[0100] Figure 8(k) shows the cumulative surface deformation observed eleven times after the earthquake;

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

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

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

[0104] Figure 8(o) shows the cumulative surface deformation observed for the twelfth time after the earthquake. Detailed Implementation

[0105] Specific implementation method one: as follows Figure 1 As shown, 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 as follows:

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

[0108] Steps 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, set thresholds for the temporal and spatial baselines. Using the preprocessed SAR imagery and the thresholds for the temporal and spatial baselines, obtain interferometric pairs and generate the original interferogram. Specifically:

[0110] Step 2.1, during the imaging time t0~t N N+1 SAR images are acquired for the same observation area. One of the SAR images is selected as the super master image. The other SAR images are registered with the super master image. The registered SAR images are then filtered and divided to obtain M interferometric pairs.

[0111] The process of filtering and grouping the registered SAR images specifically involves:

[0112] First, all the preprocessed SAR images are paired up to form interferometric pairs;

[0113] Pairing refers to any SAR image forming an interferometric pair with all other SAR images;

[0114] Then, thresholds are set for the temporal and spatial baselines. Interferometric pairs smaller than these thresholds are considered acceptable interferometric pairs (diversity results). The number of acceptable interferometric pairs is M, and M satisfies the following conditions:

[0115]

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

[0117] Step 2: Calculate the phase Φ of each SAR image in the interferometric image pair and the differential interferometric phase ΔΦ of each interferometric image pair to generate the original interferogram.

[0118] Step 3: Remove the interfering phase from the original interferogram and perform phase unwrapping to obtain the true differential interferometric phase, specifically:

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

[0120] Then, a flat-ground phase in the original interferogram is removed using a frequency shift estimation-based method for removing flat-ground effects.

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

[0122] After removing the flat terrain effect, eliminating topographic phase interference, and unwrapping the phase, the true differential interference phase is obtained.

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

[0124] Step 41: Construct the design matrix A using the true differential interferometric phase:

[0125] First, take one SAR image from each interferometric image pair obtained in step two as the main image and the other SAR image as the auxiliary image.

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

[0127] For example:

[0128]

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

[0130] The design matrix A in this step is an approximate correlation matrix, in which the specific values ​​are based on the actual interferometric image pair grouping situation, and are used to mark the position and time interval of the primary and secondary images.

[0131] Step 4.2: Construct the relationship between the phase change rate vector and the observation vector using the design matrix and the true differential interferometric phase, specifically as follows:

[0132] Step 421: Construct the phase vector and observation vector using the true differential interferometric phase, specifically as follows:

[0133] A certain coordinate is between t0 and t N The vector formed by the phases in a time series is an unknown quantity, which can be represented as:

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

[0135] Wherein, φ(t) N ) is time t N The phase value on, Φ T It is the phase vector;

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

[0137] The vector formed by the true differential interferometric phases at the same coordinate on this time series is the observation vector, expressed as:

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

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

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

[0141] First, the main images in the M image pairs are sorted in chronological order, and the main image number sequence is as follows:

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

[0143] Among them, IM j It is the original label of the j-th main image after sorting by time, j∈[1,M], and IM is the main image label sequence;

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

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

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

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

[0148] Then, based on the secondary image satisfying the temporal relation IM... j >IS j (j=1,2,…,M), the observation equation can be written as:

[0149]

[0150] in, It is the main image IM j phase, It is an auxiliary image IM j phase, It is the main image IM j At the corresponding time, It is an auxiliary image IS j At the corresponding time, ΔΦ j It is the observation equation for the j-th interferometric image pair;

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

[0152] ΔΦ=AΦ

[0153] Step 423: Construct the relationship between the phase change rate vector and the observation vector using the relationship between the phase vector and the observation vector, specifically as follows:

[0154] First, to make the obtained solution more meaningful, the problem of solving for the deformation phase can be transformed into the problem of solving for the deformation rate, with the phase change rate vector being:

[0155]

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

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

[0158]

[0159] Where k is the time label, v k It is the rate of change of the kth phase, t k It is the k-th moment;

[0160] Then, the differential interferometric phase of the j-th interferometric pair is the integral of the phase change rate over the time interval between the acquisition of the main and auxiliary images. Based on the differential interferometric phase of the j-th interferometric pair, the relationship between the phase change rate vector and the observation vector is constructed:

[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 in the j-th row and k-th column is B(j,k), and the values ​​in the other positions are 0. k+1 It is the (k+1)th time, t k It is the k-th 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 regularization parameters, specifically as follows:

[0165] For the equation BV = ΔΦ (and similarly for AΦ = ΔΦ), if the design matrix B is full rank, it can be solved directly using the least squares method to minimize the norm of the error vector. The least squares objective function is as follows:

[0166]

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

[0168]

[0169] in, It is an estimate of the deformation rate;

[0170] Since the regularization parameter λ is an empirical value, setting it too large can lead to underfitting and excessively smooth results, while setting it too small can lead to overfitting and amplified noise, thus affecting the inversion results. Therefore, this invention uses the L-curve method to optimize the regularization parameter λ.

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

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

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

[0174] Step 432, as follows Figure 2 As shown, the L-curve is used to optimize λ, and the optimized regularization parameter λ is obtained. opt Specifically:

[0175] First, obtain the curvature of a 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 a point on the L curve.

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

[0181]

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

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

[0184]

[0185] Where ||·||2 is the L2 norm;

[0186] Step 44: Solve for the objective function to obtain the optimal phase change rate, thereby completing the surface deformation inversion, specifically as follows:

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

[0188]

[0189] in, It is an estimate of the rate of phase change;

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

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

[0192]

[0193] Where U is an M×N matrix, and each column is BB T The eigenvectors are given by W, which is an N×M matrix where each column is a B vector. T The eigenvectors of B, matrices U and W are both orthogonal matrices, σ i u is a diagonal element of matrix S. i w is the diagonal element of matrix U. i Let be the diagonal element of matrix W, i be the label of the diagonal element, and S be a diagonal matrix whose diagonal elements are B. T B corresponds to the square root of the eigenvalue. Since the rank of matrix B is N-L+1, matrix S has only N-L+1 non-zero diagonal values, where L is a constant.

[0194] Set the eigenvalue threshold τ = εσ1, ε∈

[10] -3 10 -1 ], where σ1 is the largest eigenvalue, the phase change rate can be expressed as:

[0195]

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

[0197] Applying ridge regression to the low-dimensional space truncated by TSVD, i.e., performing ridge regression on the retained singular values, the solution after combining ridge regression regularization and the 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 and TSVD method, which represents the surface deformation rate and then obtains the surface accumulation type variable, thereby completing the inversion of landmark deformation information.

[0200] Example: To verify the beneficial effects of the present invention, the following experiments were conducted:

[0201] Example 1:

[0202] This embodiment simulates a multi-scale real deformation scenario, constructs a rank-deficient matrix B to simulate highly ill-conditioned data, and introduces high noise interference into the observation data to verify the anti-interference effect of the present invention.

[0203] The time series is set to a length of 40, the number of interferograms is 60, and the first 30 rows of matrix B consist 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, least squares, ridge regression, TSVD, and hybrid regularization methods were used for deformation inversion, respectively. The final results are as follows: Figure 3 As shown in Table 1, the performance of the four methods for deformation parameter inversion is compared. The hybrid regularization method can effectively and significantly improve the robustness of deformation inversion.

[0205] Table 1

[0206]

[0207] Example 2:

[0208] This embodiment selects 15 Sentinel-1 radar images with relatively uniform temporal sequence as the radar data source for temporal InSAR analysis to carry out monitoring and pattern research on surface deformation. The radar SLC images are selected according to an average time interval of 12 days. The time span of the acquired SAR images is from August 26, 2024 to February 22, 2025. The imaging mode is IW mode, vertical polarization SAR images are selected, and the orbit direction is the ascending direction. Figure 4 This is the research area for deformation inversion.

[0209] In the analysis of measured deformation data based on the Sentinel-1 satellite, earthquake events have become key nodes affecting surface deformation patterns. At 9:05 AM on January 7, 2025, a magnitude 6.8 earthquake occurred in Dingri County, Shigatse City, Tibet Autonomous Region (28.50°N, 87.45°E), with a focal depth of 10 km. To study the impact of earthquakes on surface deformation, this embodiment selects the area near the epicenter of Dingri County as the study area, connected to the epicenter region via the Pengqu River. The terrain in this area is mainly mountainous, with some flat areas consisting of cities and towns. Fifteen SAR images were freely combined to generate 105 interferometric pairs. The temporal and spatial baseline parameters of these SAR images were acquired, with a temporal baseline threshold of 40 days and a spatial baseline threshold of 200 meters. After screening, 36 interferometric pairs met the criteria. Differential interferometry was performed on these pairs to extract the deformation phase for subsequent inversion. This embodiment introduces the Copernicus 30-meter precision digital elevation model (DEM) and uses the bisection method to obtain the differential interferometric phase. After removing the flat phase and phase unwrapping, the true deformation phase is obtained, which is only related to surface deformation. Using days as the unit, a design matrix is ​​constructed according to the spatiotemporal baseline combination method, and a hybrid regularization method combining L-curve parameter optimization is used to invert the surface deformation rate and deformation. Due to the abrupt changes in deformation and deformation rate caused by the earthquake, the average deformation rates before and after the earthquake are plotted in the two SAR images from January 5th and January 17th, 2025. Figures 5(a)-5(c) As shown. The overall deformation rate within the time series is as follows. Figure 6 As shown.

[0210] Pre-earthquake surface deformation manifests as seasonal subsidence or linear deformation caused by tectonic activity; it exhibits slow accumulation characteristics related to meteorological lines, with a small deformation rate, averaging -0.8 mm / day to 0.8 mm / day; during an earthquake, crustal movement causes severe uplift and subsidence of the ground, with the average deformation rate reaching up to 20 mm / day within the observation period; after the main shock, secondary displacement occurs due to the influence of multiple aftershocks, with the average surface rate in the observation area ranging from -0.2 mm / day to 1.4 mm / day, higher than the pre-earthquake deformation rate.

[0211] The cumulative deformation of the high coherence point in the time series after inversion is as follows: Figure 7 As shown. Pre-earthquake surface deformation is seasonal deformation, with overall deformation exhibiting distinct topographic features, and cumulative deformation within 50 mm. Inversion results within the earthquake's timeframe show significant deformation changes, with cumulative deformation reaching a maximum of 276 mm. Post-earthquake cumulative surface deformation is shown below. Figures 8(a)-8(o) As shown, the overall trend is slightly subsidence. The cumulative deformation inversion data shows that the earthquake rupture zone is characterized by subsidence on the west side and uplift on the east side, which is consistent with the actual situation.

Claims

1. A SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement, characterized in that... The specific process of the method is as follows: Step 1: Acquire SAR images and preprocess them to obtain preprocessed SAR images; Step 2: Set the thresholds for the temporal and spatial baselines, obtain the interferometric pairs using the preprocessed SAR image and the thresholds for the temporal and spatial baselines, and generate the original interferogram; Step 3: Remove the interfering phase from the original interferogram and perform phase unwrapping to obtain the true differential interferometric phase; Step 4: Construct the design matrix. Using the design matrix and the true differential interferometric phase, establish the relationship between the phase change rate vector and the observation vector. Then, construct the objective function using a hybrid regularization method and the relationship between the phase change rate vector and the observation vector. Solve the objective function to obtain the optimal phase change rate vector, thereby completing the surface deformation inversion. Specifically: Step 41: Construct the design matrix A using the actual differential interferometric phase; Step 4.2: Construct the relationship between the phase change rate vector and the observation vector using the design matrix and the true differential interferometric phase, specifically as follows: Step 421: Construct the phase vector and observation vector using the true differential interferometric phase, specifically as follows: in, It is a moment Phase value on, It is the phase vector. It is the observation vector. It is the first The true differential interference phase of each interference pair; Step 422: Establish the relationship between the phase vector and the observation vector, specifically as follows: in, yes The first in One element, It is the main image At the corresponding time, It is an auxiliary image At the corresponding time, It is the main image phase, It is an auxiliary image phase, It is a secondary image label sequence. It is the main image label sequence. , It is the first one after being sorted by time. The original label of the main image. It is the first one after being sorted by time. The original label of each auxiliary image. ; Step 423: Construct the relationship between the phase change rate vector and the observation vector using the relationship between the phase vector and the observation vector, specifically as follows: in, yes The Middle One element, yes The values ​​in the matrix, where B is an M×N matrix. It is the value in row j and column k of B. It is the (k+1)th time. It is the k-th time. It is a time stamp. It is the rate of change of the kth phase. , It is the phase change rate vector. It is the rate of change of the Nth phase. It is the phase change rate label. It is the rate of change of the l-th phase. It is a moment Phase value on, It is a moment Phase value on; 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 parameters. The objective function is constructed using regularization parameters, specifically as follows: Step 431: Construct the L-curve based on the relationship between the phase change rate vector and the observation vector, specifically as follows: in, It is given The phase change rate vector at the value, It is an L-curve. It is a regularization parameter; Step 432: Use the L-curve to... Optimize to obtain optimized regularization parameters Specifically: First, obtain the curvature of a point on the L-curve, specifically: in, It is the residual norm. It is the solution norm. yes The first derivative, yes The second derivative, yes The first derivative, yes The second derivative, It is the curvature of a point on the L-curve; Then, the point corresponding to the maximum curvature The value is used as the optimized regularization parameter. : in, It is the maximum curvature of a point on the L-curve; Step 433: Set the regularization parameter Construct the objective function using regularization parameters: in, It is a norm 2; Step 44: Solve for the objective function to obtain the optimal phase change rate, thereby completing the surface deformation inversion, specifically as follows: First, when the regularization matrix is ​​the identity matrix I, the estimated phase change rate is obtained: in, It is an estimate of the rate of phase change; Then, singular value decomposition is performed on matrix B to obtain the eigenvalue matrix S, as follows: in, yes The matrix, Each column is eigenvectors, yes The matrix, Each column is eigenvectors, matrices sum matrix All are orthogonal matrices. It is a matrix diagonal elements, It is a matrix diagonal elements, These are the labels for the diagonal elements. It is a diagonal matrix. The value on the diagonal The square root of the corresponding eigenvalue, These are the diagonal elements of matrix S; Finally, the optimal phase change rate vector is obtained using the characteristic matrix S and the estimated phase change rate, specifically: in, It is the optimal phase change rate vector.

2. The SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement according to claim 1, characterized in that: The acquisition of SAR images in step one, and the preprocessing of the SAR images to obtain preprocessed SAR images, specifically involves: Step 1: Acquire SAR images, convert SAR images to single-look complex format, and then eliminate phase deviations caused by orbital errors; Steps 1 and 2: Crop the SAR image according to the study area to obtain the preprocessed SAR image.

3. The SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement according to claim 2, characterized in that: In step two, setting thresholds for the time and spatial baselines, using the preprocessed SAR image and the thresholds for the time and spatial baselines, an interferometric pair is obtained, and the original interferogram is generated. Specifically: Step 2.1, during imaging time N+1 SAR images are acquired for the same observation area. One of the SAR images is selected as the super master image. The other SAR images are registered with the super master image. The registered SAR images are then filtered and divided to obtain M interferometric pairs. Step 22: Calculate the phase of each SAR image in the interferometric image pair. The differential interference phase of each interference image pair This generates the original interferogram.

4. The SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement according to claim 3, characterized in that: The step 2.1 of filtering and grouping the registered SAR images specifically involves: First, all the preprocessed SAR images are paired up to form interferometric pairs; Then, set thresholds for the time baseline and the spatial baseline. Interference pairs that are less than the thresholds for the time baseline and the spatial baseline are considered to be qualified interference pairs. The number of qualified interference image pairs is M, and M meets the following conditions: N+1 represents the total number of SAR images in the current observation area.

5. The SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement according to claim 4, characterized in that: Step three, which involves removing the interference phase from the original interferogram and performing phase unwrapping to obtain the true differential interferometric phase, specifically involves: First, the differential interferometric phase of the original interferogram is obtained using the bisection method, and then external DEM data is introduced to remove the interference of the terrain phase. Then, a flat-ground phase in the original interferogram is removed using a frequency shift estimation-based method for removing flat-ground effects. Finally, the phase unwrapping method based on the minimum cost flow of the network model is used to unwrap the phase of the original interferogram, and the true differential interferometric phase is finally obtained.

6. The SBAS-InSAR surface deformation inversion method based on hybrid regularization improvement according to claim 5, characterized in that: The step 4.1, which involves constructing the design matrix A using the true differential interferometric phase, specifically involves: First, take one SAR image from each interferometric image pair obtained in step two as the main image and the other SAR image as the auxiliary image. Then, a 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 vectors of the design matrix A represent the positions of the interferometric pairs in an original interferogram, and the column vectors represent the SAR image at a certain moment. In matrix A, 1 represents the main image, -1 represents the auxiliary image, and other positions are set to 0.

Images