SAR Radiometric Calibration Based on Differential Geometry: From Theory to Experimentation on SAOCOM Imagery

Abstract

The modeling and simulation of topography-induced imaging distortions are crucial for consistent radiometric information exploitation in current and forthcoming SAR-based Earth observation missions with a high spatial and temporal resolution, with relevance in several applications. In this paper, for the first time, we specifically investigate the compensation of topography-induced radiometric distortions affecting SAR images acquired by the L-band Argentinian satellite SAOCOM. We adopt a recently developed calibration method relying on an analytical formulation derived in the rigorous framework of the differential geometry of surfaces. We first provide an original interpretation of the analytical formulation, thus providing further insights into the relevant area-stretching-based formalism. Then, the numerical implementation of the method is specialized to systematically process the data acquired by SAOCOM sensors; hence, the resulting sensor-specific prototype solution processor is employed in this study. Finally, experiments performed over a real scenario in the southern part of Italy, characterized by large topography variations, are presented and discussed, thus elucidating the effectiveness of the adopted method applied to SAOCOM images. The adopted effective SAR calibration strategy opens up the way to its operational use in large-scale SAOCOM data processing.

1. Introduction

There are many possible sources of radiometric distortions in satellite SAR imaging systems [1]. SAR radiometric distortions can be grouped into systemic and scene-dependent ones, both affecting the ability to properly interpret the SAR data information.

Systemic radiometric distortions (e.g., antenna pattern, range spreading, etc.) are intrinsically generated by the SAR sensor structure and can be removed, insofar as quantitative information on the sensor’s system parameters is available.

By contrast, scene-dependent radiometric distortions arise from the inherent 2D imaging process of a 3D irregular ground shape, and they are inherently local in nature [2,3,4,5,6,7,8,9,10,11]. Accordingly, the 2D imaging process of an irregular relief intrinsically implies unavoidable radiometric distortions. This effect may have a significant influence on the radiometric information conveyed by SAR images, and therefore the proper rectification of topography-induced SAR imaging distortions becomes crucial for inverse-scattering model applications and land cover identification and classification [12,13,14].

It is then clear that the availability of effective methods and tools for rectifying the local radiometric distortions affecting SAR images is of crucial importance for converting SAR images into radiometric meaningful quantities, especially if they are supported by rigorous mathematical background with clear meaning and interpretation.

Recently, we have witnessed an increasing interest in the development of SAR systems operating at the L-band. Due to its inherent penetration capabilities, L-band radiation is useful for investigating backscattering from vegetated soils, in order to retrieve accurate estimates of biophysical parameters related to forest structure and biomass. In particular, L-band radar signatures contain important information about crop classification and crop biomass, thus making L-band SAR suitable for monitoring agricultural processes and, in particular, soil moisture [13,14].

In this context, we focus on the Argentinean SAOCOM mission (Satélite Argentino de Observación con Microondas) developed by Comisión Nacional de Actividades Espaciales (CONAE). It consists of two identical Low-Earth-Orbit (LEO) satellites (SAOCOM-1A and -1B), launched on 7 October 2018 and 30 August 2020. Each satellite is equipped with an L-band synthetic aperture radar (SAR) instrument operating at a frequency of 1.275 GHz, with an average altitude of 637 km and a 16 days repetition cycle (thus, the full constellation has a temporal revisit of 8 days). SAOCOM can operate in different acquisition modalities (Stripmap, TOPSAR narrow, and TOPSAR wide swath) and polarimetric modes (single, dual, quad, and compact polarization), and it is equipped with a 3.5 m by 10 m active phased-array antenna. The main goals of the SAOCOM mission specifically address environmental resource management and natural and anthropic hazards (e.g., landslides, oil patches, volcanic eruptions, forest fires, flooding, and coastal pollution) monitoring.

SAOCOM calibration and validation activities were based on external calibration combined with internal calibration, and an ad hoc antenna model was adopted. The external calibration activities aimed at characterizing the calibration constants, which were carried out during the commissioning phase of the SAOCOM-1A satellite, are discussed in [15]. In particular, these calibration activities include absolute radiometric calibration, antenna pattern verification, and noise measurements. Furthermore, an assessment of the SAR product quality carried out during the commissioning phase of SAOCOM-1B is provided in [16]. The abovementioned calibration activities, however, focus on systemic radiometric distortions.

Regarding topography-induced SAR radiometric distortions, some calibration approaches are available in the literature [2,3,4,5,6,7,8,9,10]. More recently, an effective method for radiometric calibration derived in the rigorous framework of differential geometry of the surfaces has been developed, and it exhibits important advantages [10,17]. In particular, as discussed in [10] the method essentially requires image-domain processing based on the evaluation of a scalar function, with important implications in terms of the computational burden. However, to the best of our knowledge, the compensation of local SAR radiometric distortions induced by topography for SAOCOM images hitherto has not been investigated.

