Block-parametric approach to non-destructive control data analysis of railroad layered foundations

. The evaluation of the parameters of multi-layered foundations (ballast prisms, railroad basis, foundations of railway structures, etc.) plays an important role in ensuring the safe movement of trains. In the last decades new equipment arrived intended for non-destructive monitoring and estimating of the condition of such structures. These instruments can combine the high speed of scanning the diagnosed surface with high measurement accuracy of scanned parameters. But the use of similar equipment requires sophisticated methods for control data analysis including the building of adequate mathematical model of considered structure, efficient methods of discretization and forward problem solving, inverse problem analysis and based on such analysis surrogate model development for online data treatment, scanned by high-speed deflectometers or other scanning equipment. Presented paper devoted to the most important part of this chain of data processing procedures, namely the inverse problem analysis. The proposed approach is based on informational-probabilistic approach for inverse problem analysis, whose task is to obtain a posteriori probability density in the space of unknown parameters. The essence of the approach is the block-parametric approximation of the a priori probability density and likelihood function in the space of parameters and model data of the problem. The method allows estimating the parameters of the a priori distribution of unknown variable parameters, identifying and excluding outliers of the measured data from the created model, and constructing a posteriori estimation of the unknown parameters’ probability density with acceptable resolution. Proposed method can be used to create a new generation of equipment intended for non-destructive monitoring and estimating of the condition of the railroad basis and other roadbed structures. The appropriate software of such equipment based on the developed methods of data processing can be developed. The use of such equipment allows to operatively analyzing the state of individual areas of the railroad to decide on the need of repairing or replacing the railroad base or foundation of other elements of railroad infrastructure.


Introduction
Problem statement.Issues of non-destructive monitoring and control of the mechanical and geometric parameters of roadbed structures are one of the most important in the practice of operating a network of railroads, runways of airfields and other similar elements of transport infrastructure.The main source data supplied by modern devices for calculating such parameters are: -mechanical displacements of surface points of the testing pavement under the action of a given load (for example, car wheels, train weight etc); -values of the prone radiation reflected by the coating layers generated by the corresponding sources.The main characteristics of quality for such measuring instruments are accuracy, the ability to evaluate as many parameters of the tested coating as possible, and the high speed of taking relevant data.Devices of the first type, allowing to assess the degree of deformation of the test surface, are called deflectometers.Possible construction of such device equipped with rangefinder sensors for simultaneous scanning of the surface shape of the foundation under study is shown in the Fig. 1.Instruments of the second type include, first of all, Ground Penetrating Radars (GPR).The implementation of such devices raises the problem of creating methods that are more accurate for analyzing the diagnostic parameters of the roadbed structures to predict possible unallowable deformations and/or damage.

Review of resent researches and publications.
There is a variety of works devoted to roadbed state diagnostics.
For example, E.V. Nepomnyashchikh and K.A. Kirpichnikov in [1] set the tasks of diagnostics of the roadbed, namely: determining the parameters, that define the deformation properties of the roadbed; predicting the values of such parameters in time and their influence on the deformations occurrence of the roadbed construction; development of new principles, methods and technical means of diagnosing in order to detect "bottlenecks" in such structures.V.D. Petrenko, D.O.Yampolsky, I.O.Svyatko devoted his work [2] to the analysis of the existing modern and classical methods of numerical modeling of behavior of soil samples under the action of a static load.They proposed a new computational model, in which, in addition to the classical approach to the analysis of the stress state of the roadbed soil, the deformation state is taken into account.In the article [3] by scientists Luchko Y.Y. and Kravts, I. B., the condition of the roadbed on the railways of Ukraine is analyzed, the requirements for it and the methods of monitoring are considered.The attention is focused on the importance of monitoring the state of railway infrastructure facilities in the context of introducing high-speed traffic on the Ukrainian railways.The ground penetrating radar method (GPR) as one of the most promising methods for monitoring of the roadbed is proposed.The world and domestic experience in the use of GPR for non-destructive monitoring of the condition of the roadbed, ballast layer and artificial structures is also given [4,5].Its advantages are efficiency of work and low labor intensity.However, the methods proposed in this work provide information only about the geometrical parameters of the railroad bed base and do not give any information of its mechanical properties.
Regarding recent research, we can mention the work of M. Sysyn et al. [13], where the authors considered the issues of laboratory evaluation of railway ballast consolidation by the non-destructive testing.The relations of the residual stresses in the ballast layer to the elastic wave propagation were investigated in [14,15].Analysis of stressed state of sand-soil using ultrasound were performed in [16].
The purpose of the article.The purpose of this work is an attempt to develop such technique based on the information-probabilistic ("Bayesian") approach, which has shown its effectiveness in solving a wide range of problems in mathematical physics (there are relatively few works devoted to the application of this approach to solving inverse problems for layered foundations [6]; most of the works use the traditional deterministic ("Tikhonov's") approach).

