#### Abstract

A solution to the inverse problem for a three-layer medium with nonsmooth boundaries, representing a large class of natural subsurface structures, is developed in this paper using simulated radar data. The retrieval of the layered medium parameters is accomplished as a sequential nonlinear optimization starting from the top layer and progressively characterizing the layers below. The optimization process is achieved by an iterative technique built around the solution of the forward scattering problem. The forward scattering process is formulated by using the extended boundary condition method (EBCM) and constructing reflection and transmission matrices for each interface. These matrices are then combined into the generalized scattering matrix for the entire system, from which radar scattering coefficients are then computed. To be efficiently utilized in the inverse problem, the forward scattering model is simulated over a wide range of unknowns to obtain a complete set of subspace-based equivalent closed-form models that relate radar backscattering coefficients to the sought-for parameters including dielectric constants of each layer and separation of the layers. The inversion algorithm is implemented as a modified conjugate-gradient-based nonlinear optimization. It is shown that this technique results in accurate retrieval of surface and subsurface parameters, even in the presence of noise.

#### 1. Introduction

Many environmental, industrial, military, geological, and civil engineering applications require fast and accurate noncontact characterization of subsurface structures that have multiple non-smooth layers. These applications include estimating the depth of bodies of fresh water, through-the-wall imaging, detecting buried or sunken objects, and estimating root-zone soil moisture profiles. There have been considerable efforts to develop accurate, versatile, and computationally efficient forward and inverse scattering models for this type of structure for application to practical remote-sensing systems. Penetrating abilities of low-frequency electromagnetic radiation make it a good candidate for this group of applications, provided sufficient resolution can be achieved, and a diversity of measurements (such as multiple polarizations) is available to provide sufficient independent information for retrieval of several unknowns. In particular, multipolarization radar systems at frequencies of L band (~1.2 GHz) and lower are well suited for this purpose. Through the use of accurate and efficient forward and inverse scattering techniques, these radar measurements can be converted into estimates of the sought-after parameters of the surface and subsurface media. The retrieval of unknown parameters is accomplished through a sequential layer characterization, implemented as a nonlinear optimization starting from the top layer and progressively working through the layers below. The optimization process is accomplished by an efficient iterative technique built around the solution of the forward scattering problem. Forward scattering models can usually be classified as analytical, numerical, or methods combining analytical and numerical approaches. One of the most popular analytical methods is the small perturbation method (SPM) developed for a single rough interface by Rice in the early 1950s [1]. In recent years the method was extended to include higher-order expansions [2, 3] and to apply to two or more rough interfaces [4–7]. Since SPM is a purely analytical model, it is very computationally efficient. However, since it is based on truncation of a series, the range of validity of the model is limited to surfaces with small and slowly varying roughness. Fully numerical techniques such as methods based on finite difference time Domain (FDTD) [8] and method of moments (MOMs) [9] can in principle be applied to almost any distribution and roughness scale. When applied to random rough surfaces, many realizations are needed to fully capture surface statistics [10, 11]. In addition, the surface needs to be densely discretized to ensure accuracy of solutions. Therefore, such methods are computationally expensive and are rarely used in inversion algorithms. Methods such as the ones based on Extended Boundary Condition method (EBCM) are significantly more efficient than fully numerical methods, account for higher-order scattering, and have a significantly greater region of validity [10, 12] than the approximate analytical solutions. However, these methods are still significantly slower than analytical methods and are too computationally expensive to be practically used in inversion algorithms.

The method proposed in this work is a transformation of an EBCM-based algorithm that produces simple analytical expression relating parameters of interest (e.g., dielectric constants and layer separation) to backscattering coefficients. Once the model coefficients are obtained, all subsequent calculations are analytical and do not accumulate error. The analytical model is computationally efficient and can be used in both local and global inversion algorithms.

There has been significant progress in development of efficient inverse scattering algorithms which could be applied to real data. In [13], a fast forward solver based on the steepest descent fast multipole method is used in conjunction with “marching-on” optimization strategy. A statistical inversion technique is implemented in [14] to estimate soil moisture and surface roughness parameters from backscattering data. The performance of layer-stripping and mean-square error minimization algorithms is analyzed in [15] as applied to pulse radar data. In [16] a Distorted-Born Iterative Method (DBIM) and Local Shape Function (LSF) algorithms were used as a part of a practical nondestructive testing system.