Therefore, in this paper, the focus is on the modeling, computation, and compensation of topography-induced radiometric distortions affecting SAR imaging, with special emphasis on the use of SAOCOM data. In particular, this investigation relies on the analytical formulation of the recently developed method proposed in [10], and it is aimed at extending its physical interpretation and applicability. Therefore, the novelty of the present work is twofold, pertaining to both the strengthening of the theoretical understanding of the calibration framework originally proposed in [10] and the implementation and analysis of its application to SAOCOM data. More specifically, the three main contributions of this paper are:

(i)

Providing novel insights into the differential-geometry-based SAR image radiometric calibration method [10], thus clarifying the meaning of its general formulation and, accordingly, providing an original interpretation of the method’s analytical expression, thus shedding new light on the inherent area-stretching function-based formalism [10];

(ii)

Developing a software prototype to systematically process the data acquired by the SAOCOM sensor, thus extending the usability of the original method implementation;

(iii)

Providing a quantitative analysis of the effectiveness of the adopted methodology, through an experimental investigation conducted over a significant mountainous region using the data acquired by the recently launched SAOCOM satellite’s SAR sensor.

The paper is structured as follows. Section 2 introduces the mathematical background underlying the adopted radiometric calibration method. A novel geometrical interpretation of formalism is provided in Section 3. The experimental results obtained by using SAOCAM data are illustrated and discussed in Section 4. Section 5 concludes the paper.

2. Theoretical Background

This section provides an overview of the SAR image calibration method used in this investigation. Specifically, the adopted method is based on the analytical formulation of topography-induced radiometric distortions in SAR imaging that has been recently provided [10] by adopting the theoretical notions of the differential geometry of surfaces [18]. The relevant formulation incorporates a description of the ground surface shape in the 3D object space, by resorting to a suitable surface parametrization process, and a (sensor-centric) cylindrical coordinate system, whose longitudinal axis corresponds to the sensor flight direction. The sensor flight (azimuth) and slant-range directions define the SAR image space. Accordingly, $r$ and $a$ represent the slant-range and azimuth (continuous) coordinates in the SAR image space, respectively. The adopted cylindrical coordinate system suggests the adoption of the following parametric representation of the surface:

$$\mathsf{\phi }:\mathrm{\Omega } {\subset \mathbb{R}}^2\to {\mathbb{R}}^3, \mathsf{\phi }\left(r,a\right)=\left(r\mathrm{c}\mathrm{o}\mathrm{s}\theta ,r\mathrm{s}\mathrm{i}\mathrm{n}\theta ,a\right)$$

(1)

where $\mathrm{\Omega }$ is the parametric domain, and $r$ and $a$ are the radar (parametric) range and azimuth coordinates in the image space, respectively. As a result, a certain location $\mathbf{q}\equiv \left(r,a\right)$ in the SAR image space generally corresponds to a point $\mathbf{P}=\mathsf{\phi }\left(r,a\right)$ of the 3D object space. As a result, the ground surface $\mathrm{\Phi }=\mathsf{\phi }\left(\mathrm{\Omega }\right)$ results to be naturally parametrized in terms of the function $\theta =\theta \left(r,a\right)$, which represents the look-angle function defined in the image space (Figure 1).

Let ${\mathrm{\Omega }}_{i,j}$ be a prescribed rectangle (e.g., a pixel) in the $\left(r,a\right)$-image space, with surface area $\mathcal{A}\left({\mathrm{\Omega }}_{i,j}\right)$, where i and j are the (discrete) coordinates of the pixel along range and azimuth directions, respectively. The corresponding ground surface portion in the 3D object space is mathematically described in terms of the transformation ${\mathrm{\Phi }}_{i,j}=\mathsf{\phi }\left({\mathrm{\Omega }}_{i,j}\right)$; thus, the ground surface area is denoted by $\mathcal{A}\left({\mathrm{\Phi }}_{i,j}\right)$.

The normalized radar cross-section (NRCS) or scattering coefficient of a distributed target is defined as the average radar cross-section (RCS) per unit area on the ground [19,20,21,22]. It is traditionally denoted by ${\sigma }^0$ (sigma naught). The NRCS associated with a prescribed rectangle ${\mathrm{\Omega }}_{i,j}$ in the image domain, denoted as ${\sigma }_{ij}^0$, can be formally written as follows [10]:

$${\sigma }_{ij}^0={\beta }_{ij}^0\frac{\mathcal{A}\left({\mathrm{\Omega }}_{i,j}\right)}{\mathcal{A}\left({\mathrm{\Phi }}_{i,j}\right)}$$

(2)

where ${\beta }_{ij}^0$ is the (observable) radar brightness (also referred to as beta-naught) associated with the pixel of discrete coordinates $i$ and $j$ [23]. It is defined as the average RCS per unit area in the SAR (slant range, azimuth) image space. Note that the area associated with a SAR image pixel, $\mathcal{A}\left({\mathrm{\Omega }}_{i,j}\right)={\mathcal{A}}_{pixel}$, is a constant. To proceed with the energy-preserving reconstruction of the radiometric quantity of interest ${\sigma }_{ij}^0$ in the SAR image domain, the predictable distortive factor $\mathcal{A}\left({\mathrm{\Omega }}_{i,j}\right)/\mathcal{A}\left({\mathrm{\Phi }}_{i,j}\right)$ has to be computed and compensated for pixel by pixel. Equation (2) takes into account the local radiometric distortions in the SAR imaging system induced by the ground topography; thus, this inherent compensation ensures correct normalization and energy conservation.

