Fractal-Based Local Range Slope Estimation from Single SAR Image with Applications to SAR Despeckling and Topographic Mapping

: In this paper, we propose a range slope estimation procedure from single synthetic aperture radar (SAR) images with both methodological and applicative innovations. The retrieval algorithm is based on an analytical linearized direct model, which relates the SAR intensity data to the range local slopes and encompasses both a surface model and an electromagnetic scattering model. Scene topography is described via fractal geometry, whereas the Small Perturbation Method is adopted to represent the scattering behavior of the surface. The range slope map is then used to estimate the surface topography and the local incidence angle map. For topographic mapping applications, also referred to as shape from shading, a regularization procedure is derived to recover the azimuth local slope and reduce distortions. Then we present a new intriguing application of the inversion procedure in the ﬁeld of SAR despeckling. Proposed techniques and high-level products are tested in a wide series of experiments, where the algorithms are applied to both simulated (canonical) and actual SAR images. It is proved that the proposed range slope retrieval technique can (1) provide an estimate of the surface shape, with overall better performance w.r.t. typical models used in this ﬁeld and (2) be useful in advanced despeckling techniques.


Introduction
In recent years, theoretical analyses and experimental proofs highlighted the possibility of extracting more and more information from synthetic aperture radar (SAR) images, obtaining value-added products that help in making SAR systems more accessible and understandable to non-expert users [1]. In addition, the study of extraterrestrial bodies and the search for extrasolar planets is always of great interest, as demonstrated by recent missions such as the Cassini-Huygens mission to Saturn and the ESA planned JUICE (Jupiter Icy moons Explorer) mission [2]. In this context, the estimation of topographic parameters (e.g., local slopes) from (single) SAR images represents an important topic not only per se, but also because from the topographic content of the sensed surface geologists and astrophysicists can derive other important indications regarding both planet activity and future missions' design [3][4][5][6].
Up to now, many techniques aimed at retrieving topographic content from remotely sensed data have been developed, namely stereoscopy [7,8], interferometry [9], radarclinometry-or shape from shading (SfS)- [10][11][12][13][14][15] and polarimetry [16,17], each with its own advantages and drawbacks. Unlike the other techniques, which require (a minimum of) two complex SAR images, possibly acquired at the same time, SfS has the main feature and advantage of requiring just one intensity image. Other existing techniques, as radar and laser altimetry, require no images at all [18,19]. However, they provide one-dimensional (1-D) height profiles rather than two-dimensional (2-D) maps, thus requiring a huge number of passes and some kind of interpolation methods in order to link the different profiles and provide an estimation of the underlying topography. The multiple-pass acquisition system and the relative low accuracy make these techniques not suitable for digital elevation model (DEM) retrieval. Nonetheless, altimetry data are useful in SfS and can be exploited as a priori information.
Since the seminal work by J. van Diggelen [20], and the subsequent works of Wildey [11,21], Horn [12,13], and Frankot and Chellappa [22], SfS techniques have now reached a high level of accuracy, reliability, and robustness within the area of optical imaging, where very accurate SfS algorithms have been developed [13]. Conversely, little has been obtained from SfS techniques applied to SAR images. The main reason for this lack of results stands in the ill-posedness of the problem due to the fact that two variables, namely the local slopes on orthogonal directions, have to be retrieved from a single intensity image. Roughly speaking, SfS techniques have to deal with the retrieval of two unknowns from a single equation. In addition, a large number of physical parameters affect SAR image formation; e.g., the dielectric properties and the small-scale roughness of the observed surface. Therefore, the SfS problem is, in general, under-constrained and simplifying assumptions on the dielectric and roughness behavior of the surface are enforced in the existing literature. In particular, the dependence of image intensity on large-scale slope is mostly modeled using the same Lambertian or Lambertian-like models used for optical images, while the dependence on the physical parameters of the surface (electromagnetic parameters and small-scale roughness) are rarely modeled at all. Therefore, the widespread assumption of uniform dielectric and roughness properties for the surface is not based on a physical reasoning, but simply on the rough guess that the topography provides the main contribution to SAR image formation, at least in scenes where it is significantly present. Actually, without further insights on the nature of the involved dependencies, we are left in the presence of a conceptually unsolvable problem. Hence, the research on this topic reached a point in which the only possibility to obtain new results lies in the introduction of sound direct models for microwave scattering that explicitly depend on the physical parameters of the problem. In this case, even in situations where most of these parameters are unknown, the model can help in gaining a better understanding on some relevant properties of the scene.
In this paper, we propose a novel approach to topographic content retrieval based on fractals, for modeling both the roughness of the observed natural surface and the electromagnetic scattering. In contrast with other existing techniques ( [10,16,17,[23][24][25][26][27])-in which very simple models (or their generalization) appropriate to describe scattering in the optical region are extended in some way in the microwave region or some heuristic approaches are introduced in order to correct them-we propose a complete model-based approach. In particular, starting from a fractal surface model appropriate for natural landscapes, we introduce a direct model able to properly describe the complex electromagnetic phenomena characterizing the scattering from natural surfaces in the microwave region, where SAR systems usually operate. To this end, the problem is decomposed in two steps: first, the surface roughness is described via fractal geometry as a 2-D fractional Brownian motion (fBm) and the Small Perturbation Method (SPM) suitable for such a surface model is adopted as a scattering model to derive the normalized radar cross-section (NRCS). Then, the NRCS is linked to the SAR intensity map through an imaging model. In order to cope with the ill-posedness of the problem and simplify the inversion procedure, we derive a linearized model, which is valid in a small-slope regime and enables us to estimate the range slope. Moreover, as already shown in [28], the first order expansion of the model is shown to be independent on the azimuth slope. An absolute calibration procedure is provided as well.
In the very frequent case in which no a priori knowledge about the surface properties is available, the involved physical parameters can be assumed to be constant over the scene, on the basis of the observation that the main contribution is the topographic one. In actual cases, this hypothesis Remote Sens. 2018, 10, 1294 3 of 33 is unreal and very restrictive; however, in presence of little variation of the characteristics of the surface, the algorithm can still provide acceptable results. Anyway, the model itself does not dictate the requirement of homogeneity of surface characteristics; in fact, whenever some kind of a priori knowledge is available it can be easily injected within the model. For instance, the model can readily include information on the fractal parameters of the scene retrieved from the same SAR amplitude image using the technique proposed in [28], or any kind of information on soil characteristics estimated from multispectral sensors.
The range slope map is then exploited in both standard and innovative applications, namely topographic mapping; i.e., DEM estimation, and SAR despeckling based on a priori scattering information. Concerning the former application, while most state-of-the art techniques are based on the use of polarimetric data (and thus multiple acquisitions) to account for azimuth slopes [16] or simply neglect azimuth slopes [23], we propose an algorithm that employs a regularization procedure to recover the azimuth slopes. In [14,15] the azimuthal component of the slope is estimated through a regularization procedure based on an iterative optimization process. In the present paper we show that, by strengthening the direct model through the introduction of proper direct models for surface and scattering mechanisms, it is possible to estimate the underlying topography also with a very simple and low computational complexity inversion technique. With regard to the other parameters needed by the proposed inversion algorithm, we use reference values or we assume their a priori knowledge, as detailed in the body of the paper. Even if on the Earth's surface more convenient and accurate methods can be enforced for DEM estimation, this application still holds a prominent role in extraterrestrial bodies' exploration, where even a coarser DEM can be useful for geophysical interpretation [3][4][5][6]29].
Besides topographic mapping, it is here demonstrated that range slopes can be conveniently applied to SAR despeckling. To this aim, local slopes are exploited for the estimation of the local incidence angle, which in turn is exploited in the SB-PPB filter, a recent despeckling algorithm based on scattering [30][31][32]. Such algorithms introduce physical insights in conventional despeckling by injecting a priori scattering information about the illuminated surface in pre-existing filters. The required a priori information mainly concerns the local incidence angle, since this is the most relevant factor influencing SAR image formation [30]. In [30][31][32], a DEM gathered from independent sources (e.g., SAR interferometry, stereoscopy, lidar) is exploited to retrieve such a priori knowledge. In this paper, we demonstrate that the topographic content required in scattering-based despeckling can be conveniently estimated directly from the SAR image, thus avoiding the requirement of extra information and data.
The paper is organized as follows. Section 2.1 is devoted to derive an analytical relationship between SAR image intensity and the parameters of interest; i.e., the local slopes of the illuminated surface. The inversion of the direct model is discussed in Section 2.2, and related high-level products and applications (topographic mapping, despeckling) are defined in Section 2.3. Experimental results are presented and discussed in Section 3. First, we apply our framework on simulated data: in this case, synthesized DEMs are generated and used as ground-truth, while corresponding SAR images are simulated with a SAR signal simulator and used as input for SfS. The results presented in Section 3.1 are useful to investigate the validity of model hypotheses. For the other canonical cases, results concerning both DEM estimation and despeckling are reported. A quantitative assessment of the applications is presented. Finally, Section 4 concludes the paper with some remarks and recommendations for future work.

