Interval maximum likelihood inversion for well logging data and statistical properties of estimated parameters

During the conventional petrophysical interpretation fixed (predefined) zone parameters are applied for every interpretation zone (depth intervals). Their inclusion in the inversion process requires the extension of a likelihood function for the whole zone. This allows to define the extremum problem for fitting the parameter set to the full interval petrophysical model of the layers crossed by the well, both for the parameters associated with the depth points (e.g. porosity, saturations, rock matrix composition etc.) and for the zone parameters (e.g. formation water resistivity, cementation factor etc.). In this picture the parameters form a complex coupled and correlated system. Even the local parameters associated with different depths are coupled through the zone parameter change. In this paper, the statistical properties of the coupled parameter system are studied which fitted by the Interval Maximum Likelihood (ILM) method. The estimated values of the parameters are coupled through the likelihood function and this determines the correlation between them.


Introduction-the data system
In the process of well logging interpretation and inversion, the data are usually analyzed in predefined, disjoint depth intervals (interpretation intervals). For a given depth interval, the applied petrophysical model (petrophysical equations) and the associated predefined zone parameters (e.g. formation water resistivity, cementation exponent, shale and matrix characteristics etc.) are the same. In conventional inversion the estimation of rock parameters is performed on a depth by depth, usually as a solution to a weakly overdetermined problem, despite the high sensitivity of the model equation to zone parameters.
The focus of this study is mainly the zone parameter problem. It should be noted that interval inversion has been used by other authors to approximate the depth function of the local parameters (Dobróka, Szabó 2011;Dobróka et al. 2016).

3
If the zone parameters are included in the inversion, the fitting process for the different depth point data are coupled due to the common parameter(s). When the maximum likelihood (ML) method is used, the likelihood function should be extended over the entire interval data set (interval likelihood function).
The input data set for the inversion problem is the sampled well logging measurements M×N . Under the assumption of statistically independent measurement errors, the interval likelihood function (conditional density function of the measurement data on the parameters and petrophysical model) is as follows: where, N is the number of depth points in the interpretation interval, M is the number of logs, i is the depth index, j is the index of the measurement type (i ∈ (1..N);j ∈ (1..M)) , i ∈ ℜ K is the parameter vector (local parameters) associated with the ith depth point, which is the row vector of the parameter matrix P, K is the number of local parameters at a given depth point, p z is the vector of zone parameters. In case of zero mean, additive Gaussian noise, the interval likelihood function will be (Tarantola 1987;Szatmáry 2002): where w j is the weight associated with a given measurement type, inversely proportional to the variance of the measurement type w j = 2 2 j , σ is the normalization factor for the weights, and f j () is the model function (direct problem for the jth measurement) over the parameter space. L i () is the local likelihood function connected to the ith depth point. Taking the logarithm of interval likelihood to derive the functional to be minimized (Q) for the whole interpretation interval: The sum of squares for the interval is the sum of weighted squares (Q i ) for the local (single depth point) fitting problem. The parameters behave as a complex coupled system whose state is determined by the minimum value of Q.
For ease of handling, let us convert the elements of the measurement data system Y into the following vector form (y), ordered by depth: The measurement data vector defined above are elements of a smooth manifold (M y ). The parameter vectors ∈ ℜ N⋅K+K z to be determined are the elements of a manifold (M p ): K z is the number of zone parameters involved in the fitting (in this paper we consider the case of one zone parameter, but the statements can be generalized). The elements of the space M p are mapped by a ∶ ℜ N⋅K+K z → ℜ N⋅M vector-valued function into M f manifold, which is a subspace of M y (M f ⊂ M y ). Hence M f is parameterizable by elements of the parameter space. Given a proper log set, the relation M p and M f is bijective (which is also bijective locally per depth point). The model-dependent vector-valued function f() between the two manifolds: In the process of inversion (which means the minimizing of the functional Q), the current element of M y (the measurement data vector for the interval) is projected onto M f and, due to the homeomorphic relationship between M f and M p , also into the parameter space. This projection determines not only the parameter estimation (point estimate), but also the confidence interval of the parameters through the projection of the confidence interval associated with the measurement error distribution. This homeomorphism is the criterion for the necessary relationship between the measurement complex and the model space to ensure the unambiguity of the inversion. In the inversion, the inverse function between manifolds M f and M p is just approximated using Taylor series (mostly only up to the linear term).

