A subspace pursuit method to infer refractivity in the marine atmospheric boundary layer

A subspace pursuit method to infer refractivity in the marine atmospheric boundary layer

Marc Aurèle Gilles, Christopher Earls and David Bindel This work was supported by the Office for Naval Research through grant N00014-16-1-2077M. A. Gilles is with the Center for Applied Mathematics, Cornell University, Ithaca, NY, 14850 USA e-mail: mtg79@cornell.eduC. Earls is with the School of Civil and Environmental Engineering and the Center for Applied Mathematics, Cornell University, Ithaca, NY, 14850 USAD. Bindel is with the Department of Computer Science and the Center for Applied Mathematics, Cornell University, Ithaca, NY, 14850 USAManuscript received September 27, 2018; revised January 11, 2019.

Inferring electromagnetic propagation characteristics within the marine atmospheric boundary layer (MABL) from data in real time is crucial for modern maritime navigation and communications. The propagation of electromagnetic waves is well modeled by a partial differential equation (PDE): a Helmholtz equation. A natural way to solve the MABL characterization inverse problem is to minimize what is observed and what is predicted by the PDE. However, this optimization is difficult because it has many local minima. We propose an alternative solution that relies on the properties of the PDE but does not involve solving the full forward model. Ducted environments result in an EM field which can be decomposed into a few propagating, trapped modes. These modes are a subset of the solutions to a Sturm-Liouville eigenvalue problem. We design a new objective function that measures the distance from the observations to a subspace spanned by these eigenvectors. The resulting optimization problem is much easier than the one that arises in the standard approach, and we show how to solve the associated nonlinear eigenvalue problem efficiently, leading to a real-time method.

Ducting, Electromagnetic propagation, Eigenvalues and Eigenfunctions, Inverse Problems, Partial differential equations, Propagation Charaterization, Radar Remote Sensing.

I Introduction

The marine atmospheric boundary layer (MABL) is the part of the lower troposphere in direct contact with the ocean. This contact creates a zone of particularly high inhomogeneity due to the exchanges of heat, moisture, and momentum between the atmosphere and the ocean [1]. Within the lower MABL, the index of refraction - the speed of light in the medium relative to that in a vacuum - may change rapidly with height above the ocean surface; this causes ducting, i.e., bending of EM waves to the surface. Atmospheric ducting greatly changes the behavior of EM propagation within the MABL from what is expected in a “standard” atmosphere [2]. Ducting impacts radio communication, and it is also detrimental to maritime radars: it creates radar holes where no EM wave can travel, increases sea surface clutter, and changes the maximal operating range (illustrated in Fig. 1). It is therefore of great interest to be able to identify and characterize the presence of ducts in real time.

Fig. 1: Effects of atmospheric ducting on EM waves [3].

A variety of methods have been developed to characterize EM ducts. A few methods link refractive index profiles and weather conditions [4] to estimate the MABL characteristics. Such estimates of meteorological conditions may come from numerical weather predictions [5, 6]. Although these numerical weather predictions are successful for long-term, general characterization, their accuracy is unsatisfactory for characterizing ducting in a local sense and in real time. Another way to estimate meteorological conditions conducive to ducting is by using radiosondes or rocket sondes [7] but those methods are costly, slow to deploy and very local. Yet another way to estimate ducting conditions is by lidar [8] but this approach is sensitive to clouds, fog, and aerosols.

Other methods that rely on global positioning system (GPS) satellites have been proposed which use information about the distortion of the GPS signal to estimate refractivity profiles [9, 10]. These methods rely on the GPS being placed over the horizon with respect to the receiver, which renders this method impractical for other contexts than targets of opportunity.

In the last decade, refractivity from clutter (RFC) methods have received a lot of attention in the literature. RFC methods use a radar to estimate the refractivity profile by emitting radiation and measuring the backscattered signal from the rough ocean surface (also called clutter). RFC methods typically use either a forward model of EM propagation or a database to “predict” the clutter under some EM condition, and compare the measurements to the prediction to infer the refractivity profiles. For a review of RFC methods, see [11].

The most popular forward model in RFC applications is the parabolic equation (PE) specialization of Maxwell’s equations, briefly discussed in section II. A wealth of inverse solution methods for characterizing MABL refractivity have been developed that combine the accurate and relatively fast solvers inherited by the PE with a statistical or machine-learning method. For example, [12] uses a recursive Bayesian approach, [13] uses Kalman filters, [14] uses support vector machine, and [15] uses Bayesian Monte Carlo analysis. An example of an RFC method that learns from a database is [16]. It uses the proper orthogonal bases of collected data to form an approximate forward model to be used in an inversion aimed at characterizing the EM duct itself.

