Cross Sensor Simulation of Tomographic SAR Stacks

Abstract

This paper presents an algorithm for simulating tomographic synthetic aperture radar (SAR) data based on another stack actually gathered by a real acquisition system. Through the procedure here proposed, the simulated system can be evaluated according to its capability to image complex natural media rather than reference point targets. This feature is particularly important whenever the biophysical properties of the target of interest must be preserved and cannot be easily modeled. The system to be simulated may be different from the original one concerning resolution, off-nadir angles, bandwidth and central frequency. The algorithm here proposed handles these differences by properly taking into account the wavenumbers of the target illuminated by the real survey and requested by the simulated one. The complex images constituting the synthetic stack are associated with the effective vertical interferometric wavenumber peculiar of the geometry to be simulated, regardless of the original data. Furthermore, the three-dimensional resolution cell of the simulated tomographic system is consistent with the simulated geometry concerning size and spatial orientation. These two latter features cannot be guaranteed by simply filtering the original stack. The simulator here proposed has been used to simulate the tomographic stack expected from the forthcoming European Space Agency (ESA) BIOMASS mission. The relationship between baseline distribution and 3D focusing capability was explored; special attention has been paid to the robustness of tomographic power at being a good proxy for the above ground biomass in tropical regions.

1. Introduction

The design of any spaceborne synthetic aperture radar (SAR) mission poses several challenges, among others the definition of the orbital path or the amount of propellent on board. The vast majority of the system parameters cannot be adjusted after the launch, so their tuning cannot take advantage of the feedback provided by the first acquisitions. This issue, together with the great investment of resources, drives the attention on the reduction of hazards, namely on derisking. Usually, both technological choices and image formation can be properly modeled or simulated. End-to-end descriptions of the acquisition system are developed and can provide the expected output given any input target. Recently, several algorithms have been recently developed for this purpose. Efforts have been made to reduce the computational burden of the simulator especially whenever non-straight orbits are considered [1,2].The electromagnetic features of the structures to be simulated have been studied too, geometrical optics, physical optics and full wave approaches have been proposed [3,4]. Furthermore, ground topography and interferometric processing have been explored [3,5]. These approaches are particularly effective whenever the reference targets are simple (their response is almost known) and the processing chain, up to parameters retrieval, is consolidated. These two conditions are not met in the case of complex natural media and tomographic [6,7] measurements. Targets such as snowpacks [8], glaciers [9] and unmanaged forests [10,11,12] are currently under investigation and any model takes the risk of masking meaningful features. In particular, the analysis of dense forests is the main driver of the forthcoming ESA BIOMASS mission [13] whose launch is expected in 2022. The penetration capability of its P-band signal is expected to provide an insight of many environments as deserts and glaciers in addition to any forests (refer to [14] for further details). Acquisitions on vegetated areas, most importantly on tropical regions, will be processed to generate a map of Above Ground Biomass (AGB) at global scale. Tropical forests are very hard to reduce to a small number of significant elements, moreover the relationship between AGB and radar backscatter is not yet fully understood [15,16]. As a consequence, assessing the impact of any change in BIOMASS system parameters on its capability to retrieve the AGB amount is not trivial and requires a specific strategy. The method here proposed relies on the exploitation of real tomographic SAR images gathered on forests by airborne systems for obtaining a description of the scene; this latter being recovered by means of the back-projection [17] algorithm. Then, the 3D reconstruction of the forest is projected into a tomographic stack in accordance with the geometry of acquisition and radar parameters of the spaceborne system to be simulated. This simple approach requires special care due to the different wavenumbers that gave rise to the airborne radar echo and that are illuminated by the simulated spaceborne system. Also, the three dimensional resolution cells of the two tomographic systems differ as regards their size (dependent on the bandwidth) but also their orientation because of the different acquisition geometry. This latter issue in particular cannot be accounted for by simply filtering the real airborne data stack.

