*2.4. Aliasing and Floor Error*

Since the visibility equation is fundamentally a Fourier transform, discretization grids for regular sampling in both domains must be reciprocal to each other. The lattice depends on the overall geometry: Rectangular for U-shaped instruments [14], or hexagonal [15] for Y-shaped (MIRAS), hexagonal [16] or triangular ones. In any case, the visibility sampling coordinates (*u*, *v*) are a subset of the grid points in the fundamental period (a square for rectangular grids and a star or an hexagon for hexagonal grids). The minimum number of grid points in the fundamental period needed to include all measured samples is *N*<sup>2</sup> *<sup>T</sup>* where *NT* = 4*N*EL + 1 for a hexagonal or triangular instrument, *NT* = 3*N*EL + 1 for a Y-shaped instrument or *NT* = 2*N*EL + 1 for a rectangular instrument [8]. In all cases, *N*EL is the number of elements in one arm. Since MIRAS has *N*EL = 21, it follows that *NT* = 64 in this case (The SMOS Level 1 Operational Processor uses *NT* = 196). The number of points in the fundamental period of the corresponding (*ξ*, *η*) reciprocal grid is also *N*<sup>2</sup> *<sup>T</sup>*. If reciprocal grids are used, the elementary area Δ*ξ*Δ*η* in Equation (5) becomes then equal to 1/(*N*<sup>2</sup> *<sup>T</sup>d*2) for rectangular grids or 1/(*N*<sup>2</sup> *<sup>T</sup>d*<sup>2</sup> sin 60◦) for hexagonal grids [8] where *d* is the minimum antenna spacing normalized to the center wavelength.

Figure 2 shows the MIRAS reciprocal grids for *NT* = 64. The fundamental hexagon is depicted in blue in both domains, the green star in (*u*, *v*) includes all measured visibility points and their conjugate ones; and the extension of the (*ξ*, *η*) grid to the unit circle is drawn in gray. Both fundamental hexagons have the same number of points, in this case equal to 64<sup>2</sup> = 4096. A larger value just enlarges the (*u*, *v*) hexagon while keeping the star shape identical, and provides a thicker grid in (*ξ*, *η*) [8]. The fundamental hexagon in this domain has a fixed side, independent of the number of points, equal to 2/(3*d*) where *d* is defined in the previous paragraph. This hexagon always falls inside the unit circle if *<sup>d</sup>* <sup>&</sup>gt; 1/√3.

**Figure 2.** MIRAS reciprocal grids for visibility (**left**) and brightness temperature (**right**) using *NT* = 64. Fundamental hexagons in both domains are drawn in blue. Two areas with hermitic (*u*, *v*) points are highlighted in the left plot. Unit circles aliases are added in the right plot.

A Fourier transform needs zero padding to complete the (*u*, *v*) fundamental period. The resultant modified brightness temperature is then obtained in the (*ξ*, *η*) fundamental period. The periodic repetition of the phase produces the well known phenomenon of aliasing, resulting in that the same image is repeated at all adjacent periods. All unit circle aliases for the MIRAS case are depicted in Figure 2. The zone in which these circles do not overlap is the alias-free field of view. Out of it, the image is always contaminated with replicas of other areas of the same image. In principle, error-free imaging is only possible in the alias-free field of view unless the overlapped image content is null, as for example in the case of having a small target in the center of the field of view surrounded by a very low background. Contrarily, for wide field of view imagers, aliasing is a strong source of errors. Changing the normalized antenna spacing *d* to a value lower than the Nyquist sampling rate (1/2 for rectangular grids and 1/ <sup>√</sup><sup>3</sup> for hexagonal grids) would make the circle be inscribed within the fundamental period, so eliminating the aliases.

For the minimum MIRAS reciprocal grid of Figure 2, the total number of (*ξ*, *η*) points within the unit circle is *Np* = 8491, so this is the number of columns of the G-matrix in this case. Splitting the columns into the fundamental hexagon *GH* and the rest of the unit circle *GNH* (see Figure 2), the discretized visibility equation can be written as

$$V = \begin{bmatrix} \ G\_H & G\_{NH} \end{bmatrix} \begin{bmatrix} & T\_H \\ & T\_{NH} \end{bmatrix} = G\_H T\_H + G\_{NH} T\_{NH} \tag{8}$$

where the same nomenclature applies to *TH* and *TNH*. Clearly, inverting only the G-matrix in the fundamental hexagon *T* = *G*−<sup>1</sup> *<sup>H</sup> V* neglects the term *GNHTNH* and produces what is sometimes called "floor error" [17,18]. Contrarily to the case of aliases in Fourier inversion, this one also spreads into the alias-free field of view unless considering identical antenna patterns and no fringe washing function. In this limiting case, imaging with G-matrix is equivalent to an inverse Fourier transform and the floor error is reduced to the aforementioned aliasing error.