Direct Model
In this section, the rationale for expressing the link between SAR image intensity and the local slopes of the surface is presented. The proposed direct model is divided in three parts: first, the natural surface under study is properly modeled using fractal geometry; in particular, a fBm is Remote Sens. 2018, 10, 1294 4 of 33 used. Actually, to cope with the non-differentiability of the fBm process and the limited range of fractalness of natural surfaces, a smoothed version of the original fBm process must be introduced (physical fBm).
The second step consists in adopting a scattering model; i.e., a method for the evaluation of the field scattered from the previously modeled natural surface. In this step, a relationship between the backscattering coefficient and the geometrical (in particular, local slopes) and electromagnetic parameters of the surface is derived.
The last step consists in the determination of a model for the SAR image, linking the measured intensity map to the backscattering coefficient and to the local slopes of the surface. This complete model-based (surface, scattering and imaging models) approach is a key feature to properly devise an inversion procedure allowing to somehow circumvent the ill-posedness of the problem and physically justify information extraction from even a single SAR image. The proposed model is in sharp contrast with other existing techniques in which only the imaging models, mostly derived from optics, are used [10,16] and only sparse considerations on electromagnetic scattering are made [33].

Surface Model
Several models describing natural surfaces' roughness have been presented. Furthermore, there is increasing experimental evidence that the fractal geometry represents the most appropriate mathematical environment to describe the shape of natural surfaces [34][35][36]; in addition, extensive literature deals with fractal scaling of natural surface roughness [37][38][39]; furthermore, fractal approaches provide a description of the surface with a minimum number of independent parameters [40]. One of the reasons for this success is the ability of fractal models to properly account for the statistical scale-invariance properties (in particular, self-affinity) of natural surfaces. Analytic models reported here are able to handle a variety of natural surfaces: bare and moderately vegetated soils, as well as ocean surfaces. In particular, we adopt here a (topological) 2-D fBm stochastic process z(x,y) defined as follows [40]: where Pr{} stands for "probability", ζ is the considered height increment, z(x,y) is the surface elevation and is the distance between the two considered points of coordinates (x,y) and (x ,y ) and , i.e., the distance over which chords joining points on the surface have a surface-slope mean-square deviation equal to unity.
Fractal geometry exhibits self-affinity properties at any scale and does not allow the derivative operation at any point. Therefore, mathematical fractals cannot be directly applied in scattering problems. However, natural surfaces are observed, sensed, measured, and represented via instruments that are, for their intrinsic nature, band-limited. In other words, no actual natural surface holds property (1) at every scale, and some properties of fBm mathematical surfaces may be relaxed [40].

Scattering Model
The second step is the selection of a scattering model. Here we consider the SPM since it provides analytical closed-form equations, showing a range of validity adequate to SAR applications. Indeed, this method provides a very simple relation between fractals parameters and the backscattered field. Considering a monostatic configuration and assuming that the surface can be described as a physical fBm, the backscattering coefficient under SPM can be expressed as [41] σ 0 mn = 2π8k 4 cos 4 θ|β mn | 2 wherein m and n stand for the transmitted and received polarization (horizontal or vertical), respectively, k is the electromagnetic wavenumber of the incident field. S 0 characterizing the spectral behavior of the physical fBm surface, is expressed in [m (−2−2H) ] and related to T and H [40]; β mn , accounting for the incident and reflected fields polarization, is a function of both the complex dielectric constant ε r of the surface and the local incidence angle θ [40]. Note that, with this model, we are able to deal only with the co-polarized case; in this case, it is possible to consider the term |β mn | 2 constant with θ in the angular interval of interest. Figure 1 shows a graphical comparison between the proposed fractal scattering model and the Lambertian one. The strong non-linearity of the proposed fractal model can be appreciated; unlike Lambertian law, it is also not upper-bounded, thus forecasting an infinite energy scattered in normal incidence conditions. This behavior, clearly in contrast with the physics of scattering phenomena, is due to the approximation used to obtain the SPM solution, which is applicable only in a small slope regime. However, both the scattering models go to zero under grazing incidence conditions. The models are calibrated in order to be both equal to one for a flat surface. Both the scattering models go to zero under grazing incidence conditions; the fractal scattering model is not upper-bounded.

Imaging Model
In this section, we derive a relationship between the SAR intensity image and the slopes of the observed surface in order to apply the SfS algorithm. The following model for the intensity SAR image is used [42]: where I is the SAR intensity image and G is a calibration constant.  (4), we obtain the desired relationship between the intensity map I and the range and azimuth slopes In order to express the intensity image in terms of the local surface slopes, let us introduce the partial derivatives of the surface z(x,y) along ground-range (p) and azimuth (q) directions Hence, it is necessary to link the local incidence angle θ to the local slopes p and q; it can be easily noted that θ is the angle between the propagation unit vector and the surface normal unit vector, so that it can be expressed as: Expressing also sin θ in terms of local slopes, the backscattering coefficient σ 0 mn can be rewritten in the following form: Finally, the SAR image intensity can be linked to the local slopes as follows: Incidentally we note that the direct model presented in (9) clarifies the ill-posedness of the SfS problem (two unknowns, p and q, to be retrieved from (9)) and identifies the numerous parameters (linked to both the sensor and the surface) on which SAR imagery depends. The use of analytical models for both scattering phenomena and surface shape makes us more aware on how these parameters affect SAR image formation, emphasizing the ill-posedness of the SAR SfS problem, and more in general, of any inversion procedure from a single SAR image. At the same time, a proper modeling of the problem gives the chance of circumventing the ill-posedness, making possible the estimation of the parameters of interest. Indeed, the sensitivity analysis shown in [30] demonstrates the major contribution to the SAR image formation due to macroscopic roughness (i.e., to the local incidence angle), whereas the other parameters provide a minor contribution.

