Interferometric Calibration Based on a Constrained Evolutionary Algorithm without Ground Control Points for a Tiangong-2 Interferometric Imaging Radar Altimeter

Abstract

The interferometric imaging radar altimeter (InIRA), mounted on the Tiangong-2 space laboratory, utilizes a small incidence and a short interferometric baseline to achieve altimetry for wide swathes of ocean surface topography and inland water surface elevation. To obtain a high-precision digital elevation model (DEM), calibration of the interferometric system parameters is necessary. Because InIRA utilizes the small-incidence interference system design, serious coupling occurs between the interferometric parameters. Commonly used interferometric calibration methods tend to fall into the local optimal solution for InIRA. Because evolutionary algorithms have a stronger robustness and global search ability, they are better suited to handling the solution space structure under the coupling of complex interferometric parameters. This article establishes an interferometric calibration optimization model for InIRA by utilizing the relative flatness of the lake surface as an inequality constraint. Furthermore, an adaptive penalty coefficient constraint evolutionary algorithm is designed to solve the model. The proposed method was tested on actual InIRA data, and the results indicate that it efficiently adjusts interferometric parameters, enhancing the precision of measurements for Qinghai Lake elevation.

1. Introduction

A radar altimeter is an active remote sensing system that measures the distance between a platform and the Earth’s surface by transmitting electromagnetic pulses to the ground, receiving reflected echoes, and calculating the time delay between the transmitted and received signals [1]. Spaceborne radar altimeters are primarily used for ocean observation. They obtain the sea surface echo waveform and precisely track the half-power point of the echo waveform to determine the accurate distance from the satellite platform to the mean sea level. The height of the mean sea level can then be determined using precise orbit determination data from the satellite. Over the past thirty years, numerous altimeter satellites have been launched globally, with the accuracy of data observation increasing from meter to centimeter levels [1]. Satellites such as Seasat [2], Geosat [3], Jason-3 [4], and Sentinel-3 [5] have been launched in the past five decades for oceanography purposes.

Although traditional ocean altimeters have reached centimeter-level measurement accuracy, they are limited by narrow swath and low spatial resolution. Altimeters that use interferometry technology can achieve high spatial resolution and wide swath when measuring the sea level, which is crucial for observing mesoscale and submesoscale ocean phenomena. The Surface Water and Ocean Topography (SWOT) mission [6,7], launched in December 2022, will provide important observational data for precise measurements of global surface water and ocean surfaces. The KaRIN interferometric altimeter is an important payload on the SWOT satellite. It is a Ka-band radar interferometer, including two 5 m-long radar antennas on the left and right. The observation width of SWOT can reach 120 km, and it can obtain two-dimensional height information at the same time, which is beyond the reach of traditional altimeters.

Calibrating the interferometric system parameters with external calibration methods is necessary to enhance the accuracy of the elevation data. Several research efforts have focused on InSAR calibration [8]. The primary aim of these studies is to analyze different error parameters and develop an error model that impacts elevation data accuracy. This error model is subsequently applied to correct the digital elevation model (DEM) product and enhance DEM data accuracy. The baseline is a crucial interferometric calibration parameter among other interferometric system parameters. It affects not only the removal of the earth’s flat phase but also the conversion between terrain phase and elevation, as well as the precision of surface deformation measurements. Various methods have been developed to enhance the accuracy of the baseline and DEM products. These methods include the traditional orbit estimation method [9,10], high-precision ground control point (GCP) method [11,12,13,14], and the fringe change rate estimation method [15]. Selecting an appropriate method to estimate or correct the baseline is essential to guarantee precise elevation data.

The GCP method is the most accurate method for estimating baselines, albeit it is also the most complicated in terms of calculation and implementation. This method necessitates precise ground elevation data and the resolution of equations that combine the baseline and control points. In contrast to other baseline estimation algorithms, the GCP method offers superior baseline estimation accuracy and can fulfill the demands of high-precision terrain mapping. Nevertheless, deploying ground control points in the survey area must satisfy specific prerequisites, and there is no assurance that appropriate control points fulfilling the requirements exist globally. Therefore, investigating interferometric parameter calibration in regions without ground control points is essential. Some researchers have examined natural landscapes as a means of calibrating InSAR system parameters without ground control points [16], while others have used flat inland regions to calibrate airborne InSAR interferometric parameters [17]. However, this method must be implemented when the parameter error range is small; otherwise, it results in a consistent elevation error.

The primary method of calibrating interferometric parameters using ground control points involves solving the interferometric parameters through a sensitivity matrix [8] and optimization equation [18]. The sensitivity matrix method separates the influence of interferometric parameters on the elevation error, expressing the error caused by each interferometric parameter error source linearly. The interferometric parameter sensitivity matrix is then established, and the interferometric parameter error is solved using linear least squares. However, the approximating linearization errors between the interferometric parameters and the elevation error can result in certain errors in the results. Additionally, the sensitivity matrix’s condition number can be too large, making the inversion of the matrix unstable. Another commonly used optimization method minimizes the objective equation using traditional optimization algorithms, reducing the elevation difference between the generated elevation and the ground control point, ensuring the interferometric parameters meet the required equation’s requirements. And traditional optimization algorithms often struggle with multidimensional searches due to the strong coupling between interferometric parameters during the solution, resulting in local optimal solution.