Jacobian matrix
The Jacobian matrix (N⋅M×(N⋅K+Kz)) has fundamental importance in the relationship between the manifolds M f and M p in determining the inversion, i.e. the spatial projection to M f . Jacobian is also extended over the whole depth interval. Derivative according to local parameters define the matrix elements: Similarly the derivative according to zone parameters: The interval Jacobian matrix defined above provides a point-to-point local coordinate system (column vectors) on tangent space of M f to describe the projection of the error vectors. It plays a fundamental role in providing the projection that determines the linearized weighted least squares estimate. As is well known, at close to the optimum where the linear approximation is valid, the projection of Δy (error vector) on the parameter domain can be approximated (Tarantola 1987): where In the case of an additive error model, the estimated error vector belongs to the kernel of the operator on the right-hand side of Eq. 6.
The weights (w) in the quadratic form (Eq. 3) are also transformed to adapt for the "interval" formalism. Thus, for ease of handling can be written in the form of a diagonal weight matrix M⋅N×M⋅N with diagonal elements: The structure of the interval Jacobian matrix can be seen on Fig. 1. The coupling part related to the zone parameter is represented by the rightmost column of the matrix (vector J z ). If there are no common zone parameters in the inversion, then the fitting problem is separable for different depth points and the inversion can be performed in the conventional way. The matrix is also essential for inversion and later for estimating the covariance matrix of the parameters. Based on the known expression (Tarantola 1987;Szatmáry 2002) for the covariance matrix of the fitted parameters (applying Eq. 6): Fig. 1 The structure of interval Jacobian matrix. a Jacobian in the case of 1 zone parameter. b Jacobian without coupling: conventional approach In the transformation of above formula, the interval diagonal covariance matrix of the independent measurement data was used: The parameter σ in Eq. 9 can be estimated using the minimum of Q functional (denominator is the number of degrees of freedom): The following notation is introduced for the contraction of the J-matrix: The J-matrix and thus the matrix R for the interval can be decomposed into submatrices since the parts corresponding to the parameters of each depth point and the zone parameter are separated. The R-matrix of the local parameters associated with the oth depth point: The contribution of the zone parameter and the local parameters of the oth depth point submatrix: The R-matrix element defined by only the zone parameter: In the case of interval inversion, the structure of this matrix is also special, known as "arrowhead" matrix (Fig. 2).

Inverse of arrowhead matrix
Equation 8 shows that the covariance matrix of the parameters is determined by the inverse of the matrix R. It is also well known that the covariance matrix of the fitted function values is also based on this (Szatmáry 2002): The inverse of the matrix with the "arrowhead" structure can be expressed analytically (Jakovcevic et al 2015). Consider the following decomposition of matrix R (submatrices are defined in Eq. 12): where R 0 is the uncoupled R-type matrix of local parameters. This matrix also contains the previously defined quadratic local matrices (R 1 ,R 2 …) in a diagonal arrangement. The kth element of the vector R z is: This vector can also be decomposed into sub-vectors per depth point (Fig. 2). The inverse of the matrix R can be generated as the sum of two matrices, the first part can be written as the inverse of the uncoupled part of the matrix R, the second element can be defined in a dyadic form: The definition of vector u: The coefficient (ρ) of part 2: where I is the identity matrix. The correctness of the above inverse can be verified by the Sherman-Morrison theorem (Sherman, Morrison 1949;Saberi Najafi et al. 2014;Stanimirovic et al. 2019). This theorem can be applied to special matrices consisting of a diagonal and a dyadic part (Eq. 16). Thanks to the inverse given in the above analytical form, the covariance matrix (variances, covariances) of the fitted parameters can be studied directly (Fig. 3).

Sigma parameter
Besides the R matrix, the other important factor that determines the covariance matrix of the estimated parameters is σ: It also determines the probability distribution of the minimum value of Q (Szatmáry 2002): It can be estimated using the minimum value of Q, but its value is also affected by the number of degrees of freedom (DF) associated with the estimate. As the number of Fig. 3 Structure of arrowhead matrix inverse. The left part contains the R-inverse of conventional inversion −1 0 parameters decreases, the number of degrees of freedom increases, but the value of the squared deviation (Q) increases due to a potentially worse-fitting (or "smoother") model. For the analysis, let us take the reference case, a conventional depth-per-depth inversion to ensure the smallest weighted quadratic deviation (Q 0 ), where the parameter p z , which is later considered as the zone parameter, is also fitted per depth point (p z;i ) using Eq. 6. Then σ is the reference value (σ 0 ) for the whole interval: If the zone parameter is fitted as a constant (p z ) over the entire interval, the number of degrees of freedom increases significantly (NM-NK-1), so the variance of the parameters could be significantly reduced, but the Q min value increases.
The constant zone parameter is not optimal locally at the given depth point, so locally differs from the reference values (p zi ). This change also means a change in local parameters to "compensate", that is minimize Q i value with a correlated change in the new condition. Let's examine the change in Q i,min . The changes in Q i value can be approximated (by linear term): This inequality is persisting during the change Q i ≥ Q 0,i . The ΔQ i increment is expressed by Taylor-expansion around the reference ( Q 0i ): Plotting the quadratic form above in matrix form: This defines a quadratic hyper-surface as a function of the parameter change vector (ΔQ i (ΔP,Δp z;i )), which takes a value of zero at the origin. Again, a new minimum search, the direction of the smallest Q i changes can be determined: From here, the equation of the parameter change: Substituting this back into the equation for ΔQ i : Summarizing the increments for all depth points to get the change for the whole interval: Taking this effect into account, the value of σ is: By including the zone parameter into the inversion, the improvement in the degree of freedom is not significant.