Model Inversion
The retrieval of information from remote sensing data of the atmosphere, ocean, and land is basically an inverse problem. It is frequently an ill-posed problem, since many geophysical parameters concur to determine the quantities detected by the sensor (the observations), where different combinations of these parameters can give rise to similar observed quantities. Moreover, errors are always present in the measurements as well as in the knowledge of the causal model (the so-called direct or forward model) that relates the parameters to the observations. Due to the usual lack of a priori knowledge of the characteristics of the observed surface, the proposed inversion method assumes spatial invariance of the dielectric properties of the sensed scene, more precisely the spatial homogeneity of the dielectric constant of the imaged topography. This hypothesis is embedded in the stricter assumption of invariance of the Fresnel reflection coefficient β mn , which is made for the sake of simplicity in order to obtain a less involved inversion technique. Spatial changes of dielectric characteristics (i.e., in general of the complex relative permittivity ε r ) in SAR images can be mainly related to changes in the material that covers surface (for example water, snow, soil, rock, ice, oil, etc.) or to changes in soil moisture content (dry or wet soil).
As described previously in the paper, natural landscapes exhibit a fractal behavior varying both in space and in scale. This means that both H and T may be modeled as a function of spatial and scale variables. However, in order to deal with a simplified forward model and to exploit existing algorithms for fractal parameter estimation, e.g., [28,43], the proposed model assumes a scale-invariance of both the Hurst coefficient and topothesy; i.e., the multifractal behavior of natural surfaces is neglected.
While these problems are relevant to every SAR image and inversion technique, SfS has to deal with another important issue: being available only one intensity image, the inversion procedure has to retrieve two unknowns-p and q-from a single intensity measurement. From these considerations we can understand the successful choice of fractal geometry: it is required the knowledge of only two independent parameters, that, in turn, are independent from radar viewing geometry and sensor parameters (frequency, resolution, look angle) [43].
In this section, the proposed inversion technique is analytically derived, and its limits are analyzed. In order to simplify the inversion procedure, and then retrieve the range slope map, we assume a small-slope regime for the surface. Thanks to this hypothesis, a first-order approximation of the intensity map I, through a McLaurin series expansion with respect to p and q can be considered: where o(·) stands for small-o notation. It can be proved that ∂h ∂q p = 0 q = 0 = 0 [28], so that we obtain a linear function of the partial derivative p only. The intensity function I is, in the first-order approximation, linearly linked only to the range slope of the surface. As a consequence, The coefficients of the McLaurin series expansion, a 0 and a 1 , depend on the specific adopted direct scattering model. In particular, these coefficients are a function of the sensor look angle, which is hence an important parameter for the determination of the validity limits of the proposed linear model. Finally, we note that the obtained result highlights a key property of SAR-and, more in general, of side-looking radars-imaging behavior, showing a clear mathematical definition of a preferential imaging direction due to their particular acquisition geometry.
It can be shown that for the expression in (11) Linearization greatly simplifies the inversion procedure, allowing the estimation of the local range slopes in a very straightforward and fast way.
In order to get a valid slope map, it is necessary to devise a procedure aimed at evaluating the calibration constant G that appears in (11). In this work, we propose a very simple way to estimate G from the data; from (11) we have Since a 0 is known, (13) implies that it is possible to estimate the value of G simply by measuring from data the intensity of a flat homogeneous region. Speckle and roughness could affect the estimation, so it is necessary to consider a sufficiently wide flat area or, alternatively, a despeckling algorithm should be applied to the noisy image to obtain a good estimate of G. Obviously, the main disadvantage of this procedure is that it requires a basic knowledge of the sensed scene in order to identify a flat region in the SAR image.
If such a priori knowledge is not available, assuming the validity of the linear model in (11), the calibration constant G can be calculated as follows where µ I is the estimated mean value of the intensity map I. This procedure, based on the assumption that, according to the fractal model assumed, the global mean slope observed in a SAR image of a natural surface is zero, greatly simplifies the inversion technique, since the knowledge of a flat region present in the intensity image is no longer necessary. Furthermore, this hypothesis drastically reduces the amount of a priori knowledge required: if only the shape of the surface is of interest, no knowledge is required, apart from sensor parameters. Additionally, H can be estimated from data using a linear regression of the range cut of the image power spectra (for more details see [28]). It is noteworthy that the retrieval algorithm proposed in [28] estimates a large-scale Hurst coefficient (i.e., the macroscopic roughness), while, at least in principle, we need a microscopic characterization of the resolution cell. However, thanks to the scale-invariance hypothesis, H is assumed to be constant at any observation scale. In addition, it can be easily shown that both T and ε r do not affect topography estimation if (14) is used, since they simplify in (13). Thus, using (14): This means that, if G is evaluated according to (14), the estimated range slope depends only on the ratio between the coefficients a 0 and a 1 . According to (12), in this ratio both |β mn | 2 (in which ε r is present) and S 0 (in which T is present) cancel out. From now on, in order to simplify the notation, we perform a formal substitution, incorporating the calibration constant within the coefficients a 0 and a 1 , so that, hereafter, the linear imaging model becomes Error Source Evaluation In this section, some considerations about range slope estimation are emphasized and the main error sources are analyzed and discussed. Starting from (11), it is clear that errors on range slope estimation are caused by errors on the linear model coefficients and the calibration constant. Errors on coefficients are linked to model error; i.e., the linear approximation of the intensity image and to parameter errors; i.e., errors on geometric and electromagnetic parameters of both surface and sensor; errors on G are due to model error concerning G evaluation (see Equation (14)) and, as we will also show, to errors committed on a 0 and a 1 .
In the following, we note byx an estimator of x. Considering an additive error ∆G on G estimation we obtain:p Remote Sens. 2018, 10, 1294 9 of 33 Equation (17) means that an error ∆G cause a rescaling and a translation of the estimated range slope and, as a consequence, the estimated topography suffers from a rescaling and a tilt depending on the relative error and the coefficients ratio.
The same effects are caused by an error on model coefficients: where ∆a 0 and ∆a 1 the error on coefficients.
In a more realistic scenario, errors on both calibration constant and coefficients are present; in this case, more severe distortion effects will be experienced, thus: Indeed, (11) suggests that an error on model coefficients causes an error on calibration constant estimation, thus: where µ p is the range-slope expected value and G truth is the true calibration constant. Equation (20) establishes that an additive error on a 0 and a 1 causes a multiplicative error on G, and then an additive error given by The error ∆G takes into account both calibration constant model error ∆G m and coefficients error ∆G c . It can be noted that the two error sources are factorizable, thus Finally, Figure 2 shows at a glance the error source on range slope estimation.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 33 where µp is the range-slope expected value and Gtruth is the true calibration constant. Equation (20) establishes that an additive error on a0 and a1 causes a multiplicative error on G, and then an additive error given by The error ΔG takes into account both calibration constant model error Δ and coefficients error Δ . It can be noted that the two error sources are factorizable, thus Finally, Figure 2 shows at a glance the error source on range slope estimation.

High-Level Products
In this section, we present some meaningful value-added products derived exploiting the proposed range slope estimation method. The first, most obvious, application consists in the retrieving of topography from a single SAR image; i.e., SfS. The proposed approach is based on a two-step procedure and provides as output a raw estimate of the DEM of the imaged scene. The second application is in the context of scattering-based despeckling [30,32]. In particular, we propose

High-Level Products
In this section, we present some meaningful value-added products derived exploiting the proposed range slope estimation method. The first, most obvious, application consists in the retrieving of topography from a single SAR image; i.e., SfS. The proposed approach is based on a two-step procedure and provides as output a raw estimate of the DEM of the imaged scene. The second application is in the context of scattering-based despeckling [30,32]. In particular, we propose to use the retrieved slopes as a priori information, avoiding the necessity of an external DEM, which is usually required by this type of techniques.

