High-precision inversion method for tropospheric refractive index based on spaceborne spotlight SAR interferometry

Abstract

The application discloses a tropospheric atmospheric refractive index high-precision inversion method based on spaceborne spotlight SAR interference, and solves the resolution mismatch problem caused by unreasonable observation geometry configuration and viewing angle difference in the method by constructing a multi-target optimization mechanism and a resolution consistency constraint framework. Compared with the existing technology which relies on artificial experience setting or adopts a simple and uniform division method, equivalent multi-angle observation conditions are constructed by using spaceborne SAR spotlight imaging data through data processing, three-dimensional inversion of the tropospheric atmospheric refractive index is realized under the premise of not increasing the scale and system complexity of the constellation, and the long synthetic aperture information formed by spotlight imaging is effectively introduced into the three-dimensional tomographic inversion process through optimal sub-aperture division and coherence constraint, so that the inversion accuracy and stability are improved from the signal level, and high-precision three-dimensional reconstruction of the tropospheric atmospheric refractive index is realized, which can provide reliable atmospheric structure information support for fine monitoring and prediction of severe convective weather disasters.

Description

Technical Field

[0001] This application relates to the field of data inversion technology, and in particular to a high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry. Background Technology

[0002] In my country's natural disaster system, severe convective meteorological disasters such as rainstorms, typhoons, blizzards, and strong winds are characterized by high frequency of occurrence, wide impact range, and strong destructiveness, accounting for approximately two-thirds of all natural disasters. They cause an average of over 200 billion yuan in direct economic losses annually, seriously threatening people's lives and property and the stable operation of the economy and society. Therefore, conducting high-precision, forward-looking risk monitoring, early warning, and forecasting for severe convective meteorological disasters is of great practical significance and an urgent need.

[0003] The core of disaster prevention lies in "sensing first, then predicting," meaning that only by achieving accurate perception of atmospheric conditions can reliable predictive data and scientific decision-making support be provided for the occurrence, development, and evolution of disasters. However, constrained by observation systems and technological conditions, existing meteorological observation methods cannot simultaneously meet the needs for wide-area coverage, high spatial resolution, and all-day, all-weather, and detailed perception of atmospheric elements.

[0004] Traditional meteorological observation methods, such as ground-based automatic weather stations and radiosonde observations, have limited spatial coverage and are unable to continuously depict the spatiotemporal evolution of large-scale weather systems. Optical and infrared remote sensing methods are easily affected by cloud and rain obstruction, significantly limiting their observation capabilities under severe convective weather conditions. In contrast, spaceborne Synthetic Aperture Radar (SAR) offers advantages in all-weather, all-time imaging and, combined with interferometry techniques, has demonstrated good application potential in the inversion of meteorological elements such as precipitable water over large areas, providing a new technical approach for atmospheric perception under complex meteorological conditions.

[0005] However, existing SAR-based meteorological inversion methods mostly focus on two-dimensional parameter inversion, which is insufficient to effectively characterize the vertical layering structure of the atmosphere and its three-dimensional evolution characteristics, and thus cannot meet the higher demands for fine structural characterization and accurate forecasting of severe convective systems. Therefore, three-dimensional reconstruction of meteorological elements has become a key breakthrough direction for improving disaster monitoring and forecasting capabilities.

[0006] From the perspective of observation mechanism, the acquisition of three-dimensional atmospheric information usually relies on multi-angle observations. Existing multi-angle observations typically rely on multi-orbit, multi-satellite observations to construct multi-view information, which is costly in terms of system configuration, task scheduling, and data acquisition. Spotlight SAR, by continuously pointing the radar beam at the same ground area, can achieve high azimuth resolution imaging and has the characteristics of controllable beam pointing and adjustable observation angle. Existing technologies mostly regard spotlight SAR as a high-resolution two-dimensional imaging method, and have not fully explored the observational geometric diversity implied by the long synthetic aperture formed during spotlight imaging, making it difficult to construct a multi-angle observation structure that meets the requirements of three-dimensional tomography inversion under single-satellite, single-observation conditions. Summary of the Invention

[0007] The main objective of this application is to provide a high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry, in order to solve the problems mentioned in the background art.

[0008] To achieve the above objectives, this application provides the following technical solution: A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry is characterized by the following specific steps: S1. Data Acquisition and Purification: Acquire raw echo data from the spaceborne focused synthetic aperture radar in the target observation area, and load digital elevation model and precise orbital ephemeris data; sequentially perform deskewing, coarse focusing and Doppler parameter estimation on the raw echo data to obtain basic imaging data with a unified spatiotemporal reference. S2. Phase purification and isolation: Calculate the topographic phase component, flat terrain phase component, and ionospheric error phase component according to the imaging geometry, and remove these phase components from the interferogram; remove noise through a filtering algorithm and extract the residual interferometric phase that characterizes the tropospheric atmospheric delay. S3. Accuracy Model Construction: Based on the geometric characteristics of the long synthetic aperture in the clustering mode, the observation matrix expression for the tropospheric atmospheric refractive index tomography inversion is established; the posterior covariance matrix and tomographic grid coverage function are derived to construct a quantitative evaluation model for the inversion accuracy. S4. Geometric structure optimization: Establish an optimization function with the dual objectives of maximizing the coverage of the tomographic grid and minimizing the theoretical inversion error; use a non-dominated sorting genetic algorithm to iteratively divide the long synthetic aperture, and obtain the optimal number of sub-apertures and the corresponding center oblique angle values ​​of each sub-aperture under the condition of meeting the preset accuracy threshold. S5. Consistency Constraint: Based on the optimal sub-aperture center position obtained in step S4, calculate the distance from the center of each sub-aperture to the two sides of the overall composite aperture, and take the minimum value among all distances as the uniform sub-aperture length so that the imaging results of each sub-aperture maintain the same resolution in the azimuth direction. S6. Parallel Imaging Interferometry: Imaging processing is performed on the data of each sub-aperture separately to obtain complex images corresponding to each sub-aperture; the complex image corresponding to the central viewpoint is selected as the registration reference, and sub-pixel level registration and resampling are performed on the complex images of the remaining sub-apertures to obtain multiple sets of interferograms; S7. Construction of tomographic equations: The target atmospheric space is divided into a three-dimensional voxel grid. The path length of the radar beam through each voxel is calculated based on the orbital parameters of each sub-aperture to construct the observation matrix. A three-dimensional tomographic linear equation set between the multi-angle observation phase and atmospheric refractive index parameters is established. S8. Sparse Regularized Inversion: Construct an objective function containing data fitting terms and sparse regularization terms, and use an iterative threshold shrinkage algorithm to solve the objective function to obtain a stable solution for the atmospheric refractive index parameter. S9. Three-dimensional field reconstruction: The atmospheric refractive index parameters obtained in step S8 are mapped to a three-dimensional spatial grid. After spatial interpolation and smoothing, the three-dimensional spatial distribution data and vertical profile data of the tropospheric atmospheric refractive index are output.

[0009] Preferably, in step S1, the specific method is as follows: S1.1 Acquire the raw echo data of the satellite-borne focused synthetic aperture radar in the target observation area, simultaneously load the digital elevation model data and satellite precise orbit ephemeris data covering the observation area, unify the spatial coordinate system of the digital elevation model with the radar imaging coordinate system, establish the spatial correspondence of multi-source data, and provide terrain reference and orbital geometric parameters for subsequent phase separation. S1.2. The original echo data is sequentially deskewing to reduce the data sampling rate, coarse focusing to compress the range and azimuth signals, and the Doppler center frequency and Doppler modulation frequency parameters are estimated. The original echo data is then converted into complex imaging data with a unified spatiotemporal reference. This complex imaging data retains complete phase information and provides a data basis for sub-aperture segmentation and interferometric processing.

[0010] Preferably, in step S2, the specific method is as follows: S2.1 Based on the digital elevation model data and orbital geometric parameters obtained in step S1, calculate the terrain phase component and the flatland phase component according to the principle of radar interferometry, estimate the ionospheric phase component using the split spectrum method, and subtract the terrain phase component, flatland phase component and ionospheric phase component from the phase value of the original interferogram one by one to obtain mixed residual phase data containing tropospheric delayed phase, orbital residual phase and noise phase. S2.2. A Gaussian low-pass filter with a window size of 32 pixels × 32 pixels is used to filter the mixed residual phase data and remove high-frequency components with spatial wavelengths less than 100 meters. The low-frequency trend term of the filtered data is fitted with a second-order polynomial as the orbital residual phase. The orbital residual phase is subtracted from the filtered data to obtain the tropospheric atmospheric delay phase data, in which the phase value of each pixel is proportional to the atmospheric path integral delay in the radar slant range direction at the corresponding position.