According to [10], the area $\mathcal{A}\left({\mathrm{\Phi }}_{i,j}\right)$ of the illuminated patch on the ground ${\mathrm{\Phi }}_{i,j}$ can be written in the following form:

$$\mathcal{A}\left({\mathrm{\Phi }}_{i,j}\right)=\underset{{\mathrm{\Phi }}_{i,j}=\mathsf{\phi }\left({\mathrm{\Omega }}_{i,j}\right)}{\iint }dS=\underset{{\mathrm{\Omega }}_{i,j}}{\iint }\mu \left(r,a\right)drda$$

(3)

where the differential quantity $dS=\mu \left(r,a\right) drda$ in (3) describes the elements of the surface area in the 3D object space. According to the metric tensor of the surface, the inherently dimensionless and positive function $\mu \left(r,a\right)$ has been established in the following form [10]:

$$\mu \left(r,a\right)=\sqrt{1+{\left|r{\nabla }_{\perp }\theta \right|}^2}$$

(4)

where the vector ${\nabla }_{\perp }\theta$ represents the gradient of the look-angle function at position $\left(r,a\right)$ in the image space, with

$${\nabla }_{\perp }\theta \left(r,a\right)=\left(\begin{array}{c}\frac{\partial \theta }{\partial r}\left(r,a\right)\\ \frac{\partial \theta }{\partial a}\left(r,a\right)\end{array}\right)$$

(5)

and where the 2D operator ${\nabla }_{\perp }\equiv \left(\frac{\partial }{\partial r},\frac{\partial }{\partial a}\right)$ is given in terms of Cartesian coordinates in the image space.

Therefore, the area-stretching function (4) quantitatively expresses, in a simple and compact way, the inherent topography-induced (local) radiometric distortions, analytically derived in purely geometrical terms within the rigorous framework of the differential geometry of surfaces [18]. According to [10], Expression (4) is established in terms of the gradient magnitude of the look-angle function ($\left|{\nabla }_{\perp }\theta \right|=\sqrt{{\nabla }_{\perp }\theta {\cdot \nabla }_{\perp }\theta }$), uniquely defined in the image domain:

$$\left|{\nabla }_{\perp }\theta \right|=\sqrt{{\left(\frac{\partial \theta }{\partial r}\right)}^2+{\left(\frac{\partial \theta }{\partial a}\right)}^2}$$

(6)

where the functional dependence on the image space position $\left(r,a\right)$ of the look-angle function, its partial derivatives, and its gradient is suppressed and understood.

It is then clear that, according to the adopted formulation, the slopes of the look-angle function, rather than topographic elevation slopes [3], play a crucial role.

Furthermore, it is worth emphasizing that the area-stretching function $\mu$ arises from a metric determinant quantity [18], which is relevant to the spatial transformation between the SAR image space and the ground surface in the observed 3D object space [10]. As a result, $\mu$ describes how the image space is locally stretched as mapped on the ground surface. Conversely, factor $1/\mu$ (with $1/\mu <1)$ describes the local compression (foreshortening) of the elementary ground surface area $dS$ when it is mapped on the corresponding elementary area $drda$ in the SAR image space (Figure 1).

It should be noted that the right-hand side of (3) denotes the integration of the area-stretching function $\mu$ over the area of the elemental domain ${\mathrm{\Omega }}_{i,j}$ (e.g., a pixel) in the image space (see Figure 1). As a result, the formulation in [10] is shown to be mathematically concise and amenable to a straightforward implementation.

It is worth highlighting that radiometric calibration cannot be rigorously obtained in layover regions, and it is meaningless in shadow regions [6]. Accordingly, the calibration operation is applied only to those regions for which a full radiometric compensation is achievable [10], thus discarding layover and shadow regions for which radiometric compensation becomes suboptimal and meaningless, respectively.

A topological remark is now in order. As a matter of fact, the assumption of one-to-one correspondences between the image domain and the ground surface domain does not hold as far as layover and shadow regions are concerned [24]. This characteristic pertains to the different topological nature of layover and shadow phenomena with respect to foreshortening, regardless of the inherent discretization operated by the digital imaging system. It is worth noting that, excluding layover and shadow regions, (1) constitutes a one-to-one mapping, as discussed in [10].

Finally, we address the resulting SAR-calibration processing, which is indeed easily implemented. A conceptual scheme describing the processing in the imaging domain is depicted in Figure 2. A digital elevation model (DEM) of the terrain was used to describe the Earth’s shape, and the geodetic aspects were considered through the use of rigorous transformations. First, geodetic coordinates and DEM elevation values were transformed into 3D coordinates in a geocentric Earth-centered, Earth-fixed (ECEF) reference Cartesian system [25]. Then, according to rigorous SAR geolocation, cylindrical coordinates were derived by using the range-Doppler method [26,27], so that the look-angle function defined in the image space was reconstructed on a regular grid (first box in Figure 2). Then, the computation of the radiometric distortion can be straightforwardly carried out by evaluating the partial derivatives of a scalar function (second box in Figure 2), which were the basis for the evaluation of the ground surface area using Expressions (3)–(6) (third box in Figure 2), and finally, the resulting values were used for radiometric correction according to (2) (last box in Figure 2). More details on the numerical implementation of the method can be found in [10].