DEM Generation
To obtain a DEM, the sole knowledge of range slopes is not sufficient, but it is necessary to estimate also the azimuth slope. As a consequence, DEM estimation requires further processing, whose details are provided in the following. In particular, we propose a two-step procedure: • First, neglecting the local azimuth slope q, it is possible to solve (11) for the unique unknown p, thus obtaining an initial estimate of the height map through slope integration in the range direction. Indeed, according to (11) image amplitude is linearly related to the range component of terrain slope, at least in a first-order approximation. Hence, individual SAR range lines can be independently integrated to generate elevation profiles.

•
The raw guess of the DEM obtained in the first step does not account for azimuth slopes, so that the reconstructed DEM presents unnatural patterns along the range axis, due to independent range line integration. Therefore, the second step consists in a regularization procedure able to correct the integration errors, also taking into account the azimuth component of the slope q.
Since the intensity map is acquired in the cylindrical coordinate system of the SAR sensor, from equation (12) In order to derive the height map in azimuth-ground range, i.e., z(x,y), an azimuth-ground range transformation is required, namely r = y sin θ 0 , where, for the sake of simplicity, flat ground is assumed. Therefore, the surface can be evaluated in the azimuth-slant-range coordinate system as z(m, n) = z(m, n − 1) + p(m, n)∆y, (23) where ∆y = ∆r sin θ 0 is the ground-range resolution. Since ∆y is assumed here to be constant, foreshortening effects will be not corrected, although taken into account, being the DEM estimated in the slant range coordinate system.
It is well-known that the solution of a partial differential equation requires appropriate boundary conditions to provide a unique solution. From a numerical viewpoint, the knowledge of the surface height along an arbitrary curve, intercepting each azimuth line of the image only once, is required. This represents the minimal amount of information necessary to solve (23). In this way, we obtain a first approximation of the topography, consisting in independent range profiles. This concludes the first step.
Once slant-range integration has been performed, a coarse DEM is obtained. Low accuracy is also due to the independent integration of the range profiles, starting from the same starting point, selected neglecting azimuth slopes. In order to overcome this issue, it is necessary to introduce some regularization procedure in the azimuth direction. In fact, a satisfactory DEM can be obtained only linking the different range profiles derived from the first step along the azimuth direction. In contrast with other existing approaches that either use polarimetric data [16] or neglect azimuth slopes [23] or use complex regularization procedures [14,15], we propose a very simple approach for the regularization procedure, based on a Bayesian minimum mean squared error (MMSE) estimation.
In particular, the proposed approach is based on the definition of the azimuth slope MMSE estimator ∆z aMMSE . It can be demonstrated that∆z aMMSE can be obtained as an average (within an appropriate window W k of size K) of the azimuth increments∆z a of the first-step DEM: where the weight function w D (m, n) depends on the statistical characterization and is defined in Appendix A, along with the detailed derivation of (24). Once∆z aMMSE is found, the estimated DEM can be obtained integrating in azimutĥ whereẑ(0, n) is known from first-step integration in y.

Scattering-Based Despeckling
The proposed range slope estimation procedure can be exploited to retrieve the local incidence angle map (through (7)), which provides relevant information about the surface topography and can be used for example in radiometric slope correction algorithms [44]. Another unconventional application of such retrieved topographic information is in the frame of scattering-based SAR despeckling [30][31][32].
Speckle filtering represents a fundamental preliminary step for many remote-sensing applications: this motivates the great interest of the scientific community in the development of new despeckling techniques. In particular, in the last decade, the introduction of non-local approaches represented a breakthrough in SAR speckle filtering. However, most state-of-the-art despeckling techniques are based on statistical and/or geometrical concepts, characterized by limited physical insight. Actually, electromagnetic scattering plays a key role in SAR imaging and in speckle formation. Consequently, more physical-based approaches taking into account the electromagnetic phenomenology can be set. As discussed throughout the paper, the use of physical models is not always an easy task: in this context, the use of machine-learning and neural networks [45][46][47], which are aimed at learning empirical models from data, can represent an interesting alternative to purely model-based approaches.
Following the model-based approach, the authors recently introduced novel scattering-based despeckling algorithms, applying electromagnetic concepts [32] to state-of-the-art non-local filters, namely SB-PPB [30] (based on the PPB filter [48]) and SB-SARBM3D [31] (based on the SARBM3D filter [49]). This kind of approach requires the injection of a priori information about the surface scattering in the despeckling chain. In particular, the algorithms in [30][31][32] are based on the availability of a DEM to be used for scattering evaluation. However, even if DEMs are available on most of Earth's surface, their low resolution can significantly limit the filters' performance [50,51]. As an alternative, in this work, we explore the possibility to fruitfully exploit high-resolution SAR data to retrieve the a priori information needed in the abovementioned scattering-based despeckling algorithms. To this end, the local incidence angle, estimated via (7) assuming q = 0, is used.

Results and Discussion
In this section, experimental results validating the theoretical framework developed in the previous sections are presented. Both simulated and actual SAR images are considered. For SAR image simulation we use SARAS, a SAR raw signal and image simulator [52]. Coherently with the proposed theoretical approach, the backscattering signal is simulated using the SARAS SPM option; furthermore, typical parameters of the COSMO-SkyMed sensor are set. Unless otherwise stated, surfaces with uniform parameters are synthesized. In particular, the same electromagnetic and geometrical parameters have been used both in the direct and inverse models: relative dielectric constant ε r = 3, conductivity σ = 10 −2 S/m, T = 10 −5 m, H = 0.5. According to this setting, the coefficients a 0 and a 1 of (11) are constant over the whole image. Regarding DEM generation, range integration is performed from the mid-range to the borders. Moreover, the canonical SAR images are simulated both with and without speckle in order to better analyze the behavior of the proposed algorithm in specific cases and to isolate relevant effects.
To evaluate the accuracy of the proposed method, we consider both quantitative and qualitative reference measures, using as ground truth the simulated DEM and an available DEM in the canonical and actual cases, respectively. In particular, for qualitative evaluation we show the high-level products previously identified in the paper; i.e., the DEM, the incidence angle map, and the despeckled SAR image. Regarding despeckling, we consider the SB-PPB filter presented in [30] and compare its performances with those of the original PPB filter [48]. The provided experiments are mainly intended at highlighting the advantage of using the SB-PPB filter, which in the version proposed here does not require any additional external DEM, with respect to the standard version of the PPB filter. Therefore, the provided performance assessment leaves aside many well-performing despeckling filters and does not mean to be fully comprehensive. Default filter parameters presented in corresponding works are used for despeckling purposes. Regarding DEM generation, we present range and azimuth profiles of both the reference DEMs and the estimated ones, before (one-step DEM) and after (estimated DEM) the azimuth regularization procedure. In order to evaluate the effectiveness and benefits provided by the azimuth regularization procedure for DEM estimation, in some cases we present also the DEM estimated after the range integration.
For a quantitative evaluation, we consider relevant statistical parameters such as median value, expected value and standard deviation of the absolute error in surface, local range and azimuth slopes estimation and local incidence angle. Quantitative evaluation of the despeckling accuracy is performed through the coefficient of variation Cx, which represents a good quality index in case of heterogeneous and textured areas [53]. In the simulated cases, the coefficient of variation is used as a reference measure; i.e., it has been evaluated on the reference speckle-free image-conversely, in the actual case it has been used as a no-reference measure; i.e., it has been estimated from the single-look image [53]. Run-time RT is also measured in order to quantitatively evaluate the computational burden of the despeckling filters. All the algorithms run on a dual-core workstation equipped with 3 GHz CPUs and 8 GB RAM.
Since geocoding issues are out of the scope of this work and the estimated products are evaluated in the azimuth-slant range coordinate system, in order to take into account the non-regular sampling of the estimated DEMs in the ground range direction (due to foreshortening effects), the comparison between the output of the algorithm and the ground truth is conducted in the SAR coordinate system. Anyway, for visualization purposes the results are shown in an azimuth-ground range reference system assuming flat ground.
For evaluating the impact of a priori knowledge on DEM generation, a comparison between DEMs estimated assuming unknown and known starting points is conducted. A 31 × 2 filtering window is used in the azimuth regularization procedure. When specified, noisy SAR images are filtered via a 10 × 10 incoherent spatial multilook.

