PS deformation estimation method of navigation satellite bistatic InSAR (Interferometric Synthetic Aperture Radar) based on time coherence

Abstract

The invention belongs to the technical field of bistatic synthetic aperture radars, and particularly relates to a time-coherence-based PS deformation estimation method of a navigation satellite bistatic I nSAR, and the method comprises the specific steps: carrying out the PS point interferometric phase modeling based on the deformation phase, the atmospheric phase, the DEM error phase, the space baseline phase and the noise phase of a PS point, and carrying out the matrix processing of an interferometric phase model; the method comprises the following steps: aiming at a navigation satellite double-foundation InSAR image, combining an optical image, selecting a plurality of PS points which are relatively large in amplitude and have the number exceeding a set number of days according to an amplitude deviation method, and performing synchronous difference on phases of the plurality of PS points to obtain an interference phase; setting a speed, an acceleration and a DEM as cyclic variables, and setting a termination condition for terminating cyclic iteration; and when the termination condition is not met, compensating the interferometric phase by using the currently obtained cyclic variable, estimating and compensating the atmospheric phase and the constant phase by using the interferometric phase, and estimating the deformation speed, the acceleration and the DEM error of each PS point by using the compensated phase.

Description

Technical Field

[0001] This invention belongs to the field of bistatic synthetic aperture radar technology, specifically relating to a PS deformation estimation method for bistatic InSAR of navigation satellites based on temporal coherence. Background Technology

[0002] Bistatic InSAR (GNSS-based InBSAR) is a typical application of bistatic synthetic aperture radar, possessing the capability to monitor bistatic deformation in both time and space. It uses a navigation satellite as the transmitter, with the receiver stationary on the ground, collecting both the direct wave signal from the navigation satellite and the reflected signal from the monitored area. This not only achieves technical coupling between navigation and radar at the signal level but also integrates the advantages of differential GPS and InSAR technologies, offering benefits such as high data rate, low cost, and ease of deployment.

[0003] Currently, for the study of time-series SAR imagery, PS-InSAR technology can be used for GNSS-based InBSAR. Traditional InSAR technology, when processing complex terrain and long-term data series, is often affected by various factors such as atmospheric delay, orbital errors, and phase noise, leading to a decrease in measurement accuracy. Meanwhile, the large amount of data acquired by bistatic InSAR from navigation satellites requires effective processing and analysis methods to extract representative deformation features, providing a foundation for subsequent data analysis and interpretation, and contributing to a better understanding of the spatiotemporal evolution of surface deformation. Summary of the Invention

[0004] In view of this, the present invention proposes a PS deformation estimation method based on temporal coherence bistatic InSAR for navigation satellites. This method introduces an overall coherence coefficient based on the phase error by modeling the line-of-sight interferometric phase, and can iteratively solve the deformation rate of the PS point.

[0005] The technical solution for implementing the present invention is as follows:

[0006] A method for estimating PS deformation based on bistatic InSAR of navigation satellites with temporal coherence is described below:

[0007] Step 1: Based on the deformation phase, atmospheric phase, DEM error phase, spatial baseline phase, and noise phase of the PS point, perform PS point interferometric phase modeling, and matrix-ize the interferometric phase model;

[0008] Step 2: For the dual-ground InSAR images of navigation satellites, combined with optical images, and based on the amplitude difference method, select multiple PS points with large amplitudes that have existed for more than a set number of days, and obtain the interferometric phase by performing synchronous differential on the phases of multiple PS points.

[0009] Step 3: Set velocity, acceleration, and DEM as loop variables, and set termination conditions for the loop iteration. When the termination conditions are not met, use the currently acquired loop variables to compensate for the interference phase, and use the interference phase to estimate and compensate for the atmospheric phase and constant phase. Use the compensated phase to estimate the deformation velocity, acceleration, and DEM error of each PS point.

[0010] Optionally, the interference phase model described in this invention can be represented in matrix form as follows:

[0011] Φ=Φ01+ηR+Bh+Tv+T 2 a+E

[0012] Where η is the atmospheric phase scaling factor, R is the target receiver slant range, B represents the spatial baseline, h represents the DEM error, T represents the adjacent time baseline, v represents the average deformation rate, a represents the average deformation acceleration, E represents the error phase, and Φ0 represents the constant phase term.

[0013] Optionally, the method of compensating for the interference phase using the currently acquired cyclic variable described in this invention is as follows:

[0014] Φ (n) =Φ-(Bh (n) +Tv (n) +T 2 a (n) )

