Signal Processing and Waveform Re-Tracking for SAR Altimeters on High Mobility Platforms with Vertical Movement and Antenna Mis-Pointing

Abstract

Synthetic aperture radar (SAR) altimeters can achieve higher spatial resolution and signal-to-noise ratio (SNR) than conventional altimeters by Doppler beam sharpening or focused SAR imaging methods. To improve the estimation accuracy of waveform re-tracking, several average echo power models for SAR altimetry have been proposed in previous works. However, these models were mainly proposed for satellite altimeters and are not applicable to high-mobility platforms such as aircraft, unmanned aerial vehicles (UAVs), and missiles, which may have a large antenna mis-pointing angle and significant vertical movement. In this paper, we propose a novel semi-analytical waveform model and signal processing method for SAR altimeters with vertical movement and large antenna mis-pointing angles. A new semi-analytical expression that can be numerically computed for the flat pulse response (FSIR) is proposed. The 2D delay–Doppler map is then obtained by numerical computation of the convolution between the proposed analytical function, the probability density function, and the time/frequency point target response of the radar. A novel delay compensation method based on sinc interpolation for SAR altimeters with vertical movement is proposed to obtain the multilook echo, which can optimally handle non-integer delays and maintain signal frequency characteristics. In addition, a height estimation method based on least squares (LS) estimation is proposed. The LS estimator does not have an analytical solution, and requires iterative solving through gradient descent. We evaluate the performance of the proposed estimation strategy using simulated data for typical airborne scenarios. When the mis-pointing angles are within 10 degrees, the normalized quadratic error (NQE) of the proposed model is less than 10⁻¹⁰ and the RMSE of τ obtained by the re-tracking method fitted by the proposed model is less than 0.2 m, which indicates the high applicability of the model and accuracy of the re-tracking method.

1. Introduction

Altimeters play a vital role in the effective operation of onboard platforms. The altitude data they provide supports critical functions such as maintaining stable flight levels, and plays an integral role in implementing advanced automated features. This includes the execution of autonomous return to launch points and automatic takeoff and landing processes. For instance, state-of-the-art consumer unmanned aerial vehicles (UAVs) utilize the Internal Measurement Unit (IMU) and Global Positioning System (GPS) to measure height. However, both sensors can only estimate the flight altitude relative to their starting position, and not the actual altitude over current terrain or vegetation [1]. Using radar sensors is advantageous, as they can directly sense range and velocity and are not affected by lighting conditions or contrast [2]. By integrating high-accuracy radar altimeters, the altitude measurement capabilities of onboard platforms can be significantly enhanced, enabling them to adapt more proficiently to dynamic flight conditions. This integration allows for the execution of intricate navigational tasks with heightened precision and dependability. Consequently, the importance of airborne altimeters is irrefutable in UAV operations [1,3] and aircraft takeoff and landing [4,5,6,7].

Radar altimeters can determine the distance between the radar and the ground surface. Pulse radar altimeters function by emitting a series of short pulses towards the ground and then measuring the time it takes for these pulses to be reflected back to the altimeter. The distance from the surface is determined by the round-trip delay of radar signals. The issue of high-accuracy altimetry has been extensively explored. The rough estimation of height can be obtained through tracking processing. Re-tracking is a crucial procedure applied in processing radar altimeter data to get more accurate height parameters [8]. It is an algorithm used to refine the estimated parameters of a radar echo, such as the time delay, power, and shape, which have been distorted by various factors during the signal’s round-trip travel from the altimeter to the ground surface and back. Re-tracking methods mitigate these distortions to obtain more accurate and reliable measurements. This is achieved by adjusting the parameters of the waveform model, which is representative of the idealized or expected shape of the radar echo. The re-tracking of altimeter echoes is commonly performed using statistical inference-based algorithms such as the Least Squares (LS) method [9,10,11,12], Maximum Likelihood Estimation (MLE) algorithm [13], and Bayesian inference [14,15]. The accuracy of altitude measurements from radar altimeters can be significantly improved through re-tracking.

To improve the estimation accuracy of waveform re-tracking, a series of average echo power models for SAR altimetry have been proposed. Brown made a series of assumptions about the scattering scenario of a flat sea surface in the 1970s, suggesting the seminal three-fold convolution form of the radar altimeter echo model. This model states that the echo power equals the convolution of the FSIR, point target response (PTR), and scatterer height density function (PDF) [16]. Based on the Brown model, a series of semi-analytical and analytical models have been developed to accommodate different scenarios [13,17,18]. SAR altimetry harnesses the Doppler properties of high pulse repetition frequency-transmitted pulses to enable a higher spatial resolution through along-track aperture synthesis [19]. SAR altimetry can mitigate speckle noise compared to conventional altimetry (CA) by permitting multiple observations of the same scatterers through synthetic aperture processing. The resulting shape of the multilook echo diverges from that of the CA echo, necessitating a new altimetric signal model. Numerous studies have introduced numerical models suited to SAR altimeters [20]. If SAR altimetry adheres to the same five assumptions as Brown’s model, the convolution model of traditional radar altimeters can be directly extended to it [21]. In analytical modeling of SAR altimeters, the sinc function obtained by coherent along-track accumulation is generally regarded as the sub-beam pattern after beam sharpening. Then, the sub-beam pattern is multiplied by the antenna pattern to establish the echo model according to the CA [16]. This extends the average flat pulse response in the average echo model to a two-dimensional description of the illumination response of different Doppler frequencies of sub-beams [21,22,23]. A rectangular function is used to approximate the sub-beam pattern; a semi-analytical model of the SAR altimeter echo was proposed based on the geometric approximation of the delay–Doppler bin in [22]. Their model was later enhanced by incorporating the effect of the antenna mis-pointing angle [12]. Furthermore, a Gaussian function was employed to approximate the sub-beam pattern, and the echo model included small antenna mis-pointing angles [23]. A closed-form expression for the SAR altimeter waveform was formulated in [24] by expressing the waveform in terms of a set of basis functions. An innovative processing scheme known as fully focused SAR (FF-SAR) was proposed in [25] to maximize the altimeter’s along-track resolution to the theoretical limit of half the antenna length, proving valuable in SAR imaging and radar altimetry. A new analytical derivation of a delay–Doppler map model for FF-SAR has been proposed, including vertical wave particle motion and a more general surface representation [26]. Unfocused SAR (UF-SAR) retracking methods (e.g., the physical SAMOSA model) have been utilized to retrack fully-focused SAR (FF-SAR) waveforms in [27], ultimately improving data quality for SAR altimetry in coastal regions. A novel coastal retracking algorithm for SAR altimetry was presented in [28] that employed an adaptive interference masking scheme to mitigate interfering signals from reflective coastal targets, significantly increasing the number of valid records in the coastal zone without compromising the quality of the estimated significant wave height.

