A high-precision 3D deformation inversion method based on spatiotemporal continuity of BeiDou InSAR

Abstract

This invention discloses a high-precision three-dimensional deformation inversion method based on spatiotemporal continuity of BeiDou InSAR. By using the spatiotemporal continuity of deformation as a constraint and high-precision observation points as constraint points to correct the observation results of other points, it can achieve high-precision three-dimensional deformation inversion even when the number of satellites participating in deformation inversion is small. This solves the problem of low accuracy of three-dimensional deformation inversion caused by insufficient number of observation satellites in the BeiDou bistatic InSAR system, and improves the deformation detection accuracy and applicable scenarios of the BeiDou InSAR system.

Description

Technical Field

[0001] This invention belongs to the field of bistatic synthetic aperture radar technology, specifically relating to a three-dimensional high-precision deformation inversion method for BeiDou InSAR based on spatiotemporal continuity. Background Technology

[0002] The BeiDou-InSAR (BeiDou based Interferometric Synthetic Aperture Radar System) system can be used for three-dimensional deformation inversion. This system utilizes in-orbit BeiDou satellites as transmitters and deploys geostationary receivers on the ground to form a bistatic SAR system, such as... Figure 1 As shown, deformation monitoring is then achieved using heavy orbit SAR images. This system inherits the advantages of the BeiDou positioning system and radar system, and can achieve three-dimensional deformation measurement of the opposite scene with a single device. Compared with traditional deformation detection methods, it has advantages such as low cost and short monitoring cycle.

[0003] Achieving 3D deformation inversion requires combining observation information from multiple different angles. However, due to the varying scattering characteristics of the monitored scene at different angles, the number and distribution of satellite imagery points (PS points) differ. Consequently, the number of satellites that can observe different targets varies during multi-satellite joint processing. In 3D deformation inversion, the more effective observation angles a target has, the less noise is affected, and the higher the deformation inversion accuracy. Conversely, when the effective observation angles of a target are limited, such as when a target can only be observed by three or four satellites, its 3D accuracy will fail to meet technical requirements. Summary of the Invention

[0004] In view of this, the present invention provides a high-precision three-dimensional deformation inversion method based on spatiotemporal continuity of BeiDou InSAR. It uses the spatiotemporal continuity of deformation as a constraint and high-precision observation points as constraint points to correct the observation results of other points. Even when the number of satellites participating in deformation inversion is small, high-precision three-dimensional deformation inversion can be achieved.

[0005] The present invention provides a three-dimensional high-precision deformation inversion method based on spatiotemporal continuity of BeiDou InSAR, which includes the following steps:

[0006] Least squares estimation is used to solve the 3D deformation results of the scene under multi-satellite observations. The PS point set is divided into high-precision PS point set and low-precision PS point set according to the accuracy. The expected deformation of the whole scene is obtained by interpolation using the high-precision PS point set data. For the low-precision PS point set, the difference between the deformation estimate of the least squares estimation and the expected deformation is used to obtain the penalty function. The penalty function is used to optimize the result of the least squares estimation. Different penalty function coefficients are selected for different PS points in the low-precision PS point set to complete the constrained least squares estimation of the PS points of the whole scene, and finally obtain the high-precision 3D deformation inversion result of the whole scene.

[0007] Furthermore, the method for using least squares estimation to solve for the three-dimensional deformation results of the scene under multi-star observation is as follows:

[0008] The relationship between the observations obtained from different angles by the navigation satellite and the deformation is expressed as: Φ M×1 =H M×3 ·D 3×1 +n M×1 ,in:

[0009]

[0010]

[0011] D 3×1 =[D x D y D z ] T

[0012] n M×1 =[n1 n2…n M ] T

[0013] Φ M×1 For the observation results of M satellites, H M×3 D is the matrix of deformation measurement results. 3×1 Let n be the true shape variable matrix of the target. M×1 For the observation noise of M satellites, P s For satellite position, P E For the receiver position, P Q For the target location;

