A stereo SAR geometric calibration method based on alternating direction method of multipliers

Abstract

The application discloses a kind of stereoscopic SAR geometric calibration methods based on alternating direction multiplier method, comprising: obtaining multiple SAR images covering the same area, extracting satellite ephemeris parameters, observation slant range and Doppler center frequency, and obtaining elevation model data as prior elevation information;Based on the range-Doppler model, the range equation and Doppler equation containing the estimated parameter vector are established;Joint adjustment global optimization objective function with L2 regularization constraint is constructed;Auxiliary variables are introduced to convert the objective function into an equality constrained optimization problem, and an augmented Lagrangian function is constructed;Iterative solution is carried out using alternating direction multiplier method, and the estimated parameter vector, auxiliary variable and dual variable are updated alternately;After meeting the convergence condition, the optimal estimate of ground point three-dimensional coordinates and slant range system error is output.The application can realize high-precision and high-stability three-dimensional geometric calibration without ground control points.

Description

Technical Field

[0001] This invention belongs to the fields of radar remote sensing technology and image processing technology, specifically a stereo SAR geometric calibration method based on the alternating direction multiplier method. Background Technology

[0002] Synthetic Aperture Radar (SAR), with its all-weather, all-day high-resolution imaging capabilities, has become an important tool for Earth observation and remote sensing positioning. With the in-depth advancement of my country's major project on high-resolution Earth observation systems, the successive launches of the Remote Sensing series and the Gaofen series of SAR satellites have placed higher demands on the geometric quality and positioning accuracy of SAR images. Geometric calibration of SAR images is a prerequisite for achieving high-precision positioning, and its core lies in the accurate solution of key parameters in the Range-Doppler (RD) model.

[0003] Currently, geometric calibration techniques for SAR images are mainly classified into the following categories:

[0004] (1) Ground control point-based field calibration method: As a traditional high-precision calibration method, this method obtains the coordinates of ground control points (GCPs) by setting up corner reflectors or active calibrators in the calibration field and combining them with high-precision measurements. Although the calibration accuracy is high, it requires a lot of manpower and resources for field measurements, and is limited by the construction and maintenance of the calibration field, resulting in high costs and long cycles, making it difficult to meet the needs of high-frequency and large-scale calibration in satellite operational use.

[0005] (2) Cross-calibration method: This method uses a reference image or reference data with known geometric accuracy as the control source, extracts corresponding points through image matching, and then calculates the systematic error of the image to be calibrated. Although it avoids field measurement, its accuracy depends heavily on the accuracy and timeliness of the reference image, and requires the image to be calibrated to have a high degree of overlap and consistency with the reference image. The constraints are relatively strict, and its applicability is limited in areas with drastic terrain undulations or frequent changes in ground features.

[0006] (3) Autonomous calibration method based on multi-view images: To address the shortcomings of the above methods, researchers began to explore the ability of satellites to image the same area multiple times, i.e., autonomous calibration. The basic principle is to use multiple overlapping SAR images to establish an RD positioning model for each image, and to perform joint adjustment and solution using the slant range system error and the ground target point (GTP) coordinates as unknowns, thereby achieving geometric calibration without the need for GCP. Therefore, the autonomous calibration method has become a new direction in the development of geometric calibration technology, showing broad application prospects.

[0007] The autonomous calibration system uses multiple SAR images from different perspectives to solve for the three-dimensional coordinates of ground points by combining the distance equation and the Doppler equation. and slant range systematic error Current research on autonomous calibration methods mainly focuses on optimizing model solutions.

[0008] Traditional Least Squares Method (LSM): Early autonomous calibration methods mostly used the overall least squares method to solve the problem, incorporating unknowns such as satellite orbital state vector, system error parameters, and a large number of ground point coordinates into the same normal equation for unified solution.