Airborne altimeters frequently encounter skewed flight conditions and larger antenna mis-pointing angles, which is due to several reasons. Airborne platforms such as aircraft and drones are subject to constant changes in flight dynamics and conditions. Unlike satellites, which have a relatively stable trajectory in space, airborne vehicles often have to adjust their course and altitude to accommodate weather conditions, air traffic, and mission objectives. Such adjustments can lead to vertical movement. In addition, the maneuvers required for aircraft to follow specific survey paths or avoid obstacles can result in off-nadir pointing, leading to significant antenna mis-pointing angles. It must be emphasized that existing research has predominantly focused on altimeters in horizontal flight scenarios, and as such is not suitable for airborne scenarios. For instance, satellite altimeters such as CryoSat-2 generally consider mis-pointing angles of less than 1 degree barring specific maneuvering instances [29]. However, these angles can greatly exceed 1 degree during oblique flights due to substantial pitch and roll angles. Furthermore, mechanical vibrations and the inherent instability of aircraft flight can contribute to greater mis-pointing angles. Despite advancements in flight control systems, it remains challenging to maintain perfect alignment of the radar beam with the vertical axis [30]. All these factors are inherent to airborne operations, and present substantial challenges for high-accuracy altimetry based on echo data from airborne altimeters.

In this paper, we present a strategy for high-accuracy altimetry using airborne altimeters. We propose a generalized semi-analytical echo model for the re-tracking process of SAR altimetry with consideration of vertical movement and large antenna mis-pointing angles. This model encompasses a generalized analytical expression for the FSIR, which accounts for Doppler shift effects and distortions in illumination geometry induced by the flight path angle as well as for pronounced mis-pointing angles in oblique flight scenarios. The Gaussian approximation of the antenna mis-pointing angles is considered. The resulting delay–Doppler map (DDM) obtained by three-fold convolution is governed by six altimetric parameters: the epoch $\tau$, significant wave height, Pu, flight path angle, and on-track and cross-track mis-pointing angles. Compared to satellite altimeters, the lower altitude of airborne altimeters results in a narrower Doppler beam, leading to a decrease in the precision of traditional delay compensation methods. A novel delay compensation method based on sinc interpolation for SAR altimeters with vertical movement is proposed to obtain the multilook echo, which can optimally handle non-integer delays and maintain the signal frequency characteristics. The validity of the proposed model is ascertained through the performance of the widely adopted least squares (LS) algorithm for altimetric parameter estimation [22,31]. We propose a five-parameter estimation strategy in which the antenna mis-pointing angle and flight path angle are updated with each iteration. The effectiveness of the proposed model and algorithm is demonstrated in ideal conditions, with the results confirming our model’s ability to accurately estimate the altimetric parameters.

2. Numerical Mean Echo Power Model for SAR Altimetry

Unlike CAs, which only perform power accumulation in range resolution bins, SAR altimeters transform the echo signal into the range–Doppler domain using the slow-time Fourier transform. The proposed model provides the semi-analytic expression of the mean signal power in each range–Doppler bin. As demonstrated in Figure 1, each range–Doppler bin is associated with two ground areas; the signal power can be expressed as the summation of the scatterers’ echo power from the two corresponding areas.

In this section, we derive a mean echo power model for an SAR altimeter with antenna mis-pointing and vertical movement. The altimetry scenario is introduced in Section 2.1. We first derive the echo signal from a single point target in Section 2.2, followed by pulse compression of the signal in Section 2.3 and the along-track FFT of the compressed signal in Section 2.4. Finally, the final numerical echo model of all the scatterers on the surface is provided in Section 2.5.

2.1. Altimetric Scenario

The altimetric scenario and coordinate system are shown in Figure 1, where the radar altimeter moves in the y–z plane at constant velocity v and incident angle $\mu$. Here, the x–y plane corresponds to the flat surface illuminated by radar antenna, which can be considered flat for airborne applications, while the z axis is provided by the line from the nadir point on the ground to the radar antenna. Thus, the location of the radar at time t is provided by

$$x_r\left(t\right)=0,y_r\left(t\right)=vtcos\mu ,z_r\left(t\right)=h-vtsin\mu ,$$

(1)

where h is the altitude when $t=0$.

Due to platform’s movement, the radar antenna may not look vertically downward towards the nadir point, and its beam directions are defined by angle $\xi$ with respect to the z axis and angle $\tilde{\varphi }$ with respect to the x axis. For simplification, we consider a Gaussian approximation of the radar antenna gain in the form

$$G\left(\theta \right)=G_0{e}^{-(2/\gamma ){sin}^2\theta },$$

(2)

where $\theta$ is the angle relative to the antenna’s boreside direction and $\gamma$ determines the antenna’s beamwidth, which is provided by ${\theta }_{3\mathrm{dB}}=2{sin}^{-1}\sqrt{(\gamma /2)log2}$.

2.2. Radar Signal Model

In this section, we derive the echo signal from a single scatterer located at point $(x,y,z)$ on the flat surface. To begin with, the radar transmit waveform is provided by

$$s\left(t\right)=\sum _{n=0}^{N-1}{s}_0(t-n{T}_r)e^{j2\pi f_{0}t},$$

(3)

where $s_0\left(t\right)=1/\sqrt{T}rect(t/T)exp\left\{j2\pi f_{0}t+j\pi K{t}^2\right\}$ is the baseband waveform in each pulse, $f_0$ denotes the carrier frequency, $T_r$ denotes the pulse repetition interval (PRI), and N is the number of pulses. Denoting the distance from the radar to the scatterer at time t by $R\left(t\right)$, the echo signal from this scatterer after down-conversion can be written as

$$s_r\left(t\right)=\sum _{n=0}^{N-1}{A}_0{s}_0\left(t-n{T}_r-\frac{2R\left(t\right)}{c}\right)exp\left(-j4\pi f_0\frac{R\left(t\right)}{c}\right),$$

(4)

where c is the propagation velocity, $A_0$ is a complex amplitude factor, and

$$R\left(t\right)=\sqrt{x^2+{(y-vtcos\mu )}^2+{(h-z-vtsin\mu )}^2}.$$

Let $t=t_r+n{T}_r$, where n is pulse index and $t_r$ denotes the radar’s fast time. Thus, the received signal during the n-th PRI is

$$S_r\left(t_r,n\right)=s_r\left(t_r+n{T}_r\right)=A_0{s}_0\left(t_r-\frac{2R\left(t_r+n{T}_r\right)}{c}\right)exp\left(-j4\pi f_0\frac{R\left(t_r+n{T}_r\right)}{c}\right)$$

(5)

for $n=0,\dots ,N-1$ and $0\le t_r\le T_r$.

2.3. Pulse Compression

We obtain the range profile of the scatterers via pulse compression. The signal after pulse compression in the n-th PRI is provided by

$$\begin{array}{cc}\hfill C\left(t_r,n\right)& ={\int }_{-\infty }^{\infty }{S}_r\left(\tau ,n\right)s_0^{*}\left(\tau -t_r\right)\mathrm{d}\tau \hfill \\ \hfill & \approx A_0{W}_r\left(t_r-\frac{2R\left(t_r+n{T}_r\right)}{c}\right)exp\left(-j4\pi f_0\frac{R\left(n{T}_r\right)}{c}\right),\hfill \end{array}$$

