Long-wave sar atmospheric error correction method and device, equipment and medium

Abstract

The application discloses a long-wave SAR atmospheric error correction method, device, equipment and medium, comprising: acquiring a multi-temporal dual-polarized SAR image data set, and introducing a reference DEM for differential interference; for each polarization channel, using coherence as a weight to weight the multi-temporal differential interferogram to estimate an initial coarse terrain residual; using a bare area to fit an orbit error, and estimating a coarse terrain residual removing the orbit error; according to the differential interference phase of a single pair of interference pairs and the coarse terrain residual, calculating a composite residual phase; performing wavelet correlation analysis on the composite residual phase of the dual-polarization channel, identifying a common-mode error signal composed of atmospheric delay and orbit error, and stripping the common-mode error signal from the composite residual phase to obtain a high-frequency terrain residual compensation phase; converting the high-frequency terrain residual compensation phase to the height domain, and superimposing the high-frequency terrain residual compensation phase on the reference DEM to generate a high-precision DEM. The application can effectively separate the atmospheric delay, the orbit and the local terrain residual under sparse time sampling.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing topographic mapping technology, specifically relating to a long-wave SAR atmospheric error correction method, device, equipment, and medium. Background Technology

[0002] Forest ecosystems play a crucial role in regulating global carbon cycles and climate change. Accurately acquiring true forest understory topography (digital terrain models, DTM) and forest aboveground biomass (AGB) is fundamental for estimating global forest carbon storage and understanding ecosystem dynamics. Among numerous Earth observation technologies, synthetic aperture radar (SAR), with its all-weather, all-time capability and penetrating power, has become a core tool for forest mapping and parameter inversion. Traditional short-wavelength SAR (such as X and C bands) is limited by wavelength and cannot penetrate dense forest canopies, while long-wavelength SAR systems, represented by L and P bands, offer significant opportunities for topographic mapping in complex, penetrable environments. In particular, BIOMASS, the world's first P-band fully polarimetric spaceborne SAR satellite launched by the European Space Agency (ESA), with a wavelength of approximately 70 centimeters, possesses extremely strong canopy penetration capabilities.

[0003] However, in spaceborne repeated orbit interferometry configurations, the long revisit cycles of several days and complex environmental fluctuations often lead to degradation in interferometric quality. After introducing an external reference DEM to assist in removing known terrain phase, the actual differential interferometric phase is simplified to a mixture of local terrain residual phase (corresponding to actual elevation updates), canopy scattering phase components, orbital errors, and atmospheric delay. Since suboptimal spatiotemporal baseline configurations easily induce interferometric decoupling, and atmospheric changes and satellite orbital drift often introduce elevation errors on the order of hundreds of meters into the interferograms, effectively decoupling the orbital error phase, atmospheric delay phase, and local terrain residual phase in space constitutes a current technical bottleneck for high-precision mapping using long-wavelength repeated orbit spaceborne SAR.

[0004] Currently, the mainstream methods for error separation in interferograms typically employ a step-by-step processing strategy: estimating orbital errors through low-order polynomial fitting, followed by correcting atmospheric delay errors. However, in dense forest penetration scenarios and sparse temporal sampling (such as the BIOMASS mission), existing methods have significant drawbacks: First, traditional polynomial orbital fitting is usually applied directly to a single pair of interferograms, ignoring the interference of spatially inhomogeneous atmospheric and topographic residuals on trend fitting, which can easily lead to mathematical overfitting and phase component coupling; Second, limited by the low spatiotemporal resolution of global atmospheric models and water vapor products, model correction based on external datasets is difficult to meet the needs of high-precision DEM generation; Third, in homogeneous vegetation areas such as tropical rainforests, the phase of topographic residuals caused by differences in forest penetration depth often exhibits a trend distribution in space, making it difficult for simple spatial scale or spectral separation methods to distinguish atmospheric delay from trend-based topographic residuals. Summary of the Invention

[0005] This invention provides a method, apparatus, device, and medium for long-wave SAR atmospheric error correction, which can effectively separate atmospheric delay error, orbital error, and local terrain residual under sparse time sampling.