[0009] Existing research either utilizes symmetric geometric constraints to achieve geometric self-calibration of SAR images, or performs geometric self-calibration based on the consistency constraints of conjugate image point localization, in order to achieve accurate estimation of geometric parameters. These methods typically require a sufficient number of overlapping images or strict geometric constraints to ensure calibration quality. When the number of overlapping images is limited, the stability of the calibration results is poor. Furthermore, by normalizing the RD model and introducing a weighting strategy, model errors caused by different scales can be effectively suppressed, further improving positioning accuracy.

[0010] In summary, geometric calibration technology is evolving from the traditional model that relies on field control points to autonomous calibration based on multi-view images. Current research focuses on improving the accuracy and stability of geometric calibration under uncontrolled conditions through model solution optimization.

[0011] However, existing autonomous calibration methods still have the following technical problems:

[0012] (1) High sensitivity to initial values: Traditional local optimization methods (such as least squares method) heavily rely on the selection of initial values. When the initial value deviation is large or the observed geometry is poor, it is easy to get trapped in local minima, resulting in calibration failure or distorted results.

[0013] (2) The solution of the ill-conditioned equation is unstable: When the number of images is small (e.g., only 3 images) or the geometric configuration is weak, the observation equation is ill-conditioned, the normal equation matrix is ​​close to singular, resulting in large fluctuations in the solution results and poor robustness.

[0014] (3) Coupling of regularization and data fitting: In order to improve ill-conditioned problems, regularization terms are often introduced. However, directly merging regularization terms with residual terms for optimization increases the nonlinearity of the problem and makes it difficult to flexibly handle different types of regularization constraints. Summary of the Invention

[0015] To address the problems existing in the prior art, the present invention aims to provide a stereo SAR geometric calibration method based on the alternating direction multiplier method. This method can achieve high-precision and high-stability three-dimensional geometric calibration using multi-view SAR images without ground control points. Here, "stereo" refers to the processing method of jointly calculating the three-dimensional spatial position of ground targets based on at least two SAR images from different perspectives.

[0016] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0017] A stereo SAR geometric calibration method based on the alternating direction multiplier method includes the following steps: Step S1: Acquire multiple SAR images covering the same area, extract satellite ephemeris parameters, observation slant range, and Doppler center frequency for each image, and acquire elevation model data of the area as prior elevation information; Step S2: Based on the range-Doppler model, establish a range equation and a Doppler equation for each SAR image containing a vector of parameters to be estimated, wherein the vector of parameters to be estimated includes the three-dimensional coordinates of ground target points and the slant range systematic error; Step S3: Construct a joint adjustment global optimization objective function with L2 regularization constraints, wherein the objective function consists of a nonlinear data fitting term and a linear regularization term; Step S4: Introduce auxiliary variables with the same dimension as the vector of parameters to be estimated, and... The global optimization objective function is transformed into an optimization problem with equality constraints, and an augmented Lagrangian function is constructed. Step S5: The augmented Lagrangian function is iteratively solved using the alternating direction multiplier method. This iterative solution includes: fixing the auxiliary variable and the dual variable, minimizing the parameter vector to be estimated to obtain an updated parameter vector; fixing the updated parameter vector to be estimated and the dual variable, analytically solving the auxiliary variable to obtain an updated auxiliary variable; updating the dual variable based on the updated parameter vector to be estimated and the auxiliary variable. Step S6: When the original residual and the dual residual satisfy the preset convergence condition, the iteration stops, and the current parameter vector to be estimated is used as the optimal estimate of the ground point coordinates and slant range systematic error to complete the geometric calibration.

[0018] As a further improvement to the above technical solution:

[0019] Step S1 further includes: calculating the atmospheric delay compensation value using external meteorological data or atmospheric data, and performing atmospheric delay error compensation processing on the observation slant distance.

[0020] The distance equation and Doppler equation established in step S2 are as follows:

[0021]