Inversion algorithms can be generally classified as local or global optimization methods. Local optimization methods such as variations on Newton's method [17] and conjugate gradient methods [18] are iterative methods, which generally rely on gradients and partial derivatives to obtain the next estimate. Since numerical evaluation of derivatives is not accurate and accumulates error, these methods are best suited for situations when derivatives can be computed analytically. Local optimization methods usually converge rapidly (at least as compared to global methods) but may get trapped in a local minimum and fail to converge to the correct solution. Global optimization methods such as simulated annealing method [19] are mathematically guaranteed to converge but may take hundreds of thousands of iteration to do so. With the tremendous increase in computational power in recent years, there have been a number of works which successfully apply global optimization algorithms to solve inverse scattering problems [19, 20]. While it is feasible to use global optimization methods, local optimization techniques are still preferred from the point of view of computational efficiency. The algorithm used in this work is a local optimization algorithm that is based on a Conjugate Gradient method. Techniques that increase the chance of converging to a correct solution and improve efficiency of the algorithm are described later in the paper.

Section 2 describes the forward model, first presenting the full EBCM derivation, followed by closed form representation and sequential layer characterization. The sensitivity of the forward model to errors in ancillary parameters and to subsurface parameters is explored. Section 3 describes the inversion technique, computational efficiency of the inversion algorithm, and effect of observation parameters on inversion results. The inversion results are presented along with error analysis.

#### 2. Forward Scattering Model

The geometry of interest in this paper is a three-layer medium with a periodic interface between the first and second layers and a rough interface between the second and third layers, as shown in Figure 1. This geometry occurs frequently in natural and man-made structures (see [22]). Lakes and rivers often have this cross-section, with the periodic top surface representative of waves and the rough bottom surface representative of the lake or river beds. Plowed agricultural fields are another example of periodic-over-rough cross-section geometries, where there is often a rough interface below the top periodic surface, where transition between sandy and clay soils occurs. The first layer is typically air, that is, a homogeneous dielectric with the dielectric constant of one. The 2nd and 3rd layers are also modeled as homogeneous layers with relative dielectric constants of , , which are in general complex. Figure 1 depicts the relevant geometric parameters for the problem: amplitude and period for the periodic surface, RMS height , and correlation length for the rough surface, layer separation *d*, dielectric constants , , , incidence angle , and polarization.

The forward scattering model is based on cascading reflection and transmission matrices for each interface to form a generalized reflection, or scattering, matrix for the entire problem. Normalized radar cross-sections are then derived from the generalized reflection matrix. The method used to obtain reflection and transmission matrices for each interface is the Extended Boundary Condition Method (EBCM). Reflection and transmission matrices for each interface are given by the following (see [10] for more detail): with R and T standing for reflection and transmission matrices looking down and and looking up. The and matrices are given by the following: where and are given by the following: In (3) is a function that describes the height profile of the interface as follows:

In this method a continuous angular spectrum is represented by a finite number of propagating Floquet modes. A mode is propagating if is real, and evanescent if imaginary. The expressions for and represent integrals over the surface for one period . To properly match the modes from periodic and rough interfaces, the solutions are required to have the same period . In the case of a rough surface it is an artificial period which must be large enough (typically at least 20 ) to obtain accurate statistical representation of the surface. For the periodic surface it must equal an integer multiple of its natural period. The total reflection matrix of the entire system is given by the following: where is a phase shifting matrix and is given by the following: and through are angles of the Floquet modes. Once a total reflection matrix for the system is obtained, radar scattering coefficients can be computed as follows (see [23]): where and [] is column vector containing the incident field coefficients. The model has been validated for several special cases with SPM (see [10] for the two-rough-surface case).

