A multi-dimensional deformation and differential tropospheric delay inversion method for distributed spaceborne D-InSAR

Abstract

The present invention provides a distributed spaceborne D-InSAR multidimensional deformation and differential tropospheric delay inversion method, which adaptively achieves high-precision and high-resolution joint inversion based on the spatial characteristics of multidimensional deformation and differential atmospheric delay, thereby improving the precision of multidimensional deformation inversion. The method comprises: performing D-InSAR phase unwrapping processing on re-orbit SAR interferometry data acquired by distributed spaceborne SAR satellites; estimating the ionospheric phase error in each satellite differential interferogram using an optimal spectrum splitting method, and performing ionospheric phase error compensation on each satellite differential interferogram after phase unwrapping; constructing a multi-channel adaptive minimum mean square error inversion matrix in the wavenumber domain based on the wavenumber domain power spectrum of each parameter to be estimated and the wavenumber domain power spectrum of the phase error component and the observation equation of the distributed spaceborne SAR system; and performing joint inversion of multidimensional deformation and differential tropospheric delay using the established minimum mean square error inversion matrix based on the differential interferogram obtained after processing.

Description

Technical Field

[0001] The present invention belongs to the technical field of synthetic aperture radar, and in particular relates to a multi-dimensional deformation and differential tropospheric delay inversion method of distributed spaceborne synthetic aperture radar differential interferometry (D-InSAR). Background Art

[0002] Since its introduction into remote sensing and mapping in the 1950s, synthetic aperture radar (SAR) has played a significant role in global environmental monitoring and disaster early warning, leveraging its all-day, all-weather, and wide-area observation capabilities. Spaceborne differential interferometric SAR (D-InSAR) utilizes phase difference interferometry patterns obtained from repeated orbits of SAR satellites. Using external terrain data to remove the unchanging terrain phase, D-InSAR can reveal minute surface deformations, making it an important tool for monitoring landslides, volcanic activity, and urban subsidence.

[0003] With the significant advancement of aerospace technology, the gradual maturation of SAR satellite manufacturing processes, and the growing demand for wide-area, high-temporal-spatial Earth observation, spaceborne SAR systems are gradually entering the distributed system era. Distributed spaceborne SAR systems, consisting of several to hundreds of SAR satellites working together, will become essential for ensuring wide-area, high-temporal-spatial Earth observation. Furthermore, and more importantly, utilizing distributed spaceborne SAR systems for D-InSAR surface deformation monitoring will overcome the bottleneck of traditional single-satellite D-InSAR, which only measures deformation along the line of sight. Instead, they can comprehensively assess multidimensional deformation in disaster scenarios, enhancing disaster early warning and awareness capabilities.

[0004] However, multidimensional deformation inversion using distributed spaceborne D-InSAR faces two challenges. First, as with traditional single-satellite D-InSAR, the differential atmospheric delay perturbation phase is of the same magnitude as, or even greater than, the scene deformation phase, resulting in the deformation signal being completely masked by the atmospheric delay signal. Second, the unified smoothing filter used in the joint inversion of multidimensional deformation using multiple satellites in distributed spaceborne D-InSAR results in the inability to observe small-scale structural deformation or poor noise reduction, creating a conflict between high-resolution and high-precision inversion. Therefore, research is needed on high-precision, high-resolution-preserving multidimensional deformation inversion for distributed spaceborne D-InSAR under differential atmospheric delay signal perturbations, which has not been addressed in existing distributed spaceborne D-InSAR research. Summary of the Invention

[0005] To solve the above problems, the present invention provides a multidimensional deformation and differential tropospheric delay inversion method for distributed spaceborne D-InSAR. This inversion method can adaptively achieve high-precision and high-resolution joint inversion according to the spatial characteristics of multidimensional deformation and differential atmospheric delay, thereby improving the fineness of multidimensional deformation inversion.

[0006] The present invention is achieved through the following technical solutions.

[0007] A distributed spaceborne D-InSAR multi-dimensional deformation and differential tropospheric delay inversion method includes:

[0008] Step 1: Perform D-InSAR phase unwrapping on the heavy-orbit SAR interferometric data acquired by the distributed spaceborne SAR satellite;

[0009] Step 2: Use the optimal spectrum decomposition method to estimate the ionospheric phase error in the differential interferogram of each satellite, and perform ionospheric phase error compensation on the differential interferogram of each satellite after phase unwrapping;

[0010] Step 3: Construct a multi-channel adaptive minimum mean square error inversion matrix in the wavenumber domain based on the wavenumber domain power spectrum of each parameter to be estimated, the wavenumber domain power spectrum of the phase error component, and the observation equation of the distributed spaceborne SAR system;

[0011] Step 4: Based on the differential interferogram obtained after processing in step 2, the minimum mean square error inversion matrix established in step 3 is used to jointly invert the multidimensional deformation and differential tropospheric delay.

[0012] Beneficial effects of the present invention:

[0013] The present invention accurately models the multi-source errors of distributed spaceborne D-InSAR and, based on this, constructs a multi-channel adaptive multidimensional deformation and differential tropospheric delay joint inversion coefficient matrix under the minimum mean square error criterion, thereby achieving high-precision and high-resolution joint inversion of multidimensional deformation and differential atmospheric delay, and improving the refinement of multidimensional deformation inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 This is a flow chart of the multi-dimensional deformation and differential tropospheric delay inversion method of the distributed spaceborne D-InSAR of the present invention;

[0015] Figure 2 A diagram of a multi-dimensional deformation and differential tropospheric delay simulation scenario in a specific embodiment of the present invention;

[0016] Figure 3 This is the joint inversion result of dimensional deformation and differential tropospheric delay in a specific embodiment of the present invention. DETAILED DESCRIPTION

[0017] The present invention will be described in detail below with reference to the accompanying drawings.

[0018] like Figure 1 As shown, a multi-dimensional deformation and differential tropospheric delay inversion method of a distributed spaceborne D-InSAR according to a specific embodiment of the present invention specifically includes the following steps:

[0019] Step 1: Perform D-InSAR phase unwrapping on the heavy-orbit SAR interferometric data acquired by the distributed spaceborne SAR satellite;

[0020] In this embodiment, the D-InSAR preprocessing specifically includes the following steps:

[0021] 1.1 Use the cross-correlation method or the interferometric data registration method based on the external digital elevation model (DEM) to register the data of multiple distributed spaceborne SAR interferometric pairs;

[0022] 1.2 Conjugate multiply the primary and secondary images of each distributed SAR satellite interferometric pair registered in 1.1 to generate an interferogram;

[0023] 1.3 Using the precise orbit determination ephemeris data of each distributed SAR satellite and external DEM data, the reference terrain and flat-earth interferometric phases of different distributed satellite interferometric pairs are calculated according to the SAR range-Doppler equation.

[0024] 1.4 Conjugate multiply the interferogram of each satellite generated in 1.2 by the corresponding reference terrain and flat ground interferogram phase calculated in 1.3 to remove the terrain and flat ground phase information, thereby obtaining a differential interferogram of each distributed SAR satellite;

[0025] 1.5 Performing phase filtering on the distributed SAR satellite differential interferograms obtained in 1.4 using the Goldstein phase filtering method;

[0026] 1.6 Use the least squares phase unwrapping method to perform phase unwrapping on the interference pattern after phase filtering in 1.5.

[0027] Step 2: Use the optimal spectrum decomposition method to estimate the ionospheric phase error in the differential interferogram of each satellite, and perform ionospheric phase error compensation on the differential interferogram of each satellite after phase unwrapping;

[0028] The rationale behind this step is that, due to the high altitude of the ionosphere, the correlation between the ionospheric phase components observed by different satellites in a distributed SAR system is weak. Therefore, prior to joint inversion, ionospheric phase error estimation and compensation must be performed based on the optimal spectral splitting method. Since the spectral splitting method uses spatial averaging to estimate the ionospheric phase error, it acts as a low-pass filter. To ensure a certain spatial resolution, the spatial smoothing window of the jointly estimated interferogram is smaller than that used in the spectral splitting method. Therefore, after compensation, the residual ionospheric phase is determined by both the estimation accuracy of the spectral splitting method and the residual high-frequency ionospheric phase component.