In any case, the floor error can be mitigated by subtracting from the visibilities an estimation based on a forward a priori brightness temperature model outside the hexagon *MNH*. Using Equation (8), the equation to invert becomes

$$V - G\_{\rm NH}M\_{\rm NH} \approx G\_{\rm H}T\_{\rm H} \tag{9}$$

which leads to

$$T\_H \approx G\_H^{-1} (V - G\_{\rm NH} M\_{\rm NH}) = G\_H^{-1} V - F E \, M\_{\rm NH} \tag{10}$$

where *FE* = *G*−<sup>1</sup> *<sup>H</sup> GNH* is the floor error matrix. Even though computationally expensive, the floor error matrix is specific of the sensor and thus needs to be computed only once.

In conclusion: Image reconstruction is carried out by multiplying the inverse of the regularized G-matrix in the fundamental hexagon by the measured visibilities and substracting from the result an estimation of the floor error, which is equal to the product of the floor error matrix times a scene model outside the fundamental hexagon. Needless to say, the closer the model to the actual image, the lower the reconstruction error.

#### *2.5. Matrix Extension and Inversion*

Assuming that the regularization described in Section 2.3 has been applied and that the minimum reciprocal grids are used, the complex matrix *GH* has, in the MIRAS case, 2791 rows and 4096 columns, corresponding respectively to the (*u*, *v*) unique grid points with measured visibilities (green star of Figure 2) and the fundamental hexagon in (*ξ*, *η*). Its condition number is about 3.2. The matrix *GH* can thus be inverted using a Moore–Penrose pseudoinverse algorithm, as in Corbella et al. [8], by a conjugate-gradient method as in Camps et al. [13] or by other methods listed also in this reference.

A different approach is proposed here. First, the matrix is extended in the (*u*, *v*) domain (rows) up to the whole principal hexagon (blue dots of Figure 2) using an average antenna pattern and unit fringe washing function. Using Equation (5), the G-matrix rows corresponding to the blue dots in Figure 2 (left), that is outside the star, are computed as

$$\mathcal{G}\_{\rm lm} = \Delta\_5^{\rm z} \Delta \eta \, \frac{|\mathcal{F}\_{\rm n}(\boldsymbol{\xi}, \boldsymbol{\eta})|^2}{\sqrt{1 - \boldsymbol{\xi}^2 - \eta^2} \boldsymbol{\Omega}} e^{-j2\pi(\boldsymbol{u}\_5^{\rm z} + \boldsymbol{v}\boldsymbol{\eta})}.\tag{11}$$

The extended G-matrix becomes then square with size *N*<sup>2</sup> *<sup>T</sup>* × *<sup>N</sup>*<sup>2</sup> *<sup>T</sup>* and keeps the condition number, so it can be straightforwardly inverted with standard algorithms. Note that, if all antenna patterns were substituted by the average pattern and the fringe washing function was neglected, this extended G-matrix would be a Fourier matrix with all columns multiplied by the average antenna pattern. The product *GHT* would become then equivalent to the product of a Fourier matrix and the modified brightness temperature, as expected. The proposed method can be viewed as a modification of the Fourier inversion to include different antenna patterns. Although not specifically reported elsewhere, this inversion method has been successfully implemented in the MIRAS Testing Software [19] since its first version.

#### *2.6. Image Reconstruction*

The brightness temperature map is recovered by multiplying the calibrated visibility by the inverted extended *GH* matrix. For X or Y polarization the brightness temperature is real and *G*−<sup>1</sup> *<sup>H</sup>* is hermitic, so the first term of Equation (10) can be written as

$$\mathbf{G}\_{H}^{-1}\mathbf{V} = \left[\mathbf{G}\_{H}^{-1}\begin{vmatrix}\mathbf{G}\_{H}^{-1}\end{vmatrix}\_{A}\mathbf{G}\_{H}^{-1}\begin{vmatrix}\mathbf{G}\_{H}^{-1}\end{vmatrix}\_{B}\right]\begin{bmatrix}V\_{A}\\V\_{0}\\V\_{B}\end{bmatrix} = \mathbf{G}\_{H}^{-1}\begin{vmatrix}\mathbf{V}\_{0} + 2\Re\varepsilon\left[\mathbf{G}\_{H}^{-1}\right]\_{A}V\_{A}\end{bmatrix},\tag{12}$$

