Satellite-ground bistatic synthetic aperture radar (SAR) long-time sequence processing method

Abstract

The invention belongs to the technical field of synthetic aperture radars, and particularly relates to a satellite-ground bistatic SAR long-time sequence processing method. The method specifically comprises the following steps of: processing an interferogram: acquiring N satellite-ground bistatic SAR heavy-orbit interferograms, and performing error compensation to obtain a corrected heavy-orbit interferogram; based on the characteristic that targets in the same elevation direction have the same Doppler and the same bistatic distance, geocoding is carried out on the corrected heavy rail interferogram; performing mask processing on the geocoded interferogram, and performing phase filtering after the mask processing; pS candidate point screening and processing: screening the interferogram after phase filtering by using an amplitude deviation and coherence dual-threshold method to obtain PS candidate points; obtaining related data of the PS point and performing time sequence processing to obtain a terrain phase; and performing long-time sequence processing: calculating a deformation quantity by using the terrain phase, and calculating an average settlement rate based on the deformation quantity.

Description

Technical Field

[0001] The present invention belongs to the technical field of synthetic aperture radar, and in particular relates to a long time series processing method for satellite-ground bistatic SAR. Background Art

[0002] Deformation monitoring is a key application area for spaceborne synthetic aperture radar (SAR). By monitoring surface deformation, it is possible to provide pre-disaster warnings and damage assessments for geological hazards such as landslides, subsidence, and collapses. Traditional InSAR technology is affected by many factors when handling surface deformation monitoring, such as phase noise and decoherence caused by atmospheric disturbances, changes in surface cover (such as vegetation growth and seasonal variations), and dramatic terrain changes, as well as phase unwrapping in complex terrain conditions. These factors can prevent InSAR technology from providing accurate and reliable surface deformation information in certain situations. Time-series InSAR technology can effectively overcome the effects of spatiotemporal decoherence and separate the atmospheric delay phase from the deformation phase, thereby improving the accuracy of deformation measurements.

[0003] Most remote sensing satellites currently in orbit are polar-orbiting satellites. When performing D-InSAR, these satellites can only obtain line-of-sight deformation from ascending and descending orbits. Furthermore, the north-south components of these two angles are very small, resulting in low accuracy in the inverted north-south deformation. To address the limitations of spaceborne SAR in deformation inversion, a space-surface bistatic SAR (SS-BSAR) system can be utilized. This system uses a remote sensing satellite as an illumination source and a ground-based receiving system to receive scene echoes. Furthermore, a space-surface bistatic SAR system can utilize a single illumination source and deploy multiple receiving stations on the ground to form a multi-angle receiving system. This system can acquire multiple sets of scene data from different line-of-sight directions during a single satellite reorbit cycle, improving observation data acquisition efficiency.

[0004] Current research on satellite-ground bistatic SAR imaging and interferometry processing primarily focuses on the transmit antenna mainlobe signal, focusing on imaging, elevation inversion, and deformation inversion. There is a lack of research on using satellite-ground bistatic SAR for long-term time series processing or on using satellite-ground bistatic SAR data to obtain annual average sedimentation rates. Summary of the Invention

[0005] To solve the above problems, the present invention provides a long time series processing method for satellite-ground bistatic SAR, which can obtain the annual average sedimentation rate of the detection scene using satellite-ground bistatic SAR data.

[0006] The technical solutions for implementing the present invention are as follows:

[0007] A long time series processing method for satellite-ground bistatic SAR, the specific process is as follows:

[0008] Interferogram processing: N satellite-to-ground bistatic SAR re-orbit interferograms are acquired and error-compensated to obtain the corrected re-orbit interferograms. Based on the characteristic that targets in the same elevation direction have the same Doppler and bistatic distance, the corrected re-orbit interferograms are geocoded. Masking is performed on the geocoded interferograms, and phase filtering is performed on the masked interferograms.

[0009] PS candidate point screening and processing: The interferogram after phase filtering is screened using the amplitude deviation and coherence double threshold to obtain PS candidate points; the relevant data of the PS point is obtained and time series processed to obtain the terrain phase;

[0010] Long time series processing: using the terrain phase to calculate the deformation, and calculating the average sedimentation rate based on the deformation.

[0011] Optionally, the specific process of obtaining the corrected heavy track interferogram by error compensation according to the present invention is as follows:

[0012] Using an interferometric error fringe estimation method based on precise local phase estimation, we perform frequency domain transformation on the interferogram data along the range and azimuth directions, extracting and eliminating linear trends in the frequency domain to compensate for interferometric phase errors. This method simultaneously eliminates re-orbit interferometric phase errors caused by along-track error, cross-track error, and ephemeris matching errors.

[0013] Optionally, the specific process of the error compensation of the present invention to obtain the corrected heavy track interferogram is as follows: int (n) performing a two-dimensional Fourier transform to estimate the linear frequency range in the azimuth / range direction, and setting the initial sampling frequency θ0 and the sampling frequency interval φ0 according to the linear frequency range; performing a Chirp-Z transform on the heavy track interferogram, and taking the frequency corresponding to the transformed peak as the frequency error estimate; finally, using the frequency error estimate result to remove the phase error at the heavy track interferogram to obtain a corrected interferogram.

[0014] Optionally, the present invention sets the linear frequency range to [f min ,f max ],

[0015]

[0016] Among them, N c is the number of sampling intervals;

[0017] The frequency corresponding to the peak after transformation is used as the frequency error estimate:

[0018]

[0019] Among them, Sint,r (n) with S int,a (n) represents the Chirp-Z transform of one-dimensional re-orbit interferometry data along the range and azimuth directions, respectively;

[0020] The corrected interference pattern is s int_c ,

[0021]

[0022] Among them, s int is the uncorrected re-orbit interferometry data, y and x are column vectors containing the azimuth and range sequence numbers of the interferogram, respectively, and 1 is a column vector whose elements are all 1.

[0023] Optionally, the specific process of geocoding the corrected heavy track interferogram according to the present invention is as follows:

[0024] An external digital elevation model (DEM) is used to project each DEM point onto the imaging plane along the elevation direction, and a mapping relationship is established between the longitude and latitude coordinate system and each coordinate point on the imaging plane. Based on the above mapping relationship, the heavy track interferogram is reversely projected into the longitude and latitude coordinate system to obtain the geocoded interferogram.

[0025] Optionally, the present invention performs mask processing on the geocoded interference map and performs phase filtering on the masked image. The specific process is as follows:

[0026] The digital elevation model is used to calculate the elevation and azimuth angles of each target point in the receiver's observation scene. At the same elevation and azimuth angles, the target point closest to the receiver is set as the visible point, and the remaining target points are set as blocked points.

[0027] The interference pattern is masked, that is, the interference phase of the invisible area is set to 0, and the interference phase after masking is phase filtered.

[0028] Optionally, the present invention sets the target point closest to the receiver as the visible point and the remaining target points as the blocked points, which is specifically expressed as follows:

[0029]

[0030] Among them, Flag(k) is the discrimination result of the kth DEM sampling point, Flag = 0 corresponds to an invisible point, Flag = 1 corresponds to a visible point, σ(k) = [δ(k); η(k)], δ(k) is the pitch angle of the kth DEM sampling point observed by the receiver, η(k) is the azimuth angle of the kth DEM sampling point observed by the receiver, R R (k) is the slant distance from the receiver to the kth DEM sampling point.

[0031] Optionally, the present invention utilizes the terrain phase Calculate deformation for:

[0032]

[0033] Where β is the bistatic angle and λ is the radar wavelength.

[0034] Optionally, the annual average sedimentation rate v of the present invention is:

[0035]

[0036] Among them, G is the design function and W is the weight function.

[0037] Optionally, the design function G and weight function W of the present invention are respectively:

[0038]

[0039] Among them, t1~t n represents the time information of the observation data, σ ij represents the covariance between the i-th and j-th observations.

[0040] Optionally, the present invention selects one-dimensional re-orbit interferograms in a flat area from the N acquired satellite-ground bistatic SAR re-orbit interferograms and performs error compensation to obtain a corrected re-orbit interferogram.

[0041] Beneficial effects:

[0042] First, the present invention utilizes an interference error fringe estimation method based on local phase precise estimation. By performing frequency domain transformation on the interference pattern data along the range and azimuth directions, the linear trend in the frequency domain is extracted and eliminated, thereby ultimately compensating for the error caused by the interference phase.

[0043] Second, since the interference result corresponding to the blocked area has a low signal-to-noise ratio and poor interference quality, which will affect the subsequent phase unwrapping processing, the interference pattern after geocoding is masked and phase filtered to reduce phase noise and smooth the interference phase. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0045] Figure 1 This is the elevation inversion geometry diagram of the satellite-ground bistatic SAR;