(6)

where $W_r\left(·\right)$ is the autocorrelation function (ACF) of $s_0\left(t\right)$.

Before proceeding any further, applying approximations to the distance $R\left(t\right)$ can be beneficial for simplifying the expression of $C\left(t_r,n\right)$. For a small N and narrow antenna beam, we have

$$R\left(t\right)\approx R_0-z-\frac{v\left(ycos\mu +hsin\mu \right)}{R_0}t=R_0-z-v_{r}t,$$

(7)

where $R_0=\sqrt{x^2+y^2+h^2}$ and $v_r=v\left(ycos\mu +hsin\mu \right)/R_0$ denotes the relative velocity of the scatterer with respect to the radar. Substituting (7) into the phase history in $C\left(t_r,n\right)$ and assuming that the variation of $R\left(t\right)$ during a coherent processing interval exceeds the range resolution of radar, we have

$$C\left(t_r,n\right)\approx A_1{W}_r\left(t_r-\frac{2\left(R_0-z\right)}{c}\right)exp\left(j4\pi f_0(v_r/c)n{T}_r\right).$$

(8)

Here, the phase terms that are independent to n and $t_r$ are merged into the new factor $A_1=A_{0}exp(j4\pi f_0(-R_0+z)/c)$. From (8), the Doppler frequency is provided by

$$f_d=2v_r{f}_0/c=\frac{2v{f}_0(ycos\mu +hsin\mu )}{R_{0}c}.$$

(9)

2.4. Along-Track FFT

SAR altimeters can achieve higher along-track resolution through Doppler beam sharpening, which is accomplished using the inter-pulse Fourier transform. For convenience, we assume that N is an even number. Following discrete Fourier transformation, the signal output is

$$\begin{array}{cc}\hfill \gamma \left(t_r,k\right)& =\frac{1}{N}\sum _{n=0}^{N-1}C\left(t_r,n\right)exp\left(-j\frac{2\pi }{N}nk\right)\hfill \\ \hfill & =A_2{W}_r\left(t_r-\frac{2\left(R_0-z\right)}{c}\right)W_a(k-N{f}_d{T}_r)\hfill \end{array}$$

(10)

for $k=-N/2,\dots ,N/2-1$. Here, $A_2=A_{1}exp\left[j\left(N-1\right)\pi \left(f_d{T}_r-k/N\right)\right]$ and the point target response $W_a(·)$ in the Doppler domain is provided by

$$W_a\left(k^{\prime }\right)=\frac{1}{N}\frac{sin\left(\pi k^{\prime }\right)}{sin(\pi k^{\prime }/N)}\approx \mathrm{sinc}\left(k^{\prime }\right).$$

(11)

2.5. Backscattered Waveform Power Model

Following [6,16,24], the average backscattered power from a flat surface can be obtained by an integral of the echo signal power over all the scatterers on the surface. For a rough surface with roughness characterized by a height probability density $p\left(z\right)$, the average power in the delay–Doppler bin $\left(t_r,k\right)$ is

