Method for determining a Cramer-Rao bound of distributed spaceborne InSAR height retrieval

Abstract

The present application belongs to the technical field of synthetic aperture radar, and particularly relates to a method for determining the Cramer-Rao bound of distributed spaceborne InSAR height inversion. The present application can solve the problem that the traditional CRB calculation method cannot obtain accurate lower bound estimation under the non-Gaussian and non-stationary atmospheric condition by combining the post-CRB method with the idea of Gaussian mixture model. Since the independent calculation under the Gaussian assumption is performed on a window-by-window basis in space, the stationary and Gaussian assumption for the atmospheric phase of the whole interferogram is avoided, and the estimation process of the height inversion under the influence of the non-stationary and complex distribution atmospheric in the actual process is more close to the actual process.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar technology, and particularly relates to a method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion. Background Technology

[0002] Digital elevation models (DEMs), as important geographic information products, play a vital role in fields such as geology and meteorology. Distributed spaceborne synthetic aperture radar (SAR) systems, as a relatively new spaceborne SAR technology, possess flexible and versatile system configurations and the ability to acquire high spatiotemporal resolution data, thus offering advantages in achieving high-precision and high-resolution elevation inversion.

[0003] In elevation inversion using distributed spaceborne SAR, the atmosphere significantly impacts the accuracy of the inversion; therefore, accuracy evaluation methods that consider atmospheric characteristics are crucial. The Cramer-Rao boundary (CRB) is an important metric for evaluating parameter performance and can be used to measure the performance of unbiased estimators. Reference 1 (GUARNIERI AM, TEBALDINI S. Hybrid...) Bounds for Crustal Displacement Field Estimators in SAR Interferometry[J / OL].IEEE Signal Processing Letters, 2007, 14(12): 1012-1015.DOI: 10.1109 / LSP.2007.904705.) presents a method for evaluating line-of-sight deformation to the surface using the hybrid Cramer-Rao boundary (HCRB). This method considers the atmospheric phase screen as a random parameter and provides a feasible performance evaluation. Reference 2 (PRATS-IRAOLA P, LOPEZ-DEKKER P, DE ZAN F, et al. Performance of 3-D Surface Deformation Estimation for Simultaneous Squinted SAR Acquisitions[J / OL].IEEE Transactions on Geoscience and Remote Based on the work in Sensing, 2018, 56(4):2147-2158.DOI:10.1109 / TGRS.2017.2776140.), further research was conducted using HCRB and the formula was extended to obtain the three-dimensional deformation performance. To our current knowledge, existing CRB assessment methods that consider atmospheric influences all treat atmospheric phase as a random variable and assume that it has Gaussian properties.

[0004] However, in reality, the atmospheric phase composition is complex. Reference 3 (MULDER G, VAN LEIJEN FJ, BARKMEIJER J, et al. Estimating Single-Epoch Integrated Atmospheric Refractivity from InSAR for Assimilation in Numerical Weather Models[J / OL]. IEEE Transactions on Geoscience and Remote Sensing, 2022:1-1. DOI:10.1109 / TGRS.2022.3177041.) demonstrates that the atmosphere does not possess Gaussian properties, therefore the Gaussian assumption does not hold in CRB calculations. Furthermore, due to the variability of meteorological environments, different meteorological conditions such as wind, rain, and snow lead to different atmospheric statistical characteristics; spatially and temporally, the atmospheric phase is not a stationary process. Therefore, in order to evaluate the unbiased estimator for distributed spaceborne InSAR elevation inversion, a CRB calculation method considering the non-stationary, non-Gaussian characteristics of the atmosphere is crucial. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion. This method no longer makes Gaussian and stationary assumptions about the atmospheric phase in the entire interferogram. Instead, it adopts a posterior CRB approach in space, using a Gaussian mixture model to fit the error distribution to obtain the final CRB result.

[0006] The technical solution of this invention is:

[0007] A method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion, comprising the following steps:

[0008] Step 1: Obtain the reference digital elevation model H for the distributed spaceborne InSAR observation area, and derive the digital elevation model to be estimated based on the obtained reference digital elevation model H.

[0009] The reference digital elevation model H for the observation area can be obtained from the Space Shuttle Radar Topography (SRTM) data;

[0010] The process involves obtaining the digital elevation model to be estimated based on the obtained reference digital elevation model H. The method is as follows:

[0011] Let [H] be the pixel in the m-th row and n-th column of the reference digital elevation model H. m,n =h m,nDigital elevation model to be estimated The cell in the m-th row and n-th column is but

[0012]

[0013] Wherein, △h m,n =h m,n+1 -h m,n ; ∈ is a parameter between 0 and 1, v m,n for The process noise caused by the error between H and the mean is a process noise with zero mean and variance. Gaussian noise model, zero mean and variance All of these are related to the precision of H;

