A Robust and High-Precision Three-Step Positioning Method for an Airborne SAR Platform

: When airborne synthetic aperture radar (SAR) encounters long-time global navigation satellite system (GNSS) denial, the system cannot eliminate inertial navigation system (INS) accumulated drift. Platform positioning technology based on SAR image-matching is one of the important auxiliary navigation methods. This paper proposes a three-step positioning method for an airborne SAR platform, which can achieve the robust and high-precision estimation of platform position and velocity. Firstly, the motion model of the airborne SAR platform is established and a nonlinear overde-termined equation set of SAR Range-Doppler based on the ground-control points set obtained by SAR image-matching is constructed. Then, to overcome the ill-conditioned problem generated by the singular Jacobian matrix when solving the equations directly, a three-step robust and high-precision estimation of platform position and velocity is achieved through singular value decomposition and equation decoupling. Furthermore, the error transfer model of systematic and random platform positioning errors is derived. Finally, a set of semi-physical simulation experiments of airborne SAR is conducted to verify the effectiveness of the positioning method and the accuracy of the error model presented in this paper.


Introduction
The airborne synthetic aperture radar (SAR) provides an important source of information that can conduct detection and reconnaissance without limitations imposed by time, weather, and environmental factors [1].Typically, the INS/GNSS combination navigation system is used for SAR payload platforms; however, GNSS satellite signals are susceptible to interference in complex electromagnetic environments, making it difficult to obtain highprecision navigation information, and is even impossible to use in severe cases, referred to as "GNSS denial" [2].In such an environment, the flight platform can only rely on an inertial navigation system to obtain high-precision positioning, making it challenging to conduct detection tasks [2], and even resulting in a loss of communication.
The image-matching-assisted navigation system is one of the feasible navigation solutions under GNSS denial.Using the measured map provided by the onboard highresolution radar and comparing it with the digital map stored in the onboard computer, the geographical coordinates of several control points are obtained based on image-matching technology in the SAR image.The platform position and speed are computed according to SAR imaging geometry.This method was first applied to INS/SAR combination navigation, where positional and heading deviations are used as filter observations to correct platform position, and researchers have conducted extensive research in this area [2][3][4][5].However, when the INS measurement is not precise enough or the platform's position has already developed an irreducible drift, the SAR system can be used for positioning [6].
This method has been extensively researched in missile guidance.Using the distance-Doppler data of airborne SAR, the missile body's position is estimated; INS measures the missile's velocity and the radar altimeter measures the missile's height [7][8][9][10].The impact of Earth curvature and rotation on the deviation is subsequently considered in the ultra-long-range, and the method is further optimized based on the characteristics of airborne SAR.The point-based approach is widely applied in missile re-entry-stage guidance; however, it cannot exploit multi-points information to conduct adjustments.Under relatively mild flight attitudes, it is not the optimal choice due to its susceptibility to errors.
For the airborne SAR scenario, research scholars have proposed the location method based on the azimuth gate, using the precise slant range at the focused azimuth moment and the geographic coordinates of the control points to conduct positioning.Referring to the cosine theorem, scholars proposed to find three control points within a single azimuth gate, and the carrier's three-dimensional position can be determined accordingly [11].Scholars in [12][13][14] used the connecting line of the control points on the azimuth line to constrain the platform's position.The projection point of the platform on the horizontal plane was calculated using the slant range and altimeter to obtain the positioning of the platform at the intermediate imaging moment.However, constraining the platform's position with control points connecting lines can turn a matching error into an angular error and magnify it by distance.Scholars in [15] improved this approach by using the coarsely estimated location of each azimuth moment to accurately estimate the SAR platform's motion equation using a least-squares support vector regression.
All of the positioning methods mentioned above use the air pressure altimeter for positioning.However, the difference between the pressure altitude and the actual geographical altitude can result in significant positioning errors, and the azimuth gate positioning method inadequately utilizes Doppler information.Scholars in [16] proposed using 3D maps and flight parameters to construct simulated SAR images and compare them with real images to estimate the position by correcting flight parameters.They achieved good results.Based on the abovementioned methods, this paper proposes a robust and highly accurate three-step positioning method for an airborne SAR platform.
In Section 2, a motion model of the airborne SAR platform is established, and an equation is constructed based on the SAR scene-matching and the obtained set of ground feature points.This equation represents the nonlinear overdetermined relationship between SAR range and Doppler.During the process of solving the equations, a three-step approach is employed to achieve robust and high-precision estimations of the platform's position and velocity.This approach addresses the ill-conditioned problem caused by the singularity of the Jacobian matrix during direct solving, through techniques such as singular value decomposition and the decoupling of the equation system.Section 3 analyzes the systematic and random errors of the estimation method.In Section 4, semi-physical simulation experiments are conducted using actual SAR images and platform flight records.The influence of system parameters on estimation errors is also analyzed.Finally, a discussion and summary are provided to conclude the findings.

