Airborne Platform Three-Dimensional Positioning Method Based on Interferometric Synthetic Aperture Radar Interferogram Matching

: As the demand for precise navigation of aircraft increases in modern society, researching high-precision, high-autonomy navigation systems is both theoretically valuable and practically significant. Because the inertial navigation system (INS) has systematic and random errors, its output information diverges. Therefore, it is necessary to combine them with other navigation systems for real-time compensation and correction of these errors. The SAR matching positioning and navigation system uses synthetic aperture radar (SAR) image matching for platform positioning and compensates for the drift caused by errors in the inertial measurement unit (IMU). Images obtained by SAR are matched with digital landmark data, and the platform’s position is calculated based on the SAR imaging geometry. However, SAR matching positioning faces challenges due to seasonal variations in SAR images, the need for typical landmarks for matching, and the lack of elevation information in two-dimensional SAR image matching. This paper proposes an airborne platform positioning method based on interferometric SAR (InSAR) interferogram matching. InSAR interferograms contain terrain elevation information, are less affected by seasonal changes, and provide higher positioning accuracy and robustness. By matching real-time InSAR-processed interferograms with simulated interferograms using a digital elevation model (DEM), three-dimensional position information about the matching points has been obtained. Subsequently, a three-dimensional positioning model for the platform has bene established using the unit line-of-sight vector decomposition method. In actual flight experiments using an FMCW Ku-band Interferometric SAR system, the proposed platform positioning framework demonstrated its ability to achieve precise positioning in the absence of signals from the global navigation satellite system (GNSS).


Introduction
Radar-based positioning and navigation methods include one-dimensional terrain contour matching [1,2] and two-dimensional SAR image matching navigation [3,4].Radar, as an additional sensor, is used to address the drift problem caused by errors in the IMU, especially in situations where GNSS signals are lost or interrupted.The principle of digital terrain matching navigation is that as the aircraft flies over certain specific terrain areas along its flight path, it uses radar altimeters to measure the terrain elevation profile along the path.These real-time measurements are then correlated with pre-stored reference maps to determine the aircraft's geographical position.Existing terrain matching methods mainly involve correlating altimeter-measured terrain profiles with those in reference maps, which is a line matching approach that lacks resolution in the cross-track direction, leading to potential mismatches, and requiring areas with typical terrain slope characteristics for matching.
Apart from one-dimensional terrain matching, SAR matching navigation achieves a two-dimensional resolution.It matches real-time radar-acquired images with stored SAR images in its database to correct current positional errors.However, SAR scene matching navigation is affected by seasonal changes in terrain features, requiring typical landmarks for reliable matching performance and is challenging to map in areas with significant terrain variations, leading to poor navigation accuracy in mountainous regions.Researchers have focused on two main areas: one is the matching of heterogeneous SAR images, proposing algorithms tailored to SAR image features to improve matching precision and robustness [5,6].The other involves using the positional information of SAR matching points and SAR's slant range data to establish the relationship between the aircraft and matching points to calculate the aircraft's position [3,7,8].
InSAR uses interferometry to accurately measure terrain elevation and achieve twodimensional imaging.It effectively combines the advantages of terrain matching navigation and the SAR image matching navigation methods.InSAR interferometric matching navigation has different working modes.When the terrain undulations in the navigation area are obvious, the InSAR interferometric matching navigation mode is enabled, and the generated interferograms with terrain undulation characteristics are used for matching to achieve platform positioning.When the terrain within the navigation area is relatively flat, singlechannel SAR images can be used to complete SAR image matching navigation.Therefore, introducing InSAR technology into positioning and navigation systems offers significant advantages and broad application scenarios.The European Defence Agency (EDA) conducted preliminary analyses of the potential of InSAR-assisted navigation [9][10][11], exploring InSAR parameters and demonstrating the feasibility of InSAR-assisted navigation.They also studied the impact of position and attitude angle errors on phase and positioning errors, as well as the characteristics of digital terrain models.Subsequently, researchers established a platform position inversion model using matched interferogram offsets and proposed using interferometric phase and Doppler frequency to infer the platform's attitude angles [12,13].Further, a mountain branch point-based InSAR interferometric fringe matching algorithm was proposed to achieve precise matching between interferometric fringes [14].Later, researchers began to theoretically analyze high-speed platform positioning based on wrapped InSAR interferograms, providing simulation experiments [15].
While InSAR demonstrates advantages in accuracy and robustness for matching positioning and navigation compared to traditional SAR, there has not been a complete methodological workflow with real flight data validation for InSAR interferogram matching positioning and navigation so far.Most studies have been based on theoretical analysis and simulations.Therefore, this paper proposes a comprehensive three-dimensional platform positioning method based on InSAR interferogram matching, establishing a robust and highprecision platform positioning and navigation system.The rest of the paper is organized as follows: Section 2 presents the proposed three-dimensional positioning framework based on InSAR interferogram matching.Section 3 proposes a method for constructing reference interferograms using external source DEM.Section 4 proposes a feature matching algorithm suitable for interferometric fringe images, achieving precise matching between actual interferograms and reference interferograms.Section 5 proposes a three-dimensional positioning model for the platform based on InSAR altimetry geometry.Section 6 conducts actual InSAR airborne flight experiments, obtaining corresponding positioning results and analyzing positioning accuracy.Finally, Section 7 presents the conclusions of the paper.

Platform Positioning Framework
InSAR has the ability to generate terrain elevation models and can match these models with DEM.However, considering real-time data processing efficiency, interferograms rather than digital elevation models are used for matching.Figure 1

