Ocean Wave Inversion Based on Hybrid Along- and Cross-Track Interferometry

Abstract

The hybrid interferometric synthetic aperture radar system is a combination of an along-track configuration and cross-track configuration. Based on linear ocean wave theory, an ocean wave inversion algorithm for a hybrid interferometric synthetic aperture radar system is proposed in this work. Using the interferometric synthetic aperture radar images acquired by the TerraSAR-X and TanDEM-X satellites and the interferometric synthetic aperture radar images acquired by an airborne interferometric radar altimeter with a certain degree of squint, the profile of ocean waves and the corresponding orbital velocities were retrieved by combining the new inversion algorithm with the cross-spectra. Meanwhile, key parameters of ocean waves, such as the significant wave height, significant wave orbital velocity, propagation direction, and wavelength of the dominant waves, were also extracted from the ocean wave spectra retrieved in this study. In order to evaluate the reliability of the new inversion algorithm, the retrieved significant wave heights were compared with those provided by the European Centre for Medium-Range Weather Forecasts and measured by a Global Navigation Satellite System buoy. The results showed that for the ocean waves retrieved from the spaceborne hybrid interferometric synthetic aperture radar images, the differences between the retrieved significant wave heights of the four subareas selected in this paper and those provided by European Centre for Medium-Range Weather Forecasts were approximately 0.01, –0.17, –0.55, and –0.37 m, respectively, and for the ocean waves retrieved from the airborne interferometric radar altimeter images, the differences between the retrieved significant wave heights corresponding to the M920 and M3120 images used in this paper and those measured by the Global Navigation Satellite System buoy were approximately –0.05 and –0.09, respectively. Therefore, the method proposed in this work could retrieve the ocean wave spectra well when the velocity bunching had a small influence; however, as the nonlinear influence of the velocity bunching increased, the difference between the significant wave heights retrieved using this method and provided by the European Centre for Medium-Range Weather Forecasts also increased.

1. Introduction

As one of the most common phenomena on sea surfaces, ocean waves have an important influence on ocean energy exchange and ocean–atmosphere interaction. Therefore, theoretical research and experimental observations of ocean waves have always been an important research topic in the marine sciences. Ocean waves can be measured in many ways such as in situ measurement with buoys and remote measurement with different microwave radars. In 1996 and 2018, Hasselmann et al. and Yu et al. improved existing algorithms of ocean wave information retrieval based on synthetic aperture radar (SAR) image spectra and traditional altimeter data, respectively [1,2]. In 2021, on the basis of emerging interferometric radar altimeter (IRA) images, Jiang et al. retrieved two-dimensional ocean wave spectra, and the wind wave information was also effectively extracted [3]. In addition, in 1991 and 2022, Shemer et al. and Li et al. introduced sea surface simulation methods corresponding to an interferometric synthetic aperture radar (InSAR) system and an SAR system, respectively [4,5]. For in situ measurements, it has the advantage of high accuracy; however, its applicability to extensively measure ocean waves is seriously hampered due to the limited measurement space and high cost.

SAR images can be used for ocean wave inversion due to the fact of their high resolution [6,7]. In 1991, on the basis of hydrodynamic modulation, tilt modulation and velocity bunching modulation, Hasselmann et al. derived a nonlinear mapping that can be used to transform ocean wave spectra into SAR image spectra [8]. In 1995 and 2005, on the basis of the cross-spectra of SAR images, Engen et al. and Schulz-Stellenfleth et al. successfully retrieved ocean waves using different methods [9,10]. In 2002 and 2015, Marghany et al. and Ren et al. retrieved the significant wave height (SWH) using the azimuth cutoff of ERS-1 images and the azimuth cutoff of RADARSAT-2 images, respectively [11,12]. In 2004 and 2006, based on polarimetric SAR images, Schuler et al. and He et al. used two different methods to effectively retrieve the sea surface slope, respectively [13,14]. In 2004, He et al. successfully retrieved ocean wave spectra based on a linear polarimetric SAR image [15], and Zhang et al. validated the measurement of ocean waves with RADARSAT-2 fully polarimetric SAR images in 2010 [16]. Moreover, Wang et al. evaluated the observation of swell height of the Sentinel-1A/B wave mode in combination with buoy in situ measurement data in 2022 [17]. Unfortunately, currently, the hydrodynamic modulation and relaxation rate of the hydrodynamic modulation are not well understood [18]. In 2020, on the basis of polarized SAR (PolSAR) images, Yanmin Zhang et al. [19] derived a new nonlinear mapping that can be used to map ocean wave spectra into polarimetric SAR image spectra.