[0029] In this embodiment, the ionospheric phase error in each satellite differential interferogram is estimated using the optimal spectrum separation method, specifically:

[0030] Under the condition of spatial smoothness constraint, the window size g of the spectrum separation method is x and g y Determined by minimizing the ionospheric residual error:

[0031]

[0032] The integration domain D corresponds to the space grid size (2π / g x and 2π / g y ) and is smaller than the area of ​​the joint estimation space grid size,

[0033]

[0034] PR(k x , k y ) is represented by the wavenumber domain power spectrum of the Rino ionospheric phase, and the input parameters are the ionospheric scintillation CkL index, spectral index, outer scale of ionospheric irregularities and geographical location information of the observation scene.

[0035] Step 3: Construct a multi-channel adaptive minimum mean square error inversion matrix in the wavenumber domain based on the wavenumber domain power spectrum of each parameter to be estimated, the wavenumber domain power spectrum of the phase error component, and the observation equation of the distributed spaceborne SAR system. Specifically, the following steps are included:

[0036] 3.1 Estimating the variance of the differential tropospheric delay phase based on an invariant reference region in the scene The power spectrum of the differential tropospheric delayed phase is:

[0037]

[0038] in is the power spectrum coefficient, k x and k y is the two-dimensional wave number.

[0039] 3.2 Based on the scene deformation measurements obtained by the prior deformation model or other measurement methods (such as navigation satellites), the wavenumber domain power spectrum of the multidimensional deformation field is calculated using the periodogram; the specific formula is:

[0040]

[0041] Where d represents the deformation variable in a certain dimension obtained by a priori scene deformation model or other measurement methods.

[0042] 3.3 Estimate the wavenumber domain power spectrum of each error source; specifically including:

[0043] 3.3.1 Based on the coherence coefficient of each satellite interferogram and the multi-look number, the wavenumber domain power spectrum of the phase noise is estimated to be a white noise spectrum; the specific formula is:

[0044]

[0045] Among them, N L is the number of multi-views, γ is the coherence coefficient, K x and K y is the support domain range of the two-dimensional wavenumber domain;

[0046] 3.3.2 Calculate the maximum baseline phase error for each star based on the precise ephemeris data; the specific formula is:

[0047]

[0048] Where λ is the carrier wavelength, θ l is the viewing angle from the SAR satellite, e h and e v are the vertical baseline error and the horizontal baseline error, respectively. x and y are the positions of the scene targets. Assuming that the satellite re-orbit errors are independent and stable, then and σ h and σ v They are the precise ephemeris cross-orbit and radial orbit determination accuracies, using e h and e v The orbit error estimate of the scene is obtained by the 1σ standard deviation, and the estimation of the power spectrum of the phase wavenumber domain of the baseline error is obtained by Fourier transform in the spatial domain; the specific formula is:

[0049]

[0050] 3.3.3 Based on the external ionospheric data and the spectrum decomposition method parameters, the wavenumber domain power spectrum of the residual ionospheric phase error is calculated. The specific formula is:

[0051] P NI (k x , k y )=P R (k x , k y )|F s (k x , k y )| 2

[0052] Among them, F s is the spatial smoothing filter transfer function used in the spectrum separation method.

[0053] 3.3.4 Based on the observation geometry of distributed SAR satellites, calculate the troposphere penetration distance of the observation line of sight between any two satellites in the distributed SAR system. The specific formula is:

[0054]

[0055] Where Δx x and Δx y are the projection distances of the puncture point vector in the troposphere in two dimensions, d x and d y are the spatial distances of distributed SAR satellites in two dimensions, h s is the satellite orbit altitude, h tr is the effective height of the troposphere.

[0056] 3.3.5 Based on the observation geometry of the distributed SAR satellites and the wavenumber domain power spectrum of the differential tropospheric delay phase, the wavenumber domain power spectrum of the differential tropospheric delay decorrelation phase is estimated. The specific formula is:

[0057]