$$P_I\left(t_r,k\right)={\int \int \int }_{{\mathbb{R}}^3}p\left(z\right)\frac{{\lambda }^2{G}^2\left(\theta \right){\sigma }_0}{{\left(4\pi \right)}^3{R}_0^4{L}_p}PTR\left(t_r-2\left(R_0-z\right)/c,k-N{f}_d{T}_r\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z$$

(12)

when $z\ll h$ for $k=-N/2,\dots ,N/2-1$ and $0\le t_r\le T_r$. Here, $\lambda$ denotes the wavelength, ${\sigma }_o$ is the backscattering cross-section per unit scattering area, and $L_p$ is the two-way propagation loss. The point target response is

$$PTR(t_r,k)={\left|W_r\left(t_r\right)W_a\left(k\right)\right|}^2.$$

In (12), $G\left(\theta \right)$, $R_0$, and $f_d$ all depend on x and y, while ${\sigma }_0$ is assumed to be independent of x and y for simplify.

Generally, it is not possible to provide a closed form result for the integration in (12). Nevertheless, a semi-analytic expression can be provided via the convolution model proposed in [16], as shown in the next section.

3. Semi-Analytical Echo Model of Airborne SAR Altimeters

In this section, we provide a semi-analytical echo model of airborne SAR altimeters through an approximation that reformulates $P_I\left(t_r,k\right)$ into the following convolution:

$$P_I\left(t_r,k\right)=\frac{c}{2}{p}^{\prime }\left(t_r\right)\ast PTR(t_r,k)*FSIR(t_r,k).$$

(13)

Here, $p^{\prime }\left(t_r\right)=p(-c{t}_r/2)$ and $FSIR(t_r,k)$ is the impulse response from a flat surface, which is defined in (17). The derivation of the convolution model is first provided in Section 3.1. Then, an approximate semi-analytical expression of $FSIR\left(t_r,k\right)$ is provided in Section 3.2. The echo power after three-term convolution requires delay compensation to obtain multilook waveforms; thus, a novel delay compensation method is provided in Section 3.3. The multilook waveform model used for airborne SAR altimeters with antenna mis-pointing and vertical movement is characterized by six altimetric parameters: the significant wave height (SWH), epoch ($\tau$), amplitude ($P_u$), flight path angle ($\mu$), and on-track and cross-track mis-pointing angles (${\psi }_{ac}$ and ${\psi }_{al}$, respectively). Finally, the effects of the mis-pointing angles on the multilook echoes are analyzed in Section 3.4.

3.1. Convolution Model

To obtain the convolution model in (13), we can first consider a flat surface with a small scale of roughness, namely, $z\approx 0$, for all the scatterers. Then, the average power in the delay–Doppler bin $\left(t_r,k\right)$ can be provided by

$$P_{FS}\left(t_r,k\right)={\int \int }_{{\mathbb{R}}^2}\frac{{\lambda }^2{G}^2\left(\theta \right){\sigma }_0}{{\left(4\pi \right)}^3{R}_0^4{L}_p}PTR(t_r-2R_0/c,k-N{f}_d{T}_r)\mathrm{d}x\mathrm{d}y$$

(14)

for $k=-N/2,\dots ,N/2-1$ and $0\le t_r\le T_r$. We define $p^{\prime }\left(t_r\right)=p(-c{t}_r/2)$. It was proposed in [16] that $P_I\left(t_r,k\right)$ can be expressed as a convolution between $p^{\prime }\left(t_r\right)$ and $P_{FS}\left(t_r,k\right)$:

$$P_I\left(t_r,k\right)=\frac{c}{2}{\int }_{\mathbb{R}}{P}_{FS}\left(t_r-\tau ,k\right)p^{\prime }\left(\tau \right)\mathrm{d}\tau .$$

(15)

As it is difficult to determine the analytical expression of $P_{FS}\left(t_r,k\right)$, we can consider an ideal point target response that is a delta function in the time domain and rectangular in the Doppler domain, namely, $PT{R}^{\prime }(t_r,k)=\delta \left(t\right)I\left(k\right)$ and

$$I\left(k\right)=\left\{\begin{array}{cc}1,& \left|k\right|\le 1/2,\hfill \\ 0,& \left|k\right|>1/2.\hfill \end{array}\right.$$

(16)

For $PT{R}^{\prime }(t_r,k)$, $P_{FS}\left(t_r,k\right)$ reduces to the impulse response from the flat surface

$$FSIR\left(t_r,k\right)={\int \int }_{{\mathbb{R}}^2}\frac{{\lambda }^2{G}^2\left(\theta \right){\sigma }_0}{{\left(4\pi \right)}^3{R}_0^4{L}_p}\delta \left(t_r-2R_0/c\right)I\left(k-N{f}_d{T}_r\right)\mathrm{d}x\mathrm{d}y$$

(17)

for $k=-N/2,\dots ,N/2-1$ and $0\le t_r\le T_r$.

The convolution between $FSIR\left(t_r,k\right)$ and the actual point target response is

$$\begin{array}{cc}\hfill & {\tilde{P}}_{FS}\left(t_r,k\right)=\sum _{k^{\prime }=-\infty }^{+\infty }{\int }_{\mathbb{R}}FSIR\left(t_r-\tau ,k-k^{\prime }\left)PTR\right(\tau ,k^{\prime }\right)\mathrm{d}\tau \hfill \\ \hfill & =\sum _{k^{\prime }=-\infty }^{+\infty }{\int \int }_{{\mathbb{R}}^2}\frac{{\lambda }^2{G}^2\left(\theta \right){\sigma }_0}{{\left(4\pi \right)}^3{R}_0^4{L}_p}PTR(t_r-2R_0/c,k^{\prime })I(k-k^{\prime }-N{f}_d{T}_r)\mathrm{d}x\mathrm{d}y.\hfill \end{array}$$

(18)

Using staircase approximation on the point target function

$$PTR(t_r-2R_0/c,k)\approx \sum _{k^{\prime }=-\infty }^{+\infty }PTR(t_r-2R_0/c,k^{\prime })I(k-k^{\prime }),$$

we have ${\tilde{P}}_{FS}\left(t_r,k\right)\approx P_{FS}\left(t_r,k\right)$. Combining (15) and (18), we can finally write $P_I\left(t_r,k\right)$ as a convolution of the height probability density, point target response, and impulse response from the flat surface, as shown in (13).

The corresponding FSIR is obtained by integrating over each Doppler beam. The k-th Doppler beam for time instant t, as shown in Figure 2, is defined by an integral angle $\varphi$ varying within the range $D_{t_r,k}=\left[{\varphi }_{t_r,k},{\varphi }_{t_r,k+1}\right]\cup \left[{{\varphi }^{\prime }}_{t_r,k},{{\varphi }^{\prime }}_{t_r,k+1}\right]$, leading to

$${\varphi }_{t_r,k}=arcsin\left(y/\rho \left(t_r\right)\right)=\arcsin \left(\frac{f_d\lambda R_0+2v_{c}hsin\mu }{2v_c\rho \left(t_r\right)cos\mu }\right),$$

(19)

where $-2v\left(\rho cos\mu +hsin\mu \right)/\lambda R_0\le f_d\le -2v\left(hsin\mu -\rho cos\mu \right)/\lambda R_0$; here, the parameters $xy$ in $f_d$ are omitted for brevity. Note that in the absence of antenna mis-pointing the two delay–Doppler bins exhibit symmetry along the track, resulting in equal power contributions to the DDM. However, their contributions differ when mis-pointing is considered, resulting in

$${\varphi }_{t_r,k}=\pi -{{\varphi }^{\prime }}_{t_r,k+1}\mathrm{and}{\varphi }_{t_r,k+1}=\pi -{{\varphi }^{\prime }}_{t_r,k}.$$

(20)

In SAR altimetry, the integral area of each Doppler beam varies due to the Doppler effect. The integration interval corresponding to the k-th Doppler beam is denoted as $D_{t_r,k}$, resulting in

$$FSIR\left(t_r,k\right)=\frac{{\lambda }^2{\sigma }^0}{{\left(4\pi \right)}^3{L}_p}{\int }_{{\mathbb{R}}^{+}\times D_{t_r,k}}\frac{\delta \left(t_r-2R_0/c\right)G^2\left(\rho ,\varphi \right)}{r^4}\rho \mathrm{d}\rho \mathrm{d}\varphi .$$

(21)

Using $r=h\sqrt{1+{\epsilon }^2}$, we obtain

$$FSIR\left(t_r,k\right)=\frac{{\lambda }^2{\sigma }^0{G}_0^2}{{\left(4\pi \right)}^3{L}_p{h}^4}{\int }_{{\mathbb{R}}^{+}\times D_{t_r,k}}\frac{\delta \left(t_r-2R_0/c\right)exp\left\{-\left(\frac{4}{\gamma }{sin}^2\left[\theta \left(\rho ,\varphi \right)\right]\right)\right\}}{{\left(1+{\epsilon }^2\right)}^2}\rho \mathrm{d}\rho \mathrm{d}\varphi$$

(22)

with

$$cos\left[\theta \left(\rho ,\varphi \right)\right]=\frac{cos\left(\xi \right)+\frac{\rho }{h}sin\left(\xi \right)cos\left(\tilde{\varphi }-\varphi \right)}{\sqrt{1+\frac{{\rho }^2}{h^2}}}.$$

(23)

3.2. Semi-Analytical Expression of $FSIR\left(t_r,k\right)$

Using the change of variable $\begin{array}{c}x=\left(2h/c\right)\sqrt{1+{\epsilon }^2}\end{array}$ and the property of the Dirac distribution $h\left(t\right)={\int }_{c_1}^{+\infty }\delta \left(x-t\right)h\left(x\right)\mathrm{d}x$ for $t\ge c_1$, integrating with respect to x leads to

$$\begin{array}{c}FSIR\left(t_r,k\right)=\frac{P_u}{2\pi }{\left(1+\frac{ct}{2h}\right)}^{-3}J\left(t\right)\times \left\{{\int }_{D_{t_r,k}}F\left[\tilde{\varphi }-\varphi ,\epsilon ,\xi ,\mu \right]\mathrm{d}\varphi \right\}\end{array}$$

(24)

with

$$\begin{array}{c}F\left[\tilde{\varphi }-\varphi ,\epsilon ,\xi ,\mu \right]=exp\left\{-4\left[1-{cos}^2\left(\xi \right)/\left(1+{\epsilon }^2\right)\right]/\gamma \right.\left.+b\left(1-{sin}^2\left(\tilde{\varphi }-\varphi \right)\right)+acos\left(\tilde{\varphi }-\varphi \right)\right\}\hfill \end{array},$$

$$\epsilon =\rho /h,$$

$$P_u=\frac{{\lambda }^2{G}_0^{2}c{\sigma }^0}{4{\left(4\pi \right)}^2{L}_p{h}^3},$$

$$a=\frac{4\epsilon sin\left(2\xi \right)}{\gamma \left(1+{\epsilon }^2\right)},$$

$$b=\frac{4{\epsilon }^2{sin}^2\left(\xi \right)}{\gamma \left(1+{\epsilon }^2\right)},$$

(25)

where $J\left(·\right)$ is the step function and $D_{t,k}$ and ${\phi }_{t,k}$ are both provided in (19). Performing a variable conversion $w=\tilde{\varphi }-\varphi$ in (24) results in

$$\begin{array}{c}\mathrm{FSIR}\left(t_r,k\right)=\frac{P_u}{2\pi }{\left(1+\frac{ct}{2h}\right)}^{-3}J\left(t\right)\times exp\left\{-4\left[1-{cos}^2\left(\xi \right)/\left(1+{\epsilon }^2\right)\right]/\gamma +b/2\right\}\\ \times \left[f\left(\tilde{\varphi }-{\varphi }_{t_r,k+1},\tilde{\varphi }-{\varphi }_{t_r,k}\right)+f\left(\tilde{\varphi }-{{\varphi }^{\prime }}_{t_r,k+1},\tilde{\varphi }-{{\varphi }^{\prime }}_{t_r,k}\right)\right]\end{array}$$

(26)

with

$$f\left(w_1,w_2\right)={\int }_{u_1}^{u_2}exp\left[acos\left(w\right)+\frac{b}{2}cos\left(2w\right)\right]\mathrm{d}w,$$

(27)

where $f\left(·\right)$ can be expanded using Bessel functions as follows:

$$exp\left[acos\left(w\right)\right]=I_0\left(a\right)+2\sum _{k=1}^{\infty }{I}_k\left(a\right)cos\left(kw\right)$$

$$exp\left[\frac{b}{2}cos\left(2w\right)\right]=I_0\left(\frac{b}{2}\right)+2\sum _{k=1}^{\infty }{I}_k\left(\frac{b}{2}\right)cos\left(2kw\right),$$

(28)

where $I_k$ denotes the first-kind modified Bessel function of order k. Note that if any of these terms are approximated, the higher-order expansions of $exp\left[\frac{b}{2}cos\left(2w\right)\right]$ can lead to an unacceptable increase in model error when the mis-pointing angle exceeds 1 degree. In this paper, the complete analytical form of the Bessel function approximation is obtained without neglecting any terms. Subsequently, upon substituting (28) into (27) and simplifying, we have

$$\begin{array}{cc}\hfill f\left(w_1,w_2\right)& =\sum _{k=0}^{\infty }\sum _{j=0}^{\infty }{C}_k{C}_j{I}_k\left(a\right)I_j\left(\frac{b}{2}\right)/2\times \{w_2\left[\mathrm{sinc}\left(\left(k-2j\right)w_2\right)-\mathrm{sinc}\left(\left(k+2j\right)w_2\right)\right]\hfill \\ \hfill & -w_1\left[\mathrm{sinc}\left(\left(k-2j\right)w_1\right)-\mathrm{sinc}\left(\left(k+2j\right)w_1\right)\right]\}\hfill \end{array},$$

(29)

where $C_k,C_j=\left\{\begin{array}{c}1,k=0\hfill \\ 2,k\ge 1\hfill \end{array}\right.$. The infinite sum provided in (29) can be truncated to reduce computational complexity by maintaining a finite number m of elements based on the accuracy requirement, as follows:

$$\begin{array}{cc}\hfill f_1\left(w_1,w_2\right)& =\sum _{k=0}^m\sum _{j=0}^m{C}_k{C}_j{I}_k\left(a\right)I_j\left(\frac{b}{2}\right)/2\times \{w_2\left[\mathrm{sinc}\left(\left(k-2j\right)w_2\right)-\mathrm{sinc}\left(\left(k+2j\right)w_2\right)\right]\hfill \\ \hfill & -w_1\left[\mathrm{sinc}\left(\left(k-2j\right)w_1\right)-\mathrm{sinc}\left(\left(k+2j\right)w_1\right)\right]\}.\hfill \end{array}$$

(30)

Finally, the analytical expression of FSIR is provided by

$$\begin{array}{c}\mathrm{FSIR}\left(t_r,k\right)=\frac{P_u}{2\pi }{\left(1+\frac{ct}{2h}\right)}^{-3}J\left(t\right)\times exp\left\{-4\left[1-{cos}^2\left(\xi \right)/\left(1+{\epsilon }^2\right)\right]/\gamma +b/2\right\}\\ \times \left[f_1\left(\tilde{\varphi }-{\varphi }_{t_r,k+1},\tilde{\varphi }-{\varphi }_{t_r,k}\right)+f_1\left(\tilde{\varphi }-{{\varphi }^{\prime }}_{t_r,k+1},\tilde{\varphi }-{{\varphi }^{\prime }}_{t_r,k}\right)\right].\end{array}$$

(31)

3.3. Multilooking

The reflection power is obtained through the triple convolution in (13), where the expressions of $FSIR(t_r,k)$ are provided in (31). To obtain multilook waveforms, delay compensation can be applied to the DDM [19]. The range of each beam n needs to be offset for a particular range bin in order to correspond to the same across-track position for all beams [24]. After compensation, the range indicated for each scatterer over its entire illumination history is equal to its minimum range. The delay of each beam ${\delta }_r\left(n\right)$ is obtained by the difference between the modulus of the position vector $r_n=\sqrt{h^2+y_n^2}$ and the minimum height between the surface and the platform, as follows:

$${\delta }_r\left(n\right)=r_n\left(1/sin{\theta }_n-1\right),$$

$$cos{\theta }_n=\frac{\lambda f_n}{2v_{c}cos\mu }+\frac{h_{0}tan\mu }{r_n},$$

(32)

where the minimum height $h_0$ can be obtained by tracking algorithm. Then, the final multilook waveform can be obtained through the reflected power after delay compensation, as follows:

$$s\left(t\right)=\sum _{n=1}^{N}P\left(t-2{\delta }_r\left(n\right)/c,f_n\right).$$

(33)

The traditional method for delay compensation in SAR altimetry, as described in [19], typically employs nearest neighbor interpolation in the range domain. In this paper, a novel delay compensation method for SAR altimeters with vertical movement is proposed. The proposed method utilizes sinc interpolation, as depicted in Algorithm 1, to optimally handle non-integer delays and maintain the frequency characteristics of the signal. The resulting multilook waveform is tracked according to the maximum value to find the range gate. We compare the tracked waveform to the proposed theoretical echo model in order to determine the parameter value exhibiting the highest similarity through the corresponding re-tracking algorithm detailed in Section 4.

Algorithm 1 Delay Compensation for Airborne SAR Altimetry

Input: range r, Doppler frequency $f_d$, The reflected power $P_I\left(r,f_d\right)$, minimum height $h_0$  1: resolution of the along-track $dy=vcos\mu N\mathrm{PRI}$;  2: $\mathrm{ylimit}=\lambda f_{d,max}{r}_{max}/2vcos\mu$;  3: $y_n=-\mathrm{ylimit}+htan\mu :dy:\mathrm{ylimit}+htan\mu$;  4: for each i in r do  5:   $f_n$ according to (9);  6:   $P\left(r,f_n\right)=Interpolation\left(f_d,P\left(r,f_d\right),f_n\right)$;  7: end for  8: $r_n=\sqrt{r·r+{\left(y_n·y_n\right)}^{\mathrm{T}}}$  9: for each n in $y_n$ do 10:   $P\left(r_n,f_n\right)=Interpolation\left(r,P\left(r,f_n\right),r_n\right)$; 11: end for Output: The reflected power after delay compensation $P\left(r_n,f_n\right)$

3.4. Analysis of the Semianalytical Echo Model

The model aggregates the backscattered energies within the illuminated region to obtain a 2D FISR. The temporal dimension is determined by the propagating circle, while the frequency dimension is defined by the rectangular Doppler beam [22]. Each circle with a radius $\rho \left(t\right)=\sqrt{hct}$ and each approximate rectangular Doppler beam divides the illuminated area into distinct delay/Doppler bins. Note that the flight path angle $\mu$ introduces a Doppler frequency shift when the airborne platform flies obliquely (see Figure 3). As $\mu$ increases, all Doppler beams exhibit a negative Doppler shift $f_{d,shift}=-2v_{c}hsin\mu /\lambda r_{xy}$. Additionally, the power within each delay/Doppler bin changes, deviating from a symmetric distribution. Due to disparities in offset range among different delay–Doppler bins after multilooking, the leading edge of the multilook waveforms is correspondingly affected, as shown in Figure 4.

However, changes in the velocity direction of the airborne platform generally result in antenna mis-pointing angles. Although the mapping relationship is unchanged when considering $\xi$ and $\tilde{\varphi }$, the backscattered energy within the delay–Doppler bins is impacted. In the case where the antenna is vertically fixed to the aircraft, the two mis-pointing angles $\xi$ and $\tilde{\varphi }$ are dependent on the pitch angle ${\psi }_{\mathrm{al}}$ and the roll angle ${\psi }_{\mathrm{ac}}$ as follows:

$$\begin{array}{c}{\psi }_{\mathrm{al}}=arctan\left(tan\left(\xi \right)sin\left(\tilde{\varphi }\right)\right)\hfill \\ {\psi }_{\mathrm{ac}}=arctan\left(tan\left(\xi \right)cos\left(\tilde{\varphi }\right)\right).\hfill \end{array}$$

(34)

Figure 5 illustrates the variations in antenna gain across the various mis-pointing angles. It can be observed that the maximum antenna gain is achieved at the nadir point when the antenna mis-pointing angle is zero. As expected, the maximum antenna gain shifts along either the x or y axis as the antenna mis-pointing angle changes (see Figure 5c,d). The corresponding multilook echoes are significantly influenced by the variations in antenna gain. Next, we analyze the impact of antenna gain variation on multilook echoes.

We can evaluate the effect of antenna mis-pointing angle on both the multilook and normalized multilook echoes. As shown in Figure 6, increasing the mis-pointing angle along the track leads to a decrease in the amplitude of the multiview echo. However, this has a minimal effect on the normalized multiview echo, indicating that the mis-pointing angle along the track primarily affects the echo’s amplitude rather than its shape.

In addition, Figure 7 demonstrates that the mis-pointing angle across the track has a more pronounced impact on the waveform’s amplitude compared to the along-track mis-pointing angle, particularly affecting the trailing edge of the waveform. Note that the mis-pointing echoes obtained by the proposed model are influenced by six parameters $\left(P_u,\mathrm{SWH},\tau ,\mu ,{\psi }_{ac},{\psi }_{al}\right)$. The waveform’s shape is influenced by both the antenna mis-pointing and flight path angles, which should be considered in parameter estimation. Therefore, our focus is on analyzing the effect of the parameter vector $\left(\mathrm{SWH},\tau ,\mu ,{\psi }_{ac},{\psi }_{al}\right)$ on the normalized echo.

4. Re-Tracking Algorithm Based on Echo Model Fitting

Re-tracking processing involves parameterizing the echo function curve based on the rough tracking result. It uses a certain algorithm to fit the echo model to the actual waveform to achieve the most consistent state, thereby retrieving the parameters contained in the echo. It can be solved using the hill-climbing algorithm [32,33] or the least squares (LS) method [34]. However, in certain scenarios the paths defined by the hill-climbing algorithm do not monotonically reduce the distance to the local optimum [35]. Thus, the nonlinear least squares (LS) algorithm is widely used in altimetry parameter estimation. The LS estimator under consideration is defined as

$${\widehat{\theta }}_{LS}=\underset{\theta }{argmin}\frac{1}{2}\sum _{k=1}^K{g}_k^2\left(\theta \right),$$

(35)

where $g_k\left(\theta \right)=y_k-s_k\left(\theta \right)$ is the residual vector, $s\left(\theta \right)$ represents the multilook echoes, and $y={\left(y_1,\dots ,y_K\right)}^T$ is the echo destroyed by speckle noise. Because $g_k^2\left(\theta \right)$ is a complicated nonlinear function of SWH and $\tau$, the optimization problem does not admit a closed-form expression, and the LS estimator (35) is solved using a numerical optimization method.

The Levenberg–Marquardt (LM) algorithm is an optimization algorithm used for nonlinear least squares fitting and optimization [12,22,36]. It exhibits greater robustness and effectiveness than traditional LS methods. In this paper, a five-parameter iterative LS algorithm based on the LM algorithm is used to estimate the altimetric parameters, as shown in Figure 8.

The parameter update of the iterative LM algorithm is defined by ${\theta }^{\left(i+1\right)}={\theta }^{\left(i\right)}+e^{\left(i\right)}$, where ${\theta }^i$ is the estimate of $\theta$ at the i-th iteration. The choice of $e^i$ is based on Taylor expansion (at the first order) of g in the neighborhood of ${\theta }^i$

$$g\left({\theta }^{\left(i\right)}+e^{\left(i\right)}\right)=g\left({\theta }^{\left(i\right)}\right)+J\left({\theta }^{\left(i\right)}\right)e^{\left(i\right)}$$

(36)

where $J=\left[\partial g\left({\theta }^{\left(i\right)}\right)/\partial {\theta }_1,\dots ,\partial g\left({\theta }^{\left(i\right)}\right)/\partial {\theta }_J\right]$ is a partial differential matrix representing the gradient of g. After replacing (36) in (35), the following result is obtained:

$$G\left(\theta +e\right)\simeq L\left(e\right)=\frac{1}{2}l{\left(e\right)}^{T}l\left(e\right).$$

(37)

The descent direction e is then obtained by minimizing $L\left(e\right)$. By setting the derivative $L^{\prime }\left(e\right)=J{\left(\theta \right)}^{T}g+J{\left(\theta \right)}^{T}J\left(\theta \right)e$ to zero and adding a regularization parameter $\mu$, leading to

$$\left[J{\left(\theta \right)}^{T}J\left(\theta \right)e+\mu I_J\right]e=-J{\left(\theta \right)}^{T}g$$

(38)

where $I_J$ is the identity matrix, this algorithm employs a gradient descent method to update the parameters $\theta =\left(\mathrm{SWH},\tau ,\mu ,{\psi }_{ac},{\psi }_{al}\right)$, as follows:

$${\theta }^{\left(i+1\right)}={\theta }^{\left(i\right)}-{\left[J^{T}J+{\mu }_i{I}_J\right]}^{-1}{J}^{T}g\left({\theta }^{\left(i\right)}\right)$$

(39)

where ${\mu }_i$ controls the convergence speed of the algorithm. Note that the derivatives appearing in $\frac{\partial g\left(\theta \right)}{\partial {\theta }_J}$ can be computed numerically using the finite difference, as follows:

$$\frac{\partial g\left(\theta \right)}{\partial {\theta }_j}\simeq -\frac{s\left({\theta }_j+\mathrm{\Delta }{\theta }_j\right)-s\left({\theta }_j\right)}{\mathrm{\Delta }{\theta }_j}.$$

(40)

It should be noted that the initial value of the iterated parameters in the re-tracking algorithm is ${\theta }^0=\left({\tau }^0,\mathrm{SW}{\mathrm{H}}^0,{\mu }^0,{\psi }_{ac}^0,{\psi }_{al}^0\right)$, where ${\tau }^0$ is obtained by rough tracking. Furthermore, the real mis-pointing angles (i.e., the roll and pitch angles) are used as inputs, deviating from their actual values by approximately $1^{\circ }$. These mis-pointing angles are recorded by the IMU during airborne experiments, and the deviation in the input mis-pointing angles depends on the precision of the IMU devices. This enhances the fitness between the theoretical echo model and the actual airborne waveform, thereby reducing the possibility of tracking failure caused by the low similarity between the theoretical model and actual waveform [12,37].

5. Results

In this section, we present the simulation results obtained with the proposed model. First, the approximations used to obtain the analytical FSIR are justified in Section 5.1. Then, the accuracy of the proposed model is evaluated through comparison with the fully numerical model in Section 5.2. Finally, in Section 5.3 we verify the effectiveness of the proposed estimation algorithm by comparing it with an algorithm that does not consider mis-pointing angles.

5.1. Analysis of the FSIR Approximations

The purpose of this section is to quantitatively analyze the error of the proposed model caused by the m-order approximation used in FSIR. The multilook echoes obtained from the proposed model are analyzed to examine the effect of different antenna mis-pointing angles and varying numbers of terms m. The validation involves comparing the model (13) with the numerical calculation (21). The approximation error of the proposed model is evaluated using the normalized quadratic error (NQE), which is defined as follows:

$$\mathrm{NQE}\left(s,r\right)=\sqrt{\frac{{\sum }_{k=1}^K{\left(s_k-r_k\right)}^2}{{\sum }_{k=1}^K{r_k}^2}}.$$

(41)

The approximation of the FSIR is affected by the mis-pointing angle of the antenna and the number of terms m. Considering the use of normalized waveforms in the re-tracking algorithm, it is necessary to explain the relationship between the maximum error of the echoes and the number of terms m. NQE is used to describe the error of the maximum value of the echo $\mathrm{NQ}{\mathrm{E}}_{max}=\sqrt{{\left(s_{max}-r_{max}\right)}^2/r_{max}^2}$.

Figure 9 shows the results of the NQE of the semi-analytical model as it changes with ${\psi }_{ac}$. It can be seen that NQE increases as the mis-pointing angle ${\psi }_{ac}$ becomes larger. When considering ${\psi }_{ac}\le {18}^{\circ }$ (${\psi }_{ac}<{\theta }_{3\mathrm{dB}}/2$), $m=4$ is adequate to achieve the required level of error.

For the case in which $m=4$, the global NQE matrix of the antenna mis-pointing angle is shown in Figure 10. Note that the NQE of a noisy echo may increase in real application, which requires m to take a larger value in order to reduce the error.

5.2. Comparison with Numerical Waveforms

The proposed semi-analytical model reduces the high dimensionality of the proposed numerical model by introducing approximations into the provided impulse response under the hypothesis of larger variation of the mis-pointing angles. We performed echo simulation in typical airborne scenarios; the simulation parameters are shown in

Table 1. Simulation parameters.

Parameter	Value
Frequency	2.95 GHz
Wavelength ($\lambda$)	0.1 cm
Bandwidth (B)	100 MHz
Mean flight altitude (h)	2000 m
Pulse repetition frequency (PRF)	5000 Hz
3dB Antenna beam width (${\theta }_{3\mathrm{dB}}$)	${40}^{\circ }$
Velocity ($v_s$)	100 m/s
Pulse per burst	100 pulses

. The multilook waveform obtained by the proposed model is considered to be a DDA echo, while the multilook waveform corresponding to the numerical integral form of (14) is considered to be a numerical echo. By comparing the DDA echo and numerical echo, we aim to assess how the accuracy of the semi-analytical model decreases as the mis-pointing angles increase. In addition, the semi-analytical model with all mis-pointing angles at zero is defined as G-DDA for comparison. Through comparison with the generalized model, we can evaluate the accuracy of the semianalytical waveforms with respect to the fully numerical waveforms for different values of the mis-pointing angles and flight path angle.

In Figure 11, the multilook waveforms computed from the different model implementations are shown when $\left(\mu ,{\psi }_{ac},{\psi }_{al}\right)=\left(0^{\mathrm{o}},0^{\mathrm{o}},0^{\mathrm{o}}\right)$ and $\left(\mu ,{\psi }_{ac},{\psi }_{al}\right)=\left(18^{\mathrm{o}},0^{\mathrm{o}},18^{\mathrm{o}}\right)$. The waveforms have been normalized with respect to the peak power. From Figure 11 (left side), it can be seen that both the DDA echo and G-DDA echo are close to the numerical echo in the absence of mis-pointing. From Figure 11 (right side), it can be noticed that the shape of the G-DDA echo is far from the expected one when the pitch angle ${\psi }_{al}$ increases, while this does not occur in case of the DDA echo. In Figure 12, the multilook waveforms are computed when the angles are set as $\left(\mu ,{\psi }_{ac},{\psi }_{al}\right)=\left(10^{\mathrm{o}},10^{\mathrm{o}},10^{\mathrm{o}}\right)$ and $\left(\mu ,{\psi }_{ac},{\psi }_{al}\right)=\left(0^{\mathrm{o}},10^{\mathrm{o}},0^{\mathrm{o}}\right)$. From Figure 12, it can be noticed that while the shape of the trailing edge in numerical echo has a different decay with respect to the DDA echo and G-DDA echo when the roll increases, the DDA echo is closer to the numerical echo.

5.3. Analysis of Re-Tracking Precision

This section evaluates the performance of the proposed re-tracking algorithms for the numerical waveforms affected by speckle noise. To simulate the real echo, speckle noise is added to the multilook waveform (33) following the approach described in [12,14,22]. Consequently, the noisy echo $y\left(t\right)$ can be obtained by

$$y\left(t\right)=\sum _{k=1}^{N}P\left(t_r-2{\delta }_r\left(f_d\right)/c,k\right)q\left(t_r-2{\delta }_r\left(f_d\right)/c,k\right),$$

(42)

where $q\left(t_r-2{\delta }_r\left(f_d\right)/c,k\right)$ is gamma-distributed, multiplicatively independent, and identically distributed speckle noise expressed as $\Gamma \left(L,1/L\right)$.

The normalized echo used for parameter estimation depends on the unknown parameter vector ${\theta }_5={\left(\mathrm{SWH},\tau ,\mu ,{\psi }_{ac},{\psi }_{al}\right)}^T$. The effectiveness of the strategy can be evaluated through the simulation of multilook echoes. The root mean square error (RMSE) of the altimetric parameters is adopted to evaluate the algorithm’s performance. All results presented in this section have been averaged using NMC = 500 Monte Carlo runs. The RMSE of the i-th altimetric parameter ${\theta }_i$ is defined as

$$\mathrm{RMSE}\left({\theta }_i\right)=\sqrt{\frac{1}{N}\sum _{j=1}^N{\left[{\theta }_i-{\widehat{\theta }}_i\left(j\right)\right]}^2},i=1,\dots ,J,$$

(43)

where ${\theta }_i$ is the true parameter, ${\widehat{\theta }}_i\left(j\right)$ is the estimated parameter for the j-th waveform, and $N=500$ is the Number of Monte Carlo simulations.

To verify the effectiveness of the proposed model, a generalized delay–Doppler altimetric model is used for comparison. We replace ${\psi }_{ac}$ and ${\psi }_{al}$ with zeros and estimate the two remaining parameters ${\left(\mathrm{SWH},\tau \right)}^T$ using the Levenberg–Marquardt algorithm. The resulting estimation strategy is denoted as G-PRA (indicating the generalized altimetric model with two unknown parameters used in the parameterized re-tracking algorithm). The estimations provided by G-PRA need to be compared with those obtained with the proposed algorithm, which is denoted as PRA.

The first group of experiments corresponds to level flight in the absence of mis-pointing while considering different SWH; the other altimetric parameters are set as ${\left(P_u,\tau ,\mu ,{\psi }_{ac},{\psi }_{al}\right)}^T={\left(1,30\mathrm{gates},0^{\mathrm{o}},0^{\mathrm{o}},0^{\mathrm{o}}\right)}^T$. The parameter RMSEs obtained by the algorithm are shown in Figure 13, and the figure shows similar performance for the G-PRA and PRA algorithms. Thus, the results demonstrate that there is not much performance reduction in estimating $\tau$ and SWH in the absence of mis-pointing. Additionally, the $\mathrm{RMSE}\left(\mathrm{SWH}\right)$ is a decreasing function of SWH.

The second set of experiments aimed to investigate the performance of parameter estimation for oblique flight with varying mis-pointing angles ${\psi }_{al}$ and ${\psi }_{ac}$. The remaining altimetric parameters are set as ${\left(\mathrm{SWH},P_u,\tau ,\mu \right)}^T={\left(2\mathrm{m},1,30\mathrm{gates},6^{\mathrm{o}}\right)}^T$. Figure 14 demonstrates that the proposed PRA algorithm leads to a smaller $\mathrm{RMSE}\left(\tau \right)$ compared to G-PRA in the presence of ${\psi }_{ac}$. As ${\psi }_{ac}$ increases, the difference in $\mathrm{RMSE}\left(\tau \right)$ and $\mathrm{RMSE}\left(\mathrm{SWH}\right)$ obtained by the two algorithms becomes larger. However, the improvement is less in the presence of ${\psi }_{al}$, as shown in Figure 15. This property can be explained by the fact that the cross-track mis-pointing angle has a greater influence on the waveform shape compared to ${\psi }_{al}$. In addition, It can be observed that as the mis-pointing angle increases, the RMSE of $\tau$ obtained by the PRA algorithm initially remains at a low level and then gradually increases. This may be attributed to the non-negligible error in model approximation. The mis-pointing angle in airborne scenarios is typically not greater than 10 degrees, indicating that the estimation accuracy of this method falls within an acceptable range.

6. Conclusions

In order to perform complex navigation tasks, high-precision altimetry on highly mobility platforms is attracting growing attention. However, the waveform re-tracking methods proposed in previous works cannot be applied to highly mobility platforms, which may have a large antenna mis-pointing angle and significant vertical movement. To improve the estimation accuracy of waveform re-tracking, in this paper we propose a novel semi-analytical waveform model and signal processing method for SAR altimeters with vertical movement and large antenna mis-pointing angles. Our main conclusions are as follows:

(1)

An analytical expression that considers mis-pointing angles, circular antenna patterns, and vertical velocity is introduced for the echo model. The proposed model shows that the across-track mis-pointing angle affects the shape and the amplitude of the altimetric echo, whereas the along-track mis-pointing angle mainly affects the amplitude of the echo. When we consider $\xi \le {20}^{\circ }$ ($\xi <{\theta }_{3\mathrm{dB}}/2$), the expansion order $m=4$ is adequate to achieve the required level of error, which means that the global NQE is less than ${10}^{-10}$.

(2)

A novel delay compensation method based on sinc interpolation is proposed to obtain the multilook echo in order to optimally handle non-integer delays and maintain the signal frequency characteristics. Compared with the generalized model, the proposed semi-analytic model is closer to the fully numerical waveform at different values of mis-pointing angle ${\psi }_{ac}$ and ${\psi }_{al}$, especially when the roll angle ${\psi }_{ac}$ is large.

(3)

A five-parameter estimation strategy using the least squares procedure is proposed for SAR altimeters with vertical movement. The performance of the proposed model is evaluated using simulated data for typical airborne scenarios, analyzing the influence of different mis-pointing angles and flight path angle on parameter estimation. When the mis-pointing angles are within 10 degrees, the RMSE of $\tau$ obtained by the re-tracking method fitted by the proposed model is less than 0.2 m, and increases more slowly compared to the one fitted by the generalized model.

In February 2024, a more in-depth airborne radar altimeter experiment will be conducted in southern China. The model and re-tracking methods proposed in this paper will be validated against real data collected during this experiment. In addition, future research will focus on exploring the optimal performance of the estimation strategies under various antenna mis-pointing angles; it is worth considering a further modification to the strategy to involve fewer parameters in iterative updates.