In the current paper, we are interested in a different sampling approach: the bistatic case [17, 18, 19, 20, 21, 22, 23, 24]. An example situation could involve two separate phased arrays, one transmitting and one receiving down range, with the receiver able to sample at different heights. Our method is similar to RFC methods in that it uses a forward model to predict clutter, as it also relies on a forward model of EM propagation. However, the proposed method does not involve actually solving the associated differential equation to predict the signal. Instead, we exploit a structure that is present in the partial differential equation which governs the physics: namely the approximately low-rank structure of the field within specific parts of the domain. The low-rank structure arises because only a few eigenvectors are needed to reconstruct the PDE when it is solved through separation of variables. This allows us to design an algorithm that seeks a refractivity profile associated with eigenvectors that best fit the data. The method presented in this paper is close in spirit to the idea presented in [3] and [16]. However, while they used a basis induced by the data, we use a basis induced by the forward model.

The rest of the paper is organized as follows: in section II we state some background on the problem; in section III we motivate and describe our algorithm; in section IV we give details needed for a fast implementation; and in section V we present our numerical results.

Ii Background

Ii-a Forward problem: propagation

The physics that govern electromagnetic wave propagation are described by Maxwell’s equations. Assuming a horizontal polarization and suppressing a time dependence of the form , Maxwell’s equation can be transformed into the Helmholtz equation (cf. (1), along with useful boundary conditions (2), (3), (4)), by means of an exact earth flattening transformation [25] :


These equations are in 2D cartesian coordinates, where denotes the horizontal range, denotes the vertical altitude (the direction of invariance), , is the free-space wavenumber, is the wavelength, is the complex dielectric constant at the ocean free surface, is the radius of the earth, is the index of refraction, and denotes the electric field in horizontal polarization. This Helmholtz equation is equipped with boundary conditions. In the case of the MABL, this is achieved by imposing continuity of the tangential field components by modeling the sea surface as a locally homogeneous dielectric and specifying a surface boundary condition [25]. This surface boundary condition is implemented via the Leontovich surface impedance condition, which for horizontal polarization is expressed as (2). Equation (3) is the boundary at and represents the source, i.e. the transmitter antenna. The domain is semi-infinite in both the and direction, for these boundaries radiating boundary conditions of the form of (4) are appropriate [25].

In the particular case of the MABL, the index of refraction is often approximated to be horizontally constant [15, 14]. This assumption seems to be approximately valid for open-sea for a small region (less than  km) [26] but may not hold within coastal regions. In the argument that follows, we will make this assumption, and thus fix . In section V, we relax this assumption. That is, we allow for some horizontal change in the refractivity and attempt to characterize the mean refractive index, where the mean is taken over the downrange distance.

Ii-B Index of refraction

The literature suggests that in the MABL the refraction can be well approximated by employing a modified refractivity , defined as


where  km is the radius of the earth. is modeled as a tri-linear function represented by four coefficients [27]: is the height of the base of the duct in meters, is the thickness of the duct in meters, is the M-deficit in M-units, and is the slope of the lowest linear portion in M-unit/meter. This parameterization is typically used to represent surface-base ducts [18] for which the height of the ducts is a few tens of meters. Figure 2 displays graphically the parametrization of , and (6) shows the analytical form. The slope M-unit/m of the upper part of the refractivity profile is consistent with the mean over the whole of the United States and is not a sensitive parameter in the inversion [18].

Fig. 2: Example of a modified refractivity profile of a surface based duct. The modified refractivity profile is modeled as a tri-linear function represented by four coefficients: is the height of the base of the duct in meters, is the thickness of the duct in meters, is the M-deficit in M-units, and is the slope of the lowest linear portion in M-unit/meter.

For the rest of the paper, we will refer to the parametrization of , through as . We note that , the modified refractivity at the mean free surface of the ocean, could be included in the parametrization, but propagation of EM waves have been found to be insensitive to this parameter [12]; therefore following the example of the authors in [12], we fix it to a typical value of . We note that the method described in section III does not rely on this parametrization, and any other parametrization could be used.

Ii-C Ssfpe

The most used method for solving the PDE in (1) together with boundary conditions of the form of (2), (3), (4) is to use the split step Fourier transform for the parabolic equation method (SSFPE) [25, 28, 29, 2]. This method relies on a parabolic equation approximation of the Helmholtz equation. It is accurate, stable, relatively fast, and fairly easy to implement. Our surrogate field data is obtained by this method. In particular, our code is based on the software PETOOL, described in [29].

Ii-D Modal solution

For simple boundary conditions and geometries, the Helmholtz equation can be solved exactly by separation of variables [30]; that is,

where the eigenpair 111We use superscripts to emphasize quantities that depend on the refractivity profile parametrized by . are solutions to a Sturm-Liouville (SL) eigenvalue problem:

together with the associated boundary conditions. The functions are called eigenfunctions or modes, and the scalars are the associated eigenvalues ( are also called associated wavenumbers). Throughout, we assume that are normalized so that .