[0046] Figure 2 Geometry diagram of satellite-ground bistatic SAR re-orbit interferometry;

[0047] Figure 3 This is a schematic diagram of the satellite-ground bistatic SAR system configuration;

[0048] Figure 4 This is a schematic diagram of single-angle re-orbit interferometry of satellite-ground bistatic SAR;

[0049] Figure 5 is a geocoded, image-registered, and phase-smoothed differential interferogram;

[0050] Figure 6 The distribution map of PS candidate points screened by double threshold method;

[0051] Figure 7 is the winding phase diagram after combined filtering;

[0052] Figure 8 is the distribution map of PS candidate points after secondary screening;

[0053] Figure 9 is the unwrapped phase image with the error removed;

[0054] Figure 10 is the annual average sedimentation rate map in the study area;

[0055] Figure 11 is the histogram of the annual average sedimentation rate. DETAILED DESCRIPTION

[0056] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

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

[0058] It should be noted that various aspects of the embodiments within the scope of the appended claims are described below. It should be apparent that the aspects described herein can be embodied in a wide variety of forms, and any specific structure and / or function described herein is merely illustrative. Based on this disclosure, it should be understood by those skilled in the art that an 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 described herein can be used to implement an apparatus and / or practice a method. In addition, other structures and / or functionalities other than one or more of the aspects described herein can be used to implement this apparatus and / or practice this method.

[0059] This embodiment of the application provides a long time series processing method for satellite-ground bistatic SAR, the specific process of which is as follows:

[0060] Interferogram processing: N satellite-to-ground bistatic SAR re-orbit interferograms are acquired and error-compensated to obtain a corrected re-orbit interferogram. An elevation-direction expression is established, assuming that targets in the same elevation direction have the same Doppler and bistatic distance. The corrected re-orbit interferogram is geocoded using the elevation-direction expression. The geocoded interferogram is masked and then phase filtered.

[0061] PS candidate point screening and processing: The interferogram after phase filtering is screened using the amplitude deviation and coherence double threshold to obtain PS candidate points; the relevant data of the PS point is obtained and time series processed to obtain the terrain phase;

[0062] Long time series processing: using the terrain phase to calculate the deformation, and calculating the average sedimentation rate based on the deformation.

[0063] The above process of the present invention is described in detail below: interference pattern processing (corresponding to steps S1 to S4), PS candidate point screening and processing (corresponding to steps S5 to S8) and long time series processing (corresponding to step S9).

[0064] S1. Obtain N main lobe re-orbit interferograms of the satellite-ground bistatic SAR;

[0065] S2. Acquisition of heavy track interferogram: Error compensation is performed on the N heavy track interferograms obtained in S1 to obtain the corrected heavy track interferogram s int_c ; The specific process of this step is:

[0066] This step uses the interference error fringe estimation method based on local phase precision estimation. By performing frequency domain transformation on the interference pattern data along the range and azimuth directions, the linear trend in the frequency domain is extracted and eliminated, and finally the phase error caused by the interference is compensated. Specifically:

[0067] This method will simultaneously eliminate the re-orbit interferometry phase error caused by along-track error, cross-track error, and ephemeris matching error. The principle of high frequency resolution of Chirp-Z transform is as follows:

[0068] Assume the one-dimensional time domain signal is:

[0069] s int (n),n=0,...,N-1 (1)

[0070] Among them, s int (n) is the one-dimensional re-orbit interferometric data with interference error (along the range or azimuth direction), and N is the number of sampling points of the one-dimensional data. Then sint (n) The corresponding Chirp-Z transform expression is:

[0071]

[0072] Among them, M is the number of Chirp-Z transform sampling points, z k It is called the Z transform sampling point, and its expression is:

[0073]

[0074] in, represents the starting complex sampling point, θ0 is the initial sampling frequency, is the complex sampling interval, and φ0 is the frequency sampling interval. The frequency sampling interval φ0 can be set to improve the resolution of the corresponding spectrum.

[0075] The interference data corresponding to the relatively flat terrain in the experimental scene are selected to avoid the influence of terrain phase on fringe estimation. int (n) Perform a two-dimensional Fourier transform to roughly estimate the range of the linear frequency in the azimuth / range direction, denoted as [f min ,f max ]. According to the estimation results of the linear frequency range [f min ,f max ]Set the initial sampling frequency θ0 and the sampling frequency interval φ0:

[0076]

[0077] Among them, N c is the number of sampling intervals.

[0078] Next, the Chirp-Z transform is performed on the one-dimensional heavy track interferometry data in the flat area, and the frequency corresponding to the peak after the transformation is taken as the frequency error estimation result, which can be written as:

[0079]

[0080] Among them, S int,r (n) with S int,a (n) represents the Chirp-Z transform of the one-dimensional re-orbit interferometry data along the range and azimuth directions, respectively.

[0081] Finally, the frequency error estimation result is used to remove the phase error in the re-orbit interferometry data to obtain the corrected interferogram s int_c , which can be expressed as:

[0082]

[0083] Among them, s intis the uncorrected re-orbit interferometry data, y and x are column vectors containing the azimuth and range sequence numbers of the interferogram, respectively, and 1 is a column vector whose elements are all 1.

[0084] S3. Geocoding image registration of heavy track interferogram: Using the same elevation direction target with the same Doppler and the same bistatic distance, establish the elevation direction expression, and use the elevation direction expression to geocode the corrected heavy track interferogram; after obtaining the absolute phase of the scene terrain, it is necessary to establish h S The corresponding relationship between the scene terrain phase and the scene terrain; the specific steps include:

[0085] S31. Geocode the terrain phase and project it into a latitude and longitude coordinate system.

[0086] The geometric diagram of the satellite-ground bistatic SAR elevation inversion is shown in Figure 1. Assume that S is the satellite position, A e is the position of the receiver, S E is the position of the equivalent phase center of the SS-BSAR system, is the altitude of the imaging plane, T is the position of a target in the scene, and T0 is the position of the target T projected on the imaging plane along the elevation direction. In the SS-BSAR system, the range gradient With vector The same direction is the direction of the bisector of the transmitting and receiving line of sight angles, which can be expressed as:

[0087]

[0088] Among them, n ST and n AT Respectively represent and vector (i.e., the position vector from the satellite position to a target in the scene) and vector (i.e., the position vector from the receiver to the position vector of a target in the scene) is a unit vector in the same direction. In the SS-BSAR system, the Doppler gradient The effective velocity vector is in the same direction as the transmitter, which can be expressed as:

[0089]

[0090] Where I is the unit matrix. In SS-BSAR imaging, targets in the same elevation direction have the same Doppler frequency and the same bistatic distance, that is, T0 and T have the same Doppler frequency and bistatic distance. Let the elevation direction be Based on the above theory, it can be expressed as:

[0091]

[0092] The elevation direction expression shows that the elevation direction in the SS-BSAR scene changes with the change of the spatial position of the target point, that is, the elevation direction in the SS-BSAR scene is space-variant.

[0093] S32. Since the elevation direction expression (9) shows that the elevation direction in the SS-BSAR scene is space-varying, the heavy track interferogram can be geocoded in the following way.

[0094] Specifically, an external digital elevation model (DEM) is used to project each DEM point onto the imaging plane along the elevation direction, establishing a mapping relationship between the latitude and longitude coordinate system and the coordinate points on the imaging plane. Based on this mapping relationship, the heavy track interferogram is then back-projected into the latitude and longitude coordinate system to obtain a geocoded interferogram.

[0095] S4. Masking the geocoded interferogram and performing phase filtering on the masked interferogram. The specific steps include:

[0096] S41. Using an external DEM, calculate the elevation and azimuth angles of each target point in the receiver's observation scene. Then, at the same elevation / azimuth angle, set the target point closest to the receiver as the visible point, and the remaining target points as the blocked points, that is:

[0097]

[0098] Among them, Flag(k) is the discrimination result of the kth DEM sampling point, Flag = 0 corresponds to an invisible point, and Flag = 1 corresponds to a visible point. σ(k) = [δ(k); η(k)], δ(k) is the pitch angle of the kth DEM sampling point observed by the receiver, η(k) is the azimuth angle of the kth DEM sampling point observed by the receiver, R R (k) is the slant distance from the receiver to the kth DEM sampling point.

[0099] S42 , performing mask processing on the interference pattern, that is, setting the interference phase of the invisible area to 0. Then, performing phase filtering on the interference phase after masking to reduce phase noise and smooth the interference phase.

[0100] S5. PS candidate point screening: Use the amplitude deviation and coherence double threshold method to screen and obtain PS candidate points. The specific steps include:

[0101] S51, for the processing result of S42, for a certain resolution unit, its coherence coefficient The neighboring pixel information within a certain window containing the pixel can be selected for estimation. The formula is as follows:

[0102]

[0103] Among them, S1,n 、S 2,n Represents the first interference pattern and the second reflection pattern, n represents the pixel point, N represents the total number of pixel points, and the average value of the correlation coefficient of each pixel in the time series interference pattern is obtained, that is, the time series coherence coefficient Then, by setting a suitable threshold T γ ,Filter PS candidate points based on the coherence coefficient threshold, namely:

[0104] S52. Based on the processing results of S42, select a flat area in the satellite-ground bistatic SAR re-orbit data, calculate the average amplitude of the area, and perform amplitude comparison on all data.

[0105] S53, the amplitude deviation threshold method is based on the fact that the stability of the PS point can be expressed by the statistical characteristics of the echo phase in the time series. A suitable phase deviation threshold is selected, and pixels larger than this value are considered to have relatively stable scatterers. The expression of the amplitude deviation index is:

[0106]

[0107] Among them, σ A is the standard deviation of the time series amplitude, m A is the amplitude mean. The average value of the amplitude deviation of each pixel in the time series interference pattern is obtained, that is, the time series amplitude deviation value D A , and then by setting a suitable threshold T D ,Filter PS candidate points based on amplitude deviation threshold, namely: D A ≤T D .

[0108] S54, according to S51 and S53, the intersection of candidate points is obtained to obtain the PS candidate points of amplitude deviation and coherence double threshold method, such as Figure 6 shown.

[0109] S6. Calculate the data matrix including the time-space baseline, phase value, bistatic angle, observation angle, elevation value, amplitude deviation value, longitude and latitude values of the selected PS candidate point.

[0110] according to Figure 2 Geometry diagram of the satellite-ground bistatic SAR re-orbit interferometry, where S1 and S2 are the positions of the synthetic aperture center at the two satellite passes, respectively. ⊥ is the vertical baseline of the re-orbit interferometry. In the re-orbit SS-BSAR image after the ground is leveled, the target T phase can be expressed as:

[0111]

[0112] Furthermore, the phase after the interference between track S1 and track S2 can be expressed as:

[0113]

[0114] in, and and The unit vector in the same direction, R is approximately the slant distance from the transmitter to the target elevation projection point T0, that is, The modulus of θ is the angle between the target elevation and the horizontal plane, and α is the angle between the vertical baseline and the target elevation. From (14), we can see that the SS-BSAR re-orbit interferometric phase and the vertical baseline b of the re-orbit interferometer are ⊥ It is related to the slant distance R from the transmitter to the target elevation projection point.

[0115] SS-BSAR heavy track vertical baseline b ⊥ It cannot be directly calculated through ephemeris. To estimate the vertical baseline corresponding to the re-orbit interferometry, an area with a high correlation coefficient is selected, and the phase unwrapped re-orbit interferometry phase of this area (after geocoding) is obtained. The re-orbit baseline is estimated in combination with the external DEM, which can be expressed as:

[0116]

[0117] in, is the phase difference of the re-orbit interferometry between the selected reference points, and Δh is the difference in the DEM values of the selected reference points in the area. To improve the estimation accuracy of the vertical baseline of the re-orbit interferometry, multiple groups of reference points can be selected in this area, and the reference values of the vertical baselines of the re-orbit interferometry can be obtained based on formula (15). The reference values are then averaged to improve the estimation accuracy.

[0118] according to Figure 3 The schematic diagram of the satellite-ground bistatic SAR system configuration gives the calculation formula of the bistatic angle:

[0119]

[0120]

[0121] S7, performing time series processing on the relevant data of the PS points obtained through S1 to S6. The specific steps are as follows:

[0122] S71, phase data Combined filtering is performed, including adaptive smoothing filtering and low-pass filtering, to suppress noise and high-frequency spatial variations while retaining low-frequency phase trends.

[0123]

[0124] The combined filter is expressed as:

[0125] G=G a +G low (19)

[0126] Adaptive filter G a The composition is as follows:

[0127]

[0128] Where B is a Gaussian window used to smooth the amplitude spectrum. is the phase data. α is an exponential parameter used to adjust the weight distribution. β is a weight parameter used to control the strength of the adaptive filter. * indicates a two-bit convolution process.

[0129] Low-pass filter G low The composition is as follows:

[0130]

