Fourier-Correlation Imaging

TL;DR

This paper explores how correlating Fourier components from satellite-based electric field measurements can enable two-dimensional imaging of Earth's surface temperature with about one kilometer resolution, highlighting limitations and potential improvements.

Contribution

It introduces a method using Fourier component correlation for satellite imaging and analyzes the resolution limits with two antennas, proposing ways to enhance radiometric resolution.

Findings

Spatial resolution of about 1 km with two antennas at 100m separation.

Radiometric resolution improves with more antennas, but requires averaging multiple profiles.

Two antennas alone are insufficient for precise temperature field reconstruction.

Abstract

We investigate to what extent correlating the Fourier components at slightly shifted frequencies of the fluctuations of the electric field measured with a one-dimensional antenna array on board of a satellite flying over a plane, allows one to measure the two-dimensional brilliance temperature as function of position in the plane. We find that the achievable spatial resolution resulting from just two antennas is of the order of $h\chi$, with $\chi=c/(\Delta r \omega_0)$, both in the direction of flight of the satellite and in the direction perpendicular to it, where $\Delta r$ is the distance between the antennas, $\omega_0$ the central frequency, $h$ the height of the satellite over the plane, and $c$ the speed of light. Two antennas separated by a distance of about 100m on a satellite flying with a speed of a few km/s at a height of the order of 1000km and a central frequency of order GHz allow therefore the imaging of the brilliance temperature on the surface of Earth with a resolution of the order of one km. For a single point source, the relative radiometric resolution is of order $\sqrt{\chi}$, but for a uniform temperature field in a half plane left or right of the satellite track it is only of order $1/\chi^{3/2}$, indicating that two antennas do not suffice for a precise reconstruction of the temperature field. Several ideas are discussed how the radiometric resolution could be enhanced. In particular, having $N$ antennas all separated by at least a distance of the order of the wave-length, allows one to increase the signal-to-noise ratio by a factor of order $N$, but requires to average over $N^2$ temperature profiles obtained from as many pairs of antennas.

Figures

Equations

$${\bf E}({\bf r}_{1}+{\bf v}_{s}t,t)=-\frac{\mu_{0}}{4\pi}\int d^{3}r^{\prime\prime}\frac{1}{R(t)}\partial_{t^{\prime}}{\bf j}({\bf r}^{\prime\prime},t^{\prime})\Big{|}_{t^{\prime}=t-R(t)/c}\,,$$

$$\mathbf{j}({\mathbf{r}}^{\prime \prime },t)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }d{\omega }^{\prime }{e}^{i{\omega }^{\prime }t}\tilde{\mathbf{j}}({\mathbf{r}}^{\prime \prime },{\omega }^{\prime }),$$

$${\mathbf{E}}_{{\mathbf{r}}_1}(t)\equiv \mathbf{E}({\mathbf{r}}_1+{\mathbf{v}}_{s}t,t)$$

$${\tilde{\mathbf{E}}}_{{\mathbf{r}}_1}({\omega }_1)$$

$$=$$

$$⟨{\tilde{j}}_i({\mathbf{r}}_1^{\prime \prime },{\omega }_1){\tilde{j}}_j^{\ast }({\mathbf{r}}_2^{\prime \prime },{\omega }_2)⟩={\delta }_{ij}\frac{l_c^3}{{\tau }_c}\delta ({\mathbf{r}}_1^{\prime \prime }-{\mathbf{r}}_2^{\prime \prime })\delta ({\omega }_1-{\omega }_2)⟨|{\tilde{j}}_i({\mathbf{r}}_2^{\prime \prime },{\omega }_2){|}^2⟩,$$

$$⟨|\widetilde{j_i}(x,y,\omega ){|}^2⟩=K_{2}T(x,y),$$

$$⟨{\tilde{j}}_i({\mathbf{r}}_1^{\prime \prime },{\omega }_1){\tilde{j}}_j^{\ast }({\mathbf{r}}_2^{\prime \prime },{\omega }_2)⟩={\delta }_{ij}{K}_3\delta ({\mathbf{r}}_1^{\prime \prime }-{\mathbf{r}}_2^{\prime \prime })\delta ({\omega }_1-{\omega }_2)T(x^{\prime \prime },y^{\prime \prime })\delta (z^{\prime \prime }),$$

$$C_{ij}^F({\mathbf{r}}_1,{\mathbf{r}}_2,{\omega }_1,{\omega }_2)$$

$$C_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2,{\omega }_1,{\omega }_2)$$

