Cone Bearing Estimation Utilizing a Hybrid HMM and IFM Smoother Filter Formulation

Erick Baziw^1*, Gerald Verbeek²

¹Baziw Consulting Engineers, Vancouver, BC, Canada.
²Baziw Consulting Engineers, Tyler, Texas, USA.

DOI: 10.4236/ijg.2021.1211055                

Abstract

Cone penetration testing (CPT) is a widely used geotechnical engineering in-situ test for mapping soil profiles and assessing soil properties. In CPT, a cone on the end of a series of rods is pushed into the ground at a constant rate and resistance to the cone tip is measured (q_m). The q_m values are utilized to characterize the soil profile. Unfortunately, the measured cone tip resistance is blurred and/or averaged which can result in the distortion of the soil profile characterization and the inability to identify thin layers. This paper outlines a novel and highly effective algorithm for obtaining cone bearing estimates q_t from averaged or smoothed q_m measurements. This q_t optimal filter estimation technique is referred to as the q_tHMM-IFM algorithm and it implements a hybrid hidden Markov model and iterative forward modelling technique. The mathematical details of the q_tHMM-IFM algorithm are outlined in this paper along with the results from challenging test bed. The test bed simulations have demonstrated that the q_tHMM-IFM algorithm can derive accurate q_t values from challenging averaged q_m profiles. This allows for greater soil resolution and the identification and quantification of thin layers in a soil profile.

Keywords

Bayesian Recursive Estimation (BRE), Cone Penetration Testing (CPT), Geotechnical Site Characterization, Hidden Markov Model (HMM), Iterative Forward Modelling (IFM), Smoothing

1. Introduction

Cone Penetration Test (CPT) [1] [2] [3] [4] is extensively used in geotechnical engineering to determine the in-situ subsurface stratigraphy and to estimate geotechnical parameters of the soils present. Geotechnical engineers use CPT to characterize and quantify soil properties and groundwater conditions so that the infrastructure (e.g., bridges, roads, buildings) construction requirements can be determined. In CPT a cone penetration test rig pushes the steel cone vertically into the ground at a standard rate and data are recorded at regular intervals during penetration. The cone penetrometer has electronic sensors to measure penetration resistance at the tip (q_m) and friction in the shaft (friction sleeve) during penetration. A CPT probe equipped with a pore-water pressure sensor is called a piezo-cones (CPTU cones). CPT penetrometers with other sensors such as a seismic sensor are also used for in-situ site characterization. Figure 1 illustrates a schematic and the associated terminology of a cone penetrometer.

For piezo-cones with the filter element right behind the cone tip (i.e., the u2 position), it is standard practice to correct the recorded tip resistance for the impact of the pore pressure on the back of the cone tip. This corrected cone tip resistance is normally referred to as q_t, but in this paper, we focus on an additional correction that should be made to address the averaging that takes place when performing CPT to obtain the actual cone tip resistance values.

Boulanger and DeJong [5] outlined the distortions which occur when obtaining q_m measurements and proposed an “inverse” algorithm where the results of the distortion could be optimally removed. In their work Boulanger and DeJong incorrectly described the distortions as a convolution operation (Equations (1), (2), (10), (12), (13), and (15)). In fact, the tip-bearing distortions are an averaging process [2] [5] where cone tip values measured at a particular depth are affected by values above and below the depth of interest. This averaging or smoothing results in the inability to identify thin layers which is critical for liquefaction assessment and the reduction in soil layer resolution. The averaging/smoothing effect is subsequently described along with a proposed algorithm which combines the Bayesian recursive estimation Hidden Markov Model (HMM) filter with Iterative Forward Modelling (IFM) parameter estimation in a smoother formulation for optimal estimation of true q_t cone bearing values.

Figure 1. Schematic and terminology for cone penetrometer ( [1] ).

2. Mathematical Background

2.1. Cone Penetration Testing Model

When performing CPT the layers above and below the cone tip affect the measured tip resistance as illustrated in Figure 2.

The measured cone penetration tip resistance q_m can then be described as

$\begin{array}{l}{q}_m\left(d\right)={\displaystyle \underset{j=1}{\overset{60\times \left(\frac{C_d}{\Delta }\right)}{\sum }}{w}_c\left(j\right)\times q_t\left({\Delta }_{qt}+j\right)}+v\left(d\right)\\ {\Delta }_{qt}=\left(d-{\Delta }_{wc}\right),{\Delta }_{wc}=30\times \left(\frac{d_c}{\Delta }\right)\end{array}$ (1)