[0011] Preferably, in step S3, the specific method is as follows: S3.1 The target atmospheric space is divided into M layers along the height direction. The atmospheric refractive index in each layer is considered to be uniformly distributed. Based on the angular span formed by the long synthetic aperture of the beam-gathering mode in the azimuth direction, the slant range path length of the radar beam through each atmospheric layer is calculated at N different observation angles. The observation equation is established according to the linear relationship between phase delay and path integral refractive index, forming an N-row M-column tomographic observation matrix A. The matrix element Aij represents the propagation path length of the radar beam at the i-th observation angle in the j-th atmospheric layer. S3.2. Based on the observation matrix A and the covariance matrix Cφ of the observation phase noise established in step S3.1, derive the posterior covariance matrix Cn=(A^T·Cφ^(-1)·A)^(-1) of the atmospheric refractive index inversion result according to Bayesian estimation theory. Extract the square root of the diagonal elements of the covariance matrix Cn as the theoretical inversion standard deviation of the atmospheric refractive index of each layer. Calculate the ratio of the number of non-zero elements in the observation matrix A to the total number of elements as the tomographic grid coverage index.

[0012] Preferably, in step S4, the specific method is as follows: S4.1. The tomographic grid coverage calculated in step S3.2 is denoted as f1, and the maximum value of the theoretical inversion standard deviation of the atmospheric refractive index of each layer is denoted as f2. A bi-objective optimization function F(f1,f2) is established with the goal of maximizing f1 and minimizing f2. The value range of the number of sub-apertures K is set to 3 to 10. The value range of the oblique angle θk of the center of each sub-aperture is set to the starting angle to the ending angle of the total synthetic aperture angle span. The theoretical inversion standard deviation threshold is set to 0.5ppm as the optimization constraint. S4.2. The non-dominated sorting genetic algorithm is used to iteratively solve the optimization function F(f1,f2). The population size is set to 100, the crossover probability is 0.9, the mutation probability is 0.1, and the number of iterations is 200. In each iteration, the objective function value corresponding to different sub-aperture numbers K and central oblique angle combinations {θ1,θ2,...,θK} is calculated. The solution set that satisfies the theoretical inversion standard deviation is less than 0.5ppm is selected. The scheme with the largest tomographic grid coverage is selected from the solution set as the optimal sub-aperture configuration parameters.

[0013] Preferably, in step S5, the specific method is as follows: S5.1. Based on the K sub-aperture center oblique angles {θ1, θ2, ..., θK} obtained in step S4.2, denote the starting oblique angle of the total composite aperture as θstart and the ending oblique angle as θend. Calculate the angular distance from the center of the kth sub-aperture to the starting boundary and the angular distance from the center of the kth sub-aperture to the ending boundary. Calculate the above two distance values ​​for all sub-apertures respectively, and select the minimum value from the 2K distance values ​​as Δθmin. The formula for calculating the angular distance from the center of the kth sub-aperture to the starting boundary is as follows: Δθk_start = θk - θstart; The formula for calculating the angular distance from the center of the kth sub-aperture to the termination boundary is as follows: Δθk_end = θend - θk; S5.2. Take twice the minimum angular distance Δθmin calculated in step S5.1 as the uniform sub-aperture angular length. With the oblique angle θk of the center of each sub-aperture as the center, extend the angular range of Δθmin to both sides. Determine the starting oblique angle of the kth sub-aperture as θk-Δθmin and the ending oblique angle as θk+Δθmin. According to the inverse relationship between azimuth resolution and synthetic aperture angular length, each sub-aperture with the same angular length Lθ corresponds to the same azimuth resolution. The formula for calculating the uniform sub-aperture angle length is as follows: Lθ=2Δθmin.

[0014] Preferably, in step S6, the specific method is as follows: S6.1. Based on the starting and ending angles of each sub-aperture determined in step S5.2, extract echo data of the corresponding angle range from the complex imaging data obtained in step S1.2. Perform range pulse compression and azimuth matched filtering on the echo data of the K sub-apertures respectively to obtain K complex images {I1,I2,...,IK}. The pixel values ​​of each complex image contain amplitude information and phase information. Select the complex image whose center angle is closest to the center angle of the total synthetic aperture as the main image Im. S6.2. Using the pixel coordinate system of the main image Im as the registration reference coordinate system, calculate the geometric offset relative to the main image for the remaining K-1 complex images. Use the frequency domain registration method to estimate the sub-pixel offset parameters of each complex image in the range and azimuth directions. Perform bilinear interpolation resampling on each complex image according to the offset parameters. Multiply each resampled complex image with the main image Im pixel by pixel conjugate to obtain K-1 interferograms {Φ1,Φ2,...,ΦK-1}, where the pixel value of each interferogram is the complex phase difference.

[0015] Preferably, in step S7, the specific method is as follows: S7.1 Divide the target atmospheric space into a grid in the horizontal direction according to the pixel intervals in the range and azimuth directions, and divide it into layers in the vertical direction according to the height interval of 500 meters to form a three-dimensional voxel grid. The total number of voxels is denoted as P. Based on the oblique angle of each sub-aperture center obtained in step S4.2 and the precise orbital ephemeris data loaded in step S1.1, calculate the spatial intersection relationship between K-1 observation lines and each system. Use the ray tracing algorithm to calculate the path length lij of the i-th observation line through the j-th voxel, and construct a path length matrix L with (K-1) rows and P columns. S7.2. Expand the phase values ​​of the K-1 interferograms obtained in step S6.2 into an observation vector Φobs with a dimension of (K-1)×1 according to the pixel position. Form the atmospheric refractive index parameters in P voxels into a vector N to be determined with a dimension of P×1. Based on the relationship that the phase delay is equal to the integral of the product of the path length and the refractive index along the path, establish the observation equation Φobs=L·N. This set of equations contains K-1 observation equations and P unknown parameters. The element lij of matrix L represents the sensitivity weight of the i-th observation angle to the refractive index of the j-th voxel.

[0016] Preferably, in step S8, the specific method is as follows: S8.1. To address the underdetermined problem caused by the number of observations K-1 being much smaller than the number of unknown parameters P in the observation equation Φobs=L·N established in step S7.2, a discrete cosine transform is performed on the refractive index vector N to be determined to obtain the transform domain coefficient vector α=DCT(N). Taking advantage of the physical property that the atmospheric refractive index field changes slowly in the vertical direction, most coefficients in the transform domain are close to zero. The objective function J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1 is constructed, where the first term is the data fitting term, the second term is the L1 norm sparse regularization term, and λ is the regularization parameter. The formula for constructing the objective function is as follows: J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1; S8.2. Set the initial value of the regularization parameter λ to 0.01, set the iteration termination threshold to 10^(-6), and use the iterative threshold shrinkage algorithm to minimize the objective function J(α). Calculate the gradient and update the coefficient vector in the t-th iteration. Terminate the iteration when ||αt+1-αt|| < 10^(-6). Perform the inverse discrete cosine transform N = DCT^(-1)(α) on the final coefficient vector to obtain the atmospheric refractive index vector. The formula for calculating the gradient is as follows: J(αt)=2L^T(L·DCT^(-1)(αt)-Φobs); The formula for updating the coefficient vector is as follows: αt+1=soft(αt-μ) J(αt),λμ); Where soft is the soft thresholding function and μ is the step size parameter.

[0017] Preferably, in step S9, the specific method is as follows: S9.1. The atmospheric refractive index vector N obtained in step S8.2 is assigned P refractive index parameter values ​​to the corresponding three-dimensional voxels in the order of the voxel grid divided in step S7.1 to form discrete three-dimensional refractive index field data. For voxel positions in the observation matrix L where all column elements are zero, the average refractive index of the adjacent non-zero voxels is used to fill the position, thus establishing a complete three-dimensional voxel refractive index dataset. S9.2. The three-dimensional voxel refractive index dataset is smoothed using a three-dimensional Gaussian filter. The filter window size is set to 3×3×3 voxels. The voxel mesh is refined to twice the original resolution using bilinear interpolation along the horizontal direction. The height interval is refined to 250 meters using linear interpolation along the vertical direction. The interpolated three-dimensional refractive index distribution data is then output.

[0018] Compared with the prior art, the beneficial effects of the present invention are: 1. This method addresses the resolution mismatch caused by unreasonable observation geometry configuration and viewing angle differences by constructing a multi-objective optimization mechanism and a resolution consistency constraint framework. Compared to existing technologies that rely on manual experience or simple uniform partitioning, this scheme utilizes a quantitative accuracy evaluation model to automate the configuration of sub-aperture parameters, ensuring uniform imaging quality across all viewing angles while expanding the effective coverage of the tomographic observation grid. Furthermore, the introduced sparse regularized inversion algorithm leverages the prior characteristics of atmospheric physical distribution to alleviate the instability of the tomographic equations, enhances the inversion process's ability to suppress observation noise, ensures the numerical stability and physical rationality of the refractive index inversion results, and avoids non-physical oscillations and data drift during the solution process. Through this systematic parameter optimization and algorithmic constraints, this invention improves the robustness of the inversion process, providing robust mathematical model support for obtaining high-quality atmospheric refractive index inversion results.