[0006] To achieve the above technical objectives, the present invention adopts the following technical solution: A long-wave SAR atmospheric error correction method includes: A multi-temporal dual-polarization SAR image dataset of the study area was obtained, and differential interferometry was performed by introducing an external reference DEM to generate a multi-temporal dual-polarization differential interferogram. For the multi-temporal differential interferograms of each dual-polarization channel, the initial coarse topographic residuals are estimated by weighting the spatial baseline and interferometric coherence. Furthermore, the original coarse topographic residuals are estimated by fitting a polynomial orbital error to the bare land area. The composite residual phase is calculated based on the differential interferometric phase of a single interferometric pair and the coarse topographic residual after removing orbital errors. Dual-polar wavelet correlation analysis was performed on the composite residual phase of the dual-polarized channel to identify the common-mode error signal composed of atmospheric delay and orbital error, and it was separated from the composite residual phase to obtain the high-frequency terrain residual compensation phase. The high-frequency terrain residual compensation phase is transformed to the elevation domain and the coarse terrain residual is refined. It is then superimposed on the external reference DEM to generate the final high-precision dual-channel DEM.

[0007] Furthermore, the dual polarization includes co-polarization HH and cross-polarization HV, which are respectively introduced into different external reference DEMs for differential interferometry processing, specifically: The Copernicus DEM, which integrates X-band and optical elevation data, is used as the external reference DEM for the HH polarization channel. The SRTM DEM derived from C-band interferometry was used as the external reference DEM for the HV polarization channel.

[0008] Furthermore, the method of using spatial baseline and interferometric coherence as weights to perform weighted estimation of the initial coarse topographic residuals for multi-temporal dual-polarization differential interferograms specifically includes: Construct the following weighted least squares objective function : (1); in, For the number of multi-temporal interference pairs, The initial coarse topographic residual is to be determined; For the index of the interference pair, This is the filtered differential interference phase. Vertical wavenumber: ; in The vertical spatial baseline length, For radar operating wavelength, This is the slant distance from the radar sensor to the ground surface. The angle of incidence of radar electromagnetic waves; For the first Coherence of each interference pair The determined weights are taken as follows: for pixels with coherence higher than a preset value, the weight is... = For pixels with coherence less than a preset value, take... =0; Let the objective function Unknown initial coarse terrain residuals The partial derivatives are zero, and the joint estimation yields the initial coarse topographic residuals: .

[0009] Furthermore, the step of using bare land areas to perform polynomial orbital error fitting, and then estimating the coarse topographic residuals after removing the orbital errors, specifically includes: Based on land cover classification products, extract non-forest cover bare land areas from the images; Based on the initial coarse topographic residuals of the bare land area, a systematic track tilt error is fitted. : ; in, These are the azimuth and range coordinates of the image pixels, respectively. These are the fitting coefficients; From the initial coarse terrain residual Subtract track error This yields the coarse terrain residuals after removing track errors. : .

[0010] Furthermore, the composite residual phase is calculated based on the differential interferometric phase of a single interferometric pair and the coarse topographic residual after removing orbital errors, specifically expressed as follows: ; in, For the first The composite residual phase of each interferometric pair is composed of low-frequency atmospheric delay and orbital error independent of polarization, and high-frequency local topographic residual components that are sensitive to polarization. This is the filtered differential interference phase. Vertical wavenumber; To remove coarse terrain residuals from track errors.

[0011] Furthermore, the dual-polarization wavelet correlation analysis of the composite residual phase of the dual-polarization channel specifically includes: Step a1, firstly utilize the composite residual phase of the HH polarization channel Perform optimal wavelet decomposition scale search: Two-dimensional discrete wavelet decomposition is performed on it using a series of candidate scales, and the decomposition performance index of the low-frequency reconstructed signal relative to the original composite residual phase is calculated at each candidate scale, including root mean square error, signal-to-noise ratio and smoothness. Step a2: Calculate the absolute values ​​of the first-order differences of the three performance indicators between adjacent decomposition scales, and normalize them to obtain the normalized change of the root mean square error. Normalized change in signal-to-noise ratio and the normalized change in smoothness ; Step a3: Use the information entropy weighting method to quantify the uncertainty of each normalized change quantity, and determine... , , Corresponding weights , , This leads to the construction of a weighted mixed evaluation index. : ; Using robust nonlinear least squares method based on median absolute deviation to evaluate mixed evaluation indices With scale The added variation curve is fitted and smoothed, and the inflection point of the smoothed curve is adaptively detected by the second derivative. The scale corresponding to the inflection point is determined as the unified optimal wavelet decomposition scale. Step a4: Based on the determined optimal wavelet decomposition scale, perform complete two-dimensional discrete wavelet decomposition on the composite residual phase of the dual-polarization channel. Step a5: For each decomposition level, extract the high-frequency coefficients of each polarization channel in the horizontal, vertical, and diagonal directions; and in the above three directions, extract the composite residual phase of the HV polarization channel. The high-frequency coefficients and the composite residual phase of the HH polarization channel The high-frequency coefficients were constructed into two-dimensional scattered coordinates, and a straight line was independently fitted for each. Then, the perpendicular distance from each two-dimensional scattered coordinate to the corresponding fitted line was calculated. And construct a distance weighting factor based on Gaussian decay: ;