$$i\Phi =i\left[\tau \left(-{\omega }^{\prime }+\frac{{\omega }_2+{\omega }_1}{2}\right)+t({\omega }_2-{\omega }_1)-\frac{{\omega }^{\prime }}{c}\left(|{\mathbf{r}}_1-{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|-|{\mathbf{r}}_2+{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|\right)\right].$$

$$\int T(x^{\prime }+v_{s}t,y^{\prime })e^{i({\omega }_2-{\omega }_1)t}dt$$

$$\equiv$$

$${\tilde{T}}_{x^{\prime },y^{\prime }}({\kappa }_x)$$

$$=$$

$$C_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2,{\omega }_1,{\omega }_2)$$

$${\int }_{-\infty }^{\infty }\exp \left[\dots \right]d{\omega }^{\prime }=2\pi \delta (\tau +\frac{1}{c}(|{\mathbf{r}}_1-{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|-|{\mathbf{r}}_2+{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|))e^{i\tau \frac{{\omega }_1+{\omega }_2}{2}}.$$

$$|{\mathbf{r}}_1-{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|-|{\mathbf{r}}_2+{\mathbf{v}}_s\tau /2-{\mathbf{r}}^{\prime }|\simeq -{\hat{e}}_{\mathbf{R}-{\mathbf{r}}^{\prime }}\cdot (\Delta \mathbf{r}+{\mathbf{v}}_s\tau ).$$

$$C_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2,{\omega }_1,{\omega }_2)$$

$$C_{ij}({\mathbf{r}}_1,{\mathbf{r}}_1+\Delta r{\hat{e}}_y,{\omega }_1,{\omega }_2)=K_6{\delta }_{ij}{e}^{-i{\kappa }_x{x}_0}{\int }_{-\infty }^{\infty }\frac{d\eta }{\sqrt{{\eta }^2+{\tilde{h}}^2}}K({\kappa }_x\Delta r\sqrt{{\eta }^2+{\tilde{h}}^2},\frac{\Delta r{\omega }_c}{c}\frac{\eta }{\sqrt{{\eta }^2+{\tilde{h}}^2}})\tilde{T}({\kappa }_x,\eta ),$$

$$K(\alpha ,\beta )={\int }_{-\infty }^{\infty }d\xi \frac{e^{-i(\alpha \xi +\frac{\beta }{\sqrt{{\xi }^2+1}})}}{{\xi }^2+1},$$

$$\alpha$$

$$\beta$$

$$K(\alpha ,\beta )$$

$$K(\alpha ,0)=\pi e^{-\alpha }.$$

$$K(\alpha ,\beta )$$

$$=$$

$$K(\alpha ,\beta )\simeq \sqrt{\frac{2\pi }{\beta }}{e}^{i\pi /4}{e}^{-i\beta },$$

$$K(0,\beta )={\int }_{-\infty }^{\infty }d\xi \frac{e^{-i\frac{\beta }{\sqrt{{\xi }^2+1}}}}{{\xi }^2+1}=\pi (J_0(\beta )-i{H}_0(\beta )),$$

$$\frac{\eta }{{\eta }^2+{\tilde{h}}^2}\gg \frac{c{\kappa }_x}{{\omega }_c}=\frac{c}{v_s}\frac{\Delta \omega }{{\omega }_c},$$

$${\eta }_{1,2}=\frac{1\pm \sqrt{1-4{\delta }^2{\tilde{h}}^2}}{2\delta }\equiv {\eta }_{1,2}(\kappa ),$$

$$K(\alpha ,\beta )\simeq \sqrt{\frac{2\pi }{|\beta |}}{e}^{\mathrm{sign}(\beta )i\frac{\pi }{4}}{e}^{-i\beta }$$

$$C_{ii}({\mathbf{r}}_1,{\mathbf{r}}_1+\Delta r{\hat{e}}_y,\kappa ,{\tilde{k}}_c)$$

Full text

I Introduction

Spatial aperture synthesis is a standard technique in radio-astronomy Interferometry_Synthesis_Radio_Astronomy . It allows one to achieve the fine resolution of a large antenna by correlating time-delayed signals received from the different antennas in an antenna array. In satellite-based remote sensing, spatial aperture synthesis is a technique of choice when relatively long wave-lengths are imposed by the applications, such as the measurement of sea surface salinity or surface soil moisture. When operating in the protected L-band (1400-1427 MHz), a resolution of 10km would require already a single antenna with a size of 32 meters. Spatial aperture synthesis for passive microwave-sensing was therefore proposed to ESA Kerr2001 , and implemented for the first time in the SMOS mission in 2009 that still operates today Kerr2001b ; Kerr2010 . The satellite uses a deployable Y-shaped antenna array and provides a spatial resolution between 27-60km.

With the application-driven need for higher spatial resolution down to the order of 1km, even spatial aperture synthesis leads to forbiddingly large antenna arrays, and there is therefore an ongoing quest for finding alternative concepts (see e.g. camps_two-dimensional_2001 and references therein). Compared to stationary antenna arrays on Earth used for astronomy, one may wonder whether the motion of the satellite could be used for creating a two-dimensional (2D) artificial antenna array out of a one-dimensional (1D) moving array, oriented perpendicular to the motion of the satellite. It turns out that this is not possible when directly correlating the observed microwave fields in the time-domain: the useful phase-shift gained due to the motion of the satellite is, to first order in $v_{s}/c$ cancelled by the Doppler shift, where $v_{s}$ is the speed of the satellite and $c$ the speed of light braun_generalization_2016 .

In this paper, we examine another idea: instead of correlating the signals in the time domain, we consider the correlations between their Fourier components at slightly different frequencies. This may appear surprising at first, as, at the level of the sources, the standard model assumption is that different frequencies are entirely decorrelated. Nevertheless, a hypothetical monochromatic point source is seen by different antennas at slightly different frequencies due to the slightly different Doppler effect, and hence it makes sense to correlate different frequency components from different antennas with each other. The useful frequency differences are tiny, down to below one Hertz, and correspondingly long acquisition times are needed. However, one may hope that this opens at least in principal a new way of achieving a resolution of the order of a kilo-meter in passive microwave remote sensing in the L-band by using the motion of the satellite for reducing a 2D antenna-array to a 1D array.

We derive the principles of this “Fourier-correlation imaging” (FouCoIm) technique in detail, and calculate the achievable spatial and radiometric resolution. An emphasis is put on pushing analytical calculations as far as possible, and testing the method at the hand of simple situations, namely a single point source and a uniform temperature field. Estimation of numerical values will be done with a standard set of parameters: $h=700$ km, $v_{s}=7$ km/s, $\omega_{0}=2\pi\times 1.4$ GHz, $T=300$ K, $B=20$ MHz, and $\Delta r=100$ m. This leads to the important dimensionless parameters $\beta_{s}=v_{s}/c=2.33\cdot 10^{-5}$, $\chi=c/(\Delta r\omega_{0})=3.41\cdot 10^{-4}$, and $\tilde{h}\equiv h/\Delta r=7000$.

II Model

We assume that the fluctuating micro-wave fields measured at the position of the satellite are created by fluctuating microscopic electrical currents at the surface of Earth that are in local thermal equilibrium at absolute temperature $T(x,y)$, where $x,y$ are coordinates of a point on the surface of Earth. The entire analysis will be in terms of classical electro-dynamics. In braun_generalization_2016 we derived the expression

[TABLE]

for the time dependent electric field arising from the current fluctuations at the position of the satellite, with $R(t)=|{\bf r}_{1}+{\bf v}_{s}t-{\bf r}^{\prime\prime}|$, where ${\bf r}_{1}$ is the position of the antenna at time $t=0$, ${\bf v}_{s}$ the speed of the satellite in the Earth-fixed reference frame, $\mu_{0}$ the magnetic permeability of vacuum, and ${\bf j}({\bf r}^{\prime\prime},t)$ the current density as function of space and time. All expressions are in the Earth-fixed reference frame, which is more convenient for the present study than the satellite-fixed reference frame. It was shown in braun_generalization_2016 that (1) is the correct far-field up to relativistic corrections of the prefactor of order $\beta_{s}$ (due to the mixing of electric and magnetic fields in a moving reference frame), and the neglect of terms of order $\beta_{s}^{2}$ in the phase. Eq.(1) does contain in the phase the linear Doppler shift and relativistic effects (including time dilation) up to order $\beta_{s}$. The far-field approximation is justified for $R(t)\gg\lambda$, where $\lambda$ (of order cm in the micro-wave regime) is the wave-length of the radiation (see Chapt.9 in Jackson99 ).

We substitute the Fourier-decomposition of the current density,

[TABLE]

into (1). The question whether one should differentiate $R(t)$ with respect to $t$ was answered to the negative in braun_generalization_2016 , but it is irrelevant if we neglect changes of order $\beta_{s}$ to the prefactor. We then find the time dependent field seen by the flying antenna,

[TABLE]

with $K_{1}=-\mu_{0}/(4\pi)$. The Fourier transform of that signal is

[TABLE]

We assume that the current sources can be described by a Gaussian process, where sources at different positions or different frequencies, or with different polarizations are uncorrelated,

[TABLE]

where we have introduced for dimensional grounds a correlation length $l_{c}$ and a correlation time $\tau_{c}$, and the polarizations are indexed by $i,j$, taking values $x,y,z$. In principle the average $\langle\ldots\rangle$ is over an ensemble of realizations of the stochastic process, but we may assume ergodicity of the fluctuations, such that they can also be obtained from a sufficiently long temporal average. In practice this means that one should average over positions considered as equivalent in terms of the ensemble, i.e. the time the satellite takes to fly over a desired pixel size. For a satellite flying at a speed of order km/s and a pixel size of order km, this means a maximal averaging time of the order of a second. This does not preclude calculating Fourier transforms with finer spectral resolution from data acquired over much longer times.

We will make the assumption that only the current intensities at the surface of Earth contribute. In reality the emission seen by the satellite arises from a thin surface layer on Earth that has a finite thickness $d$ of the order of a few centimeters Kerr2010 ; Kerr2012 , depending on the soil and its humidity, and the satellite also sees the cosmic microwave background. We approximate the surface layer as a single plane located at $z^{\prime\prime}=0$, i.e. $\langle|\tilde{j}_{i}({\bf r}_{2}^{\prime\prime},\omega_{2})|^{2}\rangle=d\langle|\tilde{j}_{i}(x^{\prime\prime},y^{\prime\prime},\omega_{2})|^{2}\rangle\delta(z^{\prime\prime})$ and neglect the cosmic microwave background as its temperature is two orders of magnitude lower than that of Earth, as well as other astronomical objects.

The current intensities are related to an effective temperature $T(x,y)$ by

[TABLE]

where $K_{2}$ is a constant (see eq.(105)). Eq.(7) is valid for $\hbar\omega\ll k_{B}T$ and hence very well adapted to micro-wave emission at room temperature.

Eq.(6) together with (7) is a standard model of classical white noise currents, and appears in many places in the literature, see e.g. eq.(4.16) in sharkov_passive_2003 . The equation is an instance of the fluctuation-dissipation theorem that can be found in standard text-books on statistical physics (see e.g. Part 1, Chap. XII and Part 2, Chap.VIII. inLLStatMech80 ). In the context of thermal radiation it goes back at least to the original Russian version of rytov_theory_1959 (from 1953); see also Rytov89 . The model has also been used to study coherence effects in the thermal radiation of near-fields (see eq.(3) in carminati_near-field_1999 ). For completeness, we present the derivation of (7) in the appendix, based on Planck’s law for the energy density of an e.m. field in thermal equilibrium.

Compared to a black body, the emissivity of a real body is modified by a mode-dependent emissivity factor $B_{i}(x,y;\omega,\hat{k})$, where $\hat{k}$ is the direction of emission (from the patch on ground to the satellite), and a factor $\cos\theta(x,y,h)$ of geometrical origin that takes into account the variation of the radiation with respect to the surface normal (i.e. the projection of the area of a patch of the surface onto the plane perpendicular to the propagation direction). The temperature $T(x,y)$ is then really an effective temperature, $T_{\text{eff}}(x,y)=T_{B}(x,y)\cos\theta(x,y,h)$, where the brightness temperature $T_{B}(x,y)$ is defined as the absolute temperature a black-body would need to have in order to produce the same intensity of radiation at the frequency and in the direction considered (see Appendix X.1.1). For simplifying notations, in the following we keep writing $T(x,y)$ for short instead of $T_{\text{eff}}(x,y)$, but keep in mind its physical meaning, which, after all, is crucial for data-analysis and fitting vegetation and surface models to observational data Kerr2012 . We thus arrive at the current correlator

[TABLE]

which can be considered the statistical model underlying the imaging concept, and $K_{3}=l_{c}^{3}K_{2}d/\tau_{c}$.

III Correlation of Fourier components

For each antenna, the electric field component $E_{i,{\bf r}_{1}}(t)$ is transduced into a voltage $U_{i,{\bf r}_{1}}(t)$. We denote the frequency response of the antennas and eventual subsequent filters by the complex function $A(\omega)$, the Fourier transform of the time-dependent response function of antenna and filter. In the frequency domain we have simply $\tilde{U}_{i,{\bf r}_{1}}(\omega_{1})=A(\omega_{1})\tilde{E}_{i,{\bf r}_{1}}(\omega_{1})$. With (8) we obtain the correlation function between the voltages at two different frequency components $\omega_{1}$, $\omega_{2}$ measured at the positions of the antennas with original positions ${\bf r}_{1}$ and ${\bf r}_{2}$,

[TABLE]

where now ${\bf r}^{\prime\prime}=(x^{\prime\prime},y^{\prime\prime},0)$, and $K_{5}=K_{3}K_{1}^{2}/(4\pi^{2})$. The correlation function $C_{ij}^{F}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ is the filtered version of the original unfiltered correlations $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$. We see from (9) that the latter can be obtained from the former simply by dividing through the product of the known filter functions, as long as the latter are non-zero. Of course, outside the frequency response of the antennas and filters, the measured correlations $C_{ij}^{F}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ vanish due to the vanishing of $A(\omega)$ and do not carry any information anymore. This will ultimately limit the frequency range over which information on the brightness temperature can be extracted, or, equivalently, leads to a finite geometrical resolution even if a $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ known for all frequencies would lead to perfect resolution. However, this appears only when inverting the measured signals and will be discussed in section V.1. For the moment we assume that we have access to the unfiltered $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ through (9) for all frequencies that we need, and base the general development of the theory on $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$.

We change integration variables from $t_{1},t_{2}$ to “center-of-mass” and relative times, $t=(t_{1}+t_{2})/2$ and $\tau=(t_{2}-t_{1})$, and introduce as well a new integration variable for the spatial integration, ${\bf r}^{\prime}\equiv{\bf r}^{\prime\prime}-{\bf v}_{s}t$. This implies ${\bf r}_{1}+{\bf v}_{s}t_{1}-{\bf r}^{\prime\prime}={\bf r}_{1}-{\bf v}_{s}\tau/2-{\bf r}^{\prime}$ and ${\bf r}_{2}+{\bf v}_{s}t_{2}-{\bf r}^{\prime\prime}={\bf r}_{2}+{\bf v}_{s}\tau/2-{\bf r}^{\prime}$. The Jacobian of both transformations is equal to

1. Furthermore, from now on we take the satellite to move in $x$ direction, ${\bf v}_{s}=v_{s}\hat{e}_{x}$, where $\hat{e}_{x}$ is the unit vector in $x$ direction. This leads to $T(x^{\prime\prime},y^{\prime\prime})=T(x^{\prime}+v_{s}t,y^{\prime})$.

The total phase $\Phi$ appearing as arguments of the exponential functions under the integrals in (11) is

[TABLE]

We see that $t$ only appears as prefactor of $(\omega_{2}-\omega_{1})$ in the phase (12), and as argument ${\bf v}_{s}t$ in $T(x^{\prime}+v_{x}t,y^{\prime})$. The integral over $t$ therefore boils down to a 1D Fourier transform of the intensity of the current fluctuations in the direction of the speed of the satellite, with conjugate variable proportional to the difference $\omega_{2}-\omega_{1}$ of the frequencies of the Fourier components that we correlate. This can be made more explicit by introducing a position variable $x=v_{s}t$ along the path of the satellite. For the conjugate variable we define $\kappa_{x}=(\omega_{2}-\omega_{1})/v_{s}$. We write $\kappa_{x}$ and not $k_{x}$ in order to distinguish this “wavevector” from the usual one obtained from a single frequency and dividing by $c$. We also introduce the “center of mass frequency” $\omega_{c}\equiv(\omega_{1}+\omega_{2})/2$. It will be called “center frequency” in the following for short, but should not be confused with the central frequency $\omega_{0}$ that is the fixed frequency in the middle of the band in which the satellite operates (e.g. $2\pi\times$1.4 GHz for SMOS). With all this, we see that

[TABLE]

We have defined $T(x^{\prime}+v_{s}t,y^{\prime})\equiv T_{{\bf r}^{\prime}}(x)$, where $v_{s}t=x$ is understood, and the spatial Fourier transform $\tilde{T}_{{\bf r}^{\prime}}(\kappa_{x})$ of the temperature field $T(x,y)$ in $x$-direction. This notation makes clear that in these coordinates the temperature depends both on ${\bf r}^{\prime}$ and $t$, even though the motion of the satellite combines the two arguments in a single one, ${\bf r}^{\prime}+{\bf v}_{s}t$. We can think of $\tilde{T}_{{\bf r}^{\prime}}(\kappa_{x})$ as the Fourier image of $T({\bf r}^{\prime}+x\hat{e}_{x})$ with respect to the $x$ coordinate, calculated with a starting point ${\bf r}^{\prime}$. I.e. for all ${\bf r}^{\prime}$, we have a 1D spatial Fourier transform of the intensity of the current fluctuations where the Fourier integral is defined with origin in ${\bf r}^{\prime}$. The Fourier images obtained by translation of ${\bf r}^{\prime}$ in $x$-direction are not independent. Rather we have

[TABLE]

We are thus led to

[TABLE]

We neglect the slow dependence of $\omega^{\prime 2}$ compared to the rapid oscillations of the phase factors in (III) and pull it out of the integral as a prefactor $\omega_{0}^{2}$. We can then perform the integral over $\omega^{\prime}$, and find

[TABLE]

We introduce center-of-mass and relative coordinates for ${\bf r}_{1}$ and ${\bf r}_{2}$, ${\bf R}=({\bf r}_{1}+{\bf r}_{2})/2$ and $\Delta{\bf r}={\bf r}_{2}-{\bf r}_{1}$. We further restrict $v_{s}\tau$ to values much smaller than $|{\bf R}-{\bf r}^{\prime}\pm\Delta r|$. This implies a limitation of the integration range for $\tau$ when calculating the Fourier components, but it is a mild one. Since $|{\bf R}-{\bf r}^{\prime}\pm\Delta r|\geq h$, it is enough to have $\tau\leq h/v_{s}$, which is typically of order $100$s, and therefore gives time to resolve Fourier components down to a hundredth of a Hertz. We can then approximate to first order in $v_{s}$,

[TABLE]

Neglecting terms of order $|\Delta{\bf r}+{\bf v}_{s}\tau|/|{\bf R}-{\bf r}^{\prime}|$ and of order $\beta=v_{s}/c$ in the prefactor of the exponential, as well as a second order term of order $\beta\omega_{c}\Delta r/c$ in the phase, the integral over the Dirac $\delta$-function gives

[TABLE]

and $K_{6}=(2\pi)^{3/2}\omega_{0}^{2}K_{5}/v_{s}$. The unit vector $\hat{e}_{{\bf R}-{\bf r}^{\prime}}$ is obtained by taking the original center of mass position of the antennas at ${\bf R}=(x_{0},0,h)$, and ${\bf r}^{\prime}=(x^{\prime},y^{\prime},0)$. Eq.(20) is one of the central results of this paper. It shows that the two-frequency correlation function of the fields at different antenna positions is related linearly via a 2D integral-transformation to the brightness temperature field in the source plane, or more precisely to the Fourier transform of that field in $x$-direction. With $T(x,y)$ defined on a 2D grid, the reconstruction of $T(x,y)$ from the measured correlation function thus becomes a matrix inversion problem that has to be performed numerically in general. A crucial question is the conditioning of the inversion problem. It will be studied in more detail in a subsequent paper dedicated to a numerical approach NumPaper .

Here we give a simplified analytical treatment that allows us to obtain estimates of the spatial and radiometric resolutions, and thus provide evidence that the inversion problem is sufficiently well conditioned for the reconstruction of $T(x,y)$ from the measured $C_{i,j}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$. For this, we study the situation where the vector $\Delta{\bf r}$ from antenna 1 to antenna 2 is orientated in $y^{\prime}$ direction, ${\bf r}_{2}={\bf r}_{1}+\Delta r\hat{e}_{y}$, in which case $\Delta{\bf r}\cdot\hat{e}_{{\bf R}-{\bf r}^{\prime}}=-\Delta ry^{\prime}/\sqrt{(x^{\prime}-x_{0})^{2}+y^{\prime 2}+h^{2}}$, and $\Delta r=|\Delta{\bf r}|$ denotes the spatial separation of the two antennas.

We switch to a dimensionless representation by taking as length scale the distance $\Delta r$ between the two antennas. We will express all other lengths in this unit, and introduce the dimensionless coordinates $\xi,\eta$ by $x^{\prime}=\xi\Delta r$, $y^{\prime}=\eta\Delta r$, and $\tilde{h}\equiv h/\Delta r$. The dimensionless height $\tilde{h}$ is for the standard parameters $\tilde{h}=7\cdot 10^{3}$. Eq.(20) then reads

[TABLE]

where $\tilde{T}(\kappa_{x},\eta)\equiv\tilde{T}_{0,\eta\Delta r}(\kappa_{x})$. The 1D integral kernel

[TABLE]

which is itself defined through an integral over $\xi$. For fixed $h$, $\Delta r$, $\omega_{c}$, and $\kappa_{x}$, the integral kernel $K(\alpha,\beta)$ is a function of $\eta$ that relates the 1D Fourier transform $\tilde{T}(\kappa_{x},\eta)$ to the observed correlation function by integration over $\eta$. Suppose that the integration over $\eta$ can be inverted by finding the inverse integral kernel. Integrating the inverse kernel over with the correlation function measured as function of the center frequency $\omega_{c}$, we then obtain $\tilde{T}_{0,\Delta r\eta}(\kappa_{x})$ for all $\eta$ and the chosen $\kappa_{x}$. If this can be done for all relevant $\kappa_{x}$, we obtain for each point on the $y$ axis the Fourier transform in $x$ direction of the intensity of the brightness temperature. Taking the inverse Fourier transform in $x$-direction, we obtain the full $x$- and $y$-dependent brightness temperatures. To proceed, we first study the integral kernel in the relevant parameter regimes.

IV Properties of the integral kernel

The arguments $\alpha,\beta$ of $K$ are given by eq.(21) as

[TABLE]

By their definition, we only need $\alpha,\beta\in\mathbb{R}$. For $\alpha$ we can consider that in the end the maximum $\kappa_{x}$ should be of the order of the inverse resolution $\Delta x_{\rm min}$ required in $x$ direction. Taking $\Delta x_{\rm min}$ of the order of one km, and using the standard parameters, we get $|\alpha|_{\rm max}\geq|\kappa_{x}h|\simeq 700$. With $\eta$ varying from $-\infty\ldots\infty$ (in reality, the extension of Earth limits the integration range to a maximum value of the order $10^{7}-10^{8}$), $\beta$ reaches its maximal value $\Delta r\omega_{c}/c$ for $\eta\to\infty$. For standard parameters, $\omega_{c}=2\pi\times$1.4 GHz, $|\beta|\lesssim 30$. Both $\alpha$ and $\beta$ can be positive or negative, such that there is also a regime, where $|\beta|\gg|\alpha|$, and we will find that this is the most important one. Note that from (22) we immediately obtain the relations

[TABLE]

We therefore restrict the following discussion to $\alpha,\beta\geq 0$.

Unfortunately, the integral over $\xi$ in (22) cannot be done analytically. However, we can find approximations for different cases. Consider first $\beta=0$. Using the methods of residues, one easily finds

[TABLE]

More generally, one can obtain a useful expansion for small $\beta$ by expanding $\exp(-i\beta/\sqrt{\xi^{2}+1})$ into a power series, and then integrating term by term. We find

[TABLE]

where $K_{n}(x)$ is the modified Bessel function of the second kind of order $n$. The zeroth order result (26) is recovered by observing that $K_{1/2}(x)=\sqrt{\frac{\pi}{2x}}e^{-x}$. For small $\beta$ the series converges rapidly, and one can even improve the agreement with the numerically calculated kernel by re-exponentiating the first few terms. For example, up to second order we have a polynomial $p_{0}+p_{1}\beta+p_{2}\beta^{2}$ which we wish to write as $p_{0}\exp(a_{1}\beta+a_{2}\beta^{2})$. Expanding the exponential in powers of $\beta$ and comparing powers up to order $\beta^{2}$, one finds $a_{1}=p_{1}/p_{0}$ and $a_{2}=p_{2}/p_{0}-(p_{1}/p_{0})^{2}/2$. When plotted together with the exact result, the thus obtained approximation agrees with $K(\alpha,\beta)$ for $\alpha=2$ visibly well up to $\beta\simeq 4$, i.e. well beyond the regime $\beta\ll 1$. For the fourth order re-exponentiated form the agreement extends up to about $\beta\simeq 5$. However, the exponential decay (26) already of the zeroth order term with $\alpha$ indicates that for the values of $\alpha\simeq 10^{2}$ to $10^{3}$ the contribution to the $\eta$ integral for values such that $\beta$ is of order of or smaller than one, can be entirely neglected.

In the opposite regime of large $\beta$, an approximation based on a stationary phase approximation can be found. More precisely, one needs $\beta\gg\alpha$. In that case one can treat $e^{-i\alpha\xi}/(\xi^{2}+1)$ as a slowly varying factor compared to the rapidly oscillating $e^{-i\beta/\sqrt{\xi^{2}+1}}$. The point of stationary phase of the latter term is found at $\xi=0$ (where the function has a maximum). The second derivative of the phase at $\xi=0$ equals 1. With this we get

[TABLE]

valid for $\beta/\alpha\gg 1$. Interestingly, the integral kernel becomes independent of $\alpha$ in this regime, which is of course a consequence of the fact that the stationary phase point is at $\xi=0$, thus eliminating the factor $\alpha$ in the phase of the prefactor. We furthermore see that in this regime there is no exponential suppression of the kernel.

For $\alpha=0$, the kernel can be evaluated exactly,

[TABLE]

where $J_{0}$ is the zeroth Bessel function, and $H_{0}$ the zeroth Struve function. Their asymptotic behavior gives back (29).

In Fig.1 we plot $\alpha,\beta$ as function of $\eta$. We see that a regime $\alpha<\beta$ exists for $\kappa<\kappa_{\rm Max}$, which defines $\kappa_{\rm Max}$ (see after eq.(32) for its precise value). For $\kappa\ll\kappa_{\rm Max}$, $\alpha\ll\beta$. The regime $\alpha\sim\beta\ll 1$ is also possible, but it is restricted to a tiny $\eta$ interval, such that its contribution to the integral over $\eta$ is negligible. For $\alpha\sim\beta\gg 1$, the stationary phase points of $\xi+(\beta/\alpha)/\sqrt{\xi^{2}+1}$ become relevant. Fig.1 shows the Im- and Re-parts of the six roots of the corresponding stationary phase equation. We see that only for $(\beta/\alpha)\gtrsim 2.5$ real stationary phase points exist. Since on the other hand $\alpha\sim\beta\gg 1$ occurs for sufficiently large $\kappa$ for almost all $\eta$ only for $\beta<\alpha$ (see Fig.1), the kernel is exponentially small in this regime $\alpha\sim\beta\gg 1$. Altogether, the only relevant regime is thus $\beta\gg\alpha\gg 1$.

While the asymptotic form of the integral kernel suggests the use of the orthogonality relations of Bessel functions, inverting (20) is nevertheless non-trivial due to the more complicated dependence of $\alpha$ and $\beta$ on $\eta$. However, the above asymptotic form allows one to obtain an approximate analytical inversion of the kernel that allows for an estimation of the resulting resolution, as we will show now.

V Estimation of geometrical resolution

V.1 Approximate analytical inversion of the integral kernel

At first sight, the requirement $\beta\gg\alpha$ appears unnatural given that $\alpha$ can already of the order $10^{2}$ to $10^{3}$. And indeed, this leads to a first rather stringent condition which must be met in order for the correlation function $C$ to be non-zero. In terms of the original parameters, $\beta/\alpha=\omega_{c}\eta/(c\kappa_{x}(\eta^{2}+\tilde{h}^{2}))$. For this to be much larger than one, one needs

[TABLE]

or $\Delta\omega/\omega_{c}\ll(v_{s}/(2c\tilde{h}))$, where we have used already the maximum value $1/(2\tilde{h})$ of the function of $\eta$ on the left hand side (lhs) in (31). For the standard parameters, we find $\Delta\omega/\omega_{c}\ll 1.66\cdot 10^{-9}$. When operating at $\omega_{c}$ in the GHz regime, this means that the correlation function essentially vanishes for $\Delta\omega$ larger than a few Hertz and thus bears no more information for the measurement of the position dependent brightness temperature.

Another way of seeing this is to observe that (31) limits the integration range for $\eta$: The lhs of (31) is a function that starts of at [math] for $\eta=0$, increases linearly, reaches a maximum of $1/2\tilde{h}$ at $\eta=\tilde{h}$, and decays as $1/\eta$ for large $\eta$. Condition (31) then limits the integration range of $\eta$ to an interval $\eta_{1}\leq\eta\leq\eta_{2}$ with

[TABLE]

where $\delta\equiv c\Delta\omega/(v_{s}\omega_{c})=c\kappa_{x}/\omega_{c}=c\kappa/(\omega_{c}\Delta r)=\chi\kappa$, $\kappa\equiv\kappa_{x}\Delta r$, and $\chi=c/(\omega_{c}\Delta r)\ll 1$ (see introduction and Fig.2). A finite real integration range exists only for $\delta<1/(2\tilde{h})$, equivalent to $\kappa<\frac{\omega_{c}\Delta r}{2c\tilde{h}}\equiv\kappa_{\rm Max}$. For given $\eta$, we have $\kappa\leq\kappa_{\rm max}(\eta)\equiv\frac{\omega_{c}\Delta r}{c}\frac{\eta}{\eta^{2}+\tilde{h}^{2}}\leq\kappa_{\rm Max}$. A finite minimal value of $\kappa$ can be deduced from a maximum desired snapshot size in $x-$direction. Also the requirement $\tau<h/v_{s}$ may bound the relevant values of $\kappa$ from below, as it leads to a smallest resolvable frequency, and thus also smallest resolvable $\Delta\omega$: $\Delta\omega>2\pi/\tau\Rightarrow\kappa=\Delta r\Delta\omega/v_{s}>\Delta r\frac{2\pi}{\tau v_{s}}\equiv\kappa_{\rm min}$.

As the contributions from areas outside the allowed range $\eta_{1}\leq\eta\leq\eta_{2}$ (or, correspondingly for negative $\eta$, $-\eta_{2}\leq\eta\leq-\eta_{1}$) are exponentially suppressed, we can limit the integration range of $\eta$ to that interval for a given $\kappa_{x}$, and replace the integral kernel by its approximate form, eq.(29), extended to $\beta<0$ by (25), yielding

[TABLE]

in the allowed range, and zero elsewhere. After the substitution $\zeta=\eta/\sqrt{\eta^{2}+\tilde{h}^{2}}$, the result for $C_{ii}$ can be written as

[TABLE]

where $\tilde{T}(\kappa,\eta)\equiv\tilde{T}_{0,\eta\Delta r}(\kappa_{x})$ (with $\kappa=\kappa_{x}\Delta r$, $\tilde{x}_{0}=x_{0}/\Delta r$), and $w(\zeta_{1},\zeta_{2},\zeta)$ is a window function equal to one for $\zeta_{1}\leq\zeta\leq\zeta_{2}$ and zero elsewhere. The window functions translate in a straight forward fashion the integration range for $\eta$ into an integration range for $\zeta$. By definition, $\zeta$ ranges from $-1,\ldots,1$. So $\zeta_{1},\zeta_{2}$, lie within this interval, $-1\leq\zeta_{1},\zeta_{2}\leq 1$, and the window functions take care of restricting the argument $\zeta$ of the integrand to the intervals $\pm[\zeta_{1},\zeta_{2}]$. We have replaced $\omega_{1},\omega_{2}$ by the equivalent information $\kappa\equiv\kappa_{x}\Delta r$ (related to $\Delta\omega$) and $\tilde{k}_{c}=\Delta r\omega_{c}/c$ (related to $\omega_{c}$), and consider $i=j$ only. Given eq.(34) it is tempting to try to recover $F(\kappa,\zeta,\tilde{k}_{c})$ by Fourier transform. However, the ${\rm sign}(\tilde{k}_{c})$ functions that appear in $F(\kappa,\zeta,\tilde{k}_{c})$ prevent (34) from being a simple Fourier integral. Moreover, from the measured data we only have $C_{ij}^{F}$, the filtered version of $C_{ij}$, that is restricted to a frequency range $\omega_{1},\omega_{2}\in\pm[\omega_{0}-\pi B,\omega_{0}+\pi B]$, where $B$ is the bandwidth (20MHz in SMOS for the L-band). We assume here for simplicity a Gaussian filter and the same for both antennas. For a real filter response function $A(t)$, its Fourier transform must satisfy $A(\omega)=A^{*}(-\omega)$. Taking also $A(\omega)$ as real, we can write it as

[TABLE]

where $G(\omega;\omega_{0},b)=\exp(-(\omega-\omega_{0})^{2}/(2b^{2}))/(\sqrt{2\pi}b)$ is a normalized Gaussian centered at $\omega_{0}$ with standard deviation $b\equiv 2\pi B$. The factor $\sqrt{b}\pi^{1/4}$ assures that for $\omega_{0}\gg b$, $A(\omega)$ is normalized according to $\int_{-\infty}^{\infty}|A(\omega)|^{2}\,d\omega=1$. Under the same condition we have

[TABLE]

where $\tilde{k}_{c0}=\Delta r\omega_{0}/c=1/\chi\simeq 2932.55$, and $\tilde{b}=\Delta rb/c\simeq 41.89$. This contains the approximation of using only the “diagonal” terms in the product of $A(\omega_{1})A^{*}(\omega_{2})$, i.e. the ones with ${\rm sign}(\omega_{1})={\rm sign}(\omega_{2})$, which is justified by the fact that $\Delta\omega\ll b\ll\omega_{0}$. Fourier transforming $C_{ii}^{F}$ with respect to $\tilde{k}_{c}$ (denoted by ${\cal F}_{\tilde{k}_{c}\to\zeta}$) gives a convolution product between the FT of the Gaussians (which is $(\sqrt{2}/\tilde{b})G(\zeta;0,\sqrt{2}/\tilde{b})e^{\pm i\tilde{k}_{c0}\zeta}$) and $F(\kappa,\zeta,\tilde{k}_{c})$, and leads to

[TABLE]

where we used that the sign of $\tilde{k}_{c0}$ in (37) determines the one of $\tilde{k}_{c}$ in (35). Thus, we get back the original function $\tilde{T}(\kappa,\frac{\zeta\tilde{h}}{\sqrt{1-\zeta^{2}}})=\tilde{T}(\kappa,\eta)$, cut by the two window functions and multiplied with $1/(\sqrt{\zeta}(1-\zeta^{2}))$, convoluted with the product of a Gaussian of width $\sqrt{2}/\tilde{b}$ and a rapidly oscillating cosine function. The factor $1/(\sqrt{\zeta}(1-\zeta^{2}))$ can be tracked back to the change of variables from $\eta$ to $\zeta$ and will distort the image at the nadir and at infinity. Sources at positive or negative $\eta$ contribute differently due to the different sign of the $\pi/4$ phase shift. This arises already in (33) due to the different phase shift in the asymptotics of the Struve functions for negative or positive arguments and leads to the sum over $\sigma=\pm$. In general, an exact inversion can not be simply done by Fourier transform but needs a numerical approach. Nevertheless, we can arrive at an estimation of the resolution by considering a single point source, as then only one of the two terms in the sum over $\sigma$ in (V.1) contributes, and the factor $1/(\sqrt{\zeta}(1-\zeta^{2}))$ becomes a simple numerical factor given by the position of the source.

V.2 Single point source and geometric resolution

V.2.1 Correlation function and reconstructed image

Let the point source be at position $x^{\prime\prime}=0,y^{\prime\prime}=\eta_{s}\Delta r$ and with polarization $i$, where $\eta_{s}$ is situated in the allowed range $0\leq\eta_{1}(\kappa)\leq\eta_{s}\leq\eta_{2}(\kappa)$ for some $\kappa$ in the desired range up to the largest considered $\kappa=2\pi/p_{x}$, where $p_{x}$ is the pixel size. We thus have

[TABLE]

which together with eq.(16) yields

[TABLE]

As $\eta$ is in the allowed range, we can use the approximate analytical form of the integral kernel, eq.(33), to get from eq.(34) the correlation function

[TABLE]

where $\theta(x)$ is the Heaviside theta-function. Considering (V.1), we may define an approximative reconstructed source function suitable for sources at $\eta_{s}>0$ through

[TABLE]

where ${\cal N}$ is a normalization constant. Due to the $\kappa-$dependence of the window functions, and the $\zeta$ dependence of the integral transform as compared to a simple Fourier-transform of $\tilde{T}$, one cannot get a normalization constant independent of the source field. In particular, for the single point source, ${\cal N}$ would depend on the position of the point source. However, we use $T_{\rm rec}$ only for estimating the geometric and radiometric resolution. For the former, all prefactors are irrelevant. For the latter, we avoid the problem by calculating relative uncertainties of $\sigma(T_{\rm rec})/T_{\rm rec}$ only, where any prefactor cancels. We hence set ${\cal N}=1$ in the following.

Inverting the Fourier transform in $\kappa$ leads to

[TABLE]

This equation is valid for all sources located in the positive $y$ plane, not necessarily point sources. When we re-express the correlation function through (34) and perform the Gaussian integral over $\tilde{k}_{c}$, we find a direct approximate formal relation between the FT of the original $T(x,y)$ in the upper half plane, and its reconstructed image $T_{\rm rec}(x,y)$,

[TABLE]

Using this expression, or by inserting (42) into (37), and the resulting filtered correlation function into (44), we find the reconstructed image of the single point source

[TABLE]

where $\zeta_{s}=\eta_{s}/\sqrt{\eta_{s}^{2}+\tilde{h}^{2}}$, $\kappa_{\rm max}(\zeta_{s})\equiv\zeta_{s}\sqrt{1-\zeta_{s}^{2}}/(\chi\tilde{h})$, and ${\rm sinc}(x)\equiv\sin(\pi x)/(\pi x)$.

V.2.2 Geometric resolution

We see that the reconstructed image of the point source is a series of narrow peaks spaced by the inverse of $\tilde{k}_{c_{0}}$ due to the rapidly oscillating $\cos$-function, under an approximate Gaussian in $y$-direction centered at the position of the source with a width in $\eta$ given by $\Delta\eta=\sqrt{2}\sqrt{\eta_{s}^{2}+\tilde{h}^{2}}/\tilde{b}\geq\sqrt{2}\tilde{h}/\tilde{b}=hc/((\Delta r)^{2}\sqrt{2}\pi B)$. It reminds one of a diffraction image from a double slit, even though there the envelope is a ${\rm sinc}$-function, not a Gaussian. Nevertheless, we adapt the definition of resolution from that example, namely that the best resolution is obtained from the smallest shift that makes a peak move into the next trough. This leads to

[TABLE]

hence $\Delta y=\Delta r\pi(\eta_{s}^{2}+\tilde{h}^{2})^{3/2}/(\tilde{k}_{c_{0}}\tilde{h}^{2})$. For $y\simeq h$, this is of the order $2\sqrt{2}\pi hc/(\Delta r\omega_{c})=\chi h$. The numerical value for the standard parameters gives $\Delta y\simeq$ 2.1 km, i.e. a resolution of the order of a kilo-meter. However, for actually achieving this resolution for an extended source, one has to face two issues: i.) The reconstructed point-source image should be brought as close as possible to a single narrow peak; and ii.) one has to deal with the different phases from sources at positive or negative $\eta$. The first issue can be addressed by superposing correlation functions from pairs of antennas at different separation, and/or changing the considered central frequency. This shifts the pattern of peaks due to the cos-function, and one can engineer a rather narrow central peak (see NumPaper for details). The second issue should be absent in a numerically exact inversion of the integral kernel. The Gaussian envelope has a width $hc/(\sqrt{2}\pi B\Delta r)$ given by the inverse bandwidth, which is much larger than the width of a single peak, namely by a factor $\omega_{c}/(4\pi B)\simeq 35$ for the standard parameters.

The resolution in $x$-direction follows from the effective wave vector $\kappa_{\rm max}$ in the ${\rm sinc}$ function. It depends on the position of the source and reaches its maximum possible value $\kappa_{\rm Max}$ for $\eta_{s}=\tilde{h}$ (i.e. $y_{s}=h$). The inverse of $\kappa_{\rm Max}$ thus gives the best possible resolution in $x$-direction:

[TABLE]

We conclude that both in $x$- and $y$-direction one can expect a geometric resolution of the order $h\chi=c/(\Delta r\omega_{c})$ for sources close to $y=h$. For sources close to $y=0$, $\kappa(\eta_{s})$ goes to zero $\propto\eta_{s}$, whereas for larger $y_{s}$ the decay of $\kappa(\eta_{s})$ is $\propto 1/\eta_{s}$. The geometric resolution in $x$-direction deteriorates correspondingly. The resolution in $y$-direction, on the other hand, depends only weakly on the source-position, as $(\eta_{s}^{2}+\tilde{h}^{2})^{3/2}/(\tilde{k}_{c}\tilde{h}^{2})$ increases monotonically from $\tilde{h}$ at $y_{s}=0$ to $2\sqrt{2}$ at $y_{s}=h$, and keeps growing slowly beyond $y_{s}=h$. It is remarkable that correlating electric fields at two different frequencies can lead to a resolution that is given by the central frequency.

The definition of $\kappa_{\rm max}$ is based on the request that the stationary phase approximation (SPA) holds in the regime $\beta\gg\alpha$. In practice, the SPA is almost always better than expected, such that in the end the result $h\chi$ might be a conservative estimate of the geometric resolution.

VI Radiometric resolution

Besides the geometric resolution, the radiometric resolution (RR), i.e. the smallest difference in temperature that the system can measure for a given pixel, is the most important characteristics of the satellite imaging system. Here we calculate the RR for the idealized situations of a single point source considered above and for a uniform temperature field in the positive half-plane $y>0$.

VI.1 Fluctuations of the reconstructed temperature profiles

The idea behind the calculation of RR is that the electric field measurements yield random values, whose fluctuations and correlations reflect the thermal nature of the radiation field. Thus, if with the same field $T(x,y)$ one repeated the measurements many times, one would obtain different correlation functions in each run, and thus, after inverting the linear relationship between $C_{ij}$ and $T(x,y)$, also different reconstructed $T(x,y)$ (called $T_{\text{rec}}$ in the following) in each run. The (relative) RR is then defined as the standard deviation $\sigma(T_{\text{rec}}(x,y))$ divided by the average $T_{\text{rec}}(x,y)$ for a given position $x,y$. In general it will vary as function of $x,y$, and also depend on the temperatures at all positions, a behavior well known from standard spatial aperture synthesis. In reality, things become still a bit more complicated, as the measured signal is a superposition of the e.m. field emitted by the antenna itself (at temperature $T_{a}$), and the radiated field from the surface of Earth. However, these fields are uncorrelated, and their averaged squares just add up. For simplicity, we will neglect the noise contribution of the antennas in this first analysis, which amounts to calculating lower bounds of $\sigma(T_{\text{rec}})$.

Starting point of the calculation is the assumption that the current fluctuations ${\bf j}({\bf r}^{\prime\prime},t)$ which are at the origin of the radiated thermal field are described by a random Gaussian process, both in time and space (see sec. II). This implies immediately, that also the temporal FT $\tilde{{\bf j}}({\bf r}^{\prime\prime},\omega^{\prime})$ of the current fluctuations is a Gaussian process, now over space and frequency. Finally, the connection between $\tilde{{\bf j}}({\bf r}^{\prime\prime},\omega^{\prime})$ and $\tilde{E}_{{\bf r}_{1}}(\omega_{1})$ is linear, which implies that $\tilde{E}_{{\bf r}_{1}}(\omega_{1})$ is a Gaussian process over ${\bf r}_{1}$ and $\omega_{1}$. By the nature of this variable, it is a complex Gaussian process. One easily shows that the average of $\tilde{E}_{{\bf r}_{1}}(\omega_{1})$ equals zero (if the average of all current components is, which must be true at thermal equilibrium). The correlation function $C_{ij}$ is the (complex) covariance matrix of this Gaussian process, and all higher correlations can be expressed in terms of it.

In order to assess the fluctuations of $T_{\rm rec}$ we first define a product of Fourier coefficients of ${\bf E}$ from a single run (denoted by a $\hat{}$ ),

[TABLE]

and its corresponding filtered version $\hat{C}^{F}({\bf r}_{1},{\bf r}_{2},\kappa,\tilde{k}_{c})=\hat{C}({\bf r}_{1},{\bf r}_{2},\kappa,\tilde{k}_{c})A(\omega_{1})A^{*}(\omega_{2})$ (with $\omega_{1},\omega_{2}$ expressed in terms of $\kappa,\tilde{k}_{c}$).

The fluctuations of $T_{\rm rec}(x,y)$ are defined as $\Delta T_{\rm rec}(x,y)\equiv\langle T_{\rm rec}(x,y)^{2}\rangle-\langle T_{\rm rec}(x,y)\rangle^{2}$, where the average is over the thermal ensemble. With $T_{\rm rec}(x,y)$ from (44), one finds

[TABLE]

It can be shown (see Sec.X.1.2) that in the narrow frequency intervals considered here the Gaussian random processes given by the $\tilde{E}_{z,{\bf r}_{i}}(\omega)$ is circularly symmetric. It hence enjoys the property (see eq. 8.250 in Barrett03 ),

[TABLE]

where we have abbreviated $E_{i}\equiv\tilde{E}_{z,{\bf r}_{(i\,{\rm mod}\,2)}}(\omega_{i})$. The correlation function contained in the large parentheses of the second line of (51) becomes

[TABLE]

with $\kappa_{ij}=\Delta r(\omega_{j}-\omega_{i})/v_{s}$, $\tilde{k}_{cij}=\Delta r(\omega_{i}+\omega_{j})/(2c)$ $\forall i,j$, and where we have used $\langle\hat{C}^{F}({\bf r}_{1},{\bf r}_{2},\kappa,\tilde{k}_{c})\rangle=C^{F}({\bf r}_{1},{\bf r}_{2},\kappa,\tilde{k}_{c})$.

The fact that $C_{zz}^{F}$ and $C_{zz}^{F*}$ in (53) contain the same position arguments twice makes that we cannot evaluate it directly through eq.(21), as the coordinate transformation to dimensionless variables based on the rescaling with $\Delta r$ becomes singular. We therefore have to go back a step to eq.(20) which yields

[TABLE]

Comparing with (21) we see that this result corresponds formally to $\omega_{c}=0$ in that equation, rather than $\Delta r=0$, and (55) is recovered by using the exact result (26). We can now re-introduce dimensionless variables via the same rescaling with $\Delta r$, where, however, $\Delta r$ is still given by $\Delta r=|{\bf r}_{2}-{\bf r}_{1}|$, and ${\bf r}_{i}$ denote as before the positions of the two antennae at $t=0$, only one of which still enters as argument in $C_{zz}({\bf r}_{1},{\bf r}_{1},\kappa_{13},\tilde{k}_{13})$, respectively $C_{zz}^{F*}({\bf r}_{2},{\bf r}_{2},\kappa_{24},\tilde{k}_{c24})$. This gives

[TABLE]

We calculate $\Delta T_{\rm rec}$ for the single point source considered in Sec. V.2 at the position of the source, i.e. $\Delta T_{\rm rec}(x_{s},y_{s})=\langle T_{\rm rec}(x_{s},y_{s})^{2}\rangle-\langle T_{\rm rec}(x_{s},y_{s})\rangle^{2}$, and for a uniform temperature field.

VI.2 Single point source

For the single point source at position $(0,y_{s})$, the correlation function $C_{zz}({\bf r}_{1},{\bf r}_{1},\kappa,\tilde{k})$ becomes (see (41))

[TABLE]

Insert this into eq.(51) to find

[TABLE]

We change integration variables to $\kappa_{13},\tilde{k}_{c13},\kappa_{24},\tilde{k}_{c24}$. The Jacobian is 1. In the product of the Gaussians only the diagonal terms (i.e. with the same signs in front of $\tilde{k}_{c0}$) contribute in the relevant regime $\tilde{k}_{c0}\gg\tilde{b}$, as for opposite signs $\tilde{k}_{c12}\simeq\tilde{k}_{c34}\simeq 0$. Finally, we approximate

[TABLE]

and pull that factor out from the integral which is permissible for all ranges of variables for which the product of Gaussians is non-negligible. The integrals can then be performed, and we find for the standard deviation $\sigma(T_{\rm rec}(0,y_{s}))\equiv\sqrt{\Delta(T_{\rm rec}(0,y_{s}))}$

[TABLE]

where in the last step we have used that $\tilde{h}\gg 1$ and $\beta\zeta_{s}\ll 1$ (where once more $\beta=v_{s}/c$).

From (46) we find

[TABLE]

Combined with (60) we obtain the relative RR

[TABLE]

For $\eta_{s}=\tilde{h}$, the relative RR is of order 0.055, corresponding at $T=300$ K to $\sigma(T_{\rm rec}(0,h))\simeq 16.5$ K.

VI.3 Uniform temperature field

We now look at the second standard situation considered commonly for the determination of the radiometric resolution, namely a field of constant temperature. More precisely, we consider

[TABLE]

The restriction to sources in the upper plane is due to the fact that we still want to use eq. (44) for calculating the reconstructed temperature profile. The cut-off $\hat{y}$ arises physically from the size of the Earth and prevents a divergence of the correlation function.

From (13) we obtain

[TABLE]

The correlation function (56) becomes

[TABLE]

with $\hat{\eta}\equiv\hat{y}/\Delta r$. For $\hat{y}=R_{E}\simeq 6370$ km the radius of Earth, and $\Delta r=100$ m, $\hat{\eta}=63700$. We see that here the correlation function is perfectly diagonal in frequency, which of course reflects the lack of structure of the temperature field in $x$-direction. Hence, we can set everywhere $\kappa=0$, which greatly simplifies the analysis. The cutoff $\hat{\eta}$ in eq.(65) prevents a logarithmic divergence that arises from $1/\sqrt{\eta^{2}+\tilde{h}^{2}}\sim 1/\eta$ for $\eta\to\infty$. Eq.(65) can be extended to a temperature field that is uniform everywhere from $-\hat{y}$ to $\hat{y}$. In this case, the $\zeta$-integral starts at $-1+\epsilon$ rather than at 0. However, in this situation we cannot use (44) anymore as it is valid only for sources at positive $y$ (see the discussion after (V.1)).

Eq.(65), when inserted into (51), and with the same change of integration variables and approximation (59), leads to

[TABLE]

The reconstructed temperature field (45) is given by

[TABLE]

Unfortunately, no closed analytical form could be found for the remaining integral, and even a numerical evaluation is not straight forward, as the Gaussian yields a very narrow peak, broader, however, than the period of the $\cos$-function. But we can get an estimate of $T_{\rm rec}$ by replacing the Gaussian (normalized to an integral equal 1) with a rectangular peak of width $a\sigma$ and height $1/(a\sigma)$ centered, as the Gaussian, at $\zeta$. Here, $\sigma=\sqrt{2}/\tilde{b}$, and $a$ is a parameter of order 1. This gives

[TABLE]

A numerical evaluation of the integral is now relatively straight forward and shows a slowly varying $T_{\rm rec}(0,y)$ as function of $\zeta$ in the interval $\zeta\in[a\sigma/2,\hat{\zeta}-a\sigma/2]$, whereas outside this interval it oscillates rapidly. The slow variation arises from the factor $\sqrt{|\zeta^{\prime}|}(1-\zeta^{\prime 2})$ that distorts this approximately reconstructed image. Pulling out this slowly varying factor in order to get an analytical estimate of the order of magnitude of $T_{\rm rec}(x,y)$, we are led to

[TABLE]

for $\zeta\in[a\sigma/2,\hat{\zeta}-a\sigma/2]$. Hence, in this interval and apart from the distorting factor $1/\sqrt{|\zeta|}(1-\zeta^{2})$ identified previously, we recover a constant temperature field. The value of the reconstructed temperature depends on the precise value of $a$ as well as the ratio $\tilde{k}_{c0}/\tilde{b}$. Outside the mentioned interval, $T_{\rm rec}(x,y)$ oscillates again as function of $\zeta$, which can be understood from the fact that the box is cut-off when $\zeta$ gets within a distance $a\sigma/2$ of [math] or $\hat{\zeta}$. The sought-for order of magnitude can be estimated from the maximum value of (70) as function of $a$. As for standard parameters $\tilde{k}_{c0}/\tilde{b}\simeq 70$, we can bound the $\sin$-function by one (while still having $a\sim 1$), in which case we obtain $T_{\rm rec}(x,y)\simeq T_{0}/\tilde{k}_{c0}=T_{0}\chi$ in the mentioned $\zeta$-interval. With all this, and approximating $\sqrt{\tilde{h}^{2}+\hat{\eta}^{2}}\simeq\hat{\eta}$ in the logarithm in (67), we find the order of magnitude $\sigma(T_{\rm rec}(x,y))/T_{\rm rec}(x,y)\sim\ln(2\hat{\eta}/\tilde{h})/\chi^{3/2}$ For $\zeta\sim 1$, this is of order $\sim 10^{5}$ for standard parameters, i.e. a catastrophically large uncertainty. A small value of $\sigma(T_{\rm rec}(x,y))/T_{\rm rec}(x,y)$ is possible only if $\sqrt{|\zeta|}(1-\zeta^{2})$ is very small, but apart from the fact that one should not rely on this image-distorting factor, it could only be sufficiently small for $y$ unrealistically close to the nadir.

If one traces back the difference to the single-point source, one realizes that while $\sigma(T_{\rm rec})$ scales in both cases as $1/\sqrt{\chi}$, the difference comes from $T_{\rm rec}$ itself: for the single point source it is of order $1/\chi$, but for uniform $T$ in the upper half plane of order $\chi$, which explains a factor $1/\chi^{2}$ worse relative RR for the latter compared to the former. The factor $1/\chi$ in the single-point $T_{\rm rec}$ arises from the cut-off of the $\kappa$ integral: $\kappa_{\rm Max}$ scales as $1/\chi$, and for $x=0$ the $\kappa$-integral in (44) just gives a factor $2\kappa_{\rm Max}\sim 1/\chi$, as the correlation-function is independent of $\kappa$ in this case. On the other hand, for constant temperature in the upper half-plane, the cut-offs $\kappa_{\rm Max}$ do not play a role, as the $\delta(\kappa)$-function picks up only $\kappa=0$. This leads to the loss of one factor $1/\chi$ in $T_{\rm rec}$. The second one comes from the integration over $\zeta^{\prime}$ in (68): the rapidly oscillating $\cos$-function leads to a factor $1/\tilde{k}_{c0}=\chi$, whereas for the point-source only a single point $\zeta=\zeta_{s}$ contributes, such that the cosine is of order one.

In the light of the RR of standard radiometers that typically scales as $\sigma(T)\propto 1/\sqrt{Bt_{\rm int}}$, where $t_{\rm int}$ is the integration time, the fact that $\sigma(T_{\rm rec})$ in eqs.(60,67) is independent of the bandwidth is rather surprising. Formally, the disappearing of $\tilde{b}$ can be traced back to using the lowest order in the Laplace approximation of (58). The next order corrections are of order $\tilde{b}$, such that $\sigma(T_{\rm rec})\to\sigma(T_{\rm rec})(1+{\cal O}(\tilde{b}))$. One expects the sign of the correction to be positive, as the integrand is positive everywhere, and the lowest order approximations amounts to replacing the Gaussians by normalized Dirac-delta functions. Hence, for small but finite $\tilde{b}$, $\sigma(T_{\rm rec})$ is expected to increase with $\tilde{b}$, which is contrary to the behavior of standard radiometers. Standard radiometers are based on the van Cittert-Zernike theorem, which gives the reconstructed temperature field as Fourier transform of the observed visibilities at a fixed frequency. Different frequencies at the source are uncorrelated, and the scaling of $\sigma(T_{\rm rec})/T_{0}\propto 1/\sqrt{Bt_{\rm int}}$ just reflects averaging over a number of independent measurements that scales $\propto B\,t_{\rm int}$. In FouCoIm the information is in the correlation between different, very narrowly spaced Fourier components, and averaging over the central frequency does not lead to additional information (see also Sec.VII.2). Therefore a larger bandwidth does not improve RR.

VII Noise reduction

The bad signal-to-noise ratio for the radiometric resolution in the case of a uniform temperature field in the upper half plane, makes it essential to consider measures that lead to a noise reduction and in particular averaging schemes.

VII.1 Averaging over time

Instead of examining $T_{\rm rec}(x,y)$, we consider here for simplicity directly the fluctuations of the measured ”single shot” correlation function $\hat{C}_{ij}$. The first idea that comes to mind for reducing the fluctuations of $\hat{C}_{ij}$ is to average over the origin of the time-interval from which we construct the Fourier transform. Note that this is very different from an ensemble that one would obtain by displacing the initial position ${\bf r}_{i}$. But averaging over the origin of time only leads to an overall factor:

[TABLE]

So it is obvious that this kind of averaging is useless. This is indeed to be expected, as all available data were already used. The situation improves only slightly if the FTs are calculated from a finite stretch of data (say over a duration $\tau_{F}$). Then shifting the origin in time will include some new random data, but since we must have $\tau_{F}\gg\tau_{a}$, it is clear that we still use essentially the same data with the exception of some new data points at the edge of the interval of length $\tau_{F}$.

VII.2 Additional frequency pairs

Using only a small frequency separation of width $\Delta\omega$ about the central frequency $\omega_{c}$ which itself is allowed to vary over a large bandwidth $B$ appears to be a very wasteful use of all the pairs of frequency components $(\tilde{E}_{z,{\bf r}_{1}}(\omega_{1}),\tilde{E}_{z,{\bf r}_{2}}(\omega_{2}))$. Can we use different measured correlations $\hat{C}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ with sufficiently different $\omega_{c}=(\omega_{1}+\omega_{2})/2$ as independent data for improving the radiometric sensitivity? In order to answer this question, we need to calculate the covariance matrix $V$ between two different correlators,

[TABLE]

as well as the pseudo-covariance matrix $M$,

[TABLE]

Both matrices together determine the statistical properties of the random process $\hat{C}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$. Note that despite the fact that $E_{z,{\bf r}}(\omega)$ can be considered a circularly symmetric Gaussian process (see Appendix) over ${\bf r}$ and $\omega$ in the narrow frequency band we are interested in, the same is not true for $\hat{C}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})-\langle\hat{C}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})\rangle$ (which is not even Gaussian). We need to know whether both $V$ and $M$ essentially vanish for almost all pairs of pairs of frequencies, with the first pair $(\omega_{1},\omega_{2})$ in a first region (notably in the central narrow strip $S\equiv\omega_{2}\in[\omega_{1}-\Delta\omega,\omega_{1}+\Delta\omega]$), and the second pair $(\omega_{1}^{\prime},\omega_{2}^{\prime})$ in another region in the $(\omega_{1},\omega_{2})$ plane that we may want to consider, whereas the correlation functions $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$ and $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1}^{\prime},\omega_{2}^{\prime})$ themselves should still be non-zero. Such a situation would signal statistically independent non-vanishing correlation functions. However, we saw that only within a central narrow strip $S$ (whose width is given by $\kappa_{\rm max}(\eta)$) in the $(\omega_{1},\omega_{2})-$plane $C_{ij}$ is non-zero, and within that strip all pairs of frequencies are used for obtaining a single profile $T(x,y)$. Here we show the same thing once more by proving that for $M$ and $V$ to vanish the second pair of frequencies $\omega_{1}^{\prime},\omega_{2}^{\prime}$ must be not in $S$ — where, however, $C_{ij}({\bf r}_{1},{\bf r}_{2},\omega_{1}^{\prime},\omega_{2}^{\prime})$ vanishes.