[0019] 2. This method represents a leap from traditional two-dimensional path integral observation to three-dimensional spatial refractive index field reconstruction. Traditional detection methods are limited by the observation dimension, making it difficult to distinguish the complex heterogeneity of the troposphere in the vertical direction. This method, however, uses multi-angle tomographic inversion and multi-scale interpolation processing techniques to acquire three-dimensional distribution data and vertical profile structure, thus improving the resolution capability of atmospheric spatial detection. This improvement enriches the data dimensions of atmospheric remote sensing, providing accurate data support for fields such as weather forecasting, atmospheric environment correction, and disaster monitoring. Building upon the advantages of all-weather, all-time large-area monitoring by spaceborne radar, this invention enhances the spatiotemporal precision of tropospheric atmospheric detection, providing a technical approach for realizing atmospheric environmental perception. This three-dimensional inversion architecture breaks through the limitations of single-path observation, enabling the realistic depiction of detailed changes in atmospheric vertical stratification, thereby enhancing the application value of satellite remote sensing technology in refined meteorological detection. Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating the steps of the method described in this application. Detailed Implementation

[0021] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of this application.

[0022] The terms "first," "second," and "third" in this application are for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first," "second," or "third" may explicitly or implicitly include at least one of that feature. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified. All directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of this application are only used to explain the relative positional relationships and movements between components in a specific orientation (as shown in the figures). If the specific orientation changes, the directional indications also change accordingly. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices.

[0023] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a mutually exclusive, independent, or alternative embodiment. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0024] Example 1: Please refer to Figure 1 A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry is characterized by the following specific steps: S1. Data Acquisition and Purification: Acquire raw echo data from the spaceborne focused synthetic aperture radar in the target observation area, and load digital elevation model and precise orbital ephemeris data; sequentially perform deskewing, coarse focusing and Doppler parameter estimation on the raw echo data to obtain basic imaging data with a unified spatiotemporal reference. S2. Phase purification and isolation: Calculate the topographic phase component, flat terrain phase component, and ionospheric error phase component according to the imaging geometry, and remove these phase components from the interferogram; remove noise through a filtering algorithm and extract the residual interferometric phase that characterizes the tropospheric atmospheric delay. S3. Accuracy Model Construction: Based on the geometric characteristics of the long synthetic aperture in the clustering mode, the observation matrix expression for the tropospheric atmospheric refractive index tomography inversion is established; the posterior covariance matrix and tomographic grid coverage function are derived to construct a quantitative evaluation model for the inversion accuracy. S4. Geometric structure optimization: Establish an optimization function with the dual objectives of maximizing the coverage of the tomographic grid and minimizing the theoretical inversion error; use a non-dominated sorting genetic algorithm to iteratively divide the long synthetic aperture, and obtain the optimal number of sub-apertures and the corresponding center oblique angle values ​​of each sub-aperture under the condition of meeting the preset accuracy threshold. S5. Consistency Constraint: Based on the optimal sub-aperture center position obtained in step S4, calculate the distance from the center of each sub-aperture to the two sides of the overall composite aperture, and take the minimum value among all distances as the uniform sub-aperture length so that the imaging results of each sub-aperture maintain the same resolution in the azimuth direction. S6. Parallel Imaging Interferometry: Imaging processing is performed on the data of each sub-aperture separately to obtain complex images corresponding to each sub-aperture; the complex image corresponding to the central viewpoint is selected as the registration reference, and sub-pixel level registration and resampling are performed on the complex images of the remaining sub-apertures to obtain multiple sets of interferograms; S7. Construction of tomographic equations: The target atmospheric space is divided into a three-dimensional voxel grid. The path length of the radar beam through each voxel is calculated based on the orbital parameters of each sub-aperture to construct the observation matrix. A three-dimensional tomographic linear equation set between the multi-angle observation phase and atmospheric refractive index parameters is established. S8. Sparse Regularized Inversion: Construct an objective function containing data fitting terms and sparse regularization terms, and use an iterative threshold shrinkage algorithm to solve the objective function to obtain a stable solution for the atmospheric refractive index parameter. S9. Three-dimensional field reconstruction: The atmospheric refractive index parameters obtained in step S8 are mapped to a three-dimensional spatial grid. After spatial interpolation and smoothing, the three-dimensional spatial distribution data and vertical profile data of the tropospheric atmospheric refractive index are output.

[0025] In this embodiment: Step S1 establishes a unified spatiotemporal reference framework through the collaborative loading and systematic preprocessing of multi-source data. Deskewing reduces the data sampling rate while alleviating storage burden, improving the computational efficiency of subsequent processing by approximately tenfold. Coarse focusing enhances signal energy concentration while preserving complete phase information through range and azimuth signal compression. Doppler parameter estimation ensures the accuracy of imaging geometry by accurately calculating the center frequency and modulation frequency. These preprocessing operations work together to provide the obtained complex imaging data with a unified spatiotemporal reference, thereby avoiding the accumulation of phase errors introduced by inconsistent references during subsequent sub-aperture segmentation, and laying a reliable data foundation for high-precision interferometry and tomographic inversion.

[0026] Step S2 achieves high-purity separation of the tropospheric atmospheric delay signal by systematically eliminating non-tropospheric phase components. This step first uses a digital elevation model and orbital geometry parameters to calculate and remove topographic and flat phases, eliminating the influence of surface undulations on the observation results. Then, the split spectrum method is used to estimate and remove the ionospheric error phase, avoiding interference from the upper atmosphere. Finally, Gaussian low-pass filtering and polynomial fitting are used to remove residual orbital errors and high-frequency noise. After this series of processing steps, the signal-to-noise ratio of the extracted tropospheric delay phase is improved by approximately 15 dB, and the phase accuracy reaches within 0.1 radians, providing high-quality observational input for subsequent tomographic inversion and effectively preventing the propagation and amplification of error phases during the inversion process.

[0027] Step S3 achieves a quantitative assessment of inversion accuracy by establishing the observation matrix expression and deriving the posterior covariance matrix. Based on the geometric characteristics of the long synthetic aperture in the clustering mode, this step establishes a linear mapping relationship between the observation phase and atmospheric refractive index, clarifying the sensitivity of different observation angles to each atmospheric layer. The derived posterior covariance matrix quantifies the uncertainty of the inversion results into the theoretical standard deviation of the refractive index of each atmospheric layer. Simultaneously, the introduced tomographic grid coverage function assesses the adequacy of spatial sampling of the observation system by statistically analyzing the proportion of non-zero elements in the observation matrix. These quantitative indicators transform the abstract inversion quality into a computable mathematical expression, providing a clear objective function and criterion for subsequent sub-aperture optimization, and avoiding the uncontrollable accuracy problem caused by the reliance on empirical parameters in traditional methods.

[0028] Step S4 achieves automated optimization of sub-aperture configuration parameters through multi-objective optimization using a non-dominated sorting genetic algorithm. This step treats maximizing tomographic grid coverage and minimizing theoretical inversion error as a dual-objective optimization problem. Utilizing the global search capability of the genetic algorithm, it searches for a Pareto optimal solution set in the multi-dimensional parameter space of the number of sub-apertures and the central oblique angle. By setting a constraint that the theoretical inversion standard deviation is less than 0.5 ppm, it filters configuration schemes that meet accuracy requirements from the optimal solution set, ultimately selecting the scheme with the highest coverage as the optimal configuration. This automated optimization method increases the tomographic grid coverage to over 85%, improving inversion accuracy by approximately 30% compared to traditional uniform partitioning schemes, overcoming the limitations of empirical configuration methods that struggle to balance coverage and accuracy.

[0029] Step S5 achieves strict consistency in the azimuth resolution of the imaging results of each sub-aperture by calculating the minimum boundary distance and unifying the sub-aperture length. This step calculates the angular distance from the center of each optimized sub-aperture to the start and end boundaries of the overall synthetic aperture, and selects the minimum value from all distances as the unified sub-aperture half-length. This setting method ensures that all sub-apertures do not exceed the effective range of the overall synthetic aperture, and also ensures that each sub-aperture has the same angular span. Since azimuth resolution is inversely proportional to the angular length of the synthetic aperture, the unified angular length ensures that the azimuth resolution deviation of each sub-aperture image is less than 0.1 meters, thus solving the interference decoherence problem caused by resolution mismatch in traditional methods and maintaining the interference coherence coefficient above 0.8.

[0030] Step S6 achieves efficient generation of multiple sets of highly coherent interferograms through parallel imaging processing and sub-pixel-level precise registration. This step first applies range pulse compression and azimuth matched filtering to each sub-aperture data, improving imaging efficiency using a parallel computing architecture. Then, the central viewpoint image is selected as the registration reference, and the geometric offset of the remaining sub-aperture images is estimated using a frequency domain registration method, improving the registration accuracy to within 0.1 pixels. Finally, bilinear interpolation resampling eliminates geometric distortion, and the interferogram is generated by conjugate multiplication with the main image. The generated interferograms exhibit a mean coherence coefficient exceeding 0.85 and a phase standard deviation less than 0.15 radians. Compared to traditional single-view interferometry methods, this adds multiple independent observation dimensions, providing sufficient angular diversity information for three-dimensional tomographic inversion.