Reference Interference Fringe Generation
The preparation of reference fringe maps is a prerequisite for InSAR matching and positioning, and the accuracy of these reference maps affects both the accuracy of the matching and the positioning.Therefore, higher accuracy in the reference fringe maps can effectively improve positioning accuracy.Currently, there is no interferometric fringe database available, so we use a method of inferring interferometric fringes from DEMs to construct an interferometric fringe library [16,17].The key lies in calculating the slant range for each point in the DEM of the imaging area based on the geometry of bistatic SAR imaging and the interferometric phase model.Then, the interferogram can be obtained based on the relationship between the difference in slant range and interferometric phase.In actual airborne SAR and navigation guidance applications, it is necessary to The platform positioning and navigation framework based on InSAR interferogram matching comprises five parts: (1) First, real-time automatic registration imaging is used to generate actual interferograms, referred to as "actual interferograms".(2) Subsequently, based on the geometric principles of InSAR imaging, reference digital elevation models are inverted to simulate interferograms, known as "reference interferograms".For detailed information, please refer to Section 3. (3) Next, a feature matching algorithm is used to match the reference interferograms with the actual interferograms.This process extracts homonymous points between the two interferograms, allowing for the accurate extraction of homonymous point positions from the reference interferograms.For more details, please refer to Section 4. (4) In the following step, we developed a model to invert and analyze the platform's position.This model utilizes the concept of InSAR three-dimensional inverse positioning.Please see Section 5 for more information.(5) Finally, the platform's position and the INS positioning information are combined through a filtering process to output the platform's positioning information, achieving long-term autonomous navigation.

Reference Interference Fringe Generation
The preparation of reference fringe maps is a prerequisite for InSAR matching and positioning, and the accuracy of these reference maps affects both the accuracy of the matching and the positioning.Therefore, higher accuracy in the reference fringe maps can effectively improve positioning accuracy.Currently, there is no interferometric fringe database available, so we use a method of inferring interferometric fringes from DEMs to construct an interferometric fringe library [16,17].The key lies in calculating the slant range for each point in the DEM of the imaging area based on the geometry of bistatic SAR imaging and the interferometric phase model.Then, the interferogram can be obtained based on the relationship between the difference in slant range and interferometric phase.In actual airborne SAR and navigation guidance applications, it is necessary to consider the impact of the squint angle on the imaging area.We have proposed an algorithm for generating reference interferometric phase maps based on DEM.A flowchart for this process is shown in Figure 2.
Remote Sens. 2024, 16, x FOR PEER REVIEW 4 of 22 consider the impact of the squint angle on the imaging area.We have proposed an algorithm for generating reference interferometric phase maps based on DEM.A flowchart for this process is shown in Figure 2.

Determination of Imaging Area
The IMU typically records track information in the East-North-Up (ENU) coordinate system, while the determination of the imaging area needs to be carried out in the body coordinate system of the aircraft.Therefore, it is necessary to transform the East-North-Up coordinate system into the aircraft's body coordinate system.The coordinates ( , , ) x y z represent the body's position in the northeast-up coordinate system.By applying the rotation matrix R, the coordinates in the body coordinate system are obtained as ( , , ) After completing the coordinate system transformation, information on the position of the platform can be determined.The next step is to preliminarily determine the imaging area by combining the imaging parameters and the imaging model.The geometric relationship is shown in Figure 3.The projection point of the aircraft on the ground is denoted as shown below.In the InSAR oblique imaging model, there are offsets along the flight direction and perpendicular to the flight direction.We can use the relationship between the oblique angle and the slant range to calculate the offset for each resolution unit: where h is the elevation of the ground point.Using the projection point O of the IMU on the ground, as well as the offset along the flight direction x Δ and the offset perpendicular to the flight direction y Δ , the imaging area in the reference DEM can be calculated.

Determination of Imaging Area
The IMU typically records track information in the East-North-Up (ENU) coordinate system, while the determination of the imaging area needs to be carried out in the body coordinate system of the aircraft.Therefore, it is necessary to transform the East-North-Up coordinate system into the aircraft's body coordinate system.The coordinates (x, y, z) represent the body's position in the northeast-up coordinate system.By applying the rotation matrix R, the coordinates in the body coordinate system are obtained as (x ′ , y ′ , z ′ ).
After completing the coordinate system transformation, information on the position of the platform can be determined.The next step is to preliminarily determine the imaging area by combining the imaging parameters and the imaging model.The geometric relationship is shown in Figure 3.The projection point of the aircraft on the ground is denoted as shown below.In the InSAR oblique imaging model, there are offsets along the flight direction and perpendicular to the flight direction.We can use the relationship between the oblique angle and the slant range to calculate the offset for each resolution unit: where h is the elevation of the ground point.Using the projection point O of the IMU on the ground, as well as the offset along the flight direction ∆x and the offset perpendicular to the flight direction ∆y, the imaging area in the reference DEM can be calculated.

DEM Resampling
The calculated DEM data are located within the ground range-azimuth plane, while the imaging results of SAR are in the slant range-azimuth plane.Due to terrain undulations, the echo signal from the near-ground point reflecting to the receiver may arrive later than that from the far-ground point.Therefore, it is necessary to resample the DEM based on the size of the slant range, so that the corresponding slant ranges of the DEM are ar-

DEM Resampling
The calculated DEM data are located within the ground range-azimuth plane, while the imaging results of SAR are in the slant range-azimuth plane.Due to terrain undulations, the echo signal from the near-ground point reflecting to the receiver may arrive later than that from the far-ground point.Therefore, it is necessary to resample the DEM based on the size of the slant range, so that the corresponding slant ranges of the DEM are arranged from small to large, as shown in Figure 4:

DEM Resampling
The calculated DEM data are located within the ground range-azimuth plane, while the imaging results of SAR are in the slant range-azimuth plane.Due to terrain undulations, the echo signal from the near-ground point reflecting to the receiver may arrive later than that from the far-ground point.Therefore, it is necessary to resample the DEM based on the size of the slant range, so that the corresponding slant ranges of the DEM are arranged from small to large, as shown in Figure 4: From Figure 4, it can be observed that, due to the undulations in the DEM, the echo from the far point P2 arrives at the receiver earlier than the echo from the near point P1.As a result, the slant range R2 corresponding to the far point is smaller than the slant range R1 corresponding to the near point.After DEM resampling, the DEM in the slant rangeazimuth plane is arranged in ascending order based on the slant range.The resampled DEM can be directly used to generate interferometric fringes using the interferometric geometry model.Figure 5 presents the results of DEM resampling, where Figure 5a shows the initial DEM, Figure 5b shows the DEM after slant range resampling, and Figure 5c shows the difference after resampling.It can be observed that, due to the relatively gentle terrain in this area, the DEM difference after resampling is not significant.From Figure 4, it can be observed that, due to the undulations in the DEM, the echo from the far point P2 arrives at the receiver earlier than the echo from the near point P1.As a result, the slant range R2 corresponding to the far point is smaller than the slant range R1 corresponding to the near point.After DEM resampling, the DEM in the slant rangeazimuth plane is arranged in ascending order based on the slant range.The resampled DEM can be directly used to generate interferometric fringes using the interferometric geometry model.Figure 5 presents the results of DEM resampling, where Figure 5a shows the initial DEM, Figure 5b shows the DEM after slant range resampling, and Figure 5c shows the difference after resampling.It can be observed that, due to the relatively gentle terrain in this area, the DEM difference after resampling is not significant.

Calculate Interferometric Phase and Remove Topographic Contributions
After obtaining the resampled DEM in the imaging area, the interferometric phase of the interferogram can be calculated using the fundamental principles of InSAR to derive the DEM [18]: The above formula is a nonlinear equation to calculate the slant range difference R Δ .It can be iteratively solved using optimization methods.After meeting the accuracy conditions, the iterative results are obtained.For an airborne dual-antenna InSAR system, the interferometric phase is directly obtained from the slant range difference:

Calculate Interferometric Phase and Remove Topographic Contributions
After obtaining the resampled DEM in the imaging area, the interferometric phase of the interferogram can be calculated using the fundamental principles of InSAR to derive the DEM [18]: The above formula is a nonlinear equation to calculate the slant range difference ∆R.It can be iteratively solved using optimization methods.After meeting the accuracy conditions, the iterative results are obtained.For an airborne dual-antenna InSAR system, the interferometric phase is directly obtained from the slant range difference: Here, Q represents the interferometric mode, where Q = 1 corresponds to singlepass interferometry, and Q = 2 corresponds to repeat-pass interferometry.The obtained interferometric phase includes the topographic phase, and there is no phase unwrapping at this stage for efficient real-time navigation matching.To avoid excessive line features caused by phase wrapping, we have removed the flat earth phase, and the flat earth phase is: Here, θ = arccos(H/R), and the symbol ⌊ ⌋ denotes the floor function.The interferometric phase for matching, which includes the wrapped phase, can be expressed as:

Interference Fringe Feature Extraction and Matching
In this section, we need to achieve precise matching between actual interferometric fringes and reference interferometric fringes.In interference fringe matching, due to the influence of drift errors, flight attitude errors, and interference parameter errors, the interference fringe will produce geometric distortion and nonlinear phase errors, so the algorithm based on template matching is not applicable under such conditions.On the other hand, deep learning-based methods are not suitable since there is no relevant interference fringe database.Therefore, we adopt a feature-based image matching algorithm framework.
The classical SIFT algorithm [19] establishes a linear scale space using Gaussian filters.However, it may not preserve significant phase changes and edge information in the interferometric fringes, affecting the detection of feature points.Instead, we use nonlinear diffusion filtering to construct the scale space.When applied to images, non-linear diffusion filtering [20] represents the brightness of the image through the divergence of the flow function.Non-linear diffusion filtering is expressed by a partial differential equation as follows: Here, φ represents the interferometric phase value, t is the diffusion time, div and ∇ are the divergence and gradient operators, and c(x, y, t) is the conductivity function.By choosing an appropriate conductivity function, diffusion can be adapted to the local structure of the image.Perona and Malik proposed a way to construct the conductivity function [21]: Here, ∇φ σ is the gradient of the interferometric phase ∇φ σ after Gaussian smoothing, and the function g is defined as: Here, K is the contrast factor, controlling the diffusion level.A larger value retains less edge information.By solving the non-linear partial differential equation using the Additive Operator Splitting (AOS) algorithm, we obtain images corresponding to different times in the non-linear scale space: Here, I is the identity matrix, and A l represents the derivative along the I direction.After establishing the nonlinear scale space, all pixel points were compared with the surrounding and adjacent pixel points in the two layers above and below.After obtaining feature points, a three-dimensional quadratic function was used for fitting to obtain subpixel positioning results for feature points.The use of nonlinear diffusion filtering to extract feature points in the nonlinear scale space can suppress noise while preserving edge information, thereby increasing the quantity and accuracy of feature points.As shown in Figure 6, a comparison was made between feature points detected in the linear scale space using the SIFT algorithm and those detected using nonlinear diffusion filtering.It can be observed that the SIFT algorithm detects fewer feature points, and edge points and corners on the fringe map are not sufficiently detected under the smoothing effect of Gaussian linear space.In contrast, the feature points detected using the nonlinear diffusion filter retain edge points and corners, resulting in a higher number of detected feature points.After determining the positions of feature points, a selection process is applied using the interferometric coherence coefficient to filter out feature points with low coherence.A descriptor is then constructed to form a subset of features.Once key points are selected in the scale space, these key points, being extrema in the scale space, undergo wavelet transformation [22][23][24][25] to separate multiscale features around the key points, achieving effective feature extraction.Rectangular sub-images are established with key points as centers to describe key point information using sub-region image features, enabling mutual matching between key points in two images.Two layers of 2D discrete wavelet transforms are performed on the sub-region around key points.Since interference fringes contain a significant amount of noise in high-frequency components and texture information in low-frequency components, the first layer's low-frequency components and the second layer's low-frequency components are used for feature extraction.Simultaneously, lowfrequency components obtained through wavelet decomposition possess translation invariance, enhancing the robustness of fringe matching and making it more resistant to noise, distortions, and other interferences.
After removing the flat ground effect from the interference fringes, which reflects terrain undulation information with a gradually changing gradient, analyzing local gradient information around feature points can effectively describe the current feature points.Histogram of oriented gradient (HOG) is a technique used for texture-based image analysis [24][25][26], simplifying images by extracting gradient information.HOG extracts features that have a locally distinctive shape based on edges or gradient structures.Through block operations, HOG makes local geometric and intensity changes easy to control.If the translation or rotation of an image is much smaller than the scale of local spatial and directional units, the detected differences are minimal.Detecting HOG features in the wavelet domain of two heterogeneous images can effectively describe the texture information After determining the positions of feature points, a selection process is applied using the interferometric coherence coefficient to filter out feature points with low coherence.A descriptor is then constructed to form a subset of features.Once key points are selected in the scale space, these key points, being extrema in the scale space, undergo wavelet transformation [22][23][24][25] to separate multiscale features around the key points, achieving effective feature extraction.Rectangular sub-images are established with key points as centers to describe key point information using sub-region image features, enabling mutual matching between key points in two images.Two layers of 2D discrete wavelet transforms are performed on the sub-region around key points.Since interference fringes contain a significant amount of noise in high-frequency components and texture information in lowfrequency components, the first layer's low-frequency components and the second layer's low-frequency components are used for feature extraction.Simultaneously, low-frequency components obtained through wavelet decomposition possess translation invariance, enhancing the robustness of fringe matching and making it more resistant to noise, distortions, and other interferences.
After removing the flat ground effect from the interference fringes, which reflects terrain undulation information with a gradually changing gradient, analyzing local gradient information around feature points can effectively describe the current feature points.Histogram of oriented gradient (HOG) is a technique used for texture-based image analysis [24][25][26], simplifying images by extracting gradient information.HOG extracts features that have a locally distinctive shape based on edges or gradient structures.Through block operations, HOG makes local geometric and intensity changes easy to control.If the translation or rotation of an image is much smaller than the scale of local spatial and directional units, the detected differences are minimal.Detecting HOG features in the wavelet domain of two heterogeneous images can effectively describe the texture information around feature points in both images, making it suitable for matching between two images with rich texture features.
First, the low-frequency approximate image extracted through wavelet transformation is divided into blocks.The gradients G x G y in the x and y directions for each pixel in each block are calculated using gradient templates: For the input image I, performed convolution operations with the template, respectively, to obtain the gradients I x and I y in the x and y directions: Thus, we obtained the magnitude |G| = I x 2 + I y 2 and orientation θ = arctan(I y /I x ) of the gradient.Then, we divided 0-360 • into several ranges, counted the occurrences of gradient directions in different ranges, and formed a histogram of gradient directions.
We connected the histograms of gradient directions for each block region to create the HOG feature.This process resulted in a feature vector for each keypoint's surrounding region in the wavelet domain.The first-level low-frequency components were divided into 4 × 4 cells, and the orientation gradients were calculated for each block region to form a histogram of gradient directions.The second-level low-frequency components were divided into 8 × 8 cells, aiming to increase gradient information around key points and improve matching accuracy.Figure 7 clearly shows the process of descriptor construction.
Remote Sens. 2024, 16, x FOR PEER REVIEW 9 of 22 Thus, we obtained the magnitude Once the descriptor was constructed, the next step involved using the minimum Euclidean distance method for real-time matching of feature vectors between actual and reference interferograms.This process eliminates outlier matches, resulting in pairs of matching points.The three-dimensional position information of the matching points can be obtained by querying latitude, longitude, and height information.The algorithmic flowchart for interferometric fringe feature matching is illustrated in Figure 8, with the Once the descriptor was constructed, the next step involved using the minimum Euclidean distance method for real-time matching of feature vectors between actual and reference interferograms.This process eliminates outlier matches, resulting in pairs of matching points.The three-dimensional position information of the matching points can be obtained by querying latitude, longitude, and height information.The algorithmic flowchart for interferometric fringe feature matching is illustrated in Figure 8, with the technical process on the left and key step results on the right.

Platform Three-Dimensional Positioning Model Based on InSAR Geometry
Through the matching of the interferograms, the corresponding position and elevation information for the tie points was obtained, and we established the inversion model of the platform position as shown in Figure 9.The x-axis represents the flight direction of the platform, the z-axis represents the elevation direction, and the y-axis forms a righthanded coordinate system with the x-axis and z-axis. 1 A and 2 A represent the positions of the master and slave radar antennas, respectively, P represents the matched tie point, and R is the slant range vector from the master antenna to point P, referred to as the look vector.b is the baseline vector, where n b represents the length of the cross-track baseline and v b represents the length of the along-track baseline.Orthogonal unit vectors constitute a moving coordinate system, where v is in the flight direction, n is in the same direction as the cross-track baseline, and ŵ is determined using the right-hand rule.