where

d: the cone depth;

d_c: the cone tip diameter;

Δ: the q_t sampling rate;

q_m(d): the measured cone penetration tip resistance;

q_t(d): the true cone penetration tip resistance;

w_c(d): the q_t(d) averaging function;

v(d): additive noise, generally taken to be white with a Gaussian pdf.

In Equation (1) it is assumed that w_c averages q_t over 60 cone diameters centered at the cone tip. Boulanger and DeJong [5] outline how to calculate w_c using the equations shown below (after correcting the equation for w₁).

The cone penetration averaging function w_c for varying $q_{t,z^{\prime }}/q_{t,z^{\prime }=0}$ ratios is illustrated in Figure 3. As outlined out by Boulanger and DeJong [5], w_c is highly nonlinear and depth variant. This proves to be a significant challenge in obtaining the optimal estimates of q_t but this can be addressed by applying a BRE filter which incorporates smoothing.

Figure 2. Schematic of thin layer effect for a sand layer embedded in a clay layer [5].

Figure 3. Normalized cone penetration filter versus normalized depth from the cone tip ( $q_{t,z^{\prime }}/q_{t,z^{\prime }=0}=0.01,0.10,10$ and 100) [5].

$w_c=\frac{w_1{w}_2}{{\displaystyle \sum w_1{w}_2}}$ (2a)

$w_1=\frac{C_1}{1+\left|{\left(\frac{z^{\prime }}{{z^{\prime }}_{50}}\right)}^{m_z}\right|}$ (2b)

$w_2=\sqrt{\frac{2}{1+{\left(\frac{q_{t,z^{\prime }}}{q_{t,z^{\prime }=0}}\right)}^{m_q}}}$ (2c)

2.2. Bayesian Recursive Estimation

Bayesian Recursive Estimation (BRE) is a filtering technique based on state-space, time-domain formulations of physical problems [6] [7]. Application of this filter type requires that the dynamics of the system and measurement model, which relates the noisy measurements to the system state equations, be describable in a mathematical representation and probabilistic form that uniquely define the system behaviour. The potentially nonlinear discrete stochastic equation describing the system dynamics is defined as follows:

$x_k=f_{k-1}\left(x_{k-1},u_{k-1}\right)\leftrightarrow p\left(x_k|x_{k-1}\right)$ (3)

In Equation (3), the vector f_k is a function of the state vector x_k and the process or system noise u_k. It is assumed that Equation (3) describes a Markov process of order one. The sampled potentially nonlinear measurement equation is given as

$z_k=h_k\left(x_k,v_k\right)\leftrightarrow p\left(z_k|x_k\right)$ (4)

In Equation (4), h_k depends upon the index k, the state x_k, and the measurement noise v_k at each sampling time. The probabilistic state-space formulation described by Equation (3) and the requirement for updating the state vector estimate based upon the newly available measurements described by Equation (4) are ideally suited for the Bayesian approach to derive the optimal estimation. In this approach it is attempted to construct the posterior estimate of the state given all available measurements. In general terms, it is desired to obtain estimates of the discretized system equation states x_k, based on all available measurements up to time k (denoted as z_1:k) by constructing the posterior p(x_k|z_1:k). The posterior Probability Density Function (PDF) then allows the calculation of the conditional mean estimate of the state (E(x_k|z_1:k)).

BRE is a two step process consisting of prediction and update. In the prediction step the system equation defined by Equation (3) is used to obtain the prior PDF of the state at time k using the Chapman-Kolmogorov equation, which is given as

$p\left(x_k|z_{1:k-1}\right)={\displaystyle \int p\left(x_k|x_{k-1}\right)p\left(x_{k-1}|z_{1:k-1}\right)\text{d}{x}_{k-1}}$ (5)

The update step then computes the posterior PDF from the predicted PDF and the newly available measurement as follows:

$p\left(x_k|z_{1:k}\right)=\frac{p\left(z_k|x_k\right)p\left(x_k|z_{1:k-1}\right)}{p\left(z_k|z_{1:k-1}\right)}$ (6)

The recurrence Equations (5) and (6) form the basis for the optimal Bayesian solution. The BRE of the posterior density can generate an exact solution when the state-space equations fit into a Kalman Filter (KF) formulation or a Hidden Markov Model (HMM). Otherwise, BRE will generate an estimation numerically using Particle Filters (PF) when deriving the posterior PDF.