[0014] The objective function is: ε 2 =||Φ-H·D|| 2 , where ε represents the difference;

[0015] The estimation result of D obtained by least squares estimation is as follows:

[0016] Let the result set of multi-angle correlation be The three-dimensional deformation is obtained by least-squares estimation for each point in the point set:

[0017] Furthermore, the method of using the penalty function to optimize the least squares estimation result is as follows:

[0018] The penalty function is: Let the expected deformation variables be the objective function for optimizing the least squares estimate using the penalty function. The obtained deformation estimation results are It is a low-precision PS point set.

[0019] Furthermore, the penalty function coefficients are determined as follows:

[0020] Let the PS point set on day q-1 be The expected value of the deformation is The final deformation inversion data are The neighborhood S(A) of target point A is defined as S(A) = {B||A,B|<r}, where B is a neighboring point of A and r is the radius of the neighborhood. The standard deviation St between the actual deformation and the predicted value within the neighborhood of target point A is also defined. t q-1 (A) is:

[0021]

[0022] Let the PS point set on day q be The set of observation satellites for target point A is S. a q (A), then we have:

[0023] Step 4.1: Based on the observed satellite set S a q (A), obtain the transformation matrix H on day q. q (A);

[0024] Step 4.2, with S t q-1 (A) As the expected constrained least squares output, calculate the observation Φ for each star. q (A): Φ q (A) = H q (A)×S t q-1 (A)+n, where n is Gaussian noise with a mean of 0;

[0025] Step 4.3, let For each value of k, compute its constrained least squares solution:

[0026]

[0027] The observation error is:

[0028] k q The estimation result for (A) is: k q (A) = arg min(|err q (A)|);

[0029] Step 4.4: Perform multiple Monte Carlo experiments to modify the error. Take the penalty function coefficient that minimizes the standard deviation error between the obtained deformation inversion accuracy and the target point A at day q-1 as the penalty function coefficient of target point A.

[0030] Furthermore, the method of interpolating the high-precision PS point set data to obtain the expected deformation of the entire scene is as follows: interpolating the high-precision PS point set data using the Kriging interpolation method.

[0031] Beneficial effects:

[0032] This invention solves the problem of low accuracy in three-dimensional deformation inversion caused by insufficient number of observation satellites in the BeiDou bistatic InSAR system, and improves the deformation detection accuracy and applicable scenarios of the BeiDou InSAR system. Attached Figure Description

[0033] Figure 1 A schematic diagram of the configuration of the BeiDou satellite bistatic SAR system used in the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by the present invention.

[0034] Figure 2 This is a flowchart illustrating the three-dimensional high-precision deformation inversion method of BeiDou InSAR based on spatiotemporal continuity provided by the present invention.

[0035] Figure 3 This is a schematic diagram of the penalty function coefficient selection process in the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by the present invention.

[0036] Figure 4 The image shows the direct imaging results of the deformation scene using the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by this invention.

[0037] Figure 5 The diagram shows the change in deformation accuracy in the east direction before and after compensation using the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by this invention.

[0038] Figure 6The diagram shows the change in deformation accuracy in the north direction before and after compensation using the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by this invention.

[0039] Figure 7 This is a diagram showing the change in deformation accuracy between the past and future directions, obtained by using the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by this invention. Detailed Implementation

[0040] The following examples illustrate the invention in detail.

[0041] The core idea of ​​this invention, a high-precision 3D deformation inversion method based on spatiotemporal continuity using BeiDou InSAR, is to correct the deformation inversion results of low-precision points using deformation inversion results from high-precision PS points (points observable by multiple satellites). Specifically, least squares estimation is used to obtain the 3D deformation results under multi-satellite observations. PS point sets are then divided according to their precision. The expected deformation of the entire scene is obtained by interpolation using the high-precision PS point set data. For the low-precision PS point set, the difference between the least squares estimated deformation and the expected deformation is used to obtain a penalty function, which is then used to optimize the least squares estimation result. However, when the overall scene deformation is large, the spatial correlation of deformation between neighboring points weakens, leading to a significant difference between the interpolated expected deformation value and the actual value. Conversely, when the scene deformation is small, the expected deformation value is more accurate, and in this case, the weight of the expected value in the objective function should be increased. Therefore, for different deformation states of the scene, different penalty function coefficients are selected, and high-precision full-scene 3D deformation inversion results are obtained through least squares estimation, thereby improving the accuracy of deformation inversion and effective early warning of deformation.