The simulator here proposed has been used to estimate the spaceborne tomographic reconstruction provided by the BIOMASS system. Particular attention has been paid to the power fluctuations due to inaccurate orbital control, that is to irregular baseline sampling. This fluctuations impact on the capability of the BIOMASS system of producing measurements to be used as proxy for AGB estimates. The deterioration of the correlation with AGB of tomographic derived quantities could be observed as a function of orbital inaccuracies.

This paper is organized as follows: in Section 2 a short description of tomographic SAR imaging in the wavenumber domain is presented, highlighting the differences between airborne and spaceborne systems; in Section 3 the procedure for simulating tomographic stacks starting from real acquisitions is detailed (The simulation strategy described in this paper was used to simulate BIOMASS data on boreal and tropical forests that were delivered to European Space Agency (ESA) as a part of the SARSIM database [18]); Section 4 illustrates some features of the synthetic tomographic stack; Section 5 presents the analysis of the impact of orbital inaccuracies on AGB estimates; in Section 6 conclusions are drawn.

2. Wavenumbers in SAR Acquisitions

This section recalls the basic principles of any SAR surveys. The transmitted signal, the scene under observation and their interaction are described in two conjugated domains. The former is the direct domain, that is time or space; the latter is the frequency domain, usually referred to as wavenumber domain when spatial coordinates are transformed. This representation is particularly useful due to the possibility of interpreting SAR acquisitions in terms of enlightenment of the target’s components. Also, the simulation procedure, core of this work, can only be described by resorting to this approach.

The change of domain is achieved by computing the Fourier transform [19] of the function of interest, that is projecting it on a set of complex exponentials of the form shown by Equations (1) and (2) for signals defined over time and spatial maps respectively.

$$\begin{array}{ccc}\hfill G{e}^{j\omega t}& =& e^{j2\pi ft}\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill e^{j{\mathbf{k}}^{*}\mathbf{r}}& =& e^{j\left(kx·x+ky·y+kz·z\right)}\hfill \end{array}$$

(2)

The basis function expressed by Equation (2) has unitary magnitude and phase varying along a specific direction only. This direction is defined by the direction of the real vector $\mathbf{k}$; by moving on the plane perpendicular to $\mathbf{k}$ no changes occur. For this reason, vector $\mathbf{k}$ is referred to as wavevector, its components $kx$, $ky$ and $kz$ are known as wavenumbers. The speed of the phase rotation is determined by the magnitude instead: the longer $\mathbf{k}$ the faster the oscillations. Due to the orthogonality of two complex exponentials having different wavevectors, the function in Equation (2) defines a point in the $\{kx,ky,ky\}$ domain, that is the wavenumber or the spatial frequency domain. Signals traveling in both space and time can be described either in the frequency domain or in the wavenumber domain by freezing time or by focusing on a specific spatial coordinate; the relationship between the propagation in time and in space is ruled by the propagation speed, the speed of light in vacuum c is assumed throughout all this discussion.

The transmitted radar pulse can be described in the temporal frequency domain through a function defined on a limited spectral support: it is centered on the carrier frequency $f_0$ and width is given by the bandwidth B. In the far field approximation, this very same signal can be described in the spatial domain as a plane wave. According to this approximation, the direction of propagation, that is the direction of the wavevector $\mathbf{k}$, is defined by the line joining the sensor and the point (P in Figure 1) far away where the signal is evaluated. The extended bandwidth implies the existence of many waves of the form expressed by Equation (1), each one oscillating at a different frequency. Then, the corresponding spectral support in the wavenumber domain is a line directed as $\mathbf{k}$, with length proportional to the system bandwidth and distance from the origin linked to the carrier frequency. It follows that the link between temporal and spatial frequencies depends on the geometry of acquisition, in particular on the look angle (or off-nadir angle) specified by $\vartheta$ in Figure 1.

The interaction between wave and target generates a set of waves scattered in all directions. In this discussion only the waves traveling back to the radar are considered, that is, only the monostatic case is analyzed. The two-way travel path implies that the wavenumbers of the target responsible for the backscattered echo are twice the ones of the signal impinging on the target itself [20]. Then, the spectral components giving rise to the radar echo, lie on a line centered on $k_0$ and with length $\Delta k$; these latter quantities are defined by Equations (3) and (4).

