Application of variable projection method based on QR decomposition in surface deformation

Abstract

The present invention discloses an application method of a variable projection method based on QR decomposition in surface deformation. The fitting of surface deformation in a study area in time and space is taken as the research focus. On the basis of applying SBAS-InSAR time series analysis technology to explore the surface deformation in the study area, a complex exponential model with a separable structure is used to fit the settlement situation of the surface in spatial position. The model parameters are solved by a variable projection algorithm based on QR decomposition. The settlement profile curve is approximated with high precision in spatial position. At the same time, the surface is predicted in time in combination with the Logistic model, thereby realizing comprehensive monitoring of the surface deformation in time and space, and providing reliable technical support for surface deformation monitoring.

Description

Technical Field

[0001] The invention belongs to the technical field of satellite remote sensing, and in particular relates to a surface deformation application method of a variable projection method based on QR decomposition. Background Art

[0002] Although the traditional method of monitoring surface deformation has high accuracy, it can only monitor the deformation of a limited number of ground points and cannot monitor the overall deformation of the surface over a large area, so it has certain limitations. Differential Interferometric Synthetic Aperture Radar (D-InSAR) has the advantages of all-day and large-scale monitoring of ground deformation, breaking the limitations of traditional methods (Li Shanshan et al., 2013; Gabriele A. Ketal., 1989; Massonnet D. et al., 1993). However, the number of SAR images used in D-InSAR technology is small, there is a lack of redundant observations, and when the time span is large, it is easy to cause decoherence problems, making it difficult to obtain a continuous deformation field. To overcome this problem, Berardino et al. proposed SBAS-InSAR time series analysis technology based on D-InSAR. Through time series analysis of high-coherence point targets in a large amount of data, the surface deformation rate and time deformation sequence of the study area are obtained (Berardino P. et al., 2002). This technology is a data post-processing technology based on the least squares method, which can effectively improve the reliability of the results and can be used to monitor long-term slow deformation of the surface (Wang Chao et al., 2002).