Platform Three-Dimensional Positioning Model Based on InSAR Geometry
Through the matching of the interferograms, the corresponding position and elevation information for the tie points was obtained, and we established the inversion model of the platform position as shown in Figure 9.The x-axis represents the flight direction of the platform, the z-axis represents the elevation direction, and the y-axis forms a right-handed coordinate system with the x-axis and z-axis.A 1 and A 2 represent the positions of the master and slave radar antennas, respectively, P represents the matched tie point, and R is the slant range vector from the master antenna to point P, referred to as the look vector.b is the baseline vector, where b n represents the length of the cross-track baseline and b v represents the length of the along-track baseline.Orthogonal unit vectors constitute a moving coordinate system, where v is in the flight direction, n is in the same direction as the cross-track baseline, and ŵ is determined using the right-hand rule.In the moving coordinate system, the unit look vector can be represented as: where ˆ, ς =< > r w , and < ⋅ > represents the inner product oper- ation between vectors.Therefore, the representation of the unit look vector in the moving coordinate system can be calculated as [27][28][29]: where dc f is the Doppler centre frequency, v is the platform velocity.The positive or negative sign is determined by the side-looking direction of the radar.The unit look vector is transformed from the moving coordinate system to the carrier coordinate system using the transformation matrix T : where transformation matrix T is given by: where b θ is the angle between the baseline and the x-y plane.Therefore, the three-dimensional position vector of the platform A can be represented as: In the moving coordinate system, the unit look vector can be represented as: where µ =< r, v >, η =< r, n >, ς =< r, ŵ >, and < • > represents the inner product operation between vectors.Therefore, the representation of the unit look vector in the moving coordinate system can be calculated as [27][28][29]: where f dc is the Doppler centre frequency, v is the platform velocity.The positive or negative sign is determined by the side-looking direction of the radar.The unit look vector is transformed from the moving coordinate system to the carrier coordinate system using the transformation matrix T: r = T × rvnw (15) where transformation matrix T is given by: where θ b is the angle between the baseline and the x-y plane.Therefore, the threedimensional position vector of the platform A can be represented as:

Experiments and Results
We validated the proposed method through semi-physical experiments using data acquired by a dual-antenna Ku-band FMCW InSAR system developed by the Institute of Electronics, Chinese Academy of Sciences.The InSAR system was mounted on an unmanned aerial vehicle (UAV).The InSAR system is shown in Figure 10, and Table 1 provides the airborne InSAR flight parameters.
We validated the proposed method through semi-physical experiments u acquired by a dual-antenna Ku-band FMCW InSAR system developed by the I Electronics, Chinese Academy of Sciences.The InSAR system was mounted manned aerial vehicle (UAV).The InSAR system is shown in Figure 10, and Ta vides the airborne InSAR flight parameters.The radar system is equipped with two position and orientation syste POS610 and Mini-POS.Mini-POS is a lightweight and compact system with lo racy.We conducted matching and positioning experiments using the low-prec from Mini-POS, which has an accelerometer drift error of 0.1 mg and a gyros error of 1°/h.Our method was validated using high-precision positioning POS610.To improve the accuracy of matching and positioning, we utilized ligh and ranging (LiDAR) DEM data with a grid resolution of 1 m as a reference DE   The radar system is equipped with two position and orientation systems (POS): POS610 and Mini-POS.Mini-POS is a lightweight and compact system with lower accuracy.We conducted matching and positioning experiments using the low-precision IMU from Mini-POS, which has an accelerometer drift error of 0.1 mg and a gyroscope drift error of 1 • /h.Our method was validated using high-precision positioning data from POS610.To improve the accuracy of matching and positioning, we utilized light detection and ranging (LiDAR) DEM data with a grid resolution of 1 m as a reference DEM.

Generation of Actual Interferograms
Real-time generation of InSAR interferograms mainly involves dual-antenna SAR imaging, self-focusing, registration, interferometric processing, and ground removal.These processing techniques are well-established, and we present the interferogram processing steps and results for our data directly.Figure 13a shows the interferogram before

Generation of Actual Interferograms
Real-time generation of InSAR interferograms mainly involves dual-antenna SAR imaging, self-focusing, registration, interferometric processing, and ground removal.These processing techniques are well-established, and we present the interferogram processing steps and results for our data directly.Figure 13a shows the interferogram before registration, Figure 13b displays the interferogram after fine registration, Figure 13c illustrates the interferogram after filtering, and Figure 13d presents the interferogram after removing the flat Earth phase.

Interferogram Matching and Platform Localization Results
During the flight, due to interruptions in GNSS signals, there was no real-time correction applied to the IMU, leading to drift errors.Therefore, the purely inertial navigation solution using Mini-POS causes the position error to increase over time.As a result, using its localization information for navigation introduces position offsets, causing misalignment between actual interferograms and reference interferograms.The results obtained after applying our proposed interferogram feature matching algorithm are shown in Figure 14, where the left side depicts actual generated interferograms, and the right side shows reference interferograms.And, after using the ransac algorithm to remove outliers, the correct matching point pair was obtained.