To see this, one first shows with the help of (52) and in a few lines of calculation that

[TABLE]

We have $C_{zz}({\bf r}_{1},{\bf r}_{1},\omega_{1},\omega_{1}^{\prime})$ from (56), where now $\kappa=(\omega_{1}^{\prime}-\omega_{1})\Delta r/v_{s}$, and correspondingly for $C_{zz}({\bf r}_{2},{\bf r}_{2},\omega_{2},\omega_{2}^{\prime})$. Whether $V,M$ are large or small can be judged by comparing it to the product of the standard deviations of each factor. This corresponds to calculating Pearson’s product-moment coefficients pearson_note_1895 $V^{\rm res}\equiv\frac{V}{\sigma(\hat{C})\sigma(\hat{C^{\prime}})}$ and $M^{\rm res}\equiv\frac{M}{\sigma(\hat{C})\sigma(\hat{C^{\prime}})}$ where we define, for complex $\hat{C}$, $\sigma(\hat{C})\equiv\sqrt{\sigma^{2}(\Re\hat{C})+\sigma^{2}(\Im\hat{C})}$, and $\hat{C}\equiv\hat{C}_{zz}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$, $\hat{C^{\prime}}\equiv\hat{C}_{zz}({\bf r}_{1},{\bf r}_{2},\omega_{1}^{\prime},\omega_{2}^{\prime})$ for short. Going through the same calculation as for $V$, we find after some algebra