3. Geometrical Interpretation of the Area-Stretching Function-Based Formalism

In this section, to provide further insights into the differential-geometry-based method for describing radiometric distortions in SAR images, originally developed in [10] and presented in Section 2, we provide an original geometrical interpretation of the relevant analytical formalism.

In particular, for a better understanding of the geometrical meaning of the area-stretching function (4), it is instructive to analyze how the general formulation is specialized to some meaningful canonical cases. First, the significance of the general analytical expression is illustrated by specializing it to an azimuth-invariant topographic surface case (Figure 3), and subsequently, the idealized flat-Earth case is discussed for comparison (Figure 4).

3.1. Azimuth-Invariant Topography Elevation Case

For the sake of convenience, we consider the ideal case in which the Earth’s shape is azimuth-invariant over the area of interest. According to the considered configuration (see Figure 3), the function $\theta =\theta \left(r,a\right)$ does not depend on the azimuthal coordinate; thus, we have $\theta =\theta \left(r\right)$. In this case, the area-stretching function $\mu \left(r,a\right)$ reduces to the following expression:

$$\mu \left(r\right)=\sqrt{1+{\left(r\frac{d\theta }{dr}\right)}^2}$$

(7)

The function $\mu =\mu \left(r\right)$ in (7) is the radial-stretching factor. Equation (7) has a proper interpretation in terms of the projection on the ground of the elemental radial-distance variation, as clarified in the following.

Let $∆\theta$ be the increment of the look-angle corresponding to the (small) increments (in the imaging domain) along the slant-range and azimuthal coordinates, $∆r$ and $∆a$, respectively. The surface patch on the ground associated with the rectangular element in the $(r,a)$ imaging domain with area $∆r∆a$ is now considered.

First, we estimate the length of the arc on the ground, $∆r_g$, corresponding to a small radial change $∆r$, in the cylindrical coordinate system, as shown in Figure 3. According to the Pythagorean theorem, we assume $∆r_g$ is well approximated by its hypotenuse $\sqrt{{∆r}^2+{\left(r∆\theta \right)}^2}$. Second, we expect:

$$\frac{d{r}_g}{dr}=\underset{∆r\to 0}{\lim }\frac{∆r_g}{∆r}=\underset{∆r\to 0}{\lim }\sqrt{1+{\left(r\frac{∆\theta }{∆r}\right)}^2}=\mu \left(r\right)$$

(8)

where $∆\theta$ represents the angle subtended by $∆r_g$ (Figure 3). According to (8), the ground-range increment $∆r_g$ associated with $∆r$ can then be written to the first order of approximation in the following form:

$$∆r_g=\frac{d{r}_g}{dr}∆r=\mu \left(r\right)∆r$$

(9)

The geometric interpretation of (9) is therefore pictorially shown in Figure 3. Notice that $∆r_g=∆r_g\left(r\right)$ varies across the scene from near to the far range, according to (9). In this special case, the area of the patch on the ground, $∆A_g=∆r_g∆a$, corresponding to a prescribed rectangular element, ${∆A}_{image}=∆r∆a$, in the SAR image domain is consistently given by

$${∆A}_g=\mu \left(r\right){∆A}_{image}$$

(10)

Therefore, the area on the ground ${∆A}_g$ locally represents a stretched version of the corresponding area ${∆A}_{image}$, according to the area-stretching factor $\mu \left(r\right)$ (Figure 3).

In the analyzed case, the geometrical effect of the variability in ground-range resolution, with respect to the nominal resolution of the SAR image, can also be expressed in terms of a nonlinear stretch along the slant-range direction. As a result, this simple geometrical interpretation provides important insights into the general mathematical formulation [10] (see Section 2).

3.2. Canonical Flat-Earth Case

Two simplifying assumptions are here introduced. First, the irregular reliefs are considered absent. Second, the Earth’s shape is considered flat over the area of interest. Accordingly, with reference to the scheme in Figure 4, we have $r\cos \theta =H$, with $H$ and $r$ representing sensor height and the slant range, respectively. In the considered case, the expression of the look-angle function $\theta =\theta \left(r,a\right)$ is specialized as $\theta \left(r\right)=\arccos \left(\frac{H}{r}\right);$ thereby, in this configuration, the look-angle function does not depend on the azimuthal coordinate.

It is straightforward to specialize the expression of the area-stretching function (4) to the considered case. Differentiating $\theta \left(r\right)$ with respect to the slant-range coordinate $r$, we have:

$$\frac{d\theta }{dr}=\frac{H}{r^2}\frac{1}{\sqrt{1-{\left(\frac{H}{r}\right)}^2}}$$

(11)

By using the relation $\frac{H}{r}=\cos \theta$ in (11), after simple manipulations, we obtain:

$$r\frac{d\theta }{dr}=\frac{1}{\tan \theta }$$

(12)