Platform Positioning Model
In motion compensation for airborne SAR, the flight trajectory is often sampled and fitted as a straight line for imaging processing.The inverted flight trajectory obtained from the imaging results should also be a straight line.Therefore, a first-order polynomial with constant coefficients is used to fit the aircraft trajectory.If the coefficients can be determined, the fitted position of the platform at any time can be determined according to Equation (1).A local image coordinate system is established with the image center as the origin, the azimuth direction as the Y-axis, the range direction as the X-axis, and the elevation direction as the Z-axis.As shown in Figure 1.
In this coordinate system, the equation of motion trajectory for the flight platform is established.This trajectory is the ideal trajectory during imaging, which is the trajectory after motion compensation.
The coordinates (x 0 , y 0 , z 0 ) represent the zero-time position in the three directions, while v x , v y , v z represent the average velocity in the three directions, which can be rewritten as a matrix. where In the equation, η i represents the zero-Doppler location of the target, i is focused after imaging, and X ηi represents the coordinates of the radar at the slow time η i .
Under the motion model of the platform, the position of each feature point in the SAR image, the radar position at the imaging time, the slant range, and the Doppler frequency all satisfy the following equation system: where x ηi , y ηi , z ηi represents the coordinates of the radar at the slow time η i and (x i , y i , z i ) represents the coordinates of the control point X i .R i represents the slant range at the azimuth time of the control point and f D represents the Doppler frequency.The coordinates of control points can be obtained through heterogeneous image-matching with optical images containing high-precision positioning information.The slant range of control points can be obtained from focused images, while the center frequency of the Doppler can typically be obtained from the imaging processor.
Based on Equation (4), considering the errors caused by inaccurate matching, imaging, and DEM data, for any feature point X i (x i , y i , z i ) in the image, we have: where R i is obtained from the focused image and the ideal value of ∆R is 0. In reality, it is composed of the deviation of the matching control-point position, matching error, and slant range error.Similarly, Equation ( 5) has: where ∆ f d is composed of the estimation error of the Doppler frequency and the position error of the control points.Two equations are established for each control point.If there are n control points, the coordinates and slant range and Doppler frequency of the control points can be used to obtain a nonlinear overdetermined equation system with the number of equations being 2n.The coefficient matrix A is estimated to minimize ∆R and ∆ f d , which can be regarded as the solution of A, and then the coordinates of the radar at each time during the imaging process can be obtained.During the estimation process, we assumed that the observed measurements, including control point positions, slant ranges, and Doppler frequencies, were accurate.We iterated based on these measurements to estimate the correction amount for linear trajectory errors.However, in practical applications, these data were not entirely precise, introducing errors to the estimation results.Therefore, in Section 3, we conducted a theoretical analysis of these errors, including systematic and random components.In Section 4, we plotted error curves to visualize their characteristics.