where the subscript 0 refers to the origin (*u* = *v* = 0) and *A* and *B* are two grid point subsets (matrix columns) having opposite (*u*, *v*) signs. The ones used in the MIRAS Testing Software are shown in Figure 2 but the splitting is arbitrary. For points outside the star the visibility is ignored (see comment below about zero padding), so the corresponding columns of *G*−<sup>1</sup> *<sup>H</sup>* are not used. The last equality in the above equation holds because of the hermiticity property of visibility function. Since both *G*−<sup>1</sup> *H* 0 and *V*<sup>0</sup> are real, this equation can be written in a more compact form as

$$\mathcal{G}\_{H}^{-1}V = \mathfrak{Re}\left\{ \left[ \left. G\_{H}^{-1} \right|\_{0} \left. \mathcal{2} G\_{H}^{-1} \right|\_{A} \right] \left[ \begin{array}{c} V\_{0} \\ V\_{A} \end{array} \right] \right\}.\tag{13}$$

In MIRAS, using the minimum reciprocal grids, this operation involves the multiplication of a 4096 × 1395 complex matrix by a complex vector of 1395 elements.

If the complex polarimetric brightness temperature *Txy* is being imaged, the full *G*−<sup>1</sup> *<sup>H</sup>* matrix should be used instead of its real part, although points outside the star should also be discarded.

The second term of Equation (10) does not depend on the measurement since it is computed using a model outside the hexagon. The hermiticity property of the *X* and *Y* polarizations can also be used to reduce the size of the floor error matrix in this case.

$$FE = \mathfrak{Re}\left\{ \left. \left[ G\_H^{-1} \right|\_0 \left. 2G\_H^{-1} \right|\_A \right] \left. \begin{bmatrix} \left. G\_{NH} \right|\_0 \\\left. G\_{NH} \right|\_A \end{bmatrix} \right\} \right. \tag{14}$$

The size of the MIRAS single polarization floor error matrix using the minimum size grids is always 4096 × 4395 corresponding to the (*ξ*, *η*) points in both the fundamental hexagon and outside it respectively. This matrix is real for *Tx* and *Ty* and complex for *Txy*.

A final comment about zero-padding is worth mentioning. The above equations detail the computation of each one of the two terms of Equation (10) separately but in a consistent way. Considering the version of this equation written at the first equal sign (that is *TH* = *G*−<sup>1</sup> *<sup>H</sup>* (*V* − *GNHMNH*)), it comes out that zero padding outside the star means using in these points the product *GNHMNH*, but not zero. In practice this needs not to be done explicitly, as it suffices just to ignore the columns of *G*−<sup>1</sup> *<sup>H</sup>* outside the star.

#### *2.7. Apodization*

Due to the limited visibility coverage in the (*u*, *v*) plane, the reconstructed brightness temperature is affected by the Gibbs effect showing ripples around abrupt changes in the original scene, as for example coastlines. As it is well known from Fourier imaging, ripples can be reduced at the expense of degrading spatial resolution by using an apodization window in the original domain. In Corbella et al. [8] the apodization window was directly applied to the measured visibilities, which is correct if a Fourier inversion is used but it is at least questionable for the G-matrix technique. A more rigorous approach is to window the Fourier components of the reconstructed brightness temperature image. In this case, the apodized brightness temperature is related to the reconstructed brightness temperature *T* by

$$T\_{\text{apodizad}} = \mathcal{F}^{-1}\{\mathcal{W}\mathcal{F}\{T\}\},\tag{15}$$

where *W* is the window function, which in SMOS is always of Blackman type.

The reconstructed brightness temperature is defined in the (*ξ*, *η*) fundamental hexagon, so its Fourier components, to which the window function is applied, are defined in the (*u*, *v*) fundamental hexagon. The DFT operation involved in Equation (15) assumes implicitly that the reconstructed brightness temperature (*T*) is a periodic function in (*ξ*, *η*) replicating itself in hexagons adjacent to the fundamental one. Since in typical SMOS images the earth disk is at the bottom of the hexagon and the sky at the top, there are abrupt changes at the border of the fundamental period that may induce ripples. To mitigate them, constant temperature levels are subtracted from the sky and earth zones so as to have a zero mean image. Specifically, the constant temperature subtracted to the sky pixels is computed as the median of the recovered image in them, while the value subtracted to the Earth pixels is computed so as to cancel the Fourier component at the origin. The consequent reduction of the contrasts within the image cause the minimization of the associated ripples. The constant temperatures are added back after Fourier inversion.

This approach is strongly inspired on the incremental visibility image reconstruction method proposed in Camps et al. [13]. In that case, however, the method was applied to visibilities instead to frequency components, but the idea is the same.