Hypotheses Validation
In the development of the technique, the spatial homogeneity of electromagnetic and fractal parameters has been assumed, justified by the major dependence of SAR image intensity on local topography; i.e., on local incidence angle. In this section, we use simulations to assess the impact of the presence of spatial variations of these parameters. In particular, we analyze how variations of the parameters over the scene affect the overall performances of the proposed range slope retrieval algorithm. In order to emphasize and analyze the effects of dielectric and roughness variations over the scene and of their incomplete a priori knowledge, we present two different experiments.
In the first one, we consider a background with fixed electromagnetic (case 1) (Figure 3a-c) or roughness (case 2) (Figure 3d) parameters and generate dielectric (case 1) or fractal (case 2) variations randomly distributed both in location and intensity. In particular, we consider three possible variation patterns: pixel, small patches (20 × 20 pixels), big patches (200 × 200 pixels). In case 1, the background is a standard dry soil; i.e., ε r = 4 and σ = 10 −3 S/m. Electromagnetic characteristics of r% of scene size are then perturbed and dielectric variations are simulated through random gaussian perturbations of the parameters, in particular ε r~N (4,20) and σ~N(10 −3 , 1/2) with a correlation coefficient ρ = 0.7. Pixels with non-physical parameters values (e.g. ε r < 1 and σ < 0) are saturated. In case 2, the background (i.e., macroscopic roughness) is described by H = 0.75 and T =10 −5 m, and fractal parameters of r% of overall scene are then perturbed according to the following: H~N(0.75, 0.1) and T~Log N(−5, 2) with a correlation coefficient ρ = 0.2. Parameter values not belonging to the following ranges are saturated: 0.5 < H < 1 and 10 −10 < T < 0.1. In the inversion procedure, we assume a homogenous background with known parameters, thus neglecting the presence of inhomogeneities. Figure 3 shows the absolute error on the range slope with respect to the percentage (r) of perturbed pixels for all the configurations. The error increases with r and its behavior is independent of the perturbation pattern. This means that the error on the slope does not depend on the spatial distribution of the electromagnetic variations. Moreover, the error behavior is similar in both cases, showing that electromagnetic variations or roughness variations over the scene affect the proposed technique in the same way. Starting from this result, for the second graph (Figure 4) only the small-patches configuration has been considered. In Figure 4 the behavior of the range slope error for a fixed r (r = 0.2) with respect the standard deviation of the topography is shown. In particular, Figure 4a concerns dielectric inhomogeneity, while roughness variation effects are shown in Figure 4b. In order to evaluate the effect of topography on the range slope estimation in the presence of inhomogeneities, the range slope error is normalized to the error committed if a complete a priori knowledge about the parameters' variations is assumed; i.e., to the case in which the electromagnetic and roughness variations are correctly taken into account. As can be seen, this relative error approaches one, graphically demonstrating that for non-flat surfaces, image intensities variations due to electromagnetic or fractal variations are well masked by topography. This behavior shows that it is possible to extract the range slope contribution from a single SAR image and justifies the proposed approach, according to which, in the absence of any a priori knowledge, slopes can be estimated assuming electromagnetic and roughness homogeneity for the surface.  Figure 4) only the smallpatches configuration has been considered. In Figure 4 the behavior of the range slope error for a fixed r (r = 0.2) with respect the standard deviation of the topography is shown. In particular, Figure  4a concerns dielectric inhomogeneity, while roughness variation effects are shown in Figure 4b. In order to evaluate the effect of topography on the range slope estimation in the presence of inhomogeneities, the range slope error is normalized to the error committed if a complete a priori knowledge about the parameters' variations is assumed; i.e., to the case in which the electromagnetic and roughness variations are correctly taken into account. As can be seen, this relative error approaches one, graphically demonstrating that for non-flat surfaces, image intensities variations due to electromagnetic or fractal variations are well masked by topography. This behavior shows that it is possible to extract the range slope contribution from a single SAR image and justifies the proposed approach, according to which, in the absence of any a priori knowledge, slopes can be estimated assuming electromagnetic and roughness homogeneity for the surface.  (c) (d)   Among all the applications of the proposed technique, DEM estimation is surely the most affected by dielectric variations over the scene. Therefore, we now focus on this specific high-level product. To this aim, we generate a flat surface with variations of the complex dielectric constant. In particular, the relative permittivity is set to 4 (dry soil), except in a square where it has been set to 10 (damp soil). Moreover, the imaginary part of the complex dielectric constant has been assumed spatially variant, in particular the electrical conductivity has been assumed equal to 10 −3 S/m except in the aforementioned square where it is set to 10 −2 S/m. In this case T = 10 −5 m has been used. Due to the hypothesis of spatial homogeneity of the surface, the algorithm maps the spatial variations of the intensity induced by variations of ε r into local range slope variations, thus making gross errors in topography estimation; in particular, the errors are proportional to the range of variability of ε r . Therefore, in general, in absence of any a priori knowledge about the electromagnetic characteristics of the surface, an error on the dielectric constant will be committed and this will map into an error on the linear imaging model coefficients and, hence, on the calibration constant. However, the error will not affect slope estimates, which will be unreliable only within the bright square. Conversely, if the variations of ε r are properly taken into account, the SfS algorithm should obviously return a constant DEM. To experimentally confirm this behavior, several types of a priori knowledge are considered ( Figure 5). First, no a priori knowledge is assumed and a reference value for ε r is used to evaluate G, a 0 and a 1 (Case 1), in particular (20) is used for G evaluation. This guarantees independence and robustness of the proposed algorithm against uncertainties in parameters estimation. In order to separate the effects on the calibration constant by those on the expansion coefficient, the knowledge of the background electromagnetic properties is assumed and used to correctly evaluate the calibration constant only in Case 2, and also the coefficients in Case 3; finally, in Case 4 a complete a priori knowledge is assumed and used to correctly evaluate all parameters. As can be seen in Case 1 and 2, higher values of the relative permittivity and of the electrical conductivity determine a bright square in the intensity map, which the inversion algorithm interprets as a sharp change in local slopes. In particular, the algorithm senses a rect-like trend in range slopes so that the estimated DEM presents an inexistent ramp. In Figure 6 the mean absolute slope error against the normalized error on a dielectric constant is shown for Case 2. In particular, this performance parameter has been evaluated as follows: let us sayp(m, n) and p(m,n) is the estimated and actual range slope at pixel (m,n), with (m,n)∈R, being R the region affected by dielectric variation. The mean absolute slope error is defined as where dim(R) is the number of pixels in R.
Besides a slope error increasing with the permittivity error, the graphs clearly show that higher deformations are present in presence of an underestimation of the dielectric constant. In Case 3, thanks to the background a priori knowledge, the surface is correctly reconstructed, with the exception of the bright region, in which a ramp trend is clearly visible. However, if complete knowledge of the electromagnetic characteristics is available, the underlying topography is estimated with high accuracy.
In conclusion, this experiment shows that the spatial homogeneity of the illuminated surface should be assumed to obtain a reliable DEM, unless a deep a priori knowledge is at hand; moreover, if possible, it is desirable to use a priori knowledge regarding a wide region of the scene. Obviously, as shown in Figure 5, the more heterogeneous the electromagnetic characteristics of the surface, the more unreliable the estimated DEM.