In order to obtain interferometric parameters with a high accuracy and improve the optimization capability of the interferometric calibration algorithm, this paper proposes a novel interferometric calibration method without ground control points. This method first establishes optimization equations based on the relative flatness of the terrain and solves them using the intelligent evolutionary algorithm. Then, an adaptive penalty coefficient is introduced to dynamically adjust the optimization range of the algorithm, which adapts its size during the search process. In experiments, the proposed algorithm is validated using InIRA data. Compared with several commonly used algorithms, it demonstrates advantages in terms of accuracy.

2. InIRA Three-Dimensional Reconstruction and Elevation Error Analysis

2.1. InIRA Three-Dimensional Reconstruction Model

The InIRA, developed by the National Space Science Center, Chinese Academy of Sciences, is a wide swath near-nadir incidence cross-orbit InSAR borne on the Tiangong-2 space laboratory. The baseline length between its two Ku-band antennas is 2.3 m. To improve surveying accuracy and increase the surveying width, the InIRA incidence is designed to obtain a swath with a width of over 40 km [6]. Compared with an altimeter, InIRA can obtain a larger swath width, and compared with synthetic aperture radar (SAR) and interferometric SAR, InIRA has higher resolution. At the near-nadir incidence of InIRA, the water surface echo is dominated by specular reflection energy, making the data acquired by InIRA remarkable due to the near-nadir incidence and large incidence variation range. Traditional synthetic aperture radars with large incidences usually show low echo and signal-to-noise ratio (SNR) in water areas, appearing as dark areas in radar images. However, InIRA data present a higher SNR and coherence in the water area due to the small incidence, allowing the radar to receive reflected water surface echoes. And dual antennas are used to obtain highly coherent sea surface echoes for surveying and mapping.

The basic principle of the three-dimensional reconstruction model of InIRA is to obtain the three-dimensional coordinates of the ground target point by combining the distance equation, Doppler equation and interferometric phase equation. Figure 1 shows the geometric relationship of InIRA measurement in the geocentric coordinate system [19], where $\overrightarrow{S_1}$, $\overrightarrow{S_2}$ is the position vector of the master and slave antenna phase center, $\overrightarrow{R_1}$, $\overrightarrow{R_2}$ represents the position vector of the main and auxiliary antenna phase center to the ground point, $\overrightarrow{P}$ is the ground target point position vector and $\overrightarrow{B}$ is the baseline vector. In the triangle composed of the target point $\overrightarrow{P}$, the origin $\overrightarrow{O}$ of the geocentric space coordinate system and the phase center $S_1$ of the main antenna, the vector expression of the target point in the geocentric space rectangular coordinate system can be obtained. The three-dimensional reconstruction model is composed of the distance Equation (1), radar Doppler Equation (2) and interferometric phase Equation (3):

$$\overrightarrow{P}=\overrightarrow{S_1}+\overrightarrow{R_1}$$

(1)

$$f_{dop}=-\frac{2\overrightarrow{v}\left(\overrightarrow{S_1}-\overrightarrow{P}\right)}{\lambda R_1}$$

(2)

$$\phi =\varphi +2k\pi =\frac{2Q\pi \left(R_1-R_2\right)}{\lambda }$$

(3)

where $R_1$ is the distance from the target point to the phase center of the main antenna; $R_2$ is the distance from the auxiliary antenna to the target point; $\overrightarrow{v}$ is the satellite velocity vector; $f_{dop}$ is the Doppler center frequency; $\lambda$ is the radar wavelength; $\phi$ is the absolute interferometric phase; $\varphi$ is the wrapped phase; k is the integer cycle number of the interferometric phase; Q is the antenna mode and 2 and 1, respectively, represent the repeated orbit interferometric mode and the double antenna interferometric mode.

To obtain an analytical solution to the above equation and analyze the interferometric sensitivity of elevation to various parameters while minimizing the approximation error, the paper introduces a moving coordinate system with the phase center of the main antenna serving as the origin. The unit view vector is then orthogonally projected onto this coordinate system. Subsequently, the paper achieves the transformation of the unit view vector from the moving coordinate system to the imaging coordinate system through coordinate transformation, facilitating terrain reconstruction. Based on the InIRA geometric relationship, the target point’s coordinates are expressed as follows:

$$\overrightarrow{P}=\overrightarrow{S_1}+R_1·T_{vpq}·{\widehat{r}}_{vpq}$$

(4)