Therefore, by using (12), the area-stretching function $\mu \left(r,a\right)$ in the flat-Earth special case is simplified as

$$\mu \left(r\right)=\sqrt{1+{\left(r\frac{d\theta }{dr}\right)}^2}=\frac{1}{\sin \theta }$$

(13)

The geometric interpretation of (13) is now addressed (Figure 4). Likewise, let $∆\theta$ be the increment of the look-angle function corresponding to the (small) increments along the slant-range and azimuthal coordinates, $∆r$ and $∆a$, respectively. The increment of the look angle $∆\theta$ can be written according to (12), to the first order of approximation, in the following form:

$$∆\theta =\frac{d\theta }{dr}∆r=\frac{∆r}{r\tan \theta }$$

(14)

As sketched in Figure 4, $∆r_g$ denotes the ground-range increment associated with $∆r$, with $∆\theta$ representing the angle subtended by $∆r_g$. By substituting (13) in (9), it turns out that $∆r_g=∆r/\sin \theta$. Accordingly, from (14), $∆r_g\cos \theta =r∆\theta$ is simply obtained. This interpretation is pictorially represented in Figure 4.

The evaluation of the ground surface area, ${∆A}_g$, associated with a prescribed rectangular element in the $(r,a)$ image domain with area ${∆A}_{image}=∆r∆a$ is now in order. In the considered special case, it can be consistently evaluated as follows:

$${∆A}_g=\frac{1}{\sin \theta }{∆A}_{image}$$

(15)

It should be noted that the range-dependent function $\theta =\theta \left(r\right)$ varies across the scene from a near range, ${\theta }_{near}$, to a far range, ${\theta }_{far}$. Hence, the ground area $∆A_g$ associated with the prescribed area ${∆A}_{image}$ varies across the swath according to (15). Therefore, also in this ideal case, the associated radiometric distortion varies across the swath. It is then clear that a suitable radiometric calibration method for SAR images is generally recommended in the flat-ground case, too (see Figure 4).

The intuitive interpretation we provided highlights the importance, in terms of descriptive power and conceptual simplicity, of the differential-geometrical approach used for SAR radiometric distortion modeling and compensation [10].

4. Experimental Investigation

In order to investigate the calibration of SAOCOM imagery, we first developed a SAOCOM-oriented specialized version of the calibration processor implementing the adopted calibration method; then, a representative case study is presented in this section.

4.1. Implementation of a Prototype for SAOCOM Images Processing

In order to apply the adopted calibration method to SAOCOM data, we developed a dedicated numerical implementation, and the resulting prototype was used in this study. Specifically, the specialized prototype was obtained by suitably adapting the original processing kernel for processing SAOCOM data [10,17]. For such purpose, we focused on (1) the data ingestion operation, with the preprocessing subsystem implementing the construction of the dataset of a prescribed meta-sensor, (2) the suitable updating of the data management operation, and (3) the testing and validation of the implemented prototype. We highlight that, in the SAOCOM case, all the metadata associated with the SAR image were provided in XML files, whereas binary files contained the SAR image.

Concerning absolute radiometric calibration, it is worth noting that the SAOCOM data were obtained in a form in which the calibration constant had already been applied, as derived by relying on the reference targets at the calibration site [15].

4.2. Study Area and Dataset

The study area is located in the Calabria region (the southernmost region of the Italian Peninsula), which includes both vegetated areas with significant topographic variations and urban areas (see Figure 5).

In particular, more than 90% of the region is characterized by mountains and hills, including the Sila Mountains, the plateau of Aspromonte, and the massif of Pollino.

The Calabrian region is located between the Ionian Sea and the Tyrrhenian Sea and has almost 800 km of coastline. The reason for selecting this region is twofold. It is characterized by a complex topography. Moreover, it is also particularly important for applications, being characterized by extreme hydrological events, including floods and landslides [28], and is one of the most seismically active regions in the Mediterranean [29]. The experimental investigation was performed considering a dual-polarization SAOCOM-1A (SAO1A) L-band SAR data acquired over the selected mountainous area.

The SLC (single-look complex) data were acquired in the stripmap acquisition mode on 8 October 2021. Figure 5 shows the SAR data footprint superimposed on the Google Earth map of the study area. The nominal resolution was 4.99 m in the azimuth direction, and 5.42 m in the range direction. The main parameters of the selected SAR data are summarized in

Table 1. SAR dataset characteristics.

SAR Platform		SAO1A
Acquisition date		8 October 2021
Observation direction		Right looking
Polarization		VV + VH
Orbit direction		Ascending
Carrier frequency (GHz)		1.275
Off-nadir angle (degree)		32.2174
Sampling frequency (MHz)		30.00
Chirp bandwidth (MHz)		24.40
PRF (Hz)		1857.00
Azimuth bandwidth (Hz)		1229.94
Azimuth-pixel spacing (m)		3.74
Range pixel spacing (m)		5.00
Azimuth resolution (m)		4.99
Range resolution (m)		5.42
Azimuth lines		26749
Range samples		7935
First near	(latitude (deg), longitude (deg))	(39.087114, 16.105089)
First far		(39.967124, 15.893462)
Last near		(40.074031, 16.669855)
Last far		(39.193983, 16.871278)