[0031] Step S7, through three-dimensional voxel mesh generation and ray tracing algorithms, achieves a dimensional expansion from two-dimensional interferometry to three-dimensional atmospheric structure inversion. This step divides the target atmospheric space into regular voxel meshes at pixel intervals horizontally and at 500-meter height intervals vertically, establishing a discretized representation of the atmospheric space. Subsequently, based on the orbital parameters and observation geometry of each sub-aperture, the ray tracing algorithm calculates the path length of the radar beam through each voxel, constructing an observation matrix. Based on the physical relationship that phase delay equals the integral of the product of path length and refractive index along the path, a three-dimensional tomographic linear equation system connecting the phase of multi-angle observations and the voxel refractive index is established. This three-dimensional modeling method transforms the complex inversion problem into a standard linear algebraic problem, achieving an atmospheric vertical resolution of 500 meters, overcoming the fundamental limitation of traditional integral methods in being unable to resolve vertical structures.

[0032] Step S8 achieves a stable solution to the ill-posed tomographic equations by introducing sparsity constraints in the discrete cosine transform domain and an iterative threshold shrinkage algorithm. This step addresses the underdetermined problem caused by the number of observations being much smaller than the number of unknown parameters. Utilizing the physical property that atmospheric refractive index changes slowly in the vertical direction, most coefficients of the refractive index vector approach zero after transformation to the discrete cosine domain. Based on this sparsity prior, an objective function is constructed, including a data fitting term and an L1 norm regularization term. The regularization term suppresses high-frequency oscillations in the solution, and the iterative threshold shrinkage algorithm is used to solve this objective function, converging to a stable solution within 200 iterations. The standard deviation of the obtained refractive index parameter solution is reduced by approximately 50% compared to the unconstrained least squares method, effectively avoiding solution drift and noise amplification in ill-posed inversion, and keeping the refractive index value within a physically reasonable range of 200 ppm to 400 ppm.

[0033] Step S9, through spatial mapping, interpolation refinement, and smoothing, transforms discrete inversion parameters into a continuous three-dimensional refractive index field product. This step first maps the solved refractive index parameters to three-dimensional space according to a voxel grid. For voxels in the observation blind zone, the average value of adjacent non-zero voxels is used to fill in the spatial data, completing the spatial data. Then, a three-dimensional Gaussian filter with a window size of 3×3×3 voxels is used for smoothing to suppress local noise. Finally, bilinear interpolation is used in the horizontal direction and linear interpolation in the vertical direction, increasing the horizontal resolution to twice the original pixel interval and refining the vertical resolution to 250 meters. The output three-dimensional refractive index distribution data and vertical profile data can be directly used for meteorological analysis and atmospheric correction applications. Compared with traditional two-dimensional products, it adds complete vertical structure information, providing high-quality data products for refined research of the troposphere.

[0034] This method addresses the resolution mismatch caused by unreasonable observation geometry configuration and viewing angle differences by constructing a multi-objective optimization mechanism and a resolution consistency constraint framework. Compared to existing technologies that rely on manual experience or simple uniform partitioning, this scheme utilizes a quantitative accuracy evaluation model to automate the configuration of sub-aperture parameters, ensuring uniform imaging quality across all viewing angles while expanding the effective coverage of the tomographic observation grid. Furthermore, the introduced sparse regularized inversion algorithm leverages prior characteristics of atmospheric physical distribution to mitigate the instability of the tomographic equations, enhances the inversion process's ability to suppress observation noise, ensures the numerical stability and physical rationality of the refractive index inversion results, and avoids non-physical oscillations and data drift during the solution process. Through this systematic parameter optimization and algorithmic constraints, this invention improves the robustness of the inversion process, providing robust mathematical model support for obtaining high-quality atmospheric refractive index inversion results.

[0035] This method represents a leap from traditional two-dimensional path integral observation to three-dimensional spatial refractive index field reconstruction. Traditional detection methods are limited by the observation dimension, making it difficult to distinguish the complex heterogeneity of the troposphere in the vertical direction. This method, however, uses multi-angle tomographic inversion and multi-scale interpolation techniques to acquire three-dimensional distribution data and vertical profile structure, thus improving the resolution capability of atmospheric spatial detection. This improvement enriches the data dimensions of atmospheric remote sensing, providing accurate data support for fields such as weather forecasting, atmospheric environment correction, and disaster monitoring. Building upon the advantages of all-weather, all-time large-area monitoring by spaceborne radar, this invention enhances the spatiotemporal precision of tropospheric atmospheric detection, providing a technical approach for atmospheric environmental perception. This three-dimensional inversion architecture overcomes the limitations of single-path observation, realistically depicting the detailed changes in atmospheric vertical stratification, thereby enhancing the application value of satellite remote sensing technology in refined meteorological detection.

[0036] Example 2: Please refer to Figure 1 In step S1, the specific method is as follows: S1.1 Acquire the raw echo data of the satellite-borne focused synthetic aperture radar in the target observation area, simultaneously load the digital elevation model data and satellite precise orbit ephemeris data covering the observation area, unify the spatial coordinate system of the digital elevation model with the radar imaging coordinate system, establish the spatial correspondence of multi-source data, and provide terrain reference and orbital geometric parameters for subsequent phase separation. S1.2. The original echo data is sequentially deskewing to reduce the data sampling rate, coarse focusing to compress the range and azimuth signals, and the Doppler center frequency and Doppler modulation frequency parameters are estimated. The original echo data is then converted into complex imaging data with a unified spatiotemporal reference. This complex imaging data retains complete phase information and provides a data basis for sub-aperture segmentation and interferometric processing.

[0037] In this embodiment: Step S1.1 establishes a spatial mapping relationship between radar echoes, digital elevation models, and precise orbital data through synchronous loading of multi-source data and coordinate system transformation. This step completes the geometric alignment of data from different observation dimensions, eliminates spatial matching errors caused by reference differences, and provides terrain references and orbital geometric parameters for subsequent stripping of terrain and flat phases. Through the unified transformation of the spatial coordinate system, the geometric trajectory calculation of the radar beam's path through the atmosphere has a high-precision positioning reference, solving the problem of insufficient correlation of multi-source data in the tomographic modeling process.

[0038] Step S1.2 transforms the original echo signal into complex imaging data with spatiotemporal consistency by performing deskewing, coarse focusing, and Doppler parameter estimation. Deskewing reduces the data sampling rate and computational load of subsequent processing, while coarse focusing achieves preliminary signal energy compression while preserving complete phase information. The complex imaging data produced in this step provides a physically consistent data carrier for subsequent sub-aperture segmentation, ensuring phase continuity of sub-apertures at different observation angles during interferometric processing and reducing information loss during signal conversion.

[0039] Compared to traditional techniques that directly process raw echoes, this method introduces coordinate system transformation and refined Doppler parameter estimation at the data processing front end, enhancing the spatiotemporal correlation between radar payload, terrain benchmark, and satellite orbit. This improvement enhances the geometric fidelity of imaging data, provides standardized data input for constructing the observation matrix in subsequent tomographic inversion, thereby reducing systematic errors from signal acquisition to phase extraction and improving the initial data quality for atmospheric refractive index acquisition.

[0040] Example 3: Please refer to Figure 1 In step S2, the specific method is as follows: S2.1 Based on the digital elevation model data and orbital geometric parameters obtained in step S1, calculate the terrain phase component and the flatland phase component according to the principle of radar interferometry, estimate the ionospheric phase component using the split spectrum method, and subtract the terrain phase component, flatland phase component and ionospheric phase component from the phase value of the original interferogram one by one to obtain mixed residual phase data containing tropospheric delayed phase, orbital residual phase and noise phase. S2.2. A Gaussian low-pass filter with a window size of 32 pixels × 32 pixels is used to filter the mixed residual phase data and remove high-frequency components with spatial wavelengths less than 100 meters. The low-frequency trend term of the filtered data is fitted with a second-order polynomial as the orbital residual phase. The orbital residual phase is subtracted from the filtered data to obtain the tropospheric atmospheric delay phase data, in which the phase value of each pixel is proportional to the atmospheric path integral delay in the radar slant range direction at the corresponding position.

[0041] In this embodiment: Step S2.1, through the comprehensive application of digital elevation model and orbital geometry parameters, accurately separates the topographic and flat-ground phase components, and uses split spectrum technology to suppress ionospheric phase errors. This step subtracts multiple deterministic physical components from the original interferometric phase, reducing the complexity of the signal components. By systematically subtracting various interference terms, hybrid residual phase data focusing on tropospheric delay characteristics is obtained, providing a good signal background for subsequent phase purification and reducing the impact of surface undulations and upper-level atmospheric interference on lower-level inversion.

[0042] Step S2.2 achieves precise separation of high-frequency noise and orbital residual errors in the mixed phase data through Gaussian low-pass filtering within a specific window and second-order polynomial fitting. This step filters out random spatial disturbances and removes low-frequency orbital trend terms through modeling, obtaining clean tropospheric atmospheric delay phase data. The processed phase data establishes a correspondence with the atmospheric path integral delay, enhancing the physical interpretability of the observations, ensuring the accuracy of the input vector in 3D tomographic modeling, and providing data support for improving the robustness of atmospheric refractive index inversion.