where ${\widehat{r}}_{vpq}$ represents the unit vector of the view vector $\overrightarrow{R_1}$ in the direction of v axis, p axis and q axis; $T_{vpq}$ is the conversion matrix from the view vector in the VPQ coordinate system to the Cartesian coordinate system in the geocentric space. The VPQ coordinate system takes the phase center of the main antenna as the origin of the coordinates, the direction of the velocity vector is the v axis, the cross product direction of the velocity and the baseline is the q axis and the v axis, p axis and q axis satisfy the right-handed coordinate system relationship. The visual vector decomposition method converts the visual vector from the moving coordinate system to the geocentric space rectangular coordinate system, so as to obtain the three-dimensional position information of the target point:

$$\overrightarrow{P}=\overrightarrow{S_1}+R_1·\left[\begin{array}{ccc}\frac{\overrightarrow{v}}{v}\hfill & \frac{(\overrightarrow{v}\otimes \overrightarrow{B})\otimes \overrightarrow{v}}{|(\overrightarrow{v}\otimes \overrightarrow{B})\otimes \overrightarrow{v}|}& \left.\frac{\overrightarrow{v}\otimes \overrightarrow{B}}{|\overrightarrow{v}\otimes \overrightarrow{B}|}\right]\end{array}\right.·\left[\begin{array}{c}\frac{\lambda f_{dop}}{2v}\\ \frac{R_1}{2b_{pv}}\left[1+\frac{b^2}{R_1^2}-{\left(1-\frac{\Delta r}{R_1}\right)}^2\right]-\frac{b_v{r}_v}{b_{pv}}\\ \pm \sqrt{1-r_v^2-r_p^2}\end{array}\right]$$

(5)

where $b_v$ is the component of the baseline in the direction of velocity; $b_{pv}$ is the component of the baseline perpendicular to the direction of velocity. In InIRA, the baseline is perpendicular to the velocity direction; that is, $b_v=0$, and the squint angle is 0; that is, $f_{dop}=0$, at this angle, Equation (5) degenerates into a common radar side-looking equation.

2.2. InIRA Elevation Error Analysis

The baseline is an important parameter in the InIRA system. The length and direction of the baseline vector directly affect many factors such as image coherence, altimetry sensitivity, and altimetry error transfer coefficient. The baseline vector error is also one of the most important error sources affecting the measurement accuracy. Refer to the elevation error analysis method of TanDEM-X [20]; the baseline is decomposed into a vertical baseline $\overrightarrow{B_{\perp }}$ perpendicular to the line of sight and a parallel baseline $\overrightarrow{B_{\Vert }}$ along the line of sight. The effects of $\Delta B_{\perp }$, $\Delta B_{\Vert }$ and $\Delta \varphi$ on InIRA elevation error are analyzed below.

InIRA is a spaceborne dual-antenna interferometric system. The main antenna transmits, and the two antennas receive simultaneously. The relationship between the elevation of the ground point and the phase after removing the ground effect is:

$$h=h_{amb}\frac{\varphi }{2\pi }=\frac{\lambda R_1\sin \theta }{B_{\perp }}\frac{\varphi }{2\pi }$$

(6)

where $\varphi$ is the absolute interferometric phase after phase unwrapping, and $h_{amb}$ is the interferometric system ambiguity height, defined as:

$$h_{amb}=\frac{\lambda R_{1}sin\theta }{Bcos\left(\theta -\alpha \right)}=\frac{\lambda R_{1}sin\theta }{B_{\perp }}$$

(7)

where $\theta$ is the incidence, if there is a baseline error $\Delta B_{\perp }$ in the system, then the introduced altimetry error $\Delta h$ is:

$$\Delta h=-h\frac{\Delta B_{\perp }}{B_{\perp }+\Delta B_{\perp }}\approx -h\frac{\Delta B_{\perp }}{B_{\perp }}$$

(8)

Assuming that the terrain height is 3000 m, when $\theta =7^{\circ }$, $B_{\perp }=2.2986\mathrm{m}$, the vertical baseline error $\Delta B_{\perp }=\pm 1$ mm will cause an elevation error of 1.3 m, and the error changes slowly with the increase in incidence in the range direction.

$\Delta B_{\Vert }$ is the main error source of the interferometric system. The baseline error along the line of sight directly affects the interferometric slant range difference. Therefore, $\Delta B_{\Vert }$ and $\Delta \varphi$ are coupled with each other during the calibration process. The elevation error caused by the combination can be expressed as:

$$\Delta h=h_{amb}(\frac{\Delta B_{\Vert }}{\lambda }+\frac{\Delta \varphi }{2\pi })$$

(9)

Because the ambiguity height $h_{amb}$ changes with the range direction, $\Delta B_{\Vert }$ and $\Delta \varphi$ will cause the ground elevation to change along the range direction, which will cause the distance slope. If the distance length is s, the slope is:

$$k=\frac{\Delta h}{\Delta s}=\frac{\Delta B_{\Vert }}{B_{\perp }}+\frac{\Delta \varphi }{2\pi }·\frac{\lambda }{B_{\perp }}$$

(10)

