Back
 IJG  Vol.12 No.11 , November 2021
Cone Bearing Estimation Utilizing a Hybrid HMM and IFM Smoother Filter Formulation
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 (qm). The qm 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 qt from averaged or smoothed qm measurements. This qt optimal filter estimation technique is referred to as the qtHMM-IFM algorithm and it implements a hybrid hidden Markov model and iterative forward modelling technique. The mathematical details of the qtHMM-IFM algorithm are outlined in this paper along with the results from challenging test bed. The test bed simulations have demonstrated that the qtHMM-IFM algorithm can derive accurate qt values from challenging averaged qm profiles. This allows for greater soil resolution and the identification and quantification of thin layers in a soil profile.

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 (qm) 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 qt, 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 qm 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 qt 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 qm can then be described as

q m ( d ) = j = 1 60 × ( C d Δ ) w c ( j ) × q t ( Δ q t + j ) + v ( d ) Δ q t = ( d Δ w c ) , Δ w c = 30 × ( d c Δ ) (1)

where

d: the cone depth;

dc: the cone tip diameter;

Δ: the qt sampling rate;

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

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

wc(d): the qt(d) averaging function;

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

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

The cone penetration averaging function wc for varying q t , z / q t , z = 0 ratios is illustrated in Figure 3. As outlined out by Boulanger and DeJong [5], wc is highly nonlinear and depth variant. This proves to be a significant challenge in obtaining the optimal estimates of qt 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 / q t , z = 0 = 0.01 , 0.10 , 10 and 100) [5].

w c = w 1 w 2 w 1 w 2 (2a)

w 1 = C 1 1 + | ( z z 50 ) m z | (2b)

w 2 = 2 1 + ( q t , z q t , z = 0 ) 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 ( x k 1 , u k 1 ) p ( x k | x k 1 ) (3)

In Equation (3), the vector fk is a function of the state vector xk and the process or system noise uk. 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 ( x k , v k ) p ( z k | x k ) (4)