2.3. Hidden Markov Model (HMM) Filter

The HMM filter (also termed a grid-based filter) has a discrete state-space representation and has a finite number of states. In the HMM filter the posterior PDF is represented by the delta function approximation as follows:

$p\left(x_{k-1}|z_{1:k-1}\right)={\displaystyle {\sum }_{i=1}^{N_s}{w}_{k-1\backslash k-1}^i\delta \left(x_{k-1}-x_{k-1}^i\right)}$ (7)

where $x_{k-1}^i$ and $w_{k-1|k-1}^i$, $i=1,\cdots ,N_s$, represent the fixed discrete states and associated conditional probabilities, respectively, at time index k − 1, and N_s the number of particles utilized. The governing equations for the HMM filter are derived by substituting Equation (7) into Equations (5) and (6). This substitution results in the HMM prediction and update equations which are outlined in Table 1.

2.4. Iterative Forward Modelling

Iterative forward modeling (IFM) is a parameter estimation technique which is

Table 1. HMM governing equations.

based upon iteratively adjusting the parameters until a user specified cost function is minimized. The desired parameter estimates are defined as those which minimize the user specified cost function. The IFM technique which is utilized within the q_t estimation algorithm is the downhill simplex method (DSM) originally developed by Nelder and Mead [8]. The DSM in multidimensions has the important property of not requiring derivatives of function evaluations and it can minimize nonlinear-functions of more than one independent variable. Although it is not the most efficient optimization procedure, the DSM is versatile, robust and simple to implement. A simplex defines the most elementary geometric figure of a given dimension: a line in one dimension, the triangle in two dimensions, the tetrahedron in three, etc; therefore, in an N-dimensional space, the simplex is a geometric figure that consists of N + 1 fully interconnected vertices. For example, in determining the location of a seismic event, a three-dimensional space is searched, so the simplex is a tetrahedron with four vertices. The DSM has been used in a variety of scientific applications such as obtaining seismic source locations [9] [10] ) tomographic imaging [9], and blind seismic deconvolution [11].

The DSM starts at N + 1 vertices that form the initial simplex. The initial simplex vertices are chosen so that the simplex occupies a good portion of the solution space. In addition, it is also required that a scalar cost function be specified at each vertex of the simplex. The general idea of the minimization is to keep the minimum within the simplex during the optimization, at the same time decreasing the volume of the simplex. The DSM searches for the minimum of the costs function by taking a series of steps, each time moving a point in the simplex away from where the cost function is largest. The simplex moves in space by variously reflecting, expanding, contracting, or shrinking. The simplex size is continuously changed and mostly diminished, so that finally it is small enough to contain the minimum with the desired accuracy. The DSM incorporates the following basic steps:

1) Specify initial simplex vertices.

2) Specify the cost function at each vertex of the simplex.

3) Compare the cost function for each vertex and determine the lowest error “best” and highest error “worst” vertices.

4) Sequentially locating first the reflected, then if necessary, the expanded, and then if necessary, the contracted vertices, and calculating for each the corresponding cost function and comparing it to the worst vertex; if at any step the cost function of the new trial point is less than the value at the worst vertex; then this vertex is substituted as a vertex in place of the current worst vertex.

5) If the process in step 4 does not yield a lower error value than the previous worst, then the other vertices are shrunken towards the best vertex.

6) At each stage of shrinking, the distances between vertices are calculated and compared to a set tolerance value to check if the simplex has become sufficiently small for termination of the estimation; when the test criterion is reached, the previous best vertex becomes the solution.

7) At each stage of shrinking, the cost function values at the vertices is compared to a set minimum value to check if the error residual has become sufficiently small for termination of the estimation; when the test criterion is reached, the previous best vertex becomes the solution.

3. q_tHMM-IFM Algorithm

The q_t optimal filter estimation technique is referred to as the q_tHMM-IFM algorithm and it consist of a BRE smoother and an IFM component.

3.1. q_tHMM Algorithm Formulation

The HMM portion of the q_tHMM-IFM algorithm (so called q_tHMM algorithm) implements a BRE smoothing filter. BRE smoothing is a non-real time filter that uses all measurements available to estimate the state of a system at a certain time or depth in the q_t estimation case [6] [7]. The q_tHMM algorithm smoother consists of two parts: the forward and the backward formulation. The forward-depth filter ( ${\hat{q}}_k^F$ ) processes measurement data (q_m) above the cone tip ( $j=1$ to

$30\times \left(\frac{d_c}{\Delta }\right)$ in (1)). Next a backward-depth formulation ( ${\hat{q}}_k^B$ ) is implemented, where the filter recurses through the data below the cone tip ( $j=30\times \left(\frac{d_c}{\Delta }\right)$ to $60\times \left(\frac{d_c}{\Delta }\right)$ in (1)) starting at the final q_m value. The optimal estimate for q_t is then defined as

