*2.8. Full Polarimetric Case*

Considering a baseline formed by two dual-polarization antennas, the full polarimetric discretized visibility equation [20] can be written in terms of G-matrices as

$$V\_{xx} = \mathbf{G\_{xx}^{RR}}\mathbf{T\_x} + G\_{xx}^{CC}T\_y + G\_{xx}^{RC}T\_{xy} + G\_{xx}^{CR}T\_{yx} \tag{16}$$

$$V\_{yy} = G\_{yy}^{CC} T\_x + \mathbf{G\_{yy}^{RR}} \mathbf{T\_y} + G\_{yy}^{CR} T\_{xy} + G\_{yy}^{RC} T\_{yx} \tag{17}$$

$$V\_{xy} = G\_{xy}^{RC} T\_x + G\_{xy}^{CR} T\_y + \mathbf{G\_{xy}^{RR}} T\_{xy} + G\_{xy}^{CC} T\_{yx} \tag{18}$$

$$V\_{yx} = G\_{yx}^{CR} T\_x + G\_{yx}^{RC} T\_y + G\_{yx}^{CC} T\_{xy} + \mathbf{G\_{yx}^{RR}}^{RR} \mathbf{T\_{yx}} \tag{19}$$

where, for example, *GRC xy* denotes the G-matrix computed according to Equation (5) using the Reference (co-polar) pattern of the X polarization antenna and the Cross-polar pattern of the Y polarization antenna. The terms marked in boldface are the dominant ones, in case of antennas with negligible cross-polar patterns, as considered in Equations (1) and (2).

Averaging redundant baselines, adding hermitic points and extending the matrix to the full hexagon is here carried out for each of the sub-matrices using the procedures detailed in Sections 2.3 and 2.5. Once this is done, the combination of the four equations can be written as *V* = *GT* where now *G* is a matrix with dimension 4*N*<sup>2</sup> *<sup>T</sup>* × <sup>4</sup>*Np*, where *<sup>N</sup>*<sup>2</sup> *<sup>T</sup>* is the total number of points in the fundamental hexagon and *Np* the number of points in the unit circle (4096 and 8491 respectively for the MIRAS minimum grid). The floor error mitigation through the application of a model outside the fundamental period can be carried out in the polarimetric case in exactly the same manner as done for a single polarization. The part of the extended G-matrix inside the fundamental period has now a size of 4*N*<sup>2</sup> *<sup>T</sup>* × <sup>4</sup>*N*<sup>2</sup> *<sup>T</sup>* and, when inverted, provides as result a square matrix:

$$\begin{aligned} \, \_{1}G\_{H}^{-1} = \begin{bmatrix} \, \_{I\mathcal{G}\_{11}} & \, \_{I\mathcal{G}\_{12}} & \, \_{I\mathcal{G}\_{13}} & \, \_{I\mathcal{G}\_{14}} \\ \, \_{I\mathcal{G}\_{21}} & \, \_{I\mathcal{G}\_{22}} & \, \_{I\mathcal{G}\_{24}} & \, \_{I\mathcal{G}\_{24}} \\ \, \_{I\mathcal{G}\_{31}} & \, \_{I\mathcal{G}\_{32}} & \, \_{I\mathcal{G}\_{33}} & \, \_{I\mathcal{G}\_{34}} \\ \, \_{I\mathcal{G}\_{41}} & \, \_{I\mathcal{G}\_{42}} & \, \_{I\mathcal{G}\_{43}} & \, \_{I\mathcal{G}\_{44}} \end{bmatrix} \, \end{aligned} \tag{20}$$

where *IGij* are submatrices of size *N*<sup>2</sup> *<sup>T</sup>* × *<sup>N</sup>*<sup>2</sup> *T*.

Applying the hermiticity property to the matrices of the first two rows using the procedures of Section 2.6 the following equations are obtained:

*Tx* = *e IG*<sup>11</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>11</sup> *A Vx*<sup>0</sup> *VxA* + *IG*<sup>12</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>12</sup> *A Vy*<sup>0</sup> *VyA* + (21) + *IG*<sup>13</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>13</sup> *A Vxy*<sup>0</sup> *VxyA* + *IG*<sup>14</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>14</sup> *A Vyx*<sup>0</sup> *VyxA* , *Ty* = *e IG*<sup>21</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>21</sup> *A Vx*<sup>0</sup> *VxA* + *IG*<sup>22</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>22</sup> *A Vy*<sup>0</sup> *VyA* + (22) + *IG*<sup>23</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>23</sup> *A Vxy*<sup>0</sup> *VxyA* + *IG*<sup>24</sup> <sup>0</sup> <sup>2</sup> *IG*<sup>24</sup> *A Vyx*<sup>0</sup> *VyxA* ,