[0014] Step 2, calculate the wavenumber of distributed spaceborne InSAR, the specific method is as follows:

[0015] Obtain the interferogram of the observation area described in step 1 from distributed spaceborne InSAR. The total number of interferograms is denoted as N; in the k-th interferogram, the pixel in the m-th row and n-th column is represented as... Where k = 1, 2, ..., N; then the interference phase composition vector obtained at all pixels in the m-th row and n-th column is...

[0016] Let λ be the radar wavelength corresponding to each interferogram. k First, calculate the interferometric baseline between radars corresponding to each interferogram, and denote the k-th interferogram as X. k The corresponding interferometric baseline is B. k Next, calculate the slant distance from the center of the main track aperture to the center of the scene for each interferogram, and denote the k-th interferogram as X. k The corresponding slope distance is R k Finally, calculate the incident angle corresponding to each interferogram, and the X of the k-th interferogram. k The angle of incidence is denoted as θ k Based on the above parameters, the interferogram X is calculated. k The corresponding wave number ξ k :

[0017]

[0018] Step 3, based on the reference digital elevation model H obtained in Step 1 and the wave number ξ obtained in Step 2... k Calculate the phase difference at the m-th row and n-th column pixel in the k-th interferogram.

[0019]

[0020] Step 4, adjust the phase difference obtained in step 3. Error phase processing is performed to obtain the mean value. variance is The Gaussian distribution; specifically:

[0021] For interferograms For the pixel in the m-th row and n-th column of the k-th image, select the phase difference interferogram within a W×W window surrounding it. (W is an odd number), denoted as Then we have:

[0022]

[0023] Using Gaussian mixture model to By fitting the data, a suitable distribution is obtained. To make subsequent calculations of the expected value more convenient, based on the central limit theorem, and since the number of pixels in the W×W window is very small relative to the total number of pixels in the entire interferogram, it is assumed that the atmospheric signal within the window is stationary and the error phase approximately conforms to a Gaussian distribution. Thus, the suitable distribution is a Gaussian distribution. Furthermore, since the assumption of a Gaussian distribution within the window is made for each pixel independently, it can also represent the non-stationary characteristics of the atmospheric phase of the entire interferogram.

[0024] Suppose that after fitting, for the pixel in the m-th row and n-th column, its The distribution of the mean is variance is Gaussian distribution;

[0025] Step 5, from the digital elevation model to be estimated Starting from the first column of each row, CRB is calculated pixel by pixel. The specific calculation method is as follows:

[0026] Record the digital elevation model to be estimated Let the Fisher information content be J, and let the Fisher information content in the m-th row and n-th column be [J]. m,n =J m,n In the CRB calculation process, the parameters to be estimated are decomposed into The corresponding Fisher information content is defined as:

[0027]

[0028] For the m-th row and n-th column, we have:

[0029]

[0030]

[0031]

[0032]

[0033] The Fisher information in the m-th row and n+1-th column is obtained by recursion:

[0034]

[0035] Taking the reciprocal of the Fisher information content, we obtain the CRB at each pixel:

[0036]

[0037] Step 6: Compare the digital elevation models to be estimated To determine the accuracy of elevation inversion, based on the definition of mean square error (MSE), the mean square error (MSE) is compared with that of the digital elevation model to be estimated. Average CRB:

[0038]

[0039] Step 7: Perform topographic mapping based on the CRB obtained in Step 6.

[0040] The beneficial effects of this invention are as follows:

[0041] This invention addresses the problem of traditional CRB calculation methods failing to accurately estimate the lower bound under non-Gaussian and non-stationary atmospheric conditions by combining the posterior CRB method with a Gaussian mixture model. Because it performs independent calculations under the Gaussian assumption on a pixel-by-pixel basis in space, it avoids applying stationary and Gaussian assumptions to the atmospheric phase of the entire interferogram, thus more closely reflecting the actual process of elevation inversion under the influence of non-stationary and complex atmospheric conditions. Attached Figure Description

[0042] Figure 1 This is a schematic diagram of the re-orbit interferometry of a distributed spaceborne SAR system.

[0043] Figure 2 To compare the elevation inversion results with the actual elevation using 15 interferograms under different atmospheric stationary conditions;

[0044] Figure 3 To compare the CRB calculation results and the actual elevation inversion results under different atmospheric stability conditions using different numbers of interferograms. Detailed Implementation

[0045] The invention will now be described in detail with reference to the accompanying drawings.

[0046] Achieving distributed spaceborne SAR elevation inversion requires multiple interferometric maps, which can be obtained by multiple satellites in a distributed system performing several re-orbit interferometry operations on the same area, such as... Figure 1As shown. The specific steps of this invention are as follows:

[0047] Step 1: Obtain a reference digital elevation model H for the radar observation area, such as Space Shuttle Radar Topography (SRTM) data; Let the digital elevation model to be estimated be... The cell in the m-th row and n-th column of both is [H]. m,n =h m,n and Based on this, the state equation is obtained:

[0048]

[0049] Wherein, △h m,n =h m,n+1 -h m,n ; ∈ is a parameter between 0 and 1, v m,n for The process noise caused by the error between H and H is modeled as having zero mean and variance. Both the Gaussian noise model and the accuracy of H are related.

[0050] Step 2, calculate wavenumbers. Obtain interferograms of the same region from the distributed spaceborne SAR system. The total number of interferograms is denoted as N; the pixel in the m-th row and n-th column of the k-th interferogram can be represented as... Where k = 1, 2, ..., N; then the interference phases obtained at all pixels in the m-th row and n-th column can form a vector.

[0051] Let λ be the radar wavelength corresponding to each interferogram. k First, calculate the interferometric baseline between radars corresponding to each interferogram, and denote the k-th interferogram as X. k The corresponding interferometric baseline is B. k Next, calculate the slant distance from the center of the main track aperture to the center of the scene for each interferogram, and denote the k-th interferogram as X. k The corresponding slope distance is R k Finally, calculate the incident angle corresponding to each interferogram, and the X of the k-th interferogram. k The angle of incidence is denoted as θ k The calculation method described above can be found in reference 4 (R. Hanssen, Radar Interferometry Data Interpretation and Error Analysis. 2001.); based on the above parameters, the interferogram X is calculated. k Corresponding wave number:

[0052]

[0053] Step 3: Calculate the phase difference based on the obtained interferogram and the calculated terrain phase. The phase difference at the pixel in the m-th row and n-th column of the k-th interferogram is denoted as... but The phase difference at point can be expressed as

[0054] Step 4: Process the phase error of the interferogram. For the pixel in the m-th row and n-th column of the k-th image, select the phase difference within a W×W window around it (W is an odd number), denoted as... Then we have:

[0055]

[0056] Using Gaussian mixture model to A fitting is performed to obtain a suitable distribution. To facilitate subsequent calculations of the expected value, based on the central limit theorem, and since the number of pixels within the W×W window is relatively small compared to the total number of pixels in the entire interferogram, it can be assumed that the atmospheric signal within the window is stationary and the error phase approximately follows a Gaussian distribution. Furthermore, since the assumption of an independent Gaussian distribution within the window is made for each pixel, this can also represent the non-stationary characteristics of the atmospheric phase in the entire interferogram.

[0057] After fitting, for the pixel in the m-th row and n-th column, its The distribution of the mean is variance is The Gaussian distribution.

[0058] Step 5: Starting from the first column of each row, calculate the CRB pixel by pixel. Let J be the Fisher information content of the entire image, and let [J] be the Fisher information content of the m-th row and n-th column. m,n =J m,n In the CRB calculation process, the parameters to be estimated are decomposed into... The corresponding Fisher information content can then be defined as:

[0059]

[0060] For the m-th row and n-th column, we have:

[0061]

[0062]

[0063]

[0064]

[0065] The Fisher information in the m-th row and n+1-th column can be obtained recursively:

[0066]

[0067] Taking the reciprocal of Fisher information to obtain the CRB at each pixel:

[0068]

[0069] To compare the elevation inversion accuracy of the entire map, according to the definition of mean square error (MSE), we can compare the MSE with the average CRB of the entire map:

[0070]

[0071] Example

[0072] The following provides an implementation example with specific parameters.

[0073] In this example, we consider a three-satellite distributed spaceborne SAR system in the same orbital plane. The satellite orbital elements and platform parameters are shown in Table 1. It is assumed that the three-satellite system performs five reorbit interferometry maneuvers over six days, obtaining a total of 15 interferograms. Each interferogram is 182 pixels × 182 pixels with a pixel spacing of 1 meter. The elevation inversion scene is centered at 69°E, 30°N, and is selected from the 30-meter resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM), scaled down. Images of the simulated scene are shown below. Figure 2 .

[0074] Table 1

[0075]

[0076] First, following step 1, obtain a 30m resolution DEM of the observation area as the true elevation, then scale it down proportionally to a maximum height of 18.3m. The reference DEM is obtained by averaging the true elevation using a window.

[0077] Perform step 2 to calculate the wavenumber corresponding to each interferogram.

[0078] Perform step 3 to calculate the phase difference between the actual obtained interferogram and the calculated terrain phase.

[0079] Execute step 4 and select a window size W = 21. Then, for each pixel position, use the phase difference within the 21×21 window around it to fit the variance and mean of the pixel at each position.

[0080] Perform step 5 to calculate the final Cramero boundary for elevation inversion.