Compared with conventional SAR, InSAR has unique advantages in measuring ocean waves. When the InSAR system was originally applied in the marine technology field, it was used to measure surface currents. Since the phase obtained by the along-track InSAR (ATI-SAR) is proportional to the orbital radial velocity of the corresponding sea surface scatterer, Goldstein et al. proposed that the images of ATI-SAR can be used to retrieve the surface currents [20,21]. Marom et al. carried out two airborne ATI-SAR experiments in the Marina Beach region in 1990 and the Gulf of California in 1991. The two-dimensional ocean wave spectra retrieved based on the phase images of these two experiments were consistent with the buoy in situ measurement results [22,23]. In 1991, Lyzenga et al. pointed out that the dual or multi-antenna ATI-SAR could be used to solve the disadvantages of the azimuth wavenumber cutoff due to the velocity bunching [24]. Lyzenga and Malinas further verified this conclusion based on the measured ATI-SAR phase images in 1996 and showed that the azimuth wavenumber was directly proportional to the baseline length between two antennas [25]. In order to study the ocean wave imaging mechanism of ATI-SAR, Bao et al. derived a two-dimensional ocean wave imaging model of ATI-SAR that included the modulation of normalized radar cross-sections (NRCS) by long waves, velocity bunching modulation, and azimuth image smear caused by orbital acceleration of long waves and orbital velocity spread within the resolution unit [26]. Meanwhile, the intensity and the phase image spectra of ATI-SAR at different sea state conditions and radar configurations were also calculated based on the Monte Carlo method. Due to the large number of calculations used in the Monte Carlo method, Bao et al. derived a new nonlinear integral transform between ATI-SAR phase image spectra and ocean wave spectra in 1999 that could be used to simulate the corresponding phase image spectra from ocean wave spectra [27]. In 1999, Vachon et al. validated the nonlinear integral transform between ATI-SAR phase image spectra and ocean wave spectra derived by Bao et al. by using C-band horizontally polarized airborne ATI-SAR phase images, and they further obtained two quasi-linear integral transforms using Taylor series expansion [28]. However, the δ function was omitted in the expression of the ATI-SAR phase image established by Bao et al., resulting in the loss of the derivative term of the wave orbital radial velocity in the established nonlinear integral transform. In order to solve this problem, based on the expression of the ATI-SAR phase image, He et al. added the δ function omitted when Bao et al. established the nonlinear integral transform and established a more precise nonlinear integral transform between ATI-SAR phase image spectra and ocean wave spectra, in 2003 [29]. Based on the nonlinear mapping relationship between ATI-SAR phase spectra and ocean wave spectra, Zhang et al. developed a parametric inversion model for retrieving ocean wave direction spectra from ATI-SAR phase images in 2009, and they successfully retrieved ocean wave wavelengths, ocean wave directions, and SWHs [30].

On the basis of a cross-track InSAR (XTI-SAR) imaging mechanism of ocean waves, Bao derived a nonlinear integral transform that can be used to map ocean wave spectra into phase image spectra, in 1999 [31]. In order to simulate the traditional SAR intensity image spectra and phase image spectra, Schulz-Stellenfleth et al. established a wave imaging model on the basis of a cross-track interference mechanism in 2001 [32]. In addition, Schulz-Stellenfleth et al. also estimated the distorted (bunched) variance spectra of digital elevation models (DEMs) on the basis of the Monte Carlo method and forward simulations. It was shown that a straightforward estimation of sea surface elevation using the distorted InSAR elevation model was well consistent with the true sea surface elevation for low-amplitude swells, and the error depended on the water wave propagation direction and the coherence time of the echoes. Schulz-Stellenfleth et al. retrieved the bunched DEMs using airborne X-band horizontal polarization XTI-SAR wind wave phase images near the North Sea in the same year [33]. The results showed that when the nonlinearity caused by velocity bunching was small, the SWH and one-dimensional wavenumber spectra retrieved from the models were consistent with the in situ measurements with buoys. However, the difference between the propagation direction of the retrieved dominant waves and the buoy measurements was up to 30°, which depended on the flight direction and the flight level. In addition, since XTI-SAR images can simultaneously provide information on sea surface height and normalized radar cross-sections (NRCS), Schulz-Stellenfleth et al. proposed a new approach to measure the modulation transfer function (MTF) of real aperture radar (RAR) using cross-correlating bunched DEMs and registered SAR intensity images.

However, due to the technical limitations, most baselines of real InSAR systems contain both the cross-track and along-track configurations, which was defined as a hybrid InSAR system in this paper, and the shuttle radar topography mission (SRTM) is a typical case of a hybrid InSAR system [34]. In addition, most airborne InSAR systems use hybrid geometry which introduces squint angles to solve the effect of crosswind. Lehner et al. observed waves using an airborne hybrid InSAR system consisting of three antennas and successfully retrieved the corresponding ocean wave spectra, in 2001 [35]. Siegmund et al. successfully measured sea surface currents using an airborne hybrid InSAR system. However, since the inertial navigation data error will introduce a systematic error into the interferometric data, and the data accuracy of the different inertial navigation systems may be different, therefore, the accuracy of the measured currents is related to the accuracy of the navigation data of the platform [36,37]. In addition, in 2001, Bao et al. briefly summarized ocean wave spectra measurements using SRTM images, which mainly demonstrated the possibility of observing large-scale swell features. However, the wind waves could not be observed due to the low spatial resolutions of the SRTM images [38]. In addition, Bao et al. did not describe in detail the method for ocean wave inversion using the hybrid interferometry system.

In this paper, a spaceborne hybrid InSAR system, composed of two satellites (i.e., TerraSAR-X and TanDEM-X) and an airborne hybrid InSAR system with a small incidence angle, were used to retrieve ocean wave information including ocean wave spectra, SWH, and radial velocities. In addition, in this paper we took the SWHs provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) and the SWH measured by Global Navigation Satellite System (GNSS) buoys as the reference for evaluating the accuracy of the SWHs measured by the spaceborne hybrid InSAR system and the airborne hybrid InSAR system in the present work [39,40]. In addition, the cross-spectra of ocean waves were also retrieved to determine the propagation directions of ocean waves retrieved by the spaceborne hybrid InSAR system.

This paper is organized as follows: Section 2 introduces the basic theories of sea surface profile and radial velocities of ocean waves measurements by the hybrid InSAR system. In Section 3, the hybrid InSAR images acquired by the TerraSAR-X and TanDEM-X satellites at Pentland Firth, which is north of mainland Scotland, and the hybrid InSAR images acquired by airborne IRA from the Qingdao Xiaomaidao (XMD) sea area are introduced briefly. The data preprocessing of the hybrid InSAR images is presented in Section 4. The retrieved results, such as the ocean wave spectra, cross-spectra, and radial velocities of the ocean waves, are presented in Section 5. The conclusions of the present work are given in Section 6.

2. Algorithm of Ocean Wave Inversion