In Equation (4), hk depends upon the index k, the state xk, and the measurement noise vk 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 xk, based on all available measurements up to time k (denoted as z1:k) by constructing the posterior p(xk|z1:k). The posterior Probability Density Function (PDF) then allows the calculation of the conditional mean estimate of the state (E(xk|z1: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 ( x k | z 1 : k 1 ) = p ( x k | x k 1 ) p ( x k 1 | z 1 : k 1 ) 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 ( x k | z 1 : k ) = p ( z k | x k ) p ( x k | z 1 : k 1 ) p ( z k | z 1 : k 1 ) (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 ( x k 1 | z 1 : k 1 ) = i = 1 N s w k 1 \ k 1 i δ ( x k 1 x k 1 i ) (7)

where x k 1 i and w k 1 | k 1 i , i = 1 , , N s , represent the fixed discrete states and associated conditional probabilities, respectively, at time index k − 1, and Ns 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 qt 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. qtHMM-IFM Algorithm

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

3.1. qtHMM Algorithm Formulation

The HMM portion of the qtHMM-IFM algorithm (so called qtHMM 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 qt estimation case [6] [7]. The qtHMM algorithm smoother consists of two parts: the forward and the backward formulation. The forward-depth filter ( q ^ k F ) processes measurement data (qm) above the cone tip ( j = 1 to

30 × ( d c Δ ) in (1)). Next a backward-depth formulation ( q ^ k B ) is implemented, where the filter recurses through the data below the cone tip ( j = 30 × ( d c Δ ) to 60 × ( d c Δ ) in (1)) starting at the final qm value. The optimal estimate for qt is then defined as

q ^ k t = ( q ^ k F + q ^ k B ) / 2 (12)

Since the structure of the backward-depth qtHMM filter is similar to that of the forward-depth qtHMM filter only the forward-depth formulation is outlined. In this case a HMM filter is utilized where a bank of discrete qt values (i = 1 to N) varying from low (qtL) to high (qtH) (e.g., 0 MPa to 120 MPa) and a corresponding qt resolution qtR (e.g., 0.2 MPa) are specified. The required number of fixed grid HMM states is given as NS = (qtHqtL)/qtR. In Table 1 the notation of the states xi is mapped to qi to reflect the bank of qt values. The measurement equation given by Equation (1) is modified as outlined below for the forward-depth case:

z k i = j = 1 30 × ( d c Δ ) w c ( j ) × q k i ( Δ q t + j ) + v k Δ q t = ( d Δ w c ) , Δ w c = 30 × ( d c Δ ) (13)

The transitional probabilities (i.e., p ( x k i | x k 1 j )  or p ( q k i | q k 1 j ) ) 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 qtL to qtH. The likelihood PDF p ( z k | q k i ) in the HMM filter outlined in Table 1 is calculated based upon an assumed Gaussian measurement error as follows:

p ( z k | q k i ) = 1 2 π σ e [ q m ( d ) z k i 2 σ 2 ] (14)

where σ2 is the variance of the measurement noise. The HMM forward-depth estimated qtHMM cone tip values ( q ^ k F ) are calculated as follows:

q ^ k F = i = 1 N s w k | k i q k i (15)

3.2. qtHMM Test Bed Example

The previously outlined qtHMM 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 qt values below the cone tip (e.g., 54 unknown qt values for the case of a 10 cm2 cone tip (dc = 0.0357 m) and a sampling rate of 0.02 m (∆ = 0.02) for the forward-depth qtHMM Formulation. The same obviously applied for the backward-depth qtHMM Formulation, but in that case the 54 unknown qt values are above the cone tip. It should be noted, however, that with this qtHMM 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 qt value was not allowed to exceed 100. This is illustrated below for the case outlined in Figure 4(b) and Hsand/dc = 20.

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

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

Figure 5. Values of qm (black series) and qtHMM feed-forward qt (red series) estimates for varying qt maximum values specified for the case outlined in Figure 4(b) and Hsand/dc = 20. (a) qt maximum = 100, (b) qt maximum = 90, (c) qt maximum = 120 and (d) qt maximum = 250.

those of the true qt 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 qt 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 qt when Boulanger and DeJong used their inversion estimation algorithm for a situation where soil layers with a qt value of 12 were embedded in a uniform deposit with a qt 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 qtHMM Algorithm

IFM is incorporated into the qtHMM algorithm to address the poor test bed results. In this case, initial estimates for qt are derived utilizing IFM. Instead of attempting to estimate all the unknown qt 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 qt values are required to be estimated. In this process constant layer qt values and their corresponding depth extents are estimated for a maximum number of layers (specified by the user) within the next wc window.

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

As an example, assuming that a 10 cm2 cone is utilized for the sounding (with a diameter of 36 mm), then the extent of the wc 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 d1, d2, q2 and q3 have to be estimated, as the value of q1 is estimated with forward-depth formulation of the qtHMM-IFM algorithm where HMM filter transitional probabilities are taken into account. This is only 4 parameters as opposed to 54 qt values without the incorporation of IFM into qtHMM 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. qtHMM-IFM Test Bed Example

The performance of the qtHMM-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 qt values (light grey line in Figure 8(a)) and the resulting qm values were then calculated (black line in Figure 8(a)). Using the qtHMM-IFM algorithm the qt values were then estimated based on the qm values (black dotted line in Figure 8(a)). It shall be obvious that the algorithm performed well as the derived qt values closely matched the originally specified qt 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 qt between 84 and 96) has been lost, and that when focusing on qm 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 qtHMM-IFM algorithm is shown in Figure 9. In Figure 9 a soil profile was defined through qt values (grey line in Figure 9(a)) and the resulting qm values were then calculated

Figure 7. S Schematic outlining proposed IFM portion of the qtHMM-IFM algorithm parameters to be estimated. Parameters q1, q2, and q3 denote the cone bearing values for soil layers 1, 2 and 3, respectively. Parameters d1, d2, and d3 denote the depths of soil layers 1, 2 and 3, respectively.

(a)(b)

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

(a)(b)

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

(black line in Figure 9(a)). Using the qtHMM-IFM algorithm the qt values were then estimated based on the qm values (black dotted line in Figure 9(a)). It shall be obvious that the algorithm performed well as the derived qt values closely matched the originally specified qt 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 qm values thin soil layers are completely overlooked. The light grey ovals in Figure 9(a) clearly illustrate that by applying qm 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 qt value) can generate significantly different tip resistance readings (qm 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 qt value from recorded qm 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 qt estimation algorithm is referred to as qtHMM-IFM. Challenging test bed simulations have demonstrated that the qtHMM-IFM algorithm can derive accurate qt values from a qm 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 qtHMM-IFM algorithm on real data sets.

Cite this paper: Baziw, E. and Verbeek, G. (2021) Cone Bearing Estimation Utilizing a Hybrid HMM and IFM Smoother Filter Formulation. International Journal of Geosciences, 12, 1040-1054. doi: 10.4236/ijg.2021.1211055.
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

 
 
Top