Multi-snapshot inversion of interferometric
synthetic aperture microwave radiometry
observations for remote sensing

Max Dunitz*

0000-0002-3244-1244

Eric Anterrieu†

0000-0002-8825-9652

Ali Khazaal†

khazaala@cesbio.cnes.fr

Nemesio Rodriguez-Fernandez†

0000-0003-3796-149X

Yann Kerr†

0000-0001-6352-1717

Jean-Michel Morel*

0000-0002-6108-897X

Miguel Colom*

0000-0003-2636-0656

Abstract—The Soil Moisture and Ocean Salinity (SMOS) satellite [3] produces images of L-band brightness temperature using interferometry. The brightness temperatures that gave rise to a set of observed visibilities are inferred by solving an inverse problem based on the Van Cittert-Zernike theorem. The images are generated from visibilities one correlation period at a time, even though a given patch of Earth’s surface may be visible to the satellite for 100 consecutive periods. The recovered brightness temperature images are noisy and, due to the antenna array’s undersampling of the u-v frequency plane, aliased. A simultaneous inversion of multiple snapshots (call it “satellite-rotation synthesis” or a “global inversion”) could take advantage of these multiple looks to improve resolution, field of view, and noise level, but appeared prohibitively expensive computationally and ill-posed: due to directional emissivity, many brightness temperatures must be recovered for each pixel. However, this directional emissivity can be sparsely modeled; accordingly, for satellites like SMOS, the observation model becomes shift-invariant across an orbital segment when expressed in geodetic coordinates relative to the trace. The observation model performs a convolution on the image of parameters (uniform in geodetic coordinates) governing observed brightness as a function of satellite position, and we can use the “Fourier transform trick” introduced in [6], to enable a global inversion.

For a future L-band mission, such as the proposed L-Band Interferometer for Fluxes and Interfaces in the Environment (LIIFE) [7], [11], we present a realistic global inversion of the observation model across the orbital trace, which simultaneously unfolds and denoises the aliased snapshots. By using sparse spline models of the brightness temperature parameters, we can ensure this inversion is accurate and stable, even with a realistic simulation that accounts for the difficulty of modeling the Van Cittert-Zernike theorem (expressed in satellite-centered direction-cosine coordinates) using geodetic coordinates.

Aliasing is a fact of life for satellite missions, whose budgetary and engineering constraints limit sensor deployment. In interferometric microwave radiometry, aliasing takes the form of a characteristic folding pattern in the recovered image. The folded zone of SMOS images, situated in the extremities, is determined by the regular spacing of antennas on the MIRAS instrument. This folding complicates the

formulation of snapshot fusion or joint inversions, for the
brightness temperatures recovered outside the alias-free field
of view (i.e., outside the hexagonal shape on the left of
Fig. 1) are superpositions of brightness temperatures at
several pixels. These corrupted images are hard to fuse.

In remote-sensing applications, inversions are typically
accomplished by left-multiplying the visibilities vector by
a (regularized or constrained) pseudo-inverse of the matrix,
often called the **G** matrix in the SMOS literature, that
expresses, for each visibility, the discretized Van Cittert-
Zernike integral over the direction cosines at which the
brightness temperature image is recovered. This makes
efficient algorithms for performing multisnapshot inversions
difficult to construct: the number of rows of the **G** matrix
(visibilitys) and columns (brightness temperature pixels<fnref n="1" />)
both scale linearly with the number of snapshots. Simulta-
neously inverting two snapshots takes eight times as long.

In [6], we introduced a “Fourier transform trick” that enables efficient multisnapshot inversion. By taking the discrete Fourier transform of each visibility along the snapshot index, we can invert an image of brightness temperature parameters orbital frequency by orbital frequency. While single-snapshot imaging involves inverting a $N_{\text{visibilities}} \times N_{\text{pixels}}$ **G** matrix and a traditional multisnapshot inversion $O(N_{\text{snapshots}} \cdot N_{\text{visibilities}}) \times O(N_{\text{snapshots}} \cdot N_{\text{pixels}})$, here we have $N_{\text{snapshots}}$ matrices to invert, each of size $N_{\text{visibilities}} \times O(\sqrt{N_{\text{pixels}}})$. These skinnier, better-conditioned matrices denoise and unfold the images across the orbit.

This article incorporates splines in the computation of
the observation model integrals, so that the global inversion
may be performed with realistic acquisition geometry and
satellite elevation. We use a parametric representation of
brightness temperatures, validated using level 1C SMOS
data in Section II, to account for directional emissivity.
Section III explains the global inversion technique and