[0022] in, It is the slant distance from the satellite to the ground target point. It is the satellite's position vector. These are the coordinates of an unknown ground target point. It is the slant range systematic error. It is atmospheric delay error. It is the Doppler center frequency of the ground target point. It is the radar wavelength. It is the satellite's velocity vector.

[0023] The joint adjustment global optimization objective function with L2 regularization constraints constructed in step S3 is: :

[0024]

[0025]

[0026]

[0027] in, For the first SAR image ( The distance residuals of (=1,2,...,n). For the first Doppler residuals of SAR images For regularization parameters, Let be the vector of parameters to be estimated. , Elevation values ​​obtained from the DEM. It is the first The satellite position vector corresponding to each SAR image. It is the first The observation slant range corresponding to each SAR image. It is the first The Doppler center frequency corresponding to each SAR image. It is the first Satellite velocity vectors corresponding to SAR images.

[0028] The augmented Lagrangian function constructed in step S4 for:

[0029]

[0030] in, For scaled dual variables; To introduce the vector of parameters to be estimated Auxiliary variables with consistent dimensions; For penalty parameters; For data fitting terms, ; For regularization terms, .

[0031] The p-subproblem in step S5 for minimizing the parameter vector to be estimated is expressed as:

[0032]

[0033] This subproblem is solved using the Levenberg-Marquardt algorithm for inner-layer iterations, where... For the number of iterations, For the first The parameter vector to be estimated after the next iteration The value, For the first The value of the auxiliary variable z after the nth iteration. For the first The dual variable after the second iteration The value of .

[0034] The z-subproblem in step S5, which involves analytically solving the auxiliary variable, is represented as follows:

[0035]

[0036] Its closed-form solution is:

[0037]

[0038] in, For the first The value of the auxiliary variable z after the next iteration.

[0039] The formula for updating the dual variable in step S5 is as follows:

[0040]

[0041] in, For the first The dual variable after the second iteration The value of .

[0042] The original residual in step S6 With dual residual They are defined as follows:

[0043]

[0044]

[0045] The preset convergence condition is: and ,in This is the preset minimum tolerance threshold.

[0046] The number of SAR images shall be no less than two.

[0047] The beneficial effects of this invention are:

[0048] (1) Superior decoupling framework and high convergence stability: This invention is the first to apply the Alternating Direction Multiplier Method (ADMM) to stereo SAR geometric calibration. By splitting variables and constructing augmented Lagrangian functions, the complex nonlinear data fitting problem is decoupled from the simple linear regularization problem. This approach makes each sub-problem easier to solve, resulting in fast overall convergence and high stability. At the same time, the augmented Lagrangian form of the ADMM framework effectively improves the ill-conditioned nature of the normal equations, and stable solutions can be obtained even when the number of images is small (e.g., only 3) and the geometric structure is weak, significantly expanding the applicability of stereo SAR autonomous calibration.

[0049] (2) High solution efficiency: For z-subproblems containing L2 regularization terms, this invention utilizes the characteristics of convex quadratic programming to derive and employ an explicit analytical solution (closed solution) for updating. Compared to traditional regularization solution methods that require inner-layer iterations, this step does not require any additional iteration process, significantly reducing computational complexity and improving overall scaling efficiency.

[0050] (3) Strong robustness and avoidance of local optima: Through the alternating iteration mechanism of ADMM, this invention coordinates the optimization of prior elevation constraints (regularization term) and observation data fitting (residual term) in separate subspaces. Compared with the traditional overall optimization method that directly couples the two, this invention reduces the nonlinearity of the objective function, effectively avoids the problems of high sensitivity to initial values ​​and iteration getting trapped in local minima, and improves the robustness and accuracy of parameter estimation. Attached Figure Description

[0051] Figure 1 This is a schematic diagram of the overall process of the present invention. Detailed Implementation

[0052] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0053] A stereo SAR geometric calibration method based on the alternating direction multiplier method is proposed. This method transforms the geometric calibration problem into a least squares problem with L2 regularization constraints. Then, the alternating direction multiplier method (ADMM) is used to decompose the original problem into two subproblems: one is a nonlinear data fitting subproblem (distance equation and Doppler equation), and the other is a linear regularization subproblem. By alternately solving these two subproblems and using dual variable updates, high-precision and high-stability parameter estimation is achieved, while effectively avoiding local minima and ill-conditioned problems.