.

The ground shape description used in this investigation was derived from the digital elevation model (DEM) of the space shuttle radar topography mission (SRTM) with 1 arc-second spacing (i.e., the sampling step of the data product was approximately 30 m at the equator). The world geodetic system (WGS84) datum was used in this study [30].

Accordingly, we preliminarily converted SRTM DEM heights from EGM96 geoid to the WGS84 ellipsoid.

As in common practice, a suitable upsampling was applied to the original SRTM data, thus obtaining a DEM with pixel spacing comparable with the SAR data spacing. This operation was performed only for the benefit of the numerical implementation of the inherent mapping [10], and it did not affect the original data resolution.

It should be noted that the higher the geometric resolution of the original DEM is, the more accurate the estimation of the scattering area is. Nonetheless, high-resolution DEMs might be hard to obtain on a global scale, while SRTM-collected data cover about 80 percent of all the land on Earth.

Regarding the numerical parameters of the computation, we underline that the SAR system’s sampling parameters (see Table 1) determine the pixel size (azimuth × slant-range) in the image domain and hence the pixel size of the simulated radiometric distortion in the same domain. It should be noted that in the present SAOCOM case, the pixel size of the image (4.99 m × 3.74 m) was larger than that of the Cosmo-Skymed sensor (2.33 m × 1.33 m) used in [10]. Based on this parameter, DEM upsampling was arranged accordingly. On the other hand, we highlight that the image size (26,749 × 7935) was the main factor (Table 1) affecting the overall computational cost of the implemented procedure.

4.3. Experimental Results

In this study, the calibration method presented in Section 2 was applied for the radiometric calibration of SAOCOM data, thus specifically demonstrating its usefulness and relevance for sensors operating in the L-band as well.

The topographic pattern of the investigated area is illustrated in Figure 6. The SAR image footprint covered an area ranging between 39.0871°–40.0740° N and 15.8934°–16.8712° E.

The mean altitude of the investigated area was approximately 494 m, ranging roughly from 32 m to 2300 m. As far as scenarios with remarkable topographic reliefs are concerned, the associated large modulations severely affect the SLC image amplitude, and they can easily be compensated for in the imaging domain by specifically using the implemented solution.

In accordance with the method introduced in Section 2, the distortion simulation strategy relied on the computation of the magnitude of the gradient of the look-angle function evaluated in the image domain. First, the pattern of radiometric distortion was simulated, according to the processing scheme described in Figure 2, and then it was used to rectify the SLC image, as detailed in the following.

The look-angle function was evaluated over the swath, ranging from 30.4° to 34.5° (see Figure 7). This uniquely captures the signature of the topography as viewed by the SAR sensor; thus, it concisely conveys all the information needed for computing the correctable SAR radiometric distortions induced by topography. This is a straightforward but central observation, as it renders many complicated calculations unnecessary.

Secondly, according to the processing scheme in Figure 2, the (range-weighted) partial derivative of the look-angle function along the range and azimuth directions in the image space were computed separately.

The magnitude of the (range-weighted) partial derivative of the look-angle function, along the range and azimuth directions, are represented in Figure 8 and Figure 9, respectively.

It is clear that a steeper local look-angle function slope along the range and azimuth directions caused larger distortion. According to (4)–(6), both terms associated with the look-angle partial derivatives contribute to different degrees in determining the amount of local distortion.

As can be seen from Figure 8, the derivative of the look-angle function along the slant-range (radial) direction, $\left|\frac{\partial \theta }{\partial r}\right|$ became negligible as the grazing condition was reached [10]. Figure 9 shows the contribution to the overall radiometric distortion term associated with the magnitude of the partial derivative of the look-angle function along the azimuth direction, $\left|\frac{\partial \theta }{\partial a}\right|$, which became particularly significant over some steep regions (depicted in red in Figure 9). On the other hand, the quantity $\left|\frac{\partial \theta }{\partial a}\right|$ is negligible on the flatlands and vanishes over the sea (flat) area, as can be easily recognized in Figure 9. It should be noted that the depicted quantity is relevant to the partial derivative regardless of its sign, although the look-angle function might increase or decrease along the azimuth direction.

The local incidence angle (LIA) can be also computed according to the formalism in [10], and therefore the function ${\mathrm{\chi }}_l={\mathrm{\chi }}_l\left(r,a\right)$ evaluated in the image space is depicted in the range [0, 90°] in Figure 10.

The simulation of the radiometric distortion pattern was directly performed in terms of a scalar function, i.e., the magnitude of the gradient of the look-angle function.

As seen in the processing scheme in Figure 2, the overall radiometric distortion pattern was finally computed by using (2)–(3), and the obtained result is provided (in a suitable logarithmic scale) according to a suitable color-coded representation in Figure 11.

As can be recognized in Figure 11, the larger radiometric distortions (in red) correspond to higher values of the look-angle partial derivatives (see Figure 8 and Figure 9).

