# A New Direct Sun Correction Algorithm for the Soil Moisture and Ocean Salinity Space Mission

Ali Khazaal, François Cabot, Eric Anterrieu, Yann H. Kerr

► To cite this version:

Ali Khazaal, François Cabot, Eric Anterrieu, Yann H. Kerr. A New Direct Sun Correction Algorithm for the Soil Moisture and Ocean Salinity Space Mission. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13, pp.1164-1173. 10.1109/JSTARS.2020.2971063. hal-03254361

HAL Id: hal-03254361

https://hal.inrae.fr/hal-03254361

Submitted on 9 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A New Direct Sun Correction Algorithm for the Soil Moisture and Ocean Salinity Space Mission

Ali Khazâal<sup>ID</sup>, François Cabot, Eric Anterrieu<sup>ID</sup>, Member, IEEE, and Yann H. Kerr<sup>ID</sup>, Fellow, IEEE

**Abstract**—The interferometric measurements provided by the Soil Moisture and Ocean Salinity (SMOS) satellite are contaminated by solar radiations due to a quasi-permanent presence of the Sun in the field of view of the instrument. A correction algorithm is developed and implemented in the SMOS image reconstruction processor to remove the contribution of these radiations. Since the Sun is seen under an angle much smaller than the angular resolution of SMOS, this algorithm assimilates the Sun disc to a single source located at its center. However, strong residuals persist in the retrieved images mainly due to the presence of different distributed smaller sources of brightness in the Sun disc. Thus, the performance of the algorithm degrades rapidly as soon as the dominant source is far from the center of the disc. This behavior shows the limitation of the algorithm and its strong dependencies to the position of the Sun. In this article, we propose two other algorithms to account for the solar radiations. The first one is an iterative optimization algorithm based on estimating the position of the dominant source of brightness inside the Sun disc. The second algorithm is based on estimating the contribution of the different sources of brightness within the Sun disc using an oversampled grid by solving a constrained least square optimization problem with an explicit solution.

**Index Terms**—Constrained optimization, inverse problems, regularization, solar radiations.

## I. INTRODUCTION

THE European Space Agency (ESA) Soil Moisture and Ocean Salinity (SMOS) satellite provides global measurements over the entire earth that are used to retrieve soil moisture and sea surface salinity [1]. The instrument on board SMOS is a two-dimensional *L*-band interferometer equipped with 69 equally spaced antennas/receivers operating at the central frequency $f_o = 1.413$ GHz and located along the arms of a *Y*-shaped array [2]. The interferometric measurements, also called complex visibilities, provided by SMOS are the cross correlations of the signals collected by each pair of receivers [3] and are used to retrieve the brightness temperature distribution of the observed scenes [4]. But these signals include the contribution of

other sources of brightness than the surface emission due to their presence inside the field of view (FOV) of the satellite. These sources of brightness are referred to as “foreign sources” and their contribution are removed prior to the image reconstruction to reduce any associated reconstruction biases.

Among these sources of brightness, we can distinguish the direct Sky, Sun, and Moon radiations, as well as the reflected Sun radiations (Sun glint) [5]. This article focuses only on the impact of the direct solar radiations on the retrieved images. In the current SMOS level 1 operational data processor (L1OP), a correction algorithm is already implemented [6] where the Sun disc is assimilated to a single point located at its center but strong residuals persist in the retrieved images as we will show hereafter. The main problem is that the Sun is not a single source but rather a disc with distributed smaller sources at very high temperature as showed using ground-based radio telescope [7]. In this article, we propose two other algorithms: The first one is an iterative optimization algorithm that identifies the position of the main source, which optimizes the Sun correction performance. The second algorithm estimates the different sources of brightness within the Sun disc using an *a priori* oversampling grid around the position of the Sun and by solving a constrained least square optimization problem with an explicit solution.

## II. SMOS IMAGE RECONSTRUCTION PROCESSOR

The SMOS image reconstruction processor within the L1OP converts the interferometric measurements or visibilities $V_o$ into brightness temperature images $T$. Since SMOS is a fully polarimetric instrument as described in [8], the retrieved images are obtained in $X$, $Y$, and $XY$ polarizations [4]. For simplification reasons, the image reconstruction algorithm is recalled here without referring to the different polarizations.

Before the reconstruction, the visibilities $V_o$ are first corrected from the contribution of the physical temperature of the receivers $V_{rec}$ by using a set of measured calibrated visibilities obtained from deep Sky observations [9]. Then, the remaining visibilities are corrected from the contribution of the foreign sources of brightness present in the FOV (direct Sun, Moon, and Sky radiations, as well as the Sun glint). These sources, if not eliminated, will contaminate the retrieved brightness temperatures thus leading to erroneous retrieval of ocean salinity and soil moisture. Therefore, their contribution, namely $V_{fs}$ is removed prior to the image reconstruction in order to remove any associated reconstruction biases.

Manuscript received November 27, 2019; revised January 20, 2020; accepted January 27, 2020. Date of publication March 20, 2020; date of current version April 8, 2020. This work was supported by the European Space Agency. (Corresponding author: Ali Khazâal.)

Ali Khazâal is with RDIS Conseils, 31400 Toulouse, France, and also with the Center for the Study of the BIOSphere, Université de Toulouse, CNRS, CNES, INRA, and IRD, 31400 Toulouse, France (e-mail: ali.khazaal@cesbio.cnes.fr).

François Cabot, Eric Anterrieu, and Yann H. Kerr are with the Center for the Study of the BIOSphere, Université de Toulouse, CNRS, CNES, INRA, and IRD, 31013 Toulouse, France (e-mail: francois.cabot@cesbio.cnes.fr; eric.anterrieu@cesbio.cnes.fr; yann.kerr@cesbio.cnes.fr).

Digital Object Identifier 10.1109/JSTARS.2020.2971063
The corrected visibilities, $V = V_o - V_{rec} - V_{fs}$, are then used to retrieve $T$ through the following linear relationship:

$$V(u, v) = \mathbf{G}(u, v, \xi, \eta) \times T(\xi, \eta) \quad (1)$$

where $\mathbf{G}$ is the instrument modeling operator that depends mostly on the antenna voltage patterns of SMOS [4]; $(u, v)$ are the spatial frequencies in the Fourier domain, which are confined to a limited star-shaped region and coincide with the nodes of an hexagonally sampled grid [10]; and $(\xi, \eta)$ are the director cosine coordinates defining the instantaneous FOV of SMOS, which has the extent of the unit circle. However, the reconstructed FOV is an hexagon inside the unit circle and bounded by earth and Sky aliases since the Shannon–Nyquist sampling criterion is not respected [11]. The regions with only Sky aliases are kept since the Sky temperature is well known and its impact can be accounted for. The retrieved brightness temperature $T$ is therefore restricted to an extended alias-free field of view region (EAFFOV).

To retrieve $T$ from (1), it has been shown in [12] that this linear problem is ill-posed since the number of visibilities $V$ is always smaller than the number of pixels of $T$. A physical regularization approach was also proposed in [12] that takes advantage of the limited frequency coverage inside the star to do the reconstruction in the Fourier domain instead of the spatial domain. Indeed, since the number of frequencies inside the star-shaped region is always smaller than the number of visibilities due to the redundancies of the instrument, the reconstruction problem in the Fourier domain is well-posed, and (1) can be written as

$$V = \mathbf{J} \times \hat{T} \quad (2)$$

where $\hat{T}$ are the Fourier components of $T$ inside the star, $\mathbf{J} = \mathbf{G}\mathbf{U}^*\mathbf{Z}\mathbf{J}^{-1}$ is a well-conditioned reconstruction operator, $\mathbf{U}$ is the Fourier transform operator, and $\mathbf{Z}$ is the zero-padding operator beyond the star. The unique solution of (2) is given by

$$T_r = \mathbf{U}^*\mathbf{Z}\mathbf{J}^+V \quad (3)$$

where $\mathbf{J}^+ = (\mathbf{J}^*\mathbf{J})^{-1}\mathbf{J}^*$ is the Moore–Penrose pseudoinverse of $\mathbf{J}$.

The retrieved brightness temperatures $T_r$ suffer from important spatial biases mainly due to the aliasing in the reconstructed FOV but also to the Gibbs oscillations due to the sharp transition beyond the limited star-shaped frequency coverage [13]. To reduce the biases, the L1OP applies a differential reconstruction by removing from the visibilities $V$ the contribution $\hat{V} = \mathbf{G}\hat{T}$ of an artificial scene $\hat{T}$ as close as possible to the observed one $T$ [14]. In this case, a differential system, $\delta V = \mathbf{G}\delta T$, is solved, with $\delta V = V - \hat{V}$ and $\delta T = T - \hat{T}$, and yielding to the following solution:

$$\delta T_r = \mathbf{U}^*\mathbf{Z}\mathbf{J}^+(V - \mathbf{G}\hat{T}). \quad (4)$$

Finally, the retrieved temperature $T_r$ is obtained by adding the artificial scene $\hat{T}$ to $\delta T_r$

$$T_r = T_r + \hat{T}. \quad (5)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Sun position $(\xi_s, \eta_s)$ (small circle in yellow) inside the FOV of SMOS (circle in cyan) and its alias (small circle in red) inside the reconstructed hexagon (in red). The ellipsoid in blue corresponds to the earth as seen by the instrument. The region between the earth horizon and the unit circle corresponds to the Sky. The dot (or dashed) circles and lines correspond to the six aliases.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Retrieved brightness temperature $T_r$ in $X$ (left) and $Y$ (right) polarizations without any correction of the direct solar radiations. This example is taken from an orbit over the Pacific ocean in October 20, 2011.</figcaption></figure>

### III. DIRECT SOLAR RADIATIONS IN SMOS IMAGES

The Sun is present in almost all the SMOS images either above or below the antenna plane [15]. If the Sun is above the antenna plane, the solar radiations are captured by the front-lobes antenna patterns and their contribution to the measurements is maximum. If it is below the antenna plane, the radiations are captured by the back-lobes antenna patterns (their gain is very small with respect to those in front) with a minimum and mostly local contribution.

We also indicate that the Sun is always located within the Sky region and outside the reconstructed hexagon, as shown in Fig. 1. But, one of its six aliases is most of the time located inside the EAFFOV as also shown in Fig. 1 and thus contaminating the retrieved brightness temperatures. An example of such contamination is shown in Fig. 2 for two snapshots taken over the Pacific ocean in October 20, 2011. We notice here a bright source that contaminates the retrieved images with a brightness temperature that can reach several thousands of Kelvin. This contamination also spreads into six directions (so-called Sun tails) around the Sun position due to the hexagonal grid of SMOS. Therefore, it is necessary to remove the contribution of the direct Sun to avoid any degradation of the SMOS retrievals. In the next sections, we first describe the Sun correction algorithm already implemented
in the L1OP and then we propose two other correction algorithms. Finally, the three algorithms are compared using SMOS real data.

## IV. SINGLE-SOURCE ALGORITHM

The correction algorithm implemented in the SMOS L1OP considers the Sun as a single source and not a disk. This assumption is valid since the Sun is seen under an angle smaller than the angular resolution of SMOS. In fact, the Sun diameter at L-band is about 0.5°, whereas the angular resolution of SMOS ranges from 1.60° to 2.25°. The position of the Sun ($\xi_s, \eta_s$) in the director cosine coordinates corresponds to the center of the sun disc calculated using the method proposed by Van Flandern and Pulkkinen [16] with a precision smaller than 0.01°.

The algorithm used for now was initially introduced in [6] and then modified in [17]. The Sun is modeled as a Dirac signal and its contribution $V_{\text{Sun}}$ (part of $V_{fs}$) is calculated by estimating first its brightness temperature $T_s$ and then multiply it by the system response function vector corresponding to the Sun position $G(u, v, \xi_s, \eta_s)$