Models of layered foundations
Parameters taking into account during the modeling of layered foundations can be divided into such main groups (Fig. 2).
1. Geometric characteristics of the layered packet.These parameters characterize the number and thickness of layers, as well as the equations of the interface between the layers and the upper free surface, usually determined using splines or other approximating functions.2. Mechanical properties of packet layers.These are, first of all, elastic constants (for example, shear modulus G and Poisson's ratio ν), as well as parameters characterizing plastic, viscoelastic, etc. properties of the material of the layer, as well as the degree of development of fatigue fractures.
3. Other physical parameters of the layers' material, such as the degree of compaction, moisture content, porosity, viscosity, etc.
4. Interface options between layers.The type of interaction (tight adhesion or the possibility of slippage of the layers in accordance, for example, with the Amonton-Coulomb law, etc.) and its parameters.
5. Characteristics of the soil subgrade on which the considered packet is located.Figure 2 shows a model of the Pasternak's subgrade, characterized by two parameters: the shear modulus Gs and the stiffness coefficient of the "spring" k, although models with a smaller (Winkler's subgrade) or larger (Kerr's subgrade) number of parameters can be used [7].
6.The parameters of the distributed load, the action of which a considered layered package during the test is subjected to.They characterize a spline or other approximating function q.
Model of layered foundation as a system with distributed parameters ( m (m , m , ) =  -vector of parameters), formulated in displacements, can be expressed by the operator equation There u u(x) = is a vector-function of packet points displacements, defined in composite domain 1 2 is a vector-operator, whose components (as usual in tensor analysis, summation over repeated indices is implied) are determined by the formula The components of the stress tensor σ and small strains tensor ε are related by the constituting equations in turn, the components of the small strain tensor are expressed in terms of the components of the displacement vector u by Cauchy-Green formulas Structure of displacement set 1 H [m] determines by restriction equations: boundary conditions ( ) (n is an outer normal vector to the surface of the packet, σni=σijnj are coordinates of the stress vector normal components) and conditions of interaction between the packet layers and with the soil subgrade.In the case of rigidly bonded ("glued") layers i and k on the contact The Winkler's model of the interaction between the packet and the soil subgrade is described by the relations (we neglect the friction between the K-th layer of the packet and the subgrade) where K S is a contact surface between the K-th layer of the packet and the soil subgrade, n is a unit normal vector to contact surface K S , τ is correspondent tangent vector, S (x) k is a subgrade local stiffness coefficient.In the case of the Pasternak subgrade model, instead of the first of equations ( 7), we have (8)   where S G (x) is a local shear modulus of the subgrade, ξ : 1 2 ( , ) ξ ξ are natural curvilinear (local) coordinates, defined on the contact surface Laplace operator in the specified coordinates.
In formulated models of a layered foundations, the vector of parameters m includes discretized "versions" of functions characterizing the differential operator (2) and the vector of the right-hand sides f of equation (1) (for example, the elastic characteristics G, ν of the packet layers included in the constitutive equations (3), as well as the quantities characterizing their plastic and rheological properties), domain Ω (for example, layer boundary equations for ik S and the boundary K S between the packet and the subgrade) and a system of constraints that determine the structure of the vector space H1 (surface load vector T, friction coefficients between layers ik µ , soil subgrade characteristics S k and S G ).
where ℜ(⋅): HM→Hv[m] is a forward operator, HM is model parameters vector space.For an "inverse problem" formulation one should takes into account an operator ℵ(⋅): Hv[m]→HD, that for each v(x) defines a "model data vector" of the problem where HD is a model data space, g(⋅): HM→HD is a "forward problem operator".As a model data, one should consider a radial displacement i u of certain roadbed points Ai, scanned by each used deflectometer (see Fig. 1) where lij are direction cosines of i-th deflectometer scanning ray, Kd is number of deflectometers used.The essence of inverse problem is "to invert" the forward problem operator g(⋅): for each (may be noisy) measured by the deflectometers vector d find the corresponding model parameters vector m, that satisfies (10).
Method for solving the formulated inverse problem is proposed in [6].This method is based on informationprobabilistic (Bayesian) approach to solving the inverse problem [7].According to this approach, solution of this problem is a posterior probability density of estimated model parameters vector m.For this approach, the "prior" (or "a priori") information about estimated parameters m, expressed as a correspondent prior probability density, plays a very important role.Without such information sometimes is impossible to obtain the qualitative inverse problem solution.

Block model for a problem with distributed parameters
We divide the vectors of parameters m and model data d into blocks (Fig3): We make the following assumptions regarding the introduced parameter and model data blocks.
I. Each block of parameters  , has an a priori probability distribution independent of the other blocks with a density ρ can be one of three types (densities are given up to a constant normalizing factor).
{ } { } Here ( ) ( ) are residual functions of the parameters and data blocks, respectively.
In the last expression, the covariance matrix for the data block takes into account the uncertainties of modelization and measurements.
{ } where the residual functions of the parameter and data blocks are determined by the formulas ( ) Similarly, the model data vector relative to the Boxcar block with the index 1, 2, , If the block is not of the Boxcar type, then any vector will be considered feasible with respect to it.
IV. Modelization uncertainties (d | m) The expression for the likelihood function of the i -th data block follows from (20):  .The parameter vector m will be considered feasible relative to the i -th block of model data if (m) 0 and not feasible if (m) 0 Let us obtain the expression for the posterior probability density ( 14 The residuals for the parameters m (i) S and data d (i) S blocks are expressed depending on the type of block by equations ( 14)- (19).
Common algorithm formulation.The general algorithm for analyzing inverse problems using the described block-parametric approach may look as follows.
1. Construction of blocks of observed data and unknown parameters of the generated model.This construction is performed taking into account the plausibility of the hypothesis about the probabilistic independence of the generated blocks of variables.If it is considered that the correlation between the variables should be taken into account, then these variables should be included in one block.
2. Identification of anomalous values of the observed data ("outliers").To do this, the checked data block is declared as a block of the Longtail type (the parameter blocks can be of the Boxcar type with the largest possible value of the standard deviation), after which a point in the parameter space is found corresponding to the maximum posterior hypothesis.For this parameter vector, the residual function is then analyzed for each of the points of the checking data block; the decision on whether to consider the corresponding value as "outlier" is made on the basis of statistical criteria similar to that of Grubbs, Irwin, and variety of their modern counterparts [8].The detected outlier points are either removed from consideration or replaced with interpolants for several neighboring points.
3. Rectification of the values of uncertainty parameters for the blocks of observed data and model parameters.Incorrect values of data uncertainty parameters are usually associated with incorrect determination of measurement accuracy or incorrect tuning of scanning instruments and can significantly affect the quality of the inverse problem solution.For checking, the data and parameter blocks are assigned the Boxcar type, and for the parameter blocks, the largest possible standard deviation values are set; the corresponding uncertainty parameters for the data blocks can vary.If the minimax criterion gives an unfeasible point as a solution (i.e., the posterior probability density is identically zero), this means that the root-mean-square deviation for a block of observed data is overestimated, and it is necessary to increase it.
In the case when the uncertainties of the measured data are relatively small, it makes sense to rectify the values of the uncertainties for the model parameters' blocks.To do this, in the above scheme, the standard deviations for data blocks are fixed, and for parameter blocks, on the contrary, they vary.The sought estimates for the uncertainties of the parameter blocks are those (minimum) for which the minimax criterion gives a feasible point in the parameter space.
4. At this stage, using the updated estimates of the parameter and data blocks' uncertainties, the appropriate level of parameterization is determined.For parameter and data blocks, as a rule, the Gaussian type is used, although, in fact, it is possible to use other types of uncertainties.To determine the appropriate level of parameterization, methods similar to those proposed by A.G. Ivakhnenko "methods of group accounting of arguments" [9] can be used.The main essence of this kind of methods is to build a certain hierarchy of parameterization levels and sequential verification of the model at different levels, starting from the lowest (containing the "coarsest" parameterization).The "stopping rule" at a certain level of the hierarchy is determined using the appropriate information criteria (Akaike criterion, Bayesian information criterion, etc. [10]).Some numerical results.Let us consider an elastic two-layer packet in the form of a square with side a in  (layer numbering starts from the lowest).The materials of the package layers are considered isotropic-elastic with shear moduli 1  We will approximate the unknown interface between the layers by the quadratic spline 1 1 ( ) was introduced to provide the possibility of setting an arbitrary value of the curve's The model data vector d of the problem also consists one block (1)   (1) (1) , where yj v  are dimensionless vertical displacements of the upper surface points of the packet in nodes of uniform division h Hereinafter, by "solution of the inverse problem" we mean a point in the parameter space M H , corresponding to the maximum posterior hypothesis.Consider the problem of refining the parameters of the standard deviation of the uncertain observable data.For this, we define a certain value of the parameter vector 1 T m {0.025,0.00625,0.00625, 0.025,0,0.025} ), for which we calculate the corresponding displacement field ( , ) x y v v   .The issues of constructing grids for discretizing the indicated "direct problem", as well as tuning the parameters of a Multigrid iterative algorithm for solving it, are considered in [11], [12].We form the vector of observable data obs d , taking the displacement values y v  at the corresponding nodes of the uniform division of the packet's upper surface ( N 9 d = ) and superimposing random Gaussian "noise" to them with the standard deviation parameter The low quality of the inverse problem solution (Fig5, (b)) is concerned with the (relatively) high value of the data "noise".However, such a model makes it possible to accurately estimate the "amplitude" of such noise (the abscissa of the intersection point of the graph in Fig5, (a) with the horizontal line R 1 = ).In the case when the observed data are determined quite accurately, it is also possible to evaluate the "swing range" of the parameter values.Consider a model similar to that studied when noise is not superimposed on the data (i.e.observed data is "exact"), the data and parameter blocks are of the Boxcar type, D 0.0001 σ = and the value M σ varies.The solution of the inverse problem is again determined using the minimax criterion.Consider the problem of determining the anomalous values ("outliers") of the observed data.In the problem under consideration, we introduce a perturbation of 0.001 in the "exact" value of the observed data component x  are presented, where in the nodes of the data block for the both models considered.For the first model, the data uncertainty is assumed to be of the Gaussian type with a unit correlation matrix (uncorrelated or "white" noise) and standard deviation D 0.001 σ = .In the second model, data uncertainty is considered as a Longtail type with the same standard deviation D 0.001 σ = .For the parameters' uncertainty of both models, the Boxcar type is used with a zero mean for all components and a standard deviation of M 0.1 σ = (in fact, the case of the absence of "regularization").From the graphs it can be seen that the "Gaussian" model "smears" the magnitude of the disturbance ("outlier") at neighboring points, while the "Longtail" model allows us to accurately localize the "outlier" and estimate its value.In addition, from a comparison of Fig7, 8, (b) it can be seen that the solution of the inverse problem in the "Longtail" model, in contrast to the "Gaussian" one, retains acceptable quality.

Prospects for further research
1. Clarify the effect of "internal noise" of the numerical method for solving the "direct problem" on the quality of the "inverse problem" solution within the framework of the proposed class of models.Since the main part of the computational work is precisely the approximate solution of the "direct problem", a reasonable choice of its discretization level will allow reaching an optimal price/quality ratio (as usual, "price" means the computational work and "quality" means the precision of the inverse problem solution).Such accounting for "internal noise" in the considered class of models can be carried out (assuming the normality of the "noise") by constructing the covariation matrix, taking into account both modelization and measurements uncertainties.2. Investigate how the quality of the inverse problem solution changes depending on the depth of the packet layer.Obviously, it will worsen with depth.It is of interest to find out for what degree of "data noise" the resolution of the proposed models will be sufficient to more or less satisfactorily determine the parameters of the deep layers of the packet, as well as the subgrade.3. Determine the capabilities of modern tools for multivariate function approximation (neural networks, etc.) to operatively obtain the inverse problem solution.This issue is of particular relevance for high-speed scanning devices that allow on-the-fly display of the diagnosed layered foundation parameters.

Fig. 1 .
Fig. 1.The scanning of the roadbed by deflectometers installed on the locomotive.

Fig 2 .
Fig 2. Layered Package Layout.Numbering of the groups of parameters see in the text.

Fig 3 .
Fig 3.The block model of the inverse problem.
) of the distribution of the model parameters vector m for the considered model I-IV.We divide the sets of indices of parameter I {1feasible (relative to the model I-IV) if it is feasible for each of parameters and model data blocks included in the model.Taking into account assumptions I-IV, we obtain

Fig4.
Fig4.Loading scheme, observed data and unknown parameters of the model problem (a); dependence of the magnitude of displacement y v  on the dimensionless coordinate x  (b).
the values of the function 1 ( ) S x  at the corresponding interpolation nodes.

Fig5.
Fig5.Dependence of the value of the minimax criterion R on the value of D σ (а); approximate solution of the inverse problem (bold black line) and "exact solution" (thin red line) (b).