Parameter estimation by interval likelihood
Parameter estimation can be done for the interval inversion using Eq. 6. We often use this linearized formula to find the minimum value of Q, usually by iteration with R matrix. For interval inversion, the form prescribed in Eq. (16) is preferable. The parameter estimator operator is then split into two parts (except for the zone parameter).: The coefficient of the correction term (ρ) is decreasing with N, this also indicates that the result of the conventional inversion can be a good starting value for the interval inversion. Because of the large number of parameters, a good choice of initial value is particularly important. If the overdetermination allows, the zone parameter can be estimated per depth point from which the initial value can be estimated. From the linear approximation (Eq. 28), we can determine the initial value for p z : Minimizing the increment gives this initial value estimator: The above equations can also be used to study the conventional interpretation practice: i.e., using prior information and not inversion to arrive at the zone parameter (p a ). (27b) This approximate formula shows how the goodness of fit deteriorates if the zone parameter is not from the inversion, and thus may contain some bias. (This is essentially a manifestation of Steiner's theorem.) The formula shows that the mismatch increases with N.

Statistical properties of the estimated parameters
From Eq. (30) it follows that the parameter covariance matrix also can be separated two parts. The first term contains the local covariance matrices for the local parameters and the second term containing the effects due to coupling. The latter contains the covariances between the parameters of different depth points and the covariances between the zone parameter and the local parameters (off diagonal elements). The magnitude of this correlation is greater the more similar the rock at the two depth points (under similar conditions, the local parameters "compensate" with similar way, which causes the correlation). In other words, if the local rock parameters are similar for two depth points of an interpretation zone, then a small change in the zone parameter (generated by the global minimum search), causes similar changes in the local parameters to reduce the local sum of squares.
The second part (in Eq. 3) also modifies the diagonal elements which determines the estimated value of the parameter variances.
The coefficient ρ, which determines the importance of term 2, inversely proportional to the number of depth points (N) of the interval. This means that as N increases, the term 2 tends to zero, i.e. for large values of N the coupling effects tends to zero. Asymptotically, the different depth point parameters become uncorrelated (in case of Gaussian noise: asymptotic freedom).

Variances of estimated parameters
The variances of the parameters are the corresponding diagonal elements of the covariance matrix (Eq. 8). The contribution of the diagonal elements in matrix 2 represents a "correction" that vanishes as N increases. In the case of the zone parameter, only the part 2 contributes to the variance: Thus, of course, for large N, this variance also converges to zero (consistent estimate). The variances of the local parameters, on the other hand, do not tend to zero, but to the corresponding diagonal element in the first term (Fig. 4).
Estimation of the variance of the kth local parameter (at depth point D): The interval maximum likelihood method can therefore be used to increase the efficiency of the estimation, especially for the zone parameter, which improves steadily as the interval increases.
Equation (24) shows that the covariance between zone and local parameters has opposite sign to the covariance between the local parameters. As is can be seen on Fig. 5, by the increment of the number of depth points, the covariance structure also tends to the values obtained with the conventional method (asymptotic freedom in the case of Gaussian noise).

Conclusion
If the likelihood function extended over the entire data set of a depth interval, then the zone parameters can be included in the inversion. In terms of inversion, this links the parameter changes during the minimum search, which is reflected in the estimated covariance structure. In this study the covariance matrix has been generated in a form that separates the covariance structure associated with conventional inversion and the correction from parameter coupling. The coefficient of the correction part decreases with the number of depth points in the interpretation interval. Therefore, the correlation between the zone parameter and local parameter or local parameters in different depth point is vanishing with the number of depth point (asymptotic freedom). The variance of zone parameter also decreases by N tending to zero with large point number, while the variance of local parameters approaching the value of the conventional inversion. So, the paper has shown that the zone parameter can be efficiently estimated using ILM, which is essential when considering the sensitivity of petrophysical equations. The absolute value of correlation coefficient as a function of interval depth point number. The correlation between the zone parameter and local parameter is higher than the correlation between two local parameters (which are not belong the same depth point)