[0015] Wherein, the superscript n is the cyclic variable, B represents the spatial baseline, h represents the DEM, T represents the adjacent time baseline, v represents the average deformation rate, a represents the average deformation acceleration, and Φ is the interference phase obtained in step two.

[0016] Optionally, the method of estimating atmospheric phase and constant phase using interferometric phase in this invention specifically includes:

[0017] For the interference phase Φ (n) For each row k, the atmospheric phase is estimated using the periodogram method. and constant phase Once all rows are estimated, the atmospheric phase matrix η can be obtained. (n) and constant phase matrix

[0018]

[0019] In this diagram, the superscript 'm' represents the transmitter satellite serial number, the subscript 'l' represents the PS point serial number, and 'L' represents the number of PS points. This represents the atmospheric phase scaling factor for the m-th satellite, indicating the k-th reorbiting time. The interferometric phase of the l-th PS point of the m-th satellite represents the reorbiting time of the n-th iteration. This represents the slant range of the receiver at the l-th PS point of the m-th satellite.

[0020] Optionally, the compensated atmospheric phase Φ in step three of the present invention (n-NoAPS) for:

[0021]

[0022] Where R is the slant range of the target receiver.

[0023] Optionally, for each PS point, i.e. Φ, the present invention... (n-NoAPS) For each column, weighted based on the absolute value of the coherence coefficient, the deformation rate, acceleration, and DEM error are estimated as follows:

[0024]

[0025]

[0026] in, Φ (n-NoAPS) The phase value in the k-th row and m*l-th column. This represents the time of the k-th reorbiting and the spatial baseline of the m-th satellite. and Let represent the DEM error, average deformation rate, and average deformation acceleration of the l-th PS point of the m-th satellite, respectively, with δ being a limit value.

[0027] Optionally, the iteration termination condition described in this invention is:

[0028] ||δh (n-1) ||<T h ;||δv (n-1) ||<T v ;||δa (n-1) ||<T a or n =N MAX

[0029] Among them, T h T v And T a N represents the set termination convergence condition. MAX This indicates the maximum number of loops set.

[0030] Optionally, in the initial iteration of this invention, both the loop variable and the loop variable error are set to 0.

[0031] Optionally, the present invention performs window sliding processing on the data, setting the sliding window to 7, and using 7 days of data to complete the deformation estimation and error compensation for the intermediate day.

[0032] Optionally, the present invention sets a maximum number of iterations N. MAX =200.

[0033] Beneficial effects:

[0034] This method models the line-of-sight interferometric phase, compensates for non-ideal phase errors, introduces the overall coherence coefficient, iteratively estimates the deformation velocity and deformation acceleration at the PS point, and then uses the time baseline to obtain the deformation. This provides a prerequisite for high-precision deformation acquisition and plays an important role in the practical application of GNSS-based InBSAR. Attached Figure Description

[0035] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0036] Figure 1 This is a GNSS-based InBSAR topology configuration based on a planar scene.

[0037] Figure 2 This is the verification experimental equipment used in the embodiments of the present invention.

[0038] Figure 3 The image shows the BP imaging results of the embodiment of the present invention.

[0039] Figure 4 This describes the changes in the iterative variables during the iterative estimation process.

[0040] Figure 5 The results of the iterative estimation are: deformation velocity (a) and deformation acceleration (b). Detailed Implementation

[0041] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0042] It should be noted that, in the absence of conflict, the following embodiments and features can be combined with each other; and, based on the embodiments of this disclosure, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this disclosure.

[0043] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0044] Bistatic InSAR based on navigation satellites uses an in-orbit navigation satellite as the transmitter and a receiver fixed on the ground. An omnidirectional antenna receives the direct wave signal from the navigation satellite, achieving time-frequency synchronization, while a directional antenna receives reflected signals from the monitoring area, enabling target scene data acquisition. This not only achieves technical coupling between navigation and radar at the signal level but also integrates the advantages of differential GPS and InSAR technologies. The overall system configuration is as follows: Figure 1 As shown.

[0045] SAR images based on GNSS systems are typically considered to be composed of discrete resolution units, known as PS (Perimeter Points). Deformation monitoring is performed using PS-InSAR technology. The phase of the PS point Q can be modeled as:

[0046]

[0047] Among them, P T (k), P R (k) and P Q (k) represents the positions of the satellite, receiver, and target, respectively, and λ is the carrier wavelength. The phase error introduced by noise, etc., P E =P R The position of the receiver is indicated by k, the time is indicated by k+1, and the time of the next re-railing is indicated by k+1.