[0012] Step a6: Based on the fitted straight lines in each direction, the high-frequency coefficients of the HV polarization channel are mapped to the HH polarization channel, and the distance weighting factor is used. The high-frequency coefficients after mapping are weighted to obtain the high-frequency common-mode components; Step a7: The low-frequency coefficients of the HH polarization channel obtained by two-dimensional wavelet decomposition are spliced ​​with the high-frequency common-mode components, and the complete polarization-independent common-mode interference phase is reconstructed by inverse wavelet transform, which is the common-mode error signal between the dual polarization channels.

[0013] Furthermore, the step of converting the high-frequency terrain residual phase to the elevation domain and refining the coarse terrain residual, then superimposing it on an external reference DEM to generate a final high-precision dual-channel DEM, specifically includes: Based on the coarse terrain residuals after removing track errors Reconstruct a two-channel coarse DEM using an external reference DEM: ; in, For external reference DEM, The coarse DEM obtained from the reconstruction of each channel; High-frequency terrain residual compensation phase Converted to the elevation domain and overlaid on the coarse DEM of each channel, the final high-precision DEM of the two channels is obtained: ; in, For the first The interferometric pair ultimately yields a high-precision DEM. For the first Vertical wavenumber of each interference pair.

[0014] A long-wave SAR atmospheric error correction device, comprising: The preprocessing module is used to: acquire a multi-temporal dual-polarization SAR image dataset of the study area, and generate a multi-temporal dual-polarization differential interferogram by introducing an external reference DEM for differential interferometry processing. The weighted and orbit correction module is used to: use spatial baseline and interferometric coherence as weights to perform weighted estimation of initial coarse topographic residuals for multi-temporal dual-polarization differential interferograms; and use bare land areas to perform polynomial orbit error fitting, thereby estimating coarse topographic residuals after removing orbit errors. The multi-scale phase decoupling module is used to: perform dual-polarization wavelet correlation analysis on the composite residual phase of the dual-polarization channel, identify the common-mode error signal composed of atmospheric delay and orbital error, and separate it from the composite residual phase to obtain the high-frequency terrain residual compensation phase; The elevation compensation module is used to: convert the high-frequency terrain residual compensation phase to the elevation domain, refine the coarse terrain residual, overlay it on the external reference DEM, and generate and output the final high-precision DEM with dual channels.

[0015] An electronic device includes a memory and a processor. The memory stores a computer program, which, when executed by the processor, causes the processor to implement the long-wave SAR atmospheric error correction method described above.

[0016] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the long-wave SAR atmospheric error correction method described above.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: First, by introducing different reference DEMs (e.g., using the Copernicus DEM for the HH channel and the SRTM DEM for the HV channel), the initial phase heterogeneity between the dual-polarization channels is actively amplified. Secondly, a baseline-coherence weighted stacking strategy was adopted to perform orbit fitting in conjunction with bare land areas. Stacking weakened the influence of the atmosphere on orbit error fitting, and the extraction of bare land locations weakened the bias of vegetation height on orbit error fitting. This overcame the problems of orbital surface overfitting and elevation error aliasing in traditional single interferometric pair processing and robustly removed systematic errors. Finally, by leveraging the physical mechanism differences in penetration depth between P-band HH polarization (dominantly tree trunk-surface secondary scattering, with a deeper phase center) and HV polarization (dominantly canopy scattering, with a shallower phase center), the optimal wavelet scale for dual-polarization wavelet correlation analysis (DP-WCA) in the spatial frequency domain was determined. This adaptively separated polarization-independent common-mode low-frequency errors (atmospheric and orbital errors) from polarization-related high-frequency real terrain.

[0018] Experiments conducted in dense tropical rainforest regions demonstrate that this method can significantly refine the root mean square error of absolute elevation (RMSE) of the initial DEM from approximately 26 meters to the order of 8.65 meters, breaking through the technical bottleneck of phase coupling in topographic mapping using long-wave heavy orbit SAR systems. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating the generation of a high-precision DEM based on a multi-temporal dual-polarization processing framework according to an embodiment of the present invention.

