Accelerated algorithm for back-projection fully focused SAR altimetry

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 5209816 Accelerated Algorithm for Back-Projection Fully Focused SAR Altimetry Sergi Hernández-Burgos , Ferran Gibert , Member, IEEE, Antoni Broquetas , Life Member, IEEE, Pol Villalvilla, Alejandro Egido, Gorka Moyano, Adrián Flores de la Cruz , and Marco Fornari Abstract— Fully focused (FF) synthetic aperture radar (SAR) back-projection algorithms for radar altimetry have gained attention in the altimetry community due to their improved azimuth resolution compared to traditional methods, such as the delay/Doppler algorithm. However, the significant computational cost of back-projection has limited its use in operational applications. This article presents an accelerated back-projection algorithm that reduces the runtime of the classic back-projection by a factor of 28 on a CPU-based architecture. When implemented on a GPU-based architecture, the accelerated back-projection achieves a speedup ranging from 8 to 13 times compared to the classic back-projection on the same GPU. In comparison to the classic back-projection running on a CPU, the combination of the accelerated back-projection and GPU processing results in an overall improvement of up to 1570 times. The accelerated algorithm maintains the accuracy of classic time-domain algorithms and allows for the processing of data in the order of the sensing time under a GPU-based architecture. The algorithm has been extensively validated using transponder and open ocean data. These results indicate that operational back-projection algorithms are feasible for current radar altimetry missions, such as Sentinel-6, and upcoming missions like CRISTAL and Sentinel-3 Next Generation. Index Terms— Accelerated algorithm, accelerated backprojection, back-projection, delay/Doppler, fully focused (FF), fully focused synthetic aperture radar (FF-SAR), GPU acceleration, high-resolution radar altimetry, omega-K, open ocean, radar altimetry, real time altimetry, synthetic aperture radar (SAR), transponder. I. INTRODUCTION AT THE beginning of space-based altimetry, radar altimeters missions such as the ERS-1 (1991) [1], Jason-1 (2001) [2], and ENVISAT (2002) [3] used conventional Received 18 November 2024; revised 21 February 2025; accepted 26 March 2025. Date of publication 31 March 2025; date of current version 29 April 2025. This work was supported in part by the Industrial Doctorates Plan from the Secretary of Universities and Research of the Department of Business and Knowledge of the Government of Catalonia through Project SA.45217 with Expedient Number 094 under Grant Doctorats Industrials 2020 and in part by European Space Agency (ESA) through the Poseidon-4 Ground Processor Prototype Project under ESA Contract 4000112536/14/NL/BJ. (Corresponding author: Sergi Hernández-Burgos.) Sergi Hernández-Burgos is with isardSAT, S.L., 08005 Barcelona, Catalonia, and also with the Signal Theory and Communications Department, Universitat Politecnica de Catalunya, 08034 Barcelona, Catalonia (e-mail: [email protected]). Ferran Gibert, Pol Villalvilla, Gorka Moyano, and Adrián Flores de la Cruz are with isardSAT, S.L., 08005 Barcelona, Catalonia. Antoni Broquetas is with the Signal Theory and Communications Department, Universitat Politecnica de Catalunya, 08034 Barcelona, Catalonia. Alejandro Egido is with European Space Agency, European Space Research and Technology Centre, 2201 AZ Noordwijk, The Netherlands. Marco Fornari is with the Starion Group for European Space Agency (ESA)/ESTEC, 2201 AZ Noordwijk, The Netherlands. Digital Object Identifier 10.1109/TGRS.2025.3556544 pulse transmission at low pulse repetition frequency (PRF), enabling the observation of Earth surface with an azimuth resolution limited by the antenna footprint, typically in the order of kilometers. Later on, radar altimetry missions such as CryoSat-2 (2010) [4] and Sentinel-3 (2016) [5] were based on the transmission of short bursts of pulses at very high PRF, approximately 19 kHz in a closed-burst mode of 64 pulses, assuring coherence among them and increasing the number of independent looks obtained from a single scatterer. The high PRF together with the pulse coherence allowed to use of the delay/Doppler algorithm [6], improving the azimuth resolution up to 300 m. Newer missions, such as Sentinel-6 [7], employ nearly continuous pulse transmission around 9–10 kHz by interleaving the transmission and reception of pulses, maximizing the measurement precision while ensuring a connection between data from prior altimeter missions, such as Jason series [8]. Currently, high-resolution operational processors are based on the delay/Doppler algorithm. Egido and Smith [9] introduced the fully focused synthetic aperture radar (FF-SAR) back-projection algorithm, which further improved the azimuth resolution to the theoretical maximum, in the order of submeter. This enhancement in azimuth resolution has allowed for a reduction in noise in open ocean retracked parameters such as sea surface height (SSH) and significant wave height (SWH) [10]. FF-SAR algorithms have also facilitated more precise identification of leads and icebergs [11], the study of swells from altimetry data [12], and the estimation of inland water parameters such as water elevation [13] and water extent [14]. However, the significant high-computational effort required by the classic FF-SAR back-projection algorithm to apply all the necessary phase corrections for each along-track point on the surface has so far impeded its implementation as an operational processor. Consequently, alternative and more efficient FF-SAR algorithms have been developed. In 2019, a numerical 2-D frequency domain algorithm was presented for the CryoSat-2 mission [15]. In 2024, a closed-form 2-D frequency domain algorithm was introduced for the Sentinel6 mission [16]. Nonetheless, the frequency-domain algorithms require specific local conditions for the approximations to be effective, such as constant satellite velocity and constant PRF. Moreover, these approximations can reduce the quality performance of the results compared to time-domain algorithms such as the FF SAR back-projection. These limitations motivate the need to investigate ways to accelerate the back-projection algorithm. © 2025 The Authors. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 5209816 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 In this article, we present a variant of the back-projection algorithm called accelerated back-projection. The accelerated back-projection technique reduces the computational cost of the classic back-projection algorithm through certain approximations, without significantly compromising the performance over distributed targets. For that, we assume that the radar geometry does not change significantly with respect to short along-track distances on the surface, and we can approximate this small change using a linear approximation, reducing considerably the number of operations to focus. This linear approximation is well known in the SAR radar altimeter community, as it has been commonly employed in high-resolution algorithms such as the delay/Doppler method [6] to obtain a closed-form expression for azimuth beamforming after a Fourier Transform. This same approximation is briefly mentioned in [13, Sec. 3.2.1.6]in the context of the FF-SAR back-projection algorithm, where it is applied and noted without further analysis. Moreover, as explained in [17] for SAR imaging, an additional speeding-up approach can be achieved by splitting and averaging the synthetic aperture in small subapertures, reducing considerably the number of operations needed to focus the surface. Since maintaining the accuracy of the classic backprojection is essential, a validation analysis is conducted with real Sentinel-6 data to understand the limitations of the approximations and to ensure that the quality performance of time-domain algorithms is not compromised: First, point target scenarios such as transponders [18] and corner reflectors [19] allow us to validate the real point target response (PTR) and compare it to the theoretical PTR. Also, they allow us to identify the real limits of the accelerated back-projection. After that, tests over the open ocean have been conducted to evaluate the accuracy and precision of the accelerated algorithm when retrieving SSH, SWH, and σ0. This evaluation aims to ensure that we improve considerably the runtime performance without sacrificing the accuracy with respect to the classic backprojection algorithm. Moreover, we have conducted runtime performance tests using CPU-based processors and GPU-based processors. Our results show that we can increase the runtime performance of the classic back-projection by ×28 using CPU, and by ×1570 using GPU architecture, all while assuring good quality performance. The speeding factors achieved without compromising the geophysical retrievals pave the way for the implementation of FF-SAR products in future operational processors. II. GEOMETRIC FOUNDATIONS OF THE ALGORITHM The main acceleration factor of the proposed algorithm is based on the strong linear phase dependence that is observed along the illumination time over a bright point target when focusing not at the specific location of the point target but at a certain nearby surface in the along-track direction, in the order of meters. Fig. 1shows the azimuth phase of all the pulses within the illumination time over the transponder in Gavdos [20] for a Sentinel-6 pass on March 15, 2024, when focusing at different distances along-track from the real location of the point target. The azimuth phase is determined by obtaining the phase of each received pulse within the illumination time, specifically Fig. 1. Sentinel-6 L1A azimuth phase (rad) at the range of maximum power during 3.4 s of synthetic aperture for different focusing distances with respect to the transponder of Gavdos. (Top) Along-track corrected phase. (Bottom) Detrended phase. The slope of the linear phase is a function of the along-track distance between the point target location and the focusing point. Fig. 2. Along-track phase slope difference between the Gavdos transponder and the theoretical model. The transponder phase slope is determined through linear fitting. at the range sample with the maximum power, after applying all the classic back-projection corrections and prior to the coherent integration. As expected, when focusing all the pulses at the real point target location, the phase is constant, which allows for a full constructive summation of the signal when coherently combining all the pulses within the illumination time. However, at any point with a certain along-track distance nearby to the point target, the phase presents a constant slope. The slope value is proportional to the along-track distance from the point target. Indeed, given specific conditions, such slope can be determined in a closed-form expression. In such cases, it is possible to focus the surface points at nearby along-track distances in a faster manner, just applying a secondary phase correction proportional to the phase slope (linear term). Fig. 2shows the difference between the phase slope obtained from Fig. 1by means of a linear fitting and the closed-form expression. The maximum error is below 15 mrad, less than 0.02%. HERNÁNDEZ-BURGOS et al.: ACCELERATED ALGORITHM FOR BACK-PROJECTION FF SAR ALTIMETRY 5209816 Fig. 3. Evolution of the along-track phase from waveforms obtained around a scenario with a single point target. If we apply the back-projection corrections to a surface point before the point target, namely at point A, the remaining phase slope obtained will be positive. If the focusing point is after the point target (namely at point B), the remaining phase slope obtained will be negative. As expected, the phase is completely flat only when focusing on the surface point where the point target is located. Fig. 3illustrates a qualitative explanation in the case of an isolated point target. When applying the classic back-projection algorithm to focus on point A, the resulting azimuth phase exhibits a linear behavior, proportional to the along-track distance between point A and the target. Then, after focusing on point A, the corrected waveforms can be reused by applying an additional secondary (and linear) phase correction to retrieve the azimuth phase corresponding to the point target. Similarly, the secondary phase correction can also be applied to obtain the phase corresponding to point B. In this process, most of the back-projection operations are performed only once for point A, while focusing on additional points—such as the point target and point B—requires only a secondary phase correction relative to the distance from point A. In general, by initially focusing the collection of pulses on a subset of surface points in a coarse grid, it is possible to obtain focused waveforms for a much finer grid. This approach considerably reduces the number of operations needed to focus the entire surface, thus increasing the computational efficiency of the back-projection algorithm. The method proposed consists of two steps. 1) A first step aimed at performing the classic backprojection focusing [9] at specific surface points, namely reference surface points, which are typically separated in the order of 10 m. TABLE I SENTINEL-6 POSEIDON-4 ALTIMETER INSTRUMENT PARAMETERS 2) A second step where the stacks of corrected pulses before coherent summation from reference surface points are used to focus the waveforms at nearby surface points, namely the secondary surface points, obtaining a final focusing at points separated typically between 0.5 and 1 m for an altimeter such as Sentinel-6. Furthermore, an additional acceleration can be achieved by assuming that the phase slope is sufficiently small, so the azimuth phase remains relatively unchanged over a few pulses during the illumination time. In this case, we can average the matrix of corrected pulses from the reference surfaces into subapertures and store this reduced matrix for the secondary points. This approach further decreases the number of operations required to focus on the secondary points. III. THEORETICAL FRAMEWORK A. Transmitted Pulses Most high-range resolution radars, including altimeters, transmit chirp pulses [21], as described by the following equation: st(t)=wt(t)cosh2πfct−α 2t2i −Tp 2<t<Tp 2.(1) In this equation, wt(t)represents the pulse envelope, which is usually a rectangular window (uniform energy). Moreover, fcis the carrier frequency of the modulated signal, tis the time elapsed from the beginning of the pulse, which has a limited duration Tp, in the order of tenths of microseconds. The time variable tis known as fast time. The term αis the chirp rate, which is defined as the ratio of the pulse bandwidth Bto the pulse duration: α=B/Tp. The quadratic term in the phase signal is a representation of the linear frequency modulation of the pulse. Table Ishows the Sentinel-6 nominal values. After the transmission through a nadir-pointing antenna, the pulse travels to the Earth’s surface, where part of the energy is reflected and returned back to the sensor. B. Received Signal Fig. 4illustrates a simplified geometry of a radar altimeter. The satellite moves along the azimuth direction (y-axis), transmitting pulses with a nadir-pointing antenna in a nearly continuous way. A single scatterer is periodically reflecting the transmitted pulses during a limited illumination time Ti. The slow time ηrefers to the time relative to the position of the satellite for each emitted pulse. The order of magnitude of the 5209816 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 Fig. 4. Simplified radar geometry scheme. The satellite is moving and constantly transmitting pulses toward the nadir. Part of the energy of the transmitted pulses is reflected by the Earth’s surface and returned back to the satellite antenna. If a secondary point ysec is close enough to yref, their range migration difference R(η, yref)−R(η, ysec)can be approximated by a linear ramp. slow time is seconds, while the sampling rate is milliseconds, depending on the PRI. After reception, in modern radar altimeters such as Sentinel6, the signal is first digitized and then matched-filtered. The resulting waveforms for a scatterer located at ycan be expressed as follows [16]: Sr(η, fr,y) =wη(η, y)Wr(fr) ·exp j2πfc 2 cR(η, y)−2 c(R(η, y)−Rtrk)−fd(η, y) αfr  (2) where wηrefers to the azimuth envelope, defined by the antenna pattern. The range envelope is now defined in frequency domain Wr(fr), altered by geophysical effects. The frequency range fris the frequency variable in the range domain. The factor fd(η) =2fcvr/cindicates the alteration in the range of the phase due to the Doppler effect related to the relative velocity of the satellite vr. The tracker range Rtrk is commanded by the satellite onboard. In the case of Sentinel-6, the tracker range refers to the central point within the range window. The term R(η) is the range of the satellite to the target. The first term in the exponential is called relative range phase (RRP) and represents the phase rotation due to the range migration. The second term is the CW phase, whose frequency is proportional to the delay from the satellite to the scatterer. To focus all waveforms from the same scatterer, the back-projection algorithm compensates all phase terms in function of the range between the scatterer and the satellite, and then integrates all corrected echoes, achieving maximum coherence from a single scatterer on the surface. This has to be done for each defined point on the surface, which results in a high-computational effort. C. Radar Geometry Simplifications To determine a closed-form expression of the range migration difference between two points on the surface, we must assume that the distance between a reference along-track point yref and a secondary along-track point ysec is sufficiently small to ensure the validity of the geometric principles outlined in Section II. Under these tight conditions, we can assume some typical geometry simplifications (see [15],[16],[22]): the Earth is locally spherical with a radius Re, the satellite moves along the track with a constant velocity vs, and the difference between the satellite altitude zsat and the tracker range Rtrk is zero. D. Range Migration Difference Between Two Scatterer Points To obtain a closed-form expression of the range migration difference between scatterers, we approximate the range equation to a quadratic polynomial. It is important to remark that all phase corrections applied to a reference scatterer are made using real satellite data, and we do not need to assume any geometrical approximation, as explained in [9]. On the other hand, the azimuth phase correction applied to a secondary scatterer is designed under the quadratic approximation. Therefore, what we are really approximating is not the range equation relative to a single scatterer, but the range migration difference between two scatterers. For a point located at the along-track position y, the quadratic range equation can be defined as follows: R(η, y)=R0+vzη−y vg+v2 eq 2R0η−y vg2 (3) where R0is the range at the point of closest approach (POCA), vzis the vertical motion (orthogonal to the azimuth velocity), and veq is the equivalent speed defined as follows (see [22]): veq =√vs·vg=vs·sR0 R0+Re (4) where Reis the Earth radius and vgis the velocity of the platform projected on the ground. We can compute the range migration difference between a reference point yref and any along-track point yas a function of the azimuth distance (yref −y)in the following manner: R(η, y)=R(η, yref)+vs(yref −y) R0 η(5) where constant terms have been omitted, as they do not contribute to the azimuth phase coherence. Therefore, the range migration difference between very close scatterer points can be safely expressed as a linear function, because the nonlinear terms are negligible compared to the radar pulse wavelength. For an azimuth distance between scatterers of 100 m, using nominal values of Sentinel-6 (see Table I), the linear term of the range migration difference is in the order of 0.5 m, while nonlinear terms are in the order of micrometer and less, far below the order of the pulse wavelength, approximately 22 mm for the Sentinel-6 Ku-band. HERNÁNDEZ-BURGOS et al.: ACCELERATED ALGORITHM FOR BACK-PROJECTION FF SAR ALTIMETRY 5209816 E. Range Misalignment The most critical aspect of the simplifications made earlier is the residual misalignment of the across-track waveforms in range. Recalling the phase of the L1A waveforms from (2) θ(η, fr,y) =2πfc 2 cR(η, y)−2 c(R(η, y)−Rtrk)−fd(η, y) αfr. (6) The accelerated back-projection method corrects the RRP, which corresponds to the first term (the one multiplied by the carrier frequency fc), by applying a secondary phase correction. However, the range phase delay correction, commonly known as range cell migration correction (RCMC), is not applied to secondary surfaces. The rationale behind this omission is to avoid an additional complex correction and range compression (via FFT) for each secondary surface point. Since back-projection integrates all pulses within the illumination time to achieve proper focusing, it is crucial that all pulses align accurately at the same scatterer delay sample. Nonetheless, this residual error is unavoidable, and the only viable mitigation strategy is to ensure that the distance between the reference and the secondary point remains sufficiently small to keep the residual error within negligible limits. Within this context, if a range compression (range FFT) is applied to a secondary surface point under the assumption of uniform received energy, the resulting signal would be a delayed sinc function proportional to the range-migration difference between the reference and the secondary point, as well as the difference between their respective Dopplerdelay shifts Sr(η, fr,y)=sinc[B(t−1τ ·η)]wη(η, y) ·expj2πfc 2 cR(η, y)(7) where 1τ represents the range-migration difference, as defined in (5), plus the difference in the Doppler delay shift fd. We determine the criteria such that misaligned waveforms are still coherently integrated if the waveforms at the extremes of the synthetic aperture are within the across-track resolution δrof the waveform relative to the middle of the synthetic aperture. With these criteria, we could determine the maximum distance available from a reference point as follows: yref −y≤δr·R0 vsTi (8) where for the sake of simplicity we have omitted the Doppler delay shift fd, as for very close targets (less than 100 m), the impact of the latter on the maximum distance requirement is in the order of millimeters, thus considered negligible in this case. As evident, the maximum available azimuth distance between both points decreases as the integration time increases, given the increased range migration difference at the synthetic aperture extremes. In the case of Sentinel-6, the across-track resolution is approximately 0.42 m [16]. Taking the nominal altitude as the range of closest approach, 1336 km, the nominal speed of the Sentinel-6 satellite, 7200 m/s and using 3.4 s of illumination time (the maximum available for Fig. 5. Sentinel-6 L1A normalized power (dB) from a point target after range compression. (a) Power at the point target location. (b) Power at 40 m along-track distance from the point target location. The range misalignment is due to the range migration difference between both along-track points. the Doppler bandwidth without aliasing [16]), the maximum distance between a reference point and a secondary point to keep the waveforms aligned is around 23 m. Fig. 5presents a 2-D image of the waveforms from a Sentinel-6 pass over a transponder located in Crete [18] after aligning all waveforms. The left image (a) shows the aligned waveforms at the exact transponder location, while the right image (b) illustrates the miss-aligned waveforms at a 40-m along-track distance from the transponder. As observed, waveforms are slightly misaligned due to the range migration difference between the two along-track positions. The unavoidable miss-alignment leads to an inevitable loss of resolution in both dimensions, i.e., across-track and along-track. Assuming a symmetrical received power along the illumination time around the POCA, this miscorrection should not affect range measurement accuracy, as the positive and negative uncorrected delays cancel out during azimuth integration, effectively centering the focused pulse at the correct point. However, the received power from a scatterer point during the illumination time may not be symmetric, consequently, one of the sides of the synthetic aperture receives more power than the other extreme. As a result, a range bias increasing with integration time can be introduced when all waveforms are integrated. This effect is appreciable for very bright targets, such as transponders, and extreme back-projection configurations. For example, for 20 m of distance to the reference point and 3.4 s of integration time, the range bias is in the order of 3 mm. In the case of natural scenarios, this effect will depend on the geometry context but is expected to be much lower (below the millimeter). Thus, there is a trade-off between the spacing of reference points and the accuracy of the measurements. A detailed study of this effect is provided in the validation section. IV. METHODOLOGY As mentioned earlier, the reference surface points are focused using the classic back-projection method described in [9]. The secondary surface points are focused by reusing the set of corrected waveforms from a nearby reference point before integration, with a secondary phase correction applied, proportional to the distance between the reference and secondary points. In this section, we provide a summary of the classic back-projection method in the case of Sentinel-6, along 5209816 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 with a detailed explanation of the mathematical development behind the secondary focusing. A. Reference Surface Points: Classic Back-Projection The back-projection algorithm is explained in [9] in the case of CryoSat-2, which employs a de-ramp filter instead of matched filtering as in the Sentinel-6 case. The only difference between them is the absence of the residual video phase (RVP) term in the matched filter case. The rest of the approach is applicable. To give some context, a summary of the classic back-projection is explained. The FF-SAR back-projection algorithm for a matched-filter case system consists therefore of three main well-known steps. 1) Range Cell Migration Correction: The RCMC involves correcting the scatterer delay difference between a point located at the along-track position ypt and the position of the satellite for each pulse. The range-migration is defined as δR(η, ypt)=R(η, ypt)−R0, where R0is the range at the POCA. We can express the RCMC as follows: SRCMCη, τ, ypt=Srη, fr,ypt ·exp" 2 cδRη, ypt−fdη, ypt α!fr#. (9) After this multiplication, a range-compression is applied using a Fourier Transform. The 2-D signal in the temporal domain is Srcη, τ, ypt≈TpsincB(τ−γ0) ·expj2πfc 2 cRη, ypt(10) where γ0=(R0−Rtrk)·c/2. The remaining phase term is the RRP, which is the dominant phase term. 2) Along-Track Focusing: The along-track focusing involves correcting the RRP and integrating all the coherent waveforms. This is done for all the points defined on the surface. For a point target located at ypt, the azimuth correction is defined as follows: SATη, τ, ypt|y=Srcη, τ, ypt·exp−j2πfc 2 cR(η, y) =TpsincB(τ−γ0) ·expj2πfc 2 cRη, ypt−R(η, y)(11) where yrefers to a specific along-track point on the surface. If the evaluated point ycoincides with the position of the point target ypt, the RRP is fully corrected, and the energy after integrating all the pulses will be maximum. The 2-D PTR can be computed as the integration of all waveforms in the azimuth direction after the RRP correction, in function of the range delay τand the along-track point y χpt(τ, y)=X η∈Ti SATη, τ, ypt|y ≈TpTisincB(τ−γi)sinc2Tifcvs cR0y−ypt. (12) B. Secondary Focusing Starting from (11), we can express the azimuth-corrected signal evaluated at a reference surface point yref as follows: SATτ, η, ypt|yref=TpsincB(τ−γ0) ·expj2πfc 2 cRη, ypt−R(η, yref). (13) Substituting the range migration difference defined in (5) into the latter equation, we can express the azimuth-corrected signal as a function of the waveforms evaluated at a reference surface point multiplied by an exponential linear phase term SATτ, η, ypt|yref=TpsincB(τ−γ0) ·exp"j2πfc 2vsyref −ypt cR0 η#.(14) Thus, for any secondary along-track point y, we can express the azimuth focusing as the waveforms corrected at the reference surface point (latter equation) multiplied by a secondary phase correction, proportional to the distance between yand the reference point yref SATτ, η, ypt|yref →y=SATτ, η, ypt|yref ·expj2πfc 2vs(y−yref) cR0 η =TpsincB(τ−γ0) ·exp"j2πfc 2vsy−ypt cR0 η#.(15) In summary, we apply the RCMC and RRP corrections relative to a reference point using real satellite data to obtain SAT(τ, η, ypt|yref). Then, we apply the secondary phase correction for a secondary along-track point y, close to the reference point yref, obtaining an expression equivalent to (11) SATτ, η, ypt|yref →y≈SATη, τ, ypt|y.(16) After that, we can obtain the 2-D PTR by integrating all waveforms in the azimuth direction, yielding the same result as in (12). The principal advantage is that, by utilizing a single reference-corrected matrix (14),SAT(τ, η, ypt |yref), we can focus multiple secondary points with only an additional secondary correction (15),SAT(τ, η, ypt |yref →y). The computational cost of performing the classic back-projection for every along-track point y(11) is much larger than performing a single reference correction (14) and then apply as many secondary corrections (15) as needed. It is important to note that the linear approximation is only used for the secondary phase correction (15), which is valid for azimuth distances below ∼23 m, as previously demonstrated. C. Subaperture Averaging An additional step to enhance the computational efficiency of the FF-SAR back-projection algorithm involves averaging stacks of Nsub waveforms, previously aligned with respect to a reference scatterer and storing this reduced matrix for the secondary phase correction and azimuth integration. In this HERNÁNDEZ-BURGOS et al.: ACCELERATED ALGORITHM FOR BACK-PROJECTION FF SAR ALTIMETRY 5209816 Fig. 6. Simulated point target power loss (dB) using Sentinel-6 nominal parameters for different distances to reference points, in function of the stack size. approximation, we assume that the RRP is constant during a certain number of consecutive pulses, allowing the integration of all the pulses to remain coherent without needing to apply the RRP correction. The use of subapertures to reduce computational effort has already been studied for SAR imaging, specifically in the fast factorized back-projection algorithm [17]. The focusing process is therefore divided into two integrations: one before the extra azimuth correction for a small time of subapertures Tsub and one after the secondary phase correction for the entire integration time but with fewer points. The length of each stack should be small enough so the RRP within the subaperture remains sufficiently constant to have a coherent summation. Starting from (14), this process consists of integrating the pulses within a shifting window. The integration is followed by a reduction step, where the resulting vector is scaled proportionally to the window size SAT,Tsub τ, η, ypt|yref =Tp Tsub sincB(τ−γ0) · Nsub−1 X k=0(exp"j2πfc 2vsyref −ypt cR0·η+k·Tpri#) (17) where Tpri is the pulse repetition interval in seconds. The result of this integration is a 2-D matrix with a reduction factor proportional to the stack length. The slow-time variable ηnow has a pulse repetition interval proportional to the subaperture time Tsub. After the stack averaging and reduction, the secondary points are focused as usual, applying a secondary phase correction and integrating all the pulses along track χpt(τ, y)=X η∈TillSAT,Tsub τ, η, ypt|yref ·expj2πfc 2vs(y−yref) cR0 η.(18) If the phase term after the stack averaging were zero, the integration would be fully coherent, and we would obtain the same PTR as in (12). To assess the effect of the subaperture approximation, we must analyze the coherence loss of the signal due to the remaining RRP phase that was not corrected when applying the stack averaging. Thus, evaluating the summary of (17) as a continuous integral SAT,Tsub (τ, η, y|yref) ≈Tp Tsub sincB(τ−γ0) ·ZTsub/2 −Tsub/2 exp"j2πfc 2vsyref −ypt cR0η+η′#dη′.(19) We can express the reference-corrected matrix after secondary focusing as follows: SAT,Tsub τ, η, ypt|yref →y=sincB(τ−γ0) ·sinc"fc 2vsyref −ypt cR0 Tsub# ·expj2πfc 2vs(y−yref) cR0 η. (20) The second sinc function is responsible for the subaperture approximation. As observed, the power decays rapidly with respect to the subaperture coherence time Tsub and the along-track distance between the point target and the reference point. Moreover, it is essential to note that this term depends on the distance between the reference point and the point target location, but not on the secondary surface points (y). This implies that the focused PTR will remain unchanged, as the azimuth exponential phase is still maintained to perform the focused response. Thus, the along-track resolution remains dependent on the total illumination time Ti, irrespective of the stack length. In fact, if Tsub is set to 64 pulses and we perform a power (incoherent) integration along the aperture, we would obtain the well-known unfocused SAR altimetry response [6], [9],[15]. For the FF PTR, stack averaging directly impacts the signal-to-noise ratio (SNR), leading to a noisier PTR. Fig. 6illustrates the power decay at the point target as we increase the coherence integration time and the distance to the reference point. For a 20 m distance to the reference point, we observe a power loss of 1 dB for a coherence time of 27.5 ms, which corresponds to a stack length of Nsub =247 pulses. However, the improvement in computational efficiency saturates close to 30–40 pulses, resulting in a power loss of less than 20 mdB. This power decay is influenced by satellite and radar altimeter parameters such as platform velocity and PRF, meaning these values will vary for other radar altimeter missions, such as Sentinel-3 and CryoSat-2. D. Algorithm Implementation This section outlines the steps involved in implementing the accelerated back-projection algorithm. Fig. 7shows the flowchart of the back-projection algorithm, and Fig. 8illustrates the algorithm functions scheme for both classic and 5209816 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 Fig. 7. Algorithm flowchart of the accelerated back-projection. For each grid point, we first process a reference surface (green path) and compute its waveform. The classic back-projection is applied and the waveforms before full integration are stored in a buffer. The size of the waveforms depends on the number of pulses of the stack averaging size Nsub and the zero-padded range samples Nrzp . All secondary points associated with a single reference are focused afterward (purple paths). This is repeated for all points defined in the grid. The output is a range waveform for each point on the surface. Fig. 8. (Top) FF-SAR accelerated back-projection algorithm scheme. (Bottom) Each secondary point (in purple) is associated with the closest reference point (in red) in the surface points. accelerated back-projection focusing procedures. First, all surface points are defined by projecting the satellite coordinates into the ground. Then, each surface point is classified either as a reference point or a secondary point. Then, we iterate over each defined point, where the reference points are focused using the classic back-projection method, as described in [9] and the secondary points are focused as described in this article. Therefore, given a reference surface point, all back-projection corrections are applied, and the set of phase-corrected waveforms prior to azimuth integration is saved. The reference point is focused on applying the azimuth integration as usual. The saved matrix is then averaged and decimated by a factor of Nsub (stack averaging). The reduced matrix is used to focus all secondary points nearby to the reference surface point, applying a linear azimuth phase correction. Subsequently, the collection of corrected waveforms is coherently integrated as usual, yielding a secondary-focused look. The distance between reference along-track points is configurable, depending on user preferences, along with other parameters such as along-track spacing, azimuth integration time, and subaperture integration length. V. ALGORITHM CHARACTERIZATION To characterize the limits and performance of the accelerated back-projection, we have conducted a test with Sentinel-6 data over a transponder located in Gavdos [20] with different processing configurations. Active elements such as transponders offer a high signal-to-clutter ratio (SCR) and allow us to characterize the real PTR and compare different algorithms. Additionally, conducting repeated tests on the same transponder over time enables us to determine longterm statistics, such as the range bias, in comparison between different algorithms and measurement methods. In this section, results from the accelerated back-projection are compared with the classic back-projection algorithm. The performance of the classic back-projection algorithm over transponders and corner reflectors has been already shown in [23] (transponder) and [19] (corner-reflector). HERNÁNDEZ-BURGOS et al.: ACCELERATED ALGORITHM FOR BACK-PROJECTION FF SAR ALTIMETRY 5209816 Fig. 9. 2-D PTR normalized power in dB from the transponder located in Gavdos. Accelerated back-projection has been processed over the different distances between reference points, noted at the top of each image. For all cases, 3.4 s of synthetic aperture has been integrated. The point target is located in the worst case scenario: in the middle between two reference points. The 2-D PTR degrades as we increase the distance between reference points, as expected. Fig. 10. Sentinel-6 PTR normalized power (dB) from the transponder located in Gavdos, processed with classic and accelerated back-projection algorithm for different distances between reference surfaces and 3.4 s of integration time. (Left) Across-track cut. (Right) Along-track cut. Fig. 11. PTR resolution degradation from the transponder (Gavdos), processed with classic and accelerated back-projection algorithm for different distances between reference points and 3.4 s of integration time. The resolution degradation is relative to the reference value obtained with the classic back-projection. To perform the characterization tests, the transponder coordinates are located in the middle between two reference points, that is the farthest secondary point with respect to a reference point and, therefore, the worst case of the configuration in terms of residual effects. All tests have been performed with 3.4 s of integration time. Moreover, the distance between along-track surface points is 0.05 m and the range waveforms have been oversampled with a factor of 8, which results in a range bin spacing around 0.05 m. Any antenna pattern correction has been applied. Fig. 12. PTR power decay in dB with respect to the classic back-projection from the transponder located in Gavdos, for different distances between reference points and 3.4 s of integration time. Fig. 13. Peak sidelobe (PSL) from the transponder located in Gavdos, processed with classic and accelerated back-projection algorithm for different distances between reference points and 3.4 s of integration time. The PSL is calculated relative to the maximum power. A. PTR: Reference Point Distance Analysis Fig. 9illustrates the degradation of the focused 2-D PTR as we increase the distance between reference points, noted at the top of each image. At 40 m, the expected degradation of the 2-D PTR becomes visible. Fig. 10 showcases the cross and along-track cuts at the point of maximum power. As we increase the distance between reference points, the sinc response gets wider, absorbing energy from the sidelobes. Fig. 11 shows the across-track and along-track resolution error with respect to the classic back-projection in the function of the distance between two reference points. As we can see, the error increases rapidly as we increase the distance between reference points above 40 m. Moreover, Fig. 12 depicts the power loss at the peak of the PTR, while Fig. 13 5209816 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 63, 2025 Ferran Gibert (Member, IEEE) received the M.S. degree in aerospace engineering and the Ph.D. degree in aerospace science and technology from the Universitat Politècnica de Catalunya, Barcelona, Catalonia, in 2011 and 2016, respectively, focused on the design of thermal diagnostics experiments for the LISA Pathfinder mission. From 2016 to 2018, he was a Post-Doctoral Researcher with the University of Trento, Trento, Italy, from where he supported the LISA Pathfinder operations and data analysis phase. Since June 2018, he has been a Research and Development Engineer with isardSAT, S.L., Barcelona, being currently involved in the development of algorithms for radar altimeter processors. Antoni Broquetas (Life Member, IEEE) was born in Barcelona, Catalonia, in 1959. He received the Engineering and Doctor Engineering degrees in telecommunication engineering from the Universitat Politècnica de Catalunya (UPC), Barcelona, Catalonia, in 1985 and 1989, respectively, where his dissertation was on microwave tomography. In 1986, he was a Research Assistant with Portsmouth Polytechnic, Portsmouth, U.K., where he was involved in propagation studies. In 1987, he joined the Department of Signal Theory and Communications, School of Telecommunication Engineering, UPC. From 1998 to 2002, he was the Subdirector of Research at the Institute of Geomatics, Barcelona. Since 1999, he has been a Full Professor with UPC, where he is involved in research on radar imaging and remote sensing. From 2003 to 2006, he was the Director of the Department of Signal Theory and Communications, UPC. From 2005 to 2009, he proposed and comanaged the Knowledge To Market (K2M) Program at UPC to support university spinoffs. He has published more than 200 articles on microwave tomography, radar, ISAR, and synthetic aperture radar (SAR) systems, SAR processing, and interferometry. His present research interests include SAR calibration techniques, GMTI SAR performance analysis, GEOSAR, and mmW groundbased SAR, in support of future SAR missions and radar remote sensing innovation. Pol Villalvilla received the M.S. degree in telecommunications engineering from the Universitat Politècnica de Catalunya, Barcelona, Catalonia, in 2020. From 2021 to 2022, he completed a one-year research stay at the Earth Observation and Remote Sensing Department, Swiss Federal Institute of Technology (ETH Zürich), Zürich, Switzerland, performing signal processing techniques over measurements in alpine areas to assist Persistent Scatterer Interferometry. Since 2022, he has joined isardSAT, S.L., Barcelona, to work on the development and implementation of processing algorithms of radar altimetry for future ESA missions. Alejandro Egido received the master’s degree in telecommunications engineering from the Universidad de Zaragoza, Zaragoza, Spain, in 2007, and the Ph.D. degree in telecommunications engineering from the Universitat Politècnica de Catalunya, Barcelona, Spain, in 2013. In 2007, he joined Starlab Barcelona, Barcelona, as a Scientific Researcher. In 2015, he moved to the Laboratory for Satellite Altimetry, National Oceanic and Atmospheric Administration (NOAA), College Park, MD, USA, where he worked on developing improved SAR altimetry products and algorithms for physical oceanography and sea ice applications. In 2016, he was a NOAA’s Measurement System Engineer for the Jason altimeter missions. Since 2022, he has been an Ocean and Hydrology Earth Observation Scientist in the Earth Surface and Interior Section of European Space Agency’s Mission and Science Division, based at the European Space Research and Technology Centre, Noordwijk, The Netherlands. His main research interests include using radar altimetry to characterize open ocean, coastal, and polar zones, and leveraging global navigation satellite system signals as sources of opportunity for remote sensing applications. Dr. Egido received the NOAA’s David S. Johnson Award for his work on fully focused SAR altimetry in 2020. Gorka Moyano received the M.Sc. degree in telecommunications engineering and the M.Sc. degree in research on information and communication technologies from the Universitat Politècnica de Catalunya, Barcelona, Catalonia, in 2009. From 2010 to 2013, he worked at Indra Espacio, Barcelona, as a Radionavigation Software Engineer, where he developed GNSS monitoring systems for air traffic control and automated testing tools for the Galileo Sensor Station Core Computer. In 2013, he joined isardSAT, S.L., Barcelona. He is currently a Software Architect Engineer liaising scientific algorithms and software development. His skills are key in projects, such as Sentinel-6 P4 L1 Ground Processor Prototype (GPP) development, the Altimeter L1 Product Generation Function (PGF) processors within the Sentinel-6 PDAP, the CRISTAL IRIS L1&GR GPP, and the EPS-SG MWS/MWI/ICI L1b Local Processors reengineering. Adrián Flores de la Cruz received the M.S. and Ph.D. degrees in telecommunications engineering from the Universidad Politécnica de Cartagena, Murcia, Spain, in 2014 and 2019, respectively. From 2017 to 2019, he was with Hisdesat SA, Madrid, Spain, as a Cal/Val Engineer, where he worked on the PAZ Mission. Since March 2019, he has been a Research and Development Engineer with isardSAT, S.L., Barcelona, Catalonia, working on the development and implementation of data processing techniques with a special focus on instrumental calibration. Marco Fornari received the M.Sc. degree in telecommunications engineering from the University of Rome “La Sapienza,” Rome, Italy, in 2006. He then joined the Earth Observation Directorate, ESA-ESTEC, Noordwijk, The Netherlands, first as a Young Graduate Trainee and then as a RHEA/Starion contractor until present. Since 2007, he has been involved in several ESA altimetry missions. He was initially part of the CryoSat-2 Project as Instrument Calibration Engineer, and still supports CryoSat-2 Operations via the ESA-ESRIN Team. In 2011, he joined the Sentinel-6 Project as a Ground Segment Engineer, leading the development of the altimeter ground prototype processors and still working on their future evolutions, such as the fully focused processing. More recently, he has been supporting CRISTAL and the S3NG-Topo Projects as an altimeter instrument and data processing expert.