<footnote n="1">¹In a single-snapshot inversion, the pixels are specified in direction-cosine coordinates, and then mapped to an equal-area grid. For a multi-snapshot inversion, this is not possible, as the across-track direction cosine $\xi$ of a pixel also varies as the satellite moves along-track. To get around this, the image we recover is uniform in geodetic coordinates, but we use splines to evaluate the Van Cittert-Zernike integral.</footnote>

<footnote>*Centre Borelli, ENS Paris-Saclay, Gif-sur-Yvette, France.</footnote>

<footnote>†Centre d'Etudes Spatiales de la Biosphère (CESBIO), Toulouse, France.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1: SMOS mission level 1C data. Left: a typical brightness temperature image (H polarization). The hexagonal alias-free portion of the image in direction cosines is mapped to an equal-area grid. Right: different (folding-free) measurements of brightness temperature (in H and V polarization) of a single land pixel. Both images from the file SM_OPER_MIR_SCLF1C_20210808T105109_20210808T114423_724_001_1.DBL</figcaption></figure>

demonstrates it on realistic scenarios involving SMOS and proposed LIIFE arrays.

## I. SMOS BRIGHTNESS TEMPERATURES

The SMOS instrument produces images of sky brightness in the L-band using its onboard passive instrument Microwave Imaging Radiometer for Aperture Synthesis (MIRAS) [3], which consists of a two-dimensional, Y-shaped antenna array. The instrument-plane vector spacing in wavelengths ($u_{k,l}$, $v_{k,l}$) between each antenna pair ($k, l$) is called a baseline. From each direction of arrival ($\xi, \eta$) of a planar wave from a far-field source, this baseline determines the relative lag from the source to any antenna pair whose relative spacing is the baseline. By the Van Cittert-Zernike theorem, cross-correlations $V_{k,l}$ between antenna feeds, called visibilities, are Fourier samples of the intensity of the source [5], given here as brightness temperature:

$$V_{k,l} = \iint_{\xi^2+\eta^2 \le 1} T(\xi, \eta) A_k(\xi, \eta) \overline{A_l(\xi, \eta)} e^{-2\pi i(u_{k,l}\xi + v_{k,l}\eta)} dS(\xi, \eta), \quad (1)$$

where the solid-angle differential is given by $dS(\xi, \eta) = \frac{d\xi d\eta}{\sqrt{1-\xi^2-\eta^2}}$. Here we made the simplifying assumption that the fringe-washing function is roughly unity and spatial decorrelation effects negligible, though this need not be the case. If we extend the domain of $T$ so that it assumes 0 K at non-directions (i.e., pairs of $\xi-\eta$ coordinates for which $\xi^2 + \eta^2 > 1$), (1) can be seen as a 2D Fourier transform.

Antennas spaced regularly along a Y frame induce (with redundancy) a hexagonal grid of baselines—not (only) because their convex hull is hexagonal, but because they lie on a hexagonal Bravais lattice (see Fig. 2). After taking its convex hull and rearranging, this hexagonal shape on a hexagonal lattice can be seen as a parallelogram (see Fig. 3), for which a DFT can be defined [4], [10] using a reciprocal parallelogram on a hexagonal lattice in the $\xi-\eta$ direction-cosine plane. As a compromise between hardware cost and radiometric performance and alias-free field of view (which determines swath width and thus revisit time), antennas are spaced at 0.875 wavelengths—beyond the Nyquist rate of $\frac{1}{\sqrt{3}}$, so that folding artifacts contaminate the periphery of the image. To minimize the impact of this distortion, the satellite flies with a 32° tilt with respect to nadir.<fnref n="2" />

<figure><img src="image_2.png" /><figcaption>Fig. 2: Left: antenna locations (in wavelengths of the center frequency) of the antennas on SMOS frame. Right: induced baselines, falling on hexagonal lattice points.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3: After rearrangement, the hexagonal convex hull of the SMOS baselines in the $u-v$ plane (left) forms a parallelogram on a hexagonal lattice (center). Its DFT dual grid (right), in the $\xi-\eta$ plane, extends beyond the unit circle.</figcaption></figure>

During each correlator integrator time, visibilities are computed by correlating the feeds of each antenna pair. A brightness temperature image is reconstructed using a pseudo-inverse of the matrix $G$ that discretizes the Van Cittert-Zernike-theorem-based observation model (1), mapping vectorized temperature images $T$ to the vector $V$ of visibilities measured during a correlation period. This pseudo-inverse $G^\dagger$ can be the Moore-Penrose pseudoinverse, which is the solution of minimum Euclidean norm among all solutions to the norm-approximation problem