[0081] Figure 2The elevation inversion results and actual elevations for two stationary conditions are presented under 15 interferograms. The MSE of both conditions is close to that of the CRB, but the MSE is smaller and the accuracy is higher under the stationary atmospheric condition. Figure 3 The paper presents the CRB calculated by this method and the MSE obtained by elevation inversion under stationary and non-stationary atmospheric conditions, using different numbers of interferograms. It can be seen that the CRB calculated by this method decreases with increasing number of interferograms; the CRB under non-stationary atmospheric conditions is higher than that under stationary atmospheric conditions; and experimental results demonstrate that the MSE obtained by elevation inversion is close to the CRB in both cases, thus proving the correctness of the CRB.

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

Claims

1. A method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion, characterized in that... The steps of this method include: Step 1: Obtain the reference digital elevation model H for the distributed spaceborne InSAR observation area, and derive the digital elevation model to be estimated based on the obtained reference digital elevation model H. Step 2: Calculate the wavenumber of the interferogram for the distributed spaceborne InSAR observation area; Step 3: Calculate the phase difference of the first interferogram based on the reference digital elevation model H obtained in Step 1 and the wavenumber obtained in Step 2; Step 4: Perform error phase processing on the phase difference obtained in Step 3 to obtain a Gaussian distribution that the phase difference conforms to; Step 5, from the digital elevation model to be estimated Starting from the first column of each row, calculate CRB pixel by pixel; Step 6: Determine the digital elevation model to be estimated Average CRB value: In step 1, the reference digital elevation model H of the observation area can be obtained from the Space Shuttle Radar Topographic Mapping (SRTM) data. Based on the obtained reference digital elevation model H, the digital elevation model to be estimated is obtained. The method is as follows: Let [H] be the pixel in the m-th row and n-th column of the reference digital elevation model H. m,n =h m,n Digital elevation model to be estimated The cell in the m-th row and n-th column is but Where, Δh m,n =h m,n+1 -h m,n ; ∈ is a parameter between 0 and 1, v m,n for Process noise caused by errors in H.

2. The method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion according to claim 1, characterized in that: In step 2, the method for calculating the wavenumber of the interferogram in the distributed spaceborne InSAR observation area is as follows: Obtain the interferogram of the observation area described in step 1 from distributed spaceborne InSAR. The total number of interferograms is denoted as N; in the k-th interferogram, the pixel in the m-th row and n-th column is represented as... Where k = 1, 2, ..., N; then the interference phase composition vector obtained at all pixels in the m-th row and n-th column is... Let λ be the radar wavelength corresponding to each interferogram. k First, calculate the interferometric baseline between radars corresponding to each interferogram, and denote the k-th interferogram as X. k The corresponding interferometric baseline is B. k Next, calculate the slant distance from the center of the main track aperture to the center of the scene for each interferogram, and denote the k-th interferogram as X. k The corresponding slope distance is R k Finally, calculate the incident angle corresponding to each interferogram, and the X of the k-th interferogram. k The angle of incidence is denoted as θ k Based on the above parameters, the interferogram X is calculated. k The corresponding wave number ξ k :

3. The method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion according to claim 2, characterized in that: In step 3, calculating the phase difference of the interferogram based on the reference digital elevation model H obtained in step 1 and the wavenumber obtained in step 2 refers to calculating the phase difference of the interferogram based on the reference digital elevation model H obtained in step 1 and the wavenumber ξ obtained in step 2. k Calculate the phase difference at the m-th row and n-th column pixel in the k-th interferogram.

4. The method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion according to claim 3, characterized in that: In step 4, the phase difference obtained in step 3 is... Error phase processing is performed to obtain the mean value. variance is The Gaussian distribution method is as follows: For interferograms For the pixel in the m-th row and n-th column of the k-th image, select the phase difference interferogram within a W×W window surrounding it. W is an odd number, denoted as Then we have: Suppose that after fitting, for the pixel in the m-th row and n-th column, its The distribution of the mean is variance is The Gaussian distribution.

5. The method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion according to claim 4, characterized in that: In step 5, the digital elevation model to be estimated is used... Starting from the first column of each row, CRB is calculated pixel by pixel. The specific calculation method is as follows: Record the digital elevation model to be estimated Let the Fisher information content be J, and let the Fisher information content in the m-th row and n-th column be [J]. m,n =J m,n In the CRB calculation process, the parameters to be estimated are decomposed into The corresponding Fisher information content is defined as: For the m-th row and n-th column, we have: The Fisher information in the m-th row and n+1-th column is obtained by recursion: Taking the reciprocal of the Fisher information content, we obtain the CRB at each pixel:

6. The method for determining the Cramer-Rao boundary in distributed spaceborne InSAR elevation inversion according to claim 5, characterized in that: Topographic mapping was performed based on the obtained CRB.

Images