For example, in the case where the source is modeled through a boundary condition at (as a point source at height ) and the boundary condition is homogeneous Dirichlet at along with homogeneous Neumann at , the electric field solution is [30]:

When we consider an infinite domain and more complicated boundary conditions, such as (2) and (4), separation of variables does not provide an exact solution. In this case, we use contour integration to obtain a solution involving a linear combination of modes from the discrete part of the spectrum and an integral term from the continuous spectrum. In practice, the integral term can be neglected if we are sufficiently far from the source [30]. The modes can be further divided into two categories:

  1. Leaky modes which are not observed in range.

  2. Trapped modes which propagate in range.

As noted in [30] for most long-range propagation, only the trapped modes whose wavenumber is within a certain interval of interest are important. In our case, the interval of interest contains the admissible speeds of propagation of the modes. These admissible speeds of the propagating modes are bounded by the minimum and maximum speed induced by the refractivity in the domain. In the ducting case, we are concerned primarily with the energy emitted, propagated, and received in the MABL. Therefore, by restricting our domain of dependence to the MABL, we get heuristics bounds on the set of relevant eigenvalues. Formally: we say that a solution of the SL eigenvalue problem is one of propagating modes if and

where the maximum and minimum are taken over a domain of dependence: . In our case, we define to be the maximal height at which we consider non-standard refractivity. In our numerical experiments in section V, we take m.

Formally, we have for large (on the order of km) and :


where is a propagating mode and . In the case where we use the boundary conditions in (2) & (4) along with a trilinear refractivity profile parametrized with , the first five propagating modes are shown in Fig. 3. Fig. 4 shows the approximation of the field with different numbers of modes used, and the field obtained by using the SSFPE solution. Figure 5 shows the norm difference of and at km and for . We observe that after 4 modes, the field is well reconstructed.

In the rest of the paper, we will assume that the observations are collected at a fixed range, and at multiple fixed altitudes: that is we fix and . We denote


where is the parametrization of .

(a) mode 1
(b) mode 2
(c) mode 3
(d) mode 4
(e) mode 5
Fig. 3: Plots of the first 5 propagating vertical modes induced by a particular refraction index parametrized by .
(a) 1 mode used
(b) 2 modes used
(c) 3 modes used
(d) 4 modes used
(e) 5 modes used
(f) SSFPE solution
Fig. 4: Low order approximation of the field: for increasing numbers of retained modes ((a) - (e)) and also the true field computed by SSFPE (f)
Fig. 5: Plot of error: as a function of the number of modes retained , where km .

Iii Inverse problem: characterizing refractivity

Our problem can be described in the following way: given observations of the EM response within the MABL in the form of observed data at a fixed range, , and different heights given by the vector , identify the prevailing vertical profile of the index of refraction . The BVP described in section II provides a forward model, which given an index of refraction, , lets us compute the field at any points in the domain. We would like to solve the inverse problem: find the refractivity profile given some observations of the field. A natural approach is to seek


where is a parametrization of the refractivity profile, and is the predicted observation under the forward model, e.g., as computed by an SSFP solver.

Iii-a Analysis of the inverse problem properties

The objective function in (9) has thousands of local minima, which makes global minimization very difficult. For demonstration purposes, Fig. 6 shows a plot of a two-dimensional cut of the objective function (9), with for a fixed refractivity profile parameterization .

Fig. 6: 2-dimensional slice of objective function 9. The black dotted line indicates the value of , which by definition of the function is the global minimum.

Although the complex function behavior is a property of the solution, and not of the modal approximation, it is illuminating to reason about this behavior in terms of the modal approximation. Using the modal approximation in (8), we see that in the region of interest:


The matrix spans a basis for a subspace in which the approximation is expressed and the vector represents the coordinates of the approximation within that subspace. Armed with this representation of the forward model, we would like to determine which of the two components causes the highly oscillatory behavior observed in Fig. 6. Fig. 7 shows how the subspace changes as we vary the refractivity profile from a reference value of by plotting the largest principal angle [31] between and . Fig. 8 shows how the coordinate part of the function evolves as we change in the discrete 2-norm, i.e. it is a plot of:

Fig. 7: Principal angle between two subspaces: one constant, and the other induced by different refractivity profiles. The dotted line indicates the value where .
Fig. 8: Difference of norm of “coordinates” induced by different refractivity profiles. The dotted line indicates the value where .
Fig. 9: Objective function of the proposed method plotted in a semi-log scale, defined in equation (13). The dotted line indicates the value of the true refractivity profile.

It is clear from Fig. 7 & 8 that the basis changes smoothly, whereas the coordinate of the function causes the oscillatory behavior of the function in (9). These oscillations are explained by the terms . Indeed, is typically large (here, m); hence any small change in eigenvalue caused by a change in are heavily amplified by the multiplication by , which causes wild oscillations of the term .