$$T\_{xy} = IG\_{31}V\_x + IG\_{32}V\_y + IG\_{33}V\_{xy} + IG\_{34}V\_{yx} \tag{23}$$

in which the rows of the inverted G-matrix outside the star are always discarded. Note that, by definition, *Tyx* = *T*<sup>∗</sup> *xy*, so there is no additional equation for this term. Implementation wise it is found that this identity holds to the machine precision, which is a consistency indicator.

The floor error matrix is computed analogously to the case of single polarization measurements (Section 2.6) but using the larger G-matrix defined here. The result is a matrix four times the size of that for single polarization.

## **3. Results**

The image reconstruction methodology outlined in the previous sections is fully implemented in the MIRAS Testing Software [19]. This tool is being systematically used for scientific studies as a flexible alternative to the nominal SMOS Level 1 data products and with successful results, as reported for example in González-Gambau et al. [21].

Exploiting the processing capabilities of this software, the proposed algorithm was tested relying on a set of real MIRAS measurements. Specifically, a SMOS orbit was processed in full polarimetric mode (Section 2.8) using the all-LICEF calibration approach [12]. Figure 3 shows the retrieved geo-located images at X- and Y-polarizations corresponding to a single snapshot acquired on 24 January 2019 at 21:44 UTC, when the satellite was passing over the coast of Australia during an ascending orbit. The Blackman window is applied using the methods in Section 2.7.

**Figure 3.** Geo-located images, in the extended alias-free field of view, of the two snapshots of Figure 4.

Figure 4 shows the brightness temperatures of same snapshot in the (*ξ*, *η*) plane. Clearly, the image reconstruction algorithm is able to capture all the scene features in the whole hexagon, not only in the alias-free field of view. Ocean, land and sky areas are clearly distinguished and, at this scale, aliases impact is minor. Differences between *Tx* and *Ty* are due to the stronger *Ty* increments associated to both the sea/land and Earth/sky transitions with respect to *Tx*. The model used to cancel the floor error [*MNH* in Equation (10)] consists of a constant value in land zones (258 K for X-pol and 285 K for Y-pol), specular reflection in ocean areas using Fresnel reflection coefficients with climatology salinity and temperature [22], a constant 8.5K to account for atmospheric effects and the L-Band sky and Galaxy map available as auxiliary file in the SMOS data base.