[0043] Step S2 improves upon the limitations of traditional atmospheric phase extraction techniques, such as incomplete error separation and unsatisfactory signal-to-noise ratio, by constructing a hierarchical isolation and collaborative purification process. Compared to existing methods relying on a single correction model or empirical smoothing, this approach employs a strategy combining item-by-item elimination and multi-scale processing to enhance the resolution of tropospheric atmospheric delayed phases. This improvement strengthens the analytical efficiency of the processing flow for multi-source coupling errors under complex environments, enabling the extracted phase values ​​to reflect the delayed characteristics of atmospheric propagation paths. Through this series of purification operations, the initial model bias in the tomographic equation construction stage is reduced, providing a physically logical observational data foundation for the precise reconstruction of the atmospheric three-dimensional structure.

[0044] Example 4: Please refer to Figure 1 In step S3, the specific method is as follows: S3.1 The target atmospheric space is divided into M layers along the height direction. The atmospheric refractive index in each layer is considered to be uniformly distributed. Based on the angular span formed by the long synthetic aperture of the beam-gathering mode in the azimuth direction, the slant range path length of the radar beam through each atmospheric layer is calculated at N different observation angles. The observation equation is established according to the linear relationship between phase delay and path integral refractive index, forming an N-row M-column tomographic observation matrix A. The matrix element Aij represents the propagation path length of the radar beam at the i-th observation angle in the j-th atmospheric layer. S3.2. Based on the observation matrix A and the covariance matrix Cφ of the observation phase noise established in step S3.1, derive the posterior covariance matrix of the atmospheric refractive index inversion result according to Bayesian estimation theory, extract the square root of the diagonal elements of the covariance matrix Cn as the theoretical inversion standard deviation of the atmospheric refractive index of each layer, and calculate the ratio of the number of non-zero elements in the observation matrix A to the total number of elements as the tomographic grid coverage index. The calculation formula is as follows; Cn=(A^T·Cφ^(-1)·A)^(-1).

[0045] In this embodiment: Step S3.1 establishes a linear mapping relationship between phase observations and atmospheric refractive index parameters by performing layered discretization of the target atmospheric space and geometric calculation of the radar slant range path. This step completes the transformation from physical observation space to mathematical matrix space, clarifies the propagation characteristics of the radar beam at different altitude layers, and forms a standardized tomographic observation matrix. This modeling method provides a rigorous mathematical basis for inversion calculations, establishes the detection sensitivity of atmospheric structure from different observation perspectives, and solves the problem of inaccurate spatial mapping description in the atmospheric parameter inversion process.

[0046] Step S3.2 derives the posterior covariance matrix using Bayesian estimation theory, enabling quantitative prediction of inversion errors and scientific evaluation of tomographic grid coverage effectiveness. This step completes the quantitative modeling of inversion accuracy, assessing the sufficiency of information acquisition under the current observation geometry by calculating the theoretical standard deviation of refractive index and coverage index. This evaluation method provides clear target feedback for subsequent optimization of sub-aperture configurations, enabling pre-judgment of the reliability of inversion results and improving the controllability and expected accuracy of atmospheric refractive index detection.

[0047] Step S3 improves upon the limitations of traditional atmospheric parameter inversion methods by constructing a tomographic observation matrix and an accuracy evaluation model, addressing the shortcomings in model representation and the lack of objective evaluation criteria. Compared to conventional methods using simplified path integrals or empirical interpolation, this approach enhances the analytical accuracy of the inversion algorithm for atmospheric vertical structure through discretized modeling of atmospheric stratification characteristics and the derivation of error propagation mechanisms. This improvement enhances the effectiveness of observational resource utilization, provides quantitative criteria for the automated optimization of subsequent sub-aperture parameters, reduces computational biases caused by improper geometric configuration, ensures the mathematical rigor of the tomographic inversion process, and provides reliable model support for achieving high-precision three-dimensional refractive index reconstruction.

[0048] Example 5: Please refer to Figure 1 In step S4, the specific method is as follows: S4.1. The tomographic grid coverage calculated in step S3.2 is denoted as f1, and the maximum value of the theoretical inversion standard deviation of the atmospheric refractive index of each layer is denoted as f2. A bi-objective optimization function F(f1,f2) is established with the goal of maximizing f1 and minimizing f2. The value range of the number of sub-apertures K is set to 3 to 10. The value range of the oblique angle θk of the center of each sub-aperture is set to the starting angle to the ending angle of the total synthetic aperture angle span. The theoretical inversion standard deviation threshold is set to 0.5ppm as the optimization constraint. S4.2. The non-dominated sorting genetic algorithm is used to iteratively solve the optimization function F(f1,f2). The population size is set to 100, the crossover probability is 0.9, the mutation probability is 0.1, and the number of iterations is 200. In each iteration, the objective function value corresponding to different sub-aperture numbers K and central oblique angle combinations {θ1,θ2,...,θK} is calculated. The solution set that satisfies the theoretical inversion standard deviation is less than 0.5ppm is selected. The scheme with the largest tomographic grid coverage is selected from the solution set as the optimal sub-aperture configuration parameters.

[0049] In this embodiment: Step S4.1 constructs a mathematical model for observation geometry optimization by integrating grid coverage and theoretical inversion error indices. This step establishes a selection mechanism guided by maximizing coverage and minimizing error, and sets the search range for the number of sub-apertures and the oblique viewing angle. By introducing an accuracy threshold as a constraint, a clear evaluation boundary is provided for subsequent parameter searches, ensuring that the optimization process can improve the spatial sampling efficiency of the observation system while meeting data reliability requirements.

[0050] Step S4.2 utilizes a non-dominated sorting genetic algorithm to iteratively search the optimization function, completing the screening and determination of sub-aperture configuration schemes. This step, through population evolution and parameter recombination, identifies the optimal solution set that meets accuracy requirements in a complex multi-dimensional space. Further screening of the solution set determines configuration parameters that balance inversion error and grid coverage, eliminating uncertainties introduced by human experience and ensuring that the final determined observation geometry provides a better data acquisition perspective for subsequent atmospheric three-dimensional reconstruction.

[0051] Step S4 improves upon the traditional sub-aperture configuration method, which relies on empirical settings or simple uniform partitioning, by establishing a dual-objective optimization function and an iterative solution process. Compared to fixed-parameter partitioning, this scheme enhances the utilization efficiency of observation angle resources in tomographic modeling by jointly optimizing inversion error and spatial coverage effectiveness. This improvement strengthens the matching degree between observation geometry and atmospheric layering structure, reduces the generation of observation blind spots, and provides higher-quality geometric support for the subsequent construction of tomographic equations, thereby improving the resolution accuracy of atmospheric refractive index in the vertical direction and the stability of inversion results.

[0052] Example 6: Please refer to Figure 1 In step S5, the specific method is as follows: S5.1. Based on the K sub-aperture center oblique angles {θ1, θ2, ..., θK} obtained in step S4.2, denote the starting oblique angle of the total composite aperture as θstart and the ending oblique angle as θend. Calculate the angular distance from the center of the kth sub-aperture to the starting boundary and the angular distance from the center of the kth sub-aperture to the ending boundary. Calculate the above two distance values ​​for all sub-apertures respectively, and select the minimum value from the 2K distance values ​​as Δθmin. The formula for calculating the angular distance from the center of the kth sub-aperture to the starting boundary is as follows: Δθk_start = θk - θstart; The formula for calculating the angular distance from the center of the kth sub-aperture to the termination boundary is as follows: Δθk_end = θend - θk; S5.2. Take twice the minimum angular distance Δθmin calculated in step S5.1 as the uniform sub-aperture angular length. With the oblique angle θk of the center of each sub-aperture as the center, extend the angular range of Δθmin to both sides. Determine the starting oblique angle of the kth sub-aperture as θk-Δθmin and the ending oblique angle as θk+Δθmin. According to the inverse relationship between azimuth resolution and synthetic aperture angular length, each sub-aperture with the same angular length Lθ corresponds to the same azimuth resolution. The formula for calculating the uniform sub-aperture angle length is as follows: Lθ=2Δθmin.

[0053] In this embodiment: Step S5.1 completes the boundary constraint analysis of available angle resources by traversing and calculating the angular distance from the center of each sub-aperture to the physical boundary of the composite aperture, and determines the maximum safe radius that can be selected without exceeding the total aperture range. This step achieves a unified measurement of discretely distributed sub-apertures in geometric space, provides a numerical basis for the subsequent construction of imaging sequences with consistent resolution, ensures that all data segments participating in interferometric processing have complete original echo support, and avoids image quality degradation caused by missing edge data.

[0054] Step S5.2, by uniformly setting the angular span of the sub-apertures, completes the mapping transformation from the geometric center coordinates to the physical aperture range, establishing a sub-aperture segmentation criterion with consistent azimuth resolution. Utilizing the physical inverse relationship between angular length and azimuth resolution, this step ensures the matching of multiple sets of complex images in the spatial sampling structure, eliminating the problem of decreased interferometric coherence caused by aperture length differences. This provides a physical guarantee for the generation of high-quality multi-baseline interferograms and enhances the robustness of subsequent tomographic inversion observation data.