[TABLE]

where

[TABLE]

and hence $I(0)=\int_{-\infty}^{\infty}\frac{\tilde{T}_{0,\eta}(0)}{\sqrt{\eta^{2}+\tilde{h}^{2}}}d\eta$. This implies

[TABLE]

For $M$ we have

[TABLE]

where

[TABLE]

with $\kappa_{12^{\prime}}\equiv(\omega_{2}^{\prime}-\omega_{1})/v_{s}$, $\kappa_{1^{\prime}2}\equiv(\omega_{2}-\omega_{1}^{\prime})/v_{2}$, $\omega_{c12^{\prime}}\equiv(\omega_{1}+\omega_{2}^{\prime})/2$, and $\omega_{c1^{\prime}2}\equiv(\omega_{1}^{\prime}+\omega_{2})/2$. From the properties of the integration kernel $K$ we know that $C_{zz}({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2}^{\prime})$ vanishes iff $|\omega_{1}-\omega_{2}^{\prime}|\gtrsim(\omega_{1}+\omega_{2}^{\prime})v_{s}/(2c\tilde{h})$, and correspondingly for $C_{zz}({\bf r}_{1},{\bf r}_{2},\omega_{1}^{\prime},\omega_{2})$. Hence, for $|M|\ll 1$ and $\omega_{1,2}$ and $\omega_{1,2}^{\prime}$ all of order $\omega_{0}$, we need $|\omega_{1}-\omega_{2}^{\prime}|\gtrsim\omega_{0}v_{s}/(c\tilde{h})$ or $|\omega_{1}^{\prime}-\omega_{2}|\gtrsim\omega_{0}v_{s}/(c\tilde{h})$. Note that $\omega_{0}v_{s}/(c\tilde{h})=(1/\chi)v_{s}/h\gg v_{s}/h$.

