Abstract
L-band radars have been proposed as possible way for monitoring landslides. In this paper, we examine and solve the principal difficulties arising in modeling and processing radar data, evidencing differences with more usual SAR imaging. Numerical examples in support of the proposed processing procedure are finally provided.
1. Introduction
Amongst the different possible ways to have a continuous monitoring of areas which may be subject to landslides, L-band radars offer the possibility to penetrate foliage while still being able to get some understanding of the evolution of the scenario by means of differential imaging techniques.
Differently from a large body of the literature, where imaging is performed by means of satellite-based radars (so that the movement of the satellite allows to rely on a synthetic aperture) the research activity considered in the following concern the exploitation of a fixed ground based radar, where the only eventually available movement of the sensor is achieved by a mechanical or electronic scanning of the antenna pattern.
Such a circumstance implies a number of interesting differences, which are discussed in the following, with respect to more usual radar imaging modeling and processing techniques.
In fact, a more detailed and difficult model is required for data simulation (which is useful to “tune” imaging procedure). Moreover, data processing requires more sophisticated techniques with respect to satellite-based imaging. On the other side, location on ground of the sensor allows a very simple deployment on those areas which are judged to have a risk of landslides. Also, the assumptions on “coherence” (see below), which are needed for differential imaging, are more easily verified with respect to differential interferometric SAR techniques.
For the sake of simplicity of explanation, most concepts regarding simulation and processing are explained with reference to a simple “2D geometry,” that is, in a case where both fields and the scenario are invariant along one dimension, and the field is directed along such a direction. Such a simplifying assumption is then removed in Section 4.
In the following, Section 2 is concerned with the problem of accurately simulating scenarios of interest, and attention is paid to the need of developing models which take into account the fact that the radar antenna is supposed to be in the near zone of the scenario under test. Then, Section 3 presents the basic idea for monitoring possible deformations of the soil on the basis of a differential imaging technique. In particular, the need of solving such kind of inverse source problem is underlined, and a suitable processing technique is introduced and discussed.
Finally, the extension to the 3D case is discussed in Section 4, and results of processing (on simulated data) are presented in Section 5.
Conclusions and possible developments follow.
2. Modeling Radar Data
In order to provide whatever form of imaging technique, it is mandatory to get a deep understanding of the physical mechanism underlying the process and to get a proper accurate mathematical modeling of this latter. In fact, this will allow a better comprehension of the phenomena at hand (possibly suggesting proper imaging techniques), and it will allow an accurate simulation of radar data, which is of course essential to properly tune imaging techniques.
The basic physical mechanism is the electromagnetic scattering [1]. In general case, backscattering is due to both a kind of reflection at the interface, as well as from the secondary field which is generated within the volume of interest.
By taking into account the values of permittivity and conductivity at the frequencies of our actual interest, the predominant mechanism is “surface scattering.” In fact, as it can be deduced from Figure 1, the penetration depth, concerning the four kinds of soils described in [2], at the frequency range of interest is expected to be too small to be able to consider the mechanism of “volume scattering.”
As no exact closed form, but only approximate solutions exist for the problem of scattering from a irregular surface [1], to derive a mathematical approximate expression of the scattered field, it is assumed a surface with gentle undulations [1], whose average horizontal dimension is large compared with the wavelength. This is a general case because if one wants to consider a rough surface, it is only necessary to add a term of roughness, according to “Clapps models” [3].
It is also assumed that the total field at any point on the surface can be computed as if the incident wave is impinging upon a infinite plane tangent to the point, according to the basic assumption of the Kirch off method [4–6] or “facets model” [1]. A mathematical statement of the scattered field, formulated by Stratton and Chu [4, 5] and modified by Silver [6], is as follows: where(i) a time factor of the form is understood;(ii) ;(iii) is the unit vector in the scattered direction;(iv) is the unit vector normal to interface;(v) is the vector that scans the points on the surface;(vi) and are, respectively, the intrinsic impedance and wave number of the medium in which the scattered field is evaluated;(vii) is the range from the center of the illuminated area to the point of observation;(viii) and are the total electric and magnetic fields on the interface.
Equation (1) allows computeing the scattered field in far zone [7] as superposition of fields scattered by each point on the surface. On the other side, in the case of our interest the receiving antenna is in the near zone of the scenario at hand (which acts as a source), so that (1) is noted as adequate to accurately modeling the radar data.
So, (1) has to be reformulated considering that some parameters, like the incident angle or the distance between the sensor and the scenario, change for each point on the surface. Accordingly, it is not possible to use the usual simplification for far-field region, that is, to approximate the distance factor in the denominator by and exploit a first-order approximation for the argument of the exponential.
Another basic difference with respect to SAR techniques concerns the kind of data which is collected.
Obviously, one needs some kind of data diversity to get an image after processing, and if a 2D image is the final goal, the raw data must be a function of two coordinates at least, which is commonly achieved in SAR systems by a translation of the platform and registering delays of the backscattered signals. In our scenario, the “platform” is fixed, and we would like to exploit monochromatic signals, which can allow a very simple data processing (see Section 3).
Hence, the needed diversity is achieved by a suitable 2D mechanical or electronic scanning of the antenna, so that our data will be a 2D function of , where such a couple denotes the direction where the radar antenna is pointing to.
By assuming a scenario (and incident fields) which is invariant along one direction, and assuming an incident field parallel to such a direction, the overall model becomes scalar, and the field collected at the radar position can be written as where where(i) is the observation angle of the measuring sensor; (ii) is the wave number of the medium;(iii) is the Fresnel reflection coefficient;(iv) is the array factor of the sensor used for monitoring,(v) is the unit vector that identifies the incident direction,(vi) is the observation angle that scans the elementary areas on the surface (direction where the radar antenna is pointing to).
In (2), one can identify a factor depending on the antenna, a factor depending on the round-trip, and a factor depending on the scenario (shape of the interface and its electromagnetic characteristics).
A numerical simulator based on (2) has been developed successfully and tested in some canonical scenarios (flat and sinusoidally varying surfaces).
3. 2D Case Data Processing
3.1. The Basic Idea
In order to describe the strategy, it is convenient to rewrite the relation between the scattered electromagnetic field and the surface features, already described in (2), in a different form. In particular, the scattered electromagnetic field at the radar position when the radar antenna is pointing toward can be expressed by the following equation: where:(i) is the contribution to the overall received field due to each elementary area on the surface;(ii) is the generic distance between antenna and surface, in particular , where is a “reference distance” and is the amplitude corresponding to the surface deformation; (iii) is the array factor (in transmission and receiving case);(iv) is a function, called , that depends on surface characteristics, that is, local reflection coefficient, roughness, and local height of the surface;(v) is a term that depends on the round-trip path of the electromagnetic field.
Figure 2 shows the graphic representation of the reference scenario.
So, (4) can be rewritten as follows:
The processing used herein in order to locate possible surface deformations is based on differential interferometry [8], which works as follows. The phase of the reflectivity function contains information about the signal round-trip path. Therefore, a surface deformation causes a different length of round-trip path and consequently a phase displacement of the signal. Then, if and are reflectivity functions to two times and , one can write: where and are the maxim intensities of the two signals, the complex signals and take into account possible variations on the factor , and the exponential terms contain information about surface displacements. Then, “beating” the two signals, one achieves
Then, as long as the phase of has not changed amongst the different instant times (e.g., under a “high correlation” assumptions), the phase of expression (7) gives back which is obviously related to displacements. However, a last difficulty has to be tackled. In fact, the phase factor has an ambiguity of 2 (wrapped phase). So, a (possibly effective) phase unwrapping procedure is needed.
3.2. Extracting the Reflectivity Function
Traditional synthetic aperture radar (SAR) processing procedures heavily rely on the convolutional nature of the equation to be inverted [9]. In fact, such a circumstance is exploited to develop computationally effective codes based on fast Fourier transforms.
In our case, due to the expression of and (a different choice would have been possibly by incorporating in the reflectivity, but then a much more cumbersome unknown has to be retrieved), one cannot rely any more on convolutions, so that different processing strategies have to be devised.
As in any integral equation with a smooth kernel, problem (5) is ill posed. (According to Hadamard definition, a problem is ill posed when its solution is not unique, or does not exist or does not depend continuously on the data.)
As a consequence, only an approximate version of the true reflectivity function can be found by looking for some kind of regularized solution.
Amongst the different possibilities, the well-known truncated singular-value decomposition (TSVD) [10] can be used. As a suitable alternative, a Tikhonov regularization [11] can be exploited, which is the suggested solution in case one has to deal with a very large number of unknowns.
3.3. The Phase Unwrapping Problem
As previously described, the phase factor in (8) has an ambiguity of 2 (wrapped phase) [12, 13]. The mathematical statement of the phase unwrapping problem can be expressed by and where the integer unknown has to be determined for each point of the scenario.
Two-dimensional phase unwrapping is a classical problem encountered in several fields such as electromagnetic theory, optics and DInSAR data processing. For this reason, there is a large amount of research [14, 15] and solution approaches to both 1D and 2D phase unwrapping problem. Roughly speaking, solution approaches span from model-based procedures [14] to “total least squares” [15, 16]. Compressive sensing-based procedures have also been recently proposed [17].
4. Extension to the 3D Case
The approach previously described for the 2D geometry can be extended to the 3D case. In particular, it is possible to rewrite an equivalent approach in terms of open-circuit voltage [7] by where(i) is the open-circuit voltage induced by ;(ii) is the array factor in transmission case;(iii) is the array factor in receiving case;(iv) corresponds to the polar coordinates of the points on the surface;(v) corresponds to the polar coordinates of the observation point of the measuring sensor;(vi) is a superscript that identifies the polarization channel.
In fact, in 3D case, one can consider four different types of transmission-reception setup [18]. Notably, for three of these four different “polarization channels,” a different reflectivity function could be defined, but the signal processing required to get the unknown deformation would be exactly the same. Far from being a problem, such a circumstance could be exploited in order to improve accuracy in the deformation reconstruction and/or avoid possible ambiguities in the phase unwrapping step. On the other side, such an approach would need antennas able to work in both the two independent polarizations.
5. 2D Case Data Processing: Numerical Processor and Some Examples
Using contents as above, numerical simulators (for the 1D and 2D cases) have been developed first. Then, the different steps of the processors have also been developed in a MATLAB environment [19].
After validating the numerical simulator by examining results furnished by the codes in simple cases (flat surfaces, a single small deformation, and a sinusoidally varying, deformation), the processing chain has been tested.
The first example is aimed at displaying the reconstruction of a simple deformation, that is, a rigid shift. The working frequency is 2 GHz. The reference surface is the undulating surface shown in Figure 3. Assuming that at the time this surface is rigidly shifted of , the reconstruction is shown in Figure 4. The latter is obtained by the multiplication in (7). The unwrapping procedure used herein is the standard procedure of the MATLAB library. In correspondence to the nulls of the amplitudes of the functions and , where there is a loss of information about the phase, a preprocessing, based on technique described in [15], has also been used. As it can be seen in Figure 4, the surface deformation is correctly indentified.
A second example considers a surface deformation forming by the overlapping between the undulating surface in Figure 3 and the deformation shown in Figure 5. The working frequency is 2 GHz.
The corresponding reconstruction is shown in Figure 6. As it can be seen also in this case, the deformation is correctly identified.
As it can be seen, one is able to retrieve deformations as small as a quarter of a wavelength. In all tests, the simulated data used for processing was affected by a simulated error with SNR = 20 dB.
Similar kinds of results are achieved in the preliminary test we are performing on the corresponding 3D processor.
6. Conclusions
In this paper, approaches commonly used in SAR for imaging of deformations have been generalized to the case of ground-based, nonsynthetic radars operating in the L band.
Such a case is indeed of interest in the applications, as it allows a very simple deployment of low-cost radar systems [20], and it allows to monitor regions which could not be easily accessible from satellites or airborne-based instrumentations. Moreover, the use of relatively low frequencies allows penetration of foliage. Last but not least, different from well-known DInSAR techniques, which requires multiple passages in slightly different orbits, the sensor is herein fixed, which allows a much more easy fulfillment of the coherence property recalled after (7).
On the other side, proper modeling and processing at these frequencies has required the development of both new simulators, as well as the development of the new processing techniques which have been introduced, discussed, and tested above.
Conflict of Interests
All authors state that the research herein described is not influenced by secondary interests (such as financial gain) and that they have no conflict of interests concerning all terms used in this paper.
Acknowledgment
This work has been carried out under the framework of PON 01_01503 National Italian Project “Landslides Early Warning,” financed by the Italian Ministry of University and Research.