Figure 1 is the geometric schematic of the hybrid interferometry. The x-axis is the along-track (azimuth) direction, which is also the direction of the platform velocity, V. The y-axis is the cross-track (range) direction. P₁ and P₂ denote the positions of the two antennas. B is the baseline length, which needs to be decomposed into a cross-track baseline component (B_v) and an along-track baseline component (B_p) in the present work. α denotes the roll angle of B_v relative to the y-axis, and T is the scattering element on the sea surface.

2.1. Inversion Algorithm of the Sea Surface Profile

We can retrieve the sea surface profile using XTI-SAR images and ocean wave velocity using ATI-SAR images. The sea surface profile inversion and the radar line-of-sight velocity induced by the orbital velocity of long waves inversion based on InSAR images are given by [41,42]:

$$h=-\frac{\lambda r\sin \theta }{4\pi B_{\mathrm{XTI}}\cos \left(\theta -\alpha \right)}{\phi }_{elevation}$$

(1)

$$v=-\frac{\lambda V}{4\pi B_{\mathrm{ATI}}}{\phi }_{velocity}$$

(2)

where h denotes the profile of the simulated ocean waves; v denotes the radar line-of-sight velocity induced by the ocean waves; θ denotes the incidence angle; λ denotes the wavelength of the incidence microwave wave; B_XTI denotes the cross-track baseline; B_ATI denotes the along-track baseline; φ_elevation and φ_velocity are the interference phases caused by sea surface elevation and ocean wave radial velocity, respectively.

However, when the baselines of InSAR systems contain both along-track and cross-track components, we must use another method to acquire the ocean wave information.

Based on linear sea wave theory, the profile of ocean waves can be expressed as the superposition of a series of traveling harmonic waves:

$$h\left(r,t\right)={\displaystyle \sum _k{\zeta }_k\exp \left(i\left[kr-\omega t\right]\right)}+c.c$$

(3)

where ${\zeta }_k$ is the Fourier transform coefficient of the sea surface profile; $\omega =\sqrt{gk}$ denotes the angular frequency of the gravity waves; g and k denote the gravitational acceleration constant and the wavenumber, respectively; c.c stands for the complex conjugate. Using Equation (3), the radar line-of-sight velocity induced by the ocean waves can be expressed as:

$$v\left(r,t\right)={\displaystyle \sum _k{T}_k^v{\zeta }_k\exp \left(i\left[kr-\omega t\right]\right)}+c.c$$

(4)

where $T_k^v$ denotes the range velocity transfer function, which can be expressed as:

$$T_k^v=-\omega \left(\sin \theta \frac{k_r}{\left|k\right|}+i\cos \theta \right)$$

(5)

where k_r denotes the wavenumber component in the cross-track direction.

The interference phase acquired by the hybrid InSAR system contains both a cross-track interference phase component and an along-track interference phase component [35,36,37,38], which can be expressed as:

$$\phi ={\phi }_v+{\phi }_p=a_{0}h+b_{0}v$$

(6)

where φ is interference phase acquired by the hybrid InSAR system; φ_v and φ_p are the cross-track interference phase component and along-track interference phase component, respectively. a₀ and b₀ in Equation (6) can be given by:

$$a_0=-\frac{4\pi B_v\cos \left(\theta -\alpha \right)}{\lambda R\sin \theta }$$

(7)

$$b_0=-\frac{4\pi B_p}{\lambda V}$$

(8)

where R is the slant range from antenna to the scattering element on the sea surface. Inserting Equations (3) and (4) into Equation (6) yields:

$$\begin{array}{ll}\hfill \phi & =a_0\left\{{\displaystyle \sum _k{\zeta }_k\exp \left(i\left[kr-\omega t\right]\right)+\mathrm{c}.c}\right\}+b_0\left\{{\displaystyle \sum _k{T}_k^v{\zeta }_k\exp \left(i\left[kr-\omega t\right]\right)+c.c}\right\}\\ \hfill & ={\displaystyle \sum _k\left(a_0+b_0{T}_k^v\right){\zeta }_k\exp \left(i\left[kr-\omega t\right]\right)}+c.c\end{array}$$

(9)

According to Equation (9), we can derive that:

$$\left(a_0+b_0{T}_k^v\right){\zeta }_k=\mathrm{F}\left(\phi \right)$$

(10)

where F(φ) is the Fourier transform of interference phase. Then, the Fourier transform coefficient of the sea surface profile (i.e., ${\zeta }_k$) and the sea wave-height spectrum can be obtained as:

$${\zeta }_k=\frac{\mathrm{F}\left(\phi \right)}{a_0+b_0{T}_k^v}$$

(11)

$$P(k)=\frac{{\left|{\zeta }_k\right|}^2}{\mathrm{d}{k}_x\mathrm{d}{k}_y}$$

(12)

Using Equation (11), the sea surface profile is retrieved by taking the inverse Fourier transform of ${\zeta }_k$:

$$h{=\mathrm{F}}^{-1}\left({\zeta }_k\right){=\mathrm{F}}^{-1}\left(\frac{\mathrm{F}\left(\phi \right)}{a_0+b_0{T}_k^v}\right)$$

(13)

2.2. Orbital Velocity of Long Waves Inversion Algorithm

Using Equations (6) and (13), the along-track interference phase can be calculated by:

$${\phi }_p=\phi -{\phi }_v=\phi -a_0{\mathrm{F}}^{-1}\left(\frac{\mathrm{F}\left(\phi \right)}{a_0+b_0{T}_k^v}\right)$$

(14)

Then, the radar line-of-sight velocity induced by the ocean waves is obtained as:

$$v=\frac{\phi -a_0{\mathrm{F}}^{-1}\left(\frac{\mathrm{F}\left(\phi \right)}{a_0+b_0{T}_k^v}\right)}{b_0}$$

(15)

Meanwhile, the spectrum of the radar line-or-sight velocity is also obtained as:

$$P_v(k)=\frac{{\left|T_k^v{\zeta }_k\right|}^2}{\mathrm{d}{k}_x\mathrm{d}{k}_y}$$

(16)

3. Experimental Data

3.1. Hybrid InSAR Images

Germany launched an earth observation satellite (EOS), TerraSAR-X, in June 2007. The main payload of the satellite is an X-band radar sensor with various operation modes, which allows it to record images with different swath widths, resolutions, and polarization modes, and it provides high-quality X-band SAR images. Another satellite, TanDEM-X, was launched in June 2010, the configuration of which is almost identical to that of TerraSAR-X. The spiral formation flight of TerraSAR-X and TanDEM-X form a configurable spaceborne binary InSAR system. These two satellites cooperate with each other to perform earth topography missions to obtain high-precision global DEMs.

The two satellites, TerraSAR-X and TanDEM-X, simultaneously passed over Pentland Firth at 06:41 on 19 March 2012, which formed a spaceborne hybrid InSAR system, and we could carry out interference processing with SAR images obtained in this sea area to retrieve various marine phenomena such as sea surface current field and ocean waves. Romeiser et al. successfully retrieved the sea surface current field [43] with the SAR images acquired by TerraSAR-X and TanDEM-X, and the purpose of this paper was to use the SAR images to retrieve the ocean waves.

Figure 2a is the geographic position of Pentland Firth, where the red box is the specific location of all of the SAR images used in this paper. In the present work, we selected four subareas in the images for ocean wave inversion. The positions of the four subareas are marked with boxes in different colors in Figure 2b, and the red, yellow, green, and purple boxes represent Area 1 (A₁), Area 2 (A₂), Area 3 (A₃), and Area 4 (A₄), respectively. The detailed parameters of the SAR images are given in Table 1. Figure 3a,b show the SAR intensity images obtained by TerraSAR-X and TanDEM-X, i.e., the master image and the slave image. It can be seen from Figure 3a,b that there is a lot of speckle noise in the original SAR images. In order to retrieve high-quality ocean waves, noise reduction is needed in the process of ocean wave inversion.

Table 1. SAR image parameters of the TerraSAR-X and TanDEM-X satellites and airborne IRA.

Parameters	Spaceborne InSAR System	Airborne IRA
Center frequency	9.65 GHz	35 GHz
B	301.92 m	0.30 m
Bv	290.06 m	0.299 m
Bp	83.78 m	0.021 m
Azimuth resolution	3.30 m	0.30 m
Ground range resolution	2.26 m	0.30 m
θ	29.83–32.78°	1–15°

Table 1. SAR image parameters of the TerraSAR-X and TanDEM-X satellites and airborne IRA.

Parameters	Spaceborne InSAR System	Airborne IRA
Center frequency	9.65 GHz	35 GHz
B	301.92 m	0.30 m
Bv	290.06 m	0.299 m
Bp	83.78 m	0.021 m
Azimuth resolution	3.30 m	0.30 m
Ground range resolution	2.26 m	0.30 m
θ	29.83–32.78°	1–15°

In addition, the “Guanlan” ocean science satellite project was proposed by the National Laboratory for Marine Science and Technology of China in 2018, the aim of which is to more effectively observe submesoscale ocean phenomena and ocean waves. In the project, IRA is the main payload to be used to retrieve sea surface height (SSH) by means of the interferometric imaging technology. To evaluate the feasibility of the air prototype and the effectiveness of IRA in ocean wave observation, two Ka-band airborne IRA experiments were carried out by the No. 23 Research Institute of the Second Academy of the China Aerospace Science and Industry Corporation (CASIC 23) in the Qingdao XMD sea area in March 2019 and in the Rizhao sea area in November 2020.

Two sets of airborne hybrid InSAR images (labeled as M920 and M3120) used to retrieve ocean waves in this work were provided by CASIC 23, which was carried out in the Qingdao XMD sea area in March 2019, and it has a certain degree of squint. The location of the airborne IRA experimental area can be seen in Figure 2c. The IRA images’ geographical locations in the airborne IRA experiment are given in Figure 2d, and the parameters of the airborne IRA experiment are listed in Table 1. In addition, Figure 3a,b show the SAR intensity images of M920, which are taken as an example.

It should be noted that the InSAR images of M920 and M3120 were applied to another paper on ocean wave inversion based on pure XTI-SAR images with a small incident angle [41]. We found that the squint angle has a certain influence on ocean wave inversion. Therefore, we propose a new ocean wave inversion algorithm for hybrid InSAR images with a certain degree of squint.

3.2. Reference Data

The reference data used in this present work include the fifth-generation reanalysis data ERA5 of ECMWF and the GNSS buoy measurement data. ERA5 is the latest reanalysis data of ECMWF, which can provide global data related to ocean waves. ERA5 can provide a variety of reference SWHs for this work such as the SWH composed of wind waves and swells, SWH of total wind waves, and SWH of total swells. Due to the effects of the orbital velocity and the acceleration of the water waves, only swells can be well retrieved from TerraSAR-X and TanDEM-X images, so the validation data of SWH used in this paper is the significant height of the total swells, which can be represented by SWH_EC in this paper. Without special instructions, hereinafter, the SWH_EC denotes the significant height of total swells provided by ECMWF. Figure 4 shows the SWH_EC at the sea area around Britain at 06:41 on 19 March 2012. The data were obtained by interpolating the SWH_EC at 06:00 and 07:00 on 19 March 2012 obtained from ECMWF. We can intuitively observe the variations in SWH_EC around the sea area from Figure 4. The detailed parameters of the SWH_EC provided by ECMWF are given in Table 2. Moreover, we selected several resolution cells of SWH_EC around the TerraSAR-X and TanDEM-X images for spatial interpolation to obtain the SWH_EC corresponding to each subarea (i.e., A₁, A₂, A₃, and A₄), which can be used as a references for SWH retrieved in the four subareas.