Interferogram Matching and Platform Localization Results
During the flight, due to interruptions in GNSS signals, there was no real-time correction applied to the IMU, leading to drift errors.Therefore, the purely inertial navigation solution using Mini-POS causes the position error to increase over time.As a result, using its localization information for navigation introduces position offsets, causing misalignment between actual interferograms and reference interferograms.The results obtained after applying our proposed interferogram feature matching algorithm are shown in Figure 14, where the left side depicts actual generated interferograms, and the right side shows reference interferograms.And, after using the ransac algorithm to remove outliers, the correct matching point pair was obtained.The reference interferogram provides queryable location and elevation information for each matching point, allowing us to obtain three-dimensional position information for the matching points.Utilizing the three-dimensional localization model proposed in Section 5, we calculated the three-dimensional position of the platform.Figure 15 illustrates the differences between the positions obtained by the InSAR matching positioning system and the high-precision POS610 system at each matching point.It can be observed that the x-axis error is within 3 m, the y-axis error is within 7 m, and the z-axis error is within 2 m, resulting in a total positioning error within 8 m.The maximum error in the y-axis is attributed to the imprecise roll angle of the platform.In subsequent experiments, we plan to calibrate the interferometric parameters using the three-dimensional information of the matching points to achieve higher positioning accuracy.The minimum positioning error is observed in the z-axis, highlighting the advantages of interferometric technology in elevation measurement.
To conduct a more in-depth analysis of the data, we employed statistical indicators such as maximum error, average error, and RMSE to evaluate positioning accuracy.As shown in Table 2, the three-dimensional positioning using InSAR interferograms has a maximum error of 7.49 m, an average error of 6.50 m, and an RMSE of 6.51 m.The positioning accuracy restricted to being within 10 m demonstrates that the algorithm achieved high precision in matching navigation.It is worth noting that our experiment is a semiphysical simulation, utilizing a high-precision LiDAR DEM as a reference DEM.In the The reference interferogram provides queryable location and elevation information for each matching point, allowing us to obtain three-dimensional position information for the matching points.Utilizing the three-dimensional localization model proposed in Section 5, we calculated the three-dimensional position of the platform.Figure 15 illustrates the differences between the positions obtained by the InSAR matching positioning system and the high-precision POS610 system at each matching point.It can be observed that the x-axis error is within 3 m, the y-axis error is within 7 m, and the z-axis error is within 2 m, resulting in a total positioning error within 8 m.The maximum error in the y-axis is attributed to the imprecise roll angle of the platform.In subsequent experiments, we plan to calibrate the interferometric parameters using the three-dimensional information of the matching points to achieve higher positioning accuracy.The minimum positioning error is observed in the z-axis, highlighting the advantages of interferometric technology in elevation measurement.
shown in Table 2, the three-dimensional positioning using InSAR interferograms has a maximum error of 7.49 m, an average error of 6.50 m, and an RMSE of 6.51 m.The positioning accuracy restricted to being within 10 m demonstrates that the algorithm achieved high precision in matching navigation.It is worth noting that our experiment is a semiphysical simulation, utilizing a high-precision LiDAR DEM as a reference DEM.In the next steps, we plan to implement a specific hardware system for corresponding positioning experiments.

Discussion
In this paper, we introduced the application of interferometry techniques to the method of positioning for airborne SAR platforms, utilizing interferograms for matching and positioning.This approach combines the advantages of traditional SAR images and terrain matching positioning methods, leveraging the long-term elevation stability provided by DEMs generated through InSAR mapping, as well as the three-dimensional positioning capabilities of InSAR, thereby improving the final positioning accuracy to be within 10 m.The high precision of InSAR fringe matching is primarily due to the additional information in the elevation dimension provided by InSAR compared to SAR, enhancing the positioning capability from two dimensions to three dimensions.On the other hand, InSAR fringes reflect the undulating nature of the terrain, offering more stability compared to SAR images and thus providing higher stability in the matching process.The error of the results is explained below and some issues are discussed.

Positioning Error Analysis
The errors in InSAR interference fringe matching navigation mainly include match- To conduct a more in-depth analysis of the data, we employed statistical indicators such as maximum error, average error, and RMSE to evaluate positioning accuracy.As shown in Table 2, the three-dimensional positioning using InSAR interferograms has a maximum error of 7.49 m, an average error of 6.50 m, and an RMSE of 6.51 m.The positioning accuracy restricted to being within 10 m demonstrates that the algorithm achieved high precision in matching navigation.It is worth noting that our experiment is a semi-physical simulation, utilizing a high-precision LiDAR DEM as a reference DEM.In the next steps, we plan to implement a specific hardware system for corresponding positioning experiments.

Discussion
In this paper, we introduced the application of interferometry techniques to the method of positioning for airborne SAR platforms, utilizing interferograms for matching and positioning.This approach combines the advantages of traditional SAR images and terrain matching positioning methods, leveraging the long-term elevation stability provided by DEMs generated through InSAR mapping, as well as the three-dimensional positioning capabilities of InSAR, thereby improving the final positioning accuracy to be within 10 m.The high precision of InSAR fringe matching is primarily due to the additional information in the elevation dimension provided by InSAR compared to SAR, enhancing the positioning capability from two dimensions to three dimensions.On the other hand, InSAR fringes reflect the undulating nature of the terrain, offering more stability compared to SAR images and thus providing higher stability in the matching process.The error of the results is explained below and some issues are discussed.

Positioning Error Analysis
The errors in InSAR interference fringe matching navigation mainly include matching errors, platform attitude angle errors, and interference system parameter errors.The interference parameter errors can be eliminated after calibration, so we will not go into details here.Therefore, this section discusses the matching error and platform attitude angle error and their impact on platform positioning, and conducts an error analysis based on the actual parameters of airborne InSAR.