Conversely, at near grazing incidence (i.e., LIA approaching 90°), the look-angle partial derivatives are negligible, and accordingly, the associated local radiometric distortion becomes unitary (in linear scale). Note that over flat regions, the radiometric distortion, however, is not unitary (i.e., it is not zero on a logarithmic scale) and depends on the look-angle as well, in accordance with (13).

An important advantage of the adopted SAR image calibration method resides in its underlying compact analytical formulation, which essentially requires the calculation of the partial derivatives of the look-angle function defined in the image space.

The backscattering coefficient ${\tilde{\mathrm{\sigma }}}_{\mathrm{V}\mathrm{V}}^0$ (in dB) obtained from the entire SLC image acquired in VV polarization is depicted in Figure 12.

We specifically used the symbol ${\tilde{\mathrm{\sigma }}}_{\mathrm{V}\mathrm{V}}^0$ to emphasize that the (pseudo)backscattering coefficient was obtained without applying the compensation of topography-induced radiometric distortions. As can be clearly recognized through a visual inspection of Figure 11 and Figure 12, the simulated distortion pattern in Figure 11 follows the same structure of the ${\tilde{\mathrm{\sigma }}}_{\mathrm{V}\mathrm{V}}^0$ image affected by topography-induced distortions (Figure 12).

Finally, the computed distortion pattern was applied to rectify the imaging distortions induced by the topography (also see Figure 2), and the resulting (calibrated) backscattering coefficient ${\mathrm{\sigma }}_{\mathrm{V}\mathrm{V}}^0$ is presented in Figure 13. A visual inspection of Figure 12 and Figure 13 indicates that the local radiometric distortions were largely rectified.

To better elucidate the obtained results, we now focus on two selected subregions, denoted as R₁ and R₂, both indicated in Figure 12 (with a box with perimeter lines dotted in red).

In particular, the backscattering coefficients obtained for the selected subregion R₁, before and after the radiometric compensation are shown in Figure 14a,b, respectively. The effect of the applied compensation is supported by a visual inspection of Figure 14.

In order to understand the effect of the radiometric compensation in Figure 14 quantitatively, we depict the relevant backscattering coefficient (in dB) versus the local incidence angle (in degrees) obtained before (Figure 15a) and after (Figure 15b) the topography-induced radiometric distortion compensation was applied.

We underline that, in the scatterplots in Figure 15, the number of occurrences (in %) normalized to the total number of image pixels (4,000,000 samples) is taken into account, with the color bars representing the different colors assigned to different classes.

From an inspection of the scatterplots in Figure 15, the effect of the applied radiometric compensation on the overall dependency of the backscattering coefficient on the local incident angle is evident. For both cases (Figure 15a,b), robust linear regression was also estimated by using a least absolute deviation (or L₁-norm regression) technique, which minimizes the L₁-norm of the residuals. Specifically, in Figure 15a,b, we also depict the estimated regression lines (solid black line), with $y=-8.07-0.07x$ and $y=-12.44-0.02x$, respectively. Remarkably, we found that the slope (dB/degree) value reduced from $0.07$ (uncompensated topography-induced radiometric distortions) to $0.02$ after radiometric calibration.

Similarly, the results obtained for the selected subregion R₂ (see Figure 12), which had the same dimensions as R₁, before and after the radiometric distortion compensation was applied are depicted in Figure 16a,b, respectively. Likewise, a representation in terms of scatterplots is given in Figure 17: in particular, the estimated regression lines for Figure 17a,b are $y=7.63-0.097x$ and $y=-11.64-0.051x$, respectively.

As can be seen from Figure 14a and Figure 15a (or Figure 16a and Figure 17a), the distortions associated with the imaging of the underlying topography caused by the effect of scattering ground area indicated important alterations in the radiometric information of the SAR image, which tend to dominate for small local incident angles. Therefore, the applied calibration operation significantly reduced the trend of decreasing backscattering with increasing local incidence angle (see Figure 14b and Figure 15b, or Figure 16b and Figure 17b). It should be noted that the backscattering cross-section per unit area on the ground, derived by removing the artifacts caused by radiometric imaging distortions, remains a function of the local incident angle, i.e., ${\mathrm{\sigma }}^0={\mathrm{\sigma }}^0\left({\chi }_l\right)$, in accordance with the electromagnetic scattering theory [19,20,21,22].

4.4. Remarks on the Backscattering Angular Dependence

It is worth noting that some topographic effects (see, for instance, Figure 13, Figure 14b, or Figure 16b) persisted in the sigma-naught image (compensated for topography-induced radiometric distortions). This is because backscattering coefficients indeed are an intrinsic function of the local incidence angle, which in turn depends on the local topography, as noted previously.

As a matter of fact, the backscattering dependence on the local incidence angle ${\chi }_l$ was initially described according to Lambert’s law, by resorting to simple cosine model-based expression ${\mathrm{\sigma }}^0={\mathrm{\sigma }}_m^0{\left(\cos {\chi }_l\right)}^n$, where ${\mathrm{\sigma }}_m^0$ is the (maximum) backscattering value independent of the incidence angle effect, and $n$ is the power index characterizing surface roughness [31,32]. Few surfaces behave as truly Lambertian surfaces in nature, and accordingly, the related empirical or semiempirical models exhibit limited validity. Nonetheless, backscattering coefficients are related to several parameters describing physical (geometric and electromagnetic) characteristics of the distributed target, including local incidence angle, soil moisture, surface roughness, and dielectric constant.