[0042] The present invention provides a three-dimensional high-precision deformation inversion method based on spatiotemporal continuity of BeiDou InSAR, the process of which is as follows: Figure 2 As shown, the specific steps include:

[0043] Step 1: Use least squares estimation to obtain the three-dimensional deformation results under multi-satellite observations.

[0044] The navigation satellite obtains three-dimensional deformation results through observations from multiple angles. The relationship between the deformation and the observation at each angle is as follows:

[0045] Φ M×1 =H M×3 ·D 3×1 +n M×1 (1)

[0046] in:

[0047]

[0048]

[0049] D 3×1 =[D x D y D z ] T

[0050] n M×1 =[n1 n2…n M ] T

[0051] Φ M×1 For the observation results of M satellites, H M×3 D is the matrix of deformation measurement results. 3×1 Let n be the true shape variable matrix of the target. M×1 For the observation noise of M satellites, P s For satellite position, P E For the receiver position, P Q The target location.

[0052] The objective function is:

[0053] ε 2 =||Φ-H·D|| 2 (2)

[0054] Here, ε represents the difference. The least squares method can be used to estimate D. for:

[0055]

[0056] Let the result set of multi-angle correlation be The three-dimensional deformation is obtained by applying the least squares method to each point in the point set, as shown in the following formula:

[0057]

[0058] Step 2: Divide the PS point set according to its accuracy, and use the high-precision PS point set data for interpolation to obtain the expected deformation of the entire scene.

[0059] Specifically:

[0060] Correlated points are divided according to the number of observed stars:

[0061]

[0062] in, For a high-precision point set, This is a point set with low precision.

[0063] Using a high-precision point set, kriging interpolation is used to obtain the deformation of the entire scene, which is then used as the expected value of the deformation:

[0064]

[0065] in, This represents the estimated location of the target, z(x). i y i ) represents the known quantities of points around the target location, λ i is a coefficient.

[0066] λ i This can be obtained by solving the following system of equations:

[0067]

[0068] Where, r ij Representing point (x) i y i ) and point (x) j y j The semivariance fit value between the two is calculated as follows:

[0069] First, calculate the initial values ​​of the pairwise distances and semivariances in the observed data:

[0070]

[0071] Fitting and This allows us to calculate the semivariance value corresponding to any distance:

[0072]

[0073] Furthermore, the fitted value r of the semivariance is... ij Substituting the values, we can calculate the coefficient λ. i Then, according to formula (6), the final estimation result is obtained. By interpolating the deformation variables in the X, Y, and Z directions respectively, the expected deformation variable result for the entire scene can be obtained, as shown in the following formula:

[0074]

[0075] Step 3: For a low-precision PS point set, the difference between the least squares estimated deformation estimate and the expected deformation is used to obtain the penalty function, and the penalty function is used to optimize the least squares estimation result.

[0076] The specific process is as follows:

[0077] Using penalty function The estimation results of constrained least squares are given, and the objective function is set as follows:

[0078]

[0079] make Then formula (11) can be rewritten as:

[0080]

[0081] The estimation result for D1 is as follows:

[0082]

[0083] For low-precision PS point sets The final deformation estimation result is:

[0084]

[0085] Step 4: Select different penalty function coefficients for different PS points in the low-precision PS point set to complete the constrained least squares estimation of PS points for the entire scene, and finally obtain high-precision full-scene 3D deformation inversion results.