For determining the properties of $V$, we consider our two previous cases of sources.

Case 1: Single point source. Here we have $\tilde{T}$ independent of $\kappa$, see eq.(41), which inserted into (79) yields

[TABLE]

and hence

[TABLE]

For sources at $\eta_{s}\sim\tilde{h}$, we have therefore $|V^{\text{res}}|\ll 1$ iff $|\omega_{1}^{\prime}-\omega_{1}|>\delta\omega$ or $|\omega_{2}^{\prime}-\omega_{2}|>\delta\omega$ where $\delta\omega\equiv v_{s}/(\Delta r\sqrt{\eta_{s}^{2}+\tilde{h}^{2}})\sim v_{s}/(\Delta r\tilde{h})=v_{s}/h\sim 10^{-2}$ Hz.

Case 2: Constant temperature field in the positive upper half plane. Here, $\tilde{T}(\kappa_{x},\eta)$ is given by eq.(64). Hence, $V^{\rm res}=0$ as soon as $\omega_{1}^{\prime}\neq\omega_{1}$ or $\omega_{2}^{\prime}\neq\omega_{2}$. Of course, the $\delta(\kappa_{x})$ function in (64) arises from the complete lack of structure of the temperature profile in $x$-direction. More realistic is at least a cut-off at the size of Earth, which we take as the same as in $y$ direction. In that case one finds $\tilde{T}_{{\bf r}^{\prime}}(\kappa_{x})=(\hat{y}T_{0}/\pi){\rm sinc}\left(\frac{\kappa_{x}\hat{y}}{\pi}\right)$ and hence