Sinusoidal DEM
In this experiment, a sinusoidal DEM is generated and the corresponding SAR image is simulated with and without a Rayleigh-distributed speckle (Figure 7a,b). In Figure 8a-e the estimated DEM is shown and compared with the original one before and after azimuth regularization and with both speckle-free and speckled SAR images; numerical results are reported in Tables 1 and 2, respectively. The proposed algorithm provides a very accurate topography reconstruction, succeeding in estimating and recovering the underlying sinusoidal pattern, even if the original SAR image is affected by multiplicative noise; it is important to emphasize that no despeckling algorithm is applied. This means that, without any knowledge about the starting points used in range integration, the algorithm is able to recover the shape of the original imaged surface, even if not the exact height. The azimuth filter significantly improves the accuracy of the algorithm, especially in presence of speckle, as can be seen in Figure 8d, thus greatly reducing the linage effect due to independent integration in range direction, is clearly visible in Figure 8c. The smoothing effect due to azimuth filtering is especially visible in the azimuth cut provided in Figure 9. Therefore, the azimuth regularization procedure provides significant improvements in terms of both speckle reduction and increased correlation of range profiles.
In the speckle-free case, the effect of azimuth filtering provides no relevant enhancement (Figure 8a,b), even if after the application of the azimuth regularization procedure, the estimated DEM is a bit more accurate than that obtained after range integration only ( Table 1).
The distortion that can be noticed in the estimated DEMs is due to the variable height of the surface along the near range azimuth line used as starting line in the range integration step. Figure 8e shows the estimated DEM after azimuth filtering in absence of speckle and assuming known starting points: the waveform distortion is significantly reduced, see also Table 1. It is noteworthy that opposite properties are desired for DEM accuracy in presence or in absence of geometrical a priori knowledge: the more irregular the known azimuth line, the greater the benefits provided by the a priori knowledge. On the other hand, in the absence of a priori knowledge, the greater the maximum difference in the original topography along the first azimuth line, the more distorted the estimated elevation map in the azimuth direction. Consequently, combining a SAR sensor with altimetry data could significantly improve the performance of the algorithm and much less distortion is expected. Finally, the presented technique is also compared with the Lambertian model and results are provided in Tables 1 and 2. Thanks to the fractal models adopted, a significant improvement of DEM accuracy is obtained.
The local incidence angle maps estimated from speckle-free and multilook SAR images are shown in Figure 10a,b, respectively, and compared with the actual one shown in Figure 10c.
The incidence angle map in Figure 10b is then exploited in SB-PPB to estimate the scattering a priori information. SB-PPB and PPB are applied to the single-look SAR image and results are shown in Figure 7c,d, respectively. Quantitative performance parameters are reported in the figure caption, showing that comparable results are obtained with the two filters. However, it is noteworthy that SB-PPB implies a significant run-time reduction, thus almost halving the computational burden.
Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 33 that opposite properties are desired for DEM accuracy in presence or in absence of geometrical a priori knowledge: the more irregular the known azimuth line, the greater the benefits provided by the a priori knowledge. On the other hand, in the absence of a priori knowledge, the greater the maximum difference in the original topography along the first azimuth line, the more distorted the estimated elevation map in the azimuth direction. Consequently, combining a SAR sensor with altimetry data could significantly improve the performance of the algorithm and much less distortion is expected. Finally, the presented technique is also compared with the Lambertian model and results are provided in Table I and II. Thanks to the fractal models adopted, a significant improvement of DEM accuracy is obtained. The local incidence angle maps estimated from speckle-free and multilook SAR images are shown in Figure 10a,b, respectively, and compared with the actual one shown in Figure 10c.
The incidence angle map in Figure 10b is then exploited in SB-PPB to estimate the scattering a priori information. SB-PPB and PPB are applied to the single-look SAR image and results are shown in Figure 7c,d, respectively. Quantitative performance parameters are reported in the figure caption, showing that comparable results are obtained with the two filters. However, it is noteworthy that SB-PPB implies a significant run-time reduction, thus almost halving the computational burden. that opposite properties are desired for DEM accuracy in presence or in absence of geometrical a priori knowledge: the more irregular the known azimuth line, the greater the benefits provided by the a priori knowledge. On the other hand, in the absence of a priori knowledge, the greater the maximum difference in the original topography along the first azimuth line, the more distorted the estimated elevation map in the azimuth direction. Consequently, combining a SAR sensor with altimetry data could significantly improve the performance of the algorithm and much less distortion is expected. Finally, the presented technique is also compared with the Lambertian model and results are provided in Table I and II. Thanks to the fractal models adopted, a significant improvement of DEM accuracy is obtained. The local incidence angle maps estimated from speckle-free and multilook SAR images are shown in Figure 10a,b, respectively, and compared with the actual one shown in Figure 10c.
The incidence angle map in Figure 10b is then exploited in SB-PPB to estimate the scattering a priori information. SB-PPB and PPB are applied to the single-look SAR image and results are shown in Figure 7c,d, respectively. Quantitative performance parameters are reported in the figure caption, showing that comparable results are obtained with the two filters. However, it is noteworthy that SB-PPB implies a significant run-time reduction, thus almost halving the computational burden.

Fractal DEM
The proposed algorithm is applied to a fBm surface of controlled parameters synthesized using the Weierstrass-Mandelbrot function [40,54,55] (Figure 11). Figure 12 shows the simulated SAR images: (a) speckle-free, (b) with Rayleigh-distributed speckle (single-look).
Estimated DEMs are reported in Figure 13. In particular, in Figure 13a,b the DEMs estimated, in the speckle-free case, after azimuth filtering assuming unknown and known starting points,

