The frequency shift method removes the flat earth phase based solely on the interferogram and it does not require any additional parameters. When no flat earth effect is present, the spatial frequency of the phase fringe (i.e., the frequency spectrum peak value of the interferogram) should be zero (Ai and Li, 2009; Gatelli et al., 1994). Therefore, after performing a Fourier transform of the interferogram, the flat earth phase can be removed by searching for the frequency peak value and shifting it to zero using a circular shifting method. The two-dimensional interferogram needs to remove the flat earth phase along both the row and column directions.

Assume the master and slave images after registration are $M(m, n)$ and $S(m, n)$, respectively. The interferogram is $I(m, n) = $$M(m, n) \cdot {S^ * }(m, n) $, and the two-dimensional interferometric phase is $\varphi (m, n) = \arg (I(m, n))$. Here, the asterisk * indicates the conjugate complex, $\arg (\;)$ calculate the argument, and m and n are the row and column coordinates in the image coordinate system, respectively. Generally, the interferogram phase changes periodically along the range and azimuth directions. Therefore, processing along both the row and column directions may be necessary. Because both directions involve the same processing steps, row processing, which is the range direction in the data used in this paper, will be described next as an example.

(1) Calculating the frequency spectrum peak value of the phase fringe

First, the two-dimensional phase is transformed to the frequency domain by a one-dimensional fast Fourier transform (FFT) row by row. The values are then summed column by column, resulting in a one-dimensional array that corresponds to the rows. The position of the maximum value of this one-dimensional array corresponds to the approximate frequency spectrum peak value $\,{f_1}$ in the row direction.

To obtain a more accurate frequency, a high-order smooth interpolation such as the cubic spline method is applied to the 5 points centered on the maximum value of the one-dimensional array. Then, the 5-point interval is interpolated into 400 grid points here. The position of the maximum value of interpolated results corresponds to the exact frequency spectrum peak value $\,{f_2}$. The phase corresponding to $\,{f_2}$ is the flat earth phase, which can be described as follows:

where n indicates the column number of a resolution unit within a row, and j is the imaginary unit.

(2) Removing the flat earth phase

To remove the flat earth effect, the interferogram is multiplied by the phase shift factor corresponding to ${f_2}$:

where $\varphi '(m, n)$ indicates the interferometric phase after removing flat earth phase.

When two apertures S₁ and S₂ successively image the target point P, the phase difference ${\phi _f}$ corresponding to the distance difference $\Delta R$ is the flat earth phase, which can be described as follows:

where ${S_2}P$ indicates the distance between the aperture S₂ and the target point P, and ${S_1}P$ is the distance between the aperture S₁ and P. Based on the orbital parameters from the head file and the geographic coordinates of reference points from a geographical reference file, the geocentric coordinates of S₁, S₂ and P can be obtained as described below (Cao et al., 2013).

(1) The aperture position $S(x, y, z)$ when the target is being imaged

First, based on its row and column number (${m_P}$, ${n_P}$) in the original satellite data, the imaging time of P point is given by

where ${t_s}$ indicates the imaging time of the first resolution unit, and PRF and RSF indicate the pulse repetition frequency and the range sampling frequency, respectively.

Second, the orbit equation, $x(t)$, $y(t)$ and $z(t)$, is fitted based on the orbit information from the head file.

where x, y and z are the orbit coordinates, a_i, b_i and c_i (i = 0, 1, 2, 3) are the fitting coefficients, and t is the time variable. Equation (5) describes the change of the orbit coordinates with time. Then, the aperture position $S(x, y, z)$ as the point P is being imaged can be obtained by substituting the imaging time $({t_P})$ of point P into Eq. (5). Note that the positions of the two satellites in formation need to be calculated, respectively.

(2) The target position $P({x_P}, {y_P}, {z_P})$ in the geocentric coordinate system

First, the latitude (lat) and longitude (lon) are fitted as a function of the row number (m) and column number (n) based on the location information from the geographical reference file.

where ${\alpha _i}$ and ${\beta _i}$ (i = 0, 1, 2, 3, 4, 5) are the fitting coefficients. Then, the latitude and longitude of point P can be obtained by substituting the row and column number (${m_P}$, ${n_P}$) of point P into Eq. (6).

Second, the latitude and longitude of point P are transformed to geocentric coordinates using the following formula:

where a, b and e indicate the long axis, short axis and eccentricity of the WGS84 ellipsoid, respectively.