$$\begin{array}{c}\hfill k_0=\frac{4\pi }{c}{f}_0\end{array}$$

(3)

$$\begin{array}{c}\hfill \Delta k=\frac{4\pi }{c}B\end{array}$$

(4)

The dashed line in Figure 2 is an example of the spectral components illuminated in a highly localized target by a radar in fixed position. The length of this line is connected to the resolution that the radar system can provide along the line-of-sight. The capability of resolving target along the azimuth direction typical of SAR systems can also be interpreted in terms of wavenumbers. A set of acquisitions gathered along the platform trajectory are characterized by spectral lines whose distance from the origin and length are constant (as they depend on radar parameters) but change orientation (being the same as the line-of-sight). The spectral region illuminated by a SAR survey is shown in Figure 2 as a grey surface. The extension of this patch is linked to the resolution that the SAR system can provide; the 2D extent of this surface reflects the lack of resolution of SAR systems along the direction orthogonal to the line-of-sight and to azimuth. That is, standard SAR systems cannot resolve targets lying along the cross range direction $\zeta$ depicted in Figure 1.

The grey surface in Figure 2 is the spectral enlightenment associated with a single (or a narrow range) look angle $\vartheta$; the corresponding surface for a look angle $\vartheta +\delta \vartheta$ can be obtained by rotating this patch around the $kx$ axis by an angle $\delta \vartheta$.

The set of all these surfaces defines a volume in the $\{kx,ky,kz\}$ domain, its edges are shown as thick lines in Figure 2. The spectral spread along the cross range direction should not suggest the possibility of resolving targets in the 3D space though. As a matter of fact, this spectral spread only exists if different parts of the scene are considered together, it disappears for a single target. The very same feature can be appreciated by analyzing the spectral enlightenment along the azimuth direction; in order to do this the constraint of local analysis must be relaxed and the waves must be considered as spreading spherically from the radar rather than being plane. For any fixed azimuth position occupied by the sensor, a large spatial region is reached by the transmitted signal, its extension along x depending on the antenna beamwidth. It follows that many wavenumbers along $kx$ are illuminated at once; however, they are associated with different targets on the ground thus providing a narrow band signal (hence no resolution) for each of them. A wideband signal along $kx$ can only be produced by combining the echoes gathered by the sensor along its path, that is building a synthetic antenna. Analogously, many images of the same target are needed if resolution along the cross range is desired, each one gathered from a different look angle. This issue is discussed more in detail in the following subsection.

2.1. Wavenumbers in SAR Tomography

For the sake of simplicity, the discussion in this section hypothesizes the absence of a squint angle in the antenna pointing and in the focusing algorithm. Holding this condition, after azimuth compression this direction can be considered as decoupled from the other two. For this reason, the following analyses are carried out considering only the zero Doppler plane shown in Figure 1 rather than the whole three-dimensional space. This plane hosts many targets whose common feature is the position of the center of the synthetic antenna, this point also lies on the zero Doppler plane and it is often synthetically referred to sensor position. Whenever the area is imaged several times, each target can be associated with many sensor positions: a proper acquisition plan and processing steps ensure that every sensor position lies on the zero Doppler plane. One of them shall be chosen as reference sensor in order to define the nominal geometry (range axis, look angle and so on), this favored sensor is referred to as master whereas the others as slaves. The definition of a nominal geometry enables to parametrize the zero Doppler plane by means of a pair of orthonormal unit vectors directed as the line joining master and target (slant range) and cross range as shown in Figure 1. As opposed to a $\{y,z\}$ coordinate system, $\{r,\zeta \}$ rotates with the range coordinate however it always highlights the directions along which a single SAR image has resolution capabilities (slant range) and where it lacks (cross range).