[0020] Figure 2 This is a schematic diagram showing the location and dataset coverage of a national park test site in an embodiment of the present invention.

[0021] Figure 3This is a comparison of the effects of applying baseline-coherence weighted stacking and joint orbit correction to generate coarse terrain residuals in an embodiment of the present invention; subplot (a) is the coarse terrain residual obtained by the 0123_0126 interferometric pair, subplot (b) is the coarse terrain residual after weighted stacking and orbit correction, and subplot (c) is the DEM error statistical histogram obtained by various methods.

[0022] Figure 4 This is a comparison diagram of the overfitting caused by traditional single-pair orbit fitting in the embodiments of the present invention and the one-dimensional elevation profile of the weighted stacking strategy of the present invention. Detailed Implementation

[0023] The technical solution of the present invention will be further explained below with reference to the accompanying drawings and embodiments. The embodiments are described in detail using a typical tropical rainforest core test field (covered with dense and structurally complex high biomass forest) and BIOMASS Level-1A single-view multiple (SLC) imagery as examples.

[0024] This embodiment provides a long-wave SAR atmospheric error correction method based on weighted stacking and wavelet correlation analysis. The overall processing flow is described in [link to documentation]. Figure 1 As shown, it mainly includes four steps: Step 1: Acquire multi-temporal dual-polarization SAR images and generate initial differential interferograms.

[0025] This embodiment first obtains the coverage study area (location and data coverage details are available in [link]). Figure 2 The multi-temporal fully polarimetric BIOMASS SLC interferometer pair, with spatiotemporal baselines controlled within a 3-day revisit period over the time series span, corresponds to a vertical wavenumber range of... The average coherence ranges from 0.45 to 0.62. The multi-temporal dual-polarization dataset extracted includes co-polarization (HH) and cross-polarization (HV) channels. The interferometric pair summaries of the BIOMASS dataset used are shown in Table 1. Table 1 shows the BIOMASS dataset information used.

[0026] The naming rules in the "Interference Pair Number" column indicate the acquisition date (MMDD 2026) of the main image (left) and the secondary image (right).

[0027] To actively introduce and amplify the initial phase heterogeneity between the dual-polarization channels, this embodiment introduces two types of external reference elevation sources with characteristic differences for differential interferometry processing: For the HH channel, a Copernicus DEM that integrates high-precision X-band and optical remote sensing elevation data is introduced as an external reference DEM. ; For the HV channel, a traditional SRTM DEM derived from the C-band is used as an external reference DEM. Its penetration depth is relatively deeper than that of the Copernicus DEM.

[0028] In forest scenarios, due to differences in radar bands and sensor types, the terrain recorded by these two DEMs exhibits spatial distribution differences. This embodiment introduces a reference DEM with different acquisition times and band penetration characteristics to actively amplify the characteristic differences between the initial differential interference phases of the HH and HV channels.

[0029] The process involves using interferometric processing software for precise registration of radar images, differential multiplication to remove flat terrain effects and known terrain phases, and finally adaptive filtering and phase unwrapping. For any given... An interference pair, its unwrapped differential interference phase The model is as follows: ; in, , , These represent the orbital error phase, atmospheric delay phase, and noise phase, respectively. This is the topographic residual phase, which is updated with the true surface elevation of the relative reference DEM. The physical mapping relationship is as follows: ; Vertical wavenumber Defined as: ; In equation (3), For the first Vertical spatial baseline length of each interference pair The operating wavelength of the radar (the operating wavelength of the BIOMASS P-band is approximately...) ), This is the slant distance from the radar sensor to the ground surface. This is the incident angle of the radar electromagnetic wave.

[0030] After step 1 is completed, the generated initial differential interferometric phase contains the real terrain residual information coupled with the unknown orbital error and atmospheric delay error.

[0031] Step 2: Coarse terrain residual estimation and trajectory correction based on weighted stacking.

[0032] After obtaining the multi-temporal differential interferometric phases of each of the two channels, directly performing polynomial orbit fitting will be significantly affected by atmospheric interference, leading to fitting surface distortion and overfitting. Therefore, this embodiment introduces a spatial baseline (which affects the vertical wavenumber) as a reference. and interference coherence The joint weighted time series stacking strategy aims to effectively reduce the influence of the atmosphere on orbital error estimation by using time series averaging.