[TABLE]

The exponential factor in the integral makes again that $V^{\rm res}$ vanishes essentially if $|\omega_{i}^{\prime}-\omega_{i}|\gtrsim v_{s}/h$ for $i=1$ or $i=2$.

Comparing with the situation for $M$, we find that for both types of sources considered, $V$ vanishes much more rapidly as function of the separation of two frequencies, as there is no factor $1/\chi$ multiplying $v_{s}/h$. Hence, the request for vanishing $M$ is more restrictive.

The question of the usefulness of considering other frequency pairs can now be phrased as: Can one find pairs of frequencies $(\omega_{1}^{\prime},\omega_{2}^{\prime})$ such that $|\omega_{2}^{\prime}-\omega_{1}|\gg\Delta\omega\equiv(1/\chi)v_{s}/h$ or $|\omega_{1}^{\prime}-\omega_{2}|\gg\Delta\omega$ while still $|\omega_{2}^{\prime}-\omega_{1}^{\prime}|\lesssim\Delta\omega$, for all frequencies $\omega_{1},\omega_{2}$ with $|\omega_{2}-\omega_{1}|\lesssim\Delta\omega$ used in the reconstruction of a temperature profile from $C({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})$? For a single frequency pair $(\omega_{1}^{\prime},\omega_{2}^{\prime})$ all conditions can be easily satisfied. It is enough that both pairs $(\omega_{1},\omega_{2})$ and $(\omega_{1}^{\prime},\omega_{2}^{\prime})$ be inside the strip $S$, and at the same time far away from each other, i.e. $|\omega_{1}-\omega_{1}^{\prime}|\gg\Delta\omega$, which implies $|\omega_{2}^{\prime}-\omega_{1}|\gg\Delta\omega$ and $|\omega_{1}^{\prime}-\omega_{2}|\gg\Delta\omega$ at the same time. However, the difficulty arises from the fact that we use already all pairs $(\omega_{1},\omega_{2})$ in the full available band-width for the reconstruction of a single temperature profile. This can be seen e.g. from eq.(V.1), where we integrate over all $\tilde{k}_{c}=\Delta r(\omega_{1}+\omega_{2})/(2c)$ for recovering $\tilde{T}$. Hence, there are really no new frequency pairs that can be used for improving the signal/noise ratio of the reconstructed temperature profile.