The vector joining any two sensors is known as baseline, it can be conveniently expressed as a sum of its components along the slant range and cross range; the former is referred to as parallel baseline whereas the latter as normal baseline. The larger the baseline the further the two sensors, however only the normal component is responsible for a change in the look angle associated with the two positions. The existence of a difference in the look angles implies different directions of arrival of the wave when it impinges on the target. This, in turn, leads to different wavenumbers of the same target illuminated by the two sensors. The geometrical locus in the $\{ky,kz\}$ plane associated with the tomographic acquisition is a set of segments [21], each of them with a different tilt; they are shown in Figure 3 in case of regularly spaced normal baselines and for a nominal look angle of ${45}^{\circ }$. The spread of the spectrum provided by the various acquisitions results in the possibility of resolving targets along the cross range direction as well. The relationship between resolution and spectral width is expressed by Equation (5).

$${\rho }_{\zeta }=\frac{2\pi }{\Delta k_{\perp }}$$

(5)

$\Delta k_{\perp }$ can be expressed as a function of the central frequency and look angles as shown in Equation (6).

$$\Delta k_{\perp }=\frac{4\pi }{c}{f}_0\Delta \vartheta$$

(6)

where $\Delta \vartheta$ is the largest difference of look angles. The order of magnitude of $\Delta \vartheta$ can be appreciated by considering a typical P-band tomographic system whose parameters of interest are summarized in

Table 1. Parameters of a typical P-band tomographic system.

Parameter	Value
Central frequency	400 MHz
Cross range resolution	20 m
Number of images	7

.

By combining Equations (3), (5) and (6) a value of about $1^{\circ }$ for $\Delta \vartheta$ results. It must be kept in mind that the wavenumber region whose width is $\Delta k_{\perp }$, is sampled in the transformed domain thus leading to a periodic signal in the direct domain in case of regularly spaced baselines. The period in the cross range direction is related to the sampling step according to Equation (7).

$${\zeta }_{amb}=\frac{2\pi }{\delta k_{\perp }}$$

(7)

where $\delta k_{\perp }$ is the separation in the wavenumber domain between subsequent passes, it can be expressed as $\delta k_{\perp }=\Delta k_{\perp }/\left(N-1\right)$, N being the total number of passes.

2.2. Wavenumbers in Airborne and Spaceborne SAR Surveys

The much greater distance from sensor to target makes spaceborne SAR surveys rather different from airborne ones. Main differences are due to the significant atmospheric layer through which the signal propagates and to the similarity between near and far range. This section focuses on this latter issue and assumes that the propagation has been properly taken into account thus achieving a good azimuth compression.

Airborne systems usually flight between ∼2 km and ∼10 km of altitude thus resulting in a wide range of look angles within the swath; as a general rule, near range targets are associated with look angles of about ${15}^{\circ }$ reaching ${55}^{\circ }$ in far range, this values being mainly determined by the antenna beamwidth in the elevation direction. On the contrary, the swath width of a spaceborne system shall be set according to the need of avoiding undesired ambiguous echos and nadir returns. These constraints together with the large distances at stake lead to a very small range of look angles; as a rule of thumb, it may change by $3^{\circ }$ from near to far range. These features impact on the spectral region illuminated by the two systems in the $\{ky,kz\}$ domain. Figure 4 shows the wavenumbers excited by an airborne system (in green) and a spaceborne one (in purple). With respect to the airborne, the spaceborne system here considered has a higher carrier frequency but a much lower bandwidth. The dashed line in Figure 4 shows the wavenumbers associated with the carrier frequency as the acquisition moves from near to far range. Its distance from the origin is constant and expressed by Equation (3) whereas its rotation depends exclusively on the geometry of acquisition. According to this, a polar system of coordinates turns out to be particularly helpful in this context; the relationship between the two coordinate systems are expressed by Equations (8)–(10).

$$\begin{array}{c}\hfill \mathbf{k}={\left[\begin{array}{cc}ky& kz\end{array}\right]}^T\end{array}$$

(8)

$$\begin{array}{c}\hfill ky=\left|\mathbf{k}\right|sin\left(\vartheta \right)\end{array}$$

(9)