[0033] For the unwrapped filtered phase, construct the following weighted least squares objective function. In the time dimension, the uncorrelated random atmospheric delay component is suppressed in a least-squares sense: ; In equation (4), The total number of multi-temporal interference pairs participating in the solution. The weights are determined by coherence. For pixels with coherence higher than a preset value (0.3 in this embodiment), a weight is assigned. = For pixels with coherence less than a preset value, take... =0. The weighting factor is designed to be consistent with the coherence of the interference pair. This correlation ensures that interferometric pairs with high coherence in time series and long spatial baselines have higher confidence in terrain estimation. Let the objective function... Elevation residuals for unknown terrain The partial derivatives are zero, and the unbiased initial coarse topographic residuals are derived and calculated: ; Although weighted stacking suppresses random atmospheric delay in the time dimension, the obtained Systematic orbital tilt caused by baseline slope estimation errors still remains. To address this issue and mitigate the impact of vegetation canopy height and volume scattering on orbital surface calculation, this embodiment combines Sentinel-2 land cover classification data to remove forests, water bodies, and buildings across the entire image, extracting bare land areas with near-zero canopy cover. On these extracted bare land areas, a low-order two-dimensional polynomial is used to fit the orbital tilt error surface. : ; In equation (6), These are the azimuth and range coordinate pixel indices corresponding to the image, respectively. These are the surface coefficients obtained through least-multiplication fitting. Subtracting the track inclination surface from the initial coarse elevation residuals yields the coarse topographic residuals, which have been adjusted to remove systematic track errors. : ; Combined with the external reference DEM introduced in step 1 Calculate and reconstruct a coarse digital elevation model : ; After processing, systematic tilt errors were effectively removed (see comparison of coarse terrain residuals before and after weighted stacking and track correction). Figure 3 Since this area is mostly covered by vegetation, only the area to the right of the central river consists of relatively shallow grassland and bare land. Figure 3 The distribution of negative signals at the top and bottom and positive signals in the middle in (a) indicates overfitting during single-scene track error removal. However, through step 2, Figure 3 In (b), the topographic residuals show a uniform penetration depth (P-band has stronger penetration capability), and there are no obvious residual errors in grassland and bare land areas. Furthermore, Figure 3 The error statistics histogram shown in (c) also demonstrates its superiority. Therefore, selecting control points in bare land areas to estimate orbital errors can overcome the influence of the difference in forest penetration depth (up to 30m) between long-wave SAR and the reference DEM on orbital error estimation, resulting in more accurate coarse terrain residuals. Thus, compared to conventional single-scene orbital error removal, this method can better mitigate the influence of atmospheric delay errors and more closely approximate the actual terrain. Figure 4 ).

[0034] Step 3: Adaptive wavelet phase decoupling analysis.

[0035] The coarse topographic residual removes systematic orbital errors. Although multi-phase stacking has been performed, some atmospheric delay error remains. This embodiment returns to the single interferometric processing level for fine processing.

[0036] First, the coarse terrain residuals obtained in step 2 are... By converting back to the phase domain using formula (2), the phase is subtracted from the original unwrapped interference phase of a single interference pair to calculate the composite residual phase containing only noise and error. : ; The residual The problem involves the superposition of spatially correlated low-frequency atmospheric delay errors, orbital errors, and high-frequency topographic residuals. In this step, this embodiment leverages the heterogeneity of forest scattering physics under different polarization channels at long wavelengths (P-band). HH polarization exhibits extremely strong forest penetration depth, with dominant scattering originating from secondary scattering between the tree trunk and the ground surface; its phase center height is close to the actual ground surface. In contrast, HV polarization experiences strong volume scattering in the canopy and treetop layers, resulting in a higher phase center height, farther from the actual ground surface. Based on this high spatial heterogeneity, this embodiment introduces a dual-polarization wavelet correlation analysis (DP-WCA) algorithm in the space-frequency domain. By performing multi-resolution discrete wavelet transform (DWT) on the dual-channel composite residual phase, it decomposes it into low-frequency and high-frequency coefficients. For the composite residual phase... and Its wavelet multiresolution decomposition formula is expressed as: ; ; In equations (10) and (11), and These are the wavelet low-frequency (LF) and high-frequency (HF) coefficients obtained from the decomposition, respectively; and These are the smoothing scaling function and the mother wavelet function, respectively. The wavelet decomposition scale; and The window size for wavelet decomposition; These represent high-frequency detail signals from the horizontal, vertical, and diagonal directions, respectively. Represents spatial pixel coordinates.