Table 1 compares computational efficiencies of MOM and EBCM methods. The simulations were performed on a machine with a single core 3.0 GHz Intel CPU and a 3 GB of RAM. Special care is needed to address the differences in the type of metrics used to assess computational efficiency of different numerical models. In the case of method of moments (MOMs), unknowns represent induced current elements on a surface of a scattering target and, therefore, have to be densely discretized to preserve the accuracy of the model. Generally, it is recommended to have at least 10 unknowns per wavelength. In the case of EBCM, modes represent discrete scattering directions and, therefore, usually many fewer modes than MOM segments suffice. EBCM becomes even more computationally advantageous as compared to MOM when applied to multilayer problems since scattering matrices are developed independently for each interface and then cascaded using only a few multiplications.

##### 2.1. Creating Closed Form Representations

Even though the computational complexity of the forward model based on EBCM is significantly less than a fully numerical technique such as MOM, it is still too costly for use in most iterative inversion algorithms. This problem is especially serious for global optimization algorithms which often require hundreds of thousands of iterations to converge. Local optimization methods usually converge faster, but they require multiple computations of partial derivatives which quickly accumulates significant errors if done numerically, increasing the likelihood of being trapped in erroneous and/or local minima. To arrive at a model that is more suitable for inversion, we use the above model to derive an equivalent analytic (closed form) representation by precomputing the backscattering coefficients for a comprehensive set of parameters, followed by several function fits. A preliminary version of this equivalent-model approach was developed in [24] for scattering from vegetation canopies. More specifically, the full model is simulated for a range of parameters such as parameters related to dielectric constants (water content), surface height statistics, and layer separation. Then, the dependence on each of these parameters is sequentially modeled using much simpler, analytically differentiable functions, such as polynomials of arbitrary orders. The result is a multidimensional nonlinear, but closed-form (analytic) function. Once the proper closed-form model is developed, the subsequent evaluations of the forward model and its partial derivatives are extremely fast, and the forward model is suited for both local and global inversion techniques. The derivation of this analytic-form model has several intricacies, which are discussed later in this paper.

##### 2.2. Sequential Layer Characterization

The inversion process requires at least as many independent data points (measurements) as there are unknowns. The initial overhead cost of simulating the function for a range of values of all of the unknowns may be too large to be practical, especially considering that the lower interface is a random rough surface, requiring the scattering simulation of many realizations of the surface. Therefore, a sequential layer characterization algorithm is applied by using a multifrequency radar scenario. Even if the inverse problem is solved simultaneously for the unknowns of all layers that include complex permittivities, roughness parameters, and layer separation, at least two frequencies are required to obtain accurate results (see [25]). Assuming the medium between the first and second interfaces is uniform but lossy, the scattering problem is first simulated at a sufficiently high frequency such that the effect of the 2nd interface is negligible. With the assumption that only the top layer affects the cross-section coefficients at this frequency, the number of unknowns is greatly reduced, and the retrieval of the top-layer unknowns can be efficiently accomplished. Since the top layer is a deterministic periodic surface, there is no need for computing multiple realizations of the surface to obtain statistically representative values. With the top layer characterized, the scattering contribution of the subsurface layer can be simulated at a lower frequency. The coupling between the two interfaces is still fully represented in the solutions of the forward and inverse models through the lower-frequency radar data, but the retrieval of their properties has been effectively decoupled through this approach.

##### 2.3. Sensitivity of Forward Model to Errors in Ancillary Parameters

While the forward model contains many parameters, only a few, usually parameters related to dielectric constants (for example, volumetric moisture content for soils) and separation of layers are the parameters of interest. For a successful inversion, all other parameters, which we will call ancillary parameters, still have to be obtained as well. Generally there are three methods of obtaining these parameters. These parameters can be treated as unknowns, with the inversion algorithm producing values for these parameters along with the primary sought-for parameters. While this approach is often the only feasible way to obtain ancillary parameters, it has some major drawbacks. Each additional unknown significantly slows down the initial simulation. Generally, the minimum number of simulations required to obtain a unique solution is given by the following:
is Number of simulations, is Number of unknown parameters, and is Number of values of *m*th parameter to cover the desired range.

Increasing the number of unknowns quickly leads to prohibitively costly initial simulations, as well as to the increased likelihood of the inversion algorithm converging to a local minimum. Furthermore, a large number of unknowns requires a correspondingly large number of measurements. Another method for obtaining ancillary parameters is a-priori approximation, provided that the expected error in such approximation is small and quantifiable. However, careful sensitivity analysis must precede the approximation. This method is generally only feasible if the sensitivity of scattering coefficients to that parameter is low in the region of approximation.