$$\begin{array}{c}\hfill kz=-\left|\mathbf{k}\right|cos\left(\vartheta \right)\end{array}$$

(10)

The wavevectors pointed out by the dashed line in Figure 4 obey to the relationship expressed by Equation (11).

$${\left(\frac{4\pi }{c}{f}_0\right)}^2=k{y}^2+k{z}^2$$

(11)

A comparison between the impulse response functions (that is, the point spread functions) of airborne and spaceborne systems is shown in Figure 5 for difference range values.

3. Cross Sensor Simulation

This section illustrates a procedure for the simulation of a tomographic SAR stack starting from another one with possibly different carrier frequency, bandwidth and geometry of acquisition. Quantities associated with the first system will be labeled with a ${}_S$, the second with a ${}_A$; this follows the practical need of simulating a spaceborne system starting from an airborne one although the algorithm here illustrated can be applied to any systems pair.

The algorithm here presented relies on the exploitation of Equation (12) [22,23,24] for each azimuth position.

$$I\left(r\right)=\int \int s\left(y,z\right)g\left(r-R\left(y,z\right)\right)e^{\frac{4\pi }{\lambda }R\left(y,z\right)}dydz$$

(12)

where $s\left(y,z\right)$ represents the complex reflectivity distribution in the $\{y,z\}$ plane, $g\left(·\right)$ is the waveform envelope, r is the range coordinate and R is the distance between each point scatterer and sensor. Equation (12) models the generation of any range line of a SLC SAR image: the integral expresses the coherent sum of the scattering coefficients associated with targets placed at the same distance from the sensor, the complex exponential describes the demodulation of the target spectrum in accordance with the radar central frequency and the convolution represents the filtering due to the limited signal bandwidth. The knowledge of system parameters and acquisition geometry leads to two unknowns only in Equation (12), they are the scatterers’ distribution $s\left(y,z\right)$ and the SAR image $I\left(r\right)$. The availability of a stack of SAR images allows SAR tomography to be implemented, through which an estimation of the reflectivity map can be obtained. On the other hand, by feeding Equation (12) with any scatterers’ distribution, a synthetic SAR image can be straightforwardly produced. The generation of a synthetic stack starting from real tomographic data relies on these two steps. Firstly, the reflectivity profile is estimated by performing SAR tomography on real data and then it is used, along with the parameters of the system to be simulated, to produce the synthetic images. However, due to the difference between the two acquisition systems, an extra processing step is needed to make the reconstruction provided by the first compatible with the second.

The complex exponential in Equation (12) is responsible for a demodulation of the spectrum of the target before being filtered by $h\left(r\right)$; as a consequence, only the passband components of $s\left(y,z\right)$ contribute to $I\left(r\right)$. However the spatial direction along which this demodulation takes place is ruled by the range axis, whose direction changes from near to far range. Or rather, the wavevectors of $s\left(y,z\right)$ that contribute to $I\left(r\right)$ have constant magnitude for each value of r but they change direction in accordance with the look angle. It follows that $s\left(y,z\right)$ should lie on the spectral support defined by the acquisition system to be simulated for generating correct range lines. The spectral properties of the target are properly taken into account by the back-projection tomographic algorithm: the reconstructed profile is band-limited and pass-band. However, it lies on the spatial frequencies illuminated by the first system which might be different from the ones required by the second.

The reflectivity map returned by the back-projection algorithm when applied to real data will be referred to as $s_A\left(y,z\right)$. This function, if fed directly into Equation (12) would lead to incorrect results the more the simulated system is different from the real one. The convolution expressed in Equation (12) amounts to a product in the wavenumber domain between the Fourier transform of the system filter (whose support corresponds to the spectral enlightenment) and the reflectivity spectrum of the target. Then, for an acquisition to be properly simulated, these two spectral regions must overlap for all the range line, that is, for every look angle. The overlap must occur in the $ky,kz$ plane so the reflectivity spectrum must match the spectral enlightenment concerning both angle and central frequency. Typically, the angular agreement between the two systems is reached for a narrow interval of distances sampled by the airborne system. It follows that feeding $s_A\left(y,z\right)$ directly into Equation (12) would result in SLCs almost completely empty as a consequence of the different spectral support of $s_A\left(y,z\right)$ and the system filter. Moreover, a large discrepancy in the central frequencies (or a relatively small bandwidth) gets the two spectral regions separated regardless of the look angle. In order to guarantee the correspondence of wavenumbers, the estimated reflectivity profile $s_A\left(y,z\right)$ must be spectrally translated around the wavenumbers illuminated by the simulated system; in the remainder of this section, it is illustrated how this translation can be achieved through a spatial demodulation followed by a modulation.