[0055] The aforementioned steps, through the introduction of length unification processing under boundary constraints, improve the interferometric decoherence problem caused by azimuth resolution mismatch during multi-view imaging. Compared to the limitations of existing technologies that struggle to balance edge sub-aperture integrity with global resolution consistency, this scheme enhances the coherence characteristics between sub-aperture images at different oblique angles by performing boundary calculations and setting equivalent lengths for the synthetic aperture space. This improvement enhances the signal-to-noise ratio of the interferometric phase and reduces random phase noise introduced by resolution mismatch, thereby providing a highly consistent observation vector for the accurate construction of the three-dimensional tomographic linear equations and ensuring the accuracy of tropospheric atmospheric refractive index inversion in the vertical dimension.

[0056] Example 7: Please refer to Figure 1 In step S6, the specific method is as follows: S6.1. Based on the starting and ending angles of each sub-aperture determined in step S5.2, extract echo data of the corresponding angle range from the complex imaging data obtained in step S1.2. Perform range pulse compression and azimuth matched filtering on the echo data of the K sub-apertures respectively to obtain K complex images {I1,I2,...,IK}. The pixel values ​​of each complex image contain amplitude information and phase information. Select the complex image whose center angle is closest to the center angle of the total synthetic aperture as the main image Im. S6.2. Using the pixel coordinate system of the main image Im as the registration reference coordinate system, calculate the geometric offset relative to the main image for the remaining K-1 complex images. Use the frequency domain registration method to estimate the sub-pixel offset parameters of each complex image in the range and azimuth directions. Perform bilinear interpolation resampling on each complex image according to the offset parameters. Multiply each resampled complex image with the main image Im pixel by pixel conjugate to obtain K-1 interferograms {Φ1,Φ2,...,ΦK-1}, where the pixel value of each interferogram is the complex phase difference.

[0057] In this embodiment: Step S6.1 achieves the systematic generation of multi-view complex images through precise extraction and imaging processing of echo data from each sub-aperture. It completes pulse compression of the echo signal in the range and azimuth directions, preserving the original signal amplitude and phase information while establishing a stable reference benchmark for interferometric processing by selecting the central oblique view image as the master image. This processing method provides a highly consistent data source for subsequent multi-baseline interferometry, ensuring that phase information from different observation angles can be effectively compared within a unified imaging geometry framework, thus improving the phase fidelity during the conversion of the original echo to a complex image.

[0058] Step S6.2 utilizes frequency domain registration technology and sub-pixel-level resampling to achieve high-precision spatial alignment among multiple complex images. By calculating geometric offsets and implementing bilinear interpolation, geometric distortions and coordinate misalignments introduced by differences in observation angles are eliminated, ensuring point-to-point correspondence between pixels in the main and auxiliary images. Through conjugate multiplication, interferometric patterns reflecting spatial phase difference characteristics are successfully extracted, providing multi-dimensional observational inputs for subsequent atmospheric tomography inversion. This high-precision registration process reduces the impact of spatial decoherence noise, improves the quality of the interferometric phase, and enhances the robustness of the inversion modeling.

[0059] The aforementioned steps, through the construction of a sub-aperture parallel imaging and sub-pixel-level interferometric processing workflow, improve upon the technical limitations of traditional interferometry, such as high phase noise and poor coherence caused by excessively large viewing angle spans and insufficient registration accuracy. Compared to existing methods using conventional image registration or single-view interferometry, this scheme enhances the coupling consistency of multi-view observation sequences in spatial coordinates and phase references by performing refined frequency domain registration and resampling on multi-aperture complex images. This improvement not only expands the observation angle information required for atmospheric tomography inversion but also improves the extraction accuracy of interferometric phases through high-precision coordinate alignment, thereby providing more accurate and reliable observation vectors for the subsequent construction of three-dimensional linear equations and laying a high-quality data foundation for the precise reconstruction of the tropospheric atmospheric refractive index.

[0060] Example 8: In step S7, the specific method is as follows: S7.1 Divide the target atmospheric space into a grid in the horizontal direction according to the pixel intervals in the range and azimuth directions, and divide it into layers in the vertical direction according to the height interval of 500 meters to form a three-dimensional voxel grid. The total number of voxels is denoted as P. Based on the oblique angle of each sub-aperture center obtained in step S4.2 and the precise orbital ephemeris data loaded in step S1.1, calculate the spatial intersection relationship between K-1 observation lines and each system. Use the ray tracing algorithm to calculate the path length lij of the i-th observation line through the j-th voxel, and construct a path length matrix L with (K-1) rows and P columns. S7.2. Expand the phase values ​​of the K-1 interferograms obtained in step S6.2 into an observation vector Φobs with a dimension of (K-1)×1 according to the pixel position. Form the atmospheric refractive index parameters in P voxels into a vector N to be determined with a dimension of P×1. Based on the relationship that the phase delay is equal to the integral of the product of the path length and the refractive index along the path, establish the observation equation Φobs=L·N. This set of equations contains K-1 observation equations and P unknown parameters. The element lij of matrix L represents the sensitivity weight of the i-th observation angle to the refractive index of the j-th voxel.

[0061] In this embodiment: Step S7.1 involves dividing the target atmospheric space into three-dimensional voxel grids and performing ray tracing calculations using precise orbital ephemeris and observation oblique angles to complete the geometric path modeling of the radar beam within the discretized spatial units. This step establishes the spatial intersection relationship between the observation line of sight and the voxel units, accurately quantifying the specific path length contribution of each atmospheric layer to radar signal propagation. By constructing a path length matrix, a precise geometric mapping foundation is provided for three-dimensional tomographic inversion, solving the problem of quantitatively describing continuous atmospheric space in inversion modeling.

[0062] Step S7.2 establishes a three-dimensional tomographic linear equation system based on the path integral physics principle by vectorizing the multi-angle interferometric phase and voxel refractive index parameters. This step completes the mathematical mapping from radar interferometric phase observations to the atmospheric physical parameters to be determined, clarifying the physical meaning of the observation matrix elements as weights for refractive index sensitivity. By constructing this observation equation, the complex atmospheric refractive index inversion problem is transformed into a standard linear algebra solution problem, providing a complete mathematical framework for subsequently introducing sparse constraints and obtaining highly stable refractive index spatial solutions.

[0063] Step S7 improves upon the limitations of traditional atmospheric remote sensing by establishing a three-dimensional voxel grid and a tomographic linear equation system, overcoming the technical bottleneck of only being able to acquire total path delay and failing to distinguish vertical stratification. Compared to existing methods that use simplified path models or statistical empirical interpolation, this scheme enhances the analytical capability for atmospheric vertical structural heterogeneity through the synergistic application of ray tracing calculations and physical path integral modeling. This improvement enhances the physical fidelity of the observation equations, enabling the inversion system to capture subtle changes in the troposphere at different altitudes. Through this series of modeling operations, core algorithmic support is provided for achieving high spatial resolution three-dimensional precise reconstruction of atmospheric refractive index, thereby improving the precision of satellite remote sensing technology's perception of the atmospheric environment.

[0064] Example 9: Please refer to Figure 1 In step S8, the specific method is as follows: S8.1. To address the underdetermined problem caused by the number of observations K-1 being much smaller than the number of unknown parameters P in the observation equation Φobs=L·N established in step S7.2, a discrete cosine transform is performed on the refractive index vector N to be determined to obtain the transform domain coefficient vector α=DCT(N). Taking advantage of the physical property that the atmospheric refractive index field changes slowly in the vertical direction, most coefficients in the transform domain are close to zero. The objective function J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1 is constructed, where the first term is the data fitting term, the second term is the L1 norm sparse regularization term, and λ is the regularization parameter. The formula for constructing the objective function is as follows: J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1; S8.2. Set the initial value of the regularization parameter λ to 0.01, set the iteration termination threshold to 10^(-6), and use the iterative threshold shrinkage algorithm to minimize the objective function J(α). Calculate the gradient and update the coefficient vector in the t-th iteration. Terminate the iteration when ||αt+1-αt|| < 10^(-6). Perform the inverse discrete cosine transform N = DCT^(-1)(α) on the final coefficient vector to obtain the atmospheric refractive index vector. The formula for calculating the gradient is as follows: J(αt)=2L^T(L·DCT^(-1)(αt)-Φobs); The formula for updating the coefficient vector is as follows: αt+1=soft(αt-μ) J(αt),λμ); Where soft is the soft thresholding function and μ is the step size parameter.

[0065] In this embodiment: Step S8.1 addresses the underdetermined solution problem caused by the limited number of observations by utilizing the physical property of atmospheric refractive index variation in the vertical direction and mapping the refractive index parameter to the sparse domain through spatial transformation. This step completes the construction of the objective function, which includes data fitting terms and norm regularization terms, and introduces physical prior constraints. This modeling approach provides a mathematical means to solve the instability of the tomographic equations, establishes optimization criteria for extracting key physical features and suppressing invalid interferences in the transform domain, and solves the problem of instability in the inversion model caused by the lack of observation perspectives.