The original interferogram from the ATI complex data introduced in Section 2 is shown in Fig. 2. Note that the data are presented in image coordinates, the radar look direction (i.e., the range direction) is exactly from left to right in the figures. The phase changes periodically during the interval (–π, π] along the range and azimuth directions; this phase change is caused primarily by the flat earth effect. The flat earth phase obscures the ocean surface current information; thus it should be removed. Here, the conventional frequency shift method and the orbital parameters method are used to remove the flat earth effect, respectively. The result is shown in Fig. 3, where it can be found that, although there is phase wrapping (mainly caused by the ocean surface current), the periodic fringes have been removed by the two methods, indicating that both methods have removed the flat earth phase. Nevertheless, the methods’ performances are quite different; the orbital parameters method accurately and effectively removes the flat earth phase based on the precise orbit parameters of the TerraSAR-X satellite (Werninghaus and Buckreuss, 2010; Yoon et al., 2009); however, a large amount of residual flat earth phase remains in the frequency shift method result. In the ideal case (i.e., a zero across-track baseline), ATI-SAR detects only the ocean surface velocity field. Because the movement speed is zero over the land or island regions, the absolute value of the calibrated interferometric phase should also be zero. When no phase calibration is conducted, the phase relative value after flat earth phase removal generally has an offset from zero, but it tends to be homogeneous (i.e., the phase variation should be small). In contrast, in the actual results shown in Fig. 3a, the phase variation across different land masses and islands is small, indicating that the orbital parameters method accurately removed the flat earth phase. However, the phase variation is significant across the various islands in Fig. 3b, indicating that a large amount of residual flat earth phase remains in the result after applying the frequency shift method.

To present our argument more clearly, the unwrapped phase calibrated by averaging the land region phase to zero, which is the absolute value of the interference phase, is shown in Fig. 4. In Fig. 4a, the orbital parameters method causes the phase over land masses and islands to all tend toward zero; and the little remaining phase variations can be attributed primarily to scattered noise and topological fluctuations. However, in the frequency shift method results shown in Fig. 4b, the phases of several islands deviate significantly from zero. For example, the phase of the land at the lower left is less than zero, while that of the upper right island is greater than zero. Overall, the orbital parameters method accurately removes the flat earth phase, while the frequency shift method leaves a large amount of residual flat earth phase. In the following section, an improved method will be proposed for the frequency shift method that removes the residual flat earth phase.

The improved frequency shift method is conducted using the following steps.

(1) Divide the phase diagram into sub-blocks to calculate the frequency spectrum peak value of phase fringes.

In an actual ATI interferogram shown in Fig. 2, the spatial frequency of the phase fringes is generally not a constant; instead, it varies slowly along the range direction. In the conventional frequency shift method described above, the frequency spectrum peak value is the average spatial frequency of the phase fringe over the entire image and does not consider variations in spatial frequency with position. This omission results in obvious errors of flat earth phase removal at both ends of the range direction. To account for the variations of fringe spatial frequency, we divide the interferogram into several sub-blocks to calculate the spatial frequency of the phase fringes. Using this approach, we can not only improve the precision but also obtain its spatial variation.

First, we divide the interferogram into several sub-blocks along the range direction and calculate the frequency spectrum peak value f of each block using the first step of the conventional frequency shift method. The number of sub-blocks is governed by two contradictory factors. On one hand, in principle, the larger the number of sub-blocks is, the more distinct the spatial frequency variations are. On the other hand, if the number of sub-blocks is too large, the length of each sub-block in the range direction will be too small, which makes calculating the spatial frequency accurately a challenge. Although the periodical variations from –π to π of the original interference phase are mainly caused by the flat earth effect, as shown in Fig. 2, ocean surface currents and wave motion also affects the spatial distribution of the phase fringes—the distribution of phase fringe is especially distorted in regions with large velocity gradients. These distortions directly affect the spatial frequency estimation of the phase fringes. The smaller the sub-block size is, the more serious the influence of phase fringe distortion is. To suppress this effect, the length of each phase diagram sub-block in the range direction should be larger than several spatial periods, and abnormal spatial frequency values should be discarded. Based on the above principles, while also considering that the following second-order polynomial fitting requires a minimum of 3 sub-blocks, we divide the phase diagram into five sub-blocks.

Second, we fit the second-order functional relationship between the spatial frequency f and the position n of the center point of each sub-block in the range direction.

where a_i (i=0, 1, 2) are the fitting coefficients.

Finally, we calculate the frequency spectrum peak value in each grid point using Eq. (8) and then remove the flat earth effect through the second step of the conventional frequency shift method.