$$V_{\text{Sun}}(u, v) = G(u, v, \xi_s, \eta_s) \times T_s(\xi_s, \eta_s). \quad (6)$$

The estimation of the Sun brightness temperature $T_s$ is done using the following equation:

$$T_s = \frac{\mathcal{F}^{-1}(V_o - V_{\text{rec}})(\xi_s, \eta_s) - \bar{T}_{\text{scene}}}{\mathcal{F}^{-1}(G(u, v, \xi_s, \eta_s))(\xi_s, \eta_s)} \quad (7)$$

where the following statements hold.

1) $\mathcal{F}^{-1}(V_o - V_{\text{rec}})(\xi_s, \eta_s)$ is the inverse Fourier transform applied to $V_o - V_{\text{rec}}$ over the Sun position $(\xi_s, \eta_s)$.

2) $\bar{T}_{\text{scene}}$ is the average temperature of the observed scene. It is computed as the average of $\mathcal{F}^{-1}(V_o - V_{\text{rec}})$ over a square of 11 × 11 pixels, centered at $(\xi_s, \eta_s)$.

3) $\mathcal{F}^{-1}(G(u, v, \xi_s, \eta_s))(\xi_s, \eta_s)$ is the inverse Fourier transform of the system response function vector $G(u, v, \xi_s, \eta_s)$ calculated at the Sun position $(\xi_s, \eta_s)$.

Fig. 3 shows the retrieved brightness temperature images $T_r$ using the single-source algorithm for the two snapshots of Fig. 2 where we observe important residuals over the Sun alias. Also shown in Fig. 3 are the differences of reconstruction $\delta T_r$ with and without applying the Sun correction algorithm. These differences show that even though the Sun correction algorithm is only applied over one single pixel, its impact is affecting all the area surrounding it, mainly the tails but also creating oscillations all over the image. Therefore, any under or over correction of the Sun radiations, will propagate through the reconstruction affecting not only the position of the alias of the Sun but also the tails and to a less degree the surrounding pixels. As a consequence, all pixels of the Sun disc and the tails are for now flagged and not used by the SMOS level 2 processor in the retrieval of salinity and soil moisture.

## V. ITERATIVE ALGORITHM

The single-source algorithm is very dependent on the position of the Sun in the FOV. Two possible problems can be encountered

<figure><img src="image_4.png" /><figcaption>Fig. 3. [Top] Retrieved brightness temperature $T_{r1}$ in X (left) and Y (right) polarizations after applying the single-source Sun correction algorithm. [Bottom] Difference of brightness temperature $\delta T_r = T_r - T_{r1}$ without any correction of the direct solar radiations and with the single-source algorithm.</figcaption></figure>

in this algorithm: 1) the position of the Sun is not calculated with a sufficient precision; or 2) the dominant source of brightness inside the Sun disc is far from the center of the disc. In both cases, we need a better estimation of the Sun position, namely $(\bar{\xi}_s, \bar{\eta}_s)$.

To do that, we propose here an iterative algorithm that uses the SMOS data to estimate the new position. The idea is to minimize the standard deviation of the retrieved differential brightness temperature $\delta T_r$ of (4) over a region $\mathcal{D}$ around the Sun position

$$\langle \bar{\xi}_s, \bar{\eta}_s \rangle = \min_{(\xi_s, \eta_s) \in \mathcal{D}} \sigma \left( \mathbf{U}^* \mathbf{Z} \mathbf{J}^+ (V - \tilde{V} - V_{\text{Sun}}(\xi_s, \eta_s) - \mathbf{G} \tilde{T}) \right) \quad (8)$$

with $\mathcal{D}=\{(\xi_s, \eta_s), \text{s.t.} \sqrt{(\xi_s - \xi_s)^2 + (\eta_s - \eta_s)^2} < d\}$, where $d = 0.0103$ is the sampling rate of the hexagonal grid of SMOS. Normally, the standard deviation over such a limited zone is very close to zero if the Sun is well corrected.

To solve (8), we use the interior point algorithm of the MATLAB built in function `fmincon.m` with $(\xi_s, \eta_s)$ as an initial solution. This algorithm attempts to take a direct Newton step at the beginning of each iteration. If it cannot, it attempts a conjugate gradient step. Finally, the single-source algorithm is applied over the new estimated position $(\bar{\xi}_s, \bar{\eta}_s)$.

Fig. 4 shows the retrieved brightness temperature map $T_r$ using the iterative algorithm for the two previous snapshots of Fig. 2. The estimated position $(\bar{\xi}_s, \bar{\eta}_s)$ was found equal to $(-0.9219, 0.2886)$ for the snapshot in X polarization and $(-0.9227, 0.2896)$ for the one in Y polarization, whereas the calculated position $(\xi_s, \eta_s)$ is $(-0.9217, 0.2899)$ in X polarization and $(-0.9217, 0.2901)$ in Y polarization. The Euclidean distance between the estimated dominant source and the position of the center of the Sun is 0.0013.
<figure><img src="image_5.png" /><figcaption>Fig. 4. [Top] Retrieved brightness temperature $T_{r2}$ in X (left) and Y (right) polarizations after applying the iterative Sun correction algorithm. [Bottom] Difference of brightness temperature $\delta T_r = T_{r1} - T_{r2}$ between the single-source and the iterative Sun correction algorithms.</figcaption></figure>

Looking at Fig. 4, we notice that the residuals are almost totally corrected. The difference of brightness temperature between the single-source and the iterative Sun correction algorithms is also shown in Fig. 4. This difference spreads mainly from the Sun alias position through the tails, whereas the rest of the image is also affected by very weak oscillations. Such difference could be very significant and especially for sea surface salinity retrievals since it is very sensitive to brightness temperature [18].

On the other hand, two drawbacks of this algorithm are identified as follows.

1) Its execution time is very high compared to the single-source algorithm mainly due to the reconstruction procedure that is repeated for each iteration. We have at least a factor 100 between both algorithms making it difficult to ingest it into an operational processor.