[0054] Reference Figure 1 The specific implementation steps of the present invention are as follows:

[0055] Step S1: Acquire multi-view SAR image data and perform preprocessing.

[0056] Step S11: Select multiple SAR images with a certain degree of overlap that were observed in the same area as the data source.

[0057] Step S12: Extract satellite ephemeris parameters for each image, including satellite position vectors. and velocity vector Observation slant range R and Doppler center frequency .

[0058] Step S13: Obtain the Digital Elevation Model (DEM) data covering the area to extract prior elevation information of ground targets. .

[0059] Step S14: Addressing atmospheric delay error Atmospheric delay compensation value is calculated using external meteorological or atmospheric data. In this embodiment, the atmospheric delay compensation value It is given by the following formula:

[0060]

[0061] in, It is the angle of incidence. It is zenith delay. It consists of two parts: ionospheric zenith delay and tropospheric zenith delay Zenith delay For tropospheric zenith delay With ionospheric zenith delay sum.

[0062] Temperature, air pressure, and humidity are fundamental parameters describing tropospheric characteristics, and these three parameters are the main factors affecting tropospheric zenith delay. Knowing altitude, air pressure, and water vapor content allows for successful simulation of tropospheric zenith delay caused by dry and wet air. It is given by the following formula:

[0063]

[0064] in, For the variable height, It is the refractive index along the zenith direction. It is the refractive index. The delay is caused by dry air composed of oxygen and nitrogen. It is composed of water vapor and The delay caused by the composition of moist air It is the pressure of dry air. It's temperature. It is the compressibility of dry air. It is the pressure of moist air. It is the compressibility of moist air. It is the radar wavelength. and It is a function related to the wavelength of the radar signal, proposed by Owens. and The empirical formula is as follows:

[0065]

[0066] During SAR imaging, atmospheric data from the National Centers for Environmental Prediction (NCEP) were used to acquire pressure and temperature data, thereby ensuring the accuracy of tropospheric group delay correction.

[0067] The ionospheric zenith delay must also be considered. The speed at which radar signals travel through the ionosphere is affected by dispersion. The delay (in meters) is:

[0068]

[0069] in, The carrier frequency is TEC, and the total vertical electron content is 100 volts per minute. (TECU). TEC data is provided daily by the Center for Orbit Determination in Europe (CODE). After preprocessing and compensation, it can be considered suitable for subsequent model solutions. .

[0070] Step S2: Establish the basic geometric calibration model of SAR imagery.

[0071] For any SAR image used in calibration, based on the range-Doppler (RD) model, the range equation and Doppler equation between the target point and the antenna phase center are established. The mathematical expression of the geometric calibration basis model is:

[0072]

[0073] in, It is the slant distance from the satellite to the ground target point (i.e., the observation slant distance). It is the Doppler center frequency of the ground target point. These are the coordinates of an unknown ground target point. It is the slant range systematic error.

[0074] Step S3: Construct the joint adjustment objective function with regularization constraints.

[0075] In traditional autonomous calibration based on multi-view SAR images, the least squares method is usually used to solve the nonlinear range equation and Doppler equation using Taylor expansion. However, when the number of SAR images involved in the intersection is small, or when the spatial geometry between the images is weak (such as a small spatial intersection angle), the solution equations are prone to ill-conditioned problems, making the positioning results extremely sensitive to small perturbations in the initial parameters, and even causing iterative divergence.

[0076] To overcome the above problems, this invention introduces an L2 regularization constraint based on the classical adjustment model. The three-dimensional coordinates of the ground target points to be determined and the slant range systematic errors are combined into a parameter vector to be estimated. Therefore, the following joint adjustment global optimization objective function with regularization term is constructed. :