Direct measurement is another alternative for obtaining ancillary parameters. Since the cost of the measurement is usually proportional to the accuracy of the measurement performed, careful sensitivity analysis is needed to estimate the maximum allowable error in the ancillary parameter that would meet the accuracy requirements for the scattering measurements. For a periodic surface the surface height amplitude is usually an ancillary parameter, whereas the dielectric constant is a parameter of interest. Figure 2 shows the sensitivity of the model to errors in the amplitude. The backscattering coefficients are very sensitive to errors in the amplitude, and, therefore, this parameter has to be obtained to a high degree of accuracy. For the rough interface the ancillary parameters are usually surface roughness and correlation length . As can be seen from Figure 3, the backscattering coefficients are quite sensitive to the errors in of the lower surface. However, this sensitivity is not uniform and decreases for increased roughness. The sensitivity to errors in correlation length is relatively small and is nonuniform as well (see Figure 4). For regions where the sensitivity is the smallest, a simple a-priori approximation is enough to obtain acceptable retrieval results for other parameters of interest.

##### 2.4. Sensitivity of the Forward Model to Subsurface Parameters

There is great interest in the Earth Science community to map soil moisture on a global scale. In most cases the air/soil system is modeled as two homogeneous dielectric regions separated by a rough interface. Such models ignore the effects of subsurface interfaces that are often present in reality. Soils, particularly croplands, can often be considered to be periodic interfaces with a small roughness on top and a rough interface between the top layer of soil and the subsurface layer (for example, bedrock). It is, therefore, important to investigate the cases when the effects of the 2nd interface are significant and when they can be neglected. There are several key parameters that are especially significant: depth of the subsurface layer, difference in soil composition, and water content (dielectric contrast) between the layers, roughness (and periodic properties) of top and bottom interfaces, and loss in the middle layer.

Figure 5 shows the effect of the depth of the 2nd interface. In this case, the medium under the first interface is not very lossy (); therefore, there is a clear oscillatory pattern in the error with the envelope of the error that is gradually getting smaller (due to loss). The error can be even more significant if the roughness in the interface between the 2nd and 3rd layers is large (see Figure 7). If the losses are small and the 2nd interface is close to the top surface, the modeling errors can be too large without the inclusion of the subsurface in the scattering models. At higher frequencies, the effects of the subsurface rapidly diminish. This fact makes low frequencies very attractive when characterizing the subsurface properties. Figure 6 shows the error in predictions of backscattering cross-section as a function of layer separation for L band. At higher frequencies the effect of the subsurface is only significant when the 2nd interface is shallow.

#### 3. Inversion Algorithm

Once the numerical simulations over the desired range of parameters are complete, special care is needed to select the proper polynomial representation for the data. While it may be tempting to fit a high-order polynomial such that the residual error is exactly zero, the polynomial function may oscillate rapidly between the sample points, which represents nonphysical behavior and significantly degrades the inversion algorithm. A lower-degree fit, on the other hand, could produce a less accurate function representation with errors accumulating rapidly with each subsequent dimension. To address both of these problems, the method developed here optimally balances the quality of the closed form models and the robustness of the inversion algorithm, as described below. The initial simulated data space is broken up into subspaces of the same dimension but with only a subset of the total points along each dimension (typically 4 points). For example, for the case of 3 unknown variables, the complete data space is a cube, and the subspaces are smaller cubes that completely fill the larger space. Typically, we use polynomials of 3rd order to fit the data in each subspace, resulting in an analytical model for that subspace (Figure 8). The 3rd-order polynomials are fitted with high accuracy and generally produce nonzero 2nd order partial derivatives, which are used in a conjugate-gradient-based inversion algorithm. If a unique solution for the problem exists, it must be contained in one of the data subspaces shown in Figure 8.