2) This algorithm estimates the Sun position using SMOS visibilities. Therefore, the accuracy of the estimated position varies with the level of contamination of the visibilities by other sources of brightness and mainly radio frequency interference (RFI) highly present in the SMOS data [19].

## VI. MULTIPLE SOURCES ALGORITHM

Despite the fact that the position of the center of the Sun ($\xi_s, \eta_s$) is located very accurately, it seems that the Sun is not a point source located at its center, but rather a distribution of much smaller sources at very high temperature. This has been documented using ground-based radio telescope [7], as shown in Fig. 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Image of the Sun Stokes 1 brightness temperature (color range goes from 0 to 2.2 × 10$^6$ K) as reported by the ground-based radio telescope [7].</figcaption></figure>

The size and distances between the individual sources of the Sun within the disc are well below SMOS geometric resolution as referred above. But, the poor performance of the single-source Sun correction algorithm and the improvements obtained using the iterative algorithm have shown that it is important to localize these sources as much precise as possible. To improve the localization of the sources, the algorithm proposed hereafter uses an oversampling grid within the Sun disc around ($\xi_s, \eta_s$). As for the Sun brightness temperatures, they are estimated by solving a constrained least square optimization problem with an explicit solution.

### A. Oversampling Grid

To localize finely all sources within the Sun disc, we propose to use an oversampled grid around the computed position of the Sun ($\xi_s, \eta_s$). This fine grid preserves the hexagonal sampling and is limited to the maximum expected extent of the active Sun at L-band, which is slightly larger than the optical apparent diameter of the Sun. This extent is always smaller than the pixel spacing $d = 0.0103$ of the original ($\xi, \eta$) grid. Thus, it makes sense to set the oversampling grid support equal to $d$.

Fig. 6 shows an example of the oversampling grid ($\xi_g, \eta_g$) surrounding the Sun center ($\xi_s, \eta_s$) for an oversampling ratio $R$, which leads to $N_g = (1 + \sum_{r=1;R} 6 \times (r-1))$ subpixels including ($\xi_s, \eta_s$). The value of $R$ is fixed by the following three parameters.

1) The computational time of the overall algorithm that is an increasing function of $R$.

2) The density of the oversampled grid: We have noticed that the Euclidean distance between the center of the Sun disc ($\xi_s, \eta_s$) and the estimated ones ($\bar{\xi}_g, \bar{\eta}_g$) is in the range of $10^{-3}$ in the ($\xi, \eta$) grid. Therefore, the distance between two subpixels ($\frac{d}{R}$) should be close to this value.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Oversampling grid ($\xi_g, \eta_g$) shown in red around the Sun position (in yellow) and drawn over the standard grid (in blue). This example is shown for an oversampling ratio $\mathcal{R} = 4$, which leads to $N_g = 37$ subpixels.</figcaption></figure>

<table><caption>TABLE I
CONDITIONING NUMBER (C) OF THE SYSTEM RESPONSE FUNCTION G AND SPACING BETWEEN TWO SUBPIXELS ($\frac{d}{\mathcal{R}}$) CALCULATED FOR THE $N_g$ SUBPIXELS ($\xi_g, \eta_g$) IN FUNCTION OF THE OVERSAMPLING RATIO $\mathcal{R}$.</caption><thead><tr><th><i>R</i></th><th><i>N</i><sub>g</sub></th><th><i>d</i>/<i>R</i></th><th>C</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0.0103</td><td>1</td></tr><tr><td>2</td><td>7</td><td>0.0065</td><td>236.3</td></tr><tr><td>3</td><td>19</td><td>0.0043</td><td>1.5 × 10<sup>4</sup></td></tr><tr><td><b>4</b></td><td><b>37</b></td><td>0.0032</td><td><b>2.5 × 10<sup>6</sup></b></td></tr><tr><td>5</td><td>61</td><td>0.0026</td><td>3.2 × 10<sup>8</sup></td></tr><tr><td>6</td><td>91</td><td>0.0022</td><td>2.1 × 10<sup>11</sup></td></tr><tr><td>7</td><td>127</td><td>0.0019</td><td>3.7 × 10<sup>13</sup></td></tr><tr><td>8</td><td>169</td><td>0.0016</td><td>2.3 × 10<sup>15</sup></td></tr><tr><td>9</td><td>217</td><td>0.0014</td><td>3.6 × 10<sup>15</sup></td></tr><tr><td>10</td><td>271</td><td>0.0013</td><td>4.8 × 10<sup>15</sup></td></tr></tbody></table>

3) The conditioning of the system response matrix G calculated for the oversampled grid ($\xi_g, \eta_g$) and defined as the ratio between the maximum and the minimum singular values (SV) of G

$$C = \frac{\max(\mathbf{SV}(G))}{\min(\mathbf{SV}(G))}. \qquad (9)$$

The closer C is to 1, the more $G(u, v, \xi_g, \eta_g)$ is well conditioned.

Table I summarizes the values of the spacing between two subpixels ($\frac{d}{\mathcal{R}}$) and the conditioning number C in function of the oversampling ratio $\mathcal{R}$. A compromise between the computational time, the conditioning of G, and still having a dense grid is needed. Based on that table, we have chosen to set $\mathcal{R} = 4$ with $N_g = 37$ subpixels. The $G(u, v, \xi_g, \eta_g)$ matrix is in this case ill-conditioned and a regularization is needed to ensure a unique and stable solution.

## B. Sun Temperature Estimation

The single-source algorithm is initially applied over $V_o - V_{rec}$. For this algorithm, we propose to apply it over the differential visibilities $\delta V = V_o - V_{rec} - V_{fs} - \bar{V}$. By this, we reduce as much as possible the impact of the observed scene T in the

<figure><img src="image_8.png" /><figcaption>Fig. 7. Regions of the FOV considered in the criterion J to estimate the $N_g$ Sun temperatures $T_g$. [Blue] Polluted region by the solar radiations including the disc of the Sun alias and the six tails. [Green] Clean region around the Sun alias considered as not contaminated by the solar radiations.</figcaption></figure>