[0037] Since the phase errors caused by atmospheric delay and orbital deviation are polarization-independent components, they exhibit a trend-like characteristic in space. Therefore, in the low-frequency wavelet coefficients, HH polarization and HV polarization show extremely high common-mode correlation, requiring no additional processing during subsequent reconstruction. However, some errors still exist in the high-frequency coefficients. Therefore, this embodiment further performs correlation analysis on the high-frequency coefficients, extracting the common components through weighted extraction, and then reconstructing the weighted high-frequency coefficients with the unprocessed low-frequency coefficients to obtain the complete polarization-independent common-mode interference phase, which is the common-mode error signal between the dual-polarization channels. To accurately remove the interference of the aforementioned common-mode error signal, this embodiment constructs a selection strategy that integrates the variation characteristics of multiple statistical evaluation indicators, adaptively and robustly determining the optimal wavelet decomposition scale for wavelet correlation analysis in the space-frequency domain. Specifically, different wavelet decomposition scales are calculated respectively. Rate of change of root mean square error (RMSE) Signal-to-noise ratio (SNR) change rate and the rate of change of smoothness The information entropy weighting method is used to quantify the uncertainty of each indicator and determine its corresponding weight. , , Thus, a weighted mixed evaluation index is constructed. : ; Subsequently, based on this hybrid evaluation index With scale The added change curve is used to adaptively detect the inflection point of the curve using the Huber exponent, and the scale corresponding to the inflection point is determined as the optimal wavelet decomposition scale.

[0038] At the optimal wavelet decomposition scale, this embodiment proposes a high- and low-frequency hierarchical reconstruction strategy to adaptively separate polarization-independent common-mode interference phases. First, for the low-frequency components, given that atmospheric delay and residual orbital errors dominate the low-frequency spatial signal, this embodiment treats the wavelet low-frequency coefficients of the HH polarization channel as pure atmospheric and orbital errors. Second, for the high-frequency detail components at each layer, since they are mixed with polarization-sensitive local topographic residuals and polarization-independent high-frequency noise, this embodiment introduces multi-polarization spatial orthogonal distance regression analysis layer by layer. Specifically, it calculates the linear fitting relationship of the high-frequency coefficients of the dual-polarization channel in the horizontal, vertical, and diagonal directions, and calculates the vertical distance from the auxiliary polarization coefficient mapping point to the fitted line. And construct a distance weighting factor based on Gaussian decay: ; The fitted high-frequency coefficients are weighted using this weighting factor to accurately extract highly correlated high-frequency common-mode noise components in the dual-polarization channel. Finally, the low-frequency coefficients and the weighted high-frequency coefficients from each layer are arrayed and stitched together. The complete polarization-independent common-mode interference phase is then reconstructed using inverse wavelet transform (IWT). .

[0039] High-frequency terrain information exhibits polarization sensitivity in the high-frequency wavelet region due to the difference in the penetration phase center heights of HH and HV. By subtracting the separated common-mode low-frequency interference phase from the original single-interferometric pair composite residual phase, adaptive decoupling of atmospheric-orbit errors and terrain residuals can be achieved, thereby extracting the high-frequency terrain residual compensation phase. : ; It should be noted that the stability of DP-WCA high-frequency decoupling is physically constrained by coherence. When severe temporal decoherence occurs in local regions of the interferogram, it can lead to reduced coherence. When the wavelet coefficients fall below the theoretical physical lower bound, they will be overwhelmed by unwrapping noise. In this embodiment, a low-coherence mask or a reduction in the high-frequency compensation weight is applied to the decoherent region to prevent the introduction of synthesized high-frequency noise.

[0040] Step 4: Elevation refinement compensation and final DEM generation.

[0041] Decoupling to obtain high-frequency terrain residual compensation phase Subsequently, this embodiment uses vertical wavenumber This is then converted to the elevation domain. The refined elevation compensation is used to refine the coarse topographic residuals, and then superimposed on the reference DEM to obtain the final decoupled and refined high-precision digital elevation model. : ; In this embodiment, within a dense tropical rainforest region, the absolute accuracy of each elevation model retrieved in this invention was evaluated and verified using synchronously acquired airborne high-precision LiDAR DTM (real forest understory topography). Taking the 0123_0126 interferometric pair as an example, the quantitative evaluation results at each stage of the experiment are shown in Table 2: Table 2. Quantitative Accuracy Assessment of Digital Elevation Models (DEMs) at Each Stage Before and After Terrain Residual Compensation