[0058] in, is the distance vector of the puncture point in the equivalent troposphere observed by the line of sight of two distributed SAR satellites.

[0059] 3.4 Based on the wavenumber domain power spectrum of each parameter to be estimated and the wavenumber domain power spectrum of the phase error component and the observation equation of the distributed spaceborne SAR system, the multi-channel adaptive inversion coefficient matrix is ​​calculated; the specific formula is:

[0060]

[0061] in, is the power spectrum of the parameter to be estimated, H ij is an element in the observation matrix determined by the observation geometry of the distributed spaceborne SAR system. The power spectrum of the differential interferogram after ionospheric phase error compensation in the wavenumber domain is:

[0062]

[0063] In the case of jointly estimating the two-dimensional deformation (LOS and along-track direction of the reference satellite) and the differential tropospheric delay, it is expressed as a matrix consisting of the cross-spectral variables of the parameters to be estimated in the wavenumber domain:

[0064]

[0065] Then the sum of the wavenumber domain power spectra of each error source is:

[0066]

[0067] Step 4: Input the differential interferogram obtained after processing in step 2, and use the minimum mean square error inversion matrix established in step 3 to jointly invert the multidimensional deformation and differential tropospheric delay.

[0068] In this embodiment, the joint inversion of multidimensional deformation and differential tropospheric delay specifically includes the following steps:

[0069] 4.1 Transform W(k) obtained in step 3 into the inversion coefficient matrix w in the spatial domain;

[0070] 4.2 The differential interferogram of each satellite of the distributed SAR obtained after step 2 Respectively with the corresponding inversion coefficient matrix w i Perform spatial convolution operations;

[0071] 4.3 Group and sum the results of the spatial convolution operation of each distributed SAR satellite to obtain estimates of multidimensional deformation and differential tropospheric delay Specifically:

[0072]

[0073] in Represents the spatial convolution operation.

[0074] Embodiment 1:

[0075] In this example, the distributed SAR satellite orbit altitude is 693 km and consists of three satellites evenly spaced along the orbital axis, with only the center satellite transmitting. The radar payload parameters are shown in Table 1. The simulation generates the reference satellite's LOS, azimuth deformation, and zenith-direction differential tropospheric delay. In this example, the simulated differential interferometry data used to implement this method already includes phase noise, baseline error, residual ionospheric phase error, and differential tropospheric delay decorrelation phase.

[0076] Table 1

[0077] parameter Numerical Carrier center frequency (GHz) 5.4 Bandwidth (MHz) 50 Scene size (km) 50 Grid size (m) 500 View count 25 (distance) × 50 (direction)

[0078] Figure 3 The inversion results for multidimensional deformation are presented for two intersatellite distances: 1) 350 km (first row) and 2) 50 km (second row). The corresponding inversion root mean square errors (RMS) for the reference satellite's LOS, azimuth deformation, and zenith differential tropospheric delay are 3.9 mm, 2.8 mm, and 3.5 mm, and 4.3 mm, 3.1 mm, and 3.8 mm, respectively. The inverted multidimensional deformation and differential tropospheric delay are very close to the reference values ​​set by the simulation, demonstrating the effectiveness of this method.

[0079] Of course, the present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A method for multi-dimensional deformation and differential tropospheric delay inversion for distributed spaceborne D-InSAR, characterized in that: The following steps are involved: Step 1: Perform D-InSAR phase unwrapping on the heavy-orbit SAR interferometric data acquired by the distributed spaceborne SAR satellite; Step 2: Use the optimal spectrum decomposition method to estimate the ionospheric phase error in the differential interferogram of each satellite, and perform ionospheric phase error compensation on the differential interferogram of each satellite after phase unwrapping; Step 3: Construct a multi-channel adaptive minimum mean square error inversion matrix in the wavenumber domain based on the wavenumber domain power spectrum of each parameter to be estimated, the wavenumber domain power spectrum of the phase error component, and the observation equation of the distributed spaceborne SAR system; The step three specifically includes the following steps: 3.1 Estimate the variance and power spectrum of the differential tropospheric delay phase based on the non-deformed reference area in the scene; 3.2 Based on the scene deformation measurements obtained by the prior deformation model or other measurement methods, the wavenumber domain power spectrum of the multidimensional deformation field is calculated using the periodogram; 3.3 Estimate the wavenumber domain power spectrum of each error source; 3.4 Calculate the multi-channel adaptive inversion coefficient matrix based on the wavenumber domain power spectrum of each parameter to be estimated, the wavenumber domain power spectrum of the phase error component, and the observation equation of the distributed spaceborne SAR system; Step 4: Based on the differential interferogram obtained after processing in step 2, the minimum mean square error inversion matrix established in step 3 is used to jointly invert the multidimensional deformation and differential tropospheric delay.