[0077]

[0078]

[0079]

[0080] in, For the first SAR image ( The distance residuals of (=1,2,...,n). For the first Doppler residuals of SAR images For regularization parameters, Elevation values ​​obtained from the DEM. It is the first The satellite position vector corresponding to each SAR image. It is the first The observation slant range corresponding to each SAR image. It is the first The Doppler center frequency corresponding to each SAR image. It is the first Satellite velocity vectors corresponding to SAR images.

[0081] Step S4: Introduce an auxiliary variable z to split the variable and construct the augmented Lagrangian function.

[0082] The objective function constructed above includes both highly nonlinear data fitting terms and linear regularization penalty terms. Directly using traditional Newton's method or gradient descent to optimize the overall function can easily lead to local optima and incurs enormous computational costs. Therefore, this invention uses the ADMM framework for decoupling and simplification.

[0083] To achieve decoupling between fitting complex nonlinear data and solving with simple regularization, a parameter vector to be estimated is introduced. Auxiliary variables with consistent dimensions This transforms the original unconstrained global optimization problem into an optimization problem with equality constraints:

[0084]

[0085] in, For data fitting terms; This is a regularization term.

[0086] Constructing the augmented Lagrangian function :

[0087]

[0088] in, For the scaled dual variables (Lagrange multiplier vectors). This is the penalty parameter.

[0089] Step S5: Solve using ADMM iteratively.

[0090] Step S51: Initialize auxiliary variables Dual variables and penalty parameters The iteration count starts from 0, i.e., the 0th iteration, k=0,1,2,3,... Using the idea of ​​the ADMM algorithm, the complex augmented Lagrangian optimization problem is decomposed into alternating updates of two independent subproblems and an update step of the dual variable.

[0091] Step S52: - Sub-problem (data fitting). Fixing the current round. and Only for data Minimize the solution:

[0092]

[0093] in, For the first The parameter vector to be estimated after the next iteration The value, For the first The value of the auxiliary variable z after the nth iteration. For the first The dual variable after the second iteration The value of .

[0094] This subproblem is essentially transformed into a nonlinear least squares problem with a quadratic damping term. This invention employs the Levenberg-Marquardt (LM) algorithm for inner-layer iterative solution, based on current... Using the initial value as an example, an approximate solution is obtained through several iterations. .

[0095] Step S53: z-subproblem (regularization). Fixed update and Solving only for the auxiliary variable z:

[0096]

[0097] in, For the first The value of the auxiliary variable z after the next iteration.

[0098] Since both terms of this subproblem are quadratic convex functions of z, we can directly obtain its closed-form solution by taking its derivative and setting the first derivative to zero:

[0099]

[0100] Step S54: Update the dual variable. Update the Lagrange dual variable using the gradient ascent method. This makes the equality constraint approximate the condition:

[0101]

[0102] Step S55: Convergence Judgment. After each outer layer iteration, calculate the original residual. and dual residuals Determine if the infinite norm satisfies and ( Stop iteration when the set minimum tolerance threshold is reached; otherwise, let... Then return to step S52 to continue execution.

[0103] Step S6: Output the results.

[0104] The final vector of parameters to be estimated The value is used as the coordinate of the ground point. and systematic error The optimal estimate is obtained, and geometric calibration is completed.

[0105] Finally, it is necessary to state that the above embodiments are only used to further illustrate the technical solution of the present invention in detail, and should not be construed as limiting the scope of protection of the present invention. Any non-essential improvements and adjustments made by those skilled in the art based on the above content of the present invention shall fall within the scope of protection of the present invention.

Claims

1. A stereo SAR geometric calibration method based on the alternating direction multiplier method, characterized in that, Includes the following steps: Step S1: Acquire multiple SAR images covering the same area, extract the satellite ephemeris parameters, observation slant range and Doppler center frequency of each image, and obtain the elevation model data of the area as prior elevation information; Step S2: Based on the range-Doppler model, establish a range equation and a Doppler equation for each SAR image, which include the vector of parameters to be estimated. The vector of parameters to be estimated includes the three-dimensional coordinates of the ground target points and the slant range system error. Step S3: Construct a joint adjustment global optimization objective function with L2 regularization constraints, wherein the objective function consists of a nonlinear data fitting term and a linear regularization term; Step S4: Introduce auxiliary variables with the same dimension as the parameter vector to be estimated, transform the global optimization objective function into an optimization problem with equality constraints, and construct the augmented Lagrangian function; Step S5: The augmented Lagrangian function is iteratively solved using the alternating direction multiplier method. This iterative solution includes: By fixing the auxiliary and dual variables, the estimated parameter vector is minimized to obtain the updated estimated parameter vector. By fixing the updated parameter vector to be estimated and the dual variable, the auxiliary variable is analytically solved to obtain the updated auxiliary variable; The dual variable is updated based on the updated parameter vector to be estimated and the auxiliary variable; Step S6: When the original residual and the dual residual satisfy the preset convergence condition, stop the iteration and use the current parameter vector to be estimated as the optimal estimate of the ground point coordinates and slant distance systematic error to complete the geometric calibration.

2. The calibration method according to claim 1, characterized in that: Step S1 further includes: calculating the atmospheric delay compensation value using external meteorological data or atmospheric data, and performing atmospheric delay error compensation processing on the observation slant distance.

3. The calibration method according to claim 1, characterized in that: The distance equation and Doppler equation established in step S2 are as follows: ； in, These are the coordinates of an unknown ground target point. It is the slant distance from the satellite to the ground target point. It is the satellite's position vector. It is the slant range systematic error. It is atmospheric delay error. It is the Doppler center frequency of the ground target point. It is the radar wavelength. It is the satellite's velocity vector.

4. The calibration method according to claim 3, characterized in that: The joint adjustment global optimization objective function with L2 regularization constraints constructed in step S3 is: : ； ； ； in, For the first SAR image ( The distance residuals of (=1,2,...,n). For the first Doppler residuals of SAR images For regularization parameters, Let be the vector of parameters to be estimated. , Elevation values ​​obtained from the DEM. It is the first The satellite position vector corresponding to each SAR image. It is the first The observation slant range corresponding to each SAR image. It is the first The Doppler center frequency corresponding to each SAR image. It is the first Satellite velocity vectors corresponding to SAR images.

5. The calibration method according to claim 4, characterized in that: The augmented Lagrangian function constructed in step S4 for: ； in, For scaled dual variables; To introduce the vector of parameters to be estimated Auxiliary variables with consistent dimensions; For penalty parameters; For data fitting terms, ; For regularization terms, .

6. The calibration method according to claim 5, characterized in that: The p-subproblem in step S5 for minimizing the parameter vector to be estimated is expressed as: ； This subproblem is solved using the Levenberg-Marquardt algorithm for inner-layer iterations, where... For the number of iterations, For the first The parameter vector to be estimated after the next iteration The value, For the first The value of the auxiliary variable z after the nth iteration. For the first The dual variable after the second iteration The value of .

7. The calibration method according to claim 6, characterized in that: The z-subproblem in step S5, which involves analytically solving the auxiliary variable, is represented as follows: ； Its closed-form solution is: ； in, For the first The value of the auxiliary variable z after the next iteration.

8. The calibration method according to claim 7, characterized in that: The formula for updating the dual variable in step S5 is as follows: ； in, For the first The dual variable after the second iteration The value of .

9. The calibration method according to claim 7, characterized in that: The original residual in step S6 With dual residual They are defined as follows: ； ； The preset convergence condition is: and ,in This is the preset minimum tolerance threshold.

10. The calibration method according to claim 1, characterized in that: The number of SAR images shall be no less than two.

Images