Three-Step Solution Method 2.2.1. Gauss-Newton Iteration
According to the parameter estimation theory, the Gauss-Newton iterative method can be used to solve nonlinear over-determined equations [17].The core idea is to locally linearize the nonlinear equations and iteratively minimize the residuals.Equation ( 6) can be rewritten as: where δx 0 , δv x , δy 0 , δx y , δz 0 , δv z is the error correction term and the unknown parameter to be estimated.The first-order expansion of Equation ( 8) is: where the coefficients are: Assuming n available feature points, these equations can be expressed in matrix form: The bias estimation is obtained from the least-squares principle: Utilize the estimated δA correction parameter matrix A to recalculate the Jacobian matrix and obtain a new least-squares solution.Continue to iterate until ∆R is below the threshold or the iteration limit is reached.
Similar to the slant range equation, the Doppler equation can also be rewritten and expanded as: where f D is the Doppler frequency obtained during the imaging process, and f d is the Doppler frequency calculated from the position and velocity parameters.Assuming n available feature points, the equation can be expressed in matrix form: The bias estimation is obtained from the least-squares principle: Utilize the estimated δA correction parameter matrix A to recalculate the Jacobian matrix and obtain a new least-squares solution.Continue to iterate until ∆ f d is below the threshold or the iteration limit is reached.

Singular Value Decomposition
In the Gauss-Newton iterative method, the matrix needs to be inverted (Equations ( 13) and ( 17)); however, matrices J T r J r and J T f J f cannot ensure full rank.When the platform is imaged in a positive side view, y ηi − y i in the equation tends to zero.As a result, J r3 , J r4 is approximated to be zero, indicating matrix singularity and an ill-posed equation, making it impossible to solve the equation.The singular value decomposition truncation method is used to remove the minimum singular values from the matrix [18].Let N = J T J, then N can be decomposed as: where U is the eigenvector matrix of NN T , V is the eigenvector matrix of NN T , and Σ is the diagonal matrix composed of NN T positive singular values.Singular values less than 10 −2 in Σ are truncated, and Equation ( 13) is modified to Equation ( 19):

Three-Step Robust Solution
After the singular value decomposition removes the minimum singular values, the equation loses its estimation function for the parameters corresponding to the small singular values.When imaged in a positive side view, the parameters y 0 and v y corresponding to J r3 and J r4 cannot be corrected and estimated.The estimation depends on other equations.In the Jacobian matrix of the Doppler equation: Under the condition of a positive side view, if the platform position error is not considered: At this point, the Jacobian matrix components J f 1 J f 3 J f 5 and J f 2 J f 4 J f 6 are linearly dependent, that is, the matrix J T f J f is singular.In practice, when the platform position error is considered as a small quantity, J T f J f is approximately singular.Although the singular value decomposition can solve the equation singularity problem, the effect of truncating small singular values results in some platform state parameters cannot be effectively estimated.
To avoid the linear dependence problem in the Jacobian matrix of the Doppler equation, this paper proposed a robust and high-precision three-step positioning method for an airborne SAR platform.Based on the solution of the slant range equation, the method divides the unknown parameter vector of the Doppler equation into two groups, namely, position and velocity vectors, and solves them step by step to solve the singularity problem in the Jacobian matrix of the Doppler equation, and finally achieve high-precision platform positioning.

1.
In the first step, the six parameters are roughly estimated by the slant range equation.
In the second step, the three x 0 , y 0 , z 0 parameters are accurately estimated by the Doppler equation, namely, the position estimation at zero time. 3.
In the third step, the three v x , v y , v z parameters are accurately estimated by the Doppler equation, namely, the average velocity estimation in three directions.
In all three steps, the Gauss-Newton iterative method based on singular value decomposition is used to effectively avoid the ill-posed problem of the equation.To compare the performance difference of the solution of the positioning equation before and after the decomposition of the Doppler equation, Monte-Carlo simulations were performed 1000 times using the simulation parameters in Section 4.1, and the solution errors of each step equation were statistically analyzed.The method that does not group the Doppler equation is called the two-step method.The convergence comparison of each parameter in the two-and three-step methods is shown in Figure 2, where s1, s2, and s3 represent the three-step equation-solving operations in one iteration, and s3 is ignored in the two-step method without updates.From the experiments, it can be observed from subfigure (b) that the azimuth position y 0 of the platform cannot be effectively estimated using the two-step method, and the error remains consistent with the initial value, while the three-step method can rapidly converge.From subfigure (f) the vertical velocity V z fluctuates in the two-step method, while it stably converges in the three-step method.And from subfigures (a), (c), (d) and (e), the error convergence values of the distance and vertical position x 0 , z 0 and the distance and azimuth velocity V x , V y obtained from the two methods are basically consistent.