Fractal DEM
The proposed algorithm is applied to a fBm surface of controlled parameters synthesized using the Weierstrass-Mandelbrot function [40,54,55] (Figure 11). Figure 12 shows the simulated SAR images: (a) speckle-free, (b) with Rayleigh-distributed speckle (single-look).
Estimated DEMs are reported in Figure 13. In particular, in Figure 13a,b the DEMs estimated, in the speckle-free case, after azimuth filtering assuming unknown and known starting points, respectively, are shown. Figure 13c-f show DEMs estimated, in presence of speckle, before and after the regularization procedure, considering the single-look image (Figure 13c,d) and the spatial multilook one (Figure 13e,f). Performance parameters for DEM retrieval are presented in Tables 3-5.
Surface topography is reconstructed with sufficient accuracy even in the single-look noisy case, where the algorithm identifies most of the main topographic features, especially after azimuth filtering. Obviously, the best result is obtained in the speckle-free case with known starting points. In this case, good azimuth slope estimation is achieved (Table 3). Moreover, the first two rows of Table 3 highlight a very interesting result: the proposed azimuth regularization procedure is not able to compensate the lack of knowledge about range integration starting points. This experiment shows that speckle has a very detrimental effect on range integration, thus greatly affecting the capability of the algorithm to retrieve the 2-D underlying DEM in a satisfactory way. However, azimuth filtering is able to reduce the impact of speckle on both azimuth and range integration (Figures 14b and 15b). It is important to note that speckle also affects the estimation of the calibration constant, which is very sensitive to multiplicative noise.
Even simple spatial multilook despeckling greatly improves the performances of the algorithm, which is able to identify most topographic features, also without azimuth regularization (Figure 13e,f): good results are obtained both in range ( Figure 14c) and azimuth (Figure 15c). Indeed, a much less noisy one-step DEM is obtained owing to despeckling. Finally, we note that the Lambertian model provides a significantly overestimated topography, as shown in Figures 14a and 15a.
Reference and estimated incidence angle maps are shown in Figure 16. This experiment confirms the results obtained in the previous one: the reference incidence angle (Figure 16a) is estimated with great accuracy in the speckle-free case (Figure 16b), whereas noisier estimates are obtained in the presence of speckle.
Regarding the despeckling high-level product, in Figure 12c the result obtained applying SB-PPB, using the estimated incidence angle shown in Figure 16c, is shown. For comparison, results obtained via PPB are reported in Figure 12d: indeed, results are similar, as also indicated by the coefficient of variation reported in the figure caption. This may be also due to the relative smoothness of the simulated DEM, which enables PPB to preserve image details. However, it is noteworthy that the adaptive procedure in SB-PPB allows for a significant run-time reduction, thus halving the computational burden in the selected study case.
Remote Sens. 2018, 10, x FOR PEER REVIEW  19 of 33 respectively, are shown. Figure 13c-f show DEMs estimated, in presence of speckle, before and after the regularization procedure, considering the single-look image (Figure 13c,d) and the spatial multilook one (Figure 13e,f). Performance parameters for DEM retrieval are presented in Table 3-5. Surface topography is reconstructed with sufficient accuracy even in the single-look noisy case, where the algorithm identifies most of the main topographic features, especially after azimuth filtering. Obviously, the best result is obtained in the speckle-free case with known starting points. In this case, good azimuth slope estimation is achieved (Table 3). Moreover, the first two rows of Table  3 highlight a very interesting result: the proposed azimuth regularization procedure is not able to compensate the lack of knowledge about range integration starting points. This experiment shows that speckle has a very detrimental effect on range integration, thus greatly affecting the capability of the algorithm to retrieve the 2-D underlying DEM in a satisfactory way. However, azimuth filtering is able to reduce the impact of speckle on both azimuth and range integration (Figures 14b and 15b). It is important to note that speckle also affects the estimation of the calibration constant, which is very sensitive to multiplicative noise.
Even simple spatial multilook despeckling greatly improves the performances of the algorithm, which is able to identify most topographic features, also without azimuth regularization ( Figure  13e,f): good results are obtained both in range ( Figure 14c) and azimuth (Figure 15c). Indeed, a much less noisy one-step DEM is obtained owing to despeckling. Finally, we note that the Lambertian model provides a significantly overestimated topography, as shown in Figures 14a and 15a.
Reference and estimated incidence angle maps are shown in Figure 16. This experiment confirms the results obtained in the previous one: the reference incidence angle (Figure 16a) is estimated with great accuracy in the speckle-free case (Figure 16b), whereas noisier estimates are obtained in the presence of speckle.
Regarding the despeckling high-level product, in Figure 12c the result obtained applying SB-PPB, using the estimated incidence angle shown in Figure 16c, is shown. For comparison, results obtained via PPB are reported in Figure 12d: indeed, results are similar, as also indicated by the coefficient of variation reported in the figure caption. This may be also due to the relative smoothness of the simulated DEM, which enables PPB to preserve image details. However, it is noteworthy that the adaptive procedure in SB-PPB allows for a significant run-time reduction, thus halving the computational burden in the selected study case.

Actual Case Study
In this experiment, an actual SAR image is provided as input to the algorithm in order to visually and quantitatively evaluate the effectiveness of the proposed technique in retrieving DEMs from actual imagery. To this end, we choose a 1950 × 2430 subset of a stripmap single-look amplitude image of the Vesuvius volcano, Italy, acquired in HH channel by COSMO-SkyMed on 1 August 2011. We tried to consider an area almost free of man-made objects, an optical image of the selected area is shown in Figure 17. The radar-look angle is ϑ 0 = 35 • , the operating frequency is f = 9.6 GHz; azimuth and ground range resolution are both equal to 2.5 m, while the pixel spacing is 2.07 m and 2.06 m in azimuth and ground range, respectively. It is important to note that the variation of the radar-look angle over the scene is negligible (fractions of degree) due to its limited size, so it is assumed constant in the algorithm. The spatial multilook image (Figure 18a) is used for range slope estimation. Reference values for T, ε r and σ have been used, in particular T = 10 −4 m, ε r = 10 and σ = 10 −3 S/m; H has been estimated using the technique presented in [28], obtaining a mean value of 0.8, which has been used as input value.
The estimated range slope map is shown in Figure 18b; Figure 18c,d show the images despeckled using PPB and SB-PPB, respectively. The estimated incidence angle map in Figure 19a is exploited in SB-PPB for the evaluation of the scattering a priori information and is compared with the actual map shown in Figure 19b. SB-PPB provides a better texture preservation with respect to the PPB filter, as testified by visual inspection and by the measured coefficients of variation, reported in the figure caption. The adaptive iteration scheme in SB-PPB also reduces the computational burden, since RT = 994 s for PPB and RT = 293 s for SB-PPB.
Retrieved DEMs are shown in Figure 20a-c and compared with a 5m-resolution aero-photogrammetric DEM used as reference (Figure 20d). The crater of the Vesuvius and, beyond it, the Mt. Somma ripples are identified by the algorithm despite the fact that range integration has been performed from starting points located at the same height (Figure 20a). The effects of range-independent integration are sufficiently reduced once the azimuth filtering is applied and the Vesuvius-Somma complex is clearly visible in the estimated DEM ( Figure 20b).
As previously discussed, lack of knowledge about starting points in the range integration step may cause a significant distortion of the estimated DEM, especially in the azimuth direction. For this reason, Figure 20c shows the DEM estimated from the multilook image assuming the knowledge of the mid-azimuth line. Thanks to this a priori knowledge, the linage effect due to the independent range integration is greatly reduced, especially in the central part of the image. This accuracy improvement can be appreciated also from Figure 21, where range and azimuth cuts highlighted in Figure 20d are shown and compared with both the reference profiles and the Lambertian model. The proposed fractal model obtains better performance, especially in the azimuth direction (see also Table 6). Furthermore, if knowledge about starting points in the range integration step is available, more features are identified by the algorithm as shown in Figure 21c, where the gain in performance due to the knowledge of the starting points is clearly visible (see the red circled area in Figure 21c).
The underestimation of range slopes resulting from the proposed fractal approach requires an additional comment: the optical image in Figure 17 testifies the presence of massive vegetation all around the Mount Vesuvius cone and on Mt. Somma. In this situation volumetric scattering due to vegetation cannot be neglected, especially in strong backslope regions, where it represents the dominant scattering contribution. This volumetric contribution leads to a general increase of the backscattered signal [56] thus causing an overestimation of the calibration constant G and, consequently, an underestimation of range slopes. Volumetric scattering cannot be neglected in SfS techniques applied to SAR images of Earth (at least in X-band) and it would require appropriate modeling [56], but, in a first approximation, it can be modeled to be independent from the local incidence angle and can be treated as an additive constant. This contribution can be estimated from intensity data in a region with negligible surface scattering; i.e., from strong backslope regions of the image. For this purpose, we estimated the required additive constant from the region marked in Figure 20a. Taking into account this volumetric contribution, a significant enhancement of range slope estimation results is obtained, as shown in Figure 22. In particular, the median error on the range slope decreases to 7.42 • for the proposed approach, while it increases to 74.96 • for the Lambertian model. In conclusion, although the strong geometric distortions present in the image-foreshortening, layover on Mt. Somma, and shadowing in the crater of the Vesuvius-the algorithm provides encouraging results. consequently, an underestimation of range slopes. Volumetric scattering cannot be neglected in SfS techniques applied to SAR images of Earth (at least in X-band) and it would require appropriate modeling [56], but, in a first approximation, it can be modeled to be independent from the local incidence angle and can be treated as an additive constant. This contribution can be estimated from intensity data in a region with negligible surface scattering; i.e., from strong backslope regions of the image. For this purpose, we estimated the required additive constant from the region marked in Figure 20a. Taking into account this volumetric contribution, a significant enhancement of range slope estimation results is obtained, as shown in Figure 22. In particular, the median error on the range slope decreases to 7.42° for the proposed approach, while it increases to 74.96° for the Lambertian model. In conclusion, although the strong geometric distortions present in the imageforeshortening, layover on Mt. Somma, and shadowing in the crater of the Vesuvius-the algorithm provides encouraging results. consequently, an underestimation of range slopes. Volumetric scattering cannot be neglected in SfS techniques applied to SAR images of Earth (at least in X-band) and it would require appropriate modeling [56], but, in a first approximation, it can be modeled to be independent from the local incidence angle and can be treated as an additive constant. This contribution can be estimated from intensity data in a region with negligible surface scattering; i.e., from strong backslope regions of the image. For this purpose, we estimated the required additive constant from the region marked in Figure 20a. Taking into account this volumetric contribution, a significant enhancement of range slope estimation results is obtained, as shown in Figure 22. In particular, the median error on the range slope decreases to 7.42° for the proposed approach, while it increases to 74.96° for the Lambertian model. In conclusion, although the strong geometric distortions present in the imageforeshortening, layover on Mt. Somma, and shadowing in the crater of the Vesuvius-the algorithm provides encouraging results.