Table 2. The detailed parameters of the SWHEC provided by ECMWF.

Parameters	Parameter Values
Spatial range	Global scale
Spatial resolution	0.5 × 0.5°
Time range	1979–
Temporal resolution	1 h

Table 2. The detailed parameters of the SWHEC provided by ECMWF.

Parameters	Parameter Values
Spatial range	Global scale
Spatial resolution	0.5 × 0.5°
Time range	1979–
Temporal resolution	1 h

In addition, a GNSS buoy was deployed in the Qingdao XMD sea area to measure the ocean wave spectrum during the IRA experiment, which can be used as the reference for IRA measurements. The location of the GNSS buoy can be seen in Figure 2d, and the wind direction (WD) is also shown in Figure 2d, which can be used to determine the propagation direction of the wind waves. The GNSS buoy measurement data can provide the ocean wave parameters of the sea area during the IRA experiment such as the SWH, the significant wave orbital velocity (SWV), and the dominant ocean wave wavelength (λ_p). The in situ sea surface elevation data measured by the GNSS buoy is shown in Figure 4b, and the ocean wave spectra and orbital velocity spectra measured by the GNSS buoy are shown in Figure 4c–f. It can be seen from Figure 4c,d that there are two wave modes in the airborne IRA experiment sea area, and the two wave modes correspond to the wind waves and swells, respectively.

4. Images Preprocessing

Image preprocessing includes images registration, random phase noise suppression, and flat earth effect removal. This paper took TerraSAR-X and TanDEM-X image preprocessing as an example, and IRA image preprocessing can be referred to in [41].

4.1. Images Registration

In order to retrieve high-quality ocean waves, the registration between master and slave SAR images must be performed first. There are many existing image registration methods, and the popular methods used for SAR image registration can refer to the related literatures [44,45,46]. Among the SAR image registration methods, the correlation function method is the most basic statistical method and the basis of many registration algorithms, which is widely used in various types of image registration. In addition, it also has the characteristics of simple operation and strong robustness. At present, the registration algorithms of the only two global interference mapping tasks, SRTM and TanDEM-X interference data processing modules, were all based on the correlation measure function [47,48]. Therefore, the correlation function method was selected for SAR image registration in this paper.

Figure 5a shows the correlation coefficients after performing the registration between the TerraSAR-X image and the TanDEM-X image. Here, we took the correlation coefficients of A₄ as an example. We can see from Figure 5a that the correlation coefficient between the TerraSAR image and the TanDEM image was mainly in the range of 0.4–0.7, which is basically consistent with the correlation coefficient calculated by Romeiser [43]. In addition, Figure 5b shows the interferogram after registration. We can see from Figure 5b that the fringes of the interferogram were very sharp.

4.2. Random Phase Noise Suppression

Due to the existence of random phase noise (as shown in Figure 5b), caused by various decorrelation factors, such as baseline decorrelation, system thermal decorrelation, registration decorrelation, the correlation between the SAR master and slave images were reduced, and the accuracy of the retrieved sea surface profile and ocean wave velocity were also seriously reduced. Therefore, it was necessary to suppress the influence of phase noise first when we retrieved ocean waves, and we used multilook processing to reduce the effect of phase noise in this work. The detailed process of random phase noise suppression can be referred to in the literature [41].

In order to discuss the effect of random phase noise, we calculated the theoretical and the measured root mean square error (RMSE) of the random phase noise, which are denoted by σ_TN and σ_SAR, using Equations (17) and (18), respectively.

$${\sigma }_{\mathrm{TN}}=\sqrt{{\displaystyle \underset{-\pi +\overline{\psi }}{\overset{\pi +\overline{\psi }}{\int }}{(\psi -\overline{\psi })}^2}{p}_{\psi }^{NL}(\psi )\mathrm{d}\psi }$$

(17)

$${\sigma }_{\mathrm{SAR}}=\sqrt{{\displaystyle \underset{-\pi +\overline{\psi }}{\overset{\pi +\overline{\psi }}{\int }}{({\psi }_{\mathrm{SAR}}-\overline{\psi })}^2}{p}^{NL}({\psi }_{\mathrm{SAR}})\mathrm{d}{\psi }_{\mathrm{SAR}}}$$

(18)

where $\overline{\psi }$ denotes the mean value of the random phase, $\psi$; ${\psi }_{\mathrm{SAR}}$ is the phase measured by SAR; $p^{NL}(\psi )$ is the probability density distribution of the random phase, $\psi$; the expression of $p^{NL}(\psi )$ can be referred to the in literature [41].

The variations in σ_SAR with the looknumber (NL) are shown as the blue lines in Figure 6. In addition, the σ_TN is also shown in Figure 6 with red lines as a contrast. To remove the flat earth effect along the range direction on the interference phase, here, the spaceborne interferogram was only smoothed along the extension direction of the interference fringes. We can find from Figure 6 that with the increase in the NL, the decrease rate of the σ_SAR value was significantly slower than that of σ_TN, which was just as expected. In addition, the relative difference between σ_SAR and σ_TN increased with the NL, which means that the effect of the phase noise on the spaceborne interferogram became weaker with the increase in the NL. Moreover, we can also use Equation (19) to evaluate the relative error caused by the phase noise when we retrieved the ocean waves.

$$\mathrm{Err}=\frac{{\sigma }_{\mathrm{SAR}}-\sqrt{{\sigma }_{\mathrm{SAR}}^2-{\sigma }_{\mathrm{TN}}^2}}{\sqrt{{\sigma }_{\mathrm{SAR}}^2-{\sigma }_{\mathrm{TN}}^2}}$$

(19)

Two columns of the interference phase were selected along the extension direction of interferogram, which were used to evaluate the relative errors with Equation (19) for different NLs, and the relative errors are given in