$$\hat{T} = \arg \min_{T \in \mathbb{R}^{N_{\text{pixel}}}} \|GT - V\|_{\mathbb{R}^{N_{\text{visibilities}}}}^2, \quad (2)$$

perhaps after setting the smallest singular values to 0; a Tikhonov-regularized version thereof; or a bandlimited version thereof with a linear equality constraint that ensures the discrete Fourier transform of the recovered temperature image is annihilated at unobserved frequencies [1].

The inversion recovers brightness temperature images on a 196x196 hexagonal grid in the direction cosine plane, which are in turn mapped, based on the satellite's position and orientation, to points on an Equal Area Scalable Earth (EASE) grid [2] of Earth's surface. A level 1C brightness temperature snapshot is presented in Fig. 1. Due to aliasing and mapping geometry, only 5532 EASE cells of the 38416 hexagonal grid points are assigned valid temperatures.<fnref n="3" />

Error in discretizing (1) on the 196 × 196 direction-cosine grid is substantial, though lower than the noise level (see Fig. 4), leaving space for reparameterization of the integral.

<footnote n="2">²This tilt increases the contribution of the Sun, Moon, and sky background to the visibilities, necessitating corrections to the visibilities that remove these effects in order to ensure the invariance the global inversion approach requires.</footnote>

<footnote n="3">³3558 of the 38416 direction cosine hexagonal cells lie beyond the unit circle, and thus do not correspond to directions. The simulated brightness temperature along these pseudodirections is 0 K. The remaining 29326 cells excluded from the brightness temperature image correspond to directions that do not intersect with Earth, or those that do, but are corrupted by aliasing.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4: Left: the discretization error of each visibility computed on a representative scene with a 196 × 196 hexagonal grid, when a 1400 × 1400 quadrature grid is taken as the true value. Right: the evolution of the visibility computed of a particular (long) baseline, computed on a realistic scene with measured antenna patterns, as the grid densities by multiples of 4 sample points, from 128 × 128 to 1400 × 1400.</figcaption></figure>

## II. PARAMETRIC REPRESENTATIONS OF SMOS BRIGHTNESS TEMPERATURES

A given EASE grid cell may correspond to a recovered brightness temperature in the hexagonal direction-cosine grid dozens of times over an orbit, with dozens of Kelvin of variation. These temperatures can, in principle, vary in both of the direction-cosine coordinates at which the pixel is viewed by the satellite, but depend mostly on incidence angle, much like a Fresnel model [9], [12], [13]. When plotted against the incidence angle, these consecutive brightness temperature measurements form curves that can be approximated as the sum of a mean profile, a uniform shift in brightness temperature across all incidence angles, and an incidence-angle-dependent shift. Fig. 5 gives the mean profile and first four principal components of a PCA model fit to all pixels with fully observed data (i.e., with data in all 12 incidence-angle bins) in 8 days of L3TB data computed during the sixth reprocessing campaign. For the H polarization, the mean profile is concave and decreasing, which is characteristic of a Fresnel model,<fnref n="4" /> but brighter relative to Earth's surface temperature, due to surface roughness and multiple interfaces. The curves share a similar form but differ in mean value and decay extent; nearly all inter-pixel variance in the L3TB data set is explained by two components: a uniform shift in the brightness of the pixel and an adjustment governing the scale of its decay in incidence angle. The same analysis is performed using noisier level 1C data, after binning. The computed principal components are similar, but exhibit processing-chain artifacts (such as inexact corrections to the visibilities due to acquisition conditions). In both cases, the principal components and decline in explained-variance ratios were nearly identical for models trained on the V-polarization data. Moreover, the explained-variance ratio beyond principal component 2 is less than $10^{-10}$, signaling a high degree of regularity in these curves. Letting $R$ be the radius of a spherical approximation of Earth, $h$ the satellite's altitude, and $u_{max} = \sin^{-1}(\frac{R}{R+h})$ the incidence angle at the horizon, the first four components of this Karhunen-Loève basis of $L^2([0, u_{max}])$ nearly agree with the first four eigenfunctions of the Laplacian matrix