The reflectivity profile at radio-frequency $s_A\left(y,z\right)$ can be taken to baseband (in the $\{ky,kz\}$ domain) by compensating the phases associated with the local distances. The spatial map demodulating the reflectivity profile can be expressed as in (13).

$$F_{demod}\left(y,z\right)=e^{-j\frac{4\pi }{{\lambda }_A}{r}_A\left(y,z\right)}$$

(13)

where $r_A\left(y,z\right)$ represents the distance of each pixel of $s_A\left(y,z\right)$ from the master sensor and ${\lambda }_A$ the central wavelength used by the airborne system. The product between $s_A\left(y,z\right)$ and $F_{demod}\left(y,z\right)$ modifies the spectral support of the reflectivity function and takes it into baseband. Locally, that is for every look angle, it amounts to a translation of the corresponding spectral line along the direction defined by the spectral line itself; however, different look angles are associated with spectral lines with different orientations so that the translation always follows the radial direction. As a consequence, the shape of the spectral support of $s_A\left(y,z\right)$ is modified: from a section of a ring to a “butterfly” shaped region as shown in Figure 6. The thickness of the butterfly around the origin in the direction orthogonal to the line-of-sight is connected to the resolution of the reflectivity map along the cross-range direction that the airborne system is able to provide. Real spectra before and after the spatial demodulation are shown in Figure 7. The reflectivity profile can be taken in the surrounding of the wavenumbers illuminated by the spaceborne system by another modulation; the modulating map is shown in Equation (14).

$$F_{mod}\left(y,z\right)=e^{j\frac{4\pi }{{\lambda }_S}{r}_S\left(y,z\right)}$$

(14)

where $r_S\left(y,z\right)$ represents the distance of each pixel of $s_A\left(y,z\right)$ from the master sensor of the spaceborne system to be simulated and ${\lambda }_S$ its central wavelength. The generation of the reflectivity map whose spectral support is suited for the projection into a stack of spaceborne images, is shown in Equation (15).

$$s_S\left(y,z\right)=s_A\left(y,z\right)·F_{demod}\left(y,z\right)·F_{mod}\left(y,z\right)$$

(15)

$s_S\left(y,z\right)$ can be fed into Equation (12) to generate the tomographic stack of the system to be simulated. This latter can be explored with the back-projection tomographic algorithm and the 3D reconstruction compared with the corresponding airborne; this comparison is shown in Figure 8.

4. Features of the Simulated Stack

The correct interpretation of the simulated stack of images is not trivial; some features can be understood by exploring the processing chain outlined in Section 3. Again, the real system will be referred to as “airborne” whereas the simulated one as “spaceborne” although this layout is not mandatory. Also, this paragraph focuses on the impact of the $kz$ distribution only rather than both $ky$ and $kz$ due to the 1D nature of the issues here raised.

In general, the tomographic airborne system makes available some wavenumbers of the target whereas the spaceborne system requests others. Each image of the airborne system samples a specific $kz$ so that a sampled version of the Fourier transform of the vertical reflectivity profile of the target is collected. It follows that the reconstruction of the reflectivity profile, obtained for example through back-projection, is periodic. This periodicity implies that its Fourier transform exists exclusively in correspondence of the original $kz$ values, it goes to zero for every other value. Should the spaceborne system require a new $kz$ then some expedient is needed. A possible solution amounts to interpolating the signal in correspondence of the desired $kz$ value; this, in turn, corresponds to windowing the reconstructed reflectivity profile lying in the z domain. It follows that the reflectivity profile may be reconstructed for a limited range of elevations avoiding the replicas. This precaution saves computational time when computing the back-projection and, at the same time, allows the generation of SLC images in correspondence of new $kz$ values. The generation of data in correspondence of very large $kz$ values (outside from the sampled range) is possible only through extrapolation and hence it should be avoided as it amounts to fabricating details.