[0086] The present invention determines the penalty function coefficients of each PS point within a low-precision PS point set by utilizing the temporal continuity of scene deformation, estimating the current day's deformation based on the scene's deformation from the previous day, and then determining the penalty function coefficients for that day. The solution process is as follows: Figure 3 As shown.

[0087] Let the PS point set of day q-1 be... Let point PS be point A. According to step 2, we can obtain... The expected value of deformation at any PS point in the point set on day q-1 and the final deformation inversion result are used to obtain the standard deviation between the actual deformation and the predicted value in the neighborhood S(A) of the target point A.

[0088] Since the deformation is continuous in time, the standard deviation of the target point A on day q-1 can characterize the discreteness of the deformation on day q. Using this as the expected least squares output, and combining it with the observed satellite set for point A on day q, we can obtain the differential phase of point A observed by each satellite. The obtained differential phase of point A is then used as the input to constrained least squares, for each penalty function coefficient [0, k... max The deformation inversion accuracy under the penalty function coefficient is obtained by solving the problem. Multiple Monte Carlo experiments are conducted, and the penalty function coefficient that minimizes the standard deviation error between the obtained deformation inversion accuracy and the target point A at day q-1 is taken as the penalty function coefficient of the PS point.

[0089] By performing the above operation on each PS point in the low-precision point set, the penalty function coefficients of each PS point in the low-precision PS point set can be obtained.

[0090] Let the PS point set on day q-1 be The expected value of the deformation and the final deformation inversion data are respectively and For target point A, its neighborhood S(A) is defined as:

[0091] S(A)={B||A,B|<r} (15)

[0092] Where B represents a neighboring point of A, and r is the radius.

[0093] Based on the predicted and actual values ​​from the previous day, the standard deviation between the actual deformation and the predicted value within the vicinity of target A can be obtained:

[0094]

[0095] Since the deformation is continuous in time, the S of the previous day... t q-1 It can characterize the discreteness of the data for the next day, that is, the value of the penalty function coefficient k for the next day can be determined according to S. t q-1 That's for you to decide.

[0096] Let the PS point set on day q be For any PS point The set of satellites observed at point A can be obtained as S. a q (A) and the expected discrete case S t q-1 (A). According to S a q (A) and S t q-1 (A) Determine the penalty function coefficient k for the day. q (A)

[0097] The specific steps are as follows:

[0098] S1, based on satellite set S a q (A), obtain the transformation matrix H on day q. q (A)

[0099] S2, with S t q-1 (A) As the expected least-squares output, calculate the observations for each star:

[0100] Φ q (A) = H q (A)×S t q-1 (A)+n (17)

[0101] Where n is Gaussian noise with a mean of 0.

[0102] S3, Order For each value of k, calculate its constrained least squares solution:

[0103]

[0104] The observation error is:

[0105]

[0106] k q The estimation result for (A) is:

[0107] k q (A) = argmin(err) q (A)) (20)

[0108] S4. Adjust the error and conduct multiple Monte Carlo experiments to obtain the final k. q The estimation result of (A) is obtained. Finally, P is obtained. b q The penalty function coefficient k at each point within the range q .

[0109] This allows for constrained least squares estimation of PS points across the entire scene.

[0110] Example:

[0111] In this embodiment, a 600m × 500m slope terrain is used as the deformation scenario, such as... Figure 4 As shown, the deformation is calculated using the BeiDou InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity provided by this invention.

[0112] The deformation direction is assumed to be slope aspect, and the input data is the deformation measured by other equipment for 11 days from May 24 to June 6 (excluding May 31, June 1, and June 2).

[0113] 800 PS points were selected from the data, and the number of points that can be observed by 4-8 stars is shown in Table 1:

[0114] Table 1 Number of PS points under different numbers of observation satellites

[0115]

[0116] Using 21 points associated with 8 stars as reference points, the deformation inversion results of stars 4-7 were optimized.

[0117] The accuracy comparison results before and after compensation in the three directions of East, North, and Sky are as follows: Figure 5 , Figure 6 and Figure 7 As shown.