[0066] Step S8.2 optimizes the objective function using an iterative threshold shrinkage algorithm, completing gradient calculation, soft threshold update, and inverse transformation of transform domain coefficients to spatial domain refractive index. This step achieves stable extraction of atmospheric refractive index parameters, and the iterative convergence mechanism ensures both accuracy and computational efficiency. This solution method avoids numerical oscillations common in traditional linear solutions, ensuring that the output atmospheric refractive index conforms to physical measurements, providing reliable parameter input for the final three-dimensional field reconstruction, and improving the inversion success rate in complex environments.

[0067] Step S8 improves upon the technical shortcomings of traditional tomographic inversion, which is prone to solution drift or non-physical oscillations when the observation angle is limited, by introducing transform domain sparsity constraints and an iterative shrinkage solution mechanism. Compared to existing methods that rely on least squares or empirical smoothing, this scheme utilizes the physical prior of atmospheric vertical variation to transform the solution of the instability equation into an optimization problem with sparsity protection characteristics, enhancing the inversion algorithm's resistance to system noise and model errors. This improvement enhances the inversion stability of the spatial distribution of atmospheric refractive index, ensuring that physically reasonable three-dimensional parameter solutions can still be obtained under specific conditions. Through this systematic constraint solution process, the dependence of the inversion solution on initial values ​​is reduced, ensuring the mathematical rigor of the tomographic inversion process and providing core algorithmic support for achieving high-precision reconstruction of the three-dimensional structure of the troposphere.

[0068] Example 10: Please refer to Figure 1 In step S9, the specific method is as follows: S9.1. The atmospheric refractive index vector N obtained in step S8.2 is assigned P refractive index parameter values ​​to the corresponding three-dimensional voxels in the order of the voxel grid divided in step S7.1 to form discrete three-dimensional refractive index field data. For voxel positions in the observation matrix L where all column elements are zero, the average refractive index of the adjacent non-zero voxels is used to fill the position, thus establishing a complete three-dimensional voxel refractive index dataset. S9.2. The three-dimensional voxel refractive index dataset is smoothed using a three-dimensional Gaussian filter with a filter window size of 3×3×3 voxels. The voxel mesh is refined to twice the original resolution using bilinear interpolation along the horizontal direction, and the height interval is refined to 250 meters using linear interpolation along the vertical direction. The interpolated three-dimensional refractive index distribution data is output, the average refractive index of each height layer in the vertical direction is extracted, and the vertical profile data of refractive index variation with height is output.

[0069] In this embodiment: Step S9.1 completes the numerical backfilling from one-dimensional mathematical vectors to three-dimensional physical space by spatially mapping the inverted refractive index vectors according to a preset grid order. For the blind voxels in the observation matrix caused by the radar line of sight not being passed through, this step uses a neighborhood mean filling strategy to fill in the missing data, establishing a spatially continuous and complete three-dimensional voxel dataset. This processing method eliminates data gaps caused by discrete inversion, ensures the geometric integrity of the atmospheric refractive index field, and provides fully covered initial field data for subsequent fine interpolation and product output.

[0070] Step S9.2 utilizes three-dimensional Gaussian filtering and multi-scale interpolation techniques to denoise and smooth the refractive index field data and enhance its spatial resolution. By implementing encrypted interpolation in both the horizontal and vertical directions, this step transforms discrete system center values ​​into high-precision continuous spatial distribution data and further extracts profile information reflecting the vertical structure of the atmosphere. This processing method improves the refinement of the final output product, enabling the inversion results to directly serve downstream applications such as meteorological analysis, and establishing an atmospheric environmental perception data product with high horizontal and vertical resolution.

[0071] Step S9 improves upon the limitations of traditional atmospheric sounding methods by constructing a spatial mapping backfilling and multi-scale enhancement processing flow, addressing the issues of limited data product format, poor spatial continuity, and insufficient resolution. Compared to existing methods that only output gridded data or coarse vertical profiles, this approach enhances the smoothness and precision of the atmospheric refractive index field's spatial distribution through intelligent filling of voxel blind zones and three-dimensional filtering interpolation. This improvement not only enhances the characterization and physical consistency of the data products but also transforms atmospheric parameters from discrete solutions to high-quality continuous field products. This provides more detailed and accurate three-dimensional information support for the refined structural study of the troposphere and quantitative meteorological correction, thereby increasing the application value of spaceborne radar in atmospheric remote sensing monitoring.

[0072] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated units described above can be implemented in hardware or as software functional units. The above are merely embodiments of this application and do not limit the patent scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

[0073] The specific embodiments of the invention have been described in detail above, but they are only examples, and this application is not limited to the specific embodiments described above. For those skilled in the art, any equivalent modifications or substitutions to the invention are also within the scope of this application. Therefore, all equivalent changes, modifications, and improvements made without departing from the spirit and principles of this application should be covered within the scope of this application.

Claims

1. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry, characterized in that: The specific method is as follows: S1. Data Acquisition and Purification: Acquire raw echo data from the spaceborne focused synthetic aperture radar in the target observation area, and load digital elevation model and precise orbital ephemeris data; sequentially perform deskewing, coarse focusing and Doppler parameter estimation on the raw echo data to obtain basic imaging data with a unified spatiotemporal reference. S2. Phase purification and isolation: Calculate the topographic phase component, flat terrain phase component, and ionospheric error phase component according to the imaging geometry, and remove these phase components from the interferogram; remove noise through a filtering algorithm and extract the residual interferometric phase that characterizes the tropospheric atmospheric delay. S3. Accuracy Model Construction: Based on the geometric characteristics of the long synthetic aperture in the clustering mode, the observation matrix expression for the tropospheric atmospheric refractive index tomography inversion is established; the posterior covariance matrix and tomographic grid coverage function are derived to construct a quantitative evaluation model for the inversion accuracy. S4. Geometric structure optimization: Establish an optimization function with the dual objectives of maximizing the coverage of the tomographic grid and minimizing the theoretical inversion error; use a non-dominated sorting genetic algorithm to iteratively divide the long synthetic aperture, and obtain the optimal number of sub-apertures and the corresponding center oblique angle values ​​of each sub-aperture under the condition of meeting the preset accuracy threshold. S5. Consistency Constraint: Based on the optimal sub-aperture center position obtained in step S4, calculate the distance from the center of each sub-aperture to the two sides of the overall composite aperture, and take the minimum value among all distances as the uniform sub-aperture length so that the imaging results of each sub-aperture maintain the same resolution in the azimuth direction. S6. Parallel Imaging Interferometry: Imaging processing is performed on the data of each sub-aperture separately to obtain complex images corresponding to each sub-aperture; the complex image corresponding to the central viewpoint is selected as the registration reference, and sub-pixel level registration and resampling are performed on the complex images of the remaining sub-apertures to obtain multiple sets of interferograms; S7. Construction of tomographic equations: The target atmospheric space is divided into a three-dimensional voxel grid. The path length of the radar beam through each voxel is calculated based on the orbital parameters of each sub-aperture to construct the observation matrix. A three-dimensional tomographic linear equation set between the multi-angle observation phase and atmospheric refractive index parameters is established. S8. Sparse Regularized Inversion: Construct an objective function containing data fitting terms and sparse regularization terms, and use an iterative threshold shrinkage algorithm to solve the objective function to obtain a stable solution for the atmospheric refractive index parameter. S9. Three-dimensional field reconstruction: The atmospheric refractive index parameters obtained in step S8 are mapped to a three-dimensional spatial grid. After spatial interpolation and smoothing, the three-dimensional spatial distribution data and vertical profile data of the tropospheric atmospheric refractive index are output.

2. The high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 1, characterized in that, In step S1, the specific method is as follows: S1.1 Acquire the raw echo data of the satellite-borne focused synthetic aperture radar in the target observation area, simultaneously load the digital elevation model data and satellite precise orbit ephemeris data covering the observation area, unify the spatial coordinate system of the digital elevation model with the radar imaging coordinate system, establish the spatial correspondence of multi-source data, and provide terrain reference and orbital geometric parameters for subsequent phase separation. S1.

2. The original echo data is sequentially deskewing to reduce the data sampling rate, coarse focusing to compress the range and azimuth signals, and the Doppler center frequency and Doppler modulation frequency parameters are estimated. The original echo data is then converted into complex imaging data with a unified spatiotemporal reference. This complex imaging data retains complete phase information and provides a data basis for sub-aperture segmentation and interferometric processing.

3. The high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 2, characterized in that, In step S2, the specific method is as follows: S2.1 Based on the digital elevation model data and orbital geometric parameters obtained in step S1, calculate the terrain phase component and the flatland phase component according to the principle of radar interferometry, estimate the ionospheric phase component using the split spectrum method, and subtract the terrain phase component, flatland phase component and ionospheric phase component from the phase value of the original interferogram one by one to obtain mixed residual phase data containing tropospheric delayed phase, orbital residual phase and noise phase. S2.