2. The method according to claim 1, wherein The D-InSAR phase unwrapping process specifically includes the following steps: 1.1 Use the cross-correlation method or the interferometric data registration method based on the external digital elevation model (DEM) to register the data of multiple distributed spaceborne SAR interferometric pairs; 1.2 Conjugate multiply the primary and secondary images of each distributed SAR satellite interferometric pair registered in step 1.1 to generate an interferogram; 1.3 Using the precise orbit determination ephemeris data of each distributed SAR satellite and external DEM data, the reference terrain and flat-earth interferometric phases of different distributed satellite interferometric pairs are calculated according to the SAR range-Doppler equation. 1.4 Conjugate multiply the interferogram of each satellite generated in step 1.2 by the corresponding reference terrain and flat-ground interferogram phase calculated in step 1.3 to remove the terrain and flat-ground phase information, thereby obtaining a differential interferogram of each distributed SAR satellite; 1.5 Perform phase filtering on the distributed SAR satellite differential interferograms obtained in step 1.4 using the Goldstein phase filtering method; 1.6 Use the least squares phase unwrapping method to perform phase unwrapping on the interferogram after phase filtering in step 1.

5.

3. The method according to claim 1 or 2, wherein: The optimal spectrum splitting method is used to estimate the ionospheric phase error in the differential interferogram of each satellite. The specific formula is: Among them, under the condition of spatial smoothness constraint, the window size g of the spectrum separation method is x and g y Determined by minimizing the ionospheric residual error, the integration domain D corresponds to the region where the wave number is larger than the spatial grid size of the spectrum separation method and smaller than the spatial grid size of the joint estimation; Where, γ is the coherence coefficient, N L It is the number of multi-views; P R (k x ,k y ) is represented by the wavenumber domain power spectrum of the Rino ionospheric phase, and the input parameters are the ionospheric scintillation CkL index, spectral index, outer scale of ionospheric irregularities and geographical location information of the observation scene.

4. The method according to claim 1, wherein The wavenumber domain power spectrum of each error source is estimated in the following manner: 3.3.1 The wavenumber domain power spectrum of the phase noise is estimated to be a white noise spectrum based on the coherence coefficient of each satellite interferogram and the multi-look number; 3.3.2 Calculate the maximum baseline phase error for each star based on the precise ephemeris data; 3.3.3 Calculate the wavenumber domain power spectrum of the residual ionospheric phase error based on the external ionospheric data and the spectrum decomposition method parameters; 3.3.4 Based on the observation geometry of distributed SAR satellites, calculate the tropospheric penetration distance of the observation line of sight between any two satellites in the distributed SAR system; 3.3.5 Estimate the wavenumber domain power spectrum of the differential tropospheric delay decorrelation phase based on the observation geometry of the distributed SAR satellites and the wavenumber domain power spectrum of the differential tropospheric delay phase.

5. The method according to claim 1 or 2, wherein: The joint inversion of the multidimensional deformation and the differential tropospheric delay specifically includes the following steps: 4.1 Transform the calculation obtained in step 3 into an inversion coefficient matrix in the spatial domain; 4.2 Performing spatial convolution operations on the differential interferograms of each distributed SAR satellite obtained after processing in step 2 and the corresponding inversion coefficient matrix; 4.3 The results of the spatial convolution operation of each distributed SAR satellite are grouped and summed to obtain estimates of the multidimensional deformation and differential tropospheric delay.

Images