Table 3. The relative error (%) of the interference phase caused by phase noise with different looknumbers in the 700th and 2200th columns, respectively.

Column Number	NL = 10	NL = 20	NL = 40	NL = 60	NL = 80	NL = 100	NL = 120	NL = 140	NL = 160
700th	25.76	13.45	9.15	7.81	6.90	6.16	5.77	5.47	5.18
2200th	86.74	32.25	11.62	8.84	7.48	6.58	6.07	5.74	5.27

. In order to maintain the relative error below 10% in the present work, we chose NL = 6 × 10 (NL = 6 × 10 means that there were 6 pixels and 10 pixels along azimuth and range directions, respectively) for random phase noise suppression.

4.3. The Flat Earth Effect Removal

In order to retrieve ocean waves, the flat earth phase, which will increase the difficulty in inverting the sea surface profile and sea surface velocity, needs to be removed. We can use Equation (20) to calculate the flat earth phase, which is based on the spaceborne platforms’ positions and the ground geographical positions of the TerraSAR-X and TanDEM-X images.

$${\phi }_f=-\frac{4\pi }{\lambda }\left(\left|P_2-T\right|-\left|P_1-T\right|\right)$$

(20)

where φ_f denotes the flat earth phase; P₁ and P₂ are the antenna locations of TerraSAR-X and TanDEM-X; T is the geographical location of each resolution cell of the TerraSAR-X and TanDEM-X images.

Figure 7a,b show the interferogram after multilook processing and the flat earth phase, and we also selected A₄ as an example. By removing the flat earth phase from the interferogram, the interference phase caused by sea surface elevation and sea surface velocity was obtained as shown in Figure 7c. To clearly display the interference phase caused by the sea surface profile and the sea surface velocity, the subregion marked as a blue box in Figure 7c was selected for display as shown in Figure 7d.

5. Ocean Wave Inversion and Discussion

5.1. Sea Surface Profile Inversion

On the basis of the interference phase caused by the sea surface profile and sea surface velocity, the sea surface profile can be retrieved using Equation (13). To clearly display the ocean wave inversion results, we selected a subarea from A₁, A₂, A₃, and A₄ separately for display, as shown in Figure 8a–f, of the sea surface profiles retrieved using the airborne hybrid InSAR images, which corresponded to M920 and M3120, respectively. According to Figure 8, we can see that the sea surface profiles of ocean waves could be clearly retrieved using the spaceborne hybrid InSAR images and airborne hybrid InSAR images. It can also be seen from Figure 8 that the dominant wavelengths in A₃ and A₄ were longer than those in A₁ and A₂. The SWHs (4σ_h, σ_h is the root mean square height of the ocean waves) of the retrieved ocean waves in A₁, A₂, A₃, and A₄ were 0.79, 1.05, 0.64, and 0.67 m, respectively, which are listed in

Table 4. The ocean wave parameters retrieved from the spaceborne hybrid InSAR images and the SWHs provided by ECMWF.

Areas	A1	A2	A3	A4
λp	210 m	209 m	292 m	305 m
Direction	100°	100°	230°	230°
SWV	0.57 m/s	0.80 m/s	0.56 m/s	0.58 m/s
SWVSP	0.46 m/s	0.61 m/s	0.13 m/s	0.15 m/s
SWH	0.79 m	1.05 m	0.64 m	0.67 m
SWHSP	0.75 m	0.92 m	0.32 m	0.34 m
SWHEC	0.74 m	1.09 m	0.87 m	0.71 m

. The SWHs of the retrieved ocean waves shown in Figure 8e,f were 0.37 and 0.38 m, respectively, which are listed in

Table 5. The ocean wave parameters retrieved from the airborne hybrid InSAR images and m-measured by the GNSS buoy.

Parameters		M920	M3120
SWV	IRA	0.31 m/s	0.32 m/s
	Buoy	0.64 m/s	0.68 m/s
SWVSP	IRA	0.23 m/s	0.23 m/s
	Buoy	0.64 m/s	0.68 m/s
SWH	IRA	0.37 m	0.38 m
	Buoy	0.33 m	0.38 m
SWHSP	IRA	0.28 m	0.29 m
	Buoy	0.33 m	0.38 m
SWHC	IRA	0.39 m	0.39 m
SWHSPC	IRA	0.31 m	0.32 m
Direction	IRA	270° (Swell)	270° (Swell)
	Buoy
λp	IRA	87.27 m (Swell)	89.76 m (Swell)
	Buoy	70.69 m (Swell) 12.85 m (Wind waves)	100.98 m (Swell) 17.67 m (Wind waves)

SWH_C and SWH_SPC denote the SWHs retrieved without considering the influence of squint angle.

. On the basis of the SWH_EC provided by the ECMWF around the SAR images, the SWH_EC can linearly be interpolated to acquire the SWH_EC at the center of A₁, A₂, A₃, and A₄ as a reference, and the interpolated SWH_ECs in A₁, A₂, A₃, and A₄ were 0.74 m, 1.09 m, 0.87, and 0.71 m, respectively. The differences among the SWHs retrieved using the spaceborne hybrid InSAR images and SWH_ECs in A₁, A₂, A₃, and A₄ were 0.05, −0.04, −0.23, and −0.04 m, respectively. Before comparing the retrieved SWHs using airborne hybrid InSAR images with the SWHs measured by the GNSS buoy, it was necessary to use a low-pass filter to process ocean wave spectra measured by the GNSS buoy in order to make the ocean wave spectra measured by the GNSS buoy comparable to the ocean wave spectra retrieved using the airborne hybrid InSAR images. According to the wind direction, as shown in Figure 2d, and the resolutions of the airborne IRA images, the time window (T_LF) of the low-pass filter can be calculated by [41]:

$$T_{LF}=\sqrt{\frac{2\pi L}{g}}$$

(21)