[0118] The accuracy comparison results before and after compensation in the three directions are as follows:

[0119] Table 2 Accuracy Comparison

[0120]

[0121] As can be seen from the results in Table 2, the accuracy of the three-dimensional deformation inversion after processing by the proposed BeiDou satellite bistatic InSAR three-dimensional high-precision deformation inversion method based on spatiotemporal continuity is improved compared with that before processing, which proves the effectiveness of the present invention.

[0122] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A three-dimensional high-precision deformation inversion method based on spatiotemporal continuity of BeiDou InSAR, characterized in that, Includes the following steps: Least squares estimation is used to solve for the 3D deformation of the scene under multi-satellite observations. The PS point set is divided into high-precision PS point set and low-precision PS point set according to the accuracy. The expected deformation of the whole scene is obtained by interpolation using the high-precision PS point set data. For the low-precision PS point set, the difference between the deformation estimate obtained by least squares estimation and the expected deformation is used to obtain the penalty function. The penalty function is used to optimize the result of least squares estimation. Different penalty function coefficients are selected for different PS points in the low-precision PS point set to complete the constrained least squares estimation of the PS points of the whole scene, and finally, a high-precision 3D deformation inversion result of the whole scene is obtained. The method for solving the three-dimensional deformation results of the scene under multi-star observation using least squares estimation is as follows: The relationship between the observations obtained from different angles by the navigation satellite and the deformation is expressed as: Φ M×1 =H M×3 ·D 3×1 +n M×1 ,in: D 3×1 =[D x D y D z ] T n M×1 =[n1n2…n M ] T Φ M×1 For the observation results of M satellites, H M×3 D is the matrix of deformation measurement results. 3×1 Let n be the true shape variable matrix of the target. M×1 For the observation noise of M satellites, P s For satellite position, P E For the receiver position, P Q For the target location; The objective function is: ε 2 =||Φ-H·D|| 2 , where ε represents the difference; The estimation result of D obtained by least squares estimation is as follows: Let the result set of multi-angle correlation be The three-dimensional deformation is obtained by least-squares estimation for each point in the point set: The penalty function coefficients are determined as follows: Let the PS point set on day q-1 be The expected value of the deformation is The final deformation inversion data are The neighborhood S(A) of target point A is defined as S(A) = {B||A,B|<r}, where B is a neighboring point of A and r is the radius of the neighborhood. The standard deviation St between the actual deformation and the predicted value within the neighborhood of target point A is also defined. t q-1 (A) is: Let the PS point set on day q be The set of observation satellites for target point A is S. a q (A), then we have: Step 4.1: Based on the observed satellite set S a q (A), obtain the transformation matrix H on day q. q (A); Step 4.2, with S t q-1 (A) As the expected constrained least squares output, calculate the observation Φ for each star. q (A): Φ q (A)=H q (A)×S t q-1 (A)+n, where n is Gaussian noise with a mean of 0; Step 4.3, let For each value of k, compute its constrained least squares solution: The observation error is: k q The estimation result for (A) is: k q (A) = argmin(|err) q (A)|); Step 4.4: Perform multiple Monte Carlo experiments to modify the error. Take the penalty function coefficient that minimizes the standard deviation error between the obtained deformation inversion accuracy and the target point A at day q-1 as the penalty function coefficient of target point A.

2. The BeiDou InSAR three-dimensional high-precision deformation inversion method according to claim 1, characterized in that, The method of using a penalty function to optimize the least squares estimation result is as follows: The penalty function is: Let the expected deformation variables be the objective function for optimizing the least squares estimate using the penalty function. The obtained deformation estimation results are It is a low-precision PS point set.

3. The BeiDou InSAR three-dimensional high-precision deformation inversion method according to claim 1, characterized in that, The method for obtaining the expected deformation of the entire scene by interpolating high-precision PS point set data is as follows: interpolate the high-precision PS point set data using the Kriging interpolation method.

Images