To avoid the multimodal behavior of the objective function (9) caused by the wild oscillations in coordinates , we propose an alternative objective:

where (defined precisely in the next section) is a projector onto the low-dimensional space spanned by the propagating modes for the refractivity profile parametrized by . In contrast to (9), the new objective does not oscillate; see (9).

Iii-B Proposed inverse solution method

Our algorithm attempts to find a low dimensional subspace that best explains the measurements. To achieve this goal, we first need a measure of optimality of a subspace. Let , and be a dimensional subspace of . It is natural to define the distance between the vector and the subspace as the minimum of the distance between the vector and any vector , as in (10). A graphical representation of the distance between a subspace and a vector in the case where and is shown in figure 10.

Fig. 10: Distance between a vector and a one dimensional subspace (a line)

Let the columns of span , then the squared distance between and is defined as:


Accordingly, we define the distance between and a collection of vectors , which are the columns of as:


Equation (11) provides a way to measure the fit of some collection of measurements into a subspace spanned by the columns of .

In our case, we have (the EM observations), and (the matrix of propagating modes induced by sampled at the observation points). Therefore, to find the suspace which best fits the observation, we seek to perform the following minimization:


According to our numerical tests, the crucial consideration to make this method successful is to characterize the dimension of the subspace properly. Indeed, it can be observed that the number of modes used, even for a fixed physical setting (i.e. frequency, boundary values, range) depends heavily on the actual refractivity profile that produces the modes. Intuitively, this can be understood in the following way: if the refractivity profile represents a very strong duct, then the trapping increases and therefore more modes are trapped under the duct and propagate in range. In order to automatically choose the number of modes to consider, we propose two strategies:

Iii-B1 Algorithm 1 Minimal subspace dimension

One strategy is to search over increasingly high dimensional subspaces until we find one that fits enough of the data according to some criterion. The specific criterion we use is: the objective function must be lower than some threshold value of , which depends on the noise level222In section V we set where is the noise level.. A drawback of this method is that it requires a good estimate of the noise level. We include a regularization term of into the objective function. We choose the regularization parameter by optimizing over a representative training set. We define a regularized residual function (plotted in Fig. 9):


The algorithm is summarized in Fig. 11.

1:  Input: EM observations , threshold
4:  while   do
5:     minimize
9:  end while
10:  return  
Fig. 11: Algorithm 1: Minimal subspace dimension

Iii-B2 Algorithm 2 Filtering eigenvalues

Another way to determine the number of eigenvectors to include is to use an a priori bound on the eigenvalues to be included. As discussed in section II, we can consider the propagating modes as the eigenvectors associated with eigenvalues which fall within a physically inspired interval and include only such modes. However, using such a hard threshold for included modes within an interval has some drawbacks. One is that at all places of the resulting objective function where a mode is included or excluded, the objective function is discontinuous. The other one is that our characterization of a propagating mode is physically inspired by what we define to be the domain of dependence according to , but the true domain of dependence is infinite. Instead of a strict cut-off, we use a soft thresholding by defining a filter, , in the following way:


and are constants that relate to the a priori bounds of the eigenvalues and is a smooth interpolation. In particular, we set

The modes which have a filtered eigenvalue of 0 are excluded from the objective function. The resulting optimization problem is:


where is the diagonal matrix of eigenvalues . We describe each individual term of the objective function:

  • This term is the same as defined in (12). As discussed earlier, it models how well the data fits into the subspace spanned by the columns of .

  • This is the regularization term of . It forces the choice of coordinate in the basis spanned by to be smaller in the directions associated with eigenvalues that have a small filtered value: that is those that we think a priori are less likely to be propagating. is a regularization parameter is chosen by optimizing over a training set.

  • The norm used in this term is the sum of the absolute values of the matrix entries. This term is the regularization term on the size of the subspace. Indeed, the objective function without this term would be naturally biased towards refractivity profiles that induce a large number of propagating modes. To give intuition as to why this happens, suppose and are the refractivity profiles with propagating modes in common, but that admits an additional st mode that is not propagating for . Then, almost any data that fits will fit even better, even if the extra mode only fits model error of measurement noise. Finally, is a regularization parameter chosen by optimizing over a training set.

It is interesting to note that, unlike an objective function such as the one in (9), the algorithm described in section III-B is agnostic to information about the source or the range of the transmitter. One could imagine that this may be useful if a good model of the source is unavailable, or if the receiver wants to be a passive listener only. The method also combines different observations from different sources or different ranges at virtually no additional computational expense which can increase accuracy.

Iv Implementation

Iv-a A computational shortcut

As described in section II, the Sturm-Liouville eigenvalue problem that arises as a result of the separation of variables in the case we are interested in is posed on an infinite domain:


In our case, we have , , , . One can solve the SL eigenvalue problem with an infinite boundary condition by solving an equivalent problem on a finite domain, but with a boundary condition that is a function of the eigenvalues, for an example of a supporting derivation, see [30]. This new SL eigenvalue problem takes the form:


For the application at hand, is assumed to be linear above and therefore and can be expressed in terms of parabolic cylinder functions. One can then discretize this new SL eigenvalue problem and attempt to solve it numerically. This discretized SL eigenvalue problem gives rise to a non-linear algebraic eigenvalue problem. Nonlinear eigenvalue problems are in general expensive to solve, but since the nonlinearity is only in the boundary condition (in other words, in a single entry of the matrix), it is possible to implement fast solvers. For example, in [32] the author presents an algorithm to solve a similar problem based on Newton’s method. However, in our case, this computational expense is unnecessary as we can deal with a standard linear eigenvalue problem instead. Indeed the modes that we are interested in are the so-called trapped modes. The trapped modes are those which are “trapped” physically low within the domain, and therefore their support lie below a threshold. Note that if satisfies (19) and (20), and , , then also satisfies (21), and therefore solves the non-linear eigenvalue problem. This implies that if we solve the linear eigenvalue problem associated with the SL problem with homogeneous Dirichlet boundary conditions at , and the eigenvector’s derivative is zero at the upper boundary, then this solution also solves the nonlinear eigenvalue problem. Otherwise, we can use this eigenpair as a first guess in a Newton’s method for the non-linear problem as described in [32]. In practice, we find that this step is rarely necessary provided that is chosen high enough, and thus we only solve the Dirichlet problem to do the inference in the numerical experiments in section V in order to save the extra computational expense.

Iv-B Computation of the modes

Implementing the algorithm involves solving the Sturm-Liouville eigenvalue problem numerically. We solve this continuous problem by discretizing using finite differences; a treatment and derivation can be found in [30]. This reduces the problem to one of computing a subset of eigenvectors and eigenvalues of a tridiagonal symmetric matrix, for which optimized solvers can be used. As pointed out in [30], one should use 5 to 10 discretization per wavelength in this type of computation. In our case, this induces a discretized system of size approximately . We have observed the results of the inversion to be insensitive to the discretization size chosen.

Iv-C Inner minimization

Each function evaluation in the algorithm described above involves a minimization over . However, in that inner minimization, is fixed. Thus the inner problem is a linear least squares problem which can be solved in closed form. Furthermore, is a small matrix with a number of rows equal to the number of non-zero filtered eigenvalues (usually fewer than 10) and a number of columns equal to the number of vertical slices of sampled taken (on the order of 30). Therefore the inner minimization may be computed at the cost of a small linear solve.

Iv-D Outer minimization

The most computationally expensive part of the evaluation of the objective function is the computation of a few eigenvectors. However, the computation of first and second derivatives of the objective function is much cheaper computationally as it only involves matrix multiplications and linear solves of small matrices. This fact, coupled with the small number of local minima of the objective function, motivated the use of derivative-based local optimization method. We use MATLAB’s sequential quadratic programming method, described in [33], and perform multistart. The starting points are chosen by Latin hypercube sampling. Five starting points are used for each subspace dimension in Alg. 1, and ten starting points are used in Alg. 2.

V Numerical Experiments

V-a Error measures

Defining an error measure is crucial to perform the optimization over the parameters in Algorithm 1 and and in Algorithm 2 over a training set, as well as to evaluate the performance of our method. The error measure is defined as the relative normalized (RNL2) error between a trial and a true . First we define the normalized error by:


We define the RNL2 by dividing the normalized error by the expected value of the normalized error of two random refractivity profiles coming from a very large representative set.


The expectation is taken over the parameters for which we perform the optimization, and is approximated by averaging trials. For example, a score of signifies that the algorithm performs times better than a random guess.

V-B Experiment 1: Simulated data originating from a trilinear, horizontally constant index of refraction.

We simulate electromagnetic wave propagation by solving the partial differential equation in (1) using the SSFPE method. We set the wavelength m, use a Gaussian antenna pattern as the source, and sample the field at a fixed range km and heights of for . Therefore, each measurement is a vector of length 30, whose entry corresponds to the field at range and altitude . We simulate 5 such measurements for each test case, where different measurements are obtained by varying the height (between  m and  m) and tilt (between and degrees off of horizontal) of the antenna (represented as a Gaussian source):

where denotes the uniform distribution on . This forms a matrix of measurements of dimension . We then contaminate the data with Gaussian white noise of standard deviation . The observations are then normalized to have unit norm in order to keep the different terms of the objective function scaled relative to each other.