Platform Attitude Angle Error
Figure 16 illustrates the positioning error in a plane under the presence of yaw angle error.Since the attitude angle mainly affects the accuracy of plane positioning, an analysis of positional offsets in the azimuthal and range directions is conducted.In the diagram, the aircraft flies along the ideal trajectory.The actual trajectory and the beam illumination area when there is a certain yaw angle error, denoted as θ y , the actual are indicated in red.The platform flies along the azimuthal direction x-axis, range direction y-axis, and vertical direction z-axis, with H representing the platform's flying height.Neglecting terrain fluctuations, the actual slant range matches the ideal trajectory's slant range, denoted as R 0 .The initial look angle is θ, the actual look angle is θ 1 , and the slant angle is θ sq .The positioning offset is represented as X g .The slant range R 0 , initial range positioning distance G, X g , and θ sq are derived from the geometric relationships in the diagram.

Remote
Therefore, the positioning errors in the azimuthal direction ∆x and in the range direction ∆y are: From the equation, we can determine the influence of yaw angle error on positioning offsets.Using the airborne flight parameters provided in this paper for simulation, with a platform flying height of 1898.86 m and a yaw angle range set to [−0.2, 0.2] degrees, and considering the actual downward viewing angle, which ranges from 30 • to 60 • from near slant range to far slant range.The results are shown in Figure 17.It can be observed that the azimuthal positioning error is greatly influenced by yaw angle error, and this influence becomes more pronounced as the downward viewing angle increases.However, the impact on range direction positioning is minimal.Figure 18 depicts the conceptual illustration of positioning errors in a plane under the presence of pitch angle error.Similar to the analysis method for yaw angle, the aircraft flies along the ideal trajectory.The actual trajectory and the beam illumination area when there is a certain pitch angle error, denoted as p θ , are depicted in red.The platform flies along the azimuthal direction X-axis, range direction Y-axis, and vertical direction Z-axis, with H representing the platform's flying height.Initially, under the ideal conditions, the slant range is denoted as 0 R , while after introducing yaw angle error, it is represented as Figure 18 depicts the conceptual illustration of positioning errors in a plane under the presence of pitch angle error.Similar to the analysis method for yaw angle, the aircraft flies along the ideal trajectory.The actual trajectory and the beam illumination area when there is a certain pitch angle error, denoted as θ p , are depicted in red.The platform flies along the azimuthal direction x-axis, range direction y-axis, and vertical direction z-axis, with H representing the platform's flying height.Initially, under the ideal conditions, the slant range is denoted as R 0 , while after introducing yaw angle error, it is represented as R. The initial downward viewing angle is θ, the actual downward viewing angle is θ 1 , and the slant viewing angle is θ sq .The positioning offset is represented as X g .
The slant range R 0 , initial range positioning distance G, and X g are derived from the geometric relationships depicted in the diagram.
along the azimuthal direction X-axis, range direction Y-axis, and vertical direction Z-axis, with H representing the platform's flying height.Initially, under the ideal conditions, the slant range is denoted as 0 R , while after introducing yaw angle error, it is represented as R .The initial downward viewing angle is θ , the actual downward viewing angle is 1 θ , and the slant viewing angle is sq θ .The positioning offset is represented as g X .
Ideal  Therefore, the positioning errors in the azimuthal direction, represented by ∆x, and in the range direction, represented by ∆y, are: From the equation we can determine the influence of pitch angle error on positioning offsets.Using the airborne flight parameters provided in this paper for simulation, with a platform flying height of 1898.86 m and a yaw angle range set to [−0.2, 0.2] degrees, and considering the actual downward viewing angle, which ranges from 30 • to 60 • from the near slant range to the far slant range.The results are shown in Figure 19.It can be observed that pitch angle mainly affects azimuthal positioning error and remains relatively constant regardless of changes in the downward viewing angle.The impact on range positioning is minimal.
Therefore, the positioning errors in the azimuthal direction, represented by ∆x, and in the range direction, represented by ∆y, are:

Fringe Matching Error
Matching errors mainly stem from two aspects: the accuracy of the reference map and the accuracy of the matching algorithm.The spatial resolution and elevation accuracy of the reference map are critical factors determining its quality.A lower spatial resolution implies that each pixel in the image represents a larger actual ground area, which may result in the inability to accurately identify smaller ground features during the matching process, thereby increasing positioning errors.Additionally, if the reference map exhibits geometric distortions, such as shape distortions caused by imaging angles or differences in ground heights, it directly affects the accuracy of positioning.These geometric distortions may cause the matching algorithm to fail to correctly align the reference map with the actual observed image, introducing additional positioning errors.Moreover, inaccuracies in feature point or region extraction during the matching algorithm's process can

Fringe Matching Error
Matching errors mainly stem from two aspects: the accuracy of the reference map and the accuracy of the matching algorithm.The spatial resolution and elevation accuracy of the reference map are critical factors determining its quality.A lower spatial resolution implies that each pixel in the image represents a larger actual ground area, which may result in the inability to accurately identify smaller ground features during the matching process, thereby increasing positioning errors.Additionally, if the reference map exhibits geometric distortions, such as shape distortions caused by imaging angles or differences in ground heights, it directly affects the accuracy of positioning.These geometric distortions may cause the matching algorithm to fail to correctly align the reference map with the actual observed image, introducing additional positioning errors.Moreover, inaccuracies in feature point or region extraction during the matching algorithm's process can also lead to misalignment and matching errors.
According to the geometric positioning model, the matching error transfer relationship in the x direction, y direction, and z direction is: That is, the matching error and the reference DEM error directly affect the positioning error.Therefore, improving the accuracy of the matching algorithm and using highprecision reference maps can increase the positioning accuracy in equal proportions.