The retrieval problem is solved for each individual data subspace, and the best solution is then selected based on the magnitude of an appropriately defined cost function. The inversion algorithm used in this work is a local optimizer based on conjugate gradient method. The cost function to be minimized can be expressed as
where is a solution vector, with corresponding to the dimension of the inversion, is the *n*th component of independent observation vector, and is the *n*th estimate. The procedure starts by computing two initial vectors: and . Then, the following equation is evaluated iteratively until sufficient level of convergence or the maximum number of iterations is reached:
where and are given by

Since all of the operations are performed on closed-form functions, these operations, including multiple gradient computations, are analytic and do not accumulate errors. The algorithm is very efficient compared to most global optimization techniques and usually converges in less than 50 iterations. Since the inversion algorithm is a local optimizer, the convergence to the global minimum is not guaranteed [25]. Breaking up the problem into subspaces increases the chances of finding a global minimum, since the algorithm can only get trapped in the local minimum and not converge to a global minimum if both are present within one data subspace; the probability of this event is very small compared to the case where the search is over the entire data domain.

##### 3.1. Computational Efficiency of Inversion Algorithm

The inversion algorithm described in this work is an efficient local optimizer which is able to achieve high levels of accuracy in relatively few steps. As the data space is broken into subspaces, the problem is solved individually in every subspace. If a unique solution is to be obtained, and if the maximum allowable error is set to a sufficiently low values, the algorithm should only converge in one of the data subspaces and run out of the maximum allowable iterations in others. Since the majority of time is spent looking for the solution in subspaces which do not contain the global minimum, the maximum allowable number of iterations is a critical parameter for computational efficiency of the algorithm. The overall computational time beyond the initial overhead of obtaining the closed form solutions is approximately equal to is total run time to complete the inversion process, is runtime for a single iteration of CG within a subspace, is maximum number of iterations, and is number of subspaces.

The number of iterations it usually takes to reach convergence depends on the type of closed formed functions used to fit the data, the maximum allowable error, and the overall dimension of the problem. Table 2 summarizes the empirically found results. As expected, higher dimensions require significantly more iterations to reach convergence.

The maximum acceptable value for the cost function was set to . The user can set the parameter, the value of which is critical to the computation time of the overall inversion algorithm. If we assume there is a unique solution and there is no subspace overlap, the solution must be contained in only one subspace. Therefore, the algorithm runs times in all other subspaces without reaching convergence. In the case of nonunique solutions, or solutions at the boundaries of subspaces, fewer iterations are typically needed since the algorithm will converge faster than in more than one subspace.

##### 3.2. Effect of Observation Parameters on Inversion Results

In order for the inversion algorithm to produce a unique solution for *N*-independent variables, there needs to be at least *N*-independent observations. However, depending on the characteristics of the forward scattering model, there are still ambiguous cases when two or more combinations of parameters produce the same value for the backscattering coefficients. To illustrate this point using a one-variable case, Figure 9 shows the dependence of backscattering coefficient on the layer separation at two different frequencies. Since it is a one-variable case with the only unknown being the layer separation *d*, one observation could potentially be enough to obtain a unique solution. If the true value of layer separation is , this value corresponds to the cross-section values of −33.7 dB and −36.8 dB for 120 MHz and 150 MHz respectively. If only one observation is used, the inversion algorithm converges to two solutions for each of the frequencies (marked with asterisks), and it is impossible to pick the right solution without additional information. However, when both observations are used, the algorithm inevitably picks out the correct solution. The same principle applies to higher dimensions.

Figure 10 shows the results for one observation (one-variable case). As can be seen, the algorithm often converges to multiple solutions, and more observations are needed to remove the ambiguity. Introducing another observation that is closely spaced in frequency (20 MHz apart) and therefore not completely independent from the original observation helps eliminate a few double solutions, but the algorithm still fails to produce a unique solution every time (see Figure 11). If the second observation is replaced by one taken at a very different frequency (for example, 450 MHz), the algorithm converges to a unique solution every time (see Figure 12). Additional observations would not result in further improvement in the absence of noise, but could be helpful in practical situations when measurement noise is present. The same trend continues for 2 and 3 dimensions (see Figures 13 and 14) although it normally takes more observations for the algorithm to converge to a correct solution every time. When applied to radar remote sensing, these observations are backscattering cross-section measurements at difference frequencies, polarizations, and incidence angles. In the case of observations taken at different frequencies, the frequency separation should be large enough to ensure a necessary degree of linear independence. When observations are cross-section measurements taken at frequencies that are spaced close to one another, the algorithm typically produces multiple solutions.