For the test cases, we fix M-unit/m[17], which is consistent with the mean over the whole of the United States, and has been observed to have very little variability [34]. In order to produce unbiased test cases, we generate twenty synthetic refractivity profiles by randomly sampling the parameters in the following way:

These parameters are consistent with low altitude surface-based ducts. The realizations are shown in Fig. 12. In terms of physical domain, this means that we are observing data from to m in altitude, and are trying to invert parameters that define non-standard refractivity profiles from to m. We attempt the inverse problem of identifying the refractivity profiles from the observational data. We use the algorithms described in section III-B. For algorithm 1, we set . For algorithm 2, we set , and . These parameters were obtained by minimizing the RNL2 score over a separate training set and were found to be insensitive to the noise level. The result of Alg. 1 and 2 on the twenty refractivity profiles are shown in Fig. 12.

Fig. 12: Results for algorithm 1 and 2 on experiment 1. The solid back lines are the true profile that generated the data: , the red dashed plot are the inverted profiles obtained using algorithm 1, and the dotted blue lines are the inverted profile obtained using algorithm 2. For each case, the RNL2 score is given below the plots. The first number is the RNL2 score of algorithm 1, and the second is the score of algorithm 2.

V-C Experiment 2: Simulated data originating from a trilinear, horizontally varying index of refraction.

We present results on simulated data originating from a horizontally varying index of refraction. This experiment violates the assumption that is horizontally constant, thus one cannot expect in general to find an approximate solution of the form (7) which induces the low-rank decomposition of the solution of (1). It is this low-rank structure that is exploited by Alg. 1 and 2 and therefore we expect that their performance would degrade in this experiment.

Aside from the way that the refractivity profiles are generated, the setup is the same as experiment 1. The horizontal variation in the refractivity profiles is achieved by interpolating two trilinear refractivity profiles. That is, we set to some trilinear refractivity profile, and to some other trilinear refractivity profile. We then pointwise linearly interpolate the two refractivity profiles (that is, not the parametrizations) between , to produce refractivity profiles on for all  km. As a result, the interpolated refractivity profiles are not trilinear. We define the “true” refractivity profile as the average of the refractivity profiles between the transmitter and receiver: where is the range of the receiver. We set the range of the receiver to km. Note that this average will also not be trilinear. We generate twenty synthetic refractivity profiles by randomly sampling the parameters in the following manner:

and the parametrization of by sampling

This allows for a variation in the refractivity profile on the order of of the maximal value of each parameter along the propagation path. We then contaminate the data with Gaussian white noise of standard deviation . The result of Alg 1 and 2 on these twenty synthetic test cases are shown in Fig. 13.

Fig. 13: Results for Alg. 1 and 2 on experiment 2. The solid black lines on the left are , the solid black lines the right are , the solid black lines in the middle are , the red dashed plot are the inverted profiles obtained using Alg. 1, and the dotted blue lines are the inverted profile obtained using Alg. 2. For each case, the RNL2 score is given below the plots. The first number is the RNL2 score of Alg. 1, and the second is the score of Alg. 2.

Vi Discussion

Alg. Dataset Mean error Median error Standard deviation
1 1 0.331 0.296 0.2808
2 1 0.231 0.1982 0.129
1 2 0.441 0.370 0.2917
2 2 0.274 0.200 0.208
TABLE I: Error statistics of Alg. 1 and 2 on the datasets 1 and 2.

We observe that in all but two cases (Alg. 1 in Fig. (k)k and Alg. 2 in Fig. (a)a) the RNL2 scores are significantly below one, indicating that the algorithm performed far better than a random guess. Qualitatively, the height and structure of the inverted ducts and true ducts are in most cases similar. Overall, Alg. 2 performs better than Alg. 1 and does not require an estimate of the noise level, thus Alg. 2 is preferable to Alg. 1. The median RNL2 score indicates that the output of the Alg 2. is times more accurate than a random guess. The addition of horizontal variation does not seem to greatly affect the quality of the inference (especially for Alg. 2, where the median RNL2 score is 0.198 and 0.200 for the horizontally constant and varying case respectively).

The two cases where the algorithms produce a RNL2 score greater than one are attributed to a failure in the minimization of the objective function and could be resolved by using more starting points in the local optimization algorithm at the cost of a higher computational expense. However, we note that the ten starting points used in Alg. 2 seem sufficient in the vast majority of cases to allow the use of local optimization algorithms. Thus, we conclude that the issue of multimodality pointed out in III-A is largely resolved.

The timings shown in Table II were performed on a single Intel i7 core processor operating at 3.60GHz. We note that the timings should be interpreted as lower bounds as the algorithms were implemented in MATLAB, and were not fully optimized.

The most important feature of this method is that it is able to overcome the highly multimodal behavior associated with the physics of EM wave propagation. This allows the use a local optimization method instead of a global optimization method which is typically used in the literature (such as genetic algorithms in, for example, [3, 27, 35, 19, 23]).