[0003] In 1973, Golub and Pereyra proposed the variable projection method (VP) for the separable nonlinear least squares problem (Golub GH et al., 1973). Since then, many scholars have conducted in-depth research on this algorithm that separates variables and then solves them (Kaufman L., 1975; Ruhe A. et al., 1980; Corradi C., 1981; Dianne P. O'Leary et al., 2013). This algorithm has been widely developed and applied in the fields of numerical analysis, neural networks, electronic communications and engineering, environmental science and time series analysis. Scholars have made many attempts to use the variable projection method in different scientific research fields and have achieved relatively rich research results. The commonly used model for predicting surface subsidence based on measured data is the Logistic curve equation, and the scientificity and rationality of this model have been verified through a large number of engineering examples (Zhu Zhiduo et al., 2009). However, the Logistic model is mostly used to predict the subsidence process in the time series and cannot fit the surface state in the spatial position. Summary of the invention

[0004] The purpose of the present invention is to overcome the above-mentioned defects and provide an application method of the variable projection method based on QR decomposition in surface deformation. The fitting of the surface deformation in time and space of the study area is taken as the research focus. On the basis of applying SBAS-InSAR time series analysis technology to explore the surface deformation in the study area, a complex exponential model with a separable structure is used to fit the settlement situation of the surface in spatial position. The model parameters are solved by the variable projection algorithm based on QR decomposition, and the settlement profile curve is approximated with high precision in spatial position. At the same time, the Logistic model is combined to predict the surface in time, thereby realizing comprehensive monitoring of the surface deformation status in time and space, and providing reliable technical support for surface deformation monitoring.

[0005] In order to achieve the above-mentioned object of the invention, the present invention provides the following technical solutions:

[0006] An application of a variable projection method based on QR decomposition in surface deformation includes the following steps:

[0007] S1 uses the SBAS-InSAR time series analysis method to generate matching image pairs by setting time thresholds, thereby obtaining the original surface subsidence data;

[0008] S2 uses the original surface settlement data and the variable projection algorithm based on QR decomposition to fit the surface settlement in space. The specific method is as follows:

[0009] S2.1 Substituting the original surface settlement data into the separable exponential model to obtain the separable exponential model residual function; the separable exponential model residual function is in the form of a separable nonlinear equation system;

[0010] S2.2 Minimize the residual function of the separable exponential model under the least squares principle, and use LM iteration to obtain the nonlinear parameter estimation in the residual function of the separable exponential model

[0011] S2.3 Computing linear parameter estimates in the residual function of the separable exponential model

[0012] S2.4 Estimate the nonlinear parameters and linear parameter estimates Substituting the separable index model into the model, a solved separable index model is obtained, and a surface settlement profile and a settlement fitting diagram are drawn according to the solved separable index model;

[0013] S3 uses the original surface settlement data to fit the surface settlement in the time series based on the Logistic model. The specific method is as follows:

[0014] S3.1 Substitute the original surface settlement data into the Logistic model to obtain the Logistic model residual function; the Logistic model residual function is in the form of a nonlinear equation system;

[0015] S3.2 Minimize the residual function of the Logistic model under the least squares principle and use LM iteration to obtain the parameter estimates in the residual function of the Logistic model;

[0016] S3.3 Substitute the parameter estimates into the Logistic model to obtain the solved Logistic model, and draw a curve of surface settlement changing with time and a curve fitting diagram based on the solved Logistic model.

[0017] Furthermore, in step S1, the specific method of obtaining the original surface subsidence data using the SBAS-InSAR time series analysis method is:

[0018] S1.1 matches the remote sensing images by setting a time threshold to obtain matching image pairs;

[0019] S1.2 performing differential interference on the matching image pair to generate an interference image pair;

[0020] S1.3 performs filtering and enhancement processing on the interference pattern composed of the interference image pairs in sequence, and after phase unwrapping, estimates the surface deformation rate and residual topography according to the interference pattern, and calculates the surface deformation variable in the time series according to the surface deformation rate and residual topography;

[0021] S1.4 Remove the atmospheric phase from the surface shape variables in the time series and convert them into the geographic coordinate system to obtain the original surface settlement data.

[0022] Furthermore, in step S2.1, the separable index model is:

[0023]

[0024] The residual function of the separable exponential model is: r = y-η = y-Φβ;

[0025] Where y is the surface settlement value, η is the estimated surface settlement including the spatial position variable u, the nonlinear parameter α and the linear parameter β, and the nonlinear function matrix r is the residual of the separable exponential model, i represents the number of nonlinear functions in the model, i = 1, 2, ... n; the nonlinear parameter α includes

[0026] λ and z, α=(λ 1 ,λ 2 ,…,λ n ; z), β=(β 0 ,β 1 ,…,β n ).

[0027] Furthermore, in step S2.2, the nonlinear parameter estimation in the residual function of the separable exponential model is obtained by LM iteration. The method is:

[0028] S2.2.1 Solve the first-order partial derivative of the residual function of the separable exponential model with respect to the nonlinear parameter α to obtain the first Jacobian matrix J;

[0029] S2.2.2 Given the initial value of nonlinear parameter α 0 and the first damping factor adjustment coefficient τ, set the convergence criteria: the first iteration termination threshold ε and the first maximum number of iterations k max ;

[0030] S2.2.3 α 0 Substitute the first Jacobian matrix J and the nonlinear function matrix Φ, perform QR decomposition on the nonlinear function matrix Φ, and obtain the first variable projection function based on QR decomposition

[0031] Let the first damping factor μ = τ·max{diag(J T J)}, calculate the first iteration step length h LM =-(J T J+μI) -1 J T r, α k+1 =αk +h LM , k ≥ 0;

[0032] S2.2.4 If h LM Make Then accept α k+1 And reduce the first damping factor μ and go to the next step; otherwise, increase the first damping factor μ and recalculate h LM , until

[0033] S2.2.5 If or k>k max , then the iteration is terminated and the α at this time is output k+1 , which is the nonlinear parameter estimation Otherwise, repeat steps S2.2.3 and S2.2.4.

[0034] Furthermore, in step S2.2.3, the nonlinear function matrix Φ is subjected to QR decomposition to obtain a variable projection function based on QR decomposition: The method is:

[0035] Perform QR decomposition on Φ to obtain the decomposition formula Φ=QR, let Q=[Q 1 Q 2 ], Q 1 , Q 2 They are m×p and m×(mo) order matrices respectively, m is the number of surface settlement data, o is the rank of Φ, and R is an upper triangular matrix;

[0036] Let the linear parameter β = Φ + y, where Φ + is the generalized inverse matrix of Φ;

[0037] Combine the decomposition and β=Φ + y obtains the nonlinear function matrix Φ based on the variable projection function of QR decomposition

[0038] Furthermore, in step S2.3, the linear parameter estimation in To estimate the nonlinear parameters Substitute the generalized inverse matrix of the nonlinear function matrix obtained.

[0039] Furthermore, in step S3.1, the Logistic model is

[0040] The residual function of the logistic model is r′=yS, where y is the surface settlement value, S is the time variable t and parameter A. 1 ,A 2 ,A3 is the estimated surface subsidence, and r′ is the residual of the Logistic model.

[0041] Furthermore, in step S3.2, the method for obtaining the parameter estimation in the residual function of the Logistic model by using LM iteration is:

[0042] S3.2.1 Solve the first-order partial derivative of the residual function of the Logistic model with respect to parameter A to obtain the second Jacobian matrix J′, where parameter A is the matrix containing parameter A. 1 ,A 2 and A 3 Parameters;

[0043] S3.2.2 Given the initial value of parameter A 0 and the second damping factor adjustment coefficient τ′, set the convergence criteria: the second iteration termination threshold ε′ and the second maximum number of iterations k′ max ;

[0044] S3.2.3 A 0 Substitute the second Jacobian matrix J′ and r′ to obtain the residual sum of squares

[0045] Let the second damping factor μ′=τ′·max{diag(J′ T J′)}, calculate the second iteration step length h′ LM =-(J′ T J′+μ′I) -1 J′ T r′,A k+1 =A k +h LM , k ≥ 0;

[0046] S3.2.4 If h′ LM Make Then accept A k+1 And reduce the second damping factor μ′ and proceed to the next step; otherwise, increase the second damping factor μ′ and recalculate h′ LM , until

[0047] S3.2.5 If or k>k′ max , then the iteration is terminated and the output is A at this time k+1 , which is the estimated value of parameter A Otherwise, repeat steps S3.2.3 and S3.2.4.

[0048] Furthermore, in step S1.2, the following formula is used to perform differential interference on the matching image pair to generate an interference image pair:

[0049] δφ j (p,q)=φ(t B ,p,q)-φ(t A ,p,q);

[0050] Among them, t A ,t B is the acquisition time of two SAR images, (p, q) represents the pixel position, φ represents the phase of the corresponding pixel, δφ represents the phase difference of the corresponding pixel, and j is the time at t A and t B The interferometric image pair sequence number after the interference of two SAR images obtained at the moment.

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

[0052] (1) Based on the surface deformation data obtained by the SBAS-InSAR method, the present invention adopts the variable projection method of QR decomposition and an exponential model to realize the spatial monitoring of the surface deformation conditions. At the same time, the Logistic model is used to perform temporal fitting of the subsidence area, thereby realizing the temporal monitoring of the surface deformation conditions. The monitoring results are accurate and comprehensive.

[0053] (2) The present invention utilizes a variable projection algorithm based on QR decomposition to separate the linear and nonlinear parameters in an exponential model with a separable structure and solve them in sequence, and performs geographical location fitting of the typical subsidence areas in the study area. The solution accuracy is high, and the calculation efficiency is improved while taking into account the accuracy.

[0054] (3) The present invention uses Qilian County, Qinghai Province as the research area to verify the technical solution, and uses the LM algorithm to solve the exponential model and the logistic model, which has a high fitting accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 A map of the Qilian County area and surrounding geographical environment in Example 1 of the present invention;

[0056] Figure 2 The time series surface deformation map and the annual average deformation rate map in Example 1 of the present invention are shown in FIG. Figure 2 (1) to (6) are time series surface deformation maps. Figure 2 (7) is the annual average deformation rate diagram;

[0057] Figure 3 The final surface deformation map in Example 1 of the present invention;

[0058] Figure 4 The separable exponential model is fitted to the sedimentation profile in Example 1 of the present invention;

[0059] Figure 5 It is a settlement profile prediction diagram in Example 1 of the present invention;

[0060] Figure 6 The logistic model is used to fit the time series sedimentation value in Example 1 of the present invention. DETAILED DESCRIPTION

[0061] The following detailed description of the present invention will make the features and advantages of the present invention more clear and explicit.

[0062] The word "exemplary" is used exclusively herein to mean "serving as an example, embodiment, or illustration." Any embodiment described herein as "exemplary" is not necessarily to be construed as preferred or advantageous over other embodiments. Although various aspects of the embodiments are shown in the drawings, the drawings are not necessarily drawn to scale unless otherwise noted.

[0063] Aiming at the lack of redundant observations in D-InSAR technology and the decoherence problem caused by a large time span, the present invention applies SBAS-InSAR time series analysis technology to surface deformation monitoring, monitors surface settlement changes based on remote sensing image data of different periods, and applies a separable nonlinear model, a Logistic model and a variable projection method based on QR decomposition to settlement profile fitting, providing technical support for monitoring surface deformation.

[0064] The present invention provides an application of a variable projection method based on QR decomposition in surface deformation, comprising the following steps:

[0065] (1) Apply the SBAS-InSAR time series analysis method to surface deformation monitoring to understand the surface subsidence conditions and laws. This step selects multiple images based on remote sensing data from different periods, pairs them according to the time and space baseline thresholds to generate interference, and performs filtering and phase unwrapping, and finally calculates the surface deformation amount and deformation rate.

[0066] The specific method of this step is:

[0067] The remote sensing images are matched by setting a time threshold, and the matched image pairs are connected by lines to generate a time-space baseline connection map. The matched image pairs are differentially interfered using formula (6) to generate interference image pairs:

[0068] δφ j (p,q)=φ(t B ,p,q)-φ(t A ,p,q) (6).

[0069] In the formula, t A ,t Bis the acquisition time of two SAR images, (p, q) represents the pixel position, φ represents the phase of the corresponding pixel, δφ represents the phase difference of the corresponding pixel, and j is the time at t A and t B The interferometric image pair sequence number after interferometry of two SAR images obtained at the moment. The interferogram composed of the interferometric image pair is filtered and enhanced to improve the interferometric quality, and phase unwrapping is performed to extract surface deformation information, estimate deformation rate and residual terrain, calculate deformation in time series, estimate and remove atmospheric phase, and convert the result into geographic coordinate system to obtain original surface settlement data.

[0070] (2) Using an exponential function model with a separable structure, the variable projection algorithm based on QR decomposition is applied to solve the linear and nonlinear parameters of the model and fit the settlement profile of the surface at the spatial position;

[0071] (2.1) Substituting the original surface settlement data into the separable exponential model, we obtain the residual function shown in formula (1), forming a separable nonlinear equation system. The separable exponential model expression is:

[0072]

[0073] The residual function of the separable exponential model is: r = y-η = y-Φβ; (1);

[0074] Where y is the surface settlement value, η is the estimated surface settlement including the spatial position variable u, the nonlinear parameter α and the linear parameter β, and the nonlinear function matrix r is the residual of the separable exponential model, i represents the number of nonlinear functions in the model, i = 1, 2, ... n; the nonlinear parameter α includes the nonlinear parameters λ and z, α = (λ 1 ,λ 2 ,…,λ n ; z), β=(β 0 ,β 1 ,…,β n ).

[0075] (2.2) Minimize the residual function of the separable exponential model under the least squares principle and use LM iteration to obtain the nonlinear parameter estimation

[0076] According to the least squares principle, the solution of linear and nonlinear parameters is to find the minimum sum of squares of the residual function of the separable exponential model. and Right now:

[0077]

[0078] Decompose the matrix Φ by QR decomposition method: Φ = QR, let Q = [Q 1 Q 2 ], Q 1 , Q 2 are m×p and m×(mo)-order matrices respectively, m is the number of surface subsidence data, o is the rank of Φ, and R is an upper triangular matrix; let the linear parameter β in formula (2) = Φ + y, combined with the QR decomposition of Φ, the variable projection function is as follows:

[0079]

[0080] In the formula, Φ + is the generalized inverse matrix of Φ.

[0081] Using LM iteration to obtain nonlinear parameter estimates The specific process is as follows:

[0082] (2.2.1) Take the first-order partial derivative of the residual function with respect to the nonlinear parameter and obtain the first Jacobian matrix J;

[0083] (2.2.2) Given the initial value of the nonlinear parameter α 0 , the first damping factor adjustment coefficient τ, set the convergence criteria - the first iteration termination threshold ε and the first maximum number of iterations k max ;

[0084] (2.2.3) Substitute the initial value α 0 Get the numerical matrices J and Φ, perform QR decomposition on Φ, and then get the first variable projection function based on QR decomposition as shown in formula (3): Let the first damping factor μ = τ·max{diag(J T J)}, calculate the first iteration step length h LM =-(J T J+uI) -1 J T r, α k+1 =α k +h LM ;

[0085] (2.2.4) If h LM Make Then accept α k+1 And reduce the first damping factor and go to the next step, otherwise, increase the first damping factor and recalculate h LM , until

[0086] (2.2.5) If or k>k max, the iteration is terminated, otherwise, (2.2.3) and (2.2.4) are repeated.

[0087] (2.3) The linear parameter β is solved as follows: β = Φ + y (4);

[0088] (2.4) Substitute the linear and nonlinear parameter estimates into the separable exponential model and draw the surface subsidence profile and subsidence fitting diagram.

[0089] (3) Use the Logistic model to fit the surface settlement over time;

[0090] (3.1) Substitute the original surface subsidence data into the Logistic model to obtain the residual function of the Logistic model and form a nonlinear equation system;

[0091] The Logistic model expression is:

[0092] The residual function of the logistic model is r′=ys (5);

[0093] Among them, y is the surface settlement value, s is the time variable t and parameter A 1 ,A 2 ,A 3 is the estimated surface subsidence, and r′ is the residual of the Logistic model.

[0094] (3.2) Minimize the residual function of the Logistic model under the least squares principle, and use LM iteration to obtain parameter estimates:

[0095] (3.2.1) Solve the first-order partial derivative of the residual function of the Logistic model with respect to parameter A to obtain the second Jacobian matrix J′; the parameter A is the matrix containing parameter A 1 ,A 2 and A 3 Parameters

[0096] (3.2.2) Given the initial value of parameter A 0 and the second damping factor adjustment coefficient τ′, set the convergence criteria: the second termination threshold ε′ and the second maximum number of iterations k′ max ;

[0097] (3.2.3) Substitute the initial value A 0 Get the numerical matrices J′ and r′, and get the residual sum of squares

[0098] Let the second damping factor μ′=τ′·max{diag(J′ T J′)}, calculate the second iteration step length hL ' M =-(J′ T J′+μ′I) -1 J′ T r′,A k+1 =A k +h L ' M , k ≥ 0;

[0099] (3.2.4) If h′ LM Make Then accept A k+1 And reduce the second damping factor μ′ and proceed to the next step; otherwise, increase the second damping factor μ′ and recalculate h′ LM , until

[0100] (3.2.5) If or k>k′ max , then the iteration is terminated and the output is A at this time k+1 , which is the estimated value of parameter A Otherwise, repeat steps (3.2.3) and (3.2.4).

[0101] (3.3) Substitute the estimated parameter values ​​into the Logistic model and draw the curve of surface subsidence changing with time and the curve fitting diagram.

[0102] Embodiment 1:

[0103] 1. Obtaining deformation information using SBAS-InSAR technology

[0104] The case study takes Qilian County, Qinghai Province as the research area. The geographical environment of Qilian County and its surrounding areas is as follows: Figure 1 The SAR data used are 7 Sentinel-1 ascending orbit data, C band, VV polarization mode, image coverage range is 37°55′N~38°25′N, 99°40′E~100°30′E, and the time span is from February 19, 2019 to June 25, 2020. The time information of the image data is shown in Table 1. The data was processed based on SBAS-InSAR technology, and 6 surface deformation maps and 1 annual average deformation rate map of the study area in the time series were obtained, as shown in Figure 2 shown.

[0105] Table 1 Experimental data time information

[0106]

[0107] In the surface deformation maps and annual average deformation rate maps, uplift areas are represented by positive values ​​and subsidence areas are represented by negative values. Figure 2From the annual average deformation rate diagram in , we can see that the uplift rate of cities and surrounding areas is relatively fast, while the mountainous areas are mainly subsiding, and the surface deformation rate in most areas is concentrated between 0 and 20 mm / year, with a relatively slow subsidence rate, indicating that the surface is relatively stable. From the time series deformation diagram, we can see that the surface of the mountainous area has begun to deform significantly since December 2019, and the deformation has gradually increased over time. Figure 2 (6) Surface deformation map. The subsidence area and the uplift area in the study area are clearly divided. The urban area mainly shows uplift, while the mountainous area mainly shows subsidence. This is consistent with the conclusion drawn from the annual average deformation rate map. The maximum subsidence value in the study area exceeds 230 mm, but the subsidence values ​​in most subsidence areas are concentrated in the range of [-60,0]. In order to more intuitively display the surface subsidence in the study area, a line segment QL is selected in the study area to run through the urban area of ​​Qilian County. The geographic coordinates of the two end points of the line segment are (38.201°N, 100.231°E) and (38.162°N, 100.264°E), respectively. Figure 3 As shown in , the settlement value of each point on the line segment is extracted, the separable exponential model is used to fit and analyze the settlement of the spatial position, the settlement profile is drawn, and the Logistic model is used to fit the settlement value of a point on the line segment in the time series.

[0108] 2. Surface settlement fitting in spatial position based on variable projection method based on QR decomposition

[0109] The following separable exponential model was used to fit the sedimentation curve:

[0110]

[0111] Where α=(λ 1 ,λ 2 ,…,λ n ;z 1 ,z 2 ,…,z n ) T , This model has excellent performance in curve fitting and has a structure in which linear parameters and nonlinear parameters can be separated, so it can be solved using the variable projection algorithm. Let n=4 in the separable exponential model, then there are 9 nonlinear parameters and 5 linear parameters in the separable exponential model. The variable projection algorithm is used to fit the first 300 sets of settlement value data extracted on the line segment QL and predict the surface deformation, and the prediction results are verified using the 301st to 400th sets of data. The data are substituted into the separable exponential model to form a nonlinear equation group. The same nonlinear parameter iteration initial value of 0.5 is given. Based on the QR decomposition variable projection algorithm (VP QR ) fitting surface settlement method, singular value decomposition variable projection algorithm (VPSVD ) and full rank decomposition variable projection algorithm (VP FRD ), and use LM to iteratively solve the parameters. The linear parameters are obtained by formula (4).

[0112] The results of fitting the settlement values ​​of each point with the separable index model are as follows Figure 4 Table 2 lists the number of iterations, residual 2-norm (r 2-norm ), root mean square error (RMSE) and operation time. Figure 4 It can be seen that the surface settlement in the Qilian County area where the line segment passes is relatively small, with the maximum settlement being about 0.07 mm and some areas experiencing uplift. FRD The solution result is poor, and VP QR With VP SVD The corresponding curves of the algorithm are almost coincident, and the fitting degree with the original data is high, and the root mean square error value is on the order of 10 -6 , the deformation value of each point and the estimated error are extremely small, indicating that the parameter solution accuracy is high and the model has a good fit. Under the premise of equivalent solution accuracy, VP QR The algorithm has the shortest operation time and higher calculation efficiency. QR The algorithm has good applicability in fitting the surface subsidence profile at the spatial position.

[0113] Table 2 Iteration times, residual 2-norm (r 2-norm ), root mean square error (RMSE) and calculation time

[0114]

[0115] The prediction results of the separable index model for the settlement value of each point are as follows: Figure 5 Table 3 lists the corresponding residual 2-norm (r 2-norm ) and root mean square error (RMSE). Figure 5 It can be seen that VP FRD The prediction effect is poor, while VP QR With VP SVD The prediction curve corresponding to the algorithm has a consistent trend with the original settlement data, and its prediction accuracy is 10 -6 , the error between the deformation value and its predicted value is extremely small.

[0116] Table 3 Residual norm of model prediction (r 2-norm ) and root mean square error (RMSE)

[0117]

[0118] 3. Sedimentation Fitting of Logistic Model in Time Series

[0119] Select the 417th pixel on the line segment QL, extract the settlement value of the pixel in each period, and obtain 6 settlement data. The Logistic model is as shown in formula (5). Substitute the 6 settlement values ​​into the Logistic model, give the iterative initial value of the Logistic model parameters, use the LM iterative algorithm to solve, and fit the settlement value of each period at this point, as shown in Figure 6 As shown in the figure, it can be judged that the curve has a good fitting effect. The Logistic model is not only suitable for fitting the surface settlement in the time series, but also can be used to predict the settlement trend of ground points over time.

[0120] The present invention has been described in detail above in conjunction with specific implementations and exemplary examples, but these descriptions cannot be understood as limiting the present invention. Those skilled in the art understand that, without departing from the spirit and scope of the present invention, a variety of equivalent substitutions, modifications or improvements may be made to the technical solution of the present invention and its implementation methods, all of which fall within the scope of the present invention. The scope of protection of the present invention shall be subject to the attached claims.

[0121] The contents not described in detail in the specification of the present invention belong to the common knowledge of those skilled in the art.

Claims

1. An application of variable projection method based on QR decomposition in surface deformation, characterized in that: The following steps are involved: S1 uses the SBAS-InSAR time series analysis method to generate matching image pairs by setting time thresholds, and then obtains the original surface subsidence data; S2 uses the original surface settlement data and the variable projection algorithm based on QR decomposition to fit the surface settlement in space. The specific method is as follows: S2.1 Substituting the original surface settlement data into the separable exponential model to obtain the separable exponential model residual function; the separable exponential model residual function is in the form of a separable nonlinear equation system; S2.2 Minimize the residual function of the separable exponential model under the least squares principle, and use LM iteration to obtain the nonlinear parameter estimation in the residual function of the separable exponential model S2.3 Computing linear parameter estimates in the residual function of the separable exponential model S2.4 Estimate the nonlinear parameters and linear parameter estimates Substituting the separable index model into the model, a solved separable index model is obtained, and a surface settlement profile and a settlement fitting diagram are drawn according to the solved separable index model; S3 uses the original surface settlement data to fit the surface settlement in the time series based on the Logistic model. The specific method is as follows: S3.1 Substitute the original surface settlement data into the Logistic model to obtain the Logistic model residual function; the Logistic model residual function is in the form of a nonlinear equation system; S3.2 Minimize the residual function of the Logistic model under the least squares principle and use LM iteration to obtain the parameter estimates in the residual function of the Logistic model; S3.3 Substitute the parameter estimates into the Logistic model to obtain the solved Logistic model, and draw a curve of surface settlement changing with time and a curve fitting diagram based on the solved Logistic model.

2. The application of the variable projection method based on QR decomposition in surface deformation according to claim 1 is characterized in that: In step S1, the specific method of obtaining the original surface subsidence data using the SBAS-InSAR time series analysis method is: S1.1 matches the remote sensing images by setting a time threshold to obtain matching image pairs; S1.2 performing differential interference on the matching image pair to generate an interference image pair; S1.3 performs filtering and enhancement processing on the interference pattern composed of the interference image pairs in sequence, and after phase unwrapping, estimates the surface deformation rate and residual topography according to the interference pattern, and calculates the surface deformation variable in the time series according to the surface deformation rate and residual topography; S1.4 Remove the atmospheric phase from the surface shape variables in the time series and convert them into the geographic coordinate system to obtain the original surface settlement data.

3. The application of the variable projection method based on QR decomposition in surface deformation according to claim 1 or 2 is characterized in that: In step S2.1, the separable index model is: The residual function of the separable exponential model is: r = y-η = y-Φβ; Where y is the surface settlement value, η is the estimated surface settlement including the spatial position variable u, the nonlinear parameter α and the linear parameter β, and the nonlinear function matrix r is the residual of the separable exponential model, i represents the number of nonlinear functions in the model, i = 1, 2, ... n; the nonlinear parameter α includes λ and z, α = (λ1, λ2, ..., λ n ; z), β=(β0,β1,…,β n ).

4. The application of the variable projection method based on QR decomposition in surface deformation according to claim 3 is characterized in that: In step S2.2, the nonlinear parameter estimation in the residual function of the separable exponential model is obtained by LM iteration. The method is: S2.2.1 Solve the first-order partial derivative of the residual function of the separable exponential model with respect to the nonlinear parameter α to obtain the first Jacobian matrix J; S2.2.2 Given the initial value of the nonlinear parameter α0 and the first damping factor adjustment coefficient τ, set the convergence criteria: the first iteration termination threshold ε and the first maximum number of iterations k max ; S2.2.3 Substitute α0 into the first Jacobian matrix J and the nonlinear function matrix Φ, perform QR decomposition on the nonlinear function matrix Φ, and obtain the first variable projection function based on QR decomposition Let the first damping factor μ = τ·max{diag(J T J)}, calculate the first iteration step length h LM =-(J T J+μI) -1 J T r, α k+1 =α k +h LM , k ≥ 0; S2.2.4 If h LM Make Then accept α k+1 And reduce the first damping factor μ and proceed to the next step; Otherwise, increase the first damping factor μ and recalculate h LM , until S2.2.5 If or k>k max , then the iteration is terminated and the α at this time is output k+1 , which is the nonlinear parameter estimation Otherwise, repeat steps S2.2.3 and S2.2.

4.

5. The application of the variable projection method based on QR decomposition in surface deformation according to claim 4 is characterized in that: In step S2.2.3, the nonlinear function matrix Φ is subjected to QR decomposition to obtain a variable projection function based on QR decomposition The method is: Perform QR decomposition on Φ to obtain the decomposition formula Φ=QR, let Q=[Q1Q2], Q1 and Q2 are m×p and m×(mo)-order matrices respectively, m is the number of surface settlement data, o is the rank of Φ, and R is an upper triangular matrix; Let the linear parameter β = Φ + y, where Φ + is the generalized inverse matrix of Φ; Combine the decomposition and β=Φ + y obtains the nonlinear function matrix Φ based on the variable projection function of QR decomposition 6. The application of the variable projection method based on QR decomposition in surface deformation according to claim 5 is characterized in that: In step S2.3, the linear parameter estimation in To estimate the nonlinear parameters Substitute the generalized inverse matrix of the nonlinear function matrix obtained.

7. The application of the variable projection method based on QR decomposition in surface deformation according to claim 1 or 2, characterized in that: In step S3.1, the Logistic model is The residual function of the Logistic model is r′=ys, where y is the surface settlement value, s is the estimated surface settlement including the time variable t and parameters A1, A2, A3, and r′ is the residual of the Logistic model.

8. The application of the variable projection method based on QR decomposition in surface deformation according to claim 7 is characterized in that: In step S3.2, the method for obtaining the parameter estimation in the residual function of the Logistic model by using LM iteration is: S3.2.1 solving the first-order partial derivative of the residual function of the Logistic model with respect to parameter A to obtain the second Jacobian matrix J′, wherein parameter A is a parameter including parameters A1, A2 and A3; S3.2.2 Given the initial parameter value A0 and the second damping factor adjustment coefficient τ′, set the convergence criteria: the second iteration termination threshold ε′ and the second maximum number of iterations k′ max ; S3.2.3 Substitute A0 into the second Jacobian matrix J′ and r′ to obtain the residual sum of squares Let the second damping factor μ′=τ′·max{diag(J′ T J′)}, calculate the second iteration step length h′ LM =-(J′ T J′+μ′I) - 1 J′ T r′,A k+1 =A k +h LM , k ≥ 0; S3.2.4 If h′ LM Make Then accept A k+1 And reduce the second damping factor μ′ and proceed to the next step; otherwise, increase the second damping factor μ′ and recalculate h′ LM , until S3.2.5 If or k>k′ max , then the iteration is terminated and the output is A at this time k+1 , which is the estimated value of parameter A Otherwise, repeat steps S3.2.3 and S3.2.

4.

9. The application of the variable projection method based on QR decomposition in surface deformation according to claim 2 is characterized in that: In step S1.2, the following formula is used to perform differential interference on the matching image pair to generate an interference image pair: δφ j (p,q)=φ(t B ,p,q)-φ(t A ,p,q); Among them, t A ,t B is the acquisition time of two SAR images, (p, q) represents the pixel position, φ represents the phase of the corresponding pixel, δφ represents the phase difference of the corresponding pixel, and j is the time at t A and t B The interferometric image pair sequence number after the interference of two SAR images obtained at the moment.

Images