Conclusions
In this paper a novel model-based inversion technique for retrieving the range slope map from a single SAR image of natural surfaces has been developed. The main innovations are both conceptual and applicative: the use of fractal geometry to model both natural surfaces and scattering phenomena represents the core of the proposed theoretical framework, filling a general gap in the state of the art, where most of the works use Lambertian-like models, inadequate to properly model the complex dependencies of SAR data. This complete theoretical background allows for a systematic analysis of all the parameters influencing the SAR data formation. Therefore, the proposed approach provides a way to separate the most relevant contribution (i.e., the range slopes) from the others, thus somehow circumventing the ill-posedness of the problem.
The proposed technique candidates itself for innovative applications, beyond topography and toward SAR image analysis, elaboration, and interpretation. In this direction, starting from the range slopes map, the authors developed a new intriguing perspective in the frame of SAR despeckling. In this field, it has been demonstrated that a recently introduced scattering-based despeckling approach can be supported by the proposed approach: in particular, the a priori information about surface topography required by the filter can be easily estimated from the SAR image itself, thus avoiding the need of additional information and data.

Conclusions
In this paper a novel model-based inversion technique for retrieving the range slope map from a single SAR image of natural surfaces has been developed. The main innovations are both conceptual and applicative: the use of fractal geometry to model both natural surfaces and scattering phenomena represents the core of the proposed theoretical framework, filling a general gap in the state of the art, where most of the works use Lambertian-like models, inadequate to properly model the complex dependencies of SAR data. This complete theoretical background allows for a systematic analysis of all the parameters influencing the SAR data formation. Therefore, the proposed approach provides a way to separate the most relevant contribution (i.e., the range slopes) from the others, thus somehow circumventing the ill-posedness of the problem.
The proposed technique candidates itself for innovative applications, beyond topography and toward SAR image analysis, elaboration, and interpretation. In this direction, starting from the range slopes map, the authors developed a new intriguing perspective in the frame of SAR despeckling. In this field, it has been demonstrated that a recently introduced scattering-based despeckling approach can be supported by the proposed approach: in particular, the a priori information about surface topography required by the filter can be easily estimated from the SAR image itself, thus avoiding the need of additional information and data. The proposed theoretical framework has a direct and fascinating impact also in SfS. As shown in the paper, only with a proper modeling step, there is a chance to reconsider the position of SfS principles within the information retrieving methodologies from the SAR data. Indeed, it has been demonstrated that a proper modeling of the SfS problem ensures reliable results even by means of an extremely simple inversion procedure. In particular, in this paper, we introduced a linearization of the direct model in order to provide a straightforward evaluation of the local range slopes from intensity data. A MMSE regularization procedure has been developed to recover azimuth slopes and the elevation map of the surface. Both the direct model and the SfS inversion method have been conceptually assessed, analytically derived and an extensive series of experiments have been set up to evaluate the performance of the proposed inversion algorithm using both simulated and actual SAR images.
The main strengths of the proposed SfS technique are (1) the extremely low computational cost if compared with other information retrieval methods and (2) the small amount of a priori information required; indeed, using typical values for T and ε r , no a priori knowledge is needed to recover the relative surface elevation. However, the forward model developed allows to easily introduce a priori knowledge about the electromagnetic and/or fractal parameters of the surface, thus improving the overall performance of the SfS algorithm. The proposed procedure for the estimation of the calibration constant still requires the interpretation of SAR data and then still limits their use to expert users.
In the future, the introduction of the complete direct model in the inversion procedure will ensure a general performance improvement on range slope estimation and, in turn, on all the presented value-added products. This was envisaged in [57] with regard to the local incidence angle map. In this way, the small-slope regime hypothesis could be removed, thus widening the range of applicability of the proposed method. For what concerns DEM and incidence angle estimation, better performance will be surely experienced if a more complex iterative retrieval procedure suitable for fractal surfaces is set up, and, more in general, if a convenient set of procedures and algorithms will be developed to support its principles (like already done in interferometry). A further improvement could derive from a more appropriate integration step and a more suitable global azimuth regularization procedure suitable for fractal surfaces. In this sense, an approach similar to 2-D phase unwrapping methods could be developed.
We also note that the use of lower frequencies can be convenient whenever the assumption of electromagnetic homogeneity of the surface must be enforced, especially in the presence of vegetated areas. Before we apply the proposed SfS technique, a SAR image segmentation for the homogenous segments could be a reasonable way to ensure accurate results over wide areas. Moreover, the necessity of gaining a better insight on volumetric scattering effects (especially at low grazing angles) emerges from the actual case study.
SAR images acquired with high radar-look angle could greatly reduce geometric distortions, primarily layover and shadowing, whose irreversible effects irreparably affect the technique. Experiments have shown that also in small-slope regime foreshortening effects could be evident, impairing surface and incidence angle reconstruction.
Another key issue to deal with is speckle, since it can seriously compromise the reliability of the output provided by the proposed approach. As shown in the last experiment, despeckling greatly improves the estimation of the calibration constant from data, thus improving the overall performance of the algorithm. Therefore, in contrast with other existing SfS techniques, the proposed method does not require a preliminary speckle filtering step [58].
Author Contributions: All the authors equally contributed to the conceiving, writing, and experiments of the work.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.