##### 3.3. Inversion Results and Sensitivity to Errors in the Backscattering Coefficients

The algorithm described in this work is an efficient and robust method for retrieving parameters of interest from radar cross-section measurements. Once the initial simulation over multiple parameters is complete and subspace-based closed formed expressions are obtained, all subsequent operations such as computing 2nd order partial derivatives are analytical and therefore do not accumulate errors. In this section we show validation results, first using noise-free simulated data for a wide variety of layered medium parameters, followed by the results for the noisy case. Figures 15 and 16 show estimated parameters versus true values for a 3-variable inversion (the other variable being the roughness of the layer). In this case six observations at different frequencies were used (120 MHz, 150 MHz, 180 MHz, 400 MHz, 435 MHz, and 460 MHz).

Since all measurements include some error, it is important to analyze the robustness of the inversion algorithm when real observations containing errors are used. Figure 17 shows the same scenario as Figure 15 for layer separation retrieval, but this time 0.1 dB of Gaussian noise was added to the observations. Increasing the error to 0.5 dB (see Figure 18), we can see the increase in the error in the retrieved parameters, but, overall, the errors remain relatively small, and the algorithm does not break down. Increasing the error to 1 dB (see Figure 19), we notice the further increase in error, but the algorithms still perform considerably well. Similar analysis is then performed to analyze the error performance of the algorithm when retrieving the water content of the 3rd layer. Typically, the sensitivity of the backscattering coefficients to the moisture content in the 3rd layer is lower than to the layer separation, and therefore it is important to make sure that the inversion algorithm is robust enough to be effective. Figures 20, 21, and 22 demonstrate that the algorithm retains robustness in retrieval of the water content in the 3rd layer. Since the parameters retrieved from the top interface are incorporated in developing closed-form equations relating subsurface parameters to backscattering coefficients, it is especially critical that the retrieval of these parameters can be as accurate and robust as possible. To minimize the contribution of the subsurface, the simulation was performed at three L-band frequencies: 1000 MHz, 1100 MHz, and 1200 MHz. Figures 23 and 24 depict the retrieval of water fraction of the middle layer and amplitude of the periodic interface for the errorless case. With the data from the three measurement channels, the inversion results are very accurate and converge to a correct solution in every case. Figures 25, 26, and 27 depict the inversion results from the amplitude of the periodic interface when .1 dB, .5 dB, and 1 dB RMS error is introduced to the backscattering coefficients. As predicted the error in the retrieved amplitude increases with increased error in the backscattering coefficients but the algorithm retains robustness.

#### 4. Conclusion

An efficient and robust algorithm for retrieval of parameters of interest from scattering coefficients is described in this work. The forward model uses EBCM to obtain a reflection matrix for each interface. The Floquet modes for each interface are matched to relate parameters of interest to scattering coefficients for the entire system. The full model is simulated for a range of parameters, and the resulting data space is broken up into multiple subspaces, and closed-formed representations are determined for each of the subspaces. The conjugate gradient based inversion algorithm is applied to produce a solution in each individual subspace, and the best solution is then picked out based on the magnitude of an appropriately defined cost function. The sensitivity of the algorithm to errors in ancillary parameters and scattering coefficients is carefully analyzed. Experimental validation of the proposed models is an ongoing work, with some earlier tower-based radar system producing inconclusive results due to calibration and hardware problems (see [26]). Currently, a newly developed triband radar system is being tested showing promising preliminary results. With the expected availability of radar instruments—tower- or aircraft-based—operating simultaneously (or in tandem) at L and P bands, the technique described in this paper provides an effective tool for high-accuracy and robust retrievals for a variety of subsurface sensing applications. This technique can find applications in missions like the upcoming NASA Earth Ventures-1 Airborne Microwave Observatory of Subcanopy and Subsurface (AirMOSS), with the possibility of flying AirMOSS and the NASA/JPL Unmanned Aerial Vehicle Synthetic Aperture Radar (UAVSAR) (P-band and L-band) pods on the same aircraft.