The phase after removing flat earth effect by dividing the phase diagram into sub-blocks is shown in Fig. 5. Figure 5a shows the relative phase after removing the flat earth effect, while Fig. 5b shows the relative phase after further unwrapping. It is found that, compared with the conventional frequency shift method, the initial improvement to the method from dividing the phase diagram into sub-blocks significantly suppresses the phase spatial variation in the land regions, indicating that the improvement significantly reduces the residual flat earth phase. However, land region phase still varies to a certain extent among the various spatial positions. This remaining variation can be attributed to two factors if random phase noise is not considered. One factor is residual flat earth phase, which causes the phases of different islands over large distances to vary with the spatial position in both the range and azimuth directions (e.g., the phase of the island in the upper right of Fig. 6a is greater than zero, while the land region in the lower left is less than zero). Therefore, further removal of the residual flat earth phase is needed. The other factor is the effect of topography. In the actual TanDEM-X satellites formation, both AT baseline and XT baseline (24 m in our data) occur; thus, the interferometric phase will be affected by topography. To observe the effect of topography, the phases of the land and island regions from Fig. 5b are shown in Fig. 6. Note that the phase of the island in the upper right corner of Fig. 6a has obvious spatial variations caused primarily by topography rather than by residual flat earth phase because areas as small as this island should result in a homogeneous residual flat earth phase. Thus, topography should also be considered to further remove residual flat earth phase.

(2) Correct the residual flat earth phase by fitting the function between the phase and position of land areas.

Regardless of random phase noise, the spatial variation of the residual flat earth phase should be slow and regular. Therefore, by fitting the function of the phase of the land region to the location coordinates, the residual flat earth phase can be further removed. Note that the effect of topography fluctuations also needs to be avoided or suppressed. Because the wrapped phase may jump from π to –π, which is not convenient for processing, subsequent processing operates on the unwrapped phase (Fig. 5b).

First, we remove the high elevation region of the island, where the phase is seriously affected by the topography. Then, the phase variation in the remaining land region can be considered to be caused only by residual flat earth phase. In Fig. 6a, the phase of the island in the upper right corner is seriously affected by the island’s topography (the elevation varies from 1 m to approximately 50 m and the corresponding phase varies from approximately 1.0π to approximately 1.6π). The elevated areas above 25 m are removed using an ATI phase threshold, and the results are shown in Fig. 6b. The ATI phase threshold is selected not only to suppress the topographical impact but also to ensure that a sufficiently large area remains for phase averaging to remove the random phase noise.

Second, the three island regions (indicated by the black boxes in Fig. 6b) are selected to calibrate the phase. Because each selected area is small, the residual flat earth of each area is approximately constant. Therefore, by averaging the phase of each region and fitting the function between the average phase and the center position of each region, the residual flat earth phase can be further removed. Note that the average phase contains two parts. One is the overall offset of the phase relative value, which does not vary with the position. The other is the residual flat earth phase, which slowly varies with position. A binary linear equation is used to fit the function between the phase correction and position due to the slow change of the residual flat earth phase.

where b_i (i = 0, 1, 2) are the fitting coefficients, and m and n indicate the row and column number, respectively.

Finally, we calculate the phase correction of every grid point using Eq. (9) and subtract the results from the results of the preliminarily improved method. The calibrated phase is shown in Fig. 7.

The calibrated phase obtained through the improved frequency shift method is shown in Fig. 7. As shown, in the land and island regions, the phase tends to be zero with little fluctuation, which indicates that the improved frequency shift method accurately removes the flat earth effect. To further verify the effect of the improved method, the results of the conventional and the improved frequency shift methods are statistically compared with those of the orbital parameters method. Under conditions where accurate orbital parameters and geographic information are provided by TanDEM-X L1B data, the orbital parameters method is considered to be reasonably accurate; in other words, its errors when removing the flat earth phase can be ignored. Scatter diagrams and statistical analysis results are shown in Fig. 8. We computed the shown correlation and regression coefficients, the mean differences, and the residual root mean square (RMS) differences (i.e., the RMS differences after subtracting the mean differences). The results show that a large difference exists between the conventional frequency shift method and the orbital parameters method; this difference occurs because the conventional phase-shift method still retains a large number of unevenly distributed residual flat earth phase. In contrast, the improved frequency shift method agrees well with the orbital parameters method; the correlation coefficient and the regression coefficient are 0.984 and 1.031, respectively, which are great improvements from 0.824 and 0.531 values respectively achieved by the conventional frequency shift method. The accuracy is also greatly improved, the average difference decreases from –0.22π to –0.03π, and the root mean square error decreases from 0.44π to 0.11π. These reductions indicate that the accuracy of the improved method is similar to that of the orbital parameters method.

In view of the above analysis, the orbital parameters method may be the first choice for ATI-SAR flat earth phase removal. When sufficiently accurate orbital parameters are available, the orbital parameters method can accurately remove the flat earth phase. However, for some satellites, the orbital parameters are relatively inaccurate causing the orbital parameters method to have low accuracy when removing the flat-earth phase. In particular, when the orbital parameters are unavailable, the method cannot be used. In such cases, the improved frequency shift method proposed in this paper is an alternative for removing the flat-earth phase. The improved frequency shift method achieves approximately the same accuracy for flat earth phase removal as the orbital parameters method when the accurate parameters are available; however, the improved frequency shift method is based solely on the interferogram and does not require any additional parameters.