The aforementioned concepts are strictly valid for regular spatial sampling only; highly irregular sampling of the $kz$ axis in the airborne leads to artifacts in the simulated SLCs. Windowing the reflectivity profile amounts to a convolution between the signal and the kernel of the window, for instance a sinc. Should two subsequent sampled $kz$ be too far away, then the interpolated signal in the middle will be weaker; that is, the synthetic SLC corresponding to a $kz$ far from the sampled ones will be less powerful. Also, due to the high variability of the $kz$ distribution in airborne surveys this power loss is dependent on the spatial coordinate. The synthetic SLCs may exhibit power fluctuations due to the space varying sampling of the $kz$ axis. However, the tomographic profile estimated with the synthetic images is consistent with the one resulting from the original one up to resolution losses. Instability in the $kz$ distribution of the airborne survey leads to fluctuations in the interferometric phase, in the simulated spaceborne stack this latter is exchanged for amplitude fluctuations of the SLC images. This effect is shown in Figure 9.

It must be observed that the simulated stack is not affected by the propagation of the signal more than the real stack is. For a spaceborne simulation to be more realistic the effects of the propagation of the wave through troposphere and ionosphere should be included. However, phase screens may be straightforwardly added to the simulated SLC images; for a even more faithful simulation, they can be defocused and the atmospheric contribution applied to the synthetic raw images.

A final comment on the physical phenomena giving rise to the radar echo is in order. The target features responsible for the backscattered power change with the central frequency of the signal; for instance, shorter waves experience a stronger attenuation and are influenced by smaller structures. The possibility of “changing” the central frequency offered by the algorithm here proposed does not imply the exploration at another scale of the target. The radiometric properties together with the polarimetric signature are determined by the real dataset. The signal gathered at a certain frequency is simply mapped to another. The response of the target when observed at a certain frequency does not provide any indication of its interaction with a signal with different wavelength unless explicitly postulated by certain physical model.

5. Application to BIOMASS: Impact of Orbital Control

The concepts illustrated in the previous sections have been used to simulate the spaceborne tomography that is expected from the ESA BIOMASS mission [13]. The effect of non regularly spaced orbital paths was explored; in particular, their impact on the estimated tomographic power and its correlation to the Above Ground Biomass (AGB). The real airborne stack providing the input have been gathered by ONERA in 2009 over the Paracou area in French Guiana in the context of the TropiSAR campaign [26]; the phases of this image stack were further calibrated according to the procedure described in Reference [27]. The parameters of interest of both systems are listed in

Table 2. Parameters of the P-band airborne radar used in the ONERA TropiSAR survey over Paracou and simulated spaceborne system.

	Airborne	Spaceborne (Simulated)
Number of images	6	7
Wavelength	75 cm	69 cm
Bandwidth	125 MHz	6 MHz
${\vartheta }_{near}$	${28}^{\circ }$	${27.93}^{\circ }$
${\vartheta }_{far}$	${46}^{\circ }$	${28.08}^{\circ }$
$k{z}_{max,near}$	$0.36$ rad/m	$0.264$ rad/m
$k{z}_{max,far}$	$0.25$ rad/m	$0.262$ rad/m

; the scene has been cropped to include only the ROIs with the ground truth measurements of the AGB. The azimuth resolution of the spaceborne system was set to $12.5$ m by means of a simple low-pass filtering. A description of the site providing the ground truth measurements can be found in Reference [28] and references therein whereas the estimation of the above ground biomass is detailed in Reference [29].