As this method is able to cheaply find an estimate of the refractivity profile using local optimization, it could also be used to warm-start a different method which would typically require a global optimization method. That is, in a first step, one could use this method to find a good first guess. Then, in a second step, a local optimization method started at that initial guess could be used to minimize a multimodal, but perhaps more accurate objective function such as the ones in [27, 35, 19, 23]. This should allow for more accurate prediction while still benefiting from the lower computational cost associated with local optimization algorithms.

Algorithm Experiment Avg. run time
1 1 3 min 15 sec
1 2 2 min 59 sec
2 1 5 min 12 sec
2 2 5 min 19 sec
TABLE II: Average running time on experiments 1 and 2

Vii Conclusion

We presented a new method for characterizing the refractivity profile in the MABL which relies on the low-rank structure of the field within parts of the domain, inherited by the governing Helmholtz equation. This low-rank structure allows us to formulate the inverse problem in terms of a minimization of a new objective function. The objective function performs a projection operation of the data onto the subspace spanned by the eigenvectors that form the low-rank approximation of an electromagnetic field induced by a particular refractivity profile. Performing an optimization on this well behaved objective function allows us to accurately solve the inverse problem in around five minutes, allowing for real-time characterization of the MABL.

We conducted two numerical experiments to demonstrate the efficacy of the method. The first experiment was conducted on noisy simulated data computed with horizontally constant refractivity profiles, and the second experiment was performed on noisy simulated data computed with horizontally varying refractivity profiles. We observed that in both setups, the method is able to infer the refractivity profile using local optimization algorithms with few starting points thanks to the small number of local minima of the objective function.