estimation of the $N_g$ Sun temperatures $T_g$. To estimate $T_g$, we propose to minimize the following criterion:

$$\mathbf{J}(\xi_g, \eta_g) = \left\| \mathcal{F}^{-1}(\delta V)_p - \overline{\mathcal{F}^{-1}(\delta V)}_c - \sum_{g=1}^{N_g} T_g \mathcal{F}^{-1}(\mathbf{G}(\xi_g, \eta_g)) \right\|^2 \quad (10)$$

where $\mathcal{F}^{-1}$ represents the inverse Fourier transform; $p$ represent the pixels inside the hexagon polluted by the solar radiations including the disc of the Sun alias and the corresponding six tails, as shown in blue in Fig. 7; and $c$ represent the clean pixels around the Sun alias considered as not contaminated by the solar radiations and shown in green in Fig. 7. The term $\overline{\mathcal{F}^{-1}(\delta V)}_c$ represents the average of the residual image over the $c$ clean pixels.

As mentioned above, the $G(u, v, \xi_g, \eta_g)$ matrix used in the criterion J is ill-conditioned and should be regularized. To overcome this problem, we propose to add two constraints to J weighted by the same Lagrangian $\lambda$ given by

$$\mathbf{L}(\xi_g, \eta_g) = \left\| \sum_{g=1}^{N_g} \left( T_g - \frac{T_o}{N_g} \right) \right\|^2 + \left\| \left( \sum_{g=1}^{N_g} T_g \right) - T_o \right\|^2 \quad (11)$$

where $T_o$ is the estimated brightness temperature using the single-source algorithm. The first constraint will limit the estimated temperature of each source around an average value equal to $\frac{T_o}{N_g}$, whereas the second one will limit the sum of the estimated temperatures to $T_o$. Now, to estimate $T_g$, we propose to minimize the following constrained optimization problem:

$$\min_{(\xi_g, \eta_g)} \mathbf{J}(\xi_g, \eta_g) + \lambda \mathbf{L}(\xi_g, \eta_g). \quad (12)$$

The choice of using the same Lagrangian $\lambda$ was because of the similarity of the two constraints and more particularly to simplify the optimization. The value of $\lambda$ is chosen as the tradeoff between the deviation from the data (represented by J) and the deviation from the solution (represented by L). For a fixed value of $\lambda$, (12) is linear with respect to $T_g$ and can be written in the following
<figure><img src="image_9.png" /><figcaption>Fig. 8. L curve representing the deviation from the solution L as a function of the residual J for all values of the regularization parameter λ. The optimal value of λ was found equal to 10⁻⁶ (dot in red).</figcaption></figure>

linear form:

$$ \begin{pmatrix} \mathcal{F}^{-1}(\delta V)_p - \overline{\mathcal{F}^{-1}(\delta V)_c} \\ \frac{\lambda T_o}{N_g} \times \mathbf{O}_g \\ \lambda T_o \end{pmatrix} = \begin{pmatrix} \mathcal{F}^{-1}(\mathbf{G}(\xi_g, \eta_g)) \\ \lambda \times \mathbf{I} \\ \lambda \times \mathbf{O}_g^t \end{pmatrix} \times T_g \quad (13) $$

with I is the square identity matrix ($N_g \times N_g$) and $\mathbf{O}_g$ is a vector ($N_g \times 1$) of values equal to 1. This equation can be written as

$$ y = \mathbf{A} \times T_g \quad (14) $$

with an analytic solution in the least square sense

$$ T_g = (\mathbf{A}^t \mathbf{A})^{-1} \mathbf{A}^t \times y. \quad (15) $$

Once $T_g$ are estimated, the Sun visibilities are computed as the sum of all contributions from the subsampled sources

$$ V_{\text{Sun}} = \mathbf{G}(u, v, \xi_g, \eta_g) \times T_g(\xi_g, \eta_g). \quad (16) $$

As for the Lagrangian λ, it is chosen by plotting the L curve representing the deviation from the solution L as a function of the residual norm J for all values of the regularization parameters λ. This characteristic curve has an L-shaped dependence and the optimal value of regularization parameter can be found at the corner of the L curve. Fig. 8 shows the L curve obtained for the snapshot over the ocean in X polarization used in the previous sections. Similar curves are also obtained for all the tested snapshots and the optimal value of λ has been found to be equal to 10⁻⁶.

Fig. 9 shows the retrieved brightness temperature map $T_r$ using the multiple sources algorithm for the two snapshots of Fig. 2. We notice that the residual solar radiations are almost completely reduced and that the retrieved maps are very close to those obtained using the iterative algorithm but with a processing time slightly higher than the single-source algorithm. The difference of brightness temperature between the single-source and the multiple source Sun correction algorithms is also shown in Fig. 9. As for the iterative algorithm, we notice important differences over the Sun alias and the tails and also weak

<figure><img src="image_10.png" /><figcaption>Fig. 9. [Top] Retrieved brightness temperature $T_{r3}$ in X (left) and Y (right) polarizations after applying the iterative Sun correction algorithm. [Bottom] Difference of brightness temperature $\delta T_r = T_{r1} - T_{r3}$ between the single-source and the multiple sources Sun correction algorithms.</figcaption></figure>

oscillations elsewhere. In the next section, larger data sets of SMOS data are processed with the three algorithms in order to compare their impact on the retrieved brightness temperatures.

## VII. RESULTS AND COMPARISON

The three algorithms described above are tested and compared using SMOS interferometric measurements processed up to the level of brightness temperature. The comparison is done first over the open ocean and then over the land surfaces. To push furthermore the comparison, a global test over the entire earth surface is also shown using ten days of SMOS data.