[0131] Among them, f i represents the one-dimensional Butterworth response frequency axis, f0 represents the cutoff frequency, and μ represents the low-pass wavelength.

[0132] S72. Estimate the terrain error phase of PS point

[0133]

[0134] Among them, K is the terrain error coefficient, B is the baseline matrix, C is the constant term, and n is the noise term.

[0135] In least squares fitting, the noise term n is usually treated as a random error, and its impact is reduced by minimizing the residual sum of squares. Specifically, the noise term is included in the residuals during the fitting process, and the values of K and C are estimated by minimizing the residual sum of squares.

[0136] S73. Calculate the phase noise standard deviation of each PS point, and remove those PS points whose phase noise standard deviation is greater than the standard deviation threshold.

[0137] Calculation of phase noise standard deviation: Triangulate the phase after removing the unstable PS points to determine the geometric connection between the PS points. For each edge, calculate the phase difference between the two endpoints.

[0138]

[0139] in, It represents the phase difference between the two end points. represents the phase after removing the unstable PS point. express The end pixel index of Indicates the starting pixel index

[0140]

[0141] in, Represents the weighted average phase, wf is a weighting factor based on the time difference. The phase difference is adjusted to subtract the weighted average

[0142]

[0143] Smoothing the phase difference in time series to obtain the smoothed phase difference

[0144]

[0145] Where G is the design function, lscov represents the weighted least squares method, and m2 represents the result of the second least squares method.

[0146] Calculate the phase angle difference between the phase difference and the smoothed phase difference to estimate the phase noise. Calculate the noise standard deviation.

[0147]

[0148]

[0149] Among them, B(ifgindex) is the spatial baseline of the interference pattern, ifg var express Variance

[0150] S74. Perform spatially independent viewing angle error correction on the twisted phase to improve the quality of the interference pattern and the accuracy of the phase.

[0151]

[0152] in, represents the phase after correction, B represents the spatial baseline, C ps and K ps Represent the constant term and slope term in the phase offset respectively

[0153] Phase unwrapping is performed on the corrected selected PS pixels using a 3D cost function phase unwrapping algorithm.

[0154] S75. Estimate the spatially correlated incidence angle error (SCLA). The calculation process is as follows:

[0155] Construct the observation matrix G

[0156]

[0157] in, is the average value of the vertical baseline, is the average of the dates.

[0158] Constructing an observation model

[0159]

[0160] in, is the unwrapped phase, m is the parameter vector, including the incident angle error coefficient K ps and the main image noise term C ps

[0161]

[0162] K ps The value of m is obtained from the first column m1 of m, C ps The value of is obtained from the second column m2 of m.

[0163] From this we can get the incident angle error phase

[0164]

[0165] The main image noise term is estimated by the least squares method

[0166]

[0167] S76. Calculate the track inclination error. The track inclination error can be expressed as a linear function in the form of:

[0168]

[0169] This can be expressed in matrix form:

[0170]

[0171] Using the least squares method to solve the above system of equations, the solution is:

[0172]

[0173] In this way, a, b, and c can be obtained, and the phase of the track inclination error can be calculated.

[0174] S77, Extract linear atmospheric delay phase error. Using the linear model method, based on a single interferogram, the tropospheric phase delay It can be estimated based on the relationship between the interference phase and the ground elevation:

[0175]

[0176] Wherein, K represents the phase-elevation linear scaling factor, which is a constant obtained by fitting the relationship between the phase and elevation h of the interferogram and is estimated based on the entire interferogram. Represents the overall shift in the phase of a single interferogram.

[0177] S8, use the unwrapped phase obtained in S74 to remove the terrain error in S72, the incidence angle error (SCLA) in S75, the orbit error in S76, and the linear atmospheric delay phase in S77 to obtain the terrain phase

[0178] S9. Calculate the linear deformation phase and obtain the average sedimentation rate. The calculation steps are as follows:

[0179] S91, deformation calculation method of satellite-ground dual-base single-angle re-orbit interferometry inversion: According to Figure 4 Schematic diagram of single-angle re-orbit interferometry of satellite-ground bistatic SAR, where S1 and S2 are the positions of the synthetic aperture center of the satellite's two illumination scenes, T1 and T2 are the target positions of the satellite's two illumination scenes (deformation and position change occurred during the re-orbit), A e is the position of the receiving antenna, β is the bistatic angle, S E is the bistatic equivalent phase center. Based on the SS-BSAR re-orbit interferometry geometry, the re-orbit interferometry phase can be expressed as:

[0180]

[0181] For the slope distance difference term in (41), Simplify:

[0182]

[0183] In addition, since the deformation is usually in the centimeter level, the satellite's re-orbit baseline is usually in the hundreds of meters level, and the slant range from the satellite to the target is usually in the hundreds of kilometers level, we can get:

[0184]

[0185] Based on formula (43), the inclusion of as well as At the same time,

[0186] Similar to the formula, Taylor expansion is performed on formula (42), and Further simplified to:

[0187]

[0188] in, is with The same unit vector. The term simplifies to:

[0189]

[0190] Substituting equations (44) and (45) into equation (41), we obtain:

[0191]

[0192] Based on formula (46), the heavy orbit interference phase can be divided into The relevant deformation components and the terrain component related to the spatial baseline Respectively expressed as:

[0193]

[0194] in, It is the unit vector along the SS-BSAR equivalent line of sight. From (47), we can see that the deformation measured by the satellite-ground bistatic SAR system is the real deformation vector (three-dimensional vector ) is a scalar projected onto the bisector of the bistatic angle, i.e., the deformation direction measured by the SS-BSAR system is the bisector of the bistatic angle. In summary, the deformation inverted by the SS-BSAR single-angle re-orbit interferometry can be expressed as:

[0195]

[0196] S92. Calculate the annual average sedimentation rate v

[0197]

[0198] Among them, G is the design function, W is the weight function, is the observation value vector, i.e., the corresponding unwrapped phase after removing each error phase, and d is the deformation.

[0199] G contains the time information of the observed data and is used to describe the independent variables in the linear model. Its form is:

[0200]

[0201] W is used to represent the covariance matrix of the observed data, reflecting the correlation and uncertainty between the observed data. Its form is:

[0202]

[0203] Among them, σ ij represents the covariance between the i-th and j-th observations.

[0204] At this point, all steps are completed.

[0205] Next, an implementation example is given with specific parameters and specific data.

[0206] In this example, 28 mainlobe images of the Lutan-1 satellite-to-ground bistatic SAR were collected from November 24, 2023, to July 26, 2024, and the data from March 16, 2024, were selected as the main image. The location is Jinan City, Shandong Province, and the relevant parameters of the study area are shown in Table 1 below.

[0207] Table 1

[0208]

[0209] The date information and spatiotemporal baseline parameters of the 28 images are shown in Table 2 below.

[0210] Table 2

[0211]

[0212]

[0213] After executing step S1, 27 satellite-ground bistatic SAR differential interferograms with March 12, 2024 as the main image were obtained.

[0214] After executing steps S2 to S4, 27 differential interferograms were obtained after re-orbit interferometry error compensation, geocoding, image registration, and phase smoothing. Figure 5 shown.

[0215] After executing step S5, 42,729 PS candidate points are obtained after screening by the double threshold method of coherence coefficient and amplitude deviation, and the results are shown in Figure 6.

[0216] After executing step S6, the data matrix of the selected PS candidate point, including the time-space baseline, phase value, bistatic angle, observation angle, elevation value, amplitude deviation value, longitude and latitude values, etc., is obtained.

[0217] Execute steps S71 to S73 to perform combined filtering on the winding phase, calculate the terrain error phase of the PS candidate point, and perform secondary screening of the PS point according to the phase noise standard deviation. The results are as follows: Figure 7 shown.

[0218] After executing step S74, the unwrapped phase after spatially irrelevant viewing angle error correction and unwrapping is obtained. If Figure 8 shown.

[0219] After executing steps S75 to S77, the spatially correlated incident angle error phase, orbital inclination error phase, and linear atmospheric delay phase are obtained.

[0220] After executing step S8, the unwrapped phase with the error removed is obtained, as shown in the following example: Figure 9 shown.

[0221] After executing step S9, the average sedimentation rate map in the study area is obtained, and the results are as follows: Figure 10 As shown. Draw the average sedimentation rate histogram, the result is as follows Figure 11 shown.

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

Claims