The baseline distribution of the BIOMASS system is designed to span the $90\%$ of the critical baseline with 7 subsequent passes regularly spaced, resulting in a vertical resolution of about 20 m. Any orbital drift changes the spatial sampling of the total baseline aperture which, in turn, modifies the vertical Impulse Response Function (IRF) of the tomographic system. The vertical IRF determines the effectiveness of the system in focusing on a specific elevation and rejecting the surrounding echoes. Then, it is expected that any deviations from the nominal orbits perturb the correlation between tomographic power and AGB [15]. In order to evaluate this perturbation the spaceborne system has been simulated several times, each time with slightly different orbital paths. Then, the spaceborne tomography has been generated and the resulting power was put in relation with AGB.

The displacement of each sensor from its nominal position was simulated by means of a Guassian random variable. The displacements involve only the longitude of the orbits, being kept constant the elevation above the Earth’s surface. Also, the total baseline aperture was kept constant and equal to the nominal one, only the five innermost tracks were perturbed. Their fluctuations were modeled by a zero average Gaussian random variable whose standard deviation was increased from $0\%$ to $10\%$ of the critical baseline; for each value of this standard deviation 10 outcomes were generated.

Figure 10 shows the nominal positions of the sensors together with the outcomes; each panel refers to a different perturbation strength. The corresponding vertical IRF is shown in Figure 11; the horizontal axis refers to the elevation whereas the vertical represents the attenuation in dB, the vertical bars refer to $\pm 1\sigma$. It may be observed that driving away from the regular sampling leads to smaller rejection of the echoes coming from surrounding elevations and a larger uncertainty in the exact dB value.

For each outcome of the sensors positions, the tomographic power coming from 40 m above the ground level was estimated and then put in relation with AGB. The choice of 40 m rather than 30 m [15] was driven by the need to take the ground level in correspondence of a zero of the vertical impulse response function (see Figure 11) as suggested in Reference [30]. Many scatterplots as shown by the top panel of Figure 12 were computed and for each one the linear fit estimated. The effectiveness of the tomographic power at being a good proxy for AGB has been assessed by evaluating the correlation coefficient and the sensitivity to AGB; this latter expressed in dB/(Mg/Ha·100) so that a high value suggests more robust estimations. The behavior of the correlation and the sensitivity varying the orbital irregularities is shown in the bottom panels of Figure 12; again, the error bars express $\pm 1\sigma$. As expected, the correlation between tomographic power and AGB drops by moving toward an irregular configuration, at the same time the dispersion around the average value increases. The same happens for the sensitivity, reaching a plateau after an error of about $4\%$ of the critical baseline.

6. Conclusions

In this paper, a new procedure for simulating a tomographic SAR stack was presented; the main free parameters of the system to be simulated are: sensors positions, central frequency and bandwidth. The target of interest comes from tomographic reconstruction based on actual data acquired by a real SAR system. The possibility of simulating SAR images based on real data is particularly convenient whenever the scene is complex and all physical phenomena must be respected rather than modeled. The two systems may be tuned on different central frequencies, have different bandwidths, cover different look angles and provide different vertical resolutions. Special attention was paid to account for all these differences resulting in a simple yet effective procedure. The algorithm here proposed relies on the identification of the spatial frequencies of the target made available by the real system and required by the simulated one. This analysis was made possible by interpreting SAR tomography in the wavenumber domain. The procedure here proposed allows for the desired geometry of acquisition to be properly simulated, in terms of both $kz$ distribution and shape of the three-dimensional resolution cell. It must be observed that these properties cannot be imposed through a simple filtering of the original stack. The features of the resulting stack were discussed as well; in particular, it was shown that each image of simulated stack may be characterized by power fluctuations depending on the instability of the platform providing the original stack.

The simulator here proposed was implemented and used to simulate the tomographic stack that will be provided by the forthcoming ESA BIOMASS mission; the input data coming from the TropiSAR airborne campaign. It was possible to evaluate the relationship between tomographic power and above ground biomass as a function of the baselines irregularities.