Over open ocean, the variation of brightness temperature is expected to be quite smooth from one snapshot to another. Fig. 10 shows the spatial standard deviation in the EAFFOV, $\sigma_{T_r}(\xi, \eta)$, using the three algorithms for all snapshots within a latitude range [45°S 5°N] at boresight ($\xi = 0, \eta = 0$). The data are taken from an orbit over the Pacific Ocean in October 20, 2011, where the entire FOV is 100% ocean. We notice here the clear advantage of the iterative and multiple sources algorithms when compared to the single-source algorithm. Even though the solar radiations have not entirely been removed using the two proposed algorithms, we notice an important decrease of the residuals around the Sun alias.

A second comparison over the open ocean is to look at the latitudinal variation of $T_r$. For that, we consider two regions of pixels: 1) the circle centered at the Sun alias and of radius $d_1 = 0.0103$, and 2) the region between two circles centered at the Sun alias and of radius $d_1 = 0.0103$ and $d_2 = 0.05$. The first region shows the impact over the Sun disc whereas the second reflects the impact over the tails and the clean pixels surrounding
<figure><img src="image_11.png" /><figcaption>Fig. 10. Spatial standard deviation of 350 retrieved brightness temperature images over the ocean in X (left) and Y (right) polarizations after correcting the direct solar radiations using the: [bottom] single-source algorithm, [middle] iterative algorithm, and [top] multiple source algorithm. The data are taken from an orbit over the Pacific Ocean in January 20, 2011 and the snapshots are located within a latitude range at boresight between 45°S and 5°N.</figcaption></figure>

the Sun disc. The results are only shown here in X polarization since the impact on Y polarization is exactly the same. Fig. 11 shows the mean and standard deviation ($\overline{T_r}(\xi_d, \eta_d)$, $\sigma_{T_r}(\xi_d, \eta_d)$) for each snapshot over the two regions. From this figure, we deduce the following.

1) Over the Sun disc, the multiple sources and the iterative algorithms are close to each other and outperform the single-source algorithm. The variation of the mean $\overline{T_r}(\xi_d, \eta_d)$ and the standard deviation $\sigma_{T_r}(\xi_d, \eta_d)$ are much smoother and less noisy using the two proposed algorithms. The standard deviation is also extremely reduced.

2) Far from the Sun disc (the second region), we also notice that both the mean and standard deviation are improved using the two proposed algorithms (smoother variation, and an important decrease of the values of the standard deviation). These results show that the improvements are not restricted over the Sun disc and that the Sun tails and the clean region also benefit from these algorithms.

Moving on now to the comparison over land surfaces, Fig. 12 shows an example of the mean and standard deviation

<figure><img src="image_12.png" /><figcaption>Fig. 11. Latitudinal variation over ocean: Variation of the mean $\overline{T_r}$ (bottom) and the standard deviation $\sigma_{T_r}$ (top) of the retrieved brightness temperatures $T_r$ in X polarization using the single source (blue), the iterative (green), and the multiple sources (red) algorithms over: (a) a circle centered at the Sun alias position and for a radius $d_1 = 0.0103$ and (b) the region between two circles centered at the alias of the Sun and of radius $d_1 = 0.0103$ and $d_2 = 0.05$. The data are taken from an orbit over the Pacific Ocean in October 20, 2011 and the snapshots are located within a latitude range at boresight between 45°S and 5°N.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Latitudinal variation over land: Variation of the mean $\overline{T_r}$ (bottom) and the standard deviation $\sigma_{T_r}$ (top) of the retrieved brightness temperatures $T_r$ in X polarization using the single source (blue), the iterative (green), and the multiple sources (red) algorithms over: (a) a circle centered at the Sun alias position and for a radius $d_1 = 0.0103$ and (b) the region between two circles centered at the alias of the Sun and of radius $d_1 = 0.0103$ and $d_2 = 0.05$. The data are taken from an orbit over the USA and Mexico in January 20, 2013 and the snapshots are located within a latitude range at boresight between 25°N and 45°N.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Mean [top] and standard deviation [bottom] of $T_r$ over a disc of radius $d = 0.0103$ centered at the Sun alias inside the EAFFOV for ten days of descending orbits in X polarization after accounting for the direct solar radiations using the single source [left] and the multiple source [right] algorithms. The data were taken for the period of January 1–10, 2011. The regions in white are not contaminated by the solar radiations since the Sun is eclipsed by the earth at this period of time.</figcaption></figure>

($\overline{T_r}(\xi_d, \eta_d), \sigma_{T_r}(\xi_d, \eta_d)$) calculated over the same two regions of pixels as above for more than 150 consecutive snapshots of an orbit taken on January 20, 2013, over USA and Mexico. The results over land confirm those obtained over the ocean where the iterative and multiple sources algorithms are on the same level and outperform the single-source algorithm. These results should open the door on relaxing some of the rough filtering done at level 2 soil moisture and sea surface salinity retrievals.

Both previous tests over land and ocean were done over a limited number of snapshots taken within one orbit. Therefore, it is important to validate the proposed algorithms over larger data sets that cover the entire earth's surface. For that, we have processed ten days (January 1–10, 2011) of SMOS interferometric measurements using the multiple sources algorithm. These data are compared to those processed by the L1OP using the single-source algorithm. For each retrieved image $T_r$, we look for the Sun alias inside the EAFFOV and calculate its corresponding latitude and longitude as well as the mean and standard deviation over the disc $d_1$ described above and the results are shown in Fig. 13. We precise that the iterative algorithm, which cannot be used operationally because of its high execution time, is not considered in this test and that the results are only shown in X polarization. Looking back at Fig. 13, the following conclusions can be drawn.

1) The mean brightness temperatures suffer from important oscillations using the single-source algorithm that are

significantly reduced with the proposed algorithm leading to a much smoother variation.

2) The standard deviation map shows high residuals of the solar radiations using the single-source algorithm mainly in the northern hemisphere where the Sun is above the antenna plane at this period of time. These residuals are extremely reduced using the proposed new algorithm.

3) Over the southern hemisphere, the Sun is below the antenna plane at this period of time and the solar radiations are not accounted for in the current version of the L1OP thus the difference of contrast in the mean brightness temperature map near latitude 40°S. For the multiple source algorithm, the solar radiations are corrected efficiently even though the Sun is below the antenna plane. This map shows the need for extending the correction to both cases where the Sun is above and below the antenna plane as proposed in [15].