1. A long time series processing method for satellite-ground bistatic SAR, characterized in that: The specific process is: Interferogram processing: N satellite-ground bistatic SAR re-orbit interferograms are acquired and error-compensated to obtain the corrected re-orbit interferograms. Based on the characteristic that targets in the same elevation direction have the same Doppler and bistatic distance, the corrected re-orbit interferograms are geocoded. Perform mask processing on the geocoded interferogram and perform phase filtering on the masked interferogram; PS candidate point screening and processing: The interferogram after phase filtering is screened using the amplitude deviation and coherence double threshold to obtain PS candidate points; the relevant data of the PS point is obtained and time series processed to obtain the terrain phase; Long time series processing: using the terrain phase to calculate the deformation, and calculating the average sedimentation rate based on the deformation.

2. The long time series processing method for satellite-ground bistatic SAR according to claim 1, characterized in that: The specific process of obtaining the corrected heavy track interferogram by error compensation is as follows: An interference error fringe estimation method based on local phase precision estimation is used. By performing frequency domain transformation on the interference pattern data along the range and azimuth directions, the linear trend in the frequency domain is extracted and eliminated to compensate for the interference phase error.

3. The long time series processing method for satellite-ground bistatic SAR according to claim 2, characterized in that: The specific process of obtaining the corrected heavy track interferogram by error compensation is as follows: int (n) performing a two-dimensional Fourier transform to estimate the linear frequency range in the azimuth / range direction, and setting the initial sampling frequency θ0 and the sampling frequency interval φ0 according to the linear frequency range; performing a Chirp-Z transform on the heavy track interferogram, and taking the frequency corresponding to the transformed peak as the frequency error estimate; finally, using the frequency error estimate result to remove the phase error at the heavy track interferogram to obtain a corrected interferogram.

4. The long time series processing method for satellite-ground bistatic SAR according to claim 3, characterized in that: Assume the linear frequency range is [f min ,f max ], Among them, N c is the number of sampling intervals; The frequency corresponding to the peak after transformation is used as the frequency error estimate: Among them, S int,r (n) with S int,a (n) represents the Chirp-Z transform of one-dimensional re-orbit interferometry data along the range and azimuth directions, respectively; The corrected interference pattern is s int_c , Among them, s int is the uncorrected re-orbit interferometry data, y and x are column vectors containing the azimuth and range sequence numbers of the interferogram, respectively, and 1 is a column vector whose elements are all 1.

5. The long time series processing method for satellite-ground bistatic SAR according to claim 1, characterized in that: The specific process of geocoding the corrected heavy track interferogram is as follows: An external digital elevation model (DEM) is used to project each DEM point onto the imaging plane along the elevation direction, and a mapping relationship is established between the longitude and latitude coordinate system and each coordinate point on the imaging plane. Based on the above mapping relationship, the heavy track interferogram is reversely projected into the longitude and latitude coordinate system to obtain the geocoded interferogram.

6. The long time series processing method for satellite-ground bistatic SAR according to claim 1, characterized in that: The interferogram after geocoding is subjected to mask processing, and phase filtering is performed on the interferogram after mask processing. The specific process is as follows: The digital elevation model is used to calculate the elevation and azimuth angles of each target point in the receiver's observation scene. At the same elevation and azimuth angles, the target point closest to the receiver is set as the visible point, and the remaining target points are set as blocked points. The interference pattern is masked, that is, the interference phase of the invisible area is set to 0, and the interference phase after masking is phase filtered.

7. The long time series processing method for satellite-ground bistatic SAR according to claim 6, characterized in that: The target point closest to the receiver is set as the visible point, and the remaining target points are blocked points, which is specifically expressed as: Among them, Flag(k) is the discrimination result of the kth DEM sampling point, Flag = 0 corresponds to an invisible point, Flag = 1 corresponds to a visible point, σ(k) = [δ(k); η(k)], δ(k) is the pitch angle of the kth DEM sampling point observed by the receiver, η(k) is the azimuth angle of the kth DEM sampling point observed by the receiver, R R (k) is the slant distance from the receiver to the kth DEM sampling point.

8. The long time series processing method for satellite-ground bistatic SAR according to claim 1, characterized in that: The use of the terrain phase Calculate deformation for: Where β is the bistatic angle and λ is the radar wavelength.

9. The long time series processing method for satellite-ground bistatic SAR according to claim 8, characterized in that: The annual average sedimentation rate v is: Among them, G is the design function and W is the weight function.

10. The long time series processing method for satellite-ground bistatic SAR according to claim 1, characterized in that: From the N acquired satellite-ground bistatic SAR re-orbit interferograms, one-dimensional re-orbit interferogram data in a flat area is selected and error compensation is performed to obtain the corrected re-orbit interferogram.

Images