In the general case, the observed scene can be assimilated to an N-layered structure; thus, its backscattering coefficients might be parametrically expressed as follows [33,34]:

$${\mathrm{\sigma }}_{pq}^0=f\left({\chi }_l,\lambda ,\mathbf{\epsilon },\mathit{t},\mathit{s},\mathit{l}\right)$$

(16)

where ${\chi }_l$ is the local incidence angle,$\lambda$ is the operative wavelength, and p and q are the polarizations of the incident and scattered waves, respectively. In addition, $\mathbf{\epsilon }=[{\mathrm{\epsilon }}_1,{\mathrm{\epsilon }}_2,...,{\mathrm{\epsilon }}_N]$ is the vector of the complex dielectric constant of the different layers, $\mathit{t}$ = [${\mathrm{t}}_1$,${\mathrm{t}}_2..,{\mathrm{t}}_N]$ is the vector of the thickness of the different layers, $\mathit{s}$ = [${\mathrm{s}}_1$,${\mathrm{s}}_2..,{\mathrm{s}}_N]$ is the vector of the rough interface root mean square height, and $\mathit{l}$ = [$l_1$,$l_2..,l_N]$ is the vector of the rough interface correlation length. Resorting to a layered description turns out particularly appropriate insofar as the L-band wave-penetration capabilities are more prominent with respect to higher frequency (C or X) bands.

It is then clear that the radiometric calibration operation and model-based application have to be adopted in an integrated way to proceed with a meaningful quantitative analysis (surface properties retrieval) and interpretation of SAR data, especially when the joint analysis of the heterogeneous multichannel (i.e., multisensor or multiorbit) SAR observations is considered. In this circumstance, the backscattering coefficients, obtained by suitably compensating the SAR images for different radiometric distortions, refer to different local incident angles, and therefore their meaningful exploitation requires the use of suitable parametric scattering models.

5. Conclusions

In this paper, the compensation of topography-induced radiometric distortions was investigated using SAOCOM images. The adopted calibration method, which has been recently developed, relies on an analytical formulation derived from the rigorous framework of differential geometry [10].

This manuscript conveys novelties in several respects, with both theoretical and applicative relevance:

(1)

From a theoretical perspective, we provided an original interpretation of the analytical expressions of the formulation in [10], thus providing further insights into the area-stretching-based formalism;

(2)

The numerical implementation of the method was specialized to process SAOCOM data, with special emphasis on data ingestion operation, meta-sensor data structure assembly, and related management operations. A software prototype was conceived to systematically process the radar data acquired by SAOCOM, and the tested prototype was subsequently used in this study;

(3)

The experimental investigation was conducted by using the prototype processor for SAOCOM image calibration and was supported with illustrations and critical discussion of the obtained quantitative results, thus elucidating the effectiveness of the adopted methodology. Specifically, the experimental results were obtained by using SAOCOM data acquired over a mountainous region in the southern part of Italy;

(4)

The developed prototype provides a useful tool potentially exploitable in all remote sensing applications relying on the SAR amplitude information, thus enabling the operational use of the adopted differential-geometry-based SAR radiometric calibration method in large-scale SAOCOM data processing.

Some considerations on the comparison of the adopted method with respect to other SAR calibration algorithms are in order. Indeed, comparing the results obtained with the diverse existing radiometric calibration approaches [2,3,4,5,6,7,8,9] is not a simple task, insofar as they may include different assumptions, parameters, and specific algorithmic implementations, as noted, for instance, in [35]. Consequently, a consistent and meaningful assessment might be generally difficult to obtain with a superficial direct comparison. Moreover, this difficulty is exacerbated by the lack of SAR-independent reference results, with respect to which none of the existing approaches have been validated hitherto. Nonetheless, a theoretical comparison of the method adopted in this paper with respect to the method based on the projected local incidence angle (or projection cosine approach) [3] is provided in [10], thus demonstrating the local consistency of the two methods.

A fair comparison between different calibration methods would require the availability of an appropriate “ground-truth” reference reflectivity, which is, however, not available for actual images. In fact, working towards the development of such a benchmarking tool, possibly based on system-level simulation, is certainly interesting and an object of current research by the authors.

In future work, we also plan to realize an efficient implementation of the adopted processing strategy for SAR image calibration, which is indeed amenable to be inherently parallelized by employing shared-memory architectures and multithreading [36] or advances in GPU technology [37]. The relevant parallel prototype will boost its operational use in large-scale SAR data processing. In addition, quad-pol data processing will be addressed in future investigations [38].

Another perspective involves an extension of the proposed method to the case of heterogeneous multichannel (i.e., multisensor or multiorbit) SAR observations, in which dissimilar viewing geometries might be prone to misleading interpretation. As a result, the joint use of the information from heterogeneous multichannel SAR observations is not trivial, and therefore the full exploitation of consistently integrated SAR images affected by different distortions still presents unexplored aspects in the domain of remote sensing.