4) The regions in dark red over Europe, the Middle East, and China correspond to data contaminated by RFI.

5) Above 50°N of latitude, the Sun is eclipsed by the earth at this period of time and the Sun correction algorithm is not applied. Furthermore, the Sun alias is outside the EAFFOV thus explaining that there are no information in this region.

Finally, based on the results presented above, the ESA has implemented the multiple source Sun correction algorithm in
the latest version of the L1OP (v720) that will be used for the
next complete archive reprocessing and operational production.

VIII. CONCLUSION

Direct solar radiations contaminate the SMOS interferometric
measurements whether the Sun is below or above the antenna
plane of the satellite. A correction algorithm is already imple-
mented in the SMOS operational processor where the Sun disc
is considered as a single source located at its center since it is
seen under an angle much smaller than the angular resolution of
the satellite. We have shown that important residuals of the solar
radiations are still present in the retrieved brightness temperature
images mainly because the dominant source of brightness within
the Sun disc is not always located at the center. As a consequence,
the level 2 soil moisture and sea surface salinity processors are
for now filtering all pixels around the Sun position and tails in
order to reduce the impact of the residual radiations yielding to
less retrievals.

To overcome this issue, we have proposed two algorithms
in order to reduce furthermore the residual radiations. The first
one is an iterative algorithm that searches for the position of
the dominant source of brightness by minimizing the standard
deviation of the retrieved image around the Sun disc. We have
shown that this algorithm does reduce the residual radiations
but requires a very high computational time and therefore is not
deemed to be operational.

The second algorithm considers the Sun as a disc with multiple
sources defined over a fine grid. The brightness temperature of
the different sources are estimated using a constrained linear
optimization problem with an explicit solution. We have shown
that this algorithm is quite close from a performance point of
view to the iterative one and its computational time is slightly
higher than the single-source algorithm. We also indicate that the
multiple source algorithm is now implemented in the up-coming
version of the SMOS operational processor (L1OP v720) that
will be used for the next complete archive reprocessing and op-
erational production. The use of the new algorithm in the L1OP
will open up the question on reducing some of the restrictions at
level 2 leading to more retrievals and thus improving the quality
of sea surface salinity and soil moisture.

ACKNOWLEDGMENT

The authors would like to thank the CNES for funding this
work through the project TOSCA.

REFERENCES

[1] Y. H. Kerr et al., “The SMOS mission: New tool for monitoring key elements of the global water cycle,” in *Proc. IEEE*, vol. 98, no. 5, May 2010, pp. 666–687.

[2] Y. H. Kerr, P. Waldteufel, J.-P. Wigneron, J.-M. Martinuzzi, J. Font, and M. Berger, “Soil moisture retrieval from space: The soil moisture and ocean salinity (SMOS) mission,” *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 8, pp. 1729–1735, Aug. 2001.

[3] I. Corbella, N. Duffo, M. Vall-llossera, A. Camps, and F. Torres, “The visibility function in interferometric aperture synthesis radiometry,” *IEEE Trans. Geosci. Remote Sens.*, vol. 42, no. 8, pp. 1677–1682, Aug. 2004.

[4] A. Khazâal, D. J. Leroux, F. Cabot, P. Richaume, and E. Anterrieu, “Effect of the polarization leakage on the SMOS Image reconstruction algorithm: Validation using ocean model and in situ soil moisture data,” *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 9, pp. 4961–4971, Sep. 2015.

[5] A. Gutierrez et al., “SMOS L1 processor prototype: From digital counts to brightness temperatures,” in *Proc. IEEE Int. Geosci. Remote Sens. Symp.*, 2007, pp. 3626–3630.

[6] A. Camps et al., “Sun effects in 2-D aperture synthesis radiometry imaging and their cancelation,” *IEEE Trans. Geosci. Remote Sens.*, vol. 42, no. 10, pp. 1161–1167, Jun. 2004.

[7] G. A. Dulk and D. E. Gary, “The Sun at 1.4 GHz: Intensity and polarisation,” *Astron. Astrophys.*, vol. 124, pp. 103–107, 1983.

[8] M. Martín-Neira, S. Ribø, and A. J. Martín-Polegre, “Polarimetric mode of MIRAS,” *IEEE Trans. Geosci. Remote Sens.*, vol. 40, no. 8, pp. 1755–1768, Aug. 2002.

[9] M. Martín-Neira, M. Suess, J. Kainulainen, and F. Martín-Porqueras, “The flat target transformation,” *IEEE Trans. Geosci. Remote Sens.*, vol. 46, no. 3, pp. 613–620, Mar. 2008.

[10] M. Martín-Neira, Y. Menard, J. M. Goutoule, and U. Kraft, “MIRAS, a two dimensional aperture synthesis radiometer,” in *Proc. IEEE Int. Geosci. Remote Sensing Symp., Pasadena, CA, USA, 1994, vol. 3*, pp. 1323–1325.

[11] P. Waldteufel, J. Boutin, and Y. Kerr, “Selecting an optimal configuration for the soil moisture and ocean salinity mission,” *Radio Sci.*, vol. 38, no. 3, pp. 16–1–16–8, Jun. 2003.

[12] E. Anterrieu, “A resolving matrix approach for synthetic aperture imaging radiometers,” *IEEE Trans. Geosci. Remote Sens.*, vol. 42, no. 8, pp. 1649–1656, Aug. 2004.

[13] E. Anterrieu, “On the reduction of the reconstruction bias in synthetic aperture imaging radiometry,” *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 3, pp. 592–601, Mar. 2007.

[14] A. Khazâal, P. Richaume, F. Cabot, E. Anterrieu, A. Mialon, and Y. H. Kerr, “Improving the spatial bias correction algorithm in SMOS image reconstruction processor: Validation of soil moisture retrievals using in situ data,” *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 1, pp. 277–290, Jan. 2019.