[0042] Table 2 and related analysis results show that the initial input Copernicus DEM and SRTM DEM have a large positive elevation deviation relative to the actual forest understory due to band and canopy reflection limitations. After baseline-coherence weighted stacking and bare ground fitting in step 2, the overall systematic orbital error and most of the random atmospheric error are removed, and the RMSE of the generated coarse DEM is significantly reduced. The overall terrain trend returned to normal. Further steps, including step 3 using the DP-WCA algorithm to decouple low-frequency atmospheric remnants from a single interferometer pair in the space-frequency domain, and step 4 for high-frequency terrain detail refinement compensation, restored the detailed features and micro-topography of the DEM, achieving a final RMSE accuracy refined to [value missing]. .

[0043] Experimental data quantitatively verified the superior performance and industrial mapping value of this invention in high-precision phase decoupling and terrain restoration in long-wavelength repetitive orbit interferometry phase processing.

[0044] The above embodiments are preferred embodiments of this application. Those skilled in the art can make various changes or improvements based on them. Without departing from the overall concept of this application, these changes or improvements should fall within the scope of protection claimed in this application.

Claims

1. A method for correcting atmospheric errors in long-wave SAR, characterized in that, include: A multi-temporal dual-polarization SAR image dataset of the study area was obtained, and differential interferometry was performed by introducing an external reference DEM to generate a multi-temporal dual-polarization differential interferogram. For the multi-temporal differential interferograms of each dual-polarization channel, the spatial baseline and interferometric coherence are used as weights to estimate the initial coarse topographic residuals. Furthermore, polynomial orbital error fitting was performed using bare land areas to estimate the coarse topographic residuals after removing orbital errors; The composite residual phase is calculated based on the differential interferometric phase of a single interferometric pair and the coarse topographic residual after removing orbital errors. Dual-polar wavelet correlation analysis was performed on the composite residual phase of the dual-polarized channel to identify the common-mode error signal composed of atmospheric delay and orbital error, and it was separated from the composite residual phase to obtain the high-frequency terrain residual compensation phase. The high-frequency terrain residual compensation phase is transformed to the elevation domain and the coarse terrain residual is refined. It is then superimposed on the external reference DEM to generate the final high-precision dual-channel DEM.

2. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, The dual polarization includes co-polarization HH and cross-polarization HV, which are respectively introduced into different external reference DEMs for differential interferometry processing. Specifically: The Copernicus DEM, which integrates X-band and optical elevation data, is used as the external reference DEM for the HH polarization channel. The SRTM DEM derived from C-band interferometry was used as the external reference DEM for the HV polarization channel.

3. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, Using spatial baseline and interferometric coherence as weights, the initial coarse topographic residuals are estimated by weighting multi-temporal dual-polarization differential interferograms, specifically including: Construct the following weighted least squares objective function : （1）； in, For the number of multi-temporal interference pairs, The initial coarse topographic residual is to be determined; For the index of the interference pair, This is the filtered differential interference phase. Vertical wavenumber: ； in The vertical spatial baseline length, For radar operating wavelength, This is the slant distance from the radar sensor to the ground surface. The angle of incidence of radar electromagnetic waves; For the first Coherence of each interference pair The determined weights are taken as follows: for pixels with coherence higher than a preset value, the weight is... = For pixels with coherence less than a preset value, take... =0; Let the objective function Unknown initial coarse terrain residuals The partial derivatives are zero, and the joint estimation yields the initial coarse topographic residuals: 。 4. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, Polynomial orbital error fitting is performed using bare land areas to estimate coarse topographic residuals after removing orbital errors. Specifically, this includes: Based on land cover classification products, extract non-forest cover bare land areas from the images; Based on the initial coarse topographic residuals of the bare land area, a systematic track tilt error is fitted. : ； in, These are the azimuth and range coordinates of the image pixels, respectively. These are the fitting coefficients; From the initial coarse terrain residual Subtract track error This yields the coarse terrain residuals after removing track errors. : 。 5. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, The composite residual phase is calculated based on the differential interferometric phase of a single interferometric pair and the coarse topographic residual after removing orbital errors, and is specifically expressed as follows: ； in, For the first The composite residual phase of each interferometric pair is composed of low-frequency atmospheric delay and orbital error independent of polarization, and high-frequency local topographic residual components that are sensitive to polarization. This is the filtered differential interference phase. Vertical wavenumber; To remove coarse terrain residuals from track errors.

6. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, The dual-polarization wavelet correlation analysis of the composite residual phase of the dual-polarization channel specifically includes: Step a1, firstly utilize the composite residual phase of the HH polarization channel Perform optimal wavelet decomposition scale search: Two-dimensional discrete wavelet decomposition is performed on it using a series of candidate scales, and the decomposition performance index of the low-frequency reconstructed signal relative to the original composite residual phase is calculated at each candidate scale, including root mean square error, signal-to-noise ratio and smoothness. Step a2: Calculate the absolute values ​​of the first-order differences of the three performance indicators between adjacent decomposition scales, and normalize them to obtain the normalized change of the root mean square error. Normalized change in signal-to-noise ratio and the normalized change in smoothness ; Step a3: Use the information entropy weighting method to quantify the uncertainty of each normalized change quantity, and determine... , , Corresponding weights , , This leads to the construction of a weighted mixed evaluation index. : ； Using robust nonlinear least squares method based on median absolute deviation to evaluate mixed evaluation indices With scale The added variation curve is fitted and smoothed, and the inflection point of the smoothed curve is adaptively detected by the second derivative. The scale corresponding to the inflection point is determined as the unified optimal wavelet decomposition scale. Step a4: Based on the determined optimal wavelet decomposition scale, perform complete two-dimensional discrete wavelet decomposition on the composite residual phase of the dual-polarization channel. Step a5: For each decomposition level, extract the high-frequency coefficients of each polarization channel in the horizontal, vertical, and diagonal directions; and in the above three directions, extract the composite residual phase of the HV polarization channel. The high-frequency coefficients and the composite residual phase of the HH polarization channel The high-frequency coefficients were constructed into two-dimensional scattered coordinates, and a straight line was independently fitted for each. Then, the perpendicular distance from each two-dimensional scattered coordinate to the corresponding fitted line was calculated. And construct a distance weighting factor based on Gaussian decay: ； Step a6: Based on the fitted straight lines in each direction, the high-frequency coefficients of the HV polarization channel are mapped to the HH polarization channel, and the distance weighting factor is used. The high-frequency coefficients after mapping are weighted to obtain the high-frequency common-mode components; Step a7: The low-frequency coefficients of the HH polarization channel obtained by two-dimensional wavelet decomposition are spliced ​​with the high-frequency common-mode components, and the complete polarization-independent common-mode interference phase is reconstructed by inverse wavelet transform, which is the common-mode error signal between the dual polarization channels.

7. The long-wave SAR atmospheric error correction method according to claim 1, characterized in that, The process of converting the high-frequency terrain residual phase to the elevation domain and refining the coarse terrain residual, then superimposing it on an external reference DEM to generate a final high-precision dual-channel DEM, specifically includes: Based on the coarse terrain residuals after removing track errors Reconstruct a two-channel coarse DEM using an external reference DEM: ； in, For external reference DEM, The coarse DEM obtained from the reconstruction of each channel; High-frequency terrain residual compensation phase Converted to the elevation domain and overlaid on the coarse DEM of each channel, the final high-precision DEM of the two channels is obtained: ； in, For the first The interferometric pair ultimately yields a high-precision DEM. For the first Vertical wavenumber of each interference pair.

8. A long-wave SAR atmospheric error correction device, characterized in that, include: The preprocessing module is used to: acquire a multi-temporal dual-polarization SAR image dataset of the study area, and generate a multi-temporal dual-polarization differential interferogram by introducing an external reference DEM for differential interferometry processing. The weighted and orbit correction module is used to: use spatial baseline and interferometric coherence as weights to perform weighted estimation of initial coarse topographic residuals for multi-temporal dual-polarization differential interferograms; and use bare land areas to perform polynomial orbit error fitting, thereby estimating coarse topographic residuals after removing orbit errors. The multi-scale phase decoupling module is used to: perform dual-polarization wavelet correlation analysis on the composite residual phase of the dual-polarization channel, identify the common-mode error signal composed of atmospheric delay and orbital error, and separate it from the composite residual phase to obtain the high-frequency terrain residual compensation phase; The elevation compensation module is used to: convert the high-frequency terrain residual compensation phase to the elevation domain, refine the coarse terrain residual, overlay it on the external reference DEM, and generate and output the final high-precision DEM with dual channels.

9. An electronic device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the computer program is executed by the processor, the processor implements the long-wave SAR atmospheric error correction method as described in any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the long-wave SAR atmospheric error correction method as described in any one of claims 1 to 7.

Images