2. A Gaussian low-pass filter with a window size of 32 pixels × 32 pixels is used to filter the mixed residual phase data and remove high-frequency components with spatial wavelengths less than 100 meters. The low-frequency trend term of the filtered data is fitted with a second-order polynomial as the orbital residual phase. The orbital residual phase is subtracted from the filtered data to obtain the tropospheric atmospheric delay phase data, in which the phase value of each pixel is proportional to the atmospheric path integral delay in the radar slant range direction at the corresponding position.

4. The high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 3, characterized in that, In step S3, the specific method is as follows: S3.1 The target atmospheric space is divided into M layers along the height direction. The atmospheric refractive index in each layer is considered to be uniformly distributed. Based on the angular span formed by the long synthetic aperture of the beam-gathering mode in the azimuth direction, the slant range path length of the radar beam through each atmospheric layer is calculated at N different observation angles. The observation equation is established according to the linear relationship between phase delay and path integral refractive index, forming an N-row M-column tomographic observation matrix A. The matrix element Aij represents the propagation path length of the radar beam at the i-th observation angle in the j-th atmospheric layer. S3.

2. Based on the observation matrix A and the covariance matrix Cφ of the observation phase noise established in step S3.1, derive the posterior covariance matrix of the atmospheric refractive index inversion result according to Bayesian estimation theory, extract the square root of the diagonal elements of the covariance matrix Cn as the theoretical inversion standard deviation of the atmospheric refractive index of each layer, and calculate the ratio of the number of non-zero elements in the observation matrix A to the total number of elements as the tomographic grid coverage index. The calculation formula is as follows; Cn=(A^T·Cφ^(-1)·A)^(-1).

5. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 4, characterized in that, In step S4, the specific method is as follows: S4.

1. The tomographic grid coverage calculated in step S3.2 is denoted as f1, and the maximum value of the theoretical inversion standard deviation of the atmospheric refractive index of each layer is denoted as f2. A bi-objective optimization function F(f1,f2) is established with the goal of maximizing f1 and minimizing f2. The value range of the number of sub-apertures K is set to 3 to 10. The value range of the oblique angle θk of the center of each sub-aperture is set to the starting angle to the ending angle of the total synthetic aperture angle span. The theoretical inversion standard deviation threshold is set to 0.5ppm as the optimization constraint. S4.

2. The non-dominated sorting genetic algorithm is used to iteratively solve the optimization function F(f1,f2). The population size is set to 100, the crossover probability is 0.9, the mutation probability is 0.1, and the number of iterations is 200. In each iteration, the objective function value corresponding to different sub-aperture numbers K and central oblique angle combinations {θ1,θ2,...,θK} is calculated. The solution set that satisfies the theoretical inversion standard deviation is less than 0.5ppm is selected. The scheme with the largest tomographic grid coverage is selected from the solution set as the optimal sub-aperture configuration parameters.

6. The high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 5, characterized in that, In step S5, the specific method is as follows: S5.

1. Based on the K sub-aperture center oblique angles {θ1, θ2, ..., θK} obtained in step S4.2, denote the starting oblique angle of the total composite aperture as θstart and the ending oblique angle as θend. Calculate the angular distance from the center of the kth sub-aperture to the starting boundary and the angular distance from the center of the kth sub-aperture to the ending boundary. Calculate the above two distance values ​​for all sub-apertures respectively, and select the minimum value from the 2K distance values ​​as Δθmin. The formula for calculating the angular distance from the center of the kth sub-aperture to the starting boundary is as follows: Δθk_start = θk - θstart; The formula for calculating the angular distance from the center of the kth sub-aperture to the termination boundary is as follows: Δθk_end = θend - θk; S5.

2. Take twice the minimum angular distance Δθmin calculated in step S5.1 as the uniform sub-aperture angular length. With the oblique angle θk of the center of each sub-aperture as the center, extend the angular range of Δθmin to both sides. Determine the starting oblique angle of the kth sub-aperture as θk-Δθmin and the ending oblique angle as θk+Δθmin. According to the inverse relationship between azimuth resolution and synthetic aperture angular length, each sub-aperture with the same angular length Lθ corresponds to the same azimuth resolution. The formula for calculating the uniform sub-aperture angle length is as follows: Lθ=2Δθmin.

7. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 6, characterized in that, In step S6, the specific method is as follows: S6.

1. Based on the starting and ending angles of each sub-aperture determined in step S5.2, extract echo data of the corresponding angle range from the complex imaging data obtained in step S1.

2. Perform range pulse compression and azimuth matched filtering on the echo data of the K sub-apertures respectively to obtain K complex images {I1,I2,...,IK}. The pixel values ​​of each complex image contain amplitude information and phase information. Select the complex image whose center angle is closest to the center angle of the total synthetic aperture as the main image Im. S6.

2. Using the pixel coordinate system of the main image Im as the registration reference coordinate system, calculate the geometric offset relative to the main image for the remaining K-1 complex images. Use the frequency domain registration method to estimate the sub-pixel offset parameters of each complex image in the range and azimuth directions. Perform bilinear interpolation resampling on each complex image according to the offset parameters. Multiply each resampled complex image with the main image Im pixel by pixel conjugate to obtain K-1 interferograms {Φ1,Φ2,...,ΦK-1}, where the pixel value of each interferogram is the complex phase difference.

8. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 7, characterized in that, In step S7, the specific method is as follows: S7.1 Divide the target atmospheric space into a grid in the horizontal direction according to the pixel intervals in the range and azimuth directions, and divide it into layers in the vertical direction according to the height interval of 500 meters to form a three-dimensional voxel grid. The total number of voxels is denoted as P. Based on the oblique angle of each sub-aperture center obtained in step S4.2 and the precise orbital ephemeris data loaded in step S1.1, calculate the spatial intersection relationship between K-1 observation lines and each system. Use the ray tracing algorithm to calculate the path length lij of the i-th observation line through the j-th voxel, and construct a path length matrix L with (K-1) rows and P columns. S7.

2. Expand the phase values ​​of the K-1 interferograms obtained in step S6.2 into an observation vector Φobs with a dimension of (K-1)×1 according to the pixel position. Form the atmospheric refractive index parameters in P voxels into a vector N to be determined with a dimension of P×1. Based on the relationship that the phase delay is equal to the integral of the product of the path length and the refractive index along the path, establish the observation equation Φobs=L·N. This set of equations contains K-1 observation equations and P unknown parameters. The element lij of matrix L represents the sensitivity weight of the i-th observation angle to the refractive index of the j-th voxel.

9. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 7, characterized in that, In step S8, the specific method is as follows: S8.

1. To address the underdetermined problem caused by the number of observations K-1 being much smaller than the number of unknown parameters P in the observation equation Φobs=L·N established in step S7.2, a discrete cosine transform is performed on the refractive index vector N to be determined to obtain the transform domain coefficient vector α=DCT(N). Taking advantage of the physical property that the atmospheric refractive index field changes slowly in the vertical direction, most coefficients in the transform domain are close to zero. The objective function J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1 is constructed, where the first term is the data fitting term, the second term is the L1 norm sparse regularization term, and λ is the regularization parameter. The formula for constructing the objective function is as follows: J(α)=||Φobs-L·DCT^(-1)(α)||2^2+λ||α||1; S8.

2. Set the initial value of the regularization parameter λ to 0.01, set the iteration termination threshold to 10^(-6), and use the iterative threshold shrinkage algorithm to minimize the objective function J(α). Calculate the gradient and update the coefficient vector in the t-th iteration. Terminate the iteration when ||αt+1-αt|| < 10^(-6). Perform the inverse discrete cosine transform N = DCT^(-1)(α) on the final coefficient vector to obtain the atmospheric refractive index vector. The formula for calculating the gradient is as follows: J(αt)=2L^T(L·DCT^(-1)(αt)-Φobs); The formula for updating the coefficient vector is as follows: αt+1=soft(αt-μ) J(αt),λμ); Where soft is the soft thresholding function and μ is the step size parameter.

10. A high-precision inversion method for tropospheric atmospheric refractive index based on spaceborne focused SAR interferometry according to claim 7, characterized in that, In step S9, the specific method is as follows: S9.

1. The atmospheric refractive index vector N obtained in step S8.2 is assigned P refractive index parameter values ​​to the corresponding three-dimensional voxels in the order of the voxel grid divided in step S7.1 to form discrete three-dimensional refractive index field data. For voxel positions in the observation matrix L where all column elements are zero, the average refractive index of the adjacent non-zero voxels is used to fill the position, thus establishing a complete three-dimensional voxel refractive index dataset. S9.

2. The three-dimensional voxel refractive index dataset is smoothed using a three-dimensional Gaussian filter with a filter window size of 3×3×3 voxels. The voxel mesh is refined to twice the original resolution using bilinear interpolation along the horizontal direction, and the height interval is refined to 250 meters using linear interpolation along the vertical direction. The interpolated three-dimensional refractive index distribution data is output, the average refractive index of each height layer in the vertical direction is extracted, and the vertical profile data of refractive index variation with height is output.

Images