Error Analysis
Due to the inaccuracy of the input parameters, there may be errors between the equation solving results and the actual position of the platform.These errors can be analyzed as systematic and random errors.Systematic errors are caused by the unified bias of the parameters, including the overall slant range error ∆R t , the Doppler parameter error ∆ f d , and the overall positioning offset ∆x t , ∆y t , ∆z t of the ground-control points in three directions, which mainly affect the mean deviation of the estimation results.Random errors are caused by the random errors of the parameters, such as the random errors in the coordinates of each control point in all directions and the random errors in the imaging slant range of each point, which mainly affect the statistical variance of the estimation results.

Systematic Slant Range Error
During SAR imaging, due to atmospheric refraction, receiver delay, and other reasons, there may be a fixed slant range bias ∆R in the imaging results.According to the relative position relationship between the platform and the control points in the slant range equation, the systematic slant range error produces a positional system bias in three directions according to the projection of the downward view angle and the slant view angle, it can be intuitively seen from Figure 3, and the mathematical expression is as shown in Equation (32): ∆x p = sin(θ s ) cos(γ)∆r t ∆y p = sin(γ)∆r t ∆z p = cos(θ s ) cos(γ)∆r t (32) sin cos sin cos cos

Systematic Control-Point Positioning Error
Due to errors in image-matching, optical map geodetic coordinate errors, and other reasons, there may be a systematic error in the geodetic positioning coordinates of the control points compared to the true coordinates.Assuming that there is a systematic error ∆x t , ∆y t , ∆z t in all control points in three directions, which is equivalent to translating the image coordinate system by (−∆x t , −∆y t , −∆z t ), the relative relationship between the platform positioning results and the control points remains unchanged.Therefore, a bias of (∆x t , ∆y t , ∆z t ) also occurs in the image coordinate system, as shown in Equation (33): ∆x p = ∆x t ∆y p = ∆y t ∆z p = ∆z t (33)

Doppler Shift
In typical scenarios, the Doppler center frequency of control points is estimated during the imaging process and used for image-focusing.Our method intended to obtain this parameter from the imaging processor; at this point, it should be highly accurate.However, in special cases, this parameter may not be available and needs to be estimated through alternative methods.In such cases, the target Doppler center frequency we used may deviate from the true Doppler frequency; it affects the position and velocity estimations in the two-step estimation of the Doppler equation.The second-step equation is used to estimate the position, which produces position estimation errors.Expanding the vector form of the Doppler equation: It can be observed that the Doppler error mainly affects y p during the position estimation.
The third-step equation was used to estimate the velocity, which produces velocity estimation errors.Assuming that ∆ → v p−t is the velocity error in the direction from the platform to the control points, the Doppler error can be expressed as follows: Decomposing ∆ → v p−t into three directions, we achieve:

Random Error Analysis
In addition to systematic biases, the positioning coordinates of control points and the slant range at each point of the image generally exhibit varying degrees of random errors.The following analysis was conducted from the perspective of the three-step equation.

Error Analysis of Slant Range Equation
In the first step of solving the slant range equation, the error source ∆R was derived from two parts: the slant range random error ∆R 1 during imaging, and the slant range error ∆R 2 indirectly caused by the random error of the control points location.The two are independent.
(1) Slant range random error ∆R 1 : assuming that the slant range random error of each control point follows an independent zero-mean Gaussian distribution N 0, σ r 2 , I N represents an N × N identity matrix.
(2) Random error ∆R 2 indirectly caused by control points' positioning: linearize the slant range equation for each target point coordinate.
Assuming that the three-directional random error of each control point follows an independent zero-mean Gaussian distribution N 0, σ t 2 .
The estimated error of the slant range equation can be obtained:

Doppler Equation Position Estimation Error Analysis
In the second step of solving the Doppler equation, the error source ∆ f d1 comes from two parts: the calculation error ∆ f 1 caused by inaccurate control point positioning coordinates, and the error ∆ f 2 indirectly introduced by the inaccurate velocity estimate in the first step of estimation.The two are mutually independent.
Error ∆ f 1 indirectly caused by control points: linearize the Doppler equation for control point coordinates.
(1) Assuming that the three-directional random error of each control point follows an independent zero-mean Gaussian distribution N 0, σ t 2 , the error ∆ f 1 is as follows.
(2) Error ∆ f 2 indirectly caused by inaccurate velocity estimates: linearize the Doppler equation for the three-directional velocity: Simplify the coefficients for ease of expression: Due to the correlation of velocity deviations, expand according to the covariance definition: The cov(∆V x ), cov ∆V y in the equation represents the velocity estimate covariance from the first-step slant range equation error analysis.Finally, according to Equation (51), the estimate error of the Doppler equation for the position can be calculated:

.3. Doppler Equation Velocity Estimation Error Analysis
In the third step of the Doppler equation, the error source ∆ f d2 also comes from two parts: the calculation error ∆ f 1 caused by inaccurate control point positioning coordinates (which is the same as Equation (47)), and the error ∆ f 3 indirectly introduced by the inaccurate position estimate in the second step of the estimation.The two are mutually independent.
Linearize the platform position parameters for the first order: The covariance expression is as follows: Finally, according to Equation (55), the estimate error of the Doppler equation for the velocity can be calculated: In summary, since both the slant range equation and the Doppler equation estimate all six parameters, each of the three-step equations produce their own estimate errors, and the final estimate error depends on the smallest of the two values:

Error Propagation Curve
A simulation experiment was designed to conduct a single-factor analysis of various errors and observe the impact of various errors on the parameters of SAR platform motion in the nadir and side-looking modes.The system parameters are shown in Table 1.The control points were uniformly selected at 200 m intervals in the imaging area.Figure 4 shows the diagram.In Figure 5, due to the same error propagation coefficient in Equation (33), the threedimensional positioning error of the systematic control points introduces an equivalent positioning translation error in the three directions of the SAR platform position, and the three error propagation curves are the same straight line.
In Figure 6, the azimuth position y 0 is not affected by the systematic slant range error, while the range position x 0 is considerably affected, and the vertical position z 0 is affected the least.
In Figure 7, the positions of x 0 , z 0 , v y are not affected by the Doppler bias, while the azimuth position y 0 , the vertical velocity v z , and the range velocity v x are considerably affected, while the others are hardly affected.
Figures 8 and 9 show the error propagation curves of random errors.It can be seen from Figure 8 that the random control points positioning error has different degrees of influence on the six parameters, with the greatest impact on the vertical position z 0 and velocity v z .
According to Figure 9, the random slant range error affects the estimation of the three directions of position and the vertical velocity to varying degrees, and has little effect on the horizontal velocity.

System Parameter Effects on Positioning Error
When the SAR system operates under different imaging parameters, the same input error can result in different positioning errors.A single-factor analysis was conducted on two parameters, platform flight altitude and imaging center distance, to investigate their impacts on the positioning error.In the experiments in this section, the error settings in Table 2 were used.As shown in Figure 10, with the increase in the imaging center distance, the positioning errors in range and azimuth directions x 0 , z 0 increase rapidly, while the y 0 , v x , v y , v z error increases relatively slowly.From Figure 11, when the platform flight altitude is increased, the vertical position z 0 and velocity v z error decrease rapidly, while the range and azimuth position x 0 , y 0 errors increase slightly.This indicates that appropriately increasing the platform altitude can improve the positioning accuracy.
By keeping the distribution of ground-control points unchanged and changing the number of control points, it can be seen from Figure 12 that the positioning errors of all parameters decrease rapidly at first and then tend to be flat in a similar pattern, indicating that appropriately increasing the number of control points can improve the positioning.