The authors wish to thank Dr. Frank Ryan, of Applied Technology, Inc., for his advice and many helpful suggestions during this research work. The authors gratefully acknowledge ONR Division 331 & Dr. Steve Russell for the financial support of this work through grant N00014-16-1-2077.


  • [1] T. D. Sikora and S. Ufermann, Marine Atmospheric Boundary Layer Cellular Convection and Longitudinal Roll Vortices.   NOAA, 2004.
  • [2] I. Sirkova, “Brief review on PE method application to propagation channel modeling in sea environment,” Open Engineering, vol. 2, no. 1, Jan 2012.
  • [3] V. Fountoulakis and C. Earls, “Inverting for maritime environments using proper orthogonal bases from sparsely sampled electromagnetic propagation data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 12, p. 7166–7176, 2016.
  • [4] B. R. Bean and E. J. Dutton, Radio meteorology.   For sale by the Supt. of Docs., U.S. Govt. Print. Off., 1966.
  • [5] G. Lefurjah, R. Marshall, T. Casey, T. Haack, and D. D. F. Boyer, “Synthesis of mesoscale numerical weather prediction and empirical site-specific radar clutter models,” IET Radar, Sonar & Navigation, vol. 4, no. 6, p. 747, 2010.
  • [6] T. Haack, C. Wang, S. Garrett, A. Glazer, J. Mailhot, and R. Marshall, “Mesoscale modeling of boundary layer refractivity and atmospheric ducting,” Journal of Applied Meteorology and Climatology, vol. 49, no. 12, p. 2437–2457, 2010.
  • [7] J. R. Rowland and S. M. Babin, “Fine-scale measurements of microwave refractivity profiles with helicopter and low-cost rocket probes,” Johns Hopkins APL Tech. Dig, vol. 8, no. 4, pp. 413–417, 1987.
  • [8] A. Willitsford and C. Philbrick, “Lidar description of the evaporative duct in ocean environments,” in Optics & Photonics 2005.   International Society for Optics and Photonics, 2005, pp. 58 850G–58 850G.
  • [9] A. R. Lowry, C. Rocken, S. V. Sokolovskiy, and K. D. Anderson, “Vertical profiling of atmospheric refractivity from ground-based GPS,” Radio Science, vol. 37, no. 3, 2002.
  • [10] F. Xie, S. Syndergaard, E. R. Kursinski, and B. M. Herman, “An approach for retrieving marine boundary layer refractivity from GPS occultation data in the presence of superrefraction,” Journal of Atmospheric and Oceanic Technology, vol. 23, no. 12, pp. 1629–1644, 2006.
  • [11] A. Karimian, C. Yardim, P. Gerstoft, W. S. Hodgkiss, and A. E. Barrios, “Refractivity estimation from sea clutter: An invited review,” Radio Sci. Radio Science, vol. 46, no. 6, 2011.
  • [12] S. Vasudevan, R. H. Anderson, S. Kraut, P. Gerstoft, L. T. Rogers, and J. L. Krolik, “Recursive Bayesian electromagnetic refractivity estimation from radar sea clutter,” Radio Sci. Radio Science, vol. 42, no. 2, 2007.
  • [13] C. Yardim, P. Gerstoft, and W. S. Hodgkiss, “Atmospheric refractivity tracking from radar clutter using kalman and particle filters,” 2007 IEEE Radar Conference, 2007.
  • [14] R. Douvenot, V. Fabbro, P. Gerstoft, C. Bourlier, and J. Saillard, “A duct mapping method using least squares support vector machines,” Radio Sci. Radio Science, vol. 43, no. 6, 2008.
  • [15] C. Yardim, P. Gerstoft, and W. Hodgkiss, “Estimation of radio refractivity from radar clutter using Bayesian Monte Carlo analysis,” IEEE Transactions on Antennas and Propagation, vol. 54, no. 4, p. 1318–1327, 2006.
  • [16] V. Fountoulakis and C. Earls, “Duct heights inferred from radar sea clutter using proper orthogonal bases,” Radio Science, 2016.
  • [17] D. F. Gingras, P. Gerstoft, and N. L. Gerr, “Electromagnetic matched-field processing: Basic concepts and tropospheric simulations,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 10, pp. 1536–1545, 1997.
  • [18] P. Gerstoft, D. F. Gingras, L. T. Rogers, and W. S. Hodgkiss, “Estimation of radio refractivity structure using matched-field array processing,” IEEE Transactions on Antennas and Propagation, vol. 48, no. 3, pp. 345–356, 2000.
  • [19] S. Penton and E. Hackett, “Rough ocean surface effects on evaporative duct atmospheric refractivity inversions using genetic algorithms,” Radio Science, vol. 53, no. 6, pp. 804–819, 2018.
  • [20] J. Pozderac, J. Johnson, C. Yardim, C. Merrill, T. de Paolo, E. Terrill, F. Ryan, and P. Frederickson, “-band beacon-receiver array evaporation duct height estimation,” IEEE Transactions on Antennas and Propagation, vol. 66, no. 5, pp. 2545–2556, 2018.
  • [21] M. Wagner, P. Gerstoft, and T. Rogers, “Estimating refractivity from propagation loss in turbulent media,” Radio Science, vol. 51, no. 12, pp. 1876–1894, 2016.
  • [22] Q. Zhang and K. Yang, “Study on evaporation duct estimation from point-to-point propagation measurements,” IET Science, Measurement & Technology, vol. 12, no. 4, pp. 456–460, 2018.
  • [23] X. Zhao, “Evaporation duct height estimation and source localization from field measurements at an array of radio receivers,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 2, pp. 1020–1025, 2012.
  • [24] X.-F. Zhao, S.-X. Huang, and H.-D. Du, “Theoretical analysis and numerical experiments of variational adjoint approach for refractivity estimation,” Radio Science, vol. 46, no. 1, 2011.
  • [25] F. J. Ryan, User’s Guide for the VTRPE (Variable Terrain Radio Parabolic Equation) Computer Model, 1991.
  • [26] D. E. Kerr, Propagation of short radio waves.   McGraw-Hill, 1951.
  • [27] P. Gerstoft, L. T. Rogers, J. L. Krolik, and W. S. Hodgkiss, “Inversion for refractivity parameters from radar sea clutter,” Radio science, vol. 38, no. 3, 2003.
  • [28] G. Dockery, “Modeling electromagnetic wave propagation in the troposphere using the parabolic equation,” IEEE Trans. Antennas Propagat. IEEE Transactions on Antennas and Propagation, vol. 36, no. 10, p. 1464–1470, 1988.
  • [29] O. Ozgun, G. Apaydin, M. Kuzuoglu, and L. Sevgi, “PETOOL: Matlab-based one-way and two-way split-step parabolic equation tool for radiowave propagation over variable terrain,” Computer Physics Communications, vol. 182, no. 12, p. 2638–2654, 2011.
  • [30] F. B. Jensen, Computational ocean acoustics.   American Inst. of Physics, 1994.
  • [31] G. H. Golub and C. F. Van Loan, Matrix computations.   Johns Hopkins University Press, 1996.
  • [32] L. Kaufman, “Eigenvalue problems in fiber optic design,” SIAM. J. Matrix Anal. & Appl. SIAM Journal on Matrix Analysis and Applications, vol. 28, no. 1, p. 105–117, 2006.
  • [33] J. Nocedal and S. J. Wright, Numerical optimization.   Springer, 2006.
  • [34] L. T. Rogers, “Demonstration of an efficient boundary layer parameterization for unbiased propagation estimation,” Radio Science, vol. 33, no. 6, pp. 1599–1608, 1998.
  • [35] R. Douvenot, V. Fabbro, P. Gerstoft, C. Bourlier, and J. Saillard, “A duct mapping method using least squares support vector machines,” Radio Sci. Radio Science, vol. 43, no. 6, 2008.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description