and

$$L=L_{xf}\cos \gamma +L_{yf}\sin \gamma$$

(22)

where L_xf and L_yf are the range and azimuth scales of the two-dimensional filter used in [41], respectively, and γ is the angle between the wind direction and the range direction of the IRA image. Therefore, the calculated L and T_LF were approximately 12.21 m and 2.80 s, respectively, and the low-pass filter mainly suppressed wind waves with a wavelength of approximately 12 m and a frequency of approximately 0.36 Hz. After processing the ocean wave spectra measured by the GNSS buoy, the SWHs measured by the GNSS buoy were 0.33 and 0.38 m, respectively, and the differences between the SWHs retrieved using the airborne hybrid InSAR images and the SWHs measured by the GNSS buoy were 0.04 m and 0, respectively.

Moreover, we can also see from Figure 8 that the ocean waves in A₁ and A₂ propagated almost along the range direction, and the ocean waves in A₃ and A₄ propagated approximately along the direction of 45° or 225° with respect to the along-track direction. In order to remove the 180° directional ambiguity from the ocean wave spectra, we retrieved the cross-spectra based on the TerraSAR-X intensity images. The inversion of the cross-spectra was based on the split-look processing technique, which divides the azimuth reference signal into multiple subsignals with different center times. Therefore, split-look processing can generate a series of time series subimages, and then the cross-spectra can be obtained by performing a correlation calculation on the subimages. For the specific method of retrieving cross-spectra, please refer to the literature [49]. The cross-spectra of ocean waves in the four subareas are shown in Figure 9, in which Figure 9a–d display the real parts, and Figure 9e–h show the imaginary parts. Here, k_r represents the components of the wavenumber in the range direction, and k_a represents that in the azimuth direction. The angle of the polar coordinates is in degrees clockwise from true north. It can be seen from Figure 9e–h that the imaginary parts of the TerraSAR-X image cross-spectra had two spectral peaks, one of which was positive and the other which was negative. If the ocean wave spectra were unimodal, the ocean wave spectra were located in those spectral domains, where the imaginary part of the cross-spectra were positive [50].

According to the sea surface profiles in A₁, A₂, A₃, and A₄, and combined with the imaginary parts of the SAR images’ cross-spectra, the corresponding two-dimensional ocean wave spectra could be acquired as shown in Figure 10a–d. Based on the corresponding two-dimensional ocean wave spectra shown in Figure 10a–d, the wavelengths of the dominant waves corresponding to the spectral peaks in A₁, A₂, A₃, and A₄ could be calculated, which were approximately 210, 209, 292, and 305 m, respectively. The wavelengths of the dominant waves in A₁ and A₂ were consistent with those retrieved by Romeiser [43]. In addition, the propagation directions of the dominant waves in A₁, A₂, A₃, and A₄ are approximately 100°, 100°, 230°, and 230° clockwise versus true north, respectively. For ocean waves retrieved using the airborne hybrid InSAR images, the two-dimensional ocean wave spectra of M920 and M3120 are shown in Figure 10e,f. In order to better compare the retrieved ocean wave spectra with that measured by the GNSS buoy, the retrieved two-dimensional ocean wave spectra were converted to one-dimensional ocean wave spectra, which are shown in Figure 10g,h. It can be seen from Figure 10g,h that the retrieved swell wavelengths of M920 and M3120 were approximately 87.27 and 89.76 m, respectively; the measured swell wavelengths of M920 and M3120 were approximately 70.69 and 100.98 m, respectively; the swell wavelength differences of M920 and M3120 between the retrieved and measured were 16.58 and −11.22 m, respectively. However, due to the wind waves that propagated along the azimuth direction of M920 and M3120, the velocity bunching had a strong nonlinear influence on the wind waves imaging, and the wind waves could not be retrieved well. For the ocean wave propagation directions retrieved using the airborne hybrid InSAR images, the swells propagation directions of M920 and M3120 could be determined by combining with the geographical locations of the experimental areas, and the swells propagation directions of M920 and M3120, shown in Figure 10e,f, were all approximately 270° counterclockwise versus true north, respectively.

It is worth mentioning that there were differences between the SWHs calculated based on the ocean wave spectra (as shown in Figure 10) and those evaluated directly from the sea surface profiles in Figure 8. Here, the retrieved SWHs based on the ocean wave spectra are marked as SWH_SP and listed in Table 4 and Table 5. It can be seen from Table 4 that the values of the SWH_SP were significantly smaller than the SWHs, especially in sea areas A₃ and A₄. This was because the nonlinear influence of the velocity bunching in these two areas was more significant than that in A₁ and A₂. For the airborne hybrid InSAR system, although the velocity bunching also had a significant nonlinear influence on M920 and M3120, the IRA images of M920 and M3120 were acquired at low sea state conditions; therefore, the differences between the SWH_SP and those measured by the GNSS buoy were not as significant as that in A₃ and A₄. For the InSAR systems, the larger the angle between the ocean wave propagation direction and the range direction, the more significant the nonlinear influence of the velocity bunching [28]. We can see from Figure 8 and Figure 10 that the dominant waves in A₁ and A₂ almost propagated along the range direction, while the angles between the dominant waves’ propagation direction and the range direction in A₃ and A₄ were approximately 50°. Therefore, the nonlinear influence of the velocity bunching in A₃ and A₄ was more significant. During the process of ocean wave imaging using the SAR systems, the velocity bunching not only caused nonlinear shifting of the ocean waves but also the superposition of ocean waves at different positions. In addition, the imaging results caused by velocity bunching were different between the SAR systems and the InSAR systems. For the SAR systems, velocity bunching caused a position shift in the sea surface scattering elements along the azimuth of the SAR image, and the echo intensity produced convergence and divergence effects along the azimuth of the radar. The scattering intensities in the convergence region increased and that in the divergence region decreased. However, for the InSAR systems, the ocean waves’ shift caused by velocity bunching caused the superposition of the partial scattering fields, and the corresponding interference phase also changed indeterminately, and based on this analysis, the retrieved SWH based on the InSAR images may be less than the true value of the SWH. In addition, the nonlinear shifting of the ocean waves was caused by velocity bunching because the high-frequency components were mixed with the random noise. Therefore, the speckles in Figure 8 were not all random noise, and some of them were the drift ocean waves. When retrieving the two-dimensional ocean wave spectra, we used a low-pass filter to filter out the high-frequency components and retain only the low-frequency components. The final result is shown in Figure 10. Thus, the retrieved significant wave height based on the InSAR images should be less than the true value of the significant wave height, and the stronger the influence of the velocity bunching, the greater the difference between the retrieved value and the true value.