Figure 5 shows the X-polarization case with both terms of Equation (10) drawn separately as well as the difference between them, which is the final reconstructed image. The left panel shows the result of multiplying the calibrated visibility by the inverse of the G-matrix (first term of Equation (10), the center panel is the floor error and the right panel the corrected image. As expected, the floor error is very small in regions where the earth aliases do not enter into the hexagon, which is the so-called "extended alias-free field of view", nominally used in SMOS images (as in the geo-located images of Figure 3). Most of the artifacts of the raw image at the left are effectively removed in the corrected one. Specifically, strong improvement in the image reconstruction is found in the sky area as well as in both the bottom left and right strips. Contrarily, much lower effects are exhibited in the extended alias-free field of view, where the floor error is minimum. Later it will be shown that there is indeed a small improvement in this area. Apodization of the corrected image with a Blackman window, using the procedure of Section 2.7, provides the final result seen at the left of Figure 4.

**Figure 4.** Image reconstruction result of two individual Soil Moisture and Ocean Salinity (SMOS) snapshots over the coast of Australia. Left is for X-pol and right is for Y-pol.

**Figure 5.** Floor error correction of the X-pol snapshot of Figure 4. Left and center images are first and second terms of Equation (10) respectively, and the difference is shown at right.

The scale of Figure 4 does not allow one to assess the effect of the floor error correction in the alias-free field of view. This is only possible using a differential image, plotting the brightness temperature bias with respect to its expected value. Additionally, averaging several snapshots is desirable in order to reduce thermal noise. Both requirements can be met if the scene is limited to snapshots over the ocean, for which a very comprehensive model is available from the SMOS science community (The authors would like to thank Joseph Tenerelli (OceanDataLab, France) for providing the ocean forward model).

Figure 6 shows the bias with respect to the model of all snapshots ranging in a latitudinal range of [−40◦, −5◦] over the Pacific ocean in an ascending orbit of 28 January 2011 (chosen quite often for these kinds of analysis by the SMOS Level 1 team). In the area imaging the sky, the model is the standard SMOS Galaxy map. All four polarimetric products are included in Figure 6, brightness temperature at X and Y polarizations and real and imaginary parts of the complex brightness temperature. Basic statistics, referring to the alias-free field of view only, are shown at the bottom of each panel of Figure 6. Namely, spatial standard deviations and mean values are indicated by *σ* and *T* respectively. The relatively large negative bias in the Y-pol image may be due to poor modeling at high incident angles. In any case, these residual error images are compatible with the ones obtained with the SMOS Level-1 Operational Processor and have similar statistics.

**Figure 6.** Spatial bias computed as the difference between the reconstructed brightness temperature and a model for the ocean. Blackman window is applied.

The same images have been produced without correcting the floor error, that is using only the first term of Equation (10). Results are shown in Figure 7. In this case, the error outside the alias-free field of view increases dramatically, especially in the areas in which the earth enters the hexagon. In the alias-free field of view there is a small impact in spatial bias, quantified in the standard deviation shown at the bottom. In all cases it increases with respect to the numbers provided in Figure 6. Note that the extended alias-free field of view is quite well recovered even if no correction is applied.

**Figure 7.** Spatial bias of Figure 6 without removing the floor error.

#### **4. Discussion**

Years before SMOS launch in November 2009, methods for 2D interferometric radiometer imaging were derived by different groups with the aim of having working algorithms as soon as visibility measurements were provided by MIRAS [8,11,23]. Radioastronomy heritage showed to be not directly applicable due to different instrument layouts, as for example the reduced antenna spacing, and especially because of the larger instantaneous field of view typical of earth observation. Working with simulated data it was early discovered that after a complete forward-backward simulation of a known scene, the original brightness temperature was not perfectly recovered even in the alias-free field of view. This misfit was called "scene-dependent bias" and was attributed to an underdetermination of the mathematical problem. An efficient mitigation algorithm was proposed in Corbella et al. [24], used in Anterrieu et al. [25] and improved later in Camps et al. [13]. This consists of subtracting from the visibility measurements different contributions estimated by simulation, and inverting the resultant "differential visibilities". Some of the contributions are later added back to the reconstructed image to recover the final brightness temperature map. The most recent implementation of this idea, described in Khazâal et al. [26], is included in the version 7 of the SMOS Level 1 Operational Processor. The method, even in its simplest version, is highly effective in canceling the sky aliases so expanding the alias-free field of view—limited by the unit circle aliases—to the extended alias-free field of view, limited by the earth shape. In its latest version [26], using a sophisticated model as "artificial scene" it is able to further reduce the error in the alias-free field of view.

All these methods are based on inverting only the portion of the G-matrix inside the fundamental hexagon discarding the rest. As already pointed out in Corbella et al. [17], discarding the G-matrix outside the fundamental hexagon is responsible for the appearance of the aliases. If all antenna patters were identical, this error would become limited to only the aliasing regions. In the real case, with different antenna patterns, the error is indeed larger in these regions but spreads also in the alias-free field of view (see center panel of Figure 5). The G-matrix in the whole unit circle is always used to estimate the corresponding visibility of the model or artificial scene, so imaging differential visibilities can be interpreted as removing an estimation of the aliases.

Using the concept of extended G-matrix of Section 2.5, it is easily shown that the method of imaging differential visibilities is equivalent to removing the floor error. The reconstructed brightness temperature from differential visibilities is

$$T = G\_H^{-1}(V - GM) + M\_H = G\_H^{-1}V - (G\_H^{-1}G - lI)M,\tag{24}$$

where *U* is a matrix with the same size of *G* equal to the identity matrix for the columns inside the hexagon and zero outside *U* = [*I* 0]. In this equation, the G-matrix is extended, so having as many rows as number of points in the fundamental (*u*, *v*) hexagon, and thus *GH* is a square and invertible matrix. The first term of this equation is the same as that of Equation (10). The second term can be expanded as

$$\left(\mathbf{G}\_{H}^{-1}\begin{bmatrix}\ \mathbf{G}\_{H} & \mathbf{G}\_{\rm NH} \end{bmatrix} - \begin{bmatrix} I & \mathbf{0} \end{bmatrix}\right) \begin{bmatrix} M\_{H} \\ \ M\_{\rm NH} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{G}\_{H}^{-1}\mathbf{G}\_{\rm NH} \end{bmatrix} \begin{bmatrix} M\_{H} \\ \ M\_{\rm NH} \end{bmatrix} = \mathbf{G}\_{H}^{-1}\mathbf{G}\_{\rm NH}M\_{\rm NH} \tag{25}$$

which coincides with Equation (10). It is important to point out that this is only true if the visibility simulation of the model uses a G-matrix in the full unit circle defined in the same grid as the one used for inversion. That is, the matrix inverted is a subset of the complete one. Conceptually, a different matrix could be used since *V* = *GT* is just a mathematical model of what the instrument actually measures. Also, the equivalence shown is based on the use of the extended G-matrix concept defined in Section 2.5. Different implementations of matrix inversion in Equation (24) are not equivalent. The main advantage of using the floor error matrix defined in Equation (10) is that processing does not

require estimating visibilities with the full G-matrix and does not require adding back any artificial scene to the result. It is a correction applied directly to the reconstructed image.

The outcome of the proposed methodology is the brightness temperature map of the scene. This is not the case for the SMOS Level 1 Operational Processor which, based in the proposal of Anterrieu et al. [11], defines the Level 1B data as the spatial frequency components of the brightness temperature. This is done through the concept of J-matrix, the concatenation of the G-matrix with the Fourier operator as dully explained in Khazâal et al. [26]. This approach is a different method to regularize the system of equations that relates all measured visibility, including redundant baselines, to brightness temperature. Due to redundancies the G-matrix in this case is ill-posed but the corresponding J-matrix is well behaved and can be inverted. As pointed out in Section 2.7, since the brightness temperature is not a periodic function, when recovering it from the corresponding frequencies, ripples can appear in the limits of the fundamental hexagon. In practice, this does not happen due to the imaging of differential visibilities used in the processor. The Level 1B data consists actually of frequency components of the difference between the brightness temperature and an artificial scene. After reconstructing the map including apodization, the latter is added back. Versions prior to 7 of the SMOS Level 1 Operational Processor use as artificial scene a constant in the earth disk and a sky map for the sky. In the end the constant earth is added to all points, and this is why the nominal SMOS images at Level 1B in the full hexagon show high temperature in the sky (see Figure 8). Here it has been shown that averaging redundant baselines improves the condition of the matrix and makes it invertible, so there is no need to compute the J-matrix. As a post-processing, however, the frequency components are computed in order to apply the Blackman window, as explained in Section 2.7.

For completeness, equivalent images as those shown in Figures 4 and 6 are provided respectively in Figures 8 and 9 using SMOS Level 1B data read from the ESA SMOS Online Dissemination Service, produced by version 6.21 of the SMOS Level 1 Operational Processor.

**Figure 8.** Same snapshots as in Figure 4 but using data from the SMOS Level 1 Operational Processor version 6.21. Deimos Engenharia, Lisbon (Portugal).

**Figure 9.** *Cont*.

**Figure 9.** Same differential images as those in Figure 6 but using data from the SMOS Level 1 Operational Processor version 6.21. Deimos Engenharia, Lisbon (Portugal).

#### **5. Conclusions**

Accurate brightness temperature retrieval in 2D interferometric radiometers with a large field of view is not straightforward. The simplest approach consists of applying an inverse Fourier transform, which provides reasonable results in the alias-free field of view but large errors outside. The G-matrix approach takes into account individual differences between antenna patterns and also the fringe washing function and thus provides images of improved quality. In any case, the major error contribution is localized in the alias regions where brightness temperatures from several spatial directions overlap. For zones where the overlap is the sky the impact is low, but where aliases are produced by the earth it may become significant, affecting also the alias-free field of view. The origin of these errors is found in the contribution to the visibility of the brightness temperature from directions that fall outside the fundamental period. Subtracting an estimation of these visibilities to the measurements greatly reduces the error. Correction is carried out through the definition of the floor error matrix and relying on a brightness temperature model defined only outside the fundamental period.

Expressing the visibility equation as a linear system of equations, image reconstruction can be implemented through the inversion of the associated matrix, called G-matrix. Unfortunately, the G-matrix results are ill-conditioned and cannot be inverted without regularization. This is obtained by averaging redundant visibilities and extending the G-matrix to fill the whole hexagon in the (*u*, *v*) plane, resulting in a square and well-conditioned matrix that is easy to invert. Using this approach, combined with the floor error removal, provides high quality images in the full hexagon, and especially in the extended alias-free field of view. The method is applicable to the full polarimetric operation of SMOS, with the only difference being increasing the size of the G-matrix. Results of SMOS complex images and full polarimetric error maps over ocean demonstrate the procedure.

The inversion method presented in this paper uses the minimum number of grid points, allowing for very fast and efficient programming. It has been implemented in the MIRAS Testing Software for several years and is being successfully used by science teams for salinity and soil moisture retrievals. In any case, the procedures outlined are not intended to replace the ones used by the SMOS Level 1 Operational Processor, but they are presented to the community with the objective of reporting other means of solving the same problem with similar results.

**Author Contributions:** Conceptualization, I.C.; Investigation, I.C. and I.D.; Methodology, F.T. and N.D.; Software, I.C.; Supervision, M.M.-N.; Validation, V.G.-G.; Writing—original draft, I.C.; Writing—review & editing, V.G.-G. and M.M.-N.

**Funding:** This research was funded by the European Space Agency through SMOS P7 subcontract DME CP12 no. 2015-005 with Deimos Enginheria (Portugal) and by Ministerio de Economia, Industria y Competitividad, Gobierno de España, projects TEC2014-58582-R, TEC2017-88850-R and ESP2015-67549-C3-1-R.

**Conflicts of Interest:** The authors declare no conflict of interest.