[0048] Based on the high coherence between adjacent secondary reorbital SAR images, the interferometric phase model for point Q at PS is:

[0049]

[0050] Where φ0 represents the constant phase term, which is independent of the PS point and is equivalent to the mean of the phase of the entire interference pattern. For deformation phase (including deformation rate and deformation acceleration), Atmospheric phase, Indicates the DEM error phase. For spatial baseline phase, Indicates the noise phase.

[0051] Atmospheric phase, introduced by two echoes of the same target at different times, can be expressed as:

[0052]

[0053] Where, n atm Let Δn represent the atmospheric refractive index at the current time t and the current location, and R represent the bistatic distance difference. However, in actual data processing, Δn... atm It can be set to a fixed value.

[0054] Compared to SAR satellites, GNSS lacks strict trajectory control. At each reorbiting event, a trajectory deviation exists relative to the previous reorbiting, known as the space baseline, which can cause a space baseline phase. The interferometric phase φ introduced from the space baseline... B It can be represented as:

[0055]

[0056] Where R represents the slant range of the transmitter satellite, Ψ represents the projection matrix of the unit vector of the transmitter satellite's line-of-sight direction, and B... ⊥ Indicates the size of the projection of the spatial baseline in the celestial direction.

[0057] The baseline can be approximated as:

[0058]

[0059] Among them, R T0 B represents the slant distance from the center of the satellite synthetic aperture to the center of the scene. ⊥ This represents the size of the spatial baseline projected upwards. Therefore, by compensating for non-ideal phases, the accuracy of deformation phase extraction at a single angle can be effectively improved.

[0060] By matrixing the interferometric phase model, we obtain:

[0061] Φ=Φ01+ηR+Bh+Tv+T 2 a+E(6)

[0062] Where η is the atmospheric phase scaling factor, R is the target receiver slant range, B represents the spatial baseline, h represents the DEM error, T represents the adjacent time baseline, v represents the average deformation rate, a represents the average deformation acceleration, E represents the error phase, and Φ0 represents the constant phase term.

[0063]

[0064] Based on this, the overall coherence coefficient is introduced, which allows for iterative solutions to atmospheric phase, DEM error, deformation rate, and acceleration. The specific process is as follows:

[0065] Step 1: Based on the deformation phase, atmospheric phase, DEM error phase, spatial baseline phase, and noise phase of the PS point, perform PS point interferometric phase modeling, and matrix-ize the interferometric phase model;

[0066] Step 2: For the dual-ground InSAR images of navigation satellites, combined with optical images, and based on the amplitude difference method, select multiple PS points with large amplitudes that have existed for more than a set number of days, and obtain the interferometric phase by performing synchronous differential on the phases of multiple PS points.

[0067] Step 3: Set velocity, acceleration, and DEM as loop variables, and set termination conditions for the loop iteration. When the termination conditions are not met, use the currently acquired loop variables to compensate for the interference phase, and use the interference phase to estimate and compensate for the atmospheric phase and constant phase. Use the compensated phase to estimate the deformation velocity, acceleration, and DEM error of each PS point.

[0068] The third step of the present invention will be described in detail below.

[0069] (1) The acquired interference phase Φ is taken as the initial phase Φ (0) Let the loop variable be n, starting from n = 0:

[0070]

[0071] Where, δv (0) δa (0) and δh (0) This represents the velocity, acceleration, and DEM error when the loop variable n = 0. (0) a (0) and h (0) The values ​​represent velocity, acceleration, and DEM when the cyclic variable n = 0; the interferometric phase Φ is a K*ML matrix; K represents the number of reorbital times, L represents the number of PS points, and M represents the number of transmitter satellites. This can be represented as follows:

[0072]

[0073] (2) Repeat steps (a)-(f) until convergence:

[0074] (a) Update the loop variable and calculate velocity, acceleration, and DEM:

[0075] n = n + 1

[0076] h (n) =h (n-1) +δh (n-1)

[0077] v (n) =v (n-1) +δv (n-1)

[0078] a (n) =a (n-1) +δa (n-1) (8)

[0079] (b) When n≥2, the iteration termination condition is used for judgment. If the termination condition is met, the following step (c) is not entered. The termination condition is:

[0080] ||δh (n-1) ||<T h ;||δv (n-1) ||<T v ;||δa (n-1) ||<T a or n = N MAX (9)

[0081] Among them, T h T v And T a N represents the set termination convergence condition. MAX This indicates the maximum number of loops set.