${\hat{q}}_k^t=\left({\hat{q}}_k^F+{\hat{q}}_k^B\right)/2$ (12)

Since the structure of the backward-depth q_tHMM filter is similar to that of the forward-depth q_tHMM filter only the forward-depth formulation is outlined. In this case a HMM filter is utilized where a bank of discrete q_t values (i = 1 to N) varying from low (q_tL) to high (q_tH) (e.g., 0 MPa to 120 MPa) and a corresponding q_t resolution q_tR (e.g., 0.2 MPa) are specified. The required number of fixed grid HMM states is given as N_S = (q_tH − q_tL)/q_tR. In Table 1 the notation of the states xⁱ is mapped to qⁱ to reflect the bank of q_t values. The measurement equation given by Equation (1) is modified as outlined below for the forward-depth case:

$\begin{array}{l}{z}_k^i={\displaystyle \underset{j=1}{\overset{30\times \left(\frac{d_c}{\Delta }\right)}{\sum }}{w}_c\left(j\right)\times q_k^i\left({\Delta }_{qt}+j\right)}+v_k\\ {\Delta }_{qt}=\left(d-{\Delta }_{wc}\right),{\Delta }_{wc}=30\times \left(\frac{d_c}{\Delta }\right)\end{array}$ (13)

The transitional probabilities (i.e., $p\left(x_k^i|x_{k-1}^j\right)$  or $p\left(q_k^i|q_{k-1}^j\right)$ ) for each HMM state (i.e., discrete cone tip, q^,i) is set equal due to the fact that there is equal probability of moving from a current cone tip value to any other value between the range q_tL to q_tH. The likelihood PDF $p\left(z_k|q_k^i\right)$ in the HMM filter outlined in Table 1 is calculated based upon an assumed Gaussian measurement error as follows:

$p\left(z_k|q_k^i\right)=\frac{1}{\sqrt{2\pi \sigma }}{\text{e}}^{\left[-\frac{q_m\left(d\right)-z_k^i}{2{\sigma }^2}\right]}$ (14)

where σ² is the variance of the measurement noise. The HMM forward-depth estimated q_tHMM cone tip values ( ${\hat{q}}_k^F$ ) are calculated as follows:

${\hat{q}}_k^F={\displaystyle {\sum }_{i=1}^{N_s}{w}_{k|k}^i{q}_k^i}$ (15)

3.2. q_tHMM Test Bed Example