<figure><img src="image_5.png" /><figcaption>Fig. 5: Top: the mean profile (left) and first four principal components (right), trained on L3TB data. The mean profile is compared to a Fresnel model with a refractive index fit empirically. (Earth's surface is not a uniform planar interface and has no such global refractive index.) Bottom: the $\mathbb{R}^{12}$-normalized first four principal components of PCA fit to level 1C data (left) and $\mathbb{R}^{1200}$-normalized first four eigenvectors of the chain-graph discretization of a manifold of incidence angles (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: SMOS level 1C brightness temperature curve, as well as the best fit using 1-4 Laplace-Beltrami eigenfunctions of a discrete approximation to the compact manifold of observed incidence angles. Projection onto the first principal components of a PCA model fit to L3TB data (after subtracting the mean profile) yields similar results.</figcaption></figure>

of the chain-graph approximation to $[0, u_{max}]$, with nodes placed uniformly in the direction-cosine plane. That these two orthonormal bases of $L^2([0, u_{max}])$ nearly agree suggests the SMOS brightness temperature curves are quite smooth, since the Laplace-Beltrami eigenfunctions are sorted by Dirichlet energy (which equals the eigenvalue). Even granting that future L-band missions will produce images of higher resolution using a different instrument, reconstructed brightness temperature curves are unlikely to need more than 2-4 parameters to adequately model each pixel's change in apparent brightness as a function of the satellite's position. Indeed, the residuals between example brightness temperature curves and a weighted sum of just the first two basis functions or Laplace-Beltrami eigenfunctions can profitably be thought of as uncorrelated noise (see Fig. 7).

## III. SATELLITE-ROTATION SYNTHESIS IN GEODETIC COORDINATES RELATIVE TO THE TRACE

We model Earth's surface as a sphere of radius $R$ and the satellite's orbit as a great circle, at altitude $h$. Thanks

<footnote n="4">⁴With respect to the plane of incidence, this is the *s* polarization.</footnote>
<figure><img src="image_7.png" /></figure>

Fig. 7: Residual (left) between observed level 1C brightness temperatures of a pixel at various incidence angles and the parameterized brightness temperature curve. A Kolmogorov-Smirnov test was performed to see whether the “bag of residuals” can be profitably modeled as normal distribution of mean 0 and standard deviation equal to the sample standard deviation. In five level 1C data files, the null hypothesis was never rejected at confidence level p = 0.1. Similarly, Ljung-Box tests find no significant autocorrelation in the residual (right) at any tested lag, from 1 incidence angle bin up to half the number of incidence angle bins.

to a yaw correction, the boresight center and subsatellite point have the same relative displacement each snapshot. In these simulations, we assume the boresight center is at nadir and thus intersects Earth at the subsatellite point, though this is not necessary. The trace, that is, the curve drawn by the subsatellite point, is therefore a closed geodesic. This approximation has subpixel accuracy locally. Without loss of generality, we treat the trace as the Equator and adopt the geodetic coordinate system of [6], with $x_a$ being the across-track geodetic distance and $y_a$ the along-track geodetic coordinate relative to a fixed meridian.

As the subsatellite point traces a closed geodesic on the spherical surface $S$, in the $x_a-y_a$ plane $\mathbb{R}^2$, the set of points visible to the satellite forms the band $B = [-d, d] \times \mathbb{R}$, where $d = \cos^{-1}(\frac{R}{R+h})$ is the horizon in geodetic coordinates. We can choose the reference meridian so that, in geodetic coordinates, the subsatellite point is at $(0, y_s(t))$ in the $t$th snapshot, where $y_s(t) = \Delta y_a t$, and $\Delta y_a$ is the geodetic distance the subsatellite point moves in a correlation period.

Letting $r(t)$ be the satellite's position with respect to an inertial frame whose origin is the center of a nonrotating, spherical Earth, the direction of motion $y_{SAT}(t) = r'(t)/||r'(t)||$, negative angular momentum vector $x_{SAT}(t) = -(r(t) \times r'(t))/||r(t) \times r'(t)||$, and anticorrelated vector $z_{SAT}(t) = y(t) \times x(t) = -r(t)/||r(t)||$ define a right-handed, satellite-centered coordinate system, from which spherical coordinates $(\rho, \theta, \phi)$ and direction cosines $(\xi, \eta)$ are defined. The subsatellite point SSP, satellite SAT, and visible patch of Earth OBS define a triangle that, along with its projection on the equatorial plane (shaded in Fig. 8), defines conversions between the two systems:

$$ \rho = (R+h)\sqrt{1-\xi^2-\eta^2} - \sqrt{R^2-(\xi^2+\eta^2)(R+h)^2} \quad (\text{law of cosines}); $$

$$ y_a = \sin^{-1} \left( \frac{\rho(\xi, \eta)\eta}{\sqrt{R^2 - \rho^2(\xi, \eta)\xi^2}} \right); \quad x_a = \sin^{-1} \left( \frac{-\rho(\xi, \eta)\xi}{R} \right). $$

### A. The observation model in semi-geodetic coordinates