Semi-Physical Simulation Verification
Actual SAR images were used for verification, as shown in Figure 13.The left half shows the SAR image used and the right half shows the corresponding optical reference image.The SRAWG feature description method proposed by the author was used for heterogeneous image-matching [20].The feature points of the image were evenly distributed in the image by block partitioning [21], and the matching accuracy could reach the pixel level by using a dense feature description.The yellow dots in the figure are the control points obtained after image-matching, and after obtaining the corresponding feature point pairs in the images, we transformed the SAR image to the geographic coordinate system using a bilinear interpolation, in order to obtain the geographical coordinates of the control points.The imaging parameters of the system are shown in Table 3, and the Scenarios of semi-physical simulation are shown in Figure 14.During the process of heterogeneous image-matching, it was necessary to register the images to the geographic coordinate system and establish the correspondence between the focused images and the geographical coordinates.In this process, we employed a bilinear interpolation.
The systematic offset and random error of each error source in the simulation were set as shown in Table 4.If the platform state was directly solved using the range and Doppler equations, the solution could not converge due to the ill-posed problem of the equations.The three-step method proposed in this paper was applied to solve the platform position and velocity, and 1000 Monte-Carlo experiments were conducted to statistically analyze the error of the parameter estimation.The results were compared with the theoretical error analysis values shown in Tables 5-7.The theoretical error analysis model is described in Sections 3.1 and 3.2.The experiment showed that this method could achieve robust and high-precision estimations of SAR platform motion parameters.Under the influence of typical systematic errors, the estimation accuracy of the distance position of the SAR platform was about 6.5 m, and the estimation accuracy of the azimuth position was about 17.6 m, most of which was affected by systematic errors.The estimation accuracy of the vertical position was about 22.5 m, which was more affected by random errors than by systematic errors.As for the velocity estimation accuracy, the estimation accuracy of the range velocity was better than 0.03 m/s, mainly affected by systematic errors.The estimation accuracy of the azimuth velocity was the highest, about 0.005 m/s, and the error mainly came from random errors.The estimation accuracy of the vertical velocity was slightly lower, about 0.2 m/s.Comparing the estimation error with the theoretical error predicted in Sections 3.1 and 3.2, the two were highly consistent, further verifying the effectiveness and accuracy of the error analysis model derived in this paper.

Discussion
Based on the theoretical analysis and experimental results, the estimation in the vertical direction was the most susceptible to random errors.In the curve of random control point position error and random slant distance error transmission, the derivatives of the vertical position and velocity curves were both maximum values.From the perspectives of systematic and random errors, the former mainly affected horizontal positioning, while the latter mainly affected vertical positioning.From the perspective of system parameters, increasing the flight altitude and shortening detection distance were beneficial for improving positioning accuracy, especially in the vertical direction.In terms of control point selection, a certain number of high-reliability control points should be chosen instead of a large number of low-precision control points, as the benefits of increasing the control point quantity are diminishing.The abovementioned analysis was validated in the semi-physical simulation.
In the future work, we will construct a complete simulation model to integrate the platform positioning method into the entire process loop of navigation.We will attempt to assist INS navigation in the absence of GNSS, and further explore practical testing methods.

Figure 2 .
Figure 2. Convergence curves of parameter estimation errors.

Figures 5 -
Figures 5-7 show the error propagation curve of systematic errors.The solid line represents the position estimation error and the dashed line represents the velocity estimation error.The red, green, and blue colors represent the x, y, z directions, respectively.

Figure 5 .
Figure 5. Systematic control points positioning error propagation curve.

Figure 8 .
Figure 8. Random control points positioning error propagation curve.

Figure 10 .
Figure 10.Curve of positioning error with respect to the center distance.

Figure 11 .
Figure 11.Curve of positioning error with respect to the flight altitude.

Figure 12 .
Figure 12.Curve of positioning error with respect to the number of control points.

Figure 13 .
Figure 13.Actual SAR and corresponding optical images in the region of interest.

Table 3 .
Semi-physical simulation system parameters.

Table 5 .
Platform positioning error introduced by systematic errors.

Table 6 .
Platform Positioning Error Introduced by Random Errors.

Table 7 .
Platform Positioning Error Introduced by the Combination of Two Types of Errors.