The previously outlined q_tHMM algorithm was subjected to extensive test bed simulations. Unfortunately, it was concluded that this formulation would not work due to the significant challenges in estimating the unknown q_t values below the cone tip (e.g., 54 unknown q_t values for the case of a 10 cm² cone tip (d_c = 0.0357 m) and a sampling rate of 0.02 m (∆ = 0.02) for the forward-depth q_tHMM Formulation. The same obviously applied for the backward-depth q_tHMM Formulation, but in that case the 54 unknown q_t values are above the cone tip. It should be noted, however, that with this q_tHMM Formulation it was nearly possible to duplicate the results of Boulanger and DeJong shown in Figure 4. In this case, the forward-depth formulation was utilized and the maximum q_t value was not allowed to exceed 100. This is illustrated below for the case outlined in Figure 4(b) and H_sand/d_c = 20.

Figure 5 illustrates values of q_m and the forward-depth q_t estimates for varying maximum q_t values specified for the case outlined in Figure 4(b) and H_sand/d_c = 20. As is shown in Figure 5(a) the estimated q_t values closely match

Figure 4. Values of q_t, q_m and q_inv (estimated q_t value) for idealize profiles with interlayers of a strong soil with q_t = 100 embedded in a weaker soil with: (a) q_t = 1, (b) q_t = 10, and (c) q_t = 50. Note that q_inv almost perfectly overlays q_t [5].

Figure 5. Values of q_m (black series) and q_tHMM feed-forward q_t (red series) estimates for varying q_t maximum values specified for the case outlined in Figure 4(b) and H_sand/d_c = 20. (a) q_t maximum = 100, (b) q_t maximum = 90, (c) q_t maximum = 120 and (d) q_t maximum = 250.

those of the true q_t values (i.e. 100) if the specified maximum is 100. However, if the specified maximum value is changed to 250 the results (as illustrated in Figure 5(d)) are far from impressive. This clearly suggests that the good correlation achieved in Figure 5(a) is due to the restrictions placed upon the maximum allowable q_t value and not a demonstration of algorithm performance. It should be noted that the output is shown in Figure 5(d) is similar to the results of Boulanger and DeJong illustrated in Figure 6. This figure illustrates the significant instability in the estimates of q_t when Boulanger and DeJong used their inversion estimation algorithm for a situation where soil layers with a q_t value of 12 were embedded in a uniform deposit with a q_t value of 10. This instability caused Boulanger and DeJong to incorporate an ad-hoc smoothing filter followed by a low-pass spatial filter into their algorithm.

3.3. Incorporation of IFM into the q_tHMM Algorithm

IFM is incorporated into the q_tHMM algorithm to address the poor test bed results. In this case, initial estimates for q_t are derived utilizing IFM. Instead of attempting to estimate all the unknown q_t values (below the cone depth for the forward-depth analysis, and above the cone for the backward-depth analysis) IFM is utilized where only a fraction of the q_t values are required to be estimated. In this process constant layer q_t values and their corresponding depth extents are estimated for a maximum number of layers (specified by the user) within the next w_c window.

Figure 6. Significant instability in the estimates of q_t when using the Boulanger and DeJong inversion estimation algorithm [5].

As an example, assuming that a 10 cm² cone is utilized for the sounding (with a diameter of 36 mm), then the extent of the w_c averaging window (equal to 60 cone diameters) is approximately 2 m and the depth interval below the cone for the forward-depth analysis is approximately 1 m. Assuming that a maximum of three possible layers exist within this depth interval, then only values for d₁, d₂, q₂ and q₃ have to be estimated, as the value of q₁ is estimated with forward-depth formulation of the q_tHMM-IFM algorithm where HMM filter transitional probabilities are taken into account. This is only 4 parameters as opposed to 54 q_t values without the incorporation of IFM into q_tHMM estimation algorithm. The initial estimates derived with IFM are then used as the base line for subsequent iterations using Equation (12). The process is then repeated until a desired error residual is obtained or until a prespecified number of iterations has been reached (Figure 7).

3.4. q_tHMM-IFM Test Bed Example

The performance of the q_tHMM-IFM algorithm was evaluated by carrying out challenging test bed simulations. This section outlines two of these challenging test bed simulations. The first test bed simulation of which is illustrated in Figure 8. A soil profile was defined through q_t values (light grey line in Figure 8(a)) and the resulting q_m values were then calculated (black line in Figure 8(a)). Using the q_tHMM-IFM algorithm the q_t values were then estimated based on the q_m values (black dotted line in Figure 8(a)). It shall be obvious that the algorithm performed well as the derived q_t values closely matched the originally specified q_t values (with the percentage difference shown by the black line in Figure 8(b)). Interesting to note in this simulation is that the layering identified by the black oval (i.e., variation in q_t between 84 and 96) has been lost, and that when focusing on q_m values thin soil layers are completely overlooked. In addition inserting these thin layers also results in large differences between the measured and the true tip resistance values for the entire interval where these thin layers occur, as shown by the black dotted line in Figure 8(b).

A second test bed simulation of the performance of the q_tHMM-IFM algorithm is shown in Figure 9. In Figure 9 a soil profile was defined through q_t values (grey line in Figure 9(a)) and the resulting q_m values were then calculated

Figure 7. S Schematic outlining proposed IFM portion of the q_tHMM-IFM algorithm parameters to be estimated. Parameters q₁, q₂, and q₃ denote the cone bearing values for soil layers 1, 2 and 3, respectively. Parameters d₁, d₂, and d₃ denote the depths of soil layers 1, 2 and 3, respectively.

(a)(b)

Figure 8. TEST BED 1 (a) Specified q_t values (grey line), derived q_m values (black line) and estimated q_t values based on q_m values (black dotted line). (b) Percent differences between specified and estimated q_t values (black line) and q_m values and estimated q_t values (black dotted line).

(a)(b)

Figure 9. TEST BED 2 (a) Specified q_t values (grey line), derived q_m values (black line) and estimated q_t values based on q_m values (black dotted line). (b) Percent differences between specified and estimated q_t values (black line) and q_m values and estimated q_t values (black dotted line).

(black line in Figure 9(a)). Using the q_tHMM-IFM algorithm the q_t values were then estimated based on the q_m values (black dotted line in Figure 9(a)). It shall be obvious that the algorithm performed well as the derived q_t values closely matched the originally specified q_t values (with the percentage difference shown by the black line in Figure 9(b)). Similar to test bed 1, the layering identified by the black oval has been lost and that when focusing on q_m values thin soil layers are completely overlooked. The light grey ovals in Figure 9(a) clearly illustrate that by applying q_m values the actual tip resistance values are masked and blurred. In addition inserting these thin layers also results in large differences between the measured and the true tip resistance values for the entire interval where these thin layers occur, as shown by the back dotted line in Figure 9(b).

4. Conclusions

Cone penetrometer testing (CPT) is an effective, fast and relatively inexpensive system for determining the in-situ subsurface stratigraphy and estimating geotechnical parameters of the soils present. When performing CPT the layers above and below the cone tip affect the measured tip resistance. The extent of this issue can be significant and is dependent upon variable in-situ soil properties (i.e., it is site specific). As a result, a specific soil (with a specific q_t value) can generate significantly different tip resistance readings (q_m value) based upon the properties of the bounding soils, especially in soil profiles with thin soil layers. For this reason, it is recommended that an algorithm is implemented to generate the actual q_t value from recorded q_m values.

This paper has outlined an algorithm which utilizes a hybrid HMM and IFM filter for the purpose of obtaining CPT true cone tip bearing values from measured blurred measured values. The q_t estimation algorithm is referred to as q_tHMM-IFM. Challenging test bed simulations have demonstrated that the q_tHMM-IFM algorithm can derive accurate q_t values from a q_m profile. This allows for the identification and quantification of thin layers in a soil profile. The authors will carry out further test bed simulations and subsequently apply the q_tHMM-IFM algorithm on real data sets.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lunne, T., Robertson, P.K. and Powell, J.J.M. (1997) Cone Penetrating Testing: In Geotechnical Practice. Taylor & Francis, 1997.
[2]	Robertson, P.K. (1990) Soil Classification Using the Cone Penetration Test. Canadian Geotechnical Journal, 27, 151-158. https://doi.org/10.1139/t90-014
[3]	(2017) ASTM D6067/D6067M-17 Standard Practice for Using the Electronic Piezocone Penetrometer Tests for Environmental Site Characterization and Estimation of Hydraulic Conductivity. Soil and Rock, 4, 324-333.
[4]	Cai, G.J., Liu, L.Y., Tong and Du, G.Y. (2006) General Factors Affecting Interpretation for the Piezocone Penetration Test (CPTU) Data. Journal of Engineering Geology, 14, 632-636.
[5]	Boulanger, R.W. and DeJong, T.J. (2018) Inverse Filtering Procedure to Correct cone Penetration Data for Thin-Layer and Transition Effects. In: Hicks, M.A., Pisano, F. and Peuchen, J., Eds., Cone Penetration Testing 2018, CRC Press, London, 25-44.
[6]	Arulampalam, M.S., Maskell, S. and Clapp, T. (2002) A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking. IEEE Transactions on Signal Processing, 50, 174-188. https://doi.org/10.1109/78.978374
[7]	Baziw, E. (2007) Application of Bayesian Recursive Estimation for Seismic Signal Processing, Ph.D. Thesis, University of British Columbia, Columbia, Canada.
[8]	Nelder, J.A. and Mead, R. (1965) A Simplex Method for Function Optimization. Computing Journal, 7, 308-313. https://doi.org/10.1093/comjnl/7.4.308
[9]	Gibowicz, S.J. and Kijko, A. (1994) An Introduction to Mining Seismology. Academic Press, CA.
[10]	Baziw, E., Nedilko, B. and Weir-Jones, I. (2004) Microseismic Event Detection Kalman Filter: Derivation of the Noise Covariance Matrix and Automated First Break Determination for Accurate Source Location Estimation. Pure and Applied Geophysics, 161, 303-329. https://doi.org/10.1007/s00024-003-2443-8
[11]	Baziw, E. (2011) Incorporation of Iterative Forward Modeling into the Principle Phase Decomposition Algorithm for Accurate Source Wave and Reflection Series Estimation. IEEE Transactions on Geoscience and Remote Sensing, 49, 650-660. https://doi.org/10.1109/TGRS.2010.2058122