5.2. The Orbital Velocity of Ocean Waves Inversion

On the basis of the interference phase caused by sea surface elevation and sea surface velocity, the orbital velocity of the ocean waves can be acquired by using Equation (15). To clearly display the orbital velocities of ocean wave inversion results, we selected a subarea from A₁, A₂, A₃, and A₄ separately for display as shown in Figure 11a–d. Figure 11e,f are the orbital velocities of the ocean waves retrieved using the airborne hybrid InSAR images, which corresponded to M920 and M3120, respectively. According to the retrieved orbital velocities of the ocean waves and combining with the imaginary parts of the TerraSAR-X image cross-spectra, the corresponding two-dimensional velocity spectra can be obtained, which are shown in Figure 12a–d. Figure 12e,f show the two-dimensional velocity spectra retrieved using the airborne hybrid InSAR images, and Figure 12g,h show the corresponding one-dimensional velocity spectra. We can see from Figure 11 and Figure 12 that the orbital velocity of the ocean waves in A₂ were larger than those in A₁, A₃, and A₄. The significant wave orbital velocity (SWV) calculated based on the velocities (as shown in Figure 11) (SWV = 4σ_v, σ_v is the root mean square orbital velocity of ocean waves) and calculated based on two-dimensional velocities spectra (as shown in Figure 12a–f) (SWV_SP) are also listed in Table 4 and Table 5. We can find from Table 4 and Table 5 that the values of SWV_SP were significantly smaller than SWVs, especially in A₃, A₄, M920, and M3120. The reason for the differences is the same as that for SWH. However, since the ECMWF does not directly provide the radial velocity of ocean waves, there is no available reference orbital velocities of ocean waves to compare with the results retrieved in this paper. For the orbital velocity of ocean waves retrieved using the airborne hybrid InSAR images, the SWVs and the SWV_SP were significantly smaller than that measured by the GNSS buoy. From the significant differences, as shown in Figure 12g,h, it can be deduced that the SWVs measured by the GNSS buoy included the SWVs of swells and wind waves, and when using IRA images to retrieve the orbital velocities of ocean waves, due to the strong nonlinear influence of velocity bunching, only the orbital velocity of swells could be retrieved.

6. Conclusions

For the hybrid interferometry systems, which are a combination of along-track and cross-track configurations, we presented a new method for retrieving ocean wave information such as the sea surface profile and the orbital velocity of ocean waves. This method was derived on the basis of the linear modulation theory in which the sea surface profile can be represented as the superposition of a series of traveling wave components. On the basis of the above research, the following conclusions can be drawn:

(1)

On the basis of the spaceborne hybrid InSAR images acquired by the TerraSAR-X and TanDEM-X and the hybrid InSAR images acquired by airborne IRA, ocean waves, ocean waves spectra, orbital velocities of ocean waves, and the velocity spectra were retrieved using a new method. For the spaceborne hybrid InSAR images, since the effects of the orbital velocity and the acceleration of the water waves, only swell were well retrieved from the satellite SAR images in the present work. For the airborne hybrid InSAR images used in this paper, due to the nonlinear influence of the velocity bunching along the azimuth direction, only swells could be retrieved well;

(2)

Based on the ocean wave spectra and velocity spectra retrieved using this method, we calculated SWHs, peak wavelengths of dominant waves, and SWVs. In addition, we also obtained the propagation direction of the dominant waves in A₁, A₂, A₃, and A₄ on the basis of the retrieved cross-spectra. The wavelengths of dominant waves in A₁ and A₂ were consistent with those mentioned by Romeiser [43], and the differences between the retrieved SWHs in A₁ and A₂ and that provided by ECMWF were small, while the differences between the retrieved SWHs in A₃ and A₄ and that provided by ECMWF were large, which was mainly caused by the strong nonlinear influence of velocity bunching in A₃ and A₄. The differences between the retrieved SWHs using the M920 and M3120 images and that measured by the GNSS buoy were small, while the retrieved SWVs using the M920 and M3120 images were significantly smaller than that measured by the GNSS buoy. The significant differences were because the SWVs measured by the GNSS buoy included the SWVs of swells and wind waves, and for the retrieved SWVs, due to the strong nonlinear influence of velocity bunching, only the orbital velocity of the swells could be retrieved;

(3)

There was a close relationship between the accuracy of the retrieved ocean wave spectrum, orbital velocity spectrum, and velocity bunching. When the nonlinear influence of the velocity bunching was strong, the accuracy of the retrieved ocean wave spectrum and orbital velocity spectrum was low, and when the nonlinear influence of the velocity bunching was weak, the accuracy of the retrieved ocean wave spectrum and orbital velocity spectrum was high. Therefore, addressing the nonlinear influence of velocity bunching on ocean wave inversion is still a challenge in future work. Moreover, since the airborne hybrid interferometric images used in this paper all had strong velocity bunching nonlinearity, the accurate evaluation of ocean wave inversion based on the airborne hybrid interference images with small velocity bunching nonlinearity is also an important task in the future.