[15] A. Khazâal, E. Anterrieu, F. Cabot, and Y. H. Kerr, “Impact of direct solar radiations seen by the back-lobes antenna patterns of SMOS on the retrieved images,” *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 10, no. 7, pp. 3079–3086, Jul. 2017.

[16] T. C. Van Flandern and K. F. Pulkkinen, “Low-precision formulae for planetary positions,” *Astrophys. J. Suppl. Ser.*, vol. 41, pp. 391–411, Nov. 1979.

[17] A. Camps et al., “Sun self-estimation algorithm,” Polytechnic Univ. Catalunya, Barcelona, Spain, Tech. Rep. SMOSP3-UPC-TN-0002 v 1.0, 2007.

[18] S. H. Yueh, R. West, W. J. Wilson, F. K. Li, E. G. Njoku, and Y. Rahmat-Samii, “Error sources and feasibility for microwave remote sensing of ocean salinity,” *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 5, pp. 1049–1060, May 2001.

[19] A. Khazâal, F. Cabot, Y. Soldo, and E. Anterrieu, “A Kurtosis based approach to detect RFI in SMOS image reconstruction data processor,” *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 11, pp. 7038–7047, Nov. 2014.

<figure><img src="image_15.png" /></figure>

Ali Khazâal was born in Lebanon, in 1981. He received the Engineering degree in telecommunication and computer science from the Faculty of Engineering, Lebanese University, Beirut, Lebanon, in 2003, the M.S. degree in signal, image, acoustic, and optimization from the Institut National Polytechnique de Toulouse, Toulouse, France, in 2006, and the Ph.D. degree in signal/image processing and optimization from the University Paul Sabatier, Toulouse, France, in 2009. The subject of his thesis was the Image Reconstruction Algorithms for the SMOS mission.

He was working with the Signal, Image and Instrumentation Group, Laboratoire d'Astrophysique de Toulouse-Tarbes, Observatoire Midi-Pyrénées, University Paul Sabatier. From 2009 to 2018, he was a Research Scientist with the Centre d'Études Spatiales de la BIOsphère, Toulouse, France. In 2019, he started his own company "RDIS Conseils" specialized in the Research and Development in Spatial Imaging. His research interests include numerical analysis and signal and image processing with direct application on ESA's SMOS mission and is a member of ESA's SMOS payload calibration committee. He is also working on the future microwave missions using aperture synthesis and in particulary CNES's SMOS-HR project.
<figure><img src="image_16.png" /></figure>

**François Cabot** received the Ph.D. degree in optical sciences from the University of Paris-Sud, Orsay, France, in 1995.

Between 1995 and 2004, he was with CNES wide field of view instruments quality assessment department, working on absolute and relative calibration of CNES-operated optical sensors over natural terrestrial targets. In 2004, he joined the Center for the Study of the BIOSphere, Toulouse, France, as an SMOS System Performance Engineer. His research interests include radiative transfer both optical and microwave and remote sensing of terrestrial surfaces. He has been a Principal Investigator or Coinvestigator for various calibration studies for MSG, Terra, ENVISAT, and ADEOS-II.

<figure><img src="image_17.png" /></figure>

**Eric Anterrieu** (Member, IEEE) was born in Brive, France, in 1965. He received the Engineering and M.S. degrees in solid-state physics from the Institut National des Sciences Appliquées, Toulouse, France, in 1988, and the M.S. and Ph.D. degrees in image reconstruction in astronomy from the Paul Sabatier University, Toulouse, France, in 1989 and 1992, respectively. The subject of his thesis was the image reconstruction algorithms for multiple aperture interferometry.

Since 1993, he has been an Engineer of Research in Computer Science with the Centre National de la Recherche Scientifique, Toulouse, France. He was with the Radio and Optical Aperture Synthesis Group of LATT (UMR 5572) from 1993 to 2000, with the Signal and Image Processing Team of CERFACS (URA 1875) from 2000 to 2004, and with the Signal, Image and Instrumentation Group of IRAP (UMR 5277) from 2005 to 2016. In 2017, he joined the SMOS Team at the Centre d'Etudes Spatiales de la BIOSphère, Toulouse. His research interests include numerical analysis and image and signal processing with a particular emphasis on the SMOS mission and on the future microwave missions using aperture synthesis.

<figure><img src="image_18.png" /></figure>

**Yann H. Kerr** (Fellow, IEEE) received the Engineering degree from Ecole Nationale Supérieure de l'Aéronautique et de l'Espace, Toulouse, France, in 1977, the M.Sc. degree in electrical and electronic engineering from Glasgow University, Glasgow, U.K., in 1981, and the Ph.D. degree from Université Paul Sabatier, Toulouse, France, in 1992.

From 1980 to 1985, he was with CNES. In 1985, he joined LERTS, where he was the Director in 1993–1994. He spent 19 months with the Jet Propulsion Laboratory, Pasadena, CA, USA, in 1987–1988. He has been with the Center for the Study of the BIOSphere, Toulouse, France since 1995 (Deputy Director and Director since 2007). He has been involved with many space missions. He was an EOS Principal Investigator (interdisciplinary investigations) and Principal Investigator and precursor of the use of the SCAT over land. In 1990, he started to work on the interferometric concept applied to passive microwave earth observation and was subsequently the science lead on the MIRAS project for the European Space Agency (ESA) with Matra Marconi Space and Observatoire Midi Pyrénées. He was also a Coinvestigator on IRIS, OSIRIS, and HYDROS for NASA. He was a science advisor for MIMR and Coinvestigator on AMSR. In 1997, he first proposed the natural outcome of the previous MIRAS work with what was to become the SMOS Mission, which was eventually selected by ESA in 1999 with him as the SMOS mission Lead-Investigator and Chair of the Science Advisory Group. He is also in charge of the SMOS science activities coordination in France. He has organized all the SMOS Science workshops. His research interests include the theory and techniques for microwave and thermal infrared remote sensing of the earth, with emphasis on hydrology, water resources management, and vegetation monitoring.