The same conclusion can be arrived at more formally by calculating the correlations between temperature profiles obtained from different center frequencies. Let $T_{\text{rec}}(x,y;\omega_{0})$ be the reconstructed temperature profile given by eq.(44), where we now keep explicit the dependence on the center frequency $\omega_{0}$, hidden in that equation in the filter functions $A(\omega_{1},\omega_{0})$, see eqs.(37,38), and $C_{ii}^{F}\to\hat{C}_{ii}^{F}$ is understood, so as to get the temperature profile from a single realization of the noise process. We define the correlation function

[TABLE]

and its renormalized dimensionless version

[TABLE]

that obviously satisfies $K_{\text{rel}}(T_{\text{rec}}(\omega_{01}),T_{\text{rec}}(\omega_{01}))=1$. We have

[TABLE]

where $\kappa_{1}=(\omega_{2}-\omega_{1})\Delta r/v_{s}$, $\kappa_{2}=(\omega_{2}^{\prime}-\omega_{1}^{\prime})\Delta r/v_{s}$, $\tilde{k}_{c1}=(\omega_{1}+\omega_{2})\Delta r/(2c)$, $\tilde{k}_{c2}=(\omega_{1}^{\prime}+\omega_{2}^{\prime})\Delta r/(2c)$, $\tilde{k}_{c0}^{(i)}=\omega_{0}^{(i)}$ ($i=1,2$), $F(\tilde{k}_{c},\tilde{k}_{c0},\frac{\tilde{b}}{\sqrt{2}})=A(\omega_{1},\omega_{0})A^{*}(\omega_{2},\omega_{0})$ with $\tilde{k}_{c0}=\Delta r\omega_{0}/c$ (see eq.(36)), and we have used that $T_{\text{rec}}\in\mathbb{R}$. We evaluate $K(T_{\text{rec}}(\omega_{01}),T_{\text{rec}}(\omega_{02}))$ for the case of constant temperature in the upper half plane. Using (73), (65), and switching momentarily to integration variables $\omega_{1},\omega_{2},\omega_{1}^{\prime},\omega_{2}^{\prime}$, and then back to $\kappa_{1}\,\tilde{k}_{c1}$, we are led to

[TABLE]

The integral is clearly real, as it should. The integral over $\kappa_{1}$ leads, when integrated from $-\infty$ to $\infty$ to a divergent factor, but that factor cancels (together with the remaining prefactor $T_{0}^{2}v_{s}/(2\pi c)$) when we consider the re-scaled version of the correlation function $K_{\text{rel}}(T_{\text{rec}}(\omega_{01}),T_{\text{rec}}(\omega_{02}))$. If we set $\tilde{k}_{c0}^{(2)}=\tilde{k}_{c0}^{(1)}+\delta\tilde{k}_{c0}$, it is clear that the only remaining scale for $\delta\tilde{k}_{c0}$ is $\tilde{b}/\sqrt{2}$. The remaining integral over $\tilde{k}_{c1}$ in (90) can in fact be evaluated analytically. The result is too cumbersome to be reported here, but plotting it as function of $\delta\tilde{k}_{c0}$ shows that indeed the correlations decay only on a scale of order $\tilde{b}$. This proves that by shifting the center frequency within the available bandwidth one cannot gain independent estimates of $T_{\text{rec}}$ that would allow one to improve substantially the signal-to-noise ratio.

VII.3 Additional antennas

So far we considered only two antennas. As mentioned before, in order to obtain a reconstructed single source image with a single peak, one may sum the correlated signals from several antenna pairs. It is to be expected that this will reduce $\sigma(T_{\rm rec})/T_{\rm rec}$, but we have to figure out how far two pairs of antenna have to be separated in order to produce essentially uncorrelated correlation functions. To answer this question, we have to generalize eq.(73) to pairs of correlators at different points ${\bf r}_{1}^{\prime},{\bf r}_{2}^{\prime}$. We define

[TABLE]

where we take ${\bf r}_{i+2}={\bf r}_{i}+\rho_{i}\hat{e}_{y}$ for $i=1,2$, i.e. the antennas in the second pair are shifted by distance $\rho_{i}$ in $y$-direction compared to the corresponding ones in the first pair. From (52) we have

[TABLE]

where from (21) and (30)

[TABLE]

