Abstract
This paper studies an inverse problem that can be used for reconstructing initial wave field of a nonsteady-state wave propagation. The inverse problem is ill posed in the sense that small changes in the input data can greatly affect the solution of the problem. To address the difficulty, the problem is formulated in the form of an inference problem in an appropriately constructed stochastic model. It is shown that the stochastic inverse model enables the initial surface disturbance to be reconstructed, including its confidence intervals given the noisy measurements. The reconstruction procedure is illustrated through applications to some simulated data for two- and three-dimensional problem.
1. Introduction
Water waves are one of important subjects in ocean-related fields. They are generated on the free surface, between air and water, from various sources such as a storm at sea, moving ship in calm water, and impact from a falling object. Once they are generated, wave propagation occurs on the surface of water away from the source. The propagation of water wave is a little more complex than many other kinds of waves in nature due to its dispersion.
In this paper, we are interested in a nonsteady-state water wave propagation generated by an initial surface disturbance. The classical problem in water waves is the computation of the time and spatial evolution in a whole body of water with this given initial wave field. That is, the aim of such problems is to predict the propagation of water waves as effects of the initial free surface disturbance, which is the cause of this phenomenon. This kind of problem has been widely studied because of its importance [1–6]. We name it forward or direct problem. Although much progress has been achieved in the analysis of resulting wave flow, there are only few studies [7–11] on inverse problems, whose aim is to find the cause of the wave propagation.
In this study, our attention is focused on the inverse problem of reconstructing initial wave field from the limited measurement on the specific boundary of the fluid body for two- and three-dimensional problems. Based on the assumption of linear dispersive wave theory, the reconstructing problem can be formulated as an inverse problem. The problem requires solving the Fredholm integral equation of the first kind, which causes severe ill posedness when numerically approximating. As a result, the inverse problem (reconstruction of initial wave field) is mathematically and numerically much more challenging than the forward problem (simulation of future wave field given initial disturbance). For such problems, the solution may not exist when noise is present in the measurements. Moreover, it may also loses stability although it exists, that is, lacks continuous dependence on the data. Thus, a small perturbation of the data may yield erroneous results of the inverse solution.
To overcome the difficulties, we adopt a stochastic inverse method, which has recently attracted much research attention in diverse fields [12–15]. Stochastic inverse methods solve the inverse problem in a systematic way by modeling all parameters of interests as random variables with joint probability distributions. This randomness can be considered as parameter variability since it is related to the uncertainty of the true values. The key idea behind stochastic inverse method is to restate the original inverse problem in the form of probabilistic description of the unknown parameters, that is, stochastic inverse problem. Accordingly, the solution of the inverse problem is represented by a probability distribution of random unknowns. This probability distribution is known as the posterior probability distribution function incorporating the degree of confidence of unknown parameter given the measurement.
Two aspects of this study are worth noting. Firstly, the reconstruction procedure was derived from the limited measurement on the specific boundary of the fluid body. Secondly, a stochastic inverse method is introduced to enable the stable reconstruction procedure, since it requires solving the ill-posed problem. It also provides a way of quantifying an uncertainty in the solution, as well as smoothing the solution by expanding original solution spaces into probability spaces.
This study is organized as follows. Section 2 reviews two- and three-dimensional problems regarding a nonsteady-state water wave propagation due to initial surface disturbance. In Section 3, a procedure is described for the reconstruction of the initial wave field based on the limited measurements at a fixed spatial point. The stochastic inversion method is illustrated in Section 4 for the stable solution procedure. In Section 5, the reconstruction procedure is validated through applications to some digitally simulated data. Finally, we make concluding remarks in Section 6.
2. Water Waves due to an Initial Impulse
2.1. Mathematical Description of the Problem
It is assumed here that the motion of water waves is described by the Euler equations [1–6]. That is, the fluid is homogeneous, incompressible, and inviscid. The water depth is constant , and the surface tension is neglected. In this context, the corresponding fluid motion is governed by the Laplace equation with the linearized boundary conditions where , denote the surface elevation, gravitational acceleration, respectively.
The Euler equations can mathematically be viewed as a Cauchy problem [6], studying a partial differential equation that satisfies certain conditions which are given on a hypersurface in the domain. Given the initial spatial surface disturbance at time zero or , the above equations give a way to determine the whole behavior of fluid with respect to time and space.
2.2. Surface Waves due to Localized Impulse
If the initial wave disturbance , as shown in Figure 1, is generated by the localized sudden disturbances such as underwater explosions or sudden vertical seafloor displacement, the resulting nonsteady-state waves can be described by the axisymmetric form [1–6, 8, 9] where , is the usual Bessel function of the first kind of order zero, is the wave number, and is wave frequency. Two parameters , are related to each other through the dispersion relation . The coefficient function is determined with the initial wave disturbance , which is defined by spatial wave field when , by the Hankel transform: If the Hankel transform of the initial wave disturbance is known, the resulting wave motion can be predictable with (3).
For a two-dimensional system , the problem for the nonsteady-state waves can be simply described [1–6] by where the coefficient function defined by the Fourier transform of is
2.3. Forward and Inverse Problems
If the initial wave disturbance is specified, the coefficient function can be directly obtained by the Hankel transform (4) for the three-dimensional problem and the Fourier transform (6) for the two-dimensional problem. The wave evolution with respect to space and time variables is then simply computed by integrating (3) or (5) for the given coefficient function. The computation of the wave evolution from knowledge of the initial wave disturbance is defined as the forward problem.
On the other hand, the inverse problem of reconstructing initial wave field on the basis of the measured wave elevation on the water surface is not simple. Such inverse procedures often involve ill-posed problems whose solution does not exist or lacks continuous dependence on the given data [7–11]. Thus, the problems are mathematically and numerically much more challenging than the forward problem. In the present study, the aim is focused to reconstruct the initial wave field given a wave history data recorded at a specific spatial point.
3. Reconstruction of Initial Wave Field
3.1. Reconstructing Procedure
The inverse procedure for reconstructing initial wave field can be achieved by two steps: (i) estimating the coefficient function based on the measurement wave elevation at a fixed point, and (ii) transforming the estimated coefficient function into the initial wave field.
Formally, the first step for the reconstruction can be expressed as an operator equation: where denotes the unknown coefficient function with a finite interval and denotes the wave elevation measured at fixed-spatial point during the time interval and the operator is defined by In (8), and represent specific measurement points for 2D problem and 3D problem, respectively.
Let us denote a surface elevation measured at time by ; that is, , . The problem (7) is then approximated by the following uniformly discretized model: where are weights for a quadrature rule, are a vector representation for the coefficient function , that is, ; are an additive random noise, and the elements of are given by
Once the solution is obtained by inverting (9), the initial wave field can be reconstructed by the second step: where , are the inverse Fourier transform and the inverse Hankel transform, respectively. The accent hat symbol means estimated value since it cannot be identical with the true value.
3.2. Instability
It is worth noting that the numerical system (9) is given by discretizing the integral operator (8), which is classified as Fredholm integral operator of the first kind. It is known that the approximation of such integral operator with a regular kernel yields a highly ill-conditioned system [7–11, 16–18]. As a result, the system (9) is ill conditioned and the inverse is greatly affected by small changes in input data (measurement). More specifically, this kind of system causes the following main difficulties in computational procedures [19]: (i) very large condition number of the matrix system, (ii) no guarantee about a useful solution from the replacement into a well-conditioned matrix.
Therefore, standard methods for obtaining the inverse solution cannot be simply used. Instead, more sophisticated methods should be applied to obtain a physically meaningful solution. In this study, we will address the difficulty by formulating the problem as an inference problem in an appropriately constructed stochastic model.
4. Stochastic Model
4.1. Bayesian Formulation
The aim of stochastic modeling is to derive the probability distribution, which is called the posterior probability density function (PPDF), for the unknown parameter conditional on available data according to Bayes’ formula [18]: where is the likelihood function and is the prior density function.
Bayesian formulation solves the ill-posed problem in a systematic way by modeling all parameters of interests as random variables with joint probability distributions. The probability distribution (12) can provide the information on the unknown parameter with the degree of confidence given the measurement.
The PPDF (12) can be specified through the appropriate probabilistic modeling for each component. Since the distribution (12) is proportional to the product of its likelihood function and the prior density , the probabilistic modeling can be achieved by separately modeling these two probability distributions.
If we assume that the random error in the measurement is an additive Gaussian random noise with zero mean and standard deviation , the likelihood function can be modeled [12–15, 18] as where refers to Euclidean norm and is the number of measurements. It is worth emphasizing that other model can also be used to represent the likelihood function. But the Gaussian model (12) is the most commonly used since it is easy to handle and well fit the actual distribution.
Next, we need to consider the prior distribution reflecting the knowledge of the unknowns before the data is measured. For the present study, there is no information about prior probability. For this case, a noninformative prior such as uniform distribution can be a possible choice. However, a noninformative prior may not provide sufficient regularity for the solution of an inverse problem which is severely ill posed. In this study, the following pairwise Markov random field (MRF) model, which has been successfully used in various applications [12–15], is used for the smoothness of the unknown coefficient function : where the matrix is determined by where is the number of neighbors for the point and means that and are adjacent. The MRF is a model for a basis of contextual constraints in the prior distribution. It provides a way of synthesizing and capturing the characteristics of a stochastic process.
With the likelihood (13) and the prior models (14), the PPDF in (12) now can be expressed by
Here, it is worth noting that the expression (16) is dependent on some parameters and , which are used for modeling likelihood function and the prior density. If these parameters are given, then the desired unknown parameter can be estimated by exploring the PPDF (16). However, in most cases, the choice of these parameters is also nondeterministic. In the stochastic model, these parameters are also modeled as random variables and have their own priors. The PPDF is then formulated as the following hierarchical model: where is a pair of gamma distribution for and is a pair of inverse gamma distribution for the prior of .
4.2. Markov Chain Monte Carlo
Now the objective is to extract information on the unknown parameter from the constructed PPDF (17). The resultant posterior distribution can be analyzed using a Markov chain Monte Carlo (MCMC) sampling techniques. The PPDF (17) can be explored by the following hybrid algorithm which is designed by mixing Metropolis-Hastings steps in the Gibbs sampler [12–15]. (i) Initialize , , and ,(ii) For Sample Sample , , , , , Sample , , , , , Sample Sample if , , else Sample Sample if , else where the full conditionals are derived as where represents gamma distribution and represents inverse gamma distribution.
5. Numerical Experiments
5.1. Example I: 2D Problem with Smooth-Type
As the first example, we consider the 2D problem of reconstructing the initial wave field with a smooth-type coefficient function . To illustrate the method, the following particular form of the initial surface disturbance is chosen: The transform (4) of (19) yields the smooth-type coefficient function Figure 2 shows the true initial disturbance (19) and the exact coefficient function (20).
(a)
(b)
The forward problem is first solved with the given initial wave field. The resulting nonsteady-state waves can be simply computed through (5) with the coefficient function (20). Figure 3 shows the computed spatial-time evolution of a nonsteady-state water wave propagation generated by the initial wave disturbance (19). For the forward calculation, the integral operator is directly discretized by Simpson’s integration rule. The water depth is normalized and , are considered as time and space steps on the physical domain . To discretize the integral operator (8), the integration limit and are taken to be 5 and 0.05, respectively.
From now, we consider the inverse problem of reconstructing the initial surface disturbance from the wave history data at a specific spatial point . For illustrating purpose, a particular time history of wave elevation, , is collected at and is shown in Figure 4. To mimic a practical situation, the Gaussian random noise is generated with the density and added to the collected data. The resulting data with a random noise is considered as input data for the inverse problem.
As explained earlier, the system is severely ill posed and thus the inverse is mathematically and numerically much more challenging than the forward calculation. To examine the cause and the degree of instability of the system, the discretized system is decomposed as where ’s are the singular values ordered, , and , are orthonormal singular vectors. The preceding decomposition is well known as the singular value decomposition [16–19]. The unknown can then be obtained by the least squares estimate, that is, :
With the measurement data in Figure 4, we first obtained the inverse solution using the least squares estimate (22). It can be seen from Figure 5 that the computed inverse solutions are totally unstable and thus useless. The cause and the degree of instability in the inverse solution can be understood by analyzing the behavior of singular systems, more specifically, each component in the right-hand-side of (22). Figure 6 shows the Picard plot [19] that shows behaviors of , , and the ratio . It can be found that the coefficients become larger than most of singular values . Thus, the least squares solution for , which is the sum of the ratio , is dominated by the small singular values. This explains the reason for the unstable solution in Figure 5. The instability, frequent sign changes, may be accelerated and amplified with the increasing amount of noise in data.
With this analysis in mind, we next construct the stochastic inverse model (17) in terms of the measured wave history data. To extract useful information from the probability density, MCMC algorithm described in Section 4 is used. The initial guess for , , and are taken to be 1, 0.1, and zero vector, respectively. The pairs of parameters , are taken to be , , respectively. The total number of samples for the algorithm was taken to be 50,000 and the last 25,000 samples were used to estimate statistics for since it allows the beginning of the Markov chain. After obtaining random samples from the MCMC algorithm, it is necessary to examine the issues such as convergence and mixing of the chain to determine whether simulated set of samples can reasonably be treated as a set of random realizations from the target distribution. In this study, marginal trace plots are used to check if the chain works properly. All the components have the stationary distribution as its marginal distribution.
The true and estimated coefficient functions are shown in Figure 7. The upper and lower dotted lines denote the 95% credible interval. The credible interval quantifies the degree of uncertainties in the solution. The mean estimate of the posterior density is fairly stable and accurate compared with the results in Figure 8. Finally, the initial wave field is reconstructed using the posterior mean of the coefficient function. The reconstructed initial surface disturbance in Figure 8 represents accurate approximation to the true initial disturbance.
5.2. Example II: 3D Problem with Nonsmooth-Type
As the second example, we consider the 3D problem of reconstructing the initial wave field with a nonsmooth-type coefficient function , which is believed to be a more challenging case than the smooth type. For the purpose, we first specify the initial surface disturbance for 3D problem as where is Bessel function of the first kind of order one and . The transform of (21) yields the nonsmooth type coefficient function: Figure 9 shows the true initial disturbance (21) and the exact coefficient function (22).
(a)
(b)
In the similar manner of the previous example, the forward problem is first solved with the given initial surface disturbance. For the calculation, the water depth is normalized and , are considered as time and space steps on the physical domain . When discretizing the integral operator (8), the integration limit and are taken to be 5 and 0.05, respectively. Figure 10 shows the computed spatial-time evolution for 3D problem.
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
| (f) |
The wave history data in Figure 11, which is recorded at , is used to illustrate the reconstructing procedure. To generate noisy data, the random noise with the variance is added to the collected data. With this measurement, the coefficient function is estimated through the MCMC algorithm. The same conditions as in the first example are used for the MCMC algorithm. The posterior mean estimate of the simulated random samples is shown in Figure 12. For the purpose of comparison, the least squares estimation of the coefficient function is also illustrated. The result indicates that the posterior mean is in good agreement with the exact coefficient function while the result from least squares estimate is totally unstable. Using this posterior mean estimate, the initial wave field is reconstructed and shown in Figure 13. It can be found that the estimated result accurately approximate to the true initial disturbance.
(a)
(b)
It is also worth to compare the performances of the existing algorithms and the stochastic inverse approach illustrated in this paper to underline the advantages and disadvantages. For this purpose, Tikhonov regularization method [8, 11, 16–19], which is known as a representative deterministic inverse approach, was employed to solve the second example with the same conditions. Tikhonov regularized solution for is the minimizer of the functional where is the regularization parameter.
In deterministic inverse approach, the choice of the regularization parameter is very important because it controls the performance of the regularization. However, selecting an appropriate regularization parameter is not trivial since it depends on both the regularization method itself and the mathematical model of interest in a very complex way. Generally, the regularization parameter is selected by the heuristic method. First, the inverse problem is solved with a set of regularization parameters. The solution remains practically unchanged for a certain range of the parameter. The regularization parameter can then be chosen from this range. Figure 14 shows the Tikhonov regularized solutions to three different regularization parameters. It can be observed that the poor solution is obtained when the regularized solution to the regularization parameter is large. Furthermore, the solution loses its stability when the regularization parameter is small. Therefore, it is necessary to pay attention to the behavior of the regularized solution, when choosing the regularization parameter. It is seen by comparing with Figure 12 that the proposed stochastic inverse approach is smoother than the Tikhonov method.
| (a) |
| (b) |
| (c) |
Compared with the deterministic inverse methods for inverse problems, the proposed stochastic inverse approach has some advantages. Firstly, it provides a natural framework for quantifying the uncertainties by considering a complete probabilistic description of the inverse problem of interest. Generally, there exist many possible candidates for the solution, because of the ill-posed nature of the inverse problem. But, the deterministic inverse methods provide only a single estimate taken from these possible candidates based on the choice of the regularization parameter. Secondly, it also provides a flexible way to choose a nondeterministic parameter by formulating the hierarchical model.
However, there also exist disadvantages. In this work, a computational framework for reconstructing the initial wave field was built based on the Markov chain Monte Carlo technique. In computational mathematics, the efficient computation is important. The key behind the Markov chain is to generate samples from the target distribution based on the iteration. To ensure the chain convergence, quite large numbers of samples are required. Thus, it is expected that the solver will be computationally very expensive because the computational algorithm of the stochastic inverse model requires the direct simulation. In contrast, the deterministic inverse methods can remarkably reduce the computational time by avoiding the time-consuming iterations.
6. Conclusions
The conclusions can be summarized as follows.(1)The fluid motion has been taken to be homogeneous, incompressible, and inviscid such that the velocity potential exists and satisfies Laplace’s equation. On this basis, it has been shown that, for a given initial surface disturbance, the surface wave motion can be mathematically expressed as an integral transform, which can be used to compute the time-spatial evolution of the wave motion. (2)The inverse problem of reconstructing the initial surface disturbance was formulated with a time history data of surface elevation which is collected at a specific spatial point. The inverse problem requires solving an ill-posed system that causes numerical instability. It was shown that, for such system, a standard method such as least squares estimate fails to obtain the stable solution. To address the difficulty, the system was formulated as an inference problem in an appropriately constructed stochastic model. (3)The reconstruction procedure based on the stochastic inverse model has been evaluated through an analysis of some digitally simulated data for numerical examples related to two- and three-dimensional fluid motion. It has been shown that the proposed procedure can be used to not only reconstruct the initial surface disturbance but also detect uncertainties in solution arising from measurement error.In this study, the quality of the reconstruction has been illustrated by two examples, where the simulated data fulfilled the assumptions about the noise. Future work on this study will also include the application of real data for ensuring robustness of the method.