From the above analysis, it can be concluded that the conventional space-borne InSAR elevation error is manifested as a tilt error in the range direction, which is caused by the parallel baseline error and the interferometric phase error, and the relationship between the two is coupled. Since InIRA is different from traditional long-baseline and large-incidence spaceborne interferometric systems, the error of InIRA’s vertical baseline will also have a greater impact on the elevation.

In order to analyze and determine the actual influence of baseline vector and interferometric phase on elevation error, we use the system parameters of a scene InIRA image to analyze the sensitivity of interferometric parameters to elevation based on the elevation reconstruction model. Specifically, we introduce linear errors to the baseline length, baseline inclination, and interferometric phase at the image’s center point while keeping other interferometric parameters constant. The resulting elevation error is displayed in Figure 2.

Figure 2 reveals that the influence of baseline length B, baseline inclination $\alpha$ and absolute interferometric phase $\phi$ on elevation is roughly linear. This is due to InIRA’s smaller incidence, which leads to a linearly coupled error contribution to the elevation, making it challenging to achieve parameter correction by separating the impact of each parameter on elevation. When conventional optimization methods are used to search for a multi-dimensional solution space, coupled interferometric parameters increase the risk of local optimal solution convergence [21,22,23]. Non-convex functions often feature multiple local optimal solutions, which impede algorithms from determining the optimal solution. Moreover, these algorithms heavily rely on gradient information or higher-order derivatives that may limit their ability to escape local optimal solutions. Hence, this paper proposes the use of a constrained optimization evolutionary algorithm in the following section to overcome this issue.

3. Optimization Calibration Method Based on Inland Water Body

3.1. InIRA Optimization Calibration Model

According to the elevation error analysis of interferometric parameters in the previous section, we take the baseline length B, baseline inclination $\alpha$ and absolute interferometric phase $\phi$ as the parameters to be calibrated, which cause slope errors in the InIRA elevation, and the interferometric system quantities that do not participate in the calibration are regarded as constants; the elevation of the surface object is a nonlinear function of the system parameter variable $X=\left[B,\alpha ,\phi \right]$; and the estimation accuracy of the target point is related to the system parameter X; that is, the elevation h can be expressed as a nonlinear function of each system parameter.

The parameter error of the actual InSAR system is $\Delta X=\left(\delta B,\delta \alpha ,\delta \phi \right)$, and each error has a certain error range. Let $\delta B\in \left[B_a,B_b\right]$, $\delta \alpha \in \left[{\alpha }_a,{\alpha }_b\right]$, $\delta \phi \in \left[{\phi }_a,{\phi }_b\right]$; that is, the solution space domain R of the system parameter error is $\Delta X=\left\{\left[B_a,B_b\right],\left[{\alpha }_a,{\alpha }_b\right],\left[{\phi }_a,{\phi }_b\right]\right\}$. Assuming that the initial value of the obtained interferometric system parameter vector is $\widehat{X}$ and the real value is X, for the calibration based on control points, the elevation at the $n\left(n=1,2,\cdots ,N\right)$ control points is set to $h_n$, and the single point objective optimization objective function can be expressed as [18]:

$$F_n\left(\Delta X\right)=f_n\left(\widehat{X}+\Delta X\right)-h_n$$

(11)

where $f_n\left(·\right)$ is the elevation reconstruction formula at the n th calibration point; that is, the calibration problem of the InSAR system is transformed into the problem of adjusting the error parameters to minimize the objective function, so that the elevation value after adjusting the system parameter error is the minimum distance from the elevation value at the control point. According to the interferometric coherence, different control points are weighted. Combined with the least square method to solve the model parameters, we can obtain the optimal objective function under different weights $w_n$:

$$F\left(\Delta X\right)=\underset{n=1}{\sum ^N}{w}_n{\left[f_n\left(\widehat{X}+\Delta X\right)-h_n\right]}^2$$

(12)

Different from the traditional InSAR, the near-nadir interference makes InIRA have high SNR and high correlation on the water surface. Water surfaces exhibit the characteristic of height consistency, which is an inherent property of the data. Therefore, lake elevation information can be used for interferometric calibration, utilizing the discreteness of lake elevation as a constraint term in the calibration equation, as follows:

$$\left\{\begin{array}{c}minF\left(\Delta X\right)={\displaystyle \underset{n=1}{\sum ^N}}{w}_n{\left[f_n\left(\widehat{X}+\Delta X\right)-h_n\right]}^2\\ s.t.f_1\left(\widehat{X}+\Delta X\right)=f_2\left(\widehat{X}+\Delta X\right)=\cdots =f_M\left(\widehat{X}+\Delta X\right)\end{array}\right.$$

(13)

where N is the number of calibration points on the land elevation and M is the number of calibration points on the lake elevation. In fact, due to the influence of coherence difference and random noise, the lake elevations are not completely equal, so we use the statistical information of the lake elevation to assist with calibration, and describe the constraint as the discreteness of the lake elevation, with a penalty coefficient $\xi$ added to the Equation (19), the value of $\xi$ represents the constraint ability of the constraint term on the objective function:

$$fitness(\Delta X)=F(\Delta X)+\xi G(\Delta X)$$

(14)

where $fitness(\Delta X)$ is the fitness function after adding a penalty term and $G(\Delta X)$ is a constraint term that converts equality constraints into inequality constraints, designed as:

$$G(\Delta X)=\sum _{m=1}^M{\left[f_m(\widehat{X}+\Delta X)-\overline{f_m(\widehat{X}+\Delta X)}\right]}^2$$

(15)

3.2. Whale Optimization Algorithm (WOA)

Under the system design of InIRA with small incidence angles, the coupling between the interferometric parameters and the characteristics of the non-convex function makes it challenging for traditional optimization algorithms to find the global optimal solution [24,25]. Genetic algorithm [26] demonstrates its advantages in solving InSAR calibration problems with large incidence angles. However, when facing highly coupled optimization problems, it may encounter difficulties solving these problems due to limited search capability and challenges in parameter settings. However, insufficient search diversity in genetic algorithm can lead to the issue of premature convergence. Therefore, we use the WOA intelligent evolutionary algorithm [27] to calibrate the InIRA system parameters. WOA [28,29] is a meta-heuristic optimization algorithm that simulates the hunting behavior of humpback whales. It uses random or optimal search agents to simulate the hunting behavior and uses spirals to simulate the bubble-net predation mechanism of humpback whales, which can effectively solve complex optimization problems. WOA can perform parallel processing on multi-dimensional interferometric parameters. By searching multiple interferometric parameters at the same time, the entire search space can be better explored, thereby increasing the probability of searching for the global optimal solution. Secondly, the WOA algorithm has adaptability based on the search history. It overcomes the local optimal solution in the search space by adaptively adjusting the search step size and search direction. In the search process, the search direction is adjusted according to the position of the historical optimal solution, increasing the likelihood of finding the global optimal solution. The WOA algorithm uses a spiral search strategy to develop unexplored areas. In the search process, the whale will search near the known optimal solution, but it will also randomly explore other regions, thereby increasing the diversity of the search and helping to avoid falling into the local optimal solution.

In the WOA algorithm [28], each whale represents a potential optimal solution to the extreme value optimization problem. The whales continue to explore and hunt prey in the iterative process of the algorithm, and continue to try to improve. The WOA algorithm first initializes the search range of the interferometric calibration parameters in the feasible solution space as the position of the whale individual, and uses the interferometric optimization calibration model as the fitness function. Each time the whale population updates its position, the fitness is calculated once, and the current optimal fitness value is updated by comparison. By continuously updating the position of the whale in the solution space, the global optimal solution is finally obtained. In the process of optimization, WOA algorithm has three kinds of behavior, which are surround search, random search and spiral search.

3.2.1. Encirclement Search Mode

The population sets the optimal position according to the value of the fitness function, and the optimal solution in D dimensions is expressed as $x_p=[x_p^1,x_p^2,···,x_p^D]$, the solutions at other positions move to the vicinity of the currently known optimal position to achieve the approach to the current optimal solution. At this stage, the update formula of the interferometric parameters position is shown in Equations (16) and (17):

$$x_i^d(k+1)={x_p}^d\left(k\right)-A\times {D_{1,i}}^d$$

(16)

$$D_{1,i}^d=\left|C\times x_{\mathrm{p}}^d\left(k\right)-x_i^d\left(k\right)\right|$$

(17)

where $x_i^d$ represents the position of the i th interferometric parameter in the d-dimensional space, k represents the current iteration number, and A and C are the distance adjustment parameters [28].

3.2.2. Random Search Mode

In the process of searching the solution space, the optimal solution is randomly searched according to the position of other interferometric parameters. The mathematical model of parameter random search is shown in Equations (18) and (19):

$$x_i^d(k+1)=x_{\mathrm{rand}}^d\left(k\right)-A\times D_{2,i}^d$$

(18)

$$D_{2,i}^d=\left|C\times x_{\mathrm{rand}}^d\left(k\right)-x_i^d\left(k\right)\right|$$

(19)

where $x_{\mathrm{rand}}=\left[x_{\mathrm{rand}}^1,x_{\mathrm{rand}}^2,\cdots ,x_{\mathrm{rand}}^D\right]$ is the position vector of random interferometric parameters in the sample. At this stage, $\left|A\right|>1$.

3.2.3. Spiral Search Mode

In addition to the above two search methods, there is a spiral search method whose probability is P, and the probability of choosing the encircling or random method is $1-P$. Equations (20) and (21) are search models in spiral mode:

$$x_i^d(k+1)=x_{\mathrm{p}}^d\left(k\right)+D_{3,i}^d\times {\mathrm{e}}^{bl}\times cos\left(2\pi l\right)$$

(20)

$$D_{3,i}^d=\left|x_{\mathrm{p}}^d\left(k\right)-x_i^d\left(k\right)\right|$$

(21)

where l is a random number in [0, 1], and b is a constant parameter.

The WOA algorithm first randomly initializes a set of interferometric parameters according to the range of interferometric parameters. In each iteration, according to the values of A and p, WOA can switch between three different search methods, search and update the current optimal solution, and end the WOA algorithm by satisfying the termination criterion condition.

3.3. WOA Interferometric Calibration Algorithm Based on Adaptive Penalty Function

To adaptively adjust the random search capability of the WOA algorithm, we designed an adaptive penalty coefficient based on the constraint of lake surface elevation. The information obtained from the evolutionary process of the population is utilized as feedback to dynamically adjust the penalty coefficient. In the early stages of the search, the search range is expanded, and in the later stages, the search is focused within the feasible domain. This approach enhances the robustness of the algorithm by avoiding convergence difficulties caused by large initial errors in interferometric parameters. The penalty fitness function is constructed as follows:

$$fitness\left(x\right)=f\left(x\right)+\xi \left(t\right)G(\Delta X)$$

(22)

$\xi \left(t\right)$ is updated in each generation as follows:

$$\xi (t+1)=\left\{\begin{array}{cc}\left(1/{\beta }_1\right)\xi \left(t\right),& \mathrm{if}\mathrm{case}1\\ {\beta }_2\xi \left(t\right),& \mathrm{if}\mathrm{case}2\\ \xi \left(t\right),& \mathrm{otherwise}\end{array}\right.$$

(23)

where ${\beta }_1>{\beta }_2>1$, case 1 indicates that the optimal parameters found in the past n (defined parameter) generations are within the feasible region, indicating that the previous penalty coefficient is large enough, and $\xi$ can be appropriately reduced to reduce the penalty pressure on the infeasible solution, which is conducive to enhancing the global search ability of the algorithm around the feasible region. And case 2 indicates that the best individuals found in the past n generations are outside the feasible region, indicating that the previous penalty coefficient is too small, and it is necessary to appropriately increase the penalty for infeasible solutions to avoid excessive search outside the feasible region.

In summary, the main steps of the InIRA calibration method based on the WOA algorithm are as follows:

Step 1: Set the number of searched populations N and the maximum number of iterations $tmax$ of the algorithm, and initialize the location information. The settings of N and $tmax$ depend on the number of interferometric calibration parameters and the complexity of the optimization model. A larger population size can increase the algorithm’s search range and diversity, but it also increases computational complexity and memory consumption. Similarly, the maximum number of iterations should be sufficiently large to ensure the algorithm has enough time to conduct a global search and find optimal solutions. However, setting a very large number of iterations may lead to longer execution times. Figure 3 shows the iteration results with population sizes set to 10, 40, and 80. It can be observed that with a smaller population size, the fitness function value is higher compared to the results obtained with larger population sizes, indicating that a smaller population size limits the search capability and leads to us getting trapped in local optima. On the other hand, increasing the population size does not yield significantly better results. Therefore, we choose an appropriate population size and number of iterations. Specifically, in our experiments, we set the population number to 40 and the number of iterations to 200.

Step 2: Calculate the fitness of each search point through the interferometric calibration model, find the position of the current optimal search solution and save it.

Step 3: Calculate parameters a, p and coefficient vectors A, C. Judge whether the probability p is less than 0.5. If it is, go directly to step 4. Otherwise, use Equation (20) to update the position.

Step 4: Determine whether the absolute value of the coefficient vector A is less than 1. If so, update the optimal solution position according to Equation (16); otherwise, search for the optimal solution globally at random, and update the position according to Equation (18).

Step 5: At the end of the position update, the fitness of each interferometric parameter solution is calculated and compared with the previously retained optimal interferometric parameters. If it is better than the previous optimal solution, the new optimal solution is used for replacement. Determine whether the current optimal solution is within the feasible region and use Equation (23) to dynamically update $\xi$.

Step 6: Determine whether the current search number reaches the maximum number of iterations. If it reaches the optimal solution, the calculation is completed. Otherwise, enter the next iteration, and return to step 3.

Figure 4 is the flowchart of the InIRA interferometric calibration algorithm. The orange section represents the improvements we made to the original WOA evolutionary algorithm:

4. Experiments and Analysis

4.1. Interferometric Data Processing in the Experimental Area

The experiment uses the InIRA data borne on the Tiangong-2 space laboratory. As a near-nadir interferometric SAR, the Ku-band dual-antenna radar is used. The main system design parameters are shown in the following

Table 1. InIRA system parameters, including wavelength, baseline length, baseline inclination, nearest slant range, and incidence angle range.

$\mathit{\lambda }$ [m]	B [m]	$\mathit{\alpha }$ [${}^{\circ }$]	${\mathit{R}}_{\mathit{min}}$ [m]	Incidence Range [${}^{\circ }$]
0.0221	2.3m	5	392,081.51	3∼8

, where $R_{min}$ is the nearest slant distance, which is the shortest distance from the InIRA antenna to the imaging area. The data were acquired from the northern region of Qinghai Lake in Qinghai Province, China, in October 2017, and the specific geographical location is shown in the following Figure 5. We first perform basic interferometric processing on the two-channel SAR complex image [30,31]. Figure 6a is a single-channel SAR imaging amplitude image, in which the horizontal axis is the azimuth direction and the vertical axis is the range direction. It can be observed that the near-range direction is the inland area, and the far-range direction is the water body, which is the Qinghai Lake; Figure 6b is the interferometric phase image obtained after the interferometric processing of the registered SAR image, in which the interferometric phase is wrapped between 0 and $2\pi$. It can be observed that the interferometric fringes of the InIRA data on the water surface are similar to the interferometric fringes of the traditional InSAR on the flat ground, showing periodic characteristics. After interferometric phase filtering, Figure 6b is unwrapped by the minimum cost flow to obtain Figure 6c. Due to the short baseline design of InIRA, the interferometric phase winding is not particularly serious, and it can be unwrapped directly without removing the flat-earth phase. Figure 6d is the coherence image of the master and slave SAR images. According to the coherence, the quality of the interference in the imaging area can be judged, and the coherence is used to weight the control points in the interferometric calibration.

4.2. Experimental Design and Procedures

The initial physical baseline length of the InIRA is 2.3 m, and the baseline inclination is 5${}^{\circ }$. This will be used as the initial parameter of the calibration, and the initial elevation is obtained based on the unwrapped interferometric phase and the phase-elevation conversion relationship. In the absence of ground control points, the reference DEM is used to assist the calibration [32,33]. The reference DEM we use is the shuttle radar topography mission (SRTM) DEM with a resolution of 30 m. In our InIRA data, the near-range direction is the land part, which can be calibrated via SRTM DEM, and the far-range direction is the water part without elevation information. Using the high coherence and high SNR characteristics of the InIRA water part, this information is used to compensate for the missing elevation information, thereby realizing the calibration of the interferometric parameters. As shown in Figure 7, in the imaging area, we select multiple SRTM elevation points along the track in the inland part as the external control points. And we select uniformly distributed elevation points on the lake as constraints. After establishing the optimization calibration model of inland water body proposed in this paper, the improved WOA algorithm is used to solve the model.

At the same time, in order to increase the comparability of the experiment, we compare our proposed method with three calibration methods without ground control points: the reference DEM calibration method [18,32], the flat ground calibration method [17] and the flat-earth phase method [8]. Figure 8 shows the differences between our proposed method and the three calibration methods. The green flowchart represents the reference DEM method, which incorporates external reference DEM information to obtain corresponding external control point elevation data. It then utilizes Equation (12) for optimization to obtain the interferometric parameters. The blue flowchart represents the key steps of the flat ground calibration method, which does not involve external information but relies on the elevation characteristics of its own data. Typically, it utilizes the flatness of the terrain to construct the optimization equation and solve for the interferometric parameters. The yellow flowchart represents the flat earth phase method; it utilizes a precise correspondence between the flat earth phase and the interferometric baseline to refine the baseline accurately. Additionally, it employs polynomial fitting to mitigate errors in interferometric measurements caused by terrain deformation.

To verify the proposed method, the high-precision Sentinel-3 data [34,35] in this region is used to evaluate the generated elevation. The Sentinel-3 data are independent of the external control point information we used. They have not been used in our proposed method or in the two additional methods used for comparison. Figure 9 shows the location of the elevation checkpoints on the SAR image. The Sentinel-3 data were acquired on 5 October and 1 November 2017, which closely aligns with the InIRA data acquisition time on 6 October 2017. We sampled and screened 48 elevation points as checkpoints to quantitatively analyze the generated lake elevation and land elevation. The spatial resolution of the sentinel-3 used is 300 m, and the spatial resolution of the InIRA data at the checkpoints at the far-range is about 30 m [6]. The InIRA data were spatially averaged over a $10\times 10$ window to a comparable spatial resolution to that of Sentinel-3 data.

4.3. Results and Analysis

Figure 10 depicts the distribution of generated elevation and elevation error, while

Table 2. Quantitative evaluation of the calibration methods.

Method	Number of Checkpoints	VAR [m]	Average Error [m]	RMSE [m]
Reference DEM	48	30.86	−29.68	30.18
Flat ground	48	0.95	18.31	18.34
Flat-earth phase	48	4.76	−0.97	2.36
Proposed	48	0.95	0.31	1.01

presents the quantitative comparison results between the proposed method and the other three methods. The elevation error observed when using different methods is evaluated based on three criteria: variance (VAR), average error, and root mean square error (RMSE) [36]. The reference DEM method uses inland external control points for calibration. The calibration points cannot cover the entire range direction. Therefore, it can be observed that the elevation error exhibits a slope effect, the elevation errors near the land are smaller, and the elevation error on the lake surface is particularly serious. This is because the elevation error in the area with control points is small, while the elevation error in the lake area without control points is large. The flat ground method can eliminate the range slope to some degree, but the generated elevation will have a fixed constant deviation from the actual elevation. It can be seen from the average value that the fixed elevation error is about 18.31 m, and the error variance is small. In the flat-earth phase method, the average error is −0.97 m, and there are some points with larger errors. The flat-earth phase residual is obtained by using global polynomial fitting for the interferogram, so that the flat-earth phase may absorb some other signals, which is the cause of the error. In our proposed method, the computed baseline length is 2.3359 m, the baseline inclination is 5.0382${}^{\circ }$ and the absolute interferometric phase is precomputed by external DEM, and the calibration value is 0.041 rad. The range slope caused by interferometric parameter errors has been effectively eliminated, and the constant elevation error has returned to near zero value. The average elevation error is 0.31 m, with an RMSE of 1.01 m and a VAR of 0.95 m. According to these three evaluation indicators, our proposed method outperforms the other three methods, demonstrating the effectiveness of our calibration method.

5. Discussion

Unlike traditional InSAR measurements that focus on land surfaces with large incidence angles, InIRA can be used for both land and water body measurements. Our research area comprises both land regions and lake regions. The interferometric calibration calculation based on the adaptive constrained evolutionary algorithm we proposed is suitable for inland areas containing lakes such as InIRA images. In this case, the relative elevation flatness of the lake is used as a constraint and an adaptive penalty coefficient is used. The method we proposed is also applicable to traditional InSAR land mapping areas. At this time, when flat terrain such as plains exists, it can be used as a constraint. If it does not exist, there will be no constraint terms and the evolutionary algorithm will be used to solve it.

From Figure 6, it can be observed that we arranged calibration points on both land and the lake according to Equation (13). The number of calibration points on land is denoted as M, while the number of calibration points on the lake is denoted as N. Following the calibration point arrangement principles in interferometric calibration, calibration points on land are placed to cover the entire range direction as much as possible. If feasible, they are distributed across multiple azimuth directions to enhance stability. In our experiments, the arrangement of M calibration points on land follows the aforementioned principles.

Inspired by the calibration point layout on land, we divided the distribution of N calibration points on the lake into four categories: near-range single-range direction, far-range single-range direction, near-range multi-range direction, and far-range multi-range direction. The specific placement locations are illustrated in Figure 11. Employing these four distinct placement configurations, the proposed method in this paper is applied for processing, and the corresponding results are depicted in Figure 12. Quantitative analyses are also presented in

Table 3. Quantitative evaluation of the four different calibration point placements.

Placement	Number of Checkpoints	VAR [m]	Average Error [m]	RMSE [m]
Near-range single-range	48	30.86	−29.68	30.19
Far-range single-range	48	19.25	30.15	30.46
Near-range multi-range	48	0.95	0.31	1.01
Far-range multi-range	48	0.95	0.31	1.01

.

We devised these four distinct calibration point placements, encompassing a comparison between near-range and far-range placement, as well as the influence of the number of control points—comparing multi-range placement to single-range placement. From Figure 12, it can be observed that placing calibration points along a single slant range on the lake still leads to the presence of the slope effect in elevation errors, with a significant variance. Placing calibration points at different slant ranges on the lake noticeably resolves the slope effect, leading to a reduction in variance and the mean value becoming very close to the validation elevation values. Across the same 48 validation points, the average error drops to 0.31 m, and the RMSE reduces to 1.01 m. These results align with the previous findings from the uniform distribution of points on the lake surface. Conclusions drawn from the results in Figure 12 and Table 3 reveal that single-range calibration points on the lake surface do not positively impact outcomes, regardless of whether they are placed near the near-range or far-range end. Multi-range placement, whether near-range or far-range, proves beneficial and exhibits similar effectiveness to placing control points across the entire lake. For stability, it is advisable to cover the entire lake with as many calibration points as possible.

6. Conclusions

This study addresses the calibration problem of short-baseline interferometric SAR with near-nadir interference. Firstly, we analyze the impact of interferometric parameters on elevation errors and find through system simulation that these parameters have a linear impact on elevation within the error range, and that there is severe coupling among them. Next, we propose a calibration model based on inland water bodies without ground control points due to the high coherence of backscattering signals on water surfaces in InIRA data. The relative flatness of the lake surface is used as an inequality constraint to compensate for the missing information of the lake part without ground control points. To solve this calibration model, we introduce an intelligent evolutionary algorithm. The model’s constraint ability is adaptively adjusted during the iteration process, and the WOA algorithm’s global search ability and wide search range are utilized to avoid local optimal points, thus determining optimal interferometric parameters and achieving high-precision inland water level generation. Through comparisons with several other calibration methods, the effectiveness and advantages of our proposed approach have been demonstrated.

Our study does not consider changes in baseline with azimuth time. Therefore, the next step is to investigate the characteristics of InIRA baseline with azimuth time and develop an accurate model to enable more precise generation of inland and lake elevations.