The corresponding result for $C({\bf r}_{4},{\bf r}_{2},\omega_{2},\omega_{2})$ is obtained from the last line in (93) by replacing $\rho_{1}\to\rho_{2}$. For the uniform temperature field in the upper half plane up to a cut-off $\hat{y}$ and also a cut-off of the same value in $x$-direction, we have $T(0,\eta)=\hat{y}T_{0}/\pi$ for $0\leq y\leq\hat{y}$. No closed form was found for the remaining integral over $\zeta$, but a closed form is easily obtained if we neglect the slowly varying envelope $1/(1-\zeta^{2})$, which is legitimate for cut-offs $\hat{\eta}$ not too close to 1 and gives an idea on which length-scale $V_{r}$ will vanish. Plotting the results of the integration one finds that both real and imaginary part decay on a scale of $\rho_{i}\omega_{0}/c\sim 1$, where we have used again $\omega_{1}\simeq\omega_{2}\simeq\omega_{0}$. Hence, for the correlation functions of two pairs of antenna to decorrelate, it is enough that one antenna in one pair be at a distance of order $r\gtrsim c/\omega_{0}=\lambda/2\pi$, i.e. of the order of the central wave-length $\lambda$ with respect to at least one antenna of the other pair. For standard parameters, this is of the order 10 cm, neglecting factors of order 1 (the $2\pi$ helps of course, but for the imaginary part of $C({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{1})$ there is a comparable factor in the scale). Extending the separation of the two antennas in the original pair to $2\Delta r=200$ m, one would have place for about 2000 antennas in between. This in turn would then allow to be built correlations from $10^{6}$ pairs of antenna, where the antennas in each pair are still separated at least by $\Delta r=100$ m. Considering that averaging of $N$ temperature profiles obtained from $N$ statistically independent correlation functions improves the signal-to-noise ratio of the average temperature profile by a factor $\sqrt{N}$, we can improve the SNR by a factor $10^{3}$. If considering the prefactors of order one, 10 times more antennas can be used, the SNR could be improved by a factor $10^{4}$. But even such a large improvement is not yet sufficient to beat the low SNR of order $\chi^{3/2}\simeq 10^{-5}$. It is quite likely, however, that a displacement of an antenna also in $x$ direction by a distance of the order $\lambda/(2\pi)$ leads to a completely decorrelated correlation function. If so, one might gain another factor up to $10^{3}$ in the SNR by considering quasi-1D antenna arrangements, with a width in $x$-direction of order 10 meters. In the latter case one should then be able (after averaging temperature profiles obtained from some $10^{14}$ correlation functions from that many pairs of antennas) to achieve an SNR of $10^{2}$, and hence a RR of order of a few Kelvin. However, it is obvious that the effort for doing so is humungous, and the same geometrical and radiometric resolution might be achievable more easily with other means.

Other interesting ideas of improving the SNR involve using focussing antennas for increasing the flux, and/or exploiting higher order correlation functions as well, but these are beyond the scope of the present investigation.

VIII Discussion

We examined the fundamental feasibility of a new type of passive remote microwave-imaging of a 2D scenery with a satellite having only a 1D antenna array, arranged perpendicular to the direction of flight of the satellite. We analyzed the simplest possible configuration of only two antennas. The scheme is based on correlating Fourier components of the observed electric field fluctuations at the position of the two antennas at slightly different frequencies $\omega_{1}$ and $\omega_{2}$, and leads effectively to a mapping of the 2D brightness temperature as function of position $x,y$ to correlations as function of the center frequency $\omega_{c}=(\omega_{1}+\omega_{2})/2$ and the frequency difference $\Delta\omega=\omega_{1}-\omega_{2}$. With two antennas separated by $\Delta r$, center frequency $\omega_{c}$ and a satellite flying at height $h$, the resolution both in $x$- and $y$-direction is of the order $h\chi=hc/(\Delta r\omega_{c})$. Only very small frequency-differences lead to correlations of finite, useful magnitude. For typical intended SMOS-NEXT values they are of the order of at most 10Hz, which, however, still have to be divided by the number of points in $x$ direction that one wants to resolve within a snapshot. This implies that one must be able to measure GHz frequencies with accuracy of the order of a 1/10-1/100 Hz. The speed $v_{s}$ of the satellite only enters in the maximum frequency difference useful for correlating the signals, which is given by $\Delta\omega\lesssim(\Delta r/h)(v_{s}/c)\omega_{c}$.

In the minimal situation of two antennas, the relative radiometric resolution $\sigma(T)/T$ is, for a single point source of order $\sqrt{\chi}$, whereas for a uniform temperature field in the positive half plane $y>0$ (up to some large cut-off of the size of Earth), $\sigma(T_{\text{rec}})/T_{\text{rec}}\sim 1/\chi^{3/2}$ which for standard parameters is of order $10^{5}$. We have neglected so far the additional noise that comes from the antennas themselves, such that our results should be considered as lower bounds for $\sigma(T_{\text{rec}})/T_{\text{rec}}$. Unfortunately, this large uncertainty prevents a direct application of the method with just two antennas, and massive noise reduction is required. Some ideas are discussed in Sec.VII, where it was found that one can obtain statistically independent correlation functions by displacing one antenna by a distance of the order of $\lambda/2\pi$, where $\lambda$ is the central wave-length. Hence, the signal-to-noise ratio $T_{\text{rec}}/\sigma(T_{\text{rec}})$ can be massively increased by a factor $N$ when using the correlations from $\sim N(N-1)/2$ pairs of antennas from $N$ antennas separated all by at least a distance of order $\lambda/2\pi$. However, the computational effort and the size of the overall structure appear forbiddingly large for achieving a radiometric resolution of order of a few Kelvin with a geometrical resolution of order one kilometer.

An alternative application might be the precise localization of very strong point sources that by far dominate the more or less uniform background from Earth’s thermal emission. As long as one is not interested in a very precise measurement of the intensity of the source, one might localize it very precisely using just two widely separated antennas. These need not even by an board of the same satellite. By having two satellites with well-known distance separated by about 100km for instance, the geometrical resolution achievable in the microwave regime would be of the order of a meter in both $x-$ and $y-$ direction and with rather small computational effort, opening interesting perspectives for such applications. It should also be kept in mind that the method can be easily transferred to other types of waves, sources, and media. For example 2D ultra-sound imaging might be possible by beating the signals of just two moving microphons. Different physical systems can be easily mapped to each other by comparing the corresponding dimensionless parameters introduced in Sec.I.

IX Acknowledgments

The authors are very thankful to the CESBIO SMOS team and the CNES project management, without whom this work would not have succeeded. We have more particularly appreciated the technical support from François Cabot, Eric Anterrieu, Ali Khazaal, Guy Lesthiévant, and Linda Tomasini.

X Appendix

X.1 Current fluctuations and temperature

The connection between the intensity of the current fluctuations and the local temperature can be found e.g. in sharkov_passive_2003 ; rytov_theory_1959 ; Rytov89 ; carminati_near-field_1999 . For being self-contained and relating to the notations used in this paper, we here give a short derivation of this connection. We also show that $\hat{E}_{z,{\bf r}}(\omega)$ is, in the frequency range considered, a circularly symmetric Gaussian process.

X.1.1 Thermal radiation

We begin by recalling the energy density of electromagnetic black body radiation at frequency $\omega$, $u(\omega)=\hbar\omega\rho(\omega)f(\omega,T)$, where $\rho(\omega)=\omega^{2}/(\pi^{2}c^{3})$ is the density of states (number of modes between frequencies $\omega$ and $\omega+d\omega$ per volume), and $f(\omega,T)=1/(e^{\hbar\omega/(k_{B}T)}-1)$ the thermal Bose occupation factor, with $k_{B}$ the Boltzmann constant, and $T$ the absolute temperature of the radiation field. An infinitesimal patch on the surface at position $x,y$ with surface $dA$ and temperature $T(x,y)$ in thermal equilibrium with the radiation field in its immediate vicinity, radiates off an amount of energy per unit time and at frequency $\omega$ given by $dA\,u\,c\cos\theta$ in direction $\theta$ with respect to the the surface normal. The energy density for both polarization directions received at the position of the satellite at distance $R$ from this patch also varies $\propto\cos\theta$, and energy conservation requires

[TABLE]

Earth is rather a grey than a black body, and we therefore have to include the emissivitiy of the patch $B(x,y;\omega,\theta,\varphi)$ in the direction of the satellite given by polar and azimuthal angles. It can also depend on polarization, which we skip here for simplicity. Integration over the whole radiating surface gives the entire energy density at the position of the satellite at this frequency,

[TABLE]

In the microwave regime and temperatures $T\simeq 300K$, $\hbar\omega$ is about four orders of magnitude smaller than $k_{B}T$, such that the Bose factor becomes, to first order in $\hbar\omega/k_{B}T$, $f(\omega,T)\simeq k_{B}T/(\hbar\omega)$, with corrections of order $10^{-4}$. This simplifies $u_{s}$ to

[TABLE]

where we defined the brightness temperature $T_{B}(x,y)\equiv T(x,y)B(x,y;\omega,\hat{k})$, i.e. the absolute temperature a black body would need to have for producing the same thermal radiation intensity at the frequency and in the direction $\hat{k}$ considered, specified explicitly by the two angles $(\theta,\varphi)$. At the same time, the total energy density (integrated over all frequencies) at position ${\bf r}_{1}$ of antenna 1 is $U_{s}=\int d\omega u_{s}(\omega)=\frac{\epsilon_{0}}{2}\langle{\bf E}^{2}({\bf r}_{1})\rangle$, where the average is over the thermal ensemble, but due to ergodicity we may also average over time, $\langle\ldots\rangle_{\tau_{a}}=\frac{1}{\tau_{a}}\int_{-\tau_{a}/2}^{\tau_{a}/2}(\ldots)dt$. In the end one should take the limit $\tau_{a}\to\infty$. In fact, we may even average over both the thermal ensemble and time, i.e. $U_{s}=\frac{\epsilon_{0}}{2}\langle\langle{\bf E}^{2}({\bf r}_{1})\rangle\rangle_{\tau_{a}}$ Expressing then the electric field in terms of its Fourier transform, the time integral leads to a sinc-function,

[TABLE]

with ${\rm sinc}(x)\equiv\sin(\pi x)/(\pi x)$. For large $\tau_{a}$, the sinc-function can be replaced by $(2\pi/\tau_{a})\delta(\omega-\omega^{\prime})$, and we are then left with

[TABLE]

Therefore, the energy density per unit frequency at frequency omega is given by $u_{s}(\omega)=\frac{\epsilon_{0}}{2\tau_{a}}\langle|\tilde{{\bf E}}_{{\bf r}_{1}}(\omega)|^{2}\rangle$. Together with eq.(97) we thus have

[TABLE]

The connection to the current fluctuations is found by comparing this expression to what we obtain from eq.(11) if we do not use (7) yet. There we set $i=j$, ${\bf r}_{1}={\bf r}_{2}=(0,0,h)$, $\omega_{1}=\omega_{2}$, and $v_{s}=0$, as we are interested in the energy density in a given fixed point ${\bf r}_{1}$, identical to the original position of the antenna. This gives

[TABLE]

where $K_{4}=K_{1}^{2}l_{c}^{3}d/(4\pi^{2}\tau_{c})$, and we have already restricted the current density to the surface of Earth, i.e. assumed $\langle|\tilde{j_{i}}({\bf r}^{\prime},\omega^{\prime})|^{2}\rangle=d\langle|\tilde{j_{i}}(x,y,\omega^{\prime})|^{2}\rangle\delta(z)$. In practice, the time integrals originating from the Fourier transforms will be taken over a finite time $\tau_{F}$. Since the only time-dependence is in the exponential, the time-integrals can be done exactly, leading to

[TABLE]

For large $\tau_{F}$, this function is highly peaked at $\omega=\omega^{\prime}$ and behaves as $2\pi\tau_{F}\delta(\omega-\omega^{\prime})$, where the prefactor may be verified by integrating over. We are thus led to

[TABLE]

The thermal fluctuations of the electric field are isotropic in their intensity, such that one third of the energy is in a given polarization direction $i$, i.e. $\langle|\tilde{E}_{i,{\bf r}_{1}}(\omega)|^{2}\rangle=\frac{1}{3}\langle|\tilde{{\bf E}}_{{\bf r}_{1}}(\omega)|^{2}\rangle$. Inserting eq.(100) for the latter quantity, we are led to

[TABLE]

Comparison with eq.(103) allows one to identify

[TABLE]

with $K_{2}=32\tau_{a}\tau_{c}k_{B}/(3\tau_{F}l_{c}^{3}d\mu_{0}c)$ and $T_{\text{eff}}(x,y)\equiv T_{B}(x,y)\cos\theta(x,y,h)$. Thus, the current fluctuations are given directly by the brightness temperature (rescaled by the directional-$\cos\theta(x,y,h)$), up to a constant prefactor. As mentioned in the Introduction, we write $T$ for short for $T_{\text{eff}}$ in the rest of the article. The constant prefactor depends on the time intervals for averaging and the Fourier transforms, but in the end we will always be interested in relative radiometric resolution, i.e. $\sigma(T(x,y))/T(x,y)$, where $\sigma(T(x,y))$ denotes the standard deviation of the reconstructed temperatures over the thermal ensemble of the radiation field, such that the constant prefactor cancels out.

X.1.2 Circular symmetry

A Gaussian distribution of a complex jointly-Gaussian random vector $z=(z_{1},z_{2},\ldots,z_{n})\in\mathbb{C}^{n}$ is fully characterized by the expectation values, $E(z_{i})\,\forall i$, the covariance matrix $K=E[z\,z^{\dagger}]$, and the pseudo-covariance matrix $M=E[z\,z^{t}]$. Both matrices together specify the correlations between the four different combinations of real and imaginary parts of the $z_{i}$. The Gaussian distribution is called circularly symmetric, if $P(z)$ is invariant under the transformation $z\mapsto ze^{i\phi}$ with an arbitrary real phase $\phi$. One shows that a distribution is Gaussian symmetric if and only if $M=0$. This implies immediately that $E[z_{i}]=0\,\forall i$ Gallager08 .

The corresponding definitions and statements for complex Gaussian processes are easily obtained by replacing the discrete index $i$ in $z_{i}$ by a continuous one, e.g. a time argument, or in our case of $\hat{E}_{z,{\bf r}}(\omega)$, a 4-component real vector with a “continuous index” $\omega,{\bf r}$. In order to show that $\hat{E}_{z,{\bf r}}(\omega)$ is a circularly symmetric complex Gaussian process over $\omega,{\bf r}$, we need to prove that $0=M({\bf r}_{1},{\bf r}_{2},\omega_{1},\omega_{2})\equiv\langle\hat{E}_{z,{\bf r}_{1}}(\omega_{1})\hat{E}_{z,{\bf r}_{2}}(\omega_{2})\rangle$, at least in the narrow frequency band that we are interested in. In view of eq.(5), for this it is enough to show that $M_{J}\equiv\langle\tilde{j}_{z,{\bf r}_{1}}(\omega_{1})\tilde{j}_{z,{\bf r}_{2}}(\omega_{2})\rangle=0$. Expressed as Fourier transforms of the time-dependent current-densities, this correlator equals

[TABLE]

The physical origin of the current fluctuations are thermal fluctuations, and the condition of thermal equilibrium implies that the current correlator is invariant under global time-translation (i.e. a shift of the origin of the time axis of $t_{1}$ and $t_{2}$ by the same amount) and hence depends only on $t_{2}-t_{1}$, $\langle j_{z,{\bf r}_{1}}(t_{1})j_{z,{\bf r}_{2}}(t_{2})\rangle=f({\bf r}_{1},{\bf r}_{2},\tau)$, where $\tau=t_{2}-t_{1}$, and we will also use $t=(t_{2}+t_{1})/2$. With this we get

[TABLE]

where $\tilde{f}({\bf r}_{1},{\bf r}_{2},\omega)$ is the Fourier transform of $f({\bf r}_{1},{\bf r}_{2},t)$ with respect to time. The $\delta$-function implies that $M_{J}$ vanishes unless $\omega_{1}=-\omega_{2}$. But we are interested only in frequencies in the small interval $\omega_{2}\in[\omega_{1}-\Delta\omega,\omega_{1}+\Delta\omega]$, centered close to $\omega_{0}$ of the order of 1.4GHz. Hence, in this frequency interval we have indeed $M_{J}=0$, and the complex Gaussian process $\hat{E}_{z,{\bf r}}(\omega)$ over $\omega,{\bf r}$ can be considered as circularly symmetric and eq.(53) valid. In particular, it does not contain any correlator of the type $E\,E$, but only of the type $E\,E^{*}$.