Accurately discretizing the observation model in the $x_a-y_a$ plane requires many sample points and results in unwieldy (> 1 TB) and badly conditioned G matrices. To preserve shift-invariance, so that each visibility signal across an orbit is computed via a convolution of the image of

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

Fig. 8: Left: visual aid for conversions between satellite-centered direction cosines $(\xi, \eta)$ and spherical coordinates $(\rho, \phi, \theta)$ and geodetic coordinates relative to the trace $(x_a, y_a)$. Right: simulated trajectory of satellite superimposed on image of brightness temperature parameters $\alpha$ (the RGB values of NASA's Blue Marble image).

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

Fig. 9: Left: the simulated data, uniform in $x_a-y_a$ coordinates. Right: the same, plotted in the $\xi-y_a$ plane. The Riemann sum grid is uniform in $\xi-y_a$ coordinates, negotiating the tradeoff between geodetic coordinates (which offer interpretability and invariance) and satellite-centered direction-cosine coordinate (which offer smaller Riemann sum grids for accurate visibilities and more stable inversions).

brightness temperature parameters and fixed function, we discretize (1) in $\xi-y_a$ coordinates:

$$ V_{k,l} = \iint_{S_d(y_s)} T(\xi, y_a) A_k(\xi, y_a) A_l(\xi, y_a) e^{-2\pi i (u_k, l\xi + v_k, l\eta)} dS(\xi, y_a), $$

where the $\eta$-coordinate of a pixel with given $\xi$-coordinate and $y_a$ coordinate (relative to the subsatellite point $y_s$ is

$$ \eta(\xi, y_a) = \sin(y_a) \sqrt{\frac{R^2}{\rho^2(\xi, y_a)} - \xi^2}; $$

the visible disk $S_d(y_s)$ is the set of coordinates $(\xi, y_a) \in [-1, 1] \times [y_s - d, y_s + d]$ that satisfies

$$ |\xi| \le \frac{R}{(R+h) \cos(y_a - y_s)} \sqrt{\frac{(R+h)^2 \cos^2(y_a - y_s) - R^2}{(R+h)^2 - R^2}}; \quad (3) $$

and the solid-angle differential can be found using

$$ dS(\xi, y_a) = \frac{\partial y_a}{\partial \eta} dS(\xi, \eta) = \frac{\partial y_a}{\partial \eta} \frac{d\xi}{\sqrt{1-\xi^2-\eta^2}}; \quad \text{where} $$

$$ \frac{\partial y_a}{\partial \eta} = \frac{\rho + \eta \frac{\partial \rho}{\partial \eta} \frac{R^2}{R^2 - \rho^2 \xi^2}}{\sqrt{R^2 - \rho^2(\xi^2 + \eta^2)}}, \quad \frac{\partial \rho}{\partial \eta} = \frac{\frac{(R+h)\rho\eta}{\sqrt{1-\xi^2-\eta^2}}}{(R+h)\sqrt{1-\xi^2-\eta^2-\rho^2}}. $$

Due to sampling irregularity<fnref n="5" /> and the multitude of $\xi$ values at which a geodetic pixel is observed over an orbit, we use splines to express the Riemann sum in $\xi-y_a$ coordinates.

<footnote n="5">⁵As Earth bends away from the satellite toward the horizon, an interval in the $\xi$ coordinate coincides with a larger patch of surface; see Fig. 9.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 10: Top: parameter image α reproduced from columnwise spline fits. Bottom: error image between the first parameter and its columnwise spline fit resampled on the planisphere (mean absolute error: 0.020, for the red channel varying between 0 and 1). In both images, the 20 knots— an unusually small number, for illustration purposes—are indicated in blue on the right.</figcaption></figure>

## B. Parametric Representation of Brightness Temperature

We define our parametric representation of brightness temperature as follows:

$$T(\xi, y_a) = T_{mean}(\xi, y_a) + \sum_{c=0}^{C-1} I_c(\xi, y_a) \sigma_c T_c(\xi, y_a), \quad (4)$$

where the basis functions $T_c$ sparsely capture each point's contribution to the irradiance profile and are, along with the (unknown) brightness temperature parameters (the $I_c$ images), expressed in $\xi-y_a$ coordinates. In practice, they depend on $(\xi, y_a)$ only via the incidence angle $u$:

$$u(\xi, y_a) = \cos^{-1} \left( \frac{(R+h)^2 + \rho(\xi, y_a)^2 - R^2}{2(R+h)\rho(\xi, y_a)} \right).$$

## C. Spline Model on the Parametric Representation

We begin with an $M \times N \times C$ equirectangular image $I$ (e.g., Fig. (8), right) of the $C$ brightness temperature parameters. As this image is uniform in geodetic coordinates, we fit a spline model to each fixed-$y_a$ column, using knots spaced so that an equal number of Riemann sum grid points (green exes in Fig. 11), over a snapshot, reposes between each knot; visibility accuracy is not sensitive to the periphery. We denote the $L \times N \times C$ image of $L$-parameter columnwise spline models as $\alpha$. An equirectangular image resampled from such a model is sharp in the center and blurry near the horizon, as in Fig. 10. Natural cubic splines, B-splines, and Gaussian processes all work, with well-selected knots to ensure accurate visibilities in simulation and stable matrices in inversion. Since the splines are columnwise, adjacent columns may exhibit disagreement.<fnref n="6" /> Given any spline model, let $B_j$ be the matrix that maps the spline weights for any brightness temperature parameter channel $c \in \mathbb{Z}/C\mathbb{Z}$ to the evaluated brightness temperature parameters in the column with offset index $j$ relative to the subsatellite point. In the $t$th snapshot, the subsatellite point is in column $t$, and the brightness temperature parameters of column $t+j$ at the column's Riemann sum sample points (green exes in the $j$th column of Fig. 11) are $B_j\alpha[\cdot; t+j, c]$.

<figure><img src="image_11.png" /><figcaption>Fig. 11: Left: Riemann sum grid (green), uniform in $\xi-y_a$ coordinates, along with brightness temperature parameters, uniform in geodetic coordinates relative to the trace. The 20 spline knots given in blue. Right: spline fits to the (a priori, unknown) brightness temperature parameters in a given column.</figcaption></figure>

## D. The Fourier Transform Trick

Let $S_{\xi,\eta}^d$ be the intersection of the disk $S_d(y_s)$ and the odd-parity Riemann sum grid $[-|\frac{R}{(R+h)\Delta\xi}|, -|\frac{R}{(R+h)\Delta\xi}| + \Delta\xi, \dots, |\frac{R}{(R+h)\Delta\xi}|] \times \Delta y_a$ $[[-\frac{N y}{2}, \dots, \frac{N y}{2}]$. We define

$$Q_{k,l}^c(\xi, y_a) = \sigma_c T_c(\xi, y_a) \cdot A_k(\xi, y_a) \overline{A_l(\xi, y_a)} e^{-2\pi i (u_{k,l}\xi + v_{k,l}\eta(\xi, y_a))}, \\ \frac{\partial y_a(\xi, y_a)}{\partial \eta} / \sqrt{1 - \xi^2 - \eta(\xi, y_a)^2 \Delta\xi \Delta y_a}.$$

The term $Q_{k,l}^c(\xi, y_a)$ does not depend on what is being observed; it merely describes how the instrument is observing. It does not depend on the satellite's position and changes only when the antenna patterns are re-estimated during calibration. This allows us to write our observation model as a convolution between our image of brightness temperature parameters and this function. The terms that constitute the product $Q_{k,l}^c(\xi, y_a)$ can either be computed directly on the choice of locations of sample points (as is the case of the Jacobian or the complex exponential) or interpolated onto these points (as is the case with the well-sampled, smooth antenna radiation patterns $A_l$).<fnref n="7" />

We subtract out the contribution of the mean temperature profile $T_{mean}$ to the visibilities; the integral giving the centered visibility $\tilde{V}_{k,l}$ associated with antennas $k$ and $l$ is discretized as follows:

$$\tilde{V}_{k,l}(y_s) = \sum_{(\xi, y_a) \in S_{\xi, y_a}^d} \sum_{c=0}^{C-1} Q_{k,l}^c(\xi, y_a) I_c(\xi, y_a + y_s). \quad (5)$$

where $I_c(\xi, y_a + y_s)$ gives the (spline-modeled) weight on the $c$th brightness temperature basis function for the patch of area $\Delta\xi\Delta y_a$ centered at $(\xi, y_a + y_s)$. Taking the partial Fourier transform along the snapshot number $t$, letting $N_x^j$ give the number of Riemann sum sample points in column $j$ relative to the subsatellite point, whose sorted and vectorized

<footnote n="6">⁶Two-dimensional splines, such as thin-plate splines [14], can in principle be used, as can the Riesz representers of the bounded, linear observation functionals (1), though the latter involves intractable computations to form the $G$ matrices and the former poses (perhaps surmountable) conditioning issues.</footnote>

<footnote n="7">⁷For long baselines, the complex exponential term varies the fastest, and for more accuracy (or a smaller Riemann sum grid), we replace it and the Jacobian with their average value over the Riemann sum cell of area $\Delta\xi\Delta y_a$, minding the singularity in the Jacobian at the horizon.</footnote>
ξ coordinates are ξ<sub>j</sub>, we write

$$
\begin{align*}
\tilde{V}_{k,l}(\omega) &= \sum_{t=0}^{T-1} \tilde{V}_{k,l}(y_s(t)) e^{-2\pi i t \omega / T} \\
&= \sum_{c=0}^{C-1} \sum_{j=-\left\lfloor \frac{N_y}{2} \right\rfloor}^{\left\lfloor \frac{N_y}{2} \right\rfloor} \sum_{n=0}^{N_j^C-1} Q_{k,l}^C(\xi_j[n], j\Delta y_a) \\
&\qquad \underbrace{\mathbf{B}_j \left( \sum_{t=0}^{T-1} \alpha[n, j+t \pmod{N}, c] e^{-2\pi i t \omega / T} \right)}_{(\mathbf{B}_j \hat{\alpha}_{\text{partial}}[:; \omega, c])[:n] e^{2\pi i j \omega / T}} \\
&= \sum_{c=0}^{C-1} \sum_{j=-\left\lfloor \frac{N_y}{2} \right\rfloor}^{\left\lfloor \frac{N_y}{2} \right\rfloor} Q_{k,l}^C(j)^T \mathbf{B}_j \hat{\alpha}_{\text{partial}}[:; \omega, c] e^{2\pi i j \omega / T},
\end{align*}
$$

where $Q_{k,l}^c(j)^T$ is the length-$N_x^j$ row vector
$(Q_{k,l}^c(\xi_j[0], j\Delta y_a), \dots, Q_{k,l}^c(x_j[N_x^j-1], j\Delta y_a))$. In the
second line of the derivation above, the function

$$
\hat{\alpha}_{\text{partial}} \stackrel{\text{def}}{=} \sum_{t=0}^{T-1} \alpha[n, j+t \pmod{N}, c] e^{-2\pi i t \omega / T}
$$

is the discrete Fourier transform of channel *c* of the image
*α* along row *n* if *N* = *T*. That this does not hold in general
is the source of apparitions at the boundary of an image.

Our unknowns are now the spline basis function weights
α[:; ;, c], after taking the *T*-point DFT. We can recover it one
frequency ω ∈ Z/TZ at a time. For each ω ∈ Z/NZ, let
I(ω) be the length-LC column vector formed by stacking
the unknown spline weights on the brightness temperature
parameters at frequency ω

$$
I(\omega) = (\hat{\alpha}_{\text{partial}}[:; \omega, 0]^T, \dots, \hat{\alpha}_{\text{partial}}[:; \omega, C-1]^T)^T \quad (6)
$$

and let $V(\omega)$ be the length-$N_{\text{baseline}}$ vector of transformed visibilities at orbital frequency $\omega$:

$$
V(\omega) = \left( \hat{\bar{V}}_{0,0}(\omega), \dots, \hat{\bar{V}}_{0,N_{\text{ant}}-1}(\omega), \dots, \hat{\bar{V}}_{N_{\text{ant}}-1,0}(\omega), \dots, \hat{\bar{V}}_{N_{\text{ant}}-1,N_{\text{ant}}-1}(\omega) \right)^T
$$

The $kN_{ant} + l$th row of the $N_{baseline} \times$ LC matrix $\mathbf{G}(\omega)$ can be written

$$
\mathbf{G}(\omega)[kN_{ant} + l, \cdot] =
\begin{pmatrix}
\sum_{j=-\lfloor N_y/2 \rfloor}^{\lfloor N_y/2 \rfloor} (Q_{k,l}^0 (j)^T \mathbf{B}_j), & \dots, & \sum_{j=-\lfloor N_y/2 \rfloor}^{\lfloor N_y/2 \rfloor} (Q_{k,l}^{C-1} (j)^T \mathbf{B}_j) \\
\vdots & \ddots & \vdots
\end{pmatrix}
$$

We have partitioned row $kN_{ant} + l$ of $\mathbf{G}$, of length $CL$, into $C$ length-$L$ vectors; the $c$th such vector is $\sum_{j=-\lfloor N_y/2 \rfloor}^{\lfloor N_y/2 \rfloor} (Q_{k,l}^c(j)^T \mathbf{B}_j)$. The spline weights on the brightness temperature parameters at frequency $\omega$ can be found by solving $V(\omega) = \mathbf{G}(\omega)T(\omega)$ for each $\omega \in \mathbb{Z}/TZ$. $\mathbf{G}$ matrices have a clear singular-value cutoff, but the singular value computed on the nonzero singular values can vary over a single global inversion from $10^2$ to $10^7$, depending on the orbital frequency $\omega$. Without noise, constructions using a fixed Tikhonov parameter have slight imperfections, as in Fig. 12. We add i.i.d. complex Gaussian noise with variance $\frac{V_{DC}^2}{2B\tau}$ to visibilities simulated$^8$ using (5) and invert

<figure><img src="image_12.png" /><figcaption>Fig. 12: Left: error in brightness temperature sampled at 22.5° incidence angle of the brightness spline weights recovered from SMOS visibilities, resampled on the planisphere, using Tikhonov parameter $10^{-3}$ for each orbital frequency $\omega$. Right: error between simulated and reconstructed brightness temperatures at three distinct incidence angles, using SMOS instrument. The swath is wider and error in K substantially lower than single-snapshot inversions performed on a hexagonal grid.</figcaption></figure>

them to recover an image of spline weights, from which
the brightness temperature parameters can be reconstructed.
The brightness temperature error image of a sample orbit is
given in Fig. 12 (right).

IV. CONCLUSION

The multisnapshot inversion technique takes advantage
of the coordinate invariance and measurement redundancy
to unfold and denoise interferometric images, improving
satellite swath width and is well-suited to interferometric
processing of irregular antenna arrays [8]. While enabling
new orbits and instrument designs, the technique imposes
the constraint of look-invariance.

REFERENCES

[1] E. Anterrieu, “A resolving matrix approach for synthetic aperture imaging radiometers,” *IEEE Trans Geosci Remote Sens*, 42(8), 2004.

[2] M.J. Brodzik, B. Billingsley, T. Haran, B. Raup, and M.H. Savoie, “EASE-Grid 2.0: Incremental but significant improvements for Earth-gridded data sets,” *ISPRS Int J Geo-Inf*, 1(1), 32-45, 2012.

[3] H.M. Barré, B. Duesmann, and Y.H. Kerr, “SMOS: The mission and the system,” *IEEE Trans Geosci Remote Sens* 46(3), 587-593, 2008.

[4] A. Camps, J. Bará, I.C. Sanahuja, and F. Torres, “The processing of hexagonally sampled signals with standard rectangular techniques: Application to 2-D large aperture synthesis interferometric radiometers.” *IEEE Trans Geosci Remote Sens* 35(1), 183-190, 1997.

[5] I. Corbella et al. “The visibility function in interferometric aperture synthesis radiometry.” *IEEE Trans Geosci Remote Sens* 42(8), 2004.

[6] M. Dunitz et al., “An alternative concept for SMOS-HR: unfolding the brightness temperature map by along-the-track inversion of the Van Cittert-Zernike equation,” *IEEE Conf Antenna Meas Appl CAMA 2021*, 414-419.

[7] Y. Kerr et al. “The SMOS mission: New tool for monitoring key elements of the global water cycle.” *Proc IEEE* 98(5):666-687, 2010.

[8] P. Krzakala et al. “Irregular Layout for a Satellite’s Interferometric Array.” *IEEE J Set Top Appl Earth Obs Remote Sens*, 14, 2021.

[9] M. Leduc-Leballeur et al., “Modeling L-Band Brightness Temperature at Dome C in Antarctica and Comparison With SMOS Observations,” *IEEE Trans Geosci Remote Sens* 53(7), 4022-4032, 2015.

[10] R.M. Mersereau, “The processing of hexagonally sampled two-dimensional signals.” *Proc IEEE* 67(6), 930-949, 1979.

[11] M. Martin-Neira et al., “TriHex: Combining formation flying, general circular orbits and alias-free imaging, for high resolution L-band aperture synthesis,” *IEEE Trans Geosci and Remote Sens*, 2023.

[12] J.F. Paris. (1969). “Microwave radiometry and its application to marine meteorology and oceanography” (No. NASA-CR-103237).

[13] T. Pellarin et al., “Two-year global simulation of L-band brightness temperatures over land,” *IEEE Trans Geosci Remote Sens* 41(9), 2003.

[14] G. Wahba, “Spline interpolation and smoothing on the sphere.” *SIAM J Scien Stat Comp* 2(1), 5-16, 1981.

<footnote>$^8V_{DC}$ is the mean noise-free visibility associated with the baseline (0, 0), and the bandwidth-correlation time product $B\tau$ is $2e7$ for SMOS. The noise associated with the baseline $(-u, -v)$ is the complex conjugate of the noise associated with $(u, v)$.</footnote>