High-Precision Needs and Future Prospects
The generation of interference fringes is a basic and critical step.It is worth mentioning that in the actual navigation process, the system is a continuous iterative process.After the satellite signal is rejected, the InSAR matching navigation process is started in time, and combined filtering with the INS is performed to obtain the corrected position and attitude information.And imaging, matching, and positioning are continuously performed to reasonably control the interval, accurate position, and attitude information.In addition, even in the absence of precise positioning information to achieve high-precision imaging, image generation and interference processes can be achieved through self-focusing algorithms to obtain interference fringes.
In InSAR matching positioning, the interferometric baseline parameters have a significant impact on the generation of interferograms and the accuracy of three-dimensional positioning, leading to considerable nonlinear distortion in the phase of the reference interferogram and significant errors in the three-dimensional positioning.Therefore, calibration of the interferometric parameters is necessary to ensure high positioning accuracy.Our method involves using matching points on the reference DEM as control points and optimizing the interferometric parameters using the least squares method to ensure the accuracy of these parameters.
Our next step is to optimize the efficiency of the program and adopt parallel processing to meet the needs of real-time operation on airborne platforms, aiming to implement a complete positioning and navigation system.
illustrates a work flow of navigation based on InSAR interferogram matching.

( 4 )Figure 1 .
Figure 1.InSAR interferogram matching navigation framework.The orange boxes contain the main technical content.

Figure 1 .
Figure 1.InSAR interferogram matching navigation framework.The orange boxes contain the main technical content.

Figure 3 .
Figure 3. Determination of the reference mapping area based on the aircraft imaging geometry.

Figure 3 .
Figure 3. Determination of the reference mapping area based on the aircraft imaging geometry.

Figure 3 .
Figure 3. Determination of the reference mapping area based on the aircraft imaging geometry.

Figure 4 .
Figure 4. DEM slant range resampling schematic, with the left image showing the DEM before resampling and the right image showing the DEM after resampling.

Figure 4 .
Figure 4. DEM slant range resampling schematic, with the left image showing the DEM before resampling and the right image showing the DEM after resampling.

Figure 5 .
Figure 5. Results before and after DEM slope range resampling.(a) Original DEM.(b) DEM after slant range resampling.(c) Difference in DEM after resampling.

Figure 5 .
Figure 5. Results before and after DEM slope range resampling.(a) Original DEM.(b) DEM after slant range resampling.(c) Difference in DEM after resampling.

Figure 6 .
Figure 6.Comparison of feature point detection methods , where feature points are represented in yellow.(a) Detecting feature points using the SIFT algorithm.(b) Detecting feature points in nonlinear scale space.

Figure 6 .
Figure 6.Comparison of feature point detection methods, where feature points are represented in yellow.(a) Detecting feature points using the SIFT algorithm.(b) Detecting feature points in nonlinear scale space.

Figure 9 .
Figure 9. Three-dimensional position inversion model of the platform.

Figure 9 .
Figure 9. Three-dimensional position inversion model of the platform.

Figure 11
Figure 11 depicts the imaging area determined by the reference trajectory o for the reference DEM.The imaging area is highlighted within the red box.F shows the range image after range resampling, Figure 12b displays the compu map, Figure 12c illustrates the interferogram with the removed ground phase, a 12d presents the reference interferogram with added wrapped phase for match

Figure 11 Figure 11 .
Figure 11  depicts the imaging area determined by the reference trajectory of the IMU for the reference DEM.The imaging area is highlighted within the red box.Figure12ashows the range image after range resampling, Figure12bdisplays the computed phase map, Figure12cillustrates the interferogram with the removed ground phase, and Figure12dpresents the reference interferogram with added wrapped phase for matching.Remote Sens. 2024, 16, x FOR PEER REVIEW

Figure 11 .
Figure 11.The reference DEM determines the imaging area based on the IMU reference trajectory.The imaging area is marked within the red box area.

Figure 11 .Figure 12 .
Figure 11.The reference DEM determines the imaging area based on the IMU reference trajectory.The imaging area is marked within the red box area.

Figure 12 .
Figure 12.Processing steps for the reference interferogram.(a) The reference DEM undergoes slant range resampling.(b) Interferometric phase map.(c) Interferometric phase map after removing the flat Earth phase.(d) Interferogram with wrapped phase.
registration, Figure13bdisplays the interferogram after fine registration, Figure13cillustrates the interferogram after filtering, and Figure13dpresents the interferogram after removing the flat Earth phase.

Figure 15 .
Figure 15.Positioning error of each axis and total positioning error.(a) Positioning error for each axis, (b) Total positioning error.

Figure 15 .
Figure 15.Positioning error of each axis and total positioning error.(a) Positioning error for each axis, (b) Total positioning error.

Figure 17 .
Figure 17.Yaw angle error and positioning offset relationship.

Figure 17 .
Figure 17.Yaw angle error and positioning offset relationship.

Figure 18 .
Figure 18.Effect of pitch angle error on positioning.The slant range 0 R , initial range positioning distance G , and

Figure 18 .
Figure 18.Effect of pitch angle error on positioning.

Figure 19 .
Figure 19.The influence of pitch angle on azimuth and range positioning.

Figure 19 .
Figure 19.The influence of pitch angle on azimuth and range positioning.

Table 2 .
Positioning error statistical results of different axes.

Table 2 .
Positioning error statistical results of different axes.