[0082] (c) Compensation for interference phase introduced by DEM errors and deformation errors:

[0083] Φ (n) =Φ-(Bh (n) +Tv (n) +T 2 a (n) (10)

[0084] (d) Regarding the interference phase Φ introduced by the compensation DEM error and deformation error (n) For each row k (each row corresponds to the interferometric phase obtained from one re-regression), the atmospheric phase is estimated using the periodogram method. and constant phase Obtain the atmospheric phase matrix η (n) and constant phase matrix

[0085]

[0086] In this diagram, the superscript 'm' represents the transmitter satellite serial number, the subscript 'l' represents the PS point serial number, and 'L' represents the number of PS points. This represents the atmospheric phase scaling factor for the m-th satellite, indicating the k-th reorbiting time. The interferometric phase of the l-th PS point of the m-th satellite represents the reorbiting time of the n-th iteration. This represents the slant range of the receiver at the l-th PS point of the m-th satellite.

[0087] (e) Utilizing the atmospheric phase matrix η (n) and constant phase matrix Compensation for atmospheric phase and constant phase includes:

[0088]

[0089] Where R is the slant range of the target receiver.

[0090] (f) For each PS point, i.e. Φ (n-NoAPS) For each column, weighted based on the absolute value of the coherence coefficient, the deformation rate, acceleration, and DEM error are estimated as follows:

[0091]

[0092]

[0093] in, Φ (n-NoAPS) The phase value in the k-th row and m*l-th column. This represents the time of the k-th reorbiting and the spatial baseline of the m-th satellite. and These represent the DEM error, average deformation rate, and average deformation acceleration of the l-th PS point of the m-th satellite, respectively. δ is just a constraint value and is a known input quantity.

[0094] At this point, the deformation velocity and acceleration at point PS can be obtained. To increase the reliability of the method, the data is processed by a sliding window by inputting more data. The window size can be determined based on the spatial baseline caused by the satellite trajectory. To maximize the use of time correlation while avoiding the introduction of excessive spatial baseline phase error, the sliding window is generally set to 7 days. This allows for deformation estimation and error compensation for the intermediate days using 7 days of data, thereby improving the accuracy of the interferometric phase.

[0095] The specific implementation examples are as follows:

[0096] Step 1: Experimental Design and Data Acquisition

[0097] This embodiment was conducted in Chongqing, China (30.88892N, 109.47449E). Using BDS-2IGSO1, BDS-2IGSO3, BDS-2IGSO6, BDS-2IGSO7, and BDS-3IGSO1 as transmitters, a receiving system was established to collect raw data for 10 days, from May 4th to May 13th, 2024. The experimental equipment in this embodiment is as follows... Figure 2 As shown, in bistatic InSAR of navigation satellites, accurate acquisition of deformation phase can help predict deformation in the monitored area, providing data reference for geological disaster early warning.

[0098] Step 2: Target Phase Extraction

[0099] In this embodiment, the raw data is processed using the BP imaging method to obtain SAR image results, which are then combined with optical images. Figure 3 As shown in the left figure. Then, based on the amplitude deviation method, 50 PS points with large amplitudes and existing for more than 7 days were selected as the basis. Figure 3 As shown in the right figure, the PS points are mostly located on both sides of the ravine, mainly distributed between -600m and 300m to the east, and between 0m and 200m to the north, allowing for accurate and complete monitoring of the target area. The interferometric phase results are obtained by synchronously differentiating the phase of the PS points.

[0100] Step 3: Iterative solution of the deformation velocity at point PS

[0101] The original data from the example is selected over a sliding period of 7 days for deformation rate and acceleration estimation. The result is used as the deformation rate for the intermediate day. A maximum number of iterations N is set. MAX =200, the number of PS is 50.

[0102] Therefore, this invention uses a temporally coherent PS deformation estimation method to estimate the deformation rate of PS points, providing a guarantee for subsequent high-precision deformation inversion and laying the foundation for the practical application of bistatic InSAR based on navigation satellites.

[0103] Figure 4 The changes in iterative variables during the estimation process are shown. (a) shows the spatial baseline phase result. It can be seen that the spatial baseline phase increases with the increase of the number of re-orbiting events, indicating that the influence of the spatial baseline becomes increasingly significant over time. (b) shows the atmospheric phase coefficient. During the iteration process, the atmospheric phase coefficient did not converge. (c) shows the final atmospheric phase estimation result. (d) shows the image coherence coefficient estimation result, which is used for weighted processing.

[0104] The final estimation results are as follows Figure 5As shown, the left side represents the velocity estimation results, and the right side represents the acceleration estimation results. For the deformation results of the PS point each day, the mean square deviation is calculated in the spatial dimension as the deformation accuracy. The deformation accuracy results for a single angle are shown in the table below. The best accuracy for a single star on a single day can reach 16.14 mm, and the average deformation accuracy for 5 stars is 19.48 mm.

[0105] Table 1 Single Angle Deformation Accuracy

[0106]

[0107] 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 method for estimating PS deformation based on temporal coherence bistatic InSAR of navigation satellites, characterized in that, The specific process is as follows: Step 1: Based on the deformation phase, atmospheric phase, DEM error phase, spatial baseline phase, and noise phase of the PS point, perform PS point interferometric phase modeling, and matrix-ize the interferometric phase model; Step 2: For the dual-ground InSAR images of navigation satellites, combined with optical images, and based on the amplitude difference method, select multiple PS points with large amplitudes that have existed for more than a set number of days, and obtain the interferometric phase by performing synchronous differential on the phases of multiple PS points. Step 3: Set velocity, acceleration, and DEM as loop variables, and set the termination condition for terminating the loop iteration; When the termination condition is not met, the interferometric phase is compensated using the currently acquired cyclic variables, and the atmospheric phase and constant phase are estimated and compensated using the interferometric phase. The deformation rate, acceleration and DEM error of each PS point are estimated using the compensated phase.

2. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 1, characterized in that, The matrix representation of the interference phase model is as follows: Φ=Φ01+ηR+Bh+Tv+T 2 a+E Where η is the atmospheric phase scaling factor, R is the target receiver slant range, B represents the spatial baseline, h represents the DEM error, T represents the adjacent time baseline, v represents the average deformation rate, a represents the average deformation acceleration, E represents the error phase, and Φ0 represents the constant phase term.

3. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 1, characterized in that, The method of compensating for the interference phase using the currently acquired cyclic variable is as follows: F (n) =Φ-(Bh (n) +Tv (n) +T 2 a (n) ) Wherein, the superscript n is the cyclic variable, B represents the spatial baseline, h represents the DEM, T represents the adjacent time baseline, v represents the average deformation rate, a represents the average deformation acceleration, and Φ is the interference phase obtained in step two.

4. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 3, characterized in that, The method of estimating atmospheric phase and constant phase using interferometric phase estimation is as follows: For the interference phase Φ (n) For each row k, the atmospheric phase is estimated using the periodogram method. and constant phase Once all rows are estimated, the atmospheric phase matrix η can be obtained. (n) and constant phase matrix In this diagram, the superscript 'm' represents the transmitter satellite serial number, the subscript 'l' represents the PS point serial number, and 'L' represents the number of PS points. This represents the atmospheric phase scaling factor for the m-th satellite, indicating the k-th reorbiting time. The interferometric phase of the l-th PS point of the m-th satellite represents the reorbiting time of the n-th iteration. This represents the slant range of the receiver at the l-th PS point of the m-th satellite.

5. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence according to claim 4, characterized in that, The compensated atmospheric phase Φ in step three (n-NoAPS) for: Where R is the slant range of the target receiver.

6. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 5, characterized in that, For each PS point, i.e. Φ (n-NoAPS) For each column, weighted based on the absolute value of the coherence coefficient, the deformation rate, acceleration, and DEM error are estimated as follows: in, Φ (n-NoAPS) The phase value in the k-th row and m*l-th column. This represents the time of the k-th reorbiting and the spatial baseline of the m-th satellite. and Let represent the DEM error, average deformation rate, and average deformation acceleration of the l-th PS point of the m-th satellite, respectively, with δ being a limit value.

7. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence according to claim 1, characterized in that, The iteration termination condition is: ||δh (n-1) ||<T h ;||δv (n-1) ||<T v ;||δa (n-1) ||<T a or n=N MAX Among them, T h T v And T a N represents the set termination convergence condition. MAX This indicates the maximum number of loops set.

8. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 7, characterized in that, The maximum number of iterations N is set. MAX =200.

9. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence as described in claim 1, characterized in that, In the initial iteration, both the loop variable and the loop variable error are set to 0.

10. The PS deformation estimation method for bistatic InSAR navigation satellites based on temporal coherence according to claim 1, characterized in that, The data is processed by a sliding window, with the sliding window set to 7 days. The deformation estimation and error compensation for the middle day are completed using the data from the 7 days.

Images